Charge carrier complexes in monolayer semiconductors

.


I. INTRODUCTION
The optical properties of layered semiconductors, such as transition-metal dichalcogenides (TMDs), change as the sample thickness is reduced from bulk (B) down to a single layer (1L) [1].Indirect band gaps in B-TMDs are often observed to transition to direct band gaps in 1L [2], accompanied by the emergence of photoluminescence (PL) [3].Excitonic effects are enhanced in 1Ls relative to B-TMDs, due to reduction in electrostatic screening of the interactions between charge carriers [4].In many 1L-semiconductors, including TMDs with honeycomb lattices, spin-orbit coupling splits the conduction (CB) and valence (VB) bands at their extrema at the K and K ′ points of the Brillouin zone [5].This results in optically controllable spin and valley degrees of freedom [6][7][8].Valley polarization is retained for> 1ns [6], ideal for quantum device applications, such as quantum lightemitting diodes [9][10][11][12].Localized single-photon emitters that can be controlled by electroluminescence [9,12] are also promising for quantum photonics.
The binding energy (BE) of an exciton may be calculated from first principles by solving the Bethe-Salpeter equation (BSE) [13] on top of manybody perturbation theory calculations within the GW approximation [14][15][16], or by quantum Monte Carlo (QMC) methods [17].However, studying charge-carrier complexes, such as quintons, using these approaches is computationally expensive [18].Instead, the effectivemass approximation [19] can be used, whereby the ground-state energy is modelled by considering an electron (e) and a hole (h) interacting within a two-band model [20], and their effective masses are defined by experiment or by first principles band structure calculations.In effective-mass models of charge-carrier complexes in layered semiconductor materials (LSMs), it is crucial to take into account the two-dimensional (2d) na-ture of the electrostatic screening, as this modifies the form of the interaction between carriers [21,22].The situation for LSMs differs from III-V semiconductor heterostructures with a thickness> 1µm [23], in which the Coulomb 1/r interaction between charge carriers scales down with the permittivity of the host material [23].In LSMs, the so-called Rytova-Keldysh interaction (RKI) potential [21,22] provides a more accurate interaction between charge carriers.Refs.24-26 studied the formation of multicarrier bound states in 1L-semiconductors using QMC methods, such as path integral Monte Carlo [27] and diffusion Monte Carlo (DMC) [18] to solve the Schrödinger equation for quasiparticles interacting via the RKI potential.DMC is particularly powerful in studies of complexes with distinguishable quasiparticles [28], as it is numerically exact in this case [29].Ref. 26 used DMC to predict the stability of negative quintons in TMDs with distinguishable charge carriers (in which all three e species have different spin and/or valley degrees of freedom).These predictions were confirmed in experimental studies, which provided evidence of quintons in hBN-encapsulated 1L-WSe 2 [30][31][32][33], 1L-MoSe 2 on sapphire [34], and 1L-WSe 2 on Si/SiO 2 [35].
In Mo and W-TMDs, the VB spin-splitting is sufficiently large that the lower spin-split bands are always occupied at room temperature (RT) [5], while the CB spinsplitting is comparable with RT [5].As a result, there are effectively 4 e and 2 h species available to form chargecarrier complexes at and below RT [26].
One can distinguish dark [36], bright [36], and semidark [37] charge-carrier complexes in TMDs.In dark complexes (Fig. 1a), radiative e-h recombination is not allowed due to spin and/or momentum mismatch between the constituent e/h [38], while in bright complexes (Fig. 1b), direct radiative e-h recombination is allowed by conservation of linear and angular momentum [39].In semi-dark complexes (Fig. 1c), radiative recombination can take place following an intervalley scattering event assisted by a phonon that maintains spin, but swaps an e, e.g., from valley K ′ to K [37], accompanied by an energy shift due to the change in occupation of the upper and lower spin-split bands [37].
Due to the nature of the CB spin-splitting of Mo-TMDs [5] (Fig. 1a), bright states are energetically lower than dark [40].Hence, at low temperature T < 100K, e in exciton (X), negative trion (X − ), and biexciton (XX) complexes occupy the lower spin-split bands.X complexes therefore travel only a small distance, e.g.∼ 1µm in 1L-Mo-TMDs [41,42], before radiative recombination, which reduces the chance to bind with another chargecarrier complex [43].Furthermore, the XX PL peak may be difficult to distinguish from that of X − , due to the small energy difference∼ 10meV between their BEs [41].Ref. 34 detected XX and quintons (XX − ) in 1L-MoSe 2 by 2d coherent spectroscopy (2dCS) [44].This method can focus on a delay time∼10ps, over which XX or XX − are likely to form [45,46]. However in 1L-W-TMDs, the most energetically stable excitonic states are dark [40], so that X have longer lifetimes (∼1ps) [40] than in 1L-Mo-TMDs (∼0.5ps) [40], favouring larger than X charge-carrier complexes.We therefore focus on 1L-W-TMDs when comparing theory with experiments.(1) those with 1e in the upper spin-split CB and 2e in the lower spin-split CB, and (2) those with 1e in the upper spin-split CB and 1e in the lower spin-split CB.The CB spin splittings in 1L-Mo-and W-TMDs are∼ 3meV [5] and∼ 30meV [5], respectively.These are much less than the XX − BEs∼ 50meV [30,34,47], as reported in Table I.The fact that the XX − BE is larger than the spin splitting implies XX − complexes are thermodynamically stable at T close to 0K, even taking into account the energy required to excite 1e to the upper spin-split CB.Assuming the CB spin-orbit splitting ∆ ′ of 1L-TMDs to be≪XX − BE, E b XX − , each XX − can be treated as a twostate system [48].For k B T ≪ ∆ ′ ≪ E b XX − , with k B the Boltzmann's constant, the fraction of XX − with 1e and 2e in the upper spin-split CB, hence the PL intensity of the corresponding XX − , is [49]: for 1 e in upper spin band e −∆ ′ /(kBT ) for 2 e in upper spin band Here, we use DMC within the effective-mass approximation to calculate XX − energies in 1L-semiconductors.XX − are the largest free charge-carrier complexes in TMDs [26].We provide an interpolation formula for XX − BEs for all 1L-semiconductors as a function of e and h effective masses, permittivity of the surrounding media, and in-plane susceptibility of the 1L-semiconductor.We also use DMC to calculate the energies of chargecarrier complexes in the presence of out-of-plane magnetic and in-plane electric fields, to identify whether the behavior in external fields can be used to investigate PL peaks.We find that applying an external magnetic field helps identifying charge-carrier complexes in 1L-semiconductors which have different e and h effective masses, while electric fields can be used to identify charge-carrier complexes in all 1L-semiconductors.We explore the accuracy of the RKI potential by comparing BEs with results obtained using ab initio random-phase approximation (RPA) interaction potentials [50].We find that, within the effective-mass approximation, RKI can describe quasiparticles on length scales larger than the lattice constant.Therefore our results can be used to determine the PL spectra of excitonic charge complexes. .The e and h effective masses me and m h in terms of the free e mass m0 are taken from many-body GW calculations [14,51].We assume the materials suspended in vacuum, or encapsulated in hBN, or placed on top of a substrate such as SiO2 [35] and sapphire [34].Wherever ǫ = ǫ0, the vacuum r * is used Ref. 26 give XX and X − BEs∼19.9,29.2meV, respectively.c The experimental X − BE is∼ 30meV [35].With r * = 44.33Å and ǫ = ǫ 0 , Eq.49 of Ref. 26 gives BE∼ 30meV.d The experimental XX BEs are∼ 18.2meV [30] or∼ 20.1meV [47], and those of X − are∼ 27.1meV [30] or∼ 29.7meV [47].With r * = 48 Å and ǫ = ǫ 0 , Eqs.48,49 of Ref. 26 give XX and X − BEs∼18.9,28.1meV, respectively.
To convert the Bohr radius, flux density, electric field, and energy from e.u. to SI units, each value needs to be multiplied by a * 0 , B * , F * , and 2R * y , respectively.To convert from logarithmic e.u. to SI units, each value needs to be multiplied by √ 2r 0 , B 0 , F 0 , and E 0 , respectively.Table II summarizes all acronyms used in this paper.

