Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities

For a system strongly coupled to a heat bath, the quantum coherence of the system and the heat bath plays an important role in the system dynamics. This is particularly true in the case of non-Markovian noise. We rigorously investigate the influence of system-bath coherence by deriving the reduced hierarchal equations of motion (HEOM), not only in real time, but also in imaginary time, which represents an inverse temperature. It is shown that the HEOM in real time obtained when we include the system-bath coherence of the initial thermal equilibrium state possess the same form as those obtained from a factorized initial state. We find that the difference in behavior of systems treated in these two manners results from the difference in initial conditions of the HEOM elements, which are defined in path integral form. We also derive HEOM along the imaginary time path to obtain the thermal equilibrium state of a system strongly coupled to a non-Markovian bath. Then, we show that the steady state hierarchy elements calculated from the real-time HEOM can be expressed in terms of the hierarchy elements calculated from the imaginary-time HEOM. Moreover, we find that the imaginary-time HEOM allow us to evaluate a number of thermodynamic variables, including the free energy, entropy, internal energy, heat capacity, and susceptibility. The expectation values of the system energy and system-bath interaction energy in the thermal equilibrium state are also evaluated.

For a system strongly coupled to a heat bath, the quantum coherence of the system and the heat bath plays an important role in the system dynamics. This is particularly true in the case of non-Markovian noise. We rigorously investigate the influence of system-bath coherence by deriving the reduced hierarchal equations of motion (HEOM), not only in real time, but also in imaginary time, which represents an inverse temperature. It is shown that the HEOM in real time obtained when we include the system-bath coherence of the initial thermal equilibrium state possess the same form as those obtained from a factorized initial state. We find that the difference in behavior of systems treated in these two manners results from the difference in initial conditions of the HEOM elements, which are defined in path integral form. We also derive HEOM along the imaginary time path to obtain the thermal equilibrium state of a system strongly coupled to a non-Markovian bath. Then, we show that the steady state hierarchy elements calculated from the real-time HEOM can be expressed in terms of the hierarchy elements calculated from the imaginarytime HEOM. Moreover, we find that the imaginary-time HEOM allow us to evaluate a number of thermodynamic variables, including the free energy, entropy, internal energy, heat capacity, and susceptibility. The expectation values of the system energy and system-bath interaction energy in the thermal equilibrium state are also evaluated.

