The role of collective motion in the ultrafast charge transfer in van der Waals heterostructures

The success of van der Waals heterostructures made of graphene, metal dichalcogenides and other layered materials, hinges on the understanding of charge transfer across the interface as the foundation for new device concepts and applications. In contrast to conventional heterostructures, where a strong interfacial coupling is essential to charge transfer, recent experimental findings indicate that van der Waals heterostructues can exhibit ultrafast charge transfer despite the weak binding of these heterostructures. Here we find, using time-dependent density functional theory molecular dynamics, that the collective motion of excitons at the interface leads to plasma oscillations associated with optical excitation. By constructing a simple model of the van der Waals heterostructure, we show that there exists an unexpected criticality of the oscillations, yielding rapid charge transfer across the interface. Application to the MoS2/WS2 heterostructure yields good agreement with experiments, indicating near complete charge transfer within a timescale of 100 fs.

C harge transfer is a fundamental process governing numerous biological, chemical and physical systems. Electrons or holes transferred to different regions or functional groups can trigger chemical reactions and are integral to phenomena such as photosynthesis and water splitting 1,2 . In heterostuctures, which are a basic platform for solid-state applications such as solar cells, light-emitting diodes and highmobility transistors [3][4][5] , charge transfer controls electron-hole recombination and hence plays a crucial role in device operation 6,7 . While the structure of traditional (covalent) heterostructures are limited by the stringent requirement of lattice matching between constituent materials to avoid defect formation and large strain at the interface 8 , the severity of such issues are greatly diminished in the emerging van der Waals (vdW) heterostructures consisting of two-dimensional atomic crystals [9][10][11] . The vdW heterostructures are very promising as their components can be chosen from a variety of systems, such as graphene, hexagonal boron nitride, transition metal dichalcogenides (TMDs), phosphorene and silicene. These atomically thin materials present a set of attractive physical properties such as high thermal and electrical conductivity 12 , Klein tunnelling 13 , valleytronics [14][15][16] , strong absorption of light 17 and strong photoluminescence [17][18][19] . For vdW interfaces, one may expect that the component materials are relatively isolated from each other maintaining their individual properties because the interfacial vdW spacing is considerably wider than a covalent or ionic chemical bond, and hence electronic coupling between the layers is expected to be weak. This presumed weak coupling has motivated recent proposals for using such materials for photovoltaic and diode-like current rectification applications [20][21][22][23] . This makes the recent findings of photoluminescence and femtosecond pump-probe spectroscopy experiments 20,24 that the charge transfer in MoS 2 /WS 2 is an ultrafast phenomena (within 100 fs) quite surprising. Similarly, fast charge transfer has also been observed for the MoS 2 /MoSe 2 heterostructure 22 , suggesting that in these new heterostructures, the physics of charge transfer may fundamentally depart from traditional pictures of non-coherent charge transfer [25][26][27] . Developing a basic understanding of the mechanisms of charge transfer across the vdW hetostructure is not only of fundamental interest but essential for the successful application of these materials.
In this work, time-dependent density functional theory (TDDFT) is combined with molecular dynamics (MD) to understand the physics governing the ultrafast hole transfer from MoS 2 to WS 2 . We find that significant hole transfer takes place within a timescale of 100 fs of the initial excitation in agreement with experiments 20,24 . Although conventional wisdom suggests that throughout this process accumulated holes would lead to a dipole that would tend to hamper further hole transfer-our analysis of quantum states dynamics during the transfer process suggests precisely the opposite. Our study reveals a strong dynamic coupling, whereby charge transfer across the heterostructure significantly increases the coupling between the hole states in MoS 2 and WS 2 . It follows that after an initial induction time of relatively slow charge transfer, the charged dipole layer causes the rapid transfer of holes across the heterostructure despite the fact that the internal electric field opposes the direction of hole transfer. A model for the hole transfer shows that this mechanism is critically dependent on the magnitude of the dipole transition matrix element, which determines the coupling across the heterostructure. As a result, the existence of ultrafast charge transfer is sensitive to the atomistic details of the interface, in particular the heterostructure stacking. While the energetically favourable 2H stacking (in which the S atoms of WS 2 lie directly above Mo) is found to exhibit ultrafast charge transfer, dipole coupling in 3R and AB 0 stacking configurations 28,29 (in which the W atom of WS 2 lies directly above S or Mo from MoS 2 , respectively) is found to remain below the minimum threshold, resulting in much slower charge transfer.

