Diffusion quantum Monte Carlo study of excitonic complexes in two-dimensional transition-metal dichalcogenides

Excitonic effects play a particularly important role in the optoelectronic behavior of two-dimensional semiconductors. To facilitate the interpretation of experimental photoabsorption and photoluminescence spectra we provide (i) statistically exact diffusion quantum Monte Carlo binding-energy data for a Mott-Wannier model of (donor/acceptor-bound) excitons, trions, and biexcitons in two-dimensional semiconductors in which charges interact via the Keldysh potential, (ii) contact pair-distribution functions to allow a perturbative description of contact interactions between charge carriers, and (iii) an analysis and classification of the different types of bright trion and biexciton that can be seen in single-layer molybdenum and tungsten dichalcogenides. We investigate the stability of biexcitons in which two charge carriers are indistinguishable, finding that they are only bound when the indistinguishable particles are several times heavier than the distinguishable ones. Donor/acceptor-bound biexcitons have similar binding energies to the experimentally measured biexciton binding energies. We predict the relative positions of all stable free and bound excitonic complexes of distinguishable charge carriers in the photoluminescence spectra of WSe$_2$ and MoSe$_2$.


I. INTRODUCTION
The last decade has witnessed a remarkable surge of interest in the properties of truly two-dimensional (2D), atomically thin semiconductors. These include monolayer transition-metal dichalcogenides (TMDCs) such as MoS 2 , MoSe 2 , WS 2 , and WSe 2 , which acquire a directgap character in hexagonal monolayer form. [1][2][3][4] The di-rect gap and strong optical absorption of TMDCs suggest a range of potential optoelectronic applications, e.g., in photodetectors, photovoltaics, and light-emitting diodes. A particularly interesting aspect of monolayer TMDCs is the strong excitonic effects present in their photoabsorption and photoluminescence spectra, [5][6][7] including nonhydrogenic Rydberg spectra 8,9 and lines ascribed to trions (charged excitons) [10][11][12] and biexcitons (bound pairs of excitons). [13][14][15][16] The nonhydrogenic nature of the excitonic energy spectrum is due to lateral polarization effects in 2D crystals, which modify the form of the Coulomb interaction between charge carriers. Mott-Wannier models of 2D trions and biexcitons have been studied using quantum Monte Carlo (QMC) methods, [17][18][19][20][21][22][23][24][25] , variational methods [26][27][28] , and hyperspherical harmonics approaches, 29 and interpolation formulas linking the 2D-screened and 1/r Coulomb interaction regimes have been proposed. Here we extend these studies to provide numerically exact binding-energy data for all nonlocal screening strengths, including an analysis of limiting behavior, and we classify the types of trion and biexciton that can be observed in different TMDCs. We also investigate donor-and acceptor-bound chargecarrier complexes in TMDCs, such as donor-bound biexcitons and quintons, which have not to our knowledge been studied before.
The rest of the article is structured as follows. In Sec. II we describe the band structures of molybdenum and tungsten dichalcogenides and analyze the nature of the trions and biexcitons in these materials; furthermore, we perform a group theoretical analysis of exciton properties. In Sec. III we explain the Keldysh form of the screened Coulomb interaction between charges in 2D semiconductors, describe the ways in which chargecarrier complexes are expected to dissociate and recombine, and explain the importance of the contact pair distribution function (PDF). In Sec. IV we describe our computational methodology for solving the Mott-Wannier model of charge-carrier complexes. We present our numerical results for the binding energies and PDFs of the different complexes in Sec. V. Finally, we draw our conclusions in Sec. VI.

A. Classification of trions and biexcitons
In monolayer molybdenum and tungsten dichalcogenides the conduction-band minimum and valence-band maximum occur at the K and K points of the hexagonal Brillouin zone. Spin-orbit coupling induces a significant splitting of both the valence band and the conduction band at K and K . In molybdenum diselenides, the valence-band maximum has the same spin as the conduction-band minimum within each valley, while in tungsten dichalcogenides such states have opposite spins. 3 Figure 1(a) presents examples of the ways in which biexcitons can be formed in molybdenum and tungsten dichalcogenides. The spin-splitting of the valence band (0.15-0.5 eV) is sufficiently large that no holes in the lower spin-split valence band are expected at room temperature; however, the spin-splitting of the conduction band (∆ = 3-50 meV) is small enough that electrons can be found in the upper spin-split conduction band at room temperature. 3 An exciton, biexciton or trion is said to be either dark or semidark when the recombination of an electron and hole is forbidden by spin and momentum conservation; otherwise the complex is said to be bright. Semidark complexes are those in which recombination can in fact take place due to intervalley scattering with an accompanying energy shift. The precise photon energies depend on whether the electrons occupy the higher-or lowerenergy spin-split bands in the initial and final states. Furthermore, the intensity of a spectral line depends on the thermal occupancy of the initial state. Figures 1(b) and 1(c) present a classification of biexcitons in molybdenum and tungsten dichalcogenides with respect to the recombination energy and the intensity of the emitted photons. This intensity has the following temperature dependence: for no electrons in the upper spin-split conduction band e −∆ /(kBT ) for one electron in the upper spin-split conduction band e −2∆ /(kBT ) for two electrons in the upper spin-split conduction band , where ∆ is the spin-orbit-induced splitting of the conduction band, k B is Boltzmann's constant, and T is the temperature. A similar classification can be made for trions: see Fig. 2. In a photoluminescence experiment, we expect to see energies attributed to different kinds of biexcitons and trions and emission lines of varying intensity, as explained in Sec. III C.
The opposite spin-splittings of the conduction and valence bands in tungsten dichalcogenides results in the ground-state trions and biexcitons being dark, with the two electrons residing in opposite valleys. These dark complexes are coupled through an intervalley electronelectron scattering to their excited bright counterparts with both electrons residing in the upper spin-split conduction band. This coupling gives a finite oscillator strength to the dark ground states that is proportional to [µ bd /(2∆ )] 2 , where µ bd is the coupling matrix element between dark and bright states. As a result, the We only show the spin-splitting of the conduction band; the spinsplitting of the valence band is much larger, so that no holes in the lower spin-split valence band are expected at room temperature. 3 (b) and (c) Classification of biexciton recombination processes in molybdenum and tungsten dichalcogenides, respectively. ∆ is the band gap, while ∆ is the spin-splitting of the conduction band. E δ ≡ EXX − EX is the difference between the total energies EXX and EX of a biexciton and an exciton. ω indicates the photon energies at which peaks in photoluminescence spectra are expected to appear. XX k 1 σ 1 k 2 σ 2 k 3 σ 3 k 4 σ 4 denotes a biexciton consisting of conduction-band electrons in valleys k1 and k2 with spins σ1 and σ2 and valence-band holes in valleys k3 and k4 with spins σ3 and σ4.   Fig. 1, but for negative trions in molybdenum and tungsten dichalcogenides. E δ ≡ E X − is the total energy E X − of a negative trion. T k 1 σ 1 k 2 σ 2 k 3 σ 3 denotes a trion consisting of conduction-band electrons in valleys k1 and k2 with spins σ1 and σ2 and a valence-band hole in valley k3 with spin σ3. expected photoluminescence spectrum contains two additional lines resulting from the recombination of these "semidark" trions and biexcitons, at an energy shifted downwards by 2∆ relative to the bright complexes, and having a temperature-independent intensity.

