We introduce a simple and effective method to decompose the energy dissipation in the dynamics of open quantum systems into contributions due to individual bath components. The method is based on a vibronic extension of the Förster resonance energy transfer theory that enables quantifying the energy dissipated by specific bath degrees of freedom. Its accuracy is determined by benchmarking against mixed quantum–classical simulations that reveal that the method provides a semi-quantitative frequency-dependent decomposition of the overall dissipation. The utility of the method is illustrated by using a model donor–acceptor pair interacting to a thermal harmonic bath with different coupling strengths. The method can be used to identify the key features of a bath that leads to energy dissipation as required to develop a deep understanding of the dynamics of open quantum systems and to engineer environments with desired dissipative features.

## I. INTRODUCTION

Open quantum system dynamics refers to the evolution of a quantum subsystem interacting with an environment.^{1,2} This setup describes a vast range of chemical processes that occur in the condensed phase, in which the behavior of a few relevant molecular degrees of freedom (DOFs) of interest are influenced by a macroscopic thermal environment. An important class of environments is those that can be described as a collection of harmonic oscillators,^{3} such as photonic, vibrational, and solvation environments.^{4,5} For harmonic baths, the key quantity to describe the subsystem–bath interaction is the bath spectral density (BSD), which characterizes the profile of the subsystem–bath coupling strength for different bath frequencies. The shape of the BSD shows great diversity, as can be observed from various analytic models^{6–9} and realistic models obtained from atomistic simulations.^{10–12} The effect of the bath on the subsystem is intricately related to the BSD.

To better understand the dynamics of an open quantum system, it is highly desirable to be able to resolve the effects arising from individual bath components and, in particular, their contributions to the overall energy dissipation. In this manner, one can gain knowledge about the relaxation channels of the system and related time scales, as needed to develop a clear picture of the major factors governing the dynamics. Such investigations are essential to address the grand challenge of controlling the open quantum system dynamics through reservoir engineering.^{13–16}

However, direct calculation of dissipation routes is challenging as it requires dynamical information about the environment. In the context of electron–nuclear problems, where vibrations provide a thermal environment for the electrons, a popular strategy is to explicitly include the vibrational modes in the subsystem and follow the dynamics of the resultant electronic-vibrational (vibronic) Hilbert space.^{17–19} However, this approach can only treat a small number of vibrational modes at once, due to the exponential increase in computational cost with subsystem dimensionality. Efficient wavefunction propagation schemes such as multi-configurational time-dependent Hartree (MCTDH)^{20–22} and time-dependent density matrix renormalization group (TD-DMRG)^{23–25} can treat relatively large dimensionalities in a formally exact manner, although simulating nonzero temperatures usually requires truncation of the vibrational subspace to reduce the computational cost,^{26,27} leading to modest accuracy. This is particularly troublesome for treating the low-frequency part of the BSD, for which many levels are significantly occupied at room temperature.

Mixed quantum–classical simulation methods^{28–31} can be used to directly calculate dissipation from the positions and momenta of classical bath trajectories, by splitting the total energy of the system into contributions from the subsystem and individual bath modes. The advantages of this strategy are that the computational cost only increases linearly with the number of vibrational modes and that truncation of the vibrational Hilbert space is not needed. However, to obtain statistically meaningful results, one still needs to propagate a large number of trajectories with a small integration time step.^{19}

Therefore, there is a primary need for new strategies that can efficiently resolve the overall dissipation into contributions from individual bath components.

In this work, we develop an efficient dissipation rate equation for molecular systems coupled to harmonic environments that conveniently decomposes the overall dissipation into contributions from individual vibrational modes with a significantly low computational cost. We accomplish this by extending Förster resonance energy transfer (FRET) theory^{32} into vibronic subsystems and performing an analytical summation over the vibrational Hilbert space. The accuracy of our theory is assessed by benchmarking against the dissipation calculated by a mixed quantum–classical simulation method, specifically, the Poisson bracket mapping equation with a non-Hamiltonian modification (PBME-nH). An extensive benchmark against a numerically exact simulation method has shown that PBME-nH provides qualitatively correct results for moderate subsystem–bath couplings.^{19}

This paper is organized as follows: In Sec. II, we introduce the Hamiltonian model and theory underlying our work and derive expressions for calculating dissipation rates. Section III illustrates the basic features of the theory and provides benchmark computations. Section IV uses these expressions to analyze the physical origin of frequency dependence of dissipation in a model donor–acceptor system each coupled to a harmonic bath. Section V summarizes the main findings, and discusses the computational cost and the future prospects for the proposed strategy.

## II. THEORY

The developed theory is based on FRET, which describes the transfer of electronic excitation among weakly interacting chromophore molecules. Below, we review FRET and its basic assumptions and use this theoretical framework to capture the dissipation by individual vibrational modes. This is done by explicitly including the vibrational mode into the subsystem and constructing vibronic FRET equations. As discussed, we analytically combine microscopic energy transfer events between all pairs of vibronic states, to isolate convenient expressions for calculating dissipation.

### A. Hamiltonian

Consider *N* two-level molecules composed of their ground and first excited electronic states, each of them interacting with a local harmonic environment. We focus on the transfer of electronic excitation in the single-excitation manifold {|*A*⟩}, where |*A*⟩ describes a state for which only molecule *A* is in its excited state, while the rest are in their ground state. We write the total electron–nuclear Hamiltonian of the system as

Here,

is the bare electronic Hamiltonian where *E*_{A0} is the zero-phonon (0–0) transition energy between the ground and excited states of molecule *A*. The nuclear DOFs are described by a collection of quantum harmonic oscillators,

where the origin of the coordinates defines the minimum of the ground-state potential energy surface (PES). Here, the operators $p^Ak$ and $x^Ak$ are the mass-weighted momentum and coordinate operators of the *k*th vibrational mode coupled to molecule *A*. In turn, the vibronic coupling Hamiltonian

modifies the energy gap between the ground and excited state PESs along the vibrational coordinates {*x*_{Ak}}. The strength of the vibronic coupling along *x*_{Ak} is proportional to the displacement *d*_{Ak} between the ground and excited state PESs. The reorganization Hamiltonian,

