Quadriexciton binding energy in electron-hole bilayers

Excitonic condensation and superfluidity have recently received a renewed attention, due to the fabrication of bilayer systems in which electrons and hole are spatially separated and form stable pairs known as indirect excitons. Dichalcogenides- and graphene- based bilayers are nowadays built and investigated, giving access to systems with (i) only spin degeneracy, (ii) spin and valley degeneracy. Simulation studies performed in the last decade at $T=0$ for simple, model electron-hole bilayers, as function of inter-layer distance and in-layer carrier density, have revealed in case (i) the formation of biexcitons in a tiny region of parameter space and in case (ii) the formation of stable compounds made of 4 electrons and 4 holes (quadriexcitons) in a sizable region of parameter space. Of some interest is the relation of the properties of isolated biexcitons (quadriexcitons) and those of their finite density counterpart. In fact, the isolated biexciton has been repeatedly studied in the last years with simulations and other techniques. No simulations, instead, are available to our knowledge for the isolated quadriexciton, for which we present here results of the first quantum Monte Carlo (QMC) study. Stability with respect to the dissociation into biexcitons, and the pair correlations with varying the inter-layer distance $d$ are discussed.


Introduction
An exciton, a bound electron-hole pair in a semiconductor [1], is an elementary optical excitation or quasiparticle in solid state physics, reminiscent of positronium or hydrogen in atomic physics. Here we are interested in Mott-Wannier excitons. Interactions between two excitons can result in the formation of an excitonic molecule [2,3], a new quasi-particle known as biexciton [3]. The study of excitons and biexcitons is most naturally performed within the envelope function and effective mass approximations [4], which provide an effective Hamiltonian for interacting electron and hole quasi-particles, near energy band extrema. Whenever a band has multiple equivalent minima (maxima) the description of electrons (holes) implies a new discrete quantum number, the valley index v, which identifies a minimum (maximum). The valley index can be formally treated as a pseudospin; the number g v of equivalent valleys fixes the length τ of the pseudospin, 2τ + 1 = g v → τ = (g v − 1)/2. We can introduce a flavor (or valley-spin) index σ = (τ z , s z ) to fix the spin and pseudospin states. As electrons and holes have spin 1/2, the total number of flavors is N c = 2g v for each particle type, electron or hole. We assume that the Hamiltonian does not contain spin or pseudospin operators. It is evident that if one considers a finite electron-hole system with N e ≤ N c , N h ≤ N c , with N e (N h ) the number of electrons (holes), the ground state wave function Φ can be exactly factorized in Φ = Ψζ, with Ψ a nodeless wave function of the particles' Cartesian coordinates and ζ a wave function of pseudospin and spin coordinates; Ψ will be symmetric for any pair exchange of electrons (holes) coordinates and ζ will be antisymmetric for any pair exchange of electrons (holes) coordinates. One can of course choose Ψ to also have other symmetries possessed by the Hamiltonian, provided that this does not spoil the exchange symmetry and Ψ remains nodeless.
Long ago it was observed [5] that when there are g v equivalent minima in the electron conduction band, up to N c = 2g v electrons can occupy the same molecular orbital and therefore electron-hole complexes, where up to N c excitons are possible; the same is true [5] when N c equivalent maxima are present in the valence band. Clearly, in a system with g v = 2, such as coupled graphene bilayers [6][7][8][9][10][11], N c = 4 and quadriexcitons should be possible. We recall that the problem of an isolated biexciton has already been studied [12][13][14] with DMC and other techniques.
Before tackling below the treatment of an isolated quadriexciton (or biexciton ) at T = 0, in the simplest, model electron-hole bilayer (the paramagnetic, symmetric bilayer: m * e = m * h = m b ), we briefly summarize here the results of available QMC simulations of systems with finite in-layer carrier density n [15][16][17], which we specify via the dimensionless r s parameter, defined by πr 2 s a * B 2 = 1/n where a * B = 4πϵ 0 ϵh 2 /(m b e 2 ). The first simulations [15] were performed in the region of parameter space 0 < d/(r s a * B ) ≤ 3, 0 ≤ r s ≤ 30) with g v = 1. It was found that at large density r s ≲ 1, the plasma phase was stable at all distances d, due to the screening of the Coulomb attraction by carriers. When increasing r s , an excitonic condensate appeared for d < d X (r s ) with d X (r s ) an increasing function of r s . For d > d X (r s ), the plasma phase was stable up to an intermediate density, r s ≲ 20, while turning into a Wigner crystal for 20 ≲ r s ≲ 30. More recent simulations [16], performed in a much smaller region of the parameter space , using a much more flexible wave function and up-to-date computer resources, confirmed semiquantitatively the phase boundary between the excitonic and plasma phases as well as the condensate fractions found in [15], while finding for 4 ≲ r s ≲ 8 a new, stable biexcitonic phase for d ≲ d 2X (r s ) = 0.05a * B . Finally, the results of QMC simulations for and electron-hole bilayer with g v = 2, 0 < d/a * B ≤ 3.5 and 0 ≤ r s ≤ 8 have just been published [17]. As predicted by [5], a quadriexcitonic phase was found and in a sizeable region of the parameter space, compared with the region of stability of the biexcitonic phase found for g v = 1. The new phase appeared for r s ≳ 1 and 0 ≲ d/ ≲ d 4X (r s ), with d 4X (r s ) an increasing function of r s . At a given r s ≳ 1 and for d 4X (r s ) ≲ d ≲ d X (r s ), there was the excitonic phase, followed by the plasma phase for d ≳ d X (r s ). No trace was found of the biexcitonic phase.

