Microscopic modelling of exciton-polariton diffusion coefficients in atomically thin semiconductors

In the strong light-matter coupling regime realized e.g. by integrating semiconductors into optical microcavities, polaritons as new hybrid light-matter quasi-particles are formed. The corresponding change in the dispersion relation has a large impact on optics, dynamics and transport behaviour of semiconductors. In this work, we investigate the strong-coupling regime in hBN-encapsulated MoSe$_2$ monolayers focusing on exciton-polariton diffusion. Applying a microscopic approach based on the exciton density matrix formalism combined with the Hopfield approach, we predict a drastic increase of the diffusion coefficients by two to three orders of magnitude in the strong coupling regime. We explain this behaviour by the much larger polariton group velocity and suppressed polariton-phonon scattering channels with respect to the case of bare excitons. Our study contributes to a better microscopic understanding of polariton diffusion in atomically thin semiconductors.

In the strong light-matter coupling regime realized e.g. by integrating semiconductors into optical microcavities, polaritons as new hybrid light-matter quasi-particles are formed. The corresponding change in the dispersion relation has a large impact on optics, dynamics and transport behaviour of semiconductors. In this work, we investigate the strong-coupling regime in hBN-encapsulated MoSe2 monolayers focusing on exciton-polariton diffusion. Applying a microscopic approach based on the exciton density matrix formalism combined with the Hopfield approach, we predict a drastic increase of the diffusion coefficients by two to three orders of magnitude in the strong coupling regime. We explain this behaviour by the much larger polariton group velocity and suppressed polaritonphonon scattering channels with respect to the case of bare excitons. Our study contributes to a better microscopic understanding of polariton diffusion in atomically thin semiconductors.

