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

Jang, Seogjoo J.

doi:10.1063/5.0106546

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, the 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. A formal framework to address this issue is developed in the present work 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 second 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 of a two-level system coupled to a bath of harmonic oscillators, with test calculations focused on those due to homogeneous relaxation terms, demonstrate the feasibility and the utility of the present approach.

I. INTRODUCTION

Polaron transformation (PT)^1–31 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 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 within the limit of strong couplings to the bath. As long as physical observables of interest remain unentangled with the bath through the PT, a quantum master equation (QME) derived by projecting out the bath in the polaron picture can be used for calculating such observables. Indeed, utilizing this fact has led to a general approach called polaron-transformed QME (PQME).^14–22

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 PQMEs^14–17 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 advantageous because it invokes an over-relaxation of some slow modes, which can cause premature suppression of important coherence terms, in particular, within the second order time-local approximation. Variational PQME^18,19,21,22 that combines a variational ansatz^6,9–11 with PQME and, more recently, frozen mode PQME³² approaches have been developed to address this issue. This work develops a formulation that can put these works^18,21,32 into a broader context and can ameliorate the issue of over-relaxation by extending the PQME for partial PTs of a fairly general kind, where high frequency bath modes are transformed preferentially. It is shown here that a general closed form expression of a second order time-local PQME including all the inhomogeneous terms can be derived as in the case of the full PT,¹⁵ 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 the utility of the p-PQME approach.

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

II. THEORY

A. Hamiltonian and partial polaron transformation

Let us consider a quantum system consisting of N coupled quantum states, |j⟩s, 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 〉 〈 j | + \sum_{j, k = 1}^{N} J_{j k} | j 〉 〈 k |,

(1)

where E_j is the energy of the site local state |j⟩ and J_jk, for j ≠ k, is the electronic coupling between states |j⟩ and |k⟩, which 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 include 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⟩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_{s b},

(2)

where

H_{b} = \sum_{n} ℏ ω_{n} (b_{n}^{†} b_{n} + \frac{1}{2}),

(3)

H_{s b} = \sum_{j = 1}^{N} \sum_{n} ℏ ω_{n} g_{n, j} (b_{n} + b_{n}^{†}) | j 〉 〈 j | .

(4)

In the above expressions, ω_n and b_n $(b_{n}^{†})$ 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⟩. The nature of this bath can be characterized collectively by introducing the following bath spectral density:

J_{j j^{'}} (ω) = π ℏ \sum_{n} δ (ω - ω_{n}) ω_{n}^{2} g_{n, j} g_{n, j^{'}} .

(5)

The total Hamiltonian, Eq. (2), is a straightforward generalization of the spin-boson Hamiltonian^33,34 and has been used widely for molecular excitons^35–37 and charge transport dynamics.^34,36,38 The full information on the combined system and bath degrees of freedom, in general, requires determining the total density operator ρ(t), which is governed by the following quantum Liouville equation:

\frac{d}{d t} ρ (t) = - i L ρ (t) \equiv - \frac{i}{ℏ} [H, ρ (t)] .

(6)

Due to a 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 are 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, namely, when Eq. (5) behaves linearly or sub-linearly in the low frequency limit, the sum of Huang–Rhys factors for the bath degrees of freedom diverges.³⁴ 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(ω) with the following limiting behavior:

W_{h} (ω) = \{\begin{align} O (ω^{α}) & for ω \to 0, \\ 1 & for ω \to \infty . \end{align})

(7)

The scaling behavior for ω → 0 in the above equation, for a sufficiently large value of α, 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, this means that α ≥ 1 at least. If the bath spectral density has sub-Ohmic components, the lower bound for α 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} (ω_{n}) (b_{n}^{†} - b_{n}) | j 〉 〈 j | .

(8)

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

\begin{matrix} {\tilde{H}}_{s} & = & e^{G} H_{s} e^{- G} \\ = & \sum_{j = 1}^{N} E_{j} | j 〉 〈 j | + \sum_{j, k = 1}^{N} J_{j k} θ_{j}^{†} θ_{k} | j 〉 〈 k |, \end{matrix}

(9)

where

θ_{j} = e^{- \sum_{n} g_{n, j} W_{h} (ω_{n}) (b_{n}^{†} - b_{n})} .

(10)

On the other hand, it is straightforward to show that

\begin{align} e^{G} (H_{b} + H_{s b}) e^{- G} & = H_{b} + \sum_{j = 1}^{N} \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) (b_{n} + b_{n}^{†}) \times | j 〉 〈 j | \\ - \sum_{j = 1}^{N} \sum_{n} ℏ ω_{n} g_{n, j}^{2} W_{h} (ω_{n}) \times (2 - W_{h} (ω_{n})) | j 〉 〈 j | . \end{align}

(11)

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

\tilde{H} = e^{G} H e^{- G} = {\tilde{H}}_{s}^{p} + {\tilde{H}}_{s b} + H_{b},

(12)

where

{\tilde{H}}_{s}^{p} = \sum_{j = 1}^{N} {\tilde{E}}_{j} | j 〉 〈 j | .

(13)

In the above expression, ${\tilde{E}}_{j}$ is a partially renormalized energy for site j given by

{\tilde{E}}_{j} = E_{j} - λ_{j},

(14)

with λ_j being the corresponding reorganization energy defined as

\begin{matrix} λ_{j} & = & \sum_{n} ℏ ω_{n} g_{n, j}^{2} W_{h} (ω_{n}) (2 - W_{h} (ω_{n})) \\ = & \frac{1}{π} \int_{0}^{\infty} d ω \frac{J_{j j} (ω)}{ω} W_{h} (ω) (2 - W_{h} (ω)) . \end{matrix}

(15)

In the second equality of the above equation, the definition of the bath spectral density, Eq. (5), was used.

The second term in Eq. (12), ${\tilde{H}}_{s b}$ ⁠, is a partially renormalized system–bath interaction Hamiltonian given by

{\tilde{H}}_{s b} = \sum_{j, k = 1}^{N} J_{j k} θ_{j}^{†} θ_{k} | j 〉 〈 k | + \sum_{j = 1}^{N} \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) (b_{n} + b_{n}^{†}) | j 〉 〈 j |,

(16)

where $θ_{j}^{†} θ_{k}$ can be expressed as

θ_{j}^{†} θ_{k} = e^{\sum_{n} δ g_{n, j k} W_{h} (ω_{n}) (b_{n}^{†} - b_{n})},

(17)

with δg_n,jk = (g_n,j − g_n,k).

B. 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 Sec. II A, is provided below.

1. Quantum Liouville equation in the polaron and interaction picture

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

\tilde{H} = {\tilde{H}}_{0} + {\tilde{H}}_{1},

(18)

where

{\tilde{H}}_{0} = {\tilde{H}}_{s}^{p} + 〈 {\tilde{H}}_{s b} 〉_{b} + H_{b} \equiv {\tilde{H}}_{0, s} + H_{b},

(19)

{\tilde{H}}_{1} = {\tilde{H}}_{s b} - 〈 {\tilde{H}}_{s b} 〉_{b} .

(20)

In above expressions, $〈 {\tilde{H}}_{s b} 〉_{b} = {Tr}_{b} {{\tilde{H}}_{s b} ρ_{b}}$ with $ρ_{b} = e^{- β H_{b}} / {Tr}_{b} {e^{- β H_{b}}}$ ⁠. This term represents the average system–bath interaction in the partial polaron picture.

In Eq. (19), the system part of the new zeroth order Hamiltonian, ${\tilde{H}}_{0, s} = {\tilde{H}}_{s}^{p} + 〈 {\tilde{H}}_{s}^{c} 〉_{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 〉 〈 j | + \sum_{j, k = 1}^{N} {\tilde{J}}_{j k} | j 〉 〈 k |,

(21)

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

\begin{matrix} w_{j k} & = & 〈 θ_{j}^{†} θ_{k} 〉 = 〈 θ_{k}^{†} θ_{j} 〉 \\ = & e^{- \sum_{n} \coth (β ℏ ω_{n} / 2) δ g_{n, j k}^{2} W_{h} {(ω_{n})}^{2} / 2} . \end{matrix}

(22)

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 given that the weighting function satisfies Eq. (7).

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

{\tilde{H}}_{1} = \sum_{j, k = 1}^{N} {\tilde{B}}_{j k} | j 〉 〈 k |,

(23)

where

{\tilde{B}}_{j k} = J_{j k} (θ_{j}^{†} θ_{k} - w_{j k}) + δ_{j k} D_{j} .

(24)

In the above expressions δ_jk is the Kronecker-delta symbol and

D_{j} = \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) (b_{n} + b_{n}^{†}) .

(25)

Thus, the bath operator ${\tilde{B}}_{j k}$ given by Eq. (24) is a sum of the renormalized system–bath interaction term (relative to its average) due to partial PT (for j ≠ 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

\begin{matrix} {\tilde{H}}_{1, I} (t) & \equiv & e^{i {\tilde{H}}_{0} t / ℏ} {\tilde{H}}_{1} e^{- i {\tilde{H}}_{0} t / ℏ} \\ = & \sum_{j, k = 1}^{N} {\tilde{B}}_{j k} (t) T_{j k} (t), \end{matrix}

(26)

where

T_{j k} (t) = e^{i {\tilde{H}}_{0, s} t / ℏ} | j 〉 〈 k | e^{- i {\tilde{H}}_{0, s} t / ℏ},

(27)

\begin{matrix} {\tilde{B}}_{j k} (t) & = & e^{i H_{b} t / ℏ} {\tilde{B}}_{j k} e^{- i H_{b} t / ℏ} \\ = & J_{j k} (θ_{j}^{†} (t) θ_{k} (t) - w_{j k}) + δ_{j k} D_{j} (t) . \end{matrix}

(28)

In the above expression,

θ_{j}^{†} (t) θ_{k} (t) = e^{\sum_{n} δ g_{n, j k} W_{h} (ω_{n}) (b_{n}^{†} e^{i ω_{n} t} - b_{n} e^{- i ω_{n} t})},

(29)

D_{j} (t) = \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) (b_{n} e^{- i ω_{n} t} + b_{n}^{†} e^{i ω_{n} t}) .

(30)

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

\frac{d}{d t} {\tilde{ρ}}_{I} (t) = - i {\tilde{L}}_{1, I} (t) {\tilde{ρ}}_{I} (t) \equiv - \frac{i}{ℏ} [{\tilde{H}}_{1, I} (t), {\tilde{ρ}}_{I} (t)],

(31)

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

2. Quantum master equation for the reduced density operator

Taking trace of ${\tilde{ρ}}_{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{σ}}_{I} (t) \equiv {Tr}_{b} \{{\tilde{ρ}}_{I} (t)\} .

(32)

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{σ}}_{I} (t)$ employing the standard projection operator technique^37,39 for a well-known projection operator $P (\cdot) \equiv ρ_{b} {Tr}_{b} {(\cdot)}$ as follows:

\begin{align} \frac{d}{d t} {\tilde{σ}}_{I} (t) & = - \int_{0}^{t} d τ {Tr}_{b} \{{\tilde{L}}_{1, I} (t) e_{(+)}^{- i \int_{τ}^{t} d τ^{'} Q {\tilde{L}}_{1, I} (τ^{'})} Q {\tilde{L}}_{1, I} (τ) ρ_{b}\} \times {\tilde{σ}}_{I} (τ) \\ - i T r_{b} \{{\tilde{L}}_{1, I} (t) e_{(+)}^{- i \int_{0}^{t} d τ Q {\tilde{L}}_{1, I} (τ)} Q \tilde{ρ} (0)\}, \end{align}

(33)

where e₍₊₎ denotes the exponential operator with chronological time ordering, $Q = 1 - P$ and

Q \tilde{ρ} (0) = e^{G} ρ (0) e^{- G} - ρ_{b} {Tr}_{b} {e^{G} ρ (0) e^{- G}} .

(34)

Alternatively, replacing ${\tilde{ρ}}_{I} (τ)$ with the backpropagation of ${\tilde{ρ}}_{I} (t)$ from t to τ within the projection operator formalism,^37,40 one can obtain the following formally exact time-local equation:

\begin{align} \frac{d}{d t} {\tilde{σ}}_{I} (t) & = - \int_{0}^{t} d τ {Tr}_{b} \{{\tilde{L}}_{1, I} (t) {(1 + i Γ_{1, I} (t))}^{- 1} e_{(+)}^{- i \int_{τ}^{t} d τ^{'} Q {\tilde{L}}_{1, I} (τ^{'})}) \\ \times (Q {\tilde{L}}_{1, I} (τ) P e_{(-)}^{i \int_{τ}^{t} d τ^{'} {\tilde{L}}_{1, I} (τ^{'})} ρ_{b}\} {\tilde{σ}}_{I} (t) \\ - i P {\tilde{L}}_{1, I} (t) {(1 + i Γ_{1, I} (t))}^{- 1} e_{(+)}^{- i \int_{0}^{t} d τ Q {\tilde{L}}_{1, I} (τ)} Q \tilde{ρ} (0), \end{align}

