We employ quasiparticle path integral molecular dynamics to study how the excitonic properties of model semiconductors are altered by electron–phonon coupling. We describe ways within a path integral representation of the system to evaluate the renormalized mass, binding energy, and radiative recombination rate of excitons in the presence of a fluctuating lattice. To illustrate this approach, we consider Fröhlich-type electron–phonon interactions and employ an imaginary time influence functional to incorporate phonon-induced effects nonperturbatively. The effective mass and binding energies are compared with perturbative and variational approaches, which provide qualitatively consistent trends. We evaluate electron-hole recombination rates as mediated through both trap-assisted and bimolecular processes, developing a consistent statistical mechanical approach valid in the reaction limited regime. These calculations demonstrate how phonons screen electron–hole interactions, generically reducing exciton binding energies and increasing their radiative lifetimes.

## I. INTRODUCTION

The application and design of photovoltaic devices rely on understanding the photophysics of semiconducting materials. Recent studies into novel low dimensional and hybrid perovskite semiconductors have highlighted the need to incorporate the effects of a fluctuating lattice on the stationary behavior of excitons, moving away from the traditional perspectives in which screening is presumed to be largely determined by electronic degrees of freedom.^{1–7} While the study of electron–phonon coupling for free charges has a long history, including foundational studies on polarons,^{8–15} there is comparatively little known concerning the effects of electron–phonon coupling on excitonic properties. Motivated by observations that suggest polaronic effects play an important role in renormalizing exciton mobilities,^{16} binding energies, recombination rates,^{17,18} and photoluminescence yields,^{19} we aim to fill this knowledge gap. In this work, we explore the effects of phonons on the excitonic properties of traditional and hybrid perovskite materials using a path integral approach. Working within a Fröhlich model Hamiltonian,^{8} we evaluate numerically exactly the role of phonons, finding that they generally reduce exciton binding energies and increase radiative lifetimes.

The static properties of excitons determine the power conversion efficiencies of photovoltaics, the quantum yields for light emission, and more. Predicting these properties from molecular models is an area of active development.^{20} Most widely used approaches build upon ground state density functional theory, employing corrections from many-body physics, including the GW approximation and Bethe–Salpeter equations.^{21–23} These and related approaches^{24} have been successful for a wide range of semiconducting materials.^{25,26} However, these theories traditionally ignore dynamical effects from phonons, as including them within this framework is challenging. While historically, analytical approaches based on model Hamiltonians have been developed,^{27–29} recent efforts have focused on numerical methods to describe the effects of phonons approximately.^{30–33} While these approaches leverage powerful *ab initio* many body theories, they have been limited in the strength of the electron–phonon coupling that can be considered and the time and lengthscales approachable.

Here, we present a path integral approach for describing electron and hole quasiparticles interacting with phonons. Analogous quasiparticle path integral approaches have been utilized to describe photoinduced phase separation,^{34–36} charge trapping,^{37,38} and charge recombination^{39} as well as confinement effects.^{40} With this approach, we can derive an imaginary time influence functional, allowing us to incorporate dynamical and quantum mechanical effects of phonons within a harmonic approximation. To sample the resultant theory, we apply path integral molecular dynamics (MD) and study how electron–phonon coupling renormalizes the band mass, exciton binding energy, and electron–hole recombination rate for Fröhlich-type interactions between electron and hole quasiparticles and phonons. In the following, we elaborate theoretical details of the path integral approach applied to exciton with phonons and present how to calculate each property within this framework, followed by a discussion on the effects of phonons.

## II. THEORY

We consider a system composed of an electron and a hole interacting with a field of phonons, whose effective Hamiltonian consists of three parts,

a part due to the electronic degrees of freedom, $H\u0302eh$, a part from the lattice, $H\u0302ph$, and their interaction, $H\u0302int$. The electronic part includes the kinetic energies of an electron and a hole and a Coulomb interaction,

