# Partially polaron-transformed quantum master equation for exciton and charge
transport dynamics

Seogjoo J. Jang<EMAIL_ADDRESS>Department of Chemistry and
Biochemistry, Queens College, City University of New York, 65-30 Kissena
Boulevard, Queens, New York 11367 & PhD Programs in Chemistry and in Physics,
Graduate Center of the City University of New York  Korea Institute for
Advanced Study, Seoul 02455, South Korea

(Submitted to the Journal of Chemical Physics on March 5, 2022)

###### Abstract

Polaron-transformed quantum master equation (PQME) offers a unified framework
to describe the dynamics of quantum systems in both limits of weak and strong
couplings to environmental degrees of freedom. Thus, PQME serves as an
efficient method to describe charge and exciton transfer/transport dynamics
for a broad range of parameters in condensed or complex environments. However,
in some cases, the polaron transformation (PT) being employed in the
formulation invokes an over-relaxation of slow modes and results in premature
suppression of important coherence terms. This work develops a formal
framework to address this issue by employing a partial PT that has smaller
weights for low frequency bath modes. It is shown here that a closed form
expression of a 2nd order time-local PQME including all the inhomogeneous
terms can be derived for a general form of partial PT, although more
complicated than that for the full PT. All the expressions needed for
numerical calculation are derived in detail. Applications to a model with two
system states demonstrate the feasibility and utility of the present approach.

##  I Introduction

Polaron transformation (PT)Landau and Pekar (2008); Fröhlich (1954); Holstein
(1959a, b, 1978); Rackovsky and Silbey (1973); Jackson and Silbey (1983);
Silbey and Harris (1984); Harris and Silbey (1985a, b); Nitzan (2006); Cheng
and Silbey (2008); Jang  _et al._ (2008); Jang (2009, 2011); Nazir (2009);
McCutcheon and Nazir (2011); Yang, Devi, and Jang (2012); Pouthier (2013);
Chen  _et al._ (2011); Zhao  _et al._ (2012); Chorosajev  _et al._ (2014);
Hamm and Tsironis (2008); Lee, Moix, and Cao (2015); Xu and Cao (2016); Wang
and Zhao (2020); Balzer  _et al._ (2021) has served as both an important
conceptual framework and an efficient computational tool for describing charge
and exciton transfer/transport dynamics in various condensed and complex
media. PT creates a polaron picture where molecular vibrations and phonon
modes of environments, collectively referred to as bath here, relax
differently with respect to different site localized states. When such
responses to localized states are fast, the states “dressed” with those
responses, so called polarons, serve as effective means to represent the
contribution of the bath because they have already taken some contributions of
the bath responses into consideration up to an infinite order. Alternatively,
PT can simply be viewed as a useful unitary transformation in the combined
space of system and bath, which produces a new renormalized interaction
Hamiltonian term that remains small even in the limit of strong couplings to
the bath and can be used for the derivation of a new kind of quantum master
equation (QME). Indeed, utilizing this fact has led to a general approach
called polaron-transformed QME (PQME).Jang  _et al._ (2008); Jang (2009,
2011); Nazir (2009); McCutcheon and Nazir (2011); Yang, Devi, and Jang (2012)

The merit of PQME is its efficiency in accounting for both weak and strong
couplings to the bath, while offering reasonable interpolation between the two
limits. Earlier versions of PQMEsJang  _et al._ (2008); Jang (2009, 2011);
Nazir (2009) employed a full PT where all the bath degrees of freedom, modeled
as harmonic oscillators, are fully relaxed to site local system states.
However, in some cases, this is not appropriate because it invokes an over-
relaxation of some slow modes, which can cause premature suppression of
important coherence terms within the second order approximation. Variational
PQMEMcCutcheon and Nazir (2011); Pollock  _et al._ (2013) and, more recently,
frozen mode PQMETeh, Jin, and Cheng (2019) approaches have been developed to
address this issue. This work develops a formulation that can put these
worksMcCutcheon and Nazir (2011); Pollock  _et al._ (2013); Teh, Jin, and
Cheng (2019) into a broader context and address the issue of over-relaxation
by extending the PQME for partial PTs of fairly general kind, where high
frequency bath modes are transformed preferentially. It is shown here that a
general closed form expression of a 2nd order time-local PQME including all
the inhomogeneous terms can be derived as in the case of the full PT,Jang
(2009) although the resulting expressions are more complicated. All the terms
needed for the calculation of a second order partially polaron-transformed QME
(p-PQME) are derived in detail. Then, numerical results are provided for a
model with two system states, which demonstrate the feasibility and utility of
the p-PQME approach.

The paper is organized as follows. Section II provides the model Hamiltonian
and presents derivation of all terms involved in the second order time-local
p-PQME. Section III considers a case of two-state system and presents results
of model calculations. Section IV summarizes main results of this paper and
offers concluding remarks.

##  II Theory

###  II.1 Hamiltonian and partial polaron transformation

Let us consider a quantum system consisting of $N$ coupled quantum states,
with each representing a localized electronic excitation or charge carrying
state. The Hamiltonian of the quantum system consisting of these states in
general can be expressed as

$H_{s}=\sum_{j=1}^{N}E_{j}|j\rangle\langle
j|+\sum_{j,k=1}^{N}J_{jk}|j\rangle\langle k|,$ (1)

where $E_{j}$ is the energy of the site local state $|j\rangle$, and $J_{jk}$,
for $j\neq k$, is the electronic coupling between states $|j\rangle$ and
$|k\rangle$ that is assumed to be a real number here. The double summation in
the above equation also includes the case $j=k$, for which $J_{jk}$ is assumed
to be zero.

All other degrees of freedom constituting the total Hamiltonian are referred
to here as bath, which includes all the vibrational modes of molecules and the
polarization response of environmental degrees of freedom. Some of these may
exhibit significant anharmonic character in real molecular systems. However,
for the sake of simplicity, it is assumed here that all of them can be modeled
as coupled harmonic oscillators. In addition, all couplings of the bath to the
system Hamiltonian are assumed to be diagonal with respect to the site local
states $|j\rangle$’s and to be linear in the displacements of the bath
oscillators. Thus, the total Hamiltonian is assumed to be of the following
standard form:

$H=H_{s}+H_{b}+H_{sb},$ (2)

where

$\displaystyle
H_{b}=\sum_{n}\hbar\omega_{n}(b_{n}^{\dagger}b_{n}+\frac{1}{2}),$ (3)
$\displaystyle
H_{sb}=\sum_{j=1}^{N}\sum_{n}\hbar\omega_{n}g_{n,j}(b_{n}+b_{n}^{\dagger})|j\rangle\langle
j|.$ (4)

In the above expressions, $\omega_{n}$ and $b_{n}$ ($b_{n}^{\dagger}$) are the
frequency and the lowering (raising) operator of each normal mode constituting
the bath degrees of freedom, and $g_{n,j}$ represents the (dimensionless)
coupling strength of each mode to the state $|j\rangle$. The total
Hamiltonian, Eq. (2), is a straightforward generalization of the spin-boson
HamiltonianCaldeira and Leggett (1983); Weiss (1993) and has been used widely
for molecular excitonsKenkre and Reineker (1982); May and Kühn (2011); Jang
(2020) and charge transport dynamics.Weiss (1993); May and Kühn (2011);
Coropceanu  _et al._ (2007)

The full information on the combined system and bath degrees of freedom in
general requires determining the total density operator $\rho(t)$, which is
governed by the following quantum Liouville equation:

$\frac{d}{dt}\rho(t)=-i{\mathcal{L}}\rho(t)\equiv-\frac{i}{\hbar}[H,\rho(t)].$
(5)

Due to the large number of bath degrees of freedom typically involved, solving
the above equation exactly is difficult even for simple forms of $H_{b}$ and
$H_{sb}$ given by Eqs. (3) and (4). In the quantum master equation approach,
only the reduced system degrees of freedom is solved explicitly whereas the
effects of the bath are treated only to certain extents that are necessary.
Typically, these effects can be encoded entirely into appropriate bath
spectral densities.

For the case of Ohmic or sub-Ohmic bath spectral density, the sum of Huang-
Rhys factors for the bath degrees of freedom diverges.Weiss (1993) This
implies vanishing Debye-Waller factors, which result in premature suppression
of some coherence terms when the full PT is applied first. This issue is
significant in particular for charge transfer processes, where the bath
spectral density is typically known to contain Ohmic or even sub-Ohmic low
frequency components. A simple way to avoid such suppression is to limit the
PT to only fast enough bath modes. To this end, let us introduce a weighting
function $W_{h}(\omega)$ with the following limiting behavior:

$W_{h}(\omega)=\left\\{\begin{array}[]{cl}O(\omega^{\alpha}),&\mbox{, for
}\omega\rightarrow 0\\\ 1&\mbox{, for
}\omega\rightarrow\infty\end{array}\right.$ (6)

The scaling behavior for $\omega\rightarrow 0$ in the above equation, for
sufficiently large value of $\alpha$, can suppress the sluggish components of
the bath spectral density and thus prevents the corresponding Debye-Waller
factor from becoming zero. For the case of Ohmic bath, it turns out that any
value $\alpha>1$ results in well-behaving equation. In case the bath spectral
density has sub-Ohmic components, the lower bound for $\alpha$ should increase
accordingly.

Let us now define a generating function of a partial PT as follows:

$G=\sum_{j=1}^{N}\sum_{n}g_{n,j}W_{h}(\omega_{n})(b_{n}^{\dagger}-b_{n})|j\rangle\langle
j|$ (7)

The corresponding PT, when applied to the system Hamiltonian, Eq. (1), results
in

$\displaystyle\tilde{H}_{s}$ $\displaystyle=$ $\displaystyle e^{G}H_{s}e^{-G}$
(8) $\displaystyle=$ $\displaystyle\sum_{j=1}^{N}E_{j}|j\rangle\langle
j|+\sum_{j,k=1}^{N}J_{jk}\theta_{j}^{\dagger}\theta_{k}|j\rangle\langle k|,$

where

$\theta_{j}=e^{-\sum_{n}g_{n,j}W_{h}(\omega_{n})(b_{n}^{\dagger}-b_{n})}.$ (9)

On the other hand, it is straightforward to show that

$\displaystyle e^{G}(H_{b}+H_{sb})e^{-G}$
$\displaystyle=H_{b}+\sum_{j=1}^{N}\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))(b_{n}+b_{n}^{\dagger})|j\rangle\langle
j|$
$\displaystyle-\sum_{j=1}^{N}\sum_{n}\hbar\omega_{n}g_{n,j}^{2}W_{h}(\omega_{n})(2-W_{h}(\omega_{n}))|j\rangle\langle
j|.$ (10)

Combining the above expression with Eq. (8), one can obtain the following
expression for the transformed total Hamiltonian:

$\tilde{H}=e^{G}He^{-G}=\tilde{H}_{s}^{p}+\tilde{H}_{sb}+H_{b},$ (11)

where

$\tilde{H}_{s}^{p}=\sum_{j=1}^{N}\tilde{E}_{j}|j\rangle\langle j|,\\\ $ (12)

with a partially renormalized energy for site $j$ given by

$\tilde{E}_{j}=E_{j}-\sum_{n}\hbar\omega_{n}g_{n,j}^{2}W_{h}(\omega_{n})(2-W_{h}(\omega_{n})).$
(13)

The second term in Eq. (11), $\tilde{H}_{sb}$, is a partially renormalized
system-bath interaction Hamiltonian given by

$\displaystyle\tilde{H}_{sb}=\sum_{j,k=1}^{N}J_{jk}\theta_{j}^{\dagger}\theta_{k}|j\rangle\langle
k|\ $
$\displaystyle+\sum_{j=1}^{N}\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))(b_{n}+b_{n}^{\dagger})|j\rangle\langle
j|,$ (14)

where $\theta_{j}^{\dagger}\theta_{k}$ can be expressed as

$\theta_{j}^{\dagger}\theta_{k}=e^{\sum_{n}\delta
g_{n,jk}W_{h}(\omega_{n})(b_{n}^{\dagger}-b_{n})}\ ,$ (15)

with $\delta g_{n,jk}=(g_{n,j}-g_{n,k})$.

###  II.2 Partially polaron transformed quantum master equation (p-PQME)

A complete derivation of a QME defined in the partial polaron transformation
(p-PT) space, as defined in the previous subsection, is provided here.

####  II.2.1 Quantum Liouville equation in the polaron and interaction picture

The total Hamiltonian in the partial polaron picture, Eq. (11), can be divided
into new effective zeroth and first order terms as follows:

$\tilde{H}=\tilde{H}_{0}+\tilde{H}_{1},$ (16)

where

$\displaystyle\tilde{H}_{0}=\tilde{H}_{s}^{p}+\langle\tilde{H}_{sb}\rangle_{b}+H_{b}\equiv\tilde{H}_{0,s}+H_{b},$
(17)
$\displaystyle\tilde{H}_{1}=\tilde{H}_{sb}-\langle\tilde{H}_{sb}\rangle_{b}.$
(18)

In above expressions,
$\langle\tilde{H}_{sb}\rangle_{b}=Tr_{b}\\{\tilde{H}_{sb}\rho_{b}\\}$ with
$\rho_{b}=e^{-\beta H_{b}}/Tr_{b}\\{e^{-\beta H_{b}}\\}$. This term represents
the average system-bath interaction in the partial polaron picture.

In Eq. (17), the system part of the new zeroth order Hamiltonian,
$\tilde{H}_{0,s}=\tilde{H}^{p}_{s}+\langle\tilde{H}^{c}_{s}\rangle_{b}$,
includes the average effect of system-bath interactions in the partial polaron
picture, and can be expressed as

$\tilde{H}_{0,s}=\sum_{j=1}^{N}\tilde{E}_{j}|j\rangle\langle
j|+\sum_{j,k=1}^{N}\tilde{J}_{jk}|j\rangle\langle k|,$ (19)

where $\tilde{E}_{j}$ has been defined by Eq. (13) and
$\tilde{J}_{jk}=w_{jk}J_{jk}$ with

$\displaystyle w_{jk}$ $\displaystyle=$
$\displaystyle\langle\theta_{j}^{\dagger}\theta_{k}\rangle=\langle\theta_{k}^{\dagger}\theta_{j}\rangle$
(20) $\displaystyle=$ $\displaystyle
e^{-\sum_{n}\coth(\beta\hbar\omega_{n}/2)\delta
g_{n,jk}^{2}W_{h}(\omega_{n})^{2}/2}.$

Unlike the case of the full PT, the Debye-Waller factor $w_{jk}$ given above
is non-zero even for the Ohmic bath spectral density because of the weighting
function that satisfies Eq. (6).

Similarly, the first order term $\tilde{H}_{1}$, Eq. (18), can be expressed as

$\tilde{H}_{1}=\sum_{j,k=1}^{N}\tilde{B}_{jk}|j\rangle\langle k|,$ (21)

where

$\tilde{B}_{jk}=J_{jk}(\theta_{j}^{\dagger}\theta_{k}-w_{jk})+\delta_{jk}D_{j}.$
(22)

