Exciton ionization in multilayer transition-metal dichalcogenides

Photodetectors and solar cells based on materials with strongly bound excitons rely crucially on field-assisted exciton ionization. We study the ionization process in multilayer transition-metal dichalcogenides (TMDs) within the Mott-Wannier model incorporating fully the pronounced anisotropy of these materials. Using complex scaling, we show that the field-dependence of the ionization process is strongly dependent on orientation. Also, we find that direct and indirect excitons behave qualitatively differently as a result of opposite effective anisotropy of these states. Based on first-principles material parameters, an analysis of several important TMDs reveals WSe2 and MoSe2 to be superior for applications relying on ionization of direct and indirect excitons, respectively.


Introduction
Transition-metal dichalcogenides (TMDs) including MoS 2 , MoSe 2 , WS 2 , and WSe 2 are layered twodimensional semiconductors with unique electronic and optical properties. They are highly promising for optoelectronic applications such as photodetectors [1][2][3][4][5], solar cells [6,7], and light emitting diodes [8]. In their monolayer form, MoS 2 , MoSe 2 , WS 2 , and WSe 2 are direct bandgap semiconductors. Importantly, the low dimensionality and reduced screening leads to highly prominent exciton effects with binding energies of several hundred meVs [9,10]. Such excitons greatly modify both linear [11] and nonlinear [12] optical properties. Regarding applications, exciton binding energies significantly larger than the thermal energy at room temperature (~25 meV) increase radiative electron-hole recombination and are, therefore, beneficial for efficient light emission. In contrast, strongly bound excitons may reduce the efficiency of photodetectors and solar cells because these devices require exciton ionization in order to separate electrons and holes. The first monolayer TMD-based photodetectors showed relatively low efficiencies [1]. By improving material quality and thereby increasing carrier mobility, the responsivity was subsequently greatly improved [2,4]. These devices all operated in the parallel collection mode, in which an in-plane bias between metal contacts drives the current along the monolayer. Recently, devices based on perpendicular collection, i.e. transport between layers, have emerged as promising alternatives. In particular, both efficient and ultrafast photoresponse has been demonstrated for MoS 2 [3] and WSe 2 [5] photodetectors. Importantly, this approach allows for devices based on stacking of appropriate two-dimensional materials. Thus, contacts for carrier collection can be fabricated by encapsulation of the photoactive semiconductor between conducting graphene sheets [3,5].
The highly promising characteristics demonstrated in [3,5] were obtained for multilayer samples having thicknesses of 50 nm [3] and between 2.2 and 40 nm [5], respectively. In such samples, electron and holes are delocalized across several layers and, moreover, screening is increased compared to monolayers. Hence, in slabs thicker than the bulk exciton Bohr radius, the exciton binding energy is significantly reduced. This is expected to contribute to the high sensitivity and response rate observed experimentally. The limiting factors for the rate are poorly understood, however. In [5], the response rate was shown to increase approximately linearly with bias voltage between the graphene contacts while an inverse quadratic dependence on sample thickness was found.
This suggests that out-of-plane drift is the limiting factor for the rate [5]. It is clear, however, that a full description of the response involves several physical mechanisms. In fact, the photoresponse itself is a combined process consisting of (i) exciton creation upon photon absorption, (ii) field-assisted ionization of the exciton, (iii) transport of dissociated carriers, and (iv) collection at the contacts. Hence, the applied bias performs a twofold function of driving both exciton ionization and carrier transport.
The process of exciton ionization is a crucial step in the photoresponse. In materials with strongly bound excitons, such as TMDs, thermal ionization is inefficient and a strong external electric field is required for efficient carrier separation. As an estimate, the required field strength for efficient dissociation is given by the ratio between exciton binding energy and Bohr radius. A full description of the ionization process is needed, however, to model the dependence of the ionization rate on field strength. Moreover, in anisotropic materials, the ionization rate will depend on the direction of the applied field. In the present paper, we study the exciton ionization process in anisotropic TMDs using a modified Mott-Wannier [13] approach incorporating the material anisotropy. In this approach, bulk dielectric constants and effective masses lead to an exciton eigenvalue problem that is mathematically identical to that of a hydrogen atom embedded in an anisotropic material. The material constants are calculated from first-principles density-functional theory. To describe field-assisted ionization, an electric field is added to the Mott-Wannier equation. We apply the complex scaling technique [14] and hypergeometric resummation [15,16] to compute ionization rates and Stark shifts of excitons in typical TMDs. These are highly non-perturbative phenomena. In particular, the ionization rate cannot be described using a finite-order perturbation expansion [14]. Our work is related to a recent work [17], in which the secondorder Stark shift in phosphorene was studied. This two-dimensional material has strong in-plane anisotropy. Also, in TMDs, the optical Stark effect, i.e. non-perturbative effects of strong laser excitation, has been demonstrated [18,19]. In addition, a recent study by some of the present authors has described exciton dissociation in monolayer MoS 2 , demonstrating a pronounced dependence on screening by the surroundings [20]. Very recently [21], experimental Stark shifts due to static perpendicular fields in MoS 2 multilayers with thicknesses between one and five layers were reported. The weak thickness dependence of the observed shift clearly indicates that exciton effects are important. So far, however, exciton ionization due to static fields in multilayer TMDs has not been discussed in any quantitative theoretical work.