B. Binding energies
We define the X, X − , and XX BEs as: where the complexes are defined in Table II.In the absence of external fields, E e = E h = 0. We define the de-excitonization energy of XX − as: and the electron affinity of XX as: By comparing Eq.7 with 5,6 for 1L-TMDs, the energy difference between bright X and XX − PL peaks is of XX − complexes in all 1L-semiconductors with all the possible values for r * /a * 0 = {0, 0.5, 1, 2, 4, 6, 8, ∞} and σ = {0, 0.1, 0.2, . . ., 1, 1.5, 4, 9, ∞}, where σ = m e /m h is the mass ratio.We fit: , where {a ij }, {b i }, {c i }, and {d i } are fitting parameters, x = σ/(σ + 1) = m e / (m e + m h ) is a rescaled mass ratio and y = r * /(r * + a * 0 ) is a rescaled in-plane susceptibility parameter.The fitting function goes as the square root of the mass at extreme mass ratios (σ = 0 and σ = ∞), as required by the Born-Oppenheimer approximation [58].We fit ] so that the asymptotic behavior at r * → ∞ obtained using the logarithmic interaction can be included in the fit.The error in the fitted E DE XX − is< 5% at each data point.The statistical error bars on the DMC E DE XX − data points are much smaller than the error in the fit.We therefore use an unweighted least-squares fit [59].We provide a program [60] that can be used to evaluate E DE XX − and XX − BE for any 1Lsemiconductor, for which effective masses, r * , and dielectric constant of the environment are the inputs.
Using the BE fits in Eqs.48,49 of Ref.26 together with Eqs.6,7,8, we calculate the XX − BEs in Fig. 2. Above the yellow line in Fig. 2c, the XX − BE =E DE  XX − (all 1L-TMDs fall in this region).Below the yellow line, the XX − BE is equal to the XX electron affinity.
Table I  The size of a charge-carrier complex in a 1L-TMD can be defined by r 0 .This is∼ 8 Å in TMDs listed in Table I, because their e and h masses and screening lengths are around the same order.Hence, we suggest that encapsulation in> 1nm ML-hBN can be described by the permittivity ǫ = 4ǫ 0 [61][62][63][64] of bulk hBN.
For ML-hBN-encapsulated TMDs we test 2 approaches to compare our results with experiments.1) We fix ǫ = ǫ 0 and determine r * by fitting Eqs.48,49 of Ref.26 to the experimental X − , XX BEs in Refs.30, 34, 35, and 47.This is reasonable because, at distances larger than the layer-layer separation, the Keldysh interaction of Eq.15 for a ML is of the same form as for a 1L [65], but with r * being the sum of r * for the different layers [65].
For r ≪ r * , only ǫr * appears in the logarithmic approximation to the Keldysh interaction in Eq.16, apart from a constant contribution to the total energy, which cancels out of E DE XX − .Hence, it is preferable to fix ǫ, and treat r * as the independent parameter.The XX − BEs calculated with this approach for 1L-WSe 2 and 1L-WS 2 encapsulated in ML-hBN agree with the experiments in Refs.30 and 47, differing by at most∼ 2meV, as for Table I.
2) We use ab initio vacuum r * for 1L-TMDs.To describe hBN encapsulation we use ǫ = 4ǫ 0 , consistent with Refs.61-64.This gives BEs∼5-18meV smaller than Refs.30and 47.This difference could either due to the phenomenological parameters (r * and ǫ), obtained by ab initio methods and used in the Mott-Wannier-Keldysh model of Eq.14, or be a result of neglecting intervalley scattering [66] and contact (exchange) interactions.
The substrate effect on the BE of a charge-carrier complex can be described by: where ǫ substrate is the bulk permittivity of the substrate.In Ref.  [34].Substrate-induced roughness can also cause inhomogeneity in the electronic structure and extra carrier scattering [69].This affects PL, leading to inhomogeneous broadening [70,71], which makes it difficult to identify charge-carrier complexes [34].In Ref. 34, PL spectra were not recorded as a function of excitation power.However, power-dependent measurements help assign the PL peaks to X − , XX, XX − , because they show, respectively, sublinear, quadratic, and superlinear dependence with excitation power.Due to the complexity in defining r * and ǫ for 1L-TMDs encapsulated in hBN or placed on a substrate, we use the first approach, where we fix ǫ = ǫ 0 and determine r * by fitting theoretical X − and XX BEs to available experiments, so to define the XX − BEs.

