Theory of Exciton-Phonon Coupling

TL;DR

This paper develops a theoretical framework to understand how excitons interact with phonons in low-dimensional materials, affecting optical properties and exciton lifetimes, with a focus on deriving a practical expression for the exciton-phonon self-energy.

Contribution

The paper introduces a new, more accurate expression for the exciton-phonon self-energy that accounts for exciton binding and electron-hole correlations, improving upon previous models.

Findings

The derived expression captures temperature-dependent shifts in optical transition energies.

The model predicts exciton lifetime broadening due to phonon interactions.

Qualitative differences from previous models highlight the importance of electron-hole correlations.

Abstract

The effects of the electron-phonon interaction on optical excitations can be understood in terms of exciton-phonon coupling, and require a careful treatment in low-dimensional materials with strongly bound excitons or strong electron-hole interaction in general. Through phonon absorption and emission processes, the optically accessible excitons are scattered into otherwise optically dark finite-momentum exciton states. We derive a practical expression for the phonon-induced term of the exciton self-energy (denoted as the exciton-phonon self-energy) that gives the temperature dependence of the optical transition energies and their lifetime broadening resulting from the exciton's interaction with the phonons. We illustrate this theory on a two-dimensional model, and show that our expression for the exciton-phonon self-energy differs qualitatively from previous expressions found in the literature that neglect the exciton binding or electron-hole correlations.

Equations

$$\big{\{}c_{i},c_{{i^{\prime}}}^{\dagger}\big{\}}=c_{i}c_{{i^{\prime}}}^{\dagger}+c_{{i^{\prime}}}^{\dagger}c_{i}=\delta_{i{i^{\prime}}}$$

$$G_{i{i}^{\prime }}^0(\omega )=\frac{1}{\omega -{\epsilon }_i\pm i\eta }{\delta }_{i{i}^{\prime }}$$

$$L_{vc{v}^{\prime }{c}^{\prime }}^0(\omega )=\frac{1}{\omega -({\epsilon }_c-{\epsilon }_v)+i\eta }{\delta }_{c{c}^{\prime }}{\delta }_{v{v}^{\prime }}$$

$$H_{vc,v^{\prime }{c}^{\prime }}^{2p}=({\epsilon }_c-{\epsilon }_v){\delta }_{c{c}^{\prime }}{\delta }_{v{v}^{\prime }}+K_{vc,v^{\prime }{c}^{\prime }}$$

$$L_{vc{v}^{\prime }{c}^{\prime }}(\omega )=\sum _S\frac{A_{vc}^S{A}_{v^{\prime }{c}^{\prime }}^{S\ast }}{\omega -{\Omega }_S+i\eta }.$$

$$L_{S{S}^{\prime }}(\omega )=\sum _{vc,v^{\prime }{c}^{\prime }}{A}_{vc}^{S\ast }{A}_{v^{\prime }{c}^{\prime }}^{S^{\prime }}{L}_{vc{v}^{\prime }{c}^{\prime }}(\omega )=\frac{1}{\omega -{\Omega }_S+i\eta }{\delta }_{S{S}^{\prime }}.$$

$$\displaystyle\epsilon^{\prime\prime}(\omega)=\frac{4\pi^{2}e^{2}}{\omega^{2}}\sum_{S}\big{|}\sum_{vc}A_{vc}^{S}\bra{v}\mathbf{e}\cdot\mathbf{v}\ket{c}\big{|}^{2}\delta(\omega-\Omega_{S})$$

$$M_{\kappa }{\omega }_{\lambda }^2{U}_{\kappa j}^{\lambda }=\sum _{{\kappa }^{\prime },j^{\prime }}{\Phi }_{j{j}^{\prime }}^{\kappa {\kappa }^{\prime }}(\mathbf{q})U_{{\kappa }^{\prime }{j}^{\prime }}^{\lambda }$$

$$D_{\lambda }(\omega )=\frac{1}{\omega -{\omega }_{\lambda }+i\eta }-\frac{1}{\omega +{\omega }_{\lambda }+i\eta }.$$

$$V_{ep}=V_{ep}^{(1)}+V_{ep}^{(2)}=\sum _{i{i}^{\prime }\lambda }{g}_{i{i}^{\prime }\lambda }{A}_{\lambda }{c}_i^{†}{c}_{i^{\prime }}+\sum _{i{i}^{\prime }\lambda {\lambda }^{\prime }}{g}_{i{i}^{\prime }\lambda {\lambda }^{\prime }}^{(2)}{A}_{\lambda }{A}_{{\lambda }^{\prime }}{c}_i^{†}{c}_{i^{\prime }}$$

$$g_{i{i}^{\prime }\lambda }(\mathbf{k},\mathbf{q})=\sqrt{\frac{\hslash }{2M{\omega }_{\mathbf{q}\lambda }}}⟨{\varphi }_{i\mathbf{k}+\mathbf{q}}|{\nabla }_{\mathbf{q}\lambda }{V}^{SCF}(\mathbf{r})|{\varphi }_{\mathbf{k}{i}^{\prime }}⟩$$

$$g_{i{i}^{\prime }\lambda {\lambda }^{\prime }}^{(2)}(\mathbf{k},\mathbf{q},{\mathbf{q}}^{\prime })=$$

$$\times ⟨{\varphi }_{i\mathbf{k}+\mathbf{q}+{\mathbf{q}}^{\prime }}|\frac{1}{2}{\nabla }_{\mathbf{q}\lambda }{\nabla }_{{\mathbf{q}}^{\prime }{\lambda }^{\prime }}{V}^{SCF}(\mathbf{r})|{\varphi }_{i^{\prime }\mathbf{k}}⟩.$$

$${\nabla }_{\mathbf{q}\lambda }{V}^{SCF}(\mathbf{r})=\sum _{l\kappa j}{U}_{\kappa j}^{\lambda }(\mathbf{q})e^{i\mathbf{q}\cdot {\mathbf{R}}_l}{\nabla }_{l\kappa j}{V}^{SCF}(\mathbf{r}).$$

$$G_{i{i}^{\prime }}(\tau )$$

$$\times {⟨T_{\tau }{c}_i(\tau )c_{i^{\prime }}^{†}(0)V_{ep}({\tau }_1)\dots V_{ep}({\tau }_n)⟩}_0$$

$$c_i(\tau )=e^{\tau {\hat{H}}_0}{c}_i(0)e^{-\tau {\hat{H}}_0}$$

$$G_{i{i}^{\prime }}(\omega ,T)$$

$${\Sigma }_{i{i}^{\prime }}(\omega ,T)={\Sigma }_{i{i}^{\prime }}^{\text{FM}}(\omega ,T)+{\Sigma }_{i{i}^{\prime }}^{\text{DW}}(T).$$

$${\Sigma }_{i{i}^{\prime }}^{\text{FM}}(i{\omega }_n)=$$

