Introduction

Excitons, bound states formed by one electron and one hole, dominate the photophysics of a wide class of novel functional materials including transition metal dichalcogenides1,2,3,4, perovskites5,6,7, Xenes8,9,10 and related van der Waals heterostructures11,12,13. The microscopic understanding of the competing processes governing the exciton dynamics is therefore of crucial importance for technological applications in optoelectronics14,15,16, photovoltaics17,18,19 and photocatalysis20,21,22.

The non-equilibrium properties of excitons are largely influenced by the electron-phonon (e-ph) coupling g23, that has been shown to determine excitonic relaxation lifetimes24,25,26, optical linewidths1,27, energy renormalizations28, dephasing rates29,30, and diffusion dynamics31,32. The theoretical description of the exciton dynamics takes advantage from a very useful quantity called exciton-phonon (X-ph) coupling33. Its expression has been independently derived by a number of authors33,34,35,36,37 and it involves the projection of the screenede-ph coupling \({\mathbb{g}}\equiv {\varepsilon }^{-1}g\) (with ε the dielectric constant of the material) between two excitonic states. Accurate calculations of the X-ph coupling have been performed and excellent agreement between theory and experiments has been reported for phonon-assisted photoluminescence spectra34,38, exciton linewidths39,40, valley depolarization time41, and polaronic redshifts in absorption spectra42.

The X-ph coupling does not exhaust the possible interactions between excitons and lattice vibrations. Transient optical spectroscopies have recently revealed periodic modulations of the excitonic resonances43,44,45,46,47, pointing to a strong coupling between excitons and coherent phonons (X-cph), i.e. waves of atomic vibrations extending over a macroscopic spatial range. In spite of the great interest in the topic, a microscopic theory of the X-cph interaction is still lacking.

In this work we show that the X-cph coupling \({G}_{\nu }^{\lambda }\) between the lowest bright exciton λ generated after resonant pumping and a coherent phonon mode ν is

$${G}_{\nu }^{\lambda }=\sum\limits_{{{{\bf{k}}}}vc{c}^{{\prime} }}{g}_{c{c}^{{\prime} }}^{\nu }({{{\bf{k}}}}){A}_{{{{\bf{k}}}}vc}^{\lambda * }{A}_{{{{\bf{k}}}}v{c}^{{\prime} }}^{\lambda }-\sum\limits_{{{{\bf{k}}}}cv{v}^{{\prime} }}{g}_{v{v}^{{\prime} }}^{\nu }({{{\bf{k}}}}){A}_{{{{\bf{k}}}}{v}^{{\prime} }c}^{\lambda * }{A}_{{{{\bf{k}}}}vc}^{\lambda },$$
(1)

where \({A}_{{{{\bf{k}}}}vc}^{\lambda }\) is the exciton wavefunction, \({g}_{ij}^{\nu }({{{\bf{k}}}})={g}_{ij}^{\nu }({{{\bf{k}}}},{{{\bf{q}}}}=0)\) is the bare e-ph coupling for an electron with momentum k to be scattered from band i to band j by the phonon mode ν (momentum transfer q = 0), and the sum runs over all valence (v) and conduction (c) bands. We further show that for resonant pumping the X-cph coupling is responsible to change the exciton energy in time according to (up to a phase)

$$\delta {E}^{\lambda }(\tau )=n\sum\limits_{\nu }\frac{| {G}_{\nu }^{\lambda }{| }^{2}}{\hslash {\omega }_{\nu }}\times \left[\frac{{n}_{\nu }}{n}\cos {\omega }_{\nu }\tau -1\right],$$
(2)

where ων ≡ ωνq=0 is the frequency of the ν-th optical mode at the Γ point, n is the excitation density (number of conduction electrons per unit cell) and nν ≤ n are “effective” densities depending on the duration TP of the pump pulse. The inequality is saturated (nν = n) only for TP shorter than the optical periods Tν = 2π/ων (displacive excitations48,49,50).

The interaction with coherent phonons does not involve a scattering between different excitons, its main effect being a polaronic-like red-shift \(-n{\sum }_{\nu }\frac{| {G}_{\nu }^{\lambda }{| }^{2}}{\hslash {\omega }_{\nu }}\) and monochromatic oscillations or beatings (depending on the number of coupled optical modes) of the exciton energy. The polaronic-like red-shift adds up to the one due to (incoherent) phonons40, which involves the square of the X-ph coupling and it is proportional to the number of phonons rather than to the excitation density n. Notably the X-cph coupling in Eq. (1) is formally identical to the X-ph coupling evaluated at the Γ point for the diagonal scattering λ → λ33,34,35,36,37, with the crucial difference that the bare g instead of the screened \({\mathbb{g}}\) appears in the former. The difference arises from the diagrammatic origin of these couplings, the X-cph one emerging from the Ehrenfest self-energy51, see below. From Eq. (2) we infer that the amplitude of the coherent oscillations in a time-resolved optical spectrum is proportional to the excitation density and the square modulus of the X-cph coupling. This prediction is confirmed through real-time (RT) simulations of transient absorption in resonantly pumped MoS2 monolayer. Our results provide a readily applicable formula for a quantitative understanding of the coherent dynamics in excitonic materials.

Results

Coupling between excitons and coherent phonons

We consider a semiconductor hosting bound excitons, and denote by λ the lowest optically bright exciton. The exciton state is described by the coherent superposition of electron-hole pairs

