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.
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 PQMEs14–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 PQME18,19,21,22 that combines a variational ansatz6,9–11 with PQME and, more recently, frozen mode PQME32 approaches have been developed to address this issue. This work develops a formulation that can put these works18,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.
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 Ej is the energy of the site local state |j⟩ and Jjk, 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 Jjk 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:
In the above expressions, ωn and bn are the frequency and the lowering (raising) operator of each normal mode constituting the bath degrees of freedom, and gn,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 Hamiltonian33,34 and has been used widely for molecular excitons35–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 Hb and Hsb 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 Wh(ω) 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
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:
In the above expression, 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), , is a partially renormalized system–bath interaction Hamiltonian given by
where can be expressed as
with δgn,jk = (gn,j − gn,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:
In above expressions, with . 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, , includes the average effect of system–bath interactions in the partial polaron picture and can be expressed as
where has been defined by Eq. (14) and with
Unlike the case of the full PT, the Debye–Waller factor wjk 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 , Eq. (20), can be expressed as
In the above expressions δjk is the Kronecker-delta symbol and
Thus, the bath operator 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 , for which exact time evolution can be implemented numerically, let us now consider the dynamics in the interaction picture of . First, in this interaction picture becomes
In the above expression,
In the interaction picture with respect to , the partially polaron-transformed total density operator becomes , which evolves according to the following time dependent quantum Liouville equation:
where the second equality serves as the definition of .
2. Quantum master equation for the reduced density operator
Taking trace of 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 employing the standard projection operator technique37,39 for a well-known projection operator as follows:
where e(+) denotes the exponential operator with chronological time ordering, and
Alternatively, replacing with the backpropagation of 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
In the above expression, and 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 with , 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 work44 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 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 (with subscript b omitted) represents averaging over the equilibrium bath density operator, ρb.
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 wjk defined by Eq. (22) can be expressed as follows:
If Wh(ω) = 1, and become zero and the expressions for and wjk reduce to those for the original second order PQME based on the full PT. On the other hand, for Wh(ω) = 0, and become zero and the expression for 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 Wh(ω) is a step function. For this case, Wh(ω)(1 − Wh(ω)) = 0 and, thus, 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 Eq. (27) with t = 0 has been used and
with the convention that wjk = 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 and , Eqs. (28) and (57), respectively. Details of this calculation are provided in Appendix B, and the resulting expression can be summed up as
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 Wh(ω) = 1 and those for the regular second order QME for Wh(ω) = 0. For the case where Wh(ω)(1 − Wh(ω)) = 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 . Let us denote the pth eigenstate and eigenvalue of as |φp⟩ and . Then,
The transformation matrix U between the site localized states |j⟩s and the eigenstates |φp⟩s can be defined such that Ujp = ⟨j|φp⟩. Then,
This transformation can be used to express defined by Eq. (27) in the basis of |φp⟩s as follows:
where . Let us also introduce Spq(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 in the basis of |φp⟩s. For a more compact expression, let us introduce
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.
which can also be obtained from Eq. (67) by replacing Sp′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:
Combining the above expression with bath correlation functions introduced in Appendix C, one can express the second order inhomogeneous term as follows:
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 = J12 = J21 and only one Debye–Waller factor, w = w12 = w21. Then, introducing 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 D1 and D2 have been defined by Eq. (25). Note also that, for the present case,
and by definition. The eigenvalues of 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 . The corresponding eigenstates are given by
Thus, for the present case, U11 = U22 = cos ξ and U21 = −U12 = sin ξ.
Let us also assume that the spectral densities for sites 1 and 2 are the same, which we denote as . Namely, . In addition, let us also assume that . Then, and . For all other indices, and .
Then, all the bath correlation functions that enter the p-PQME can be represented by the following three functions:
For all other components, , , and .
In the site basis, Eq. (40) for the present case can be expressed as follows:
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 , as detailed in Sec. II C. For this, given by Eq. (72) needs to be calculated, which in turn requires calculation of , , and 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 ’s defined by Eq. (68) for the present case reduce to one of the following two functions:
All other terms of are zero by definition. Similarly, all of nonzero ’s and ’s can be specified by
All other terms with indices different from the above are zero. Thus, with , , and as defined above, all of constituting Eq. (72) for the present case can be calculated.
The first order inhomogeneous term, Eq. (74), involves additional functions fjk,k′(t) and hj,k′(t), which can also be simplified for the present model. Let us define
Then, it is easy to show that f21,2(t) = f(t), f21,1(t) = f12,2(t) = f*(t). On the other hand, f11,1(t) = f22,1(t) = f11,2(t) = f22,2(t) = 1. For hj,k′(t), only specification of the following function is necessary:
For other cases, h1,2(t) = h2,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 Wh(ω). 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 ’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:
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, E1 = E2 and E1 − E2 = 200 cm−1 will be considered.
There are infinite number of possible choices available for Wh(ω) that satisfies the requirement of Eq. (7). One simple but flexible choice is the following function known as (cumulative) Weibull distribution45 or the complement of a stretched exponential function,
For small ω, . For large α, Wh(ω) 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 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 and , defined, respectively, by Eqs. (89) and (90), converge to finite values.
On the other hand, for the case of , 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 , 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 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 to converge to a finite value in the steady state limit.
For the case of , as shown in Fig. 3, both real and imaginary parts decay to zero quickly already for α = 1, resulting in well behaving . Summing up the results shown in Figs. 1–3, while the choice of α = 1 is acceptable as far as the Deby–Waller factor, , and are concerned, it is not appropriate due to resulting slow decay of . 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 E1 = E2 for different values of ωh with the choice of α = 2, where Wh(ω) 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 E1 − E2 = 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 Wh(ω) 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 Wh(ω) 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.
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.
Conflict of Interest
The author has no conflicts to disclose.
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).
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 with ρb in Eq. (40), namely, its thermal average is expressed as follows:
In obtaining the second equality of the above equation, the identities , , and ⟨Dj(t)⟩ = ⟨Dj′(τ)⟩ = 0 have been used.
The first term in the second equality of Eq. (A1) has the same form as that for the full PT15 except for the additional factor Wh(ωn) multiplied with δgn,jk and δgn,j′k′, respectively. Thus, following the same procedure as in previous work,15 it can be shown to be
where 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 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:
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 . For this, and need to be calculated. For the first of these terms, the following identity can be used:
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 δgn,kj = −δgn,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 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,
The above identity implies that the first term of Eq. (B1) (without Jjk) 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 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:
Further calculation of each of the terms above is straightforward as described below.
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 . Combining all of these, one can show that
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:
The second bath average term on the right-hand side of Eq. (C9) can be shown to be
defined by Eq. (C4) can be calculated in a similar manner but can, in fact, be calculated using its relation to 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: ; Mj′,kj(τ − t)* = Mj′,kj(t − τ); ; and .
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 can be combined with the above identity. When Eqs. (C22)–(C25) are inserted into Eq. (C21), the term for n = n′ cancel 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 .