{'fraction_non_alphanumeric': 0.04637769020915429, 'fraction_numerical': 0.02503789026977872, 'mean_word_length': 4.107120743034056, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electron-Hole Liquid in Monolayer Transition Metal Dichalcogenide Heterostructures

Monolayer films of transition metal dichalcogenides (in particular, MoS2, MoSe2, WS2, and WSe2) can be considered as ideal systems for the studies of high-temperature electron-hole liquids. The quasi-two-dimensional nature of electrons and holes ensures their stronger interaction as compared to that in bulk semiconductors. The screening of the Coulomb interaction in monolayer heterostructures is significantly reduced, since it is determined by the permittivities of the environment (e.g., vacuum and substrate), which are much lower than those characteristic of the films of transition metal dichalcogenides. The multivalley structure of the energy spectrum of charge carriers in transition metal dichalcogenides significantly reduces the kinetic energy, resulting in the increase in the equilibrium density and binding energy of the electron-hole liquid. The binding energy of the electron-hole liquid and its equilibrium density are determined. It is shown that the two-dimensional Coulomb potential should be used in the calculations for the electron-hole liquid.


INTRODUCTION
An increased interest in the studies of graphene as a promising material for nanoelectronics [1] has led to the development of novel two-dimensional (2D) materials, such as monolayers of hexagonal boron nitride, black phosphorus, and transition metal dichalcogenides (TMDs) [2]. Recently, vertical (van der Waals) heterostructures, where various 2D materials are combined in a given sequence, have been actively studied [3].
Monomolecular layers of TMDs with the chemical formula M X 2 , where M is a transition metal and X is a chalcogen, are of a particular interest. Among them, semiconductors with Group VI metal atoms (M = Mo, W) and S, Se, and Te chalcogens are the most studied. Bulk layered TMDs (e.g., MoS 2 , WS 2 , MoSe 2 , and WSe 2 ) have an indirect energy gap E g ∼ 1 eV [4,5], whereas the corresponding monomolecular layers are direct-gap semiconductors with E g about 2 eV [6].
Many bulk TMD samples were obtained as early as the 1960s [7]. Even at that time, their electron properties were already actively investigated [8,9]. In particular, some TMDs (M = Nb, Ta, Ti, Mo; X = S, Se) exhibit a low-temperature transition to a superconducting state. The structure, synthesis, properties, and applications of TMDs are described in detail in the recent review [10].
The optical properties of monomolecular TMD layers are mainly determined by excitons and trions. The exciton binding energy Ex in TMDs is of the order of hundreds of meV (e.g., in MoS 2 monolayers, E x = 420 meV [11]), and the corresponding energy for trions is equal to tens of meV [6].
These facts suggest that the arrays involving TMD monolayers are ideal systems for the studies of a hightemperature electron-hole liquid (EHL). The energy of one electron-hole pair in the EHL is |E EHL | ∼ E x and the critical temperature of the gas-liquid phase transition is T c ∼ 0.1 |E EHL | [12][13][14][15][16][17]; therefore, it can be expected that the EHL will be observed in TMD monolayers even at room temperature. In MoS 2 monolayers, a high-temperature strongly coupled EHL with T c ≃ 500 K has been already observed [18].
In this work, we investigate the possibility of forming the EHL in monolayers of multivalley semiconductors [9,19]. We consider a thin film of a model multivalley semiconductor deposited in vacuum on a dielectric substrate. We assume that the semiconductor has a fairly wide band gap E g ≫ |E EHL | and use the single-band approximation. The semiconductor has a large equal number of equivalent electron ν e and hole ν h valleys, ν e = ν h = ν ≫ 1 with the effective masses m e and m h of the electron and hole, respectively. The multivalley structure can be due to the presence of sev-eral molecular monolayers in the film.
In [20], it was shown that the interaction energy of charge carriers belonging to different valleys dominates in such a system at ν ≫ 1. The equilibrium density n EHL of the EHL and the energy E EHL corresponding to this density increase drastically. This increase in the density justifies the use of the random phase approximation to calculate the correlation energy.

MODEL
We study the 2D electron-hole system described by the Hamiltonian [21,22] Here, a psk (a † psk ) and b psl (b † psl ) are fermion annihilation (creation) operators for an electron and a hole having the crystal momentum p and spin projection s, which are located in kth and lth valleys, respectively. The dispersion relations of electrons and holes are ε e sk (p) = where σ = m e /m h and m = m e m h /(m e + m h ) is the reduced mass of the electron and hole. Usually, σ ≤ 1. The Coulomb interaction in the films of a finite thickness is described by the Keldysh potential [23,24] V (q) = 2π e 2 |q|(1 + r 0 |q|) .
Here, e 2 = e 2 /ǫ eff , where ǫ eff = (ǫ 1 + ǫ 2 )/2 is the effective permittivity of the media adjacent to the film (e.g., ǫ 1 = 1 is the permittivity of vacuum and ǫ 2 is the permittivity of the substrate), r 0 = d/2δ is the screening length, δ = ǫ eff /ǫ (ǫ is the permittivity of the substrate material), and d is the film thickness. For the monolayer film (d → 0), Eq. (3) leads to the following conventional expression for the Coulomb interaction of charge carriers in the 2D system: In terms of macroscopic electrodynamics, the screening of the Coulomb interaction of charge carriers is determined by the permittivities of media surrounding the film, since the electric field lines turn out to be outside the film. The introduction of the permittivity ǫ for monolayer TMD films, as well as for graphene [1], has no physical meaning.
The nonzero term r 0 |q| in the denominator of Eq. (3) is introduced in order to explain the significant deviation of the energies of several first exciton levels from the Rydberg series [6].
First, we use potential (4) (the beginning of Section 3, Subsections 3.1 and 3.2). Then, to verify the validity of using potential (4), as well as to compare our results to experimental data, we carry out calculations using the Keldysh potential (Subsection 3.3). The calculation results for the monolayer MoS2 film with both potentials are compared in Section 5.

GROUND STATE ENERGY
The ground state energy per electron-hole pair in the 2D EHL can be written as [21,22] where is the average kinetic energy, where n 2D is the 2D electron and hole density; is the exchange energy; and E cor is the correlation energy defined below. Here, r s = ν/πn 2D is the dimensionless interparticle distance. The Fermi wave vector is q F = 2πn 2D /ν = √ 2/r s . Here and further on, we use the system of units in which the binding energy and radius of the 2D exciton are equal to unity: E x = 2m e 4 / 2 = 1 and a x = 2 /2m e 2 = 1.
The main problem in determining the ground state energy of the EHL is the calculation of the correlation energy. In the simplest case of a single-valley semiconductor, it was calculated in [21,22] by the Nosieres-Pines method. It was shown that, in contrast to the three-dimensional case, the 2D EHL appears to be more favorable in energy than the exciton gas even in the isotropic case. In this situation, the main contribution to E cor comes from the momentum transfer exceeding the Fermi momentum.
The calculations of the EHL energy have been recently reported in [25]. The correlation energy of the electron gas in narrow-gap multivalley and layered semiconductors is calculated in [26,27]. The Wannier-Mott excitons in heterostructures composed of narrow-gap semiconductors are addressed in [28].
In [29], the EHL in double quantum wells with spatially separated electrons and holes in multivalley semiconductors is studied. The energy and equilibrium density of the EHL are calculated at various distances between the electron and hole layers. The procedure for calculating the correlation energy of the 2D EHL at the spatial separation of electrons and holes is described in [30].

Calculation of the Correlation Energy at a Finite Number of Valleys
The correlation energy can be represented in the form of an integral over the momentum transfer [21,22,31,32] At small (compared to q F ) q values, the function I(q) is calculated in the random phase approximation, whereas at large q values, it is determined by the sum of diagrams in the second order of perturbation theory in terms of the interaction energy. For any value of σ, the expansion of I(q) at small q is rather lengthy. We show here the result for the particular case of equal masses of electrons and holes (σ = 1) In the intermediate range q 1 ≤ q ≤ q 2 , the function I(q) is approximated by a segment of the tangent, as in [21,22]. At q ≪ 1, we integrate expansion (9) from zero to the matching point q 1 ≈ q 0 (q 0 is the point corresponding to the minimum of the function I(q)), whereas the asymptotic expression at q ≫ 1 is integrated from the matching point q 2 to infinity. Adding the contribution to the integral from the intermediate range (I(q 1 ) + I(q 2 )) (q 2 − q 1 )/2, we find where . At ν ≤ 3 and 1 r s 2,   10) and (11), (red stars) E cor calculated by Eq. (10) with the numerical solution of the equation for q 0 , and (black circles) E cor calculated taking into account corrections of the first and second orders with respect to the deviation of q 1 from q 0 (q 0 was determined numerically). The lower green horizontal straight line corresponds to estimate (12) and purple triangles are E gs values calculated by Eqs. (5)- (7) and (10). The inset shows the q 0 values (blue squares) calculated by Eq. (11) and (red stars) calculated numerically. Note that q 0 → 2 √ 2 at ν ≫ 1.
Here, we take into account that, at ν ≫ 1, the position of the minimum of the function I(q) appreciably deviates from unity (q 0 → 2 √ 2). The comparison of the dependence of the correlation energy on the number of valleys calculated by Eq. (10) with that obtained taking into account the first and second order corrections with respect to the deviation of q 1 from q 0 is shown in Fig. 1 for the case of σ = 1 and n 2D = 1/π (r s = √ ν). It is noteworthy that the corrections to E cor only slightly affect the results (red asterisks almost coincide with black points). At large ν values, the correlation energy tends to estimate (12). The same figure shows the numerically calculated dependence of the energy of the ground state E gs on ν.