$$\left\vert \lambda \right\rangle =\sum\limits_{{{{\bf{k}}}}vc}{A}_{{{{\bf{k}}}}vc}^{\lambda }{\hat{d}}_{{{{\bf{k}}}}c}^{{\dagger} }{\hat{d}}_{{{{\bf{k}}}}v}\left\vert {\Phi }_{0}\right\rangle ,$$
(3)

where the operator \({\hat{d}}_{{{{\bf{k}}}}i}^{({\dagger} )}\) annihilates (creates) an electron with momentum k in band i = {v, c}, and \(\left\vert {\Phi }_{0}\right\rangle\) denotes the ground-state (filled valence bands and empty conduction bands). The exciton wavefunction \({A}_{{{{\bf{k}}}}vc}^{\lambda }\) is normalized to unity, \({\sum }_{{{{\bf{k}}}}vc}| {A}_{{{{\bf{k}}}}vc}^{\lambda }{| }^{2}=1\), and is obtained by solving the Bethe-Salpeter equation (BSE) HAλ = EλAλ where

$${H}_{{{{\bf{k}}}}vc,{{{{\bf{k}}}}}^{{\prime} }{v}^{{\prime} }{c}^{{\prime} }}=({\epsilon }_{{{{\bf{k}}}}c}-{\epsilon }_{{{{\bf{k}}}}v}){\delta }_{{{{\bf{k}}}}{{{{\bf{k}}}}}^{{\prime} }}{\delta }_{c{c}^{{\prime} }}{\delta }_{v{v}^{{\prime} }}+{K}_{{{{\bf{k}}}}vc,{{{{\bf{k}}}}}^{{\prime} }{v}^{{\prime} }{c}^{{\prime} }}.$$
(4)

Here ϵkv and ϵkc are valence and conduction band dispersions and K is the BSE kernel in the Hartree plus statically screened exchange (HSEX) approximation.

For weak resonant pumping with the exciton energy Eλ the many-body quantum state at time t is well described by \(\left\vert \Psi (t)\right\rangle =\left\vert {\Phi }_{0}\right\rangle +\sqrt{N(t)}\,{e}^{{{{\rm{i}}}}{E}^{\lambda }t}\left\vert \lambda \right\rangle\), where N(t) is the total number of conduction electrons at time t. The weak pumping assumption implies that the excitation density n(t) = N(t)/Nk ≪ 1, Nk being the number of k-points in the first Brillouin zone. The change in the one-particle density matrix to lowest order in n is then given by

$$\delta {\rho }_{{{{\bf{k}}}}c{c}^{{\prime} }}(t)=\langle \Psi (t)| {\hat{d}}_{{{{\bf{k}}}}{c}^{{\prime} }}^{{\dagger} }{\hat{d}}_{{{{\bf{k}}}}c}| \Psi (t)\rangle =N(t)\sum\limits_{v}{A}_{{{{\bf{k}}}}v{c}^{{\prime} }}^{\lambda * }{A}_{{{{\bf{k}}}}vc}^{\lambda },$$
(5a)
$$\delta {\rho }_{{{{\bf{k}}}}v{v}^{{\prime} }}(t)=-\langle \Psi (t)| {\hat{d}}_{{{{\bf{k}}}}v}{\hat{d}}_{{{{\bf{k}}}}{v}^{{\prime} }}^{{\dagger} }| \Psi (t)\rangle =-N(t)\sum\limits_{c}{A}_{{{{\bf{k}}}}vc}^{\lambda * }{A}_{{{{\bf{k}}}}{v}^{{\prime} }c}^{\lambda }.$$
(5b)

Notice that for times subsequent to the end of the pump, i.e., t > TP, N(t) = N is independent of time52,53, and so are the density matrix elements for each k. This peculiarity arises from pumping resonantly with the lowest excitonic state.

The bare e-ph coupling governs the dynamics of coherent phonons. Let xν be the displacement of the ν-th optical mode. To lowest order in n and for times t > TP it satisfies the equation of motion (EOM)54

$$M{\ddot{x}}_{\nu }(t)=-M{\omega }_{\nu }^{2}{x}_{\nu }(t)-\frac{1}{{x}_{0\nu }{N}_{k}}\sum\limits_{{{{\bf{k}}}}ij}{g}_{ij}^{\nu }({{{\bf{k}}}})\delta {\rho }_{{{{\bf{k}}}}ji}(t),$$
(6)

with M the mass of the unit cell and \({x}_{0\nu }=\sqrt{\hslash /(M{\omega }_{\nu })}\). For resonant pumping we can solve Eq. (6) with δρ from Eqs. (5). Indeed the phononic feedback on the electronic density matrix yields a correction of order \({{{\mathcal{O}}}}({n}^{2}{g}^{2})\) that can be neglected. This allows us to rewrite the EOM as

$${\ddot{x}}_{\nu }(t)=-{\omega }_{\nu }^{2}{x}_{\nu }(t)-n(t)\frac{{G}_{\nu }^{\lambda }}{\hslash }{\omega }_{\nu }\,{x}_{0\nu },$$
(7)

where the X-cph coupling \({G}_{\nu }^{\lambda }\) has been defined in Eq. (1). The appearance of the bare e-ph coupling g in \({G}_{\nu }^{\lambda }\) implies that the interaction between excitons and coherent phonons can be substantially larger than the interaction between excitons and phonons55.