B. Group theoretical analysis of excitons
Exciton wave functions can be classified according to the irreducible representation (irrep) of the point-group symmetry of the TMDC crystal, D 3h . As the states in the two valleys are degenerate, one can treat the two valleys simultaneously by using the extended group D 3h = D 3h + tD 3h + t 2 D 3h , where t denotes translation by a lattice vector. The character table of the extended group is given in Table IX. The total exciton wave function X is given in general by the product of three components: the spatial envelope function Φ, the Bloch or lattice wave functions of the electron and hole U k , and the spin part χ: The representations of the wave functions by irreps consist of the direct product of the individual irreps corresponding to the three components: The tightly bound ground-state excitons are characterized by a maximally symmetrized envelope function corresponding to the identity irrep Γ Φ = A + 1 . Therefore the representations of the exciton states are determined by the irreps of the lattice and spin parts.
The conduction-and valence-band Bloch states transform according to the 2D irreps E 1 and E 2 , respectively. Using the product table, Table X, the lattice part of the exciton wave function transforms as where the 2D irrep E + corresponds to intravalley excitons in the K and K valleys, and E + 3 corresponds to intervalley excitons, which are dark due to momentum conservation. In the following, we will consider the E + intravalley excitons only.
The exciton spin part consists of two spin-1 /2 particles corresponding to the spinor 2D irrep D1 /2 . The direct product of the two spinors can be decomposed into the crystal point group irreps as Hence the total exciton representation is given by The E + irrep corresponds to the vector representation, and therefore the two E + irreps correspond to excitons coupled to in-plane polarized light. The z coordinate transforms as the A − 2 irrep, and therefore the A − 2 exciton is coupled to out-of-plane polarized light, which involves a spin-flip process in recombination. 30 In the case of tungsten dichalcogenides, the A − 2 exciton is the ground-state exciton, and results in photon emission at an energy that is lower than the excited bright exciton by the spin-orbit splitting of the conduction band ∆ . The A − 1 and E − excitons are not coupled to light. A summary of the classification of exciton states is given in Table I using a notation similar to that used in Figs. 1 and 2. Finally, we note that the spin-flip transition resulting in the emission of out-of-plane polarized light corresponding to the A − 2 exciton is also relevant for ground-state trions and biexcitons in tungsten dichalcogenides, resulting in trion or biexciton emission at a photon energy shifted downwards relative to the excited bright states by ∆ . I. Classification of exciton states into irreps of D 3h and the polarization ( and z for in-plane and out-of-plane, respectively) of the electric field to which the excitons are coupled.

A. Screened Coulomb interaction between charge carriers
We model the charge carriers in a 2D semiconductor using a Mott-Wannier model, in which small numbers of quasielectrons and quasiholes are treated within the band effective mass approximation and interact via an appropriately screened Coulomb interaction. The band effective masses for different 2D semiconductors are assumed to be 2D-isotropic, and are discussed in Sec. V B 2. However, unlike quasi-2D electron(-hole) systems in GaAs/InAs heterostructures, the form of the Coulomb interaction is profoundly affected by the 2D nature of single-layer TMDCs, as we will now discuss.
Consider a charge density ρ(x, y)δ(z) in the z = 0 plane of the 2D material, embedded in an isotropic medium of permittivity . The resulting electric displacement field is where φ is the electrostatic potential, P ⊥ (x, y) is the inplane polarization, and κ is the in-plane susceptibility of the material. By using Gauss's law, ∇ · D = ρδ(z), we obtain After taking the Fourier transform, denoting the wavevector in the (x, y) plane by q and the wavenumber in the z direction by k, we find However φ(q, z = 0) = 1 2π φ(q, k) dk Rearranging, we find the in-plane electric potential to be φ(q, z = 0) = ρ(q) q(2 + qκ) .
Therefore the electrostatic potential energy between charges q i and q j in a 2D semiconductor is where r * ≡ κ/(2 ). After taking the Fourier transform, the potential energy can be written as where r is the separation of the particles and where H n (x) is a Struve function and Y n (x) is a Bessel function of the second kind. This result was first derived by Keldysh,31 and we refer to the interaction of Eq. (12) as the Keldysh interaction. At long range (r r * ) this potential becomes a Coulomb interaction: while at short range (r r * ) it is approximately logarithmic: where γ is Euler's constant. We refer to the interaction potential of Eq. (14) as the logarithmic interaction. The Keldysh interaction is plotted in Fig. 3, along with the Coulomb (r * = 0) and logarithmic (r * → ∞) approximations.  (11). The inset shows the percentage error in different approximations [Eqs. (13), (14), and (15)] to the Keldysh interaction of Eq. (12).
The following approximation to Eq. (12) was introduced in Ref. 32: This form of potential was used in the diffusion quantum Monte Carlo (DMC) study of Ref. 21. It is also plotted in Fig. 3, where it can be seen that the error in Eq. (15) is as large as several percent in the region r ≈ r * . We compare DMC results obtained using Eqs. (12) and (15) in Sec. IV F. Finally, the Mott-Wannier-Keldysh Schrödinger equation for a set of charged quasiparticles in a 2D semiconductor is where m i and q i are the band effective mass and charge of particle i, r ij is the separation of particles i and j, and E is the energy eigenvalue. Now consider the situation in which the 2D semiconductor has a dielectric medium of permittivity a above it and a dielectric medium of permittivity b below it, as would be the case for a 2D semiconductor deposited on a substrate. In general this is a more complicated problem than the situation described above. However, if we take ≡ ( a + b )/2 in the expressions above, the correction to the electrostatic energy of Eq. (11) is second order in a − b . Hence the Keldysh interaction remains valid when the permittivity is chosen to be the average of the permittivities of the media on either side of the 2D semiconductor, provided these permittivities are similar.