where $p\u0302$ and $x\u0302$ are momentum and position operators of a quantum particle, *ɛ*_{r} is the dielectric constant in units of vacuum permittivity *ɛ*_{0}, and the subscripts *e* and *h* indicate the electron and hole. The masses *m*_{e} and *m*_{h} are taken as their corresponding band masses using an effective mass approximation. This simplification can be relaxed by parameterizing more elaborate kinetic energy functions with position dependent masses.^{40}

The lattice is described by a collection of harmonic modes

where $p\u0302k$ and $q\u0302k$ are the mass weighted momentum and coordinate of a phonon at wavevector **k**, respectively. Without loss of generality, we will take the frequency of the oscillators to be constant *ω*_{k} = *ω* and equal to the longitudinal optical mode. While previous work has illustrated the importance of including additional modes or their wave-vector dependence in specific materials,^{33,39} we neglect these effects here in order to benchmark the approach to a simplified model. For generalizations, see Appendix A.

We adopt a Fröhlich-type interaction^{8,41} between the charges and the phonons, where a charged particle interacts linearly with the polarization field produced by a lattice vibration,

where $q\u0302k$ corresponds to the polar displacement field that the charge can be coupled to along the **k** direction. The strength of the coupling is set by a material specific constant

where *i* indicates either the electron or hole, *V* is the volume of the system, *ℏ* is Planck’s constant divided by 2*π*, and *α*_{i} = *α* is a dimensionless Fröhlich coupling constant.^{41} To study this system, we employ a path integral formalism,^{42,43} which allows us to describe the correlated behavior of the electron, hole, and phonons quantum mechanically and on an equal footing.^{44} The partition function of the system can be written as

where the path action $S$ is defined as

with the imaginary time variable *τ* and *β*^{−1} = *k*_{B}*T*. The Hamiltonian indexed by *τ* represents the classical counterpart of Eqs. (2)–(4) at given *τ*, where $p\u0302$ and $x\u0302$ are replaced by **p**_{τ} and **x**_{τ}, respectively, for each quasiparticle.

Considering that the phonons act as a Gaussian field coupled linearly to the charge density, the phonon variables {**q**_{k}} can be integrated out,^{39,42} yielding

with $Zph$ being the partition function of phonons without the charge. The resulting effective Hamiltonian at a given *τ* can be written as a sum of four pieces,

with

where $\alpha ij=\alpha i\alpha j$, $mij=mimj$, and *σ*_{ij} is +1/−1 for the same/opposite charges. We have used the inverse Fourier representation of **k**^{−2} with the equation for *C*_{i/j} given in Eq. (5). In order to use these results computationally, we discretized the imaginary time interval into *n* slices. In this discrete formulation, the effective Hamiltonian $H=Heh+Heff$ becomes

where we denote the last sum $HC$ and

where *i*, *j* ∈ {*e*, *h*} and *t*, *s* ∈ [1, *n*]. The number of timeslices, or beads, is a convergence parameter that needs to be taken large for accuracy.^{45}

In the path integral framework, Eq. (11) implies that electron and hole quasiparticles are represented as classical ring polymers^{46} consisting of *n* identical beads, where adjacent beads are harmonically coupled and beads with the same index from the electron and hole interact through a fraction of Coulomb potential. Additionally from Eq. (12), the effective energy induced by phonons depends on the positions of two different imaginary times, represented as the interaction between beads, and the stiffness of phonons sets the decaying imaginary timescale. In this way, we employ an imaginary time influence functional formalism, which is shown schematically shown in Fig. 1. Two effects from the phonons are clear from this picture. First, individual charges are localized by the phonons due to the induced attractive self-interaction. Second, the electron–hole interaction is weakened due to the induced repulsion, which is a reflection of phonon screening. The implications of these two effects are explored below.

## III. SIMULATION DETAILS

To study the utility and efficiency of this approach, we consider models motivated by CdS and MAPbI_{3}, whose material properties are in different regimes of band mass and electron–phonon coupling strength. CdS is a traditional II–IV semiconductor, where the material specific Fröhlich coupling constant *α* is small (1.3) and the band mass of the hole is much heavier than the mass of electron. For MAPbI_{3}, a lead–halide perovskite, the electron and the hole have nearly the same band mass and the coupling strength is intermediate *α* = 2.87. For the purpose of this paper, since we aim to use the path integral method to study the general effects of electron–phonon coupling, we will ignore the anharmonic corrections from the lattice, which can be important in determining the optoelectronic properties of these materials.^{6,17–19,39} Parameters used in simulations are summarized in Table I.

