Theory of Exciton-Phonon Coupling

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.

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 1-7 . 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 [8][9][10][11][12][13][14] . This dynamics is governed by several processes: the interaction of excitons with defects, the exciton-exciton interaction, and the excitonphonon 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 [15][16][17] .
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 [18][19][20][21][22][23][24][25] . 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) [26][27][28] , 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 [29][30][31] . Going beyond DFT, electron-electron correlation effects to the electron self-energy may further be included at the GW level from ab initio 32 us- * gabriel.antonius@uqtr.ca ing the GW PT method 33 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 [34][35][36] . 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 37 . This scheme has been used to compute the absorption spectrum of h-BN and MoS 2 at finite temperature [37][38][39] . 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 [40][41][42] 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 [43][44][45][46] . 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 47 . Another effect of phonons on the exciton binding energies comes from the lattice screening, and this has been recently computed from first principles 48,49 In this work, we develop a general theory of the excitonphonon 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 excitonphonon 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 twodimensional systems. The main findings are summarized in Sec. IV. Several mathematical details of the derivation can be found in supplemental materials 50 .

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 fixedions Hamiltonian for the electrons H 0 =T k + V SCF (r), wherê T k is the kinetic energy operator, and V SCF (r) is the selfconsistent 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 φ i (r) and energies ε i , where the label i comprises a band index (n i ), a wavevector (k i ), and eventually a spin index (σ i ).
These quantities are used to construct the time-ordered Green's function, defined in the one-particle basis and in time as where c † i and c i are the electron creation and annihilation operators,T t is Wick's time-ordering operator, and "0" 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 and In terms of frequencies, the oneparticle Green's function writes as where ±η is an real infinitesimal number with the same sign as ε 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 manybody 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 electronphonon coupling elements. Such a Hamiltonian, for example, may be constructed from the self-energy computed in the GW formalism 33,51,52 .

A. 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. 24.
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 where v and c refer to the labeling of occupied (valence band) and unoccupied (conduction band) states which involve both band and k-point indices. We call L 0 the bare electron-hole propagator. It is computed with the fixed-ions Hamiltonian, but the superscript "0" 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, where K vc,v ′ c ′ 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 Ω S and electron-hole coefficients A S vc for each exciton S. The resonant part of the bare exciton propagator in the quasiparticle basis can thus be written as 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 (q S ) which is the wavevector difference (k c − k v ) between the unoccupied and occupied singleparticle 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 q S ≈ 0 are optically accessible. We will see that finite-momentum excitons (q S = 0) are important to describe scattering events by phonons.
Since the A S vc 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 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 where v is the velocity operator, 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.

B. Lattice dynamics and electron-phonon interaction
The lattice dynamics can be obtained from a self-consistent calculation of the dynamical matrix of the crystal (Φ), as detailed in Ref. 27. 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: where l labels a unit cell of the crystal with lattice vector R l , κ labels an atom in the unit cell, and j labels a Cartesian direction. The phonon frequencies ω λ and polarization vectors U λ κ j , are obtained by diagonalizing the Fourier transform of Φ as where M κ is the mass of an atom. The label λ for a phonon mode comprises a branch index and a wavevector (q or q λ ). The DFPT method allows for the self-consistent calculation of the first-order derivative of the local potential with respect to atomic displacements, ∇ lκ j V SCF (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 53 .
The phonon propagator is defined in time as qλ +a qλ is the sum of a phonon creation and annihilation operator. In frequency space, the phonon propagator writes as 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. 29 and 54. Expanding the Hamiltonian up to second order in the atomic displacements, the electron-phonon coupling potential writes as with the first-order electron-phonon coupling matrix elements and the second-order matrix elements The derivative of the potential with respect to a phonon mode is defined as The wavevectors (k, q) are written explicitly in Eq. (11) and (12), but will be omitted in the remainder of the paper to lighten the notation. Also, to make the units explicit, we usē h = 1 and the mass M, which normalizes the phonon eigenvectors according to ∑ κ M κ |U λ κ | 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 N q , where N q is the number of wavevectors used to sample the Brillouin zone.
The quantity g ii ′ λ is the electron-phonon coupling matrix element between the one-particle states i and i ′ via the phonon mode λ . Because of crystal momentum selection rule, for g ii ′ λ to be non-zero, we must have k i = k i ′ + q λ + G where G is a reciprocal lattice vector that is non-zero for an Umklapp process. For g (2) ii ′ λ λ ′ , at the lowest order of perturbation theory, we must have k i = k i ′ and q λ ′ = −q λ for the wavevectors, and λ ′ = λ for the branch indices. In practice, the second-order electron-phonon coupling matrix elements g ii ′ λ λ 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 55,56 .

C. 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 τ = it in the interval [−β , β ] with β = 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 and we adopt the convention that in each expectation value . . . , only distinct and connected diagrams should be retained 57 . The time-dependent operators are expressed in the interaction picture, that is with the mean field HamiltonianĤ 0 = ∑ i ε i c † i c i . Equation (14) can be cast into a Dyson equation for G, depicted in Fig. 2, as which defines the one-particle electron-phonon self-energy Σ.
The perturbative expansion of G with Eq. (14) yields powers of V The FM term is dynamic (i.e. complex and frequencydependent) and stems from the first-order electron-phonon coupling potential. Its analytic expression is where N B (ω) is the Bose-Einstein distribution, and f (ω) 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ω n → ω + iη where η 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 frequencyindependent), and it stems from the second-order electronphonon coupling potential. It writes The self-energy allows for the mixing of electronic states through the off-diagonal components of the self-energy (i = i ′ ). This mixing is only possible among states with the same crystal momentum (k i = k i ′ ). If the bands are well-separated in energy, one may use only the diagonal elements of the selfenergy to obtain the renormalization of the bands as and in general, ∆ε i (0) = 0. Correspondingly, the inverse lifetime of the electronic state is 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 58 .

