Real-Time and Imaginary-Time Quantum Hierarchal Fokker-Planck Equations

We consider a quantum mechanical system represented in phase space (referred to hereafter as"Wigner space"), coupled to a harmonic oscillator bath. We derive quantum hierarchal Fokker-Planck (QHFP) equations not only in real time, but also in imaginary time, which represents an inverse temperature. This is an extension of a previous work, in which we studied a spin-boson system, to a Brownian system. It is shown that the QHFP in real time obtained from a correlated thermal equilibrium state of the total system possess the same form as those obtained from a factorized initial state. A modified terminator for the hierarchal equations of motion is introduced to treat the non-Markovian case more efficiently. Using the imaginary-time QHFP, numerous thermodynamic quantities, including the free energy, entropy, internal energy, heat capacity, and susceptibility can be evaluated for any potential. These equations allow us to treat non-Markovian, non-perturbative system-bath interactions at finite temperature. Through numerical integration of the real-time QHFP for a harmonic system, we obtain the equilibrium distributions, the auto-correlation function, and the first- and second-order response functions. These results are compared with analytically exact results for the same quantities. This provides a critical test of the formalism for a non-factorized thermal state, and elucidates the roles of fluctuation, dissipation, non-Markovian effects, and system-bath coherence. Employing numerical solutions of the imaginary-time QHFP, we demonstrate the capability of this method to obtain thermodynamic quantities for any potential surface. It is shown that both types of QHFP equations can produce numerical results of any desired accuracy.

previous work, in which we studied a spin-boson system, to a Brownian system. It is shown that the QHFP in real time obtained from a correlated thermal equilibrium state of the total system possess the same form as those obtained from a factorized initial state. A modified terminator for the hierarchal equations of motion is introduced to treat the non-Markovian case more efficiently. Using the imaginary-time QHFP, numerous thermodynamic quantities, including the free energy, entropy, internal energy, heat capacity, and susceptibility can be evaluated for any potential. These equations allow us to treat non-Markovian, non-perturbative system-bath interactions at finite temperature. Through numerical integration of the real-time QHFP for a harmonic system, we obtain the equilibrium distributions, the auto-correlation function, and the first-and second-order response functions. These results are compared with analytically exact results for the same quantities. This provides a critical test of the formalism for a non-factorized thermal state, and elucidates the roles of fluctuation, dissipation, non-Markovian effects, and system-bath coherence. Employing numerical solutions of the imaginary-time QHFP, we demonstrate the capability of this method to obtain thermodynamic quantities for any potential surface. It is shown that both types of QHFP equations can produce numerical results of any desired accuracy. The FORTRAN source codes that we developed, which allow for the treatment of Wigner space dynamics with any potential form, (TanimuranFP15 and ImTanimuranFP15) are provided as supplementary materials. a) Electronic mail: tanimura@kuchem.kyoto-u.ac.jp.