where $\lambda A=\u2211k(\omega Ak2dAk2/2)$ is the reorganization energy of molecule *A*, accounts for the difference between the vertical and zero-phonon transition energies induced by the vibronic coupling [Eq. (4)]. Finally,

is the Coulombic coupling Hamiltonian that describes the inter-molecular coupling. Here, *V*_{AB} = *V*_{BA} is the electronic coupling between molecules *A* and *B*, which is assumed to be independent on the bath coordinates, and H.c. denotes the Hermitian conjugate.

For future reference, we introduce the BSD,

where

is the dimensionless Huang–Rhys factor that quantifies the strength of vibronic coupling. We also introduce the single-mode reorganization energy *λ*_{Ak},

which satisfies ∑_{k}λ_{Ak} = λ_{A}. Combining Eqs. (7) and (9) gives the relation between the BSD and the reorganization energy,

### B. Förster resonance energy transfer

The total Hamiltonian [Eq. (1)] is treated in the framework of perturbation theory by writing $\u0124=\u01240+\u01241$, where

FRET theory supposes that: (1) The inter-molecular coupling *V*_{AB} is much smaller than either the excitation energy gap |*E*_{A0} − *E*_{B0}| or the standard deviation of the vibronic coupling $\u27e8\u0124e-n2\u27e91/2$ so that $\u01241$ can be treated perturbatively. (2) Markovian bath, so the nuclear DOFs instantly relax to the equilibrium. (3) A local nuclear thermal equilibrium defined by the Hamiltonian $\u2009A|\u01240|A\u2009$. That is, the vibrational modes coupled to molecule *A* relax to the excited state equilibrium, while the other modes relax to the ground-state equilibrium. This reflects that FRET describes each molecule as a separate entity, which is valid for small *V*_{AB}.

By applying Nakajima–Mori–Zwanzig projection operator techniques under these assumptions and keeping terms up to second-order in *V*_{AB},^{33} one arrives at

where *P*_{A}(*t*) is the population of the state |*A*⟩ at time *t* and *K*_{BA} is the rate constant for the population transfer from |*A*⟩ to |*B*⟩. Note that Eq. (12) focuses on the populations and neglects any effects of coherence between chromophores, which is valid when the effect of $\u01241$ is small. The FRET rate constants are evaluated according to

where $FA(t)$ and $AB(t)$ are defined by

Here, *g*_{A}(*t*) is the line broadening function,

