Crystal Phases of Charged Interlayer Excitons in van der Waals Heterostructures

Throughout the years, strongly correlated coherent states of excitons have been the subject of intense theoretical and experimental studies. This topic has recently boomed due to new emerging quantum materials such as van der Waals (vdW) bound atomically thin layers of transition metal dichalcogenides (TMDs). We analyze the collective properties of charged interlayer excitons observed recently in bilayer TMD heterostructures. We predict new strongly correlated phases - crystal and Wigner crystal - that can be selectively realized with TMD bilayers of properly chosen electron-hole effective masses by just varying their interlayer separation distance. Our results open up new avenues for nonlinear coherent control, charge transport and spinoptronics applications with quantum vdW heterostuctures.

For bilayer TMD heterostructures, controlled optical and electrical generation of IEs and charged IEs (CIEs, also known as trions formed by indirect excitons [33]) has lately been achieved [23,24]. Their in-plane propagation through the sample was adjusted by the excitation power and perpendicular electrostatic field. These experiments exhibit a unique potential of TMD bilayers for achieving precise control over compound quantum particles of both bosonic (IE) and fermionic (CIE) nature. The CIEs offer even more flexibility in this respect as they have both net charge and permanent dipole moment as well as non-zero spin ( Fig. 1), to allow for electrical tunability and optical spin manipulation in charge transport and spinoptronics experiments with quasi-2D vdW heterostructures.
Here, we consider the collective properties of the negative and positive CIEs starting with their binding energies in bilayer quasi-2D semiconductor heterostructures. We derive the general analytical expressions as functions of the electron-hole effective mass ratio and interlayer separation distance to explain the experimental evidence earlier reported for the negative CIE to have a greater binding energy than that of the positive CIE [23].
Our analysis of the pairwise interactions between the CIEs, as sketched in Fig. 1, exhibits two scenarios for crystallization phase transitions in the collective multiparticle CIE system. They are the crystallization of the unlike-charge CIEs and the Wigner crystallization of the like-charge CIEs, which can be selectively realized in practice by choosing bilayers with appropriate electron-hole effective mass ratio and interlayer separation in addition to the standard technique of electrostatic doping. We conclude that this strongly correlated multiexciton phenomenon of CIE crystallization can be realized in layered van der Waals heterostructures such as double bilayer graphene and bilayer TMD systems [23,31], to open up new avenues for nonlinear coherent optical control and spinoptronics applications with charged interlayer excitons.

A. The binding energy
The compound structure of the CIE complexes of interest is sketched in Fig. 1. We use the configuration space approach [34] to derive the binding energy expressions for the CIEs as functions of their electron-hole effective mass ratio σ = m e /m h and interlayer separation distance d. This approach was recently proven to be efficient as applied to quasi-1D [35] and quasi-2D bilayer semiconductors [33] where it offers easily tractable analytical solutions to reveal universal relations between the binding energy of the complex of interest and that of the 1D-exciton or that of the indirect (interlayer) exciton [36], respectively. The method itself was originally pioneered by Landau [37], Gor'kov and Pitaevski [38], Holstein and Herring [39] in their studies of molecular binding and magnetism.
The negative X − (positive X + ) trion complex in Fig. 1 can be viewed as two equivalent IEs sharing the same hole (electron). The CIE bound state then forms due to the exchange under-barrier tunneling between the equivalent configurations of the electron-hole system in the configuration space of the two independent relative electron-hole motion coordinates representing the two equivalent IEs that are separated by the center-of-mass-to-center-ofmass distance ∆ρ. The binding strength is controlled by the exchange tunneling rate integral J X ± (∆ρ). The CIE binding energy is with ∆ρ X ± to be determined from an appropriate variational procedure to maximize the tunneling rate, which corresponds to the CIE ground state. This approach gives an upper bound for the (negative) ground state binding energy of an exciton complex of interest [33][34][35]. It captures essential kinematics of the formation of the complex and helps understand the general physical principles to underlie its stability.
Using the configuration space method for solving the CIE ground state binding energy problem, we obtain (see Methods) are the interlayer separation dependent constants coming from the indirect (interlayer) exciton wave function [36], and the upper or lower term should be taken in the curly brackets for the positive or negative CIE, respectively. Here the 3D "atomic units" are used [36][37][38][39], with distance and energy measured in the units of exciton Bohr radius a * B = 0.529Å ε/µ and exciton Rydberg energy Ry * =h 2 /(2µ m 0 a * 2 B ) = e 2 /(2εa * B ) = 13.6 eV µ/ε 2 , respectively, ε represents the effective average dielectric constant of the bilayer heterostructure and µ = m e /(λ m 0 ) stands for the exciton reduced effective mass (in the units of free electron mass m 0 ) with λ = 1+m e /m h = 1+σ.
The image-charge effects are neglected [36]. To properly take into account the screening effect for the charges forming the CIEs as sketched in Fig. 1, we used the Keldysh-Rytova (KR) interaction potential energy (see Refs. [40]) approximated by elementary functions in the form (atomic units) proposed for atomically thin layers in Ref. [41], to represent the effective electrostatic potential energy for like charges in monolayers. Here, ρ is the in-plane intercharge distance and r 0 = 2πχ 2D is the screening length parameter with χ 2D being the in-plane polarizability of 2D material [41,42]. For unlike charges the interlayer electrostatic potential energy is taken in the standard screened Coulomb form V C (r) = −1/r with r = ρ 2 + d 2 (atomic units).
The function J X ± (∆ρ) in Eq. (2) Substituting this in Eq. (1), one obtains the positive and negative CIE binding energies of interest.
Figure 2 (a) shows the binding energies E X + and E X − calculated from Eqs. (1), (2) and (4) with σ = 1 and 0.7 as functions of d and r 0 . For σ = 1 they coincide [33]. For σ = 0.7 the positive-negative CIE binding energy splitting is seen to occur in the entire domain of parameters used. Figure 2 (b) shows the crosscuts of Fig. 2 (a) for r 0 = 0 and 0.1 to exhibit the remarkable features of the screening and binding energy splitting effects. The screening of like charges in the CIE complex is seen to increase its binding energy. The X ± trion energy splitting at short d is such that |E X −| > |E X +|, which agrees with and thus explains the measurements reported recently for (h-BN)-encapsulated MoSe 2 -WSe 2 bilayer heterostructures [23]. As d increases the crossover occurs to give |E X −| < |E X +| with |E X − | quickly going down to zero, which is also seen in Fig. 2 (a). On closer inspection of Eqs. (2) and (4) it can be seen though that |E X −| and |E X +| swap places for σ > 1 (not shown here), thereby offering an extra functionality for properly fabricated vdW heterostructures [43,44]. Equation links the electron-hole mass ratio σ and interlayer separation d at which the crossover occurs. For σ = 1 it turns into an identity [33]. The three lines in Fig. 2 (c) present the nontrivial solution to this equation, σ(d), for three different r 0 values. The screening is seen to shrink the |E X −| > |E X +| domain and expand the |E X −| < |E X +| domain (above and below the solution line, respectively). Since the greater binding energy increases the formation probability, these domains are also those to preferentially form the X − and X + trion, respectively, while the constraint E X − = E X + defines the line of equal X ± formation probabilities. Thus by varying d for a properly chosen TMD bilayer composition with known σ, one can selectively control intrinsic positive/negative CIE formation in an undoped heterostructure as opposed to the electrostatic doping technique.

