Thermalization of photoexcited carriers in two-dimensional transition metal dichalcogenides \\ and internal quantum efficiency of van der Waals heterostructures

Van der Waals semiconductor heterostructures could be a platform to harness hot photoexcited carriers in the next generation of optoelectronic and photovoltaic devices. The internal quantum efficiency of hot-carrier devices is determined by the relation between photocarrier extraction and thermalization rates. Using \textit{ab-initio} methods we show that the photocarrier thermalization time in single-layer transition metal dichalcogenides strongly depends on the peculiarities of the phonon spectrum and the electronic spin-orbit coupling. In detail, the lifted spin degeneracy in the valence band suppresses the hole scattering on acoustic phonons, slowing down the thermalization of holes by one order of magnitude as compared to electrons. Moreover, the hole thermalization time behaves differently in MoS$_2$ and WSe$_2$ because spin-orbit interactions differ in these seemingly similar materials. We predict that the internal quantum efficiency of a tunneling van der Waals semiconductor heterostructure depends qualitatively on whether MoS$_2$ or WSe$_2$ is used.


I. INTRODUCTION
Utilizing high-energy carriers in photovoltaic devices could improve light-to-energy conversion efficiencies [1]. Despite recent progress with hot-carrier generation and injection in plasmonic nanostructures [2], the conventional semiconductor solar cells still demonstrate a superior efficiency without need for nanoscale fabrication. Alternative approaches to the problem either facilitate hot-electron transfer, e.g. from a chemically modified surface of lead selenide to titanium oxide [3], or extend photocarrier lifetimes, e.g. in some perovskites [4,5]. Van der Waals (vdW) semiconductor heterostructures formed from atomically thin two-dimensional (2D) crystals [6] might represent a suitable platform to take advantage of both phenomena, i.e. ultrafast photocarrier extraction and slow photocarrier thermalization.
The possibility to assemble vdW heterostructures layer by layer allows to tune electron transfer across interfaces in the out-of-plane dimension. This interlayer charge transfer directly competes with intralayer photocarrier thermalization. The scheme has first been realized in a graphene-boronnitride-graphene vdW heterostructure [7], where the interlayer tunneling time ranges from 1 fs to 1 ps depending on bias voltage. Since the photocarrier thermalization time in graphene spans the interval between 10 fs and 10 ps depending on carrier concentration and excitation energy [8][9][10][11][12][13], there is a parameter range within which high-energy photocarrier extraction is feasible. Substituting graphene by a 2D transition metal dichalcogenide (TMDC) in a stack [14][15][16][17] opens new perspectives thanks to the possibility to create a diode configuration [18] and to realize stronger light-matter interactions [19]. There are already several reports on the fabrication of TMDC-based vdW heterostructures and their optoelectronic properties [20][21][22]. The high-energy photocarriers may be filtered by means of a boron-nitride layer, constituting a barrier for thermalized electrons and holes with low energy [21,22]. Since the photocarrier thermalization time is sensititive to the electron and phonon spectra of each semiconductor, it is important to compare this thermalization time with the interlayer tunneling time: The high-energy photocarrier transport may FIG. 1. Photocarrier excitation, thermalization, and tunneling in a semiconductor-dielectric-semiconductor vdW heterostructure. A weak bias (not shown) is applied to create a current across the layers. The internal quantum efficiency is determined by the ratio between thermalization and tunneling rates. Here, we consider transport through the valence bands because of the longer thermalization time for photoexcited holes. not be feasible, if the former is too short. The relevant processes are sketched in Fig. 1. Note that we do not distinguish between the terms "hot" or "high-energy" carriers, but use them interchangably.
Free photoexcited carrier thermalization in 2D TMDCs, arising for excitations above the band gap, has been studied much less comprehensively [23,30,[34][35][36]40]. Ab-initio studies have shown that the carrier-carrier scattering is efficient far away from the band edges due to a quadratic increase of the available phase space [41,42]. For a clean sample with no defects, the interaction with phonons thus represents the most important nonradiative channel for thermalization and cooling of photoexcited carriers in the vicintiy of the valence and conduction band edges [25,40,43] and requires a deeper understanding.
In what follows, we use an approach already applied successfully to bulk and 2D materials [13,42,44,45]. For instance we have demonstrated in our recent work that phonon emission by photocarriers in graphene is strongly suppressed at low excess energies of about 100 meV [13]. This phenomenon occurs because of high optical phonon frequencies due to strong carbon-carbon bonding and has been observed experimentally as a thermalization bottleneck [10]. In the present work we study the influence of spin-orbit coupling (SOC) on photocarrier thermalization in single-layer TMDCs and additionally relate intralayer thermalization rates to interlayer tunneling rates in vdW heterostructures, important for future optoelectronic devices with improved internal quantum efficiency (IQE). For this purpose we make use of density functional theory (DFT) and density functional perturbation theory (DFPT) to determine the electronic and phononic spectra of two 2D TMDC representatives, namely MoS 2 and WSe 2 . Furthermore, we compute band-and momentumdependent scattering rates τ −1 nk for a given electronic band n and wave vector k. Applying the Boltzmann equation in the relaxation time approximation (RTA) we extract the total thermalization time τ th as a function of excess energy and temperature. In our parameter-free ab-initio modeling, we take relativistic SOC and all the optical and acoustical phonon branches in the first Brillouin zone (BZ) into account to provide a reliable description of electron-phonon scattering events. We show that photocarrier thermalization is slowed down by SOCinduced band splitting near valence band maxima and conduction band minima. Similar to interlayer transport, the thermalization occurs faster for photocarriers excited farther away from the band edges. We find, however, that the hole thermalization rate increases with excess energy in MoS 2 much more rapidly than in WSe 2 . We therefore expect that the IQE of the tunneling device shown in Fig. 1 will strongly depend on whether MoS 2 or WSe 2 is employed.