Results
Real-time TDDFT-MD. To investigate the photo-induced charge transfer dynamics at the MoS 2 /WS 2 interface, we built a structural model and setup similar to the experimental situation 20,24 , see   Fig. 1a, to choose an initial occupation of electrons associated with the optical excitation, and simultaneously iterate the combined electron and ionic dynamics using time domain TDDFT-MD (see the Methods section). The appropriate initial excitation is determined through examination of the band structure shown in Fig. 1b. We find that the gap is indirect with the valence band maximum (VBM), in which the wave function is delocalized to both WS 2 and MoS 2 , at G and the conduction band minimum (CBM), in which the wave function is localized to MoS 2 , at the K-point of the Brillouin zone, consistent with previous findings 17, 30 . In the vicinity of the K-point, shown in the blue rectangle and enlarged in Fig. 1c, a type-II alignment develops where the local VBM, which we denote |WS 2 i presents a charge density almost entirely localized to WS 2 . Conversely, the state immediately below the VBM, the (VBM-1) state, is almost entirely localized to MoS 2 and is denoted by |MoS 2 i in Fig. 1c. As optical absorption via the indirect transition (VBM-CBM) or involving charge transfer across the interface (VBM(at K)-CBM) is negligible 24,30 , the predominate optical absorption is associated with the direct excitation of electrons from the (VBM-1) state to the CBM in the vicinity of K, as indicated by the red arrow in Fig. 1c. The subsequent hole transfer observed in the experiment 20,24 results from the type-II alignment, whereby the holes can substantially reduce their energy by transferring from |MoS 2 i to |WS 2 i, as indicated by the blue arrow in Fig. 1c. Hence, in the TDDFT-MD simulation, two electrons were excited from the |MoS 2 i state to the CBM and our subsequent analysis focuses on the hole dynamics. Figure 2a-c shows the time evolution of the hole (h) occupation r h , calculated by projecting the time-evolved hole wave function f(t) onto |MoS 2 i and |WS 2 i, that is, r h (MoS 2 ) ¼ |hMoS 2 |f(t)i| 2 and r h (WS 2 ) ¼ |hWS 2 |f(t)i| 2 , for the case of 2H stacking at the three different ionic temperatures, T ¼ 0, 77 and 300 K, respectively. Let us first consider the case of 77 K. Figure 2b shows that r h oscillates by periodically filling and emptying |MoS 2 i and |WS 2 i states at a period of approximately 30 fs, albeit with a damped amplitude. The decay of the r h (MoS 2 ) is significant, indicating that the average hole position is transitioning in time, from being primarily localized on |MoS 2 i to reside on either side of the interface, with nearly equal probability. As the energy of the hole is lower in |WS 2 i than in |MoS 2 i, as shown in Fig. 1c, this change in expectation value for the hole position is associated with an energy relaxation process-wherein excitation energy of the hole is being dissipated by phonons throughout the  (2), where the parameters used are indicated in Table 1. The model well produces the oscillations seen in the transition. The same effect, although even more pronounced, is also observed at 300 K, as shown in Fig. 2c. To investigate which phonons are associated with the decay process, we determine the energy of each normal mode of the system as a function of time, as shown for several time snapshots throughout the 77 K simulation in Fig. 3a-d. Here it can be seen that the A 1g Raman active modes and phonons associated with the longitudinal acoustic branch with non-zero wavevectors pick up substantial energy in the decay process. In particular, the strong coupling of the two A 1g phonon modes, involving out-of-plane vibrations of S atoms shown in Fig. 3e, is due to the large associated modification of the distance between the vdW layers. It follows that these modes present a large response to a transverse electric field associated with charge transfer as evidenced by the substantial peaks of the A 1g modes in the Raman spectra of MoS 2 and WS 2 (refs 31,32). While the longitudinal acoustic phonons also show a marked increase in energy at approximately 120 fs, we note that this may be due to subsequent phonon-phonon relaxation of the higher energy A 1g modes.
To better understand the hole dynamics and the role of phonons in the charge transfer process, we decouple the ionic motion from the hole dynamics by performing a TDDFT simulation where the ionic positions are frozen throughout the simulation, as shown in Fig. 2a. The large charge sloshing and the lack of amplitude dampening here are further evidence that it is indeed the electron-phonon coupling that is responsible for the dissipation. In addition, it confirms that the observed hole oscillation is not the result of ionic motion. Interestingly, the oscillation is not sinusoidal as one would expect from a superposition of the two states, |MoS 2 i and |WS 2 i, but rather it can be characterized by broad mesas with sharp peaks. This indicates that the coupling between |MoS 2 i and |WS 2 i, which governs the characteristic frequency of oscillation, gains time dependence throughout the simulation; initially being quite small and becoming more substantial near the peaks. Furthermore this dynamic coupling is purely an electronic phenomenon and is not associated with the ionic vibration, since these were kept fixed throughout this simulation. As such, the mechanism for coherent charge oscillations here are distinct from those proposed for organic photovoltaic blend 33 and artificial light-harvesting system 34 , wherein the charge oscillations are associated with the frequency of a vibrational mode of the system. To understand the mechanism giving rise to the nonlinear charge oscillations, timedependent coupling and ultimately the ultrafast charge transfer in the vdW heterostructure, we now further examine the electron dynamics of the system. origin of this oscillating interfacial dipole is hole transfer from |MoS 2 i to |WS 2 i, here it should be noted the importance of the collective behaviour of the holes during transport. To mimic the coherent excitation associated with a pump laser setup, 20 the initial conditions of the TDDFT-MD simulation contain two holes with identical characteristics. Throughout the course of simulation, they maintain identical characteristics with both single-particle hole wavefunctions evolving identically in time. Furthermore, the periodic boundary conditions represent an entire sheet of charges, and as a result charge transfer from MoS 2 to WS 2 is associated with the collective motion of holes across the interface with a resulting macroscopic dipole similar to the parallel plate capacitor schematically shown in Fig. 4b. As this macroscopic dipole is responsible for modulating the hole dynamics, tests indicate that for slightly different initial hole character, the results do not alter qualitatively. While the E-field produced in this manner opposes hole transfer from MoS 2 to WS 2 , it significantly enhances the coupling between the two states and hence alters the characteristic timescale of the oscillations. This phenomenon can be understood in the context of an elementary model where a single hole exists in a two-state system (either in MoS 2 or WS 2 ) in the presence of an external electric field representing the hole's interaction with the dipole generated by the other holes in the system. Assuming the geometry shown in Fig. 4b, the Hamiltonian associated with the uniform E-field in the z-direction, perpendicular to the interface, where E 0 is the initial magnitude of the z-component of the electric field associated with charge rearrangement due to the optical excitation and s t ½ E is the time-dependent field, which is dynamically generated due to hole transport from MoS 2 to WS 2 . Using the ground-state eigenfunctions |WS 2 i and |MoSi as basis and choosing the origin between the two slabs, the Hamiltonian of the two-state system in the external field becomes, where E WS2 and E MoS2 are the hole energies of the ground-state eigenfunctions of the heterostructure, d is the effective distance between the two layered materials, seen here as two plates of the capacitor, d ¼ (hWS 2 |z|WS 2 i À hMoS 2 |z|MoS 2 i), and M z is the magnitude of the z-component of the dipole transition matrix element between the heterostructures, hWS 2 |z|MoS 2 i. Here it can be seen that the electric field contributes to both the diagonal and off-diagonal matrix elements. While the diagonal ones reduce the energy difference between the two levels that drive hole transfer, the off-diagonal elements increase the dipole coupling between the MoS 2 and WS 2 states. In principle, all of the parameters in equation (2) can be calculated with first principles methods, however, as the 2 Â 2 model Hamiltonian is significantly reduced from the full Hamiltonian of the actual system, several parameters are obtained by fitting to the TDDFT-MD results.
The good agreement between the model Hamiltonian with the TDDFT-MD results using fitted values for d, E and E 0 establishes the validity of the approach, as shown in Fig. 2d. Notably, instead of a sinusoidal oscillation of r h , the model correctly predicts the sharp peaks and broad mesas observed in the TDDFT-MD simulations. The ability for such a simple model to capture the nonlinear behaviour present in the TDDFT-MD results and the closeness of the fit parameters to those calculated via first principles (Table 1) validate the physical arguments used to derive it.
Investigating how the model describes the physics of charge transfer as a function of M z , we find that criticality exists and that the hole dynamics undergoes a qualitative change depending on the magnitude of the dipole-induced coupling. Keeping the initial electric field associated with the excitation fixed at 6.9 mV Å À 1 , the maximum amount of dynamic hole transfer as a function of the magnitude of the dipole coupling M z is shown in Fig. 5a. For small coupling, this serves to mix the two states leading to small sinusoidal oscillations associated with being in a linear superposition of states. However, after the coupling reaches a critical value, at approximately 1.05 Å, the magnitude of this oscillation  While the qualitative features of the charge oscillations are well reproduced for the parameter free model, the quantitative agreement between Fig. 2a,d are obtained by using the fitted parameter values in the table above. Note that the values for M z and EMoS 2 À EWS 2 À E0d are not fit but instead taken directly from the ground-state wavefunctions and the single particle energies of the TDDFT-MD simulation, respectively. The calculated value of E=E0 shown is the average of those values obtained from bulk MoS 2 and bulk WS 2 of 7.55 and 6.92, respectively. The largest deviation between the fit and calculated values is obtained for the effective distance, d ¼ hWS2 j z j WS 2 i À hMoS2 j z j MoS 2 i, between the slabs. Here we note that the Coulombic attraction between the electron and hole in the dynamics may reduce this distance as observed by the fitted value being 25% shorter than that obtained from the ground-state wavefunctions. changes discontinuously from below 5% to above 75% of the hole population. The emergence of this discontinuity is a result of the inherent non-linearity of the Hamiltonian. Through numerical investigation, we find that the effective potential that governs the dynamics of r h is bimodal. For small values of M z , the system does not have enough energy to overcome the barrier separating the two wells and the system dynamics are trapped, yielding only small values of r h and minimal charge transfer. As M z is increased, however, the barrier separating the two wells decreases, until the critical value in which the system dynamics can overcome the barrier leading to much larger oscillations in r h . Physically, beyond this critical value, despite the initial small coupling, the collective behaviour of the holes in the system leads to feedback so that the degree of hole transfer is enough to steadily increase the coupling leading to even more charge transfer. This critical dependence on M z highlights the importance of the interlayer spacing (and vdW interaction) on the ultrafast charge-transfer mechanism. Note that the existence of substantial dipole coupling, M z ¼ hMoS 2 |z|WS 2 i, relies on the spatial overlap of the |MoS 2 i and |WS 2 i wavefunctions. Hence, the nature and strength of vdW interaction at the interface, which determines the distance between layers and hence the spatial overlap of the wavefunctions, becomes an essential part of the rapid charge transfer. Furthermore, as electron-phonon coupling is proportional to the magnitude of the dipole oscillation, the fast energy dissipation to phonons and ultimately the ultrafast charge transfer can be viewed as a direct result of this marked charge oscillation resulting from the prerequisite critical dipole coupling of the heterostructure. To extend the model to incorporate interaction with ionic motion, we include a mean-field interaction with a heat bath (see Methods section). Fitting the interaction parameter to a single receiving mode at 400 cm À 1 leads to a qualitative agreement with the finite temperature TDDFT-MD results, as shown in Fig. 2e,f. Furthermore, this fitting allows us to calculate charge dynamics as a continuous function of temperature in Fig. 5b and improve on the TDDFT-MD results by considering the quantum mechanical nature of the ions in Fig. 5c revealing that near-complete charge transfer takes place within 100 fs even as T approaches 0 K. Figure 5b,c and almost the entirety of the results presented herein refer to the case of 2H stacking between MoS 2 and WS 2 , which we have predicted to exhibit ultrafast hole transfer. In contrast, for the 3R and AB 0 stacking configurations, we find that no hole transfer takes place from MoS 2 to WS 2 over the entire simulation period of 120 fs. This result is due to the critical role played by the magnitude of the dipole transition matrix element in the onset of the charge oscillations. For both the 3R and AB 0 stackings, we find M z to be far below this critical threshold, both being approximately 10 À 5 Å. This sensitivity to stacking is also seen in the transfer integral, where the electronic coupling of the 2H case is found to be 132 meV, many orders of magnitude larger than either the 3R or AB 0 cases. This sensitivity on the atomistic details of the interface can be traced back to the C 3v symmetry of the wavefunctions associated with the VBM of the isolated layers of |MoS 2 i and |WS 2 i, which are shown in Fig. 6a,b, respectively. Regions of positive and negative wave function are highlighted by the red and blue triangles, respectively. Looking at the case of AB 0 stacking, schematically shown in Fig. 6c, it can be seen that due to the shift in the unit cell, there is very little spatial overlap in the wave function of MoS 2 (represented by the open triangles) and that of WS 2 (represented by the shaded triangles). This can also be seen in the case of 3R stacking, shown in Fig. 6d, where the small overlap that does exist tends to cancel. Conversely, for 2H stacking shown in Fig. 6e, the wavefunctions line up entirely, yielding maximal spatial overlap and an associated much stronger coupling.