Hamiltonians and Wave Functions
Within the envelope-function effective mass approximation [4], the Hamiltonian of a two-dimensional bilayer reads (1) This Hamiltonian, apart from the assumption of isotropic masses, is quite general: v and v ′ run over the valleys and i and i ′ run over the number of electrons or holes in a given valley. A prime in the sum excludes terms with both v = v ′ and i = i ′ . The paramagnetic, symmetric bilayer is the one in which m * h = m * e = m b , all valleys present have the same electron (hole) population and in each valley, there is the same number of spin-up and spin-down electrons (holes). In the rest of this paper, we measure distances using the effective Bohr radius a * B = 4πϵ 0 ϵh 2 /(m b e 2 ) and energies in effective Rydbergs Ry * =h 2 /(2m b a * B 2 ), so that the Hamiltonian becomes (2) As we argued above for N e = N h = N c ≡ N, the ground-state wave function Ψ 0 NX (R e , R h ), with R e (R h ) the electrons (holes) coordinates, (i) must be the solution of (ii) it must be symmetric for any pair exchange of electrons (holes) coordinates; (iii) in the symmetric bilayer, it must be symmetric for the electron-hole exchange, i.e., Ψ 0 [12].
For the biexciton and the quadriexciton, we tackle the problem in Equation (3) resorting to variational and diffusion Monte Carlo [18][19][20] (VMC and DMC) as implemented in our own code. For a given N, at each d, an optimal trial function Ψ T NX (R e , R h ; c) is determined by minimizing the variational energy with respect to a number of optimizable parameters c [21][22][23]. We then compute the VMC estimates of the properties of interest using a Monte Carlo integration with |Ψ T NX | 2 as the importance function; the DMC estimates are obtained by employing the optimized Ψ T NX as the guiding function. We now turn to the explicit form of Ψ T NX .

Polyexciton Wave Function Ψ T NX
Here, we restrict to polyexcitons made by N c electrons and N c holes. Given the symmetry requests on the wave function it is natural to write Ψ T NX (R e , R h ) as a symmetrized Jastrow factor. We start from the unsymmetrized form embodying two-body pseudopotentials among all particles. Above the species index, µ = (t, σ) combines the particle type (t = e, h) and the flavor σ = (τ z , s z ); moreover, the primed sum for µ ′ = µ contains only the terms with i µ ̸ = j µ . Evidently Let us inspect the three terms above. As there is only one electron per flavor, only interflavor terms (σ ̸ = σ ′ ) survive in and in fact, the sum over i e,σ and j e,σ ′ in the equation above has just one term; thus, where we chose the interflavor pseudopotentials as being all equal, u e,σ;e,σ ′ (r) = ϕ(r). Since we are considering a paramagnetic, symmetric bilayer, it immediately follows that We note that with the choices made, the product J ee (R e )J hh (R h ) is symmetric under any electron-hole pair exchange as well as under any electrons-holes exchange. Let us turn now to J eh (R e , R h ). We have Above, we chose one pseudopotential for electron-hole pairs with the same flavor, u e,σ;h,σ (r) = ψ(r), and another for electron-hole pairs with different flavor, u e,σ;h,σ ′ (r) = ψ(r). The possibility ofψ(r) being different from ψ(r) should allow the fragmentation of the polyexcitons in smaller components, namely, excitons [12,14]. This choice, however, breaks the symmetry under an electron-electron pair exchange, as is clear by inspecting Equation (9). Thus, we restore the symmetry by taking the symmetrized version of whereP e permutes the electron coordinates. Finally, the trial wave function is nodeless and satisfies all the required symmetry properties: the symmetry under the pair exchange of electrons (holes) and the symmetry under the exchange of electrons and holes. We note that in the unsymmetrized electron-hole Jastrow of Equation (9), the pseudopotential ψ(r) describes intraexciton correlations, whileψ(r) describes interexciton correlations.
We take the pseudopotentials of the Padé form. For the electron-electron (hole-hole) pair: For the electron-hole pair: Some of the parameters c 1 -c 7 are fixed by exact conditions; the others are determined through an energy minimization as follows: The parameter c 1 is fixed by the electronelectron (hole-hole) Kato cusp conditions [24]. The parameter c 3 is fixed by the electron-hole cusp condition at d = 0, while for d > 0, there is no electron-hole cusp and so we set c 3 = 0 [12]. We require c 2 , c 5 , c 7 > 0 to avoid divergences and c 4 , c 6 < 0 to describe the wave function decay when the electron and the hole are far apart. The above wave function should describe separated excitons when either c 4 or c 6 goes to zero.

