Phonon wind and drag of excitons in monolayer semiconductors

TL;DR

This paper develops a theoretical framework for understanding how non-equilibrium phonons influence exciton transport in monolayer semiconductors, revealing mechanisms like phonon wind and Seebeck effect that shape exciton distribution.

Contribution

It introduces a novel theory of exciton-phonon interactions under non-equilibrium conditions, explaining exciton drag and halo formation in monolayer materials.

Findings

Exciton drag by non-equilibrium phonons can create halo-like exciton distributions.

Ballistic and diffusive phonon regimes both contribute to exciton transport.

Phonon wind effects can be significant in localized excitation scenarios.

Abstract

We study theoretically the non-equilibrium exciton transport in monolayer transition metal dichalcogenides. We consider the situation where excitons interact with non-equilibrium phonons, e.g., under the conditions of localized excitation where a ``hot spot'' in formed. We develop the theory of the exciton drag by the phonons and analyze in detail the regimes of diffusive propagation of phonons and ballistic propagation of phonons where the phonon wind is formed. We demonstrate that a halo-like spatial distribution of excitons akin observed in [Phys. Rev. Lett. 120, 207401 (2018)] can be formed as a result of the exciton drag by non-equilibrium phonons or Seebeck effect.

Equations

$$\begin{gathered} \frac{\partial N_{\mathbf{q}}}{\partial t}+{\mathbf{V}}_{\mathbf{q}}\cdot {\nabla }_{\mathbf{r}}{N}_{\mathbf{q}}+\frac{N_{\mathbf{q}}-{\bar{N}}_q}{{\tau }_p^{ph}} \\ =-\frac{N_{\mathbf{q}}}{{\tau }_0}+Q_{ph-exc}\{N_{\mathbf{q}}\}+{\dot{N}}_{\mathbf{q}}(\mathbf{r},t). \end{gathered}$$

$$\begin{gathered} \frac{\partial f_{\mathbf{k}}}{\partial t}+{\mathbf{v}}_{\mathbf{k}}\cdot {\nabla }_{\mathbf{r}}{f}_{\mathbf{k}}+\frac{f_{\mathbf{k}}-{\bar{f}}_k}{{\tau }_p} \\ =-\frac{f_{\mathbf{k}}}{{\tau }_d}+Q_{exc-ph}\{f_{\mathbf{k}}\}+{\dot{f}}_{\mathbf{k}}(\mathbf{r},t). \end{gathered}$$

$$\begin{gathered} Q_{ph-exc}\{N_{\mathbf{q}}\} \\ =\sum _{\mathbf{k}}\left[W_{\mathbf{k}}^{em}(1+N_{\mathbf{q}})f_{\mathbf{k}}-W_{\mathbf{k}}^{abs}{N}_{\mathbf{q}}{f}_{\mathbf{k}}\right], \end{gathered}$$

$$W_{\mathbf{k}}^{em/abs}=\frac{2\pi }{\hslash }|M_q{|}^2\delta (E_{\mathbf{k}}-E_{\mathbf{k}-\mathbf{q}}\mp \hslash sq)$$

$$Q_{ph-exc}\{N_{\mathbf{q}}\}={\mathcal{W}}_{\mathbf{q}}{N}_{\mathbf{q}},$$

$${\mathcal{W}}_{\mathbf{q}}=\frac{2\pi |M_q{|}^2}{\hslash }\sum _{\mathbf{k}}\delta (E_{\mathbf{k}}-E_{\mathbf{k}-\mathbf{q}}-\hslash sq)(f_{\mathbf{k}}-f_{\mathbf{k}-\mathbf{q}}).$$

$${\mathcal{W}}_{\mathbf{q}}=\frac{2\pi |M_q{|}^2}{\hslash }\sum _{\mathbf{k}}\delta (E_{\mathbf{k}}-E_{\mathbf{k}-\mathbf{q}})(\hslash sq{\bar{f}}_k^{\prime }+\delta f_{\mathbf{k}}-\delta f_{\mathbf{k}-\mathbf{q}}),$$

$$\begin{gathered} Q_{exc-ph}\{f_{\mathbf{k}}\} \\ =\sum _{\mathbf{q}}{N}_{\mathbf{q}}\left[f_{\mathbf{k}+\mathbf{q}}{W}_{\mathbf{k}+\mathbf{q}}^{em}+f_{\mathbf{k}-\mathbf{q}}{W}_{\mathbf{k}-\mathbf{q}}^{abs}-f_{\mathbf{k}}(W_{\mathbf{k}}^{em}+W_{\mathbf{k}}^{abs})\right] \\ =\frac{2\pi }{\hslash }\sum _{\mathbf{q}}|M_q{|}^2(f_{\mathbf{k}+\mathbf{q}}-f_{\mathbf{k}})[N_{\mathbf{q}}\delta (E_{\mathbf{k}+\mathbf{q}}-E_{\mathbf{k}}-\hslash sq) \\ +N_{-\mathbf{q}}\delta (E_{\mathbf{k}+\mathbf{q}}-E_{\mathbf{k}}+\hslash sq)]. \end{gathered}$$

$$\begin{gathered} Q_{exc-ph}\{f_{\mathbf{k}}\} \\ =2\pi s{\bar{f}}_k^{\prime }\sum _{\mathbf{q}}q|M_q{|}^2(N_{\mathbf{q}}-N_{-\mathbf{q}})\delta (E_{\mathbf{k}+\mathbf{q}}-E_{\mathbf{k}}). \end{gathered}$$

$$N_{\mathbf{q}}(\mathbf{r},t)=\int d{\mathbf{r}}^{\prime }\int d{t}^{\prime }{\mathcal{G}}_{\mathbf{q}}(\mathbf{r}-{\mathbf{r}}^{\prime },t-t^{\prime }){\dot{N}}_{\mathbf{q}}({\mathbf{r}}^{\prime },t^{\prime }),$$

$$\left[-\mathrm{i}\omega +\frac{1}{\tau }+\mathrm{i}({\mathbf{V}}_{\mathbf{q}}\cdot \mathbf{Q})\right]{\mathcal{G}}_{\mathbf{q}}(\mathbf{Q},\omega )=\frac{{\bar{\mathcal{G}}}_q(\mathbf{Q},\omega )}{{\tau }_p^{ph}}+1$$