B. Unlike-charge trion crystallization
For undoped structures of two monolayers with σ = 1 as well as for those with σ = 1 fabricated to hit the E X − = E X + line, both X − and X + trions are equally likely to form under intense external irradiation at not too high temperatures T < |E X ±|/k B . This results in an overall neutral two-component many-particle mixture of X − and X + trions. The aggregate state of a many-particle system is defined by its Helmholtz free energy consisting of the total energy term and the entropy term. The entropy term becomes dominant at high T to favor configurations with greater randomness. At not too high T the total energy termthe sum of kinetic, potential and binding energies of individual particles -overcomes the entropy term so that an ordered state is favored, with the order-disorder transition being predominantly determined by the interparticle pairwise interaction potential energy [45].
The long-range Coulomb interaction of the pair of CIEs (trions) is strengthened by their permanent dipole moments directed perpendicular to the plane of the structure. Their actual exact interaction potential depends on the relative orientation of the triangles formed by the three charges in a trion complex. The exact potential includes nine terms to couple electrons and holes in two complexes by means of the V eff (R, ∆ρ X ± , r 0 ) and V C (R, ∆ρ X ± , d) potentials, where R is the trion center-to-center distance (see Methods).  minimum of the potential U in Fig. 3 (a) to represent the chain period and the unlike-charge trion coupling constant, respectively, one obtains T (N) Here, V (2R min ) ≈ 0 stands for the repulsive interaction coupling constant of the like-charge trions whose sublattice period is twice greater than the period of the chain.
The top and bottom panels in Fig. 3 (b) present the exact d-dependences of R min and U min calculated for the U-potential surface shown in Fig. 3 (a). Their approximate expressions can be relatively easily found analytically by seeking the U-potential minimum under the conditions r 0 , d < 1 and ∆ρ X ± > 1 consistent with Eq. (4). This leads to R min ≈ (r ee + r hh )/2 and U min ≈ −1/d + 1/r ee + 1/r hh , where r ee = (λ/σ)∆ρ X − and r hh = λ∆ρ X + are the interelectron and interhole distances in the negative and positive CIE, respectively. These expressions are seen to reproduce the numerical calculations quite well, within the approximations used, to demonstrate the fast drop of |U min | (and T cN for the unlike-charge trion crystallization transition, accordingly) with R min slowly rising as the interlayer separation d in the heterostructure increases.

