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.

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 models6–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 methods28–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) theory32 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.

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.

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

Ĥ=Ĥe+Ĥn+Ĥe-n+Ĥλ+ĤCoul.
(1)

Here,

Ĥe=A=1NEA0|AA|
(2)

is the bare electronic Hamiltonian where EA0 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,

Ĥn=A=1Nkp^Ak22+ωAk22x^Ak2,
(3)

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 kth vibrational mode coupled to molecule A. In turn, the vibronic coupling Hamiltonian

Ĥe-n=A=1N|AA|kωAk2dAkx^Ak
(4)

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

Ĥλ=A=1NλA|AA|,
(5)

where λA=k(ω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,

ĤCoul=A=1NB<AVAB|AB|+H.c.
(6)

is the Coulombic coupling Hamiltonian that describes the inter-molecular coupling. Here, VAB = VBA 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,

JA(ω)=kωAk2sAkδ(ωωAk),
(7)

where

sAk=ωAkdAk22
(8)

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

λAk=ωAk2dAk22=ωAksAk,
(9)

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

λA=0JA(ω)ωdω.
(10)

The total Hamiltonian [Eq. (1)] is treated in the framework of perturbation theory by writing Ĥ=Ĥ0+Ĥ1, where

Ĥ0=Ĥe+Ĥn+Ĥe-n+Ĥλ,Ĥ1=ĤCoul.
(11)

FRET theory supposes that: (1) The inter-molecular coupling VAB is much smaller than either the excitation energy gap |EA0EB0| or the standard deviation of the vibronic coupling Ĥe-n21/2 so that Ĥ1 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 A|Ĥ0|A. 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 VAB.

By applying Nakajima–Mori–Zwanzig projection operator techniques under these assumptions and keeping terms up to second-order in VAB,33 one arrives at

ṖA(t)=BA[KBAPA(t)+KABPB(t)],
(12)

where PA(t) is the population of the state |A⟩ at time t and KBA 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 Ĥ1 is small. The FRET rate constants are evaluated according to

KBA=2|VBA|22 Re0FA*(t)AB(t)dt,
(13)

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

FA(t)expit(EA0λA)gA*(t),AB(t)expit(EB0+λB)gB(t).
(14)

Here, gA(t) is the line broadening function,

gA(t)=10JA(ω)cothβω21cos(ωt)ω2+isin(ωt)ωtω2dω,
(15)

where β = 1/kBT and kB is Boltzmann’s constant. The physical meaning of Eq. (13) becomes more apparent if we express it with the Fourier-transformed quantities,

FA(ω)=12πFA(t)eiωtdt,
(16a)
AB(ω)=12πAB(t)eiωtdt,
(16b)

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*(t) and A(t)=A*(t), F(ω) and A(ω) are real quantities. Rewriting Eq. (13) by using Eq. (16) leads to

KBA=2π|VBA|22FA(ω)AB(ω)dω.
(17)

Equation (17) shows that the FRET rate constant KAB is proportional to the overlap integral between the donor fluorescence line shape FA(ω) and the acceptor absorption line shape AB(ω).

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

Ĥ=Ĥe+ĤCoul+AkĤAk,
(18)

where, in the right-hand side of Eq. (18), Ĥe+ĤCoul is the electronic contribution, while A,kĤAk is the vibrational contribution. The Hamiltonian component ĤAk for an individual vibrational mode of index Ak is

ĤAk=|AA|ĥAk,e+BA|BB|ĥAk,g,
(19)

where

ĥAk,g=p^Ak22+ωAk22x^Ak2,
(20a)
ĥAk,e=p^Ak22+ωAk22(x^Ak+dAk)2
(20b)

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 {|Ie⟩} of ĥAk,e and all other electronic states {|B⟩} with the ground-state vibrational eigenstates {|Jg⟩} of ĥAk,g. According to the second and third assumptions of FRET, the RAVM is always in the local thermal equilibrium of the nuclear PESs A|ĤAk|A=ĥAk,e or B|ĤAk|B=ĥAk,g. Therefore, an advantage of using the {|A⟩|Ie⟩} and {|B⟩|Jg⟩} 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⟩|Ie⟩} and {|B, J⟩ ≡ |B⟩|Jg⟩}. 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)],

Ĥ=Ĥe+Ĥn+Ĥe-n+Ĥλ+ĤCoul,
(21)

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

Ĥe=EA0+ωAvI+12|A,IA,I|+BAEB0+ωAvJ+12|B,JB,J|,
(22)

where we have added the vibrational eigenenergies of the RAVM to EA0 and EB0. The quantity

Ĥn=kvp^Ak22+ωAk22x^Ak2+BAkp^Bk22+ωBk22x^Bk2
(23)

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

Ĥe-n=I|A,IA,I|kvωAk2dAkx^Ak+BAJ|B,JB,J|kωBk2dBkx^Bk
(24)

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, Ĥe-n still adopts a direct product form as in Eq. (4). The reorganization Hamiltonian is now