$${\mathcal{G}}_{\mathbf{q}}(\mathbf{Q},\omega )=\tau \frac{{\bar{\mathcal{G}}}_q(\mathbf{Q},\omega )/{\tau }_p^{ph}+1}{1-\mathrm{i}\omega \tau +iQl\cos \phi },$$

$${\bar{\mathcal{G}}}_{\mathbf{q}}(\mathbf{Q},\omega )=\frac{P_{Q,\omega }\tau }{1-P_{Q,\omega }\frac{\tau }{{\tau }_p^{ph}}},P_{Q,\omega }=\frac{1}{\sqrt{(1-\mathrm{i}\omega \tau {)}^2+Q^2{l}^2}}.$$

$$\begin{gathered} {\mathcal{G}}_{\mathbf{q}}(\mathbf{Q},\omega )= \\ \frac{\tau \sqrt{(1-\mathrm{i}\omega \tau {)}^2+Q^2{l}^2}}{(\sqrt{(1-\mathrm{i}\omega \tau {)}^2+Q^2{l}^2}-\tau /{\tau }_p^{ph})(1-\mathrm{i}\omega \tau +\mathrm{i}Ql\cos \phi )}. \end{gathered}$$

$${\dot{N}}_{\mathbf{q}}(\mathbf{r},t)=\pi r_0^2{\mathcal{N}}_q\delta (\mathbf{r})\delta (t).$$

$$\begin{gathered} N_{\mathbf{q}}(\mathbf{r},t)=\pi r_0^2{\mathcal{N}}_q{\mathcal{G}}_{\mathbf{q}}(\mathbf{r},t) \\ =\pi r_0^2{\mathcal{N}}_q\int \frac{d^{2}Q}{(2\pi {)}^2}\int \frac{d\omega }{2\pi }{\mathcal{G}}_{\mathbf{q}}(\mathbf{Q},\omega ){\text{e}}^{i(\mathbf{Q}\cdot \mathbf{r}-\omega t)}. \end{gathered}$$

$$N_{\mathbf{q}}^{\text{ball}}(\mathbf{r},t)=\pi r_0^2{\mathcal{N}}_q\delta \left(\mathbf{r}-st\frac{\mathbf{q}}{q}\right){\text{e}}^{-t/{\tau }_0}.$$

$${\mathcal{G}}_{\mathbf{q}}(\mathbf{Q},\omega )\approx \frac{1-i{\mathbf{V}}_{\mathbf{q}}\cdot \mathbf{Q}{\tau }_p^{ph}}{{\tau }_0^{-1}-i\omega +{\mathcal{D}}_{ph}{Q}^2},$$

$${\mathcal{D}}_{ph}=\frac{l^2}{2\tau }=\frac{s^2{\tau }_p^{ph}}{2}.$$

$$N_{\mathbf{q}}^{\text{diff}}(\mathbf{r},t)=\left(1-{\tau }_p^{ph}{\mathbf{V}}_{\mathbf{q}}\cdot {\nabla }_{\mathbf{r}}\right){\mathcal{N}}_q{r}_0^2\frac{{\text{e}}^{-r^2/4{\mathcal{D}}_{ph}t-t/{\tau }_0}}{4{\mathcal{D}}_{ph}t}.$$

$$N_{\mathbf{q}}^{\text{diff}}(\mathbf{r},t)=\left(1-{\tau }_p^{ph}{\mathbf{V}}_{\mathbf{q}}\cdot {\nabla }_{\mathbf{r}}\right)\frac{k_{B}T(\mathbf{r},t)}{\hslash sq},$$

$$\frac{\partial T(\mathbf{r},t)}{\partial t}={\mathcal{D}}_{ph}\Delta T(\mathbf{r},t)-\frac{T(\mathbf{r},t)}{{\tau }_0}+\dot{T}(\mathbf{r},t),$$

$$Q_{exc-ph}\{f_{\mathbf{k}}\}=-{\bar{f}}_k^{\prime }({\mathbf{v}}_{\mathbf{k}}\cdot {\mathcal{F}}_{\mathbf{k}}),$$

$$\begin{gathered} ({\mathbf{v}}_{\mathbf{k}}\cdot {\mathcal{F}}_{\mathbf{k}})=-\frac{ms}{\pi {\hslash }^2}\underset{0}{\overset{2\pi }{\int }}d{\phi }_{\mathbf{q}}q|M_q{|}^2\Theta (2k-q) \\ \times (N_{\mathbf{q}}-N_{-\mathbf{q}}). \end{gathered}$$

$$\delta (E_{\mathbf{k}+\mathbf{q}}-E_{\mathbf{k}})=\frac{2m}{{\hslash }^{2}q}\delta [q+2k\cos ({\phi }_{\mathbf{q}}-{\phi }_{\mathbf{k}})].$$

$$n=g\int \frac{d\mathbf{k}}{(2\pi {)}^2}{f}_{\mathbf{k}},\mathbf{j}=g\int \frac{d\mathbf{k}}{(2\pi {)}^2}{\mathbf{v}}_{\mathbf{k}}{f}_{\mathbf{k}}.$$

$$\frac{\partial n}{\partial t}+\mathrm{div}\mathbf{j}+\frac{n}{{\tau }_d}=\dot{n}(\mathbf{r},t),$$

$$\dot{n}(\mathbf{r},t)=g\int \frac{d\mathbf{k}}{(2\pi {)}^2}{\dot{f}}_{\mathbf{k}}(\mathbf{r},t)$$

$${\bar{f}}_k\propto \frac{n(\mathbf{r},t)}{T_x(\mathbf{r},t)}\exp \left(-\frac{E_k}{T_x(\mathbf{r},t)}\right),$$

$$\mathbf{j}=-D_x{\nabla }_{\mathbf{r}}n-\eta n{\nabla }_{\mathbf{r}}{T}_x+\frac{{\tau }_p}{m}\mathbf{F}(\mathbf{r},t)n.$$

Full text

I Introduction

Two-dimensional semiconductors based on transition metal dichalcogenides form a family of materials whose optical properties are dominated by excitons, electron-hole pairs bound by the Coulomb interaction Mak et al. (2010); Splendiani et al. (2010); Chernikov et al. (2014); Wang et al. (2018); Durnev and Glazov (2018). Sharp excitonic resonances with high oscillator strengths at the room temperature open up wide prospects for light-matter coupling studies and applications in the fields of photonics and polaritonics Schneider et al. (2018).