In the above expression, $\delta_{jk}$ is the Kronecker-delta symbol and

$D_{j}=\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))(b_{n}+b_{n}^{\dagger}).$
(23)

Thus, the bath operator $\tilde{B}_{jk}$ given by Eq. (22) is a sum of the
renormalized system-bath interaction term (relative to its average) due to
partial PT (for $j\neq k$) and of the remaining linear interaction term (for
$j=k$) for portions of bath modes that have not been transformed.

Having defined $\tilde{H}_{0}$, for which exact time evolution can be
implemented numerically, let us now consider the dynamics in the interaction
picture of $\tilde{H}_{0}$. First, $\tilde{H}_{1}$ in this interaction picture
becomes

$\displaystyle\tilde{H}_{1,I}(t)\equiv
e^{i\tilde{H}_{0}t/\hbar}\tilde{H}_{1}e^{-i\tilde{H}_{0}t/\hbar}$
$\displaystyle=\sum_{j,k=1}^{N}\tilde{B}_{jk}(t){\mathcal{T}}_{jk}(t),$ (24)

where

$\displaystyle{\mathcal{T}}_{jk}(t)=e^{i\tilde{H}_{0,s}t/\hbar}|j\rangle\langle
k|e^{-i\tilde{H}_{0,s}t/\hbar},$ (25)
$\displaystyle\tilde{B}_{jk}(t)=e^{iH_{b}t/\hbar}\tilde{B}_{jk}e^{-iH_{b}t/\hbar}$
$\displaystyle=J_{jk}(\theta_{j}^{\dagger}(t)\theta_{k}(t)-w_{jk})+\delta_{jk}D_{j}(t).$
(26)

In the above expression,

$\displaystyle\theta_{j}^{\dagger}(t)\theta_{k}(t)=e^{\sum_{n}\delta
g_{n,jk}W_{h}(\omega_{n})(b_{n}^{\dagger}e^{i\omega_{n}t}-b_{n}e^{-i\omega_{n}t})},$
(27) $\displaystyle
D_{j}(t)=\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))$
$\displaystyle\hskip
72.26999pt\times(b_{n}e^{-i\omega_{n}t}+b_{n}^{\dagger}e^{i\omega_{n}t}).$
(28)

In the interaction picture with respect to $\tilde{H}_{0}$, the partially
polaron-transformed total density operator becomes
$\tilde{\rho}_{I}(t)=e^{i\tilde{\mathcal{L}}_{0}t}\tilde{\rho}(t)$, which
evolves according to the following time dependent quantum Liouville equation:

$\frac{d}{dt}\tilde{\rho}_{I}(t)=-i\tilde{\mathcal{L}}_{1,I}(t)\tilde{\rho}_{I}(t)\equiv-\frac{i}{\hbar}[\tilde{H}_{1,I}(t),\tilde{\rho}_{I}(t)]\
,$ (29)

where the second equality serves as the definition of
$\tilde{\mathcal{L}}_{1,I}(t)$.

####  II.2.2 Quantum master equation for reduced density operator

Taking trace of $\tilde{\rho}_{I}(t)$ over the bath degrees of freedom leads
to the following interaction-picture reduced system density operator defined
in the p-PT system-bath space:

$\tilde{\sigma}_{I}(t)\equiv Tr_{b}\left\\{\tilde{\rho}_{I}(t)\right\\}.$ (30)

While the above reduced density operator still retains full information on the
system degrees of freedom in the p-PT space, it is important to note that the
trace operation makes it impossible to retrieve the full information on the
system prior to the application of p-PT. On the other hand, properties
diagonal in the site basis, which are not affected by p-PT, remain intact.

A formally exact time evolution equation with time-convolution can be obtained
for $\tilde{\sigma}_{I}(t)$ employing a standard projection operator
techniqueJang, Cao, and Silbey (2002); Jang (2020) and a projection operator
${\mathcal{P}}(\cdot)\equiv\rho_{b}Tr_{b}\\{(\cdot)\\}$ as follows:

$\displaystyle\frac{d}{dt}\tilde{\sigma}_{I}(t)=-\int_{0}^{t}d\tau
Tr_{b}\Big{\\{}\tilde{\mathcal{L}}_{1,I}(t)e_{(+)}^{-i\int_{\tau}^{t}d\tau^{\prime}{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau^{\prime})}$
$\displaystyle\hskip
108.405pt\times{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau)\rho_{b}\Big{\\}}\tilde{\sigma}_{I}(\tau)$
$\displaystyle\hskip 36.135pt-
iTr_{b}\\{\tilde{\mathcal{L}}_{1,I}(t)e_{(+)}^{-i\int_{0}^{t}d\tau{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau)}{\mathcal{Q}}\tilde{\rho}(0)\\},$
(31)

where ${\mathcal{Q}}=1-{\mathcal{P}}$ and

${\mathcal{Q}}\tilde{\rho}(0)=e^{G}\rho(0)e^{-G}-\rho_{b}Tr_{b}\\{e^{G}\rho(0)e^{-G}\\}.$
(32)

Alternatively, replacing $\tilde{\rho}_{I}(\tau)$ with the back propagation of
$\tilde{\rho}_{I}(t)$ from $t$ to $\tau$ within the projection operator
formalism,Shibata and Arimitsu (1980); Jang (2020) one can obtain the
following formally exact time-local equation:

$\displaystyle\frac{d}{dt}\tilde{\sigma}_{I}(t)=-\int_{0}^{t}d\tau
Tr_{b}\Big{\\{}\tilde{\mathcal{L}}_{1,I}(t)(1+i\Gamma_{1,I}(t))^{-1}$
$\displaystyle\hskip 14.45377pt\times
e_{(+)}^{-i\int_{\tau}^{t}d\tau^{\prime}{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau^{\prime})}{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau){\mathcal{P}}e_{(-)}^{i\int_{\tau}^{t}d\tau^{\prime}\tilde{\mathcal{L}}_{1,I}(\tau^{\prime})}\rho_{b}\Big{\\}}\tilde{\sigma}_{I}(t)$
$\displaystyle\hskip
7.22743pt-i{\mathcal{P}}\tilde{\mathcal{L}}_{1,I}(t)(1+i\Gamma_{1,I}(t))^{-1}e_{(+)}^{-i\int_{0}^{t}d\tau{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau)}{\mathcal{Q}}\tilde{\rho}(0),$
(33)

where

$\Gamma_{1,I}(t)=\int_{0}^{t}d\tau
e_{(+)}^{-i\int_{\tau}^{t}d\tau^{\prime}{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau^{\prime})}{\mathcal{Q}}\tilde{\mathcal{L}}_{1,I}(\tau){\mathcal{P}}e_{(-)}^{i\int_{\tau}^{t}d\tau^{\prime}\tilde{\mathcal{L}}_{1,I}(\tau^{\prime})}.$
(34)

When approximated up to the second order, it is straightforward to show that
Eq. (33) simplifies to

$\frac{d}{dt}\tilde{\sigma}_{I}(t)=-{\mathcal{R}}(t)\tilde{\sigma}_{I}(t)+{\mathcal{I}}(t),$
(35)

where

$\displaystyle{\mathcal{R}}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}d\tau
Tr_{b}\\{\tilde{\mathcal{L}}_{1,I}(t)\tilde{\mathcal{L}}_{1,I}(\tau)\rho_{b}\\},$
(36) $\displaystyle{\mathcal{I}}(t)$ $\displaystyle=$
$\displaystyle{\mathcal{I}}^{(1)}(t)+{\mathcal{I}}^{(2)}(t)$ (37)
$\displaystyle=$ $\displaystyle-
iTr_{b}\\{\tilde{\mathcal{L}}_{1,I}(t){\mathcal{Q}}\tilde{\rho}(0)\\}$
$\displaystyle-\int_{0}^{t}d\tau
Tr_{b}\\{\tilde{\mathcal{L}}_{1,I}(t)\tilde{\mathcal{L}}_{1,I}(\tau){\mathcal{Q}}\tilde{\rho}(0)\\}.$

In the above expression, ${\mathcal{I}}^{(1)}(t)$ and ${\mathcal{I}}^{(2)}(t)$
represent the first and second order inhomogeneous terms of the time evolution
equation. Note that Eq. (35) can also be obtained from Eq. (31), by simply
replacing $\tilde{\sigma}_{I}(\tau)$ with $\tilde{\sigma}_{I}(t)$, which does
not affect the accuracy at the second order level.

Equation (35) is the 2nd order time local p-PQME expressed in Liouville space.
In all previous works and in the present paper, this time local form has been
chosen due to its convenience. However, a time nonlocal 2nd order expressions
can also be derived directly from Eq. (31), and its performance compared to
the time-local form needs to be understood better through actual numerical
studies. Many numerical tests so far seem to indicate that the performance of
the time-local form is better than that of the time non-local form for exciton
and charge transfer dynamics near room temperature. However, a recent workLai
and Geva (2021) provides examples of the opposite case where the performance
of time-nonlocal QME is more satisfactory than the time-local form. Therefore,
further tests and comparative calculations are necessary to make more
comprehensive assessment of the two approaches. With this point clear, the
rest of this section provides detailed expressions for the relaxation
superoperator, Eq. (36), and the inhomogeneous term, Eq. (37), in the Hilbert
space.

####  II.2.3 Homogeneous terms of the 2nd order time local p-PQME

The Hilbert space expression for ${\mathcal{R}}(t)\tilde{\sigma}_{I}(t)$ in
Eq. (35) can be obtained by employing Eqs. (24)-(26) in Eq. (29) and by taking
advantage of the cyclic invariance of the trace operation with respect to the
bath degrees of freedom. The resulting expression is as follows:

$\displaystyle{\mathcal{R}}(t)\tilde{\sigma}_{I}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar^{2}}\sum_{j,k=1}^{N}\sum_{j^{\prime},k^{\prime}=1}^{N}\int_{0}^{t}d\tau\langle\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\rangle$
(38) $\displaystyle\hskip
57.81621pt\times[{\mathcal{T}}_{jk}(t),{\mathcal{T}}_{j^{\prime}k^{\prime}}(\tau)\tilde{\sigma}_{I}(t)]$
$\displaystyle+{\rm H.c.},$

where ${\rm H.c}$ refers to the Hermitian conjugate of all previous terms and
$\langle\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\rangle$ (with
subscript $b$ omitted) represents averaging with respect to the equilibrium
bath density operator, $\rho_{b}$.

Appendix A describes calculation of all the terms constituting
$\langle\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\rangle$. When
the resulting expressions, Eqs. (104), (115), (120), and (LABEL:eq:dd-t), are
used in Eq. (103), it can be expressed as

$\displaystyle\langle\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\rangle=\tilde{J}_{jk}\tilde{J}_{j^{\prime}k^{\prime}}\left(e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t-\tau)}-1\right)$
$\displaystyle\hskip
36.135pt+\delta_{jk}\tilde{J}_{j^{\prime}k^{\prime}}{\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t-\tau)+\delta_{j^{\prime}k^{\prime}}\tilde{J}_{jk}{\mathcal{M}}_{j^{\prime},kj}(t-\tau)$
$\displaystyle\hskip
36.135pt+\delta_{jk}\delta_{j^{\prime}k^{\prime}}{\mathcal{C}}_{jj^{\prime}}(t-\tau),$
(39)

where

$\displaystyle{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)=\sum_{n}\delta
g_{n,jk}\delta g_{n,j^{\prime}k^{\prime}}W_{h}(\omega_{n})^{2}$
$\displaystyle\hskip
21.68121pt\times\left(\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\cos(\omega_{n}t)-i\sin(\omega_{n}t)\right),$
(40)
$\displaystyle{\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)=\sum_{n}\hbar\omega_{n}g_{n,j}\delta
g_{n,j^{\prime}k^{\prime}}(1-W_{h}(\omega_{n}))W_{h}(\omega_{n})$
$\displaystyle\hskip
21.68121pt\times\left(\cos(\omega_{n}t)-i\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\sin(\omega_{n}t)\right),$
(41)
$\displaystyle{\mathcal{C}}_{jj^{\prime}}(t)=\sum_{n}\hbar^{2}\omega_{n}^{2}g_{n,j}g_{n,j^{\prime}}(1-W_{h}(\omega_{n}))^{2}$
$\displaystyle\hskip
21.68121pt\times\left(\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\cos(\omega_{n}t)-i\sin(\omega_{n}t)\right).$
(42)

Similar expressions as above have also been derived in the context of
variational PQME.Pollock  _et al._ (2013) Note that the three bath correlation
functions defined above satisfy the following symmetry properties:

$\displaystyle{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)$ $\displaystyle=$
$\displaystyle{\mathcal{K}}_{kj,k^{\prime}j^{\prime}}(t)$ (43)
$\displaystyle=$
$\displaystyle-{\mathcal{K}}_{jk,k^{\prime}j^{\prime}}(t)=-{\mathcal{K}}_{kj,j^{\prime}k^{\prime}}(t),$
$\displaystyle{\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)$ $\displaystyle=$
$\displaystyle-{\mathcal{M}}_{j,k^{\prime}j^{\prime}}(t),$ (44)
$\displaystyle{\mathcal{C}}_{jj^{\prime}}(t)$ $\displaystyle=$
$\displaystyle{\mathcal{C}}_{j^{\prime}j}(t).$ (45)

For more compact expressions for the three correlation functions, Eqs.
(40)-(42), it is useful to introduce the following bath spectral density:

${\mathcal{J}}_{jj^{\prime}}(\omega)=\pi\hbar\sum_{n}\delta(\omega-\omega_{n})\omega_{n}^{2}g_{n,j}g_{n,j^{\prime}},$
(46)

and also the following auxiliary bath spectral densities:

$\displaystyle{\mathcal{J}}^{(1)}_{j,j^{\prime}k^{\prime}}(\omega)=\pi\hbar\sum_{n}\delta(\omega-\omega_{n})\omega_{n}^{2}g_{n,j}\delta
g_{n,j^{\prime}k^{\prime}}$
$\displaystyle={\mathcal{J}}_{jj^{\prime}}(\omega)-{\mathcal{J}}_{jk^{\prime}}(\omega),$
(47)
$\displaystyle{\mathcal{J}}^{(2)}_{jk,j^{\prime}k^{\prime}}(\omega)=\pi\hbar\sum_{n}\delta(\omega-\omega_{n})\omega_{n}^{2}\delta
g_{n,jk}\delta g_{n,j^{\prime}k^{\prime}}$
$\displaystyle={\mathcal{J}}_{jj^{\prime}}(\omega)+{\mathcal{J}}_{kk^{\prime}}(\omega)-{\mathcal{J}}_{jk^{\prime}}(\omega)-{\mathcal{J}}_{kj^{\prime}}(\omega).$
(48)

Then, it is straightforward to show that Eqs. (40)-(42) can be expressed as

