Effective mass path integral simulations of quasiparticles in condensed phases

Electronic excitations in semiconducting materials form the foundation of many areas of materials and energy sciences, including solar energy harvesting and conversion and nanoelectronics. It is often advantageous to describe processes involving such excitations within the language of quasiparticles, e.g., electrons and holes with effective masses or excitons.^1–3 Due to the complexity of these descriptions, the theory and simulation of quasiparticles are often limited to coarse-grained and continuum approaches or modeling small, highly symmetric systems in quantum mechanical detail.^4–11 While these approaches are responsible for important advances in our understanding of exciton physics and nanotechnology, it is difficult for the existing methodologies to describe the sources of disorder (defects), and large system sizes are often needed to properly model the effects of charge carriers on their surroundings and vice versa.

One promising approach for simulating quantum particles in complex environments uses the path integral (PI) representation of quantum mechanics in which a quantum particle, such as an electron, can be represented as a classical ring polymer. Given pseudopotentials to describe the interactions between a quantum particle and its (often classical) environment, one can perform molecular dynamics (MD) simulations of large systems for long times. However, it is difficult for straightforward PIMD simulations to describe phenomena such as quantum confinement that manifest as a result of anisotropic electronic band structures, in addition to quasiparticles such as holes. In this work, we describe an approach that uses effective mass theory^12–17 to incorporate these aspects of anisotropic electronic band structures of materials into the path integral representation of quantum mechanics^18–22 in order to model quantum charge carriers and their excited states (e.g., excitons) in complex, atomistic environments. Such an approach can be considered as an intermediate-scale coarse-graining, which incorporates some effects of the interactions into the effective masses, leaving simpler interactions to be described by pseudopotentials. We demonstrate the utility of this effective mass path integral (EMPI) approach by modeling an exciton in crystalline potassium chloride and electron trapping in a defective monolayer of molybdenum disulfide (MoS₂).

We consider a quantum particle described by a Hamiltonian $H=T+V$⁠, which consists of a kinetic term, $T$⁠, and a potential term describing its interactions with the environment, $V$⁠. The partition function for this particle can be written as

By applying the (symmetric) Trotter factorization,^18–20,23 we can write $Z$ in a form that physically corresponds to a discretization of the (cyclic) quantum path of the particle in imaginary (or Euclidean) time,

In practice, we use a finite number of discretizations, P, and the partition function and equilibrium ensemble averages are exact in the limit P → ∞, corresponding to a continuous path. The potential term, $V(r)$⁠, can be readily evaluated, and so we focus on rewriting the kinetic part of the partition function. This arises from matrix elements of the form

which connects each discrete step in the imaginary time path of the particle, e.g., step i to j.

We now work within effective mass theory (EMT) to include some aspects of the electronic structure in our model through the above matrix elements. EMT prescribes an effective mass m^* to charge carriers, which reflects the modification of their masses, from that of a free electron, due to interactions with other charge carriers and the static nuclei.^12–16 In order to model quantum particles in a classical bath using the path integral isomorphism, we substitute the free masses of the quantum particles with those determined from EMT to include the influence of the electronic response to the environment in effective classical models. For highly symmetric materials, isotropic parabolic fits of the band structure may be sufficient such that a scalar effective mass m^* can be assigned to the particle.^24,25 However, lower symmetry materials require that the 3 × 3 mass-tensor

can be computed, leading to an inverse effective mass-tensorm⁻¹ for each quantum particle that can, in principle, be anisotropic. Note that the effective mass-tensor is symmetric, $mαγ−1=mγα−1$⁠. We also drop the star for notational clarity.

In the context of EMT, the kinetic energy of the quantum particle is now

which reduces to $T=p2/2m$ in the limit of an isotropic, diagonal mass matrix. We can then follow the typical evaluation of the partition function in the path integral isomorphism, but now with Eq. (5) for the quantum kinetic energy. The desired matrix element can be readily evaluated through Gaussian integration²⁶ to yield

where m_αγ is the α, γ element of the inverse of the matrix m⁻¹. With this expression for the matrix elements, the partition function is given by

where $HP$ is the isomorphic Hamiltonian of a classical ring polymer with harmonic bonds between neighboring beads of the polymer. This isomorphic Hamiltonian is

where

and the matrix elements m_αγ can be determined by inverting m⁻¹. The harmonic bonds between neighboring beads of the ring polymer generally are not spherically symmetric but involve different spring constants, κ_αγ, along each direction, as well as coupling between the displacements in the Cartesian components, as demonstrated below.

