Calculation of the equilibrium state of an open quantum system interacting with a bath remains a challenge to this day, mostly due to a huge number of bath degrees of freedom. Here, we present an analytical expression for the reduced density operator in terms of an effective Hamiltonian for a high temperature case. Comparing with numerically exact results, we show that our theory is accurate for slow baths and up to intermediate system–bath coupling strengths. Our results demonstrate that the equilibrium state does not depend on the shape of spectral density in the slow bath regime. The key quantity in our theory is the effective coupling between the states, which depends exponentially on the ratio of the reorganization energy to temperature and, thus, has opposite temperature dependence than could be expected from the small polaron transformation.

The equilibrium state of open quantum systems has attracted significant scientific attention from the early days of quantum theory1 and remains of considerable interest in the present day.2–7 This issue is relevant in a variety of contexts, from initial states of emission spectra in molecular systems,8,9 multichromophoric Förster transfer rates,10–12 nonequilibrium steady states in quantum transport problems13 or quantum heat engines,14 and decoherence problems in quantum information science15 to quantum corrections of thermodynamic laws.16 Unfortunately, analysis of the equilibrium state is hindered by an enormous amount of degrees of freedom of the bath.

When system–bath coupling is negligible, the equilibrium reduced density operator ρ^eq of the system is represented by the canonical Boltzmann distribution of the bare quantum system. On the other hand, this is no longer the case if the coupling between the system and the bath becomes stronger. These issues have been investigated both analytically and numerically,3–5,17 and there exist straightforward though numerically somewhat expensive methods to calculate ρ^eq exactly.3,6 Unfortunately, no simple analytical formulas are currently available in the literature. This is rectified in this letter, where we derive an analytical expression for the reduced equilibrium density operator that is valid for high temperatures and slow baths. Our result is intuitively understandable, providing a mental picture for this complex issue.

Previously, we have demonstrated non-trivial spectroscopic manifestations of the non-canonical effects and showed numerically that the effects of the bath upon the equilibrium of the system can be accounted for by introducing effective coupling between system states, Jeff.2 It was shown for a molecular dimer that the effective coupling could be well described by the relation

Jeff12=Jexpaλ1+λ2kBT,
(1)

with J denoting the original bare coupling, λ1,2 being the reorganization energies of the system states, and a being a numerical constant, the value of which was determined by the fitting procedure to be ∼0.15. Here, we provide a theoretical justification for these results.

Let us start by defining the total Hamiltonian, which is comprised of the system, bath, and system–bath interaction parts, respectively,

Ĥ=ĤS+ĤB+ĤSB.
(2)

The Hamiltonian of the system can be expressed as ĤS = Ĥε + ĤJ, with

Ĥε=nNεn|nn|
(3)

denoting the energies of the states (or sites) and

ĤJ=m,nmnJmn|mn|
(4)

denoting the couplings between the states. N numbers the states of the system.

The bath consists of an infinite set of harmonic oscillators

ĤB=jωj2p^j2+x^j2,
(5)

which interact linearly with system states,

ĤSB=nNjωjdnjx^j|nn|=nNF^n|nn|.
(6)

Here, ωj, x^j, p^j are the frequency, dimensionless position, and dimensionless momentum of the jth bath oscillator, while dnj is a dimensionless coupling between the jth bath oscillator and state n. We have also introduced bath operators F^n=jωjdnjx^j. Note that we set = 1 throughout the paper. System–bath interaction induces fluctuations of system state energies that can be parameterized by the spectral densities

Cmnω=πjωj2dmjdnj2δωωjδω+ωj.
(7)

For simplicity, we assume no correlations between the fluctuations of different sites, Cmnω=δmnCnω. The overall system–bath coupling strengths are described by the reorganization energies

λn=1π0Cnωωdω.
(8)

For convenience, we also define the reorganization energy operator, Λ^=nNλn|nn|. Calculations of the reduced equilibrium density operator require the knowledge of the imaginary time bath correlation function

Knτ2τ1=TrBF^nIBτ2F^nIBτ1ρ^Beq,
(9)

with

F^nIBτ=eτĤBF^neτĤB
(10)

denoting the operator F^m in the imaginary time interaction picture with respect to ĤB and ρ^Beq being the equilibrium bath density operator. TrB denotes the trace over the bath degrees of freedom.

The total equilibrium density operator of the whole system-plus-bath supersystem should correspond to the Boltzmann distribution