Parameter (unit) . | CdS . | MAPbI_{3}
. |
---|---|---|

Electron band mass m_{e} (m_{0}) | 0.19^{47} | 0.20^{48} |

Hole band mass m_{h} (m_{0}) | 0.80^{47} | 0.20^{48} |

Optical frequency ω (THz) | 9.14^{49} | 7.53^{50} |

Dielectric constant ɛ_{r} (ɛ_{0}) | 5.7^{47} | 6.1^{48} |

Bandgap E_{gap} (eV) | 2.58^{47} | 1.64^{51} |

Fröhlich constant α_{e} | 1.3 | 2.87 |

Parameter (unit) . | CdS . | MAPbI_{3}
. |
---|---|---|

Electron band mass m_{e} (m_{0}) | 0.19^{47} | 0.20^{48} |

Hole band mass m_{h} (m_{0}) | 0.80^{47} | 0.20^{48} |

Optical frequency ω (THz) | 9.14^{49} | 7.53^{50} |

Dielectric constant ɛ_{r} (ɛ_{0}) | 5.7^{47} | 6.1^{48} |

Bandgap E_{gap} (eV) | 2.58^{47} | 1.64^{51} |

Fröhlich constant α_{e} | 1.3 | 2.87 |

For both materials, we use MD simulations to sample the effective actions with fictious masses for the beads kept at 1 amu, and we study the renormalization of the effective mass, exciton binding energy, and recombination rate due to electron–phonon coupling. For the effective mass calculations, $Heffee$ defined in Eq. (12) is used as the system Hamiltonian, and for the other two properties, we run simulations of an electron–hole pair described as two ring polymers where one has a unit negative charge for an electron and the other has a unit positive charge for a hole with the Hamiltonian $H=Heh+Heff$ given in Eqs. (11) and (12).

To avoid the divergence in the 1/|**x**| term between attractive beads, a pseudopotential is used, where 1/|**x**| is replaced by $(rc2+x2)\u22121/2$ and *r*_{c} is chosen to reproduce the bandgap of each material.^{52} Simulations are run in an ensemble with constant volume, particle number, and temperature using a Langevin thermostat, where the total momentum averages to zero with an integration time step 1.0 fs and at room temperature unless explicitly specified, using the LAMMPS^{53} package. In Secs. IV–VI, we present a way to compute each property within the path integral framework, followed by the discussion on the effects from phonons, where *α* serves as a control parameter for the interaction strength between charges and phonons.

## IV. EFFECTIVE MASS

To validate the path integral framework, we first study the effective mass for which significant previous analysis has been undertaken.^{43,55} Since the presence of a charge induces a distortion to the lattice, motion of the charge requires moving the corresponding distortion field and renormalizes the mass of the charge, making it heavier. Within the path integral framework, since the effective mass is a momentum-dependent quantity, an inverse effective mass can be computed with an open-chain ring polymer in the low temperature limit as^{55,56}

where ⟨⋯⟩ denotes the ensemble average and Δ*r* is the distance between the first and the last beads, schematically shown in the inset of Fig. 2.

