Interlayer excitons in van der Waals heterostructures: Binding energy, Stark shift, and field-induced dissociation

Photoexcited intralayer excitons in van der Waals heterostructures (vdWHs) with type-II band alignment have been observed to tunnel into interlayer excitons on ultrafast timescales. Such interlayer excitons have sufficiently long lifetimes that inducing dissociation with external in-plane electric fields becomes an attractive option of improving efficiency of photocurrent devices. In the present paper, we calculate interlayer exciton binding energies, Stark shifts, and dissociation rates for six different transition metal dichalcogenide (TMD) vdWHs using a numerical procedure based on exterior complex scaling (ECS). We utilize an analytical bilayer Keldysh potential describing the interaction between the electron-hole pair, and validate its accuracy by comparing to the full multilayer Poisson equation. Based on this model, we obtain an analytical weak-field expression for the exciton dissociation rate. The heterostructures analysed are MoS2/MoSe2, MoS2/WS2, MoS2/WSe2, MoSe2/WSe2, WS2/MoSe2, and WS2/WSe2 in various dielectric environments. For weak electric fields, we find that WS2/WSe2 supports the fastest dissociation rates among the six structures. We, furthermore, observe that exciton dissociation rates in vdWHs are significantly larger than in their monolayer counterparts.


interlayer excitons
A bilayer vdWH supports two distinct types of excitons. The electron and hole may either be localized within the same layer, or they may reside in different layers. These two cases are referred to as intra-and interlayer excitons, respectively. The vdWHs considered in the present paper are bilayers with type-II band alignment. We name the structures TMD1/TMD2 so that the conduction band minimum and valence band maximum reside in the first and second layer, respectively. For the structures considered, the energy won from band offsets by the electron and hole residing in the first and second layer, respectively, is larger than the loss in exciton binding energy. The many-body excitonic ground state is therefore an interlayer exciton. Intralayer excitons in either layer are thus excited states, and will therefore tunnel into the ground state. Direct photoexcitation of interlayer excitons by resonance photons is not as likely as for intralayer excitons due to the weak overlap of the electron and hole wave functions. However, resonantly excited intralayer excitons in either layer quickly transition into interlayer excitons 35 that, therefore, become very important for many properties of vdWHs. To describe excitonic effects from first principles, one must turn to the many-body Bethe-Salpeter equation 41,42 . Solving it is a computationally demanding task even for simple structures. Fortunately, under well-defined approximations, the many-body problem can be simplified to the Wannier equation 43,44 , essentially reducing it to a Schrödinger-type problem with a hydrogenic Hamiltonian. The Wannier model has indeed been shown to yield a sufficiently accurate description of many excitonic properties 29,[45][46][47][48] . In terms of the relative coordinate r = r e − r h of the electron-hole pair, it reads (atomic units are used throughout) h is the reduced exciton mass, V is the screened Coulomb interaction between the electron and hole, and E E is the electric field. The electric field is taken to point along the x-axis throughout the paper, i.e. = E E e x  . As the valence band maximum and conduction band minimum at the K point are primarily composed of the d orbitals of the metal atoms 49,50 , electrons and holes will, to a good approximation, reside in the middle of their respective layers. We are therefore able to freeze their out-of-plane motion, which effectively makes solving Eq. (1) a two dimensional problem. Figure 1 shows an illustration of an interlayer exciton in a bilayer vdWH subjected to three different field strengths: zero (a), weak (b), and strong (c). When an electric field is present, the electron and hole will be pulled in opposite directions. For weak electric fields, the probability of dissociating the exciton is low. It will therefore become polarized, but most likely recombine rather than dissociate. In strong electric fields, however, field induced dissociation becomes likely, and dissociation rates may become extremely large as the field strength increases.
To describe the interaction V between the electron-hole pair, we fix the electron and hole to the middle of their respective sheets. The vertical separation between the electron and hole = − for intra-and interlayer excitons, respectively, where d 1 and d 2 are the thicknesses of the two TMD sheets (see Fig. 1). The van der Waals heterostructure is then modeled as the four layer system in Fig. 2. Interlayer exciton in a bilayer van der Waals heterostructure with zero external field (a), a weak inplane field polarizing the exciton (b), and a strong in-plane field dissociating the exciton (c).

