{'fraction_non_alphanumeric': 0.0358729074137342, 'fraction_numerical': 0.012299282541851725, 'mean_word_length': 4.348255070345331, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
To the problem of electron-hole bound state in transition-metal dichalcogenides

The interacting electron and hole in transition-metal dichalcogenides is considered. For investigation of the interaction between electron and hole we obtain the Bethe-Salpeter equation for two interacting Dirac particles. The dependence of a few lowest binding energies of electron and hole on the interaction constant for different potentials is found. We demonstrate that the behavior of the potential at small distances significantly affects on the values of the binding energies. For small interaction constant we have developed the perturbative method of the binding energy calculation. For the largre interaction constant the binding energies are found numerically. The critical values of the interaction constant for the Coulomb potential and exponential potential are found.

The interacting electron and hole in transition-metal dichalcogenides is considered. For investigation of the interaction between electron and hole we obtain the Bethe-Salpeter equation for two interacting Dirac particles. The dependence of a few lowest binding energies of electron and hole on the interaction constant for different potentials is found. We demonstrate that the behavior of the potential at small distances significantly affects on the values of the binding energies. For small interaction constant we have developed the perturbative method of the binding energy calculation. For the largre interaction constant the binding energies are found numerically. The critical values of the interaction constant for the Coulomb potential and exponential potential are found.

I. INTRODUCTION
In 2004 the famous article devoted to the study of the first two-dimensional material appeared 1 . This article opens a new direction of experimental and theoretical investigations of two-dimensional materials. After creation of the graphene many different two-dimensional and quasi two-dimensional materials were investigated both experimentally and theoretically. The interesting family of the materials is the transition-metal dichalcogenides. Some of the materials are the semiconductors which have the value of the band gap of the order of 2 eV , see eg. 2 . The investigation of the energies of the exciton states in such materials is an interesting problem, from both experimental and theoretical point of view. The experimental investigation of the exciton spectrum and comparison of the spectrum with the one predicted by the theory gives the understanding of the nature of the electron-electron interaction in the considered materials. It allows us to predict the properties of the similar materials. From theoretical point of view the problem of electron-hole interaction is interesting due to the non-triviality of including the interaction to the model. There were many attempts to include the interaction to the system, see e.g. [4][5][6][7][8] . In the Ref. 6 to obtain the Hamiltonian of the system of interacting electron and hole the authors perform the transformation of the Hamiltonian for electron-hole system without interaction to the block diagonalized form, Then the authors expand obtained Hamiltonian in the assumption that the kinetic energy is much less than the band gap. After that the Coulomb interaction was added to the expanded Hamiltonian. In Ref. 7 the Hamiltonian of the system was chosen as the sum of three terms. First term corresponds to the free particle Hamiltonian with reduced mass of electron and hole. Second term corresponds to the interaction potential, and the third therm corresponds to the corrections related to the Berry curvature. In Ref. 4 to consider the exciton spectrum and to investigate the critical behavior of the system the two-dimensional modified Dirac equation was considered. In the modified equation the momentum operator is doubled and the Coulomb field is added. In Ref. 5 to describe the bound state of electron and hole the three types of equations were considered. First equation contains sum of free particle Hamiltonian with reduced mass and interaction potential. The second term is related to the Berry curvature. Second type of the equation corresponds to the expansion of the Hamiltonian which is sum of two Hamiltonian for free Dirac particles and interaction potential in the case when the binding energy is much less than band gap. The third equation is some modification of the Bethe-Salpeter equation (BSE).
Different approaches to the problem demonstrate the difficulty of the including the interaction between electrons to the system. The difficulty is related with the appearance of the electron-hole excitations in the intermediate states 13,14 . The problem was solved in the framework of the quantum field theory. The existence of the bound states of the electron and hole manifests itself in the form of singularities in the two particle scattering amplitude. The equation which describes the singularities of the amplitude is the BSE. Therefore, to find the exciton spectrum it is necessary to obtain the BSE for the system of interacting electron and hole. The BSE for the semiconductors was obtained in Ref. 9 . In Ref. 10 the authors solved numerically the BSE and found the exciton spectrum. To obtain the spectrum the authors took into account the random phase approximation for the potential and the corrections related with the self-energy operator. However, in Ref. 10 the one particle Green's function corresponds to a particle with parabolic dependence of the energy on the momentum. In transition-metal dichalcogenides the one particle excitations are described by the two-dimensional Dirac equation. Therefore, the wave function should have two component and the dependence of the energy on the momentum is not parabolic. The approximation performed in Ref. 10 means the smallness of the kinetic energy of the particle in comparison with band gap. For excitons the kinetic and interaction energies are of the order of band gap, therefore such approximation is not correct. Thus, further study of the BSE for the interacting electron and hole in the transition-metal dichalcogenides is necessary.
In the present paper we obtain the BSE for the system of intravalley interacting electron and hole in the transitionmetal dichalcogenides in the leading order in the interaction potential. We investigate the solutions of the obtained equation for different values of the interaction constant and for different types of the potentials. For the Coulomb potential we find the analogue of the non-relativistic approximation for the equation which can be used in the case of small interaction constant. We demonstrate that this equation differs from that obtained in Ref. 5 . We also find numerically the critical value of the interaction constant, i.e. the value of the interaction constant at which the ground state energy of the interacting electron and hole reach the valence band. To investigate the behavior of the solutions in the vicinity of the critical values of the interaction constant, we consider the localized potential which does not contain any singularities.
The paper is organized as follows. In the Section II we describe the model, obtain the BSE for the electron and hole. In the Section III we consider the approximation of the BSE at small interaction constant. We demonstrate that standard method of the expansion of the BSE for the Coulomb potential leads to the appearance of the non-integrable corrections to the interaction potential III A. In the Subsection III B the correct expansion of the BSE is obtained. In the Section IV the exact numerical solution for the BSE is obtained. In this section we compare the exact solution and the one obtained using the perturbation theory. We also consider the critical values of the potential. In the Conclusion we discuss the obtained results.

II. MODEL
The starting point of our consideration is the one particle Hamiltonian suggested in 3 : where (σ x , σ y , σ z ) are the Pauli matrices, p = (p x , p y ) is the momentum operator, s z is the Pauli matrix for spin, τ is the valley index, the parameter λ denote the spin-orbit coupling parameter. In Eq. (1) we introduced the Fermi velocity parameter v F = at/ , a is the lattice constant, t is the effective hopping integral, is the Planck constant. The values of the parameters ∆, λ, t, and a can be found in the Table 1 of Ref. 3 . Below we set = v F = 1. Let us consider the eigenvalues ǫ τ,sz of the Hamiltonian (1) for different parameters τ and s z . We have So, for the valley corresponding to τ = 1 the gap between conduction and valence band is equal to ∆ − λ. The conductance band has lower bound which corresponds to the ǫ τ,sz correspond to the anti-particles (holes). One can check that the ratio λ/∆ ≪ 1, therefore below in the present paper we will consider the case λ = 0. In this case the masses m of the quasiparticles equal to ∆/2.
To find the energy spectrum of excitons (bound state of an electron and hole) it is necessary to introduce the interaction in the Hamiltonian (1).
In Refs. 4-8 the starting point of consideration is the expansion of the Hamiltonian in the vicinity of the electronhole pair energy which equals to 2m. Then the authors add the interaction potential V (r), and some terms related with Berry curvature, see, e.g., Refs. 5,7 . Such consideration is similar to the non-relativistic expansion in quantum electrodynamics. It is well known that this expansion is valid when the ratio (2m − E exc )/m ≪ 1, where E exc is the electron-hole pair energy. However, as can be seen from experimental data 11,12 , the ratio (2m−E exc )/m is of the order of unity for the ground state for such materials as W S 2 and W Se 2 . It means that the effective interaction constant between quasiparticles is of the order of unity. For such interaction parameter the analogue of the non-relativistic expansion is not applicable because the electron-hole excitations in the intermediate state can significantly change the interaction, see 13,14 . Therefore applicability of this expansion is disputable.
In quantum field theory the method of the energy calculation for the bound states is developed. It based on the investigation of the singularities in the two-particle scattering amplitude. The BSE 14 describes singularities in the two particle scattering amplitude which related to the bound states of electron and hole (electron and positron in quantum electrodynamics), therefore we should obtain the equation for the electron-hole system in the transition metal dichalcogenides.
For simplicity, below we consider the case τ = 1 and neglect spin-orbit interaction. In this case the Hamiltonian (1) has the form:Ĥ where σ = (σ x , σ y ), p = (p x , p y ). In the Hamiltonian (4) we omit spin index, since it diagonal in the spin variable.
Using the technique described in Ref. 14 we derive the BSE for interacting electron and hole in the leading order in the interaction potential: where V (q) is the Fourier transform of the electron-electron interaction potential V (r), G i,m a (ε, p) is the one-particle Green's function: where index a enumerates electron and hole. The imaginary part of the Green's function corresponds the Fermi energy equals to zero (center of the band gap). The function Ψ depends on the two energies ǫ 1 and ǫ 2 , but the wave function of an electron-hole pair should depend only on one energy parameter. So the quantity Ψ can not be treated as the wave function. To obtain the equation for the wave function we use standard prescription 14 . We change the variables ε 1 = E/2 + Ω, ε 2 = E/2 − Ω and perform the integration over Ω in the both sides of Eq. (5). As result we obtain: Here we introduce the wave function of the electron-hole pair: Below we consider the case p 1 = −p 2 = p. Performing the integration over Ω in Eq. (7), we finally obtain the equation for the wave functionψ(E, p, −p) = ψ(E, p): where we omit the indexes for σ matrices, the operators Λ ±± have the form: Here ω(p) = p 2 + m 2 . One can check that the equation (8) is similar to the BSE for electron-positron pair in quantum electrodynamics 14 . This is not surprising since the Hamiltonian (4) is the two-dimensional analogue of the Dirac Hamiltonian of the electron. In the coordinate space the Eq. (8) has the form: is the momentum operator, and Q ±± (R) is a Fourier transform of Λ ±± (p), The Eq. (11) has the conserving density ρ: Where ψ ±± (r) = dR Q ±± (r − R)ψ(t, R). The density ρ can be negative; therefore, we can not treat it as the probability density, but we can consider it as the quantity which proportional to the charge density, see Ref. 14 .

III. SMALL INTERACTION CONSTANT APPROXIMATION
Let us demonstrate that the BSE (8) can be transformed to the Schrödinger equation in the case of small interaction constant. It is convenient to present the function ψ(E, p) in the following form: where |Σ, Σ z is the eigenvectors of the operators (σ e +σ h ) 2 and (σ z e +σ z h ): We substitute the function ψ in the form (14) to the Eq. (8) and obtain the system of equations: where p ± = |p|e ±iφp , φ p is the angle of the vector p, One can see, that the component g(p) is equal to zero in the case of the nonzero energy E; therefore, below we omit the component g from the consideration. The expansion of the system (15)-(17) at small interaction constant is similar to the non-relativistic expansion of the Dirac equation in quantum electrodynamics 13 . For small interaction constant (the smallness of the potential will be considered latter) we imply that the bound states energies E obey the condition |E − 2m| ≪ m. For such energies we imply that the characteristic values of the momentum and potential obey the conditions 13 : In the inequality (21) we also imply that the main contribution to the integrals in Eqs. (15)-(17) comes from the momentum scales p ≪ m. It means that we perform the calculation under the assumption that the function f decreases fast enough and the contribution of the region p m to the integral is small. The substitution of the energy E in the form to the Eqs. (15)- (17) gives Here we retain only the terms of the leading and text-to-leading order in the parameter |p|/m. The substitution of the function d(p) in the form (25) to the Eq. (23) gives the Schödinger equation in the momentum space for the particle with reduced mass: Performing the Fourier transform of Eq. (26) we obtain: So, we show that at small interaction constant the BSE is similar to the Schödinger equation for the particle which has the mass equals to the reduced mass. The reason of the another sign before the potential in Eqs. (26), (28) is that the potential V (r) is the electron-electron interaction potential. The condition of weakness of the potential or smallness of the interaction constant means that the lowest binding energy of the electron and hole should be much less than the band gap (|κ| ≪ m).
Let us consider the solution of the Eq. (27) in the case of the Coulomb potential: where e is the electron charge, ǫ is the dielectric constant. For the Coulomb potential the eigenvalues κ of the Eq.
(27) have the form: where n is the radial quantum number, l is the angular momentum quantum number, α = e 2 /ǫ is the integration constant. The applicability condition |κ| ≪ ∆ ∼ m means that the parameter α should be much less than unity (α ≪ 1). However, the experimental results for binding energies demonstrate that the parameter α is of the order of 0.5 for WS 2 (∆ = 2.41eV, κ = 0.32eV), and 0.6 for WSe 2 (∆ = 2.02eV, κ = 0.37eV), see 11 . It means that for such values of the interaction constant it is necessary to consider the exact equation or, at least, the higher order corrections in the parameter α should be taken into account. The exact solution of the BSE for the ground state and for a few levels with l = 0 is calculated numerically in the Section IV. In the present section we find the first correction to the Eq. (27) and to the eigenvalue (30) for l = 0, and n = 0, 1.

A. Singularities in the Hamiltonian
Let us demonstrate that the naive expansion leads to the appearance of the non-integrable operators in the Hamiltonian (28). To find the next-to-leading order correction in α we keep the following terms in the system (15)-(17), and obtain: The substitution of the the function h in the form (32) to Eq. (33), and then the substitution of the result for d to the Eq. (31) gives the following equation: Here we keep only the terms which have the necessary accuracy. One can see that the operator in the right-hand side of this equation is not hermitian. The reason for that is the function f (p) can not be treated as the wave function because the normalization condition (13) for the function f (p) differs from the normalization condition for the wave function, see 13 . So, to obtain the Hamiltonian we should find the relation between the function f (p) and the Schrödinger wave function. The normalization condition (13) in the momentum space in the leading and next-to-leading order in the parameter κ/m has the form: Therefore, performing the change of the function as we obtain that the density ρ is expressed through the function ψ Sh as where we keep only the terms which have the necessary accuracy. The last expression coincides with normalization condition for the Schrödinger wave function. Finally, we substitute the function f in the form (36) to the Eq. (34), then perform the Fourier transformation and obtain the Shrödinger equation in the coordinate space: Here ∇ = ( ∂ ∂x , ∂ ∂y ), [a × b] is the vector product. In the Eq. (39) we also recover the dimension. One can see that the operator in the left-hand side of Eq. (39) is the hermitian operator, the function ψ Sh obeys the correct normalization condition. Therefore, we can treat the function ψ Sh as the wave function and the equation (39) as the Schrödinger equation. Note that the equation (39) coincides with the equation (11) of Ref. 5 , which was obtained from the nonrelativistic expansion of the the Schrödinger equation for two Dirac particles. So, the leading and next-to-leading orders of the non-relativistic expansion for the BSE and Schrödinger equation for two interacting Dirac particles coincides. The difference between the expansions of these equations appears in the higher orders corrections in the interaction constant.
Let us consider the correction δκ which is related with the Hamiltonian (39) to the ground state energy κ 0,0 for Coulomb potential, see Eq. (30): where ψ 0 (r) is the solution of the equation (27) for the case n = 0 and l = 0: where a B = 1/(mα) is the Bohr radius. Since the wave function ψ 0 (r) equals to constant at r = 0, the integral (40) diverges at small distances as: It means that the correction can not be calculated using this expansion. Note that the similar correction which appears in the non-relativistic expansion of the three dimensional Dirac equation is finite 13 . The reason of the appearance of divergence is that in our case electrons propagate in two spatial dimensions whereas electric field propagates in three dimension. Therefore, in our case the interaction potential diverges as 1/r at r → 0, see (29), whereas the Coulomb potential in two spatial dimensions diverges only as log r.
The divergence of the integral (40) indicates that the correction comes from the distances r ≪ a B . So, the assumption that the main contribution to the integrals in Eqs. (15)-(17) comes from the momentum scales p ≪ m is incorrect for the correction under consideration. Therefore, to find the correction of the order of α 4 we should take into account that the interaction potential differs from the Coulomb potential at the distances r ∼ 1/m. So, we should consider the distances p ∼ m in the system (15)-(17) more carefully.
B. The correction of the order of α 4 To find the correction of the order of α 4 for the s-states (l = 0) we substitute the functions f (p), h(p), and d(p) in the form to the Eqs. (15) -(17), and obtain: To find the Eqs.
then we substitute the energy E in the form (22), and obtain following equations for the functions a and b: The operatorṼ obeys the conditionṼ (p, k) =Ṽ (k, p), therefore the operators in the right-hand side of the Eqs. (48), (49) are hermitian one. Also, the potentialṼ decreases faster than potential V (p) at large momentum p. This property allows us to calculate the correction of the order of α 4 to the energy κ 0 .
To calculate the correction we expand the equations (48) and (49). In the necessary order the function b(p) = − p 2 4m 2 a(p).
We substitute this relation to Eq. (48), retain the terms of the order of p 4 /m 2 , and obtain the equation for the function a(p): The Eq. (51) can be solved using standard methods of the perturbation theory. To find the binding energy in the case of the Coulobm potential (29) we present κ and a in the form then substitute these representations in Eq. (51), and obtain: where the functionã So, the ground state energy in the leading and next-to-leading order in the parameter α has the form: The first term in brackets coincides with the energy (30), the second term is the correction. One can see that the α 4 corrections contain the logarithmic term log α which comes from the momentum region a −1 B ≪ p ≪ m. It means that it is necessary to take into account the corrections to the interaction potential related with electron-hole exitations in the intermediate states. The calculation of the correction for the first excited state (2s-state) gives One can see, that it also contain the logarithmic term log α. The corrections of the order of α 4 to the other states can be calculated using described procedure. So, the BSE (44)-(46) allows us to find the correction of the order of α 4 , whereas the equation (39) does not give the possibility to calculate the correction, since it does not take into account the difference of the potential form the Coulomb potential at the distances r ∼ 1/m. The results (58) and (59) are shown in Figs. 1 and 4 by the dasheddotted line. One can see that the corrections are in a good agreement with the exact numetical solution of the BSE for the parameters α 0.5. For the case α > 0.5 it is necessary to take into account the α 6 terms.
Strictly speaking, there are another correction of the order of α 4 which we do not take into account. These corrections are related with the Uehling potential, the cross-box diagram contribution and annihilation diagram contribution. Here we do not consider these corrections and concentrate our attention only on the BSE in the leading order in the interaction potential. The contribution of the mentioned corrections will be consider in the following papers. In the previous section we obtain that the correction of the order of α 4 is not small therefore we solve the BSE exactly for the case of the Coulomb potential and for the states with zero angular momentum l = 0. For this case it is convenient to present the function ψ(p) in the following form: The substitution of the function ψ(r) in the form (60) to the Eq. (8) gives the following system of the equations One can see that the Eq. (63) is trivial, therefore we can consider only Eqs. (61) and (62). Below we solve these equations numerically for two potentials: the Coulomb potential Eq. (29), and the localized potential without singularity at the point r = 0: where U 0 is the depth of the potential well, R is the width of the potential well. In the momentum space the potential has the form: A. Coulomb potential We start our consideration with the Coulomb potential. The Coulomb potential decreases only as 1/q at large q, therefore to increase the accuracy of the numerical calculation we extract the asymptotic behavior of the functions F (p) and H(p) analytically at large momentum q ≫ m. The substitution of the functions F (p) and H(p) in the form gives the following equations for the β and β 1 in the leading order in the parameter α: In the leading order in the parameter α we obtain the solutions of the equations: Therefore, for the numerical calculations we substitute the functions F (p) and H(p) to the Eqs. (61) and (62) in the form: whereF (p) andH(p) tends to constant at p ≫ m. The solution E/2m for different values of the parameer α is presented in Figs. 1 and 3. In Fig. 1 the dependence of the ground state energy in the units 2m on the interaction constant is plotted. One can see that the exact solution of the BSE significantly differs from the result (30) only for α > 0.7. However, for α > 0.8 the difference is sufficient. Moreover, the critical value of α c (E(α c ) = −2m) for BSE differs from the critical value of the interaction constant for the Schödinger equation. The ground state energy E for BSE achieves the valence band (E = −2m) at α = 1 ± 0.05. So, with our calculation precision, we obtain that the critical value of the interacting constant α c = 1. The uncertainty is related with the dependence of the result on the sampling step when calculating the integrals. The results E/(2m) for α < 0.9 is weakly depend on the discretization step of the numerical integration. For α > α c the energy level of the ground state disappears from the discrete spectrum, see Fig. 3. In Fig. 3 we plot the dependence E(α) for the ground state and the first excited state with n = 1, l = 0. One can see that when the α is approaching to unity the slope of the function E(α) goes to −∞.
For α > 1 the state with n = 1, l = 0 becomes the ground state. Such behavior is similar to that of electron in the Coulomb field in the three dimensional electrodynamics.

B. Exponential potential
To investigate the behavior of the energy levels in the vicinity of critical value of the interacting constant in detail we consider the toy potential (64). In Fig. 4 the dependence E/2m on U 0 /m for three lowest binding energies of  BSE is presented. One can see that when the interaction constant is increasing the energy level is decreasing then it cross the line E = 0 and then riches the valence band E = −2m atŨ =Ũ c . For the exponential potential (64) the slope of the function E(Ũ ) does not go to −∞ whenŨ in the vicinity ofŨ c , because of the potential does not have singularities at r → 0. Such behavior resembles the behavior of the electron in the potential well in three dimensional quantum electrodynamics.

V. CONCLUSION
In the present paper we have investigated the electron-hole interaction in the transition-metal dichalcogenides using the BSE obtained in the leading order in the potential. We have obtained the equation for the four component wave function of the exciton. We have investigated the dependence of a few lowest binding energies on the interaction constant. For small interaction constant the equation (51) and analytical expression for the ground state energy have been found up to the terms of the order of α 4 . We have demonstrated that the expansion of BSE in the leading order in the interaction constant coincides with expansion of the Schrödinger equation for the two interacting Dirac particles. However, in the next-to-leading order in the interaction constant the equations are different. The expansion of the Schrödinger equation leads to the appearance of the operators with non-integrable singularities. The appearance of the singularities means that even for small α at small distances the problem is "relativistic". It indicates that the distances r ≪ a B ∼ m/α give the contribution to correction for the energy levels of exciton. The correct expansion of the BSE does not contain the singular operators. Using the equation we have calculated the correction to the binding energies and demonstrated that it contains the term α 2 log α which indeed comes from the distances (1/m) ≪ r ≪ a B .
We have investigated the binding energies for the Coulomb potential and for the exponential potential (64). We have obtained that for the Coulomb potential the value α = 1 is critical. When the parameter α approaches to unity the energy of the ground state goes to zero. If the value of the parameter α is slightly greater than 1 this state disappears from the discrete spectrum. This picture is close to that in quantum electrodynamics, see 13 . The reason of such behavior is the singularity of the Coulomb potential. To confirm this statement we found the spectrum for the potential without singularities. We showed that for the increasing interaction constant the bound state energies decreases from E = 2m smoothly to E = 0, and then to E = −2m. The energy level disappears when it reaches the valence band (E = −2m).
We have obtained that when the interaction constant α is close to 0.5-0.6 the binding energies slightly differ from that obtained in the leading order, see Eq. (30). So, the BSE obtained in the leading order does not describe the experimental results. However, the investigation of the BSE shows that the structure of the potential at small distances should affect energy levels significantly. Therefore, to describe the experimental data it is necessary to include the Uehling type diagrams, and cross-box type of the diagrams in the BSE. These types of the corrections will be considered in the following papers.