Excitonic units
The energies of complexes interacting via the Keldysh or Coulomb interactions are given in terms of the exciton Rydberg, R * y = µe 4 /[2(4π ) 2 2 ], and lengths are given in terms of the exciton Bohr radius, a * 0 = 4π 2 /(µe 2 ), where µ = m e m h /(m e + m h ) is the reduced mass of electron-hole pairs, with m e and m h being the electron and hole masses, respectively.
Letr i = r i /a * 0 . Then Eq. (16) can be written as where E = E/R * y . Note that µ/m i only depends on the electron-hole mass ratio σ ≡ m e /m h . Hence for a fixed value of r * /a * 0 , the dimensionless energy eigenvalues E only depend on the mass ratio, not on the absolute masses. Furthermore, for an exciton we may write the Schrödinger equation in terms of the difference coordinate r eh as so that for a given value of r * /a * 0 , the dimensionless exciton energy eigenvalues E X are also independent of the mass ratio. For the case of the Coulomb interaction (r * = 0), the dimensionless ground-state energy of an isolated exciton is E X = −4, irrespective of the mass of the electron or the hole. The binding energies in excitonic Rydbergs of donor-bound trions, biexcitons, and donor-bound biexcitons only depend on r * /a * 0 and the electron-hole mass ratio σ. Unfortunately, the energies of the different complexes go to zero in these units in the limit that r * → ∞, and so a separate set of units is required for the case of the logarithmic interaction, as discussed in Sec. III B 2.

Logarithmic interaction
For the limit r * → ∞, where the interaction is of logarithmic form, we use the dimensionless units introduced in Ref. 19. The Schrödinger equation for a charge car-rier complex with the logarithmic approximation to the interaction [Eq. (14)] is and Defining dimensionless coordinatesr i = r i /r 0 and a dimensionless energy E = E/E 0 , the Schrödinger equation can be written as The only dependence of the dimensionless energy E of the complex on r * is through the pairwise additive constant Note that for an exciton or donor atom −1 for a trion or donor-bound exciton −2 for a biexciton or donor-bound trion −2 for a donor-bound biexciton (n+−n−) 2 −n+−n− 2 for a complex of n + charges +e and n − charges −e .
Hence the additive constant C cancels out of the binding energies of the different charge-carrier complexes defined in Sec. III C.
For an isolated exciton, we may write the Schrödinger equation in terms of the difference coordinate r eh and reduced mass, giving The only dependence of the dimensionless energy eigenvalue E X on the mass ratio and r * comes from the constant term log(r 0 /r * ) in the Hamiltonian. Hence we may write the ground-state dimensionless energy as where E X0 = 0.41057747491(7) was evaluated by a finite-element method (see Sec. V A).

C. Binding energies and spectra of charge-carrier complexes
We define the binding energies D 0 X , and E b D + XX of a trion, biexciton, donor-bound exciton, donor-bound trion, and donor-bound biexciton, respectively, as follows: where The energy difference between the exciton peak in a photoluminescence experiment and the peak corresponding to a particular complex is equal to the energy required to separate a single exciton from that complex. Thus the energy difference between the exciton peak and the trion peak is E X − E X − = E b X − ; the energy difference between the exciton peak and the biexciton peak is 2E X − E XX = E b XX ; and the energy difference between the exciton peak and the donor-bound trion peak is On the other hand, the energy difference between the exciton peak and the donor-bound exciton peak is E X − E D + X = E b D + X + E X − E D 0 , and the energy difference between the exciton peak and the donor-bound biexciton peak is Some of these peaks are shown in Fig. 4. In addition there are expected to be offsets to the peak positions due to the spin-splitting of the conduction bands of TMDCs, as described in Sec. III.
In Sec. V F we report DMC binding energies for quintons and other large charge-carrier complexes in tungsten and molybdenum dichalcogenides. In each of these cases the binding energy is defined to be the energy required to remove an exciton from the complex; this is the binding energy with respect to dissociation into the most energetically competitive products.

D. Contact and exchange interactions between charge carriers
The Mott-Wannier model of a charge-carrier complex is valid provided the complex extends over many unit cells of the underlying crystal. However, when charge carriers are present at the same point in space there is an energy contribution due to local exchange and correlation effects. 25 Although the excitons in TMDCs are Mott-Wannier-like, their wave functions only extend over a small number of primitive unit cells, so that local exchange and correlation effects are expected to be significant. We may represent this effect within a Mott-Wannier model by introducing additional pairwise contact interaction potentials. For example, for a biexciton the Hamiltonian should include an additional term of the form where A ee , A hh , and A eh are constants and r ee , r hh , and r eihj are the electron-electron separation, the hole-hole separation, and the separation of electron i and hole j, respectively. Evaluating A ee , A hh , and A eh by ab initio calculations is challenging, and so we leave them as free parameters to be determined in experiments or subsequent ab initio calculations. If we evaluate the expectation value of this contact interaction then we find that the first-order perturbative correction to the total energy can be written as A eh g eh XX (0) + A ee g ee XX (0) + A hh g hh XX (0), where the electron-electron, hole-hole pair, and electronhole PDFs are respectively. We report contact PDF data within the Mott-Wannier model.
In addition to the role of the contact PDF in evaluating perturbative corrections due to contact interactions, the PDF and contact PDF contain a wealth of physical information. The exciton recombination rate of a chargecarrier complex is proportional to the electron-hole contact PDF. Furthermore, the PDF gives a very direct indication of the spatial size and shape of a charge-carrier complex.
The contact PDF also plays a role in the intervalley scattering of carriers. As the intervalley scattering involves a large momentum transfer of the order of the inverse lattice constant, the interaction is short range and can be modelled by a contact interaction with both carriers in the same position. In particular, the electronelectron contact PDF for the semidark trion and biexciton in tungsten-based TMDCs determines the coupling strength of the dark and bright states as µ bd ∝ g ee (0) and hence determines the recombination rates of the semidark states. 33 3 EX is the total energy of an exciton. The lines show the frequency relative to the bright exciton peak arising due to the recombination of a single electron-hole pair in each complex: see Sec. III C. E.g., the D 0 X line shows the frequency relative to the exciton peak of the process D 0 X → D 0 + γ. The trion and biexciton peaks labelled "SD" arise from semidark complexes, and are offset by 2∆ , as explained in Sec. II A; the exciton peak labelled "dark" arises from the process described in Sec. II B; the other peaks arise from bright complexes. Donor-and acceptor-bound exciton peaks are shown with very low intensity due to the marginal stability of these complexes.