I. INTRODUCTION
A Brownian oscillator (BO) model, which consists of a primary system coupled to a harmonic oscillator bath, is a versatile model that has been used to investigate fundamental problems in physics, chemistry and biology. [1][2][3][4][5][6][7][8] The key feature of the Brownian model is that it describes irreversible dynamics through which the system evolves toward the thermal equilibrium state at finite temperature. This feature arises from interaction with the heat bath, which exhibits the canonical distribution at temperature T . To make the heat bath an unlimited heat source that possesses infinite heat capacity, the number of heat bath oscillators is effectively made infinitely large by replacing the spectral distribution of the system-oscillator coupling, J(ω), which was originally defined as the discretized distribution J(ω) = c 2 j δ(ω − ω j ) (where c j is the coupling strength between the system and the jth bath oscillator with frequency ω j ), with a continuous distribution, for example, J(ω) ∝ ω.
Because the time-evolution of the total system is described by the Schrödinger equation, the total energy is conserved and the dynamics are reversible. In the reduced description of the system obtained by tracing over the bath degrees of freedom using such methods as the path integral method 1 or the projection operator method, 4,7 however, the energy is no longer conserved, and its dynamics are irreversible, because the reduced system is merely a part of the total system. Heat bath effects arise in the reduced dynamics as fluctuation and dissipation in the reduced main system. These satisfy the classical or quantum version of the fluctuation-dissipation theorem. The reduced system evolves in an irreversible manner toward the thermal equilibrium state, in which the energy supplied by the fluctuations and the energy lost through dissipation are balanced, while the bath temperature does not change, because its heat capacity is infinite.
With the above described features, the Brownian model exhibits wide applicability, despite its simplicity. This is because the influence of the environment can in many cases be approximated by a Gaussian process, due to the cumulative effect of the large number of weak environmental interactions, in which case the ordinary central limit theorem is applicable, 9 while the distribution function of the harmonic oscillator bath itself also exhibits a Gaussian distribution. By adjusting the form of the spectral distribution, the properties of the bath can be adjusted to represent a variety of environments consisting of, for example, solid state materials, solvates, and protein molecules. This model has been used to solve various problems of practical interest, in particular to investigate tunneling processes, 2,3,10 chemical reaction, 11,12 non-adiabatic transition, 13,14 quantum device systems, 15 ratchet rectification, 16,17 to evaluate the efficiency of SQUID rings, 18,19 and to analyze the line shapes in laser spectra. 20,21 While the Brownian model itself is fairly simple, it is somewhat difficult to apply in the quantum mechanical case not only analytically but also numerically, due to the infinite number of bath degrees of freedom. Analytically exact solutions of Green's function for the BO Hamiltonian have been obtained only in the cases of a harmonic oscillator, 5,6 a free particle, 22 and a free rotator, 23 using the path integral approach. Several approximate approaches have been developed to facilitate application of the Brownian model to more complicated systems. These approaches involve variational methods to study polarons 24,25 and the optical response of an anharmonic oscillator 26 using a damped oscillator as a trial function, an instanton method for estimating the tunneling rate using instantaneously jumping paths between tunneling wells, 2,3 a WKB method for evaluating the density matrix along a classical minimal action path, 10,11 and diagrammatic expansion methods to study the anharmonicity of potentials and the nonlinearity of the system-bath coupling. [27][28][29] The analytical expressions obtained in these studies are helpful to gain insight into the role of dissipative environments in the dynamics of systems, but they do not allow us to study situations investigated in modern experiments that are usually described by complex potentials driven with time-dependent external forces.
A great deal of effort has been dedicated to numerically calculating the time evolution of BO systems under external perturbations. Widely used approaches employ a reduced equation of motion that can be derived from the quantum Liouville equation with the full Hamiltonian by reducing the heat bath degrees of freedom. To obtain reduced equations of motion in a compact form, one usually employs the Markovian assumption, in which the correlation time is very short in comparison to the characteristic time of the system dynamics. In this case, the noise can be regarded as white. The quantum Langevin equation and the quantum Fokker-Planck equation have been derived with the projection operator method and the path integral method, for example. [30][31][32][33][34][35] In the classical case, the Langevin equation 36 and the Fokker-Planck (or Kramers) equation 37,38 have proved to be useful in the treatment of transport problems, and they have even been included in algorisms employed in molecular dynamics simulations. However, the applicability of the quantum forms of these equations is very limited, because they cannot be derived in a quantum mechanical framework without severe approximations and/or assumptions. For example, in the treatment of the quantum Langevin equation expressed in operator form, it is generally assumed that the antisymmetric correlation function of the noise is very short (the Markovian assumption) and positive. A similar Markovian assumption has been used in the treatment of the quantum Fokker-Planck equation. But in order for these assumptions to be valid, the heat bath must be at a sufficiently high temperature, in which case most of the important quantum dynamical effects play a minor role. This implys that the Markovian assumption is incompatible with obtaining a quantum mechanical description of dissipative dynamics at low temperature. 39 An Ohmic spectral distribution is generally assumed to realize Markovian noise. As we show in Appendix B, however, even if the dissipation process is Markovian, the fluctuation process may not be, because it must satisfy the fluctuation-dissipation theorem. 8 For this reason, if we apply the equation of motion under Markovian assumption to low temperature systems, then the positivity of the probability distributions of the reduced system cannot be maintained. 40 This is a fundamental limitation, known as the "positivity problem," which is particularly significant for the quantum master equation. [41][42][43][44][45][46][47][48] If the system is not time dependent and if the system Hamiltonian and the system-bath interaction commute, the time-convolutionless (TCL) master equation becomes exact. [49][50][51] For the time-dependent case and/or non-commuting case, however, this master equation is valid only to second order with respect to the system-bath interaction, and the positivity condition is again broken. 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 in order to put the equation of motion in the Lindblad form. However, the RWA alters the thermal equilibrium state and the dynamics of the reduced system. These changes are particularly large in the case of a strong system-bath coupling and at low temperature. Moreover, in a typical quantum transport problem, the system is described by continuous energy states, and the energy levels of the heat bath and the system overlap. For this reason, the RWA cannot be used. Treatments of these kinds are therefore not sufficient to construct fully quantum mechanical descriptions of broad validity.
Path integral Monte Carlo simulations do not have any of the limitations of the approaches discussed above, but this approach is computationally intensive, because the number of paths to be evaluated grows rapidly with time, while sampling fails, due to the phase cancellation of wave functions. [52][53][54] Much effort has been made to overcome these problems and extend the applicability of this method. [55][56][57][58][59][60][61][62][63][64] Because this approach can easily incorporate the semiclassical approximation for the bath, it may be advantageous in the study of polyatomic systems treated in multi-dimensional coordinates, but applications to this point incorporating full quantum mechanical dynamics have been limited to relatively small systems without time-dependent external force.
Wave function based methodologies for the full Hamiltonian have been developed in order to avoid the reduced description of the system. The multi-configurational time-dependent Hartree (MCTDH) approach 65-72 employs time-dependent basis sets to represent the total wave function. Then, a variational principle is applied to derive the optimal equation of motion in order to reduce the bath degrees of freedom. This approach can be used to treat nonlinear system-bath coupling and anharmonic bath modes. 68 However, the number of bath modes must be increased until convergence is reached. This implies that the study of long time behavior requires more basis sets, which makes the calculation more difficult. In the effective-mode approach, the heat bath degrees of freedom are mapped to a linearly coupled harmonic oscillator chain. Then, the dynamics of the system are described by the wave function of the system with a finite number of chained oscillators using a truncation scheme [73][74][75] or by utilizing the density matrix renormalization group method. 76 Strictly speaking, the time evolution obtained with the wave function based approach describes time-reversible processes and thus, within this approach, there exists no thermal equilibrium state. However, in practice, this kind of approach has wider applicability than the reduced equation of motion. At this stage, the results obtained from these approaches have been limited to relatively simple systems. In particular, the inclusion of time dependent external forces is not as straightforward in these approaches as in the case of reduced equation of motion, because the energy of the total system changes due to the presence of an external force if the perturbation is strong, and hence the optimal basis set may also be changed.
The reduced hierarchal equations of motion (HEOM), which are derived by differentiating the reduced density matrix elements defined by path integral, are reduced equations of motion that can describe the dynamics of the system for non-perturbative and non-Markovian system-bath interactions with any desired accuracy under strong time-dependent perturbations at finite temperature. 8 In this formalism, the effects of higher-order non-Markovian system-bath interactions are mapped into the hierarchal elements of the reduced density matrix. In their original formulation, these equations of motion were limited to the case in which the spectral distribution function takes the Drude form (i.e., the Ohmic form with a Lorentzian cut-off) and the bath temperature is high. 77 However, with the inclusion of low temperature corrections terms, this temperature limitation has been eliminated. [78][79][80][81] In addition, with the extension of the dimension of the hierarchy, in its present form, this approach is capable of treating a great variety of spectral distribution functions. [82][83][84][85][86][87][88][89][90] This formalism is valuable because it can be used to treat not only strong system-bath coupling but also quantum coherence between the system and bath, which is essential to study a system subject to a time-dependent external force 8 and nonlinear response functions. [91][92][93] The system-bath coherence becomes particularly important if the bath interaction is regarded as non-Markovian, as found from femtosecond nonlinear optical measurements, which are carried out on time scales that are much shorter than the noise correlation time of environmental molecules. 20 For a Brownian system, the reduced hierarchal equations of motion are expressed in the Wigner space representation. [94][95][96][97][98][99][100][101][102][103][104][105][106][107] In the Markovian limit, these equations of motion reduce to the Caldeira-Leggett quantum Fokker-Planck equation, 30,31 and in the classical limit, they reduce to the classical Fokker-Planck (Kramers) equation. 37,38 Recently, the author derived the HEOM not only in real time, but also in imaginary time, which represents an inverse temperature, starting from correlated initial conditions for a system described by discretized energy states. 108 Reduction of these HEOM to a system represented in Wigner space is not straightforward, because they involve derivatives with respect to the position and momentum that require a careful treatment with regards to the order of time slices in the path integral formalism. In this paper, we present the derivation of real-and imaginary-time HEOM in Wigner space and demonstrate the validity of these equations.
The organization of the paper is as follows. In Sec. II, we present a model Hamiltonian and its influence functional with correlated initial conditions. In Sec. III, we derive the realtime quantum hierarchal Fokker-Planck (the real-time QHFP) equations using the influence functional given in Sec. II. In Sec. IV, we derive the imaginary-time quantum hierarchal Fokker-Planck (the imaginary-time QHFP) equations, which are convenient for evaluating thermodynamic quantities of the system. In Sec. V, the validity of our approach is demon-strated through numerical integration of the real-and imaginary-time QHFP equations for a harmonic system and comparing the calculated results with the exact results obtained from analytical calculations. Section VI is devoted to concluding remarks.