Equation (7) admits an analytic solution for any time-dependent excitation density n(t). Let us discuss some physically relevant limiting cases. For displacive excitations48,49,50, i.e., if the duration TP of the pump pulse is much shorter than the phononic period Tν = 2π/ων, then we can set n(t) = n and solve Eq. (7) with initial conditions \({x}_{\nu }({T}_{{{{\rm{P}}}}})={\dot{x}}_{\nu }({T}_{{{{\rm{P}}}}})=0\). The solution is

$${x}_{\nu }(t)={x}_{0\nu }\frac{n{G}_{\nu }^{\lambda }}{\hslash {\omega }_{\nu }}[\cos {\omega }_{\nu }(t-{T}_{{{{\rm{P}}}}})-1].$$
(8)

Equation (8) shows that fast resonant pumping produces coherent oscillations of the nuclear displacements with amplitude \({a}_{\nu }=| {{\mathsf{x}}}_{\nu }| \equiv \frac{{x}_{0\nu }n| {G}_{\nu }^{\lambda }| }{\hslash {\omega }_{\nu }}\). In addition, the oscillations occur around an average displaced position \({{\mathsf{x}}}_{\nu }=-{{{\rm{sign}}}}({G}_{\nu }^{\lambda })| {{\mathsf{x}}}_{\nu }|\), that can be positive or negative, depending on the sign of the X-cph coupling. If, on the other hand, the pump duration is comparable or larger than the phonon period then the average displaced position does not change while the amplitude of the oscillations reduces: \({a}_{\nu }=({n}_{\nu }/n)| {{\mathsf{x}}}_{\nu }|\) with nν < n. The ��effective” densities nν = nν(TP) are decreasing functions of TP which approach n for TP ≪ Tν and zero for TP ≫ Tν. In fact, in the limit of very slow pumping the displacement \({x}_{\nu }(t)\,\approx \,{{\mathsf{x}}}_{\nu }\) becomes independent of time, consistently with the fact that the system is adiabatically driven toward a nonequilibrium steady-state.

Coherent modulation of the exciton energy

We now address how coherent phonons modify the exciton energy Eλ as measured in a transient optical experiments43,44,45,46,47. We are interested in probing the system at a time τ subsequent the phonon-induced dephasing of the electronic polarization29,56,57 (typical time-scales are a few hundreds of femtoseconds). In this depolarized regime the inter-band elements of the density matrix vanish58, i.e., ρkvc = 0, and the system is in a nonequilibrium steady-state characterized by a density matrix ρ = ρeq + δρ, with ρeq the ground state density matrix and δρ as in Eqs. (5). The quasi-particle Hamiltonian then reads \({h}_{{{{\bf{k}}}}}^{{{{\rm{qp}}}}}(\tau )={h}_{{{{\bf{k}}}}}^{{{{\rm{HSEX}}}}}[\rho ]+{h}_{{{{\bf{k}}}}}^{{{{\rm{Eh}}}}}(\tau )\), where the first term is the HSEX Hamiltonian (which is a functional of the time-independent density matrix ρ) and the second term is the time-dependent Ehrenfest potential given by \({h}_{{{{\bf{k}}}}ij}^{{{{\rm{Eh}}}}}(\tau )={\sum }_{\nu }{g}_{ij}^{\nu }({{{\bf{k}}}})\frac{{x}_{\nu }(\tau )}{{x}_{0\nu }}\). The Ehrenfest potential arises from the photo-excited coherent motion of the nuclei and acts on the electrons with an intensity proportional to the bare e-ph couplings g51.

For probe pulses shorter than the optical phonon periods the transient absorption spectrum can be obtained using the adiabatic approach of ref. 59. The absorption energies at time τ are the eigenvalues of the non-equilibrium BSE with Hamiltonian H(τ) = H + δH(τ) with δH(τ) = δHEh + δK. The change δK of the BSE kernel K = δhHSEX/δρ in Eq. (4) is due to the pump-induced change of the density matrix59. The term δHEh does instead accounts for the renormalization of the quasi-particle energies caused by the Ehrenfest potential. It is responsible for the temporal modulation of the excitonic energies and reads

$$\delta {H}_{{{{\bf{k}}}}vc,{{{{\bf{k}}}}}^{{\prime} }{v}^{{\prime} }{c}^{{\prime} }}^{{{{\rm{Eh}}}}}(\tau )={\delta }_{{{{\bf{k}}}}{{{{\bf{k}}}}}^{{\prime} }}\sum\limits_{\nu }\frac{{x}_{\nu }(\tau )}{{x}_{0\nu }}\left[{\delta }_{v{v}^{{\prime} }}{g}_{c{c}^{{\prime} }}^{\nu }({{{\bf{k}}}})-{\delta }_{c{c}^{{\prime} }}{g}_{{v}^{{\prime} }v}^{\nu }({{{\bf{k}}}})\right].$$
(9)

For small excitation densities n, both δHEh and δK can be treated perturbatively, yielding the τ-dependent modification of the \({\lambda }^{{\prime} }\) excitonic energy

