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

TL;DR

This study uses diffusion quantum Monte Carlo to accurately compute binding energies of excitonic complexes in 2D transition-metal dichalcogenides, providing insights into their stability and spectral signatures.

Contribution

It offers statistically exact binding energy data and classification of excitonic complexes in 2D semiconductors, enhancing interpretation of experimental spectra.

Findings

Biexcitons are only bound when identical particles are much heavier.

Donor/acceptor-bound biexcitons have energies similar to experimental values.

Predicted spectral positions of excitonic complexes in WSe₂ and MoSe₂.

Abstract

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$.

Tables

Table 1. Table 1: Classification of exciton states into irreps of D 3 ​ h ′′ superscript subscript 𝐷 3 ℎ ′′ D_{3h}^{\prime\prime} and the polarization ( ∥ parallel-to \parallel and z 𝑧 z for in-plane and out-of-plane, respectively) of the electric field to which the excitons are coupled.

Irrep	Excitons	Field
$E^{+}$	$X_{K\uparrow }^{K\uparrow }$, $X_{K^{\prime }\downarrow }^{K^{\prime }\downarrow }$	${𝐄}_{\parallel }$
$E^{+}$	$X_{K\downarrow }^{K\uparrow }$, $X_{K^{\prime }\uparrow }^{K^{\prime }\downarrow }$
$A_2^{-}$	X (dark): $\left\{X_{K\uparrow }^{K\downarrow }-X_{K^{\prime }\downarrow }^{K^{\prime }\uparrow }\right\}$	${𝐄}_z$
$E^{-}$	$X_{K\uparrow }^{K\downarrow }$, $X_{K^{\prime }\downarrow }^{K^{\prime }\uparrow }$	–
$A_1^{-}$	$\left\{X_{K\downarrow }^{K\uparrow }+X_{K^{\prime }\uparrow }^{K^{\prime }\downarrow }\right\}$

Table 2. Table 2: Effect of changing the number of free parameters in η 𝜂 \eta [i.e., the values of N η hh superscript subscript 𝑁 𝜂 hh N_{\eta}^{\rm hh} and N η eh superscript subscript 𝑁 𝜂 eh N_{\eta}^{\rm eh} in Eqs. ( 42 ) and ( 43 )] on the VMC ground-state energy ( E VMC subscript 𝐸 VMC E_{\rm VMC} ), VMC energy variance, and DMC energy ( E DMC subscript 𝐸 DMC E_{\rm DMC} ) for biexcitons with indistinguishable holes interacting via the logarithmic interaction. The mass ratio is σ = 0.1 𝜎 0.1 \sigma=0.1 and the reduced mass is μ = 0.5 ​ m 0 𝜇 0.5 subscript 𝑚 0 \mu=0.5m_{0} , where m 0 subscript 𝑚 0 m_{0} is the bare electron mass. In each case r ∗ = r 0 subscript 𝑟 ∗ subscript 𝑟 0 r_{\ast}=r_{0} .

$N_{\eta }^{\mathrm{hh}}=N_{\eta }^{\mathrm{eh}}$	$E_{\mathrm{VMC}}$ ($E_0$)		Variance ($E_0^2$)		$E_{\mathrm{DMC}}$ ($E_0$)
2	$0.$	$7604\left(3\right)$	$0.$	$00920$	$0.$	$7585\left(2\right)$
3	$0.$	$7602\left(3\right)$	$0.$	$00908$	$0.$	$7584\left(2\right)$
4	$0.$	$7605\left(3\right)$	$0.$	$00914$	$0.$	$7579\left(2\right)$
6	$0.$	$7606\left(3\right)$	$0.$	$00927$	$0.$	$7580\left(2\right)$

Table 3. Table 3: Effect of changing the number of free parameters in η 𝜂 \eta [i.e., the values of N η hh superscript subscript 𝑁 𝜂 hh N_{\eta}^{\rm hh} and N η eh superscript subscript 𝑁 𝜂 eh N_{\eta}^{\rm eh} in Eqs. ( 42 ) and ( 43 )] on the VMC ground-state energy ( E VMC subscript 𝐸 VMC E_{\rm VMC} ), VMC energy variance, and DMC energy ( E DMC subscript 𝐸 DMC E_{\rm DMC} ) of biexcitons with indistinguishable holes interacting via the Keldysh interaction, with an electron-hole mass ratio of σ = 0.1 𝜎 0.1 \sigma=0.1 . r ∗ = 0 subscript 𝑟 ∗ 0 r_{\ast}=0 corresponds to the Coulomb interaction.

$r_{\ast }/a_0^{\ast }$	$N_{\eta }^{\mathrm{hh}}=N_{\eta }^{\mathrm{eh}}$	$E_{\mathrm{VMC}}$ ($R_y^{\ast }$)		Variance ($R_y^{\ast }{}_{}{}^2$)		$E_{\mathrm{DMC}}$ ($R_y^{\ast }$)
0	2	$-8.$	$608\left(1\right)$	$0.$	$1709$	$-8.$	$6100\left(4\right)$
0	3	$-8.$	$608\left(1\right)$	$0.$	$1658$	$-8.$	$6112\left(4\right)$
0	4	$-8.$	$606\left(1\right)$	$0.$	$1718$	$-8.$	$6108\left(4\right)$
0	6	$-8.$	$608\left(1\right)$	$0.$	$1719$	$-8.$	$6108\left(4\right)$
8	2	$-0.$	$6304\left(1\right)$	$0.$	$0008$	$-0.$	$6308\left(1\right)$
8	3	$-0.$	$63020\left(8\right)$	$0.$	$0007$	$-0.$	$6306\left(2\right)$
8	4	$-0.$	$6301\left(1\right)$	$0.$	$0007$	$-0.$	$6308\left(1\right)$
8	6	$-0.$	$63024\left(8\right)$	$0.$	$0007$	$-0.$	$6309\left(1\right)$