II. REDUCED HIERACHAL EQUATIONS OF MOTION FROM CORRELATED INITIAL CONDITIONS
We consider the 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 Brownian Hamiltonian expressed as 1-8 whereĤ is the Hamiltonian for the system with mass m and potential U(q) described by the momentump and positionq. The bath degrees of freedom are treated as an ensemble of harmonic oscillators, and the momentum, position, mass, and frequency of the jth bath oscillator are given byp j ,x j , m j and ω j , respectively. In the conventional Brownian model, the systembath interaction is represented by a bilinear function of the system and bath coordinates as H I = −q j α jxj . Brownian models employing this bilinear interaction have been studied with various approaches. [4][5][6][7] In this paper also we restrict our investigation to this bilinear form to simplify the derivation, but we note that extension to the non-bilinear case is possible. 8,[98][99][100][101][102] To maintain translational symmetry in the case of U(q) = 0, required to describe the motion of a free Brownian particle, we include the counter-term j α 2 The heat bath we consider is characterized by the spectral distribution function defined by J(ω) ≡ j ( α 2 j /2m j ω j )δ(ω − ω j ) and the inverse temperature, β ≡ 1/k B T , where k B is the Boltzmann constant. The path integral used here to derive the HEOM is expressed in terms of an influence functional with correlated initial conditions. The influence functional that we employ, F CI [t, β ], is calculated by taking the trace over the heat bath degrees of freedom, starting from the thermal equilibrium state of the total Hamiltonian. The calculation of the influence functional for a heat bath consisting of harmonic oscillators is analogous to that of the generating functional for a Brownian oscillator system if we regard the system operator in the system-bath interactionq as an external force acting on the bath. [109][110][111] As shown in Appendix A, the reduced density matrix elements of the system with correlated initial conditions can be expressed as where S A [q; t] is the action for the Hamiltonian of the system, Eq. (2), given by andρ eq [q; β ] is the initial thermal equilibrium state, with the heat bath defined by Eq.(A12).
We assume that the spectral density J(ω) has an Ohmic form with a Lorentzian cut-off and write 8 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. This spectral distribution approaches the Ohmic distribution, J(ω) ≈ mζω/π, for large γ. In Appendix B, we present several profiles of fluctuation and dissipation terms for the Drude distribution to illustrate the origin of the positivity problem in the Markovian master and Redfield equations.
With J(ω) given by Eq. (5), the influence functional with correlated initial conditions is expressed as 108 where, for the Matsubara frequency ν k ≡ 2πk/β , we have defined and for k ≥ 1, and where C k ≡ ν 2 k /(ν 2 k + ω 2 c ) is the correction factor that counteracts the overestimation of the contribution of higher-order Matsubara frequencies approximated by the delta function with cut-off number, K, introduced in Appendix B for the characteristic frequency of the system, ω c . This modification improves the convergence of hierarchies at lower temperature. We now introduce the hierarchal elements that play an essential role in our formalism: where F (n) for nonnegative integers n, j 1 , . . . , j K . From the above definition, the first hierarchal element and the reduced density matrix given by Eq.(3) are identical: ρ(q, q ′ ; t) = ρ (0) 0,...,0 (q, q ′ ; t). As shown in Appendix C, we then have the following equations of motion: where and Φ(q, q ′ ), Θ k (q, q ′ ), and Ξ ′ (q, q ′ ) are defined by Eqs. (7), (11), and (13) by making the replacements q(t) → q and q ′ (t) → q ′ . In the HEOM formalism, only the first element has a physical meaning and the other elements ρ (n) j 1 ,...,j K (q, q ′ ; t) are introduced in numerical calculations in order to treat the non-perturbative and non-Markovian system-bath interaction. We can evaluate ρ (0) 0,...,0 (q, q ′ ; t) through numerical integration of the above equations.
We next explain the truncation scheme that we use for the hierarchical equations, which is different from the scheme used in previous studies. [102][103][104][105] First, we choose the number of Matsubara frequencies to be included in the HEOM, K, such that it satisfies K ≫ ω 0 /ν 1 .
Then, we introduce the scaled integer K γ as K γ ≡ int(Kν 1 /γ) for ν 1 > γ and K γ ≡ K for ν 1 ≤ γ, which allows us to make calculations in the highly non-Markovian case more efficiently. The index for the hierarchy, denoted by n, for a given value of γ, then runs from 0 to K γ . The total number of hierarchy members to be included in the calculations is then given by N ≡ (K γ + K + 1)!/(K + 1)!/K γ !. For the case K k=1 j k > K, we truncate the hierarchal equations by replacing Eq.(16) with In practice, we can simply set ρ (n) j 1 ,...,j K (q, q ′ ; t) = 0 instead of employing the above equation, because ρ (n) j 1 ,...,j K (q, q ′ ; t) decays to zero as t becomes large. 79 For the K γ and K γ + 1 members of the hierarchy, we have the following relation, valid to order δt: This asymptotic relation allows us to obtain the terminator for given γ in the form This equation reduces to the quantum Fokker-Planck equation in the Markovian limit, i.e., the Ohmic distribution (γ → ∞) with the high-temperature limit. 30,94 While the termsΘ andΨ k from the correlated initial state do not appear in Eqs. (16) and (21), they define the hierarchal elements for the correlated initial equilibrium state. 108 To demonstrate this point, we consider the initial states of the density operators, obtained by setting t = 0 in Eq. (14): whereρ (m) are the equilibrium hierarchal elements. Here, Z A , Z tot , and Z B are the partition functions of the system, total system, and bath, respectively, related as It is important to note that the steady state of ρ n j 1 ,...,j K (t) for n > 0 in Eq.(16) is slightly shifted from the initial thermal equilibrium state as a result of the influence of the sum in Eq. (22). However, because ρ  Fig. 9, the correlations responsible for the correlated initial state are represented by green arcs, and the static correlations are represented by red arcs. After the time evolution, the elements ρ (n) j 1 ,...,j K (q, q ′ ; t) describe the dynamical correlation, represented by the blue arcs and lines in Fig. 9.