Recently, an increased interest to the transport properties of excitons has appeared Mouri et al. (2014); Yuan et al. (2017); Cadiz et al. (2018); Kulig et al. (2018). It is motivated by the possibility to access fundamental parameters of exciton dynamics such as momentum scattering time and diffusion coefficient as well as by the possibilities of applications of two-dimensional materials for the energy harvesting. In recent work Kulig et al. (2018) different regimes of exciton transport have been revealed in monolayer WS2: At low exciton densities the classical diffusion of excitons has been observed. An increase in the density of excitons enables the Auger process Abakumov et al. (1991); Sun et al. (2014); Kumar et al. (2014); Mouri et al. (2014); Moody et al. (2016); Manca et al. (2017); Han et al. (2018) giving rise to an efficient non-radiative exciton-exciton annihilation. It changes the shape of the exciton distribution profile in the real space and gives rise to an apparent increase of the observed diffusion coefficient. Finally, at even higher exciton densities the Gaussian-like distribution of the excitons evolves into long-lived halo shaped profiles with the dip in the middle. The effect is attributed to the memory effects resulting, e.g., from the heating of excitons due to the efficient Auger process Kulig et al. (2018).

The energy of overheated excitons can be dissipated by emission of phonons. It gives rise to a formation of a “hot spot” in the sample where the average kinetic energy of some of the phonon modes exceeds by far the lattice temperature outside of spot. These nonequilibrium phonons start to propagate out of the hot spot bringing away extra energy and momentum. As a result of the exciton-phonon interaction the excitons tend to drift away from the excitation spot as well. In this scenario phonons can drag excitons just like non-equilibrium phonons push electrons in metals and semiconductors Gurevich (1946), see Ref. Gurevich and Mashkevich (1989) for review. Furthermore, if phonons propagate ballistically, the phonon “wind” is formed Keldysh (1976); Zinov’ev et al. (1983); Kozub (1988); Bulatov and Tikhodeev (1992) resulting in unusual exciton propagation regimes somewhat resembling superfluid flow Link and Baym (1992); Bulatov and Tikhodeev (1992); Kopelevich et al. (1996); Tikhodeev et al. (1998).

While the phonon drag of excitons has been studied in detail for bulk materials and conventional quasi-two-dimensional systems, its specific features in atomically-thin transition metal dichalcogenides have not been explored. In particular, while in bulk semiconductors and quantum wells the excitons mainly interact with three-dimensional phonons because a weak acoustic contrast is insufficient to strongly confine the phonons in the quantum well plane, in monolayer-thin materials and in van der Waals heterostructures the phonons are strongly confined in the two-dimensional layer. It results in the enhancement of the phonon density of states. This together with the fact that the exciton radii in these monolayer materials are quite small results in rather strong the exciton-phonon interaction in two-dimensional materials Slobodeniuk and Basko (2016); Vuong et al. (2017); Selig et al. (2016); Christiansen et al. (2017); Shree et al. (2018); Song:2013uq ; PhysRevB.92.125431 ; PhysRevLett.122.217401 .

Here we develop the kinetic theory of the exciton transport in two-dimensional transition metal dichalcogenides in the case where excitons interact with non-equilibrium phonons. In Sec. II we formulate the problem and present the kinetic equations for the exciton and phonon distribution functions. Next, in Sec. III we reduce the kinetic equations to an effective drift-diffusion model for excitons. Section IV contains analytical and numerical results on exciton propagation under the conditions of phonon drag and phonon wind. Brief conclusion is presented in Sec. V

II Kinetic theory

II.1 Model

We consider a monolayer of transition metal dichalcogenide semiconductor and assume that at the initial time moment $t=0$ it is excited by a short and tightly-focused optical pulse with the photon energy exceeding the exciton resonance energy similarly to the experimental situation of Ref. Kulig et al. (2018). The absorption of the pulse results in the formation of electron-hole pairs and excitons which lose energy by phonon emission. We assume that initially the exciton density is high enough so that the Auger-like exciton-exciton annihilation process takes place as well. It results in the non-radiative “bi-molecular” recombination where one exciton recombines with its energy and momentum being transferred to another exciton which ends up in a highly-excited state Takeshima (1975); Kavoulakis and Baym (1996); Wang et al. (2006); Sun et al. (2014); Kulig et al. (2018); Manca et al. (2017); Han et al. (2018). It starts losing energy emitting more phonons. Importantly, high energy phonons with almost flat dispersion that can include both optical phonons at the Brillouin zone center and the zone-edge modes serve as a major source of energy loss of non-equilibrium quasi-particles Gantmakher and Levinson (1987). As a result, their velocity is close to zero and such phonons accumulate. They form a hot spot in monolayer semiconductor. Eventually, owing to lattice anharmonicity optical phonons and zone edge acoustic phonons decay to the acoustic phonons with small wavevectors and almost linear dispersion. The latter propagate out of the hot spot bringing away energy and momentum and produce the drag of the excitons. Another potential experimental option is to generate phonons by additional intense laser beam or electric current pulse, in which case one can control phonon distribution independently Moskalenko et al. (1994, 1995); Akimov et al. (2006). Schematically, the system under study is shown in Fig. 1(a) and the process of momentum transfer from non-equilibrium acoustic phonons with linear dispersion to excitons is shown in Fig. 1(b)

Thus, in order to describe the exciton dragging by non-equilibrium phonons it is sufficient to consider the acoustic phonons only taking into account the processes described above as a source of non-equilibrium acoustic phonon population. In what follows we describe the acoustic phonons by the distribution function $N_{\bm{q}}(\bm{r},t)$, where $\bm{q}$ is the wavevector of the phonon, $\bm{r}$ is the position, and $t$ is time. Similarly, the excitons are described by the distribution function $f_{\bm{k}}(\bm{r},t)$ with $\bm{k}$ and $\bm{r}$ being the exciton wavevector and its coordinate, respectively. Below in Sec. II.2 we present the set of coupled kinetic equations for the exciton and phonon distribution functions.

II.2 Kinetic equations

The kinetic equation for the phonon distribution $N_{\bm{q}}(\bm{r},t)$ describing phonon redistribution in the real and momentum space has the form