Calculation of the Correlation Energy in the Limit of a Large Number of Valleys
At ν ≫ 1, when n 2D satisfies inequalities [20,21] 1 ≪ q F ≪ n the correlation energy is given by the expression Here, where is the polarization operator in the zeroth approximation with respect to the interaction constant at high momentum transfers (q ≫ q F ) and frequencies (ω ≫ E F ). The dispersion relations are the same as in Eq. (2). Relation (14) can be easily expressed in terms of the dimensionless units q = (4πn 2D λ) 1/3 ξ and ω = (4πn 2D λ) 2/3 ζ. The correlation energy at an arbitrary ratio σ of the electron and hole masses can be written as where Figure 2 shows the numerically calculated function A(σ). It is convenient to fit it by the expression A(σ) ≈ 0.23 σ 2/3 e −4σ − 0.098σ 3 + 0.378σ 2 − 0.442σ + 4.932.
For σ = 1, we have where Γ(x) is the Euler gamma function. This value is close to the factor in estimate (12). Note that estimate (12) in the limit ν → ∞ remains a lower bound for (12). The constant A(1) is smaller than the factor in estimate (12) because the asymptotic expression (15) incompletely includes the contribution of low momenta and frequencies.

Calculation of the Ground State Energy
Using the Keldysh Potential The average kinetic energy is also given by Eq. (6). The exchange energy is expressed in the form (ρ 0 = r 0 q F ) where J(ρ 0 ) =  At characteristic densities n 2D ∼ 10 13 − 10 14 cm −2 , the dimensionless parameter is ρ 0 ≃ 2 − 9 at ν = 2. Therefore, it is interesting to calculate numerically the function J(ρ 0 ) in the range 0 < ρ 0 < 10. In the limit of large ρ 0 values, it tends to 8/3 (see Fig. 3) We determine the correlation energy by the method described in Subsection 3.1. Now, instead of the function I(q), the integral in Eq. (8) contains the function In the intermediate range of wave vectors q, the function I(q) can also be fitted by a straight line segment. Integrating I(q) over q, we obtain an expression for the correlation energy E K cor . It is rather lengthy and is not presented here. However, for large ν values, when ρ 0 q 1 ≪ 1 and ρ 0 q 2 ≪ 1 (q 1 and q 2 are the matching points and q 1 ≈ q 0 , where q 0 is the point of minimum of the function ), this expression can be expanded in powers of ρ 0 as where E cor is given by Eq. (10), whereas the correction in the approximation linear in ρ 0 reads