$$\displaystyle\times\Big{[}\frac{N_{B}(\omega_{\lambda})+f(\varepsilon_{i^{\prime\prime}})}{i\omega_{n}-\varepsilon_{i^{\prime\prime}}+\omega_{\lambda}}+\frac{N_{B}(\omega_{\lambda})+1-f(\varepsilon_{i^{\prime\prime}})}{i\omega_{n}-\varepsilon_{i^{\prime\prime}}-\omega_{\lambda}}\Big{]}$$

$$\displaystyle\Sigma^{\text{DW}}_{ii^{\prime}}=\sum_{\lambda}g^{(2)}_{ii^{\prime}\lambda\lambda}\big{[}2N_{B}(\omega_{\lambda})+1\big{]}.$$

$${\epsilon }_i+\Delta {\epsilon }_i(T)={\epsilon }_i+\text{Re}{\Sigma }_{ii}({\epsilon }_i,T)$$

$${\tau }_i^{-1}(T)=2|\text{Im}{\Sigma }_{ii}({\epsilon }_i,T)|.$$

$${\Omega }_S+\Delta {\Omega }_S(T)$$

$${\tau }_S^{-1}(T)=2|\text{Im}{\Xi }_{SS}({\Omega }_S,T)|$$

$${\Lambda }_{vc,v^{\prime }{c}^{\prime }}(\tau )$$

$$\times$$

$${\Lambda }_{S{S}^{\prime }}(\omega ,T)=L_{S{S}^{\prime }}(\omega )+\sum _{S_1{S}_2}{L}_{S{S}_1}(\omega ){\Xi }_{S_1{S}_2}(\omega ,T){\Lambda }_{S_2{S}^{\prime }}(\omega ,T).$$

$${\Lambda }_{vc,v^{\prime }{c}^{\prime }}^0(\tau )$$

$$\times$$

$${\Lambda }_{vc,v^{\prime }{c}^{\prime }}^0(\omega ,T)$$

$$+$$

Full text

A wealth of two-dimensional and nanocrystalline materials with interesting optical properties have been studied in recent years, including transition-metal dichalcogenides, layered heterostructures and halide perovskites Bhimanapati et al. (2015); Mueller and Malic (2018); Schaibley et al. (2016); Novoselov et al. (2016); Kovalenko et al. (2017); Shi et al. (2018); Chen and Shi (2017). In these systems, the optical excitations lead to the formation of strongly bound electron-hole pair states known as excitons. Many of their useful opto-electronic properties (e.g., photocurrent generation, single-photon emission, etc.) depend on the scattering dynamics and diffusion of the excitons Shi et al. (2013); Yuan et al. (2017); Kozawa et al. (2014); Ruppert et al. (2017); Raja et al. (2018); Li et al. (2019a); Chi et al. (2018). This dynamics is governed by several processes: the interaction of excitons with defects, the exciton-exciton interaction, and the exciton-phonon interaction. In particular, exciton-phonon coupling effects can be identified by their distinctive temperature dependence, whether in the exciton mean free path, lifetimes or emission spectra Selig et al. (2016); Gong et al. (2018); Lv et al. (2021).

The exciton-phonon coupling mechanism originates from a combined action of the electron-hole and the electron-phonon interactions, both of which can be described from first principles. On the one hand, the electron-hole interaction underlies the formation of excitons, and can be addressed with the Bethe-Salpeter equation (BSE) within the ab initio $GW$-BSE method Strinati (1984, 2008); Albrecht et al. (1998); Benedict et al. (1998); Rohlfing and Louie (1998a, b, 2000); Onida et al. (2002). This method solves the interacting two-particle problem for an electron and a hole, and yields the exciton energies and wavefunctions, which allow to predict the optical absorption spectra of materials The electron-phonon interaction, on the other hand, has largely been studied within density functional perturbation theory (DFPT) Gonze (1997); Gonze and Lee (1997); Baroni et al. (2001), which provides an ab initio description of the phonon energy spectrum and coupling potential. This framework has been used to study the effect of phonons on the band structure and carrier mobility as a function of temperature Giustino (2017); Miglio et al. (2020); Poncé et al. (2020). Going beyond DFT, electron-electron correlation effects to the electron self-energy may further be included at the $GW$ level from ab initio Hybertsen and Louie (1986) using the $GW$PT method Li et al. (2019b) in computing the electron-phonon interaction.

Describing the dynamics of photoexcited states from first principles is a challenging task. Simulations of hot electrons were achieved by retaining the electron-phonon interaction only, with the rationale that electrons far from the band edges would scatter freely without forming bound excitons Bernardi et al. (2014); Molina-Sánchez et al. (2017); Caruso (2021). It is necessary, however, to include the electron-hole interaction to predict the lifetime of absorption and emission states when they originate from bound excitons.

An early attempt to compute the temperature-dependent broadening and renormalization of exciton states was based on a one-particle picture of the electron-phonon coupling Marini (2008). This scheme has been used to compute the absorption spectrum of h-BN and MoS2 at finite temperature Marini (2008); Qiu et al. (2013); Molina-Sánchez et al. (2016). In this approach, the electron-phonon renormalization and broadening of the band structure is computed before solving the BSE. This method does not, however, describe correctly the process where excitons scatter into finite-momentum bound states, which is necessary to enforce energy conservation. Alternatively, the supercell BSE technique Zacharias et al. (2015); Zacharias (2016); Zacharias and Giustino (2020) does account for the phonon-mediated interaction between optical excitons and a limited number of finite-momentum excitons commensurate to the supercell size. It does so only within a static approximations, which is valid for non-polar materials. This approach predicts the exciton energy renormalization as a function of temperature, but makes no prediction on the scattering lifetime of the excitons. Recent methods formally achieved a proper description of exciton dynamics with exciton-phonon scattering amplitude deduced from Fermi’s golden rule using exciton-phonon coupling matrix elements Brem et al. (2018, 2019); Christiansen et al. (2019); Chen et al. (2020). This approach enforces energy conservation, and is consistent with the theory presented in this paper, as well as other methods derived from many-body perturbation theory Cudazzo (2020). Another effect of phonons on the exciton binding energies comes from the lattice screening, and this has been recently computed from first principles Filip et al. (2021); Adamska and Umari (2021)

In this work, we develop a general theory of the exciton-phonon coupling that is amenable to first-principles calculation. The central object is the exciton-phonon self-energy (i.e., the contribution to the exciton self-energy due to exciton-phonon interaction; there are of course contributions due to other excitations in a system), which yields the energy renormalization of the exciton states, as well as their scattering lifetime due to phonons. We apply this theory to a tight-binding model in two dimensions and discuss how it differs from other methods. This paper is organized as follows. Section I reviews the theory of electron-hole and electron-phonon interactions. Section II presents the extension of the one-particle theory to the exciton-phonon coupling, and the main equations for the phonon-induced temperature-dependent exciton lifetimes and energies. In Sec. III, we apply this scheme to a 2D tight-binding model, and discuss the consequences of electron-hole interactions on the scattering dynamics of two-dimensional systems. The main findings are summarized in Sec. IV. Several mathematical details of the derivation can be found in supplemental materials SM .