$$\begin{array}{ll}{E}^{{\lambda }^{{\prime} }}(\tau )\,=\,{E}^{{\lambda }^{{\prime} }}+\sum\limits_{{{{\bf{k}}}}{{{{\bf{k}}}}}^{{\prime} }c{c}^{{\prime} }v{v}^{{\prime} }}{A}_{{{{\bf{k}}}}vc}^{{\lambda }^{{\prime} }* }[\delta {H}_{{{{\bf{k}}}}vc,{{{{\bf{k}}}}}^{{\prime} }{v}^{{\prime} }{c}^{{\prime} }}(\tau )]{A}_{{{{{\bf{k}}}}}^{{\prime} }{v}^{{\prime} }{c}^{{\prime} }}^{{\lambda }^{{\prime} }}\\ \qquad\quad\equiv \,{E}^{{\lambda }^{{\prime} }}+\delta {E}_{K}^{{\lambda }^{{\prime} }}+\delta {E}^{{\lambda }^{{\prime} }}(\tau ).\end{array}$$
(10)

The constant shift \(\delta {E}_{K}^{{\lambda }^{{\prime} }}\) is due to δK, and it accounts for the renormalization caused by electronic correlation effects like band-gap renormalization60,61,62, Pauli blocking61,63,64 and excited-state screening63,64,65. The explicit evaluation of \(\delta {E}_{K}^{{\lambda }^{{\prime} }}\) has been the subject of several papers and is not discussed further here. The time-dependent shift \(\delta {E}^{{\lambda }^{{\prime} }}(\tau )\) origins from δHEh. Using the expression for xν(τ) in Eq. (8), derived for TP ≪ Tν, we obtain that the modulation of the \({\lambda }^{{\prime} }\) exciton energy due to resonant pumping with the lowest exciton reads

$$\delta {E}^{{\lambda }^{{\prime} }}(\tau )=n\sum\limits_{\nu }\frac{{G}_{\nu }^{\lambda }{G}_{\nu }^{{\lambda }^{{\prime} }}}{\hslash {\omega }_{\nu }}\times \left[\frac{{n}_{\nu }}{n}\cos {\omega }_{\nu }\tau -1\right],$$
(11)

that for \({\lambda }^{{\prime} }=\lambda\) reduces to Eq. (2).

We conclude that the coupling with coherent phonons in resonantly excited materials is responsible for oscillations of the excitonic peak versus the impinging time τ of the probe43,44,45,46,47. The oscillation periods are dictated by the phonon frequencies and, for the lowest-energy exciton, the amplitude is proportional to the square of the X-cph coupling. Our result also predicts a likelihood of observing beatings when more phonons with different frequencies are strongly coupled to the exciton. This finding aligns with the outcomes of recent experiments in MoSe266.

Transient absorption in monolayer MoS2

We validate our analytical findings through numerical simulations of the transient absorbance in monolayer MoS2, see Fig. 1a. We use a tight-binding description of the band structure67 and a Rytova–Keldysh potential for the screened interaction68, and consider only the out-of-plane \({A}_{1}^{{\prime} }\) normal mode. Recent experiments suggest that the photoexcited dynamics is dominated by such mode43, see Fig. 1b. For the bare e-ph coupling we use \(g\,\approx \,\varepsilon {\mathbb{g}}\) with dielectric constant ε = 1 – screening at vanishing momentum transfer is negligible in two-dimensional materials11,68. The screened e-ph couplings \({{\mathbb{g}}}_{c{c}^{{\prime} }}^{{A}_{1}^{{\prime} }}\) and \({{\mathbb{g}}}_{v{v}^{{\prime} }}^{{A}_{1}^{{\prime} }}\) have been evaluated with the Quantum Espresso package69, and are displayed in Fig. 2. We propagate in time the electronic density matrix ρk(t) in the HSEX+Ehrenfest approximation using the CHEERS code70 (see Methods).

Fig. 1: Electronic, excitonic and vibrational properties of MoS2 monolayers.
figure 1

a Crystal structure of a MoS2 monolayer. b Geometrical representation of the \({A}_{1}^{{\prime} }\) normal mode. c Electronic band structure. The yellow box at the K-point highlights the splitting due to spin-orbit interaction. d Pictorial view of the bands involved in the formation of A and B excitons.

Full size image
Fig. 2: Electron-phonon couplings.
figure 2

Color plot of the coupling between the \({A}_{1}^{{\prime} }\) optical phonon and electrons in the conduction band shown in panel a and electrons in the valence band shown in panel b.

Full size image