Anisotropic exciton model
In a truly two-dimensional material, the nonlocal wave vector dependence of the dielectric constant is important for a correct description of screening [10]. In a three-dimensional material, however, dispersion is much less pronounced and the dielectric constant can be assumed independent of wave vector, i.e. approximated by the long-range limit e  ( ) q 0  [22]. In uniaxially anisotropic three-dimensional materials such as multilayer TMDs, the dielectric constant is then given by a constant tensor e e e e = ( ) diag , , x x z  with separate elements for the out-of-plane (z-axis) and in-plane (x-axis) values. Similarly, the effective masses for these directions differ. For direct excitons, electrons and holes are located at the K point. For indirect excitons, holes reside at the Γ point whereas electrons are located at the conduction band minimum Σ roughly midway between Γ and K. In fact, the in-plane effective mass depends on rotation angle in the (x, y) plane. This relatively weak in-plane dependence will be ignored, however. Thus, below we take the effective hole masses computed for the  , ,    As illustrated in figure 1, these scaling factors are the appropriate ones for the cases of perpendicular collection (  along z) and parallel collection (  along x), respectively. In the effective exciton units, the consequences of anisotropy are entirely governed by the single parameter k. Hence, if k = 0, the material is effectively isotropic. In contrast, k > 0 means that, in scaled coordinates, the Coulomb attraction is less sensitive to the z coordinate. In fact, as k approaches unity, the z-dependence of the interaction vanishes. Hence, for k > 0 excitons will tend to delocalize perpendicular to the layers. If k < 0, delocalization within individual layers is enhanced. These trends are illustrated by the schematic insets in figure 2. We stress that this picture is only valid in scaled coordinates. As shown below, the sign of k differs for direct and indirect excitons. Hence, the large effective masses for the z-direction for direct excitons mean that k > 0 for these states in TMDs. In contrast, indirect excitons have nearly isotropic effective masses and it is mainly the dielectric constants that lead to anisotropy. The fact that e e > x z consequently means that k < 0 for indirect excitons.
To solve the eigenvalue problem, we introduce cylindrical polar coordinates r q { } z , , and expand in a Laguerre-type basis The parameters k and q can be optimized to minimize the exciton ground state energy. We find however, that = = k q 2 is very nearly the optimal choice in all cases. For = z,  cylindrical symmetry is preserved around the z-direction and for the exciton ground state = L 0 is sufficient in the expansion above. On the other hand, for = x, Hence, the size of the basis is 1800 and 2400 for the two cases, respectively.
In the limits k = 0 and k = 1, the Wannier equation describes three-dimensional and two-dimensional excitons, respectively. We note, however, that for k = 1 the exciton state is actually completely delocalized along the z-direction and the 2D behavior pertains to the in-plane character only. This follows from the behavior of the potential x y 2 2 1 2 2 2 1 2 as the limit k = 1 is approached. Hence, the in-plane and out-ofplane motions decouple completely leading to delocalization along z. The 2D character of the exciton in this limit is, however, rather different from excitons in monolayer TMDs because of differences in screening, which has a pronounced nonlocal wave vector dependence for monolayers, as discussed above. In figure 2, we show the binding energy E 0 of the ground state as a function of the anisotropy. As expected, the limits for k = 0 and respectively. Importantly, though, for k = 1 the ground state is, in fact, the lower limit of a continuum reflecting the delocalized z-dependence. In contrast, the state is fully localized for k = 0. In figure 2, we also include the perturbation result valid for small k, i.e. for materials that are only weakly anisotropic, see appendix and [25].

