Magnetoexcitons in transition-metal dichalcogenides monolayers, bilayers, and van der Waals heterostructures

We study direct and indirect magnetoexcitons in Rydberg states in monolayers and heterostructures of transition-metal dichalcogenices (TMDCs) in an external magnetic field, applied perpendicular to the monolayer or heterostructures. We calculate binding energies of magnetoexcitons for the Rydberg states 1$s$, 2$s$, 3$s$, and 4$s$ by numerical integration of the Schr\"{o}dinger equation using the Rytova-Keldysh potential for direct magnetoexcitons and both the Rytova-Keldysh and Coulomb potentials for indirect magnetoexcitons. Latter allows understanding the role of screening in TMDCs heterostructures. We report the magnetic field energy contribution to the binding energies and diamagnetic coefficients (DMCs) for direct and indirect magnetoexcitons. The tunability of the energy contribution of direct and indirect magnetoexcitons by the magnetic field is demonstrated. It is shown that binding energies and DMCs of indirect magnetoexcitons can be manipulated by the number of hBN layers. Therefore, our study raises the possibility of controlling the binding energies of direct and indirect magnetostrictions in TMDC monolayers, bilayers and heterostructures using magnetic field and opens an additional degree of freedom to tailor the binding energies and DMCs for heterostructures by varying the number of hBN sheets between TMDC layers. The calculations of the binding energies and DMCs of indirect magnetoexcitons in TMDC heterostructures are novel and can be compared with the experimental results when they will be available.