$\displaystyle{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)=\frac{1}{\pi\hbar}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}^{(2)}_{jk,j^{\prime}k^{\prime}}(\omega)}{\omega^{2}}W_{h}(\omega)^{2}$
$\displaystyle\hskip
36.135pt\times\left(\coth\left(\frac{\beta\hbar\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right),$ (49)
$\displaystyle{\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)=\frac{1}{\pi}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}^{(1)}_{j,j^{\prime}k^{\prime}}(\omega)}{\omega}(1-W_{h}(\omega))W_{h}(\omega)$
$\displaystyle\hskip 36.135pt\times\left(\cos(\omega
t)-i\coth\left(\frac{\beta\hbar\omega}{2}\right)\sin(\omega t)\right),$ (50)
$\displaystyle{\mathcal{C}}_{jj^{\prime}}(t)=\frac{\hbar}{\pi}\int_{0}^{\infty}d\omega{\mathcal{J}}_{jj^{\prime}}(\omega)(1-W_{h}(\omega))^{2}$
$\displaystyle\hskip
36.135pt\times\left(\coth\left(\frac{\beta\hbar\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right).$ (51)

Note also that $w_{jk}$ defined by Eq. (20) can be expressed as follows:

$\displaystyle w_{jk}=e^{-{\mathcal{K}}_{jk,jk}(0)/2}$
$\displaystyle=\exp\left\\{-\frac{1}{2\pi\hbar}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}^{(2)}_{jk,jk}(\omega)}{\omega^{2}}W_{h}(\omega)^{2}\coth\left(\frac{\beta\hbar\omega}{2}\right)\right\\}$

If $W_{h}(\omega)=1$, ${\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)$ and
${\mathcal{C}}_{jj^{\prime}}(t)$ become zero and the expressions for
${\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)$ and $w_{jk}$ reduce to those for
the original 2nd order PQME based on the full PT. On the other hand, for
$W_{h}(\omega)=0$, ${\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)$ and
${\mathcal{C}}_{jj^{\prime}}(t)$ become zero and the expression for
${\mathcal{C}}_{jj^{\prime}}(t)$ reduce to that for a conventional 2nd order
time-local QME (without PT). In this sense, the above expressions can be
viewed as general ones incorporating the two limiting cases. As the next
simplest case, let us consider the case where $W_{h}(\omega)$ is a step
function. For this, $W_{h}(\omega)(1-W_{h}(\omega))=0$ and thus
${\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)$ becomes zero. As a result, the
relaxation superoperator becomes a simple sum of those due to p-PT part and
untransformed linear coupling terms.

####  II.2.4 Inhomogeneous terms of p-PQME

For the calculation of inhomogeneous terms, let us assume that the initial
untransformed total density operator is given by

$\rho(0)=\sigma(0)\rho_{b},$ (53)

where

$\sigma(0)=\sum_{j,k=1}^{N}\sigma_{jk}(0)|j\rangle\langle k|.$ (54)

Then,

${\mathcal{Q}}\tilde{\rho}(0)=\sum_{j,k=1}^{N}\sigma_{jk}(0){\mathcal{T}}_{jk}(0)\delta\tilde{\rho}_{b,jk},$
(55)

where Eq. (25) with $t=0$ has been used and

$\delta\tilde{\rho}_{b,jk}=\theta^{\dagger}_{j}\rho_{b}\theta_{k}-w_{jk}\rho_{b},$
(56)

with the convention that $w_{jk}=1$ for $j=k$. Therefore, the first order
inhomogeneous term in Eq. (37) can be expressed as

$\displaystyle{\mathcal{I}}^{(1)}(t)=-iTr_{b}\\{\tilde{\mathcal{L}}_{1,I}(t){\mathcal{Q}}\tilde{\rho}(0)\\}$
$\displaystyle\hskip
36.135pt=-\frac{i}{\hbar}\sum_{j,k=1}^{N}\sum_{j^{\prime},k^{\prime}=1}^{N}Tr_{b}\left\\{\tilde{B}_{jk}(t)\delta\tilde{\rho}_{b,j^{\prime}k^{\prime}}\right\\}$
$\displaystyle\hskip
86.72377pt\times\sigma_{j^{\prime}k^{\prime}}(0)[{\mathcal{T}}_{jk}(t),{\mathcal{T}}_{j^{\prime}k^{\prime}}(0)].$
(57)

In the above expression, the trace over the bath can be calculated explicitly
employing the definitions for $\tilde{B}_{jk}(t)$ and
$\delta\tilde{\rho}_{b,jk}$, Eqs. (26) and (56), respectively. Details of this
calculation are provided in Appendix B, and the resulting expression can be
summed up as

$\displaystyle
Tr_{b}\\{\tilde{B}_{jk}(t)\delta\tilde{\rho}_{b,j^{\prime}k^{\prime}}\\}$
$\displaystyle=w_{j^{\prime}k^{\prime}}\Big{\\{}\tilde{J}_{jk}\left(e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)}f_{jk,k^{\prime}}(t)-1\right)$
$\displaystyle\hskip
36.135pt+\delta_{jk}\left({\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)+h_{j,k^{\prime}}(t)\right)\Big{\\}},$
(58)

where

$\displaystyle
f_{jk,k^{\prime}}(t)=\exp\left\\{2i\sum_{n}g_{n,k^{\prime}}\delta
g_{n,jk}W_{h}(\omega_{n})^{2}\sin(\omega_{n}t)\right\\}$
$\displaystyle=\exp\left\\{\frac{2i}{\pi\hbar}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}^{(1)}_{k^{\prime},jk}(\omega)}{\omega^{2}}W_{h}(\omega)^{2}\sin(\omega
t)\right\\},$ (59) $\displaystyle
h_{j,k^{\prime}}(t)=2\sum_{n}\hbar\omega_{n}g_{n,j}g_{n,k^{\prime}}(1-W_{h}(\omega_{n}))W_{h}(\omega_{n})\cos(\omega_{n}t)$
$\displaystyle=\frac{2}{\pi}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}_{jk^{\prime}}(\omega)}{\omega}(1-W_{h}(\omega))W_{h}(\omega)\cos(\omega
t).$ (60)

Similarly, the second order inhomogeneous term in Eq. (37) can be expressed as

$\displaystyle{\mathcal{I}}^{(2)}(t)$ $\displaystyle=$
$\displaystyle-\int_{0}^{t}d\tau
Tr_{b}\left\\{\tilde{\mathcal{L}}_{1,I}(t)\tilde{\mathcal{L}}_{1,I}(\tau){\mathcal{Q}}\tilde{\rho}(0)\right\\}$
$\displaystyle=$
$\displaystyle-\frac{1}{\hbar^{2}}\sum_{j,k=1}^{N}\sum_{j^{\prime},k^{\prime}=1}^{N}\sum_{j^{\prime\prime},k^{\prime\prime}=1}^{N}$
$\displaystyle\int_{0}^{t}d\tau
Tr_{b}\left\\{\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\delta\tilde{\rho}_{b,j^{\prime\prime}k^{\prime\prime}}\right\\}$
$\displaystyle\times\sigma_{j^{\prime\prime}k^{\prime\prime}}(0)[{\mathcal{T}}_{jk}(t),{\mathcal{T}}_{j^{\prime}k^{\prime}}(\tau){\mathcal{T}}_{j^{\prime\prime}k^{\prime\prime}}(0)]+{\rm
H.c.},$

where H.c. refers to Hermitian conjugates of all previous terms. Detailed
expressions for the trace of bath operators in the above expression are
derived in Appendix C.

As in the case of the relaxation superoperator, the expressions for the
inhomogeneous terms shown above reduce to those for the 2nd order PQME for
$W_{h}(\omega)=1$ and those for the regular 2nd order QME for
$W_{h}(\omega)=0$. For the case where $W_{h}(\omega)(1-W_{h}(\omega))=0$, they
become independent sums of those due to PT and due to untransformed linear
system-bath couplings.

###  II.3 Representation in the basis of renormalized system eigenstates

For both better conceptual understanding and efficient numerical calculation,
it is convenient to consider the dynamics in the basis of eigenstates of
$\tilde{H}_{0,s}$. Let us denote the $p$th eigenstate and eigenvalue of
$\tilde{H}_{0,s}$ as $|\varphi_{p}\rangle$ and ${\mathcal{E}}_{p}$. Then,

$\tilde{H}_{0,s}=\sum_{p=1}^{N}{\mathcal{E}}_{p}|\varphi_{p}\rangle\langle\varphi_{p}|\
.$ (62)

The transformation matrix $U$ between the site localized states $|j\rangle$’s
and the eigenstates $|\varphi_{p}\rangle$’s can be defined such that
$U_{jp}=\langle j|\varphi_{p}\rangle$. Then,

$|j\rangle=\sum_{p=1}^{N}U_{jp}^{*}|\varphi_{p}\rangle\ .$ (63)

This transformation can be used to express ${\mathcal{T}}_{jk}(t)$ defined by
Eq. (25) in the basis of $|\varphi_{p}\rangle$’s as follows:

${\mathcal{T}}_{jk}(t)=\sum_{p,q=1}^{N}U_{jp}^{*}U_{kq}e^{i\delta{\mathcal{E}}_{pq}t/\hbar}|\varphi_{p}\rangle\langle\varphi_{q}|\
,$ (64)

where $\delta{\mathcal{E}}_{pq}={\mathcal{E}}_{p}-{\mathcal{E}}_{q}$. Let us
also introduce $S_{pq}(t)$’s such that

$\tilde{\sigma}_{I}(t)=\sum_{p,q=1}^{N}S_{pq}(t)|\varphi_{p}\rangle\langle\varphi_{q}|.$
(65)

Then, with some arrangement of dummy summation indices, it is straightforward
to show that

$\displaystyle[{\mathcal{T}}_{jk}(t),{\mathcal{T}}_{j^{\prime}k^{\prime}}(\tau)\tilde{\sigma}_{I}(t)]$
$\displaystyle=\sum_{p,q=1}^{N}|\varphi_{p}\rangle\langle\varphi_{q}|\sum_{p^{\prime},q^{\prime}=1}^{N}\Big{(}\delta_{q^{\prime}q}\sum_{r=1}^{N}U_{jp}^{*}U_{kr}U_{j^{\prime}r}U_{k^{\prime}p^{\prime}}$
$\displaystyle\hskip 108.405pt\times
e^{i\delta{\mathcal{E}}_{pp^{\prime}}t/\hbar}e^{i\delta{\mathcal{E}}_{p^{\prime}r}(t-\tau)/\hbar}$
$\displaystyle\hskip 72.26999pt-
U_{jq^{\prime}}^{*}U_{kq}U_{j^{\prime}p}^{*}U_{k^{\prime}p^{\prime}}e^{i(\delta{\mathcal{E}}_{pp^{\prime}}-\delta{\mathcal{E}}_{qq^{\prime}})t/\hbar}$
$\displaystyle\hskip 108.405pt\times
e^{i\delta{\mathcal{E}}_{p^{\prime}p}(t-\tau)/\hbar}\Big{)}S_{p^{\prime}q^{\prime}}(t).$
(66)

The above expression, in combination with Eq. (39), can be used to express
${\mathcal{R}}(t)\tilde{\sigma}_{I}(t)$ in the basis of
$|\varphi_{p}\rangle$’s. For more compact expression, let us introduce

$\displaystyle{\mathcal{W}}_{jk,j^{\prime}k^{\prime}}^{pq}(t)=\tilde{J}_{jk}\tilde{J}_{j^{\prime}k^{\prime}}\int_{0}^{t}d\tau\
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}$ $\displaystyle\hskip
86.72377pt\times\Big{(}e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t-\tau)}-1\Big{)},$
(67)
$\displaystyle{\mathcal{Y}}_{j,j^{\prime}k^{\prime}}^{pq}(t)=\tilde{J}_{j^{\prime}k^{\prime}}\int_{0}^{t}d\tau
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}{\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t-\tau),$
(68) $\displaystyle{\mathcal{X}}_{jj^{\prime}}^{pq}(t)=\int_{0}^{t}d\tau
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}{\mathcal{C}}_{jj^{\prime}}(t-\tau).$
(69)

Then,

${\mathcal{R}}(t)\tilde{\sigma}_{I}(t)=\sum_{p,q=1}^{N}\sum_{p^{\prime},q^{\prime}=1}^{N}|\varphi_{p}\rangle\langle\varphi_{q}|{\mathcal{R}}_{pq}^{p^{\prime}q^{\prime}}(t)S_{p^{\prime}q^{\prime}}(t),$
(70)

where

$\displaystyle{\mathcal{R}}^{p^{\prime}q^{\prime}}_{pq}(t)=\frac{1}{\hbar^{2}}\sum_{j,k=1}^{N}\sum_{j^{\prime},k^{\prime}=1}^{N}$
$\displaystyle\hskip
36.135pt\times\Bigg{\\{}\delta_{q^{\prime}q}\sum_{r=1}^{N}U_{jp}^{*}U_{kr}U_{j^{\prime}r}^{*}U_{k^{\prime}p^{\prime}}e^{i\delta{\mathcal{E}}_{pp^{\prime}}t/\hbar}$
$\displaystyle\hskip
72.26999pt\times\left({\mathcal{W}}_{jk,j^{\prime}k^{\prime}}^{p^{\prime}r}(t)+\delta_{jk}{\mathcal{Y}}_{j,j^{\prime}k^{\prime}}^{p^{\prime}r}(t)\right.$
$\displaystyle\hskip
72.26999pt\left.+\delta_{j^{\prime}k^{\prime}}{\mathcal{Y}}_{j^{\prime},kj}^{p^{\prime}r}(t)+\delta_{jk}\delta_{j^{\prime}k^{\prime}}{\mathcal{X}}_{jj^{\prime}}^{p^{\prime}r}(t)\right)$
$\displaystyle\hskip 50.58878pt-
U_{jq^{\prime}}^{*}U_{kq}U_{j^{\prime}p}^{*}U_{k^{\prime}p^{\prime}}e^{i(\delta{\mathcal{E}}_{pp^{\prime}}-\delta{\mathcal{E}}_{qq^{\prime}})t/\hbar}$
$\displaystyle\hskip
72.26999pt\times\left({\mathcal{W}}_{jk,j^{\prime}k^{\prime}}^{p^{\prime}p}(t)+\delta_{jk}{\mathcal{Y}}_{j,j^{\prime}k^{\prime}}^{p^{\prime}p}(t)\right.$
$\displaystyle\hskip
72.26999pt\left.+\delta_{j^{\prime}k^{\prime}}{\mathcal{Y}}_{j^{\prime},kj}^{p^{\prime}p}(t)+\delta_{jk}\delta_{j^{\prime}k^{\prime}}{\mathcal{X}}_{jj^{\prime}}^{p^{\prime}p}(t)\right)\Bigg{\\}}$
$\displaystyle\hskip 36.135pt+[{\rm c.c.},p\leftrightarrow
q,p^{\prime}\leftrightarrow q^{\prime}].$ (71)

In the above expression, the last line represents complex conjugates along
with the interchange of indices, $p\leftrightarrow q$ and
$p^{\prime}\leftrightarrow q^{\prime}$, of all previous terms.

For the calculation of the first order inhomogeneous term, the commutator of
system operators in Eq. (57) can be calculated in a manner similar to Eq.
(66). The resulting expression is

