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^{32} 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,^{15} 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

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:

where

In the above expressions, *ω*_{n} and *b*_{n}$(bn\u2020)$ 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:

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:

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.^{34} 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:

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:

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

where

On the other hand, it is straightforward to show that

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

where

In the above expression, $E\u0303j$ is a partially renormalized energy for site *j* given by

with *λ*_{j} being the corresponding reorganization energy defined as

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), $H\u0303sb$, is a partially renormalized system–bath interaction Hamiltonian given by

where $\theta j\u2020\theta k$ can be expressed as

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:

where

In above expressions, $\u3008H\u0303sb\u3009b=Trb{H\u0303sb\rho b}$ with $\rho b=e\u2212\beta Hb/Trb{e\u2212\beta Hb}$. 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, $H\u03030,s=H\u0303sp+\u3008H\u0303sc\u3009b$, includes the average effect of system–bath interactions in the partial polaron picture and can be expressed as

where $E\u0303j$ has been defined by Eq. (14) and $J\u0303jk=wjkJjk$ with

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 $H\u03031$, Eq. (20), can be expressed as

where

In the above expressions *δ*_{jk} is the Kronecker-delta symbol and

Thus, the bath operator $B\u0303jk$ 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 $H\u03030$, for which exact time evolution can be implemented numerically, let us now consider the dynamics in the interaction picture of $H\u03030$. First, $H\u03031$ in this interaction picture becomes

where

In the above expression,

In the interaction picture with respect to $H\u03030$, the partially polaron-transformed total density operator becomes $\rho \u0303I(t)=eiL\u03030t\rho \u0303(t)$, which evolves according to the following time dependent quantum Liouville equation:

where the second equality serves as the definition of $L\u03031,I(t)$.

#### 2. Quantum master equation for the reduced density operator