Table 4. Table 4: 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 italic-ϵ subscript italic-ϵ 0 \epsilon=\epsilon_{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 m e subscript 𝑚 e m_{\rm e} , m h subscript 𝑚 h m_{\rm h} , and r ∗ subscript 𝑟 ∗ r_{\ast} shown in bold in Table 5 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.

	Ex. tot. en. $E_X$ (meV)		Biexciton binding energy $E_{\mathrm{XX}}^b$ (meV)
TMDC	Eq. (47)	VBerkelbach et al. (2013)	Eq. (48)	DMCMayers et al. (2015)	PIMCKylänpää and Komsa (2015)	Exp.	SVZhang et al. (2015a); Kidd et al. (2016)	HHKezerashvili and Tsiklauri (2016)
MoS2	$-546.5$	$-540$	$23.5$	$22.7\left(3\right)$	$22.7\left(5\right)$	$70$,Mai et al. (2014a) $40$–$60$Sie et al. (2015)	$22.5$	$22.1$
MoSe2	$-543.3$	$-470$	$22.9$	$17.7\left(3\right)$	$19.3\left(5\right)$	$\simeq 20$Hao et al. (2017)	$18.4$	$17.9$
MoTe2	$-385.7$		$13.4$		$14.4\left(4\right)$
WS2	$-511.9$	$-500$	$23.4$	$23.3\left(3\right)$	$23.9\left(5\right)$	$65$,Plechinger et al. (2015) $65$–$70$,Paradisanos et al. (2017) $45$Shang et al. (2015)	$23.6$	$23.1$
WSe2	$-461.1$	$-450$	$20.0$	$20.2\left(3\right)$	$20.7\left(5\right)$	$52$You et al. (2015)	$20.2$	$19.8$
WTe2	$-454.3$		$18.9$

Table 5. Table 5: Comparison of electron and hole effective masses and r ∗ subscript 𝑟 ∗ r_{\ast} values obtained by different ab initio methods in the literature [many-body G ​ W 𝐺 𝑊 GW calculations, and density functional theory either in the local density approximation (LDA) or using the Perdew-Burke-Ernzerhof (PBE) or Heyd-Scuseria-Ernzerhof (HSE06) exchange-correlation functionals]. Where multiple results are available using a given method, we have taken the average of the published results. The reported r ∗ subscript 𝑟 ∗ r_{\ast} values and effective masses are not necessarily obtained using the same method. We assume that the materials are suspended in vacuum, i.e., that ϵ = ϵ 0 italic-ϵ subscript italic-ϵ 0 \epsilon=\epsilon_{0} . The effective masses are reported in units of the bare electron mass m 0 subscript 𝑚 0 m_{0} . The binding energies E XX b subscript superscript 𝐸 b XX E^{\rm b}_{\rm XX} are calculated using Eq. ( 48 ). The effective masses and r ∗ subscript 𝑟 ∗ r_{\ast} values shown in bold are used to evaluate the binding energies reported in Tables 4 and 7 .

	Effective masses
TMDC	Method	$m_e/m_0$		$m_h/m_0$		$r_{\ast }$ (Å)		$E_{\mathrm{XX}}^b$ (meV)
MoS2	$GW$Cheiwchanchamnangij and Lambrecht (2012)	$\mathrm{𝟎}$	$.35$	$\mathrm{𝟎}$	$.428$	$\mathrm{𝟑𝟖}$	$.62$Cheiwchanchamnangij and Lambrecht (2012)	$23$	$.5$
	$G_0{W}_0$Ramasubramaniam (2012)	$0$	$.60$	$0$	$.54$	$38$	$.62$Cheiwchanchamnangij and Lambrecht (2012)	$24$	$.2$
	LDAKadantsev and Hawrylak (2012); Jin et al. (2014); Chang et al. (2014); Scalise et al. (2014); Yu et al. (2015); Kormányos et al. (2015)	$0$	$.495$	$0$	$.576$	$36$	$.28$Kumar and Ahluwalia (2012)	$25$	$.5$
	PBEWickramaratne et al. (2014); Dong et al. (2014); Yun et al. (2012); Kormányos et al. (2015); Kylänpää and Komsa (2015)	$0$	$.470$	$0$	$.575$	$44$	$.69$Kylänpää and Komsa (2015)	$21$	$.1$
	HSE06Peelaers and Van de Walle (2012)	$0$	$.37$	$0$	$.44$	$38$	$.62$Cheiwchanchamnangij and Lambrecht (2012)	$23$	$.6$
MoSe2	$GW$Shi et al. (2013)	$\mathrm{𝟎}$	$.38$	$\mathrm{𝟎}$	$.44$	$51$	$.71$Berkelbach et al. (2013)	$18$	$.0$
	$G_0{W}_0$Ramasubramaniam (2012)	$0$	$.70$	$0$	$.55$	$51$	$.71$Berkelbach et al. (2013)	$18$	$.7$
	LDAJin et al. (2014); Chang et al. (2014); Kormányos et al. (2015)	$0$	$.59$	$0$	$.686$	$\mathrm{𝟑𝟗}$	$.79$Kumar and Ahluwalia (2012)	$23$	$.7$
	PBEWickramaratne et al. (2014); Liu et al. (2011); Kormányos et al. (2015); Kylänpää and Komsa (2015)	$0$	$.546$	$0$	$.643$	$53$	$.16$Kylänpää and Komsa (2015)	$18$	$.1$
MoTe2	$G_0{W}_0$Ramasubramaniam (2012)	$\mathrm{𝟎}$	$.69$	$\mathrm{𝟎}$	$.66$	$\mathrm{𝟕𝟑}$	$.61$Kylänpää and Komsa (2015)	$13$	$.4$
	LDAChang et al. (2014)	$0$	$.64$	$0$	$.78$	$73$	$.61$Kylänpää and Komsa (2015)	$15$	$.5$
	PBELiu et al. (2011)	$0$	$.575$	$0$	$.702$	$73$	$.61$Kylänpää and Komsa (2015)	$13$	$.4$
WS2	$GW$Shi et al. (2013)	$\mathrm{𝟎}$	$.27$	$\mathrm{𝟎}$	$.32$	$\mathrm{𝟑𝟕}$	$.89$Berkelbach et al. (2013)	$23$	$.4$
	$G_0{W}_0$Ramasubramaniam (2012)	$0$	$.44$	$0$	$.45$	$37$	$.89$Berkelbach et al. (2013)	$24$	$.1$
	LDAJin et al. (2014); Chang et al. (2014); Kormányos et al. (2015)	$0$	$.312$	$0$	$.422$	$32$	$.42$Kumar and Ahluwalia (2012)	$27$	$.7$
	PBEWickramaratne et al. (2014); Liu et al. (2011); Kormányos et al. (2015); Kylänpää and Komsa (2015)	$0$	$.328$	$0$	$.402$	$40$	$.17$Kylänpää and Komsa (2015)	$22$	$.6$
WSe2	$GW$Shi et al. (2013)	$\mathrm{𝟎}$	$.29$	$\mathrm{𝟎}$	$.34$	$\mathrm{𝟒𝟓}$	$.11$Berkelbach et al. (2013)	$20$	$.0$
	$G_0{W}_0$Ramasubramaniam (2012)	$0$	$.53$	$0$	$.52$	$45$	$.11$Berkelbach et al. (2013)	$20$	$.8$
	LDAJin et al. (2014); Chang et al. (2014); Kormányos et al. (2015)	$0$	$.36$	$0$	$.476$	$34$	$.72$Kumar and Ahluwalia (2012)	$26$	$.2$
	PBEWickramaratne et al. (2014); Kormányos et al. (2015); Kylänpää and Komsa (2015)	$0$	$.342$	$0$	$.428$	$47$	$.57$Kylänpää and Komsa (2015)	$19$	$.4$
WTe2	LDAKormányos et al. (2015)	$\mathrm{𝟎}$	$.325$	$\mathrm{𝟎}$	$.460$	$\mathrm{𝟒𝟗}$	$.56$Kumar and Ahluwalia (2012); Lu et al. (2016)	$18$	$.9$
	PBELv et al. (2015)	$0$	$.307$	$0$	$.51$	$49$	$.56$Kumar and Ahluwalia (2012); Lu et al. (2016)	$19$	$.3$

Table 6. Table 6: 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 5 . m 0 subscript 𝑚 0 m_{0} is the bare electron mass.

TMDC	$\frac{\partial E_X}{\partial m_e}$		$\frac{\partial E_X}{\partial m_h}$		$\frac{\partial E_X}{\partial r_{\ast }}$		$\frac{\partial E_{\mathrm{XX}}^b}{\partial m_e}$		$\frac{\partial E_{\mathrm{XX}}^b}{\partial m_h}$		$\frac{\partial E_{\mathrm{XX}}^b}{\partial r_{\ast }}$		$\frac{\partial E_{X^{-}}^b}{\partial m_e}$		$\frac{\partial E_{X^{-}}^b}{\partial m_h}$		$\frac{\partial E_{X^{-}}^b}{\partial r_{\ast }}$
	(meV/$m_0$)		(meV/$m_0$)		(meV/Å)		(meV/$m_0$)		(meV/$m_0$)		(meV/Å)		(meV/$m_0$)		(meV/$m_0$)		(meV/Å)
MoS2	$-240$		$-160$		$10$		$1$	$.2$	$5$	$.4$	$-0$	$.56$	$10$		$9$	$.2$	$-0$	$.72$
MoSe2	$-210$		$-160$		$9$	$.9$	$1$	$.4$	$4$	$.3$	$-0$	$.53$	$8$	$.8$	$8$	$.7$	$-0$	$.70$
WS2	$-310$		$-220$		$9$	$.5$	$2$	$.1$	$6$	$.7$	$-0$	$.56$	$13$		$12$		$-0$	$.70$
WSe2	$-240$		$-180$		$7$	$.3$	$1$	$.6$	$5$	$.1$	$-0$	$.41$	$10$		$9$	$.6$	$-0$	$.52$

Table 7. Table 7: As Table 4 (using the r ∗ subscript 𝑟 ∗ r_{\ast} values and effective masses shown in bold in Table 5 ), 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), acceptor-bound 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.3 .

	Energy (meV)		Binding energy (meV)