I. INTRODUCTION
Quantum open systems have been a subject of fundamental interest for many years.
Problems in this category include those of understanding how the irreversibility of time appears in system dynamics, why macroscopic systems can be treated with classical mechanics instead of quantum mechanics, how wave functions collapse as a result of measurements done with macroscopic instruments, and why and how quantum systems approach a thermal equilibrium state through interaction with their environments. [1][2][3][4] Theories of quantum open systems have also been used to construct models of practical interest, in particular to account for line shapes in EPR, NMR 5,6 and laser spectra, 7 to evaluate chemical reaction rates 8 and electron and charge transfer rates 9,10 in chemical physics, and to explore the lifetimes of quantum entanglement states in quantum information theory. 11 The phenomena mentioned above arise from the unavoidable interaction of a system with its environment. In the quantum mechanical case, dissipative systems are often modeled as main systems coupled to heat-bath degrees of freedom at finite temperature. This coupling gives rise to thermal fluctuations and dissipation that drive the systems toward the thermal equilibrium state. The heat-bath degrees of freedom are then reduced using such methods as the projection operator method and the path integral method.
The projection operator approach is effective if the interaction between the system and the bath is weak. If one further assumes that the correlation time of the noise arising from the interaction with the bath is very short (the Markovian assumption), equations of motion for the density matrix elements can be derived, and these can be solved numerically. The most commonly used equations of this kind are the quantum master equations 3,4 and the Redfield equation. 5 It has been proven, however, that quantum master equations and the Redfield equation do not satisfy the necessary positivity condition. [12][13][14][15] Careful analyses addressing this problem have been carried out by several researchers. 16,17 As a method to preserve positivity, the rotating wave approximation (RWA), which modifies the interaction between the system and the heat bath, has been applied. 18,19 However, this may alter the thermal equilibrium state as well as the dynamics of the original total Hamiltonian. Then, building on these results, it was shown that the violation of the positivity condition results from the Markovian assumption. Specifically, it was found that even if the dissipation process is Markovian, the fluctuation process may not be, because it must satisfy the fluctuation-dissipation theorem. 20 The time convolution-less (TCL) master equation has a wider range of applicability than the quantum master equations and Redfield equation, because it allows the system-bath interaction to be non-perturbative and fluctuation and dissipation to be non-Markovian. 21,22 In order for a non-perturbative treatment to be possible, however, the system Hamiltonian and the bath interactions of the TCL equation must commute. 23,24 Thus, the TCL equation cannot be used to treat many important problems involving molecules, atoms, and spins driven by a time-dependent laser or magnetic field. In addition, because of the factorized nature of the system-bath interaction, in the case of a non-commuting excitation, the TCL equation cannot be used to calculate nonlinear response functions of the system operator involved in the optical multidimensional spectrum. 25,26 Path integral Monte Carlo simulations do not have the limitations of any of the approaches discussed above, and for this reason, they are capable of incorporating imaginary path integrals and unfactorized initial conditions more easily, but this approach is computationally heavy, because the number of paths to be evaluated grows rapidly with time, while sampling fails due to the phase cancellation of wave functions. [27][28][29] Much effort has been made to extend the applicability of this method. [30][31][32][33][34][35] Because this approach can easily incorporate the semi-classical approximation in the bath, it may have an advantage in the study of polyatomic systems treated in multi-dimensional coordinates, but applications to this point incorporating full quantum dynamics have been limited to relatively small systems.
Many of the above-mentioned limitations can be overcome with the hierarchal equations of motion (HEOM) for the reduced density matrix, which are derived by differentiating the reduced density matrix elements defined by path integrals. 20 This approach was introduced to investigate the connection between the phenomenological stochastic Liouville equation and the dynamical Hamiltonian theory, and was originally limited to the case in which the spectral distribution function takes the Drude form (i.e., the Ohmic form with a Lorentzian cutoff) and the bath temperature is high. 36 However, with the inclusion of low temperature corrections terms, this temperature limitation has been eliminated. [37][38][39][40] In addition, with the extension of the dimension of the hierarchy, this approach is capable of treating a great variety of spectral distribution functions. [41][42][43][44][45][46][47][48] This formalism is valuable because it can treat not only strong system-bath coupling but also quantum coherence between the system and bath, which is essential to calculate nonlinear response functions. The system-bath coherence becomes particularly important if the bath interaction is regarded as non-Markovian, as was found from nonlinear optical measurements in the late 1980s, when laser technology reached the femto-second time scale, which is much shorter than the noise correlation time of environmental molecules. 7 The HEOM approach has been used to study such problems, which include multi-dimensional spectroscopy. [49][50][51][52][53][54][55][56][57][58][59][60] Recently, it was shown that system-bath coherence also plays an important role in calculations of quantum measures involving concurrence 11 and non-Markovianity 61 under multiple external perturbations. [62][63][64] Because the HEOM approach is computationally heavy, a variety of methods have been developed to study dissipative dynamics in realistic situations. [65][66][67][68][69][70][71][72] It has been applied to the study of multi-dimensional vibrational spectroscopies, 52-55 photosynthetic antenna systems, 58-60,73-76 fermion systems, 77-79 quantum ratchets, 80 resonant tunneling diodes, 81,82 and dissociation of tightly bounded electron-hole pairs. 83 While the applicability of the HEOM approach continues to expand, the basic nature of the hierarchy elements has not been thoroughly explored. The purpose of this paper is to investigate the role of correlated initial equilibrium states in the HEOM formalism. Until this time, the HEOM have been derived by assuming a factorized initial state, , at inverse temperature β, whereĤ A andĤ B are the system and bath Hamiltonians, respectively, while the true thermal equilibrium state of the system is given byρ tot = exp[−β(Ĥ A +Ĥ I +Ĥ B )], whereĤ I is the system-bath interaction. The difference between the factorized and correlated initial states becomes large for strongĤ I . Analysis based on an analytic solution of a Brownian oscillator system indicates that even if we start from a factorized initial state, the system reaches the true equilibrium state, tr B {exp[−β(Ĥ A +Ĥ I +Ĥ B )]}, through transient phenomena arising from the factorized initial state, for example phenomena known as initial sweeping. [84][85][86] With the HEOM approach, we have run the HEOM program until all of the hierarchy elements reach the steady state and then used these elements as the initial conditions of the correlated thermal equilibrium state. The accuracy of this method has been confirmed by analyzing multi-dimensional spectra obtained with it. 53 Nevertheless, it would be interesting to derive the HEOM starting from a correlated initial thermal state in order to obtain an analytically derived expression for the system-bath coherence in the HEOM formalism. Moreover, with a simple generalization, we can also derive the HEOM in imaginary time, which corresponds to the inverse temperature. We show that the imaginary-time HEOM is convenient for obtaining correlated thermal equilibrium states and the thermodynamic variables of the reduced system. The organization of the paper is as follows. In Sec. II we present a model Hamiltonian and its influence functional with correlated initial states. In Sec. III, we derive the HEOM from the density matrix elements using the influence functional formalism with the correlated initial states given in Sec. II. In Sec. IV, we derive the imaginary-time HEOM, which is convenient for evaluating correlated thermal equilibrium states and the thermodynamic quantities of the system. In Sec. V, to confirm the validity and numerical efficiency of our approach, we report the results of numerical integrations of the HEOM carried out over real time and imaginary time for a spin-boson system and compare their results. Thermodynamic variables and expectation values for the spin-boson system are also calculated as a demonstration. Section VI is devoted to concluding remarks.