ρ^Teq=1ZeβĤ,
(11)

with β=kBT1 and Z=TreβĤ being the partition function. The reduced density operator of the system only is obtained by taking the trace over the bath degrees of freedom,

ρ^eq=1ZTrBeβĤ.
(12)

Unfortunately, calculations of this quantity are very difficult. In this letter, however, we derive a simple analytical expression for ρ^eq, requiring only the assumption of high temperature and, thus, keeping correction terms up to the order of β2.

This is accomplished as follows. First, we rewrite the relevant exponential operator as

eβĤ=eβĤε+ĤSB+ĤBexp+0βdτĤJI0τ.
(13)

Here, exp+ denotes the time-ordered exponential18 and I0 denotes the imaginary time interaction picture with respect to Ĥ0 = Ĥε + ĤSB + ĤB. The equilibrium reduced density operator can, therefore, be written as

ρ^eq=ZSZBZeβĤεZSTrB1ZBeβĤSB+ĤBexp+0βdτĤJI0τ.
(14)

The trace is calculated by expanding the time ordered exponential and keeping only the terms of the order β2. After lengthy but straightforward derivations, we arrive at (see the supplementary material for details)

ρ^eqZSZBZeβĤεZSÎS+βΛ^βĤJ12β2Ĥε,ĤJ+12β2ĤJ2+12β2Λ^213β2ĤJΛ^13β2Λ^ĤJ.
(15)

Here, ÎS is the identity operator in the system space. The key step is to notice that this expansion can be resummed into a double exponential function. Thus, we finally obtain

ρ^eq=1ZeffeβĤeff,
(16)

with an effective Hamiltonian

Ĥeff=ĤεΛ^+e16βΛ^ĤJe16βΛ^.
(17)

The normalization factor is defined as Zeff=TreβĤeff. Resummation techniques have been recently applied for a variety of problems in the theory of open quantum systems.19–21 

Equations (16) and (17) constitute the main result of this work. They imply that the equilibrium distribution is achieved according to energies εn − λn and couplings

Jmneff=Je16βλm+λn.
(18)

This result provides justification for our previously obtained expression for Jeff [Ref. 2 and Eq. (1)], with coefficient 16 being very close the numerically obtained value of 0.15.

Our expression for Jeff paints the following physical picture. With increasing system–bath coupling strength λ, excitation sinks deeper in the potential well corresponding to each site; thus, interaction between the sites becomes irrelevant. On the other hand, with increasing temperature, excitation can move more freely on the potential surface and can once again “sense” neighboring sites.

To validate our results, we performed ρ^eq calculations using the numerically exact stochastic path integral (SPI) approach of Ref. 3, which is based on the unraveling of the influence functional. From the obtained expression for the effective Hamiltonian [Eq. (17)], it follows that ρ^eq does not depend on the exact form of the spectral density. To confirm this, we performed calculations using three forms of spectral density, namely, the Debye spectral density

CDebyeω=2λγωω2+γ2,
(19)

where parameter γ is the bath relaxation rate, the Ohmic spectral density

COhmicω=λπωωcexpωωc,
(20)

and the super-Ohmic spectral density

CsuperOhmicω=λπ2ωωc3expωωc,
(21)

where in the latter two cases, parameter ωc is the cut-off frequency.

For simplicity, we performed simulations for a dimer system, with Δε = ε1ε2 = 100 cm−1 and J = 100 cm−1. We also assumed equal reorganization energies for both sites, λ1 = λ2 = λ, and performed calculations for several values of λ. In case of the Debye spectral density, we set γ−1 = 100 fs, and for both Ohmic and super-Ohmic spectral densities, we set ωc = 20 cm−1. These values correspond to the slow bath regime (see below). The temperature for calculations was set to 300 K. Interestingly, calculated equilibrium populations in the site basis remained virtually constant even for λ values up to 1000 cm−1. Evaluations of Eq. (16) predict that in case of equal reorganization energies, site basis equilibrium populations are not affected by the non-canonical effects, and this was supported by our numerical calculations. On the other hand, coherences showed significant differences. This is demonstrated in Fig. 1, where we plot the dependence of ρ12eq values on a dimensionless parameter 2λ/kBT. As this parameter increases, the true value of ρ12eq goes further away from the value corresponding to the Boltzmann distribution. As predicted by our theory, the results do not depend on the form of spectral density in the slow bath regime. The theoretical values agree extremely well with numerical calculations up to 2λ/kBT3 (λ ≲ 300 cm−1 at T = 300 K) and reasonably well up to 2λ/kBT6 (λ ≲ 600 cm−1 at T = 300 K). Application of our theory for larger systems is discussed in the supplementary material.