III. REAL-TIME QUANTUM HIERARCHAL FOKKER-PLANCK EQUATIONS
We now introduce the Wigner distribution function, which is the quantum analog of the classical distribution function in phase space. For the density matrix element ρ (n) The Wigner representation of the reduced density matrix defined in Eq. (3), W (p, q; t), and the first member of the hierarchal elements are then identical: W (p, q; t) = W (0) 0,...,0 (p, q; t). The Wigner distribution function is a real function, in contrast to the complex density matrix. In terms of the Wigner distribution, the quantum Liouvillian takes the form 115 where U W (p, q) is given by The quantum Liouvillian can also be expressed as 113,114 While the above expression is easier to integrate in the case that the potential is nearly harmonic, the expression in Eq.(25) is numerically stable, and it can be applied with any form of potential, including an unbounded potential.
Using the Wigner distribution and quantum Liouvillian, the equations of motion appearing in Eq. (16) can be expressed in the form of quantum hierarchal Fokker-Planck (QHFP) equations in real time as andΞ As in the case of the energy eigenstate representation, 108 the above equations are identical to the equations derived from factorized initial conditions. [102][103][104][105] The above equations are then truncated by using the modified "terminators" expressed in the Wigner representation. As explained in Sec. II, the number of Matsubara frequencies to be included in the calculation, K, is chosen to satisfy K ≫ ω c /ν 1 . The upper limit for the number of hierarchy members for given γ is then chosen to be K γ ≡ int(Kν 1 /γ) for ν 1 > γ and K γ ≡ K for ν 1 ≤ γ. Then, for the case K k=1 j k > K, we truncate the hierarchal equations by replacing Eq.
while, for the case n = K γ we employ We can evaluate W (n) j 1 ,...,j K (p, q; t) through numerical integration of the above equations. While only the first element W (p, q; t) ≡ W (0) 0,0,··· ,0 (p, q; t) has a physical meaning and the other elements W (n) j 1 ,...,j K (p, q; t) are initially introduced to avoid the explicit treatment of the inherent memory effects, it turns out, however, that these elements allow us to take into account the system-bath coherence, 8 entanglement 88,116,117 and expectation values that include the bath operators as Ĥ I ≡ − q α jxj . 108 The HEOM consist of an infinite number of equations, but they can be evaluated with the desired accuracy by depicting the asymptotic behavior of the hierarchal elements for different K and using this to determine whether or not there are sufficiently many members in the hierarchy. Essentially, the error introduced by the truncation to be negligibly small when K is sufficiently large.
The correlated initial equilibrium state defined by Eq. (22) is expressed in the Wigner representation accordingly. The correlated initial equilibrium state can be set in the HEOM formalism by running the HEOM program until all of the hierarchy elements reach the steady state and then use these elements as the initial state, 8 or by integrating the imaginary-time HEOM that we discuss in the next section. 108 In practice, the former approach is simpler, because it requires the real-time HEOM only. This approach has been used to set the correlated initial conditions of the HEOM derived from factorized initial conditions that are identical to those used with the present HEOM.
The HEOM in Wigner space is ideal for studying quantum transport systems, because it allows the treatment of continuous systems, utilizing open boundary conditions and periodic boundary conditions. 103,104 In addition, the formalism can accommodate the inclusion of an arbitrary time-dependent external field. [95][96][97]105 In the Markovian limit, γ → ∞, which is taken after the high temperature limit, yielding the condition β γ ≪ 1, we have the quantum Fokker-Planck equation 3,10 which is identical to the quantum master equation without the RWA. 8 Because we assume that the relation β γ ≪ 1 is maintained while taking the limit γ → ∞, this equation cannot be applied to low-temperature systems, in which quantum effects play a major role.
As in the case of the master equation without the RWA, the positivity of the population distribution, P (q) = dpW (p, q; t), cannot be maintained if we apply this equation in the low temperature case.
The classical HEOM can be derived by taking → 0. 94,95 The Wigner distribution function reduces to the classical one in this limit. The classical equation of motion is helpful, because knowing the classical limit allows us to identify the purely quantum mechanical effects. 95,102,105 IV. IMAGINARY-TIME QUANTUM HIERARCHAL FOKKER-PLANCK