The dynamics is initiated by a linearly polarized pump pulse of duration TP = 20 fs (much shorter than the optical period \({T}_{\nu = {A}_{1}^{{\prime} }}\simeq 82\) fs) and energy ℏωP = 1.9 eV, which is in resonace with the A exciton, see Fig. 1c, d. The pump pulse does therefore induce a displacive excitation and we expect \({n}_{\nu = {A}_{1}^{{\prime} }}\simeq n\). The excitation density for t > TP is found to be n = 0.0013, corresponding to a carrier density per unit area \({n}_{{{{\rm{x}}}}}=n/{A}_{{{{{\rm{MoS}}}}}_{2}}=1.5\times 1{0}^{12}\,{{{{\rm{cm}}}}}^{-2}\) (\({A}_{{{{{\rm{MoS}}}}}_{2}}=8.8\) Å2 being the area of the unit cell). About 400 fs after the photoexcitation the system enters the depolarized regime with ρkvc ≈ 0 and \({\rho }_{{{{\bf{k}}}}c{c}^{{\prime} }}\) and \({\rho }_{{{{\bf{k}}}}v{v}^{{\prime} }}\) independent of time. In Fig. 3a we plot the occupations of spin-up electrons in the conduction band, i.e., nkc↑(t) = ρkc↑c↑(t), at t = 700 fs. In the spin-up sector the A and B excitons are located at the K and \({K}^{{\prime} }\) valleys respectively, see again Fig. 1c, d, whereas in the spin-down sector this configuration reverses4. In Fig. 3 we compare nkc↑(t) (panel b) with the analytic value nkc↑ ≡ δρkc↑c↑ from Eq. (5a) (panel c), which is proportional to the square modulus of the exciton wavefunction. The agreement is fairly good, thus confirming that resonant pumping generates a finite density of A excitons. The small discrepancy is due to the finite duration of the pump, which prevents achieving a perfect resonant condition. As a result, a minimal contamination arises from the generation of a small density of B excitons, specifically observed in Fig. 3a at the \({K}^{{\prime} }\) valley.

Fig. 3: Photo-excited carrier occupations.
figure 3

a Plot of the density of spin-up electrons in the conduction band after 700 fs (depolarized regime), showing a predominance of the A exciton (K-point) over the B exciton (\({K}^{{\prime} }{{{\rm{point}}}}\)). b Magnification of panel a around the K-point. c Analytic density of spin-up electrons in the conduction band around the K-point as obtained from Eq. (5a).

Full size image

To obtain the time-resolved absorption spectrum, we further imping the system with a broadband probe pulse at different delay times τ. We then calculate for all times t the change in the density matrix Δρ(t, τ) = ρP+p(t, τ) − ρP(t), resulting from a simulation with pump and probe and a simulation with only the pump59,71. The transient absorption spectrum is calculated as59,71

$${\mathfrak{S}}(\omega ,\tau )=-2\omega {{{\rm{Im}}}}[{{{\bf{e}}}}(\omega ,\tau )\cdot {{{\bf{d}}}}(\omega ,\tau )],$$
(12)

where e(ω, τ) is the Fourier transform of the probe field and d(ω, τ) is the Fourier transform of the probe-induced dipole moment d(t, τ). The latter is related to the change of the density matrix via \({{{\bf{d}}}}(t,\tau )={\sum }_{{{{\bf{k}}}}}{{{\rm{Tr}}}}[{{{{\bf{d}}}}}_{{{{\bf{k}}}}}\Delta {\rho }_{{{{\bf{k}}}}}(t,\tau )]\), where dk is the dipole transition matrix. Equation (12) generalizes the equilibrium absorption spectrum to out-of-equilibrium conditions. In the absence of pump fields (equilibrium condition) we have \({{{\bf{d}}}}(\omega ,\tau )={{{\rm{Tr}}}}[{{{\boldsymbol{\chi }}}}(\omega )\cdot {{{\bf{e}}}}(\omega )]\) independent of τ, with χ(ω) the dipole-dipole response function, and \({\mathfrak{S}}(\omega ,\tau )\) reduces to equilibrium absorption spectrum59,72.

We have conducted simulations for several equidistant delay times τn = nΔτ, encompassing approximately 100 distinct values spanning a relatively long time-window (~1 ps). The dependence of the spectrum on τ originates from the time-dependent HSEX+Ehrenfest potential: the specific timing at which the probe interacts with the system determines the value of hqp(τ), resulting in different frequencies for Δρkvc(t, τ).

In Fig. 4a we display the logarthimic plot of \({\mathfrak{S}}(\omega ,\tau )\) for energies ℏω up to a few meV below the exciton energy EA ≃ 1.9 eV. The spectrum changes wildly until the pump-induced polarization is sizable (coherent regime), i.e., up to τ ≃ 400 fs. Both the position and the amplitude of the excitonic peak oscillate with a period of TA = 2π/EA ~ 2.2 fs. In our simulations, however, the system is probed every Δτ = 7.75 fs, preventing us to resolve the faster time-scale TA. Instead, an artificial period T ~ 31 fs [corresponding to the matching condition T ~ pTA ~ qΔτ with p and q prime numbers] is observed.

The position EA(τ) of the A exciton energy, calculated as the maximum \({\max }_{\omega }{\mathfrak{S}}(\omega ,\tau )\), is shown in Fig. 4b (blue dots). The exciton is initially blue-shifted with respect to its equilibrium value. This blue-shift undergoes relaxation within about 200 fs, eventually transitioning into a stable red-shift \(\delta {E}_{K}^{{{{\rm{A}}}}}\approx 2\) meV due to electronic mechanisms (change of the BSE kernel), see Eq. (10). The decay time of this transition is governed by the polarization lifetime which, in turns, depends strongly on the excitation density and it can vary from hundreds of femtoseconds to picoseconds. Our observations align closely with the work by Calati et al.73, where a similar blue-to-red transition has been reported in resonantly excited WS2. Entering the depolarized (or incoherent) regime the wild oscillations gradually attenuate, and the intensity of the signal begins to oscillate at the phonon frequency \({\omega }_{{A}_{1}^{{\prime} }}\), see Fig. 4a. For τ > 400 fs, the exciton position δEA(τ) features oscillations at the very same frequency around the red-shifted value \({E}^{{{{\rm{A}}}}}+\delta {E}_{K}^{{{{\rm{A}}}}}\), in qualitative agreement with Eq. (10) and Eq. (2).