where *β* = 1/*k*_{B}*T* and *k*_{B} is Boltzmann’s constant. The physical meaning of Eq. (13) becomes more apparent if we express it with the Fourier-transformed quantities,

which are actually the spectra for the fluorescence of molecule *A* and the absorption of molecule *B*, whose areas are normalized to unity. Due to the relations $F(t)=F*(\u2212t)$ and $A(t)=A*(\u2212t)$, $F(\omega )$ and $A(\omega )$ are real quantities. Rewriting Eq. (13) by using Eq. (16) leads to

Equation (17) shows that the FRET rate constant *K*_{AB} is proportional to the overlap integral between the donor fluorescence line shape $FA(\omega )$ and the acceptor absorption line shape $AB(\omega )$.

### C. vFRET: Extension of FRET to vibronic dynamics

To capture the dissipation into individual vibrational modes, it is necessary to simulate correlated quantum dynamics between electronic and vibrational DOFs. To do so, we take advantage of the structure of FRET and extend its applicability to capture vibronic dynamics. We call the resulting theory vibronic FRET (vFRET).

We start by decomposing the full Hamiltonian [Eq. (1)] into

where, in the right-hand side of Eq. (18), $\u0124e+\u0124Coul$ is the electronic contribution, while $\u2211A,k\u0124Ak$ is the vibrational contribution. The Hamiltonian component $\u0124Ak$ for an individual vibrational mode of index *Ak* is

where

are the ground and excited state PESs of the vibrational mode *Ak*.

We investigate the energy flow into a given vibrational mode by re-assigning it from the bath to the subsystem and consider the dynamics in the vibronic subsystem formed by the electronic DOFs and the re-assigned vibrational mode (RAVM). We assume that the RAVM couples to molecule *A* and denote its index as $v$. As a complete basis for the vibronic subsystem, we pair |*A*⟩ with the excited state vibrational eigenstates {|*I*_{e}⟩} of $\u0125Ak,e$ and all other electronic states {|*B*⟩} with the ground-state vibrational eigenstates {|*J*_{g}⟩} of $\u0125Ak,g$. According to the second and third assumptions of FRET, the RAVM is always in the local thermal equilibrium of the nuclear PESs $\u2009A|\u0124Ak|A\u2009=\u0125Ak,e$ or $\u2009B|\u0124Ak|B\u2009=\u0125Ak,g$. Therefore, an advantage of using the {|*A*⟩|*I*_{e}⟩} and {|*B*⟩|*J*_{g}⟩} bases is that there is no coherence so that the vibronic population transfer can be described in the framework of FRET. From now on, we denote {|*A*, *I*⟩ ≡ |*A*⟩|*I*_{e}⟩} and {|*B*, *J*⟩ ≡ |*B*⟩|*J*_{g}⟩}. Figure 1 schematically illustrates the population transfer between the subsystem basis states of the electronic and vibronic descriptions. Accordingly, we re-write the total Hamiltonian [Eq. (1)],

where each Hamiltonian component is the vibronic generalization of the corresponding component without prime in Eq. (1). To construct them, we take into account that the RAVM is now re-assigned to the subsystem and, hence, not included in the bath any more. That is, the electronic Hamiltonian is now

where we have added the vibrational eigenenergies of the RAVM to *E*_{A0} and *E*_{B0}. The quantity

defines the ground-state PES without the RAVM. In turn,

is the vibronic coupling Hamiltonian without the RAVM. Here, we have taken into account that all vibrational modes are independent of each other, and therefore, $\u0124e-n\u2032$ still adopts a direct product form as in Eq. (4). The reorganization Hamiltonian is now

where we have replaced *λ*_{A} in Eq. (5) by *λ*_{A−} ≡ *λ*_{A} − *λ*_{Av}, to take into account that the contribution of the RAVM needs to be removed. The coupling between the vibronic states is

The inter-chromophore couplings in Eq. (26) are now scaled by the overlap between the vibrational eigenstates of the RAVM. The equivalence between Eqs. (1) and (21) can be verified by evaluating the Hamiltonian matrix elements in the {|*A*, *I*⟩, |*B*, *J*⟩} basis by using both representations.

To calculate the dissipation by RAVM, we need to determine its energy. For this, we construct expressions that govern the vibronic population transfer by applying FRET to Eqs. (21)–(26). This can be achieved by modifying the electronic expressions [Eqs. (12)–(14)] by using the vibronic quantities in Eqs. (21)–(26). We start from the rate equations for vibronic populations, which are generalizations of Eq. (12) to the vibronic description,

Note that we are separately treating {|*A*, *I*⟩} and {|*B*, *J*⟩} in Eq. (27) to distinguish the ground and excited vibrational quantum states for the RAVM. We also took into account that *K*_{AI,AI′} = *K*_{BJ,BJ′} = 0, as there is no coupling between such vibronic state pairs according to Eq. (26). This reflects that the bath is always in the thermal equilibrium, and therefore, there is no direct population transfer within the vibrational manifold of a single molecule. We classify the nonzero vibronic rate constants in Eq. (27) into three types. The rate constants {*K*_{BJ,AI}} and {*K*_{AI,BJ}} govern the outward and inward population transfer with respect to molecule *A* that is coupled to the RAVM. These two types of vibronic population transfers change the PES of the RAVM and induce vibrational relaxation of the RAVM. On the other hand, {*K*_{BJ,CJ′}} keeps the RAVM at the ground-state PES and, hence, does not involve any vibrational relaxation of RAVM.

The expressions for vibronic rate constants in Eq. (27),

can be constructed from the electronic rate constant [Eq. (13)] by replacing the electronic coupling *V*_{AB} by the scaled couplings in Eq. (26) and the fluorescence and absorption line shapes $F(t)$ and $A(t)$ by their vibronic analogs determined by

Equation (29) is deduced by generalizing the electronic line shapes [Eq. (14)] by (i) converting the electronic state energies to the vibronic state energies [Eq. (22)] and (ii) removing the contribution of the RAVM from the bath-related quantities in Eqs. (29a) and (29b). As a result, *λ*_{A} and *g*_{A}(*t*) in Eq. (14) are replaced by *λ*_{A−} and *g*_{A−}(*t*) defined by

where *g*_{Av}(*t*) is the line broadening function of the RAVM,

evaluated by plugging the spectral density of the RAVM $JAv(\omega )=\u210f\omega Av2sAv\delta (\omega \u2212\omega Av)$ [Eq. (7)] into Eq. (15). Meanwhile, *λ*_{B} and *g*_{B}(*t*) in Eqs. (29c) and (29d) remain unaffected because the exclusion of the RAVM from the bath only changes the bath structure of molecule *A*, as can be seen from Eqs. (23) and (24).

Equations (27)–(29) form vFRET theory. We note that these expressions can be rigorously obtained by applying the projection operator technique to Eqs. (21)–(26), and our approach is just a convenient shortcut for the derivation. The physical consistency of vFRET is demonstrated in Appendix A by showing that the net electronic population transfer in FRET and vFRET coincide. Indeed, as the electronic and vibronic descriptions in Fig. 1 are just two different representations of the same Hamiltonian, their physical behaviors must be the same. However, do note that to make the FRET and vFRET identical, it is necessary to impose the Markov approximation on the RAVM. That is, to insist that the vibronic populations always satisfy

where $w$_{I} and $w$_{J} are the Boltzmann weights for the vibrational eigenstates of the RAVM,

### D. Dissipation rate equation

With the rate equations for the vibronic state populations [Eqs. (27)–(29)], it is now possible to derive the dissipation rate equation for the RAVM. We define the energy contained in the RAVM as the weighted average of the vibrational eigenenergies,

as there are no coherences between the vibronic states {|*A*, *I*⟩} and {|*B*, *J*⟩} in vFRET due to the second and third assumptions of FRET. Meanwhile, the populations of the vibronic states must satisfy Eqs. (32) and (33) at every instance. Plugging these equations into Eq. (34) gives

which is constant over time and, hence, does not capture the dissipation. This is because the Markov approximation [Eq. (32)] insists that any extra vibrational energy gained by dissipation must be quenched instantly. Nevertheless, we can still determine the amount of dissipated energy by differentiating Eq. (34) with respect to time and expressing the time derivatives of vibronic populations by using Eq. (27). That is,

Equation (36) corresponds to the vibrational energy gained by the RAVM per unit time and, therefore, is equivalent to the dissipation rate. This gained energy is immediately lost by imposing Eq. (32), as discussed in Sec. II C.

The amount of dissipation by the RAVM can now be calculated by integrating Eq. (36) over time. However, this procedure involves summations over infinite number of vibronic states and, hence, is not very efficient. To develop much more practical expressions, we express the vibronic populations and rate constants in Eq. (36) by using Eqs. (28) and (32). This gives

To simplify Eq. (37), we need to perform analytical summations over vibrational quantum numbers *I* and *J*. As we are summing line shapes weighted by the thermal population and vibrational overlap, we are motivated to consider each of the sums as a spectral line shape dressed by a thermal phonon. In the formulation presented in Appendix B, we show that $FA(\omega )$ and $AA(\omega )$ are related to the line shapes without the RAVM, $FA\u2212(\omega )$ and $AA\u2212(\omega )$, through the following convolution relations:^{34}

where $FA\u2212(\omega )$ and $AA\u2212(\omega )$ are constructed from Eqs. (14) and (16) by making substitutions *λ*_{A} → *λ*_{A−} and *g*_{A}(*t*) → *g*_{A−}(*t*),

and $LAv(\omega )$ is the phonon sideband of the RAVM,

Combining Eqs. 29(a) and 29(b) with Eqs. (38)–(40) gives us the relation between the electronic and vibronic line shapes,

To use Eq. (41), however, we first need to remove *I* − *J* from the prefactors of Eq. (37). This can be accomplished by using time derivatives of $FAI*(t)ABJ(t)$ and $FBJ*(t)AAI(t)$, as both of them have *I* − *J* in the exponential. As a result, we can derive the identities

We plug Eq. (42) in Eq. (37) and reallocate the dependence on *I* and *J* solely to the vibronic line shapes of molecule *A* by using

Then, to the resulting expression, we can now use Eq. (41) to sum $FA(I\u2212J)(t)$ and $AA(I\u2212J)(t)$. Furthermore, explicitly carrying out the time derivative in the prefactor by using Eq. (14) makes most terms cancel out and yields

Here, we have used the relations *λ*_{A} − *λ*_{A−} = *λ*_{Av} and *g*_{A}(*t*) − *g*_{A−}(*t*) = *g*_{Av}(*t*) [Eq. (30)]. The prefactor on the right-hand side can be further simplified by using Eqs. (9) and (31),

where $IBA(\omega )$ and $IAB(\omega )$ are

with the dissipation rate constants $KAvBA$ and $KAvAB$ given by

To generalize Eqs. (48) and (49) to treat continuous BSDs, based on Eq. (10), we make substitutions *ω*_{Av} → *ω* and $\lambda Av\u2192JA(\omega )\omega d\omega $ in Eq. (49) and define “dissipative spectral densities” $JABA(\omega )$ and $JAAB(\omega )$ as follows:

The dissipative spectral densities quantify the ability of the bath to induce dissipation and are determined as products of the electronic coupling *V*_{AB} = *V*_{BA}, [Eq. (47)], and the bath reorganization energy per unit frequency *J*_{A}(*ω*)/*ω*. Although the dissipative spectral densities are local in the frequency domain, they actually encode the full structure of the environment as $I(\omega )$ is evaluated from the spectral line shapes ${F(t)}$ and ${A(t)}$. By substituting $KAvBA$ and $KAvAB$ in Eq. (49) with $JABA(\omega )\u2002d\omega $ and $JABA(\omega )\u2002d\omega $, we can express the rate of dissipation in the frequency window [*ω*, *ω*+ *dω*] at time *t* as

Equation (51) is the main result of this work, which enables us to calculate the time-dependent dissipation by specific regions of BSD. Note that all the elements in Eq. (51) are expressed purely in terms of electronic quantities, ${FA(t)}$, ${AA(t)}$, and {*P*_{A}(*t*)}, so the calculations can be done without directly accessing the vibrational manifold. Because Eq. (51) contains {*P*_{A}(*t*)}, calculating dissipation requires knowledge of the electronic dynamics. In practice, both the electronic rate constants {*K*_{AB}} [Eq. (13)] and dissipative spectral densities ${JABA(\omega ),\u2009JAAB(\omega )}$ [Eq. (50)] are evaluated at the start of the simulation. Then, {*P*_{A}(*t*)} and ${DA(\omega ,t)}$ are propagated simultaneously by using Eqs. (12) and (51).

## III. ACCURACY AND FEATURES OF vFRET

To assess the utility of Eq. (51) to investigate dissipation, we benchmark it against PBME-nH,^{35} a mixed quantum–classical simulation method recently extended to calculate dissipation by individual vibrational modes.^{19} PBME-nH is based on real-time evolution of the nuclear phase space density and, thus, allows us to test the assumptions behind vFRET theory such as the Markovian approximation and the choice of thermal equilibrium for the nuclear DOFs. Furthermore, PBME-nH provides consistent results regardless of whether RAVM is included either in the subsystem or bath (Fig. 1),^{36} which justifies our calculation of dissipation without an explicit treatment of the RAVM. Despite these advantages, PBME-nH is still an approximate method, so we further check its accuracy by comparing the electronic dynamics to the results from the numerically exact hierarchical equations of motion (HEOM) method.^{37,38}

We note that the formulation for PBME-nH dissipation in Ref. 19 did not include the contribution from the mode reorganization energy $\lambda Ak=\omega Ak2dAk2/2$ [see the right-hand side of Eq. (20b)]. This makes the vibrational eigenenergies for the ground and excited state PESs different and introduces spurious effects in the dissipation as the electronic populations change. Therefore, the PBME-nH dissipation requires a slight modification as detailed in the supplementary material.

Section III A describes the simulation procedure and the model donor–acceptor complex employed to illustrate the theory. Section III B presents the electronic dynamics of the model and benchmarks FRET and PBME-nH against HEOM simulations. Then, in Sec. III C, we compare the dissipation calculated by vFRET and PBME-nH and numerically demonstrate that our vFRET theory satisfies the detailed balance condition.

### A. Simulation procedure

As a minimal model that allows simple interpretations, we consider a donor–acceptor complex whose excitation energy difference is *E*_{D0} − *E*_{A0} = 400 cm^{−1} and the electronic couplings are *V*_{DA} = *V*_{AD} = 50 cm^{−1}. As |*E*_{D0} − *E*_{A0}| ≫ *V*_{DA}, the system satisfies the first assumption of FRET that requires relatively weak electronic coupling (Sec. II B). Each molecule is coupled to a Drude–Lorentz spectral density,

where *λ*_{ph} is the reorganization energy and *ω*_{c} is the cutoff frequency. In this work, we vary *λ*_{ph} while keeping *ω*_{c} = 0.02 rad fs^{−1} (*ℏω*_{c} ∼ 106 cm^{−1}) constant. The temperature of the bath is *T* = 300 K. The electronic excitation is assumed to be initially localized in the donor molecule.

For FRET and vFRET simulations, the time integrals needed for calculating the electronic rate constants [Eq. (13)] and ${I(\omega )}$ [Eq. (47)] were evaluated by using the midpoint method with a grid size of 0.5 fs and an upper limit of integration of *t* = 50 ps. Electronic populations [Eq. (12)] and dissipated energies [Eq. (51)] were propagated by using the fourth-order Runge–Kutta method with a 0.5 fs time step. For PBME-nH simulations, the initial conditions for the trajectories were sampled according to Ref. 19 and propagated using the integration scheme described in Ref. 39. The propagation time steps for PBME-nH were 0.5 fs for simulating electronic dynamics and 0.05 fs for calculating the dissipation. The number of PBME-nH trajectories was 10^{5} for simulating the electronic dynamics, 10^{6} for the dissipation calculation with *λ*_{ph} = 10 cm^{−1}, and 10^{7} for all other dissipation calculations. The HEOM simulations were based on an efficient formulation recently developed by Ikeda and Scholes,^{40} with the three-term Padé series expansion of the Bose–Einstein distribution function^{41} and the hierarchy depth of 7. The integration time step was 0.1 fs.

### B. Electronic population dynamics

Before analyzing the dissipation, we look into the electronic dynamics and assess the accuracy of FRET and PBME-nH by comparing the results to numerically exact HEOM simulations. We note that HEOM cannot be used toward decomposing the dissipation, which is the reason why we only use HEOM to benchmark the accuracy of the electronic dynamics.

As a simulation method based on classical nuclei, PBME-nH is well-known to suffer from zero-point energy (ZPE) leak of the phase space density, especially for high-frequency vibrational modes whose ZPE is much larger than the thermal energy (*ℏω* ≫ 2 kT).^{19} To reduce the influence of ZPE leak on the dissipation, we have set *J*(*ω*) beyond *ℏω* = 800 cm^{−1} identically to zero for all PBME-nH simulations. The remaining part of the BSD was discretized into 2000 harmonic oscillator modes by using the discretization scheme in Ref. 42, which recovered 91.6% of the original reorganization energy of the pristine spectral density. For consistency, we also applied the same BSD cutoff to FRET simulations and calculated the reorganization energy of the modified BSD by using Eq. (10) instead of *λ*_{ph}. We have confirmed that the cutoff had only negligible effect on the electronic dynamics simulated by FRET and PBME-nH.

Figure 2(a) displays the time-dependent electronic populations of the donor and acceptor, calculated by FRET and PBME-nH with *λ*_{ph} = 50 cm^{−1}. Due to the Markovian approximation behind the FRET theory, the donor and acceptor populations can be precisely modeled by exponential profiles,

Here, *P*_{D}(∞) and *P*_{A}(∞) are the equilibrium donor and acceptor populations that we take from the end of the simulation, and *γ* is the effective population transfer rate. It is also known that the equilibrium populations also satisfy the detailed balance condition,^{33}

and, hence, is not affected by *λ*_{ph}.

In Fig. 3, we summarize the electronic dynamics simulated by FRET, PBME-nH, and HEOM for various *λ*_{ph}. Figure 3(a) presents the equilibrium acceptor population *P*_{A}(∞), which shows that *P*_{A}(∞) from FRET simulations show excellent agreement with HEOM results and remain constant for all values of *λ*_{ph}, as predicted by Eq. (54). By contrast, *P*_{A}(∞) calculated by PBME-nH matches HEOM results for small values of *λ*_{ph}, but increases as *λ*_{ph} becomes larger and eventually exceeds unity for *λ*_{ph} > 400 cm^{−1}. It is known that such unphysical behavior arises from the ZPE leak among electronic DOFs, which is a common problem for most simulation methods that employ the semi-classical description of electronic DOFs^{39,43,44} as in PBME-nH.

Figure 3(b) plots the effective transfer rate *γ* extracted from the time-dependent electronic populations. Because the relaxation time of the bath (2*π*/*ω*_{c} = 50 fs) is much quicker than the picosecond time scale of the electronic dynamics, the electronic dynamics is in the Markovian regime. As a result, almost perfect exponential fits were obtained even for HEOM and PBME-nH, which are not based on the Markovian approximation. One can observe that FRET slightly overestimates *γ* for all values of *λ*_{ph}, while PBME-nH underestimates *γ* for *λ*_{ph} > 50 cm^{−1}. Such a deterioration of PBME-nH for large bath reorganization energies has been already observed in a few studies.^{19,35} Nevertheless, both FRET and PBME-nH correctly predict that *γ* increases until *λ*_{ph} = 150 cm^{−1} − 200 cm^{−1} and decreases afterward, which is a typical manifestation of the so-called environment-assisted quantum transport.^{45,46} Such a behavior can be readily explained based on spectral line shapes,^{47} which directly affects the electronic rate constants [Eq. (13)].