[TABLE]

Here ${\bm{V}}_{\bm{q}}=s\bm{q}/q$ is the phonon velocity in the state with the wavevector $\bm{q}$ with $s$ being the speed of sound, $\tau_{0}$ is the phonon lifetime, $\tau_{p}^{ph}$ is isotropization time of the phonon distribution function (phonon momentum relaxation time), and the overline denotes averaging over the directions of $\bm{q}$. The right hand side of Eq. (1) describes the processes where the total number of phonons changes, ${Q_{ph-exc}\{N_{\bm{q}}\}}$ is the phonon-exciton collision integral and ${\dot{N}_{\bm{q}}}(\bm{r},t)$ is the generation rate in the course of the exciton cooling and Auger recombination. Here and in what follows we distinguish between the processes of acoustic phonon generation in the course of exciton energy relaxation via optical and zone-edge phonons, described by the term $\dot{N}$, and the interaction between the excitons and low-energy acoustic phonons, described by the term $Q_{ph-exc}$. Equation (1) can be derived from the continuity equation for the distribution function taking into account the incoherent processes (phonon generation, annihilation and scattering) as collision integrals, whose form is specified below.

The kinetic equation for the exciton distribution function has a form similar to Eq. (1)

[TABLE]

Here ${\dot{f}_{\bm{k}}}(\bm{r},t)$ is the exciton generation rate due to optical excitation, $\bm{v}_{\bm{k}}=\hbar\bm{k}/m$ is the exciton velocity, $m$ is its translational motion mass, $\tau_{p}$ and $\tau_{d}$ are the momentum relaxation and population decay time (lifetime) of exciton, respectively, ${Q_{exc-ph}\{f_{\bm{k}}\}}$ is the exciton-phonon collision integral. Being interested in the exciton transport in the presence of non-equilibrium phonons we also disregard any nonlinear effects resulting from exciton-exciton interaction such as exciton-exciton scattering, energy renormalization as well as the Auger recombination, see Sec. IV for brief discussion.

Kinetic equations (1) and (2) are coupled via the collision integrals ${Q_{ph-exc}\{N_{\bm{q}}\}}$ and ${Q_{exc-ph}\{f_{\bm{k}}\}}$. These integrals describe the processes of phonon emission and absorption due to exciton-phonon interaction. Hereafter we assume that the excitons are non-degenerate, $f_{\bm{k}}\ll 1$ and the collision integral operator ${Q_{ph-exc}\{N_{\bm{q}}\}}$ acting on the phonon distribution function has a form

[TABLE]

where the rates of the phonon emission, $W^{em}_{\bm{k}}$, and absorption, $W^{abs}_{\bm{k}}$, are given by the Fermi golden rule

[TABLE]

with $M_{q}$ being the matrix element of the exciton-phonon interaction Shree et al. (2018) and $E_{\bm{k}}=\hbar^{2}k^{2}/2m$ is the exciton dispersion. For simplicity and better clarity we disregard here the complex band structure involving several valleys and also spin degrees of freedom: Given similar masses and phonon coupling strengths across different valleys, it is further justified by inefficient intervalley and spin-flip scattering processes with low energy and momentum acoustic phonons.

In what follows we assume that the phonon occupancies are high, $N_{\bm{q}}\gg 1$, in agreement with the situation we are describing. This condition is fulfilled since the energies of acoustic phonons $\hbar sq$ interacting with excitons is typically much smaller than exciton energy Gantmakher and Levinson (1987); Shree et al. (2018), thus, under experimental conditions it is much smaller than the temperature both of the lattice and of the excitons. In this regime it is sufficient to consider only stimulated processes and the collision integral can be further simplified as

[TABLE]

where

[TABLE]

The exciton interaction with acoustic phonons is just weakly inelastic, see Refs. Gantmakher and Levinson (1987); Shree et al. (2018) for discussion. Thus, we have

[TABLE]

where $\bar{f}_{k}^{\prime}=d\bar{f}/dE_{k}$. Equations (5) and (6) demonstrate that exciton-phonon collisions induce (generally anisotropic) correction to the phonon lifetime, it will be ignored in what follows.

Under the same assumptions as above (i.e., low exciton occupancy, $f_{\bm{k}}\ll 1$, and high phonon occupancy, $N_{\bm{q}}\gg 1$) we arrive at the following form of the exciton-phonon collision integral:

[TABLE]

In the elastic approximation, where $\hbar sq$ in the energy conservation law is omitted, all the effect of exciton-phonon collisions is reduced to the isotropization of the exciton distribution function. In what follows we assume that this contribution to $\tau_{p}^{-1}$ is already included in the momentum relaxation rate in the kinetic equation (2).

Thus, in order to describe the exciton drag by the phonons we need to take into account non-elasticity in Eq. (7) in the lowest non-vanishing order. To that end we arrive at the collision integral in the simple form

[TABLE]

In derivation of Eq. (8) we assumed that the anisotropic in the momentum space part of the exciton distribution function is much smaller that the isotropic one and neglected it. We stress that the right hand side of the collision integral $Q_{exc-ph}\{f_{\bm{k}}\}$ Eq. (8) is non-zero due to the fact that the phonons propagate out of the hot spot, Fig. 1. Correspondingly, the spatial gradient of $N_{q}$ produces the anisotropic part of the phonon distribution function $N_{\bm{q}}-N_{-\bm{q}}\neq 0$.

In the following section we reduce the kinetic equations to the drift-diffusion equation for excitons which is suitable for description of the phonon drag.

III Exciton drift-diffusion model

III.1 Phonon distribution function

Hereafter we assume that propagating excitons provide minor effect on phonon distribution function. This is reasonable because the phonons are generated while the major part of the exciton population relaxes in energy and decays due to the Auger process. The collision integral (5) is proportional to the number of excitons and is small. Consequently, we neglect the $Q_{ph-exc}$ in the right hand side of Eq. (1). Hence, the solution of this equation is expressed via the Green’s function $\mathcal{G}_{\bm{q}}(\bm{r},t)$ as

[TABLE]

The Green’s function Fourier-transform satisfies the equation

[TABLE]

with $1/\tau=1/\tau_{0}+1/\tau_{p}^{ph}$. Its solution reads

[TABLE]