Inspecting the amplitude of the oscillations, we also appreciate an excellent quantitative agreement, see Fig. 4c. The red curve is the plot of Eq. (2) with \({n}_{{A}_{1}^{{\prime} }}=n\) (displacive excitation) and n = 0.0013 the excitation density created by the pump. The value of the X-cph coupling is \({G}_{{A}_{1}^{{\prime} }}^{{{{\rm{A}}}}}\approx 50\) meV, implying that resonant pumping with the A exciton brings the MoS2 monolayer toward a nonequilibrium steady-state characterized by a vertical compression of the sulfur atoms (\({{\mathsf{x}}}_{{A}_{1}^{{\prime} }} < 0\)).

As the analytic formula in Eq. (2) has been derived in the limit of low excitation density and short pump duration (displacive regime), we have explored the extent to which its applicability is affected when the photoexcitation deviates from these conditions, see Fig. 5. To facilitate comparisons a τ-slice of Fig. 4c is shown in panel a. In panel b we consider longer pump pulses with duration TP = 50 fs, while still maintaing a small excitation density nx = 4 × 1011cm−2. In this case a correction due to the effective density \({n}_{{A}_{1}^{{\prime} }}\) must be introduced [see Eqs. (2) and (8)] since the photoexcitation is no longer displacive. The value \({n}_{{A}_{1}^{{\prime} }}/n={a}_{{A}_{1}^{{\prime} }}/| {{\mathsf{x}}}_{{A}_{1}^{{\prime} }}| =0.71\) is extracted by taking the ratio between the amplitude of the oscillations \({a}_{{A}_{1}^{{\prime} }}\) of the \({A}_{1}^{{\prime} }\) mode as obtained from the numerical simulation and the analytic value \(| {{\mathsf{x}}}_{{A}_{1}^{{\prime} }}| =\frac{{x}_{0{A}_{1}^{{\prime} }}n| {G}_{{A}_{1}^{{\prime} }}^{\lambda }| }{\hslash {\omega }_{{A}_{1}^{{\prime} }}}\). Remarkably, the validity of the analytic formula extends well beyond the initially defined applicability region. In panel c we consider a very long pump TP = 100 fs and a relatively high excitation density nx = 2.3 × 1012cm−2. Applying the correction \({n}_{{A}_{1}^{{\prime} }}/n=0.10\), the agreement between the simulated and predicted excitonic oscillations is still reasonably good. Further increasing the excitation density, nx = 5.6 × 1012cm−2, the agreement begins to deteriorate, even within the displacive regime, see panel d, where \({n}_{{A}_{1}^{{\prime} }}/n=1\).

Fig. 4: Exciton dynamics due to coherent phonons.
figure 4

a Color plot of the logarithmic transient absorption spectrum \(\ln \left[\frac{{\mathfrak{S}}(\omega ,\tau )-{{\mathfrak{S}}}_{\min }}{{{\mathfrak{S}}}_{\max }-{{\mathfrak{S}}}_{\min }}\right]\) versus energy (vertical axis) and delay (horizontal axis). Here \({{\mathfrak{S}}}_{\min }\) and \({{\mathfrak{S}}}_{\max }\) denote the minimum and maximum values of \({\mathfrak{S}}(\omega ,\tau )\) in the displayed region. Blue dots indicate the position of the maxima of \({\mathfrak{S}}(\omega ,\tau )\) for different delays. The electronic red-shift \(\delta {E}_{K}^{{{{\rm{A}}}}}\) and the phonon-induced coherent modulation δEA(τ) are also indicated. b Temporal evolution of the exciton energy EA(τ) as a function of the pump-probe delay τ. The reference equilibrium value EA(0) = 1.9 eV is shown as a dashed line, while blue and red backgrounds are used visualize blue- and red-shift regions respectively. c Comparison of the position of the exciton peak as obtained from the real-time (RT) simulations of panel a (blue dots) and the analytic formula in Eq. (2) (red curve).

Full size image
Fig. 5: Pump-dependent dynamics of exciton oscillations.
figure 5

Comparison of the oscillations of the exciton peak δEA from RT simulations (blue dots) with the analytic result of Eq. (2) (solid curves) for different pump dutations TP and excitation densities nx. Panel a: TP = 20 fs, nx = 1.5 × 1012cm−2 (\({n}_{{A}_{1}^{{\prime} }}=n\)); Panel b: TP = 50 fs, nx = 4.0 × 1011cm−2 (\({n}_{{A}_{1}^{{\prime} }}=0.71n\)); Panel c: TP = 100 fs, nx = 2.3 × 1012cm−2 (\({n}_{{A}_{1}^{{\prime} }}=0.1n\)); Panel d: TP = 20 fs, nx = 5.6 × 1012cm−2 (\({n}_{{A}_{1}^{{\prime} }}=n\)).

Full size image

Discussion

We present the many-body theory of the interaction between resonantly pumped excitons and coherent phonons, and derive the expression of the X-cph coupling. The magnitude of this coupling can exceed that of the X-ph coupling in three dimensions, given that the bare (screened) e-ph matrix elements are used to calculate the former (latter). We also derive a readily applicable formula to quantify the modulations of the excitonic energies in resonantly pumped materials43,44,45,46,47. The accuracy of the formula is demonstrated through comparisons with model simulations of transient absorption in a MoS2 monolayer. Our findings provide a straightforward means to directly extract the value of the X-cph coupling, along with accurately determining the magnitude of the excitation density.