$\displaystyle[{\mathcal{T}}_{jk}(t),{\mathcal{T}}_{j^{\prime}k^{\prime}}(0)]$
$\displaystyle=$
$\displaystyle\sum_{p,q=1}^{N}|{\varphi}_{q}\rangle\langle\varphi_{q}|$ (72)
$\displaystyle\times\sum_{r}\left(U_{jp}^{*}U_{kr}U_{j^{\prime}r}^{*}U_{k^{\prime}q}e^{i\delta{\mathcal{E}}_{pr}t/\hbar}\right.$
$\displaystyle\left.-U_{jr}^{*}U_{kq}U_{j^{\prime}p}^{*}U_{k^{\prime}r}e^{i\delta{\mathcal{E}}_{rq}t/\hbar}\right),$

which can also be obtained from Eq. (66) by replacing
$S_{p^{\prime}q^{\prime}}(t)$ with $\delta_{p^{\prime}q^{\prime}}$ and
assuming $\tau=0$. Combining the above expression with Eq. (58), one can
express the first order inhomogeneous term, Eq. (57), as follows:

$\displaystyle{\mathcal{I}}^{(1)}(t)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\sum_{p,q=1}^{N}|\varphi_{p}\rangle\langle\varphi_{q}|\sum_{r=1}^{N}\sum_{j,k=1}^{N}\sum_{j^{\prime},k^{\prime}=1}^{N}\sigma_{j^{\prime}k^{\prime}}(0)w_{j^{\prime}k^{\prime}}$
(73)
$\displaystyle\times\left(U_{jp}^{*}U_{kr}U_{j^{\prime}r}^{*}U_{k^{\prime}q}e^{i\delta{\mathcal{E}}_{pr}t/\hbar}\right.$
$\displaystyle\hskip
14.45377pt\left.-U_{jr}^{*}U_{kq}U_{j^{\prime}p}^{*}U_{k^{\prime}r}e^{i\delta{\mathcal{E}}_{rq}t/\hbar}\right)$
$\displaystyle\times\Big{\\{}\tilde{J}_{jk}(e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)}f_{jk,k^{\prime}}(t)-1)$
$\displaystyle\hskip
14.45377pt+\delta_{jk}\left(M_{j,j^{\prime}k^{\prime}}(t)+h_{j,k^{\prime}}(t)\right)\Big{\\}}.$

Similar expressions for the second order inhomogeneous term, Eq.
(LABEL:eq:ihom2), can be obtained by replacing $S_{p^{\prime}q^{\prime}}(t)$
in the commutator of Eq. (66) with
$U_{j^{\prime\prime}p^{\prime}}^{*}U_{k^{\prime\prime}q^{\prime}}$. The
resulting expression is as follows:

$\displaystyle[{\mathcal{T}}_{jk}(t),{\mathcal{T}}_{j^{\prime}k^{\prime}}(\tau){\mathcal{T}}_{j^{\prime\prime}k^{\prime\prime}}(0)]$
$\displaystyle=\sum_{p,q=1}^{N}|\varphi_{p}\rangle\langle\varphi_{q}|\sum_{p^{\prime},q^{\prime}=1}^{N}\Big{(}\delta_{q^{\prime}q}\sum_{r=1}^{N}U_{jp}^{*}U_{kr}U_{j^{\prime}r}U_{k^{\prime}p^{\prime}}$
$\displaystyle\hskip 72.26999pt\times
e^{i\delta{\mathcal{E}}_{pp^{\prime}}t/\hbar}e^{i\delta{\mathcal{E}}_{p^{\prime}r}(t-\tau)/\hbar}$
$\displaystyle\hskip 36.135pt-
U_{jq^{\prime}}^{*}U_{kq}U_{j^{\prime}p}^{*}U_{k^{\prime}p^{\prime}}e^{i(\delta{\mathcal{E}}_{pp^{\prime}}-\delta{\mathcal{E}}_{qq^{\prime}})t/\hbar}$
$\displaystyle\hskip 72.26999pt\times
e^{i\delta{\mathcal{E}}_{p^{\prime}p}(t-\tau)/\hbar}\Big{)}U_{j^{\prime\prime}p^{\prime}}^{*}U_{k^{\prime\prime}q^{\prime}}.$
(74)

Combining the above expression with bath correlation functions introduced in
Appendix C, one can express the second order inhomogeneous term as follows:

$\displaystyle{\mathcal{I}}^{(2)}(t)$ $\displaystyle=$
$\displaystyle-\sum_{j,k=1}^{N}\sum_{j^{\prime},k^{\prime}=1}^{N}\sum_{j^{\prime\prime},k^{\prime\prime}=1}^{N}\sigma_{j^{\prime\prime}k^{\prime\prime}}(0)\sum_{p,q=1}^{N}|\varphi_{p}\rangle\langle\varphi_{q}|$
(75)
$\displaystyle\times\sum_{p^{\prime},q^{\prime}=1}^{N}\sum_{r=1}^{N}U_{j^{\prime\prime}p^{\prime}}^{*}U_{k^{\prime\prime}q^{\prime}}e^{i\delta{\mathcal{E}}_{pp^{\prime}}t/\hbar}$
$\displaystyle\times\Big{(}\delta_{q^{\prime}q}U_{jp}^{*}U_{kr}U_{j^{\prime}r}U_{k^{\prime}p^{\prime}}$
$\displaystyle\hskip
14.45377pt-\delta_{rp}U_{jq^{\prime}}^{*}U_{kq}U_{j^{\prime}p}^{*}U_{k^{\prime}p}e^{-i\delta{\mathcal{E}}_{qq^{\prime}}t/\hbar}\Big{)}$
$\displaystyle\hskip
14.45377pt\times\Big{(}J_{jk}J_{j^{\prime}k^{\prime}}\tilde{F}_{jk,j^{\prime}k^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t;\delta{\mathcal{E}}_{p^{\prime}r})$
$\displaystyle\hskip
28.90755pt+\delta_{jk}J_{j^{\prime}k^{\prime}}\tilde{H}_{j,j^{\prime}k^{\prime}}^{(1),j^{\prime\prime}k^{\prime\prime}}(t;\delta{\mathcal{E}}_{p^{\prime}r})$
$\displaystyle\hskip
28.90755pt+J_{jk}\delta_{j^{\prime}k^{\prime}}\tilde{H}_{jk,j^{\prime}}^{(2),j^{\prime\prime}k^{\prime\prime}}(t;\delta{\mathcal{E}}_{p^{\prime}r})$
$\displaystyle\hskip
28.90755pt+\delta_{jk}\delta_{j^{\prime}k^{\prime}}\tilde{L}_{j,j^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t;\delta{\mathcal{E}}_{p^{\prime}r})\Big{)}$
$\displaystyle+{\rm H.c.},$

where the definitions of Eqs. (158)-(161) in Appendix C have been used.

##  III Two state system with independent bath

###  III.1 Model and general expressions

Let us consider the simplest case where there are two site local system states
($N=2$), with only one electronic coupling $J=J_{12}=J_{21}$ and only one
Debye-Waller factor, $w=w_{12}=w_{21}$. Then, introducing
$\theta=\theta_{1}^{\dagger}\theta_{2}$ for the present case, the renormalized
zeroth order system Hamiltonian and the first order term of the Hamiltonian
defined by Eq. (21) reduce to

$\displaystyle\tilde{H}_{0,s}$ $\displaystyle=$
$\displaystyle\tilde{E}_{1}|1\rangle\langle 1|+\tilde{E}_{2}|2\rangle\langle
2|+Jw(|1\rangle\langle 2|+|2\rangle\langle 1|),$ (76)
$\displaystyle\tilde{H}_{1}$ $\displaystyle=$ $\displaystyle
D_{1}|1\rangle\langle 1|+D_{2}|2\rangle\langle 2|$ (77) $\displaystyle+$
$\displaystyle J(\theta-w)|1\rangle\langle
2|+J(\theta^{\dagger}-w)|2\rangle\langle 1|,$

where $w=\langle\theta\rangle=\langle\theta^{\dagger}\rangle$, and $D_{1}$ and
$D_{2}$ have been defined by Eq. (23). Note also that, for the present case,

$\tilde{J}_{12}=\tilde{J}_{21}=wJ,$ (78)

and $\tilde{J}_{11}=\tilde{J}_{22}=0$ by definition. The eigenvalues of
$\tilde{H}_{0,s}$ are given by

${\mathcal{E}}_{1,2}=\frac{\tilde{E}_{1}+\tilde{E}_{2}}{2}\pm\frac{\tilde{E}_{1}-\tilde{E}_{2}}{2}\sec(2\eta),\\\
$ (79)

with $1$ denoting the $+$ sign and $2$ denoting $-$ sign on the right hand
side, and $\eta=\tan^{-1}(2Jw/(\tilde{E}_{1}-\tilde{E}_{2}))/2$. The
corresponding eigenstates are given by

$\displaystyle|\varphi_{1}\rangle=\cos\eta|1\rangle+\sin\eta|2\rangle,$ (80)
$\displaystyle|\varphi_{2}\rangle=-\sin\eta|1\rangle+\cos\eta|2\rangle.$ (81)

Thus, for the present case, $U_{11}=U_{22}=\cos\eta$ and
$U_{21}=-U_{12}=\sin\eta$.

Let us also assume that the spectral densities for sites 1 and 2 are the same,
which we denote as ${\mathcal{J}}(\omega)$. Namely,
${\mathcal{J}}_{11}(\omega)={\mathcal{J}}_{22}(\omega)={\mathcal{J}}(\omega)$.
In addition, let us also assume that
${\mathcal{J}}_{12}(\omega)={\mathcal{J}}_{21}(\omega)=0$. Then,
${\mathcal{J}}^{(1)}_{1,12}(\omega)={\mathcal{J}}^{(1)}_{2,21}(\omega)={\mathcal{J}}(\omega)$
and
${\mathcal{J}}^{(2)}_{12,12}(\omega)={\mathcal{J}}^{(2)}_{21,21}(\omega)=2{\mathcal{J}}(\omega)$.
For all other indices,
${\mathcal{J}}^{(1)}_{j,j^{\prime}k^{\prime}}(\omega)=0$ and
${\mathcal{J}}^{(2)}_{jk,j^{\prime}k^{\prime}}(\omega)=0$.

Then, all the bath correlation functions that enter the p-PQME can be
represented by the following three functions:

$\displaystyle{\mathcal{K}}(t)\equiv{\mathcal{K}}_{12,12}(t)=\frac{2}{\pi\hbar}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}(\omega)}{\omega^{2}}W_{h}(\omega)^{2}$
$\displaystyle\hskip
36.135pt\times\left(\coth\left(\frac{\beta\hbar\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right)$ $\displaystyle\hskip
21.68121pt={\mathcal{K}}_{21,21}(t)=-{\mathcal{K}}_{12,21}(t)=-{\mathcal{K}}_{21,12}(t),$
(82)
$\displaystyle{\mathcal{M}}(t)\equiv{\mathcal{M}}_{1,12}(t)=\frac{1}{\pi}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}(\omega)}{\omega}W_{h}(\omega)(1-W_{h}(\omega))$
$\displaystyle\hskip 36.135pt\times\left(\cos(\omega
t)-i\coth\left(\frac{\beta\hbar\omega}{2}\right)\sin(\omega t)\right)$
$\displaystyle\hskip
21.68121pt={\mathcal{M}}_{2,21}(t)=-{\mathcal{M}}_{1,21}(t)=-{\mathcal{M}}_{2,12}(t),$
(83)
$\displaystyle{\mathcal{C}}(t)\equiv{\mathcal{C}}_{11}(t)=\frac{\hbar}{\pi}\int_{0}^{\infty}d\omega{\mathcal{J}}(\omega)(1-W_{h}(\omega))^{2}$
$\displaystyle\hskip
36.135pt\times\left(\coth\left(\frac{\beta\hbar\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right)$ $\displaystyle\hskip
21.68121pt={\mathcal{C}}_{22}(t).$ (84)

For all other components, ${\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(\omega)=0$,
${\mathcal{M}}_{j,j^{\prime}k^{\prime}}(\omega)=0$, and
${\mathcal{C}}_{jk}(\omega)=0$.

In the site basis, Eq. (38) for the present case can be expressed as follows:

$\displaystyle{\mathcal{R}}(t)\tilde{\sigma}_{I}(t)=\int_{0}^{t}d\tau\left\\{C(t-\tau)\left([{\mathcal{T}}_{11}(t),{\mathcal{T}}_{11}(\tau)\tilde{\sigma}_{I}(t)]\right.\right.$
$\displaystyle\left.\hskip
93.95122pt+[{\mathcal{T}}_{22}(t),{\mathcal{T}}_{22}(\tau)\tilde{\sigma}_{I}(t)]\right)$
$\displaystyle\hskip
7.22743pt+Jw{\mathcal{M}}(t-\tau)\left([\Delta{\mathcal{T}}_{d}(t),\Delta{\mathcal{T}}_{c}(\tau)\tilde{\sigma}_{I}(t)]\right.$
$\displaystyle\hskip
65.04256pt\left.-[\Delta{\mathcal{T}}_{c}(t),\Delta{\mathcal{T}}_{d}(\tau)\tilde{\sigma}_{I}(t)]\right)$
$\displaystyle\hskip
7.22743pt+J^{2}w^{2}\left(e^{-{\mathcal{K}}(t-\tau)}-1\right)\left([{\mathcal{T}}_{12}(t),{\mathcal{T}}_{12}(\tau)\tilde{\sigma}_{I}(t)]\right.$
$\displaystyle\hskip
101.17755pt\left.+[{\mathcal{T}}_{21}(t),{\mathcal{T}}_{21}(\tau)\tilde{\sigma}_{I}(t)]\right)$
$\displaystyle\hskip
7.22743pt+J^{2}w^{2}\left(e^{{\mathcal{K}}(t-\tau)}-1\right)\left([{\mathcal{T}}_{12}(t),{\mathcal{T}}_{21}(\tau)\tilde{\sigma}_{I}(t)]\right.$
$\displaystyle\hskip
101.17755pt\left.\left.+[{\mathcal{T}}_{21}(t),{\mathcal{T}}_{12}(\tau)\tilde{\sigma}_{I}(t)]\right)\right\\}$
$\displaystyle+{\rm H.c.},$ (85)

where

$\displaystyle\Delta{\mathcal{T}}_{d}(t)={\mathcal{T}}_{22}(t)-{\mathcal{T}}_{11}(t),$
(86)
$\displaystyle\Delta{\mathcal{T}}_{c}(t)={\mathcal{T}}_{21}(t)-{\mathcal{T}}_{12}(t).$
(87)

Equation (85) above provides clear insights into the effects of p-PT and is
useful for understanding the steady state behavior of the population dynamics.
However, for general numerical calculation, it is convenient to use
expressions in the basis of eigenstates of $\tilde{H}_{0,s}$, as detailed in
Sec. IIC. For this, ${\mathcal{R}}^{p^{\prime}q^{\prime}}_{pq}(t)$ given by
Eq. (71) needs to be calculated, which in turn requires calculation of
${\mathcal{W}}_{jk,j^{\prime}k^{\prime}}^{pq}(t)$,
${\mathcal{Y}}_{j,j^{\prime}k^{\prime}}^{pq}(t)$, and
${\mathcal{X}}_{jj^{\prime}}^{pq}(t)$ defined by Eqs. (67)-(69). For the
present model of two state systems coupled to independent baths, most of these
are zero except for few functions, as defined below.

First, all of nonzero ${\mathcal{W}}_{jk,j^{\prime}k^{\prime}}^{pq}(t)$’s
defined by Eq. (67) for the present case reduce to one of the following two
functions:

$\displaystyle{\mathcal{W}}_{-}^{pq}(t)\equiv{\mathcal{W}}_{12,12}^{pq}(t)={\mathcal{W}}_{21,21}^{pq}(t)$
$\displaystyle\hskip 14.45377pt=J^{2}w^{2}\int_{0}^{t}d\tau
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}\left(e^{-{\mathcal{K}}(t-\tau)}-1\right),$
(88)
$\displaystyle{\mathcal{W}}_{+}^{pq}(t)\equiv{\mathcal{W}}_{12,21}^{pq}(t)={\mathcal{W}}_{21,12}^{pq}(t)$
$\displaystyle\hskip 14.45377pt=J^{2}w^{2}\int_{0}^{t}d\tau
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}\left(e^{{\mathcal{K}}(t-\tau)}-1\right).$
(89)

All other terms of ${\mathcal{W}}_{jk,j^{\prime}k^{\prime}}^{pq}(t)$ are zero
by definition. Similarly, all of nonzero
${\mathcal{Y}}_{j,j^{\prime}k^{\prime}}^{pq}(t)$’s and
${\mathcal{X}}_{jj^{\prime}}^{pq}(t)$’s can be specified by

$\displaystyle{\mathcal{Y}}^{pq}(t)\equiv{\mathcal{Y}}_{1,12}^{pq}(t)={\mathcal{Y}}_{2,21}^{pq}(t)$
$\displaystyle\hskip 36.135pt=Jw\int_{0}^{t}d\tau
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}{\mathcal{M}}(t-\tau)$
$\displaystyle\hskip
36.135pt=-{\mathcal{Y}}_{1,21}^{pq}(t)=-{\mathcal{Y}}_{2,12}^{pq}(t),$ (90)
$\displaystyle{\mathcal{X}}^{pq}(t)\equiv{\mathcal{X}}^{pq}_{11}(t)={\mathcal{X}}^{pq}_{22}(t)$
$\displaystyle\hskip 36.135pt=\int_{0}^{t}d\tau
e^{i\delta{\mathcal{E}}_{pq}(t-\tau)/\hbar}{\mathcal{C}}(t-\tau).$ (91)