Exciton Wave Function
The Schrödinger equation for the ground state of the isolated exciton, in the center-ofmass reference frame, reads with r = |r e − r h | the in-plane distance between the electron and the hole. At d = 0, Equation (15) has the simple hydrogenic normalized solution Φ 0 However, at d > 0, there is no closed form for the solution, which must be computed numerically. To this end, we used the simple and accurate Numerov algorithm [25,26], which is especially suited for atomic-like problems.

Results
In Table 1, we report the ground state energy per particle of the exciton, biexciton and quadriexciton as a function of the interlayer distance d. The energies of the biexciton and quadriexciton were obtained by DMC simulations using as guiding function the optimized trial function of Equation (11). The optimization with respect to the free parameters was performed by minimizing the variational energy using the linear method [21,22] and checking, in selected cases, that the obtained minimum agreed with the one found by the improved stochastic reconfiguration method [23]. DMC simulations are affected by the walker population's bias (finite number of walkers N w ) and finite time-step bias. However, for small systems, such as those studied here, these biases can easily be made negligible. To this end we used a large number of walkers, N w = 1760, and performed an extrapolation to the zero time step. A good optimization and bias reduction are especially important when estimating pair correlation functions. Table 1. Energy per particle E X /2 (exciton), E 2X /4 (biexciton) and E 4X /8 (quadriexciton) for various distances d. Excitonic energies were estimated using the Numerov algorithm [25,26]. The energies for the excitonic complexes were obtained from DMC simulations with N w = 1760 walkers and the time-step bias removed by an extrapolation to the zero time step.  Table 1, it is evident that for d = 0, we have E X /2 > E 2X /4 > E 4X /8 and that with an increasing d, the energies of all three "phases" increase, become closer and could presumably cross or merge. So, studying, for instance, the stability of the biexciton with respect to the separation into two excitons requires the extrapolation of the available biexciton energies at larger distances d. In fact, a more efficient manner to study the stability of polyexcitons is to define a binding energy with respect to the separation into fragments. The binding energy of a polyexciton is defined as the energy released upon the formation of a bound state with respect to the initially isolated constituents (fragments). Therefore, for a biexciton, we can define and on the same footing, the quadriexciton binding energy with respect to two isolated biexcitons is defined as: In all the figures with DMC energies, the statistical error are reported, though they are not always visible.