II. EXCITON-PHONON SELF-ENERGY: SELF-ENERGY OF EXCITONS FROM COUPLING TO PHONONS
We seek similar expressions for the exciton energy (Ω S ) and lifetime (τ 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 (Ξ SS ) to be discussed below, namely, for the exciton energies, and 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 (Λ) defined as where τ + is a time infinitesimally larger than τ, and the interaction termĤ 1 (τ) is the sum of the static electron-hole interaction kernelK and the electron-phonon coupling potential V ep . Equation (24) reduces to the BSE ifV ep = 0. The exciton-phonon self-energy connects the interacting exciton propagator Λ to the bare exciton propagator L (without electron-phonon interaction) in a Dyson-like equation: (25) 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.

A. Independent electron-hole-pair-phonon self-energy
We define Λ 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 The corresponding IEHPP self-energy Ξ 0 is defined through the Dyson equation which is depicted in Fig. 3. A detailed derivation of Ξ 0 is provided in supplemental materials 50 . 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 typi- The different contributions to the IEHPP self-energy are grouped in three terms, according to the topology of the diagrams depicted in Fig. 4 Ξ The Fan-Migdal term is further split in two contributions: the dynamic term (FMd) and the static term (FMs), 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 Performing the convolution of the bare electron-hole propagator L 0 with the phonon propagator D, we obtain the analytic expression and one can safely assume that The static FM term has the same topology as the dynamic FM term. It is given by and its analytic expression is where the symbol S means that the terms in brackets should be symmetrized with the substitution v, c ↔ v ′ , c ′ 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 The analytic expression for this term is 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 Finally, the Debye-Waller contribution to Ξ 0 is simply the second-order interaction of the electron and the hole with the phonon modes, giving Consider a non-interacting electron-hole excited state with energy ε c − ε v . We can show that the diagonal component of the IEHPP self-energy for this state is 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 The broadening of a non-interacting electron-hole transition is thus the sum of the broadenings of the one-particle states.

B. Exciton-phonon self-energy
The exciton-phonon self-energy Ξ is obtained by expanding the interacting exciton propagator Λ in the bare exciton basis as 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 oneparticle 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 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 59 we can take The bare exciton propagator is diagonal in the exciton basis and is formally defined as 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 with the two-particle HamiltonianĤ 2p = ∑ S Ω S c † S c S . The perturbation terms V ep = V and Here, as defined above, A λ is a sum of a phonon creation and annihilation operator, (not to be confused with the electronhole coefficients A S vc ). The first-and second-order excitonphonon coupling matrix elements are the transition amplitude between two exciton states mediated by the electron-phonon coupling potential, e.g.
Using Eq. (40) and (1) to simplify the expectation values and discarding terms like δ SS ′ ∑ g vvλ which produce disconnected diagrams, we arrive at the expressions and g (2) Momentum conservation dictates that q S = q S ′ + q λ for g SS ′ λ to be nonzero, while q S ′ = q S for g (2) SS ′ λ λ . It is worth writing explicitly all the wavevectors in these expressions. For an exciton S ′ with center-of-mass momentum q S ′ = Q and an exciton S with center-of-mass momentum q S = Q + q, the first-order matrix element writes as where A SQ vck is the coefficient of the exciton envelope function in k-space for a hole at k and an electron at k + Q. We note that (Q, 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 ′ . 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 Q are Note that the second-order matrix elements do not transfer momentum between the exciton states, even when q = 0. We may now compute the self-energy from Eq. (39). Some of the diagrams using quasiparticle (electron/hole) basis for the exciton-phonon self-energy up to second order in electronphonon coupling are illustrated in Fig. 5. The Debye-Waller diagram coming from the second-order matrix element is 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 Ξ DW SS ′ = Ξ 0 DW SS ′ . The first-order coupling matrix elements yield the Fan-Migdal term. In the perturbative expansion of Λ, the first term of Eq. (39) produces only the dynamic part of the FM term, since only this part can be expressed in terms of g SS ′ λ . It is given by with the analytic expression 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. (53) reduces to Eq. (31).
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 Ξ.
All other diagrams of Λ in which the scattered state cannot be projected onto the exciton basis are contained in Λ 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.(33) so that with Ξ FMs SS ′ = Ξ 0 FMs SS ′ . 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 Λ 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 Ξ FMd , which we call the completion term, written as Ξ 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 Ξ C from the expression of Ξ 0 , where the intermediate states v ′′ and c ′′ 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 excitonphonon self-energy reads Thus, Ξ contains all the diagrams of Ξ 0 but with addition of the electron-hole interaction diagrams, as shown in Fig. 5. Note that Ξ C SS ′ 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 (Ω S ) we may consider that only the dynamic Fan-Migdal term contributes to the imaginary part of the self-energy, which is The temperature renormalization of the exciton energies and their lifetime can finally be computed according to Eqs. (22) and (23).

C. 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 (Ω S , T ) ≪ N B (ω λ , T ) and write for the diagonal part of the self energy 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: This approach is equivalent to the one used in previous studies [37][38][39] . 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 Ξ 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. (31), the imaginary part of Ξ UE S writes as Compare this expression with Eq. (56), 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 (ω λ , T ) → 0. For the lowest bound exciton, from Eq. (56), 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 60 . 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 61 , which mimics the spin-dependent band dispersion of transition metal dichalcogenides near the K and K ′ valleys, and yields realistic exciton binding energies. The main equations are reported below, and we refer the reader to Ref. 61 for further details of this model. A. One-particle and two-particle Hamiltonians The one-particle tight-binding Hamiltonian is H = ∑ nσ H nσ with where c † nσ R and c nσ R are creation and destruction operators for an electron in the n th orbital on the lattice site R with spin σ .
The first term describes the hopping between a site and one of its six closest neighbours, where δ δ δ 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 parametert n according to t nσ δ δ δ = t n + 4iσt n sin(K · δ δ δ ) (61) where σ = ±1/2 and K is one corner of the 2D hexagonal Brillouin zone. The second term of Eq. (60) is diagonal in lattice site, with ε n being the on-site parameters. The solutions for the one-particle energies at wavevector k are 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 ′ valleys, with band effective masses m * and a splitting ∆ between different spin bands. The on-site and hopping parameters are thus with t = 2h 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 ∆ = 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 ′ 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 welldefined center-of-mass momentum Q, and the BSE Hamiltonian writes σ R creates an electron at R e = R with spin σ , and a hole at R h = 0 with spin −σ , meaning that no spin-flip occurs in this process. The BSE hopping parameter carries the information on the exciton momentum, and writes T σ Qδ δ δ = t cσ δ δ δ e −iQ·δ δ δ /2 − t vσ δ δ δ e iQ·δ δ δ /2 .
Together, the on-site terms that involve ε c and ε 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 where ε is the dielectric constant of the medium and ∆v 0 is the difference between the screened Coulomb interaction and the bare exchange integral at R = 0, which we set to ∆v 0 = 1.6 eV. Note that, for this model, the spin degrees of freedom do not mix, and at Γ, 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 σ will thus be omitted. Due to the sparse nature of a Hamiltonian with nearestneighbour hopping, the BSE can be solved efficiently in real space for all Q-vectors. Solving the BSE yields the exciton energies Ω QS and wavefunctions in real space A S (R, Q). Their Fourier transform give the electron-hole coefficients A S (k, 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).

B. Exciton-phonon coupling
We model the lattice vibration spectrum as a single dispersionless phonon band with frequency ω 0 . The electronphonon 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 ′ = g vv ′ = 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 ω 0 = 50 meV and g = 250 meV.
The IEHPP self-energy for a hole at k and an electron at k + Γ is diagonal in k indices and is given by where P + (T ) = N B (ω 0 , T ) and P − (T ) = N B (ω 0 , T ) + 1 correspond to phonon absorption and emission channels, respectively, and both channels must be added. In the exciton basis, Ξ 0 is non-diagonal in the indices S, S ′ , but remains diagonal in the wavevector Q. Here, we consider only the optical excitons, and we use Ξ ΓS to denote the diagonal elements of Ξ. The IEHPP self-energy in the exciton basis is and we have introduced the symbols Ξ 0 Γk,q and Ξ 0 ΓS,q to denote an individual q-point's contribution to the IEHPP self-energy. The exciton-phonon coupling matrix elements are and the exciton self-energy is given by 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 ′ > N). We express the completion term in the form 50 where we define the partial sum and the corresponding summation over all exciton states can be obtained with the sum rule The ratioζ ΓS (q)/ζ ΓS (q) thus describes how much the first N excitons probe the space of available states at wavevector q for the exciton ΓS to couple with.
C. Results and discussion Figure 7(a-b) presents the q-space convergence of the imaginary part of the exciton-phonon self-energy. The use of a finite value for the infinitesimal parameter η eases the convergence with respect to the number of 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 ImΞ S with a converged q-grid thus indicates the magnitude of the error. A value of η = 10 meV yields an error smaller than 3 meV for the self-energy, and we use this value for the following computations of Ξ. For the second optical exciton state (S = 2), we conclude that a 48 × 48 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 Ξ at all frequencies. Near the S = 2 exciton energy (Ω S = 2.4 eV), the imaginary part of Ξ 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 Γ and K allow the optical excitons to scatter and diffuse, confering a finite value to the imaginary part of Ξ. 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 Unlike the full exciton-phonon self-energy, this expression only requires the computation of the exciton's wave function for the state ΓS. Figure 9, compares the real and imaginary parts of Ξ UE to those of Ξ. The uncorrelated exciton approximation overestimates the inverse lifetime (broadening) by an order or magnitude. For the phonon-induced energy shift ∆Ω 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 electronhole and electron-phonon interactions are switched on separately (L, Λ 0 ) or simultaneously (Λ). In all cases, the imaginary part of the self-energy of the particles due to electronelectron interaction is neglected. An artificial broadening parameter (η = 5 meV) is used to represent the bare exciton propagator (L), which would otherwise be infinitely sharp around the exciton energy Ω S because of the neglect of the lifetime of the particles due to electron-electron interaction . The function L 0 (ω) illustrates the spectral decomposition of L into electron-hole pairs. The spiky features in L 0 are an artifact of the finite k-points sampling used, and the function can be made smooth by using a larger broadening parameter η. Turning on the electron-phonon interaction broadens these features, as can be seen by comparing L 0 with Λ 0 . The exciton propagator in the presence of the phonons field (Λ) is red-shifted with respect to the bare exciton propagator (L), and a satellite peak appears above the exciton peak in Λ.

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 twoband 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 [62][63][64]

I. QUASIPARTICLE PROPAGATORS AND FREQUENCY SUMMATIONS
This section presents the Matsubara formalism used and a few important Matsubara summations that lead to the analytic expressions of the different self-energies presented in the manuscript.
The transform between imaginary time and imaginary frequency for a quantity A (either a propagator or a self-energy) is where β = 1/k B T , and the Matsubara frequencies are defined as Throughout the paper, we use iω n and iω l for bosonic frequencies, and iω m for fermionic frequencies. We will recover all the retarded quantities on the real frequency axis using a real positive infinitesimal number η and the analytic continuation In this formalism, the one-particle propagator is while the phonon propagator is And the bare exciton propagator is For the phonon propagator, the summation over bosonic Matsubara frequencies yields where N B (ω) is the Bose-Einstein distribution defined as For the one-particle propagator, the summation over fermionic Matsubara frequencies yields where f (ω) is the Fermi-Dirac distribution defined as assuming that ω is measured from the chemical potential µ. The bare electron-hole propagator is expressed as the convolution of two one-particle propagators Upon computing the one-particle electron-phonon self-energy, we encouter a summation of the while, for the IEHPP self-energy Ξ 0 , we encouter summations of the type The last approximation results from the large-band-gap approximation, which is explained in appendix II.

II. LARGE BAND GAP APPROXIMATION
In order to arrive at an expression for all the contributions to Ξ 0 , we will take advantage of the large-band-gap approximation for the one-particle occupation numbers, namely For materials with band gaps larger than 1 eV, this approximation holds for temperatures up to 10, 000 K. As a result, the bare electron-hole propagator writes And we have the identity Using this result, we can show that 19) which has been symmetrized in (v, v ′ ). The third line of Eq. (S. 19) gives a contribution that is exactly zero if it is summed along a frequency-independent term, such as the DW self-energy. Its contribution is also zero when evaluated at iω n = ε c − ε v , and is vanishing in the neighborhood of this value. It will thus be neglected. To simplify the notation, we introduce the symmetrization symbol S, which, for a quantity A vc,v ′ c ′ , means Quantities with two or three indices are obtaind by a contraction of indices and symmetrize as Hence, we arrived at the approximation A similar reasonning gives Another result from the large-band-gap approximation is Finally, we assume that all phonon frequencies are much smaller than the band gap and the exciton energies. This approximation allow us to write and for any pair of valence and conduction states with energies ε v and ε c ), any exciton energy Ω S , and any phonon energy ω λ .