EQUATIONS
The equilibrium reduced density matrix has been evaluated with several approaches. 118,119 By applying the methodology developed in Ref. 108 we obtain the imaginary-time QHFP equations as ∂ ∂τW where the factorsc k are expressed asc 0 = mζγ/β andc k = 2mζγ 2 /β(γ + ν k ), for k ≥ 1. We setW [m:l] k 1 ,...,k m (τ ) = 0 for higher-order elements in hierarchy denoted by m to truncate. The Euclidean Liouvillian is expressed as withŪ for the potential, U ′ (q) = U(q) + mζγq 2 /2, including the counter-term. This can also be expressed in differential form as If the anharmonicity of the potential is small, the above expression is useful. The initial By integrating Eq. (36) from τ = 0 to τ = β , we can evaluate the equilibrium distribution functionW (p, q; β ).
Once we obtain the equilibrium distribution, we can calculate the partition function employing the relation This allows us to calculate the Helmholtz free energy, F A = − ln(Z A )/β, the entropy, /∂β, and the heat capacity, for any potential. If the system is subject to an external force ∆f (p,q), where f (p,q) is any function of the momentum and position,p andq, we can also calculate the susceptibility, It should be noted that even if the potential is a function of time, we can calculate thermodynamic quantities as functions of time through Z A (β ; t), assuming that the system reaches the thermal equilibrium state faster than the change of the potential.

V. NUMERICAL RESULTS
In principle, the HEOM provide an asymptotic approach that allows us to calculate various physical quantities with any desired accuracy by adjusting the number of hierarchal elements. Here, we demonstrate the applicability and validity of the real-time and imaginarytime QHFP equations, by presenting the results obtained from numerical integrations of Eqs.(28)- (33) and Eqs. (36)- (38). For this purpose, we consider the harmonic potential However, here we consider the non-commuting case.
Below we also present our results for calculations of thermodynamic quantities obtained from the imaginary-time QHFP and compared them with analytical results.
A. Steady state distribution: Static system-bath coherence and mixed state For a harmonic system, the equilibrium distribution in the Wigner representation is analytically expressed as whereN ≡ 2π p 2 q 2 is the normalization factor and q 2 and p 2 are the mean squares of the position and momentum, respectively.
The Wigner distribution for an isolated oscillator is written W eq A (p, q). For the Hamiltonian Eq.(41), we have 113 and  (43) and (44) is set as the temporally initial state at time t = 0. After integrating the real-time QHFP and the TCL Redfield equations for a sufficiently long time (t = 100), the distribution reaches the steady state. In the real-time QHFP case, the obtained steady state is identical within numerical error to the thermal equilibrium state W eq BO (p, q) with q 2 and p 2 given by Eqs. (45) and (46), while those from the TCL Redfield equations are similar to the original factorized initial state.
This implies that the TCL Redfield equation cannot take into account the system-bath correlation properly.
The Wigner distribution for a harmonic Brownian system is denoted by W eq BO (p, q). In this case, we have 5-7 and with δΓ 2 (ω) ≡ ζγ 2 ω/(γ 2 + ω 2 ). In the Wigner representation, the thermal equilibrium state under the factorized assumption, exp(−βĤ A ) exp(−βĤ B ), is denoted by W eq A (p, q), while the true thermal equilibrium state of the reduced density operator, is denoted by W eq BO (p, q); the difference between the two distributions arises from the static system-bath coherence and represents the non-factorized effect of the thermal equilibrium state.
To obtain the thermal equilibrium state from the real-time QHFP, we integrated Eqs. In Fig.1(a), we display W (0) 0,...,0 (p, q; t) for the factorized initial state at t = 0 given by W eq A (p, q) (blue curves) and the steady state distribution at t = 100 (red curves) obtained from the real-time QHFP calculation. We found that even if we start from the factorized initial state, the steady state solution is the true thermal equilibrium state, denoted by W eq BO (p, q). This indicates that the real-time QHFP has the capability to produce the thermal equilibrium state with a static system-bath correlation through the fluctuation and dissipation terms. In the TCL Redfield equation cases, Figs. 1(b) and (c), the calculated steady states (red curves) are similar to the factorized initial states (blue curves). The peak intensity of the TCL result in the case without the RWA is slightly higher than that in the case with the RWA, because the ground and first excited populations in the former case are ρ 00 (t) = 1.083 and ρ 11 (t) = −0.090, due to the breakdown of the positivity condition, which is a physical requirement for the reduced equations of motion necessary for the pop-ulation state of the density matrix to be positive. [41][42][43][44][45][46][47][48] Other than this difference, the TCL distributions are similar to the factorized distribution. This indicates that the TCL Redfield equation cannot take into account the static system-bath coherence, because the TCL theory in the present case is valid only to second order in the system-bath coupling.
As explained in Appendix A, the effects of the system-bath coherence consist of the imaginary-time (static) part and complex-time (correlated) part represented by the red and green arcs in Fig. 9, respectively. From the equilibrium distribution, we can only observe the effects of the static part. To elucidate the correlated part, we need to calculate the nonlinear response function, as will be discussed in Sec. V-C.
In order to calculate the above functions using an equation of motion approach, we employ the following forms: 8,120 and whereq ×Â ≡qÂ −Âq,q •Â ≡qÂ +Âq,Ĝ(t)Â ≡ e −iĤtott/ Â e iĤtott/ for any operatorÂ, and ρ eq tot = e −βĤtot /Z tot with Z tot = tr{ρ eq tot }. In the reduced equation of motion approach, the density matrix is replaced by a reduced one. In the QHFP case,ρ eq tot is replaced by the hierarchy member W

Auto-correlation function: Fluctuation and temperature effects
First we study the temperature dependence of the auto-correlation function for the fixed coupling strength ζ = 1 and the inverse noise correlation time γ = 1. In Fig. 2  rium state, and the discrepancy between the TCL results and exact results arises from the equilibration process discussed in Sec. V-A.