EQUILIBRIUM DENSITY AND ENERGY OF THE ELECTRON-HOLE LIQUID
The equilibrium density of the EHL n EHL is found by the minimization of the ground state energy. Substituting the correlation energy given by Eq. (16) for a multivalley system into Eq. (5) and taking the derivative with respect to n 2D , we obtain the following equation for n EHL : To solve Eq. (20), we note that the absolute value of the exchange energy given by Eq. (7) at ν ≫ 1 is less than the absolute values of the kinetic and correlation energies. Therefore, we can first neglect the second term in Eq. (20) and then find the exchange energy correction to the equilibrium density: In our opinion, multilayer multivalley systems are promising (the number of layers is l and the distance between them is c). If d = lc a x , we can neglect the second term in the denominator of Eq. (3). In this case, Eq. (4) can be used. The effective number of valleys ν eff = lν of the system increases significantly.
The equilibrium density and energy of the EHL calculated using Eqs. (10) and (16) for the correlation energy in the case of 10 separated TMD monolayers (ν eff = 20) are shown in Fig. 4 as a function of the ratio of electron and hole masses.

COMPARISON OF THE CALCULATIONS WITH EXPERIMENTAL DATA
Let us compare the results calculated by Eqs. (5)- (7) and (10) with the experimental data for n EHL , |E EHL |, and T c obtained for the monolayer MoS 2 film [18]. With a good accuracy, we can assume that m e ≈ m h [33]. If ν = 2, we have |E EHL | = 700 meV, n EHL = 10 14 cm −2 and T c ≃ 800 K. The corresponding experimental data are |E EHL | = 480 meV, n EHL = 4×10 13 cm −2 , and T c ≃ 500 K [18]. We can give two explanations of this discrepancy.
The first explanation. The number of valleys decreases. Stresses arising in a monolayer film can lift the degeneracy of valleys [16,17]. In addition, the lifting of spin degeneracy of charge carriers is also equivalent to halving the number of valleys. This is possible due to the large spin-orbit splitting of the valence band, ∆ vb ≈ 148 meV [25]. For ν = 1, we obtain good agreement with the experiment: |E EHL | = 450 meV, n EHL = 3.3 × 10 13 cm −2 , and T c = 520 K.
The second explanation. It is necessary to use the Keldysh potential, which contains the fitting parameter r 0 . At r 0 = 0.7Å and ν = 2, the calculated EHL energy |E EHL | = 480 meV is in the best agreement with its experimental value, but the density of the EHL is overestimated: n EHL = 5.4 × 10 13 cm −2 .
In the quantitative description of the positions of exciton lines in the photoluminescence spectrum of the MoS 2 monolayer, we find r 0 = 41.47Å [34]. The discrepancy in the values of r 0 is large because the calculations of excitons and EHL involve quite different ladder and loop diagrams, respectively.
We favor the first explanation.

CONCLUSIONS
To summarize, we have obtained analytical and numerical results for the binding energy of the EHL and its equilibrium density in 2D systems with TMD monolayers at an arbitrary number of valleys.
We have also calculated the characteristics of the EHL using the Keldysh potential. It turned out that using only one parameter does not simultaneously match the binding energy of the EHL and its equilibrium density with experimental results. This suggests a quite limited range of applicability for the Keldysh potential in these calculations.
The difference between the theoretical and experimental results is due to the insufficient accuracy of the used parameters of the monolayer heterostructures, to the use of experimental E x and a x values as the units of measure, and to the possible inhomogeneity of the sample.