C. Other charge-carrier complexes
We investigate doubly charged complexes, including XX 2− (4e and 2h) and D 0 hh (one positive donor ion, 1e, 2h), in which all charge carriers are distinguishable.Optimizing wave functions with pairwise and three-body correlations by variational Monte Carlo (VMC) energy minimization [72] does not result in bound-state wave functions.If we constrain the wave function to be bound, and then perform DMC, the resulting energy confirms that complexes are unbound.Thus, doubly charged complexes are unstable for all relevant material parameters.
In Ref. 26 we considered what is the largest stable charge-carrier complex that can occur in 1L-TMDs.We showed that XX with two indistinguishable e is unstable in 1L-TMDs, because of the resulting antisymmetry of the spatial wave function.We concluded that, in bound complexes featuring only singly charged dopant ions and charge carriers, all charge carriers must be distinguishable.Our results show that a charge-carrier complex can feature at most one dopant ion.Because of the band structure (see Fig. 1a), 1L-TMDs can have 4e species and 2h species.This suggests that the largest stable cluster will have a positive donor ion, 4 distinguishable e, and 2 distinguishable h.We get bound-state wave functions describing the donor-bound double-negative XX (D − XX).These seven-body complexes are predicted to be stable in 1L-WS 2 , 1L-WSe 2 , 1L-MoS 2 , 1L-MoSe 2 in vacuum and air.The DMC-calculated BEs with respect to the most energetically favorable products [donor-bound negative X (D 0 X − )+free X] are in Table IV.Because the dominant decay products include an X, the BE gives the PL peak position of the D − XX complex relative to the X line, possible in samples containing donor defects.