Exciton ionization
As mentioned above, field-ionization is a non-perturbative phenomenon that is not captured in any finite-order perturbation expansion in field strength. In the presence of the field, the bound state energies turn into complex resonances. Writing i 2we may decompose such complex eigenvalues into their real and imaginary parts. The interpretation is then that the real part D represents the Stark energy whereas the imaginary part provides the ionization rate G. This rate can be understood as the rate of tunneling from the bound state into the dissociated state. There are two commonly sought routes to obtaining G, both of which we will pursue below. The first is the purely numerical approach of complex scaling [14,20] that immediately provides both real and imaginary parts of the resonance. Alternatively, a semi-analytical result can be obtained through analytical continuation and resummation of a low-order perturbation series [15,16,26].
In a direct bandgap material, the excitation and ionization processes are conceptually simple. Hence, for excitation by photon energies close to the fundamental exciton, the initial excitation quickly relaxes to the lowest direct exciton. From here, the state then either recombines (radiatively or non-radiatively) or becomes ionized. In contrast, an indirect bandgap material is more complicated. Assuming again excitation by low energy photons and ignoring phonon-assisted processes, the initial excitation is the direct exciton as the transition to the indirect one is optically forbidden by momentum conservation. Hence, if ionization happens sufficiently rapidly, the dissociating species is still the direct exciton. Conversely, ionization after thermal relaxation to the indirect exciton means that this species is the relevant one. It follows that a rapid ionization rate, i.e. ionization in a strong field, pertains to the direct exciton. Slow ionization in weak fields, on the other hand, will predominantly happen from the indirect exciton. In the present work, we will study ionization of both species. As our focus is on strong field ionization, however, we will mainly illustrate results relevant for direct excitons, i.e. the k > 0 regime.
If a weak electric field is applied, the resonance is approximately x z x z , 0 , 2 2 where subscripts x and z indicate parallel and perpendicular collection, respectively. The coefficient of the second order term determines the exciton polarizability a x z , through a = - ( ) E .
x z x z , , More generally, if a perturbation expansion in the electric field is applied, the ground state energy is written  direction. Using a finite-field approach, we solve the Wannier equation equation (3) in the presence of a small electric field of magnitude = - 10 3 applying the localized bases introduced above. By subtracting the field-free ground state energy, we obtain the direct exciton k > ( ) 0 polarizabilities shown in figure 3. As shown in the figure, the behavior for highly anisotropic materials is strongly dependent on the direction of the applied field. The two curves start from a common value of 9/2 in agreement with the exact value for the hydrogen atom [27,28]. For values of k approaching unity, however, their behaviors are markedly different. Thus, in the perpendicular case, the polarizability diverges reflecting the completely delocalized state along the z-direction. In contrast, in the parallel case, the polarizability remains finite and eventually reaches a value of / » 21 128 0.1641 in the 2D limit [29]. Again, in the plot, we include the approximate results found by linearizing in k as derived in the appendix, i.e. / / a k » -9 2 439 120 x and / / a k » -9 2 101 60. z It is seen that these describe the exact behavior surprisingly well even for anisotropies as large as k~0.8. By solving for several distinct (small) field strengths, the finite-field approach can be extended to provide higher order terms in the field expansion. In this way it can be shown that linearization in k remains reasonably accurate for higher order terms as well.
In fact, the full asymptotic series å = ¥ ( )  E n x z n n 0 , 2 2 is highly divergent for all finite values of the field strength. However, by appropriate resummation of a finite series, physically meaningful quantities can be obtained from it. Traditionally, Padé approximants have been applied to this end [26]. These typically require partial expansions of high order to be successful. We have recently shown [15,16] that hypergeometric resummation provides a very efficient alternative. Here, only the first five non-vanishing terms   n 0 4of the expansion are needed to find a highly accurate result. We have previously applied this result to the three-dimensional [15] and low-dimensional [16] hydrogen problems, which are mathematically very similar to the Wannier problem considered here. In the appendix, the required series are provided for both collection modes and below we demonstrate how these may be used to accurately describe both Stark energy and exciton ionization.
The complex scaling method [14] is based on a coordinate scaling  j r re .   By analytically continuing into the complex plane, the parameter j can be taken purely imaginary, i.e. j q = i with q real. In turn, the complex scaled Wannier equation reads as For finite field and rotation angle q, the eigenstates are square integrable and the eigenvalues complex, as explained above. In fact [14], the eigenvalue is independent of the value of q as long as a finite value is adopted. In practise, a small q-dependence is observed whenever expansion in a finite basis is applied. We have found, however, that this dependence is negligible if bases of the sizes discussed above are used. Below, all results will be for q = 0.4.
In figure 4, the solid lines show the complex scaling results for four values of the anisotropy ranging from none k = ( ) 0 to substantial (k = -1.0 and k = ) 0.8 . Moreover, both collection modes are analyzed. At vanishing field strength, the Stark energy, i.e. the real part of the resonance, agrees with the results in figure 2. As the field is increased, the energy initially decreases quadratically with a prefactor given by the polarizability in figure 3. However, beyond a field strength of approximately 0.1, a significantly softer behavior is found. At a similar field strength, the ionization rate increases dramatically before becoming approximately linear in the field. In weak fields, the ionization rate is dominated by an exponential G~-( )  c exp behavior. All of these features agree with the analogous findings for atomic hydrogen [15]. It is observed that increased anisotropy leads to reduced ionization, mainly as a result of the increased exciton binding energy. Moreover, in the perpendicular collection mode, the ionization rate in relatively large fields is only weakly sensitive to the anisotropy, i.e. all curves share roughly the same slope. In contrast, for the parallel case, the slope decreases markedly as k is increased. This difference relates to the delocalization behavior of the states as shown in the insets of figure 2. Thus, as k increases, the states become increasingly delocalized in the perpendicular direction. This, taken by itself, increases ionization in a perpendicular field but not in a parallel one. Hence, delocalization partially counteracts the increased exciton binding for the perpendicular collection mode.
The complex scaling results are based on diagonalization of relatively large matrices on a fine grid of electric field strengths. Such a fine grid is required in order to track the evolution of a particular state as the field is increased. This approach is therefore computationally demanding and it is of interest to compare with the much simpler hypergeometric resummation approach. As explained above, the method takes the first five terms in the asymptotic expansion of the field dependence as input. In the appendix, the required expansions are provided. We apply these within the hypergeometric resummation approach used in [16] and thereby obtain Stark energies and ionization rates at practically no computational cost. These are shown by the dots in figure 4. For vanishing anisotropy, the agreement is essentially perfect, as expected [15,16]. Remarkably, however, the agreement is excellent even for large values of k | |. This demonstrates the power of the hypergeometric resummation approach. We stress that, using the series provided in the appendix, accurate calculation of resonances for any given material (within the class studied here) can be made at a fraction of the computational cost compared to the full complex scaling.