In the limit of a diagonal m⁻¹, the ring polymer Hamiltonian becomes

In this formulation, the harmonic springs between neighboring beads of the ring polymer do not have spatially isotropic spring constants but are instead given by κ_αα = P²m_αα/β²ℏ², where α refers to a spatial coordinate. This is particularly important in systems with reduced dimensionality. For example, m_zz ≫ m_xx ≈ m_yy in monolayer MoS₂, as discussed below, resulting in ring polymers that are confined essentially to two dimensions.

To illustrate a specific case where off-diagonal coupling is present, we explicitly consider an effective mass-tensor with $mxz−1=mzx−1=myz−1=mzy−1=0$⁠; all other elements are non-zero. In this case,

with the spring constants

and

where $μxy=1−mxxmyy/mxy2$⁠. The analogous Hamiltonian when $mxz−1$ or $myz−1$ is the only non-zero off-diagonal element can be readily obtained by permuting the relevant indices in Eq. (10). Inclusion of a single off-diagonal element (xy) in the effective mass-tensor leads to a coupling between the x- and y-directions. The presence of the off-diagonal coupling additionally renormalizes the effective xx and yy spring constants by a factor of μ_xy. This may be expected from conservation of energy. Because some of the quantum kinetic energy is transferred into the coupling between the x- and y-directions, the xx and yy components of the kinetic energy must be correspondingly reduced to account for this energy transfer into the additional degree of freedom. Hence, the contribution from the off-diagonal coupling is opposite in sign to the diagonal coupling, and the prefactor is twice the magnitude of that for a single diagonal term (half from xx and half from yy).

We now demonstrate the utility of our approach in the context of simulating an exciton in an alkali-halide crystal. We model an electron–hole pair at a constant temperature of T = 300 K, and the electron and hole ring polymer each have P = 1024 beads. These MD simulations were performed using the LAMMPS software package²⁷ with Nosé–Hoover chains to maintain a constant temperature²⁸ and a Parrinello–Rahman barostat to maintain zero pressure,^29–31 allowing the crystal to relax to the presence of the exciton. The ring polymers were massively thermostatted to ensure proper sampling.^32,33 We use the Tosi–Fumi model for KCl.³⁴ Electrostatic interactions between charges of the same sign are described with a standard Coulomb potential, while those between charges of opposite sign are described by a Shaw pseudopotential,³⁵ with short-ranged cutoffs of 1.96 Å, 1.75 Å, and 1.69 Å for electron–K⁺, hole–Cl⁻, and electron–hole interactions. The electron–K⁺ cutoff is that used by Parrinello and Rahman in their seminal study of F-centers in KCl,³⁶ the hole–Cl⁻ cutoff corresponds to the Cl van der Waals radius, and the electron–hole cutoff is crudely chosen to yield the bandgap in the single bead limit; the bandgap is the energy difference between infinite separation and a perfectly overlapping electron and hole. Further refinement of the latter cutoff can be performed to match the exciton binding energy, for example, but we reserve this for future work, and note that the cutoff used here yields reasonable predictions. Long-ranged electrostatic interactions were evaluated using the particle–particle–particle mesh Ewald method.³⁷ 

Band structure calculations were performed using the GPAW software package,³⁸ in combination with the atomic simulation environment (ASE),³⁹ and employed the PBE density functional approximation⁴⁰ with a plane-wave cutoff of 1200 eV and a 12 × 12 × 12 k-point mesh. Although the bandgap is not properly described at this level of theory,^41,42 the curvature of the bands near the gap is likely adequate, and we expect the effective mass for these materials to be somewhat insensitive to the choice of semi-local functional. Note that high throughput calculations of more complex materials have also used density functional theory (DFT) effective masses successfully in order to correlate transport and other properties from the predicted band structures.^43,44 For highly accurate band structures and effective masses, approaches beyond semi-local DFT are required.^45–47 Effective masses herein were determined using the effective mass calculator (EMC) program using a five-point stencil to evaluate the second derivatives.⁴⁸ We find that m⁻¹ is approximately diagonal such that the masses of the electron and hole are approximately isotropic, m_αγ ≈ δ_αγm^*, and equal to $me*=0.45me$ and $mh*=5.2me$⁠, respectively, where m_e is the bare mass of an electron.