Despite the observed deviation of PBME-nH from exact HEOM results, the method still provides a qualitatively correct picture about how the dissipated energy is distributed among individual vibrational modes that can be employed to assess the utility of vFRET.

### C. Dissipation dynamics

#### 1. Time profile

We now examine the dissipation. To do so, we define

which quantifies the dissipation by *J*_{B}(*ω*) per unit bath frequency accumulated until time *t* by either donor (*B* = *D*) or acceptor (*B* = *A*). By inserting the time profile of the electronic populations [Eq. (53)] in the rate of dissipation [Eq. (51)], arranging the expression for $DB(\omega ,t)$ by using the fact that the net dissipation must vanish at equilibrium [$DB(\omega ,\u221e)=0$], and finally integrating the resulting expression by using Eq. (55), consequently,

For our model with *J*_{D}(*ω*) = *J*_{A}(*ω*), Eq. (56) predicts that the dissipation into the donor and acceptor will be identical as $JDAD(\omega )=JAAD(\omega )$ [Eq. (50)].

Figure 2(b) shows the total energy dissipated into the donor bath modes in the [40, 50] cm^{−1} frequency window calculated by integrating $ED(\omega ,t)$ [Eq. (55)] in this range. The dissipated energy for vFRET exhibits an exponential profile and reaches equilibrium with the electronic dynamics [Fig. 2(a)], as expected by Eq. (56). However, we also observe that the net dissipation for PBME-nH does not vanish, even after the electronic dynamics reaches equilibrium. This artifact arises from the ZPE leakage from the high-frequency modes to the low-frequency modes (not shown). To remove the artificial drifts from the PBME-nH dissipation, we correct it according to