Methods

Real-time HSEX + Eherenfest method

We describe in detail the method that we use to propagate in time the electronic one-particle density matrix ρkij together with the phonon displacement xν(t) and momentum pν(t). Here k denotes the cristal momentum, i, j the spin-band indices that can be valence (v) or conduction (c), while ν denotes the phonon mode. Our method is based on the Hartree plus statically screened exchange (HSEX) plus Ehrenfest approximation of many-body theory. The coupled equations of motion to be propagated in time read51,54

$$\begin{array}{rcl}i\frac{d}{dt}{\rho }_{{{{\bf{k}}}}ij}(t)&=&{[{h}_{{{{\bf{k}}}}}^{{{{\rm{qp}}}}}(t),{\rho }_{{{{\bf{k}}}}}(t)]}_{ij}+{I}_{{{{\bf{k}}}}ij}^{{{{\rm{coll}}}}}(t),\\ {p}_{\nu }(t)&=&-{\omega }_{\nu }{x}_{\nu }(t)-\sum\limits_{{{{\bf{k}}}},ij}{g}_{ij}^{\nu }({{{\bf{k}}}},0){\rho }_{{{{\bf{k}}}}ji}(t)\\ \frac{d}{dt}{x}_{\nu }(t)&=&{\omega }_{\nu }{p}_{\nu }(t).\end{array}$$
(13)

In Eq. (13) hqp is the quasi-particle hamiltonian given by

$$\begin{array}{rcl}{h}_{{{{\bf{k}}}}ij}^{{{{\rm{qp}}}}}(t)&=&{h}_{{{{\bf{k}}}}ij}^{{{{\rm{HSEX}}}}}(t)+{h}_{{{{\bf{k}}}}ij}^{{{{\rm{Eh}}}}}(t)\\ &=&{\delta }_{ij}{\epsilon }_{{{{\bf{k}}}}i}+\sum\limits_{{{{{\bf{k}}}}}_{1}mn}({V}_{imnj}^{{{{\bf{0}}}}{{{\bf{k}}}}{{{{\bf{k}}}}}_{1}}-{V}_{imjn}^{({{{\bf{k}}}}-{{{{\bf{k}}}}}_{1}){{{\bf{k}}}}{{{{\bf{k}}}}}_{1}})\delta {\rho }_{{{{{\bf{k}}}}}_{1}nm}(t)+{{{\bf{E}}}}(t)\cdot {{{{\bf{d}}}}}_{{{{\bf{k}}}}ij}+\sum\limits_{\nu }{x}_{\nu }(t){g}_{ij}^{\nu }({{{\bf{k}}}})\end{array}$$
(14)

This hamiltonian accounts for (i) the electron kinetic energy (via the band structure ϵki), (ii) the interaction with the external laser pulses E(t) (via the dipole matrix elements dkij), (iii) the HSEX self-energy (via the electron-electron interaction \({V}_{imjn}^{{{{\bf{q}}}}{{{\bf{k}}}}{{{{\bf{k}}}}}^{{\prime} }}\)), and (iv) the Ehrenfest potential describing the coupling of electrons with the nuclear displacements (via the bare electron-phonon coupling \({g}_{ij}^{\nu }({{{\bf{k}}}})\equiv {g}_{ij}^{\nu }({{{\bf{k}}}},0)\) at vanishing momentum transfer). The presence of the HSEX and Ehrenfest potentials guarantees that strong excitonic effects and the coupling of excitons to coherent phonons are included in the photoexcited dynamics. The SEX interaction is modeled by the Rytova–Keldysh potential68,74 (see also below) which has demonstrated remarkable accuracy in reproducing excitonic binding energies in 2D materials2,63,75. The evaluation of the Ehrenfest potential and the effective forces acting on the nuclei require only the intra-band electron-phonon couplings \({g}_{ii}^{\nu }({{{\bf{k}}}})\) at q = 051. The couplings g can be evaluated according to

$${g}_{ii}^{\nu }({{{\bf{k}}}})=\sqrt{\frac{\hslash }{2M{\omega }_{\nu }}}\frac{\partial {E}_{{{{\bf{k}}}}i}}{\partial {x}_{\nu }},$$
(15)

where Eki is the Kohn–Sham energy of the i-th band at momentum k, ων = ων(q = 0) is the frequency of mode ν at the Γ point, and uν denotes the ion displacement along the phonon mode ν with momentum q = 0.

In MoS2 monolayers the mostly coupled phonon is the out-of-plane \({A}_{1}^{{\prime} }\) optical mode43. Therefore, only this specific normal mode has been included in our simulations. The last term in the equation of motion for ρ is the collision integral Icoll(t), that incorporates dynamical correlation effects responsible for, e.g., inelastic quasi-particle scattering, decoherence and dynamical screening58,76,77,78. At low and moderate excitation densities the dominant contribution to the collision integral stems from electron-phonon scattering which drives the system toward a depolarized regime characterized by the disappearance of valence-conduction elements ρk,vc of the density matrix53,54,58. For this reason all matrix elements of the collision integral are neglected, except for the valence-conduction ones that we approximate as \({I}_{{{{\bf{k}}}},vc}^{{{{\rm{coll}}}}}(t)=-2i\gamma {\rho }_{{{{\bf{k}}}},cv}(t)\). Here ℏ/γ = 250 fs is the dephasing time of the electronic polarization experimentally observed in MoS2 at low temperature57. The three coupled equations of motion in Eq. (13) are numerically integrated by using a 4th order Runge-Kutta solver with a time step Δt = 0.1 fs as implemented in the CHEERS code70.