A. Quantum Monte Carlo modelling of excitonic complexes
Our total-energy and PDF calculations were carried out using the variational quantum Monte Carlo (VMC) and DMC approaches. 34,35 The ground-state wave function for a set of interacting, distinguishable particles is nodeless; hence the fixed-node DMC algorithm is exact for all the systems studied in this work with the exception of biexcitons with indistinguishable holes. We used a numerical representation of the potential of Eq. (12) that is accurate to at least eight significant figures. Trial wave functions were optimized using VMC with variance minimization 36,37 and energy minimization. 38 The DMC calculations were performed using time steps in the ratio 1 : 4 with the corresponding target configuration populations being in the ratio 4 : 1. Afterwards, the energies were extrapolated linearly to zero time step and hence, simultaneously, to infinite population. To perform all our calculations, the casino code was used. 39 QMC methods have previously been used to study 2D trions with nonlocal screening 19,21 and the Coulomb interaction 22 and 2D biexcitons with the Coulomb interaction (including indirect biexcitons in coupled-quantumwell heterostructures) 17,18,40,41 and in TMDCs with nonlocal screening. 21,23 In a recent work some of the present authors have investigated the binding energies of trions and biexcitons using DMC for a range of susceptibility parameters r * and effective masses, and have represented the DMC data using simple interpolation formulas. 25 It was shown that for the applicable range of r * values, 2D semiconductors are expected to show larger trion binding energies than biexciton binding energies, in contrast to the situation in quasi-2D systems such as GaAs/InAs quantum wells. Here we extend this work to include extreme cases and donor/acceptor-bound carrier complexes.

B. Wave functions for complexes of distinguishable charge carriers
Our trial wave functions for complexes of distinguishable charge carriers were of the Jastrow form Ψ = exp[J(R)], where R is the vector of all the particle coordinates. The Jastrow exponent J(R) included a pairwise sum of terms of the form 42 for the Keldysh and logarithmic interactions, where r is interparticle distance, c 1 , c 2 ≤ 0, and c 3 ≥ 0 are optimizable parameters, and for distinguishable pairs of particles of charge q i and q j and mass m i and m j . Different constants c i are used for each type of particle pair. This form satisfies the analog of the Kato cusp conditions, 43,44 i.e., it ensures that the local energy Ψ −1Ĥ Ψ is nondivergent at coalescence points, whereĤ is the Hamiltonian operator.
Where the interaction between the charge carriers was of Coulomb form, we used pairwise terms of the form in the Jastrow exponent, where c 1 ≤ 0 and c 2 ≥ 0 are optimizable parameters, and for distinguishable pairs of particles of mass m i and m j and charge q i and q j . This form satisfies the Kato cusp conditions. 43,44 Donor ions and other infinitely heavy particles were fixed point charges in our calculations. In this case u ex2D provided a one-body Jastrow term between the free particles and the fixed particles that satisfies the Kato cusp conditions. In addition, cuspless one-body, two-body and three-body polynomial terms truncated at finite range were used in our Jastrow factor. 45,46

C. Wave functions for biexcitons with indistinguishable holes
For biexcitons with indistinguishable holes we used the trial wave function where J is of the form described in Sec. IV B. For indistinguishable particles of mass m and charge q interacting via the logarithmic or Keldysh interactions, Eq.
(37) must be replaced by Γ = −q 2 m/(8a * 0 µe 2 r * ), while for indistinguishable pairs of particles interacting via the Coulomb interaction, Eq. (39) must be replaced by where η hh and η eh are smoothly truncated polynomials, with optimizable expansion coefficients, and r hh and r eihj are the hole-hole and electron-hole relative positions, respectively. Equation (41) is effectively a backflow 47,48 transformation: Ψ = exp(J)x hh introduces the correct nodal topology for the state that we want to consider and Eq. (41) maps the particle coordinates {r} to quasiparticle coordinates {r } without changing the nodal topology. In Eq. (41), and are smoothly truncated polynomials with optimizable parameters {a n } and {b n }. L is a cutoff length, N hh η and N eh η determine the amount of variational freedom, C = 3 to ensure smooth behavior at the cutoffs, and Θ denotes the Heaviside function. We require b 1 = Cb 0 /L to ensure that η does not affect the Kato cusp conditions, which are enforced by the Jastrow factor. We optimized the free parameters in our antisymmetric wave function using energy minimization. 38 For different values of N hh η and N eh η in Eqs. (42) and (43), we compare the VMC ground-state energy, variance, and DMC energy of biexcitons with indistinguishable electrons interacting via the logarithmic interaction in Table II. Analogous results for biexcitons interacting via the Keldysh interaction at finite r * are shown in Table III. Our results show that increasing N hh η and N eh η slightly decreases the variances; nevertheless, the VMC and DMC energies are independent of the number of free parameters when N hh η , N eh η ≥ 2 to within our statistical error bars. We have used N hh η = N eh η = 3 in our production calculations.
Biexcitons with distinguishable electrons and indistinguishable holes can trivially be mapped onto biexcitons with indistinguishable electrons and distinguishable holes by charge conjugation.

D. Time-step and population-control errors
We chose our DMC time steps such that the rootmean-square distance diffused by each particle in a single  (42) and (43)] on the VMC ground-state energy (EVMC), VMC energy variance, and DMC energy (EDMC) for biexcitons with indistinguishable holes interacting via the logarithmic interaction. The mass ratio is σ = 0.1 and the reduced mass is µ = 0.5m0, where m0 is the bare electron mass. In each case r * = r0.  (42) and (43)] on the VMC ground-state energy (EVMC), VMC energy variance, and DMC energy (EDMC) of biexcitons with indistinguishable holes interacting via the Keldysh interaction, with an electron-hole mass ratio of σ = 0.1. r * = 0 corresponds to the Coulomb interaction.
time step was much less than r 0 for the logarithmic interaction, much less than a * 0 for the Coulomb interaction, and much less than min{r 0 , a * 0 } for the Keldysh interaction at finite r * . In Fig. 5 we plot the DMC total energy of a biexciton with distinguishable particles against time step. The figure confirms that the linear extrapolation scheme described in Sec. IV A largely eliminates the effects of time-step bias, provided the time steps used are sufficiently small. For the logarithmic interaction with σ = 1 and r * = r 0 , the time step should evidently be rather less than 0.04 /E 0 .