D. Accuracy of the Rytova-Keldysh interaction
RKI arises from the approximation that the in-plane susceptibility of a material is a constant [21].Here, we investigate the RKI accuracy by using an alternative approach based on ab initio calculations for 1L-MoS 2 .Realistic dielectric functions exhibit spatial dependencies which differ from the Coulomb interaction at short range (r ≪ r * ).At long range (r ≫ r * ), in-plane screening becomes irrelevant, and all physical dielectric functions behave as the Coulomb interaction, as explained in Methods.Given that the binding of excitonic complexes occurs on length scales larger than the lattice spacings (Table I), where screening effects are most prominent (Fig. 3), an investigation into their effects on charge-carrier binding is warranted.Ref.50 parameterized a dielectric permittivity ǫ(q) for 1L-MoS 2 via RPA applied to Kohn-Sham orbitals from density functional theory calculations to study charged defects.We refer to the real-space interaction formed from ǫ(q) as the RPA interaction (RPAI), and compare it to RKI in Fig. 3.
We use cusp conditions [18] to prevent the wave function of charge carriers to diverge around particle coalescence points [26].We use the same trial-wave-function form as our calculations with RKI, see Methods for details.As a test, Fig. 4   Fig. 3 indicates that at distances> a 1L-MoS 2 , RPAI follows the same form as RKI (and, ultimately, Coulomb) interactions.However, at distances∼ a 1L-MoS 2 , RKI no longer overlaps RPAI, hence cannot describe the interaction between quasiparticles.Fig. 3 shows that for r ≪ 0.5a 1L-MoS 2 , RPAI reduces to an unscreened 1/r potential, while the RKI behavior is that of a logarithmic divergence [22].However, within the effective-mass approximation, we can only describe quasiparticles on length scales> a 1L-MoS 2 , as shown in Table I and Fig. 3, whose associated Bloch wave packets are localized in momentum space, with well-defined effective mass.
The RPAI BEs of charge-carrier complexes are in Table V. Removing the bare Coulomb interaction at distances< a 1L-MoS 2 is necessary to obtain results in agreement with previous experimental [73,74] and theoretical [75] works.For r c < a 1L-MoS 2 , we truncate the RPAI to a constant v(r r c ) = v(r c ).The precise value of r c is not particularly important for the BE calculation of charge-carrier complexes, as we observe a weak BE dependence on this parameter, see Table V.Table V indicates that there is no need to use an expression for the electrostatic interaction between charge carriers in LSMs more sophisticated than RKI when evaluating BEs of trions, biexcitons, and quintons.As explained in Methods, any errors in the Mott-Wannier-Keldysh model of charge-carrier complexes for an isolated 1L are either due to the parameters (effective masses, r * , environment permittivity), or to a more fundamental breakdown of the effective-mass approximation.Intervalley scattering may play an important role in the complexes' BEs [50], while exchange effects could be relevant in highly localized complexes [76].