Ĥλ=IλA|A,IA,I|+BAJλB|B,JB,J|,
(25)

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

ĤCoul=B<AI,JVABIe|Jg|A,IB,J|+BACAC<BJ,JVBCδJJ|B,JC,J|+H.c..
(26)

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.

FIG. 1.

The population transfer from molecule A to molecule B depicted in (a) the electronic description and (b) the vibronic description formed by explicit treatment of re-assigned vibrational mode (RAVM).

FIG. 1.

The population transfer from molecule A to molecule B depicted in (a) the electronic description and (b) the vibronic description formed by explicit treatment of re-assigned vibrational mode (RAVM).

Close modal

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,

ṖAI(t)=BAJ[KBJ,AIPAI(t)+KAI,BJPBJ(t)],
(27a)
ṖBJ(t)=I[KAI,BJPBJ(t)+KBJ,AIPAI(t)]+CACBJ[KCJ,BJPBJ(t)+KBJ,CJPCJ(t)].
(27b)

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 KAI,AI = KBJ,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 {KBJ,AI} and {KAI,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, {KBJ,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),

KBJ,AI=2|VBA|22|Jg|Ie|2 Re0FAI*(t)ABJ(t)dt,KAI,BJ=2|VAB|22|Ie|Jg|2 Re0FBJ*(t)AAI(t)dt,KCJ,BJ=2|VBC|22δJJ Re0FBJ*(t)ACJ(t)dt,
(28)

can be constructed from the electronic rate constant [Eq. (13)] by replacing the electronic coupling VAB 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

FAI(t)=expitEA0+ωAvI+12λAgA*(t),
(29a)
AAI(t)=expitEA0+ωAvI+12+λAgA(t),
(29b)
FBJ(t)=expitEB0+ωAvJ+12λBgB*(t),
(29c)
ABJ(t)=expitEB0+ωAvJ+12+λBgB(t).
(29d)

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 gA(t) in Eq. (14) are replaced by λA and gA(t) defined by

gA(t)gA(t)gAv(t),
(30)

where gAv(t) is the line broadening function of the RAVM,

gAv(t)=sAvcothβωAv2[1cos(ωAvt)]+i[sin(ωAvt)ωAvt],
(31)

evaluated by plugging the spectral density of the RAVM JAv(ω)=ωAv2sAvδ(ωωAv) [Eq. (7)] into Eq. (15). Meanwhile, λB and gB(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

PAI(t)=wIPA(t),PBJ(t)=wJPB(t),
(32)

where wI and wJ are the Boltzmann weights for the vibrational eigenstates of the RAVM,

wI=1ZexpβωAvI+12,wJ=1ZexpβωAvJ+12,
(33)

with the partition function Z=[2sinh(βωAv/2)]1. Equations (32) and (33) are not automatically satisfied by Eqs. (27) and (28), and the vibronic populations within each molecule must be redistributed according to Eqs. (32) and (33) after every integration time step.

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,

EAv(t)=IωAvI+12PAI(t)+BAJωAvJ+12PBJ(t),
(34)

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

EAv(t)=IwIωAvI+12,
(35)

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,

ĖAv(t)=BAI,JωAv(IJ)KBJ,AIPAI(t)+KAI,BJPBJ(t)+BACAC<BJ,JωAv(JJ)[KCJ,BJPBJ(t)+KBJ,CJPCJ(t)].
(36)

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

ĖAv(t)=2|VBA|22BAPA(t)I,JwI(IJ)ωAv|Jg|Ie|2×Re0FAI*(t)ABJ(t)dt+2|VAB|22BAPB(t)×I,JwJ(IJ)ωAv|Ie|Jg|2Re0  FBJ*(t)AAI(t)dt.
(37)

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(ω) and AA(ω) are related to the line shapes without the RAVM, FA(ω) and AA(ω), through the following convolution relations:34 

FA(ω)=dωFA(ω)LAv(ωω),
(38a)
AA(ω)=dωAA(ω)LAv(ωω),
(38b)

where FA(ω) and AA(ω) are constructed from Eqs. (14) and (16) by making substitutions λAλA and gA(t) → gA(t),

FA(ω)=12πexpit(EA0λA)gA*(t)eiωtdt,AA(ω)=12πexpit(EA0+λA)gA(t)eiωtdt,
(39)

and LAv(ω) is the phonon sideband of the RAVM,

LAv(ω)=I,JwJ|Ie|Jg|2δ(ω(IJ)ωAv).
(40)

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

I,JwI|Jg|Ie|2FA(IJ)(t)=FA(t)eiωAvt2,I,JwJ|Ie|Jg|2AA(IJ)(t)=AA(t)eiωAvt2.
(41)

To use Eq. (41), however, we first need to remove IJ 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 IJ in the exponential. As a result, we can derive the identities

(IJ)ωAvFAI*(t)ABJ(t)  =(EA0λAEB0λB)iġA(t)+ġB(t)+ddt   ×FAI*(t)ABJ(t),
(42a)
(IJ)ωAvFBJ*(t)AAI(t)  =(EA0+λAEB0+λB)+iġA(t)+ġB(t)+ddt   ×FBJ*(t)AAI(t).
(42b)

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

FAI*(t)ABJ(t)=FA(IJ)*(t)AB(t)eiωAvt2,FBJ*(t)AAI(t)=FB*(t)AA(IJ)(t)eiωAvt2.
(43)

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

I,JwI(IJ)ωAv|Jg|Ie|2FAI*(t)ABJ(t)  =λAviġAv(t)FA*(t)AB(t),I,JwJ(IJ)ωAv|Ie|Jg|2FBJ*(t)AAI(t)  =λAviġAv(t)FB*(t)AA(t).
(44)

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

λAviġAv(t)=λAvcos(ωt)icothβω2sin(ωt).
(45)

Now, from Eqs. (44) and (45), we have

I,JwI(IJ)ωAv|Jg|Ie|2 Re0FAI*(t)ABJ(t)dt  =IBA(ωAv)λAv,I,JwJ(IJ)ωAv|Ie|Jg|2 Re0FBJ*(t)AAI(t)dt  =IAB(ωAv)λAv,
(46)

where IBA(ω) and IAB(ω) are

IBA(ω)=Re 0FA*(t)AB(t)×cos(ωt)icothβω2sin(ωt)dt,IAB(ω)=Re 0FB*(t)AA(t)×cos(ωt)icothβω2sin(ωt)dt.
(47)

Plugging Eq. (46) into Eq. (37) yields

ĖAv(t)=BA[KAvBAPA(t)+KAvABPB(t)],
(48)

with the dissipation rate constants KAvBA and KAvAB given by

KAvBA=2|VBA|22IBA(ωAv)λAv,KAvAB=2|VAB|22IAB(ωAv)λAv.
(49)

To generalize Eqs. (48) and (49) to treat continuous BSDs, based on Eq. (10), we make substitutions ωAvω and λAvJA(ω)ωdω in Eq. (49) and define “dissipative spectral densities” JABA(ω) and JAAB(ω) as follows:

JABA(ω)2|VBA|22IBA(ω)JA(ω)ω,JAAB(ω)2|VAB|22IAB(ω)JA(ω)ω.
(50)

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

DA(ω,t)dω=BAJABA(ω)PA(t)+JAAB(ω)PB(t)dω.
(51)

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 {PA(t)}, so the calculations can be done without directly accessing the vibrational manifold. Because Eq. (51) contains {PA(t)}, calculating dissipation requires knowledge of the electronic dynamics. In practice, both the electronic rate constants {KAB} [Eq. (13)] and dissipative spectral densities {JABA(ω),JAAB(ω)} [Eq. (50)] are evaluated at the start of the simulation. Then, {PA(t)} and {DA(ω,t)} are propagated simultaneously by using Eqs. (12) and (51).

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 λAk=ω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.

As a minimal model that allows simple interpretations, we consider a donor–acceptor complex whose excitation energy difference is ED0EA0 = 400 cm−1 and the electronic couplings are VDA = VAD = 50 cm−1. As |ED0EA0| ≫ VDA, 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,

JD(ω)=JA(ω)=2λphπωcωω2+ωc2,
(52)

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(ω)} [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 105 for simulating the electronic dynamics, 106 for the dissipation calculation with λph = 10 cm−1, and 107 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 function41 and the hierarchy depth of 7. The integration time step was 0.1 fs.

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,

PD(t)=PD()+1PD()eγt,PA(t)=PA()(1eγt).
(53)

Here, PD(∞) and PA(∞) 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 

PD()PA()=eβ(ED0EA0),
(54)

and, hence, is not affected by λph.

FIG. 2.

Electronic and dissipation dynamics. (a) Electronic populations of the donor (solid lines) and acceptor (dashed lines) calculated by FRET and PBME-nH. (b) Energy dissipated into the donor vibrational modes in the [40, 50] cm−1 window (solid line) calculated by vFRET and PBME-nH. The dashed line shows correction due to the spurious linear drift in PBME-nH. The circle on the vertical axis marks the corrected dissipation at the long time limit, while Eq. (57) generalizes this idea to arbitrary t.

FIG. 2.

Electronic and dissipation dynamics. (a) Electronic populations of the donor (solid lines) and acceptor (dashed lines) calculated by FRET and PBME-nH. (b) Energy dissipated into the donor vibrational modes in the [40, 50] cm−1 window (solid line) calculated by vFRET and PBME-nH. The dashed line shows correction due to the spurious linear drift in PBME-nH. The circle on the vertical axis marks the corrected dissipation at the long time limit, while Eq. (57) generalizes this idea to arbitrary t.

Close modal

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 PA(∞), which shows that PA(∞) from FRET simulations show excellent agreement with HEOM results and remain constant for all values of λph, as predicted by Eq. (54). By contrast, PA(∞) 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 DOFs39,43,44 as in PBME-nH.

FIG. 3.

The electronic dynamics calculated by HEOM, FRET, and PBME-nH simulations while varying the bath reorganization energy λph from 10 cm−1 to 500 cm−1. (a) The equilibrium electronic population for the acceptor molecule PA(∞) and (b) the population transfer rate γ extracted by fitting the time-dependent population to the exponential model [Eq. (53)].

FIG. 3.

The electronic dynamics calculated by HEOM, FRET, and PBME-nH simulations while varying the bath reorganization energy λph from 10 cm−1 to 500 cm−1. (a) The equilibrium electronic population for the acceptor molecule PA(∞) and (b) the population transfer rate γ extracted by fitting the time-dependent population to the exponential model [Eq. (53)].

Close modal

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.

1. Time profile

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

EB(ω,t)=0tDB(ω,t)dt,
(55)

which quantifies the dissipation by JB(ω) 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(ω,t) by using the fact that the net dissipation must vanish at equilibrium [DB(ω,)=0], and finally integrating the resulting expression by using Eq. (55), consequently,

EB(ω,t)=JBAD(ω)γ(1eγt).
(56)

For our model with JD(ω) = JA(ω), Eq. (56) predicts that the dissipation into the donor and acceptor will be identical as JDAD(ω)=JAAD(ω) [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(ω,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

ẼB(ω,t)=EB(ω,t)tlimτĖB(ω,τ).
(57)

Equation (57) assumes that the contribution from the drift linearly increases during the entire simulation period. All PBME-nH dissipation results discussed below are corrected by Eq. (57).

2. Frequency-dependent decomposition

In this section, we decompose the total dissipation in the frequency domain, ẼD(ω,)+ẼA(ω,). Figure 4 shows ẼD(ω,)+ẼA(ω,) 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.

FIG. 4.

Frequency-dependent decompositions of dissipation in the frequency domain by (a) vFRET and (b) PBME-nH, for different values of λph.

FIG. 4.

Frequency-dependent decompositions of dissipation in the frequency domain by (a) vFRET and (b) PBME-nH, for different values of λph.

Close modal

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 VDA, thus slowing down the electronic dynamics (Fig. 5). As expected, while decreasing VDA 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.

FIG. 5.

Frequency-dependent decompositions of dissipation for λph = 200 cm−1 with (a) VDA = 50 cm−1 and (b) VDA = 10 cm−1, calculated by vFRET and PBME-nH.

FIG. 5.

Frequency-dependent decompositions of dissipation for λph = 200 cm−1 with (a) VDA = 50 cm−1 and (b) VDA = 10 cm−1, calculated by vFRET and PBME-nH.

Close modal

Finally, the λph = 10 cm−1 results also reveal that the peaks in ẼD(ω,)+ẼA(ω,) in Fig. 4 are located at slightly different frequencies for vFRET and PBME-nH. The origin of such a difference is also addressed in Sec. IV.

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 ṖD()=ṖA()=0 and DD(ω,)=DA(ω,)=0. By plugging these conditions into Eq. (51), we get

PA()PD()=JBAD(ω)JBDA(ω).
(58)

Equation (58) shows that JBAD(ω) and JBDA(ω) 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(ω)=JAAD(ω) and JDDA(ω)=JADA(ω) (Sec. III C 2). Because the equilibrium electronic populations satisfy detailed balance,33 all the ratios in Eq. (58) must be equal to exp[β(ED0EA0)] = 6.810 143 83… for our simulation conditions.

Figure 6 displays JBAD(ω)/JBDA(ω) calculated for four different values of λph. We observe that JBAD(ω)/JBDA(ω) is practically equal to exp[β(ED0EA0)] 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.

FIG. 6.

Ratios between the dissipative spectral densities [JBAD(ω)/JBDA(ω)], evaluated for different values of λph. The black line corresponds to the Boltzmann population ratio eβ(ED0EA0) that represents the detailed balance condition.

FIG. 6.

Ratios between the dissipative spectral densities [JBAD(ω)/JBDA(ω)], evaluated for different values of λph. The black line corresponds to the Boltzmann population ratio eβ(ED0EA0) that represents the detailed balance condition.

Close modal

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(ω) and JBDA(ω) [Eq. (50)] that are directly proportional to EB(ω,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 VDA as it directly affects the speed of FRET [Eq. (13)]. (ii) The “dissipative potential” I(ω) [Eq. (47)], which quantifies the ability of the bath mode of a certain frequency to participate in the dissipation. (iii) The reorganization energy density JB(ω)/ω 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 VDA is fixed in our model, we consider how JB(ω) is affected by I(ω) and JB(ω)/ω.

Figure 7(a) displays JBAD(ω) and JBDA(ω). 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(ω) as JBAD(ω) and JBDA(ω) only differ by a constant prefactor [Fig. 6(b) and Eq. (58)]. The frequency dependence of JBAD(ω) 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(ω) and JB(ω)/ω.

FIG. 7.

Analysis of the dissipative spectral densities in the frequency domain, for λph = 10 cm−1 and λph = 500 cm−1. (a) The dissipative spectral densities JBAD(ω) and JBDA(ω), (b) the dissipative potential IAD(ω) and IDA(ω), and (c) the reorganization energy density JB(ω)/ω.

FIG. 7.

Analysis of the dissipative spectral densities in the frequency domain, for λph = 10 cm−1 and λph = 500 cm−1. (a) The dissipative spectral densities JBAD(ω) and JBDA(ω), (b) the dissipative potential IAD(ω) and IDA(ω), and (c) the reorganization energy density JB(ω)/ω.

Close modal

The frequency dependence of the dissipative potentials IAD(ω) and IDA(ω) is plotted in Fig. 7(b). When λph = 10 cm−1, IAD(ω) peaks at ω ∼ 400 cm−1, which coincides with the donor–acceptor energy gap ED0EA0. As VDA ≪ |ED0EA0| in our system, ED0EA0 = 400 cm−1 is close to the eigenenergy gap of Ĥe+ĤCoul where vibronic resonance occurs, which is (ED0EA0)2+4|VDA|2412cm1. We, therefore, attribute the maximum of IAD(ω) to the vibronic resonance. On the other hand, when λph = 500 cm−1, the maximum at ω ∼ 400 cm−1 has now disappeared, and IAD(ω) 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 JB(ω)/ω [Fig. 7(c)] shows a much simpler frequency dependence compared to IAD(ω). For λph = 10 cm−1, JB(ω)/ω 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(ω) by combining the effects of both IAD(ω) and JB(ω)/ω. When λph = 10 cm−1, the vibronic resonance is so prevalent that JBAD(ω) shows a similar behavior as IAD(ω), while the additional weighting by JB(ω)/ω does not significantly affect the overall trend. Nevertheless, the concentration of JB(ω)/ω in the low-frequency region creates a minor maximum of JBAD(ω) 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(ω) is primarily determined by JB(ω)/ω as IAD(ω) 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 ED0EA0 and the vibronic resonance frequency (ED0EA0)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 VDA 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.

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 105–106 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.

See the supplementary material for the modified PBME-nH dissipation that includes the contribution of the RAVM reorganization energy and J(ω) and I(ω) for λph = 50 cm−1 and λph = 200 cm−1.

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.

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

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

I,JKBJ,AIPAI(t)=KBAPA(t),
(A1a)
I,JKAI,BJPBJ(t)=KABPB(t),
(A1b)
I,JKCJ,BJPBJ(t)=KBCPB(t),
(A1c)

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),

I,JKBJ,AIPAI(t)=2|VAB|22I,J|Jg|Ie|2PAI(t)× Re0FA(IJ)*(t)AB(t)eiωAvt2dt=2|VAB|22PA(t) Re0FA*(t)AB(t)dt=KBAPA(t),
(A2a)
I,JKAI,BJPBJ(t)=2|VAB|22I,J|Ie|Jg|2PBJ(t)× Re0FB*(t)AA(IJ)(t)eiωAvt2dt=2|VAB|22PB(t) Re0FB*(t)AA(t)dt=KABPB(t).
(A2b)

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,

J,JKCJ,BJPBJ(t)=J,JKCBδJJwJPB(t)=KCBPB(t).
(A3)

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,

Ĥ0=Eg|gg|+Ee|ee|,
(B1)

where |g⟩ and |e⟩ denote the ground and excited electronic states, and Eg and Ee 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 Ĥ0 is monochromatic,

A0(ω)=k|e|μ^eg|g|2δ(ωΩeg),
(B2)

where μ^eg is the transition dipole operator, Ωeg = EeEg 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 Ĥ0,

Ĥ1=Ĥ0+Ĥv1,
(B3)

where Ĥvn is

Ĥvn=|gg|p^n22+ωn22x^n2+|ee|p^n22+ωn22(x^n+dn)2,
(B4)

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 dn is the displacement between the ground and excited state PESs along xn. To construct the absorption line shape of Ĥ1, we consider its eigenbases {|e, I1⟩ ≡ |e⟩|I1,e⟩} and {|g, J1⟩ ≡ |g⟩|J1,g⟩}, where |I1,e⟩ and |J1,g⟩ are the vibrational eigenstates of the excited and ground-state PESs of the vibrational mode. Then, the vibronic transition energy from |g, J1⟩ to |e, I1⟩ is

Ee,I1Eg,J1=Ωeg+(I1J1)ω1.
(B5)

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

A1(ω)=k|e|μ^eg|g|2I1,J1wJ1|I1,e|J1,g|2  ×δ(ωΩeg(I1J1)ω1),
(B6)

where we have used the fact that μ^eg is an operator in the electronic subspace so that e,I1|μ^eg|g,J1=e|μ^eg|gI1,e|J1,g. It is now straightforward to show that

A1(ω)=dωA0(ω)L1(ωω),
(B7)

where

L1(ω)=I1,J1wJ1|I1,e|J1,g|2δ(ω(I1J1)ω1)
(B8)

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,

An(ω)=dωAn1(ω)Ln(ωω),
(B9)

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

Ĥn=Ĥ0+k=1nĤvk,
(B10)

and Ln(ω) is the phonon sideband of the nth vibrational mode,

Ln(ω)=In,JnwJn|In,e|Jn,g|2δ(ω(InJn)ωn).
(B11)

The physical meaning of Eqs. (B9)(B11) can be summarized as follows: The addition of Ĥvn to Ĥn1 dresses An1(ω), and Eqs. (B9) and (B11) state that the effect of the dressing can be implemented in the line shape by convoluting An1(ω) with Ln(ω).

We can now deduce Eqs. (38)(40) from Eqs. (B9) and (B11) by considering the RAVM as the nth mode.

1.
U.
Weiss
,
Quantum Dissipative Systems
, Series in Modern Condensed Matter Physics (
World Scientific
,
1999
).
2.
H.
Breuer
and
F.
Petruccione
,
The Theory of Open Quantum Systems
(
Oxford University Press
,
2002
).
3.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
, “
Dynamics of the dissipative two-state system
,”
Rev. Mod. Phys.
59
,
1
85
(
1987
).
4.
F.
Fassioli
,
R.
Dinshaw
,
P. C.
Arpin
, and
G. D.
Scholes
, “
Photosynthetic light harvesting: Excitons and coherence
,”
J. R. Soc., Interface
11
,
20130901
(
2014
).
5.
N. J.
Hestand
and
F. C.
Spano
, “
Expanded theory of H- and J-molecular aggregates: The effects of vibronic coupling and intermolecular charge transfer
,”
Chem. Rev.
118
,
7069
7163
(
2018
).
6.
A.
Ishizaki
and
G. R.
Fleming
, “
Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature
,”
Proc. Natl. Acad. Sci. U. S. A.
106
,
17255
17260
(
2009
).
7.
E. Y.
Wilner
,
H.
Wang
,
M.
Thoss
, and
E.
Rabani
, “
Sub-ohmic to super-ohmic crossover behavior in nonequilibrium quantum systems with electron–phonon interactions
,”
Phys. Rev. B
92
,
195143
(
2015
).
8.
P.
Nalbach
,
C. A.
Mujica-Martinez
, and
M.
Thorwart
, “
Vibronically coherent speed-up of the excitation energy transfer in the Fenna–Matthews–Olson complex
,”
Phys. Rev. E
91
,
022706
(
2015
).
9.
C.
Duan
,
C.-Y.
Hsieh
,
J.
Liu
,
J.
Wu
, and
J.
Cao
, “
Unusual transport properties with noncommutative system-bath coupling operators
,”
J. Phys. Chem. Lett.
11
,
4080
4085
(
2020
).
10.
M. K.
Lee
,
K. B.
Bravaya
, and
D. F.
Coker
, “
First-principles models for biological light-harvesting: Phycobiliprotein complexes from cryptophyte algae
,”
J. Am. Chem. Soc.
139
,
7803
7814
(
2017
).
11.
C. W.
Kim
,
B.
Choi
, and
Y. M.
Rhee
, “
Excited state energy fluctuations in the Fenna–Matthews–Olson complex from molecular dynamics simulations with interpolated chromophore potentials
,”
Phys. Chem. Chem. Phys.
20
,
3310
3319
(
2018
).
12.
S.
Maity
,
B. M.
Bold
,
J. D.
Prajapati
,
M.
Sokolov
,
T.
Kubař
,
M.
Elstner
, and
U.
Kleinekathöfer
, “
DFTB/MM molecular dynamics simulations of the FMO light-harvesting complex
,”
J. Phys. Chem. Lett.
11
,
8660
8667
(
2020
).
13.
D.
Kienzler
,
H.-Y.
Lo
,
B.
Keitch
,
L.
de Clercq
,
F.
Leupold
,
F.
Lindenfelser
,
M.
Marinelli
,
V.
Negnevitsky
, and
J. P.
Home
, “
Quantum harmonic oscillator state synthesis by reservoir engineering
,”
Science
347
,
53
56
(
2014
).
14.
L.-Y.
Hsu
,
W.
Ding
, and
G. C.
Schatz
, “
Plasmon-coupled resonance energy transfer
,”
J. Phys. Chem. Lett.
8
,
2357
2367
(
2017
).
15.
J. A.
Campos-Gonzalez-Angulo
,
R. F.
Ribeiro
, and
J.
Yuen-Zhou
, “
Resonant catalysis of thermally activated chemical reactions with vibrational polaritons
,”
Nat. Commun.
10
,
4685
(
2019
).
16.
K.
Ng
,
M.
Webster
,
W. P.
Carbery
,
N.
Visaveliya
,
P.
Gaikwad
,
S. J.
Jang
,
I.
Kretzschmar
, and
D. M.
Eisele
, “
Frenkel excitons in heat-stressed supramolecular nanocomposites enabled by tunable cage-like scaffolding
,”
Nat. Chem.
12
,
1157
1164
(
2020
).
17.
E. J.
O’Reilly
and
A.
Olaya-Castro
, “
Non-classicality of the molecular vibrations assisting exciton energy transfer at room temperature
,”
Nat. Commun.
5
,
3012
(
2014
).
18.
V. I.
Novoderezhkin
,
E.
Romero
,
J.
Prior
, and
R.
van Grondelle
, “
Exciton-vibrational resonance and dynamics of charge separation in the photosystem II reaction center
,”
Phys. Chem. Chem. Phys.
19
,
5195
5208
(
2017
).
19.
C. W.
Kim
and
Y. M.
Rhee
, “
Toward monitoring the dissipative vibrational energy flows in open quantum systems by mixed quantum-classical simulations
,”
J. Chem. Phys.
152
,
244109
(
2020
).
20.
H.-D.
Meyer
,
U.
Manthe
, and
L. S.
Cederbaum
, “
The multi-configurational time-dependent Hartree approach
,”
Chem. Phys. Lett.
165
,
73
78
(
1990
).
21.
G. A.
Worth
, “
Accurate wave packet propagation for large molecular systems: The multiconfiguration time-dependent Hartree (MCTDH) method with selected configurations
,”
J. Chem. Phys.
112
,
8322
8329
(
2000
).
22.
H.
Wang
, “
Multilayer multiconfiguration time-dependent Hartree theory
,”
J. Phys. Chem. A
119
,
7951
7965
(
2015
).
23.
M. A.
Cazalilla
and
J. B.
Marston
, “
Time-dependent density-matrix renormalization group: A systematic method for the study of quantum many-body out-of-equilibrium systems
,”
Phys. Rev. Lett.
88
,
256403
(
2002
).
24.
S. R.
White
and
A. E.
Feiguin
, “
Real-time evolution using the density matrix renormalization group
,”
Phys. Rev. Lett.
93
,
076401
(
2004
).
25.
J.
Haegeman
,
C.
Lubich
,
I.
Oseledets
,
B.
Vandereycken
, and
F.
Verstraete
, “
Unifying time evolution and optimization with matrix product states
,”
Phys. Rev. B
94
,
165116
(
2016
).
26.
J.
Ren
,
Z.
Shuai
, and
G.
Kin-Lic Chan
, “
Time-dependent density matrix renormalization group algorithms for nearly exact absorption and fluorescence spectra of molecular aggregates at both zero and finite temperature
,”
J. Chem. Theory Comput.
14
,
5027
5039
(
2018
).
27.
T.
Jiang
,
W.
Li
,
J.
Ren
, and
Z.
Shuai
, “
Finite temperature dynamical density matrix renormalization group for spectroscopy in frequency domain
,”
J. Phys. Chem. Lett.
11
,
3761
3768
(
2020
).
28.
H.
Kim
,
A.
Nassimi
, and
R.
Kapral
, “
Quantum-classical Liouville dynamics in the mapping basis
,”
J. Chem. Phys.
129
,
084102
(
2008
).
29.
C.-Y.
Hsieh
and
R.
Kapral
, “
Nonadiabatic dynamics in open quantum-classical systems: Forward–backward trajectory solution
,”
J. Chem. Phys.
137
,
22A507
(
2012
).
30.
J.
Liu
and
G.
Hanna
, “
Efficient and deterministic propagation of mixed quantum-classical Liouville dynamics
,”
J. Phys. Chem. Lett.
9
,
3928
3933
(
2018
).
31.
A.
Kelly
,
N.
Brackbill
, and
T. E.
Markland
, “
Accurate nonadiabatic quantum dynamics on the cheap: Making the most of mean field theory with master equations
,”
J. Chem. Phys.
142
,
094110
(
2015
).
32.
T.
Förster
, “
Energiewanderung und fluoreszenz
,”
Naturwissenschaften
33
,
166
175
(
1946
).
33.
M.
Yang
and
G. R.
Fleming
, “
Influence of phonons on exciton transfer dynamics: Comparison of the Redfield, Förster, and modified Redfield equations
,”
Chem. Phys.
275
,
355
372
(
2002
).
34.
J. M.
Hayes
,
J. K.
Gillie
,
D.
Tang
, and
G. J.
Small
, “
Theory for spectral hole burning of the primary electron donor state of photosynthetic reaction centers
,”
Biochim. Biophys. Acta Bioenerg.
932
,
287
305
(
1988
).
35.
H. W.
Kim
and
Y. M.
Rhee
, “
Improving long time behavior of Poisson bracket mapping equation: A non-Hamiltonian approach
,”
J. Chem. Phys.
140
,
184106
(
2014
).
36.
C. W.
Kim
,
W.-G.
Lee
,
I.
Kim
, and
Y. M.
Rhee
, “
Effect of underdamped vibration on excitation energy transfer: Direct comparison between two different partitioning schemes
,”
J. Phys. Chem. A
123
,
1186
1197
(
2019
).
37.
Y.
Tanimura
, “
Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath
,”
Phys. Rev. A
41
,
6676
6687
(
1990
).
38.
Y.
Tanimura
, “
Stochastic Liouville, Langevin, Fokker–Planck, and master equation approaches to quantum dissipative systems
,”
J. Phys. Soc. Jpn.
75
,
082001
(
2006
).
39.
H. W.
Kim
,
W.-G.
Lee
, and
Y. M.
Rhee
, “
Improving long time behavior of Poisson bracket mapping equation: A mapping variable scaling approach
,”
J. Chem. Phys.
141
,
124107
(
2014
).
40.
T.
Ikeda
and
G. D.
Scholes
, “
Generalization of the hierarchical equations of motion theory for efficient calculations with arbitrary correlation functions
,”
J. Chem. Phys.
152
,
204101
(
2020
).
41.
J.
Hu
,
R.-X.
Xu
, and
Y.
Yan
, “
Communication: Padé spectrum decomposition of Fermi function and Bose function
,”
J. Chem. Phys.
133
,
101106
(
2010
).
42.
H.
Wang
,
X.
Song
,
D.
Chandler
, and
W. H.
Miller
, “
Semiclassical study of electronically nonadiabatic dynamics in the condensed-phase: Spin-boson problem with Debye spectral density
,”
J. Chem. Phys.
110
,
4828
4840
(
1999
).
43.
U.
Müller
and
G.
Stock
, “
Flow of zero-point energy and exploration of phase space in classical simulations of quantum relaxation cynamics. II. Application to nonadiabatic processes
,”
J. Chem. Phys.
111
,
77
88
(
1999
).
44.
G.
Tao
and
W. H.
Miller
, “
Semiclassical description of electronic excitation population transfer in a model photosynthetic system
,”
J. Phys. Chem. Lett.
1
,
891
894
(
2010
).
45.
P.
Rebentrost
,
M.
Mohseni
,
I.
Kassal
,
S.
Lloyd
, and
A.
Aspuru-Guzik
, “
Environment-assisted quantum transport
,”
New J. Phys.
11
,
033003
(
2009
).
46.
E.
Zerah-Harush
and
Y.
Dubi
, “
Universal origin for environment-assisted quantum transport in exciton transfer networks
,”
J. Phys. Chem. Lett.
9
,
1689
1695
(
2018
).
47.
F.
Caruso
,
A. W.
Chin
,
A.
Datta
,
S. F.
Huelga
, and
M. B.
Plenio
, “
Highly efficient energy excitation transfer in light-harvesting complexes: The fundamental role of noise-assisted transport
,”
J. Chem. Phys.
131
,
105106
(
2009
).
48.
A.
Montoya-Castillo
,
T. C.
Berkelbach
, and
D. R.
Reichman
, “
Extending the applicability of Redfield theories into highly non-Markovian regimes
,”
J. Chem. Phys.
143
,
194108
(
2015
).
49.
A. G.
Redfield
, “
On the theory of relaxation processes
,”
IBM J. Res. Dev.
1
,
19
31
(
1957
).
50.
W. M.
Zhang
,
T.
Meier
,
V.
Chernyak
, and
S.
Mukamel
, “
Exciton-migration and three-pulse femtosecond optical spectroscopies of photosynthetic antenna complexes
,”
J. Chem. Phys.
108
,
7763
7774
(
1998
).
51.
S.
Jang
,
M. D.
Newton
, and
R. J.
Silbey
, “
Multichromophoric Förster resonance energy transfer
,”
Phys. Rev. Lett.
92
,
218301
(
2004
).
52.
S.
Jang
,
S.
Hoyer
,
G.
Fleming
, and
K. B.
Whaley
, “
Generalized master equation with non-Markovian multichromophoric Förster resonance energy transfer for modular exciton densities
,”
Phys. Rev. Lett.
113
,
188102
(
2014
).
53.
Y.-H.
Hwang-Fu
,
W.
Chen
, and
Y.-C.
Cheng
, “
A coherent modified Redfield theory for excitation energy transfer in molecular aggregates
,”
Chem. Phys.
447
,
46
53
(
2015
).
54.
Y.
Chang
and
Y.-C.
Cheng
, “
On the accuracy of coherent modified Redfield theory in simulating excitation energy transfer dynamics
,”
J. Chem. Phys.
142
,
034109
(
2015
).

Supplementary Material