III. MORE DETAILED DERIVATION OF Ξ 0
We proceed to a detailed derivation of the IEHPP self-energy Ξ 0 by expanding Λ 0 from Eq. (??). To the lowest order in the perturbation, the self-energy will be given by Hence, each term in Λ 0 beyond L 0 contributes to Ξ 0 .
The first set of terms are called the Debye-Waller terms, since they involve the second-order electron-phonon coupling potential. They write We recover Eq. (??) for Σ DW by making use of Eq. (S.8). The next set of terms are called the Fan-Migdal (FM) terms. They write We recover Eq. (??) for Σ FM by making use of Eq. (S.13). The last set of terms are called the phonon exchange terms. They write Using Eq. (S.19) and Eq. (S.22), the Debye-Waller contribution to Ξ 0 is Within the large-bandgap-approximation, the Debye-Waller term is simply To obtain the FM term contribution, we insert Equations (S.31) into (S.28), and use (S.22) and (S.23), giving We split the FM term into two contributions: the dynamical FM term (FMd) and the static FM term (FMs). We define the dynamical Fan term as 37) and the static FM term as We use Eq. (??) to express the one-particle self-energy as a convolution of the Green's functions of an electron and a phonon, and use Eq. (S.12) to write and × g v ′ c ′′ λ g * vc ′′ λ δ cc ′ Sδ vv ′′ + g cv ′′ λ g * c ′ v ′′ λ δ vv ′ Sδ cc ′′ (S.40) To obtain the phonon exchange term contribution, we insert Equations (S.33) into (S.28), and use (S.24) and (S.25), giving Ξ 0 X vcv ′ c ′ (iω n ) = g c ′ cλ g * v ′ vλ 1 β ∑ l L 0 vc ′ ,vc ′ (iω n + iω l ) + L 0 v ′ c,v ′ c (iω n + iω l ) D λ (iω l ) (S.41) Performing the imaginary frequency summation, we obtain Eq. (??).