C. Like-charge trion Wigner crystallization
In heterostructures of two monolayers with σ = 1 separated by an interlayer distance not tures, either X − or X + trions are most likely to form under intense irradiation. As can be seen from Fig. 2, for σ < 1 the domains |E X −| > |E X +| and |E X −| < |E X +| are located at smaller and greater d to form like-charge trions -negative and positive, respectively, as long as their binding energy absolute values exceed the thermal fluctuation energy at a given T .
An ensemble of repulsively interacting particles (or quasiparticles, structureless or compound) forms a Wigner lattice when its average potential interaction energy exceeds average kinetic energy, V / K = Γ 0 > 1. This was previously shown for systems such as 2D electron gas [46], cold polar molecules [47], and indirect excitons [4]. For like-charge trions in Fig. 1 (b), the Coulomb repulsion at large R is strengthened at shorter R by the dipole-dipole repulsion of their collinear permanent dipole moments (to result in the pairwise interaction potential V illustrated in Fig. 3), while the total kinetic energy is additionally contributed by the rotational term K (r) X ± =h 2 l(l + 1)/2I X ± with l = 0, 1, 2, ... being the orbital quantum number and I X ± = m h,e r hh,ee /2 representing the moment of inertia for CIE rotation about its permanent dipole moment direction. The low-T statistical averaging over l leads to the characteristic rotational motion "freezing" temperature T (r) X ± =h 2 /k B I X ± (see, e.g., Ref. [48]). By direct analogy with the hydrogen molecular ion problem this can be rewritten as T (r) X + ≈ σ|E X +|/k B and T (r) X − ≈ |E X −|/k B σ (see, e.g, Ref. [49]), indicating the rotational degrees of freedom to be frozen out (at least for the case of σ being close to unity typical of TMDs, in particular [50,51]) as long as the CIEs are stable against the thermal fluctuations.
With no rotational term contribution, it is straightforward to get a qualitative picture of the like-charge trion Wigner crystallization by performing an analysis analogous to that done in Ref. [46] for the 2D electron gas. With slight modifications to include the dipole repulsion in the interparticle interaction potential V and to replace the electron mass by the CIE mass in the translational kinetic energy K, the expressions for the zero-T critical density n c and for the critical temperature T (W) c of the Wigner crystallization phase transition take the form (see Methods) The quantities n cX ± and T (W) cX ± are shown on the top and bottom of Fig. 3 (c) as functions of d and σ (< 1), respectively, for moderate Γ 0 values [46]. As d increases so does V once the dipole repulsion becomes appreciable. With constant Γ 0 this leads to the K increase and n cX ± rise, accordingly. The latter is slightly lower for the negative CIE due to its smaller K because of the smaller mass than that of the positive CIE. Lowering σ generally lowers the CIE mass thus decreasing its K whereby T (W) cX ± increases. These are the general trends featured in Fig 3 (c).