(35)

where e₍₋₎ denotes the exponential operator with anti-chronological time ordering and

Γ_{1, I} (t) = \int_{0}^{t} d τ e_{(+)}^{- i \int_{τ}^{t} d τ^{'} Q {\tilde{L}}_{1, I} (τ^{'})} Q {\tilde{L}}_{1, I} (τ) P e_{(-)}^{i \int_{τ}^{t} d τ^{'} {\tilde{L}}_{1, I} (τ^{'})} .

(36)

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

\frac{d}{d t} {\tilde{σ}}_{I} (t) = - R (t) {\tilde{σ}}_{I} (t) + I (t),

(37)

where

R (t) = \int_{0}^{t} d τ {Tr}_{b} {{\tilde{L}}_{1, I} (t) {\tilde{L}}_{1, I} (τ) ρ_{b}},

(38)

\begin{align} I (t) & = I^{(1)} (t) + I^{(2)} (t) \\ = - i {Tr}_{b} {{\tilde{L}}_{1, I} (t) Q \tilde{ρ} (0)} \\ - \int_{0}^{t} d τ {Tr}_{b} {{\tilde{L}}_{1, I} (t) {\tilde{L}}_{1, I} (τ) Q \tilde{ρ} (0)} . \end{align}

(39)

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

Equation (37) is the second order time local p-PQME expressed in Liouville space. In all previous studies and in the present paper, this time local form has been chosen due to its convenience. However, a time nonlocal second order expression can also be derived directly from Eq. (33), and its performance compared with 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, considering that the second order time-nonlocal PQME is equivalent to the non-interacting blip approximation,^9,41–43 which accounts for a significant contribution of coherent dynamics even for the Ohmic bath, it is likely that the time non-local equation becomes more reliable as the bath becomes sluggish. In addition, a recent work⁴⁴ also provides examples of the 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 a more comprehensive assessment of the two approaches. With this point clear, the rest of this section provides detailed expressions for the relaxation superoperator, Eq. (38), and the inhomogeneous term, Eq. (39), in the Hilbert space.

3. Homogeneous terms of the second order time local p-PQME

The Hilbert space expression for $R (t) {\tilde{σ}}_{I} (t)$ in Eq. (37) can be obtained by employing Eqs. (26)–(28) in Eq. (31) 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:

\begin{matrix} R (t) {\tilde{σ}}_{I} (t) & = & \frac{1}{ℏ^{2}} \sum_{j, k = 1}^{N} \sum_{j^{'}, k^{'} = 1}^{N} \int_{0}^{t} d τ 〈 {\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) 〉 \\ \times [T_{j k} (t), T_{j^{'} k^{'}} (τ) {\tilde{σ}}_{I} (t)] + H . c ., \end{matrix}

(40)

where H.c refers to Hermitian conjugates of all previous terms and $〈 {\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) 〉$ (with subscript b omitted) represents averaging over the equilibrium bath density operator, ρ_b.

Appendix A describes calculation of all the terms constituting $〈 {\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) 〉$ ⁠. When the resulting expressions, Eqs. (A2), (A15), (A21), and (A22), are used in Eq. (A1), it can be expressed as

\begin{matrix} 〈 {\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) 〉 & = & {\tilde{J}}_{j k} {\tilde{J}}_{j^{'} k^{'}} (e^{- K_{j k, j^{'} k^{'}} (t - τ)} - 1) + δ_{j k} {\tilde{J}}_{j^{'} k^{'}} M_{j, j^{'} k^{'}} (t - τ) \\ + δ_{j^{'} k^{'}} {\tilde{J}}_{j k} M_{j^{'}, k j} (t - τ) + δ_{j k} δ_{j^{'} k^{'}} C_{j j^{'}} (t - τ), \end{matrix}

(41)

where

\begin{matrix} K_{j k, j^{'} k^{'}} (t) & = & \sum_{n} δ g_{n, j k} δ g_{n, j^{'} k^{'}} W_{h} {(ω_{n})}^{2} \\ \times (\coth (\frac{β ℏ ω_{n}}{2}) \cos (ω_{n} t) - i \sin (ω_{n} t)), \end{matrix}

(42)

\begin{matrix} M_{j, j^{'} k^{'}} (t) & = & \sum_{n} ℏ ω_{n} g_{n, j} δ g_{n, j^{'} k^{'}} (1 - W_{h} (ω_{n})) W_{h} (ω_{n}) \\ \times (\cos (ω_{n} t) - i \coth (\frac{β ℏ ω_{n}}{2}) \sin (ω_{n} t)), \end{matrix}

(43)

\begin{matrix} C_{j j^{'}} (t) & = & \sum_{n} ℏ^{2} ω_{n}^{2} g_{n, j} g_{n, j^{'}} {(1 - W_{h} (ω_{n}))}^{2} \\ \times (\coth (\frac{β ℏ ω_{n}}{2}) \cos (ω_{n} t) - i \sin (ω_{n} t)) . \end{matrix}

(44)

Similar expressions as above have also been derived in the context of variational PQME.²¹ Note that the three bath correlation functions defined above satisfy the following symmetry properties:

\begin{matrix} K_{j k, j^{'} k^{'}} (t) & = & K_{k j, k^{'} j^{'}} (t) \\ = & - K_{j k, k^{'} j^{'}} (t) = - K_{k j, j^{'} k^{'}} (t), \end{matrix}

(45)

M_{j, j^{'} k^{'}} (t) = - M_{j, k^{'} j^{'}} (t),

(46)

C_{j j^{'}} (t) = C_{j^{'} j} (t) .

(47)

The three correlation functions, Eqs. (42)–(44), can all be expressed in terms of the bath spectral density, Eq. (5). For more compact expressions of these, let us introduce the following auxiliary bath spectral densities,

\begin{matrix} J_{j, j^{'} k^{'}}^{(1)} (ω) & = & π ℏ \sum_{n} δ (ω - ω_{n}) ω_{n}^{2} g_{n, j} δ g_{n, j^{'} k^{'}} \\ = & J_{j j^{'}} (ω) - J_{j k^{'}} (ω), \end{matrix}

(48)

\begin{matrix} J_{j k, j^{'} k^{'}}^{(2)} (ω) & = & π ℏ \sum_{n} δ (ω - ω_{n}) ω_{n}^{2} δ g_{n, j k} δ g_{n, j^{'} k^{'}} \\ = & J_{j j^{'}} (ω) + J_{k k^{'}} (ω) - J_{j k^{'}} (ω) - J_{k j^{'}} (ω) . \end{matrix}

(49)

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

\begin{matrix} K_{j k, j^{'} k^{'}} (t) & = & \frac{1}{π ℏ} \int_{0}^{\infty} d ω \frac{J_{j k, j^{'} k^{'}}^{(2)} (ω)}{ω^{2}} W_{h} {(ω)}^{2} \\ \times (\coth (\frac{β ℏ ω}{2}) \cos (ω t) - i \sin (ω t)), \end{matrix}

(50)

\begin{matrix} M_{j, j^{'} k^{'}} (t) & = & \frac{1}{π} \int_{0}^{\infty} d ω \frac{J_{j, j^{'} k^{'}}^{(1)} (ω)}{ω} (1 - W_{h} (ω)) W_{h} (ω) \\ \times (\cos (ω t) - i \coth (\frac{β ℏ ω}{2}) \sin (ω t)), \end{matrix}

(51)

\begin{matrix} C_{j j^{'}} (t) & = & \frac{ℏ}{π} \int_{0}^{\infty} d ω J_{j j^{'}} (ω) {(1 - W_{h} (ω))}^{2} \\ \times (\coth (\frac{β ℏ ω}{2}) \cos (ω t) - i \sin (ω t)) . \end{matrix}

(52)

Note also that w_jk defined by Eq. (22) can be expressed as follows:

\begin{matrix} w_{j k} & = & e^{- K_{j k, j k} (0) / 2} \\ = & \exp \{- \frac{1}{2 π ℏ} \int_{0}^{\infty} d ω \frac{J_{j k, j k}^{(2)} (ω)}{ω^{2}} W_{h} {(ω)}^{2} \coth (\frac{β ℏ ω}{2})\} . \end{matrix}

(53)