where $l=s\tau$, $\varphi$ is the angle between $\bm{Q}$ and $\bm{q}$, and the angular-average Green’s function can be readily expressed as

[TABLE]

This yields the following closed form expression for the Green’s function:

[TABLE]

We assume that at $t=0$ the ${\cal N}_{q}$ phonons with the wavevector absolute value $q$ were generated in the small spot of the area $\pi r_{0}^{2}$ at $\bm{r}=0$, Fig. 1(a). Accordingly, we take the generation rate ${\dot{N}_{\bm{q}}}(\bm{r},t)$ in the form 111Generalization of the results to the large hot spot size is trivial and beyond the scope of the present paper. Basic results remain the same, see Appendix .

[TABLE]

In this case we get:

[TABLE]

In what follows we consider two important limits for phonon propagation. First one is the case of ballistic propagation where the momentum relaxation of phonons in unimportant, $\tau_{0}/\tau_{p}^{ph}\ll 1$. In this regime we have $\mathcal{G}_{\bm{q}}(\bm{Q},\omega)\approx(\tau_{0}^{-1}-\mathrm{i}\omega+iQs\cos{\varphi})^{-1}$, and

[TABLE]

Equation (15) clearly shows that the phonons propagate along straight lines from the origin to the point $\bm{r}$ with the fixed speed $s$ and $\bm{q}\parallel\bm{r}$.

Second important regime of phonon transport is the diffusive propagation, where the momentum scattering time $\tau_{p}^{ph}\ll\tau_{0}$. In the diffusion regime there are three small parameters: $\omega\tau_{p}^{ph},Ql,\tau_{p}^{ph}/\tau_{0}\ll 1$ and the distribution of phonons is practically isotropic in the momentum space. Since we are interested in the transfer of phonon momentum to excitons we need to account for its small anisotropy in the first order. As a result we have for the phonons Green’s function

[TABLE]

where we introduced the phonon diffusion coefficient

[TABLE]

We stress that Eqs. (16) and (17) are valid in the regime of slow, diffusive, dynamics of phonon distribution, i.e., on the frequency scales $\omega\ll 1/\tau_{p}^{ph}$ and on the wavevector scales $Q\ll 1/l$, while the product $\omega\tau_{0}$ can be arbitrary. That is why it is sufficient to take into account only the lowest non-vanishing terms in $\omega$ and $Q$ in the denominator and recover the diffusion pole in the form $(-i\omega+{\mathcal{D}_{ph}}Q^{2})^{-1}$. In this regime the short time dynamics of the phonon distribution function can be included by appropriate modification of the initial condition (13). Finally, for the phonon distribution we have

[TABLE]

It is possible to recast Eq. (18) in somewhat different notations which are convenient for what follows. To that end we represent the wavevector average phonon distribution as $\bar{N}_{\bm{q}}^{\text{diff}}(\bm{r},t)=k_{B}T(\bm{r},t)/(\hbar sq)$, where $k_{B}$ is the Boltzmann constant and $T(\bm{r})$ is the local (i.e., coordinate dependent) temperature of the small energy acoustic phonon subsystem 222Such description is valid for thermalized phonons which form Planck distribution due to anharmonic processes. Also we assume that $k_{B}T\gg\hbar sr_{0}^{-1}$, where $r_{0}$ is the hot spot size.. Hence,

[TABLE]

where $T(\bm{r})$ satisfies standard heat conduction equation

[TABLE]

with $\dot{T}(\bm{r},t)$ being the energy generation rate in the units of $k_{B}^{-1}$. Equations (19) and (20) allow one to calculate phonon flux for an intense local heating of the sample. We emphasize here that $T(\bm{r})$ is not necessary equals to the local lattice temperature and may not represent the heating of the material in the conventional sense.

III.2 Drift-diffusion equation for excitons

Let us now analyze the impact of exciton-phonon interaction on the exciton distribution. To that end we note that the collision integral (8) can be recast as

[TABLE]

where the scalar product of the exciton velocity and the effective force can be presented in the form

[TABLE]

Here the phonon wavevector absolute value in the integral and the angle $\varphi_{\bm{q}}$ are interrelated as $q=-2k\cos{(\varphi_{\bm{q}}-\varphi_{\bm{k}})}$ by virtue of the energy conservation law

[TABLE]

In what follows we assume that the excitons propagate diffusively in the sample, i.e., we are interested in the exciton dynamics on the time scales which exceed by far the exciton momentum scattering time $\tau_{p}$. Therefore, under assumption $\tau_{d}\gg\tau_{p}$, the kinetic equation (2) with the collision integral in the form of Eq. (21) can be replaced by the effective drift diffusion equations. In order to derive these equations we introduce the exciton density $n(\bm{r},t)$ and the exciton flux density $\bm{j}(\bm{r},t)$ according to

[TABLE]

Here $g$ is the spin and valley degeneracy factor which takes into account the band structure of the transition metal dichalcogenides monolayers (within simplest approximation where the intervalley excitons are neglected, $g=4$ or $2$ depending whether both bright and spin dark excitons or only one species are involved). Naturally, the integration in Eqs. (23) picks the zeroth and first angular harmonics of the exciton distribution which determine, respectively, the particle density and flux.

Integrating Eq. (2) over the wavevectors and taking into account the exciton-phonon interaction does not change the number of excitons, we arrive at the continuity equation for the exciton density

[TABLE]

where

[TABLE]

is the exciton generation rate. The equation for the flux density $\bm{j}$ can be derived in the same way by multiplying Eq. (2) by $\bm{v}_{\bm{k}}$ and integrating over $\bm{k}$. In order to simplify the resulting expressions we assume that $\tau_{p}$ is independent of exciton energy and that the excitons on average have a Boltzmann distribution

[TABLE]

with the effective exciton temperature $T_{x}(\bm{r},t)$ which, in general, is not equal to the acoustic phonon subsystem temperature $T(\bm{r},t)$ and depends on the position and time. Calculation of the thermalization rate is beyond the scope of the paper. It can be carried out following Refs. Ivchenko and Takunov (1988); Golub et al. (1996); Selig et al. (2018). As a result we have

[TABLE]

Here the first term is responsible for the exciton diffusion with

[TABLE]