A single exciton introduced into an otherwise perfect alkali-halide crystal can self-trap and create lattice defects.^7,49–54 This self-trapping results in a structure resembling a closely separated F-center–H-center pair, where the latter corresponds to a hole bridging two anions, Cl$2−$⁠, and similar states, e.g., Cl$32−$⁠. The self-trapped exciton is not spherical, as one might expect using continuum theories of excitons in condensed phases. Instead, the exciton is expected to be dipolar, with the electron and hole separated by some average distance R_eh > 0 even in the bound, excitonic state.^7,51,52

The formation of this self-trapped exciton state in solid KCl is illustrated by the snapshot in Fig. 1(a). The electron–hole pair distribution function in Fig. 1(b) demonstrates that the dipolar exciton consists of an electron and hole separated by R_eh ≈ 2 Å. This is further supported by the electron–hole spatial distribution function shown in the inset of Fig. 1(b). This is in very good agreement with previous interpretations of experiments and detailed theoretical calculations that also predict a dipolar exciton with R_eh ≈ 2 Å.⁵¹ 

While the light electron is delocalized over many ions (but localized with respect to a free electron with the same effective mass), the heavy hole is highly localized and bridges Cl⁻ ions, as illustrated by the snapshot in Fig. 1(c). The Cl–Cl pair distribution functions, g(r), of the KCl crystal in the absence and presence of the exciton are compared in Fig. 1(d), averaged over all Cl⁻ in the system. In the presence of the exciton, a peak at close Cl–Cl distances appears in g(r), consistent with the formation of Cl$2−$-like structures predicted in more detailed quantum calculations and experiments.^7,51–53 We additionally note that the hole–anion interaction potential we employ is spherically symmetric. Including directionality into the hole–anion interactions, e.g., through the use of multisite ion models,⁵⁵ for example, may lead to even better descriptions of H-centers.

The EMPI formalism can also describe the effects of reduced dimensionality. For example, monolayer transition metal dichalcogenides, such as MoS₂, are two-dimensional materials that exhibit quantum confinement in the two-dimensional plane of the lattice.^5,6,56–58 In this case, the effective masses suggest that the quasiparticles are essentially confined to the xy-plane, m_xx = m_yy ≈ 0.562m_e ≪ m_zz ≈ 960m_e. The large z-component of the effective mass-tensor manifests as a high spring constant κ_zz that confines the quasiparticle significantly in the z-direction while allowing the particle to spread in the two-dimensional (xy) plane, as illustrated by the snapshots in Figs. 2(a) and 2(b). Note that the out-of-plane effective mass (m_zz) will go to infinity as the width of the semiconductor tends to zero, i.e., a purely two-dimensional system.

We demonstrate the utility of EMPIMD simulations for systems with reduced dimensionality by studying the binding of an excess electron in monolayer MoS₂ to a sulfur vacancy, which carries an effective positive charge. Monolayer MoS₂ is modeled using the Stillinger–Weber potential designed to capture the structure and vibrational properties of the monolayer,^59,60 which readily enables modeling of defects. The interactions between the electron and the S atoms are modeled using a Coulomb potential, and those between the electron and the Mo atoms are modeled using a Shaw pseudopotential³⁵ with a short-range cutoff of 0.18 Å. The partial charges on the Mo and S atoms, used for interactions with the electron, are those from the work of Sresht et al.⁶¹ Simulations are performed using LAMMPS²⁷ with appropriate Nosé–Hoover thermostatting to maintain a constant temperature of 300 K. The anisotropic spring constants of the electron ring polymer are incorporated using PLUMED.⁶² The electron–S-vacancy potential of mean force is calculated using umbrella sampling⁶³ combined with UWHAM,^64,65 where we biased the two-dimensional (xy) distance between the electron and the location of the vacancy using harmonic potentials in PLUMED.⁶² 

Spatially resolved photoluminescence (PL) spectroscopy has discovered that excitonic hotspots appear at the location of sulfur vacancies in monolayer MoS₂.^66–69 The higher intensity peaks in the vicinity of sulfur vacancies suggest that excess electrons in doped MoS₂ are bound to these vacancies. Increasing the concentration of sulfur vacancies results in the appearance of a new, lower energy peak in the PL spectra, further suggesting the validity of this interpretation. Our EMPIMD simulations further support this picture of strong electron–S-vacancy interactions in monolayer MoS₂, leading to the formation of a trap state.

To quantify the interactions between the electron and a sulfur vacancy, the free energy as a function of electron–vacancy distance is shown in Fig. 2(c). We find a binding free energy of ∼30k_BT at 300 K, in agreement with the range of energies predicted by kinetic modeling of spectroscopic measurements.⁷⁰ This binding free energy is also in good agreement with the energy difference between the trap state and the conduction band predicted by density functional theory (DFT) calculations.⁷¹ 