Discussion
The TDDFT-MD modelling was performed for the most stable 2H stacking. However, we note that experimental studies of the excited-state hole dynamics have been carried out on far less-controlled samples than the perfect stacking since it is based on mechanical pealing and transfer 20 . The low energy of the 2H stacking suggests that despite the nominal random stacking in experiment, either the 2H stacking dominates or the large sample size guarantees that a sufficiently large proportion of isolated regions feature 2H stacking, thereby ensuring fast charge transfer.
The model describing charge transfer in vertical MoS 2 /WS 2 heterostructure is general and should be applicable to any weakly bound heterostructures where coupling is dominated by the induced electric field. Furthermore, as charge sloshing couples strongly to G-point phonons (A 1g ) associated with modulation of the vdW gap region, ultrafast charge transport across vdW heterostructures may be the rule rather than the exception. As the initiation of this mechanism for transport is critically dependent on the dipole coupling matrix element, this relationship deserves further study. In particular, our findings prompt a number of fundamental questions: What are the root causes of the discontinuity? What are the dynamics at the critical point? As the collective motion of holes is responsible for the criticality, can this mechanism be interpreted as the instability of an excitonic plasma at the interface? While the timescale of the predicted charge oscillations may be difficult to resolve with current experimental techniques, the mechanism for ultrafast charge transport presented here and the critical role of M z should be confirmed either by careful preparation of heterostructures with different stacking configurations or through an assortment of vdW heterostructures that feature different intrinsic M z . Similar to the discontinuity reproduced by the new model, it is expected that vdW heterostructures with intrinsic M z coupling larger than M crit will exhibit ultrafast charge transport while those with smaller intrinsic coupling, such as 3R and AB 0 stacked MoS 2 / WS 2 , will exhibit charge transport on a timescale several order of magnitudes larger. In summary, a combined ab initio TDDFT-MD simulation and analytic modelling reveals the nature of hole transfer as a result of carriers dynamics at the MoS 2 /WS 2 vdW heterostructures. Despite the weak interlayer vdW binding, the collective motion of the holes is found to lead to strong dynamic coupling due to the transverse electric field associated with charge transfer. In contrast to classical systems, the build-up of an interfacial dipole can serve as an enhancement factor for such an ultrafast transfer in quantum coherent systems. This ultrafast transfer is well described by a simple model Hamiltonian in which, provided the existence of a critical dipole transition matrix element, the collective motion leads to dramatic non-linear charge oscillations that couple strongly to the phonons at the G-point. We anticipate that these findings will trigger further fundamental studies and new device concepts over a broad range of two-dimensional layered materials and systems.