I Electron-hole and electron-phonon interactions

As a starting point for the treatment of electron-hole and electron-phonon interactions, we consider a mean-field fixed-ions Hamiltonian for the electrons $H_{0}=\hat{T}_{k}+V^{SCF}(\mathbf{r})$, where $\hat{T}_{k}$ is the kinetic energy operator, and $V^{SCF}(\mathbf{r})$ is the self-consistent field potential. Within the density functional theory (DFT) framework, $V^{SCF}$ includes the potential of the ions, the Hartree potential and the exchange-correlation potential. Solving the one-particle Hamiltonian yields the set of unperturbed wavefunctions $\phi_{i}(\mathbf{r})$ and energies $\varepsilon_{i}$, where the label $i$ comprises a band index ($n_{i}$), a wavevector ($\mathbf{k}_{i}$), and eventually a spin index ($\sigma_{i}$).

These quantities are used to construct the time-ordered Green’s function, defined in the one-particle basis and in time as $G^{0}_{i{i^{\prime}}}(t)\!=\!-i\braket{\hat{T}_{t}c_{i}(t)c_{i^{\prime}}^{\dagger}(0)}_{0}$, where $c_{i}^{\dagger}$ and $c_{i}$ are the electron creation and annihilation operators, $\hat{T}_{t}$ is Wick’s time-ordering operator, and ”[math]” indicates here that the expectation value is taken over a ground state that is not perturbed by the phonons. The creation and annihilation operators follow the commutation relations

[TABLE]

and $\{c_{i},c_{{i^{\prime}}}\}=\{c_{i}^{\dagger},c_{{i^{\prime}}}^{\dagger}\}=0$. In terms of frequencies, the one-particle Green’s function writes as

[TABLE]

where $\pm\eta$ is an real infinitesimal number with the same sign as $\varepsilon_{i}$, the eigenvalue measured with respect to the chemical potential.

The electronic energies and wavefunctions that define the starting point need not be obtained from DFT; they may also be obtained from a model hamiltonian or from a many-body scheme. It is required however to have an effective Hamiltonian that depends implicitly on the atomic coordinates and that can be differentiated with respect to these coordinates to obtain the dynamical matrix and the electron-phonon coupling elements. Such a Hamiltonian, for example, may be constructed from the self-energy computed in the $GW$ formalism van Schilfgaarde et al. (2006); Antonius et al. (2014); Li et al. (2019b).

I.1 Electron-hole interaction

A class of neutral excitations of an insulating system is composed of an electron being promoted into the conduction bands and leaving a hole in the valence bands. If the Coulomb attraction between the electron and the hole is sufficiently strong, they may form a bound exciton, that is, a bound state whose excitation energy is smaller than the fundamental band gap. The procedure to compute the exciton spectrum from the BSE is described in Ref. Rohlfing and Louie, 2000.

The starting point to describe excitons is the set of independent (or non-interacting) electron-hole pairs, typically with quasiparticle energies from a $GW$ calculation in the ab initio $GW$-BSE approach. The corresponding propagator for these fictitious excitations is

[TABLE]

where $v$ and $c$ refer to the labeling of occupied (valence band) and unoccupied (conduction band) states which involve both band and $\mathbf{k}$-point indices. We call $L^{0}$ the bare electron-hole propagator. It is computed with the fixed-ions Hamiltonian, but the superscript ”[math]” refers to the fact that the electron and hole propagate freely without interacting with one another.

The BSE relates the bare exciton propagator $L$ to the bare electron-hole propagator $L^{0}$ through the kernel $K$, as depicted in Fig. 1. The BSE kernel $K$ is composed of an attractive screened Coulomb interaction between the electron and the hole, and a repulsive Coulomb exchange term. In practice, the BSE is solved by diagonalizing an effective two-particle Hamiltonian,

[TABLE]

where $K_{vc,v^{\prime}c^{\prime}}$ is the static version of the BSE kernel, and the first term is the sum of the quasiparticle energies of the electron and hole. This Hamiltonian yields the exciton energies $\Omega_{S}$ and electron-hole coefficients $A_{vc}^{S}$ for each exciton $S$. The resonant part of the bare exciton propagator in the quasiparticle basis can thus be written as

[TABLE]

Like the one-particle state indices, the label $S$ comprises a set of discrete quantum numbers as well as the center-of-mass momentum of the exciton ($\mathbf{q}_{S}$) which is the wavevector difference ($\mathbf{k}_{c}-\mathbf{k}_{v}$) between the unoccupied and occupied single-particle orbitals that form the exciton. (The periodicity of the crystal dictates that all free electron-hole pairs forming the exciton have the same wavevector difference modulo a reciprocal lattice vector.) Because of the small momentum carried by photons, only excitons with $\mathbf{q}_{S}\approx 0$ are optically accessible. We will see that finite-momentum excitons ($\mathbf{q}_{S}\neq 0$) are important to describe scattering events by phonons.

Since the $A_{vc}^{S}$ are eigenfunctions of the two-particle Hamiltonian, they form a complete orthonormal basis for the space spanned by the set of valence and conduction bands used to construct $H^{2p}$. We can use these coefficients to transfer between the $vc$ basis and the $S$ basis. For example, the bare exciton propagator writes

[TABLE]

The absorption spectrum at zero temperature without electron-phonon interaction can be constructed from the exciton energies and electron-hole coefficients. It takes the form

[TABLE]

where $\mathbf{v}$ is the velocity operator, $\mathbf{e}$ is the photon’s polarization vector, and the summation over exciton states is restricted to zero-momentum excitons. The delta function in Eq. (7) is typically represented as a Lorentzian function with a certain broadening. In most past calculations, this broadening was chosen empirically to reproduce the available experimental data. The absorption line broadening is in fact related to the lifetime of the excitons and, as we show, can be computed from first principles.

I.2 Lattice dynamics and electron-phonon interaction

The lattice dynamics can be obtained from a self-consistent calculation of the dynamical matrix of the crystal ($\Phi$), as detailed in Ref. Gonze and Lee, 1997. In real space, the dynamical matrix (or force matrix) corresponds to the second-order derivative of the total energy with respect to the displacement of two atoms: $\Phi^{\kappa\kappa^{\prime}}_{jj^{\prime}}(\mathbf{R}_{l}-\mathbf{R}_{l^{\prime}})=\nabla_{l\kappa j}\nabla_{l^{\prime}\kappa^{\prime}j^{\prime}}E$, where $l$ labels a unit cell of the crystal with lattice vector $\mathbf{R}_{l}$, $\kappa$ labels an atom in the unit cell, and $j$ labels a Cartesian direction. The phonon frequencies $\omega_{\lambda}$ and polarization vectors $U^{\lambda}_{\kappa j}$, are obtained by diagonalizing the Fourier transform of $\Phi$ as

[TABLE]