D. Estimates for the effects discussed
We consider the case of the CIE formation in TMD homobilayers (both monolayers of the same material) encapsulated in bulk hexagonal boron nitride (hBN), a popular practical realization one encounters in a wide range of experiments [23,24,32,55]. Heterobilayers (two different TMD monolayers) offer many more CIE formation possibilities and therefore preferably should be analyzed individually. For the quantitative description of the effects predicted, our model requires the knowledge of the exciton reduced effective mass µ, the electron-hole effective masses m e,h associated with it, the effective average dielectric constant ε of the system, and the screening length parameter r 0 = 2πχ 2D with χ 2D being a spatially dispersive (and so nonlocal, i.e. in-plane distance-dependent) polarizability function [41,52].
We use µ, m e and m h reported recently from the first-principles calculations of the TMDmonolayer electronic structure [44]. The effective dielectric permittivity ε can be evaluated by the Maxwell-Garnett method [53], which in our case prescribes to use the weighted average of the hBN and TMD static permittivities, whereby for the hBN-monolayer number much greater than two we obtain ε = 5.87 (bulk hBN permittivity averaged over all three directions [54]). Finally, the r 0 parameter can be obtained based on the fundamental energy minimum principle [56], whereby the (negative) binding energy of a CIE complex must contribute the most in order for the CIE ensemble to be at a local minimum of its total energy in equilibrium. The r 0 parameter can therefore be found as the maximum point of the CIE binding energy absolute value |E X ± (σ, d, r 0 )| taken with both σ = m e /m h and d fixed.
We note that by its definition the KR potential screening length r 0 refers to in-plane charges which are the like-charge carriers to form the CIE in our case. These carriers are separated by distances at least of the order of 2a * B -much greater than those of the order of a * B one typically encounters in the exciton case. Therefore, being determined by greater distances, our r 0 due to its inherent nonlocality may very well be different from the values previously reported theoretically and experimentally for excitons in TMD monolayers [42,55]. With µ, m e,h and ε found as described above, we first calculate the exciton Bohr radius a * B and Rydberg energy Ry * for each case individually. Then, with known σ = m e /m h , a * B and Ry * we obtain the binding energy surfaces |E X ± (σ, d, r 0 )| in physical units from Eqs. (1), (2) and (4), determine their maximum points r X ± 0 for a particular fixed d, and compute the CIE binding energy absolute values |E X ± (σ, d, r X ± 0 )|. We do this for the interlayer distances d = 3, 4, 5 and 6Å (typical of van der Waals coupling) for each homobilayer type in order to be able to see the tendencies for the X + and X − trion formation as d increases. As an example, the left and right top panels in Fig. 4 show the X + and X − binding energy surfaces and their fixed-d crosscuts, respectively, for the MoSe 2 homobilayer. The vertical dashed lines on the right panel trace the |E X ± | maxima and their respective r X ± 0 distances. The CIE parameters thus obtained are tabulated at the bottom of Fig. 4 for all four homobilayers selected. We note the general consistency of our |E X ± | obtained both with numerical simulation data reported previously for the MoS 2 /WS 2 heterobilayer embedded in hBN (18/28 meV for the X + /X − trion [15]) and with the latest experimental observations on the MoSe 2 /WSe 2 heterobilayer system (10/15 meV for the X + /X − trion [23] and 28 meV for the X − trion [24], respectively). Highlighted greenish in the table are the largest differences between the positive and negative trion binding energies in MoSe 2 due to a significant m e and m h difference yielding σ = 0.8, which makes this homobilayer energetically favorable for the positive CIE Wigner crystallization for the interlayer distances d ranging between 3 and 5Å. As d increases from 3 to 6Å, for all types of bilayers tabulated, both |E X + | and |E X − | quickly decrease and get closer together while still remaining significant in magnitude, to make the normal unlike-charge trion crystallization energetically favorable. In the case of MoSe 2 , this implies a crossover from the Wigner crystal phase of the positive trions to the normal crystal phase of the unlike-charge trions.
A similar crossover from the Wigner crystallization of the negative trions to the normal crystallization of the unlike-charge trions, although not as pronounced as for MoSe 2 , might also be the case for WS 2 and WSe 2 according to our data tabulated. For MoS 2 , on the contrary, only the normal unlike-charge trion crystallization is energetically favorable as |E X + | and |E X − | there are about the same over the entire range of the interlayer distances d presented.
Using a * B and Ry * obtained as scaling units, it is quite straightforward to estimate the critical parameters for many-particle CIE systems in TMD homobilayers tabulated in cX ± ≈ 6 meV (to give T (W) cX ± ≈ 70 K). The fact of k B T (N) c being much greater than our |E X ±| tabulated tells that the dipoleordered normal 1D-crystal phase is the actual ground state of the many-particle unlikecharge trion system. The obtained n cX ± and T (W) cX ± are, respectively, close to and exceed those reported experimentally for IEs [23,31], suggesting that the Wigner crystallized CIE phase can be realized in properly fabricated vdW heterostructures with the twofold overbalance of negative [as in Fig. 1 (b)] or positive charge carriers. Crystallized exciton photoemission features can be found in Ref. [9].
In summary, we study the properties of charged interlayer excitons in highly excited vdW heterostructures -a compound fermion system with the permanent dipole moment observed recently in TMD bilayers [10,23]. We predict the existence of new strongly correlated collective CIE states, the long-range ordered phases of the excited heterostructure -the crystal phase and the Wigner crystal phase. We evaluate the critical temperatures and density for the formation of such many-particle cooperative compound fermion states. We demonstrate that they can be selectively realized with bilayers of properly chosen electron- respectively, as shown in Fig. 5 (a) for the negative trion case [34]. For such a quantum system the effective configuration space can be represented by the two independent in-plane projections ρ 1 and ρ 2 of the relative e-h coordinates (relative to the center of mass) of each of the IEs, whereby the X ± ground-state Hamiltonian takes the following form [33] The "atomic units" are used with distance and energy measured in the units of exciton Bohr radius a * B = 0.529Å ε/µ and Rydberg energy Ry * =h 2 /(2µ m 0 a * 2 B ) = e 2 /(2εa * B ) = 13.6 eV µ/ε 2 , respectively [36][37][38][39], µ = m e /(λ m 0 ) with λ = 1 + σ stands for the exciton reduced effective mass (in the units of free electron mass m 0 ), σ = m e /m h is the electron-to-hole effective mass ratio, and ε represents the effective average dielectric constant of the entire bilayer structure [36]. The image-charge effects are neglected.
The first two lines in Eq. (6) describe the kinetic and potential energy, respectively, for the two non-interacting IEs. Their individual e-h attractive Coulomb potentials screened, generically of the form (atomic units) with ρ being the in-plane intercharge distance, are symmetrized to account for the presence of the neighbor a distance ∆ρ away as seen from the ρ 1 -and ρ 2 -coordinate systems assigned to originate at the respective IE centers-of-mass and treated independently; see Fig. 5 (a). The last line is the interexciton exchange Coulomb interaction (or the like-charge Coulomb repulsion potential inside the trion) -h-h for X + and e-e for X − , respectively. We use the repulsive KR interaction potential to represent this interaction (atomic units) in order to properly take into account the screening effect for the like charges confined to the same monolayer [40]. Here, N 0 and H 0 are the 0th order Neumann and Struve functions, respectively, r 0 is the screening length defined in Eq. (3) for a 2D material [42], and ǫ 1,2 are the dielectric permittivities of its surroundings. To facilitate the analytical calculations, we approximate Eq. (8) by its accurate alternative (3) written in terms of elementary functions as discussed and proposed for atomically thin layers in Ref. [41].
For the CIE complex of two identical configurations with an extra charge attached to the left or right IE, the total wave function must be either symmetric or antisymmetric with respect to their interchange due to the conservation of parity. This can generally be achieved with coordinate wave functions of the form where φ IX (ρ 1 , ρ 2 ) = ψ IX (ρ 1 , d) ψ IX (ρ 2 , d) with ψ IX being the IE wave function. This involves the two terms localized at ρ 1 = ρ 2 = 0 and ρ 1 = −ρ 2 = ∆ρ, respectively, to represent the two equivalent configurations in terms of the two independent relative e-h coordinates ρ 1 and ρ 2 as shown in Fig. 5 (a). Since the total wave function of the quantum ground state must be nodeless [37], for large ∆ρ ≫ 1 the ground-state wave function of two IEs (two bosons) must be symmetric in coordinates to hold with Ψ g in Eq. (9). At shorter ∆ρ > ∼ 1 it can be multiplied by an even function of coordinates to be found from the Hamiltonian (6) in the manner similar to that developed in the past for the hydrogen molecule and molecular ion in seminal works by Landau, Gor'kov, Pitaevski, Holstein and Herring [37][38][39] and more recently by one of us for biexcitons and trions in quasi-1D/2D semiconductors [33][34][35].
Assuming further that for both configurations their respective IEs are in the spin-singlet states as dictated by the hyperfine interactions of their unlike-charge spin-1/2 fermionic (electron and hole) constituents [49], one arrives at the CIE complex featuring the ground state with two identical like-charge collinear-spin fermions in the same layer, which are thereby forced both by the Coulomb repulsion and by the Pauli exclusion principle to avoid each other at short ∆ρ < 1. Such a CIE complex is therefore only possible to form due to the asymptotic Coulomb exchange coupling at ∆ρ > ∼ 1, the domain our theory applies for. Figure 5 (b) shows a diagonal vertical crosscut of the potential energy surface (bottom) as given for X − by the second line of Eq. (6) in the two-coordinate configuration space (ρ 1 , ρ 2 ). On the main diagonal, this surface has two symmetrical minima separated by the potential barrier. The minima represent the two equivalent isolated IE states (top) given by the solution to the ground-state eigenvalue problem defined by the first two lines of the Hamiltonian (6). This solution is the product of the two ground-state IE wave functions.
The interlayer (or indirect) exciton eigenvalue problem was previously studied by Leavitt and Little [36]. Their ground-state energy E IX and the wave-function ψ IX are as follows (atomic units) where with N = 4/ 1 + 4 √ d + 8d (1 + √ d ) as per the normalization ∞ 0 dρ ρ |ψ IX (ρ, d)| 2 = 1. As described at large in Refs. [33,34], we start the CIE binding energy calculation with the (ρ 1 , ρ 2 )-configuration space transformation to the new coordinates as follows This transformation places the origin and both axes of the new coordinate system (x, y) as shown in Fig. 5 (b) -in the middle of the potential barrier that separates the two potential wells representing the two equivalent isolated IE states -to capture the maximal tunnel flow J X ± (∆ρ) between the two indistinguishable IE configurations. An approximate solution to the Schrödinger equation with the Hamiltonian (6) can be constructed using Eq. (11).
By converting Eq. (11) to the (x, y)-space per Eq. (12), we define the product wave function to describe the motion with the energy E IX inside the potential well centered at ρ 1 = ρ 2 = 0 (or x = −∆ρ/ √ 2, y = 0), while being exponentially damped outside. In just the same way, the function φ IX (−x, y) describes the motion with the same energy inside the well centered at ρ 1 = −ρ 2 = ∆ρ for the X − case shown in Fig. 5 (b) and at ρ 2 = −ρ 1 = ∆ρ for the X + case (both corresponding to x = ∆ρ/ √ 2, y = 0). Both of these functions are properly normalized to unity within their respective potential wells. Both of them are even in x and y with respect to their respective well-center positions, whereby ∂φ IX (∓∆ρ/ √ 2, y)/∂x = ∂φ IX (x, 0)/∂y = 0.
When the small probability of the underbarrier tunneling is taken into account, the energy level E IX splits into E IX − J X ± (∆ρ) and E IX + J X ± (∆ρ). Then, the correct zeroapproximation wave functions corresponding to these levels are and since φ IX (x, y)φ IX (−x, y) is vanishingly small everywhere, they are normalized so that the integrals of their squares over both wells are unity. This suggests that the actual eigenfunctions of the eigenvalues E g,u can be written as where Here the kinetic and potential energy terms are as followŝ In general, the function ψ X ± (x, y) is supposed to preserve the parity and the behavior of the function φ IX (x, y), to only depart noticeably from φ IX (x, y) in the very tail area x ∼ y ∼ 0 under the potential barrier and to overlap with ψ X ± (−x, y) in there; see Fig. 5 (b).
The overlap enables the tunnel exchange between the two indistinguishable configurations represented by φ IX (x, y) pinned to the potential well centered at ρ 1 = ρ 2 = 0 (x = −∆ρ/ √ 2, y = 0) and by φ IX (−x, y) pinned to the other potential well at ρ 1 = −ρ 2 = ∆ρ or ρ 2 = −ρ 1 = ∆ρ (x = ∆ρ/ √ 2, y = 0) for X − and X + , respectively. Under these restrictive assumptions about ψ X ± (x, y) in Eq. (14), it is possible to write down the two Schrödinger equations as follows whereT and U are those of Eq. (16). We multiply from the left the former by ψ g (x, y) and the latter by ψ X ± (x, y), subtract one from another, and integrate over x from −∞ to 0 and over y from −∞ to +∞. This includes the potential well positioned at x = −∆ρ/ √ 2, y = 0, and we find In here, with T of Eq. (16) it can be seen that its last term might only be significant at or close to x = −∆ρ/ √ 2, y = 0, but the partial derivatives of relevance are zero there, and so this term can be dropped for smallness over the entire integration domain. What remains can be integrated by parts. Bearing in mind that ψ g (0, y) = √ 2 ψ X ± (0, y), ∂ψ g (0, y)/∂x = 0 and all the functions involved as well as their derivatives must vanish at infinity, this after numerus cancelations gives From here, with just a tiny adjustment for practical application purposes, we obtain the tunnel exchange splitting integral in Eq. (1) of the following final form Here, we take into account the fast exponential drop-off of the integrand away from the y = 0-plane, whereby the integration limits can be shrunken to only include the physically significant cross-section region, see Fig. 5 (b), that controls the under-barrier tunnel probability flow -a positive quantity we wish to stress by taking the absolute value of. Such a tunnel exchange coupling binds the three-particle system to form a stable CIE state.