Methods
Real-time TDDFT-MD. The electron and ion dynamics are simulated within the TDDFT formalism 35 coupled with MD, within the local density approximation (LDA) as implemented in the SIESTA code [36][37][38] . We employ Troullier-Martins norm-conserving pseudopotentials and Ceperley-Alder 39 exchange-correlation functional. Local basis set with double-z polarized orbitals is used. The real-space grid is equivalent to a plane-wave cutoff energy of 350 Ry. A 3 ffiffi ffi 3 p Â 3 ffiffi ffi 3 p supercell containing 162 atoms with G-point for the Brillouin zone integration is used. Note that the K-point of primitive cell is folded to the G point in this supercell. In contrast to the well-known underbinding character in the PBE functional, the LDA reproduces the experimental structural parameters quite well for the vdW systems, as shown in Table 2. Although vdW forces are not explicitly included in the LDA formalism, LDA is known to yield a good description of forces between layers in TMDs as indicated by the quantitative agreement of the layer-to-layer vibration modes in these systems, as compared with experimental Raman data 31,[40][41][42][43] . Furthermore, LDA predicts the preference of the 2H stacking, similar to the vdW density functionals (which are built on top of semi-classical GGA functionals such as PBE) 29 . We consider 3R, AB 0 and 2H stacking patterns for the MoS 2 /WS 2 heterostructures 28,29 , but we focus our discussion on the 2H stacking patterns in this study. Testing for the case of 2H stacking indicates that TDDFT-MD calculations using the LDA and PBE þ vdW functionals yield quite similar results for the hole transfer dynamics across the heterostructure. We use optimized lattice constants 3.16 Å and all the atoms are fully relaxed until the residual forces are o0.025 eV Å À 1 . In the dynamics calculations, we use NVE ensemble, a time step of 24.19 attoseconds, and the Ehrenfest approximation for ionic motion. Initial ion positions and velocities are obtained from ground-state MD simulations at the target ionic temperature.
Phonon decomposition. The eigenvectors, v n ai , and eigenvalues, o n , of the normal modes of the system were determined by diagonalizing the dynamical matrix as implemented in SIESTA. Cartesian coordinates of the atomic trajectories during the TDDFT-MD simulation are recast in terms of the normal mode coordinates, q n , by projecting the ionic displacements from equilibrium onto the eigenvectors of the dynamical matrix. In this way, the energy of each normal mode is determined throughout the simulation, according to Model calculation. The model calculations are performed by numerically integrating the Liouville equation, r Á ¼ H;r ½ i' , to obtain the density matrix as a function of time. To incorporate the dissipation to the phonon modes observed in the TDDFT-MD calculations at elevated temperatures in the model, the interaction of the system with a heat bath is considered in an average way. In every step in the time evolution of r, we enforce a Kraus transformation 44 where r t þ dt ð Þ¼ 1 À P ð Þr 11 r 12 ffiffiffiffiffiffiffiffiffiffi ffi 1 À P p r 21 ffiffiffiffiffiffiffiffiffiffi ffi 1 À P p r 22 þ Pr 11 ; ð3Þ with P being the probability that the event of phonon emission with the hole transitioning to a lower energy state (that is, |MoS 2 i|ni bath -|WS 2 in þ 1i bath ) takes place in the time dt. Note here that the bath is Markovian and only dissipates energy, we do not consider events which raise the energy of the hole. From Fermi's golden rule, the probability of such a transition is directly proportional to the square of the amplitude of the receiving mode, allowing for the relation of the transition probabilities for different temperatures. Here, however, it is important to note that while the TDDFT-MD calculation is not on the Born-Oppenheimer surface, the ions are still treated as classical objects, and hence the energies of the normal modes in the TDDFT-MD calculations follow a classical equipartition theorem with E ¼ kT, instead of Bose-Einstein statistics. Hence, for the case of classical ions, assuming a single phonon decay channel associated with a frequency o, we can simultaneously fit both the dynamics at 77 K and 300 K with a single coupling parameter, g, where the temperature-dependent transition probability becomes While the amplitude distributions of the normal modes associated with the fully quantum mechanical system are not well represented in the TDDFT-MD, the underlying coupling is well described by the Ehrenfest approximation and hence it is expected that the fitted coupling parameter, g, associated with the decay process is still valid. Hence, we can improve the TDDFT-MD simulations to more adequately represent the fully quantum mechanical case by replacing the amplitudes in equation (4) by those correctly determined from the Bose-Einstein distribution yielding, By correcting for this shortcoming, we can go beyond the TDDFT-MD to more accurately calculate charge transfer dynamics at low T, as shown in Fig. 5c. In this case the temperature dependence is much less pronounced and we see nearly complete charge transfer in a 100 fs timescale even as we approach T ¼ 0 K. This relative insensitivity to temperature is a result of the decay process dominated by high-frequency phonons, which even at room temperature have only very modest phonon occupation (for instance at 300 K, nB0.16). By directly comparing Fig. 5b,c, we see that only as we approach room temperature do the classical and quantum pictures of ionic motion yield correspondence. The experimental results are taken from reference (46). The theoretical results are this work and are calculated within the local density approximation.