being the exciton diffusion coefficient, the second term accounts for the Seebeck effect Ashcroft and Mermin (1976); Askerov (1994), see also Ref. 2019arXiv190602084P for advanced calculations of the Seebeck coefficient for excitons, with

[TABLE]

and the last term accounts for the force $\bm{F}(\bm{r},t)$ acting from the phonons on the excitons (phonon wind and drag effects),

[TABLE]

Note that if $\bm{\mathcal{F}}_{\bm{k}}$ is wavevector independent then $\bm{F}(\bm{r},t)=\bm{\mathcal{F}}(\bm{r},t)$. The set of Eqs. (24) and (27) can be combined into the single drift-diffusion equation for the exciton density

[TABLE]

The first term in parentheses describes the exciton diffusion and two remaining terms describe the drift due to the Seebeck effect and interaction with phonons, respectively. Equation (31) can be recast in a more compact form using the explicit expressions for $D_{x}$ and $\eta$, Eqs. (28) and (29):

[TABLE]

Equation (31) should be also supplemented by a similar equation for the exciton temperature $T_{x}(\bm{r},t)$.

The model formulated above describes the exciton diffusion and drift due the Seebeck effect on overheated excitons and exciton-phonon interaction. In fact, similar description is valid if excitons are interacting with other particles or quasi-particles in the two-dimensional crystal. For example, the excitons can be dragged by the fluxes of electrons and holes, which can be generated in the sample, e.g., due to the hot exciton dissociation. In this case, $\bm{F}$ has a meaning of the momentum transfer rate from the electrons and holes to the excitons. We also note that non-linear effects due to exciton-exciton interactions, such as Auger recombination and renormalization of the exciton energy can be straightforwardly introduced into Eq. (31). Similar drift-diffusion description of excitons can be also applicable provided that exciton-exciton collisions Note:xx are more efficient as compared with the exciton-phonon and exciton-impurity scattering by analogy with the hydrodynamical regime of electron transport, see Refs. Gurzhi_1968 ; PhysRevLett.117.166601 and references therein.

IV Results

In this section we present the results of the exciton propagation modelling in the framework of the drift-diffusion model, Eq. (31). First, we determine the effective force $\bm{F}(\bm{r},t)$ acting from the phonons on the excitons and, second, we present the results of analytical and numerical calculations.

IV.1 Phonon-induced driving forces

The exciton-phonon coupling in transition metal dichalcogenide monolayers is mainly governed by the deformation potential interaction Christiansen et al. (2017); Shree et al. (2018). In the long wavelength limit relevant for our work, where the product $qa_{B}\ll 1$, with $a_{B}$ being the exciton Bohr radius, the phonon-induced deformation of the crystalline lattice is homogeneous on the scale of the exciton size. Thus, the energy shift of the exciton state is simply provided by the variation of the band gap energy and only longitudinal acoustic phonons contribute to the effect. Accordingly, we present the matrix element $M_{q}$ in Eqs. (4) and (8) as

[TABLE]

where $\Xi$ is a constant related to the conduction band, $D_{c}$, and valence band, $D_{v}$, deformation potentials as Shree et al. (2018)

[TABLE]

with $\mathcal{S}$ being the normalization area and $\varrho$ being the mass density of the two-dimensional crystal.

For diffusive phonon propagation where $\tau_{p}^{ph}\ll\tau_{0}$, Eq. (19) yields

[TABLE]

and by virtue of Eqs. (21), (22) and (30) we have

[TABLE]

In the last equation we introduced the characteristic exciton-phonon scattering time as

[TABLE]

Equation (35) has a clear physical sense: $-k_{B}\bm{\nabla}T$ is the momentum per unit of time carried by the phonons and the factor $\tau/\tau_{x}$ (note that in diffusive regime $\tau=\tau_{p}^{ph}$) gives the fraction of phonon momentum transferred to exciton in the course of exciton-phonon interaction. Such regime of the exciton drift where phonons are diffusive is denoted as the phonon drag regime.

It is instructive to compare the phonon drag force in Eq. (35) and the effective force due to the Seebeck effect caused by the exciton temperature gradient, second term in Eq. (27), see also second term in the brackets in Eq. (31). Under our assumptions the corresponding exciton Seebeck force is simply given by

[TABLE]

Provided that the temperature gradients of excitons and phonons are the same, the ratio $|\bm{F}_{\rm drag}|/|\bm{F}_{S}|=\tau/\tau_{x}<1$, because total scattering rate of phonons is larger than their rate of collisions with excitons. Generally, the temperatures of excitons and phonons as well as their gradients can differ, particularly, in the situation where the phonons are generated by additional light pulse. The study of competition between the Seebeck effect and phonon drag is beyond the scope of this work, we just stress that both effects produce additive contributions to the exciton flux.

We recall that according to Eqs. (18) and (19) the temperature distribution under pulsed excitation has a Gaussian form and the phonon drag force in Eq. (35) acquires a form

[TABLE]

We recall that $T_{0}$ is the effective temperature in the hot spot and $r_{0}$ is the hot spot radius, see Fig. 1. It is assumed that $r_{0}\ll r$, $\sqrt{\mathcal{D}_{ph}t}$ and $\sqrt{D_{x}t}$ at the relevant space and time scales.

For ballistic phonon propagation where $\tau_{p}^{ph}\gg\tau_{0}$, Eq. (15) yields

[TABLE]

Further, to obtain simple analytical expressions, we again approximate $\mathcal{N}_{q}$ as $\mathcal{N}_{q}=k_{B}T_{0}/\hbar sq$ with $T_{0}$ being the effective temperature of the hot spot and $r_{0}$ being its radius, Fig. 1(a). As a result, we obtain

[TABLE]

It is instructive to generalize Eq. (40) to a situation where the phonons are generated in the hot spot during a finite time $\tau_{\epsilon}$. To that end we replace $\delta(t)$ by $\tau_{\epsilon}^{-1}\exp{(-t/\tau_{\epsilon})}$ in Eq. (13). After some algebra we obtain instead of Eq. (40) the following expression for the force

[TABLE]