where $M_{\kappa}$ is the mass of an atom. The label $\lambda$ for a phonon mode comprises a branch index and a wavevector ($\mathbf{q}$ or $\mathbf{q}_{\lambda}$). The DFPT method allows for the self-consistent calculation of the first-order derivative of the local potential with respect to atomic displacements, $\nabla_{l\kappa j}V^{SCF}(\mathbf{r})$. Thanks to the Hellman-Feynman theorem and the $2n+1$ theorem, only the first-order derivatives of the potential and density need to be computed self-consistently in order to construct the dynamical matrix Gonze and Vigneron (1989).

The phonon propagator is defined in time as $D_{\mathbf{q}\lambda}(t)\!=\!-i\braket{\hat{T}_{t}A_{\mathbf{q}\lambda}(t)A_{-\mathbf{q}\lambda}(0)}_{0}$, where $A_{\mathbf{q}\lambda}\!=\!a^{\dagger}_{\mathbf{q}\lambda}\!+\!a_{\mathbf{q}\lambda}$ is the sum of a phonon creation and annihilation operator. In frequency space, the phonon propagator writes as

[TABLE]

The electron-phonon interaction stems from the perturbation to the fixed-atoms Hamiltonian created by the phonons. A thorough discussion of the electron-phonon interaction in the context of DFT is presented in Ref. Marini et al., 2015; Giustino, 2017. Expanding the Hamiltonian up to second order in the atomic displacements, the electron-phonon coupling potential writes as

[TABLE]

with the first-order electron-phonon coupling matrix elements

[TABLE]

and the second-order matrix elements

[TABLE]

The derivative of the potential with respect to a phonon mode is defined as

[TABLE]

The wavevectors ($\mathbf{k},\mathbf{q}$) are written explicitly in Eq. (11) and (I.2), but will be omitted in the remainder of the paper to lighten the notation. Also, to make the units explicit, we use $\hbar=1$ and the mass $M$, which normalizes the phonon eigenvectors according to $\sum_{\kappa}M_{\kappa}|U^{\lambda}_{\kappa}|^{2}=M$. It is useful to assign $M$ the value of the average mass of the atoms of the unit cell, so that the phonon eigenvectors are on the order of unity, and the factor $1/M$ serves as an expansion parameter. Note that each summation over phonon modes in Eq. (10) is implicitly normalized by $\sqrt{N_{\mathbf{q}}}$, where $N_{\mathbf{q}}$ is the number of wavevectors used to sample the Brillouin zone.

The quantity $g_{ii^{\prime}\lambda}$ is the electron-phonon coupling matrix element between the one-particle states $i$ and $i^{\prime}$ via the phonon mode $\lambda$. Because of crystal momentum selection rule, for $g_{ii^{\prime}\lambda}$ to be non-zero, we must have $\mathbf{k}_{i}=\mathbf{k}_{i^{\prime}}+\mathbf{q}_{\lambda}+\mathbf{G}$ where $\mathbf{G}$ is a reciprocal lattice vector that is non-zero for an Umklapp process. For $g^{(2)}_{ii^{\prime}\lambda\lambda^{\prime}}$, at the lowest order of perturbation theory, we must have $\mathbf{k}_{i}=\mathbf{k}_{i^{\prime}}$ and $\mathbf{q}_{\lambda^{\prime}}=-\mathbf{q}_{\lambda}$ for the wavevectors, and $\lambda^{\prime}=\lambda$ for the branch indices. In practice, the second-order electron-phonon coupling matrix elements $g^{(2)}_{ii^{\prime}\lambda\lambda}$ are never computed explicitly. Their contribution to the self-energy is approximated in terms of the first-order electron-phonon coupling matrix elements by making use of the accoustic sum rule and the rigid-ion approximation Allen and Heine (1976); Poncé et al. (2014).

I.3 One-particle propagators

We employ the Matsubara formalism to treat the propagators and self-energies at finite temperature. The propagators are defined on the imaginary time axis $\tau=it$ in the interval $[-\beta,\beta]$ with $\beta=1/k_{B}T$, and are made periodic outside of this range. The one-particle propagator $G$ is composed of odd (fermionic) Matsubara frequencies, while the even frequencies compose all bosonic propagators: electron-hole ($L^{0}$), exciton ($L$), and phonon ($D$).

The interacting Green’s function $G$ can be expanded in powers of the perturbation as

[TABLE]

and we adopt the convention that in each expectation value ${\langle\dots\rangle}$, only distinct and connected diagrams should be retained Mahan (2000). The time-dependent operators are expressed in the interaction picture, that is

[TABLE]

with the mean field Hamiltonian ${\hat{H}_{0}=\sum_{i}\varepsilon_{i}c^{\dagger}_{i}c_{i}}$. Equation (I.3) can be cast into a Dyson equation for $G$, depicted in Fig. 2, as

[TABLE]

which defines the one-particle electron-phonon self-energy $\Sigma$.

The perturbative expansion of $G$ with Eq. (I.3) yields powers of $V_{ep}^{(1)}$ and $V_{ep}^{(2)}$, where $V_{ep}^{(1)}$ appears an even number of times in each term. The Migdal theorem states that this expansion can be truncated to the lowest power of $1/M$. Since $V_{ep}^{(1)}$ is proportional to $1/\sqrt{M}$ and $V_{ep}^{(2)}$ is proportional to $1/M$, two terms remain in the self-energy after truncation: the Fan-Migdal (FM) term and the Debye-Waller (DW) term, written as

[TABLE]

The FM term is dynamic (i.e. complex and frequency-dependent) and stems from the first-order electron-phonon coupling potential. Its analytic expression is

[TABLE]

where $N_{B}(\omega)$ is the Bose-Einstein distribution, and $f(\omega)$ is the Fermi-Dirac distribution, both of which depend implicitly on temperature. The self-energy is evaluated on the real frequency axis with the analytic continuation $i\omega_{n}\rightarrow\omega+i\eta$ where $\eta$ is an infinitesimal real number that is taken positive for the retarded Green’s function (for unoccupied states or electrons) and negative for the advanced Green’s function (for occupied states or holes).

The Debye-Waller term is static (real and frequency-independent), and it stems from the second-order electron-phonon coupling potential. It writes

[TABLE]

The self-energy allows for the mixing of electronic states through the off-diagonal components of the self-energy ($i\neq i^{\prime}$). This mixing is only possible among states with the same crystal momentum ($\mathbf{k}_{i}=\mathbf{k}_{i^{\prime}}$). If the bands are well-separated in energy, one may use only the diagonal elements of the self-energy to obtain the renormalization of the bands as

[TABLE]

and in general, $\Delta\varepsilon_{i}(0)\neq 0$. Correspondingly, the inverse lifetime of the electronic state is

[TABLE]

In Eq. (20) and (21), the temperature renormalization and lifetime is computed in the on-the-mass-shell limit, that is, by evaluating the self-energy at the bare energy, rather than the renormalized energy. This is the preferred approach for phonon perturbations, as it gives results in good agreement with a self-consistent calculation of the self-energy Brown-Altvater et al. (2020).