A. Ab-initio theory for electronic and phononic properties
We determine electronic and phononic properties of MoS 2 and WSe 2 monolayers using DFT within the local density approximation (LDA) as implemented in QUANTUM ESPRESSO [46]. Six outermost electrons of each transition metal (Mo, W) and chalcogen (S, Se) are treated explicitly as valence electrons, while the remaining core electrons are included through norm-conserving Troullier-Martins pseudopotentials with relativistic corrections [47]. We employ a plane-wave basis set with a kinetic energy cutoff of 70 Ry and a charge density cutoff of 280 Ry. Unit cells of MoS 2 and WSe 2 monolayers are optimized with the help of the Broyden-Fletcher-Goldfarb-Shanno algorithm, neglecting SOC until the net force on atoms is less than 10 −6 Ry/a.u. and total energy changes are below 10 −8 Ry. The monolayers are placed in a cell with a vacuum that separates periodic images by 18Å to avoid artificial interactions in the out-of-plane direction. After geometry opimization we calculate the ground state density with and without SOC, sampling the BZ with a 45 × 45 × 1 Γ-centered k grid. We evaluate the phonon dispersion of MoS 2 and WSe 2 monolayers through DFPT [48] with and without SOC, employing a 12 × 12 × 1 q grid to obtain phonon dynamical matrices. Next we construct localized Wannier functions from planewave eigenfunctions. By using the Wannier function interpolation scheme, we obtain electronic eigenenergies, dynamical matrices and electron-phonon matrix elements on desired grids in the BZ [49][50][51].
After structural relaxation the in-plane lattice constants |a| for MoS 2 and WSe 2 monolayers turn out to be 3.16Å and 3.27Å, respectively, which is in good agreement with experiment and previous ab-initio calculations [52,53]. The electronic band structure is displayed in Fig. 2(a,b). We find a direct quasiparticle band gap of 1.74 eV at the K point in the BZ for MoS 2 , while WSe 2 exhibits an indirect gap of 1.44 eV with the valence band maximum at the K point and the conduction band minimum located on the Γ-K direction. As DFT is known to underestimate band gaps, we apply a rigid shift to the unoccoupied conduction bands to arrive at values of 2.82 and 2.42 eV for MoS 2 and WSe 2 , respectively, which are consistent with the GW approximation [54,55]. This shift modifies the band gap values but does not affect the time evolution of the hot carriers to be discussed later. This assumption is valid because the electronic band gap is much larger than the highest phonon energy, and electrons and holes thus ther-malize independently. Fig. 2(a,b) furthermore shows that the Wannier-function-interpolated band structure matches exactly with the band structure determined with the plane wave basis, confirming the high quality of the localized basis set. Due to three atoms in the unit cell MoS 2 and WSe 2 feature three acoustical and six optical modes of vibration, as it is visible from the phononic band structure in Fig. 2(c,d). The highest frequencies of 60 and 37 meV for MoS 2 and WSe 2 monolayers occur at the Γ point, respectively, and these values agree with previously reported ones [56,57]. An energy gap separates acoustical and optical branches, which will be exploited in the analysis of thermalization times further below.
Having determined electronic and phononic band structures, we evaluate the electronic self-energy Σ nk (T ) in the lowest order of the electron-phonon interaction for band n and wave vector k at temperature T using the EPW code [51]. It is expressed as Here, ω pq is the energy of the phonon of branch p at wave vector q, ε F = 0 is the Fermi energy, f is the Bose function, Ω BZ is the volume of the BZ, and η = 20 meV is a small broadening parameter. The electron-phonon matrix elements are defined as [51] describing transitions between the Kohn-Sham states |nk and |mk + q mediated by the phonon pq. Here, ∂ pq V is the derivative of the self-consistent Kohn-Sham potential with respect to displacements of nuclei along the phonon mode pq. The first term in the brackets of Eq. (1) arises from the absorption of phonons and the second one from their emission.
To compute Σ nk (T ), we interpolate electronic and vibrational states on fine 300 × 300 × 1 k and q grids, which we find sufficient to accurately map the first BZ and to converge the integral over q.
Our model makes the following assumptions: (i) Changes induced by the electron-phonon coupling in the electronic wavefunctions and phonon dynamical matrices are small and can be neglected [51], (ii) similarly renormalization of phonon frequencies due to anharmonic effects is ignored [58], (iii) electron and phonon baths are at the same temperature T all the time. The electron-phonon scattering time for a carrier in the state |nk is then given by Conversely, Im[Σ nk (T )] is proportional to the electronphonon scattering rate τ nk (T ) −1 .