If W_h(ω) = 1, $M_{j, j^{'} k^{'}} (t)$ and $C_{j j^{'}} (t)$ become zero and the expressions for $K_{j k, j^{'} k^{'}} (t)$ and w_jk reduce to those for the original second order PQME based on the full PT. On the other hand, for W_h(ω) = 0, $K_{j k, j^{'} k^{'}} (t)$ and $C_{j j^{'}} (t)$ become zero and the expression for $C_{j j^{'}} (t)$ reduce to that for a conventional second 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(ω) is a step function. For this case, W_h(ω)(1 − W_h(ω)) = 0 and, thus, $M_{j, j^{'} k^{'}} (t)$ becomes zero. As a result, the relaxation superoperator becomes a simple sum of those due to the p-PT part and untransformed linear coupling terms.

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

ρ (0) = σ (0) ρ_{b},

(54)

where

σ (0) = \sum_{j, k = 1}^{N} σ_{j k} (0) | j 〉 〈 k | .

(55)

Then,

Q \tilde{ρ} (0) = \sum_{j, k = 1}^{N} σ_{j k} (0) T_{j k} (0) δ {\tilde{ρ}}_{b, j k},

(56)

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

δ {\tilde{ρ}}_{b, j k} = θ_{j}^{†} ρ_{b} θ_{k} - w_{j k} ρ_{b},

(57)

with the convention that w_jk = 1 for j = k. Therefore, the first order inhomogeneous term in Eq. (39) can be expressed as

\begin{align} I^{(1)} (t) & = - i T r_{b} {{\tilde{L}}_{1, I} (t) Q \tilde{ρ} (0)} \\ = - \frac{i}{ℏ} \sum_{j, k = 1}^{N} \sum_{j^{'}, k^{'} = 1}^{N} {Tr}_{b} \{{\tilde{B}}_{j k} (t) δ {\tilde{ρ}}_{b, j^{'} k^{'}}\} σ_{j^{'} k^{'}} (0) [T_{j k} (t), T_{j^{'} k^{'}} (0)] . \end{align}

(58)

In the above expression, the trace over the bath can be calculated explicitly employing the definitions for ${\tilde{B}}_{j k} (t)$ and $δ {\tilde{ρ}}_{b, j k}$ ⁠, Eqs. (28) and (57), respectively. Details of this calculation are provided in Appendix B, and the resulting expression can be summed up as

\begin{matrix} {Tr}_{b} {{\tilde{B}}_{j k} (t) δ {\tilde{ρ}}_{b, j^{'} k^{'}}} & = & w_{j^{'} k^{'}} \{{\tilde{J}}_{j k} (e^{- K_{j k, j^{'} k^{'}} (t)} f_{j k, k^{'}} (t) - 1)) \\ + (δ_{j k} (M_{j, j^{'} k^{'}} (t) + h_{j, k^{'}} (t))\}, \end{matrix}

(59)

where

\begin{matrix} f_{j k, k^{'}} (t) & = & \exp \{2 i \sum_{n} g_{n, k^{'}} δ g_{n, j k} W_{h} {(ω_{n})}^{2} \sin (ω_{n} t)\} \\ = & \exp \{\frac{2 i}{π ℏ} \int_{0}^{\infty} d ω \frac{J_{k^{'}, j k}^{(1)} (ω)}{ω^{2}} W_{h} {(ω)}^{2} \sin (ω t)\}, \end{matrix}

(60)

\begin{align} h_{j, k^{'}} (t) & = 2 \sum_{n} ℏ ω_{n} g_{n, j} g_{n, k^{'}} (1 - W_{h} (ω_{n})) W_{h} (ω_{n}) \cos (ω_{n} t) \\ = \frac{2}{π} \int_{0}^{\infty} d ω \frac{J_{j k^{'}} (ω)}{ω} (1 - W_{h} (ω)) W_{h} (ω) \cos (ω t) . \end{align}

(61)

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

\begin{matrix} I^{(2)} (t) & = & - \int_{0}^{t} d τ {Tr}_{b} \{{\tilde{L}}_{1, I} (t) {\tilde{L}}_{1, I} (τ) Q \tilde{ρ} (0)\} \\ = & - \frac{1}{ℏ^{2}} \sum_{j, k = 1}^{N} \sum_{j^{'}, k^{'} = 1}^{N} \sum_{j^{″}, k^{″} = 1}^{N} \\ \times \int_{0}^{t} d τ {Tr}_{b} \{{\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) δ {\tilde{ρ}}_{b, j^{″} k^{″}}\} \\ \times σ_{j^{″} k^{″}} (0) [T_{j k} (t), T_{j^{'} k^{'}} (τ) T_{j^{″} k^{″}} (0)] + H . c ., \end{matrix}

(62)

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 second order PQME for W_h(ω) = 1 and those for the regular second order QME for W_h(ω) = 0. For the case where W_h(ω)(1 − W_h(ω)) = 0, they become independent sums of those due to PT and due to untransformed linear system–bath couplings.

C. 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 pth eigenstate and eigenvalue of ${\tilde{H}}_{0, s}$ as |φ_p⟩ and $E_{p}$ ⁠. Then,

{\tilde{H}}_{0, s} = \sum_{p = 1}^{N} E_{p} | φ_{p} 〉 〈 φ_{p} | .

(63)

The transformation matrix U between the site localized states |j⟩s and the eigenstates |φ_p⟩s can be defined such that U_jp = ⟨j|φ_p⟩. Then,

| j 〉 = \sum_{p = 1}^{N} U_{j p}^{*} | φ_{p} 〉 .

(64)

This transformation can be used to express $T_{j k} (t)$ defined by Eq. (27) in the basis of |φ_p⟩s as follows:

T_{j k} (t) = \sum_{p, q = 1}^{N} U_{j p}^{*} U_{k q} e^{i δ E_{p q} t / ℏ} | φ_{p} 〉 〈 φ_{q} |,

(65)

where $δ E_{p q} = E_{p} - E_{q}$ ⁠. Let us also introduce S_pq(t)s such that

{\tilde{σ}}_{I} (t) = \sum_{p, q = 1}^{N} S_{p q} (t) | φ_{p} 〉 〈 φ_{q} | .

(66)

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

[T_{j k} (t), T_{j^{'} k^{'}} (τ) {\tilde{σ}}_{I} (t)] = \sum_{p, q = 1}^{N} | φ_{p} 〉 〈 φ_{q} | \sum_{p^{'}, q^{'} = 1}^{N} (δ_{q^{'} q} \sum_{r = 1}^{N} U_{j p}^{*} U_{k r} U_{j^{'} r} U_{k^{'} p^{'}} e^{i δ E_{p p^{'}} t / ℏ} e^{i δ E_{p^{'} r} (t - τ) / ℏ} - U_{j q^{'}}^{*} U_{k q} U_{j^{'} p}^{*} U_{k^{'} p^{'}} e^{i (δ E_{p p^{'}} - δ E_{q q^{'}}) t / ℏ} e^{i δ E_{p^{'} p} (t - τ) / ℏ}) S_{p^{'} q^{'}} (t) .

(67)

The above expression, in combination with Eq. (41), can be used to express $R (t) {\tilde{σ}}_{I} (t)$ in the basis of |φ_p⟩s. For a more compact expression, let us introduce

W_{j k, j^{'} k^{'}}^{p q} (t) = {\tilde{J}}_{j k} {\tilde{J}}_{j^{'} k^{'}} \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} (e^{- K_{j k, j^{'} k^{'}} (t - τ)} - 1),

(68)

Y_{j, j^{'} k^{'}}^{p q} (t) = {\tilde{J}}_{j^{'} k^{'}} \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} M_{j, j^{'} k^{'}} (t - τ),

(69)

X_{j j^{'}}^{p q} (t) = \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} C_{j j^{'}} (t - τ) .

(70)

Then,

R (t) {\tilde{σ}}_{I} (t) = \sum_{p, q = 1}^{N} \sum_{p^{'}, q^{'} = 1}^{N} | φ_{p} 〉 〈 φ_{q} | R_{p q}^{p^{'} q^{'}} (t) S_{p^{'} q^{'}} (t),

(71)

where

\begin{matrix} R_{p q}^{p^{'} q^{'}} (t) & = & \frac{1}{ℏ^{2}} \sum_{j, k = 1}^{N} \sum_{j^{'}, k^{'} = 1}^{N} \{δ_{q^{'} q} \sum_{r = 1}^{N} U_{j p}^{*} U_{k r} U_{j^{'} r}^{*} U_{k^{'} p^{'}} e^{i δ E_{p p^{'}} t / ℏ}) (W_{j k, j^{'} k^{'}}^{p^{'} r} (t) + δ_{j k} Y_{j, j^{'} k^{'}}^{p^{'} r} (t) + δ_{j^{'} k^{'}} Y_{j^{'}, k j}^{p^{'} r} (t) + δ_{j k} δ_{j^{'} k^{'}} X_{j j^{'}}^{p^{'} r} (t)) \\ - U_{j q^{'}}^{*} U_{k q} U_{j^{'} p}^{*} U_{k^{'} p^{'}} e^{i (δ E_{p p^{'}} - δ E_{q q^{'}}) t / ℏ} (W_{j k, j^{'} k^{'}}^{p^{'} p} (t) + δ_{j k} Y_{j, j^{'} k^{'}}^{p^{'} p} (t) + (δ_{j^{'} k^{'}} Y_{j^{'}, k j}^{p^{'} p} (t) + δ_{j k} δ_{j^{'} k^{'}} X_{j j^{'}}^{p^{'} p} (t))\} + [c . c ., p \leftrightarrow q, p^{'} \leftrightarrow q^{'}] . \end{matrix}

(72)

In the above expression, the last line represents complex conjugates along with the interchange of indices, p ↔ q and p′ ↔ q′, of all previous terms.

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

\begin{matrix} [T_{j k} (t), T_{j^{'} k^{'}} (0)] & = & \sum_{p, q = 1}^{N} | φ_{q} 〉 〈 φ_{q} | \sum_{r} (U_{j p}^{*} U_{k r} U_{j^{'} r}^{*} U_{k^{'} q} e^{i δ E_{p r} t / ℏ}) \\ (- U_{j r}^{*} U_{k q} U_{j^{'} p}^{*} U_{k^{'} r} e^{i δ E_{r q} t / ℏ}), \end{matrix}

(73)

which can also be obtained from Eq. (67) by replacing S_p′q′(t) with δ_p′q′ and assuming τ = 0. Combining the above expression with Eq. (59), one can express the first order inhomogeneous term, Eq. (58), as follows:

\begin{matrix} I^{(1)} (t) & = & - \frac{i}{ℏ} \sum_{p, q = 1}^{N} | φ_{p} 〉 〈 φ_{q} | \sum_{r = 1}^{N} \sum_{j, k = 1}^{N} \sum_{j^{'}, k^{'} = 1}^{N} σ_{j^{'} k^{'}} (0) w_{j^{'} k^{'}} \\ \times (U_{j p}^{*} U_{k r} U_{j^{'} r}^{*} U_{k^{'} q} e^{i δ E_{p r} t / ℏ} - U_{j r}^{*} U_{k q} U_{j^{'} p}^{*} U_{k^{'} r} e^{i δ E_{r q} t / ℏ}) \\ \times \{{\tilde{J}}_{j k} (e^{- K_{j k, j^{'} k^{'}} (t)} f_{j k, k^{'}} (t) - 1)) \\ + (δ_{j k} (M_{j, j^{'} k^{'}} (t) + h_{j, k^{'}} (t))\} . \end{matrix}

(74)

Similar expressions for the second order inhomogeneous term, Eq. (62), can be obtained by replacing S_p′q′(t) in Eq. (67) with $U_{j^{''} p^{'}}^{*} U_{k^{''} q^{'}}$ ⁠. The resulting expression is as follows:

\begin{matrix} [T_{j k} (t), T_{j^{'} k^{'}} (τ) T_{j^{″} k^{″}} (0)] & = & \sum_{p, q = 1}^{N} | φ_{p} 〉 〈 φ_{q} | \sum_{p^{'}, q^{'} = 1}^{N} (δ_{q^{'} q} \sum_{r = 1}^{N} U_{j p}^{*} U_{k r} U_{j^{'} r} U_{k^{'} p^{'}}) e^{i δ E_{p p^{'}} t / ℏ} e^{i δ E_{p^{'} r} (t - τ) / ℏ} \\ - U_{j q^{'}}^{*} U_{k q} U_{j^{'} p}^{*} U_{k^{'} p^{'}} e^{i (δ E_{p p^{'}} - δ E_{q q^{'}}) t / ℏ} (e^{i δ E_{p^{'} p} (t - τ) / ℏ}) U_{j^{″} p^{'}}^{*} U_{k^{″} q^{'}} . \end{matrix}

(75)

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

\begin{matrix} I^{(2)} (t) & = & - \sum_{j, k = 1}^{N} \sum_{j^{'}, k^{'} = 1}^{N} \sum_{j^{″}, k^{″} = 1}^{N} σ_{j^{″} k^{″}} (0) \sum_{p, q = 1}^{N} | φ_{p} 〉 〈 φ_{q} | \sum_{p^{'}, q^{'} = 1}^{N} \sum_{r = 1}^{N} U_{j^{″} p^{'}}^{*} U_{k^{″} q^{'}} e^{i δ E_{p p^{'}} t / ℏ} (δ_{q^{'} q} U_{j p}^{*} U_{k r} U_{j^{'} r} U_{k^{'} p^{'}}) - (δ_{r p} U_{j q^{'}}^{*} U_{k q} U_{j^{'} p}^{*} U_{k^{'} p} e^{- i δ E_{q q^{'}} t / ℏ}) \\ \times (J_{j k} J_{j^{'} k^{'}} {\tilde{F}}_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t; δ E_{p^{'} r}) + δ_{j k} J_{j^{'} k^{'}} {\tilde{H}}_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t; δ E_{p^{'} r}) + J_{j k} δ_{j^{'} k^{'}} {\tilde{H}}_{j k, j^{'}}^{(2), j^{″} k^{″}} (t; δ E_{p^{'} r})) + (δ_{j k} δ_{j^{'} k^{'}} {\tilde{L}}_{j, j^{'}}^{j^{″} k^{″}} (t; δ E_{p^{'} r})) + H . c ., \end{matrix}

(76)

where the definitions of Eqs. (C28)–(C31) in Appendix C have been used.

III. TWO-LEVEL SYSTEM WITH INDEPENDENT BATH

A. 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₁₂ = J₂₁ and only one Debye–Waller factor, w = w₁₂ = w₂₁. Then, introducing $θ = θ_{1}^{†} θ_{2}$ for the present case, the renormalized zeroth order system Hamiltonian and the first order term of the Hamiltonian defined by Eq. (23) reduce to

{\tilde{H}}_{0, s} = {\tilde{E}}_{1} | 1 〉 〈 1 | + {\tilde{E}}_{2} | 2 〉 〈 2 | + J w (| 1 〉 〈 2 | + | 2 〉 〈 1 |),

(77)

\begin{align} {\tilde{H}}_{1} & = D_{1} | 1 〉 〈 1 | + D_{2} | 2 〉 〈 2 | \\ + J (θ - w) | 1 〉 〈 2 | + J (θ^{†} - w) | 2 〉 〈 1 |, \end{align}

(78)

where w = ⟨θ⟩ = ⟨θ^†⟩ and D₁ and D₂ have been defined by Eq. (25). Note also that, for the present case,

{\tilde{J}}_{12} = {\tilde{J}}_{21} = w J,

(79)

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

E_{1, 2} = \frac{{\tilde{E}}_{1} + {\tilde{E}}_{2}}{2} \pm \frac{{\tilde{E}}_{1} - {\tilde{E}}_{2}}{2} \sec (2 ξ),

(80)

with subscripts 1 and 2 on the left hand side denoting and denoting the + sign, respectively − signs on the right-hand side, and $ξ = \tan^{- 1} (2 J w / ({\tilde{E}}_{1} - {\tilde{E}}_{2})) / 2$ ⁠. The corresponding eigenstates are given by

| φ_{1} 〉 = \cos ξ | 1 〉 + \sin ξ | 2 〉,

(81)

| φ_{2} 〉 = - \sin ξ | 1 〉 + \cos ξ | 2 〉 .

(82)

Thus, for the present case, U₁₁ = U₂₂ = cos ξ and U₂₁ = −U₁₂ = sin ξ.

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

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

\begin{matrix} K (t) & \equiv & K_{12, 12} (t) = \frac{2}{π ℏ} \int_{0}^{\infty} d ω \frac{J (ω)}{ω^{2}} W_{h} {(ω)}^{2} \\ \times (\coth (\frac{β ℏ ω}{2}) \cos (ω t) - i \sin (ω t)) \\ = & K_{21, 21} (t) = - K_{12, 21} (t) = - K_{21, 12} (t), \end{matrix}

(83)

\begin{matrix} M (t) & \equiv & M_{1, 12} (t) = \frac{1}{π} \int_{0}^{\infty} d ω \frac{J (ω)}{ω} W_{h} (ω) (1 - W_{h} (ω)) \\ \times (\cos (ω t) - i \coth (\frac{β ℏ ω}{2}) \sin (ω t)) \\ = & M_{2, 21} (t) = - M_{1, 21} (t) = - M_{2, 12} (t), \end{matrix}

(84)

\begin{matrix} C (t) & \equiv & C_{11} (t) = \frac{ℏ}{π} \int_{0}^{\infty} d ω J (ω) {(1 - W_{h} (ω))}^{2} \\ \times (\coth (\frac{β ℏ ω}{2}) \cos (ω t) - i \sin (ω t)) \\ = & C_{22} (t) . \end{matrix}

(85)

For all other components, $K_{j k, j^{'} k^{'}} (ω) = 0$ ⁠, $M_{j, j^{'} k^{'}} (ω) = 0$ ⁠, and $C_{j k} (ω) = 0$ ⁠.

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

\begin{align} R (t) {\tilde{σ}}_{I} (t) & = \int_{0}^{t} d τ \{C (t - τ) ([T_{11} (t), T_{11} (τ) {\tilde{σ}}_{I} (t)] + [T_{22} (t), T_{22} (τ) {\tilde{σ}}_{I} (t)]) + J w M (t - τ) ([Δ T_{d} (t), Δ T_{c} (τ) {\tilde{σ}}_{I} (t)])) \\ ((- [Δ T_{c} (t), Δ T_{d} (τ) {\tilde{σ}}_{I} (t)]) + J^{2} w^{2} (e^{- K (t - τ)} - 1) ([T_{12} (t), T_{12} (τ) {\tilde{σ}}_{I} (t)] + [T_{21} (t), T_{21} (τ) {\tilde{σ}}_{I} (t)])) \\ (+ J^{2} w^{2} (e^{K (t - τ)} - 1) ([T_{12} (t), T_{21} (τ) {\tilde{σ}}_{I} (t)] + [T_{21} (t), T_{12} (τ) {\tilde{σ}}_{I} (t)])\} + H . c ., \end{align}

(86)

where

Δ T_{d} (t) = T_{22} (t) - T_{11} (t),

(87)

Δ T_{c} (t) = T_{21} (t) - T_{12} (t) .

(88)

Equation (86) 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 defined in the basis of eigenstates of ${\tilde{H}}_{0, s}$ ⁠, as detailed in Sec. II C. For this, $R_{p q}^{p^{'} q^{'}} (t)$ given by Eq. (72) needs to be calculated, which in turn requires calculation of $W_{j k, j^{'} k^{'}}^{p q} (t)$ ⁠, $Y_{j, j^{'} k^{'}}^{p q} (t)$ ⁠, and $X_{j j^{'}}^{p q} (t)$ defined by Eqs. (68)–(70). For the present model of two state systems coupled to independent baths, most of these are zero except for a few functions, as defined below.

First, all of nonzero $W_{j k, j^{'} k^{'}}^{p q} (t)$ ’s defined by Eq. (68) for the present case reduce to one of the following two functions:

\begin{matrix} W_{-}^{p q} (t) & \equiv & W_{12, 12}^{p q} (t) = W_{21, 21}^{p q} (t) \\ = & J^{2} w^{2} \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} (e^{- K (t - τ)} - 1), \end{matrix}

(89)

\begin{matrix} W_{+}^{p q} (t) & \equiv & W_{12, 21}^{p q} (t) = W_{21, 12}^{p q} (t) \\ = & J^{2} w^{2} \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} (e^{K (t - τ)} - 1) . \end{matrix}

(90)

All other terms of $W_{j k, j^{'} k^{'}}^{p q} (t)$ are zero by definition. Similarly, all of nonzero $Y_{j, j^{'} k^{'}}^{p q} (t)$ ’s and $X_{j j^{'}}^{p q} (t)$ ’s can be specified by

\begin{matrix} Y^{p q} (t) & \equiv & Y_{1, 12}^{p q} (t) = Y_{2, 21}^{p q} (t) \\ = & J w \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} M (t - τ) \\ = & - Y_{1, 21}^{p q} (t) = - Y_{2, 12}^{p q} (t), \end{matrix}

(91)

\begin{matrix} X^{p q} (t) & \equiv & X_{11}^{p q} (t) = X_{22}^{p q} (t) \\ = & \int_{0}^{t} d τ e^{i δ E_{p q} (t - τ) / ℏ} C (t - τ) . \end{matrix}

(92)

All other terms with indices different from the above are zero. Thus, with $W_{\pm}^{p q} (t)$ ⁠, $Y^{p q} (t)$ ⁠, and $X^{p q} (t)$ as defined above, all of $R_{p q}^{p^{'} q^{'}} (t)$ constituting Eq. (72) for the present case can be calculated.

The first order inhomogeneous term, Eq. (74), involves additional functions f_jk,k′(t) and h_j,k′(t), which can also be simplified for the present model. Let us define

\begin{matrix} f (t) & \equiv & \exp \{\frac{2 i}{π ℏ} \int_{0}^{\infty} d ω \frac{J (ω)}{ω^{2}} W_{h} {(ω)}^{2} \sin (ω t)\} \\ = & f_{12, 1} (t) . \end{matrix}

(93)

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′(t), only specification of the following function is necessary:

\begin{matrix} h (t) & \equiv & \frac{2}{π} \int_{0}^{\infty} d ω \frac{J (ω)}{ω} (1 - W_{h} (ω)) W_{h} (ω) \cos (ω t) \\ = & h_{1, 1} (t) = h_{2, 2} (t) . \end{matrix}

(94)

For other cases, h_1,2(t) = h_2,1(t) = 0.

The second order inhomogeneous term, Eq. (76), consists of much more terms that involve time integration of four different types of time correlation functions as described in Appendix C. Although the calculation of all these expressions is straightforward, actual numerical implementation is nontrivial and will be the subject of future work.

B. Model calculations without inhomogeneous terms for Ohmic spectral density

Numerical calculations neglecting the contribution of inhomogeneous terms are provided here. The objective of these calculations is to demonstrate the feasibility of using the p-PQME expressions derived in previous subsections and to investigate the effects of the weighting function W_h(ω). Thus, the results to be presented do not yet have the quantitative accuracy for all time due to the neglect of inhomogeneous terms. However, they still provide important qualitative information and reliable steady state limits for the cases where the inhomogeneous terms vanish.

First, it is useful to further simplify the expressions for $R_{p q}^{p^{'} q^{'}} (t)$ ’s. Summing up all possible 16 cases of j, k, j′, k′ = 1, 2 in Eq. (72), it is straightforward to show that it can be simplified as follows:

\begin{matrix} R_{p q}^{p^{'} q^{'}} (t) & = & \frac{1}{ℏ^{2}} \{δ_{q^{'} q} e^{i δ E_{p p^{'}} t / ℏ} \sum_{r = 1}^{2} (A_{p r r p^{'}} X^{p^{'} r} (t) + (B_{p r}^{(1)} B_{r p^{'}}^{(2)} - B_{p r}^{(2)} B_{r p^{'}}^{(1)}) Y^{p^{'} r} (t) + C_{p r r p^{'}}^{(1)} W_{-}^{p^{'} r} (t) + C_{p r r p^{'}}^{(2)} W_{+}^{p^{'} r} (t)) - e^{i (δ E_{p p^{'}} - δ E_{q q^{'}}) t / ℏ}) \\ \times ((A_{q^{'} q p p^{'}} X^{p^{'} p} (t) + (B_{q^{'} q}^{(1)} B_{p p^{'}}^{(2)} - B_{q^{'} q}^{(2)} B_{p p^{'}}^{(1)}) Y^{p^{'} p} (t) + C_{q^{'} q p p^{'}}^{(1)} W_{-}^{p^{'} p} (t) + C_{q^{'} q p p^{'}}^{(2)} W_{+}^{p^{'} p} (t))\} + [C . C ., p \leftrightarrow q, p^{'} \leftrightarrow q^{'}], \end{matrix}

(95)

where

A_{pqrs} = U_{1 p} U_{1 q} U_{1 r} U_{1 s} + U_{2 p} U_{2 q} U_{2 r} U_{2 s},

(96)

B_{p q}^{(1)} = U_{1 p} U_{1 q} - U_{2 p} U_{2 q},

(97)

B_{p q}^{(2)} = U_{1 p} U_{2 q} - U_{2 p} U_{1 q},

(98)

C_{pqrs}^{(1)} = U_{1 p} U_{2 q} U_{1 r} U_{2 s} + U_{2 p} U_{1 q} U_{2 r} U_{1 s},

(99)

C_{pqrs}^{(2)} = U_{1 p} U_{2 q} U_{2 r} U_{1 s} + U_{2 p} U_{1 q} U_{1 r} U_{2 s} .

(100)

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

J (ω) = π ℏ η ω e^{- ω / ω_{c}} .

(101)

It is assumed that the two parameters of the above spectral density are set to η = 1 and ℏω_c = 200 cm⁻¹. In addition, the temperature and electronic coupling are set to T = 300 K and J = 300 cm⁻¹. For this choice of parameters, the bath dynamics, system electronic coupling, and thermal energy are all comparable. Thus, the dynamics are expected to be the border-line case between incoherent and coherent quantum dynamics. Two cases of relative site energies, E₁ = E₂ and E₁ − E₂ = 200 cm⁻¹ will be considered.

There are infinite number of possible choices available for W_h(ω) that satisfies the requirement of Eq. (7). One simple but flexible choice is the following function known as (cumulative) Weibull distribution⁴⁵ or the complement of a stretched exponential function,

W_{h} (ω) = 1 - e^{- {(ω / ω_{h})}^{α}} .

(102)

For small ω, $W_{h} (ω) \approx {(ω / ω_{h})}^{α}$ ⁠. For large α, W_h(ω) approaches the step function at ω = ω_h.

While choosing a value α in Eq. (102) such that α ≥ 1 is sufficient for ensuring that the Debye–Waller factor does not vanish, it is not clear whether it also guarantees a well-behaving time evolution equation. In order to check this, it is important to examine how the three time correlation functions defined by Eqs. (83)–(85) behave for different choices of α. Figure 1 shows the real and imaginary parts of $K (t)$ defined by Eq. (83) for two different values of α = 1 and 2 with the choice of ω_h = ω_c. In both cases, the real and imaginary parts decay to zero quickly enough to make both $W_{-}^{p q} (t)$ and $W_{+}^{p q} (t)$ ⁠, defined, respectively, by Eqs. (89) and (90), converge to finite values.

FIG. 1.

View large Download slide

Real (upper panel) and imaginary (Lower panel) parts of $K (t)$ for ω_h/ω_c = 1 and for two cases of α = 1 and 2.

On the other hand, for the case of $M (t)$ ⁠, it turns out that the choice of α = 1 does not result in a stable time evolution equation, in general. Figure 2 shows the real and imaginary parts of $M (t)$ ⁠, also for α = 1 and 2 with the choice of ω_h = ω_c. For α = 1, while the real part decays to zero quickly, the imaginary part is seen to decay very slowly. In fact, due to the slow decay, the $Y^{p q} (t)$ defined by Eq. (91) diverges, in general, in this case. Test calculations for other values of α show that such divergence persists up to α = 3/2. On the other hand, for the case of α = 2 shown in Fig. 2, it is clear that the imaginary part decays to zero quickly, ensuring for $Y^{p q} (t)$ to converge to a finite value in the steady state limit.

FIG. 2.

View large Download slide

Real (upper panel) and imaginary (lower panel) parts of $M (t) / (ℏ ω_{c})$ for ω_h/ω_c = 1 and for two cases of α = 1 and 2.

For the case of $C (t)$ ⁠, as shown in Fig. 3, both real and imaginary parts decay to zero quickly already for α = 1, resulting in well behaving $X^{p q} (t)$ ⁠. Summing up the results shown in Figs. 1–3, while the choice of α = 1 is acceptable as far as the Deby–Waller factor, $K (t)$ ⁠, and $C (t)$ are concerned, it is not appropriate due to resulting slow decay of $M (t)$ ⁠. On the other hand, the choice of α = 2 results in all well-behaving and convergent functions that constitute the relaxation operator. This is also true for all the terms involved in the inhomogeneous terms as well.

FIG. 3.

View large Download slide

Real (upper panel) and imaginary (lower panel) parts of $C (t) / (ℏ^{2} ω_{c}^{2})$ for ω_h/ω_c = 1 and for two cases of α = 1 and 2.

Figure 4 shows the results for E₁ = E₂ for different values of ω_h with the choice of α = 2, where W_h(ω) becomes a complement of a Gaussian function. The result for the smallest value of ω_h among those shown (ω_h/ω_c = 0.5) is close to the limit of full PQME, without coherence, whereas that for the largest value of ω_h among those shown (ω_h/ω_c = 4) is close to the full second order time-local QME, which has maximum coherence. For ω_h/ω_c = 1, the time dependent population exhibits an intermediate character between the two limits.

FIG. 4.

View large Download slide

Populations of the donor (site 1) for different values of ω_h for E₁ = E₂ and α = 2 in W_h(ω), Eq. (102). Other parameters of the model are as follows: η = 1, ω_c = 200 cm⁻¹, T = 300 K, and J = 300 cm⁻¹.

Figure 5 shows the results for an asymmetric case, where E₁ − E₂ = 200 cm⁻¹. All other parameters, including α = 2, remain the same as those for Fig. 4. In this case, the steady state limits of the population as well as the coherence pattern vary with ω_h. The variation of the steady state limit with ω_h reflects a different extent of system–bath entanglement depending on the extent of polaron transformation. The smaller the value of ω_h, the closer the steady states are to the original localized states 1 and 2, for which the energy gap becomes the maximum.

FIG. 5.

View large Download slide

Populations of the donor (site 1) for E₁ − E₂ = 200 cm⁻¹ and α = 2 in W_h(ω), Eq. (102). All other parameters are the same as in Fig. 4.

IV. CONCLUDING REMARKS

This work has provided a general framework to overcome a known issue of the original second 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. (72), (74), and (76), can be applied for any kind 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. It is important to note that the expressions provided here are applicable even to the case where the same bath mode is partly transformed, with the remaining part untransformed, and is, thus, more general than the case where the bath is divided into two disjoint transformed and untransformed groups.

Numerical tests for a simple two level system coupled to an Ohmic bath demonstrate that appropriate specification of the weighting function W_h(ω) can tune the extent of coherence and the extent of system–bath entanglement in the steady state limit. This adds a new dimension of flexibility in incorporating the PT approach into a QME calculation. The flexibility in choosing the weighting function W_h(ω) in all the expressions derived for the p-PQME presented here leaves open various possibilities of adapting or improving the methodology. For example, the 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, the subject of forthcoming work, 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.

ACKNOWLEDGMENTS

This work was mainly supported by the National Science Foundation (Grant No. CHE-1900170). The author also acknowledges partial support from the U.S. Department of Energy, Office of Sciences, Office of Basic Energy Sciences (Grant No. DE-SC0021413), and support from the Korea Institute for Advanced Study (KIAS) through its KIAS Scholar program.

AUTHOR DECLARATIONS

Conflict of Interest

The author has no conflicts to disclose.

Author Contributions

Seogjoo J. Jang: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).

DATA AVAILABILITY

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

APPENDIX A: DERIVATION OF Eq. (41)

The trace of ${\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ)$ with ρ_b in Eq. (40), namely, its thermal average is expressed as follows:

\begin{align} 〈 {\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) 〉 & = 〈 θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 - w_{j k} 〈 θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 - w_{j^{'} k^{'}} 〈 θ_{j}^{†} (t) θ_{k} (t) 〉 + w_{j k} w_{j^{'} k^{'}} \\ + δ_{j k} (〈 D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 - 〈 D_{j} (t) 〉 w_{j^{'} k^{'}}) + δ_{j^{'} k^{'}} (〈 θ_{j}^{†} (t) θ_{k} (t) D_{j^{'}} (τ) 〉 - w_{j k} 〈 D_{j^{'}} (τ) 〉) + δ_{j k} δ_{j^{'} k^{'}} 〈 D_{j} (t) D_{j^{'}} (τ) 〉 \\ = 〈 θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 - w_{j k} w_{j^{'} k^{'}} + δ_{j k} 〈 D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 + δ_{j^{'} k^{'}} 〈 θ_{j}^{†} (t) θ_{k} (t) D_{j^{'}} (τ) 〉 + δ_{j k} δ_{j^{'} k^{'}} 〈 D_{j} (t) D_{j^{'}} (τ) 〉 . \end{align}

(A1)

In obtaining the second equality of the above equation, the identities $〈 θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 = w_{j^{'} k^{'}}$ ⁠, $〈 θ_{j}^{†} (t) θ_{k} (t) 〉 = w_{j k}$ ⁠, and ⟨D_j(t)⟩ = ⟨D_j′(τ)⟩ = 0 have been used.

The first term in the second equality of Eq. (A1) has the same form as that for the full PT¹⁵ except for the additional factor W_h(ω_n) multiplied with δg_n,jk and δg_n,j′k′, respectively. Thus, following the same procedure as in previous work,¹⁵ it can be shown to be

〈 θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 = w_{j k} w_{j^{'} k^{'}} e^{- K_{j k, j^{'} k^{'}} (t - τ)},

(A2)

where $e^{K_{j k, j^{'} k^{'}} (t - τ)}$ in the exponent is defined by Eq. (42).

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

\begin{align} 〈 D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 & = \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) ⟨(b_{n} e^{- i ω_{n} t} + b_{n}^{†} e^{i ω_{n} t})) \\ \times (\prod_{n^{'}} e^{δ g_{n^{'}, j^{'} k^{'}} W_{h} (ω_{n^{'}}) (b_{n^{'}}^{†} e^{i ω_{n^{'}} τ} - b_{n^{'}} e^{- i ω_{n^{'}} τ})}⟩ . \end{align}

(A3)

In the above expression, the product of bath modes (indexed by n′) can be averaged independently for n′ ≠ n, resulting in

\begin{align} ⟨\prod_{n^{'} \neq n} e^{δ g_{n^{'}, j^{'} k^{'}} W_{h} (ω_{n^{'}}) (b_{n^{'}}^{†} e^{i ω_{n^{'}} τ} - b_{n^{'}} e^{- i ω_{n^{'}} τ})}⟩ \\ = \prod_{n^{'} \neq n} e^{- \coth (β ℏ ω_{n^{'}} / 2) δ g_{n^{'}, j^{'} k^{'}}^{2} W_{h} {(ω_{n^{'}})}^{2} / 2} . \end{align}

(A4)

On the other hand, the term for n′ = n in Eq. (A3) needs to be calculated together with the linear term as follows:

\begin{align} 〈 (b_{n} e^{- i ω_{n} t} + b_{n}^{†} e^{i ω_{n} t}) e^{δ g_{n, j^{'} k^{'}} W_{h} (ω_{n}) (b_{n}^{†} e^{i ω_{n} τ} - b_{n} e^{- i ω_{n} τ})} 〉 \\ = e^{- i ω_{n} t} 〈 b_{n} e^{γ^{*} b_{n}^{†} - γ b_{n}} 〉 + e^{i ω_{n} t} 〈 b_{n}^{†} e^{γ^{*} b_{n}^{†} - γ b_{n}} 〉, \end{align}

(A5)

where γ is defined by

γ = δ g_{n, j^{'} k^{'}} W_{h} (ω_{n}) e^{- i ω_{n} τ},

(A6)

and γ* is the complex conjugate of γ. Using the fact that $e^{γ^{*} b_{n}^{†} - γ b_{n}} = e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} e^{| γ |^{2} / 2}$ and the following identity,

\begin{align} 〈 b_{n} e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 & = - \frac{\partial}{\partial γ} 〈 e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 \\ = \frac{γ^{*}}{1 - e^{- β ℏ ω_{n}}} e^{- γ γ^{*} / (1 - e^{- β ℏ ω_{n}})}, \end{align}

(A7)

it is straightforward to show that

〈 b_{n} e^{γ^{*} b_{n}^{†} - γ b_{n}} 〉 = \frac{γ^{*}}{1 - e^{- β ℏ ω_{n}}} e^{- \coth (β ℏ ω_{n} / 2) | γ |^{2} / 2} .

(A8)

Similarly, the bath average in the second term of Eq. (A5) can be expressed as

〈 b_{n}^{†} e^{γ^{*} b_{n}^{†} - γ b_{n}} 〉 = 〈 b_{n}^{†} e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 e^{| γ |^{2} / 2} .

(A9)

Using the following identity,

b_{n}^{†} e^{- γ b_{n}} = e^{- γ b_{n}} b_{n}^{†} + γ e^{- γ b_{n}},

(A10)

the average over the bath in Eq. (A9) can be expressed as

〈 b_{n}^{†} e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 = 〈 e^{- γ b_{n}} b_{n}^{†} e^{γ^{*} b_{n}^{†}} 〉 + γ 〈 e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 .

(A11)

Now, employing the following identity:

\begin{align} 〈 e^{- γ b_{n}} b_{n}^{†} e^{γ^{*} b_{n}^{†}} 〉 & = \frac{\partial}{\partial γ^{*}} 〈 e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 \\ = - \frac{γ}{1 - e^{- β ℏ ω_{n}}} e^{- γ γ^{*} / (1 - e^{- β ℏ ω_{n}})}, \end{align}

(A12)

in Eq. (A11) and combining the resulting expression with Eq. (A9), we obtain the following expression:

〈 b_{n}^{†} e^{γ^{*} b_{n}^{†} - γ b_{n}} 〉 = - \frac{γ e^{- β ℏ ω_{n}}}{1 - e^{- β ℏ ω_{n}}} e^{- \coth (β ℏ ω_{n} / 2) | γ |^{2} / 2} .

(A13)

Combining Eqs. (A8) and (A13) with the definition of Eq. (A6) leads to

\begin{align} 〈 (b_{n} e^{- i ω_{n} t} + b_{n}^{†} e^{i ω_{n} t}) e^{δ g_{n, j^{'} k^{'}} W_{h} (ω_{n}) (b_{n}^{†} e^{i ω_{n} τ} - b_{n} e^{- i ω_{n} τ})} 〉 \\ = δ g_{n, j^{'} k^{'}} W_{h} (ω_{n}) e^{- \coth (β ℏ ω_{n} / 2) δ g_{n, j^{'} k^{'}}^{2} W_{h} {(ω_{n})}^{2} / 2} \\ \times (\cos (ω_{n} (t - τ)) - i \coth (\frac{β ℏ ω_{n}}{2}) \sin (ω_{n} (t - τ))) . \end{align}

(A14)

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

\begin{align} 〈 D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) 〉 \\ = w_{j^{'} k^{'}} \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) δ g_{n, j^{'} k^{'}} W_{h} (ω_{n}) \\ \times (\cos (ω_{n} (t - τ)) - i \coth (\frac{β ℏ ω_{n}}{2}) \sin (ω_{n} (t - τ))) . \end{align}

(A15)

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

\begin{align} 〈 θ_{j}^{†} (t) θ_{k} (t) D_{j^{'}} (τ) 〉 & = \sum_{n} ℏ ω_{n} g_{n, j^{'}} (1 - W_{h} (ω_{n})) \\ \times ⟨\prod_{n^{'}} e^{δ g_{n^{'}, j k} W_{h} (ω_{n^{'}}) (b_{n^{'}}^{†} e^{i ω_{n^{'}} t} - b_{n^{'}} e^{- i ω_{n^{'}} t})}) \\ \times ((b_{n} e^{- i ω_{n} τ} + b_{n}^{†} e^{i ω_{n} τ})⟩ . \end{align}

(A16)

The bath average term in the above expression can be calculated in a manner similar to the third term of Eq. (A1), which has been described above, but with a different definition of $γ = δ g_{n, j k} W_{h} (ω_{n}) e^{- i ω_{n} t}$ ⁠. For this, $〈 e^{γ^{*} b_{n}^{†} - γ b_{n}} b_{n} 〉$ and $〈 e^{γ^{*} b_{n}^{†} - γ b_{n}} b_{n}^{†} 〉$ need to be calculated. For the first of these terms, the following identity can be used:

e^{γ^{*} b_{n}^{†}} b_{n} = (b_{n} - γ^{*}) e^{γ^{*} b_{n}^{†}} .

(A17)

Thus,

\begin{align} 〈 e^{γ^{*} b_{n}^{†} - γ b_{n}} b_{n} 〉 & = 〈 e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} b_{n} 〉 e^{| γ |^{2} / 2} \\ = 〈 e^{- γ b_{n}} (b_{n} e^{γ^{*} b_{n}^{†}} - γ^{*} e^{γ b_{n}^{†}}) 〉 e^{| γ |^{2} / 2} \\ = (〈 b_{n} e^{- γ b_{n}} e^{γ^{*} b_{n}^{†}} 〉 - γ^{*} 〈 e^{- γ b_{n}} e^{γ b_{n}^{†}} 〉) e^{| γ |^{2} / 2} . \end{align}

(A18)

Employing Eq. (A7), one can simplify the above expression as follows:

〈 e^{γ^{*} b_{n}^{†} - γ b_{n}} b_{n} 〉 = γ^{*} \frac{e^{- β ℏ ω_{n}}}{1 - e^{- β ℏ ω_{n}}} e^{- | γ |^{2} \coth (β ℏ ω_{n} / 2)} .

(A19)

Combining the above identity with Eq. (A12) followed by further calculation leads to the following expression:

\begin{align} 〈 e^{δ g_{n, j k} W_{h} (ω_{n}) (b_{n}^{†} e^{i ω_{n} t} - b_{n} e^{- i ω_{n} t})} (b_{n} e^{- i ω_{n} τ} + b_{n}^{†} e^{i ω_{n} τ}) 〉 \\ = δ g_{n, k j} W_{h} (ω_{n}) e^{- \coth (β ℏ ω_{n} / 2) δ g_{n, j k}^{2} W_{h} {(ω_{n})}^{2} / 2} \\ \times (\cos (ω_{n} (t - τ)) - i \coth (\frac{β ℏ ω_{n}}{2}) \sin (ω_{n} (t - τ))), \end{align}

(A20)

where note the use of the fact that δg_n,kj = −δg_n,jk on the right-hand side of the above equation. As a result, Eq. (A16) can be expressed as

\begin{align} 〈 θ_{j}^{†} (t) θ_{k} (t) D_{j^{'}} (τ) 〉 \\ = w_{j k} \sum_{n} ℏ ω_{n} g_{n, j^{'}} (1 - W_{h} (ω_{n})) δ g_{n, k j} W_{h} (ω_{n}) \\ \times (\cos (ω_{n} (t - τ)) - i \coth (\frac{β ℏ ω_{n}}{2}) \sin (ω_{n} (t - τ))) . \end{align}

(A21)

The calculation of the last term in the second equality of Eq. (A1) is straightforward. The resulting expression is as follows:

\begin{align} 〈 D_{j} (t) D_{j^{'}} (τ) 〉 & = \sum_{n} ℏ^{2} ω_{n}^{2} g_{n, j} g_{n, j^{'}} {(1 - W_{h} (ω_{n}))}^{2} \\ \times (〈 b_{n} b_{n}^{†} 〉 e^{- i ω_{n} (t - τ)} + 〈 b_{n}^{†} b_{n} 〉 e^{i ω_{n} (t - τ)}) \\ = \sum_{n} ℏ^{2} ω_{n}^{2} g_{n, j} g_{n, j^{'}} {(1 - W_{h} (ω_{n}))}^{2} \\ \times (\coth (\frac{β ℏ ω_{n}}{2}) \cos (ω_{n} (t - τ)) - i \sin (ω_{n} (t - τ))) . \end{align}

(A22)

In the second equality of the above equation, the fact that $〈 b_{n}^{†} b_{n} 〉 = 〈 b_{n} b_{n}^{†} 〉 - 1 = e^{- β ℏ ω_{n}} / (1 - e^{- β ℏ ω_{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. (58), can be expressed as follows:

\begin{align} {Tr}_{b} {{\tilde{B}}_{j k} (t) δ {\tilde{ρ}}_{b, j^{'} k^{'}}} & = J_{j k} {Tr}_{b} \{θ_{k^{'}} θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} ρ_{b}\} \\ - J_{j k} w_{j k} w_{j^{'} k^{'}} + δ_{j k} {Tr}_{b} \{θ_{k^{'}} D_{j} (t) θ_{j^{'}}^{†} ρ_{b}\} . \end{align}

(B1)

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

\begin{align} θ_{k^{'}} θ_{j}^{†} (t) θ_{k} (t) & = θ_{k^{'}} θ_{j}^{†} (t) θ_{k} (t) θ_{k^{'}}^{†} θ_{k^{'}} \\ = e^{\sum_{n} δ g_{n, j k} W_{h} (ω_{n}) (θ_{k^{'}} b_{n}^{†} θ_{k^{'}}^{†} e^{i ω_{n} t} - θ_{k^{'}} b_{n} θ_{k^{'}}^{†} e^{- i ω_{n} t})} θ_{k^{'}} . \end{align}

(B2)

In the exponent on the right-hand side of the above expression,

\begin{align} θ_{k^{'}} b_{n}^{†} θ_{k^{'}}^{†} & = b_{n}^{†} - W_{h} (ω_{n}) g_{n, k^{'}} [b_{n}^{†} - b_{n}, b_{n}^{†}] \\ = b_{n}^{†} + W_{h} (ω_{n}) g_{n, k^{'}}, \end{align}

(B3)

\begin{align} θ_{k^{'}} b_{n} θ_{k^{'}}^{†} & = b_{n} - W_{h} (ω_{n}) g_{n, k^{'}} [b_{n}^{†} - b_{n}, b_{n}] \\ = b_{n} + W_{h} (ω_{n}) g_{n, k^{'}} . \end{align}

(B4)

Therefore,

\begin{align} \sum_{n} δ g_{n, j k} W_{h} (ω_{n}) (θ_{k^{'}} b_{n}^{†} θ_{k^{'}}^{†} e^{i ω_{n} t} - θ_{k^{'}} b_{n} θ_{k^{'}}^{†} e^{- i ω_{n} t}) \\ = \sum_{n} δ g_{n, j k} W_{h} (ω_{n}) (b_{n}^{†} e^{i ω_{n} t} - b_{n} e^{- i ω_{n} t}) \\ + 2 i \sum_{n} δ g_{n, j k} g_{n, k^{'}} W_{h} {(ω_{n})}^{2} \sin (ω_{n} t) . \end{align}

(B5)

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

\begin{align} θ_{k^{'}} θ_{j}^{†} (t) θ_{k} (t) & = θ_{k^{'}} θ_{j}^{†} (t) θ_{k} (t) θ_{k^{'}}^{†} θ_{k^{'}} \\ = e^{\sum_{n} δ g_{n, j k} W_{h} (ω_{n}) (θ_{k^{'}} b_{n}^{†} θ_{k^{'}}^{†} e^{i ω_{n} t} - θ_{k^{'}} b_{n} θ_{k^{'}}^{†} e^{- i ω_{n} t})} θ_{k^{'}} \\ = f_{j k, k^{'}} (t) θ_{j}^{†} (t) θ_{k} (t) θ_{k^{'}} . \end{align}

(B6)

The above identity implies that the first term of Eq. (B1) (without J_jk) can be simplified to

\begin{align} {Tr}_{b} \{θ_{k^{'}} θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} ρ_{b}\} & = f_{j k, k^{'}} (t) 〈 θ_{j}^{†} (t) θ_{k} (t) θ_{k^{'}} θ_{j^{'}}^{†} 〉 \\ = f_{j k, k^{'}} (t) w_{j k} w_{j^{'} k^{'}} e^{- K_{j k, j^{'} k^{'}} (t)}, \end{align}

(B7)

where Eq. (A2) has been used.

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

\begin{align} θ_{k^{'}} D_{j} (t) & = θ_{k^{'}} D_{j} (t) θ_{k^{'}}^{†} θ_{k^{'}} = \sum_{n} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) \\ \times (θ_{k^{'}} b_{n} θ_{k^{'}}^{†} e^{- i ω_{n} t} + θ_{k^{'}} b_{n}^{†} θ_{k^{'}}^{†} e^{i ω_{n} t}) θ_{k^{'}} \\ = D_{j} (t) θ_{k^{'}} + 2 θ_{k^{'}} \sum_{n} ℏ ω_{n} g_{n, j} g_{n, k^{'}} \\ \times (1 - W_{h} (ω_{n})) W_{h} (ω_{n}) \cos (ω_{n} t), \end{align}

(B8)

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

{Tr}_{b} \{θ_{k^{'}} D_{j} (t) θ_{j^{'}}^{†} ρ_{b}\} = (M_{j, j^{'} k^{'}} (t) + h_{j, k^{'}} (t)) w_{j^{'} k^{'}},

(B9)

where the identity of Eq. (A15) for τ = 0, the definition of Eq. (43), and the definition of Eq. (61) have been used. Inserting Eqs. (B7) and (B9) into Eq. (B1), one can then obtain the following expression:

\begin{align} {Tr}_{b} {{\tilde{B}}_{j k} (t) δ {\tilde{ρ}}_{b, j^{'} k^{'}}} & = J_{j k} w_{j k} (f_{j k, k^{'}} (t) e^{- K_{j k, j^{'} k^{'}} (t)} - 1) w_{j^{'} k^{'}} \\ + δ_{j k} (M_{j, j^{'} k^{'}} (t) + h_{j, k^{'}} (t)) w_{j^{'} k^{'}} . \end{align}

(B10)

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

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. (62), can be expressed as follows:

\begin{align} {Tr}_{b} \{{\tilde{B}}_{j k} (t) {\tilde{B}}_{j^{'} k^{'}} (τ) δ {\tilde{ρ}}_{b, j^{″} k^{″}}\} & = {Tr}_{b} \{(J_{j k} (θ_{j}^{†} (t) θ_{k} (t) - w_{j k}) + δ_{j k} D_{j} (t)) (J_{j^{'} k^{'}} (θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) - w_{j^{'} k^{'}}) + δ_{j^{'} k^{'}} D_{j^{'}} (τ)) (θ_{j^{″}}^{†} ρ_{b} θ_{k^{″}} - w_{j^{″} k^{″}} ρ_{b})\} \\ = J_{j k} J_{j^{'} k^{'}} F_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t, τ) + δ_{j k} J_{j^{'} k^{'}} H_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t, τ) + J_{j k} δ_{j^{'} k^{'}} H_{j k, j^{'}}^{(2), j^{″} k^{″}} (t, τ) + δ_{j k} δ_{j^{'} k^{'}} L_{j, j^{'}}^{j^{″} k^{″}} (t, τ), \end{align}

(C1)

where

F_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t, τ) = {Tr}_{b} \{(θ_{j}^{†} (t) θ_{k} (t) - w_{j k}) (θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) - w_{j^{'} k^{'}}) (θ_{j^{″}}^{†} ρ_{b} θ_{k^{″}} - w_{j^{″} k^{″}} ρ_{b})\},

(C2)

H_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t, τ) = {Tr}_{b} \{D_{j} (t) (θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) - w_{j^{'} k^{'}}) (θ_{j^{″}}^{†} ρ_{b} θ_{k^{″}} - w_{j^{″} k^{″}} ρ_{b})\},

(C3)

H_{j k, j^{'}}^{(2), j^{″} k^{″}} (t, τ) = {Tr}_{b} \{(θ_{j}^{†} (t) θ_{k} (t) - w_{j k}) D_{j^{'}} (τ) (θ_{j^{″}}^{†} ρ_{b} θ_{k^{″}} - w_{j^{″} k^{″}} ρ_{b})\},

(C4)

L_{j, j^{'}}^{j^{″} k^{″}} (t, τ) = {Tr}_{b} \{D_{j} (t) D_{j^{'}} (τ) (θ_{j^{″}}^{†} ρ_{b} θ_{k^{″}} - w_{j^{″} k^{″}} ρ_{b})\} .

(C5)

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

The first bath term in Eq. (C1), $F_{j k, j^{'} k^{'}}^{j^{''} k^{''}} (t, τ)$ defined by Eq. (C2), can be expanded further and expressed as follows:

\begin{align} F_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t, τ) & = 〈 θ_{k^{″}} θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{j^{″}} 〉 \\ - w_{j k} 〈 θ_{k^{″}} θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{j^{″}}^{†} 〉 \\ - w_{j^{'} k^{'}} 〈 θ_{k^{″}} θ_{j}^{†} (t) θ_{k} (t) θ_{j^{″}}^{†} 〉 \\ + w_{j k} w_{j^{'} k^{'}} 〈 θ_{k^{″}} θ_{j^{″}}^{†} 〉 - w_{j^{″} k^{″}} ⟨(θ_{j}^{†} (t) θ_{k} (t) - w_{j k})) \\ (\times (θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) - w_{j^{'} k^{'}})⟩ . \end{align}

(C6)

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

\begin{align} 〈 θ_{k^{″}} θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{j^{″}} 〉 \\ = f_{j k, k^{″}} (t) f_{j^{'} k^{'}, k^{″}} (τ) 〈 θ_{j}^{†} (t) θ_{k} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{k^{″}} θ_{j^{″}} 〉 \\ = f_{j k, k^{″}} (t) f_{j k^{'}, k^{″}} (τ) w_{j k} w_{j^{'} k^{'}} w_{j^{″} k^{″}} \\ \times e^{- K_{j k, j^{″} k^{″}} (t) - K_{j^{'} k^{'}, j^{″} k^{″}} (τ) - K_{j k, j^{'} k^{'}} (t - τ)} . \end{align}

(C7)

For the second and third terms of Eq. (C6), Eq. (B7) can be used. The last term of Eq. (C6) corresponds to the first term that appears in the evaluation of $R (t)$ ⁠. Combining all of these, one can show that

\begin{align} F_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t, τ) & = w_{j k} w_{j^{'} k^{'}} w_{j^{″} k^{″}} \{e^{- K_{j k, j^{'} k^{'}} (t - τ)}) \\ \times (f_{j k, k^{″}} (t) f_{j^{'} k^{'}, k^{″}} (τ) e^{- K_{j k, j^{″} k^{″}} (t) - K_{j^{'} k^{'}, j^{″} k^{″}} (τ)} - 1) \\ (- f_{j k, k^{″}} (t) e^{- K_{j k, j^{″} k^{″}} (t)} - f_{j^{'} k^{'}, k^{″}} (τ) e^{- K_{j^{'} k^{'}, j^{″} k^{″}} (τ)} + 2\} . \end{align}

(C8)

The bath term contributing to the second term of Eq. (C1), i.e., $H_{j k, j^{'} k^{'}}^{(1), j^{''} k^{''}} (t, τ)$ defined by Eq. (C3), is expressed as follows:

\begin{align} H_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t, τ) & = 〈 θ_{k^{″}} D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{i^{″}}^{†} 〉 \\ - w_{j^{'} k^{'}} 〈 θ_{k^{″}} D_{j} (t) θ_{j^{″}}^{†} 〉 \\ - w_{j^{″} k^{″}} 〈 D_{j} (t) (θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) - w_{j^{'} k^{'}}) 〉 . \end{align}

(C9)

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

\begin{align} 〈 θ_{k^{″}} D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{j^{″}}^{†} 〉 \\ = f_{j^{'} k^{'}, k^{″}} (τ) 〈 D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{k^{″}} θ_{j^{″}}^{†} 〉 \\ + f_{j^{'} k^{'}, k^{″}} (τ) h_{j, k^{″}} (t) 〈 θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{k^{″}} θ_{j^{″}}^{†} 〉, \end{align}

(C10)

where the following identities have been used,

θ_{k^{″}} D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) = θ_{k^{″}} D_{j} (t) θ_{k^{″}}^{†} θ_{k^{″}} θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{k^{″}}^{†} θ_{k^{″}},

(C11)

θ_{k^{″}} D_{j} (t) θ_{k^{″}}^{†} = D_{j} (t) + h_{j, k^{″}} (t),

(C12)

θ_{k^{″}} θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{k^{″}}^{†} = f_{j^{'} k^{'}, k^{″}} (τ) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ),

(C13)

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

\begin{matrix} 〈 D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{k^{″}} θ_{j^{″}}^{†} 〉 & = & (M_{j, j^{'} k^{'}} (t - τ) + M_{j, j^{″} k^{″}} (t)) \\ \times w_{j^{'} k^{'}} w_{j^{″} k^{″}} e^{- K_{j^{'} k^{'}, j^{″} k^{″}} (τ)} . \end{matrix}

(C14)

Combining this with the identity given by Eq. (A2) (with the replacement of t → τ and τ → 0), one can show that Eq. (C10) can be expressed as follows:

\begin{matrix} 〈 θ_{k^{″}} D_{j} (t) θ_{j^{'}}^{†} (τ) θ_{k^{'}} (τ) θ_{j^{″}}^{†} 〉 & = & f_{j^{'} k^{'}, k^{″}} (τ) (M_{j, j^{'} k^{'}} (t - τ)) \\ + (M_{j, j^{″} k^{″}} (t) + h_{j, k^{″}} (t)) w_{j^{'} k^{'}} \\ \times w_{j^{″} k^{″}} e^{- K_{j^{'} k^{'}, j^{″} k^{″}} (τ)} . \end{matrix}

(C15)

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

〈 θ_{k^{″}} D_{j} (t) θ_{j^{″}}^{†} 〉 = (M_{j, j^{″} k^{″}} (t) + h_{j, k^{″}} (t)) w_{j^{″} k^{″}},

(C16)

where Eqs. (A21) and (C12) along with the definition of Eq. (43) have been used.

For the last term on the right-hand side of Eqs. (C9), Eq. (A15) can be employed. As a result, Eq. (C9) can be expressed as

\begin{matrix} H_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t, τ) & = & (f_{j^{'} k^{'}, k^{″}} (τ) e^{- K_{j^{'} k^{'}, j^{″} k^{″}} (τ)} - 1) \\ \times w_{j^{'} k^{'}} w_{j^{″} k^{″}} (M_{j, j^{'} k^{'}} (t - τ) + M_{j, j^{″} k^{″}} (t) + h_{j, k^{″}} (t)) . \end{matrix}

(C17)

$H_{j k, j^{'}}^{(2), j^{''} k^{''}} (t, τ)$ defined by Eq. (C4) can be calculated in a similar manner but can, in fact, be calculated using its relation to $H_{j^{'}, k j}^{(1), k^{''} j^{''}} (τ, t)$ as follows:

\begin{matrix} H_{j k, j^{'}}^{(2), j^{″} k^{″}} (t, τ) & = & H_{j^{'}, k j}^{(1), k^{″} j^{″}} {(τ, t)}^{*} \\ = & {(f_{k j, j^{″}} (t) e^{- K_{k j, k^{″} j^{″}} (t)} - 1)}^{*} w_{j k} w_{j^{″} k^{″}} \\ \times {(M_{j^{'}, k j} (τ - t) + M_{j^{'}, k^{″} j^{″}} (τ) + h_{j^{'}, j^{″}} (τ))}^{*} \\ = & (f_{j k, k^{″}} (t) e^{- K_{j k, j^{″} k^{″}} (t)} - 1) w_{j k} w_{j^{″} k^{″}} \\ \times (M_{j^{'}, k j} (t - τ) + M_{j^{'}, k^{″} j^{″}} (- τ) + h_{j^{'}, j^{″}} (τ)) . \end{matrix}

(C18)

The first equality of the above equation can be confirmed from the definitions, Eqs. (C3) and (C4), and the third equality results from the following identities: ${(f_{k j, j^{''}} (t) e^{- K_{k j, k^{''} j^{''}} (t)})}^{*} = f_{j k, k^{''}} (t) e^{- K_{j k, j^{''} k^{''}} (t)}$ ⁠; M_j′,kj(τ − t)* = M_j′,kj(t − τ); $M_{j^{'}, k^{''} j^{''}} {(τ)}^{*} = M_{j^{'}, k^{''} j^{''}} (- τ)$ ⁠; and $h_{j^{'}, j^{''}} {(τ)}^{*} = h_{j^{'}, j^{''}} (τ)$ ⁠.

Finally, Eq. (C5) can be expressed as

L_{j, j^{'}}^{j^{″} k^{″}} (t, τ) = 〈 θ_{k^{″}} D_{j} (t) D_{j^{'}} (τ) θ_{j^{″}}^{†} 〉 - w_{j^{″} k^{″}} 〈 D_{j} (t) D_{j^{'}} (τ) 〉 .

(C19)

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

\begin{matrix} 〈 θ_{k^{″}} D_{j} (t) D_{j^{'}} (τ) θ_{j^{″}}^{†} 〉 & = & 〈 D_{j} (t) D_{j^{'}} (τ) θ_{k^{″}} θ_{j^{″}}^{†} 〉 + h_{j, k^{″}} (t) 〈 D_{j^{'}} (τ) \\ \times θ_{k^{″}} θ_{j^{″}}^{†} 〉 + h_{j^{'}, k^{″}} (τ) 〈 D_{j} (t) θ_{k^{″}} θ_{j^{″}}^{†} 〉 \\ + h_{j, k^{″}} (t) h_{j^{'}, k^{″}} (τ) w_{j^{″} k^{″}}, \end{matrix}

(C20)

where the first term can be calculated as follows:

\begin{matrix} 〈 D_{j} (t) D_{j^{'}} (τ) θ_{k^{″}} θ_{j^{″}}^{†} ρ_{b} 〉 & = & \sum_{n} \sum_{n^{'}} ℏ ω_{n} g_{n, j} (1 - W_{h} (ω_{n})) ℏ ω_{n^{'}} g_{n^{'}, j^{'}} \\ \times (1 - W_{h} (ω_{n^{'}})) \{e^{- i (ω_{n} t + ω_{n^{'}} τ)}) \\ \times 〈 b_{n} b_{n^{'}} θ_{k^{″}} θ_{j^{″}}^{†} 〉 + e^{- i (ω_{n} t - ω_{n^{'}} τ)} \\ \times 〈 b_{n} b_{n^{'}}^{†} θ_{k^{″}} θ_{j^{″}}^{†} 〉 + e^{i (ω_{n} t - ω_{n^{'}} τ)} \\ \times 〈 b_{n}^{†} b_{n^{'}} θ_{k^{″}} θ_{j^{″}}^{†} 〉 + e^{- i (ω_{n} t + ω_{n^{'}} τ)} \\ \times (〈 b_{n}^{†} b_{n^{'}}^{†} θ_{k^{″}} θ_{j^{″}}^{†} 〉\} . \end{matrix}

(C21)

In the above expression, the four averages over the bath can be calculated explicitly and can be expressed as

〈 b_{n} b_{n^{'}} θ_{k^{″}} θ_{j^{″}}^{†} 〉 = w_{j^{″} k^{″}} \frac{δ g_{n, j^{″} k^{″}} W_{h} (ω_{n})}{(1 - e^{- β ℏ ω_{n}})} \frac{δ g_{n^{'}, j^{″} k^{″}} W_{h} (ω_{n^{'}})}{(1 - e^{- β ℏ ω_{n^{'}}})},

(C22)

\begin{matrix} 〈 b_{n} b_{n^{'}}^{†} θ_{k^{″}} θ_{j^{″}}^{†} 〉 = w_{j^{″} k^{″}} (\frac{1}{1 - e^{- β ℏ ω_{n}}} δ_{n n^{'}}) \\ - (\frac{δ g_{n, j^{″} k^{″}} W_{h} (ω_{n})}{(1 - e^{- β ℏ ω_{n}})} \frac{δ g_{n^{'}, j^{″} k^{″}} W_{h} (ω_{n^{'}})}{(1 - e^{- β ℏ ω_{n^{'}}})} e^{- β ℏ ω_{n^{'}}}), \end{matrix}

(C23)

\begin{matrix} 〈 b_{n}^{†} b_{n^{'}} θ_{k^{″}} θ_{j^{″}}^{†} 〉 & = & w_{j^{″} k^{″}} (\frac{e^{- β ℏ ω_{n}}}{1 - e^{- β ℏ ω_{n}}} δ_{n n^{'}} - \frac{δ g_{n, j^{″} k^{″}} W_{h} (ω_{n})}{(1 - e^{- β ℏ ω_{n}})}) \\ \times (e^{- β ℏ ω_{n}} \frac{δ g_{n^{'}, j^{″} k^{″}} W_{h} (ω_{n^{'}})}{(1 - e^{- β ℏ ω_{n^{'}}})}), \end{matrix}

(C24)

\begin{align} 〈 b_{n}^{†} b_{n^{'}}^{†} θ_{k^{″}} θ_{j^{″}}^{†} 〉 & = w_{j^{″} k^{″}} \frac{δ g_{n, j^{″} k^{″}} W_{h} (ω_{n})}{(1 - e^{- β ℏ ω_{n}})} e^{- β ℏ ω_{n}} \\ \times \frac{δ g_{n^{'}, j^{″} k^{″}} W_{h} (ω_{n^{'}})}{(1 - e^{- β ℏ ω_{n^{'}}})} e^{- β ℏ ω_{n^{'}}} . \end{align}

(C25)

It is worth noting that the additional factor for n = n′ in Eq. (C23) comes from the first term of the following identity,

\begin{align} 〈 b_{n} b_{n}^{†} e^{γ^{*} b_{n}^{†} - γ b_{n}} 〉 & = (\frac{1}{(1 - e^{- β ℏ ω_{n}})} - | γ |^{2} \frac{e^{- β ℏ ω_{n}}}{{(1 - e^{- β ℏ ω_{n}})}^{2}}) \\ \times e^{- \coth (β ℏ ω_{n} / 2) | γ |^{2} / 2} . \end{align}

(C26)

For the case of (C24), the fact that $b_{n}^{†} b_{n} = b_{n} b_{n}^{†} - 1$ can be combined with the above identity. When Eqs. (C22)–(C25) are inserted into Eq. (C21), the term for n = n′ cancel $- w_{j^{''} k^{''}} 〈 D_{j} (t) D_{j^{'}} (τ) 〉$ in Eq. (C19) and the double summation for n and n′ can be factored as follows:

\begin{align} L_{j, j^{'}}^{j^{″} k^{″}} (t, τ) & = w_{j^{″} k^{″}} (M_{j, j^{″} k^{″}} (t) + h_{j, k^{″}} (t)) (M_{j^{'}, j^{″} k^{″}} (τ)) + (h_{j^{'}, k^{″}} (τ)) . \end{align}

(C27)

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

{\tilde{F}}_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t; E) = \int_{0}^{t} d τ e^{i E (t - τ) / ℏ} F_{j k, j^{'} k^{'}}^{j^{″} k^{″}} (t, τ),

(C28)

{\tilde{H}}_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t; E) = \int_{0}^{t} d τ e^{i E (t - τ) / ℏ} H_{j, j^{'} k^{'}}^{(1), j^{″} k^{″}} (t, τ),

(C29)

{\tilde{H}}_{j k, j^{'}}^{(2), j^{″} k^{″}} (t; E) = \int_{0}^{t} d τ e^{i E (t - τ) / ℏ} H_{j k, j^{'}}^{(2), j^{″} k^{″}} (t, τ),

(C30)

{\tilde{L}}_{j, j^{'}}^{j^{″} k^{″}} (t; E) = \int_{0}^{t} d τ e^{i E (t - τ) / ℏ} L_{j, j^{'}}^{j^{″} k^{″}} (t, τ) .

(C31)

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

REFERENCES

1.

L. D.

Landau

and

S. I.

Pekar

, “

Effective mass of a polaron

,”

Ukr. J. Phys.

53

,

71

–

74

(

2008

).

Google Scholar

2.

H.

Fröhlich

, “

Electrons in lattice fields

,”

Adv. Phys.

3

,

325

(

1954

).

https://doi.org/10.1080/00018735400101213

Google Scholar

Crossref

3.

T.

Holstein

, “

Studies of polaron motion: Part 1. The molecular-crystal model

,”

Ann. Phys.

8

,

325

(

1959

).

https://doi.org/10.1016/0003-4916(59)90002-8

Google Scholar

Crossref

4.

T.

Holstein

, “

Studies of polaron motion: Part 2. The ‘small’ polaron

,”

Ann. Phys.

8

,

343

(

1959

).

https://doi.org/10.1016/0003-4916(59)90003-x

Google Scholar

Crossref

5.

T.

Holstein

, “

Quantal occurrence-probability treatment of small-polaron hopping

,”

Philos. Mag. B

37

,

49

(

1978

).

https://doi.org/10.1080/13642817808245306

Google Scholar

Crossref

6.

V. J.

Emery

and

A.

Luther

, “

Low-temperature properties of the Kondo Hamiltonian

,”

Phys. Rev. B

9

,

215

–

225

(

1974

).

https://doi.org/10.1103/physrevb.9.215

Google Scholar

Crossref

7.

S.

Rackovsky

and

R.

Silbey

, “

Electronic energy transfer in impure solids I. Two molecules embedded in a lattice

,”

Mol. Phys.

25

,

61

(

1973

).

https://doi.org/10.1080/00268977300100081

Google Scholar

Crossref

8.

B.

Jackson

and

R.

Silbey

, “

On the calculation of transfer rate between impurity states in solids

,”

J. Chem. Phys.

78

,

4193

(

1983

).

https://doi.org/10.1063/1.445096

Google Scholar

Crossref

9.

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

).

https://doi.org/10.1063/1.447055

Google Scholar

Crossref

10.

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

(

1985

).

https://doi.org/10.1063/1.449469

Google Scholar

Crossref

11.

R.

Silbey

and

R. A.

Harris

, “

Tunnelingof molecules in low-temperature media: An elementary description

,”

J. Phys. Chem.

93

,

7062

(

1989

).

https://doi.org/10.1021/j100357a010

Google Scholar

Crossref

12.

A.

Nitzan

,

Chemical Dynamics in Condensed Phases

(

Oxford University Press

,

Oxford

,

2006

).

Google Scholar

Crossref

13.

Y.-C.

Cheng

and

R. J.

Silbey

, “

A unified theory for charge-carrier transport in organic crystals

,”

J. Chem. Phys.

128

,

114713

(

2008

).

https://doi.org/10.1063/1.2894840

Google Scholar

Crossref

PubMed

14.

S.

Jang

,

Y.-C.

Cheng

,

D. R.

Reichman

, and

J. D.

Eaves

, “

Theory of coherent resonance energy transfer

,”

J. Chem. Phys.

129

,

101104

(

2008

).

https://doi.org/10.1063/1.2977974

Google Scholar

Crossref

PubMed

15.

S.

Jang

, “

Theory of coherent resonance energy transfer for coherent initial condition

,”

J. Chem. Phys.

131

,

164101

(

2009

).

https://doi.org/10.1063/1.3247899

Google Scholar

Crossref

PubMed

16.

S.

Jang

, “

Theory of multichromophoric coherent resonance energy transfer: A polaronic quantum master equation approach

,”

J. Chem. Phys.

135

,

034105

(

2011

).

https://doi.org/10.1063/1.3608914

Google Scholar

Crossref

PubMed

17.

A.

Nazir

, “

Correlation-dependent coherent to incoherent transitions in resonant energy transfer dynamics

,”

Phys. Rev. Lett.

103

,

146404

(

2009

).

https://doi.org/10.1103/physrevlett.103.146404

Google Scholar

Crossref

PubMed

18.

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

).

https://doi.org/10.1063/1.3636081

Google Scholar

Crossref

PubMed

19.

E. N.

Zimanyi

and

R. J.

Silbey

, “

Theoretical description of quantum effects in multi-chromophoric aggregates

,”

Philos. Trans. R. Soc., A

370

,

3620

(

2012

).

https://doi.org/10.1098/rsta.2011.0204

Google Scholar

Crossref

20.

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

).

https://doi.org/10.1063/1.4732309

Google Scholar

Crossref

PubMed

21.

F. A.

Pollock

,

D. P. S.

McCutcheon

,

B. W.

Lovett

,

E. M.

Gauger

, and

A.

Nazir

, “

A multi-site variational master equation approach to dissipative energy transfer

,”

New J. Phys.

15

,

075018

(

2013

).

https://doi.org/10.1088/1367-2630/15/7/075018

Google Scholar

Crossref

22.

A.

Nazir

and

D. P. S.

McCutcheon

, “

Modelling exciton-phonon interactions in optically driven quantum dots

,”

J. Phys.: Condens. Matter

28

,

103002

(

2016

).

https://doi.org/10.1088/0953-8984/28/10/103002

Google Scholar

Crossref

PubMed

23.

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

).

https://doi.org/10.1063/1.4789017

Google Scholar

Crossref

PubMed

24.

D.

Chen

,

J.

Ye

,

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

).

https://doi.org/10.1021/jp1093492

Google Scholar

Crossref

PubMed

25.

Y.

Zhao

,

B.

Luo

,

Y.

Zhang

, and

J.

Ye

, “

Dynamics of a Holstein polaron with off-diagonal coupling

,”

J. Chem. Phys.

137

,

084113

(

2012

).

https://doi.org/10.1063/1.4748140

Google Scholar

Crossref

PubMed

26.

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

).

https://doi.org/10.1063/1.4884275

Google Scholar

Crossref

PubMed

27.

P.

Hamm

and

G. P.

Tsironis

, “

Barrier crossing to the small Holstein polaron regime

,”

Phys. Rev. B

78

,

092301

(

2008

).

https://doi.org/10.1103/physrevb.78.092301

Google Scholar

Crossref

28.

C. K.

Lee

,

J.

Moix

, and

J.

Cao

, “

Coherent quantum transport in disordered systems: A unified polaron treatment of hopping and band-like transport

,”

J. Chem. Phys.

142

,

164103

(

2015

).

https://doi.org/10.1063/1.4918736

Google Scholar

Crossref

PubMed

29.

D.

Xu

and

J.

Cao

, “

Non-canonical distribution and non-equilibrium transport beyond weak system-bath coupling regime: A polaron transformation approach

,”

Front. Phys.

11

,

110308

(

2016

).

https://doi.org/10.1007/s11467-016-0540-2

Google Scholar

Crossref

30.

Y.-C.

Wang

and

Y.

Zhao

, “

Variational polaron transformation approach toward the calculation of thermopower in organic crystals

,”

Phys. Rev. B

101

,

075205

(

2020

).

https://doi.org/10.1103/physrevb.101.075205

Google Scholar

Crossref

31.

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

).

https://doi.org/10.1039/d0sc04116e

Google Scholar

Crossref

32.

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

).

https://doi.org/10.1063/1.5096287

Google Scholar

Crossref

PubMed

33.

A. O.

Caldeira

and

A. J.

Leggett

, “

Quantum tunnelling in a dissipative system

,”

Ann. Phys.

149

,

374

(

1983

).

https://doi.org/10.1016/0003-4916(83)90202-6

Google Scholar

Crossref

34.

U.

Weiss

,

Quantum Dissipative Systems

, Series in Modern Condensed Matter Physics Vol. 2 (

World Scientific

,

Singapore

,

1993

).

Google Scholar

Crossref

35.

V. M.

Kenkre

and

P.

Reineker

,

Exciton Dynamics in Molecular Crystals and Aggregates

(

Springer

,

Berlin

,

1982

).

Google Scholar

Crossref

36.

V.

May

and

O.

Kühn

,

Charge and Energy Transfer Dynamics in Molecular Systems

(

Wiley VCH

,

Weinheim

,

2011

).

Google Scholar

Crossref

37.

S. J.

Jang

,

Dynamics of Molecular Excitons

, Nanophotonics Series (

Elsevier

,

Amsterdam

,

2020

).

Google Scholar

Crossref

38.

V.

Coropceanu

,

J.

Cornil

,

D. A.

da Silva Filho

,

Y.

Olivier

,

R.

Silbey

, and

J.-L.

Brédas

, “

Charge transport in organic semiconductors

,”

Chem. Rev.

107

,

926

–

952

(

2007

).

https://doi.org/10.1021/cr050140x

Google Scholar

Crossref

PubMed

39.

S.

Jang

,

J.

Cao

, and

R. J.

Silbey

, “

Fourth order quantum master equation and its Markovian bath limit

,”

J. Chem. Phys.

116

,

2705

(

2002

).

https://doi.org/10.1063/1.1445105

Google Scholar

Crossref

40.

F.

Shibata

and

T.

Arimitsu

, “

Expansion formulas in nonequilibrium statistical mechanics

,”

J. Phys. Soc. Jpn.

49

,

891

(

1980

).

https://doi.org/10.1143/jpsj.49.891

Google Scholar

Crossref

41.

A. J.

Leggett

,

S.

Chakravarty

,

A. T.

Dorsey

,

M. P. A.

Fisher

,

A.

Garg

, and

W.

Zwerger

, “

Dynamics of the dissipative two-state system

,”

Rev. Mod. Phys.

59

,

1

–

85

(

1987

).

https://doi.org/10.1103/revmodphys.59.1

Google Scholar

Crossref

42.

C.

Aslangul

,

N.

Pottier

, and

D.

Saint-James

, “

Quantum ohmic dissipation: Transition from coherent to incoherent dynamics

,”

Phys. Lett. A

110

,

249

(

1985

).

https://doi.org/10.1016/0375-9601(85)90091-x

Google Scholar

Crossref

43.

H.

Dekker

, “

Noninteracting-blip approximation for a two-level system coupled to heat bath

,”

Phys. Rev. A

35

,

1436

–

1437

(

1987

).

https://doi.org/10.1103/physreva.35.1436

Google Scholar

Crossref

44.

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

).

https://doi.org/10.1063/5.0069313

Google Scholar

Crossref

PubMed

45.

W.

Weibull

, “

A statistical distribution function of wide applicability

,”

ASME J. Appl. Mech.

18

,

293

–

297

(

1951

).

https://doi.org/10.1115/1.4010337

Google Scholar

Crossref

2022

Author(s)

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

I. INTRODUCTION

II. THEORY

A. Hamiltonian and partial polaron transformation

B. Partially polaron transformed quantum master equation (p-PQME)

1. Quantum Liouville equation in the polaron and interaction picture

2. Quantum master equation for the reduced density operator

3. Homogeneous terms of the second order time local p-PQME

4. Inhomogeneous terms of p-PQME

C. Representation in the basis of renormalized system eigenstates

III. TWO-LEVEL SYSTEM WITH INDEPENDENT BATH

A. Model and general expressions

B. Model calculations without inhomogeneous terms for Ohmic spectral density

IV. CONCLUDING REMARKS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: DERIVATION OF Eq. (41)

APPENDIX B: EVALUATION OF THE BATH CORRELATION FUNCTION IN THE FIRST ORDER INHOMOGENEOUS TERM

APPENDIX C: EVALUATION OF THE BATH CORRELATION FUNCTION IN THE SECOND ORDER INHOMOGENEOUS TERM

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

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

I. INTRODUCTION

II. THEORY

A. Hamiltonian and partial polaron transformation

B. Partially polaron transformed quantum master equation (p-PQME)

1. Quantum Liouville equation in the polaron and interaction picture

2. Quantum master equation for the reduced density operator

3. Homogeneous terms of the second order time local p-PQME

4. Inhomogeneous terms of p-PQME

C. Representation in the basis of renormalized system eigenstates

III. TWO-LEVEL SYSTEM WITH INDEPENDENT BATH

A. Model and general expressions

B. Model calculations without inhomogeneous terms for Ohmic spectral density

IV. CONCLUDING REMARKS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: DERIVATION OF Eq. (41)

APPENDIX B: EVALUATION OF THE BATH CORRELATION FUNCTION IN THE FIRST ORDER INHOMOGENEOUS TERM

APPENDIX C: EVALUATION OF THE BATH CORRELATION FUNCTION IN THE SECOND ORDER INHOMOGENEOUS TERM

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only