Equation (41) demonstrates that, neglecting phonon decay ($\tau_{0}\to\infty$), the momentum flux by the ballistic phonons decays as inverse distance due to phonon propagation Keldysh (1976) as it follows from the momentum conservation, Fig. 1(b). Finite lifetime of phonons gives rise to additional $\propto\exp{(-r/s\tau_{0})}$ decay with $s\tau_{0}$ being the mean distance which phonon propagates before it decays. The factor $\Theta\left(st-r\right)\exp{[({r-st})/(s\tau_{\epsilon})]}$ accounts for the retardation while phonons arrive from the hot spot. The regime of exciton drift due to the interaction with ballistic phonons is denoted hereafter as a phonon wind.

IV.2 Solution of drift-diffusion equation

We start with the simplified analysis of the exciton drift-diffusion equation (31) and identify the most important limits. We also confirm our analytical results by the numerical solution of the drift-diffusion equation.

It is noteworthy that in the diffusive regime of exciton propagation the exciton momentum is lost at a very short time scale $\tau_{p}$. That is why as soon as the force field $\bm{F}$ vanishes, i.e., due to the phonon diffusion or finite phonon lifetime, the exciton distribution evolves according to the diffusion equation, spreads in the space and decays due to the finite exciton lifetime $\tau_{d}$. By contrast, while substantial number of phonons are present, a competition between the drift and diffusion takes place. The diffusive flux of excitons $\propto D_{x}|\bm{\nabla}n|\sim D_{x}n/a$, where $a$ is the exciton spot size, Fig. 1(a). The drift flux of excitons is given by $\tau_{p}Fn/m$. Thus, the drift of excitons dominates provided

[TABLE]

The latter condition follows from the exciton diffusion coefficient (28). This condition can be naturally fulfilled in the phonon wind regime. The situation can be different in the case of the Seebeck or the phonon drag effect. Indeed, taking into account that the force is given by the temperature gradient [see Eqs. (35) and (37)] and making crude estimates with $\tau_{p}^{ph}\sim\tau_{x}$ we obtain that the drift and diffusion provide comparable contributions to the exciton dynamics if $T\lesssim T_{x}$.

In order to obtain the analytical result we assume that the drift dominates, i.e., that Eq. (42) is fulfilled. In this case the diffusive term in Eq. (31) can be neglected, we neglect also the Seebeck effect putting $\bm{\nabla}T_{x}$ to zero. The remaining equation can be solved by the method of characteristics. Namely, we find the exciton dynamics from the second Newton’s law which at $t\gg\tau_{p}$ reads

[TABLE]

Here $v_{dr}$ is the exciton drift velocity in the force field $\bm{F}$. The trajectories $\bm{r}(\bm{\rho}_{0},t)$ where $\bm{\rho}_{0}$ is the initial position of the exciton at $t=0$ provide implicit dependence of $\bm{\rho}_{0}$ on $\bm{r}$ and $t$ which makes it possible to express the exciton distribution via the initial one. Assuming axial symmetry of the problem we obtain, see Appendix A for details

[TABLE]

Let us first focus on the phonon wind regime where the force field $\bm{F}(\bm{r},t)$ is given by Eq. (41). It is noteworthy that exciton drift velocity, Eq. (43), cannot exceed the speed of sound $s$. Indeed, the exciton-phonon scattering tends to equalize the velocities of the particles while the exciton momentum scattering processes result in the velocity dissipation. We assume in what follows that $v_{dr}\ll s$ in agreement with the assumptions made at derivation of Eq. (31); the effects where $v_{dr}$ is comparable with the speed of sound are beyond the scope of the present work, see, e.g., Ref. Bulatov and Tikhodeev (1992) where saturation effects were studied. Thus, in Eq. (41) one can neglect $r$ as compared with $st$ and represent $\bm{F}_{\rm wind}$ in a very simple form

[TABLE]

where the parameter $U(r,t)=U_{0}\exp{(-t/\tau_{\epsilon}-r/s\tau_{0})}$ with $U_{0}\propto k_{B}T_{0}$ describes the hot spot efficiency. Equation (45) resembles the force field produced by a Coulomb center in a two-dimensional geometry. This is indeed consistent with the coordinate dependence of the momentum flux produced by the ballistic non-equilibrium phonons Keldysh (1976). Equation (43) together with Eq. (45) can be easily integrated with the result

[TABLE]

Particularly, for the excitons starting from the $\rho_{0}=0$ the dynamics is given by the following asymptotes depending on the relations between the parameters:

[TABLE]

Physically, at $t\to 0$ the excitons are strongly accelerated by the phonon wind since the force field (45) diverges as $1/r$ and giving rise to a formally infinite velocity in Eq. (43). It results in $r\propto\sqrt{t}$ dependence, as demonstrated by Eq. (47a). As time goes by, the action of the phonon wind diminishes and excitons stop to be driven, $r=\mbox{const}$, either because the phonons vanish [short $\tau_{\epsilon}$, Eq. (47b)], or because excitons eventually left the area covered by the phonon wind [long $\tau_{\epsilon}$, Eqs. (47c) and (47d)]. Thus, at long time scales the exciton dynamics is controlled by the diffusion and recombination processes.

This regime of exciton propagation is illustrated in Fig. 2. In our numerical calculations we took the initial distribution of excitons in the form

[TABLE]

Here $N_{x}=2\pi\int_{0}^{\infty}n_{0}(\rho)\rho d\rho$ is the total number of excitons in the excitation spot and $a$ is the effective spot radius, $a^{2}=2\pi\int_{0}^{\infty}n_{0}(\rho)\rho^{3}d\rho$ (half-width at half-maximum is $a\sqrt{\ln{2}}$). The exciton distribution was calculated by solving the drift-diffusion equation (31) neglecting the Seebeck effect, i.e., putting $\eta=0$; the detailed analysis of the Seebeck effect on the exciton propagation is given in Ref. 2019arXiv190602084P . We have also replaced the generation term $\dot{n}$ in the right-hand side with the initial condition (48) assuming instantaneous generation of excitons. For the phonon wind mechanism we take the force field in the form of Eq. (45).