II Exciton-Phonon self-energy: Self-energy of excitons from coupling to phonons

We seek similar expressions for the exciton energy ($\Omega_{S}$) and lifetime ($\tau_{S}$) due to coupling to phonons that would allow one to compute absorption spectra at zero and finite temperature through Eq. (7). The phonon-induced corrections will be given by the diagonal components of the exciton-phonon self-energy ($\Xi_{SS}$) to be discussed below, namely,

[TABLE]

for the exciton energies, and

[TABLE]

for the inverse lifetime of the excitons, that is, the absorption line broadening.

The expression for the exciton-phonon self-energy is obtained by considering the interacting exciton propagator ($\Lambda$) defined as

[TABLE]

where $\tau^{+}$ is a time infinitesimally larger than $\tau$, and the interaction term $\hat{H}_{1}(\tau)$ is the sum of the static electron-hole interaction kernel $\hat{K}$ and the electron-phonon coupling potential $\hat{V}_{ep}$. Equation (II) reduces to the BSE if $\hat{V}_{ep}=0$.

The exciton-phonon self-energy connects the interacting exciton propagator $\Lambda$ to the bare exciton propagator $L$ (without electron-phonon interaction) in a Dyson-like equation:

[TABLE]

In order to express the exciton-phonon self-energy in an analytic form, it is useful to first consider the case without the electron-hole interaction.

II.1 Independent electron-hole-pair-phonon self-energy

We define $\Lambda^{0}$ as the electron-hole propagator with the electron-phonon interaction, but without the electron-hole interaction, calling it the independent electron-hole-pair-phonon (IEHPP) propagator, given by

[TABLE]

The corresponding IEHPP self-energy $\Xi^{0}$ is defined through the Dyson equation

[TABLE]

which is depicted in Fig. 3.

A detailed derivation of $\Xi^{0}$ is provided in supplemental materials SM . Throughout this derivation, we assume the existence of a band gap that does not allow for a significant density of thermal carriers at the temperature of interest, which is typically room temperature, where $E_{g}\gg k_{B}T$.

The different contributions to the IEHPP self-energy are grouped in three terms, according to the topology of the diagrams depicted in Fig. 4

[TABLE]

The Fan-Migdal term is further split in two contributions: the dynamic term (FMd) and the static term (FMs),

[TABLE]

The dynamic FM term describes the propagation of an electron-hole pair being temporarily knocked into a different electron-hole state while absorbing/emitting a phonon. Either the hole scatters into another valence state of different momentum, or the electron scatters into another conduction state. Its diagrammatic expression is given by

[TABLE]

Performing the convolution of the bare electron-hole propagator $L^{0}$ with the phonon propagator $D$, we obtain the analytic expression

[TABLE]

and one can safely assume that $N_{B}(\varepsilon_{c}-\varepsilon_{v},T)\ll N_{B}(\omega_{\lambda},T)$. The static FM term has the same topology as the dynamic FM term. It is given by

[TABLE]

and its analytic expression is

[TABLE]

where the symbol $\mathfrak{S}$ means that the terms in brackets should be symmetrized with the substitution $v,c\leftrightarrow v^{\prime},c^{\prime}$ and a factor of $1/2$. The static FM term includes transitions that cannot be described as an intermediate electron-hole pair, such as the hole being coupled to a conduction band state or the electron being coupled to a valence band state. These transitions are only virtual, in the sense that they do not conserve energy and do not contribute to the imaginary part of the self-energy.

The next set of diagrams are the phonon exchange (X) term, defined as

[TABLE]

The analytic expression for this term is

[TABLE]

In this process, the electron emits a phonon that is being absorbed later on by the hole, or vice-versa. This term is exclusively off-diagonal, since the electron and the hole exchange momentum and are being scattered into different states. It will be non-zero only when ${\mathbf{k}_{v^{\prime}}=\mathbf{k}_{v}+\mathbf{q}_{\lambda}+\mathbf{G}}$ and ${\mathbf{k}_{c^{\prime}}=\mathbf{k}_{c}+\mathbf{q}_{\lambda}+\mathbf{G}}$.

Finally, the Debye-Waller contribution to $\Xi^{0}$ is simply the second-order interaction of the electron and the hole with the phonon modes, giving

[TABLE]

Consider a non-interacting electron-hole excited state with energy $\varepsilon_{c}-\varepsilon_{v}$. We can show that the diagonal component of the IEHPP self-energy for this state is

[TABLE]

This is the expected result. Without the electron-hole interaction, the corrections to the optical excitations are simply given by the electron-phonon interaction corrections to the one-particle energies that make up the transitions. Since the imaginary part of the self-energy has opposite signs for electrons and holes, we also have that

[TABLE]

The broadening of a non-interacting electron-hole transition is thus the sum of the broadenings of the one-particle states.

II.2 Exciton-phonon self-energy

The exciton-phonon self-energy $\Xi$ is obtained by expanding the interacting exciton propagator $\Lambda$ in the bare exciton basis as

[TABLE]

where the last term regroups all the ”non-excitonic” diagrams and will be discussed further below. In this formulation, the electron-hole interaction is already included in $L$ (the bare exciton propagator without considering coupling of the exciton to phonons), and one needs to expand the perturbation in powers of $V_{ep}$ only. The solution is thus analogous to the one-particle case, i.e., considering an exciton as a single particle interacting with the phonons. The physical motivation for this formulation is that the two kind of processes occur on different time scales. The electron-hole interaction is much faster than the electron-phonon interaction, since the plasmon frequency is much larger than the typical phonon frequency.

We define the exciton annihilation operator as

[TABLE]

with the coefficient $A^{S}$ from solving the BSE, and we treat these excitations as bosons. This treatment is not exact, since the Bose commutation relations are not exactly fulfilled by the operators in Eq. (40). However, in the low exciton density limit Laikhtman (2007) we can take

[TABLE]

The bare exciton propagator is diagonal in the exciton basis and is formally defined as

[TABLE]

where the subscript $2p0$ refers to the ground state of a Hamiltonian that includes electron-hole interactions through the BSE kernel but does not include electron-phonon interactions. In this picture, the time dependence of the operators is given by

[TABLE]

with the two-particle Hamiltonian ${\hat{H}^{2p}=\sum_{S}\Omega_{S}c^{\dagger}_{S}c_{S}}$. The perturbation terms $V_{ep}=V^{(1)}_{ep}+V^{(2)}_{ep}$ appearing in Eq. (II.2) may be expressed in the exciton basis as

[TABLE]

and

[TABLE]

Here, as defined above, $A_{\lambda}$ is a sum of a phonon creation and annihilation operator, (not to be confused with the electron-hole coefficients $A^{S}_{vc}$). The first- and second-order exciton-phonon coupling matrix elements are the transition amplitude between two exciton states mediated by the electron-phonon coupling potential, e.g.

[TABLE]