I. INTRODUCTION
Semiconductors integrated into microcavities can be subject to a strong light-matter coupling [1,2]. In this regime, the coupling to light is stronger than the (difference of) lifetime of the involved particles, e.g. cavity modes and excitons [1, [3][4][5]. As a direct consequence, polaritons as dual light-matter quasi-particles are formed and considerably change the bandstructure, optics, dynamics and transport properties of semiconductors. A promising class of materials for the strong-coupling regime are monolayers of transition metal dichalcogenides (TMDs) [6][7][8][9][10]. They exhibit a large oscillator strength and exciton binding energies in the range of a few hundreds of meV governing the optoelectronic properties of these materials [11][12][13][14][15] and have been hence widely investigated in planar microcavities [6,9,16]. In the strong-coupling regime, the coupling of photons with confined excitons give rise to the formation of exciton-polaritons. The interaction leads to an avoided crossing at momenta where the cavity mode crosses the exciton dispersion resulting in a Rabi splitting, cf. Fig. 1(a) [16,17]. Hence, the polariton dispersion is characterized by two distinct energy branches denoted as the upper and the lower polariton (UP, LP). The weight of the photonic or excitonic character of polariton branches is expressed by the Hopfield coefficients C ± [18], cf. Fig. 1(b). * beatriz.ferreira@chalmers.se The quasi-bosonic polaritons show intriguing effects, such as Bose-Einstein condensation [10,19,20], (super)fluidity [21,22], topological effects [23], and promising applications from lasing [24], to integrated circuits [25] and quantum computing [26]. Exciton-polaritons maintain characteristics of both photons and excitons, such as a small effective mass, which makes them attractive also for transport purposes [27][28][29]. Previous studies have shown that TMD monolayers with their rich exciton landscape, including dark and bright exciton states [30][31][32], exhibit an interesting spatio-temporal exciton dynamics resulting in an intriguing exciton diffusion behaviour. This includes non-classical diffusion [33], transient negative diffusion [34], accelerated hot-exciton diffusion [35] or formation of spatial rings (halos) [36][37][38] and unconventional exciton funneling effects [39]. In view of their light component, polaritons show an interesting transport behaviour resulting in a fast propagation in the ballistic regime [27,29]. The polariton diffusion has already been observed e.g. in ZnSe and GaAs films [40,41]. To best of our knowledge, there have been no microscopic studies on exciton-polariton diffusion in atomically thin TMDs yet. In our work, we close this gap and present a microscopic study of polariton diffusion in an exemplary TMD monolayer integrated into an optical cavity. We investigate the change in the group velocity and polariton-phonon scattering rates as crucial ingredients determining the diffusion coefficient. Based on our microscopic approach, we predict a polariton diffusion coefficient up to three orders of magnitude higher compared to bare exciton diffusion.

II. THEORY
First, we describe the theoretical approach allowing us to microscopically investigate polariton diffusion coefficients in TMD monolayers. Exciton energies and wavefunctions in TMD monolayers are obtained by solving the Wannier equation [42][43][44] including DFT input on the characteristics of the electronic bandstructure [45]. In this work, we focus on an hBN-encapsulated MoSe 2 monolayer, which we find to be a direct semiconductor with the bright KK excitons (electron and hole located at the K point) as energetically lowest states [31]. In contrast, tungsten-based TMDs are known to be indirect semiconductors with momentum-dark excitons as energetically lowest states [31,46,47]. Since the latter are not directly affected by light-matter coupling, we expect smaller polariton-induced changes in the diffusion coefficient for tungsten-based materials. The starting point is the many-particle Hamilton operator in excitonic picture, which reads in second quantization [46,48] (1) Here,X † Q (X Q ),ĉ † Q (ĉ Q ), andb † q (b q ) are the exciton, photon, and phonon creation (annihilation) operators, respectively. The first two terms in the Hamiltonian describe the energy of excitons and photons with E X Q and E C Q , respectively. The third term corresponds to the exciton-light interaction mediated by the exciton-photon coupling matrix element g Q , where photons need to have the same in-plane momentum Q as excitons to fulfill the momentum conservation. In general, the out-of plane component influences the cavity energy and exciton-photon coupling. Here, we assume the existence of one resonant photon mode (i.e. E X 0 = E C 0 ) and we study its impact on the polariton diffusion [40]. Finally, the last term in the Hamiltonian describes the excitonphonon interaction [44], where the coupling strength is determined by the exciton-phonon matrix element D q . We use DFT input parameters for phonon dispersion and the strength of the electron-phonon coupling [49]. In particular, we consider three optical modes (LO, TO and A 1 ) with the energies E TO = 36.1 meV, E LO = 36.6 meV and E A1 = 30.3 meV, and for the two acoustic modes we consider a sound velocity of v s = 4.1 × 10 −3 nm/fs [49]. Now, we investigate the strong-coupling regime, where the exciton-photon coupling strength g Q is larger than (the difference of) cavity and nonradiative exciton decay rates [1]. This allows to form polaritons as new eigenmodes that can be obtained by applying a Hopfield transformation of the Hamilton operator yielding [1, 18] Here,Ŷ † (Ŷ ) denotes the polariton creation (annihilation) operator and E i Q the energy of the lower and upper polariton branch (i =LP,UP), cf. Fig.  1(a). In many cases, the decay rate of the cavity mode is much larger than the dissipation rate of excitons allowing us to find an analytic expression for the polariton energies For a vanishing excitonphonon coupling g Q , the lower/upper polariton branch approach the cavity and exciton energy (thin yellow and grey lines in Fig. 1(a)), respectively. Throughout this work, we focus on the resonant case, i.e. E X 0 = E C 0 . In the strong coupling regime with a large g Q , an avoided crossing occurs and two polariton branches are formed. Their separation corresponds to the Rabi splitting Ω R = 2g 0 . The two polariton branches can be visualized in optical spectra for large-enough coupling g [4,5,50].
The polariton states consist of a coherent mixture of excitons and photons with the in-plane momentum Q, thus the polariton operators can be expressed asŶ with the Hopfield coefficients h LP The second term in the Hamilton operator from Eq. (2) provides the polariton-phonon interaction. Since only the excitonic part of the polariton couples with phonons, the appearing matrix el-ementD is related to the exciton-phonon coupling [51]. Using the secondorder Born-Markov approximation [52,53], we determine the polariton-phonon scattering rate with the phonon momentum q = Q − Q, the Bose-Einstein distribution n q , and the phonon energy ω q . In our work, we partially include some effects beyond the so-called completed-collision limit [54] introducing a Lorentzian function instead of a Delta function. This is similar to the damping introduced in the second-order Born approximation including higher-order effects leading e.g. to a collisional broadening [55,56]. A microscopic calculation of the broadening γ is beyond the scope of this work. Through the manuscript, we use a value of 0.1 meV, which provides a low estimation of the presented scattering rates. As further discussed below, larger values of γ do not change the qualitative behaviour, but lead to a quantitative increase of the scattering of polaritons with acoustic modes due to higher-order non-resonant contributions.
Having determined the expression for polariton energies and polariton-phonon scattering rates, we have all necessary ingredients to investigate the polariton diffusion. Polaritons in general show a peculiar transient spatio-temporal dynamics resulting in a large ballistic propagation even at room temperature [29], as well as a non-linear transport behaviour [21,22]. While a full spatiotemporal polariton dynamics is beyond the scope of this work, we focus the investigation on polariton diffusion coefficients in steady-state limit. In TMD monolayers, excitons show -after an initial unconventional diffusion [34,35]-a regular steady-state diffusion behaviour, i.e. exhibiting a linear increase of the square width of the spatial distribution as a function of time [36,[57][58][59][60]. The rate of this increase is given by the diffusion coefficient D [35,61]. Such a regime appears when a local thermalized distribution is reached and when the scattering processes are fast enough leading to a quick thermalization compared to the transport timescale [61].
Assuming fast exciton-phonon scattering and quasithermalized local distributions [61], an analytic expression for the diffusion coefficient can be obtained with the polariton group velocity v i Q and occupation probability f i Q . The latter is assumed to be a thermalized Boltzmann distribution f i Q ∝ e −E i Q /k B T in the low-density limit (with k B as the Boltzmann constant and T as the temperature). Note that this approximation can in general lead to an overestimation of the actual occupation of quasi-photonic polariton states in the upper polariton branch. However, in the considered regime of resonant exciton and cavity energies, not too high temperatures and relatively large Rabi splittings, these effects are found to be small. According to Eq. (7), the crucial quantities determining the polariton diffusion are the polariton group velocity v i Q , the polariton-phonon scattering rate Γ i Q and the occupation of polariton states f i Q .
The many-particle mechanisms behind the diffusion can differ considerably when moving from TMD monolayers to TMD bulk materials. In the monolayer case, the reduced screening leads to large excitonic binding energies. As a consequence, the diffusion is typically dominated by excitons. Nevertheless, contributions from the faster diffusing electron-hole plasma can still appear for substrates with a large dielectric constant, as observed for hBNencapsulated TMDs at higher temperatures [62]. Bulk materials are expected to have smaller excitonic effects and thus higher diffusion coefficients in the range of 10 cm 2 /s, as observed e.g. for MoS 2 [63] and MoTe 2 [64]. In the next section, we discuss each of these microscopic quantities before we analyze the polariton diffusion coefficients and their temperature dependence.

A. Polariton group velocity and occupation
Now, we investigate the change of the excitonic band structure in the presence of a strong coupling regime. Figure 1(a) illustrates the polariton dispersion for three values of the Rabi splitting representing different exciton-photon coupling strengths. The latter depends on the oscillator strength of the material and the characteristics of the optical cavity. Exciton-photon coupling induces a Rabi splitting Ω R = E UP 0 − E LP 0 = 2g 0 and the formation of an upper and a lower polariton branch. We investigate the polariton dispersion for typical Rabi splitting values of Ω R =10, 25, 50 meV [9], which are larger than the non-radiative exciton linewidth (typically a few meV at low temperatures for hBN-encapsulation TMDs [44,65,66]) and the cavity linewidth (ranging from the meV [16] down to the µeV range [27,67]), thus allowing for strongcoupling regime [1].
At larger momenta, the LP and UP branches merge with the exciton and photon dispersion, respectively. The larger Ω R , the higher are the momentum values at which this occurs, cf. Fig.  1(a). Polaritons are coherent superpositions of excitonic and photonic states with the Hopfield coefficients defining the weights of the single constituents, cf. Fig. 1(b). The coefficient |C +,Q | 2 gives the exciton content of the lower polariton and the photon content of the upper polariton, i.e. for |C +,Q | 2 =1 the LP state |LP, Q =Ŷ LP † Q |0 coincides with the exciton state |X, Q and the UP state |UP, Q corresponds to the photon state |C, Q . In Fig. 2, we consider the polariton groupvelocity v i Q = −1 dE i Q /dQ, in particular focusing on the lower polariton branch, since the upper one has a negligible occupation and thus its impact on diffusion will be limited. We consider the two cases of Rabi splitting Ω R =25, 50 meV in comparison with the excitonic group velocity v X Q = Q/M X , where M X is the exciton mass. We find that the polariton group velocity is approximately 4 to 5 orders of magnitude larger than the excitonic one for small momenta within the light cone, cf. Fig. 2. Due to the rapidly changing polariton dispersion, we find group velocities in the range of 10 µm/ps, thus principally opening the possibility of ballistic polariton propagation for 10s µm. This has recently indeed been observed in a space-and angle-resolved photoluminescence experiments on a WS 2 monolayer in a distributed Bragg reflector cavity [29]. In addition to the remarkable magnitude difference, the group velocity for polaritons has also a qualitatively different momentum dependence. It shows a maximum in correspondence to the inflection point in the lower polariton branch and decreases toward the excitonic velocity for momenta of several µm −1 .
In a nutshell, two subset of states with a considerably different group velocity coexist in a cavity: The fast ones located within the light cone and the slow ones coinciding with conventional excitons. However, a large group velocity alone is not enough to boost the increase of diffusion coefficients, but these states also need to be occupied. This enters through the momentum-dependent Boltzmann distribution in Eq. (7). To illustrate this, we overlay the occupation of the lower polariton state on the line displaying its group velocity (reddish colours denote large occupation), cf. Fig. 2. While the excitonic occupation is momentum-independent in the considered range of momenta (cf. the thin orange line in Fig. 2), strong variations are observed for polaritons. At 20 K, the occupation of the states at larger momenta is decreased by two orders of magnitude with respect to the exciton case for both considered Rabi splittings (f LP /f X ≈ 10 −5 at Q ≈4µm −1 ), cf. orange vs black colour in Fig.  2. The curvature of the polariton branch induces a significant decrease of the occupation of the slow quasi-excitonic states at large momenta, as the energetically lower states at Q ≈0 are more efficiently populated. Increasing the temperature, the population of the former starts to increase, in particular for the smaller Rabi splitting of 25 meV, cf. Fig. 2(b). Regarding the behaviour at smaller momenta, we see that at 20 K the occupation of states with the maximum group velocity is negligible for Ω R =50 meV (black colour at approximately Q ≈ 1.3µm −1 ). However, when increasing the temperature to 40 K, we find a considerable occupation even at these states indicating the possibility of a strongly accelerated polariton diffusion.

B. Polariton-phonon scattering
Besides the group velocity and occupation of polariton states, polariton-phonon scattering plays an important role for the diffusion coefficient (cf. Eq. 7). Here, we neglect scattering with defects/disorder [29] or intervalley scattering with KK excitons. Figure 3 illustrates the polariton-phonon scattering rate at 40 K (cf. Eq. (6) for LP and UP branches around the light cone). Note that only the scattering into LP states out the light cone is efficient due to the limited number of receiving partner states available within the light cone as well as due to the negligibly small Hopfield coefficients h UP X for large-momenta UP states. This implies that the receiving LP state is quasi-excitonic, cf. Fig. 1, i.e. the associated coefficient fulfills h LP X ≈ 1 and the related scattering coefficient are proportional to  Fig. 3(b).
As a result, the scattering with phonons is driven by the excitonic component of the emitting polariton. One would naively expect larger scattering rates for LP states, as here the excitonic constituent is dominant reflected by h LP X 2 ≡ |C + | 2 ≈ 1. Surprisingly, our microscopic calculations of scattering rates show a much more efficient scattering for the UP branch, cf. Fig. 3. This can be traced back to the number of available scattering states fulfilling the momentum and energy conservation.
First, we discuss the LP scattering rates shown in Fig. 3(a). For states around Q=0, we find two orders of magnitude smaller polariton-phonon scattering compared to the exciton case (thin black line), while the Hopfield coefficient |C + | 2 = 0.5 would only imply a decrease by a factor of two. The reason for the dramatic decrease is related to the change in the dispersion relation in the strong coupling regime. The energy of acoustic phonons v s Q is almost flat compared to polaritons (cf. orange and red line in Fig. 3(c), respectively). As a consequence, when low-momentum polaritons absorb acoustic phonons, they are not able to find a resonant scattering partner. This is only possible if they are very close to the exciton energy E X 0 .The increase of the scattering rate at larger momenta can be traced back to non-resonant scattering with acoustic phonons (cf. the dashed line in Fig. 3(a) excluding acoustic phonons). This is due to the width of the Lorentzian in Eq. 6, whose origin can be related to higher-order scattering contributions inducing a softening of the energy selection rules. Note that the scattering rates show a quantitative dependence on this phenomenologically introduced width parameter, however the qualitative behaviour remains unaffected.
Next, we discuss the UP scattering rates illustrated in Fig. 3(b). At very small momenta, the scattering with acoustic modes is much more efficient and reaches a value that is approximately only two times smaller than for excitons, as expected by the Hopfield coefficient ( Fig. 1(b)). Since the UP branch is higher in energy compared to excitons, it is possible to find resonant scattering partners, cf. the crossing between the phonon dispersion (top orange line) and LP energy in Fig. 3 (c)). As a result, the scattering via acoustic phonons is efficient. At larger momenta, we observe the appearance of pronounced A and B resonances reflecting the emission of optical phonons. To better understand their origin, in Fig. 3(d) we plot the UP dispersion and the optical phonon energies (with respect to the exciton energy E X 0 ) and find crossing points exactly at the position of the A and B peaks in the scattering rate. Note that the different weight of these peaks for different Rabi splitting is due to the Hopfield coefficients. For ΩR=50 meV, the A peak appears at Q ≈ 1µm −1 , where |C + | 2 ≈ 0.4, while at 10 meV it appears at Q ≈ 1.5µ m −1 with |C + | 2 ≈ 0.04 resulting in a much smaller scattering efficiency. The UP scattering rates are dominated by resonant scattering, hence the width of the Lorentzian plays a minor role.
In a nutshell, the polariton-phonon scat- tering is strongly affected by the polariton dispersion resulting in suppressed scattering with acoustic phonons for LP and an enhanced emission of optical phonons for UP states. Note that the phononinduced scattering also contributes to the polariton linewidth in optical spectra [16,51]. Here, angleresolved spectroscopy reaching states with non-zero momenta [6] could principally allow to measure the predicted strongly momentum-dependent scattering rates via a change in the polariton linewidth, in particular for the UP branch.

C. Polariton diffusion
We have now discussed all key ingredients to evaluate and understand the polariton diffusion. Figure 4(a) shows the diffusion coefficient as a function of temperature and Rabi splitting. Based on our microscopic approach, we predict polariton dif-fusion coefficients that are two to three orders of magnitude larger than the ones from the bare exciton. This can be explained by: (i) the polariton dispersion exhibiting huge group velocities, (ii) effective occupation of fast polaritonic states, and (iii) reduced scattering with phonons of the occupied lower polariton states. These features concern, however, only small-momentum polaritons, while the diffusion coefficients depend also on large-momentum states. The latter are unaffected by point (i) and (iii), as the polariton dispersion and scattering rates correspond to the excitonic values at large momenta (Figs. 2,  3). Only the relative population remains affected by the presence of lower lying polariton states at small momenta. As a result, the polariton diffusion coefficient is the result of a non-trivial interplay between the few very fast states within the light-cone and excitonic-like states outside of the cone. Interestingly, for Ω R = 50 meV we observe a maximum in the polariton diffusion at around 40K. This can be traced back to the occupation discussed in Fig.  2: At 40 K, fast polariton states with a maximum group velocity at approximately Q 1.3µm −1 are efficiently populated resulting in a maximum diffusion. Further increasing the temperature occupies states at higher momenta (and a smaller group velocity) inducing a decrease of D.
In Fig. 4(b), we show the temperature dependence of the polariton diffusion coefficient at three fixed value for the Rabi splitting in comparison with the excitonic value (thin grey line). For increasing temperature, the polariton diffusion decreases towards the bare exciton diffusion (with D in the range of a few cm 2 /s) and this occurs faster for smaller Rabi splittings Ω R . At higher temperatures, the amount of occupied slower quasi-excitonic states outside the light cone becomes larger. Comparing the total diffusion coefficient with the contribution stemming only from the lower polariton states (thin lines), we find that the LP contribution is dominant for Ω R = 25 and 50 meV, as the UP states are only marginally occupied for the considered low temperatures. For the lower Rabi splitting of Ω R =10 meV, we find that the total and the LP diffusion start to deviate at higher temperatures indicating the increasing weight of the UP diffusion. To illustrate the impact of the polariton-phonon scattering on the diffusion, we calculate in a gedanken experiment the polariton diffusion assum-ing exciton scattering rates, cf. the dashed line in Fig.4(c). We observe a significant decrease in the polariton diffusion by more than one order of magnitude indicating the important role of polaritonphonon scattering. Note however that the polariton diffusion coefficient still remains considerably higher than the exciton one (thin grey line) reflecting the strong impact of the polariton group velocity. Interestingly, the peak at 40 K shown in the full polaritonic case disappears as the excitonic scattering rates decrease the distinction between fast and slow polariton states in view of their weak momentum dependence ( Fig. 3(a)).

IV. CONCLUSION
We have microscopically calculated the polariton diffusion coefficient for an hBN-encapsulated MoSe 2 monolayer embedded in an optical cavity. We predict a drastically enhanced polariton diffusion that is two to three orders of magnitude higher compared to excitons. We show that this accelerated polariton diffusion can be traced back to the increased group velocity of polaritons and the decreased scattering rates with phonons reflecting the dual light-matter character of polaritons.

ACKNOWLEDGMENTS
This project has received funding support from the DFG via SFB 1083 (project B9), the European Union's Horizon 2020 Research and Innovation programme under grant agreement no. 881603 (Graphene Flagship) and from the Knut and Alice Wallenberg Foundation via the Grant KAW 2019.0140. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC). We thank Jamie Fitzgerald (Chalmers) for valuable discussions that helped us in developing this work.