Linear response function: Dissipation and non-perturbative effects
As can be seen from Eq.(50), R (1) [ω] is temperature independent. Therefore, this function is convenient to study the non-perturbative effects of the system-bath coupling, ζ, and non-Markovian effects for slow modulation, controlled by the parameter γ, apart from the temperature effects. In Fig. 3 Fig. 3(a). For the strong coupling case considered in Fig. 3 (c), both the QHFP and analytical results exhibit a peak near ω 0 = 0.2. This peak arises from the strong coupling between the harmonic mode and the low frequency bath mode characterized by γ 2 ω/(γ 2 + ω 2 ) and only appears in the simultaneous non-Markovian (γ ≤ ω 0 ) and non-perturbative (ζ ≫ ω 0 ) case. 121 The existence of this peak, which we call a "non-Markovian bosonic peak," is a good indication of the applicability of non-perturbative and non-Markovian theories.
Because the TCL Redfield theory is valid only to second order in the system-bath coupling, the TCL results cannot reproduce this peak. Moreover, the spectrum calculated from the TCL Redfield equation without the RWA in the strong coupling case, shown in Fig.   3(c), is not positive for ω ≈ 5, due to the breakdown of the positivity condition. Despite this problem, however, the difference between the TCL results with and without the RWA is minor. This is because the spurious behavior caused by the positivity problem is suppressed in the non-Markovian treatment of the reduced dynamics, as explained in Appendix B.

Linear response function: Noise correlation and non-Markovian effects
We next discuss the non-Markovian effects in the Brownian oscillator system. It should be noted that when the inverse noise correlation time, γ, is decreased, the effective coupling strength becomes stronger, even if we fix ζ, because the bath can interact with the system multiple times when the correlation time is long. In order to study the pure non-Markovian effects, here we employ an effective coupling strength ζ ef f ≈ δΓ 2 (ω 0 ) = ζγ 2 ω 0 /(γ 2 + ω 2 0 ) 95 and fix it while varying γ.
In Fig. 4, we plot R (1) [ω] in the weak coupling regime corresponding to Fig. 3(a). Hereafter, we do not consider the TCL Redfield equation with the RWA, because the difference between the TCL results with and without the RWA is minor. While all of the peak profiles are similar if we fix ζ ef f , the peak position shifts slightly in the high-frequency direction, because a change of γ results in a change of δΩ 2 (ω 0 ). As the exact results and the QHFP results in Fig. 4 indicate, there is no clear indication of non-Markovian dynamics in this weak coupling regime, once we have normalized the effective coupling strength.
While the TCL Redfield results are close to the exact results in the fast modulation (weak non-Markovian) case depicted in Fig. 4(a), they differ significantly in the slow modulation (strong non-Markovian) case considered in Fig. 4(c). This is because the perturbative Thus the TCL-Redfield result without the RWA becomes negative for ω > 4.
It is seen that while the QHFP results always coincide with the exact results, the discrepancy between the TCL Redfield and exact results is large in the slow modulation (strong non-Markovian) case, due to the non-perturbative nature of the interactions. Specifically, the lack of a non-Markovian bosonic peak becomes apparent even at this intermediate coupling

C. Nonlinear response function: Dynamical system-bath coherence
As explained in Appendix C, the system-bath interaction induces static effects arising in imaginary time and dynamic effects arising in real time and complex time. While the static effects can be obtained from the thermal equilibrium distribution, as illustrated in Sec.
V-A, we have to study the nonlinear response function to elucidate the dynamic effects. It should be noted that, in addition to their inability to treat strongly non-Markovian dynamics, the conventional reduced equation of motion approaches involving the TCL Redfield equation have a severe limitation in studying systems subject to time-dependent external forces because their description of the damping kernels is based on energy eigenstates. 93 The capability of an approach to treat external forces can also be examined by calculating nonlinear response functions, because nonlinear measurements can capture the effects of multiple interactions through time-dependent external forces. Here, we calculate the second-order nonlinear response function of the position given by This is an observable in two-dimensional THz-Raman spectroscopy system. 122,123 Note that, because of the Gaussian integral involved in the expectation value ( · · · = tr{· · · exp(−βĤ tot )}), the contribution from the lowest-order response, [[q(t 1 + t 2 ),q(t 1 )],q] / 2 , vanishes. 8,123 In the harmonic case, there is also a contribution from R T RT (t 1 , t 2 ) = − [[q(t 1 + t 2 ),q 2 (t 1 )],q] / 2 , which corresponds to an observable in 2D THz-Raman-THz spectroscopy system. We find that to explore the system-bath coherence, Eq.(53) is suitable, as we show below. This response function in the harmonic Brownian case can be calculated analytically as 124 where C(t) is obtained from the Fourier transform of Eq. (49). To apply the Liouville operator formalism, we rewrite Eq.(53) as Using the above expression, we calculated R T T R (t 1 , t 2 ) for various values of t 1 and t 2 by extending the method employing Eqs. (51) and (52). 8,120 The response functions evaluated from Eqs. (54) and (55) are then Fourier transformed to obtain two-dimensional spectra, In Fig. 6, we plot 2D spectra in the frequency domain obtained from (a) the analytically exact approach, (b) the QHFP approach, and (c) the TCL Redfield without the RWA approach under the same physical conditions as in Fig. 5 (b). We find that while the analytically exact and QHFP results exhibit peaks at (ω 1 , ω 2 ) = (0, 1) and (ω 1 , ω 2 ) = (1, 1), the TCL approach cannot reproduce them. As shown in a study of multi-dimensional spectroscopy, in order to have these peaks, the dynamical system-bath coherence subject to the second interaction at time t 1 must be maintained throughout the time evolution described by Eq. (55). 93 In the TCL case, however, the time evolution is described in terms of the reduced operator tr B {Ĝ(t)}, derived from the factorization assumption with While the exact dynamics maintain the coherence during the period of length t 1 + t 2 expressed by C(t 1 + t 2 ), the TCL approach cannot maintain this coherence. In contrast to the Redfield approach, because the HEOM approach can store this coherence in the hierarchal members, it is capable of treating a nonlinear response function.
Because many modern experiments utilize the nonlinear response of a system, which is measured by applying a variety of time-dependent external forces, the capability to calculate the nonlinear response function is important. The validity of the HEOM approach has been demonstrated for systems subject to time-dependent external forces. 82,[95][96][97]105 In addition to the HEOM approach, the path integral approach has also been shown to have this capability. 60