Biexciton Binding Energy
Meyertholen and Fogler [13] showed that when increasing the interlayer distance, the binding energy of a biexciton with respect to two excitons vanished at a finite critical distance d c and for d ≲ d c , it had the behavior where D is a positive constant of order one, γ = 0.577... is the Euler-Mascheroni constant and D/(d c − d) ≫ 1. Clearly the last condition implies E B ≪ E 0 . Their derivation is built on the fact that at at large d, the biexciton is weakly bound and the exciton-exciton interaction brings into play at a large interexciton distance the dipole-dipole interaction, present for d > 0. Clearly the larger the d, the larger the exciton dipole and the greater the importance of such an interaction. In brief, they neglected what happens at a small interexcitonic distance and wrote a Schrödinger equation for two interacting excitons valid at an intermediate distance (first region) and at a large distance (second region): in the first region the binding energy was neglected with respect to the dipole-dipole interaction and in the second, the dipole-dipole interaction was neglected with respect to the binding energy. Imposing the continuity of the logarithmic derivative of the wave function when going from the first to the second region, they obtained an expression involving the binding energy of the biexciton and a characteristic energy of the problem, E 0 (d). Further manipulations finally yielded Equation (18). On the basis of their analysis, they proposed a fit of the binding energy in the form with d c , D, D 1 fitting parameters.
In Figure 1, we report our DMC results for the binding energy E B (2X, d) of the biexciton, together with the results of the stochastic variational method (SVM)of [13]. In Figure 2, we report the same data of Figure 1 in a form suitable for the fit according to Equation (20). One effect of using the "lens" provided by the logarithm is to emphasize the region of distances d near d c as well as the need of extrapolating to d c . It is also clear from Figure 2 that the SVM and DMC results are in good agreement, with the DMC results covering a larger range, which should be important for the fit. Minor differences are found in the values of d c provided by the fits using the two data sets. We should also stress that in the range of distances covered in Figures 1 and 2, E 0 (d) is at least two orders of magnitude larger than E B , so that the requirement E B ≪ E 0 is fully fulfilled.     [13]. We compare results from our DMC simulation and results from Ref. [13].

Quadriexciton Binding Energy
In principle, one may study the quadriexciton stability following the treatment for the biexciton [13], summarized in some detail in the previous section, also assuming that for the quadriexciton, the binding energy vanishes at a finite critical distance d c . Therefore, at distances d ≲ d c , the quadriexciton is weakly bound and the biexciton-biexciton interaction (dipole-dipole interaction), comes into play. The Schrödinger equation for the relative motion of the two interacting biexcitons is similar to the one for the two interacting excitons, with the following differences. With respect to the exciton-exciton case, the reduced mass is twice larger and the dipole-dipole coupling is four times larger. This yields, when imposing the continuity of the logarithmic derivative of the wave function at large distances (see the biexciton case above), an expression involving the ratio of the binding energy of the quadriexciton E B (4X, d) and a characteristic energy of the problem which turns out to be with E 0 (d) given in Equation (19). However, the self-consistency of the procedure would require E B (4X, d) ≪ E 0 (4X, d), a condition which is not fulfilled in the present case, as is clear from Figure 3, where our DMC results for E B (4X, d) are displayed together with E 0 (4X, d). We therefore fitted functions different from the one used for the biexciton and given in Equation (20). The first choice was  The results of the fits to our DMC energies are displayed in Figure 4 in terms of the logarithm of E B (4X, d), to exploit the "lens" effect. The two different recipes (the second somewhat inspired by the biexciton fit function) gave compatible results as far as the d c value was concerned. However, the second choice reproduced the DMC energies in a larger interval of distances, i.e., including also energies that were not included in the fitting.

Quadriexciton Pair Correlation Functions
In a finite system, pair correlation functions can be defined as appropriate "center-ofmass" averages of the two-body density. Let us consider electron-electron correlations. As we have four electrons (one per flavor), we can start from the two-body electronelectron density ρ σσ ′ ee (r 1 , r 2 ) = ⟨δ(r 1 − r e,σ )δ(r 2 − r e,σ ′ )⟩, where ⟨...⟩ indicates a ground state average and only the different-flavor case is present. We define new vectors as R = (r 1 + r 2 )/2 and r = r 1 − r 2 , so that r 1 = R + r/2 and r 2 = R − r/2. Substituting in Equation (24), we get and integrating over the "center-of-mass" vector R, we obtain a pair correlation function g σσ ′ ee (r) = ⟨δ(r − r e,σ + r e,σ ′ )⟩.
A few observations are in order. First, g σσ ′ ee (r) is a dimensional quantity; it has the dimensions of a density. However, for our finite system, there is no sensible density to divide by to get a dimensionless pair correlation function, as in extended systems. The integration of g σσ ′ ee (r) with respect to r over the whole space is one. One usually finds convenient to take the angles' average of g σσ ′ ee (r) to get a g σσ ′ ee (r). In a similar manner, one may define electron-hole pair correlation functions as g σσ ′ eh (r) = ⟨δ(r − r e,σ + r h,σ ′ )⟩, (27) having in this case both same-flavor σ ′ = σ and different-flavor σ ′ ̸ = σ functions. We should mention that the above definitions for the pair correlation functions coincide with that of the intracule density in chemistry as defined in [27].
In Figure 5, we show the electron-hole pair correlation functions, at various distances d, averaged on the azimuthal angle. Clearly, the contact value g σσ ′ eh (r = 0) is a decreasing function of d. At the largest d shown (d = 0.65), there is still an appreciable pileup of probability at the origin. In Figure 6, we show the electron-electron pair correlation functions. A peak at distances on the scale of a few Bohr radii is present, due to the quadriexciton binding; similarly to what happens to the contact values of g σσ ′ eh (r), such a peak is a decreasing function of d. Our results were similar to those in [14] for the biexciton. One should keep in mind that different units were used in the two studies, so our pair correlation functions should be multiplied by four, before comparing with those of the biexciton case in [14]. Finally, as a symmetry check, we verified that the electronhole pair correlation functions with the same or different flavor coincided; similarly, the electron-electron pair correlation functions did not depend on the pair of different flavors chosen.