(a) The Trion Wave Function
We seek the function ψ X ± (x, y) of Eq. (17) in the following form Here, the unknown function S X ± (x, y) is to be chosen so that S X ± (x = −∆ρ/ √ 2, y) = 0 to fulfill the condition ψ X ± (−∆ρ/ √ 2, y) = φ IX (−∆ρ/ √ 2, y) as per Eq. (14), while also being smooth and slowly varying in the domain |x|, |y| < ∆ρ/ √ 2 under the barrier, whereby its second derivatives should be negligible. Additionally, as was mentioned above, for our threeparticle X ± complexes the equivalency of the two IEs sharing the same hole (or electron) implies their identity and leads to the fact of the like-charge carriers having collinear spins.
The Coulomb repulsion strengthened by the Pauli exclusion principle forces them to avoid each other at short interexciton center-of-mass-to-center-of-mass distance ∆ρ < 1, making it possible for a stable CIE complex to only form at ∆ρ > ∼ 1, which is why 1/∆ρ can be used as a smallness parameter in analytical calculations.
For the negative CIE, plugging Eq. (18) into the Schrödinger equation with the Hamiltonian (15), (16), to the first non-vanishing order in 1/∆ρ one obtains where the second-order derivatives of S X − are neglected. To find the analytical solution to this differential equation in the domain of interest |x|, |y| < ∆ρ/ √ 2, we use V eff of Eq. (3) to replce V KR in the right-hand side of Eq. (19). The solution to fulfill the boundary condition S X − (−∆ρ/ √ 2, y) = 0 is then given by ∆ρ). (20) To calculate the integral I(x, ∆ρ) here, we first use the unit step function to write followed by the change of variable τ = ( √ 2 t − σ∆ρ)/λ to obtain Of three terms here, only the third is seen to provide the solution in the domain x < σ∆ρ/ √ 2 that includes the region x ∼ 0 of interest to us. With V eff of Eq. (3), this term can be easily calculated analytically using integration by parts. One obtains with s = (σ∆ρ − √ 2 x)/λr 0 and p = ∆ρ/r 0 . A close inspection of this expression reveals that since s < p , the first summand is predominant there and the other two are negligible for all 1 < s < p regardless of how big s and p individually are. After dropping the negligible terms, Eq. (20) in the domain of interest takes the final form as follows For the positive CIE, plugging Eq. (18) into the Schrödinger equation with the Hamiltonian (15), (16) yields to the first non-vanishing order in 1/∆ρ the equation as follows It is easy to see that this equation can be obtained from Eq. (19) by the simple replacement 1/λ ↔ σ/λ. Its solution in the domain of interest can then be obtained by applying this replacement to Eq. (21). This gives