TMDC	$E_{D^0}$	$E_{A^0}$	$E_{D^{+}X}^b$	$E_{A^{-}X}^b$	$E_{D^{0}X}^b$	$E_{A^{0}X}^b$	$E_{D^{+}\mathrm{XX}}^b$	$E_{A^{-}\mathrm{XX}}^b$
MoS2	$-638.5$	$-670.1$	$7.2$	$2.7$	$32.4$	$31.7$	$51.8$	$48.0$
MoSe2	$-636.8$	$-659.3$	$6.5$	$3.2$	$31.5$	$31.1$	$50.8$	$48.1$
MoTe2	$-447.8$	$-443.9$	$3.9$	$4.6$	$17.9$	$18.0$	$32.4$	$32.9$
WS2	$-606.8$	$-633.6$	$6.2$	$2.7$	$32.0$	$31.5$	$48.6$	$45.3$
WSe2	$-542.3$	$-563.6$	$5.5$	$2.6$	$27.5$	$27.1$	$43.4$	$40.8$
WTe2	$-519.0$	$-562.1$	$7.9$	$1.6$	$26.1$	$25.3$	$44.0$	$38.8$

Table 8. Table 8: Binding energies of larger charge-carrier complexes in different TMDCs. Binding energies are presented for: quintons (XX - ), which dissociate into excitons (X) and negative trions (X - ); donor-bound double-negative excitons (D - X), which dissociate into excitons (X) and negative donor ions (D - ); and donor-bound quintons (D 0 XX), which dissociate into excitons (X) and donor-bound trions (D 0 X). The binding energies were evaluated using the effective mass and in-plane permittivity parameters reported in bold in Table 5 . The fitting functions of Eqs. ( 49 ) and ( 51 ) were used to evaluate the energies of negative donor ions and donor-bound trions.

	Binding energy (meV)
TMDC	XX-	D-X	D0XX
MoS2	$58.6\left(6\right)$	$84.4\left(4\right)$	$61.6\left(6\right)$
MoSe2	$57.0\left(4\right)$	$57.9\left(2\right)$	$56.9\left(9\right)$
MoTe2	$33.8\left(3\right)$
WS2	$57.4\left(3\right)$	$59.2\left(4\right)$	$58.2\left(6\right)$
WSe2	$52.5\left(7\right)$	$51.3\left(4\right)$	$51\left(1\right)$
WTe2	$47.5\left(3\right)$

Table 9. Table 9: Character table for the irreps of the extended point group D 3 ​ h ′′ superscript subscript 𝐷 3 ℎ ′′ D_{3h}^{\prime\prime} . The superscript ( ± plus-or-minus \pm ) denotes the transformation under the σ h subscript 𝜎 ℎ \sigma_{h} operation and ′ denotes representations with nontrivial transformation under translation.

$D_{3h}^{\prime \prime }$	$E$	$2t$	$2C_3$	$2t{C}_3$	$2t^2{C}_3$	$9t{\sigma }_v$	${\sigma }_h$	$2t{\sigma }_h$	$2S_3$	$2t{S}_3$	$2t^2{S}_3$	$9t{C}_2^{\prime }$
$A_1^{+}$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$
$A_2^{+}$	$1$	$1$	$1$	$1$	$1$	$-1$	$1$	$1$	$1$	$1$	$1$	$-1$
$A_1^{-}$	$1$	$1$	$1$	$1$	$1$	$-1$	$-1$	$-1$	$-1$	$-1$	$-1$	$1$
$A_2^{-}$	$1$	$1$	$1$	$1$	$1$	$1$	$-1$	$-1$	$-1$	$-1$	$-1$	$-1$
$E^{+}$	$2$	$2$	$-1$	$-1$	$-1$	$0$	$2$	$2$	$-1$	$-1$	$-1$	$0$
$E^{-}$	$2$	$2$	$-1$	$-1$	$-1$	$0$	$-2$	$-2$	$1$	$1$	$1$	$0$
$E_1^{\prime +}$	$2$	$-1$	$-1$	$2$	$-1$	$0$	$2$	$-1$	$-1$	$2$	$-1$	$0$
$E_1^{\prime -}$	$2$	$-1$	$-1$	$2$	$-1$	$0$	$-2$	$1$	$1$	$-2$	$1$	$0$
$E_2^{\prime +}$	$2$	$-1$	$2$	$-1$	$-1$	$0$	$2$	$-1$	$2$	$-1$	$-1$	$0$
$E_2^{\prime -}$	$2$	$-1$	$2$	$-1$	$-1$	$0$	$-2$	$1$	$-2$	$1$	$1$	$0$
$E_3^{\prime +}$	$2$	$-1$	$-1$	$-1$	$2$	$0$	$2$	$-1$	$-1$	$-1$	$2$	$0$
$E_3^{\prime -}$	$2$	$-1$	$-1$	$-1$	$2$	$0$	$-2$	$1$	$1$	$1$	$-2$	$0$

Table 10. Table 10: Product table of the extended point group D 3 ​ h ′′ superscript subscript 𝐷 3 ℎ ′′ D_{3h}^{\prime\prime} . The relevant classification of the irreps according to C s = { I , σ h } subscript 𝐶 𝑠 𝐼 subscript 𝜎 ℎ C_{s}=\{I,\sigma_{h}\} , denoted by the superscript ± plus-or-minus \pm in Table 9 , is included for a given product by using + ⁣ ⊗ ⁣ + ⁣ = ⁣ + , + ⁣ ⊗ ⁣ − ⁣ = ⁣ − tensor-product tensor-product +\otimes+=+,+\otimes-=- and − ⁣ ⊗ ⁣ − ⁣ = ⁣ + tensor-product -\otimes-=+ and noting that all irreps in a given direct sum have the same C s subscript 𝐶 𝑠 C_{s} classification.

$\otimes$	$A_1$	$A_2$	$E$	$E_1^{\prime }$	$E_2^{\prime }$	$E_3^{\prime }$
$A_1$	$A_1$	$A_2$	$E$	$E_1^{\prime }$	$E_2^{\prime }$	$E_3^{\prime }$
$A_2$	$A_2$	$A_1$	$E$	$E_1^{\prime }$	$E_2^{\prime }$	$E_3^{\prime }$
$E$	$E$	$E$	$A_1\oplus A_2\oplus E$	$E_2^{\prime }\oplus E_3^{\prime }$	$E_1^{\prime }\oplus E_3^{\prime }$	$E_1^{\prime }\oplus E_2^{\prime }$
$E_1^{\prime }$	$E_1^{\prime }$	$E_1^{\prime }$	$E_2^{\prime }\oplus E_3^{\prime }$	$A_1\oplus A_2\oplus E_1^{\prime }$	$E\oplus E_3^{\prime }$	$E\oplus E_2^{\prime }$
$E_2^{\prime }$	$E_2^{\prime }$	$E_2^{\prime }$	$E_1^{\prime }\oplus E_3^{\prime }$	$E\oplus E_3^{\prime }$	$A_1\oplus A_2\oplus E_2^{\prime }$	$E\oplus E_1^{\prime }$
$E_3^{\prime }$	$E_3^{\prime }$	$E_3^{\prime }$	$E_1^{\prime }\oplus E_2^{\prime }$	$E\oplus E_2^{\prime }$	$E\oplus E_1^{\prime }$	$A_1\oplus A_2\oplus E_3^{\prime }$

Equations