#### 2. Frequency-dependent decomposition

In this section, we decompose the total dissipation in the frequency domain, $E\u0303D(\omega ,\u221e)+E\u0303A(\omega ,\u221e)$. Figure 4 shows $E\u0303D(\omega ,\u221e)+E\u0303A(\omega ,\u221e)$ for four different values of *λ*_{ph}, calculated by vFRET and PBME-nH. Consider the vFRET results in Fig. 4(a) first. When *λ*_{ph} = 10 cm^{−1}, the dissipation is concentrated around *ω* = 380 cm^{−1}. By contrast, for larger *λ*_{ph}, the dissipation is concentrated in the low-frequency region. A discussion of the physical origin of this behavior is included in Sec. IV.

The same trend is also observed for the PBME-nH results [Fig. 4(b)], confirming the validity of vFRET theory. A notable difference between vFRET and PBME-nH is the low-frequency region, where PBME-nH results show a sharp dip, while vFRET results do not. This implies that the low-frequency modes actually do not participate in the dissipation as actively as expected by vFRET. This can be explained by considering that the low-frequency modes are dynamically frozen and, therefore, take longer to undergo relaxation. This interpretation is consistent with a previous study,^{48} which showed a dramatic increase in the accuracy of the Markovian quantum master equation when the low-frequency modes were separated from the bath and treated as static noise. Such a quasi-static nature of the low-frequency modes causes violation of the second assumption of FRET about Markovianity (Sec. II B), which extends to a wider frequency range as the electronic dynamics becomes faster. Indeed, Fig. 4(b) shows that the low-frequency dip becomes broader as *λ*_{ph} increases up to 200 cm^{−1}, for which the population transfer is the fastest among the four values of *λ*_{ph} [Fig. 3(b)]. We also expect the dip will become narrower for *λ*_{ph} = 500 cm^{−1} compared to *λ*_{ph} = 200 cm^{−1}, although confirming this from Fig. 4(b) is difficult because of the statistical error of PBME-nH. Our observation demonstrates that even when the electronic dynamics is Markovian, the vibronic dynamics that leads to dissipation may not be, especially if the movement of the nuclear DOFs is slow.