II. INFLUENCE FUNCTIONAL WITH CORRELATED INITIAL STATES
We consider a situation in which the system interacts with a heat bath that gives rise to dissipation and fluctuation in the system. To illustrate this, let us consider a Hamiltonian expressed asĤ whereĤ A ≡ H A (â + ,â − ) is the Hamiltonian of the system, denoted by A, defined by the creation and annihilation operatorsâ + andâ − . The bath degrees of freedom are treated as an ensemble of harmonic oscillators, with the momentum, position, mass, and frequency of the jth bath oscillator given byp j , x j , m j and ω j , respectively. The system-bath interaction is given bŷ whereV (â + ,â − ) is the system part of the interaction, and α j is the coupling constant between the system and the jth oscillator.
The heat bath can be characterized by the spectral distribution function, defined by and the inverse temperature, β ≡ 1/k B T , where k B is the Boltzmann constant. Note that if a + andâ − respectively represent the creation and annihilation operators of spin states, the above Hamiltonian is the spin-boson Hamiltonian, 2,3 which has been studied with various approaches. [5][6][7][27][28][29] Now, let us introduce the fermion coherent state |φ , which satisfiesâ − |φ = φ|φ and φ|â + = φ † φ|, where φ and φ † are Grassmann numbers (G-numbers). [36][37][38] Note that here we consider a two-level system, but extension to a multi-level system is also straightforward. 42,43 In practice, we can treat a G-number system in the same manner as a c-number system, as long as we maintain the time order of the operators in the integral. In the path integral representation, the time propagator of the wave function (the Feynman propagator) for the total system is expressed as where N is the normalization constant, D[φ † (τ )φ(τ )] represents a functional integral over a set of Grassmann variables, and D[x(τ )] ≡ Π j D[x j (τ )] denotes path integrals over the bath oscillator coordinates, with mẋ 2 ≡ j m jẋ Here, the action for the system's Hamiltonian, The thermal equilibrium state can also be expressed in the path integral representation by making the replacement iτ / → τ ′ in Eq.(5). We thereby obtain where Z tot is the normalization constant, the G-numbers {φ † ,φ} form the coherent repre- The total density matrix elements is then given by The heat-bath degrees of freedom can be eliminated by integrating over the bath coordinates The reduced density operator is then expressed in the coherent representation of the G-numbers as [36][37][38] where Here, the influence functional for correlated initial states. 85 Employing the counter path illustrated in Fig. 1, we can express the influence functional in the path integral representation as (see where tively. The path integral used here to derive the HEOM is expressed in terms of an influence functional. The calculation of the influence functional for a heat bath consisting of harmonic The counter path for the influence functional given in Eq. (11).
oscillators is analogous to that of the generating functional for a Brownian oscillator system if we regard the system operator in the system-bath interactionV as an external force acting on the bath. [87][88][89][90] Then, the influence functional can be calculated analytically and is where the influence phase is expressed as (see Here, C ′ represents the counter path for s ′ that follows s ′′ along C under the condition s ′′ > s ′ and After dividing the contour of the integral in Eq. (14) into C 1 , C 2 and C 3 , we have (see where we define The functionals V × (t) and V • (t) represent the commutator and anticommutator ofV . This form of the influence functional has been used to analytically study quantum Brownian systems. 85 The first term in Eq.(16) represents a commonly used influence functional derived from the factorized initial conditions. 1-3 The collective bath oscillator coordinate,X ≡ j α jxj , is regarded as a driving force for the system through the interaction −VX. The timedependent kernels are then represented by iL 1 (t) ≡ [X(t),X] and L 2 (t) ≡ X (t)X + XX(t) /2, respectively, whereX(t) is the Heisenberg representation of the operatorX. 20 The function L 2 (t) is analogous to the classical correlation function of the bath induced noise X(t) and corresponds to the fluctuations. The dissipation corresponding toC(t) = − dtL 1 (t) is related to L 2 (t) through the quantum version of the fluctuation-dissipation theorem, L 2 [ω] = ω coth(β ω/2)C[ω]/2, which insures that the system evolves toward the thermal equilibrium state for finite temperatures. 4 The second term in Eq.(16) consists of the cross-terms between the real-time and imaginary-time integrals that describe the correlation between the initial equilibrium state and the dynamical state at time t. This term represents the contribution of the correlated initial conditions and can be regarded as the non-Markovian effects with respect to both real and imaginary times. The last term describes the influence of the heat bath on the thermal equilibrium state of the system. In the following sections, in order to derive the imaginary-time HOEM, we consider the full equilibrium state of the system, ρ eq [φ † ,φ; β ], by including the last term in ρ eq 0 [φ † ,φ; β ]. For 0 < τ < β , we have 85,88,89 where ν k ≡ 2kπ/ β are the Matsubara frequencies. We thus have andL For the case τ = 0, we use the definition given in Eq.(15) to obtain

III. REDUCED HIERARCHAL EQUATIONS OF MOTION IN REAL TIME
We assume that the spectral density J(ω) has an Ohmic form with a Lorentzian cutoff and write 20 where the constant γ represents the width of the spectral distribution of the collective bath modes and is the reciprocal of the correlation time of the noise induced by the bath.
The parameter η is the system-bath coupling strength, which represents the magnitude of damping.
With Eq.(24) for 0 < τ < β , we obtain where and c ′′ 0 = ηγ/β. At t = 0, the above equation reduces tō where and Here, we choose K so as to satisfy ν k = 2πK/(β ) ≫ ω c , where ω c represents the characteristic frequency of the system. Under this condition we can apply the approximation ν k e −ν k |t| ≃ δ(t) (for j ≥ K + 1) with negligible error at the desired temperature, 1/β. We define the equilibrium distribution function of the system under the influence of the heat bath through the replacement of the last term ofΦ[V; t, β ] (as expressed in Eq. (16)) We then obtain The influence functional F [V; t, β ] is redefined through this replacement as where and for k ≥ 1, and Note that in the high temperature limit, β γ ≪ 1, the noise correlation function reduces to where F (n) for nonnegative integers n, j 1 , . . . , j K . Among theρ 0,...,0 (t) =ρ(t) has a physical meaning, and the others are introduced for computational purposes. Differentiating ρ (n) j 1 ,...,j K (φ † , φ ′ ; t) with respect to t, we obtain the following hierarchy of equations in operator form: whereĤ × is the Liouvillian ofĤ A , and the relaxation operatorsΘ andΨ k are obtained through the replacement V × (t) →V × and V • (t) →V • in Eqs. (33) and (35), whereÔ ×f ≡ Of −fÔ andÔ •f ≡Ôf +fÔ for any operand operatorÔ andf , and The above expression is identical to the HEOM with a factorized initial state and can be truncated in the same manner as in the factorized case for large N ≡ n + Σ K k=1 j k ≫ ω c /min(γ, µ 1 ), where ω c is the characteristic frequency of the system. 20,38 If we add the counter term to the Hamiltonian (1), we have an additional term in Eq. (41). 53 While the terms from the correlated initial stateΘ andΨ k do not appear in Eq.(40), they define the hierarchy elements for the correlated initial equilibrium state. To demonstrate this point, we consider the initial states of the density operators, obtained by setting t = 0 in Eqs. (38) and (39): Here, and we have ρ eq [φ † ,φ; β ] = Z Bρ [φ † ,φ; β ]. This defines the correlated equilibrium initial conditions of Eq. (40). In the next section, we derive the equations of motion to evaluate these hierarchy elements.

IV. REDUCED HIERARCHAL EQUATIONS OF MOTION IN IMAGINARY TIME: CORRELATED THERMAL EQUILIBRIUM STATE
The thermal equilibrium stateρ[φ † ,φ; τ ] at time t = 0 and inverse temperature τ can be obtained by considering the imaginary-time derivative of the reduced density matrix elements given in Eq. (31). This is expressed as Note that the first product in Eq. To illustrate the structure of the hierarchy given in Eq. (44), here we write out the equations up to second order. The hierarchy starts from the zeroth-order equation, which is that for the thermal equilibrium state density matrix: Then, the first order consists of two equations, ∂ ∂τρ ∂ ∂τρ ρ (0) 0,...,j k =1,0,...,j k ′ =1,0,...,0 (0)= The elementsρ (0) 0,...,j k =2,0,...,0 (0) are obtained by setting k = k ′ in Eq. (56). In practice, in order to evaluate the HEOM elements in the case of correlated initial conditions from the imaginary-time HEOM, the cutoff, K ′ , must be comparable to the cutoff of used for the real-time HEOM given in Eq. (40), K. If we only need the equilibrium distribution, Eq.(52), however, we may choose K ′ even slightly smaller than K/2.
The equilibrium reduced density matrix has been evaluated from various approaches. 93,94 Equation (44) We chose the system parameters as ω 0 = 1 and ∆ = 0 or 1, and the bath parameters as β = 0.5 ∼ 5, η = 0 ∼ 2, and γ = 0.5 for the system-bath interactionV =σ x . We truncated the hierarchy by settingρ [K ′ +1:l] k 1 ,...,k K ′ +1 (τ ) = 0 for K ′ = 6 in the imaginary-time HEOM, while we truncated by settingρ Note that we must chose K ′ ≈ K in order to accurately calculate the real-time HEOM elements for the correlated initial conditions from the imaginary-time HEOM. Then, in order to obtain a better accuracy for deeper hierarchy elements in the imaginary-time HEOM, we used a small time step in the numerical integrations. For this reason, the computational costs for the real-time and imaginary-time HEOM were comparable. However, if we merely needed the equilibrium distributionρ [0:0] (β ) to one percent accuracy, we could use a smaller cuttoff K ′ and/or a larger time step for the imaginary-time HEOM and thereby reduce the computational costs to less than 1% of that for the results reported here.

B. Partition functions and thermodynamic variables
Although with the real-time HEOM, we can calculate only the probability distribution, with the imaginary-time HEOM we are able to calculate thermodynamic variables via the partition function of the reduced system, Z A = tr A {ρ [0:0] (β )}. Note that the total partition function can be expressed as Z tot = Z A Z B , where the partition function of the bath is given by Because we consider an infinite number of oscillators, however, the partition function of the bath cannot be determined. For this reason, we consider the system part, Z A , only. We [(β ω 0 /2) tanh(β ω 0 /2) -ln(2 cosh(β ω 0 /2))], and C 0 A = k B (β ω 0 ) 2 /4 cosh 2 (β ω 0 /2). Also, note that the susceptibility for finite ∆ is expressed as χ 0 A = ∆ tanh β ω 2 0 + ∆ 2 )/2 /2. The superscript "0" on these quantities indicates that these are calculated using the conventional statistical physics approach, which is equivalent to assuming a factorized thermal equilibrium state.
As seen in Fig. 2, in both cases of the spin-boson and factorized spin system, the entropy and internal energy decrease with the inverse temperature, while the heat capacities of both systems exhibit maxima at inverse temperatures near β = 2, where the thermal excitation energy becomes comparable to the excitation energy. The entropy in the spin-boson case is larger than that in the factorized case at lower temperatures because the spin-boson system involves more degrees of freedom, due to the presence of the system-bath interaction. It is also seen that the internal energy is systematically lower in the spin-boson case than in the factorized case. This indicates that the bath absorbs some of the system energy through the interaction. The degree to which the system energy is absorbed by the bath increases as β approaches the thermal excitation energy of the system and, as a result, the heat capacity of the spin-boson system becomes smaller than that of the factorized spin system near the peak position at β = 2. Compared with the other thermodynamic variables, the difference between the susceptibilities in the two cases is small. This is because the systembath interaction has the same form as the magnetic excitation, and the effects of ∆ are suppressed by the strong system-bath interaction.
It is important to note here that those states regarded as the thermal equilibrium states in the two cases compared above are different. In the conventional treatment, the thermal equilibrium state of the system corresponds to the case of a factorized partition function, while in the present treatment of the spin-boson system, we consider the thermal equilibrium state of the total system. Although the difference between the equilibrium thermodynamic quantities for the spin-boson system and the factorized spin system are rather minor in the static case considered in Fig. 2, the difference becomes significant when we study the dynamics of the system, because in this case, the positivity condition is often violated in the conventional treatment. This may indicate that treatments based on the canonical distribution are inherently incompatible with dynamical states.

C. Auxiliary hierarchy elements and expectation values
By utilizing the hierarchy elements, we can calculate expectation values of the system and bath. For example, the expectation value of the system energy, Ĥ A , is obtained from Eq. (52) as Using the first element of the hierarchy, the expectation value of the system-bath interaction, Ĥ I = tr{V c jxj exp[−βĤ tot ]}, is evaluated as   but the bath part of the contribution is much larger than the system part, because the bath contains many degrees of freedom.

VI. CONCLUDING REMARKS
In this paper, we derived the real-time and imaginary-time HEOM starting from the influence functional formalism with a correlated thermal initial state. It was shown that the thermal equilibrium state calculated from the imaginary-time HEOM is equivalent to the steady state solution of the real-time HEOM. Because the imaginary-time HEOM is defined in terms of integrals carried out over the definite time interval from τ ′ = 0 to τ ′ = β and because the elements of the imaginary-time HEOM are real, we were able to calculate the hierarchy elements more easily in this case than in the case of the real-time HEOM. Moreover, using the imaginary-time HEOM, we were able to calculate the partition function, and from this, we could directly obtain several thermodynamic quantities, namely, the free energy, entropy, internal energy, heat capacity, and susceptibility of the system in the dissipative environment. The expectation values of not only the system energy but also the systembath interaction energy were also evaluated from the hierarchy elements obtained from the real-and imaginary-time HEOM. We found that for the purpose of studying equilibrium properties, rather than dynamical behavior, the imaginary-time HEOM is more usueful than the real-time HEOM.
For the bath Hamiltonian appearing in Eq.(2) with the interaction −V α j x j , Eq.(A19) is expressed as Eq. (15). By tracing out p and r, we obtain the influence functional for correlated initial conditions, F [Ṽ C ; t, β ] = exp Φ [Ṽ C ; t, β ] .