For our simulations in MoS2 we select as active space the two highest valence and the two lowest conduction spin-bands. Therefore the spin-band index i can take the 4 values i = {(v↑), (v↓), (c↑), (c↓)}. In this way spin-orbit effects responsible for the energy splitting between A and B excitons are included.

From a knowledge of the electronic density matrix we can evaluate the momentum resolved carrier populations nki(t) according to

$${n}_{{{{\bf{k}}}}i}(t)={\rho }_{{{{\bf{k}}}}ii}(t).$$
(16)

Using the off-diagonal elements of ρ we can also calculate the probe-induced dipole moment according to

$${{{\bf{d}}}}(t)=\sum\limits_{{{{\bf{k}}}}}{{{{\bf{d}}}}}_{{{{\bf{k}}}}}(t),\quad {{{{\bf{d}}}}}_{{{{\bf{k}}}}}(t)=\sum\limits_{ij}{\rho }_{{{{\bf{k}}}}ij}(t){{{{\bf{d}}}}}_{{{{\bf{k}}}}ji}$$
(17)

Modelization of the MoS2 monolayer

We calculate the spin-orbit dependent band structure of a MoS2 monolayer from the tight-binding parametrization provided in ref. 67 that gives low-energy bands ϵki in very good agreement with DFT calculations. The conduction bands are rigidly shifted-up by 0.6 eV, to match the quasiparticle bandgap measured in ARPES experiments79. For the RT simulations we consider the two highest valence bands and the two lowest conduction bands as our active space. The Coulomb integral \({V}_{imjn}^{{{{\bf{q}}}}{{{\bf{k}}}}{{{{\bf{k}}}}}^{{\prime} }}\) accounting for the scattering amplitude of two electrons in bands j and n with quasimomenta \({{{{\bf{k}}}}}^{{\prime} }+{{{\bf{q}}}}\) and k–q to end up in bands m and i with quasimomenta \({{{{\bf{k}}}}}^{{\prime} }\) and k respectively is calculated according to \({V}_{imjn}^{{{{\bf{q}}}}{{{\bf{k}}}}{{{{\bf{k}}}}}^{{\prime} }}={v}_{q}({{{{\bf{U}}}}}_{{{{\bf{k}}}}i}^{{\dagger} }\cdot {{{{\bf{U}}}}}_{{{{\bf{k}}}}-{{{\bf{q}}}}n})({{{{\bf{U}}}}}_{{{{{\bf{k}}}}}^{{\prime} }m}^{{\dagger} }\cdot {{{{\bf{U}}}}}_{{{{{\bf{k}}}}}^{{\prime} }+{{{\bf{q}}}}j})\)80, where Uki are the eigenvectors of the Bloch Hamiltonian Hk and vq is the Rytova–Keldysh potential68,74 in momentum space, i.e.

$${v}_{q}=\frac{2\pi }{\epsilon q(1+{r}_{0}q)}.$$
(18)

In the above equation q = ∣q∣, r0 = 33.875 Å/ϵ80, and ϵ = 2.5 is the dielectric constant of a typical substrate (the topstrate is supposed to be air). The regularization of the diverging value of vq at the Γ point is expressed as81,82\(v(0)\to \frac{1}{\Omega }{\int}_{\Omega }d{{{\bf{q}}}}v(q)\). Here Ω represents a 2D domain around q = 0 of linear dimension determined by the discretization step of the first Brillouin zone. In this study, the first Brillouin zone has been discretized using a C6v-symmetric grid comprising 3072 k-points. In semiconductor materials the Coulomb integrals that do not conserve the number of valence and conduction electrons are typically negligible83 and have been excluded from our calculations. Finally the dipole matrix elements are evaluated according to63,84

$${{{{\bf{d}}}}}_{{{{\bf{k}}}}ij}=\frac{1}{i}\frac{1}{{\epsilon }_{{{{\bf{k}}}}i}-{\epsilon }_{{{{\bf{k}}}}j}}{{{{\bf{U}}}}}_{{{{\bf{k}}}}i}^{{\dagger} }\cdot {\partial }_{{{{\bf{k}}}}}{H}_{{{{\bf{k}}}}}\cdot {{{{\bf{U}}}}}_{{{{\bf{k}}}}j}.$$
(19)

Concerning the phonon input, the electron-phonon coupling in Eq. (15) are obtained by displacing manually the ions along the phonon eigenvector directions (with equilibrium lattice constant a = 3.12 Å) and by calculating the corresponding difference in the band structure. This procedure is very accurate for the out-of-plane optical mode \({A}_{1}^{{\prime} }\)85. The Kohn–Sham energies Eki for a given displaced configuration are calculated with Quantum ESPRESSO using LDA in a 24 × 24 grid. The phonon frequencies ων are evaluated within the same grid.