(b) The Tunnel Exchange Coupling Integral
It is noteworthy that both Eq. (21) and Eq. (23) are fully consistent with the result reported for σ = 1 previously [33]. The functions ψ X ± one obtains by plugging these equations into Eq. (18) can be used to evaluate the tunnel exchange coupling integrals J X ± in Eq. (17).
The differentiation therein can be conveniently done using the following easy-to-prove rule: Here ∂/∂(x, y) stands for either ∂/∂x or ∂/∂y. With this, after simplifications and elementary integration over y one obtains J X ± (∆ρ) in the form of Eq. (2) in the main text.
Seeking the extremum for J X ± (∆ρ) must only include the leading term in small 1/∆ρ to be consistent with the procedure of finding S X ± described above. Taking the derivative of J X ± over ∆ρ, equating it to zero, and solving the polynomial equation obtained to the first infinitesimal order in 1/∆ρ, results in ∆ρ X ± in the form of Eq. (4) in the main text.
(c) Remarks on the Interlayer Coulomb Interaction Potential The electrostatic interaction potential energies (7) and (8) we use in our analysis can be shown to consistently originate from the general solution to the electrostatic boundaryvalue problem that includes two coupled parallel monolayers. Such a solution was recently obtained by one of us (with coathors) as a byproduct in the bilayer optical probing experiment analysis (see Ref. [57], Appendix A). A bilayer system was considered to consist of the two parallel monolayers with individual 2D-polarizabilities χ ′ 2D and χ ′′ 2D (in our notations) that are separated by a distance d and surrounded by a dielectric medium of the static permittivity ε, with a point charge sitting at the origin of the cylindrical coordinate system placed in the bottom layer. In order to find the electrostatic interaction potential energy in the whole space, the Poisson's equation was solved in the Fourier space in the way similar to that reported in Ref. [41]. In the 2D-coordinate space, the solution obtained yields the electrostatic unlike-and like-charge interaction energies of interest as follows (atomic units) and r ′′ 0 = 2πχ ′′ 2D are the respective screening parameters for the individual monolayers. Due to the presence of the second layer, these equations do not seem to look similar to the solitary-monolayer KR potential case. However, setting d = ∞ to take the top layer away makes the former zero, while the latter integrates to yield the KR potential energy (8) with the effective screening length r 0 = r ′ 0 just as it should be. A close inspection of Eq. (24) reveals that due to the oscillatory behavior of the 0th order Bessel function J 0 (x) for all x > 1, only q < ∼ 1/ρ contribute the most to the integrals there.
In our case, ρ ≈ ∆ρ X ± as can be seen from Fig. 5 (a). Then, in the domain 1/∆ρ X ± < 1 we work within, only wave vectors q < ∼ 1/ρ ≈ 1/∆ρ X ± < 1 contribute the most to both integrals in Eq. (24), so that qd < ∼ d/∆ρ X ± < d < 1 in both integrals for all d we used in this work. This can also be seen from Fig. 6 we obtained using ∆ρ X ± of Eq. (4). Therefore, it is legitimate to neglect q 2 -terms under the integrals in Eq. (24). This gives , and the second integral turns into the KR potential energy (8) with the screening length Additionally, as per previous computational studies of monolayer TMDs [42], the monolayer screening length can be accurately represented by c(ε ⊥ −1)/2(ǫ 1 + ǫ 2 ), where c and ε ⊥ are the bulk TMD out-of-plane translation period and in-plane dielectric permittivity, respectively. For a TMD bilayer embedded in hBN with ε = 5.87 (averaged over all three directions [54]), which is the case for a variety of experiments [23,32,55], the typical parameters are c ≈ 12−13Å, ε ⊥ ≈ 14−17, ǫ 1 = (2ε ⊥ + ε )/3 with ε ≈ ε ⊥ /2 [42,54] and ǫ 2 = ε (or vice versa), to yield r 0 ≈ c(ε ⊥ − 1)/(5ε ⊥ /6 + ε)a * −1 B < 1 as a * B in TMDs is consistently greater than 1 nm both by our data (see Fig. 4) and also by others [32,42,55]. Then, we obtain r 0 /ρ ≈ r 0 /∆ρ X ± ≪ 1. With this in mind the denominator of the first integral above can be expanded in rapidly convergent binomial series, whereby after the term-by-term integration the interlayer electrostatic interaction energy takes the form Here, the second term in parentheses comes out as the 2nd (not the 1st as one would expect!) order of smallness since d/ρ ≈ d/∆ρ X ± < 1 as demonstrated in Fig. 6, and so it can be safely dropped along with the rest of higher infinitesimal order terms, whereby one arrives at the interlayer Coulomb interaction (7) we used in our calculations throughout this work. Note also that, even more generally, this series expansion can be seen to be uniformly suitable for all ρ ≥ 0, including ρ ∼ 0 as well, in which case the second term in parentheses comes out as As can be seen from the two special cases shown in Fig. 7 (a) and (b), the long-range Coulomb interaction of the pair of CIEs (trions) depends on the relative orientation of the triangles formed by the three charges in a trion complex. The exact interaction potential includes nine terms to couple the electrons and holes in the two spatially separated complexes.
To simulate the actual potential energy surfaces we use the Coulomb interaction coupling of Eq. (7) for the (unlike) charges located in the distinct monolayers and the KR interaction coupling of Eq. (8) for the (like) charges confined to the same monolayer. The explicit coupling parameter dependence is given by the functions V KR (R, r 1 , r 2 , r 0 ) and V C (R, r 1 , r 2 , d) specified below, where R is the trion-trion center-of-mass-to-center-of-mass distance and r 1,2 are the distances between the like charges in the first and second trion of the interacting trion pair. In general, r 1 = r 2 for the unlike-charge trion-trion coupling and r 1 = r 2 for the like-charge trion-trion coupling as sketched in Fig. 7 (a) and (b). Using the standard triangle similarity theorems, these distances come out as λ∆ρ X + and (λ/σ)∆ρ X − for the positive and negative trion, respectively. This is why for unlike-charge trion pairs, r 1 can only be equal to r 2 if σ = 1 (or m e = m h ).