E. Complexes in uniform magnetic fields
For an out-of-plane external magnetic field of flux density B = (0, 0, B), where B is a positive constant, we can write the Hamiltonian as: where A i = −r i × B/2 = (−y i , x i , 0)B/2 is the magnetic vector potential for particle i in the Coulomb gauge (so that ∇ i • A i = 0) [58].We neglect the charge carriers' intrinsic magnetic dipole moment energy in the external magnetic field, because this contribution cancels out.Substituting A i into Eq.10, the term ) provides a quadratic confining potential for the particles in the complex.This cannot be regarded as a perturbation for the (otherwise free) center-of-mass (CoM) motion, because there is a quantitative difference between a bound state wave function in a quadratic potential and free motion in zero potential, no matter how small the quadratic coefficient [28].The zero-point energy of the CoM motion in the confining potential results in a linear [O(B)] contribution to the total energy, as given in Eq.12.The term also weakly perturbs the relative motion within the complex, giving a quadratic [O(B 2 )] contribution to the energy.We thus include the ) term in our QMC calculations.The linear (i q i /m i )A i • ∇ i term in Eq.10 breaks time-reversal symmetry as it is imaginary [77].It only adds to the energy in second-order perturbation theory, giving another O(B 2 ) contribution.This vanishes when we use a variational Ansatz consisting of a real trial wave function.We therefore neglect it.
The ground-state energies of isolated e/h are E e = eB/(2m e ) and E h = eB/(2m h ), in the presence of a magnetic field [58].More generally, if a bound complex of N e e and N h h moves in a magnetic field, from Eq.10 the quadratic confining potential is: where R is the CoM position.The total mass of the complex is N e m e + N h m h .Hence, we obtain the CoM zero-point energy of a charge complex as: If m e = m h ≡ m then E CoM = eB/(2m), independent of N e , N h .For a bound complex, our results show that the magnetic field can always be made sufficiently weak so that the external potential is slowly varying on the length scale of the complex (i.e./(eB) < a * 0 ).Hence, Eq.12 is the leading-order contribution to the free chargecarrier complex energy in a magnetic field.(f) FIG. 5. Theoretical BEs of (a) X, (c) X − , and (e) XX as a function of perpendicular magnetic field for 1L-WSe2 in vacuum.We use the ab initio mass and r * parameters of Table I Experimental BEs of (b) X, (d) X − , (f) XX for 1L-WSe2 encapsulated in hBN, compared with DMC ones using ǫ = ǫ0 and r * = 48 Å and the fit to Eq.21 Fig. 5 plots the DMC X, X − , XX BEs for 1L-WSe 2 in vacuum, in the presence of an out-of-plane magnetic field, using RKI.m e , m h , and r * are taken from Table I.Our results are in agreement with Ref. 78.The CoM contribution of Eq.12 is a good approximation to calculate the X, X − , XX BEs in magnetic fields< 8T, because it is exact up to linear order in magnetic field within the effective-mass approximation.For the X BE in magnetic fields> 8T, we use Eq.21, derived in Methods.The fitted C in Eq.21 is 0.557 for X in 1L-WSe 2 .
Figs5b,d,f compare our DMC BEs with measurements for hBN-encapsulated 1L-WSe 2 .The sample is produced by exfoliating flux zone grown B-WSe 2 [79], then encapsulating it with ML-hBN (10nm bottom and 3nm top) using an all-dry technique [80,81].Measurements are done in a closed-cycle cryostat (Attocube Attodry 1000) at 4K with superconducting magnets allowing out-of-plane magnetic fields up to 8T.CW excitation is provided with a diode laser at 658nm, close to the 1L-WSe 2 optical band gap [82].Polarization-resolved excitation and collection pass through a confocal microscope with the sample in reflection geometry.The PL signal is sent to a liquid-N 2cooled spectrometer (Princeton).We assume r * = 48 Å and ǫ = ǫ 0 , as discussed in Sec.II B. The theoretical and experimental BEs differ< 0.3meV over the 0-8K temperature range.The O(B) magnetic-field dependence is only via the effective masses and N e , N h , via the CoM energy, Eq.12.The fact that the theoretical and experimental magnetic-field trends in Fig. 5 agree well demonstrates that the approximation with ab initio effective masses is accurate.The main challenge is to obtain a sufficiently accurate interaction between charge carriers.The BE O(B) term is the same for all complexes, in the limit m e = m h .For most 1L-TMDs, m e and m h are similar, Table I, implying that the magnetic-field dependence cannot be used to distinguish carrier complexes.Table VI has DMC and experimental X, X − , XX BEs for 1L-WSe 2 in the presence of an out-ofplane external magnetic field, as for Fig. 5.The variation of BEs of different charge complexes is the same< 8T.