For simplicity, a Feynman unit system where *k*_{B} = *ℏ* = *ω* = *m* = 1 is used in this calculation. Given that Eq. (13) is valid at low temperature, we tune the temperature and pseudopotential parameters and set *k*_{B}*T*/*ℏω* = 0.02 and *r*_{c} = 0.0707 that are low and small enough for the convergence of our results. Figure 2 describes the inverse effective mass of electron and hole quasiparticles of CdS at different *α*, where each point is the value averaged over 20 ensembles from simulations performed in constant volume, particle number, and temperature. We find that the lighter mass and the stronger coupling require larger number of beads to converge. The results from the path integral approach are consistent with the value predicted from first order perturbation theory for small *α*,^{15}

## V. EXCITON BINDING ENERGY

Within the path integral framework, we capture the full correlation energy between the electron and hole quasiparticles, allowing us to compute accurate exciton binding energies. Here, we show how electron–phonon coupling alters the binding energy. We compute the exciton binding energy from the average energy of exciton,^{57}

resulting in two pieces, the average kinetic energy ⟨*E*⟩_{K} and the average potential energy ⟨*E*⟩_{P}. For the kinetic energy, since the relevant terms produced by Eq. (15) diverge as *n* → ∞, we use a virial estimator,^{45} which is known as an efficient way to estimate the kinetic energy in path integral simulations to avoid the large fluctuations from the subtraction of two diverging terms. Using the derivative of potential energy, the average kinetic energy can be written as

and the average potential energy becomes

where $\u27e8Heffij\u27e9\u2032$ is defined as $\u27e8Heffij\u27e9$ given by Eq. (12) with an additional factor of *βℏω*|*t* − *s*|/*n* inside the summations with respect to *t* and *s*.

The exciton binding energy is the energy threshold for an optical absorption between conduction and valence bands and, thus, traditionally reported in the low temperature limit. To evaluate it, we compute the average energy difference between the exciton and separately the electron and hole at different temperatures and extrapolate to zero temperature, *E*_{B} = lim_{T→0}⟨*E*⟩_{ex} − ⟨*E*⟩_{e} − ⟨*E*⟩_{h}. The subscript ex denotes an average with both the electron and hole, while *e*/*h* refers to calculations of the self-energies of the electron/hole, where the quasiparticles interact only with surrounding phonons.

Figures 3(a) and 3(b) show the average energy difference at different temperatures for CdS and MAPbI_{3}. In both cases, we consider *α* = 0 and compare to *α* > 0 and find that *n* = 200 is large enough to converge the result. The extrapolated exciton binding energies are summarized in Table II. We find that exciton energy becomes lower at high temperature, reflective of the higher population of phonons to stabilize the free charges. In addition, as the electron–phonon interaction becomes stronger, the exciton binding energy becomes smaller, implying that the phonons screen the effective electron–hole interaction.

CdS α | 0.0 | 0.53 | 1.5 | 3.0 |

E_{B} (meV) | 70.8 | 65.0 | 54.3 | 37.4 |

MAPbI_{3} α | 0.0 | 2.87 | ⋯ | ⋯ |

E_{B} (meV) | 40.6 | 25.8 | ⋯ | ⋯ |

CdS α | 0.0 | 0.53 | 1.5 | 3.0 |

E_{B} (meV) | 70.8 | 65.0 | 54.3 | 37.4 |

MAPbI_{3} α | 0.0 | 2.87 | ⋯ | ⋯ |

E_{B} (meV) | 40.6 | 25.8 | ⋯ | ⋯ |

For the results without phonons, *α* = 0, we compare our path integral simulations with binding energies from the Wannier–Mott exciton^{58} that is equivalent to our system.^{59} The exciton binding energy is computable exactly and given by $EBH=\mu e4/2(4\pi \epsilon r)2\u210f2$, with *μ* = *m*_{e}*m*_{h}/(*m*_{e} + *m*_{h}) being the reduced mass of an exciton. Calculated hydrogenic binding energies are 64.3 and 36.6 meV for CdS and MAPbI_{3}, consistent within a few meV with our values in both materials, showing the robustness of the calculation with the path integral framework. In the presence of electron–phonon interaction with a nonzero value of *α*, we compare the binding energies from the corresponding coupling strength, *α* = 1.3 for CdS, and *α* = 2.87 for MAPbI_{3}, with the prediction from perturbation theory and the Pollmann–Buttner theory.^{28} For CdS, we interpolated the results in Table II. The approximated differences in binding energies from the first order perturbation theory, $\Delta EBF=EB,\alpha =0\u2212EB,\alpha \u22600\u2248\u22122\alpha \u210f\omega $, are 15.7 and 28.5 meV for CdS and MAPbI_{3}, respectively, which are higher than our results, 14.4 meV for CdS and 14.8 meV for MAPbI_{3}. The Pollmann–Buttner theory results from a canonical transformation of the original Hamiltonian, which in the weak electron–phonon coupling limit provides an effective potential between the electron and hole. Binding energies estimated by solving the Schrödinger equation using the Pollmann–Buttner potential between the electron and hole in Appendix B are 41.3 and 16.7 meV for CdS and MAPbI_{3}, respectively, lower than the results from the path integral simulations. We suspect that the differences are attributed to the neglect of charge density relaxation due to hybridization with the phonons in perturbation theory and the Pollmann–Buttner theory. This is likely a larger effect in MAPI_{3} because the equal masses of the electron and hole render phonon screening at short distances a higher order process. Considering the fact that the anharmonicity from the lattice is not taken into account, the path integral estimates of the binding energies are in reasonable agreement with typical experimental values of 28^{60} and 16 meV^{61} for CdS and MAPI_{3}, respectively.