Using Eq. (40) and (1) to simplify the expectation values and discarding terms like $\delta_{SS^{\prime}}\sum g_{vv\lambda}$ which produce disconnected diagrams, we arrive at the expressions

[TABLE]

and

[TABLE]

Momentum conservation dictates that $\mathbf{q}_{S}=\mathbf{q}_{S^{\prime}}+\mathbf{q}_{\lambda}$ for $g_{SS^{\prime}\lambda}$ to be nonzero, while $\mathbf{q}_{S^{\prime}}=\mathbf{q}_{S}$ for $g^{(2)}_{SS^{\prime}\lambda\lambda}$. It is worth writing explicitly all the wavevectors in these expressions. For an exciton $S^{\prime}$ with center-of-mass momentum $\mathbf{q}_{S^{\prime}}=\mathbf{Q}$ and an exciton $S$ with center-of-mass momentum $\mathbf{q}_{S}=\mathbf{Q}+\mathbf{q}$, the first-order matrix element writes as

[TABLE]

where $A^{S\mathbf{Q}}_{vc\mathbf{k}}$ is the coefficient of the exciton envelope function in k-space for a hole at $\mathbf{k}$ and an electron at $\mathbf{k}+\mathbf{Q}$. We note that $(\mathbf{Q},\mathbf{q})$ are not Fourier components of the matrix element; they are momentum information in fact that is contained in the quantum numbers $S$ and $S^{\prime}$. Here we use a somewhat redundant notations to draw out the physics of center-of-mass momentum conservation of the excitons. The second-order matrix elements between two exciton states with the same momentum $\mathbf{Q}$ are

[TABLE]

Note that the second-order matrix elements do not transfer momentum between the exciton states, even when $\mathbf{q}\neq 0$.

We may now compute the self-energy from Eq. (II.2). Some of the diagrams using quasiparticle (electron/hole) basis for the exciton-phonon self-energy up to second order in electron-phonon coupling are illustrated in Fig. 5. The Debye-Waller diagram coming from the second-order matrix element is

[TABLE]

This is the same result as for the case without electron-hole interactions, meaning that, using the basis transformation rule of Eq. (6), we find $\Xi^{\text{DW}}_{SS^{\prime}}=\Xi^{0^{\text{DW}}}_{SS^{\prime}}$.

The first-order coupling matrix elements yield the Fan-Migdal term. In the perturbative expansion of $\Lambda$, the first term of Eq. (II.2) produces only the dynamic part of the FM term, since only this part can be expressed in terms of $g_{SS^{\prime}\lambda}$. It is given by

[TABLE]

with the analytic expression

[TABLE]

We can verify that in the limit where the electron-hole interaction vanishes, that is, when each exciton is composed of a single electronic transition, Eq. (II.2) reduces to Eq. (II.1).

In the exciton basis, the phonon exchange term doesn’t need to be added; it is included in the dynamic FM term. This term arises from the cross terms between valence-to-valence and conduction-to-conduction bands coupling when squaring Eq. (47). Since the exciton states mix several electronic transitions across the Brillouin zone, the phonon exchange term does contribute to the diagonal elements of $\Xi$.

All other diagrams of $\Lambda$ in which the scattered state cannot be projected onto the exciton basis are contained in $\Lambda^{\text{non-excitonic}}$. One such diagram is obtained from the static FM term. Therefore, we complete the interacting FM self-energy by adding the static FM term of Eq.(II.1) so that

[TABLE]

with $\Xi^{\text{FMs}}_{SS^{\prime}}=\Xi^{0^{\text{FMs}}}_{SS^{\prime}}$. Using the non interacting static FM term here is an approximation. Since the static FM term is composed of virtual transitions across the band gap, the relative error induced by this approximation is at most on the order of $E^{S}_{b}/E_{g}$, where $E^{S}_{b}$ is the binding energy of the exciton.

As formulated, in practical implementation, there will be an additional term in the self-energy, which contributes to $\Lambda^{\text{non-excitonic}}$. Most of the time, the basis of valence and conduction bands used to expand the exciton wavefunctions is not complete; only a subset of the bands are used. Moreover, one typically doesn’t compute all the possible excitons states, but only the first few solutions of the BSE with the lowest eigenvalues. Therefore, there will be a missing contribution to $\Xi^{\text{FMd}}$, which we call the completion term, written as $\Xi^{\text{C}}$. Since the missing contribution involves high-energy states, this term can be computed by taking these states as free electron-hole pairs. One can construct $\Xi^{\text{C}}$ from the expression of $\Xi^{0}$, where the intermediate states $v^{\prime\prime}$ and $c^{\prime\prime}$ are projected outside of the basis of excitons that were computed. An example of this procedure is given in Sec. III.

Collecting all terms, the final expression for the exciton-phonon self-energy reads

[TABLE]

Thus, $\Xi$ contains all the diagrams of $\Xi^{0}$ but with addition of the electron-hole interaction diagrams, as shown in Fig. 5. Note that $\Xi^{\text{C}}_{SS^{\prime}}$ is complex and frequency-dependent, but its imaginary part is non-zero only at frequencies far from the bare exciton energy, and its real part varies smoothly near the exciton energy. Hence, near the exciton energy ($\Omega_{S}$) we may consider that only the dynamic Fan-Migdal term contributes to the imaginary part of the self-energy, which is

[TABLE]

The temperature renormalization of the exciton energies and their lifetime can finally be computed according to Eqs. (22) and (23).

II.3 Approximate expressions

Since the excitation energies are usually much larger than the phonon frequencies, it is a safe approximation to use the fact that $N_{B}(\Omega_{S},T)\ll N_{B}(\omega_{\lambda},T)$ and write for the diagonal part of the self energy

[TABLE]

Let us compare the full exciton-phonon self-energy with approximate expressions found in literature. We first define the ”uncorrelated exciton” (UE) approximation to the exciton-phonon self-energy:

[TABLE]

This approach is equivalent to the one used in previous studies Marini (2008); Qiu et al. (2013); Molina-Sánchez et al. (2016). Several approximations have been made between Eqs. (55) and (58). First, the UE self-energy is constructed using the independent electron-hole-pair-phonon propagator (IEHPP). Then, only the diagonal elements of $\Xi^{0}$ are used when transforming from the one-particle basis to the exciton basis. Finally, the self-energy is being evaluated at the non interacting transition energies. The uncorrelated exciton approximation reduces tremendously the computational effort needed to obtain the self-energy, since the knowledge of the finite-momentum exciton states is no longer required. However, it does not yield accurate results for the lifetime of bound excitons, as we show in the next section.

In this approximation, the uncorrelated exciton is just a wavepacket composed of electron-hole pairs that do not interact with each other and hence there is no binding energy. This becomes apparent when looking at the imaginary part of the self-energy. Using Eq. (II.1), the imaginary part of $\Xi^{\text{UE}}_{S}$ writes as

[TABLE]