FIG. 1.

Real part of the equilibrium density operator element ρ12. The solid line denotes results from Eqs. (16) and (17), while points denote exact numerical SPI calculations with different spectral densities. The dashed line corresponds to the expected value based on the Boltzmann distribution for the bare system.

FIG. 1.

Real part of the equilibrium density operator element ρ12. The solid line denotes results from Eqs. (16) and (17), while points denote exact numerical SPI calculations with different spectral densities. The dashed line corresponds to the expected value based on the Boltzmann distribution for the bare system.

Close modal

We have investigated the relationship between the accuracy of our theory and the relaxation time scale of the bath, τrel, which can be well approximated by the inverse of the frequency corresponding to the maximum of the spectral density function. For the Debye spectral density, we have τrel = 1/γ, for the Ohmic spectral density, we have τrel = 1/ωc, and for the super-Ohmic spectral density, we have τrel=1/3ωc. On the basis of extensive calculations (not shown), we conclude that our expression for the equilibrium density operator is accurate for τrel/β ≳ 1. Moreover, numerical investigations show that the slower baths result only in a decrease of factor 1/6 in Eq. (18), but the overall form of the expression for Ĥeff remains the same. This is illustrated in Fig. 2 for the Ohmic spectral density. For slow baths (T = 300 K, ωc = 100 cm−1, resulting in τrel/β ≈ 2), the effective resonance coupling, calculated from the equilibrium density operator, follows the prediction of Eq. (2). If the bath is made significantly faster (ωc = 800 cm−1, resulting in τrel/β ≈ 0.3), we have the same exponential relationship, but with different coefficient instead of 1/6. If we significantly increase the temperature for this fast bath (T = 800 K, resulting in τrel/β ≈ 0.7), the dependence of Jeff again comes very close to that predicted by Eq. (2) and would agree completely for even higher temperatures. Therefore, we can conclude that fast, or Markovian, baths do not induce corrections to the equilibrated state of an open quantum system, while the effects are clearly visible for slow baths. This observation is in line with attempts to relate Markovianity with system–bath entanglement.22 

FIG. 2.

Effective resonance coupling Jeff. The solid line denotes results from Eqs. (16) and (17), while points denote Jeff values extracted from exact numerical SPI calculations of ρ^eq (see the supplementary material for details) with different parameters for the Ohmic spectral density.

FIG. 2.

Effective resonance coupling Jeff. The solid line denotes results from Eqs. (16) and (17), while points denote Jeff values extracted from exact numerical SPI calculations of ρ^eq (see the supplementary material for details) with different parameters for the Ohmic spectral density.

Close modal

Despite the simplicity of our final expressions, in the slow bath regime, they have considerably larger range of validity than the results based on the second-order perturbation theory23 with respect to system–bath coupling (the relevant expressions are listed in the supplementary material), even though the latter is more expensive numerically. This is illustrated in Fig. 3, where it can be seen that the perturbation theory approach agrees well with numerically exact results only for 2λ/kBT2, whereas our theory is considerably more accurate throughout the demonstrated parameter regime. The accuracy of our approach can be traced to the resumation of the perturbative series (see the supplementary material).

FIG. 3.

Real part of the equilibrium density operator element ρ12. Red squares denote exact SPI numerical calculation with Debye spectral density with parameters being the same as in Fig. 1. The black solid line denotes results from Eqs. (16) and (17). The solid light blue line denotes results from the second order perturbation theory (PT) with respect to system–bath coupling. The dashed line corresponds to the expected value based on the Boltzmann distribution.

FIG. 3.

Real part of the equilibrium density operator element ρ12. Red squares denote exact SPI numerical calculation with Debye spectral density with parameters being the same as in Fig. 1. The black solid line denotes results from Eqs. (16) and (17). The solid light blue line denotes results from the second order perturbation theory (PT) with respect to system–bath coupling. The dashed line corresponds to the expected value based on the Boltzmann distribution.

Close modal