TMDCs are semiconductors that have a chemical formula MX 2 where M is a transition metal atom (Mo or W) and X is a chalcogen atom (S, Se, or Te). The comprehensive review of excitonic complexes in TMDCs monolayers is given in Refs. [14,25,26]. Excitons in monolayers have a direct energy gap. However, TMDCs heterostructures made of monolayers have indirect gaps. There are two possible arrangements of double-layer heterostructures: bilayer TMDC-TMDC [19,20,[27][28][29][30][31][32][33][34][35][36][37] and TMDC-insulator-TMDC [17,35,38]. For both arrangements, there are two types of heterostructures: type-I and type-II. Type-I heterostructure has a direct band gap alignment, and type-II has a staggered band alignment [19,20,28,30,31]. However, stacking order [39,40] or application of strain [41] can be used to induce an indirect gap into a direct gap. In addition, the application of the perpendicular electric field controlled by voltage modifies the band structure, so it is advantageous for electrons and holes to reside in different TMDC layers [17,19,42] and makes indirect excitons more preferable than direct excitons. However, the application of the electric field can lead to dissociation of indirect excitons [43]. Nevertheless, the extended lifetime of indirect excitons [28,29,34,44] compared to direct excitons presents great interest for the design of electronic devices.
Historically, the electromagnetic interaction between the electron and hole was described by Coulomb potential (V C ) for both direct and indirect excitons [45,[47][48][49][81][82][83][84][85][86]. However, as has been shown, the Rytova-Keldysh (RK) potential (V RK ) [87,88] is the appropriate potential to describe the interaction of electron and hole in two-dimensional (2D) configuration space since the RK potential takes into account the effects of the dielectric environment and the 2D confinement. As a result, currently, both the RK [26, 51, 52, 54-56, 62, 65, 69, 89-92] and Coulomb potentials are used to describe the direct exciton. For indirect excitons, the Coulomb potential is used. However, in a few cases, the RK potential is used as well [43,93]. We do the same thing: we use both the RK and Coulomb potentials for indirect excitons to investigate the importance of using the accurate potential that describes an interaction between the hole and electron.
Magnetoexcitons in TMDCs monolayers have brought the considerable recent interest with respect to magnetic field tuning than other 2D materials, since they preserve time-reversal symmetry with excitons formed at K and K points at the boundary of the Brillouin zone, which restricts valley polarization. On one side, the contribution of the magnetic field to the Rydberg states is small, however the tunability of the binding energy of magnetoexcitons brings unique potential for controlling novel device applications in optoelectronics. On the other side, as we demonstrate in this work, the realization of van der Waals heterostructures consisting of stacked 2D layers, where indirect magnetoexcitons formed, can be significantly controlled by engineering numbers of hBN monolayer, in addition to be tuned by the external magnetic field.
In this paper, we study the effects of the external magnetic field on the binding energies of 1s, 2s, 3s, and 4s Rydberg states of direct and indirect magnetoexcitons in TMDCs. A and B direct magnetoexcitons are considered in freestanding (FS) and encapsulated by the hexagonal boron nitride (hBN) MX 2 monolayers, and A and B indirect magnetoexcitons are considered in bilayer MX 2 -MX 2 and MX 2 -hBN-MX 2 . In our study to find eigenfunctions and eigenenergies for indirect magnetoexcitons, a two-particle Schrödinger equation is solved numerically using the Rytova-Keldysh and Coulomb potentials to calculate binding energies and diamagnetic coefficients. We demonstrate strong sensitivity of the energy contributions from the external magnetic field on the type of the potential. The first time diamagnetic coefficients (DMCs) for bilayer MX 2 -MX 2 and MX 2 -hBN-MX 2 heterostructures are calculated and reported. Our study raises the possibility to control the binding energies of direct and indirect magnetoexcitons in monolayer MX 2 , bilayer MX 2 -MX 2 , and MX 2 -hBN-MX 2 heterostructures using a magnetic field and opens an additional degree of freedom to tailor the binding energies and DMCs for heterostructures by varying the number of hBN sheets between TMDC layers.
The paper is organized in the following way. In Sec. II we introduce the theoretical formalism for the description of Mott-Wannier magnetoexcitons in freestanding and encapsulated MX 2 monolayers, and in TMDC heterostructures and discuss electrostatic interactions that form direct and indirect magnetoexcitons. In Sec. III we focus on the analysis and discussion of the results. In particular, contributions from the external magnetic field to binding energies of magnetoexcitons in MX 2 monolayers, MX 2 -MX 2 bilayers, and MX 2 -hBN-MX 2 heterostructures are presented and discussed in Secs. III A and III B, respectively. Sec. III C is devoted to the calculations of diamagnetic coefficients for direct and indirect magnetoexcitons. Finally, conclusion follows in Sec. IV.

II. THEORETICAL FORMALISM
In this section, we provide the theoretical formalism for describing the Mott-Wannier magnetoexciton in 2D materials and present the energy contribution from the external magnetic field to the Rydberg states binding energy of magnetoexcitons. Electrostatic interactions that due to a non-local dielectric screening of an electron-hole interaction strongly modify the electrostatic Coulomb potential and form direct magnetoexcitons in TMDC monolayers and indirect magnitoexcitons in MX 2 -MX 2 bilayers and MX 2 -hBN-MX 2 van der Waals heterostructure are presented.

A. Mott-Wannier magnetoexciton in 2D materials
Our approach is an effective mass model of excitons and we are using this model for the approximation of the solid state Hamiltonian for an interacting electron and hole. Therefore, our theory is based on two-body excitonic Hamiltonians in the effective mass approximation with screened interactions appropriate for TMDC materials. So, let's introduce the equation for the description of Mott-Wannier excitons in the external magnetic field. We employ the effective mass model for two charged point-like particles in two dimensions. In general, electrostatically-bound electrons and holes in the external magnetic field form magnetoexcitons. To find the eigenfunctions and eigenenergies of a 2D magnetoexciton in TMDCs monolayer, bilayer MX 2 -MX 2 , and MX 2 -hBN-MX 2 heterostructure in the external magnetic field, we write the Schrödinger equation for an interacting electron and hole [83] where e is the charge of the electron, the indices e and h are referring to the electron and hole, respectively, r e and r h are 2D coordinates, m e and m h are the masses of charged particles, A(r e/h ) = B × r e/h /2 is a gauge vector potential, and V (r e , r h ) is the potential of interaction between the electron and hole confined in 2D space. In threedimensional (3D) homogeneous dielectric environments, the electron-hole interaction is described by the Coulomb potential. However, in 2D monolayer this interaction has to be modified because of the reduced dimensionality and screening effects. The original derivation of two charged particles interaction in 2D space was given by Rytova [87] and a decade later was independently obtained by Keldysh [88]. The potential is called the Rytova-Keldysh (RK) potential. For almost a decade, the RK potential has been used to describe electrostatic interaction between charge carriers of few-body complexes in TMDCs, phosphorene, and Xenes monolayers (See [26] and references herein).
Following Refs. [45, 47-49, 83, 86], in Eq. (1) we introduce the coordinate of the center-of-mass R = mere+m h r h M , where M = m e + m h is the total mass of the system, and the relative motion coordinate r = r e − r h and consider the magnetic field pointing in z-direction that is perpendicular to the monolayer where the exciton is located, B = B(0, 0, B). After performing the standard procedure for the coordinate transformation to the center-of-mass, Eq. (1) becomes: where γ = m h −me m h +me and µ = mem h me+m h is the reduced mass. Equation (2) is written for the case when masses of the electron and hole are different: m h = m e , which is the case in TMDCs. The term (B × r) · ∂ ∂r = B · L = 0 since we consider Rydberg states: 1s, 2s, 3s, and 4s, for which l = 0, m l = 0.
Following Refs. [48,83,86], we introduce an operatorP : It is easy to check thatP commutes with the Hamiltonian in Eq. (1), therefore, it has the same eigenfunction as Eq.
(1). Thus, one can write the wave function for the exciton in the magnetic field as [48,86]: where the notationρ 0 = 1 eB 2 (B × P ) is introduced. After substituting the wave function ψ(R, r) in Eq. (2), the Schrödinger equation for the relative motion of the electron and hole reads: Finally, after separating the angular variable in (5), the equation with zero center-of-mass momentum reads [84,86]: Equation (6) describes the Mott-Wannier magnetoexciton in Rydberg optical states in 2D materials. This equation has a long history of the solution in the case of the electron-hole Coulomb interaction [45,[47][48][49][81][82][83][84][85][86]. However, we solve Eq. (6) using V RK for direct magnetoexcitons. For indirect magnetoexcitons, both V RK and V C potentials are used to investigate the importance of using correct potential describing electron-hole interactions. Note that Eq. (6) does not explicitly contain any spin-or valley-dependent Zeeman terms. To find binding energies, we numerically solve Eq. (6) by using the code implemented in Ref. [75] which was successfully modified and used to calculate binding energies of magnetoexcitons in TMDCs monolayers [71] and Xenes heterostructures [94]. The method is based on using the finite element method implemented in Wolfram Mathematica in the NDEigensystem function. The code was modified in a way that the Schrödinger equation explicitly contains e 2 8µ (B × r) 2 term. To check the code, we use the input parameters from respective papers listed below and calculate the binding energies of direct and indirect excitons. The code reproduces theoretical binding energies of excitons in TMDCs monolayers reported in Refs. [50,69,95], obtained in the framework of the stochastic variational, the path integral Monte-Carlo and variational methods, respectively, within 5%. The experimental binding energies reported in Refs. [53][54][55] are reproduced within 7%. Related to the bilayer system composed from two different TMDC monolayers, we reproduce theoretical binding energies, obtained using a numerical method based on exterior complex scaling given in Ref. [43], within 3% when the parameter κ is varied between 1 and 5. The theoretical binding energies reported in Ref. [93], obtained using quantum Monte Carlo method, are reproduced within 10%. It is worth noting that in the above papers for the bilayer system, the RK potential is used to described interactions between the electron and hole.
The energy contribution from the magnetic field and DMCs have been considered so far only for direct excitons in TMDCs monolayers where the RK potential is used. In Ref. [54] the external magnetic field is treated as a small perturbation. In Ref. [53] Eq. (6) is solved on a grid, i.e. an unknown function is introduced in order to split the second-order differential equation into two first-order equations. In Ref. [55] Eq. (6) is solved numerically. In Ref.
[69] the stochastic variational method with a correlated Gaussian basis is used to calculated binding energies. They provide values of binding energies of 1s state at different values of the magnetic field. We reproduce their values when we use their parameters in our code.
By calculating the magnetoexciton energy of Rydberg states at different values of the magnetic field, we can find the energy contribution from the magnetic field to the binding energy of a magnetoexciton in the following way: In Eq. (7) E 0 is the exciton binding energy when the magnetic field is absent and is calculated with respect to the two-body threshold [69,[95][96][97] while E(B) is the magnetoexciton energy at some value of the magnetic field. In literature, the Zeeman and diamagnetic shifts of excitons in TMDCs monolayers are treated in the same way as ZM and DM shifts of excitons in quantum wells and quantum dots [97][98][99][100][101][102][103][104][105]. The first step is to write a magnetoexciton energy using Taylor series [97]: In Eq. (8) the second and third terms are the ZM and DM shifts, respectively. When the magnetic field's energy contribution is small compared to the binding energy in the absence of the magnetic field, we can use the first three terms of Taylor series to describe the magnetoexciton binding energy [65,97,101,103]. In other words, the following As stated before, the treatment of the ZM and DM shifts of magnetoexcitons in TMDCs monolayers follows the same procedure applied for quantum dots and quantum wells. The second step in describing ZM and DM shifts is to define them by connecting shifts to terms of Taylor series. The Zeeman shift is defined as [57,60,62,64,65,101,103,104], where g is the effective g factor and µ B is Bohr magneton. Experimentally, the ZM shift is defined as energy difference between [51, 53-55, 61, 62, 64, 65, 102, 103]. The diamagnetic shift is defined as 53-55, 65, 97-105], where r 2 is the expectation value of r 2 over the exciton envelope wave function. In our work, we follow notation used for magnetoexciton diamagnetic coefficient in TMDCs monolayers i.e. γ 2 ≡ σ [51, 53-55, 60, 64]. According to Refs. [51, 53-55, 64, 97-105] the DMC for a carrier in a semiconductor is defined as σ = e 2 8µ r 2 . Therefore, by calculating DMCs, exciton radius, reduced mass, and dielectric properties of a material can be obtained [51,[53][54][55]97]. Experimentally, the DM shift is defined as σB 2 = E(K)+E(K ) 2 [54,57,65,101,103,105]. In Fig. 1 we show the schematic formation of bright (optically allowed) direct intravalley A and B excitons in K and K valleys under the right-and left-polarized light. We call the excitons that are formed by charge carriers with parallel spins at the conduction and valence bands separated by a small gap as direct A excitons. The excitons formed by charge carriers with parallel spins at the conduction and valence bands separated by a large gap are called direct B excitons [50,51,106]. There exist two additional exciton types, when charge carriers have antiparallel spins, but these excitons are optically forbidden and called spin-forbidden dark excitons [107,108]. In the bilayer and MX 2 -hBN-MX 2 structures, the conduction band minimum and valence band maximum reside in two different layers forming the indirect interlayer exciton. There are two possible stacking orders in TMDCs: AA and AB [109][110][111][112]. In the case of AA stacking interlayer intravalley exciton is formed by the hole from the valence band in K valley of the layer 1 and the electron from the conduction band in K valley of the layer 2 [112]. Therefore, µ of direct and indirect excitons are the same. In the case of AB stacking, the lower layer is 180 • in plane rotation of the upper layer [113]. So, band spins of K valley of the layer 1 are flipped compared to corresponding band spins of K valley of the layer 2. Therefore, the bright interlayer intravalley exciton is formed by the hole with m B in K valley of the layer 1 and the electron with m A in K valley of the layer 2 [114][115][116][117].

B. Magnetoexcitons in TMDCs monolayer
A non-local dielectric screening of an electron-hole interaction strongly modifies the electrostatic Coulomb potential. It leads to a non-hydrogenic Rydberg series of exciting magnetoexciton states. The Rytova-Keldysh potential has become a prevalent description of the electrostatic interaction of charge carriers in 2D systems. The RK potential is central potential, and the interaction between the electron and hole for direct excitons in an encapsulated TMDC monolayer has the form [87,88]: where r ≡ r eh = r e − r h is the relative coordinate between the electron and hole. In Eq. (9), κ = ( 1 + 2 )/2 describes the surrounding dielectric environment, 1 and 2 are the dielectric constants below and above the monolayer, H 0 and Y 0 are the Struve and Bessel functions of the second kind, respectively, and ρ 0 = 2πχ 2D /κ [95] is the screening length, where χ 2D is the 2D polarizability of the monolayer, which is given by χ 2D = h /4π [88], where is the bulk dielectric constant of the electron/hole containing monolayer and h is the thickness of a TMDC monolayer. At r < ρ 0 the RK potential has logarithmical dependence and diverges logarithmically at the origin. At long-range distances, r > ρ 0 , it retains the Coulomb potential behavior. Using the RK potential for description of excitons in 2D monolayer, one finds a strong dependence of the binding energy on whether the monolayer is suspended in air (Fig. 2a), encapsulated by hBN (Fig. 2b), or in heterostructures (Figs. 2c and 2d). It is worth mentioning that the RK potential is obtained for a monolayer. However, a TMDC monolayer consists of three atomic layers: chalcogen-metal-chalcogen sheets. To address this problem, a new potential is derived in Ref. [118] which considers the three atomic sheets that compose a monolayer of TMDC and can explain the non-hydrogenic Rydberg series of excitons in TMDC monolayer. The Rytova-Keldysh potential can be recovered when considering the strict 2D limit, and in our study, we are using the RK potential.

C. Indirect magnetoexcitons in TMDCs heterostructures
Van der Waals attractive forces can hold together layered TMDC materials. Let us now consider the formation of indirect magnetoexcitons in the heterostructure molded by the same type of TMDC monolayers: a bilayer MX 2 -MX 2 heterostructure, shown in Fig. 2c and van der Waals MX 2 -hBN-MX 2 heterostructure, where two TMDC monolayers are separated by N layers of hBN monolayers as the ones used in most of the experiments, presented in Fig. 2d. In such structures, the electron and hole move in the planes of different layers, form an indirect exciton, and have restricted motion between the layers due to the dielectric barrier. The latter is a two-body problem in the restricted 3D space, and the electrostatic electron-hole interaction V (r) could form a bound state, i.e., the indirect magnetoexciton. Therefore, to determine the binding energy of the magnetoexciton one must solve a two-body problem in the restricted 3D space due to the restriction of the motion in z-direction. The relative separation r between the electron and hole can be written in cylindrical coordinates as r = ρρ + Dẑ, whereρ andẑ are unit vectors. Writing r in cylindrical coordinates allows us to treat the case of direct excitons in a TMDC monolayer and spatially indirect excitons in bilayer MX 2 -MX 2 and MX 2 -hBN-MX 2 heterostructures on equal footing. If we set D = 0 in the latter expression, it becomes a purely 2D equation, with ρ representing the separation between the electron and hole sharing the same plane. For spatially indirect excitons, the relative distance between the electron and hole is r = ρ 2 + D 2 , where D is the distance between the middle of two TMDC layers assuming that the electron and hole reside in the middle of their respective sheets [43].
Thus, for the description of indirect Mott-Wanner magnetoexcitons one can use Eq. (6) with interactions for the RK potential, and for the Coulomb potential. In Eq. (10) ρ 0 = ρ 0 and ρ (2) 0 are the screening lengths of the first and second monolayer, respectively. Polarizability of both layers needs to be taken into account since the total screening length has increased. If polarizabilities of both layers are not taken into account then binding energies are higher [43]. Our calculations show that if polariziblity for only one layer is considered then the binding energies of 1s state increase by about 25-30%. Equations (10) and (11) describe the interaction between the electron located in one and hole in the other parallel TMDC monolayers. Therefore, one can obtain the eigenfunctions and eigenenergies of magnetoexcitons by solving Eq. (6) using the potential (9) for direct magnetoexcitons, or using either potential (10) or (11) for indirect magnetoexcitons. For indirect magnetoexcitons, we perform calculations using both the RK and Coulomb potentials. This allows a better understanding of the importance of the screening effect in MX 2 -MX 2 and MX 2 -hBN-MX 2 heterostructures.
It is worth mentioning that the RK potential was originally formulated as an explicitly 2D description of the Coulomb interaction. Nevertheless, there have been recent attempts to apply the RK potential to indirect excitons in van der Waals heterostructures of 2D materials such as the TMDCs, phosphorene, and Xenes [17,18,23,75,94,119,120]. The logic behind considering the RK potential for indirect excitons follows from two considerations: i. the dielectric environment is still inhomogeneous, just as in the case of the direct exciton -when the interlayer separation D is smaller than, or comparable to, the RK potential screening length ρ 0 and the excitonic gyration radius r 2 , the electron-hole interaction potential must account for both the TMDC monolayers and the interlayer dielectric, and ii. as the interlayer separation D becomes larger than ρ 0 , the total separation between the electron and hole, r = ρ 2 + D 2 , becomes much larger than ρ 0 , and, therefore, the RK potential converges towards the Coulomb potential. Let us emphasize that we are not claiming definitively that the RK potential provides the most accurate description of the spatially indirect exciton, hence, we present the extensive comparison with the corresponding results obtained using the Coulomb potential.
Above, we consider ideal cases when encapsulated TMDC monolayer, bilayer MX 2 -MX 2 , and MX 2 -hBN-MX 2 heterostructures are fabricated so that each layered material is stacked without any vacuum or air interlayer gap. In reality, between the layers always exists non-vanishing interlayer gap. In each layer the field lines of the Coulomb interaction are screened by the adjacent material, which reduces the single-particle band gap as well as exciton binding energies [121]. In Ref. [121] the authors demonstrate and give a quantitative understanding that the binding energy of excitons and electronic and optical properties are sensitive to the interlayer distances on the atomic scale.
Throughout this paper, we consider the separation between the electron and hole residing in two TMDC monolayers of the MX 2 -hBN-MX 2 van der Waals heterostructures in steps of calibrated thickness, l hBN = 0.333 nm, corresponding to the thickness of one hBN monolayer: D = h + N l hBN , where h is the TMDC monolayer thickness and N is the number of hBN monolayers. For a bilayer MX 2 -MX 2 heterostructure N = 0 and following Ref. [43] the electron and hole reside in the middle of their respective sheets and r = ρ 2 + h 2 .
The diamagnetic coefficients for the bilayer MX 2 -MX 2 and MX 2 -hBN-MX 2 heterostructures are reported for the first time. We are considering the further knob to tailor the binding energies and diamagnetic coefficients for the MX 2 -hBN-MX 2 heterostructures by varying the number of hBN sheets between TMDC layers and present corresponding calculations.
In our calculations, we vary the magnetic field in the increment of 1 T in the range from 0 T to 30 T and use input parameters given in Table A.1 presented in Appendix A. In Table A.1, the values of A high and A low correspond to the parameters found in the literature, which maximize and minimize the A exciton binding energy, respectively. It is highly likely that each parameter's true value for a given material falls somewhere within the given range. Therefore, the true magnitude of the calculated quantities for A excitons studied in this paper lies somewhere between the calculated values. Values of B correspond to parameters given in the literature for B excitons.

A. Contribution from the external magnetic field to binding energies of magnetoexcitons in a monolayer
The results of the energy contribution from the magnetic field to the binding energy of Rydberg states of direct A low magnetoexcitons in freestanding and encapsulated TMDC monolayers are reported in Fig. 3. The comparative analysis of the results presented in Fig. 3 shows the following: i. the energy contribution from the magnetic field to the binding energy of the direct magnetoexcitons in FS monolayers of TMDC materials is always less than in encapsulated monolayers, and the difference increases with the increase of the magnetic field; ii. the direct magnetoexcitons in FS monolayers are bound in 1s, 2s, 3s, and 4s states within all considered range of the magnetic field, while magnetoexcitons in encapsulated monolayers dissociated in the states 3s and 4s when the magnetic field is greater than some particular values. These conclusions are related to the direct A low magnetoexcitons. Interestingly enough that these conclusions stay the same for direct A high and B magnetoexcitons. However, the ∆E is systematically smaller for the A high , while ∆E for B magnetoexcitons lie between ∆E for the A low and A high magnetoexcitons. Analysis of the results shows that the energy contribution from the external magnetic field to the magnetoexciton binding energy depends on material parameters. In fact, the binding energy of magnetoexcitons is a function of the reduced mass µ and polarizability χ 2D : the binding energy of a magnetoexciton is bigger for the larger µ and smaller χ 2D , while, when µ is smaller and χ 2D is larger, ∆E is bigger. This pattern qualitatively coincides with previously reported calculations for magnetoexcitons in WSe 2 and MoSe 2 monolayers [71].
Results of our calculations show three distinct features: i. the energy contribution from the magnetic field to magnetoexcitons in the FS and encapsulated monolayers for the state 1s is more than one order of magnitude smaller than for the states 2s, 3s and 4s; ii. the magnetoexcitons in FS monolayers are bound in states 1s, 2s, 3s and 4s when the magnetic field varies up to 30 T; iii. magnetoexcitons in encapsulated monolayers are dissociated in states 3s and 4s at some values of the magnetic field while staying bound in states 1s and 2s at higher magnetic field values.
Let us address the a dissociation of the magnetoexcitons in encapsulated monolayers in states 3s and 4s at some values of the magnetic field. To understand the dissociation behavior of magnetoexcitons in 3s and 4s states in TMDC monolayer, following Refs. [122][123][124] we examine the characteristic behavior of the wavefunction and the total interaction potential of the electron-hole system in the external magnetic field. In Fig. 4 is shown the behavior of wavefunctions for different states and the corresponding total potential V RK + e 2 8µ B 2 r 2 . The total potential of the system, V RK (r) + e 2 8µ B 2 r 2 , and the wavefunctions, Φ 1s (r), Φ 3s (r) and Φ 4s (r) are plotted as a function of r. The total potential is given at two different values of the magnetic field when the dissociation of the 3s and 4s states occurs. As can be seen from Fig. 4, the total potential becomes positive as the magnetic field increases. The the wavefunction for the bound 1s state is localized, while the wavefunctions for the 3s and 4s states are delocalized at B = 16 T and B = 7 T, respectively.  4. Dependencies of the 1s, 3s, 4s states wavefunctions and total potential energy on the electron-hole distance r. The graphs show the behavior of the total potential and wave functions of states 1s, 3s, and 4s of A low magnetoexciton in WSe2 monolayer encapsulated by hBN at different values of the magnetic field. The magnetoexciton dissociates in states 3s and 4s at 16 T and 7 T, respectively. The graphs demonstrate the behavior of Φ3s(r) and Φ4s(r) that dissociate at 16 T and 7 T, respectively, with respect to the bound Φ1s(r) and the corresponding VRK (r) + e 2 8µ B 2 r 2 potential. B. Contribution from the external magnetic field to binding energies of magnetoexcitons in bilayer MX2-MX2 and MX2-hBN-MX2 heterostructures First, we consider bilayer system. We compare the contributions to binding energies of indirect magnetoexcitons due to the external magnetic field calculated using the RK and Coulomb potentials. In Fig. 5, as a representative case, are shown the energy contributions from the magnetic field to the binding energies of indirect magnetoexcitons in bilayer MoS 2 -MoS 2 for Rydberg states 1s (5a), 2s (5b), 3s (5c), and 4s (5d). The binding energies are calculated using V RK and V C potentials. Data are plotted for A high , A low , and B magnetoexcitons. The comparison shows that i. ∆E RK > ∆E C and the difference increases with the magnetic field increase; ii. results obtained with both potentials satisfy the following inequality: ∆E A high < ∆E B < ∆E A low . Most importantly, it is worth to mention that the contribution from the magnetic field is systematically bigger for the bilayer MoS 2 -MoS 2 than for the FS MoS 2 monolayer for all states when calculations are performed with V RK . In other words, ∆E A low for indirect magnetoexcitons in MoS 2 -MoS 2 bilayer is greater than ∆E A low for direct magnetoexcitons in FS MoS 2 monolayer. Similar patterns are observed for other bilayers constituent from MoSe 2 , WS 2 and WSe 2 monolayers, respectively. It is interesting to note that for the bilayer MX 2 -MX 2 system the energy contribution from the magnetic field for the state 1s is more than one order of magnitude smaller than for the states 2s, 3s and 4s and the states do not dissociate in the range of the varying magnetic field.
Next, we consider van der Waals heterostructure. Since the binding energies for excitons in TMDCs monolayers are reported in numerous publications, for example in Refs. [50, 51, 53-55, 69, 95]. As the first step, we calculate the binding energies of indirect excitons in MX 2 -hBN-MX 2 heterostructures using RK and Coulomb potentials in the absence of the magnetic field. The results of the calculations for 1s Rydberg state are presented in Table A.2, Appendix B. At the next step, we consider indirect magnetoexcitons in MoS 2 -hBN-MoS 2 heterostructure. In Fig. 6 are presented the results of calculations for the dependence of the ratio ∆E RK /∆E C on the magnetic field and number of the hBN layers separating two parallel TMDC layers. Calculations are performed for A high magnetoexcitons. The ratio weakly depends on the magnetic field and converges towards 1 when the number of hBN layers increases. The RK potential always gives the larger contribution than the Coulomb potential. Such a result is understandable because the potential V RK ( ρ 2 + D 2 ) converges to V C ρ 2 + D 2 when D increases. From the known asymptotic properties of the Struve and Bessel functions [125,126], it is easy to show that lim in Fig. 7. For 1s and 2s states the energy contribution increases with the increase of the external magnetic field and number of the hBN layers. Indirect magnetoexcitons dissociate in state 3s when B > 11 − 15 T, and in 4s state, when B > 5 T. In Fig. 7a, it is shown the energy contribution for the B exciton always falls between ∆E A high and ∆E A low : ∆E A high < ∆E B < ∆E A low . The same holds for the states 2s, 3s and 4s, but not shown in Figs. 7b, 7c, and 7d. It is worth mentioning that Fig. 7 shows a representative case for indirect magnetoexcitons in MX 2 -hBN-MX 2 heterostructure. The same kind patterns and quantitative dependencies we obtained for the WSe 2 -hBN-WSe 2 heterostructures shown in Appendix C, Fig. 8. Data for WS 2 -hBN-WS 2 and MoSe 2 -hBN-MoSe 2 are not shown. The differences occur only in the energy contribution magnitudes, which varies within up to 40% for MoSe 2 -hBN-MoSe 2 for the maximum value, and indirect magnetoexcitons do not dissociate in the 3s state and dissociate at larger values of the magnetic field in 4s state. The tangerine-based heterostructures have the same kind of patterns as are observed for the molybdenum-based heterostructures. In the same Appendix C, we present in Fig. 9, as a representative case, the comparison of the energy contribution to the binding energy of indirect magnetoexcitons for MoSe 2 -hBN-MoSe 2 heterostructure obtained by solving the Schrödinger equation with the V RK and V C potentials. The analysis shows that the V RK potential gives higher contributions to the binding energy than the Coulomb potential for all states, and these contributions increase with the increase of the magnetic field, and as the number of hBN layers increase contributions converge as is also shown in Fig. 6. The same qualitative patterns are observed in other heterostructures.

C. Diamagnetic coefficients
The diamagnetic coefficients σ A high/low for A and σ B for B direct magnetoexcitons in encapsulated MX 2 monolayers obtained in the framework of our approach are reported in Ref. [71]. Here, we report DMCs for direct magnetoexcitons in FS MX 2 monolayers and indirect magnetoexcitons in MX 2 -MX 2 bilayers and MX 2 -hBN-MX 2 heterostructures in Tables I -III and Table A.3, Appendix D.
The DMCs of direct magnetoexcitons in FS monolayers along with the DMCs for magnetoexcitons in encapsulated WSe 2 , WS 2 , MoSe 2 , and MoS 2 are presented in Table I. Calculations are performed using the V RK potential. The DMCs σ A high < σ B < σ A low for all states. For the 1s state σ hBN > σ FS and the ratio σ hBN /σ FS varies from 1.6 to 2.6, depending on monolayer parameters. For the 2s state this ratio significantly exceeds 2.6. Therefore, the screening effect in encapsulated monolayers substantially increases the value of DMCs. In addition, as can be seen from Table  I in states 3s and 4s E 0 ∼ |E(B) − E 0 | and Eq. (8) cannot be applied to extract σ.
Results of calculations of DMCs for the molybdenum-and the tungsten-based van der Waals heterostructures are presented in Table III and Table A.3, Appendix D, respectively. We present results for σ A high and σ A low for indirect A magnetoexcitons obtained with parameters for each material that give the highest and lowest binding energies, σ A high/low for A excitons and σ B for B magnetoexcitons, respectively. For each type of indirect magnetoexciton, we report two sets of σ: from solution of the Schrödinger equation with V RK and V C potentials, respectively. The DMCs for MoS 2 -hBN-MoS 2 are always a bit higher than for MoSe 2 -hBN-MoSe 2 heterostructure. The latter is related to the small difference of the gaps between valence and conduction bands for MoS 2 and MoSe 2 monolayers. As follows from Table A.3, Appendix D, the same conclusion can be extended to the DMCs for WS 2 -hBN-WS 2 and WSe 2 -hBN-WSe 2 heterostructures. However, the DMCs for the tungsten-based heterostructures are significantly larger than for the molybdenum-based heterostructures because the larger difference of the gaps between valence and conduction bands in the tungsten-and molybdenum-based monolayers [14].
The data analysis of Tables III and A.3 shows the following important features for both S 2 -and Se 2 -based heterostructures: i. the DMCs of the A and B magnetoexcitons increase when the number of hBN layers increase; ii. the DMCs obtained using Rytova-Keldysh potential are always bigger than one obtained for the Coulomb potential; iii. for the states 2s, 3s and 4s E 0 ∼ |E(B)−E 0 | and Eq. (8) cannot be applied to extract σ A high/low and σ B from data; iv. for both the RK and Coulomb potentials σ A high < σ B < σ A low . The latter indicates that DMCs are sensitive to the exciton reduced mass and polarizability of material: the smaller reduced mass of exciton and larger polarizability lead to the bigger DMC value. A distinct feature for the diamagnetic coefficients in bilayer and van der Waals molybdenum-based heterostructures is that σ is always smaller than for the tungsten-based heterostructures due to the difference of the valence and conduction bands in MoX 2 and WX 2 monolayers. The other distinct feature is that the diamagnetic coefficients can be extracted for the 1s, 2s, 3s, and 4s states for the bilayers, while the DMCs exist in the van der Waals heterostructure with up to six hBN layers for 1s and three hBN layers, for 1s and 2s states, respectively. The linear dependence of ∆E on B 2 is invalidated in 3s and 4s states for both the Rytova-Keldysh and Coulomb potentials.

IV. CONCLUSION
In the present paper, we have studied energy contributions from the external magnetic field to the binding energies of A and B magnetoexcitons in 1s, 2s, 3s, and 4s Rydberg states in WSe 2 , WS 2 , MoSe 2 , and MoS 2 freestanding and encapsulated monolayers, bilayer MX 2 -MX 2 , and MX 2 -hBN-MX 2 heterostructures. The first time diamagnetic coefficients for the bilayer MX 2 -MX 2 and MX 2 -hBN-MX 2 heterostructures are calculated and reported. We consider the additional degree of freedom to tailor the binding energies and diamagnetic coefficients for the MX 2 -hBN-MX 2 heterostructures by varying the number of hBN sheets between TMDC layers. The study is performed for direct and indirect A and B excitons. For the A magnetoexcitons we used two sets of parameters given in literature that provide the highest (A high exciton) and lowest (A low exciton) binding energy for excitons. For indirect magnetoexcitons in bilayer and van der Waals heterostructures the dependence of energy contribution from the external magnetic field to the binding energies and DMCs are studied using Rytova-Keldysh and Coulomb potentials.
Our study shows that the energy contributions ∆E for direct and indirect magnetoexcitons and DMCs are sensitive to the reduced mass of exciton and the screening distance ρ 0 : the smaller reduced mass and larger ρ 0 lead to the bigger ∆E and DMCs.
The energy contributions ∆E for direct magnetoexcitons increase with the increase of the external magnetic field, and ∆E is more significant for 2s, 3s and 4s states. The DMCs for direct magnetoexcitons in encapsulated monolayers are bigger than in FS monolayers, σ hBN > σ FS , and the ratio σ hBN /σ FS varies from 1.6 to 3.6 in states 1s and 2s depending on monolayer parameters and the state of magnetoexcitons. Therefore, the screening effect in encapsulated monolayers significantly increases the value of DMCs.
For our investigation, a two-particle Schrödinger equation for indirect excitons is solved using the Rytova-Keldysh and Coulomb potentials to calculate binding energies and diamagnetic coefficients. We demonstrate strong sensitivity of the energy contributions from the external magnetic field on the type of the potential: ∆E RK > ∆E C and the difference increases with the increase of the magnetic field. For all structures the energy contribution have the following order: ∆E A high < ∆E B < ∆E A low . The DMCs for MX 2 -hBN-MX 2 heterostructures obtained using Rytova-Keldysh potential are always bigger than one obtained for the Coulomb potential and the DMCs of the A and B magnetoexcitons increase with the increase of the number of hBN layers. In the considered range of the magnetic field, the direct and indirect magnetoexcitons have stable 1s and 2s Rydberg states. The 3s and 4s states magnetoexcitons dissociate at some values of the magnetic field depending on parameters of material. The 3s and 4s states indirect magnetoexcitons bound by the Coulomb potential dissociate at a higher magnetic field than magnetoexcitons bound by the RK potential.
Related to bilayer MX 2 -MX 2 : the magnetic field's contribution to the magnetoexcitons binding energies is systematically bigger for the bilayer than for the FS MX 2 monolayer for all states when the Rytova-Keldysh potential is used in calculations. The same conclusion is extended to DMCs.
A distinct feature for the diamagnetic coefficients in bilayer and van der Waals molybdenum-based heterostructures is that σ is always smaller for the tungsten-based heterostructures due to the difference of the valence and conduction bands in MoX 2 and WX 2 monolayers. The other distinct feature is that within the linear regression model with R 2 = 0.9998 the diamagnetic coefficients can be extracted for the 1s, 2s, 3s, and 4s states for the bilayers, while the DMCs exist in the van der Waals heterostructure with up to six hBN layers for 1s and three hBN layers, for 1s and 2s states, respectively. The linear dependence of ∆E on B 2 is invalidated in 3s, and 4s states for both the Rytova-Keldysh and Coulomb potentials.
Finally, our results raise the possibility of controlling the binding energies of direct and indirect magnitoexcitons in monolayer MX 2 , bilayer MX 2 -MX 2 and MX 2 -hBN-MX 2 heterostructures using the external magnetic field and open the additional degree of freedom to tailor the binding energies and diamagnetic coefficients for the MX 2 -hBN-MX 2 heterostructures by varying the number of hBN sheets between TMDC layers.  A and B excitons. For the A exciton for each material two values of µ and χ2D are given. The values of A high and A low correspond to the parameters found in the literature, which maximize and minimize the A exciton binding energy, respectively The average dielectric constant, κ = ε 1 +ε 2 2 = 4.89, is taken from Ref. [17]. κ is the same for all materials. µ is in units of electron mass, m0.The 2D polarizability χ2D and TMDC monolayer thickness h are given inÅ.
A high 0.27 [  In Fig. 8 the energy contributions from the magnetic field to the binding energies of indirect magnetoexcitons in tangerine-based WSe 2 -hBN-WSe 2 heterostructure are presented. The comparison with the results with the indirect magnetoexcitons in MoS 2 -hBN-MoS 2 heterostructure presented in Fig. 7 shows that the differences occur only in the energy contribution magnitudes, which varies within up to 40% for MoSe 2 -hBN-MoSe 2 for the maximum value, and indirect magnetoexcitons do not dissociate in the 3s state and dissociate at larger values of the magnetic field in 4s state.  Appendix D: Diamagnetic coefficients