Compare this expression with Eq. (II.2), where the scattering occurs between exciton states. This very drastic approximation in past studies creates a qualitative difference in the physics of the exciton in the low temperature limit, where ${N_{B}(\omega_{\lambda},T)\rightarrow 0}$. For the lowest bound exciton, from Eq. (II.2), no states are available for scattering through phonon emission. Hence, the electron-phonon interaction does not contribute to its inverse lifetime. The physical lifetime of this exciton comes from other processes, such as radiative recombination or scattering by impurities Spataru et al. (2005). In the UE approximation however, scattering events still occur because the binding energy of the initial and final states are ignored, and the lowest bound exciton has an unphysical lifetime at zero temperature.

For the temperature-dependent energy shift of the excitonic peaks, the accuracy of the UE approximation is unknown at this point. In the following section, we gain some intuition on this matter by computing the self-energy in a model system.

III Application to a model system

To illustrate the theory, we present a two-band model in two dimensions, and compute the energy and lifetime of optical excitations as a function of temperature. We use the triangular lattice exciton model proposed by D. Gunlycke and F. Tseng Gunlycke and Tseng (2016), which mimics the spin-dependent band dispersion of transition metal dichalcogenides near the $K$ and $K^{\prime}$ valleys, and yields realistic exciton binding energies. The main equations are reported below, and we refer the reader to Ref. Gunlycke and Tseng, 2016 for further details of this model.

III.1 One-particle and two-particle Hamiltonians

The one-particle tight-binding Hamiltonian is $H=\sum_{n\sigma}H_{n\sigma}$ with

[TABLE]

where $c^{\dagger}_{n\sigma\mathbf{R}}$ and $c_{n\sigma\mathbf{R}}$ are creation and destruction operators for an electron in the  $n^{\text{th}}$ orbital on the lattice site $\mathbf{R}$ with spin $\sigma$.

The first term describes the hopping between a site and one of its six closest neighbours, where $\bm{\delta}$ is the lattice vector from one site to the other. The hopping amplitude is composed of a spatial amplitude $t_{n}$ and a spin-orbit coupling parameter $\tilde{t}_{n}$ according to

[TABLE]

where $\sigma=\pm 1/2$ and $\mathbf{K}$ is one corner of the 2D hexagonal Brillouin zone. The second term of Eq. (60) is diagonal in lattice site, with $\varepsilon_{n}$ being the on-site parameters. The solutions for the one-particle energies at wavevector $\mathbf{k}$ are

[TABLE]

In order for this Hamiltonian to reproduce the main features of typical transition metal dichalcogenides for the valence (v) and conduction (c) band, the parameters are chosen by imposing that the band structure has a band gap $E_{g}$ located at the $K$ and $K^{\prime}$ valleys, with band effective masses $m^{*}$ and a splitting $\Delta$ between different spin bands. The on-site and hopping parameters are thus

[TABLE]

with $t={2\hbar^{2}}/{3m^{*}a^{2}}$ and $a$ is the lattice constant. The resulting band structure is shown in Fig. 6(a) with the parameters $m^{*}=0.49$, $a=3.13$ Bohr, $E_{g}=2.5$ eV, and $\Delta=425$ meV. Each spin channel admits one valence band and one conduction band. The two spin channels are degenerate for the conduction band, while for the valence band, the splitting of the up and down spins due to spin-orbit coupling reaches a maximum at the $K$ and $K^{\prime}$ points in the Brillouin zone.

The BSE Hamiltonian, which yields the exciton bands, is made of two terms: the kinetic energy, and the Coulomb interaction between the electron and the hole. Translational symmetry implies that the exciton states possess a well-defined center-of-mass momentum $\mathbf{Q}$, and the BSE Hamiltonian writes

[TABLE]

where $b^{\dagger}_{\sigma\mathbf{R}}$ creates an electron at $\mathbf{R}_{e}=\mathbf{R}$ with spin $\sigma$, and a hole at $\mathbf{R}_{h}=0$ with spin $-\sigma$, meaning that no spin-flip occurs in this process. The BSE hopping parameter carries the information on the exciton momentum, and writes

[TABLE]

Together, the on-site terms that involve $\varepsilon_{c}$ and $\varepsilon_{v}$ and the hopping term form the kinetic energy component of the excitons. In this model, the Coulomb interaction term is simply taken to be

[TABLE]

where $\epsilon$ is the dielectric constant of the medium and $\Delta v_{0}$ is the difference between the screened Coulomb interaction and the bare exchange integral at $\mathbf{R}=0$, which we set to $\Delta v_{0}=1.6$ eV. Note that, for this model, the spin degrees of freedom do not mix, and at $\Gamma$, both spin channels are equivalent and can be interpreted as singlet exciton states. In what follows, we will assume that the phonons are non-magnetic and cannot flip the spin. The index $\sigma$ will thus be omitted.

Due to the sparse nature of a Hamiltonian with nearest-neighbour hopping, the BSE can be solved efficiently in real space for all $\mathbf{Q}$-vectors. Solving the BSE yields the exciton energies $\Omega_{\mathbf{Q}S}$ and wavefunctions in real space $A^{S}(\mathbf{R},\mathbf{Q})$. Their Fourier transform give the electron-hole coefficients $A^{S}(\mathbf{k},\mathbf{Q})$. Note that, within this model, it is not necessary to specify the one-particle band indices for these coefficients, since there is only a single valence band and conduction band.

The exciton band structure is shown in Fig 6(b). The binding energy of the first four optical excitons is between 100 and 500 meV, and their wavefunction in real space is depicted in Fig. 6(c-f).

III.2 Exciton-phonon coupling

We model the lattice vibration spectrum as a single dispersionless phonon band with frequency $\omega_{0}$. The electron-phonon coupling strength is also taken to be independent of the phonon wavevector, and we use a single parameter $g$ to represent intra-band coupling ($g_{cc^{\prime}}\!=\!g_{vv^{\prime}}\!=\!g$ ; $g_{cv}\!=\!0$). For simplicity of notation, we write the intra-band coupling constants as $g_{c}$ and $g_{v}$. In the present calculation, we set $\omega_{0}\!=\!50$ meV and $g\!=\!250$ meV.

The IEHPP self-energy for a hole at $\mathbf{k}$ and an electron at $\mathbf{k}+\Gamma$ is diagonal in $\mathbf{k}$ indices and is given by

[TABLE]

where $P_{+}(T)=N_{B}(\omega_{0},T)$ and $P_{-}(T)=N_{B}(\omega_{0},T)+1$ correspond to phonon absorption and emission channels, respectively, and both channels must be added. In the exciton basis, $\Xi^{0}$ is non-diagonal in the indices $S,S^{\prime}$, but remains diagonal in the wavevector $\mathbf{Q}$. Here, we consider only the optical excitons, and we use $\Xi_{\Gamma S}$ to denote the diagonal elements of $\Xi$. The IEHPP self-energy in the exciton basis is

[TABLE]