Our results imply that the eigenbasis of the bare system Hamiltonian is no longer the eigenbasis of the equilibrated system. The true eigenbasis changes in time due to the influence of bath degrees of freedom. This process is equivalent to the basis rotation picture used in Ref. 4. Since for molecular systems, the eigenbasis of the bare system is denoted as the exciton basis, following our previous works,2,24,25 we denote the true eigenbasis as the excitonic polaron basis. To quantify the change of basis numerically, we have calculated the amount of mixing between the two bases, θ=p1|e12, where |e1⟩ denotes the lowest energy excitonic state and |p1⟩ denotes the lowest energy eigenstate. The dependence of mixing has been calculated for the previously defined dimer system at 300 K and is shown in Fig. 4. Clearly, the polaronic effects are subtle, but they tend to increase when the dimensionless parameter 2λ/kBT increases and that might be enough to break some specific symmetries in the system, as demonstrated in Ref. 2. Using the theoretical effective Hamiltonian of Eq. (17) slightly overestimates these effects. As suggested in Ref. 25, the change, or rotation, of the relevant basis during the equilibration might lead to fast transfer processes, that can be observed spectroscopically. The polaronic effects can be expected to become more significant in case the reorganization energies of system states differ substantially, due to the ĤεĤεΛ^ change in the effective Hamiltonian of Eq. (17).

FIG. 4.

Mixing of excitonic and excitonic polaron states θ=p1|e12. The solid line denotes calculations using the theoretical Hamiltonian of Eq. (17), while points denote exact numerical calculations using the equilibrium density operator calculated using the approach of Ref. 3.

FIG. 4.

Mixing of excitonic and excitonic polaron states θ=p1|e12. The solid line denotes calculations using the theoretical Hamiltonian of Eq. (17), while points denote exact numerical calculations using the equilibrium density operator calculated using the approach of Ref. 3.

Close modal

The theoretical approach developed here has straightforward practical implications. The obtained analytical solution for the equilibrium density operator could be used as a benchmark for approximate theories, for both long time limit of real time dynamics or for imaginary time equilibrium solutions, absolving the need to have exact numerical results. This type of strategy for testing accuracy of theoretical approaches for open quantum systems has been employed previously.7,26,27 Our approach is also a very inexpensive way to calculate the reduced equilibrium density operator that is a prerequisite for simulations of emission spectra by the hybrid cumulant expansion method proposed in Ref. 9.

At this point, we would like to contrast our result for effective coupling with that coming from the small polaron transformation.28–33 Indeed, when applying the small polaron transformation (also in the limiting case of the variational polaron transformation), one obtains that the coupling becomes dressed with bath variables and the thermal averaging results in bath-renormalized coupling

Jpolaron12=J12exp0dωcothβω2C1ω+C2ω2πω2.
(22)

Even though the origins of Jpolaron and Jeff are quite different, quite often Jpolaron is taken as describing the equilibrium state of the system. For this reason, it is worthwhile to compare it with our theory. There are two essential differences between the expression for Jpolaron and our result of Eq. (18). First, the polaron result shows that the bath-renormalized coupling becomes zero for spectral densities that exhibit behavior Cωω for small frequencies due to the divergence of the integral in the exponential. This suggests a sharp transition in the dynamical and equilibrium behavior of the system between cases of Ohmic and super-Ohmic spectral densities. On the other hand, our results show that the equilibrium state is completely insensitive to the exact form of the spectral density if the bath can be assumed to be slow. Second, Eq. (22) gives completely different temperature behavior. At high temperatures, cothβω/22/βω; thus, Jpolaron decreases at higher temperature. As follows from Eq. (18), however, Jeff decreases at lower temperature. We suggest that in some cases, Jpolaron might be misrepresented as being related to the equilibrium properties of the system when, instead, it is just a parameter that is useful to decrease the perturbation term. Indeed, polaron transformation leads to an accurate equilibrium state for some parameter regimes.4 

In conclusion, we have derived an accurate expression for the equilibrium reduced density operator of an open quantum system that linearly interacts with a harmonic oscillator bath at high temperatures. The equilibrium state can be described in terms of an effective Hamiltonian. Our results, however, should hold also for cases when the system–bath interaction is non-linear, provided that the bath is slow and the interaction term can be written as ĤSB=nF^n|nn| with TrBF^n=0, because the assumption of linearity was not used in our derivation. By comparison with numerically exact results, we have demonstrated that our theory is accurate for slow baths and up to intermediate values of reorganization energies. In these conditions, the equilibrium state was shown to be the same regardless of the exact form of spectral density. The effective coupling between the quantum states has contrary temperature dependence to what could be expected from the polaron transformation. We have also shown that for slow baths, our result has a significantly larger range of validity than the usual second order perturbation theory with respect to the system–bath coupling. We expect that our insights will prove relevant for analysis of equilibrium fluorescence spectra or non-equilibrium steady states, characterization of the system–bath entanglement and Markovianity, and interpretation of dynamical spectroscopy results.