## VI. ELECTRON–HOLE RECOMBINATION RATE

We now investigate the electron–hole recombination rate that typically dominates the lifetime of charge carriers in bulk semiconducting materials.^{20} It has been generally accepted that charge carrier recombination can be divided into three different mechanisms.^{17,61} The first is trap-assisted recombination where one charge carrier is trapped by a defect or an impurity and then recombination occurs from the trap state. The second is due to bimolecular recombination where an electron in the conduction band is combined with a hole in the valence band, which is a predominant radiative pathway for direct semiconductors under standard operating conditions.^{62} Finally, Auger recombination in which the electron and hole are recombined through an energy transfer to other charge carriers or phonons is yet a higher order process. While Auger recombination can be important in nanocrystals, it is typically negligible in bulk materials.^{63,64}

All the mechanisms described above contribute to the recombination process but depending on the concentration of charge carriers, the dominant pathway varies. At low density, trap-assisted recombination will dominate; at high density, Auger recombination will matter. In this work, we assume that the contribution from Auger recombination is small and consider the first and second order recombination processes. For a charge neutral system, the electron and hole concentrations are equal *ρ*_{e} = *ρ*_{h}, and the rate equation is given by

where *k*_{m/b} is the rate constant for the trap-assisted/bimolecular recombination process with *k*_{tot} being the overall rate constant.

Independent of the mechanism, the radiative recombination rate of an electron–hole pair within the path integral framework is given by^{65}

originating from Fermi’s golden rule for spontaneous emission under our effective mass approximation, where *c* is the speed of light and *E*_{gap} is the bandgap energy. The other parameters for CdS and MAPbI_{3} are summarized in Table I.

The rate in Eq. (19) is proportional to a ratio of partition functions whose subscripts 0 and *c* indicate the standard thermal trace and the trace for a *radiating path* where two separate imaginary time paths for quasiparticles are combined at a common imaginary time schematically shown in Fig. 4. This ratio is identical to a thermally averaged overlap integral between the electron and hole densities. The partition functions are related through

where $\Delta S$ is the change in path action. It is convenient for sampling purposes to rewrite this average using a conditional probability representation,

so that the ratio in the dilute limit can be written as

where **R** is the vector between electron and hole ring polymers, represented by *βF*(*R*) = −ln *P*(*R*), a free energy for changing their distance. The difference in path actions is equal to

for going between the thermal and radiating paths. The first two terms correspond to $S0$, and $Sc$ is given by the last term, and the index of beads for the combined coordinates {**x**_{c}} is schematically described in Fig. 4. In the following, we present the details on simulations and discussions on trap-assisted and bimolecular recombination rates and combine these two to estimate a total rate constant. In both, we assume that recombination is not limited by diffusion of the charge carriers so that a local equilibrium distribution is established for the relative positions of electrons and holes. For both, we will evaluate the rate at room temperature *T* = 298 K.

### A. Trap-assisted recombination rate

For trap-assisted recombination, we need to describe the trapping of a charge as well as the subsequent recombination of electron and hole quasiparticles. To describe this process, we assume that the trapped charge achieves a steady state population, and the rate is given by the likelihood of finding a trapped charge times the rate to recombine that trapped charge with an incoming charge,