Discussion
We studied a system of four electrons and four holes (quadriexciton) in a symmetric, paramagnetic bilayer, where indirect excitons are formed, focusing on (i) the unbinding of the quadriexciton into two biexcitons and (ii) the pair correlation functions. To this end, we used state-of-the-art QMC simulations in which both quadriexcitons and biexcitons were studied. We already analyzed the various properties of biexcitons and quadriexcitons in the foregoing. Here, we simply discuss what are the implications of our results and analysis in relation to the phase diagram of the symmetric, paramagnetic bilayer at a finite density [16,17].
According to our analysis and the one in [13], an isolated biexciton unbinds into two excitons at d c = 0.84a * B . Evidently, one might speculate that such result could be of some relevance to the condensed phase, though only at a very small density. The phase diagram in [16] covers the density range 0 ≤ r s ≤ 8 and reveals a biexcitonic phase only for 4 ≲ r s ≲ 8 and d ≲ d 2X (r s ) = 0.05a * B . One may therefore conclude that at such densities the electronic screening is very effective in destabilizing the biexcitonic phase and favoring its melting into excitons, to the point that the biexcitonic phase almost disappears. The situation appears qualitatively different in the symmetric, paramagnetic bilayer with valley degeneracy [17]. In this system, a robust quadriexcitonic phase is present in the same density range 0 ≤ r s ≤ 8. It first appears around r s = 1.5 at d = 0 but with increasing r s , the quadriexciton-exciton boundary increases up to d ≃ 0.65a * B at r s = 8. This d value is just 12% smaller than the d = 0.74a * B at which the isolated quadriexciton melts into two biexcitons. One may then observe, comparing Figures 1 and 3, that the binding energy of the quadriexciton is about one order of magnitude larger than that of the biexciton for 0.5a * B ≤ d ≤ 0.65a * B . Is that the reason of the much greater stability of the quadriexciton in the condensed phase, as compared with the biexciton? We plan to further study the properties of polyexcitons in the near future to answer this and other questions.
Our goal in this paper was to provide benchmark simulation results for the simplest model electron-hole bilayer when valley degeneracy was added. Though we were not trying to accurately model real devices, we could nevertheless link our results to fabricated coupled graphene bilayers, by estimating the effective Bohr radius a * B in actual devices and hence d exp /a * B , where d exp is the actual inter-bilayer distance. This procedure relied on crude assumptions on the form of the interparticle interactions, as well as on the chosen values of the dielectric constants and band masses. We gave three examples, using the isotropic interactions in Equation (1). In [11], coupled graphene bilayers with d exp = 2.2 nm were studied. Assuming m b = 0.04m e and ϵ = 3 or ϵ = 14 (the smallest or the largest of the dielectric constants of the constituent materials), we obtained d exp /a * B = 0.55 or d exp /a * B = 0.12. In [7], coupled graphene bilayers with 5 nm≤ d exp ≤ 12 nm were investigated. Choosing m b = 0.04m e and ϵ = 3, the authors obtained 1.3 ≤ d exp /a * B ≤ 3.0. Finally, in [8], where devices with the same band mass and dielectric constant as in [7] were studied, with 2nm ≤ d exp ≤ 5nm, the authors obtained 0.50 ≤ d exp /a * B ≤ 1.3. Comparing the estimates of d exc /a * B for the three sets of experiments [7,8,11] with the results for our model electron-hole bilayer, one would conclude that only the devices studied in [8,11] would be compatible with the observation of an isolated quadriexciton.
Author Contributions: Conceptualization, G.S. and S.D.P.; methodology, G.S. and S.D.P.; investigation, S.D.P. and C.M. All the authors discussed the details of the work and contributed to the writing of the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.