(a) The Pairwise Interaction Potentials for Unlike-Charge IEs
In this case, two most likely relative orientations are supported by symmetry for a pair of triangle-shaped complexes in bilayer structures we deal with. They are the coplanar and parallel biplanar orientation. Their side and top views are shown in Fig. 7 (a) and in the bottom-left inset of Fig. 8, respectively. For the former, counting e-h couplings in Fig. 7 (a) counterclockwise from top left, the total interaction potential energy U takes the form For the latter, from the inset in Fig. 8 the total interaction potential W comes out as The calculated interaction potentials U and W are presented in Fig. 3 (a) of the main text and in Fig. 8 herewith, respectively. The former is seen to be over an order of magnitude more attractive than the latter in the same parameter range, which is why the W interaction potential energy is neglected in the analysis we report about in the main text.

(b) The Pairwise Interaction Potentials for Like-Charge IEs
In this case, both coplanar and parallel biplanar relative orientations of the triangleshaped complexes are strongly repulsive and, in general, are different for positively and negatively charged trion pairs. The side view of the coplanar orientation of two positive trions is shown in Fig. 7 (b). The top view of their parallel biplanar orientation can be obtained from the sketch in Fig. 8 by setting r 1 = r 2 and relabeling e ↔ h in one of the trions. For the former, counting e-h couplings in Fig. 7 (b) counterclockwise from top left, the total interaction potential energy V takes the form For the latter, the total interaction potential energyV can be obtained from Eq. (26) by setting r 1 = r 2 and simultaneously swapping V KR ↔ V C and R 2 ↔ R 2 + d 2 . This gives For a negatively charged trion pair, r 2 should be replaced with r 1 in both of these equations.
A close inspection of Eqs. (27) and (28) reveals their very similar repulsive behavior and in fact their coincidence when R is greatly different from r 2 (both greater and less than).
The calculated interaction potential V of Eq. (27) is presented in Fig. 3 (a) of the main text.