where *P*_{trap} is the probability of an electron to be trapped. This approximation is valid in the dilute limit, provided trapping is reversible. For both CdS and MAPbI_{3}, we consider the trapping of an electron with a positively charged point defect. The point defect is described by a Coulomb potential acting between the defect and charge, determined by the corresponding dielectric constant *ɛ*_{r} and *r*_{c} that is set to recover the reported trapping energy, 1.75^{66}(0.3^{67}) eV, for CdS (MAPbI_{3}).

In equilibrium, the probability of finding an electron trapped by an isotropic point defect is given by

where *ρ*_{d} is the density of defect sites in the lattice where it is possible to trap an electron, $Rd*$ is a cutoff distance for defining the trapped state, and *F*_{d}(*R*_{d}) is the potential of mean force between the point defect and the electron, assuming both are dilute. To evaluate *P*_{trap}, we employ umbrella sampling^{68} and the Weighted Histogram Analysis Method (WHAM)^{69,70} by adding a bias potential $V(Rd)=0.5ksp(Rd\u2212Req)2$ on the distance between the defect and the centroid of the electron ring polymer with *k*_{sp} = 0.2 kcal mol^{−1} Å^{−2} and ${Req}={3,6,\u2026,195A\u030a}$. This allows us to determine the potential of mean force,

with *R*_{d} being the distance between the defect **x**_{d} and the centroid of the electron ring polymer $xec$.

The potential of mean force between an electron and the defect is shown in Fig. 5(a). For both CdS and MAPbI_{3}, the potential is monotonic. The binding free energy of the electron to the defect is much less than the bare potential energy of the pseudo-potential, reflecting the charge delocalization. This delocalization is evident in Fig. 5(b), which describes the charge distribution of the electron as a function of the distance between a point defect, *r*_{d−el} = |**x**_{d} − **x**_{e}|. For a free electron at a large distance, the electron has a spatial extent defined by the radius of gyration of the imaginary timeslices, which is related to the thermal wavelength *λ* at given temperature and mass, $Rg(\beta ,m)=\u210f\beta /2m$. For CdS and MAPbI_{3}, these free particle sizes are both nearly 20 $A\u030a$. Upon trapping to a defect, the electrons of both CdS and MAPbI_{3} become slightly more localized, with characteristic sizes of 15 and 17 $A\u030a$, respectively.

Given an equilibrium concentration of trapped electrons, the recombination rate is then evaluated by computing the likelihood of finding a hole in the vicinity of the electron and then measuring the conditional overlap in their densities, as described in Eq. (19). To evaluate both, we run MD simulations with a hole ring polymer and an electron ring polymer trapped into the point defect, where Eqs. (11) and (12) are used for two ring polymers. To sample the trajectory efficiently, we add the same harmonic potentials described above along the distance between two centroids of ring polymers and additional harmonic potentials with *k*_{sp} = 0.5 kcal mol^{−1} Å^{−2} and *R*_{eq} = 0.0 on the distance between the point defect and the centroid of the electron ring polymer to hold an electron near the defect. This potential is unweighted analogously with the WHAM to yield an unbiased distribution.

The resultant trap-assisted rate constants under different electron–phonon coupling strengths with *ρ*_{d} = 10^{18} cm^{−3} are summarized in Table III. We find that the interaction with phonons reduces the rate constant in both materials although values with finite coupling strength in CdS are not significantly distinct. This reduction results from the smaller likelihood of finding a hole in the vicinity of the electron, a manifestation of dynamical screening from the phonons. This is explored more directly for bimolecular recombination below.

CdS α | 0.0 | 0.53 | 1.5 | 3.0 |

k_{m} (μs^{−1}) | 0.93 | 0.37 | 0.36 | 0.44 |

k_{b}ρ_{e} (ns^{−1}) | 0.101 | 0.082 | 0.032 | 0.025 |

τ_{tot} (ns) | 9.79 | 12.2 | 31.3 | 39.6 |