Scientific RepoRtS |
(2020) 10:5537 | https://doi.org/10.1038/s41598-020-62431-y www.nature.com/scientificreports www.nature.com/scientificreports/ Here, the dielectric function ε z ( ) is taken to be piecewise constant given by ε a , ε 1 , ε 2 , and ε b in the superstrate, upper TMD sheet, lower TMD sheet, and substrate, respectively. By Fourier decomposing V and solving the multilayer Poisson equation for the Fourier components (see App. A), one may obtain V expressed as an integral in momentum space , as expected. For q → ∞, however, the dielectric function describing intralayer excitons tends to the dielectric constant of the layer to which they are confined, i.e. to ε 1 and ε 2 for intralayer excitons in the first and second layer, respectively. On the other hand, the function for interlayer excitons tends to the average dielectric constant of the two layers ε ε + ( )/2 1 2 . The complete interlayer dielectric function is obtained in App. A and is given by Eq. (14). Appendix A also explains how to obtain the intralayer function. The full potential may readily be obtained for real r by using standard numerical integration techniques in Eq. (2). However, in the present paper, we seek to calculate dissociation rates by using exterior complex scaling (ECS) [51][52][53][54] . This implies rotating the radial coordinate into the complex plane outside a radius R by an angle φ, i.e.
It is a simple task to show that the integrand in Eq. (2) will become an exponentially increasing oscillating function when φ > + r d R /sin by using the integral representation for the Bessel function 55 . This makes Eq. (2) extremely difficult to handle while using ECS. Fortunately, the numerical solution is very accurately approximated by a bilayer Keldysh (BLK) potential. The Keldysh potential has been used extensively to describe excitons in monolayer TMDs. The monolayer Keldysh (MLK) interaction is given by 56,57  where H 0 is the zeroth order Struve function 55 , Y 0 is the zeroth order Bessel function of the second kind 55 , is the average dielectric constant of the surrounding media, and r 0 is the screening length proportional to the polarizability of the sheet 47 . This potential diverges logarithmically at the origin. On the other hand, the potential describing interlayer excitons is finite at the origin due to the vertical separation between the electron-hole pair.
To obtain the bilayer Keldysh potential, we substitute → + r r d 2 2 and → + r r r 0 0 (1) 0 (2) into Eq. (4) which yields Here, r 0 (1) and r 0 (2) are the screening lengths of the first and second monolayer, respectively. The first substitution accounts for the possible vertical separation between electrons and holes in bilayer structures. The second www.nature.com/scientificreports www.nature.com/scientificreports/ substitution accounts for the increased thickness of the structure, as the total screening length is proportional to it 56,57 . It should be noted that making the first substitution without the second leads to a potential that is far too strongly binding 58 . Note, further, that to obtain the interaction for intralayer excitons in a bilayer structure, we use the BLK potential with d = 0. The screening lengths used in the present paper are ab initio values from ref. 21 . In Fig. 3, we compare the bilayer Keldysh potential to the full potential obtained by solving the multilayer Poisson equation. Panel (a) shows the potential for the representative WS 2 /WSe 2 case. Evidently, a good agreement is found. Both the intra-and interlayer potentials (exact and approximate) behave as −1/κr for r ≫ d, as the vertical separation becomes negligible in this region. For small r, on the other hand, the intralayer potential diverges logarithmically, while the interlayer potential behaves as a + br 2 , where a and b are constants. The quadratic form is readily understood by expanding the Bessel function in Eq. (2) to second order, i.e.
Integrating using the second order expansion shows that the integral diverges for d = 0. For d > 0, on the other hand, the exponential function leads to converging integrals even for r = 0. The quadratic behaviour has inspired the use of a harmonic oscillator model, obtained by expanding a potential similar to Eq. (5) to second order, to analyse interlayer excitons 59 . It should be noted, however, that the binding energies predicted by a second order expansion of Eq. (5) are in poor agreement with those obtained using the full potential. In contrast, the binding energies obtained with the BLK potential are in good agreement with the full potential results, as seen in panel (b) of Fig. 3. In fact, they never deviate more than 6% for the cases considered. The binding energies are, furthermore, in excellent agreement with those found in literature. As an example, the binding energy of 1s interlayer excitons in WS 2 /WSe 2 on a diamond substrate (ε a = 1 and ε b = 5) was measured recently to be 126 ± 7 meV 35 , and our model yields 125 meV. Furthermore, Ovesen et al. 60 found the binding energy of interlayer excitons in MoSe 2 /WSe 2 in free space to be 246 meV using a model similar to our full potential, where we find 260 meV with the BLK model. It should be mentioned that the model used in the present paper predicts slightly lower free-space binding energies than some of the ab initio methods [61][62][63][64] . Interlayer exciton binding energies for the six TMD bilayer combinations with type-II band alignment 65,66 are summarized in Table 1 for various dielectric environments. The effective masses used are obtained from ref. 50 , and the TMD widths have been taken as half the vertical lattice constants found in ref. 67 .