The scaling of the free energy with distance is in good agreement with that expected for Coulomb interactions within an isolated two-dimensional (2D) semiconductor,^6,72

where d is the electron–vacancy in-plane distance, H₀ is the Struve function, Y₀ is the Bessel function of the second kind, and d₀ is a screening length. We compare V_2D(d) with the simulated free energies in Fig. 2(c) for 30 Å ≤ d₀ ≤ 40 Å, as estimated in Ref. 6. We also show the three-dimensional Coulomb potential V_3D(d) = −A/d, where A was determined by fitting the simulation data from 2.5 Å to 20 Å. The simulated free energies follow the two-dimensional scaling, V_2D, for all but the smallest electron–vacancy distances and not V_3D, as expected for 2D materials. Note, however, that this only verifies that the distance scaling of the free energy is consistent with reduced-dimensionality effects on electrostatics. Quantitative agreement between simulation results and theoretical predictions⁶ that account for both the effective charge state of the vacancy and electronic screening within MoS₂ may require inclusion of spatially varying dielectric constants within the semiconductor in our simulations, as well as a reparameterization of the pseudopotentials for this environment.

Our approach additionally enables the characterization of the effects of the sulfur vacancy on the electron. For example, examination of the imaginary time root mean squared displacement,^73,74 $R(τ)=rxy(τ)−rxy(0)21/2$⁠, where r_xy indicates that distances in-plane were considered in the calculation. Comparison of $R(τ)$ for the trapped electron with that expected for a free electron in MoS₂, $Rfree(τ)=2λ2(t/βℏ)(1−t/βℏ)1/2$⁠, where λ² = βℏ²/m and m = m_xx = m_yy, indicates that binding to the sulfur vacancy traps the excess electron, as shown in Fig. 2(d). The value of $R(βℏ/2)≈5$ Å yields an estimate for the effective size of the electron in good agreement with that determined for the trap state from DFT calculations.⁷¹ Moreover, writing $R(τ)$ in the basis of electronic eigenstates,^73,74

where we have considered the 2D case assuming x and y are equivalent by symmetry, Z is the partition function, $xnm=nxm$⁠, and E_n is the energy of eigenstate n, further suggests that the electron is in a deep trap state. In order for $R(τ)$ to be independent of τ, as is the case for τ ≫ 0, Eq. (15) must be dominated by a single, localized eigenstate with a significant gap to the first excited state. This ground state that dominates the behavior of the vacancy-bound electron is the trap state. In contrast, the spatially extended states sampled by the free electron correspond to the thermally accessible excited electronic states. As shown in Fig. 2(d), an “unbound” electron, here obtained by simulating an electron in a defect-free monolayer, yields $R(τ)$ consistent with that for a free effective mass electron, with a slight renormalization due to dynamic electron–phonon coupling. Fitting to $Rfree(τ)$ yields a heavier in-plane effective mass of m ≈ 0.7m_e, as expected for electronic coupling to phononic distortions of the nuclear sites.³ These results highlight the utility of EMPIMD simulations in quantifying the thermodynamics of quasiparticle interactions in anisotropic materials.

In this Communication, we have presented a formulation of effective mass path integral simulations. This approach incorporates effective mass theory into the path integral description of electrons and holes for their simulation in condensed phases, extending previous PI-based methods using scalar effective masses to materials with anisotropic electronic structures. We expect that this approach will find wide use in a variety of applications, including the simulation of exciton and charge carrier structure and dynamics in materials.

We note that our implementation of EMPIMD simulations described above has not exploited the many significant advances made in the context of path integral simulations in recent years. However, the EMPIMD method is readily amenable to such approaches, including novel integration schemes,^75–78 thermostats,^79–81 and nearly all other developments in the context of path integral and ring-polymer techniques.^82–88 Of additional importance for the future of EMPI models is the development of accurate electron–environment and hole–environment pseudopotentials, beyond the simplistic, spherically symmetric charge–charge pseudopotentials used here,^89,90 as well as the inclusion of modified electrostatics due to dielectric screening in low-dimensional materials.^4,6 These effective interaction potentials will impact the accuracy of the predictions made by EMPI approaches. Finally, we note that this framework can also be readily used within the ring-polymer MD approximation and similar approaches^23,91–93 to model the quantum dynamics of quasiparticles in complex environments, and future work will focus on extending these approaches to EMPI simulations.