To further test our interpretation about the discrepancy in the low-frequency part, we decreased the electronic coupling *V*_{DA}, thus slowing down the electronic dynamics (Fig. 5). As expected, while decreasing *V*_{DA} had essentially no effect on the vFRET result, it made the low-frequency dip in PBME-nH results [Fig. 5(a)] nearly disappear in Fig. 5(b), leading to almost perfect agreement with vFRET.

#### 3. Detailed balance

We now numerically verify that the dissipation calculated by the vFRET theory satisfies detailed balance. By doing so, we probe the consistency between the electronic and dissipation dynamics.

We first identify the conditions met by the electronic populations and dissipation rates at equilibrium. At equilibrium, both the net population transfer and dissipation vanish so that $P\u0307D(\u221e)=P\u0307A(\u221e)=0$ and $DD(\omega ,\u221e)=DA(\omega ,\u221e)=0$. By plugging these conditions into Eq. (51), we get

Equation (58) shows that $JBAD(\omega )$ and $JBDA(\omega )$ must have opposite signs, as the electronic populations are positive. Note that we do not distinguish the dissipative spectral densities of the donor and acceptor because $JDAD(\omega )=JAAD(\omega )$ and $JDDA(\omega )=JADA(\omega )$ (Sec. III C 2). Because the equilibrium electronic populations satisfy detailed balance,^{33} all the ratios in Eq. (58) must be equal to exp[*β*(*E*_{D0} − *E*_{A0})] = 6.810 143 83… for our simulation conditions.

Figure 6 displays $\u2212JBAD(\omega )/JBDA(\omega )$ calculated for four different values of *λ*_{ph}. We observe that $\u2212JBAD(\omega )/JBDA(\omega )$ is practically equal to exp[*β*(*E*_{D0} − *E*_{A0})] for all values of *ω* and *λ*_{ph}, which confirms that the dissipation calculated by vFRET satisfies the detailed balance. An analytical proof of the detailed balance remains outstanding.

## IV. ANALYSIS OF THE DISSIPATION RATE

In Sec. III C, we have shown that vFRET can accurately capture the frequency dependence of the dissipation. We now illustrate how to use vFRET to interpret the effectiveness of different components of the bath to exert dissipation. For this purpose, we scrutinize the dissipative spectral densities $JBAD(\omega )$ and $JBDA(\omega )$ [Eq. (50)] that are directly proportional to $EB(\omega ,t)$ [Eq. (56)]. For definitiveness, we analyze cases with small (*λ*_{ph} = 10 cm^{−1}) and large (*λ*_{ph} = 500 cm^{−1}) reorganization energies and identify the main factors that determine the dissipation rate constants. The results for intermediate reorganization energies, *λ*_{ph} = 50 cm^{−1} and *λ*_{ph} = 200 cm^{−1}, can be found in Fig. S1 in the supplementary material.