B. Time-evolution of excited charge carriers
We determine the time evolution of the electronic occupation f nk (t, T ) using the Boltzmann equation in the RTA which has the solution Here, f th nk (T ) is the thermalized Fermi distribution at t → ∞, where the quasi-Fermi level is chosen to reproduce the number of carriers initially present. The description through the Boltzmann equation (4) is valid only, if the number of carriers excited at t = 0 represents a weak perturbation. RTA relaxes an excited charge carrier directly to the thermalized state, omitting all intermediate relaxation steps between the initial nonequilibrium and final relaxed occupation.
We start with a hot-carrier occupation f nk (0, T ) that consists of the sum of a Fermi-Dirac distribution f (0) nk (T ) at temperature T and a Gaussian peak centered at energy +ζ or −ζ for electrons or holes, respectively, see also Fig. 1, as The Gaussian distribution uses a small energy broadening σ = 8.47 meV. Since valence and conduction bands feature an energy-dependent density of states (DOS), λ e and λ h are adjusted such that the same number of photoexcited carriers is present at any excess energy ζ.
We define the thermalization time τ th as with the energy-, time-and temperature-dependent population In our numerical calculations, we approximate the delta function in Eq. (8) by a narrow Gaussian with a width of 20 meV. We focus on electronic excitations far enough above the quasi-Fermi levels such that the finite width does not affect τ th .