and we have introduced the symbols $\Xi^{0}_{\Gamma\mathbf{k},\mathbf{q}}$ and $\Xi^{0}_{\Gamma S,\mathbf{q}}$ to denote an individual $\mathbf{q}$-point’s contribution to the IEHPP self-energy. The exciton-phonon coupling matrix elements are

[TABLE]

and the exciton self-energy is given by

[TABLE]

In general, the completion term for the self-energy contains the contribution of the valence and conduction bands excluded from the two-particle basis, as well as the contribution of the exciton states that are not explicitly computed. Within this model, all the electron bands are included in the basis, but only the first $N$ exciton bands are computed. The completion term thus corresponds to the contribution of the remaining exciton states ($S^{\prime}>N$). We express the completion term in the form SM

[TABLE]

where we define the partial sum

[TABLE]

and the corresponding summation over all exciton states can be obtained with the sum rule

[TABLE]

The ratio $\tilde{\zeta}_{\Gamma S}(\mathbf{q})/\zeta_{\Gamma S}(\mathbf{q})$ thus describes how much the first $N$ excitons probe the space of available states at wavevector $\mathbf{q}$ for the exciton $\Gamma S$ to couple with.

III.3 Results and discussion

Figure 7(a-b) presents the $\mathbf{q}$-space convergence of the imaginary part of the exciton-phonon self-energy. The use of a finite value for the infinitesimal parameter $\eta$ eases the convergence with respect to the number of $\mathbf{q}$-points, but also introduces an arbitrarily small error in the lifetimes. For the first optical exciton ($S=1$ in the exciton labeling), the inverse lifetime due to phonons is known to be zero at $T\!=\!0$, since this exciton cannot scatter into a lower energy states by phonon emission. The finite value obtained for $\text{Im}\Xi_{S}$ with a converged $\mathbf{q}$-grid thus indicates the magnitude of the error. A value of $\eta=10$ meV yields an error smaller than $3$ meV for the self-energy, and we use this value for the following computations of $\Xi$. For the second optical exciton state ($S=2$), we conclude that a $48\times 48$ $\mathbf{q}$-grid is well converged. The temperature-dependent energy shift and the spectral width for the seven lowest exciton states are presented in Fig. 7c. The lowest exciton state, having the largest binding energy, consequently has a longer lifetime than the others at all temperatures.

Figure 8 presents the frequency-dependent self-energy for the second lowest bound exciton. The height of the shaded area represents the completion term, which accounts for a large fraction of the real part of $\Xi$ at all frequencies. Near the $S\!=\!2$ exciton energy ($\Omega_{S}\!=\!2.4$ eV), the imaginary part of $\Xi^{0}$ doesn’t have any structure, since, without electron-hole interactions, the electron-hole pairs lie above the band gap. However, the presence of bound excitons with finite crystal momentum near $\Gamma$ and $K$ allow the optical excitons to scatter and diffuse, confering a finite value to the imaginary part of $\Xi$. Just as the electron-hole interaction binds the excitons below the band gap, it also moves the spectral weight of the self-energy towards lower frequencies.

Let us now evaluate the accuracy of previously-used approximate expressions for the real and imaginary parts of the self-energy. The uncorrelated exciton approximation of Eq. (58) corresponds to writing

[TABLE]

Unlike the full exciton-phonon self-energy, this expression only requires the computation of the exciton’s wave function for the state $\Gamma S$. Figure 9, compares the real and imaginary parts of $\Xi^{\text{UE}}$ to those of $\Xi$. The uncorrelated exciton approximation overestimates the inverse lifetime (broadening) by an order or magnitude. For the phonon-induced energy shift $\Delta\Omega_{S}$, the overestimation of the negative shift (renormalization) for the seven lowest optical excitons is about 20% to 40%.

The combined effect of electron-hole and electron-phonon interactions is summarized in Fig. 10, which shows the exciton propagator for the lowest bound exciton as the electron-hole and electron-phonon interactions are switched on separately ($L$, $\Lambda^{0}$) or simultaneously ($\Lambda$). In all cases, the imaginary part of the self-energy of the particles due to electron-electron interaction is neglected. An artificial broadening parameter ($\eta\!=\!5$ meV) is used to represent the bare exciton propagator ($L$), which would otherwise be infinitely sharp around the exciton energy $\Omega_{S}$ because of the neglect of the lifetime of the particles due to electron-electron interaction . The function $L^{0}(\omega)$ illustrates the spectral decomposition of $L$ into electron-hole pairs. The spiky features in $L^{0}$ are an artifact of the finite $\mathbf{k}$-points sampling used, and the function can be made smooth by using a larger broadening parameter $\eta$. Turning on the electron-phonon interaction broadens these features, as can be seen by comparing $L^{0}$ with $\Lambda^{0}$. The exciton propagator in the presence of the phonons field ($\Lambda$) is red-shifted with respect to the bare exciton propagator ($L$), and a satellite peak appears above the exciton peak in $\Lambda$.

IV Conclusion

In summary, we have derived a rigorous expression for the exciton-phonon coupling self-energy to lowest order in the electron-phonon interaction and in the limit of low exciton density. Through the exciton-phonon coupling matrix elements, the optically accessible excitons may scatter into optically dark finite-momentum exciton states, resulting in an energy renormalization and a finite lifetime for the optical excitations. Our expressions takes into account the electron-hole interaction, and improve upon approximate expressions found in the literature by naturally enforcing energy conservation.

We implemented this theory on a two-dimensional two-band model and computed the temperature-dependent energy shift and lifetimes of the optical excitons. This model allowed us to compare our exciton-phonon self-energy with an approximate expression which we call the uncorrelated exciton (UE) approximation. We showed that the previously used approximation overestimates the inverse lifetime by an order of magnitude, making this approximation unreliable. We conclude that, in physical systems with strong electron-hole interaction such as low-dimensional materials, it is necessary to use the exact exciton-phonon coupling theory to compute accurately the lifetime of optical excitations. The UE approximation also overestimates the temperature-dependent shift of the energy of the excitons by 20-40% in our two-band model. We interpret this result as an upper bound for the error induced by this approximation in realistic systems.

The scheme developed in this paper readily applies to the study of exciton diffusion dynamics. A main challenge in applying this theory is the computation of finite-momentum excitons energies and wavefunctions. While such calculation has been demonstrated, a full sampling of the Brillouin zone remains computationally expensive and would benefit from interpolation techniques Qiu et al. (2015); Berghäuser and Malic (2014); Selig et al. (2018).

Acknowledgements.

G.A. acknowledges valuable discussions with Jonah B. Haber, Sivan Refaely-Abramson, Diana Y. Qiu, and Felipe H. da Jornada. This research was supported by the National Science Foundation under Grant No. DMR-1926004 (conceptual and formalism development), by the Center for Computational Study of Excited-State Phenomena in Energy Materials (C2SEPEM) funded by the U. S. Department of Energy, Office of Basic Energy Sciences, under Contract No. DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory (code development), and by the Natural Sciences and Engineering Research Council of Canada under Grant No. RGPIN-2019-07149 (computation).