E. PDF calculations
The PDFs defined in Sec. III D were evaluated by binning the interparticle distances sampled in VMC and DMC calculations. The errors in the VMC and DMC PDFs are linear in the error in the trial wave function; however, the error in the extrapolated estimate (twice the DMC estimate minus the VMC estimate) is quadratic in the error in the trial wave function. 49 Our reported PDFs were obtained by extrapolated estimation.
Contact PDF data have been calculated by extrapolating electron-hole and electron-electron PDFs to zero separation for each r * value and mass ratio considered. To perform the extrapolation we fitted exp[g(r)] to our PDF data at short range, 50 wherẽ g(r) = a 0 + 2Γ r 2 log(r) + a 2 r 2 + a 3 r 3 + · · · + a 6 r 6 (44) for the Keldysh and logarithmic interactions and g(r) = a 0 + 2Γr + a 2 r 2 + · · · + a 6 r 6 (45) for the Coulomb interaction (r * = 0), where Γ and Γ are defined in Eqs. (37) and (39) and a 0 , a 2 , . . . , a 6 and a 0 , a 2 , . . . , a 6 are fitting parameters. These forms satisfy (the analog of) the Kimball cusp conditions. 51 The model functions were fitted to our PDF data at small r, with the data being weighted by 2πr.

F. Sensitivity of binding energy to the form of screened interaction
We have investigated whether the approximation to the Keldysh interaction given in Eq. (15), which has been used in previous QMC studies of excitonic complexes, 21 leads to significant errors. For an exciton with r * = a * 0 /2, the DMC total energies are E X = −1.5358899(2)R * y and −1.4668074(3)R * y with the Keldysh interaction [Eq. (12)] and the approximate Keldysh interaction [Eq. (15)], respectively. This is a difference of about 4.5%, which is small but non-negligible. The DMC binding energies of trions with r * = a * 0 /2 and mass ratio σ = 1 using the exact and approximate Keldysh interactions are 0.1377(4)R * y and 0.1335(3)R * y , respectively, so the error in the binding energy due to the approximate Keldysh interaction is about 3%. Since these errors are easily avoidable, we have used the exact Keldysh interaction in our production calculations.

A. Excitons
The exciton ground-state energy is presented in Fig. 8. Our DMC data are in agreement with the results of finiteelement calculations as implemented in the Mathematica software. 52 In excitonic units, the energy of an exciton is independent of the effective masses: see Sec. III B. In the Coulomb limit, one recovers the well-known excitonic energy of −4R * y . We can determine the behavior of the energy near the Coulomb limit by evaluating the firstorder perturbative correction where ∆v = v Keldysh − v Coulomb is the difference between the Keldysh potential of Eq. (12) and the Coulomb potential of Eq. (13), and the expectation value is taken with respect to the exact ground-state wave function for the Coulomb interaction Ψ = exp(−2r/a * 0 ). The correction is shown in Fig. 8 as a green line. We have numerically evaluated the dimensionless constant E X0 in Eq. (26) to be E X0 = 0.41057739(7) using DMC and E X0 = 0.41057747491(7) using the finiteelement method. These results confirm the expected accuracy of the DMC method. The logarithmic-limit behavior from Eq. (26) is also shown in Fig. 8 (red line) and matches the DMC data near r * → ∞. The difference ∆E X /E 0 between the exciton energies in units of E 0 with the Keldysh and logarithmic interactions at large r * was calculated numerically. Using the optimized ground-state wave function for the logarithmic interaction, we used VMC to evaluate the first-order perturbative approximation ∆E X /E 0 ≈ v Keldysh − v logarithmic . The results are presented in Fig. 9 and show that the leading-order error in the exciton energy due to the logarithmic interaction goes as a * 0 /r * . We fitted the function with a 5 = −29 + 2E X0 − a 1 − a 2 − a 3 − a 4 − log 2 to our DMC exciton energy data, where y = r * /(a * 0 +r * ) and the remaining {a i } and {b i } are six free fitting parameters. The fractional error in the fit of Eq. (47) to our DMC data is everywhere less than 0.5%.
Contact PDFs were extracted as described in Sec. IV E. An example of a fit to Monte Carlo-sampled PDF data is shown in Fig. 10, and our contact PDF results are shown in Fig. 11(a). In all our plots of contact PDFs the statistical error bars from the Monte Carlo calculation are smaller than the symbols. Unlike the DMC mixed estimate of the energy, the extrapolated estimate of the PDF depends on the stochastically optimized trial wave function and hence in some cases slight noise in the g(0) data is visible.
In the Supplemental Material we provide a program for evaluating our fit to the total energy of an exciton [Eq. (47)], as well as fits to the binding energies of biexcitons, trions, donor-bound excitons, donor-bound trions, and donor-bound biexcitons. 53 In addition, the program reports fits to contact PDFs for the different clusters.