C. Interlayer charge transport
Having the thermalization time at hand, we need a reference timescale to see whether high-energy photoexcited carriers can contribute to interlayer electron transport in vdW heterostructures, as sketched in Fig. 1. An accurate calculation of the interlayer transport time is a challenging task, because the carrier motion across interfaces is subject to multiple uncontrolled effects, e.g. interfacial roughness and impurity scattering. These effects may vary even within the same batch of samples, let alone the use of different 2D materials. In what follows, we employ a concept based on the transmission coefficient for carriers tunneling through a barrier and the uncertainty relation between the photocarrier excess energy and lifetime [7]. The advantage of the approach is that we can straightforwardly relate the transmission probability of the photoexcited carriers to a measurable quantity, namely the IQE.
Our starting point is the well-known transmission probability of carriers through a rectangular barrier of width d and height ∆ 0 counted from the respective band edge, see Fig. 1. We assume the tunneling regime [7] so that the photocarrier excess energy is always below the barrier ζ < ∆ 0 . Note that according to the conventions used in this work ζ and ∆ 0 are both positive, even for holes. We assume that the in-plane photocarrier momentum is not conserved due to interfacal disorder, rendering the problem effectively one-dimensional. We must, however, change the physical meaning of the incident wave vector k z . It does not describe propagating waves anymore but is related to the out-of-plane momentum uncertainty ∆p z / , which is of the order of the size of the first Brillouin zone for 2D conductors [59]. The approximations can formally be summarized as [7,60] In the last line, we have used the relations κ = 2m(∆ 0 − ζ) and 2 k 2 z /m ≈ v z ∆p z . Since the structure is aperiodic along the out-of-plane direction, the effective mass m of the quasiparticle moving across the interface equals the free electron mass. Importantly, there is a relation between the quasiparticle velocity, momentum uncertainty, and lifetime [59,60]: v z ∆p z ≈ 1/τ th . Note that we have assumed here that the lifetime is determined by the thermalization time. Hence, the photocarrier transmission probability can be expressed as T (ζ) ≈ τ th /τ tun , where Equation (13) is somewhat similar to the transport time formula for a triangular barrier, employed in Ref. [7] to describe electron tunneling in the Fowler-Nordheim regime. In our case, we assume a low bias and a thin barrier such that the voltage does not appear in τ tun explicitly, but it may influence ∆ 0 . As we need τ tun solely for comparison with τ th by the order of magnitude, we estimate it roughly for ∆ 0 ≈ 1 eV and d ≈ 1 nm. The interlayer transport time then ranges between 100 fs to 10 ps, depending on the excess energy. To gain a substantial photocurrent, the interlayer transport must not be too slow as compared with thermalization. From the experimental point of view, the ratio τ th /τ tun is nothing else but the IQE in the tunneling limit τ th τ tun . In the following section we will further analyze τ th , τ tun and the ratio τ th /τ tun .