See the supplementary material for details of the derivation of the expression for the reduced density operator, analysis of the resummation approach, expression of Jeff from ρ^eq, and summary of the second order perturbation theory for the equilibrium reduced density operator.

Computations were performed using the resources of the High Performance Computing Center “HPC Sauletekis” at the Faculty of Physics, Vilnius University.

1.
E.
Wigner
, “
On the quantum correction for thermodynamic equilibrium
,”
Phys. Rev.
40
,
749
759
(
1932
).
2.
A.
Gelzinis
,
D.
Abramavicius
, and
L.
Valkunas
, “
Non-Markovian effects in time-resolved fluorescence spectrum of molecular aggregates: Tracing polaron formation
,”
Phys. Rev. B
84
,
245430
(
2011
).
3.
J. M.
Moix
,
Y.
Zhao
, and
J.
Cao
, “
Equilibrium-reduced density matrix formulation: Influence of noise, disorder, and temperature on localization in excitonic systems
,”
Phys. Rev. B
85
,
115412
(
2012
).
4.
C.
Kong Lee
,
J.
Cao
, and
J.
Gong
, “
Noncanonical statistics of a spin-boson model: Theory and exact Monte Carlo simulations
,”
Phys. Rev. E
86
,
021109
(
2012
).
5.
J.
Thingna
,
J.-S.
Wang
, and
H.
Peter
, “
Generalized Gibbs state with modified Redfield solution: Exact agreement up to second order
,”
J. Chem. Phys.
136
,
194110
(
2012
).
6.
Y.
Tanimura
, “
Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities
,”
J. Chem. Phys.
141
,
044114
(
2014
).
7.
J.
Iles-Smith
,
N.
Lambert
, and
A.
Nazir
, “
Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems
,”
Phys. Rev. A
90
,
032114
(
2014
).
8.
P.
Kumar
and
S.
Jang
, “
Emission lineshapes of the B850 band of light-harvesting 2 (LH2) complex in purple bacteria: A second order time-nonlocal quantum master equation approach
,”
J. Chem. Phys.
138
,
135101
(
2013
).
9.
J.
Ma
,
J.
Moix
, and
J.
Cao
, “
Förster resonance energy transfer, absorption and emission spectra in multichromophoric systems. II. Hybrid cumulant expansion
,”
J. Chem. Phys.
142
,
094107
(
2015
).
10.
S.
Jang
,
M. D.
Newton
, and
R. J.
Silbey
, “
Multichromophoric Förster resonance energy transfer
,”
Phys. Rev. Lett.
92
,
218301
(
2004
).
11.
L.
Banchi
,
G.
Costagliola
,
A.
Ishizaki
, and
P.
Giorda
, “
An analytical continuation approach for evaluating emission lineshapes of molecular aggregates and the adequacy of multichromophoric Förster theory
,”
J. Chem. Phys.
138
,
184107
(
2013
).
12.
J.
Ma
and
J.
Cao
, “
Förster resonance energy transfer, absorption and emission spectra in multichromophoric systems. I. Full cumulant expansions and system-bath entanglement
,”
J. Chem. Phys.
142
,
094106
(
2015
).
13.
J.
Thingna
,
J.-S.
Wang
, and
H.
Peter
, “
Reduced density matrix for nonequilibrium steady states: A modified Redfield solution approach
,”
Phys. Rev. E
88
,
052127
(
2013
).
14.
H. T.
Quan
,
P.
Zhang
, and
C. P.
Sun
, “
Quantum heat engine with multilevel quantum systems
,”
Phys. Rev. E
72
,
056110
(
2005
).
15.
W. H.
Zurek
, “
Decoherence, einselection, and the quantum origins of the classical
,”
Rev. Mod. Phys.
75
,
715
775
(
2003
).
16.
M.
Perarnau-Llobet
,
H.
Wilming
,
A.
Riera
,
R.
Gallego
, and
J.
Eisert
, “
Strong coupling corrections in quantum thermodynamics
,”
Phys. Rev. Lett.
120
,
120602
(
2018
).
17.
Y.
Subaşi
,
C. H.
Fleming
,
J. M.
Taylor
, and
B. L.
Hu
, “
Equilibrium states of open quantum systems in the strong coupling regime
,”
Phys. Rev. E
86
,
061132
(
2012
).
18.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
New York
,
1995
).
19.
Z.
Gong
,
Z.
Tang
,
S.
Mukamel
,
J.
Cao
, and
J.
Wu
, “
A continued fraction resummation form of bath relaxation effect in the spin-boson model
,”
J. Chem. Phys.
142
,
084103
(
2015
).
20.
T. C.
Berkelbach
,
H.-T.
Chen
, and
D. R.
Reichman
, “
On the accuracy of the Padé-resummed master equation approach to dissipative quantum dynamics
,”
J. Chem. Phys.
144
,
154106
(
2016
).
21.
A.
Erpenbeck
,
C.
Hertlein
,
C.
Schinabeck
, and
M.
Thoss
, “
Extending the hierarchical quantum master equation approach to low temperatures and realistic band structures
,”
J. Chem. Phys.
149
,
064106
(
2018
).
22.
D.
Chruściński
and
S.
Maniscalco
, “
Degree of non-Markovianity of quantum evolution
,”
Phys. Rev. Lett.
112
,
120404
(
2014
).
23.
E.
Geva
,
E.
Rosenman
, and
D.
Tannor
, “
On the second-order corrections to the quantum canonical equilibrium density matrix
,”
J. Chem. Phys.
113
,
1380
1390
(
2000
).
24.
V.
Chorošajev
,
A.
Gelzinis
,
L.
Valkunas
, and
D.
Abramavicius
, “
Dynamics of exciton-polaron transition in molecular assemblies: The variational approach
,”
J. Chem. Phys.
140
,
244108
(
2014
).
25.
J.
Pan
,
A.
Gelzinis
,
V.
Chorošajev
,
M.
Vengris
,
S. S.
Senlik
,
J.-R.
Shen
,
L.
Valkunas
,
D.
Abramavicius
, and
J. P.
Ogilvie
, “
Ultrafast energy transfer within the photosystem II core complex
,”
Phys. Chem. Chem. Phys.
19
,
15356
15367
(
2017
).
26.
C.
Kong Lee
,
J.
Moix
, and
J.
Cao
, “
Accuracy of second order perturbation theory in the polaron and variational polaron frames
,”
J. Chem. Phys.
136
,
204120
(
2012
).
27.
N.
Bellonzi
,
A.
Jain
, and
J. E.
Subotnik
, “
An assessment of mean-field mixed semiclassical approaches: Equilibrium populations and algorithm stability
,”
J. Chem. Phys.
144
,
154110
(
2016
).
28.
R.
Silbey
and
R. A.
Harris
, “
Variational calculation of the dynamics of a two level system interacting with a bath
,”
J. Chem. Phys.
80
,
2615
2617
(
1984
).
29.
S.
Jang
,
Y.-C.
Cheng
,
D. R.
Reichman
, and
J. D.
Eaves
, “
Theory of coherent resonance energy transfer
,”
J. Chem. Phys.
129
,
101104
(
2008
).
30.
S.
Jang
, “
Theory of multichromophoric coherent resonance energy transfer: A polaronic quantum master equation approach
,”
J. Chem. Phys.
135
,
034105
(
2011
).
31.
A.
Kolli
,
A.
Nazir
, and
A.
Olaya-Castro
, “
Electronic excitation dynamics in multichromophoric systems described via a polaronrepresentation master equation
,”
J. Chem. Phys.
135
,
154112
(
2011
).
32.
F. A.
Pollock
,
D. P. S.
McCutcheon
,
B. W.
Lovett
,
E. M.
Gauger
, and
A.
Nazir
, “
A multi-site variational master equation approach to dissipative energy transfer
,”
New J. Phys.
15
,
075018
(
2013
).
33.
M.
Qin
,
H. Z.
Shen
,
X. L.
Zhao
, and
X. X.
Yi
, “
Effects of system-bath coupling on a photosynthetic heat engine: A polaron master-equation approach
,”
Phys. Rev. A
96
,
012125
(
2017
).

Supplementary Material