Taking trace of $\rho \u0303I(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:

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 $\sigma \u0303I(t)$ employing the standard projection operator technique^{37,39} for a well-known projection operator $P(\u22c5)\u2261\rho bTrb{(\u22c5)}$ as follows:

where *e*_{(+)} denotes the exponential operator with chronological time ordering, $Q=1\u2212P$ and

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

where *e*_{(−)} denotes the exponential operator with anti-chronological time ordering and

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

where

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 $\sigma \u0303I(\tau )$ with $\sigma \u0303I(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^{44} 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)\sigma \u0303I(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:

where H.c refers to Hermitian conjugates of all previous terms and $\u3008B\u0303jk(t)B\u0303j\u2032k\u2032(\tau )\u3009$ (with subscript *b* omitted) represents averaging over the equilibrium bath density operator, *ρ*_{b}.

Appendix A describes calculation of all the terms constituting $\u3008B\u0303jk(t)B\u0303j\u2032k\u2032(\tau )\u3009$. When the resulting expressions, Eqs. (A2), (A15), (A21), and (A22), are used in Eq. (A1), it can be expressed as

where

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

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,

Note also that *w*_{jk} defined by Eq. (22) can be expressed as follows:

If *W*_{h}(*ω*) = 1, $Mj,j\u2032k\u2032(t)$ and $Cjj\u2032(t)$ become zero and the expressions for $Kjk,j\u2032k\u2032(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, $Kjk,j\u2032k\u2032(t)$ and $Cjj\u2032(t)$ become zero and the expression for $Cjj\u2032(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, $Mj,j\u2032k\u2032(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

where

Then,

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

with the convention that *w*_{jk} = 1 for *j* = *k*. Therefore, the first order inhomogeneous term in Eq. (39) can be expressed as

In the above expression, the trace over the bath can be calculated explicitly employing the definitions for $B\u0303jk(t)$ and $\delta \rho \u0303b,jk$, Eqs. (28) and (57), respectively. Details of this calculation are provided in Appendix B, and the resulting expression can be summed up as

where

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

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 $H\u03030,s$. Let us denote the *p*th eigenstate and eigenvalue of $H\u03030,s$ as |*φ*_{p}⟩ and $Ep$. Then,

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,

This transformation can be used to express $Tjk(t)$ defined by Eq. (27) in the basis of |*φ*_{p}⟩s as follows:

where $\delta Epq=Ep\u2212Eq$. Let us also introduce *S*_{pq}(*t*)s such that

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

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

Then,

where

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

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:

Similar expressions for the second order inhomogeneous term, Eq. (62), can be obtained by replacing *S*_{p′q′}(*t*) in Eq. (67) with $Uj\u2032\u2032p\u2032*Uk\u2032\u2032q\u2032$. The resulting expression is as follows:

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

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*_{12} = *J*_{21} and only one Debye–Waller factor, *w* = *w*_{12} = *w*_{21}. Then, introducing $\theta =\theta 1\u2020\theta 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

where *w* = ⟨*θ*⟩ = ⟨*θ*^{†}⟩ and *D*_{1} and *D*_{2} have been defined by Eq. (25). Note also that, for the present case,

and $J\u030311=J\u030322=0$ by definition. The eigenvalues of $H\u03030,s$ are given by

with subscripts 1 and 2 on the left hand side denoting and denoting the + sign, respectively − signs on the right-hand side, and $\xi =tan\u22121(2Jw/(E\u03031\u2212E\u03032))/2$. The corresponding eigenstates are given by

Thus, for the present case, *U*_{11} = *U*_{22} = cos *ξ* and *U*_{21} = −*U*_{12} = sin *ξ*.

Let us also assume that the spectral densities for sites 1 and 2 are the same, which we denote as $J(\omega )$. Namely, $J11(\omega )=J22(\omega )=J(\omega )$. In addition, let us also assume that $J12(\omega )=J21(\omega )=0$. Then, $J1,12(1)(\omega )=J2,21(1)(\omega )=J(\omega )$ and $J12,12(2)(\omega )=J21,21(2)(\omega )=2J(\omega )$. For all other indices, $Jj,j\u2032k\u2032(1)(\omega )=0$ and $Jjk,j\u2032k\u2032(2)(\omega )=0$.

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

For all other components, $Kjk,j\u2032k\u2032(\omega )=0$, $Mj,j\u2032k\u2032(\omega )=0$, and $Cjk(\omega )=0$.

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

where

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 $H\u03030,s$, as detailed in Sec. II C. For this, $Rpqp\u2032q\u2032(t)$ given by Eq. (72) needs to be calculated, which in turn requires calculation of $Wjk,j\u2032k\u2032pq(t)$, $Yj,j\u2032k\u2032pq(t)$, and $Xjj\u2032pq(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 $Wjk,j\u2032k\u2032pq(t)$’s defined by Eq. (68) for the present case reduce to one of the following two functions:

All other terms of $Wjk,j\u2032k\u2032pq(t)$ are zero by definition. Similarly, all of nonzero $Yj,j\u2032k\u2032pq(t)$’s and $Xjj\u2032pq(t)$’s can be specified by

All other terms with indices different from the above are zero. Thus, with $W\xb1pq(t)$, $Ypq(t)$, and $Xpq(t)$ as defined above, all of $Rpqp\u2032q\u2032(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

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:

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 $Rpqp\u2032q\u2032(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:

where

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

It is assumed that the two parameters of the above spectral density are set to *η* = 1 and *ℏω*_{c} = 200 cm^{−1}. In addition, the temperature and electronic coupling are set to *T* = 300 K and *J* = 300 cm^{−1}. 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*_{1} = *E*_{2} and *E*_{1} − *E*_{2} = 200 cm^{−1} 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^{45} or the complement of a stretched exponential function,

For small *ω*, $Wh(\omega )\u2248(\omega /\omega h)\alpha $. 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\u2212pq(t)$ and $W+pq(t)$, defined, respectively, by Eqs. (89) and (90), converge to finite values.

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 $Ypq(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 $Ypq(t)$ to converge to a finite value in the steady state limit.

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 $Xpq(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.

Figure 4 shows the results for *E*_{1} = *E*_{2} 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.

Figure 5 shows the results for an asymmetric case, where *E*_{1} − *E*_{2} = 200 cm^{−1}. 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.

## 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 $B\u0303jk(t)B\u0303j\u2032k\u2032(\tau )$ with *ρ*_{b} in Eq. (40), namely, its thermal average is expressed as follows:

In obtaining the second equality of the above equation, the identities $\u3008\theta j\u2032\u2020(\tau )\theta k\u2032(\tau )\u3009=wj\u2032k\u2032$, $\u3008\theta j\u2020(t)\theta k(t)\u3009=wjk$, 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^{15} 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,^{15} it can be shown to be

where $eKjk,j\u2032k\u2032(t\u2212\tau )$ in the exponent is defined by Eq. (42).

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

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

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

where *γ* is defined by

and *γ** is the complex conjugate of *γ*. Using the fact that $e\gamma *bn\u2020\u2212\gamma bn=e\u2212\gamma bne\gamma *bn\u2020e|\gamma |2/2$ and the following identity,

it is straightforward to show that

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

Using the following identity,

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

Now, employing the following identity:

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

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

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

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 $\gamma =\delta gn,jkWh(\omega n)e\u2212i\omega nt$. For this, $\u3008e\gamma *bn\u2020\u2212\gamma bnbn\u3009$ and $\u3008e\gamma *bn\u2020\u2212\gamma bnbn\u2020\u3009$ need to be calculated. For the first of these terms, the following identity can be used:

Thus,

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

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

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

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

In the second equality of the above equation, the fact that $\u3008bn\u2020bn\u3009=\u3008bnbn\u2020\u3009\u22121=e\u2212\beta \u210f\omega n/(1\u2212e\u2212\beta \u210f\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. (58), can be expressed as follows:

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

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

Therefore,

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

The above identity implies that the first term of Eq. (B1) (without *J*_{jk}) can be simplified to

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

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

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:

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:

where

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

The first bath term in Eq. (C1), $Fjk,j\u2032k\u2032j\u2032\u2032k\u2032\u2032(t,\tau )$ defined by Eq. (C2), can be expanded further and expressed as follows:

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

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

The bath term contributing to the second term of Eq. (C1), i.e., $Hjk,j\u2032k\u2032(1),j\u2032\u2032k\u2032\u2032(t,\tau )$ defined by Eq. (C3), is expressed as follows:

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

where the following identities have been used,

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:

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:

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

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

$Hjk,j\u2032(2),j\u2032\u2032k\u2032\u2032(t,\tau )$ defined by Eq. (C4) can be calculated in a similar manner but can, in fact, be calculated using its relation to $Hj\u2032,kj(1),k\u2032\u2032j\u2032\u2032(\tau ,t)$ as follows:

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: $fkj,j\u2032\u2032(t)e\u2212Kkj,k\u2032\u2032j\u2032\u2032(t)*=fjk,k\u2032\u2032(t)e\u2212Kjk,j\u2032\u2032k\u2032\u2032(t)$; *M*_{j′,kj}(*τ* − *t*)* = *M*_{j′,kj}(*t* − *τ*); $Mj\u2032,k\u2032\u2032j\u2032\u2032(\tau )*=Mj\u2032,k\u2032\u2032j\u2032\u2032(\u2212\tau )$; and $hj\u2032,j\u2032\u2032(\tau )*=hj\u2032,j\u2032\u2032(\tau )$.

Finally, Eq. (C5) can be expressed as

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

where the first term can be calculated as follows:

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

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

For the case of (C24), the fact that $bn\u2020bn=bnbn\u2020\u22121$ can be combined with the above identity. When Eqs. (C22)–(C25) are inserted into Eq. (C21), the term for *n* = *n*′ cancel $\u2212wj\u2032\u2032k\u2032\u2032\u3008Dj(t)Dj\u2032(\tau )\u3009$ in Eq. (C19) and the double summation for *n* and *n*′ can be factored as follows:

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

The above time correlation functions are useful for expressing the second order inhomogeneous terms in the basis of eigenstates of $H\u0303s$.