According to Eq. (50), the vFRET theory clearly highlights the three main factors that determine the dissipation rate among individual modes: (i) The electronic coupling *V*_{DA} as it directly affects the speed of FRET [Eq. (13)]. (ii) The “dissipative potential” $I(\omega )$ [Eq. (47)], which quantifies the ability of the bath mode of a certain frequency to participate in the dissipation. (iii) The reorganization energy density *J*_{B}(*ω*)/*ω* quantifying the strength of vibronic coupling in the frequency domain. The overall contribution of a bath component in the dissipation is determined by the product of these three quantities. We now explain the trend observed in Fig. 4(a) based on these three factors. Because *V*_{DA} is fixed in our model, we consider how $JB(\omega )$ is affected by $I(\omega )$ and *J*_{B}(*ω*)/*ω*.

Figure 7(a) displays $JBAD(\omega )$ and $JBDA(\omega )$. The dissipative spectral density for *λ*_{ph} = 10 cm^{−1} is about 20 times smaller than those for *λ* = 500 cm^{−1}, reflecting the slower relaxation due to the smaller spectral overlap between donor and acceptor line shapes. We focus on $JBAD(\omega )$ as $JBAD(\omega )$ and $JBDA(\omega )$ only differ by a constant prefactor [Fig. 6(b) and Eq. (58)]. The frequency dependence of $JBAD(\omega )$ is very different for the two reorganization energies. When *λ*_{ph} = 10 cm^{−1}, the major maximum is located at *ω* = 380 cm^{−1}, while another maximum faintly occurs at *ω* = 0 cm^{−1}. By contrast, for *λ*_{ph} = 500 cm^{−1}, *ω* = 0 cm^{−1} is the only maximum. The origin of such a stark contrast can be explained by separately examining $IAD(\omega )$ and *J*_{B}(*ω*)/*ω*.

The frequency dependence of the dissipative potentials $IAD(\omega )$ and $IDA(\omega )$ is plotted in Fig. 7(b). When *λ*_{ph} = 10 cm^{−1}, $IAD(\omega )$ peaks at *ω* ∼ 400 cm^{−1}, which coincides with the donor–acceptor energy gap *E*_{D0} − *E*_{A0}. As *V*_{DA} ≪ |*E*_{D0} − *E*_{A0}| in our system, *E*_{D0} − *E*_{A0} = 400 cm^{−1} is close to the eigenenergy gap of $\u0124e+\u0124Coul$ where vibronic resonance occurs, which is $(ED0\u2212EA0)2+4|VDA|2\u223c412\u2009cm\u22121$. We, therefore, attribute the maximum of $IAD(\omega )$ to the vibronic resonance. On the other hand, when *λ*_{ph} = 500 cm^{−1}, the maximum at *ω* ∼ 400 cm^{−1} has now disappeared, and $IAD(\omega )$ shows only a moderate variation throughout the entire frequency range. This is because the vibronic resonance has been quenched by the strong quasi-static disorder arising from the low-frequency part of the spectral density.^{36}

Finally, the reorganization energy density *J*_{B}(*ω*)/*ω* [Fig. 7(c)] shows a much simpler frequency dependence compared to $IAD(\omega )$. For *λ*_{ph} = 10 cm^{−1}, *J*_{B}(*ω*)/*ω* has a maximum at *ω* = 0 cm^{−1} and monotonically decreases as *ω* increases. The result for *λ*_{ph} = 500 cm^{−1} is just a scaled-up version of the *λ*_{ph} = 10 cm^{−1} result by a constant.

We are now ready to analyze the trend in $JBAD(\omega )$ by combining the effects of both $IAD(\omega )$ and *J*_{B}(*ω*)/*ω*. When *λ*_{ph} = 10 cm^{−1}, the vibronic resonance is so prevalent that $JBAD(\omega )$ shows a similar behavior as $IAD(\omega )$, while the additional weighting by *J*_{B}(*ω*)/*ω* does not significantly affect the overall trend. Nevertheless, the concentration of *J*_{B}(*ω*)/*ω* in the low-frequency region creates a minor maximum of $JBAD(\omega )$ at *ω* = 0 cm^{−1} and also slightly shifts the main peak from 400 cm^{−1} [Fig. 7(b)] to 380 cm^{−1} [Fig. 7(a)]. Such a shift can be interpreted as the Lamb shift associated with the bath reorganization. On the other hand, when *λ*_{ph} = 500 cm^{−1}, $JBAD(\omega )$ is primarily determined by *J*_{B}(*ω*)/*ω* as $IAD(\omega )$ only exhibits small dependence on the bath frequency.

Based on the discussions above, we can now also conclude that it is the difference between the donor–acceptor energy gap *E*_{D0} − *E*_{A0} and the vibronic resonance frequency $(ED0\u2212EA0)2+4|VDA|2$ that makes the discrepancy of ∼12 cm^{−1} between the peak locations for vFRET and PBME-nH (Fig. 4). This reflects that FRET only treats localized excitation (Sec. II B) and, therefore, neglects the contribution of *V*_{DA} in vibronic resonance.

By applying our theory of dissipation to a donor–acceptor system with varying reorganization energies, we directly visualized that increasing the subsystem–bath coupling strength *λ*_{ph} changes the dominant dissipation pathway from the resonant, single-phonon processes to multi-phonon processes, as implied by earlier studies.^{33} Overall, vFRET offers a general method of quantifying dissipation through individual bath components, which allows a detailed investigation on the dissipative dynamics that is often difficult for typical quantum dynamics simulation methods.

## V. FINAL REMARKS

In this paper, we developed a theory that can efficiently calculate the dissipation into individual bath components, by extending FRET theory into vibronic description. To do so, we first re-assigned an individual vibrational mode to the subsystem and constructed vibronic analogs of the key equations of the FRET theory. Then, analytical summations of vibronic energy transfer yielded compact expressions for calculating the dissipation without directly accessing the vibrational manifold. Our theory was benchmarked against a mixed quantum–classical simulation method by using a model donor–acceptor complex, which showed that our vFRET theory provides semi-quantitative frequency-dependent decomposition of the overall dissipation. We also numerically showed that the dissipation satisfies the detailed balance. Finally, based on the dissipative spectral densities [Eq. (50)], we identified the three factors determining the dissipation rate and explained how the dissipation channel changes as the vibronic coupling increases.

A few additional remarks are in order. Our vFRET theory is a useful tool that lets us calculate the dissipation with a significantly smaller computational cost compared to existing methods. For example, in our work, vFRET was cheaper than PBME-nH by a factor of 10^{5}–10^{6} in terms of total central processing unit (CPU) time. Moreover, the electronic rate constants and dissipative spectral densities [Eqs. (13) and (49)] need to be evaluated only once before the time propagation, as their values remain constant throughout the entire simulation. As a result, the time propagation in vFRET usually takes only a small portion of the total computational cost so that extending the propagation time does not pose significant additional computational burden. This is in contrast to other quantum dynamics simulation methods whose total cost increases linearly (PBME-nH and HEOM) or even exponentially (stochastic methods) with time.

The vFRET theory can also serve as a starting point for developing more sophisticated theories of dissipation by applying the approach of this work to similar Markovian theories such as Redfield theory,^{49} modified Redfield theory,^{50} or multi-chromophore version of FRET theory.^{51} Even non-Markovian effects may be captured by non-Markovian generalizations of FRET theory,^{52} modified Redfield theory,^{53,54} or hybridizing vFRET with Ehrenfest dynamics.^{48} With such extensions, we will be able to study many interesting aspects of dissipation, such as its dependence on the temperature or the functional form of *J*(*ω*). In addition, although vFRET predicts that the donor and acceptor identically contribute to dissipation, there may be some degree of asymmetry that can be only captured by more advanced formulations. Development of more accurate benchmarks, possibly based on the state-of-the-art mixed quantum classical methods,^{30,31} will benefit assessing the applicability and limitations of the newly developed theories.

We envision that the vFRET theory and its possible extensions will enable quantitative investigation on dissipation in a wider variety of molecular systems, with larger spatial and temporal scales than those previously studied. The insights gained from such a detailed level of studies will provide deeper understandings of the connection between the bath structure and dissipation in various molecular systems such as photosynthetic complexes, J- and H-aggregates, charge transfer complexes, electroluminescent materials, and photovoltaic devices. We hope pursuing this direction will enable us to control and design the open quantum system dynamics by properly engineering the subsystem–bath coupling, ultimately contributing to developing new molecular systems with enhanced properties.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the modified PBME-nH dissipation that includes the contribution of the RAVM reorganization energy and $J(\omega )$ and $I(\omega )$ for *λ*_{ph} = 50 cm^{−1} and *λ*_{ph} = 200 cm^{−1}.

## ACKNOWLEDGMENTS

The authors thank Professor Frank Huo for fruitful discussions about interpreting the results, the University of Rochester for support of this work through a PumpPrimer II award, and the National Science Foundation under Grant No. CHE-1553939.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: CONSISTENCY BETWEEN ELECTRONIC AND VIBRONIC POPULATION TRANSFER

In this appendix, we demonstrate the physical equivalence between electronic and vibronic descriptions by showing that the sum of vibronic population transfer over all vibrational state pairs (*I*, *J*) recovers the electronic population transfer. This is expressed as

for all three types of vibronic rate constants in Eq. (28). As in Sec. II D, analytical summations over *I* and *J* for Eqs. (A1a) and (A1b) can be done by using Eqs. (41) and (43),

Finally, proving Eq. (A1c) is rather trivial compared to the previous two equations, due to the orthonormality between the ground-state vibrational eigenstates of the RAVM,

### APPENDIX B: SPECTRAL LINE SHAPES DRESSED BY A THERMAL PHONON

In this appendix, we derive the expressions for spectral line shapes dressed by a thermal phonon [Eqs. (38)–(40)], which is a generalization of the result of Ref. 34 to nonzero-temperature cases. This is done by starting from a purely electronic Hamiltonian and adding the vibrational modes one by one, instead of directly calculating spectral line shapes from the BSD by using Eq. (15). Specifically, we prove the relation for the absorption line shape, Eq. (38b), while the proof for the fluorescence line shape [Eq. (38a)] can be constructed in a similar manner.

Consider a two-state molecular Hamiltonian that is not coupled to any vibrational modes,

where |*g*⟩ and |*e*⟩ denote the ground and excited electronic states, and *E*_{g} and *E*_{e} their electronic energies. Suppose the molecule is initially at |*g*⟩. As there is only one possible transition between the ground and excited states, the absorption line shape of $\u01240$ is monochromatic,

where $\mu ^eg$ is the transition dipole operator, *ℏ*Ω_{eg} = *E*_{e} − *E*_{g} is the energy gap between the ground and excited states, and *k* is the proportionality constant.

We now add the Hamiltonian of a vibrational mode to $\u01240$,

where $\u0124vn$ is

with the vibrational mode index *n*. Here, $p^n$ and $x^n$ are the position and momentum operators of the vibrational mode, *ω*_{n} is the characteristic frequency of the vibrational mode, and *d*_{n} is the displacement between the ground and excited state PESs along *x*_{n}. To construct the absorption line shape of $\u01241$, we consider its eigenbases {|*e*, *I*_{1}⟩ ≡ |*e*⟩|*I*_{1,e}⟩} and {|*g*, *J*_{1}⟩ ≡ |*g*⟩|*J*_{1,g}⟩}, where |*I*_{1,e}⟩ and |*J*_{1,g}⟩ are the vibrational eigenstates of the excited and ground-state PESs of the vibrational mode. Then, the vibronic transition energy from |*g*, *J*_{1}⟩ to |*e*, *I*_{1}⟩ is

Let us assume that the molecule is in thermal equilibrium on the ground-state PES, which is described by the density matrix $\u2211J1wJ1|g,J1\u2009\u2009g,J1|$ where the vibrational state populations ${wJ1}$ obey Boltzmann distribution [Eq. (33)]. With these considerations, the absorption line shape of $\u01241$ is constructed as a thermally weighted sum of individual vibronic transitions whose frequencies are determined by Eq. (B5),

where we have used the fact that $\mu ^eg$ is an operator in the electronic subspace so that $\u2009e,I1|\mu ^eg|g,J1\u2009=\u2009e|\mu ^eg|g\u2009\u27e8I1,e|J1,g\u27e9$. It is now straightforward to show that

where

is the phonon sideband of the vibrational mode. If we continue in this direction by adding more vibrational modes, we can generalize Eqs. (B7) and (B8) to an arbitrary number of vibrational modes,

where $An(\omega )$ is the absorption line shape of a molecule coupled to *n* vibrational modes, represented by the Hamiltonian

and $Ln(\omega )$ is the phonon sideband of the *n*th vibrational mode,

The physical meaning of Eqs. (B9)–(B11) can be summarized as follows: The addition of $\u0124vn$ to $\u0124n\u22121$ dresses $An\u22121(\omega )$, and Eqs. (B9) and (B11) state that the effect of the dressing can be implemented in the line shape by convoluting $An\u22121(\omega )$ with $Ln(\omega )$.