$$I(T)\sim\left\{\begin{array}[]{ll}\color[rgb]{0,0,1}{\text{const}.}&\text{for no electrons in the upper spin-split conduction band}\\ {\color[rgb]{0,0.5,0}e^{-\Delta^{\prime}/(k_{\rm B}T)}}&\text{for one electron in the upper spin-split conduction band}\\ \color[rgb]{1,0,0}{e^{-2\Delta^{\prime}/(k_{\rm B}T)}}&\text{for two electrons in the upper spin-split conduction band}\end{array}\right.,$$

$$X=\Phi ({\mathbf{r}}_{\mathrm{e}},{\mathbf{r}}_{\mathrm{h}})\otimes U_k({\mathbf{r}}_{\mathrm{e}},{\mathbf{r}}_{\mathrm{h}})\otimes \chi (s_{\mathrm{e}},s_{\mathrm{h}}).$$

$$E_1^{\prime +}\otimes E_2^{\prime +}=E^{+}\oplus E_3^{\prime +},$$

$$D_{\text{nicefrac}12}\otimes D_{\text{nicefrac}12}=A_1^{+}\oplus A_2^{+}\oplus E^{-}.$$

$$E_{\mathrm{lattice}}^{+}\otimes (A_1^{+}\oplus A_2^{+}\oplus E^{-}{)}_{\mathrm{spin}}=2E^{+}\oplus A_1^{-}\oplus A_2^{-}\oplus E^{-}.$$

$$ϵ{\nabla }^2\varphi =-\rho \delta (z)-\kappa [{\nabla }^2\varphi (x,y,0)]\delta (z).$$

$$\varphi (\mathbf{q},k)=\frac{\rho (\mathbf{q})-\kappa q^2\varphi (\mathbf{q},z=0)}{ϵ(q^2+k^2)}.$$

$$\varphi (\mathbf{q},z=0)$$

$$\varphi (\mathbf{q},z=0)=\frac{\rho (\mathbf{q})}{q(2ϵ+q\kappa )}.$$

$$v(q)=\frac{q_i{q}_j}{2ϵq(1+r_{\ast }q)},$$

$$v(r)=\frac{q_i{q}_j}{4\pi ϵr_{\ast }}V\left(\frac{r}{r_{\ast }}\right),$$

$$V\left(r/r_{\ast }\right)=\frac{\pi }{2}\left[H_0\left(\frac{r}{r_{\ast }}\right)-Y_0\left(\frac{r}{r_{\ast }}\right)\right],$$

$$V(r/r_{\ast })\approx r_{\ast }/r,$$

$$V(r/r_{\ast })\approx \left[\log \left(2r_{\ast }/r\right)-\gamma \right]=\log \left(\frac{2r_{\ast }}{\exp (\gamma )r}\right),$$

$$V(r/r_{\ast })\approx -\log \left(\frac{r/r_{\ast }}{1+r/r_{\ast }}\right)-\left[\gamma -\log (2)\right]e^{-r/r_{\ast }}.$$

$$\left[-\sum _i\frac{{\hslash }^2}{2m_i}{\nabla }_i^2+\sum _{i>j}\frac{q_i{q}_j}{4\pi ϵr_{\ast }}V\left(\frac{r_{ij}}{r_{\ast }}\right)\right]\psi =E\psi ,$$

$$\left[-\sum _i\frac{\mu }{m_i}{\tilde{\nabla }}_i^2+\sum _{i>j}\frac{2q_i{q}_j{a}_0^{\ast }}{e^2{r}_{\ast }}V\left(\frac{{\tilde{r}}_{ij}{a}_0^{\ast }}{r_{\ast }}\right)\right]\psi =\mathcal{E}\psi ,$$

$$\left[-{\tilde{\nabla }}_{\mathrm{eh}}^2-\frac{2a_0^{\ast }}{r_{\ast }}V\left(\frac{{\tilde{r}}_{\mathrm{eh}}{a}_0^{\ast }}{r_{\ast }}\right)\right]\psi ={\mathcal{E}}_{\mathrm{X}}\psi ,$$

$$\left[-\sum _i\frac{{\hslash }^2}{2m_i}{\nabla }_i^2+\sum _{i>j}\frac{q_i{q}_j}{4\pi ϵr_{\ast }}\log \left(\frac{2r_{\ast }}{\exp (\gamma )r}\right)\right]\psi =E\psi .$$

$$r_0=\sqrt{\frac{4\pi ϵr_{\ast }{\hslash }^2}{2e^2\mu }}$$

$$E_0=\frac{e^2}{4\pi ϵr_{\ast }}.$$

$$\left\{-\sum _i\frac{\mu }{m_i}{\tilde{\nabla }}_i^2-\sum _{i>j}\frac{q_i{q}_j}{e^2}\left[\log ({\tilde{r}}_{ij})+\log \left(r_0/r_{\ast }\right)+\gamma -\log (2)\right]\right\}\psi =\mathcal{E}\psi .$$

$$\mathcal{C}=-\sum _{i>j}\frac{q_i{q}_j}{e^2}\log \left(r_0/r_{\ast }\right).$$

$$\sum_{i>j}\frac{q_{i}q_{j}}{e^{2}}=\left\{\begin{array}[]{ll}-1&\text{for an exciton or donor atom}\\ -1&\text{for a trion or donor-bound exciton}\\ -2&\text{for a biexciton or donor-bound trion}\\ -2&\text{for a donor-bound biexciton}\\ \frac{(n_{+}-n_{-})^{2}-n_{+}-n_{-}}{2}&\text{for a complex of $n_{+}$ charges $+e$ and $n_{-}$ charges $-e$}\end{array}\right..$$

$$\left[-{\tilde{\nabla }}_{\mathrm{eh}}^2+\log ({\tilde{r}}_{\mathrm{eh}})+\log \left(r_0/r_{\ast }\right)+\gamma -\log (2)\right]\psi ={\mathcal{E}}_{\mathrm{X}}\psi .$$

$${\mathcal{E}}_{\mathrm{X}}={\mathcal{E}}_{X0}-\log \left(r_{\ast }/r_0\right),$$

$$E_{{\mathrm{X}}^{-}}^{\mathrm{b}}$$

$$E_{\mathrm{XX}}^{\mathrm{b}}$$

$$E_{{\mathrm{D}}^{+}\mathrm{X}}^{\mathrm{b}}$$

$$E_{{\mathrm{D}}^0\mathrm{X}}^{\mathrm{b}}$$

$$E_{{\mathrm{D}}^{+}\mathrm{XX}}^{\mathrm{b}}$$

$$A^{\mathrm{ee}}\delta ({\mathbf{r}}_{\mathrm{ee}})+A^{\mathrm{hh}}\delta ({\mathbf{r}}_{\mathrm{hh}})+A^{\mathrm{eh}}\sum _{i=1}^2\sum _{j=1}^2\delta ({\mathbf{r}}_{{\mathrm{e}}_i{\mathrm{h}}_j}),$$

Full text

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 MoS2, MoSe2, WS2, and WSe2, which acquire a direct-gap character in hexagonal monolayer form.Mak et al. (2010); Splendiani et al. (2010); Kormányos et al. (2015); Ross et al. (2014) The direct 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,Ramasubramaniam (2012); Cheiwchanchamnangij and Lambrecht (2012); Qiu et al. (2013) including nonhydrogenic Rydberg spectraChernikov et al. (2014); Ye et al. (2014) and lines ascribed to trions (charged excitons)Ross et al. (2013); Mak et al. (2013); Lin et al. (2014) and biexcitons (bound pairs of excitons).Mai et al. (2014a); Shang et al. (2015); You et al. (2015); Plechinger et al. (2015) 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,Bressanini et al. (1998); Tan et al. (2005); Ganchev et al. (2015); Velizhanin and Saxena (2015); Mayers et al. (2015); Spink et al. (2016); Kylänpää and Komsa (2015); Chaves et al. (2016); Szyniszewski et al. (2017), variational methodsBerkelbach et al. (2013); Zhang et al. (2015a); Kidd et al. (2016), and hyperspherical harmonics approaches,Kezerashvili and Tsiklauri (2016) 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 charge-carrier 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 charge-carrier 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.

II Bright and dark biexcitons and trions in molybdenum and tungsten

dichalcogenides

II.1 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^{\prime}$ 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^{\prime}$. 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.Kormányos et al. (2015) 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 ($\Delta^{\prime}=3$–50 meV) is small enough that electrons can be found in the upper spin-split conduction band at room temperature.Kormányos et al. (2015)

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 lower-energy 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:

[TABLE]

where $\Delta^{\prime}$ is the spin-orbit-induced splitting of the conduction band, $k_{\rm 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.3.

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 electron-electron 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 $[\mu_{\rm bd}/(2\Delta^{\prime})]^{2}$, where $\mu_{\rm bd}$ is the coupling matrix element between dark and bright states. As a result, the expected photoluminescence spectrum contains two additional lines resulting from the recombination of these “semidark” trions and biexcitons, at an energy shifted downwards by $2\Delta^{\prime}$ relative to the bright complexes, and having a temperature-independent intensity.

II.2 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}^{\prime\prime}=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 9.

The total exciton wave function $X$ is given in general by the product of three components: the spatial envelope function $\Phi$, the Bloch or lattice wave functions of the electron and hole $U_{k}$, and the spin part $\chi$:

[TABLE]

The representations of the wave functions by irreps consist of the direct product of the individual irreps corresponding to the three components: $\Gamma_{X}=\Gamma_{\Phi}\otimes\Gamma_{U}\otimes\Gamma_{\chi}$. The tightly bound ground-state excitons are characterized by a maximally symmetrized envelope function corresponding to the identity irrep $\Gamma_{\Phi}=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}^{\prime}$ and $E_{2}^{\prime}$, respectively. Using the product table, Table 10, the lattice part of the exciton wave function transforms as

[TABLE]

where the 2D irrep $E^{+}$ corresponds to intravalley excitons in the $K$ and $K^{\prime}$ valleys, and $E^{\prime+}_{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-$\nicefrac{{1}}{{2}}$ particles corresponding to the spinor 2D irrep $D_{\nicefrac{{1}}{{2}}}$. The direct product of the two spinors can be decomposed into the crystal point group irreps as

[TABLE]

Hence the total exciton representation is given by

[TABLE]

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.Wang et al. (2017) 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 $\Delta^{\prime}$. The $A_{1}^{-}$ and $E^{-}$ excitons are not coupled to light. A summary of the classification of exciton states is given in Table 1 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 $\Delta^{\prime}$.

III Charge-carrier complexes in 2D

semiconductors

III.1 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.2.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 $\rho(x,y)\delta(z)$ in the $z=0$ plane of the 2D material, embedded in an isotropic medium of permittivity $\epsilon$. The resulting electric displacement field is ${\bf D}=-\epsilon\nabla\phi+{\bf P}_{\bot}(x,y)\delta(z)=-\epsilon\nabla\phi-\kappa[\nabla\phi(x,y,0)]\delta(z)$, where $\phi$ is the electrostatic potential, ${\bf P}_{\bot}(x,y)$ is the in-plane polarization, and $\kappa$ is the in-plane susceptibility of the material. By using Gauss’s law, $\nabla\cdot{\bf D}=\rho\delta(z)$, we obtain

[TABLE]

After taking the Fourier transform, denoting the wavevector in the $(x,y)$ plane by ${\bf q}$ and the wavenumber in the $z$ direction by $k$, we find

[TABLE]

However

[TABLE]

Rearranging, we find the in-plane electric potential to be

[TABLE]

Therefore the electrostatic potential energy between charges $q_{i}$ and $q_{j}$ in a 2D semiconductor is

[TABLE]

where $r_{\ast}\equiv\kappa/(2\epsilon)$. After taking the Fourier transform, the potential energy can be written as

[TABLE]

where $r$ is the separation of the particles and

[TABLE]

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,Keldysh (1979) and we refer to the interaction of Eq. (12) as the Keldysh interaction. At long range ($r\gg r_{\ast}$) this potential becomes a Coulomb interaction:

[TABLE]

while at short range ($r\ll r_{\ast}$) it is approximately logarithmic:

[TABLE]

where $\gamma$ 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_{\ast}=0$) and logarithmic ($r_{\ast}\to\infty$) approximations.

The following approximation to Eq. (12) was introduced in Ref. Cudazzo et al., 2011:

[TABLE]

This form of potential was used in the diffusion quantum Monte Carlo (DMC) study of Ref. Mayers et al., 2015. 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\approx r_{\ast}$. We compare DMC results obtained using Eqs. (12) and (15) in Sec. IV.6.

Finally, the Mott-Wannier-Keldysh Schrödinger equation for a set of charged quasiparticles in a 2D semiconductor is

[TABLE]

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 $\epsilon_{\rm a}$ above it and a dielectric medium of permittivity $\epsilon_{\rm 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 $\epsilon\equiv(\epsilon_{\rm a}+\epsilon_{\rm b})/2$ in the expressions above, the correction to the electrostatic energy of Eq. (11) is second order in $\epsilon_{\rm a}-\epsilon_{\rm b}$. Hence the Keldysh interaction remains valid when the permittivity $\epsilon$ is chosen to be the average of the permittivities of the media on either side of the 2D semiconductor, provided these permittivities are similar.

III.2 Units and scaling

III.2.1 Excitonic units

The energies of complexes interacting via the Keldysh or Coulomb interactions are given in terms of the exciton Rydberg, $R_{\rm y}^{\ast}=\mu e^{4}/[2(4\pi\epsilon)^{2}\hbar^{2}]$, and lengths are given in terms of the exciton Bohr radius, $a_{0}^{\ast}=4\pi\epsilon\hbar^{2}/(\mu e^{2})$, where $\mu=m_{\rm e}m_{\rm h}/(m_{\rm e}+m_{\rm h})$ is the reduced mass of electron-hole pairs, with $m_{\rm e}$ and $m_{\rm h}$ being the electron and hole masses, respectively.

Let $\tilde{\bf r}_{i}={\bf r}_{i}/a_{0}^{\ast}$. Then Eq. (16) can be written as

[TABLE]

where ${\cal E}=E/R_{\rm y}^{\ast}$. Note that $\mu/m_{i}$ only depends on the electron-hole mass ratio $\sigma\equiv m_{\rm e}/m_{\rm h}$. Hence for a fixed value of $r_{\ast}/a_{0}^{\ast}$, the dimensionless energy eigenvalues ${\cal 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 ${\bf r}_{\rm eh}$ as

[TABLE]

so that for a given value of $r_{\ast}/a_{0}^{\ast}$, the dimensionless exciton energy eigenvalues ${\cal E}_{\rm X}$ are also independent of the mass ratio. For the case of the Coulomb interaction ($r_{\ast}=0$), the dimensionless ground-state energy of an isolated exciton is ${\cal E}_{\rm 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_{\ast}/a_{0}^{\ast}$ and the electron-hole mass ratio $\sigma$. Unfortunately, the energies of the different complexes go to zero in these units in the limit that $r_{\ast}\to\infty$, and so a separate set of units is required for the case of the logarithmic interaction, as discussed in Sec. III.2.2.

III.2.2 Logarithmic interaction

For the limit $r_{\ast}\to\infty$, where the interaction is of logarithmic form, we use the dimensionless units introduced in Ref. Ganchev et al., 2015. The Schrödinger equation for a charge carrier complex with the logarithmic approximation to the interaction [Eq. (14)] is

[TABLE]

Let

[TABLE]

and

[TABLE]

Defining dimensionless coordinates $\tilde{\bf r}_{i}={\bf r}_{i}/r_{0}$ and a dimensionless energy ${\cal E}=E/E_{0}$, the Schrödinger equation can be written as

[TABLE]

The only dependence of the dimensionless energy ${\cal E}$ of the complex on $r_{\ast}$ is through the pairwise additive constant

[TABLE]

Note that

[TABLE]

Hence the additive constant ${\cal C}$ cancels out of the binding energies of the different charge-carrier complexes defined in Sec. III.3.

For an isolated exciton, we may write the Schrödinger equation in terms of the difference coordinate $r_{\rm eh}$ and reduced mass, giving

[TABLE]

The only dependence of the dimensionless energy eigenvalue ${\cal E}_{\rm X}$ on the mass ratio and $r_{\ast}$ comes from the constant term $\log(r_{0}/r_{\ast})$ in the Hamiltonian. Hence we may write the ground-state dimensionless energy as

[TABLE]

where ${\cal E}_{\rm X0}=0.41057747491(7)$ was evaluated by a finite-element method (see Sec. V.1).

III.3 Binding energies and spectra of charge-carrier

complexes

We define the binding energies $E_{\rm X^{-}}^{\rm b}$, $E_{\rm XX}^{\rm b}$, $E_{\rm D^{+}X}^{\rm b}$, $E_{\rm D^{0}X}^{\rm b}$, and $E_{\rm D^{+}XX}^{\rm b}$ of a trion, biexciton, donor-bound exciton, donor-bound trion, and donor-bound biexciton, respectively, as follows:

[TABLE]

where $E_{\rm X}$, $E_{\rm X^{-}}$, $E_{\rm XX}$, $E_{\rm D^{0}}$, $E_{\rm D^{+}X}$, $E_{\rm D^{0}X}$, and $E_{\rm D^{+}XX}$ are the ground-state total energies of an exciton, trion, biexciton, donor atom, donor-bound exciton, donor-bound trion, and donor-bound biexciton, respectively. These are the binding energies with respect to dissociation into the most energetically competitive species. With the exceptions of the donor-bound exciton (D+X) and donor-bound biexciton (D+XX), each of the complexes dissociates into an exciton (X) plus one other complex. Binding energies of charge-conjugated complexes (such as positive trions, acceptor-bound trions, and acceptor-bound biexcitons) are defined in an analogous fashion. Note that, under the definitions of Eqs. (27)–(31), a binding energy is positive for a bound complex.

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_{\rm X}-E_{\rm X^{-}}=E_{\rm X^{-}}^{\rm b}$; the energy difference between the exciton peak and the biexciton peak is $2E_{\rm X}-E_{\rm XX}=E_{\rm XX}^{\rm b}$; and the energy difference between the exciton peak and the donor-bound trion peak is $E_{\rm X}+E_{\rm D^{0}}-E_{\rm D^{0}X}=E_{\rm D^{0}X}^{\rm b}$. On the other hand, the energy difference between the exciton peak and the donor-bound exciton peak is $E_{\rm X}-E_{\rm D^{+}X}=E_{\rm D^{+}X}^{\rm b}+E_{\rm X}-E_{\rm D^{0}}$, and the energy difference between the exciton peak and the donor-bound biexciton peak is $E_{\rm X}+E_{\rm D^{+}X}-E_{\rm D^{+}XX}=E_{\rm D^{+}XX}^{\rm b}+E_{\rm D^{0}X}^{\rm b}-E_{\rm D^{+}X}^{\rm b}$. 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.6 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.

III.4 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.Szyniszewski et al. (2017) 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

[TABLE]

where $A^{\rm ee}$, $A^{\rm hh}$, and $A^{\rm eh}$ are constants and ${\bf r}_{\rm ee}$, ${\bf r}_{\rm hh}$, and ${\bf r}_{{\rm e}_{i}{\rm h}_{j}}$ are the electron-electron separation, the hole-hole separation, and the separation of electron $i$ and hole $j$, respectively. Evaluating $A^{\rm ee}$, $A^{\rm hh}$, and $A^{\rm 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^{\rm eh}g_{\rm XX}^{\rm eh}({\bf 0})+A^{\rm ee}g_{\rm XX}^{\rm ee}({\bf 0})+A^{\rm hh}g_{\rm XX}^{\rm hh}({\bf 0})$, where the electron-electron, hole-hole pair, and electron-hole PDFs are

[TABLE]

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 charge-carrier 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 electron-electron contact PDF for the semidark trion and biexciton in tungsten-based TMDCs determines the coupling strength of the dark and bright states as $\mu_{\rm bd}\propto g^{\rm ee}(0)$ and hence determines the recombination rates of the semidark states.Danovich et al. (2017)

IV Computational methodology

IV.1 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.Ceperley and Alder (1980); Foulkes et al. (2001) 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 minimizationUmrigar et al. (1988); Drummond and Needs (2005) and energy minimization.Umrigar et al. (2007) 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.Needs et al. (2010)

QMC methods have previously been used to study 2D trions with nonlocal screeningGanchev et al. (2015); Mayers et al. (2015) and the Coulomb interactionSpink et al. (2016) and 2D biexcitons with the Coulomb interaction (including indirect biexcitons in coupled-quantum-well heterostructures)Bressanini et al. (1998); Tan et al. (2005); Lee et al. (2009); Maezono et al. (2013) and in TMDCs with nonlocal screening.Mayers et al. (2015); Kylänpää and Komsa (2015) 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_{\ast}$ and effective masses, and have represented the DMC data using simple interpolation formulas.Szyniszewski et al. (2017) It was shown that for the applicable range of $r_{\ast}$ 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.

IV.2 Wave functions for complexes of distinguishable charge

carriers

Our trial wave functions for complexes of distinguishable charge carriers were of the Jastrow form $\Psi=\exp[J({\bf R})]$, where ${\bf R}$ is the vector of all the particle coordinates. The Jastrow exponent $J({\bf R})$ included a pairwise sum of terms of the form111For biexcitons with distinguishable particles, we used the Jastrow form introduced in the Supplemental Material to Ref. Ganchev et al., 2015.

[TABLE]

for the Keldysh and logarithmic interactions, where $r$ is interparticle distance, $c_{1}$, $c_{2}\leq 0$, and $c_{3}\geq 0$ are optimizable parameters, and

[TABLE]

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,Kato (1957); Pack and Byers Brown (1966) i.e., it ensures that the local energy $\Psi^{-1}\hat{H}\Psi$ is nondivergent at coalescence points, where $\hat{H}$ is the Hamiltonian operator.

Where the interaction between the charge carriers was of Coulomb form, we used pairwise terms of the form

[TABLE]

in the Jastrow exponent, where $c_{1}\leq 0$ and $c_{2}\geq 0$ are optimizable parameters, and

[TABLE]

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.Kato (1957); Pack and Byers Brown (1966)

Donor ions and other infinitely heavy particles were fixed point charges in our calculations. In this case $u_{\rm 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.Drummond et al. (2004); López Ríos et al. (2012)

IV.3 Wave functions for biexcitons with indistinguishable holes

For biexcitons with indistinguishable holes we used the trial wave function

[TABLE]

where $J$ is of the form described in Sec. IV.2. For indistinguishable particles of mass $m$ and charge $q$ interacting via the logarithmic or Keldysh interactions, Eq. (37) must be replaced by $\Gamma^{\prime}=-q^{2}m/(8a_{0}^{\ast}\mu e^{2}r_{\ast})$, while for indistinguishable pairs of particles interacting via the Coulomb interaction, Eq. (39) must be replaced by $\Gamma=q^{2}m/(2a_{0}^{\ast}\mu e^{2})$. $x_{\rm hh}^{\prime}$ is the $x$-component of

[TABLE]

where $\eta_{\rm hh}$ and $\eta_{\rm eh}$ are smoothly truncated polynomials, with optimizable expansion coefficients, and ${\bf r}_{\rm hh}$ and ${\bf r}_{{\rm e}_{i}{\rm h}_{j}}$ are the hole-hole and electron-hole relative positions, respectively. Equation (41) is effectively a backflowFeynman (1954); Feynman and Cohen (1956) transformation: $\Psi=\exp(J)x_{\rm hh}$ introduces the correct nodal topology for the state that we want to consider and Eq. (41) maps the particle coordinates $\{{\bf r}\}$ to quasiparticle coordinates $\{{\bf r}^{\prime}\}$ without changing the nodal topology. In Eq. (41),

[TABLE]

and

[TABLE]

are smoothly truncated polynomials with optimizable parameters $\{a_{n}\}$ and $\{b_{n}\}$. $L$ is a cutoff length, $N_{\eta}^{\rm hh}$ and $N_{\eta}^{\rm eh}$ determine the amount of variational freedom, $C=3$ to ensure smooth behavior at the cutoffs, and $\Theta$ denotes the Heaviside function. We require $b_{1}=Cb_{0}/L$ to ensure that $\eta$ 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.Umrigar et al. (2007)

For different values of $N_{\eta}^{\rm hh}$ and $N_{\eta}^{\rm 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 2. Analogous results for biexcitons interacting via the Keldysh interaction at finite $r_{\ast}$ are shown in Table 3. Our results show that increasing $N_{\eta}^{\rm hh}$ and $N_{\eta}^{\rm eh}$ slightly decreases the variances; nevertheless, the VMC and DMC energies are independent of the number of free parameters when $N_{\eta}^{\rm hh},N_{\eta}^{\rm eh}\geq 2$ to within our statistical error bars. We have used $N_{\eta}^{\rm hh}=N_{\eta}^{\rm 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.

IV.4 Time-step and population-control errors

We chose our DMC time steps such that the root-mean-square distance diffused by each particle in a single time step was much less than $r_{0}$ for the logarithmic interaction, much less than $a_{0}^{\ast}$ for the Coulomb interaction, and much less than $\min\{r_{0},a^{*}_{0}\}$ for the Keldysh interaction at finite $r_{\ast}$. 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.1 largely eliminates the effects of time-step bias, provided the time steps used are sufficiently small. For the logarithmic interaction with $\sigma=1$ and $r_{\ast}=r_{0}$, the time step should evidently be rather less than $0.04\hbar/E_{0}$.

Figures 6 and 7 show similar time-step tests performed for a negative trion and a donor-bound biexciton with the Keldysh interaction. For $r_{\ast}\geq 0.25a_{0}^{\ast}$, one should use time steps of less than $0.01\hbar/R_{\rm y}^{\ast}$ to be in the linear time-step bias regime, while for $r_{\ast}<0.25a_{0}^{\ast}$, time steps of less than $0.0025\hbar/R_{\rm y}^{\ast}$ are required.

IV.5 PDF calculations

The PDFs defined in Sec. III.4 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.Ceperley and Kalos (1979) 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_{\ast}$ value and mass ratio considered. To perform the extrapolation we fitted $\exp[\tilde{g}(r)]$ to our PDF data at short range,Drummond et al. (2011) where

[TABLE]

for the Keldysh and logarithmic interactions and

[TABLE]

for the Coulomb interaction ($r_{\ast}=0$), where $\Gamma^{\prime}$ and $\Gamma$ are defined in Eqs. (37) and (39) and $a_{0}^{\prime}$, $a_{2}^{\prime}$, …, $a_{6}^{\prime}$ and $a_{0}$, $a_{2}$, …, $a_{6}$ are fitting parameters. These forms satisfy (the analog of) the Kimball cusp conditions.Kimball (1973) The model functions were fitted to our PDF data at small $r$, with the data being weighted by $2\pi r$.

IV.6 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,Mayers et al. (2015) leads to significant errors. For an exciton with $r_{\ast}=a_{0}^{\ast}/2$, the DMC total energies are $E_{\rm X}=-1.5358899(2)R_{\rm y}^{\ast}$ and $-1.4668074(3)R_{\rm y}^{\ast}$ 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_{\ast}=a_{0}^{\ast}/2$ and mass ratio $\sigma=1$ using the exact and approximate Keldysh interactions are $0.1377(4)R_{\rm y}^{\ast}$ and $0.1335(3)R_{\rm y}^{\ast}$, 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.

V Numerical results

V.1 Excitons

The exciton ground-state energy is presented in Fig. 8. Our DMC data are in agreement with the results of finite-element calculations as implemented in the Mathematica software.Mat (2016) In excitonic units, the energy of an exciton is independent of the effective masses: see Sec. III.2. In the Coulomb limit, one recovers the well-known excitonic energy of $-4R_{\rm y}^{\ast}$. We can determine the behavior of the energy near the Coulomb limit by evaluating the first-order perturbative correction

[TABLE]

where $\Delta v=v_{\text{Keldysh}}-v_{\text{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 $\Psi=\exp(-2r/a_{0}^{\ast})$. The correction is shown in Fig. 8 as a green line.

We have numerically evaluated the dimensionless constant $\mathcal{E}_{\text{X0}}$ in Eq. (26) to be $\mathcal{E}_{\text{X0}}=0.41057739(7)$ using DMC and $\mathcal{E}_{\text{X0}}=0.41057747491(7)$ using the finite-element 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_{\ast}\to\infty$. The difference $\Delta E_{\rm X}/E_{0}$ between the exciton energies in units of $E_{0}$ with the Keldysh and logarithmic interactions at large $r_{\ast}$ was calculated numerically. Using the optimized ground-state wave function for the logarithmic interaction, we used VMC to evaluate the first-order perturbative approximation $\Delta E_{\rm X}/E_{0}\approx\langle v_{\text{Keldysh}}-v_{\text{logarithmic}}\rangle$. 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 $\sqrt{a_{0}^{\ast}/r_{\ast}}$.

We fitted the function

[TABLE]

with $a_{5}=-29+2\mathcal{E}_{\text{X}0}-a_{1}-a_{2}-a_{3}-a_{4}-\log 2$ to our DMC exciton energy data, where $y=r_{\ast}/(a_{0}^{\ast}+r_{\ast})$ 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.5. 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.sup (2016) In addition, the program reports fits to contact PDFs for the different clusters.

V.2 Biexcitons

V.2.1 Binding energies

We compare the stability of biexcitons with distinguishable and indistinguishable holes in the limit of the Coulomb interaction ($r_{\ast}=0$) in Fig. 12(a) and at $r_{\ast}=8a^{*}_{0}$ in Fig. 12(b). We find that biexcitons with indistinguishable holes are unbound for $\sigma\gtrsim 0.3$, while biexcitons consisting of distinguishable particles are bound at all mass ratios. The binding energies at $\sigma=0$ are obtained using the Born-Oppenheimer potentials as a function of heavy-hole separation $r$ plotted in Fig. 13. We fitted $U(r)=\alpha+\beta\sqrt{r}+\gamma r+\delta r^{2}$, where $\alpha$, $\beta$, $\gamma$, and $\delta$ 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)=\zeta+\eta\exp(-r/d)+\kappa\log(r)$ to our data, where $\zeta$, $\eta$, $d$, and $\kappa$ are fitting parameters. The Born-Oppenheimer approximation in Fig. 12(b) for heavy holes is in agreement with our DMC calculations at small $\sigma$. Analogous results obtained with the logarithmic interaction are shown in Fig. 12(c). For $\sigma\gtrsim 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=\sigma/(1+\sigma)$ and rescaled in-plane susceptibility $y=r_{\ast}/(a^{*}_{0}+r_{\ast})$. Our results are in agreement with path-integral Monte Carlo (PIMC) data at finite $r_{\ast}$, as shown in Fig. 15.Velizhanin and Saxena (2015) However, the PIMC data obtained by Velizhanin and Saxena have much larger statistical errors and they quoted a previous DMC resultStébé and Ainane (1989) at $r_{\ast}=0$ due to the infeasibility of PIMC in this case. The function

[TABLE]

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 $\sigma\to 0$, as derived in App. B.1, and is also invariant under charge conjugation ($m_{\rm e}\leftrightarrow m_{\rm 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. Szyniszewski et al., 2017. 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.sup (2016) Binding-energy results in the limit of large $r_{\ast}$, where the interaction is of logarithmic form, are given in Sec. V.7.

In Table 4, 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,Mayers et al. (2015) and also with PIMC calculations.Kylänpää and Komsa (2015) The small differences between DMC results in the literature must be due to the use of different effective masses, etc. Unfortunately, the theoretical biexciton binding energies are up to three times smaller than those reported in experimental works.Sie et al. (2015); Mai et al. (2014b); Plechinger et al. (2015); You et al. (2015) 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,Mayers et al. (2015); Szyniszewski et al. (2017) 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 semiconductors could only arise for one or more of the following four reasons: (i) the 2D Mott-Wannier model with the Keldysh interaction between charge carriers is incorrect or incomplete; (ii) the parameters (band effective masses and $r_{\ast}$ values) used in the model are incorrect; (iii) the exciton that remains after exciton recombination in a biexciton is not in its ground state;Zhang et al. (2015a) or (iv) the experimental spectra have been misinterpreted or the peaks have been misclassified.

As explained in Sec. III.4, there should be an additional contact interaction between charge carriers; however, the Mott-Wannier model with the Keldysh interaction apparently provides a good descriptionMayers et al. (2015); Szyniszewski et al. (2017) 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.2.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.You et al. (2015); Plechinger et al. (2015) 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.6) 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.

V.2.2 Sensitivity of binding energies to effective masses and in-plane

dielectric susceptibility

In Table 5 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 6. The energies depend relatively weakly on the in-plane permittivity $r_{\ast}$; the errors arising from the uncertainty in the effective mass almost certainly dominate errors arising from the uncertainty in $r_{\ast}$. 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 5 and/or the $r_{\ast}$ value would have to be an order of magnitude smaller. While there is still appreciable uncertainty in the ab initio effective mass and $r_{\ast}$ 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.

V.2.3 PDFs

In Fig. 16, we show the PDFs of biexcitons with distinguishable particles interacting via the logarithmic interaction for two different mass ratios, $\sigma=0.4$ and $\sigma=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 electron-electron contact $(r=0)$ PDFs for a biexciton. Notice that $g_{\rm XX}^{\rm eh}\approx 2g^{\rm eh}_{\rm X}$. Fits to the contact PDFs can be evaluated using the program supplied as Supplemental Material.sup (2016)

V.3 Trions

The binding energies of negative trions are presented in Fig. 18. We have fitted the function

[TABLE]

where $x=\sigma/(1+\sigma)$, $y=r_{\ast}/\left(r_{\ast}+a_{0}^{\ast}\right)$, 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 4. 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 6, 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_{\ast}$, where the interaction is of logarithmic form, are given in Sec. V.7.

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.sup (2016)

V.4 Donor/acceptor-bound excitons

We present the binding energies of donor-bound excitons in Fig. 19. For $\sigma\gtrsim 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 $\sigma\to\infty$, 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_{\ast}$ 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:

[TABLE]

In this expression $x=\sigma/(1+\sigma)$ and $y=r_{\ast}/(a_{0}^{\ast}+r_{\ast})$, while the $\{a_{ij}\}$ are fitting parameters. Our fitting function can be evaluated using the program in the Supplemental Material.sup (2016) We summarize our theoretical predictions for the binding energies of donor/acceptor-bound excitons in various TMDCs in Table 7. Binding-energies in the limit of large $r_{\ast}$, where the interaction is of logarithmic form, are given in Sec. V.7.

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 $\sigma\to\infty$, as expected, because the light hole becomes unbound in this limit. Contact PDFs can be evaluated using the program supplied as Supplemental Material.sup (2016)

V.5 Donor/acceptor-bound trions

Figure 21 presents the binding energies of donor-bound trions. We have devised the following 30-parameter fitting formula:

[TABLE]

which includes the correct divergence as $\sigma\to\infty$ and appropriate square-root behavior for the heavy-hole limit $\sigma\to 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).sup (2016) Binding-energy results in the limit of large $r_{\ast}$, where the interaction is of logarithmic form, are given in Sec. V.7.

Table 7 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.sup (2016)

V.6 Donor/acceptor-bound biexcitons

Donor-bound biexciton binding energies have also been calculated and are presented in Fig. 23. A 38-parameter fitting formula similar to that of a donor-bound trion [Eq. (51)] was used:

[TABLE]

where $x=\sigma/(1+\sigma)$ and $y=r_{\ast}/(a_{0}^{\ast}+r_{\ast})$, 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.sup (2016) We summarize our theoretical predictions for the binding energies of donor/acceptor-bound biexcitons in various TMDCs in Table 7. Binding-energy results in the limit of large $r_{\ast}$, where the interaction is of logarithmic form, are given in Sec. V.7.

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 ($\sigma\to 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_{\ast}/a_{0}^{\ast}=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 potential-energy surface.

Closely related to donor-bound biexcitons are five-carrier 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.222The most energetically competitive products of the dissociation of a quinton are either a trion plus an exciton, or a biexciton plus a free charge carrier; the energy difference between the two possible sets of products is given by the difference of the trion and biexciton binding energies. For all the 2D semiconductors considered here, the trion binding energy is larger than the biexciton binding energy,Szyniszewski et al. (2017) so that an exciton and a trion are the most energetically competitive species into which a quinton may dissociate. Other possible large complexes are donor-bound double-negative excitons (D-X, i.e., D+e-e-e-h+), donor-bound quintons (D0XX, 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 8 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 8. As with donor-bound biexcitons, the energies required to remove excitons from the larger complexes such as quintons are relatively close to the experimentally observed biexciton binding energies.

Figure 26 presents the contact PDFs of a donor-bound biexciton, which are also reported in the Supplemental Material.sup (2016)

V.7 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 excitonic units to logarithmic units we multiply the fitting function by $R_{\rm y}^{\ast}/E_{0}=r_{\ast}/(2a_{0}^{\ast})=y/(2-2y)$ and take the limit that $r_{\ast}\to\infty$, i.e., that $y\to 1$. For complexes that have been studied previously, our results are in good agreement with earlier exact calculations.Ganchev et al. (2015)

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 responsible for the observed peak ascribed to biexcitons.

Appendix A Character and product tables for the $D_{3h}^{\prime\prime}$ symmetry group

Character and product tables for the $D_{3h}^{\prime\prime}$ symmetry group of TMDCs are given in Tables 9 and 10.

Appendix B Molecular analogies and behavior of binding energies at extreme

mass ratios

B.1 Biexcitons

In the limit that the hole mass is large, a biexciton resembles a 2D H2 molecule, and we may use the Born-Oppenheimer approximation.Spink et al. (2016) The leading-order biexciton total energy is given by the minimum of the Born-Oppenheimer potential energy $U(r)$, where $r$ is the exciton-exciton separation, plus the harmonic zero-point energy of the exciton-exciton vibrations. Higher-order corrections to the energy arise from vibrational anharmonicity.

Consider a biexciton in which the charge carriers interact via the Keldysh interaction. Let ${\cal U}(r/a_{0})$ be the Born-Oppenheimer potential energy in Rydberg units for the case that the electron mass $m_{\rm e}$ is finite but the hole mass $m_{\rm h}$ is infinite. Then, at finite electron and hole masses, the Born-Oppenheimer potential is $U(r)={\cal U}(r/a_{0})R_{\rm y}$ where $R_{\rm y}=m_{\rm e}e^{4}/[2(4\pi\epsilon)^{2}\hbar^{2}]=(m_{\rm e}/\mu)R_{\rm y}^{\ast}$ is the Rydberg and $a_{0}=4\pi\epsilon\hbar^{2}/(m_{\rm e}e^{2})=(\mu/m_{\rm e})a_{0}^{\ast}$ is the Bohr radius. Note that ${\cal U}(r/a_{0})$ does not depend on the electron or hole mass.

Near the minimum of the potential $r_{\rm min}$ we may write

[TABLE]

where $(m_{\rm e}+m_{\rm h})/2$ is the reduced mass of the two excitons, $U_{\rm min}={\cal U}_{\rm min}R_{\rm y}$ is the minimum of potential, $U^{\prime\prime}_{\rm min}={\cal U}^{\prime\prime}_{\rm min}R_{\rm y}/a_{0}^{2}$ is the second derivative of the potential at the minimum, and

[TABLE]

The resulting ground-state energy in the harmonic approximation is

[TABLE]

where we have used $m_{\rm h}\gg m_{\rm e}$ in the last step. This suggests that a suitable fitting function for the binding energy of a biexciton with small $\sigma\equiv m_{\rm e}/m_{\rm h}$ is a polynomial in powers of $\sqrt{\sigma}$. 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.

B.2 Negative trions

In the limit of heavy holes ($\sigma\to 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 $\sigma$.

In the limit of heavy electrons ($\sigma\to\infty$), 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 square-root of the mass ratio $\sigma$.

B.3 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 $\sigma$.

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.

B.4 Donor-bound trions

The Born-Oppenheimer potential energy curve of a donor-bound trion with a heavy hole is the same as that of a biexciton, but this time the reduced mass is simply equal to the exciton mass. The binding energy varies as the square root of the mass ratio $\sigma$.

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 (${\rm 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_{\rm D^{0}X}\approx E_{\rm D^{-}}+E_{\rm A^{0}}$ in this limit, where $E_{\rm 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_{\rm X}\approx E_{\rm 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 $E_{\rm D^{0}X}^{\rm b}=E_{\rm D^{0}}+E_{\rm X}-E_{\rm D^{0}X}\approx E_{\rm D^{0}}-E_{\rm D^{-}}$, 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 ($\sigma\to\infty$).

B.5 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,Martin et al. (1961) and it is known that the protons in H${}_{3}^{+}$ form an equilateral triangle. In Sec. V.6 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 electrons and two light holes. The binding energy of a donor-bound biexciton in the limit that the holes (h${}^{+}_{\rm light}$) are much lighter than the electrons (e${}^{-}_{\rm heavy}$) is

[TABLE]

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.

Acknowledgements.

We acknowledge support from the European Research Council synergy grant Hetero2D, the United Kingdom’s Engineering and Physical Sciences Research Council (EPSRC), and the European Graphene Flagship project. M.S. was funded by the EPSRC doctoral training centre “NoWNANO” (grant no. EP/G03737X/1) and R.J.H. was funded by the EPSRC doctoral training centre “Graphene NOWNANO” (grant no. EP/L01548X/1). Computer resources were provided by Lancaster University’s High-End Computing Cluster. This work made use of the facilities of N8 HPC provided and funded by the N8 consortium and EPSRC (grant no. EP/K000225/1). R.M. is grateful for financial support from MEXT-KAKENHI grant 17H05478, support by FLAGSHIP2020, MEXT for computational resources, and project nos. hp170269 and hp170220 at the K-computer. We acknowledge useful discussions with I. Aleiner, A. Tartakovski, M. Potemski, and T. Heinz.