Binding energies
We compare the stability of biexcitons with distinguishable and indistinguishable holes in the limit of the Coulomb interaction (r * = 0) in Fig. 12(a) and at r * = 8a * 0 in Fig. 12(b). We find that biexcitons with indistinguishable holes are unbound for σ 0.3, while biexcitons consisting of distinguishable particles are bound at all mass ratios. The binding energies at σ = 0 are obtained using the Born-Oppenheimer potentials as a function of heavy-hole separation r plotted in Fig. 13. We fitted U (r) = α + β √ r + γr + δr 2 , where α, β, γ, and δ are fitting parameters, to our DMC data to find the minimum and the curvature about the minimum of the Born-Oppenheimer potential. For the logarithmic interaction we fitted U (r) = ζ + η exp(−r/d) + κ log(r) to our data, where ζ, η, d, and κ are fitting parameters. The Born-Oppenheimer approximation in Fig. 12(b) for heavy holes is in agreement with our DMC calculations at small σ. Analogous results obtained with the logarithmic interaction are shown in Fig. 12(c). For σ 0.2, only biexcitons with distinguishable holes are stable. Hence it is only at extreme mass ratios, where exchange effects between the heavy particles are negligible, that biexcitons with indistinguishable particles are stable. Figure 14 shows DMC binding energies for biexcitons with distinguishable particles interacting via the Keldysh interaction as a function of x = σ/(1 + σ) and rescaled in-plane susceptibility y = r * /(a * 0 + r * ). Our results are in agreement with path-integral Monte Carlo (PIMC) data at finite r * , as shown in Fig. 15. 20 However, the PIMC data obtained by Velizhanin and Saxena have much larger statistical errors and they quoted a previous DMC result 54 at r * = 0 due to the infeasibility of PIMC in this case. The function containing 17 fitting parameters {a ij } and {b ij }, was fitted to our DMC binding-energy data, giving a fractional error of less than 1.5% everywhere. This choice of fitting function exhibits the correct behavior as σ → 0, as derived in App. B 1, and is also invariant under charge conjugation (m e ↔ m h ). Equation (48) accurately reproduces the DMC biexciton binding energies over the whole space of possible susceptibility and mass-ratio parameters, unlike the simple fitting functions reported in Ref. 25. The latter are by construction only valid in the currently experimentally relevant region and, because of the relative simplicity of the fitting function, give significantly larger fractional errors (up to 5%) than Eq. (48). The fitted binding energy can be evaluated using the program supplied in the Supplemental Material. 53 Binding-energy results in the limit of large r * , where the interaction is of logarithmic form, are given in Sec. V G.  In Table IV, we compare the DMC binding energies of biexcitons in monolayer TMDCs with experiment and with previous theoretical works. Our DMC binding energies are in good agreement with previous DMC binding energies where available, 21 and also with PIMC calculations. 23 The small differences between DMC results in the literature must be due to the use of different effective masses, etc. Unfortunately, the theoretical biex-  citon binding energies are up to three times smaller than those reported in experimental works. 15,16,55,56 There is also a striking, qualitative disagreement with the experimental works regarding the trion and biexciton binding energies: the Mott-Wannier model with the Keldysh interaction predicts that the trion has a larger binding energy than the biexciton, 21,25 while the experimental studies report that the biexciton peak occurs at lower energies than the trion peak in photoluminescence spectra (i.e., that the biexciton has a larger binding energy). The theoretical results are reported for a free-standing monolayer; any screening by the substrate and environment would further exacerbate the disagreement with experiment.
The ground-state wave function of a system of distinguishable particles is nodeless, and so DMC provides exact solutions to Mott-Wannier models of excitonic complexes. Hence the disagreement with experiment regarding the binding energies of biexcitons in 2D semiconduc-   IV. Total energies of excitons (X) and binding energies of biexcitons (XX) and trions (X − and X + ) with distinguishable particles for different monolayer TMDCs suspended in vacuum ( = 0). We compare our results with values reported in the literature obtained by DMC, PIMC, hyperspherical harmonics (HH), stochastic variational (SV), and variational (V) methods. For each complex we use the values of me, m h , and r * shown in bold in Table V to evaluate the fits of Eqs. (48) and (49). Note that the ditellurides adopt a 2H stacking arrangement in bulk and few-layer samples, which may complicate comparison with experiment.  biexciton is not in its ground state; 27 or (iv) the experimental spectra have been misinterpreted or the peaks have been misclassified.

DMC biexciton binding energy
As explained in Sec. III D, there should be an additional contact interaction between charge carriers; however, the Mott-Wannier model with the Keldysh interaction apparently provides a good description 21,25 of the energies of excitons and trions, and there is no obvious reason to believe that contact interactions should be more important in a biexciton than in a trion or exciton. Moreover, it is unlikely that the contact interactions could be responsible for the threefold difference between the theoretical and experimental biexciton binding energies.
The second possibility is that the Mott-Wannier model is in principle correct, but the band effective masses and in-plane susceptibilities used in the model are incorrect. These are taken from ab initio calculations, which might not provide a sufficiently accurate description of the electronic band structure. However, as shown in Sec. V B 2, the different mass ratios and in-plane susceptibilities reported in the literature do not significantly affect the binding energy; in fact the mass ratios and in-plane susceptibilities would need to be in error by more than an order of magnitude to explain the difference with experiment. Finally, if inappropriate model parameters are responsible for the disagreement with experiment regarding the biexciton binding energy, it is not clear why the Mott-Wannier model with the same parameters apparently provides a good description of excitons and trions.
We believe that the exciton that remains after exciton recombination in a biexciton is unlikely to be in an excited state, because the parent biexciton is in its nodeless ground state, which strongly overlaps with the product of the ground states of the two daughter excitons.
The misclassification of the experimental results may offer at least a partial explanation of the disagreement.
By considering the behavior of the photoluminescence emission intensity, it has been argued that the observed peaks do indeed correspond to trions and biexcitons. 15,16 However, another possibility is that they could correspond to charge-carrier complexes involving donor or acceptor ions. In particular, the energies required to remove excitons from donor-bound biexcitons (see Sec. V F) are similar to the experimentally observed "biexciton" binding energies. If donor-bound biexcitons are responsible for the experimentally observed "biexciton peak" then we might expect the intensity of the peak to depend strongly on the doping of the sample. It is possible that other large charge-carrier complexes could also contribute to the spectra.
None of these options offers an entirely satisfactory explanation of the discrepancy. Further experimental and theoretical modeling work is required in order to understand the excitonic properties of 2D semiconductors.

Sensitivity of binding energies to effective masses and in-plane dielectric susceptibility
In Table V we compare the DMC binding energies of biexcitons with distinguishable particles for a variety of effective masses and in-plane screening lengths obtained by different first-principles methods. Since a range of masses are reported in the literature, we have taken the average of the reported masses that were supposedly obtained using the same method. The different model parameters in the literature lead to a spread of about 1 meV in the theoretical binding energies.
The sensitivities of the exciton total energy and the trion and biexciton binding energies to the model parameters are reported in Table VI. The energies depend relatively weakly on the in-plane permittivity r * ; the errors arising from the uncertainty in the effective mass almost certainly dominate errors arising from the uncertainty in r * . The sensitivity of the exciton energy to the effective masses is an order of magnitude larger than the sensitivity of the trion binding energy, which is in turn an order of magnitude larger than the sensitivity of the biexciton binding energy. To account for the 30-40 meV disagreement with experiment over the biexciton binding energy the effective masses would have to be more than an order of magnitude larger than the ab initio values reported in Table V and/or the r * value would have to be an order of magnitude smaller. While there is still appreciable uncertainty in the ab initio effective mass and r * values, it seems very unlikely that both density functional theory and many-body GW calculations would be in error by more than an order of magnitude.

PDFs
In Fig. 16, we show the PDFs of biexcitons with distinguishable particles interacting via the logarithmic inter- action for two different mass ratios, σ = 0.4 and σ = 1. The long-range biexciton wave function is relatively independent of the mass ratio. However, at short range the electron-hole PDF shows a peak near the separation that corresponds to the minimum of the Born-Oppenheimer potential-energy surface, which gets more pronounced at extreme mass ratios. As expected, the physical size of the biexciton is a low multiple of r 0 . Figure 17 presents the electron-hole and electronelectron contact (r = 0) PDFs for a biexciton. Notice that g eh XX ≈ 2g eh X . Fits to the contact PDFs can be evaluated using the program supplied as Supplemental Material. 53 TABLE VI. Sensitivity of binding energies to the three parameters that characterize the Mott-Wannier-Keldysh model of excitonic complexes in 2D semiconductors suspended in vacuum. The derivatives are evaluated using the effective mass and in-plane permittivity parameters reported in bold for different TMDCs in Table V. m0 is the bare electron mass.