IV. COMPLETION TERM WITHIN THE TWO-BAND MODEL
In the present calculation, the wavevector basis allows for a maximum of N k exciton states, but the first N exciton states are actually computed, The exciton-phonon self-energy therefore includes a completion term to account for the missing transitions. Let |k, q denote the state with one hole at wavevector k and one electron in the conduction band at wavevector k + q, such that k, q|S, q = A S (k, q). The electron-phonon coupling operator is V ep (q) = ∑ k g c |k, q k, Γ| − g v |k − q, q k, Γ| (S.42) We complete the self-energy using the equivalent completion phase space of the non-interacting self-energy, that is In the exciton basis, the non-interacting self-energy is non-diagonal, but we consider a diagonal element (S, S) Ξ 0 SS (Γ, ω) = ∑ k ′ ,q g SS ′ (Γ, q)g * SS ′′ (Γ, q)A S ′ * (k ′ , q)A S ′′ (k ′ , q) ω − (ε k ′ +qc − ε k ′ v ) ± ω 0 + iη P ± (T ) Given that ∑ k A S ′ * (k ′ , q)A S ′′ (k ′ , q) = δ S ′ ,S ′′ (S. 43) in the limit of weak exciton binding, Ξ 0 SS must correspond to Ξ SS , namely The completion term is the missing terms, when only the first N excitons are included in the basis, that is |g SS ′ (Q, q)| 2 P ± (T ) ω − Ω S ′ ± ω 0 + iη We first approximate the contribution of the S ′ states as those of electron-hole pairs |g SS ′ (Q, q)| 2 P ± (T ) ω − Ω S ′ ± ω 0 + iη Next, define a q-dependent effective energy parameter E e f f (S, q) that allows to write the completion term as |g SS ′ (Q, q)| 2 ω − E e f f (S, q) + iη (S. 44) In the limit of complete basis (N = N k ), each pair of terms in Eq. (S.44) must cancel perfectly. In