F. Complexes in uniform electric fields
A bias voltage ∆V applied to a 1L-LSM results in an in-plane electric field.Its precise form depends on device geometry.Here, we assume a uniform electric field F = −∆V /d, where d is the distance between terminals, for simplicity.F will perturb the energies of chargecarrier complexes in the CoM frame.We therefore investigate the effects of F on BEs by including an additional term − i q i F x i in the Hamiltonian, where x i is the x coordinate of particle i. Fig. 6 plots the X BE shift as a function of electric field strengths for 1L-MoS 2 , 1L-MoSe 2 , 1L-WS 2 , 1L-WSe 2 , in vacuum and encapsulated by hBN, using the ab initio parameters in Table I and ǫ = 4ǫ 0 .In each case, X BE goes as the square of the inplane electric field, as expected for a linearly polarizable exciton [83].Thus, the total energy of an isolated neutral complex of polarizability α in a uniform F is: where E F =0 is the energy of the complex in the absence of external fields.The variation of energy with electric field strength remains quadratic up to at least∼ 50mV nm −1 .> 50 mV nm −1 we find that optimizing wave functions by VMC energy minimization does not result in bound-state wave functions.If the parameters in the wave function are fixed such that a bound state is forced, the resulting DMC calculations are unstable.It is possible that some, or all, complexes remain bound at these larger electric fields, and our QMC calculations become unstable, due to the choice of trial wave function.The form we use is isotropic, so it does not allow the complex to polarize in VMC.Polarization arises at DMC level.The XX and donor-atom BEs vary linearly with F 2 , Fig. 6.However, while the donor-atom BEs increase with F 2 , the XX BEs decrease.For a 4-particle complex, alignment of charges in the direction of the applied field places like charges closer together, and reduces BE with respect to dissociation into two-particle complexes.Trion BEs also vary linearly with F 2 , Fig. 7.However, QMC calculations become unstable at much lower F .This is reflected in the higher polarizabilities for trions than for neutral complexes, Table VII.
The predicted BE shifts of each of the complexes are in Table VIII for 1L-TMDs, both in vacuum and encapsulated by hBN, subject to F = 50mV nm −1 , beyond which VMC energy minimization does not result in bound-state wave functions.The shifts in the peaks of the trions are so large that, at the very least, they should be experimentally distinguished from the neutral complexes when an electric field is applied.Identification of a positive from a negative trion may be possible in some materials/environments, but not all.For neutral complexes, the differences of a few meV in BE shifts suggest they are unlikely to be experimentally identified by their peak shifts under an electric field.

III. CONCLUSIONS
We used DMC to calculate XX − BEs in 1L-LSMs within the effective-mass approximation, using the RKI potential.A program available online [60] can be used to evaluate interpolated XX − BEs given e and h effective masses, in-plane susceptibility, and environment permittivity for a desired 1L-LSM.The BEs of charge-carrier complexes in 1L-LSMs in vacuum from RKI are in excellent agreement with those obtained using interaction potentials taken from ab initio RPA, suggesting RKI is a reliable interaction potential to describe screened interaction between charge carriers in 1L-LSMs.
We also considered the effect of external out-of-plane magnetic fields and in-plane electric fields on BEs of charge-carrier complexes in 1L-LSMs.The resulting BE changes are linear in magnetic fields and quadratic in electric fields up to 10T and 50mV nm −1 .
We measured X, X − , XX BEs for hBN-encapsulated 1L-WSe 2 up to 8T, where the BEs vary linearly with magnetic field, and found them to be in good agreement with the effective-mass approximation using ab initio effective masses.These BE shifts could in principle be used to identify complexes in PL experiments, provided m * e and m * h are different.In practice, m * e and m * h in 1L-TMDs are too similar to distinguish complexes in external magnetic fields.In-plane electric fields should shift the BE peaks in proportion to the field strength and allow for identification of charged from neutral complexes.
We derived BEs of charge-carrier complexes in 1L-TMDs by solving the interacting quantum few-body problem for each complex, working within the effectivemass approximation, with a RKI potential between charge carriers.The BE magnetic-field dependence agrees with experiments on a sub-meV energy scale.Since this only involves m * e and m * h , and not the parameters describing the screened interaction, the approximation with ab initio effective masses is highly accurate.
Efforts to improve the quantitative accuracy of BE calculations should therefore focus on the description of substrate and environmental screening, and on the inclusion of contact interactions and intervalley scattering.

A. Effective-mass approximation
All our calculations are performed within the effectivemass approximation.For charge-carrier complexes in 1L-LSMs in the absence of external fields, we solve the Mott-Wannier-Keldysh Schrödinger equation [26]: 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, E is the energy eigenvalue, ǫ is the absolute permittivity of the surrounding medium, and r * ≡ κ/(2ǫ), where κ is the inplane susceptibility.In Eq.14 the electrostatic interaction potential V , known as RKI [21,22,84], is given by [26]: where H n (x) is a Struve function [85] and Y n (x) is a Bessel function of the second kind [85].At long range (r ≫ r * ) the potential in Eq.15 is a Coulomb interaction V (r/r * ) ∼ r * /r; at short range (r ≪ r * ), logarithmic: where γ ∼ 0.57721 is Euler's constant [85].We do not include contact interactions between charge carriers due to exchange and correlation effects that occur when they are localized on the same site [76], since these partially cancel out of BEs for complexes larger than X.