D. Thermal equilibrium state and thermodynamic quantities
We finally examine the imaginary QHFP equation by considering our results obtained through numerical integration of Eq.(36) from τ = 0 to β using the harmonic potential to compare W eq BO (p, q) presented in V-A and the partition function Z A . The number of Matsubara frequencies used in the imaginary-time QHFP is K = 4. The mesh size was optimized for the Euclidean Liouvillian, and we chose n p = 60-120 and n q = 120-240.
Because the distribution is spread relatively widely in the higher temperature case, we employed a coarser mesh in that case.  In Fig.7, we display solution of the imaginary-time QHFP, Eq.(36) with β = 1 for several values of τ . Because the damping kernels in the imaginary-time QHFP are defined by the Matsubara frequency at β , the solutions τ < β do not correspond to the equilibrium distribution at temperature τ . While the initial distribution is flat, the distribution approaches a Gaussian form due to the Euclidean and the damping operators. At τ = β , the solution coincides with the analytical solution given in Eq. (42) with Eqs. (45) and (46).
While the equilibrium distribution can also be obtained from the real-time QHFP, as shown in Sec. V-A, the thermodynamic quantities can only be calculated from the imaginary-time QHFP. We next demonstrate this point. In the BO case, the partition function can also be evaluated analytically in terms of the Matsubara frequencies as 5-7 .
We should note that, the normalization constant of the real-time QHFP isN = 2π p 2 q 2 , whereas that of the imaginary-time QHFP is Z A obtained from Eq. (40). Because Z BO involves a temperature dependent factor other thanN , we cannot calculate the partition function using the real-time HEOM approach.
To obtain thermodynamic quantities, we first repeated the integration of the imaginarytime QHFP from β = 0.025 to 3.05 with step size ∆β = 0.025 to derive Z A . Then, we calculated thermodynamic quantities through Z A . In Fig. 8, we compare the partition function given by Eq. (40) (brown curve) and that obtained from Eq. (56) (dotted curve).
As expected, the imaginary-time QHFP results coincide with the exact results. For the purpose of demonstration, we also plot the entropy, S A , the internal energy, U A , and the heat capacity, C A , calculated with the imaginary-time QHFP. The behavior in the high temperature regime is very different from that in the spin-Boson case, 108 because the BO model has an infinite number of excited states.

VI. CONCLUDING REMARKS
In this paper, we presented real-time and imaginary-time QHFP equations derived using the influence functional formalism with correlated initial conditions. While we found that the QHFP equations in real time possess the same form as those obtained from a factorized initial state, we introduced a modified terminator in order to facilitate the more efficient spectroscopy, can also be taken into account in the QHFP formalism. 8,[98][99][100][101][102] Extension to multi-potential surfaces is also possible. 96,97 Because this formalism treats the quantum and classical systems with any form of potential from the same point of view, it allows identification of purely quantum mechanical effects through comparison of classical and quantum results in the Wigner distribution. 95,102,105 We showed that the thermal equilibrium state obtained from the imaginary-time QHFP is equivalent to the steady state solution of the real-time QHFP. Because the imaginary-time QHFP is defined in terms of integrals carried out over the definite time interval, we were able to calculate the equilibrium state more easily in this case than in the case of the real-time QHFP. Moreover, using the imaginary-time QHFP, 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, and heat capacity of the system in the dissipative environment. Numerical integration of the real-time and imaginary-time QHFP equations is computationally intensive. Nevertheless, we were able to study the dynamics of onedimensional potential systems using personal computers. [102][103][104][105] Great effort has been made to reduce the computational intensiveness of algorithms used to implement the real-time HEOM approach. For example, the hierarchy has been optimized for numerical calculations, [125][126][127][128][129][130][131] and a graphic processing unit (GPU) 132 and parallel computers 133

ACKNOWLEDGMENTS
Financial support from a Grant-in-Aid for Scientific Research (A26248005) from the Japan Society for the Promotion of Science is acknowledged.