III. RESULTS
We start by discussing the imaginary part of the self-energy, which is proportional to the electron-phonon scattering rate τ −1 nk . Figure 3 shows Im[Σ nk ] for MoS 2 and WSe 2 with and without SOC as a function of energy at different temperatures. The energy dependence of Im[Σ nk ] follows that of the electronic DOS. This can be understood from Eq. (1), where the self-energy for state |nk is determined by a sum over all electronic states |mk + q . Considering that phonon energies are limited to below 60 meV for both materials, the sum es- sentially constitutes an integral over electronic states in the vicinity of nk , which is proportional to the electronic DOS. SOC reshapes the electronic band structure by splitting spindegenerate states. On the rather large energy scales shown in Fig. 3, comprising excess energies of up to 3 eV, we find that this can increase the imaginary part of the self-energy as compared to the case without SOC.
More relevant to us are the effects of SOC around the band edges. Comparing scattering times with and without SOC in Fig. 4, we observe that the τ nk are strongly increased by SOC for holes of both MoS 2 and WSe 2 and electrons of WSe 2 . The electron states of MoS 2 , on the other hand, are basically unaffected. The large increase of τ nk correlates with a large spinorbit splitting of the corresponding valence and conduction band states in Fig. 2. Let us point out that similar calculations of τ nk have already been presented by Ciccarino et al. [41]. The authors have assigned the increased scattering times to the suppression of intraband scattering and spin-valley locking. In the following we extend that study by calculating thermalization times with the Boltzmann equation, by distinuishing the roles of acoustical and optical phonons in the scattering processes, and by comparing thermalization to tunneling times in vdW heterostructures.
Typical densities of carriers excited in pump-probe experiments lie between 10 11 to 10 13 cm −2 [23,28,61]. We adjust our free parameters λ e , λ h of Eq. (6) such that we generate non-equilibrium populations corresponding to densities of 7 × 10 12 cm −2 carriers in the valence and conduction bands and use the same number for all the initial state preparations throughout the paper. The time evolution of the occupation is then calculated using Eq. (5) with the thermalized Fermi-Dirac distribution f th nk centered at a quasi-Fermi level such that the density of 7 × 10 12 cm −2 carriers is also present after thermalization. This assumption of a conserved number of carriers is valid, because radiative recombination times of electronhole pairs in MoS 2 and WSe 2 are on the order of a few ps [25,28,39]. They are thus much longer than the thermalization times that we will report in the following. Fig. 5(a) shows the evolution of the hot-carrier population for MoS 2 at T = 0 K without SOC. In this case thermalization times of electrons and holes are comparable. Upon introducing the SOC, it can be seen in Fig. 5(b) that the thermalization time of holes is increased, whereas it remains almost unchanged for electrons. WSe 2 follows the same trend as MoS 2 , see Fig. 5(c,d). We attribute this behavior to the strong spinorbit splitting of the valence band in both MoS 2 and WSe 2 . The finite energy separation between optical and acoustical branches, observed in Fig. 2, makes it possible to distinguish them in Im[Σ nk ] = Im[Σ nk ] ac + Im[Σ nk ] op and consequently in τ −1 th = τ −1 th,ac + τ −1 th,op . In Figs. 6 and 7 we show thermalization times determined by disregarding or considering SOC, respectively. Beside the total thermalization time τ th (•) we display τ th,ac ( ) and τ th,op (×), calculated by taking into account only acoustical or optical phonons. In Fig. 6 it can be seen that thermalization of electrons and holes at low temperatures near the band edges (see the low excess energy ζ = 0.08 eV) is typically dominated by low-energy acoustical phonons. The electrons of MoS 2 constitute an exception, since their thermalization is governed by optical phonons. With increasing temperature acoustical phonons keep their dominant influence on τ th or start to dominante the electron-phonon scattering at T 100 K for the electrons of MoS 2 . Away from the band edges (see the high excess energy ζ = 0.8 eV) a rather diverse picture arises, where both acoustical and optical contributions can define thermalization. Once SOC is introduced in Fig. 7, the contribution of acoustical phonons is strongly suppressed for low ζ and low T , while that of the optical phonons remains nearly unaffected. For this reason the thermalization of both electrons and holes at low temperatures is fully goverened by optical phonons. With increasing temperature acoustical phonons play an increasingly important role and at T 100 K they define the decreasing behavior of τ th . Altogehter, the effects of SOC slow down thermalization near the band edges.
We summarize the total thermalization time τ th in Fig. 8 as a function of temperature for different excess energies. The dashed curves represent the thermalization time without SOC, and the solid curves include SOC. While SOC strongly increases the thermalization time at low T near the band edges, τ th may also decrease (see for instance ζ = 0.8 eV in Fig. 8(d)) depending upon the band structure and electronphonon couplings. For a fixed excitation energy the thermalization time generally decreases with increasing temperature, because more phonons are accessible to scatter on electrons and holes. The same trend has been reported in pumpprobe experiments performed on five-layer MoS 2 [23], where the drop of τ th has been observed for electrons at around T = 300 K. Fig. 8(a,b) suggests that it occurs around 100 K for MoS 2 in our theory. It should however be emphasized that we study perfectly crystalline, free-standing layers in vacuum. The presence of a substrate, defects and interlayer interactions thus complicates a comparison with the experimental results. A qualitative agreement is nevertheless achieved.
As visible in Fig. 8, free carriers near the band edges feature longer relaxation times as compared with those at higher FIG. 7. Same as Fig. 6 but including SOC. excess energies. The reason is the lower electronic DOS available for scattering around the band edges. Caused by multiple valleys in the conduction band, see Figs. 2 and 3, electrons thermalize faster than holes in both MoS 2 and WSe 2 . Due to the heavier elements composing WSe 2 , the SOC-related splitting at both the valence and conduction band edges is larger than in MoS 2 . This alters the thermalization time of both electrons and holes. However, electrons still thermalize one order of magnitude faster than holes around band edges. We note that our calculated thermalization times of electrons in MoS 2 for ζ = 0.05 eV amount to 63 and 36 fs at 10 and 100 K, respectively, which is of a similar magnitude as found in Ref. [36], where τ th ≈ 30 fs at T ≈ 50 K.
Since our calculations predict that holes thermalize more slowly than electrons in both MoS 2 and WSe 2 , we focus on hole transport to discuss the IQE of the tunneling device presented in Fig. 1. The heterostructure that we have in mind thus consists of two single layers of 2D TMDCs at the outside, which are separated by a 1 nm thick dielectric film that establishes a tunneling barrier ∆ 0 of 1 eV to the valence band maximum of the single layers. Figure 9 shows how the photoexcited hole thermalization time τ th compares to the tunneling time τ tun of Eq. (13) as a function of the excess energy ζ. Surprisingly, τ th follows the same trend as τ tun for WSe 2 . The ratio τ th /τ tun thus depends only weakly on ζ, as visible in the inset of Fig. 9. This suggests that the IQE of the vdW heterostructure stays rather constant with regard to increasing photocarrier excess energy. This contrasts to the case of MoS 2 , where τ th drops much faster than τ tun when ζ increases. Hence, we expect that an analogous MoS 2 -based tunneling device will demonstrate a pronounced decay in the IQE once the photocarriers are excited away from the band edges. Ultimately, our calculations trace back this effect to the larger SOC-induced valence band splitting of WSe 2 as compared to MoS 2 in these otherwise similar materials, see Fig. 2.

IV. CONCLUSIONS
We have combined DFT with many-body perturbation theory to calculate scattering times of photoexcited carriers in MoS 2 and WSe 2 monolayers arising from electron-phonon interaction. Our model takes all of the phonon branches and their dispersion into account in the scattering processes and highlights the crucial influence of SOC on the thermalization time. In the relevant region near the band edges we find that the inclusion of SOC generally increases the thermalization times by up to an order of magnitude. Our analysis assigns this effect to a suppression of acoustical phonon scattering, while thermalization by optical phonon remains basically unaffected. In both monolayers electrons thermalize almost one order of magnitude faster than holes. We have additionally estimated the tunneling time of a TMDC-based heterostructure using an analytical model. The ratio of thermalization and tunneling times decreases weakly with photocarrier excess energy in WSe 2 , whereas it drops quickly for MoS 2 when moving away from the valence band maximum. Our calculations hence suggest that tunneling devices based on WSe 2 monolayers will have a higher internal quantum efficiency than MoS 2based systems.