Application to TMDs
To convert any result of the Wannier approach to physical quantities we require specific values of the anisotropic dielectric constants and effective masses. To this end, we have performed first-principles calculations for the important TMDs MoS 2 , MoSe 2 , WS 2 , and WSe 2 using the Random Phase Approximation and including spinorbit interaction. All calculations were performed with the projector augmented wave electronic structure code GPAW [30,31] and include local-field effects in the dielectric response. Full details on the calculations can be found in [32]. The first-principles results for direct and indirect excitons are listed in table 1 and 2, respectively. For direct excitons, the anisotropy manifests itself in highly different out-of-plane and in-plane effective masses. As seen in table 1, these masses differ by nearly an order on magnitude. Similarly, the out-of-plane dielectric constants are less than half the in-plane values. For both sulphides MoS 2 and WS 2 , these values translate into an anisotropy of k » 0.8. Thus, even if the two materials are obviously dissimilar, their effective anisotropy is essentially the same. The dissimilarity reveals itself in the field scaling factors, however, which are about twice as large in MoS 2 as in WS 2 as a result of the much larger effective masses of the former. This means that any normalized field strength in figure 4 translates into a much smaller physical field in WS 2 relative to MoS 2 . An approximately similar ratio is found for the field scaling of MoSe 2 relative to WSe 2. To complete the conversion and translate the ionization rate G into physical units, the frequency scale / ⁎  Ha is applied. The converted results for direct and indirect excitons obtained from this procedure are shown in figures 5 and 6, respectively.
In the effective exciton units, a field of unity magnitude corresponds to a physical field of x z x z x z , , 0 , i.e. the ratio between effective Hartree and Bohr radius for the particular field direction. From tables 1 and 2 it is seen that the direct and indirect exciton binding energies constitute about 75% and 40% of ⁎ Ha , respectively. Hence, estimating the field strength required for efficient ionization as / | | ⁎ E a x z 0 , is equivalent to unity field strength 1 in exciton units within a factor of approximately two. From   figure 4 it is seen that at such a field strength, ionization is well into the linear regime. Hence, using , as an estimate for the ionization threshold is partly supported by our results. In fact, the transition from exponential to linear field-dependence in figure 4 occurs around 0.1, which may be considered a better measure for the threshold.
For both direct and indirect excitons, the magnitude of the ionization rate reflects the exciton binding energy. Hence, there is a simple correlation between the binding energies in tables 1 and 2, on the one hand, and the ordering of the curves in figures 5 and 6 on the other, with strongly bound excitons leading to suppressed ionization. The larger differences for the direct exciton binding energies lead to more spread-out curves as compared to the indirect ones. Overall, the best candidates for efficient ionization of direct and indirect excitons are WSe 2 and MoSe 2 , respectively. For the direct excitons (figure 5) the WSe 2 rate is significantly higher than the second best material, i.e. MoSe 2 and WS 2 for perpendicular and parallel collection, respectively. In comparison, the differences for indirect excitons are modest.
We end this section with a brief comparison to the experimental photoresponse rate observed in [5]. In that work, the photoresponse rate was determined using two-pulse excitation with an adjustable delay. Hence, the extraction time of the photogenerated carriers produced by the first pulse is probed by the second pulse. The extraction time itself is a measure of at least three processes occurring in series: (i) exciton ionization, (ii) drift to the TMD/graphene interface, and (iii) transfer across the interface. Hence, importantly, the measurement is dominated by the slowest among these processes and the experimental rate cannot necessarily be identified with the actual ionization rate. The fastest response~⋅ -2 10 s 11 1 was seen for a WSe 2 device having a slab thickness of 2.2 nm. This thickness is larger than the perpendicular exciton Bohr radius (1.1 and 1.4 nm for direct and indirect excitons, respectively) and we therefore expect the bulk picture in the present work to be applicable. Bias voltages up to 1.2 V were applied in the measurement and the response rate typically increased with applied bias. The precise value of the internal electric field is subject to some uncertainty and depends on the geometrical capacitance and interface charges (see [5], supplementary information). It is clear, however, that the theoretical ionization rates computed above exceed the measured rate by a significant factor. In fact, for field strengths of m -100 V m 1 and above, the converted results in figures 5 and 6 for perpendicular collection are found to be at least ⋅ for direct and indirect excitons, respectively. This is obviously significantly larger than the experimental response rate. Hence, in agreement with [5], we conclude that exciton ionization is not the limiting process in the observed photoresponse at this field strength. Rather, drift of the carriers to the contacts after ionization is the probable cause of the reduction. Generally, carrier drift in the perpendicular direction is much slower than in-plane drift. For instance, the ratio between parallel and perpendicular carrier mobility for the related material MoSe 2 can be as large as 10 3 [33]. The limiting role of carrier drift is further evidenced by the inverse square dependence of the rate on sample thickness [5]. Finally, a transfer time of 1 ps across the interface between graphene and WS 2 was recently reported [34]. Assuming a similar rate for WSe 2 , this implies an upper limit of~-10 s 12 1 for the measured photoresponse rate.

Summary
In summary, we have considered the process of field-assisted exciton ionization in multilayer TMDs. By solving an anisotropic Mott-Wannier equation we are able to extract exciton binding energies. Combined with complex scaling and resummation techniques, we subsequently find exciton Stark shifts and ionization rates as a function of external field strength. When applied to sulphides (MoS 2 and WS 2 ) and selenides (MoSe 2 and WSe 2 ), our results show that, in scaled exciton units, all of these materials behave similarly. However, after conversion into physical units, significant differences emerge. Thus, direct excitons in the tungsten compound WSe 2 are found to ionize in substantially smaller fields as a result of the smaller exciton binding energy. Similarly, for the indirect exciton, the ionization rate of MoSe 2 is the highest among the compounds considered. The linearization in k means that the above results cannot be immediately applied in the experimentally relevant range for TMDs k | |  0.5. We therefore fit the full k -dependence to a polynomial using (i) the exact known limits for k = 0 and k = 1 [16], (2)