C. Like-charge trion Wigner crystallization parameters
An ensemble of repulsively interacting particles (or quasiparticles, structureless or compound) forms a Wigner lattice when its average potential interaction energy exceeds average kinetic energy, V / K = Γ 0 > 1 (see, e.g., Ref. [46]). For like-charge trions in Fig. 7 (b), the Coulomb repulsion at large R (≫ r 2 ) is strengthened at shorter R by the dipole-dipole repulsion of their collinear permanent dipole moments directed perpendicular to the heterostructure plane. These are the two major terms of the power series expansion in r 2 /R (< 1) of the repulsive pairwise interaction potential V presented in Fig. 3 (a) of the main text.
With rotational kinetic energy neglected for the reasons explained in the main text, the like-charge trion critical density n cX ± and temperature T (W) cX ± can be obtained by drawing an analogy to the 2D electron gas system [46] to include the extra dipole-dipole repulsion term.

(a) The Critical Density
With the commonly used notations preserved, we go on with using the atomic units introduced previously. For trion-trion separation distances R greater than the size of the trion (R ≫ r 1,2 ), the first order power series expansion of the average repulsive trion-trion interaction potential takes the form where n = 1/πR 2 is the trion surface density. Our trions are compound fermions with the occupation number where β = 1/k B T , E k =h 2 k 2 /2M, M = M X ± andμ being the trion total mass and chemical potential, respectively. At zero T this turns into a unit-step function to give n in Eq. (29) in the form where S is the surface area and k F is the trion Fermi-momentum. The average kinetic energy per particle can then be written as to result, with V of Eq. (29), in where g stands for the ratio of the electron-hole reduced mass to the trion total mass Introducing the new variable t = d √ πn turns Eq. (33) into a quadratic equation with two roots as follows of which only one, t 2 , stays finite as d goes down to zero. This root leads to and reproduces the result of Ref. [46] for d → 0 and g ± = 1.

(b) The Critical Temperature
For arbitrary nonzero T , using Eq. (30) with the new variable x =hk β/2M, the trion surface density (31) can be written in a parametric form as follows Similarly, the average kinetic energy per particle of Eq. (32) takes the form After the power series expansions of their respective denominators, these integrals can further be represented in terms of the gamma and polylogarithm functions following the rule In the classical limit (high T and/or low density; see, e.g., Ref. [56]), one has e β(E k −μ) ≫ 1, so that the occupation number (30) takes the form n k = e −β(E k −µ) = ze −βE k to simplify n in Eq. (36) as follows whereby the kinetic energy per particle of Eq. (37) takes the form as expected from the energy equipartition theorem of classical statistical mechanics.
In this equation, to make it consistent with the approximation Eq. (29) is valid within, one has to discard the terms with powers of d higher than d 2 . The quadratic equation thus obtained gives two roots for n(T ), one of which is manifestly negative and so to be discarded.
Equating the other root to n cX ± of Eq. (35) gives the constraint for the critical temperature.
Solving it for T subject to keeping powers of d no greater than d 2 , leads to k B T (W) cX ± = 4 g ± Γ 2 0 (in the units of Ry * ) with g ± (σ) given by Eq. (34). This agrees with Ref. [46] for g ± = 1. in-depth assessments of exciton and trion interaction potentials. Y.E.L. provided expertise in exciton many-particle correlations and crystallization phenomena. All authors discussed the results and commented on the ways to best represent them in the manuscript.