All other terms with indices different from above are zero. Thus, with
${\mathcal{W}}_{\pm}^{pq}(t)$, ${\mathcal{Y}}^{pq}(t)$, and
${\mathcal{X}}^{pq}(t)$ as defined above, all of
${\mathcal{R}}_{pq}^{p^{\prime}q^{\prime}}(t)$ constituting Eq. (71) for the
present case can be calculated.

The first order inhomogeneous term, Eq. (73), involves additional functions
$f_{jk,k^{\prime}}(t)$ and $h_{j,k^{\prime}}(t)$, which can also be simplified
for the present model. Let us define

$\displaystyle f(t)$ $\displaystyle\equiv$
$\displaystyle\exp\left\\{\frac{2i}{\pi\hbar}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}(\omega)}{\omega^{2}}W_{h}(\omega)^{2}\sin(\omega
t)\right\\}$ (92) $\displaystyle=$ $\displaystyle f_{12,1}(t).$

Then, it is easy to show that $f_{21,2}(t)=f(t)$,
$f_{21,1}(t)=f_{12,2}(t)=f^{*}(t)$. On the other hand,
$f_{11,1}(t)=f_{22,1}(t)=f_{11,2}(t)=f_{22,2}(t)=1$. For
$h_{j,k^{\prime}}(t)$, only specification of the following function is
necessary.

$\displaystyle h(t)$ $\displaystyle\equiv$
$\displaystyle\frac{2}{\pi}\int_{0}^{\infty}d\omega\frac{{\mathcal{J}}(\omega)}{\omega}(1-W_{h}(\omega))W_{h}(\omega)\cos(\omega
t)$ (93) $\displaystyle=$ $\displaystyle h_{1,1}(t)=h_{2,2}(t).$

For other cases, $h_{1,2}(t)=h_{2,1}(t)=0$. Calculation of the second order
inhomogeneous term, Eq. (75), although complicated, is straightforward and the
expressions defined in Appendix C can be used along with definitions of all
the constituting functions defined in the present subsection.

###  III.2 Model calculations without inhomogeneous terms for Ohmic spectral
density

In order to demonstrate the feasibility of using the p-PQME expressions
derived here and the effects of the weighting function $W_{h}(\omega$), model
calculations without including the transient inhomogeneous terms are provided
here. For this, it is useful to note that
${\mathcal{R}}_{pq}^{p^{\prime}q^{\prime}}(t)$ for the present case can be
simplified further. Summing up all possible 16 cases of
$j,k,j^{\prime},k^{\prime}=1,2$ in Eq. (71), it is straightforward to show
that it can be simplified as follows:

$\displaystyle{\mathcal{R}}_{pq}^{p^{\prime}q^{\prime}}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar^{2}}\Bigg{\\{}\delta_{q^{\prime}q}e^{i\delta{\mathcal{E}}_{pp^{\prime}}t/\hbar}\sum_{r=1}^{2}\left({\mathcal{A}}_{prrp^{\prime}}{\mathcal{X}}^{p^{\prime}r}(t)\right.$
(94) $\displaystyle\hskip
21.68121pt+({\mathcal{B}}^{(1)}_{pr}{\mathcal{B}}^{(2)}_{rp^{\prime}}-{\mathcal{B}}_{pr}^{(2)}{\mathcal{B}}_{rp^{\prime}}^{(1)}){\mathcal{Y}}^{p^{\prime}r}(t)$
$\displaystyle\hskip
21.68121pt\left.+{\mathcal{C}}_{prrp^{\prime}}^{(1)}{\mathcal{W}}_{-}^{p^{\prime}r}(t)+{\mathcal{C}}_{prrp^{\prime}}^{(2)}{\mathcal{W}}_{+}^{p^{\prime}r}(t)\right)$
$\displaystyle\hskip
7.22743pt-e^{i(\delta{\mathcal{E}}_{pp^{\prime}}-\delta{\mathcal{E}}_{qq^{\prime}})t/\hbar}\left({\mathcal{A}}_{q^{\prime}qpp^{\prime}}{\mathcal{X}}^{p^{\prime}p}(t)\right.$
$\displaystyle\hskip
21.68121pt+({\mathcal{B}}_{q^{\prime}q}^{(1)}{\mathcal{B}}_{pp^{\prime}}^{(2)}-{\mathcal{B}}_{q^{\prime}q}^{(2)}{\mathcal{B}}_{pp^{\prime}}^{(1)}){\mathcal{Y}}^{p^{\prime}p}(t)$
$\displaystyle\hskip
21.68121pt\left.+{\mathcal{C}}_{q^{\prime}qpp^{\prime}}^{(1)}{\mathcal{W}}_{-}^{p^{\prime}p}(t)+{\mathcal{C}}_{q^{\prime}qpp^{\prime}}^{(2)}{\mathcal{W}}_{+}^{p^{\prime}p}(t)\right)\Bigg{\\}}$
$\displaystyle\hskip 36.135pt+[{\rm C.C.},p\leftrightarrow
q,p^{\prime}\leftrightarrow q^{\prime}],$

where

$\displaystyle{\mathcal{A}}_{pqrs}=U_{1p}U_{1q}U_{1r}U_{1s}+U_{2p}U_{2q}U_{2r}U_{2s},$
(95) $\displaystyle{\mathcal{B}}^{(1)}_{pq}=U_{1p}U_{1q}-U_{2p}U_{2q},$ (96)
$\displaystyle{\mathcal{B}}^{(2)}_{pq}=U_{1p}U_{2q}-U_{2p}U_{1q},$ (97)
$\displaystyle{\mathcal{C}}^{(1)}_{pqrs}=U_{1p}U_{2q}U_{1r}U_{2s}+U_{2p}U_{1q}U_{2r}U_{1s},$
(98)
$\displaystyle{\mathcal{C}}^{(2)}_{pqrs}=U_{1p}U_{2q}U_{2r}U_{1s}+U_{2p}U_{1q}U_{1r}U_{2s}.$
(99)

Figure 1: Populations of donor (site $1$) for the case where a step function
that becomes one at $\omega=\omega_{h}$ is used for $W_{h}(\omega)$. The
parameters of the model are as follows: $E_{1}=E_{2}$, $\eta=1$,
$\omega_{c}=200\ {\rm cm^{-1}}$, $T=300\ {\rm K}$, and $J=300\ {\rm cm^{-1}}$.

The bath spectral density being considered is the following Ohmic form with
exponential cutoff:

${\mathcal{J}}(\omega)=\pi\hbar\omega e^{-\omega/\omega_{c}}.$ (100)

More specifically, the following set of parameters are used for the
calculation: $E_{1}=E_{2}$, $\eta=1$, $\omega_{c}=200\ {\rm cm^{-1}}$, $T=300\
{\rm K}$, and $J=300\ {\rm cm^{-1}}$.

There are infinite number of possible choices available for $W_{h}(\omega)$
that satisfies the requirement of Eq. (6). First, let us consider the simplest
case where $W_{h}(\omega)$ is given by the following step function:

$W_{h}(\omega)=\left\\{\begin{array}[]{cl}0,&,\ \omega<\omega_{h}\\\ 1&,\
\omega\geq\omega_{h}\end{array}\right.$ (101)

For this case, ${\mathcal{M}}(t)$ vanishes and the QME becomes a simple sum of
PT part and untransformed part. Figure 1 shows results of population dynamics
for three different choices $\omega_{h}/\omega_{c}$. The result for
$\omega_{h}/\omega_{c}=0.1$ is close to the limit of full PT, and that for
$\omega_{h}/\omega_{c}=10$ is close to the limit of conventional 2nd order
time local QME. As expected, the population for $\omega_{h}/\omega_{c}=1$
exhibits an intermediate character.

In order to observe the effects of ${\mathcal{M}}(t)$, let us consider a more
general weighting function given by

$W_{h}(\omega)=1-e^{-(\omega/\omega_{h})^{\alpha}}.$ (102)

Figure 2 shows results for different values of $\alpha$ for the case where
$\omega_{h}=\omega_{c}$. It is interesting to note that the population
dynamics becomes more coherent as $\alpha$ increases. For $\alpha=4$, the
extent of coherence is even bigger than the case of the step function, which
indicates that ${\mathcal{M}}(t)$ makes positive contribution to coherence in
this case. In this sense, the role of ${\mathcal{M}}(t)$ seems important in
reproducing correct behavior for a given choice of the weighting function.
However, one important practical issue that needs careful examination is
possible divergence or violation of positivity due to ${\mathcal{M}}(t)$.
Based on preliminary numerical studies, it appears that ${\mathcal{M}(t)}$
results in divergence of dynamics for $\alpha\leq 1$ and significant non-
positivity for too large value of $\omega_{h}$ for a given $\alpha>1$ in Eq.
(102). For clear theoretical understanding of this behavior and more
quantitative assessment, more comprehensive numerical study also including the
inhomogeneous terms is necessary.

Figure 2: Populations of donor (site $1$) for the case where $W_{h}(\omega)$
is given by Eq. (102). Results for three different values of $\alpha$ are
compared with that of step function as the weighting function. All other
parameters are the same as in Fig. 1.

##  IV Conclusion

This work has provided a general framework to overcome a known issue of the
original 2nd order PQME, namely, premature over-relaxation of the sluggish
bath, by deriving full expressions for the second order time-local p-PQME. The
main results of this work, represented by Eqs. (71), (73), and (75), can be
applied for any kinds of bath spectral densities and initial system states but
will be particularly useful for the cases where the bath spectral densities
are Ohmic or sub-Ohmic.

Indeed, numerical tests for a simple two level system coupled to an Ohmic bath
demonstrate that appropriate choice of the weighting function $W_{h}(\omega)$
can make the population dynamics recover some level of coherence through
appropriate choice of the functional form of $W_{h}(\omega)$ and its
parameters, adding a new dimension of flexibility in incorporating the PT
approach into a QME calculation. The flexibility in choosing the weighting
function $W_{h}(\omega)$ in the expression for the p-PQME presented here
leaves open various possibilities of adapting or improving the methodology.
For example, variational theorem can be used for its optimization.
Alternatively, benchmarking against numerically exact computational results
followed by a Machine Learning based optimization can potentially lead to an
optimized second order p-PQME that can best approximate exact dynamics. To
this end, full calculations including all the inhomogeneous terms and
benchmarking against a broad range of exact numerical results will be
necessary.

The formulation developed here also will be useful for further extension of PT
based QME approaches. For example, extension to the cases with time dependent
Hamiltonian for driven quantum systems is straightforward. The formulations
and theoretical identities employed in this work will also be useful for the
development of new PT based approaches for general anharmonic bath and for the
formulation of time dependent PT approach.

###### Acknowledgements.

This work was supported by the National Science Foundation (CHE-1900170). The
author also acknowledges support from Korea Institute for Advanced Study
(KIAS) through its KIAS Scholar program.
AUTHOR DECLARATIONS
Conflict of Interest
The author has no conflicts to disclose.
DATA AVAILABILITY
Most data that support the findings of this article are contained in this
article. Additional data are available from the corresponding author upon
reasonable request.

##  Appendix A Derivation of Eq. (39)

The trace of $\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)$ with
$\rho_{b}$ in Eq. (38), namely, its thermal average is expressed as follows:

$\displaystyle\langle\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\rangle=\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle$
$\displaystyle\hskip 21.68121pt-
w_{jk}\langle\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle-
w_{j^{\prime}k^{\prime}}\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\rangle+w_{jk}w_{j^{\prime}k^{\prime}}$
$\displaystyle\hskip 21.68121pt+\delta_{jk}\left(\langle
D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle-\langle
D_{j}(t)\rangle w_{j^{\prime}k^{\prime}}\right)$ $\displaystyle\hskip
21.68121pt+\delta_{j^{\prime}k^{\prime}}\left(\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)D_{j^{\prime}}(\tau)\rangle-
w_{jk}\langle D_{j^{\prime}}(\tau)\rangle\right)$ $\displaystyle\hskip
21.68121pt+\delta_{jk}\delta_{j^{\prime}k^{\prime}}\langle
D_{j}(t)D_{j^{\prime}}(\tau)\rangle$
$\displaystyle=\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle-
w_{jk}w_{j^{\prime}k^{\prime}}$ $\displaystyle\hskip
21.68121pt+\delta_{jk}\langle
D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle$
$\displaystyle\hskip
21.68121pt+\delta_{j^{\prime}k^{\prime}}\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)D_{j^{\prime}}(\tau)\rangle$
$\displaystyle\hskip
21.68121pt+\delta_{jk}\delta_{j^{\prime}k^{\prime}}\langle
D_{j}(t)D_{j^{\prime}}(\tau)\rangle.$ (103)