Figure 2(a) shows the exciton distribution in the presence of wind (red solid line) and in the absence of wind (black dashed line) at a fixed time $t=1$ ns. The phonon wind markedly changes the distribution giving rise to the dip in the center which evolves in time as shown in the inset and also in more detail in Fig. 2(b) where the false color plot of exciton distribution is shown. We also illustrate the propagation of excitons which started at the origin (dashed lines) calculated after Eq. (46). It illustrates rapid expansion of the cloud and formation of a ring-shaped pattern [inset in Fig. 2(b)] at short time scales, Eq. (47a), and saturation at longer times, Eq. (47b). The halo-like shape of the exciton cloud is a result of the drift dominating over the exciton diffusion: the particles leave the hot spot faster than their distribution smoothens by the diffusion process. The sensitivity of exciton density profiles to the parameters of the phonon system, $\tau_{\epsilon}$ and $\tau_{0}$, is briefly analyzed in Appendix C. Taking shorter values of these times still results in the halo formation, but the size of the ring is respectively smaller, in agreement with approximate analytical expressions (47b) and (47c).

An interplay between the drift and diffusive behavior can be observed also in the phonon drag regime as illustrated in Fig. 3. This calculation has been carried out similarly to the presented above, but the force field was taken in the form of Eq. (38). Panels (a) and (b) demonstrate the exciton distribution evolution for two values of initial temperature gradient which correspond to relatively low and relatively high drift velocity of excitons, respectively. As expected, the significant drift of excitons and the halo formation are possible only provided that the drift velocity is large enough, so that the product $v_{dr}\tau_{d}\gtrsim\sqrt{D_{x}\tau_{d}}$. Similar behavior is also possible in the regime of Seebeck effect, where the exciton drift is produced by the gradient of their temperature. The detailed study of the Seebeck effect requires self-consistent solution of the exciton drift-diffusion equation and the equation for the exciton temperature. This effect is beyond the scope of the present work.

Above we demonstrated the scenario for non-diffusive propagation of excitons: energy relaxation of excitons results in the formation of highly non-equilibrium acoustic phonons which, in turn, drag excitons out of the excitation spot. The spatial profile of the exciton density acquires a non-monotonic halo-like shape with the dip in the middle and increased exciton density on the periphery. Such profiles of the exciton density have been recently observed in Ref. Kulig et al. (2018). They were tentatively attributed to the memory effects such as the heating of the exciton gas and subsequent variation of the Auger recombination rate. The non-radiative exciton-exciton annihilation additionally decreases the number of excitons in the middle of the excitation spot. Here, the memory comes from the overheated phonon subsystem or from the elevated temperature of the excitons (in the case of the Seebeck effect induced drift). The full quantitative description of the experimental data requires also inclusion of the Auger recombination of excitons (which, according to our estimations, does not substantially qualitatively affect the profiles shown in Figs. 2 and 3) as well as the analysis of the phonon transport conditions. Particularly, at the room temperature the phonon lifetimes are expected to be sufficiently short ruling out the phonon wind effect. By contrast, at low temperatures the phonon wind could dominate the exciton drift. In this work, we abstain from the detailed comparison of the predictions with experimental results of Ref. Kulig et al. (2018), since a combination of several factors may affect the halo-like profile formation in the experiments at the room temperature: In particular, the combination of the wind and drag effects (on short- and long-timescales), possibly, exciton-exciton or exciton-free carriers interaction, as well as the Seebeck effect (discussed as an origin of the halo-like profiles in Ref. 2019arXiv190602084P ) could result in the observed exciton density profiles. Further experiments, e.g., aimed at studies of exciton propagation as a function of temperature would be helpful to finally establish the origin of the halo effect.

V Conclusion

We have developed analytical theory of exciton drift and diffusion in the presence of non-equilibrium phonons in two-dimensional transition metal dichalcogenides. We demonstrate that the flow of phonons can drag excitons out of the excitation spot giving rise to the halo-like shape of the exciton spatial profile. Different regimes of the drift are identified: for ballistic phonons the excitons are affected by the phonon wind, while for the diffusive phonons the exciton drift can be caused by the effective force resulting from the phonon temperature gradient.

Acknowledgements.

The author is grateful to L.E. Golub, A. Chernikov, E. Malic, and S.G. Tikhodeev for valuable discussions. Partial support from Russian Science Foundation (project # 19-12-00051) is acknowledged.

Appendix A Analytical solution of the drift equation

Here we present analytical solution of Eq. (31) at $D_{x}=0$ and $\tau_{d}\to\infty$. This drift equation for the density $n$ in the presence of the central force $\bm{f}$ has the following form:

[TABLE]

Then the drift equation (49) assumes the form

[TABLE]

We introduce the solution of the characteristic equation [second Newton’s law, Eq. (43)] with the initial condition ${r(t=0)=\rho_{0}}$:

[TABLE]

The following relation between the derivatives of $\rho_{0}^{2}$ takes place:

[TABLE]

We now prove that the solution of Eq. (50) has the form of Eq. (44) with $\tau_{d}\to\infty$:

[TABLE]

Indeed the left hand side has the form

[TABLE]

whereas the right hand side is:

[TABLE]

Evaluating the derivatives we see that the left and right hand sides coincide proving Eq. (53).

Appendix B Interplay of the drift and diffusion

It is instructive to provide the detailed comparison of the analytical (at $D_{x}=0$) and full numerical solutions of drift-diffusion equation (31) in the regime of the phonon wind. We take in Eq. (45)

[TABLE]

In order to produce the analytical solution we first find the exciton trajectories in the force field:

[TABLE]

Equations (56) or (57) provide implicit dependence of $\rho_{0}$ on $r$ and $t$, i.e., the initial position of exciton to reach at the time $t$ the position $\rho$. Explicitly, it reads

[TABLE]

Figure 4 illustrates good agreement of the approximate analytical and exact numerical solution of the drift-diffusion equation.

Appendix C Halo sensitivity to the system parameters

In Fig. 5 we analyze the sensitivity of the exciton density profiles in the phonon wind regime to the key parameters of the phonon system, which are largely unknown at this stage: hot spot generation time $\tau_{\epsilon}$ and phonon lifetime $\tau_{0}$. We took the same parameters of the system as we used to calculate Fig. 2 but used shorter value of $\tau_{\epsilon}$ in Fig. 5(a) and shorter value of the product $s\tau_{0}$ in Fig. 5(b). In both cases the halo formation is clearly seen, but the halo radius is smaller than that in Fig. 2(b). This is consistent with approximate analytical expressions (47b) and (47c), respectively. Minor differences from the analytical expressions results from the approximations used to derive Eqs. (47a), while the evolution of the halo peaks is well described by Eq. (46) (see dashed lines in Fig. 5).