C. Trions
The binding energies of negative trions are presented in Fig. 18. We have fitted the function where x = σ/(1+σ), y = r * / (r * + a * 0 ), and the {a ij } and {b ij } are fitting parameters, to the DMC trion binding energies. Equation (49) satisfies the limiting behavior described in App. B 2, has 31 free fitting parameters, and the fractional error in the fit to our DMC data is everywhere less than 1%. Positive trion binding energies can be obtained by charge conjugating the corresponding negative trion. The program included in the Supplemental Material can be used to evaluate Eq. (49). The resulting trion binding energies for various TMDCs are shown in Table IV. It can be seen that, in contrast to the biexciton binding energies, the trion binding energies are in excellent agreement with the available experimental results. As shown in Table VI, trion binding energies are significantly more sensitive to the effective mass values than biexciton binding energies; nevertheless, the ab initio effective masses would need to be in error by an implausibly large amount to change the trion binding energies by more than a few meV. Binding-energy results in the limit of large r * , where the interaction is of logarithmic form, are given in Sec. V G. Figures 11(a) and 11(b) present the electron-hole and electron-electron contact PDFs of trions. The fitting functions can be found in the program supplied as Supplemental Material. 53

D. Donor/acceptor-bound excitons
We present the binding energies of donor-bound excitons in Fig. 19. For σ 1, the binding energy is close to zero. In this region, the calculations were especially difficult, since the complex tends to unbind very easily. Therefore, during the wave function optimization, the cutoff lengths for the Jastrow factor were fixed at small values, to force the complex to be bound. In the limit σ → ∞, the complex is expected to be unbound (see App. B 3), which is consistent with our results. Indeed, over a broad range of large electron-hole mass ratios and large r * values, the DMC binding energy of the donor-bound exciton is either zero or extremely small, such that the binding energy cannot easily be resolved in DMC calculations. The following 50-parameter fitting formula has a fractional error that is mostly less than 2% in fits to our DMC data: (50) In this expression x = σ/(1 + σ) and y = r * /(a * 0 + r * ), while the {a ij } are fitting parameters. Our fitting function can be evaluated using the program in the Supplemental Material. 53 We summarize our theoretical predictions for the binding energies of donor/acceptorbound excitons in various TMDCs in Table VII. Bindingenergies in the limit of large r * , where the interaction is of logarithmic form, are given in Sec. V G.
We have also calculated the electron-hole contact PDFs of donor-bound excitons, which are presented in Fig. 20. Our results confirm that the contact PDFs decrease to zero as σ → ∞, as expected, because the light hole becomes unbound in this limit. Contact PDFs can be evaluated using the program supplied as Supplemental Material. 53 E. Donor/acceptor-bound trions Figure 21 presents the binding energies of donor-bound trions. We have devised the following 30-parameter fitting formula: which includes the correct divergence as σ → ∞ and appropriate square-root behavior for the heavy-hole limit σ → 0 (see App. B 4). The {a ij }, {b i }, and {c ij } are fitting parameters. The fractional error in the fit to our DMC data is less than 3%. The program in our Supplemental Material can be used to evaluate Eq. (51). 53 Binding-energy results in the limit of large r * , where the interaction is of logarithmic form, are given in Sec. V G. Table VII reports theoretical binding energies for donor-bound trions with biexciton energies for several real materials. The binding energy of a donor-bound trion is slightly larger than the binding energy of a free biexciton. This leads us to expect two lines close together in the absorption/emission spectra of TMDCs, one corresponding to biexcitons, and another at slightly larger energy corresponding to donor-bound trions.
Contact PDFs for donor-bound trions have been extracted from our QMC data and are presented in Fig. 22 and in the Supplemental Material. 53  Table V), but for donor atoms (D 0 ), acceptor atoms (A 0 ), donor-bound excitons (D + X), acceptor-bound excitons (A − X), donor-bound trions (D 0 X), acceptorbound trions (A 0 X), donor-bound biexcitons (D + XX), and acceptor-bound biexcitons (A − XX). The binding-energy results are our theoretical predictions using Eqs. (50), (51), and (52), while the energies of donor and acceptor atoms are calculated using Eq. (47) with infinite hole and electron masses, respectively. Note that the binding energy is defined with respect to dissociation into the most energetically favorable products, which do not always include an exciton: see the definitions in Sec. III C.

Energy (meV)
Binding energy (meV) where x = σ/(1 + σ) and y = r * /(a * 0 + r * ), while the {a ij }, {b i }, and {c ij } are fitting parameters. This gives a fractional error of less than 3% everywhere when fitted to our DMC data. Equation (52) can be evaluated using the program supplied as Supplemental Material. 53 We summarize our theoretical predictions for the binding energies of donor/acceptor-bound biexcitons in various TMDCs in Table VII. Binding-energy results in the limit of large r * , where the interaction is of logarithmic form, are given in Sec. V G.
The behavior of a donor-bound biexciton in the limit of heavy electrons is discussed in App. B 5. In the limit of heavy holes (σ → 0), this complex consists of three fixed positive particles and two light electrons and thus the question arises of how the three fixed, positive charges are positioned with respect to each other. The most natural position that three positive particles would assume is an equilateral triangle. To check if this assumption is correct we first determined how the Born-Oppenheimer potential energy changes if we distribute the three positive charges in the corners of equilateral triangle and then vary the triangle side. Figure 24 shows the case of r * /a * 0 = 1 as an example. After finding the side length that minimizes the Born-Oppenheimer potential energy, we changed the position of one of the positive particles (fixing the remaining two) and again observed the effect on the Born-Oppenheimer potential energy. Figure 25 presents the results, which clearly show that the equilateral triangle is a local minimum of the Born-Oppenheimer potentialenergy surface.
Closely related to donor-bound biexcitons are fivecarrier complexes known as charged biexcitons or quintons (XX − , i.e., e − e − e − h + h + ). In molybdenum and tungsten dichalcogenides these consist of two distinguishable holes with opposite spin and valley indices, and three distinguishable electrons that differ in either their spin or their valley indices: see Fig. 1(a). The binding energy of a quinton is defined as the energy required to split it into a free exciton and a free trion. 81 Other possible large complexes are donor-bound doublenegative excitons (D − X, i.e., D + e − e − e − h + ), donorbound quintons (D 0 XX, i.e., D + e − e − e − h + h + ), and even donor-bound double-negative biexcitons (D − XX, i.e., D + e − e − e − e − h + h + ). For molybdenum and tungsten dichalcogenides there are no further possibilities: we have exhausted the possible neutral or singly charged complexes that can be constructed from up to four distinguishable electrons, up to two distinguishable holes, and zero or one donor ions. Any larger charge-carrier complexes in molybdenum or tungsten dichalcogenides inevitably either include indistinguishable particles or involve the much larger energies required to excite holes  in the lower spin-split valence bands. In Table VIII we present our DMC binding-energy results for quintons and the other large complexes. Donor-bound double-negative biexcitons appear to be unstable to dissociation into free excitons plus donor-bound double-negative excitons, and hence are not included in Table VIII. As with donorbound biexcitons, the energies required to remove excitons from the larger complexes such as quintons are