field induced dissociation
When the exciton is subjected to an electrostatic field, it may be dissociated. This is realized in the Wannier model by the energy eigenvalue obtaining a non-vanishing imaginary part [27][28][29]68 . The field induced dissociation rate is then given by Γ =− E 2 Im , where the imaginary part can be obtained efficiently by utilizing the ECS technique [51][52][53][54] . As mentioned briefly in the previous section, this technique consists of rotating the radial coordinate into the complex plane outside a radius R (see Eq. (3)). The partitioning of the radial domain is efficiently dealt with by resolving the radial part of the eigenstate in a finite element basis consisting of Legendre polynomials, and the angular part in a cosine basis. We have previously used the same numerical procedure to calculate dissociation rates for excitons in various monolayer TMDs 30 , and we refer the interested reader to that paper for the technical details of the method. The field induced dissociation rates and Stark shifts for interlayer excitons in the six van der Waals heterostructures are shown in Fig. 4. As is evident, the structures support excitons that behave very  48 . The reason that these large polarizabilities are observed for interlayer excitons in bilayer vdWHs is the increased screening and the vertical separation of the electron and hole. Both reduce the binding energy and they are therefore much easier to polarize.
Turning to the dissociation rates, it is clear that the encapsulating media significantly alter how quickly interlayer excitons are dissociated. The same behaviour was observed for their monolayer counterparts 30 , and one therefore has a large degree of freedom in device design. For extremely weak fields, the dissociation rates are very low, but they grow rapidly with an increasing field strength. As an example, the dissociation rate of interlayer excitons in freely suspended MoS 2 /WS 2 is around Γ ≈ 1.7 × 10 4 s −1 already at = μ E 10 V/ m, and only Γ ≈ 5.3 × 10 −38 s −1 and Γ ≈ 2.7 × 10 −33 s −1 for monolayer MoS 2 and WS 2 , respectively 30 . It should be noted that dissociation rates of intralayer excitons are only important if they are comparable to the rate at which intralayer excitons tunnel over to interlayer excitons. In a recent experiment on MoS 2 /WS 2 structures, the holes of photoexcited excitons in the MoS 2 layer of this structure were observed to tunnel into the WS 2 layer within 50 fs 31 .   Table 1. Interlayer exciton binding energy | | E 0 in meV and the field-independent material front factor Γ 0 of the dissociation rate in atomic units for the six vdWHs in various dielectric environments κ ε . The reduced interlayer exciton mass μ is indicated for each heterostructure. = μ E 10 V/ m, we find Γ ≈ 2.1 × 10 −3 s −1 and Γ ≈ 2.9 × 10 4 s −1 , respectively. The large difference between these two rates can be traced back to the reduced masses, which are μ = 0.2513 and μ = 0.1543, respectively. The time it takes to dissociate these excitons with the given field may be approximated as τ = 1/Γ ≈ 476 s and τ ≈ 3.5 × 10 −5 s, respectively. The intralayer excitons have therefore clearly transitioned to interlayer excitons before they are dissociated by the field. Interlayer tunneling rates are likely affected by material parameters as well as the surrounding media. Assuming, however, similar time scales as observed for MoS 2 /WS 2 Γ Tunnel ≈ 10 13 s −1 , interlayer dissociation rates are the limiting factor in field induced exciton dissociation for weak to moderate fields. For the largest fields in Fig. 4, the competition between tunneling and dissociation will be important for an accurate description. Due to the risk of dielectric break-down, such large fields are best avoided in devices, however. The high interlayer dissociation rates suggest that using carefully chosen bilayer TMDs in photocurrent devices is much more attractive than their monolayer counterparts. Moreover, proper encapsulation will further improve device performance. As hBN is a very common material used to encapsulate samples, we show the interlayer dissociation rates for WS 2 /WSe 2 that is either (i) freely suspended, (ii) placed on an hBN substrate, or (iii) encapsulated by hBN in Fig. 5. Evidently, hBN surroundings increase the dissociation rates by several orders of magnitude, and for weak fields in particular.
Recently an analytical weak-field approximation was obtained for exciton dissociation rates in monolayer TMDs 30 . The derivation was made using weak-field asymptotic theory 70 , and the fact that the MLK potential has a simple asymptotic form. As mentioned in the previous section, the BLK potential has exactly the same asymptotic behaviour. The weak-field approximation for interlayer excitons therefore has exactly the same form, albeit with a different field-independent front factor. We arrive at is a field-independent material constant. The parameters needed to use Eq. (7) are presented in Table 1, where Γ 0 has been computed by the integral procedure in 30 . Panel (a) of Fig. 6 shows the interlayer dissociation rates as functions of environment screening κ for a field of 5 V/μm computed with the numerical procedure (dots) and the weak-field approximation (solid lines). Evidently a good agreement is found for such a low field. Nevertheless, as the inset shows, the weak-field approximation quickly begins to overestimate the dissociation rate for larger fields. This was also found for monolayer TMDs 30 . It is clear that WS 2 /WSe 2 supports excitons with the largest dissociation rates at = μ E 5 V/ m, making it an interesting choice in device design. Note that the fully numerical procedure breaks down if the dissociation rate becomes extremely small, and, hence, we are unable to use it to obtain dissociation rates for MoS 2 /MoSe 2 and MoS 2 /WS 2 in surroundings with (very) low screening for this field strength. It is therefore advantageous to have a formula such as Eq. (7) when the fields become sufficiently weak.
In panel (b) of Fig. 6, we show the binding energies for the same excitons as functions of κ. They clearly follow similar trends as the dissociation rates do, suggesting that binding energy has a significant impact on dissociation rates. This is to be expected, as strongly bound excitons are harder to pull apart. It should, however, be noted that the binding energy does not uniquely determine the dissociation rate. As an example, the interlayer excitons in MoS 2 /MoSe 2 and MoS 2 /WS 2 have identical binding energies for surrounding media with κ ≈ 2.8. Nevertheless, the dissociation rates are clearly different. In fact, several crossings are observed in the binding energies as a function of κ whereas only one is found in the dissociation rates at 5 V/μm. The origin is the different reduced exciton masses. For no structure does Γ 0 nor μ k / contain any crossing for κ ∈ [1; 5] . However, μ k 3 crosses for MoS 2 / WSe 2 and MoSe 2 /WSe 2 . This complicated behaviour suggests that, to obtain accurate estimates of interlayer exciton dissociation rates, one must turn to either the weak-field approximation or make a full numerical calculation.