In obtaining the second equality of the above equation, the identities that
$\langle\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle=w_{j^{\prime}k^{\prime}}$,
$\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\rangle=w_{jk}$, and $\langle
D_{j}(t)\rangle=\langle D_{j^{\prime}}(\tau)\rangle=0$ have been used.

The first term in the second equality of Eq. (103) has the same form as that
for the full PTJang (2009) except for the additional factor
$W_{h}(\omega_{n})$ multiplied with $\delta g_{n,jk}$ and $\delta
g_{n,j^{\prime}k^{\prime}}$, respectively. Thus, following the same procedure
as in previous work,Jang (2009) it can be shown to be

$\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle=w_{jk}w_{j^{\prime}k^{\prime}}e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t-\tau)},$
(104)

where $e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t-\tau)}$ is defined by Eq.
(40).

The third term in the second equality of Eq. (103) can be expressed as

$\displaystyle\langle
D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle$
$\displaystyle=\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))\Big{\langle}(b_{n}e^{-i\omega_{n}t}+b_{n}^{\dagger}e^{i\omega_{n}t})$
$\displaystyle\hskip 14.45377pt\times\prod_{n^{\prime}}e^{\delta
g_{n^{\prime},j^{\prime}k^{\prime}}W_{h}(\omega_{n^{\prime}})(b_{n^{\prime}}^{\dagger}e^{i\omega_{n^{\prime}}\tau}-b_{n^{\prime}}e^{-i\omega_{n^{\prime}}\tau})}\Big{\rangle}.$
(105)

In the above expression, the product of bath modes (indexed by $n^{\prime}$)
can be averaged independently for $n^{\prime}\neq n$, resulting in

$\displaystyle\left\langle\prod_{n^{\prime}\neq n}e^{\delta
g_{n^{\prime},j^{\prime}k^{\prime}}W_{h}(\omega_{n^{\prime}})(b_{n^{\prime}}^{\dagger}e^{i\omega_{n^{\prime}}\tau}-b_{n^{\prime}}e^{-i\omega_{n^{\prime}}\tau})}\right\rangle$
$\displaystyle=\prod_{n^{\prime}\neq
n}e^{-\coth(\beta\hbar\omega_{n^{\prime}}/2)\delta
g_{n^{\prime},j^{\prime}k^{\prime}}^{2}W_{h}(\omega_{n^{\prime}})^{2}/2}.$
(106)

On the other hand, the term for $n^{\prime}=n$ in Eq. (105) need to be
calculated together with the linear term as follows:

$\displaystyle\langle(b_{n}e^{-i\omega_{n}t}+b_{n}^{\dagger}e^{i\omega_{n}t})e^{\delta
g_{n,j^{\prime}k^{\prime}}W_{h}(\omega_{n})(b_{n}^{\dagger}e^{i\omega_{n}\tau}-b_{n}e^{-i\omega_{n}\tau})}\rangle$
$\displaystyle=e^{-i\omega_{n}t}\langle
b_{n}e^{\gamma^{*}b_{n}^{\dagger}-\gamma b_{n}}\rangle+e^{i\omega_{n}t}\langle
b_{n}^{\dagger}e^{\gamma^{*}b_{n}^{\dagger}-\gamma b_{n}}\rangle,$ (107)

where $\gamma$ is defined by

$\gamma=\delta
g_{n,j^{\prime}k^{\prime}}W_{h}(\omega_{n})e^{-i\omega_{n}\tau},$ (108)

and $\gamma^{*}$ is the complex conjugate of $\gamma$. Using the fact that
$e^{\gamma^{*}b_{n}^{\dagger}-\gamma b_{n}}=e^{-\gamma
b_{n}}e^{\gamma^{*}b_{n}^{\dagger}}e^{|\gamma|^{2}/2}$ and the following
identity:

$\displaystyle\langle b_{n}e^{-\gamma
b_{n}}e^{\gamma^{*}b_{n}^{\dagger}}\rangle=-\frac{\partial}{\partial\gamma}\langle
e^{-\gamma b_{n}}e^{\gamma^{*}b_{n}^{\dagger}}\rangle$ $\displaystyle\hskip
36.135pt=\frac{\gamma^{*}}{1-e^{-\beta\hbar\omega_{n}}}e^{-\gamma\gamma^{*}/(1-e^{-\beta\hbar\omega_{n}})},$
(109)

it is straightforward to show that

$\langle b_{n}e^{\gamma^{*}b_{n}^{\dagger}-\gamma
b_{n}}\rangle=\frac{\gamma^{*}}{1-e^{-\beta\hbar\omega_{n}}}e^{-\coth(\beta\hbar\omega_{n}/2)|\gamma|^{2}/2}.$
(110)

Similarly, using the following identities:

$\displaystyle\langle e^{-\gamma
b_{n}}b_{n}^{\dagger}e^{\gamma^{*}b_{n}^{\dagger}}\rangle=\frac{\partial}{\partial\gamma^{*}}\langle
e^{-\gamma b_{n}}e^{\gamma^{*}b_{n}^{\dagger}}\rangle$ $\displaystyle\hskip
36.135pt=-\frac{\gamma}{1-e^{-\beta\hbar\omega_{n}}}e^{-\gamma\gamma^{*}/(1-e^{-\beta\hbar\omega_{n}})},$
(111) $\displaystyle b_{n}^{\dagger}e^{-\gamma b_{n}}-e^{-\gamma
b_{n}}b_{n}^{\dagger}=\gamma e^{-\gamma b_{n}},$ (112)

it is also straightforward to show that

$\langle b_{n}^{\dagger}e^{\gamma^{*}b_{n}^{\dagger}-\gamma
b_{n}}\rangle=-\frac{\gamma
e^{-\beta\hbar\omega_{n}}}{1-e^{-\beta\hbar\omega_{n}}}e^{-\coth(\beta\hbar\omega_{n}/2)|\gamma|^{2}/2}.$
(113)

Combining Eqs. (110) and (113) with the definition of Eq. (108) leads to

$\displaystyle\langle(b_{n}e^{-i\omega_{n}t}+b_{n}^{\dagger}e^{i\omega_{n}t})e^{\delta
g_{n,j^{\prime}k^{\prime}}W_{h}(\omega_{n})(b_{n}^{\dagger}e^{i\omega_{n}\tau}-b_{n}e^{-i\omega_{n}\tau})}\rangle$
$\displaystyle=\delta
g_{n,j^{\prime}k^{\prime}}W_{h}(\omega_{n})e^{-\coth(\beta\hbar\omega_{n}/2)\delta
g_{n,j^{\prime}k^{\prime}}^{2}W_{h}(\omega_{n})^{2}/2}$
$\displaystyle\times\left(\cos(\omega_{n}(t-\tau))-i\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\sin(\omega_{n}(t-\tau))\right).$

Employing this identity in Eq. (105), one can find that

$\displaystyle\langle
D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\rangle$
$\displaystyle=w_{j^{\prime}k^{\prime}}\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))\delta
g_{n,j^{\prime}k^{\prime}}W_{h}(\omega_{n})$ $\displaystyle\hskip
14.45377pt\times\left(\cos(\omega_{n}(t-\tau))-i\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\sin(\omega_{n}(t-\tau))\right).$
(115)

The fourth term in the second equality of Eq. (103) can be expressed as

$\displaystyle\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)D_{j^{\prime}}(\tau)\rangle$
$\displaystyle=\sum_{n}\hbar\omega_{n}g_{n,j^{\prime}}(1-W_{h}(\omega_{n}))$
$\displaystyle\hskip 14.45377pt\times\Big{\langle}\prod_{n^{\prime}}e^{\delta
g_{n^{\prime},jk}W_{h}(\omega_{n^{\prime}})(b_{n^{\prime}}^{\dagger}e^{i\omega_{n^{\prime}}t}-b_{n^{\prime}}e^{-i\omega_{n^{\prime}}t})}$
$\displaystyle\hskip
21.68121pt\times(b_{n}e^{-i\omega_{n}\tau}+b_{n}^{\dagger}e^{i\omega_{n}\tau})\Big{\rangle}.$
(116)

The thermal average in the above expression can be calculated in a manner
similar to the third term in the second equality of Eq. (103) as described
above, but with a different definition of $\gamma=\delta
g_{n,jk}W_{h}(\omega_{n})e^{-i\omega_{n}t}$. Employing the following identity:

$e^{\gamma^{*}b_{n}^{\dagger}}b_{n}=(b_{n}-\gamma^{*})e^{\gamma^{*}b_{n}^{\dagger}},$
(117)

along with Eq. (109), one can easily show that

$\langle e^{-\gamma
b_{n}}e^{\gamma^{*}b_{n}^{\dagger}}b_{n}\rangle=\gamma^{*}\frac{e^{-\beta\hbar\omega_{n}}}{1-e^{-\beta\hbar\omega_{n}}}e^{-\gamma\gamma^{*}/(1-e^{-\beta\hbar\omega_{n}})}.$
(118)

Combining the above identity with Eq. (111) leads to the following expression:

$\displaystyle\langle e^{\delta
g_{n,jk}W_{h}(\omega_{n})(b_{n}^{\dagger}e^{i\omega_{n}t}-b_{n}e^{-i\omega_{n}t})}(b_{n}e^{-i\omega_{n}\tau}+b_{n}^{\dagger}e^{i\omega_{n}\tau})\rangle$
$\displaystyle=\delta
g_{n,kj}W_{h}(\omega_{n})e^{-\coth(\beta\hbar\omega_{n}/2)\delta
g_{n,jk}^{2}W_{h}(\omega_{n})^{2}/2}$
$\displaystyle\times\left(\cos(\omega_{n}(t-\tau))-i\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\sin(\omega_{n}(t-\tau))\right),$

where note the appearance of $\delta g_{n,kj}=-\delta g_{n,jk}$ on the
righthand side of the above equation. As a result, Eq. (116) can be expressed
as

$\displaystyle\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)D_{j^{\prime}}(\tau)\rangle$
$\displaystyle=w_{jk}\sum_{n}\hbar\omega_{n}g_{n,j^{\prime}}(1-W_{h}(\omega_{n}))\delta
g_{n,kj}W_{h}(\omega_{n})$
$\displaystyle\times\left(\cos(\omega_{n}(t-\tau))-i\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\sin(\omega_{n}(t-\tau))\right).$
(120)

Calculation of the last term in the second equality of Eq. (103) is also
straightforward. The resulting expression is as follows:

$\displaystyle\langle
D_{j}(t)D_{j^{\prime}}(\tau)\rangle=\sum_{n}\hbar^{2}\omega_{n}^{2}g_{n,j}g_{n,j^{\prime}}(1-W_{h}(\omega_{n}))^{2}$
$\displaystyle\hskip 36.135pt\times(\langle b_{n}b_{n}^{\dagger}\rangle
e^{-i\omega_{n}(t-\tau)}+\langle b_{n}^{\dagger}b_{n}\rangle
e^{i\omega_{n}(t-\tau)})$ $\displaystyle=$
$\displaystyle\sum_{n}\hbar^{2}\omega_{n}^{2}g_{n,j}g_{n,j^{\prime}}(1-W_{h}(\omega_{n}))^{2}$
$\displaystyle\times\left(\coth\left(\frac{\beta\hbar\omega_{n}}{2}\right)\cos(\omega_{n}(t-\tau))-i\sin(\omega_{n}(t-\tau))\right).$

In the second equality of the above equation, the fact that $\langle
b_{n}^{\dagger}b_{n}\rangle=\langle
b_{n}b_{n}^{\dagger}\rangle-1=e^{-\beta\hbar\omega_{n}}/(1-e^{-\beta\hbar\omega_{n}})$
has been used.

##  Appendix B Evaluation of the bath correlation function in the first order
inhomogeneous term

The bath portion of the first order inhomogeneous term, Eq. (57), can be
expressed as follows:

$\displaystyle
Tr_{b}\\{\tilde{B}_{jk}(t)\delta\tilde{\rho}_{b,j^{\prime}k^{\prime}}\\}=J_{jk}Tr_{b}\left\\{\theta_{k^{\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}\rho_{b}\right\\}$
$\displaystyle\hskip 36.135pt-
J_{jk}w_{jk}w_{j^{\prime}k^{\prime}}+\delta_{jk}Tr_{b}\left\\{\theta_{k^{\prime}}D_{j}(t)\theta_{j^{\prime}}^{\dagger}\rho_{b}\right\\}.$
(122)

In the first term on the righthand side of the above expression, product of
the first three operators within the trace operation can be expressed as

$\displaystyle\theta_{k^{\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)=\theta_{k^{\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{k^{\prime}}^{\dagger}\theta_{k^{\prime}}$
$\displaystyle=e^{\sum_{n}\delta
g_{n,jk}W_{h}(\omega_{n})(\theta_{k^{\prime}}b_{n}^{\dagger}\theta_{k^{\prime}}^{\dagger}e^{i\omega_{n}t}-\theta_{k^{\prime}}b_{n}\theta_{k^{\prime}}^{\dagger}e^{-i\omega_{n}t})}\theta_{k^{\prime}}$
(123)

In the exponent on the righthand side of the above expression,

$\displaystyle\theta_{k^{\prime}}b_{n}^{\dagger}\theta_{k^{\prime}}^{\dagger}$
$\displaystyle=$ $\displaystyle
b_{n}^{\dagger}-W_{h}(\omega_{n})g_{n,k^{\prime}}[b_{n}^{\dagger}-b_{n},b_{n}^{\dagger}]$
(124) $\displaystyle=$ $\displaystyle
b_{n}^{\dagger}+W_{h}(\omega_{n})g_{n,k^{\prime}},$
$\displaystyle\theta_{k^{\prime}}b_{n}\theta_{k^{\prime}}^{\dagger}$
$\displaystyle=$ $\displaystyle
b_{n}-W_{h}(\omega_{n})g_{n,k^{\prime}}[b_{n}^{\dagger}-b_{n},b_{n}]$ (125)
$\displaystyle=$ $\displaystyle b_{n}+W_{h}(\omega_{n})g_{n,k^{\prime}}.$

Therefore,

$\displaystyle\sum_{n}\delta
g_{n,jk}W_{h}(\omega_{n})(\theta_{k^{\prime}}b_{n}^{\dagger}\theta_{k^{\prime}}^{\dagger}e^{i\omega_{n}t}-\theta_{k^{\prime}}b_{n}\theta_{k^{\prime}}^{\dagger}e^{-i\omega_{n}t})$
$\displaystyle=\sum_{n}\delta
g_{n,jk}W_{h}(\omega_{n})(b_{n}^{\dagger}e^{i\omega_{n}t}-b_{n}e^{-i\omega_{n}t})$
$\displaystyle+2i\sum_{n}\delta
g_{n,jk}g_{n,k^{\prime}}W_{h}(\omega_{n})^{2}\sin(\omega_{n}t).$ (126)

Employing the above identity and using the definition of Eq. (59), one can
show that Eq. (123) can be expressed as

$\displaystyle\theta_{k^{\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)=\theta_{k^{\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{k^{\prime}}^{\dagger}\theta_{k^{\prime}}$
$\displaystyle=e^{\sum_{n}\delta
g_{n,jk}W_{h}(\omega_{n})(\theta_{k^{\prime}}b_{n}^{\dagger}\theta_{k^{\prime}}^{\dagger}e^{i\omega_{n}t}-\theta_{k^{\prime}}b_{n}\theta_{k^{\prime}}^{\dagger}e^{-i\omega_{n}t})}\theta_{k^{\prime}}$
$\displaystyle=f_{jk,k^{\prime}}(t)\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{k^{\prime}}.$
(127)

The above identity implies that Eq. (122) can be simplified to

$\displaystyle
Tr_{b}\left\\{\theta_{k^{\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}\rho_{b}\right\\}$
$\displaystyle=f_{jk,k^{\prime}}(t)\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{k^{\prime}}\theta_{j^{\prime}}^{\dagger}\rangle$
$\displaystyle=f_{jk,k^{\prime}}(t)w_{jk}w_{j^{\prime}k^{\prime}}e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)},$
(128)

where Eq. (104) has been used.

On the other hand, in the last term of Eq. (122), product of the first two
operators within the trace operation can be expressed as

$\displaystyle\theta_{k^{\prime}}D_{j}(t)=\theta_{k^{\prime}}D_{j}(t)\theta_{k^{\prime}}^{\dagger}\theta_{k^{\prime}}$
$\displaystyle=\sum_{n}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))$
$\displaystyle\hskip
36.135pt\times(\theta_{k^{\prime}}b_{n}\theta_{k^{\prime}}^{\dagger}e^{-i\omega_{n}t}+\theta_{k^{\prime}}b_{n}^{\dagger}\theta_{k^{\prime}}^{\dagger}e^{i\omega_{n}t})\theta_{k^{\prime}}$
$\displaystyle=D_{j}(t)\theta_{k^{\prime}}+2\theta_{k^{\prime}}\sum_{n}\hbar\omega_{n}g_{n,j}g_{n,k^{\prime}}$
$\displaystyle\hskip
50.58878pt\times(1-W_{h}(\omega_{n}))W_{h}(\omega_{n})\cos(\omega_{n}t),$
(129)

where Eqs. (124) and (125) have been used. The above identity leads to the
following expression for the trace of the bath operators in the last term of
Eq. (122):

$Tr_{b}\left\\{\theta_{k^{\prime}}D_{j}(t)\theta_{j^{\prime}}^{\dagger}\rho_{b}\right\\}={\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)+h_{j,k^{\prime}}(t)w_{j^{\prime}k^{\prime}}\
,$ (130)

where Eq. (115) for $\tau=0$, the definition of Eq. (41), and the definition
of Eq. (60) have been used. Inserting Eqs. (128) and (130) into Eq. (122), one
can then obtain the following expression:

$\displaystyle
Tr_{b}\\{\tilde{B}_{jk}(t)\delta\tilde{\rho}_{b,j^{\prime}k^{\prime}}\\}$
$\displaystyle=J_{jk}w_{jk}(f_{jk,k^{\prime}}(t)e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t)}-1)w_{j^{\prime}k^{\prime}}$
$\displaystyle+\delta_{jk}({\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t)+h_{j,k^{\prime}}(t))w_{j^{\prime}k^{\prime}}.$
(131)

It is easy to confirm that this expression is equivalent to Eq. (58).

##  Appendix C Evaluation of the bath correlation function in the second order
inhomogeneous term

The trace over the bath in the second order inhomogeneous term, Eq.
(LABEL:eq:ihom2), can be expressed as follows:

$\displaystyle
Tr_{b}\left\\{\tilde{B}_{jk}(t)\tilde{B}_{j^{\prime}k^{\prime}}(\tau)\delta\tilde{\rho}_{b,j^{\prime\prime}k^{\prime\prime}}\right\\}$
$\displaystyle=Tr_{b}\left\\{\left(J_{jk}(\theta_{j}^{\dagger}(t)\theta_{k}(t)-w_{jk})+\delta_{jk}D_{j}(t)\right)\right.$
$\displaystyle\times\left(J_{j^{\prime}k^{\prime}}(\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)-w_{j^{\prime}k^{\prime}})+\delta_{j^{\prime}k^{\prime}}D_{j^{\prime}}(\tau)\right)$
$\displaystyle\times\left.\left(\theta_{j^{\prime\prime}}^{\dagger}\rho_{b}\theta_{k^{\prime\prime}}-w_{j^{\prime\prime}k^{\prime\prime}}\rho_{b}\right)\right\\}$
$\displaystyle=J_{jk}J_{j^{\prime}k^{\prime}}F^{j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}k^{\prime}}(t,\tau)+\delta_{jk}J_{j^{\prime}k^{\prime}}H^{(1),j^{\prime\prime}k^{\prime\prime}}_{j,j^{\prime}k^{\prime}}(t,\tau)$
$\displaystyle+J_{jk}\delta_{j^{\prime}k^{\prime}}H^{(2),j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}}(\tau,t)^{*}+\delta_{jk}\delta_{j^{\prime}k^{\prime}}L^{j^{\prime\prime}k^{\prime\prime}}_{j,j^{\prime}}(t,\tau),$
(132)

where

$\displaystyle
F^{j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}k^{\prime}}(t,\tau)$
$\displaystyle=Tr_{b}\left\\{(\theta_{j}^{\dagger}(t)\theta_{k}(t)-w_{jk})(\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)-w_{j^{\prime}k^{\prime}})\right.$
$\displaystyle\hskip
72.26999pt\left.\times(\theta_{j^{\prime\prime}}^{\dagger}\rho_{b}\theta_{k^{\prime\prime}}-w_{j^{\prime\prime}k^{\prime\prime}}\rho_{b})\right\\},$
(133) $\displaystyle
H_{j,j^{\prime}k^{\prime}}^{(1),j^{\prime\prime}k^{\prime\prime}}(t,\tau)=Tr_{b}\left\\{D_{j}(t)(\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)-w_{j^{\prime}k^{\prime}})\right.$
$\displaystyle\hskip
72.26999pt\left.\times(\theta_{j^{\prime\prime}}^{\dagger}\rho_{b}\theta_{k^{\prime\prime}}-w_{j^{\prime\prime}k^{\prime\prime}}\rho_{b})\right\\},$
(134) $\displaystyle
H_{jk,j^{\prime}}^{(2),j^{\prime\prime}k^{\prime\prime}}(\tau,t)^{*}=Tr_{b}\left\\{(\theta_{j}^{\dagger}(t)\theta_{k}(t)-w_{jk})D_{j^{\prime}}(\tau)\right.$
$\displaystyle\hskip
72.26999pt\left.\times(\theta_{j^{\prime\prime}}^{\dagger}\rho_{b}\theta_{k^{\prime\prime}}-w_{j^{\prime\prime}k^{\prime\prime}}\rho_{b})\right\\},$
(135) $\displaystyle
L_{j,j^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t,\tau)=Tr_{b}\left\\{D_{j}(t)D_{j^{\prime}}(\tau)\right.$
$\displaystyle\hskip
72.26999pt\left.\times(\theta_{j^{\prime\prime}}^{\dagger}\rho_{b}\theta_{k^{\prime\prime}}-w_{j^{\prime\prime}k^{\prime\prime}}\rho_{b})\right\\}.$
(136)

Further calculation of each of the terms above is straightforward as described
below.

The first bath term in Eq. (132),
$F^{j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}k^{\prime}}(t,\tau)$
defined by Eq. (133), can be expanded further and expressed as follows:

$\displaystyle
F^{j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}k^{\prime}}(t,\tau)=\langle\theta_{k^{\prime\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}^{\dagger}(\tau)\theta_{j^{\prime\prime}}\rangle$
$\displaystyle\hskip 21.68121pt-
w_{jk}\langle\theta_{k^{\prime\prime}}\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip 21.68121pt-
w_{j^{\prime}k^{\prime}}\langle\theta_{k^{\prime\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip
21.68121pt+w_{jk}w_{j^{\prime}k^{\prime}}\langle\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip 21.68121pt-
w_{j^{\prime\prime}k^{\prime\prime}}\langle(\theta_{j}^{\dagger}(t)\theta_{k}(t)-w_{jk})(\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}^{\dagger}(\tau)-w_{j^{\prime}k^{\prime}})\rangle.$
(137)

The first term in the above expression can be calculated as follows:

$\displaystyle\langle\theta_{k^{\prime\prime}}\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}^{\dagger}(\tau)\theta_{j^{\prime\prime}}\rangle$
$\displaystyle=f_{jk,k^{\prime\prime}}(t)f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)\langle\theta_{j}^{\dagger}(t)\theta_{k}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}\rangle$
$\displaystyle=f_{jk,k^{\prime\prime}}(t)f_{jk^{\prime},k^{\prime\prime}}(\tau)w_{jk}w_{j^{\prime}k^{\prime}}w_{j^{\prime\prime}k^{\prime\prime}}$
$\displaystyle\hskip 36.135pt\times
e^{-{\mathcal{K}}_{jk,j^{\prime\prime}k^{\prime\prime}}(t)-{\mathcal{K}}_{j^{\prime}k^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t-\tau)}.$
(138)

For the second and third terms of Eq. (137), Eq. (128) can be used. The last
term of Eq. (137) corresponds to the first term that appears in the evaluation
of ${\mathcal{R}}(t)$. Combining all of these, one can show that

$\displaystyle
F^{j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}k^{\prime}}(t,\tau)=w_{jk}w_{j^{\prime}k^{\prime}}w_{j^{\prime\prime}k^{\prime\prime}}\left\\{e^{-{\mathcal{K}}_{jk,j^{\prime}k^{\prime}}(t-\tau)}\right.$
$\displaystyle\
\times\left(f_{jk,k^{\prime\prime}}(t)f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)e^{-{\mathcal{K}}_{jk,j^{\prime\prime}k^{\prime\prime}}(t)-{\mathcal{K}}_{j^{\prime}k^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)}-1\right)$
$\displaystyle\left.\
-f_{jk,k^{\prime\prime}}(t)e^{-{\mathcal{K}}_{jk,j^{\prime\prime}k^{\prime\prime}}(t)}-f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)e^{-{\mathcal{K}}_{j^{\prime}k^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)}+2\right\\}.$

The second bath term in Eq. (132),
$H^{(1),j^{\prime\prime}k^{\prime\prime}}_{jk,j^{\prime}k^{\prime}}(t,\tau)$
defined by Eq. (134), can be expanded further and expressed as follows:

$\displaystyle
H^{(1),j^{\prime\prime}k^{\prime\prime}}_{j,j^{\prime}k^{\prime}}(t,\tau)=\langle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{j"}^{\dagger}\rangle$
$\displaystyle\hskip 21.68121pt-
w_{j^{\prime}k^{\prime}}\langle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip 21.68121pt-w_{j^{\prime\prime}k^{\prime\prime}}\langle
D_{j}(t)(\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)-w_{j^{\prime}k^{\prime}})\rangle.$
(140)

The first term in the above expression can be shown to be

$\displaystyle\langle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle=f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)\langle
D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle+f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)h_{j,k^{\prime\prime}}(t)\langle\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle,$
(141)

where the following identities have been used,

$\displaystyle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)$
$\displaystyle=\theta_{k^{\prime\prime}}D_{j}(t)\theta_{k^{\prime\prime}}^{\dagger}\theta_{k^{\prime\prime}}\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{k^{\prime\prime}}^{\dagger}\theta_{k^{\prime\prime}},$
(142)
$\displaystyle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{k^{\prime\prime}}^{\dagger}=D_{j}(t)+h_{j,k^{\prime\prime}}(t),$
(143)
$\displaystyle\theta_{k^{\prime\prime}}\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{k^{\prime\prime}}^{\dagger}=f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau),$
(144)

along with the definitions of Eqs. (59) and (60). In Eq. (141), the first bath
average on the right hand side can be calculated employing identities similar
to those leading to Eq. (120). The resulting expression is as follows:

$\displaystyle\langle
D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle=({\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t-\tau)+{\mathcal{M}}_{j,j^{\prime\prime}k^{\prime\prime}}(t))w_{j^{\prime}k^{\prime}}w_{j^{\prime\prime}k^{\prime\prime}}e^{-{\mathcal{K}}_{j^{\prime}k^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)}.$

Combining this with the identity given by Eq. (104) (with the replacement of
$t\rightarrow\tau$ and $\tau\rightarrow 0$), one can show that Eq. (141) can
be expressed as follows:

$\displaystyle\langle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{j^{\prime}}^{\dagger}(\tau)\theta_{k^{\prime}}(\tau)\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle=f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)\left({\mathcal{M}}_{j,j^{\prime}k^{\prime}}(t-\tau)+{\mathcal{M}}_{j,j^{\prime\prime}k^{\prime\prime}}(t)+h_{j,k^{\prime\prime}}(t)\right)$
$\displaystyle\hskip 21.68121pt\times
w_{j^{\prime}k^{\prime}}w_{j^{\prime\prime}k^{\prime\prime}}e^{-{\mathcal{K}}_{j^{\prime}k^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)}.$
(146)

The second bath average term on the right hand side of Eq. (140) can be shown
to be

$\langle\theta_{k^{\prime\prime}}D_{j}(t)\theta_{j^{\prime\prime}}^{\dagger}\rangle=({\mathcal{M}}_{j,j^{\prime\prime}k^{\prime\prime}}(t)+h_{j,k^{\prime\prime}}(t))w_{j^{\prime\prime}k^{\prime\prime}},$
(147)

where Eqs. (120) and (143) along with the definition of Eq. (41) have been
used.

For the last term on the righthand side of Eq. (140), Eq. (115) can be
employed. As a result, Eq. (140) can be expressed as

$\displaystyle
H^{(1),j^{\prime\prime}k^{\prime\prime}}_{j,j^{\prime}k^{\prime}}(t,\tau)=\left(f_{j^{\prime}k^{\prime},k^{\prime\prime}}(\tau)e^{-{\mathcal{K}}_{j^{\prime}k^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)}-1\right)$
$\displaystyle\times
w_{j^{\prime}k^{\prime}}w_{j^{\prime\prime}k^{\prime\prime}}\left(M_{j,j^{\prime}k^{\prime}}(t-\tau)+M_{j,j^{\prime\prime}k^{\prime\prime}}(t)+h_{j,k^{\prime\prime}}(t)\right).$
(148)

Following a similar step, Eq. (135) can be shown to have the following
expression:

$\displaystyle
H^{(2),j^{\prime\prime}k^{\prime\prime}}_{j^{\prime},jk}(t,\tau)^{*}=\left(f_{jk,k^{\prime\prime}}(t)e^{-{\mathcal{K}}_{jk,j^{\prime\prime}k^{\prime\prime}}(t)}-1\right)$
$\displaystyle\times
w_{jk}w_{j^{\prime\prime}k^{\prime\prime}}\left(M_{j^{\prime},kj}(t-\tau)+M_{j^{\prime},k^{\prime\prime}j^{\prime\prime}}^{*}(\tau)+h_{j^{\prime},j^{\prime\prime}}(\tau)\right).$
(149)