G. Complexes with the logarithmic interaction
We have also studied complexes of distinguishable particles interacting with the purely logarithmic form of Eq. (14). The binding energies are presented in Fig. 27. The lines shown in Fig. 27 were obtained using Eqs. (47), (48), (49), (50), (51), and (52). To convert from exci- tonic units to logarithmic units we multiply the fitting function by R * y /E 0 = r * /(2a * 0 ) = y/(2 − 2y) and take the limit that r * → ∞, i.e., that y → 1. For complexes that have been studied previously, our results are in good agreement with earlier exact calculations. 19

VI. CONCLUSIONS
In summary, we have discussed the different types of biexciton and trion that can be observed in molybdenum and tungsten dichalcogenides. Furthermore, we have presented statistically exact DMC binding-energy data for biexcitons, trions, donor/acceptor-bound trions, and donor/acceptor-bound biexcitons in 2D semiconductors, including an analysis of extreme mass ratios. We have shown that biexcitons with indistinguishable charge carriers are unstable at experimentally relevant electron-hole mass ratios. Our calculations have used the effective interaction between charge carriers arising from screening effects in such materials. We have also presented contact PDF data that allow the investigation of additional contact interaction energies between charge carriers in 2D semiconductors within first-order perturbation theory. Our work provides a complete reference for the interpretation of spectral lines in photoabsorption and photoluminescence experiments on monolayer TMDCs in terms of a model of charge carriers moving within the effective mass approximation.
A broad range of theoretical works on 2D biexciton binding energies show excellent quantitative agreement with each other, but an enormous, threefold disagreement with experiment. By contrast, for trions there is good agreement between theory and experiment. We have considered and discounted various possible deficiencies in the theoretical models of charge-carrier complexes. We believe that the most likely explanation for the disagreement with experiment is a misinterpretation or misclassification of experimental optical spectra. In particular, we note that the energies require to remove excitons from donor-bound biexcitons are similar to the binding energies of experimentally observed biexcitons, suggesting that larger charge-carrier complexes could be respon-     Table IX, is included for a given product by using + ⊗ + = +, + ⊗ − = − and − ⊗ − = + and noting that all irreps in a given direct sum have the same Cs classification. ω = 2U min m h + m e = 2R y U min a 2 0 (m h + m e ) .

(B2)
The resulting ground-state energy in the harmonic approximation is where we have used m h m e in the last step. This suggests that a suitable fitting function for the binding energy of a biexciton with small σ ≡ m e /m h is a polynomial in powers of √ σ. Similar conclusions hold for the case where the interaction between the charge carriers is logarithmic.
In the limit of heavy holes, the total energies of biexcitons with distinguishable and indistinguishable holes are identical, because exchange effects become negligible as the heavy holes localize. Hence a biexciton with indistinguishable holes must be bound when the hole mass is sufficiently large. Likewise, a biexciton with indistinguishable electrons has the same total energy as a biexciton with distinguishable electrons in the limit that the electron mass is large.

Negative trions
In the limit of heavy holes (σ → 0), a negative trion resembles a 2D H − ion. The leading-order correction to the energy of an infinite-hole-mass negative trion is therefore due to the reduced-mass and mass-polarization perturbative corrections encountered in atomic physics, each of which gives a contribution to the energy that is linear in the electron-hole mass ratio σ.
In the limit of heavy electrons (σ → ∞), a negative trion resembles a charge-conjugated 2D H + 2 ion, and hence one can use the Born-Oppenheimer and harmonic approximations, as was done in App. B 1. The binding energy near the extreme mass limit varies as the squareroot of the mass ratio σ.

Donor-bound excitons
A donor-bound exciton in the limit of heavy holes is a charge conjugate of a negative trion with heavy electrons, and therefore will have a binding energy that varies as the square root of the mass ratio σ.
In a donor-bound exciton with heavy electrons, the positive donor ion and the heavy electron overlap, so the light hole is unbound. Therefore the binding energy in this limit is zero. equal to the exciton mass. The binding energy varies as the square root of the mass ratio σ. Now consider a donor-bound trion with two heavy electrons and a light hole. If the hole is very much lighter than the electrons then the hole will be extremely delocalized and will see the positive donor ion and two electrons (D − ) as a fixed, negative point charge; the system therefore resembles an acceptor atom in which the hole is bound to a fixed, negative point charge. Hence E D 0 X ≈ E D − + E A 0 in this limit, where E A 0 is the energy of an acceptor atom. In addition, if the electron mass is very much larger than the hole mass, the exciton ground-state energy is E X ≈ E A 0 . The binding energy of a donor-bound trion in the limit that the hole is much lighter than the electron mass is therefore which is the electron affinity of a donor atom. Note that the electron affinity of a donor atom is equal to the binding energy of a negative trion in the limit of large hole mass.
The exciton Rydberg goes to zero in the limit that the hole mass goes to zero; hence the binding energy of a donor-bound trion in excitonic units goes to infinity as the hole-to-electron mass ratio goes to zero (σ → ∞).

Donor-bound biexcitons
A donor-bound biexciton with two heavy holes resembles a trihydrogen cation (H + 3 ). This molecular ion is an important component of the interstellar medium, 82 and it is known that the protons in H + 3 form an equilateral triangle. In Sec. V F we verify that 2D donor-bound biexcitons with heavy holes also adopt an equilateral triangular structure, and we calculate the bond length by minimizing the Born-Oppenheimer potential energy.
Consider a donor-bound biexciton with two heavy elec-trons and two light holes. The binding energy of a donorbound biexciton in the limit that the holes (h + light ) are much lighter than the electrons (e − heavy ) is which is the hole affinity of an acceptor atom (in the limit of large electron mass, D + e − e − acts like a fixed negative point charge). Note that the hole affinity of an acceptor atom is equal to the binding energy of a positive trion in the limit of large electron mass.