B. QMC calculations
We use VMC [72] and DMC [29,86] to calculate the total energies of complexes of charge carriers in 1L-LSMs.We use the RKI potential in Eq.15 or, for the short range (r ≪ r * ) limit, the logarithmic interaction of Eq.16.Our trial wave functions for complexes of distinguishable charge carriers are of the Jastrow form [18], which includes a pairwise sum of terms depending on the distances between charge carriers, as for Ref. 26.Trial wave functions are optimized within VMC by minimizing first the energy variance [87,88], then the energy expectation [72].Our fixed-node DMC energies are exact solutions to the Mott-Wannier-Keldysh model of Eq.14.DMC calculations use time steps in the ratio 1 : 4, with the corresponding target configuration populations in the ratio 4 : 1.The resulting energies are extrapolated linearly to zero time step and to infinite population.QMC calculations are done in the casino code [18].
C. Fitting function for BE as a function of magnetic field We consider a complex of N e and N h e and h interacting via the logarithmic approximation to the Keldysh interaction in the presence of a uniform magnetic field B = (0, 0, B).Let B = B/B 0 , mi = m i /µ, qi = q i /e, ri = r i /( √ 2r 0 ), and r * = r * /( √ 2r 0 ) be magnetic field, mass, charge, and position of particle i.The screening length and the Hamiltonian Ĥ = Ĥ/E 0 in logarithmic e.u. are as defined in Sec.II A. We thus get: where we neglect the term (i q i /m i )A i • ∇ i in Eq.10 that breaks time-reversal symmetry.The energy eigenvalue Ẽ = E/E 0 is therefore the sum of a function f (σ, B), where σ = m e /m h , and an additive constant c(r * ) = i>j qi qj ln (r * ).For X in the absence of an external magnetic field, ẼX B=0 = 0.41057747491(7) − ln( √ 2) − ln(r * ) was calculated in Ref. 26.For B such that the magnetic confinement energy is larger than the log interaction, the interaction − i>j qi qj ln(e γ rij /2) is negligible compared with the magnetic confinement energy of each particle.The dimensionless total energy is the sum of the zero-point energies of the individual particles in the quadratic potential plus the constant c(r * ).Hence, at large B ≫ 1: since Ẽ B=0 ∼ c(r * ), when r * is large (r * ≫ 1).For small B ≪ 1, we use the CoM zero-point energy approximation, Eq.12, in which we assume the quadratic potential varies on the scale of the complex.Then: The total energies for X with m e = m h are calculated using the finite-element method (FEM) implemented in Mathematica [89].The results are converged by increasing the region size and decreasing the maximum cell size in order to achieve at least six digits of precision.This leads to errors comparable errors to QMC (see Sec.IV B).Subtracting the large-B, Eq.19, from the energy shift of X due to external magnetic fields, results in the logarithmic-like behavior in Fig. 8.This suggests the following formula for the energy shift of a generic charge-carrier complex due to external magnetic field: where C = C(r * , σ) is independent of r * in the r ≪ r * limit in which the logarithmic interaction is valid.We use the least-squares method to fit the FEM results for the BE shift of X with equal m * e and m * h , Fig. 8, for several values of susceptibility, and extract the fitting parameter C X for each r * .We use a polynomial fit to get the dependence of C X on susceptibility, Fig. 8: C X = −0.1020(22)+ 0.546 (9)x + 0.194(6)x 2 , ( where x = r * / (1 + r * ).Since most 1L-TMDs have effective mass ratios close to 1, Table I, we neglect the massratio dependence of C. Our fit is only valid for r * 0.5.Although the fit from Eq.21 is derived for the logarithmic interaction, it does fit well our DMC results in Fig. 5 for the full Keldysh interaction for experimentally relevant values, as shown by the red curves in Fig. 5.

FIG. 1 .
FIG. 1.(a) Upper spin-split VB and spin-split CB for 1L-MoSe2 and 1L-WSe2.The spin-split VB is> 150meV[5], so we only show the upper VB. (b,c) Classification of quinton recombination processes in 1L-Mo and 1L-W-TMDs.E δ = E XX − − E X − is the difference between the total energies E XX − and E X − of XX − and X − .ω indicates the photon energies at which XX − peaks in PL spectra are expected.(XX − ) k 1 σ 1 k 2 σ 2 k 3 σ 3 k 4 σ 4 k 5 σ 5 denotes a quinton consisting of CB e in valleys k1, k2, and k3 with spins σ1, σ2, σ3 and VB h in valleys k4 and k5, with spins σ4, σ5.E.g., the quintons in (a) are both (XX − ) K↑K↓K ′ ↓ K↑K ′ ↓ .Unlike Figs.1,2 of Ref.26, we only show complexes with distinguishable charge carriers, because they are stable and should be experimentally observable.

Fig. 1
Fig.1 classifies XX − in 1L-Mo-and W-TMDs with respect to recombination energy and T -dependence of the emitted photons' intensity.There are two XX − types:(1) those with 1e in the upper spin-split CB and 2e in the lower spin-split CB, and (2) those with 1e in the upper spin-split CB and 1e in the lower spin-split CB.The CB spin splittings in 1L-Mo-and W-TMDs are∼ 3meV[5] and∼ 30meV[5], respectively.These are much less than the XX − BEs∼ 50meV[30,34,47], as reported in Table I.The fact that the XX − BE is larger than the spin splitting implies XX − complexes are thermodynamically stable at T close to 0K, even taking into account the energy required to excite 1e to the upper spin-split CB.Assuming the CB spin-orbit splitting ∆ ′ of 1L-TMDs to be≪XX − BE, E b XX − , each XX − can be treated as a twostate system[48].For k B T ≪ ∆ ′ ≪ E b XX − , with k B the Boltzmann's constant, the fraction of XX − with 1e and 2e in the upper spin-split CB, hence the PL intensity of the corresponding XX − , is[49]:

FIG. 6 .
FIG. 6. DMC BE shift for (a) X, (b) XX, (c) donor atoms as a function of F 2 for different 1L-TMDs in vacuum and encapsulated in hBN.Error bars in (a,c) are smaller than the symbols.The solid and dashed lines are BEs determined by the polarizabilities in Table VII for 1L-TMDs in vacuum and encapsulated by hBN.The vertical dotted lines correspond to F = 50mV nm −1 , beyond which VMC energy minimization does not result in bound-state wave functions.

FIG. 7 .
FIG. 7. DMC BE shift for (a) X − and (b) X + as a function of F 2 for different 1L-TMDs in vacuum and encapsulated in hBN.Where error bars are not visible they are smaller than the symbols.The solid and dashed lines show BEs determined by the polarizabilities in Table VII (r * ) ,

FIG. 8 .
FIG. 8. (a) Dependence of CX on susceptibility.The markers show fitted CX.The line is a quadratic fit to the points as for Eq.22.(b)Shift in energy of X with equal m * e and m * h due to external magnetic field, after subtracting the large-B behavior from Eq.19, for several r * .Markers indicate the finite-element method results, while lines show the fit of Eq.21.

TABLE II .
List of acronymsSince the most stable dissociated complexes have the lowest ground-state energies, the XX − BE is the minimum of E DE XX − and E EA XX for a given r * and effective mass: 30sts the XX − BEs from DMC and the fit of Eq.8 for 1L-WS2 , 1L-WSe 2 , 1L-WTe 2 , 1L-MoS 2 , 1L-MoSe 2 , 1L-MoTe 2 .To measure X − , XX, XX − BEs, in Ref.30we used continuous wave (CW) PL at 4K for 1L-WSe 2 encapsulated between 2 10nm (bottom) and 3nm (top) ML-hBN on Si/SiO 2 .Refs.30, 35, and 47 used different experimental conditions for 1L-WSe 2 , but they all produced similar results, Table I, across various techniques and substrates.

TABLE IV .
Theoretical E b D − XX for 1L-TMDs in vacuum, with ab initio masses and r * from TableI

TABLE V .
BEs of charge-carrier complexes in 1L-MoS2 calculated using different interaction potentials.R-RPA is the rounded RPAI, with rc in brackets.RKI values from Ref.26.

TABLE VI .
DMC and experimental X, X − , XX BEs in meV for 1L-WSe2 with an out-of-plane external magnetic field.

TABLE VII .
Theoretical in-plane polarizabilities of X, XX, D 0 , X − , X + in 1L-TMDs, in vacuum and hBN encapsulated

TABLE VIII .
Calculated BE shifts of X, XX, D 0 , X − , X + using Eq.13 and polarizabilities in Table VII for 1L-TMDs, both in vacuum and encapsulated by hBN, for F = 50mV nm −1 .Not all complexes are bound at F = 50 mV nm−1 me + N h / mh N e me + N h mh −