Finally, Eq. (136) can be expressed as

$\displaystyle L^{j^{\prime\prime}k^{\prime\prime}}_{j,j^{\prime}}(t,\tau)$
$\displaystyle=$
$\displaystyle\langle\theta_{k^{\prime\prime}}D_{j}(t)D_{j^{\prime}}(\tau)\theta_{j^{\prime\prime}}^{\dagger}\rangle$
(150) $\displaystyle-w_{j^{\prime\prime}k^{\prime\prime}}\langle
D_{j}(t)D_{j^{\prime}}(\tau)\rangle\ .$

In the above expression, the first term on the right hand side can be
expressed as

$\displaystyle\langle\theta_{k^{\prime\prime}}D_{j}(t)D_{j^{\prime}}(\tau)\theta_{j^{\prime\prime}}^{\dagger}\rangle=\langle
D_{j}(t)D_{j^{\prime}}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip 14.45377pt+h_{j,k^{\prime\prime}}(t)\langle
D_{j^{\prime}}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip 14.45377pt+h_{j^{\prime},k^{\prime\prime}}(\tau)\langle
D_{j}(t)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip
14.45377pt+h_{j,k^{\prime\prime}}(t)h_{j^{\prime},k^{\prime\prime}}(\tau)w_{j^{\prime\prime}k^{\prime\prime}},$
(151)

where the first term can be calculated as follows:

$\displaystyle\langle
D_{j}(t)D_{j^{\prime}}(\tau)\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rho_{b}\rangle$
$\displaystyle=\sum_{n}\sum_{n^{\prime}}\hbar\omega_{n}g_{n,j}(1-W_{h}(\omega_{n}))\hbar\omega_{n^{\prime}}g_{n^{\prime},j^{\prime}}(1-W_{h}(\omega_{n^{\prime}}))$
$\displaystyle\hskip
28.90755pt\times\left\\{e^{-i(\omega_{n}t+\omega_{n^{\prime}}\tau)}\langle
b_{n}b_{n^{\prime}}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle\right.$
$\displaystyle\hskip
43.36243pt+e^{-i(\omega_{n}t-\omega_{n^{\prime}}\tau)}\langle
b_{n}b_{n^{\prime}}^{\dagger}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip
43.36243pt+e^{i(\omega_{n}t-\omega_{n^{\prime}}\tau)}\langle
b_{n}^{\dagger}b_{n^{\prime}}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle$
$\displaystyle\hskip
43.36243pt\left.+e^{-i(\omega_{n}t+\omega_{n^{\prime}}\tau)}\langle
b_{n}^{\dagger}b_{n^{\prime}}^{\dagger}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle\right\\}.$
(152)

In the above expression, the four thermal averages can be calculated
explicitly and can be shown to be

$\displaystyle\langle
b_{n}b_{n^{\prime}}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle=w_{j^{\prime\prime}k^{\prime\prime}}\frac{\delta
g_{n,j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n})}{(1-e^{-\beta\hbar\omega_{n}})}$
$\displaystyle\hskip 36.135pt\times\frac{\delta
g_{n^{\prime},j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n^{\prime}})}{(1-e^{-\beta\hbar\omega_{n^{\prime}}})},$
(153) $\displaystyle\langle
b_{n}b_{n^{\prime}}^{\dagger}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle=w_{j^{\prime\prime}k^{\prime\prime}}\left(\frac{1}{1-e^{-\beta\hbar\omega_{n}}}\delta_{nn^{\prime}}\right.$
$\displaystyle\left.-\frac{\delta
g_{n,j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n})}{(1-e^{-\beta\hbar\omega_{n}})}\frac{\delta
g_{n^{\prime},j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n^{\prime}})}{(1-e^{-\beta\hbar\omega_{n^{\prime}}})}e^{-\beta\hbar\omega_{n^{\prime}}}\right),$
(154) $\displaystyle\langle
b_{n}^{\dagger}b_{n^{\prime}}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle=w_{j^{\prime\prime}k^{\prime\prime}}\left(\frac{e^{-\beta\hbar\omega_{n}}}{1-e^{-\beta\hbar\omega_{n}}}\delta_{nn^{\prime}}\right.$
$\displaystyle\left.-\frac{\delta
g_{n,j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n})}{(1-e^{-\beta\hbar\omega_{n}})}e^{-\beta\hbar\omega_{n}}\frac{\delta
g_{n^{\prime},j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n^{\prime}})}{(1-e^{-\beta\hbar\omega_{n^{\prime}}})}\right),$
(155) $\displaystyle\langle
b_{n}^{\dagger}b_{n^{\prime}}^{\dagger}\theta_{k^{\prime\prime}}\theta_{j^{\prime\prime}}^{\dagger}\rangle=w_{j^{\prime\prime}k^{\prime\prime}}\frac{\delta
g_{n,j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n})}{(1-e^{-\beta\hbar\omega_{n}})}e^{-\beta\hbar\omega_{n}}$
$\displaystyle\hskip 36.135pt\times\frac{\delta
g_{n^{\prime},j^{\prime\prime}k^{\prime\prime}}W_{h}(\omega_{n^{\prime}})}{(1-e^{-\beta\hbar\omega_{n^{\prime}}})}e^{-\beta\hbar\omega_{n^{\prime}}}.$
(156)

Inserting the above expressions into Eq. (152), one can shown that

$\displaystyle
L^{j^{\prime\prime}k^{\prime\prime}}_{j,j^{\prime},j^{\prime\prime}k^{\prime\prime}}(t,\tau)$
$\displaystyle=$ $\displaystyle
w_{j^{\prime\prime}k^{\prime\prime}}({\mathcal{M}}_{j,j^{\prime\prime}k^{\prime\prime}}(t)+h_{j,k^{\prime\prime}}(t))$
(157)
$\displaystyle\times({\mathcal{M}}_{j^{\prime},j^{\prime\prime}k^{\prime\prime}}(\tau)+h_{j^{\prime},k^{\prime\prime}}(\tau)).$

For all the bath correlation functions calculated above, let us define the
following time integrals:

$\displaystyle\tilde{F}_{jk,j^{\prime}k^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t;{\mathcal{E}})=\int_{0}^{t}d\tau
e^{i{\mathcal{E}}(t-\tau)/\hbar}F_{jk,j^{\prime}k^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t,\tau),$
(158)
$\displaystyle\tilde{H}_{j,j^{\prime}k^{\prime}}^{(1),j^{\prime\prime}k^{\prime\prime}}(t;{\mathcal{E}})=\int_{0}^{t}d\tau
e^{i{\mathcal{E}}(t-\tau)/\hbar}H_{j,j^{\prime}k^{\prime}}^{(1),j^{\prime\prime}k^{\prime\prime}}(t,\tau),$
(159)
$\displaystyle\tilde{H}_{jk,j^{\prime}}^{(2),j^{\prime\prime}k^{\prime\prime}}(t;{\mathcal{E}})=\int_{0}^{t}d\tau
e^{i{\mathcal{E}}(t-\tau)/\hbar}H_{jk,j^{\prime}}^{(2),j^{\prime\prime}k^{\prime\prime}}(t,\tau),$
(160)
$\displaystyle\tilde{L}_{j,j^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t;{\mathcal{E}})=\int_{0}^{t}d\tau
e^{i{\mathcal{E}}(t-\tau)/\hbar}L_{j,j^{\prime}}^{j^{\prime\prime}k^{\prime\prime}}(t,\tau).$
(161)

The above time correlation functions are useful for expressing the second
order inhomogeneous terms in the basis of eigenstates of $\tilde{H}_{s}$.

REFERENCES

## References

  * Landau and Pekar (2008) L. D. Landau and S. I. Pekar, “Effective mass of a polaron,” Ukr. J. Phys. 53, 71–74 (2008).
  * Fröhlich (1954) H. Fröhlich, “Electrons in lattice fields,” Adv. Phys. 3, 325 (1954).
  * Holstein (1959a) T. Holstein, “Studies of polaron motion: Part 1. the molecular-crystal model,” Ann. Phys. 8, 325 (1959a).
  * Holstein (1959b) T. Holstein, “Studies of polaron motion: Part 2. the “small” polaron,” Ann. Phys. 8, 343 (1959b).
  * Holstein (1978) T. R. Holstein, Philos. Mag. B 37, 49 (1978).
  * Rackovsky and Silbey (1973) S. Rackovsky and R. Silbey, “Electronic energy transfer in impure solids i. two molecules embeded in a lattice,” Mol. Phys. 25, 61 (1973).
  * Jackson and Silbey (1983) B. Jackson and R. Silbey, ““On the calculation of transfer rate between impurity states in solids”,” J. Chem. Phys. 78, 4193 (1983).
  * Silbey and Harris (1984) 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 (1984).
  * Harris and Silbey (1985a) R. A. Harris and R. Silbey, “Variational calculation of the tunneling system interacting with a heat bath. ii. dynamics of an asymmetric tunneling system,” J. Chem. Phys. 83, 1069 (1985a).
  * Harris and Silbey (1985b) R. A. Harris and R. Silbey, J. Chem. Phys. 83, 1069 (1985b).
  * Nitzan (2006) A. Nitzan,  _Chemical Dynamics in Condensed Phases_ (Oxford University Press, Oxford, 2006).
  * Cheng and Silbey (2008) Y. C. Cheng and R. J. Silbey, J. Chem. Phys. 128, 114713 (2008).
  * Jang  _et al._ (2008) S. Jang, Y.-C. Cheng, D. R. Reichman, and J. D. Eaves, “Theory of coherent resonance energy transfer,” J. Chem. Phys. 129, 101104 (2008).
  * Jang (2009) S. Jang, “Theory of coherent resonance energy transfer for coherent initial condition,” J. Chem. Phys. 131, 164101 (2009).
  * Jang (2011) S. Jang, “Theory of multichromophoric coherent resonance energy transfer: A polaronic quantum master equation approach,” J. Chem. Phys. 135, 034105 (2011).
  * Nazir (2009) A. Nazir, “Correlation-dependent coherent to incoherent transitions in resonant energy transfer dynamics,” Phys. Rev. Lett. 103, 146404 (2009).
  * McCutcheon and Nazir (2011) D. P. S. McCutcheon and A. Nazir, “Consistent treatment of coherent and incoherent energy transfer dynamics using a variational master equation,” J. Chem. Phys. 135, 114501 (2011).
  * Yang, Devi, and Jang (2012) L. Yang, M. Devi, and S. Jang, “Polaronic quantum master equation theory of inelastic and coherent resonance energy transfer for soft systems,” J. Chem. Phys. 137, 024101 (2012).
  * Pouthier (2013) V. Pouthier, “The reduced dynamics of an exciton coupled to a phonon bath: A new approach combining the lang-firsov transformation and the perturbation theory,” J. Chem. Phys. 138, 044108 (2013).
  * Chen  _et al._ (2011) D. Chen, J. He, H. Zhang, and Y. Zhao, “On the munn silbey approach to polaron transport with off-diagonal coupling and temperature-dependent canonical transformations,” J. Phys. Chem. B 115, 5312 – 5321 (2011).
  * Zhao  _et al._ (2012) Y. Zhao, B. Luo, Y. Zhang, and J.Ye, “Dynamics of a holstein polaron with off-diagonal coupling,” J. Chem. Phys. 137, 084113 (2012).
  * Chorosajev  _et al._ (2014) V. Chorosajev, A. Gelzinis, L. Valkunas, and D. Abramavicius, “Dynamics of exciton-polaron transition in molecular assemblies: The variational approach,” J. Chem. Phys. 140, 244108 (2014).
  * Hamm and Tsironis (2008) P. Hamm and G. P. Tsironis, “Barrier crossing to the small holstein polaron regime,” Phys. Rev. B 78, 092301 (2008).
  * Lee, Moix, and Cao (2015) C. K. Lee, J. Moix, and J. Cao, J. Chem. Phys. 142, 164103 (2015).
  * Xu and Cao (2016) D. Xu and J. Cao, “Non-canonical distribution and non-equilibrium transport beyond weak system-bath coupling regime: A polaron transformation approach,” Fron. Phys. 11, 110308 (2016).
  * Wang and Zhao (2020) Y.-C. Wang and Y. Zhao, “Variational polaron transformation approach toward the calculation of thermopower in organic crystals,” Phys. Rev. B 101, 075205 (2020).
  * Balzer  _et al._ (2021) D. Balzer, T. J. A. M. Smolders, D. Blyth, S. N. Hood, and I. Kassal, “Delocalized kinetic monte carlo for simulating delocalization-enhanced charge and exciton transport in disordered materials,” Chem. Sci. 12, 2276–2285 (2021).
  * Pollock  _et al._ (2013) F. A. Pollock, D. P. S. McCutcheon, B. W. Lovett, E. M. Gauger, and A. Nazir, “A multi-site variatonal master equation approach to dissipative energy transfer,” New J. Phys. 15, 075018 (2013).
  * Teh, Jin, and Cheng (2019) H.-H. Teh, B.-Y. Jin, and Y.-C. Cheng, “Frozen-mode small polaron quantum master equation with variational bound for excitation energy transfer in molecular aggregates,” J. Chem. Phys. 150, 224110 (2019).
  * Caldeira and Leggett (1983) A. O. Caldeira and A. J. Leggett, Ann. Phys. 149, 374 (1983).
  * Weiss (1993) U. Weiss,  _Series in Modern Condensed Matter Physics Vol. 2 : Quantum Dissipative Systems_ (World Scientific, Singapore, 1993).
  * Kenkre and Reineker (1982) V. M. Kenkre and P. Reineker,  _Exciton Dynamics in Molecular Crystals and Aggregates_ (Springer, Berlin, 1982).
  * May and Kühn (2011) V. May and O. Kühn,  _Charge and Energy Transfer Dynamics in Molecular Systems_ (Wiley-VCH, Weinheim, Germany, 2011).
  * Jang (2020) S. J. Jang,  _Dynamics of Molecular Excitons (Nanophotonics Series)_ (Elsevier, Amsterdam, 2020).
  * Coropceanu  _et al._ (2007) V. Coropceanu, J. Cornil, D. A. da Silva Fihlo, Y. Olivier, R. Silbey, and J.-L. Brédas, “Charge transport in organic semiconductors,” Chem. Rev. 107, 926–952 (2007).
  * Jang, Cao, and Silbey (2002) S. Jang, J. Cao, and R. J. Silbey, “Fourth order quantum master equation and its markovian bath limit,” J. Chem. Phys. 116, 2705 (2002).
  * Shibata and Arimitsu (1980) F. Shibata and T. Arimitsu, “Expansion formulas in nonequilibrium statistical mechanics,” J. Phys. Soc. Jpn 49, 891 (1980).
  * Lai and Geva (2021) Y. Lai and E. Geva, “On simulating the dynamics of electronic populations and coherences via quantum master equations based on treating off-diagonal electronic coupling terms as a small perturbation,” J. Chem. Phys. 155, 204101 (2021).