Appendix A: Influence functional with correlated initial conditions
The reduced density matrix elements of the system are obtained in path integral form as where Z B is the partition function of the bath and the Euclidean action is given bȳ The influence functional for the correlated initial state expressed in terms of the influence phase is given by 108 where q × (t) ≡ q(t) − q ′ (t) and q • (t) ≡ q(t) + q ′ (t). Using the spectral density, J(ω), we rewrite these functions for 0 < τ < β as In the case τ = 0, we have L(t) ≡ iL 1 (t) + L 2 (t) with and in the case t = 0, we haveL(τ ) ≡ L(iτ ) with where the quantities ν k ≡ 2πk/β are the Matsubara frequencies. For later convenience, we also introduce the canonical correlation and express the counter-term of the potential using B(0). Using the relations the influence functional can be rewritten as where we have included the bath part in the initial thermal state of the systemρ eq A [q; β ] as 9. The roles of the system-bath interactions illustrated schematically. The solid black lines represent the wave function of the system A along the complex counter-path (see Fig. 1 in Ref. 108), and the arcs and blue lines represent the system-bath interactions. Because the bath degrees of freedom have been reduced, the bath interactions connect the wave function of the system A at multiple complex times. The blue arcs and lines correspond to the fluctuation and dissipation processes described by the terms containing B(t ′′ − t ′ ) and L 2 (t ′′ − t ′ ) in Eq.(A11), while the red arcs represent the static thermal system-bath correlation described byL(τ ′′ − τ ′ ) in Eq.(A12). The green arcs represent the correlation in complex time described by L(t ′ + τ ′ ), which leads to the correlated initial conditions. The contributions arising from the factorized initial conditions or correlated initial conditions consist of two parts. One is a static contribution represented by the term containing the imaginary-time integrals ofL (τ ′′ − τ ′ ) in Eq. (A12). Because of this term, the thermal equilibrium state of the system is not the equilibrium state of the system alone (pure state), but that of the combination of the system and bath (mixed state). The other is the correlated state contribution, represented by the term containing the complex time integrals of L(t ′ + iτ ′ ) in Eq.(A11). The second contribution involves the effects of the dynamical correlation and is negligible when the Markovian assumption is applied, while the first contribution always plays a significant role. It is important to note that, in addition to the fluctuation and dissipation denoted by L 2 (t) and B(t), respectively, there is a dynamical contribution from the correlated initial conditions. The role of the system-bath interaction is illustrated in terms of each contribution in Fig. 9. where , and c ′′ k = 2mζγ 3 /β(γ 2 − ν 2 k ) for k ≤ 0. At t = 0, the above equation reduces tō where ν 0 = 0,c 0 = c ′′ 0 , andc k = c ′ k + c ′′ k for 1 ≤ k, while at τ = 0, we have and As shown in Fig. 10, the noise correlation, L 2 (t), becomes negative at low temperature.
This results from the contribution of the terms with ν k = 2πk/β in the region of small t.
This behavior is characteristic of quantum noise. 8 We note that the characteristic time scale over which we have L 2 (t) < 0 is determined by the temperature and is not influenced by the spectral distribution J(ω). Thus, the validity of the Markovian (or δ(t)-correlated) noise assumption is limited in the quantum case to the high temperature regime. Approaches employing the Markovian master equation and the Redfield equation, which are usually applied to systems possessing discretized energy states, ignore or simplify such non-Markovian contributions of the fluctuation, and this is the reason that the positivity condition of the population states is broken. [41][42][43][44][45][46][47][48] As a method to resolve this problem, the rotating wave approximation (RWA) (also known as the "secular approximation") is often employed, but a system treated under this approximation will not satisfy the fluctuation-dissipation theorem, and thus the use of such an approximation may introduce significant error in the thermal equilibrium state and in the time evolution of the system toward equilibrium. Because the origin of the positivity problem lies in the unphysical Markovian assumption for the fluctuation term, the situation is better in the non-Markovian case, even within the framework of the Redfield equation without the RWA, as discussed in Sec. V. In the classical limit, with tending to zero, L 2 (t) is always positive.
While conventional approaches employing reduced equations of motion eliminate the bath degrees of freedom completely, the HEOM approach retains information with regard to the system-bath coherence in the hierarchy elements. Because of this feature, the HEOM approach can treat the reduced dynamics in a non-perturbative, non-Markovian manner.
To obtain a more compact form for the HEOM, we use the following approximate form for L 2 (t), given in Eq. (B4): 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| ≃ C k δ(t) (for k ≥ K + 1) with negligible error at the desired is the correction factor that compensates for the overestimation of L 2 (t) in the approximation at very low temperature for small cut-off K.
The accuracy of this approximation is verified on basis of the asymptotic behavior of L 2 (t) as a function of K. Then, the HEOM can be obtained by considering the time derivative of Eq. (3). 8 When the temperature becomes high (i.e. for β γ ≪ 1), the noise correlation function reduces to L 2 (t) ≃ mζγe −γ|t| /β, and hence the noise modulates the system as a Gaussian-Markovian stochastic process. 94,95 Appendix C: Derivation of the HEOM in configuration space In the present appendix, we construct the equation of motion for ρ (n) j 1 j 2 ···j K (q, q ′ ; t). In order to obtain differential equations in time, we consider the reduced density matrix elements at t + δt, ρ (n) where A is the normalization constant for the integrals over y and y ′ , and we set q = q(t) + y and q ′ = q ′ (t) + y ′ with q(t+ δt) = q and q ′ (t+ δt) = q ′ . We then expand ρ (n) j 1 j 2 ···j K (q, q ′ ; t+ δt) in terms of δt up to first order. Because y and y ′ also depend on δt, we have to expand the above equations in terms of y and y ′ . In the following, we expand the components separately.
We next consider the expansion of the factor {· · · } n in Eq.(C1), first in terms of y and y ′ , and then in terms of δt. For the expansion in y and y ′ up to second-order, we have This term reduces to ρ (n−1) j 1 ,...,j K (q, q ′ ; t), and therefore the contribution of the above to the relevant order in δt can be expressed as For the expansion of {· · · } n in terms of δt, we have We can expand the factors {· · · } j k in Eq.(C1) similarly to the {· · · } n factor. We obtain Using the definition of the hierarchy elements Eqs. (14) and (15), we obtain −nγρ (n) j 1 ,...,j K (q, q ′ ; t) + nΘ 0 (t)ρ (n−1) j 1 ,...,j K (q, q ′ ; t) δt (C13) and −j k ν k ρ (n) j 1 ,...,j K (q, q ′ ; t) − j k ν k Θ k (t)ρ

Appendix D: Time-Convolutionless (TCL) Redfield Equation
The TCL Redfield equation is the reduced equation of motion in the case of non-Markovian noise whose damping kernels are expressed in a time-convolutionless form. [49][50][51] The TCL Redfield equation is exact if the system Hamiltonian,Ĥ A , is time-independent and ifĤ A commutes with the bath interaction. However, for the BO model considered in this paper, defined by Eq.(1), the system Hamiltonian does not commute with the bath interaction.
The rotating wave approximation (RWA) is expressed asqx j = 2 /mω ′ 0 (â + +â − )(b − j + b + j ) ≈â +b− j +â −b+ j , whereâ ± andb ± j are the creation and annihilation operators of the BO oscillator and the jth bath oscillator, respectively. For j|â + |k = 0 and k|â − |j = 0, with j = k = j + 1, the interaction tensor in RWA form is given bȳ In the Wigner representation, the eigenstate elements of the density matrix are expressed as W jk (p, q) = 1 2π The total distribution is then given by W (p, q; t) = M j,k=1 ρ jk (t)W jk (p, q), where M is the number of energy eigenstates employed to solve the TCL Redfield equation.
The factorized initial state is expressed as where Z ′ A = j exp(−βE j ). By comparing the steady state solution of the TCL Redfield equation in the Wigner representation with the analytical solution of the BO model given by W eq BO (p, q) with Eqs. (45) and (46), we can check the accuracy of the steady-state distribution in the TCL formalism.