conclusion
In the present paper, we have studied binding energies, Stark shifts, and dissociation rates of interlayer excitons in van der Waals heterostructures (vdWHs). The structures analysed are the six bilayer vdWHs with type-II band alignment arising from combinations of MoS 2 , MoSe 2 , WS 2 , and WSe 2 . The bilayer excitons are described using an analytical bilayer Keldysh potential, which we have verified for accuracy by comparing to the full solution of the multilayer Poisson equation. Exciton binding energies, Stark shifts, and dissociation rates can therefore readily be calculated using the analytical potential. We find interlayer exciton binding energies ranging from 256 to 280 meV for freely suspended structures and from 78 to 92 meV for heavily screened structures, making them stable at room temperature. Furthermore, both the polarizabilities and dissociation rates found for these excitons are much larger than their monolayer counterparts. For example, interlayer excitons in freely suspended MoS 2 / WS 2 are found to dissociate at a rate of Γ ≈ 1.7 × 10 4 s −1 in a field strength of 10 V/μm whereas monolayer MoS 2 and WS 2 have Γ ≈ 5.3 × 10 −38 s −1 and Γ ≈ 2.7 × 10 −33 s −1 , respectively.
For moderate field strengths, intralayer exciton dissociation rates are significantly lower than the rate at which such excitons tunnel into interlayer excitons. For this reason, interlayer exciton dissociation rates are the limiting factor in generation of photocurrents at weak to moderate fields. Since optically excited excitons in one of the layers tunnel to long-lived interlayer excitons on ultrafast timescales, bilayer vdWHs with favourable band offsets may potentially serve as building blocks in efficient photocurrent devices. Finally, the numerically exact dissociation rates are compared to an analytical weak-field dissociation formula obtained from weak-field asymptotic theory. A good agreement is found in the weak-field limit, and Eq. (7) therefore serves as a useful formula to quickly estimate field induced dissociation rates of excitons in bilayer vdWHs for weak electric fields. in the respective regions. The Fourier components can then be found analytically by solving the system of equations that arises from the boundary conditions. To describe charges confined to distinct layers, we fix the electron and hole to z = d 1 /2 and = − ′ z d /2 2 , respectively. This leads to the interlayer exciton potential z z q 0 is the bare interaction. The effective dielectric function is given by eff with ε γ ε ε γ ε ε γ ε ε γ for q → 0 and to ε ε + ( )/2 1 2 for q → ∞, as expected. The dielectric function describing charges confined to the same layer may be obtained in a similar manner by placing both charges at z = − d 2 /2. The interaction in real space can then be obtained as the inverse Fourier transform ∫ ε − = .
∞ − − ′ V r e J qr q dq ( ) ( ) ( ) (18) z z q 0 0 eff For the dielectric constants, we have used the static in-plane dielectric constants calculated from first principles in ref. 71 .