MAPbI_{3} α | 0.0 | 2.87 | ⋯ | ⋯ |

k_{m} (μs^{−1}) | 0.115 | 0.058 | ⋯ | ⋯ |

k_{b}ρ_{e} (ns^{−1}) | 0.067 | 0.022 | ⋯ | ⋯ |

τ_{tot} (ns) | 15.0 | 45.6 | ⋯ | ⋯ |

CdS α | 0.0 | 0.53 | 1.5 | 3.0 |

k_{m} (μs^{−1}) | 0.93 | 0.37 | 0.36 | 0.44 |

k_{b}ρ_{e} (ns^{−1}) | 0.101 | 0.082 | 0.032 | 0.025 |

τ_{tot} (ns) | 9.79 | 12.2 | 31.3 | 39.6 |

MAPbI_{3} α | 0.0 | 2.87 | ⋯ | ⋯ |

k_{m} (μs^{−1}) | 0.115 | 0.058 | ⋯ | ⋯ |

k_{b}ρ_{e} (ns^{−1}) | 0.067 | 0.022 | ⋯ | ⋯ |

τ_{tot} (ns) | 15.0 | 45.6 | ⋯ | ⋯ |

### B. Bimolecular recombination rate

Bimolecular recombination is studied through the same method described above with MD simulations of electron and hole ring polymers. The bimolecular recombination rate, *k*_{b}*ρ*_{e} = *k*_{r}, requires us to evaluate the potential of mean force for localizing the two charges and their subsequent overlap density. The potential of mean force, *F*(*R*_{eh}), between the centroids of electron $xec$ and hole $xhc$ ring polymers,

is computable from umbrella sampling using the same procedure as that for *R*_{d}. The resulting function is shown in Fig. 6 for both CdS and MAPbI_{3} as functions of electron–phonon coupling, *α*. The monotonic free energies display systematic destabilization of the electron–hole pair with increasing *α*, consistent with the reduction in exciton binding energies. Similarly, the minimum is more shallow for MAPbI_{3} than for CdS.

In addition to the decreased likelihood of finding electron and hole pairs together, with increasing *α* the charge density distribution is broadened. Shown in Fig. 7 are the probability distributions of the distance between the electron and hole beads, *r*_{eh} = |**x**_{e} − **x**_{h}|, at each strength of interaction with phonons. In both materials, the stronger the phonon interaction is, the larger the average bead–bead distance becomes, implying that phonons make an effective electron–hole interaction weaker through screening. The comparison between results from path integral simulation with the probability distribution predicted from the Wannier–Mott exciton using a hydrogen model shown in purple lines in Fig. 7 implies the importance of capturing the fluctuations of quasiparticles at room temperature.

The ratio of path partition functions $Zc/Z0$ for the bimolecular rate constant can be computed and extrapolated to the large *n* limit in order to extract a converged overlap element, as shown in Appendix C. The calculated bimolecular recombination rates for both materials with typical charge density *ρ*_{e} = 10^{17} cm^{−3} are summarized in Table III as a function of *α*. The rates are found for both materials to decrease significantly over the range of electron–phonon coupling strength considered. MAPbI_{3} is found to have a longer charge carrier lifetime than the CdS, which is due to the enhanced screening of the former.

### C. Total recombination rate

Combining trap-assisted and bimolecular rates, the total recombination rates defined in Eq. (18) are summarized in Table III with a typical charge carrier density *ρ*_{e} = 10^{17} cm^{−3} and the defect density *ρ*_{d} = 10^{18} cm^{−3} for both materials. We find that the radiative recombination is predominantly determined by the bimolecular process. In both materials, electron–phonon coupling generally decreases the recombination rate, resulting in the increase in the lifetime of charge carriers. The value obtained for CdS is in very good agreement with that observed from photoluminescence lifetime measurements on large spherical nanocrystals, 13 ns, but underestimates the lifetime reported for MAPI_{3}, 70–100 ns.^{61} The latter disagreement can be attributed to the neglect of anharmonic effects accounted by the **k**-dependence of the correlation function in the optical mode, which can be considered previously^{39} and results in an effective electron–hole repulsion that has not been considered here.

