Cavity-Tuned Exciton Dynamics in Transition Metal Dichalcogenides Monolayers

A fully quantum, numerically accurate methodology is presented for the simulation of the exciton dynamics and time-resolved fluorescence of cavity-tuned two-dimensional (2D) materials at finite temperatures. This approach was specifically applied to a monolayer WSe2 system. Our methodology enabled us to identify the dynamical and spectroscopic signatures of polaronic and polaritonic effects and to elucidate their characteristic timescales across a range of exciton–cavity couplings. The approach employed can be extended to simulation of various cavity-tuned 2D materials, specifically for exploring finite temperature nonlinear spectroscopic signals.


Introduction
Monolayer transition metal dichalcogenides (TMDs) are celebrated for their broken inversion symmetry, resulting in a non-zero Berry curvature and robust valley-spin interactions.Such valley-spin interactions enable selective valley polarization through optical or electromagnetic means [1,2].The reduction in the layer number results in the transformation of TMDs from indirect to direct band gap semiconductors [3,4].This layer-number-dependent symmetry variation also allows for the targeted manipulation of valley properties, establishing TMDs as a pivotal platform for valleytronic device development.In tungsten-based monolayers, low-temperature photoluminescence studies have shown that dark exciton states [5][6][7][8], which are momentum-forbidden and significantly influence the optical and electronic behavior of TMDs [9][10][11][12], exhibit lower energy levels than their bright counterparts.
Additionally, the interaction of two-dimensional (2D) materials with light, crucial for devices such as nanolasers [13,14], can be enhanced by coupling these materials with optical cavities or through intrinsic polaritonic resonances [15,16].The strong interaction between light and matter within optical cavities leads to the creation of quasiparticles termed polaritons.In organic microcavities, when the vibronic coupling between excitons and phonons is comparable in strength to the exciton-photon coupling, multiple polariton branches have been experimentally observed [17,18] alongside phonon-assisted relaxation processes [19,20].Recently, these phenomena have been directly utilized to investigate polaritonic nonlinearities at the single-photon level [21].These entities underpin a wide array of both quantum and classical behaviors across various materials and spectral regions, such as polariton lasing [22], polariton-induced blockade [23,24], and the Bose-Einstein condensation of polaritons [25].Recent developments in 2D TMDs interfaced with photonic architectures have paved the way for detailed studies of exciton-polaritons on an atomic level.The valley-selected property of TMDs excitons, when utilized in hybrid photonic systems, has enabled the creation of valley-polarized exciton-polaritons [26].Moreover, excitons located in the K and K' valleys are interconnected through time reversal symmetry.Consequently, removing valley degeneracy in these systems might facilitate the study of topological exciton-polaritons that exhibit broken time reversal symmetry at visible frequencies [27].Another compelling dimension of exciton behavior in 2D TMDs emerges when two monolayers are angularly misaligned and superposed, creating bilayer heterostructures [28,29].The resulting moiré patterns lead to moiré excitons within individual layers, which show a pronounced increase in nonlinearity due to exciton blockade effects.This significant interaction between moiré excitons and photonic modes in cavities can produce polaritons with marked nonlinearity [30], forming an excellent platform for studying highly interactive exciton-polaritons.Such integration facilitates the development of compact, on-chip devices [31] by coupling excitons in TMDs with photonic modes within planar nanocavities [32,33], paving the way for novel light-matter interaction platforms.
Photoluminescence (PL) spectroscopy remains a fundamental technique for probing quantum phenomena and complex interactions in TMD monolayers [34].At a low temperature, PL reveals the presence of bound excitonic complexes, such as biexcitons and trions [5,[35][36][37].The significant role of phonon-assisted recombination in shaping the PL emission characteristics can be highlighted, particularly within cavity-controlled materials, leading to asymmetrical line shapes and enhanced sidebands-a phenomenon termed cavity-coupled PL emission [31].Time-resolved fluorescence (TRF) spectroscopy [38][39][40] further allows the examination of these interactions over brief intervals, providing deeper insights into exciton dynamics and relaxation processes.
Despite the utility of conventional theoretical approaches such as quantum master equations [31] and cluster expansion methods [5], these often fall short in scenarios involving strong exciton-photon or exciton-phonon interactions.However, the multiple Davydov Ansatz (mDA) technique has recently emerged as a robust, numerically precise method [41,42], applied to a diverse array of multidimensional challenges, including analyzing Landau-Zener transitions [43,44], ultrafast dynamics at conical intersections (CIs) [45,46], and hole-magnon dynamics in an antiferromagnet [47,48].The mDA method is a variational method to solve the many-body time-dependent Schrödinger equation, where the wave function for bosonic degrees of freedom is expanded into multiple coherent states for each boson mode in a multispecies boson bath of arbitrary spectral density functions.If the multiplicity is high enough, the so-obtained solutions converge to an exact solution.Here, we employ the mDA framework, complemented by thermofield dynamics (TFDs) representation [49,50], to offer a comprehensive, microscopically accurate simulation method for excitonic dynamics and spectroscopic responses of cavity-tuned single-layer WSe2 at finite temperature and multiple exciton-cavity (EC) coupling strengths.
In Section 2, we establish a microscopic model of the WSe 2 monolayer embedded in a microcavity and introduce the methodologies adopted.Section 3 is devoted to a comprehensive discussion of calculated exciton populations and TRF spectra for a range of EC coupling strengths from weak to strong.Conclusions are offered in Section 4.

Microscopic Model and Methodologies
PL and TRF profiles of WSe 2 monolayers are predominantly determined by KK, KQ, and KK ′ excitons, of which the first K represents hole momentum in the valence band, and the second symbol indicates the electron momentum in the conduction band, corresponding to momenta at 1Γ, 2Λ, and 3K respectively.Excitons of higher energy than the above states have a negligible influence on the dynamic properties and optical signals at room temperature and below.The phonon-mediated PL is attributed to the KQ and KK ′ intervalley dark states, which is in contrast to the KK intravalley excitons that are optically bright and yield direct PL emission.The contrasting dynamics of direct versus indirect PL are detailed in Figure 1.Following this qualitative picture, the exciton-polariton dynamics of the cavity-modified WSe 2 monolayer can be described by a Hamiltonian containing three excitonic states, a cavity mode, and several dominant phonon modes [5,31,51,52]: Here, operators , respectively, exciton states (iQ ∥ ), circularly polarized photon states σ + and phonon modes α with momentum q ∥ , and frequency Ω α,q ∥ .E (iQ ∥ ) is the excitonic dispersion relation, where the symbol ∥ denotes the momentum's in-plane component.The strength of coupling between the exciton and cavity photon modes is quantified by the coefficient M σ + .The tensor D ij αq ∥ , which maps the coupling of excitons with phonons, can be computed through the density functional theory (DFT) formalism (mean effective deformation potential approximation) [7,53]: where Φ i k is the excitonic wave function in the momentum space, which can be obtained by solving the Wannier equation.g c/vkij αq ∥ are the electron-phonon interaction coefficients within the conduction and valence bands [54,55], where α and q ∥ define the phonon modes scattered between the excitonic band index i and j.
are the electron-phonon matrix elements in conduction/valence band involving integrals over the unit cell (UC).M encapsulates the total mass of atoms within a unit cell, while |i, k⟩ describes the Bloch eigenstates ψ i,k (r) characterized by wave vector k and an excitonic band index i.The potential perturbation, ∆V α q ∥ , is derived using density functional perturbation theory.The first-order deformation potential, or the acoustic deformation potential, is expressed as ⟨j, k + q ∥ |∆V α q ∥ |i, k⟩/|q ∥ |.It can be shown that acoustic phonons at the Γ do not contribute to the electron-phonon coupling.The zero-order deformation potential, also known as the optical deformation potential, is defined as ⟨j, k + q ∥ |∆V α q ∥ |i, k⟩ [54].Detailed parameters for this potential are documented in Ref. [56].In this work, we neglect cavity energy losses, which can be realized in fabricated photonic crystal nanobeam (PhCnB) cavities with a notably high-quality factor [57].For instance, the exciton-photon coupling of a thin excitonic material positioned at the center of a high-quality, symmetric Fabry-Pérot cavity near the cavity resonance can be fine-tuned by adjusting the mirror reflectivity and the cavity length [58][59][60].
To treat finite temperature effects with the TFD method, we introduce the extended Hamiltonian [61][62][63]: where is the phonon Hamiltonian acting in the fictitious "tilde" vibrational space.Having performed the thermal Bogoliubov transformation, we arrive at the TFD Hamiltonian [49,50]: in which is the Bogoliubov operator, and are the mixing angles, which account for the influence of the temperature on the excitonphonon couplings.The Hamiltonian of Equation ( 6) commutes with the number operator and conserves the number of excitons.In this work, we consider the dynamics within the singly excited excitonic manifold, where ⟨N ex ⟩ = 1.In this manifold, the solution of the time-dependent Schrödinger equation with the TFD Hamiltonian of Equation ( 6) can be represented in terms of the mDA wave function of multiplicity M: Here, the indices i = 0, 1, 2, 3 denote the photon state, KK, KQ, and KK ′ excitons, respectively.B mi (t) represents the time-dependent exciton amplitude, where the indices m stand for the mth coherent state superposition (with a total of M superpositions).The operators b † l (b l ) and b † l (b l ) are the physical and tilde phonon's creation (annihilation) operators, respectively.The variables u ml (t) and ũml (t) capture displacements of the physical and tilde phonons in mth coherent state.For the calculation of exciton dynamics, the initial population is solely in the photon state, and the initial elements of u ml (t) and ũml (t) are random numbers of the order 10 −4 .The parameters B mi (t), u ml (t), and ũml (t) are determined through the equations of motion obtained following the variational principle [41,46,64], as detailed in Appendix A. All observables studied in the present work were evaluated in |D M 2 (t)⟩, as described in Appendix B.
To establish a realistic model, ab initio input parameters entering the Hamiltonian include the electronic band structure [65], phonon dispersion [56], dielectric constants [66] and electronphonon coupling elements [56].In our calculations, we assumed that E (1Γ) = 1.724 eV [67,68], E (2Λ) = 1.69 eV, and E (3K) = 1.678 eV, where the spectral exciton separations between KK, KQ, and KK ′ excitons were taken from ref. [5].The exciton energies and wave functions were derived in the effective mass approximation by solving the Wannier equation [69][70][71].This approach leads to an exciton dispersion, which is comparable to that from DFT calculations based on the Bethe-Salpeter equation [72][73][74].To amplify cavityprompted phenomena, a single cavity mode [5,75] was taken to be in resonance with the luminescent KK exciton, where hω σ + = 1.724 eV.The dielectric constant [5,66] was set to be 4.5 for the hBN-encapsulated TMDs monolayers.For the exciton-cavity couplings, we only consider the interactions between the bright KK exciton and photons, as there is no transition dipole moment of momentum-dark states.Within the photonic crystal cavities, the coupling between excitons and photons is modifiable within a range from 4 meV to 14 meV, as documented by Rosser et al. [31]; for our calculations, M σ + = 0, 4, 8, 12 meV [76][77][78] were used.Table 1 lists the employed phonon frequencies [5,56], illustrating that both longitudinal acoustic (LA) and transverse acoustic (TA) phonon branches exhibit a linear dispersion nearing the long-wavelength limit, aligning their frequencies at zero when q = 0.The optical phonon branches with significant exciton interactions include the homopolar (A 1 ) modes, characterized by out-of-plane vibrations, and the in-plane longitudinal (LO) and transversal (TO) modes.
The Hamiltonian of Equation ( 6) incorporates 23 phonon modes, three excitonic states (KK, KQ, KK ′ ), and a single photonic mode.The number of physical phonon modes (i.e., 23) surpasses 13, as outlined in Table 1, which is due to the inclusion of ± modes b α,±q ∥ with α = KQ, KK ′ .The total number of phonon modes (i.e., 46) comes about by the duplication of the number of physical phonons owing to the inclusion of the tilde modes bα,q ∥ , crucial for addressing temperature effects within the TFD framework.The model omits holes near the K point due to their KK ′ symmetry and excludes spin-forbidden dark states, which are irrelevant for the ultrafast dynamics under study.This model effectively captures coherent interactions among excitons, phonons, and photons, and it incorporates environmental dephasing as proposed, e.g., in Ref. [5].
The impact of the temperature on the exciton dynamics and TRF spectra has been comprehensively investigated in Ref. [55].Here, we set the temperature to 75 K and focused on a detailed study of how the strength of the EC coupling affects the dynamical and spectroscopic observables.All our numerical results are proven to be convergent for the multiplicity M = 48 of the mDA wave function of Equation (10).

Exciton Populations
Figure 2 provides a comprehensive analysis of how varying ECs influence the dynamics of specific photonic and excitonic modes.In the absence of EC coupling in panel (a), the photon mode population persists, reflecting the lack of energy exchange between photons and excitons.For EC = 4 meV (red lines), the photon mode population decreases, and the excitonic mode population increases almost linearly with time, exhibiting the short-time ballistic transport on the timescale of ∼400 fs.As the EC coupling strength increases, additional oscillations and nonmonotonic behaviors emerge, as depicted in panels (b)-(d).This is indicative of exciton-polariton formation resulting from intensified exciton-photon interactions within the cavity.It is the value of EC = 8 meV (yellow lines) that marks the changeover to the oscillatory strong EC coupling regime.The onset of this regime is manifested through the nonmonotonic evolution of the cavity mode population in panel (a), which reaches a minimum (almost zero) at t ≈ 350 fs (when EC = 8 meV) and t ≈ 200 fs (when EC = 12 meV) followed by the increase at longer times.Irrespective of the EC coupling strength, the excitonic states are populated sequentially [55] at short times: the KK state is populated first, followed by the KQ and KK ′ states.Enhanced exciton-cavity coupling accelerates the transfer of photon population to excitons, as evidenced by the purple lines for EC = 12 meV, reaching their minima (a) and maxima (b)-(d) more swiftly.At an EC of 12 meV both bright KK and dark KQ excitons acquire rapid population growth initially at 100 fs and 200 fs, respectively, followed by noticeable oscillations.Remarkably, an increased EC consistently sustains a higher (though nonmonotonic as a function of time) population level in the dark KK ′ exciton due to phonon-assisted transfer when the polaritonic state is formed.This observation, which demonstrates robustness of the cavityinduced population of the KK ′ excitons within a wide range of EC coupling strengths, may become important for practical applications.

TRF
To delve deeper into the dynamics of ultrafast exciton-polariton formation in WSe 2 monolayers, we computed the TRF spectra S(ω, t) presuming instantaneous excitation of the system by the pump pulse.These spectra, which describe the rate of emission of photons of frequency ω at time t [64,79,80], are evaluated with the mDA wave function, as described in Appendix B. For obtaining the spectral shape of S(ω, t), we used the overdamped harmonic oscillator lineshape function of Equation (A9) with the Stokes shift parameter λ = 5 meV and the memory rate parameter Λ = 30 meV.
A 3D view of the TRF spectra S(ω, t) is presented in Figure 3 for different ECs from zero (a) through 4 eV (b) and 8 eV (c) to 12 eV (d).The spectra are characterized by a high initial (at t = 0) peak revealing the bright KK exciton (direct PL), the intensity of which drops substantially on the timescale of several dozens femtoseconds.Rapid fluorescence attenuation of the bright KK exciton is culminated in panel (d), which is characterized by an approximately 90% loss in intensity within the initial 100 fs.This behavior is indicative of swift internal conversion processes at the (cavity-induced) CIs [45,46,81].In contrast, the absence of a cavity (EC = 0, panel (a)) sees the fluorescence intensity of direct PL at about 50% in the first 100 fs.From panels (a) to (d), the tendency of fast fluorescence decay with the EC coupling agrees with the enhanced population of the dark states (KQ and KK ′ ) indicated in Figure 2c,d.
In the time domain, S(ω, t) exhibits oscillatory behavior, which mirrors cavity-mediated exciton-phonon wave packet motion.This motion is especially pronounced for EC = 0 (panel a).However, it is not strictly periodic, since the spectral maxima are separated by the intermittent time interval of about ∼80-100 fs.Interestingly, these intervals are substantially shorter than the period corresponding to the fastest phonon mode (2π/Ω LO,Λ = 127 fs.This is a clear indication of the exciton-vibrational-photon coupling, which results in Rabi-like oscillations.As EC coupling strengths increase, the wave packet evolution retains oscillatory features, which, however, become progressively more erratic and spread over a wider range in the frequency domain. A detailed examination of the early stages of the spectral evolution is provided by the inset of Figure 3, which elucidates the wave packet dynamics influenced by the intrinsic and cavity-induced CIs.The spectral features primarily cluster around ω = 1.724 eV, aligning with the bright KK exciton energy.Notably, this observation confirms that the lower-energy KQ and KK ′ excitons, which are intrinsically dark, do not substantially contribute to the emission at these early times (cf.Ref. [82]).This scenario also illustrates the influence of polaritonic effects in Figure 3d: the fundamental excitonic KK state splits into two distinct polaritonic states, with a separation roughly equal to 2M σ − = 24 meV, thereby widening the spectrum.Additionally, interactions with phonon modes (polaron effects) further broaden the spectrum, extending it toward the blue side.Inversely, decreasing ECs causes a narrowing of S(ω, t) along the ω axis, which produces a more symmetric spectral shape and leads to a slower fluorescence quenching.For obtaining a detailed view of the spectral features, it is worthwhile to inspect S(ω, t) at specific values of t.TRF spectra at t = 0, which can be interpreted as absorption spectra [79], are shown in panel (a) of Figure 4, accompanied by relaxed fluorescence spectra at t = 1000 fs in panel (b).At t = 0, the spectra for all EC couplings exhibit similar shapes.Yet, larger EC coupling strengths reduce the intensity of the main peak associated with direct PL, as well as inducing a blue shift and spectral broadening.The red wing peak at ω ≈ E − E bright = −43 meV in Figure 4a is related to indirect PL, while the blue wing shoulder located around ω ≈ E − E bright = 30 meV represents a vibronic peak.This implies a modification (enhancement) of the Huang-Rhys factor [83], and a similar spectral redistribution is observed in surface plasmon polaritons [84].The observed enhancement of the vibronic peak intensity, as outlined in Ref. [85], stems from the closeto-resonance conditions.If the phonon frequency Ω α,q ∥ (25-30 meV) from the energy gap between the prominent KK peak and the adjacent higher energy peak matches the separation ∆ between the red wing peak and the central peak, the resonance condition is satisfied.This relationship is quantified as ∆ ≈ kΩ α,q ∥ , where k = 1, 2, . ... As EC couplings increase from 0 to 12 meV, ∆ nearly doubles the frequency of the (Λ and K) phonon modes.
The relaxed TRF spectrum at t = 1000 fs in Figure 4b shows drastic changes in lineshape as a function of the EC coupling strength.The intensity of the central direct PL peak (EC = 0) diminishes significantly due to the EC-enhanced population transfer to the dark excitons.In the strong coupling regime (EC = 8 meV and EC = 12 meV), the peak splits into two sub-peaks separated by 25 meV, as shown in the inset.These polaritonic features are not seen in the absorption spectrum.Therefore, TRF spectra deliver valuable information on the photoinduced processes in WSe 2 monolayers, which cannot be extracted from the absorption spectra.Figure 5 presents the TRF spectra S(ω, t) at intermediate times for four EC coupling strengths.For EC = 0 [panel (a)], the main dynamical events in the evolution of S(ω, t) are over before t = 250 fs (see Figure 3a).Hence, the spectra at 500 fs and 1000 fs in Figure 5a are almost identical, suggesting that the steady state fluorescence is reached quicker for weaker EC coupling.Only the intensity of the KK peak goes out slowly and remains much higher than that of the adjacent vibronic peaks.Additionally, the KK peak position stays unchanged in panels (a) and (b), contrasting with the significant shifts of the peak positions in panels (c) and (d).These shifts can be attributed to the polaritonic effects, which are manifested in significantly richer peak structures that spread over a substantially broader spectral domain.This indicates that the widths of relaxed TFR spectra can be used for estimating the EC coupling strengths.Qualitatively, the spectra in the strong EC coupling regime (panel d) are similar to those calculated in Ref. [31] using the Lindblad master equation for cavity-assisted TMDs monolayers.Nonetheless, the inter-peak separations and spectral widths in our spectra are significantly larger than those in Ref. [31] due to the stronger exciton-phonon and exciton-photon interactions.

Conclusions
By integrating the mDA method [41,42] with the TFD framework [49,50], we performed numerically accurate simulations of the exciton dynamics and TRF spectra of cavity-modulated WSe 2 monolayers at 75 K for EC coupling strengths ranging from 0 (no polaritonic effects) to 12 meV (strong polaritonic effects).The efficient computational method employed in the present work can be applied to simulations of various cavitytuned 2D materials at finite temperatures.More specifically, our approach can be extended to possibly account for the inversion of the band ordering (i.e., the lower polariton mode pushed below the WSe 2 dark state) and phonon-assisted polariton relaxation [86].Furthermore, anti-Stokes PL [87] activated in nanocavity-integrated WSe 2 monolayers can be simulated by incorporating resonant excitation of a dark exciton [87] in our Hamiltonian.
Our approach identified distinct spectroscopic signatures and characteristic timescales of polaronic and polaritonic effects in WSe 2 monolayers at various EC interactions.Notably, our investigation highlighted the crucial role of multidimensional CIs in governing the dynamics of strongly coupled excitons, photons, and phonons.It comes as no surprise that the EC coupling is an important parameter governing energy transfer in TMDs.What is more significant for applications and less evident is that we found that several aspects of the polaritonic transport in WSe 2 monolayers do not require finetuning of EC interactions.For example, increasing EC coupling strengths enhances the total population of the dark interval excitons KQ and KK ′ on the timescale of up to 200 fs, and it amplifies populations of the KK ′ excitons at times longer than 400 fs.This robustness in the cavity-mediated energy transport may be instrumental for engineering 2D materials.

Appendix A. Equations of Motion for the mDA Parameters
The time evolution of the parameters B mi (t), u ml (t), and ũmq (t) specifying the mDA wave function (10) are determined using the variational method [41,46,64]: where the Lagrangian L is defined as The multi-D 2 Ansatz converges to reproduce the accurate wave function if the multiplicity M is large enough.The multiplicity M = 48 adopted in our simulations yields the converged results.
The Hamiltonian structured within the framework of the multi-D 2 ansatz is expressed as il specifies the diagonal exciton-phonon coupling matrix elements for the exciton state |i⟩ interacting with phonon mode l.The terms D (2) l and D (3) l represent the excitonphonon couplings between KK exciton states and KQ (KK ′ ) exciton states with phonon mode l, respectively.Consequently, the governing equations for B m ′ 0 are formulated as B m ′ 1 , B m ′ 2 , and B m ′ 3 , which are formulated as Here µ i are the transition dipole moment coupling coefficients, where µ 0 = µ 1 = 1, µ 2 = µ 3 = 0. B m ′ i 1 i 0 (t) indicates the probability at time t for the exciton to occupy the state |i 1 ⟩ relative to the initial state |i 0 ⟩ with a given multiplicity M, while u m ′ qi 0 (t) and ũm ′ li 0 (t) denote displacements of the 'physical' and 'tilde' phonon modes.

Figure 1 .
Figure 1.A graphical representation of the direct and phonon-assisted PL emission processes.The momentum-dark excitons located at the K − Q and K − K ′ points undergo radiative decay via absorbing or emitting a phonon, thus contributing to the indirect PL emissions.The inset in the upper middle part illustrates the distribution of conduction electrons across various valleys within the first Brillouin zone.

Figure 4 .
Figure 4. TRF spectra S(ω, t) at t = 0 (a) and 1000 fs (b) for four EC coupling strengths, as indicated in the panels.