## VII. CONCLUSIONS

In summary, we have shown how a path integral approach can be used to study the excitonic properties in the presence of dynamical phonons. We have presented ways to compute the renormalization of the binding energy and recombination rate and validated these results in limiting regimes. While we have considered a simple model for an exciton in a polar lattice as being coupled through a Fröhlich interaction with a single optical mode, the influence functional approach employed is general and can easily be extended to many modes, arbitrary linear coupling forms, and, in principle, parameterized through *ab initio* methods. Additionally, while we have considered phonon effects on a single exciton in this work, this can be applied to multiple excitons and further combined with confining potentials.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract No. DEAC02-05-CH11231 within the Physical Chemistry of Inorganic Nanostructures Program (Grant No. KC3103). This research used the resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility. Y.P. acknowledges the Kwanjeong Educational Foundation. D.T.L. acknowledges the Alfred P. Sloan Foundation.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Yoonjae Park**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (lead); Writing – original draft (equal); Writing – review & editing (equal). **David T. Limmer**: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Resources (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: GENERALIZED INFLUENCE FUNCTIONAL

In the main text, we have focused on a simplified model whereby the electron and hole are coupled to a single optical phonon through a linear Fröhlich coupling. For more complex lattice models where multiple dispersive modes are relevant, the influence functional formalism employed can be generalized. For many modes, the classical Hamiltonian entering into the path action can be written as

where **q**_{k,τ} is the classical displacement of the **k** mode at imaginary time *τ* and *ω*_{k} is its corresponding frequency. For a generalized linear coupling between the charge density and the lattice of the form

where *C*_{e/h,k} are generalized coupling coefficients, the phonons can still be integrated out. This yields an effective potential between the electron and hole of the form

where *χ*_{k}(*τ* − *τ*′) = ⟨**q**_{k}(0)**q**_{k}(*τ* − *τ*′)⟩ is the imaginary time correlation function of mode **q**_{k}. In Fourier space, the correlation function is given by

a sum of poles. In the classical limit, *βℏω*_{k} → 0, this influence functional returns a Coulomb potential screened by a wavevector dependent dielectric susceptibility.^{39}

### APPENDIX B: POLLMANN–BUTTNER THEORY

In order to compare our exciton binding energies to the Pollman–Buttner theory,^{28} we parameterize an effective potential of the form

where Δ*m* = *m*_{h} − *m*_{e}, 1/*ɛ** = 1/*ɛ*_{r} − 1/*ɛ*_{s}, $Re/h=\u210f/2me/h\omega $ with a static dielectric constant *ɛ*_{s}, which is taken to be 8.9 and 24.1 for CdS and MAPbI_{3}, respectively. We solve the time independent Schrödinger equation using this potential in the relative coordinate system for the electron and hole, **r**_{eh}. This takes the form

where *E*_{n} and *ϕ*_{n}(**r**_{eh}) are the associated eigenvalues and eigenvectors. Solving this equation for the ground state with zero angular momentum simplifies this to putting the equation on a real space grid on *r*_{eh}, which yields the exciton binding energies reported in the main text.

### APPENDIX C: OVERLAP DENSITY CONVERGENCE

The bimolecular rate constant *k*_{b} is determined by the ratio of partition functions $Zc/Z0$. The error from the path integral scales by the number of discretization factor as *n*^{−2}, so the convergence of the radiative rate can be extrapolated by fitting the ratios of partition functions at finite *n* as a function of *n*^{−2}. Extrapolations for CdS and MAPbI_{3} are plotted in Figs. 8(a) and 8(b) under different coupling strengths. Since the path integral formalism becomes exact in the limit of large number beads, the rate constant is defined as the value extrapolated to 1/*n*^{2} → 0 limit.

## REFERENCES

_{2}

_{2}and WSe

_{2}

_{3}NH

_{3}PbI

_{3}: Effect of spin-orbit interaction, semicore electrons, and self-consistency

_{3}NH

_{3}PbI

_{3}from multiphonon Fröhlich coupling