Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM)

An open quantum system refers to a system that is further coupled to a bath system consisting of surrounding radiation fields, atoms, molecules, or proteins. The bath system is typically modeled by an infinite number of harmonic oscillators. This system–bath model can describe the time-irreversible dynamics through which the system evolves toward a thermal equilibrium state at finite temperature. In nuclear magnetic resonance and atomic spectroscopy, dynamics can be studied easily by using simple quantum master equations under the assumption that the system–bath interaction is weak (perturbative approximation) and the bath fluctuations are very fast (Markovian approximation). However, such approximations cannot be applied in chemical physics and biochemical physics problems, where environmental materials are complex and strongly coupled with environments. The hierarchical equations of motion (HEOM) can describe the numerically “exact” dynamics of a reduced system under nonperturbative and non-Markovian system–bath interactions, which has been verified on the basis of exact analytical solutions (non-Markovian tests) with any desired numerical accuracy. The HEOM theory has been used to treat systems of practical interest, in particular, to account for various linear and nonlinear spectra in molecular and solid state materials, to evaluate charge and exciton transfer rates in biological systems, to simulate resonant tunneling and quantum ratchet processes in nanodevices, and to explore quantum entanglement states in quantum information theories. This article presents an overview of the HEOM theory, focusing on its theoretical background and applications, to help further the development of the study of open quantum dynamics. Published under license by AIP Publishing. https://doi.org/10.1063/5.0011599., s


I. INTRODUCTION
Time irreversibility is not a problem to be solved, but a reality to be experienced. This is true for the physical, chemical, and biological phenomena that we encounter throughout our lives. Theoretically, in particular in the quantum case, realization of time irreversibility is difficult because the fundamental kinetic equations, including the Schrödinger equation and the Dirac equation, ensure that the dynamics are reversible in time. "Open quantum dynamics" refers to the dynamics of a system that is coupled to a bath system, for example, a surrounding radiation field or an atomic or molecular environment. The bath system is typically modeled by an infinite number of harmonic oscillators. This system-bath model can describe the time irreversibility of the dynamics toward the thermal equilibrium state in which the energy supplied by fluctuations and the energy lost through dissipation are balanced, while the bath temperature does not change because its heat capacity is infinite.
Historically, studies of open quantum dynamics were motivated by practical consideration, such as line shape analysis in nuclear magnetic resonance (NMR) 1-3 and maser and laser spectra, 4,5 or by philosophical interest in thermodynamic behavior in a quantum regime as an aspect of nonequilibrium statistical physics. 6,7 In the former context, theories have been developed to construct models to describe chemical reactions, [8][9][10] electron and charge transfer rates, [11][12][13] nonadiabatic transitions, 14,15 quantum device systems, 16 ratchet rectification, 17,18 and superconducting quantum interference device (SQUID) rings 19,20 and to facilitate the analysis of linear and nonlinear laser spectra of molecules in condensed phases. 21,22 Almost independently, theories in the philosophical category have been developed mainly on the basis of the Brownian oscillator (BO) model with the aim of understanding how time irreversibility 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 made with macroscopic instruments, and why and how quantum systems approach a thermal equilibrium state through interaction with their environments. [23][24][25][26][27][28][29][30][31][32][33][34] While the possibility of treating such problems through analytical approaches is limited in the case of a Brownian-based system, researchers, in particular, in atomic spectroscopy 35 and NMR spectroscopy 36 have used a reduced equation of motion for density matrix elements, assuming that the system-bath interaction is weak (perturbative treatment) and that the correlation time of the bath fluctuations is very short (Markovian assumption). The most commonly used approaches for this kind of problem are the Redfield equation 2 and the quantum master equation. 4,5 It has been shown, however, that these equations do not satisfy the necessary positivity condition without the imposition of a rotating wave approximation (RWA) (or a secular approximation). 23,24,37,38 Because such approximations modify the form of the system-bath interaction, [39][40][41] the thermal equilibrium state and the dynamics of the original total Hamiltonian are altered. [42][43][44] In addition to these Markovian equations, phenomenological stochastic approaches have often been employed to account for the non-Markovian dephasing effects of line shape. 3,45,46 These approaches, however, ignore the effects of dissipation and are formally equivalent to assuming that the bath temperature is infinite. 42 Since the 1980s, when it became possible to investigate the ultrafast dynamics of molecular motion using nonlinear laser spectroscopic techniques such as pump-probe and photon echo measurements, the importance of the non-Markovian effects of the environment has been realized in which the noise correlation time of the environment is comparable to the time scale of the system dynamics. 21 Moreover, in many chemical physics and biochemical physics problems, the environment is complex and strongly coupled to the system at finite temperature. Thus, in addition to the Markovian approximation, other approximations also become invalid, including the rotating wave approximation, the factorization assumption, and perturbative expansions. 42 Thus, a great deal of effort has been dedicated to studying the problems of open quantum dynamics with a nonperturbative and non-Markovian system-bath interaction.
The modified Redfield 47,48 and time-convolution-less (TCL) Redfield equations [49][50][51] are reduced equations of motion that can be derived from the quantum Liouville equation by reducing the heat bath degrees of freedom. Although these approximate approaches are handy for studying problems of open quantum dynamics, their range of validity is limited. For example, they are not applicable to a system subject to a time-dependent external force because the energy eigenstates of the system incorporated in these formalisms are altered in time by the external perturbation.
The pseudomode approach 52,53 and the reaction coordinate mapping approach [54][55][56] consider equations of motion that utilize a kind of effective mode (EM) whose dynamics are described by the Markovian master equation. While these approaches have wider applicability than the conventional reduced equation of motion approaches, the description of long-time behavior that they provide may not be accurate, in particular at very low temperatures, because the Markovian master equation cannot predict the correct thermal equilibrium state. This limitation can be relaxed by introducing a pseudo-Matsubara mode, as in the case of the hierarchical equations of motion (HEOM) formalism. 57 Several variational approaches, such as the multiconfigurational time-dependent Hartree (MCTDH) approach [58][59][60] and the time-dependent Davydov ansatz (TDDA), [61][62][63] and asymptotic approaches, such as the effective-mode (EM) approach, 64,65 the density matrix renormalization group (DMRG), 66,67 and the timeevolving density matrix using the orthogonal polynomials algorithm (TEDOPA), 68,69 have been developed on the basis of a wave-function formalism for the total system. The variational approaches can be used to treat nonlinear system-bath coupling, anharmonic bath modes, 60 and a variety of Hamiltonians, for example, the Holstein Hamiltonian; [61][62][63] however, because the bath is described as a finite number of oscillators, the number of bath modes must be increased until convergence is realized to obtain the accurate results. This implies that the study of long-time behavior using these approaches requires an intensive computational effort, whereas a reduced equation of motion approach requires a numerical effort that scales only linearly with the simulation time. Strictly speaking, these wavefunction-based approaches can describe only time-reversible processes, and thus, within these approaches, there exists no thermal equilibrium state. Moreover, the inclusion of time-dependent external forces is not as straightforward in these approaches because the energy of the total system changes owing to the presence of an external force if the perturbation is strong. However, in practice, this kind of approach has wider applicability than any reduced equation of motion approach.
As theories of open quantum dynamics, on the basis of the formally exact path integral or the formally exact reduced equation of motion formalism, numerically "exact" approaches have been developed that are not subject to the limitations of the approaches discussed above. Here, numerically "exact" indicates the ability to calculate the dynamical and thermal aspects of a reduced system with any desired accuracy that can be clearly verified through non-Markovian tests on the basis of exact analytical solutions.
The path-integral approaches, most notably path-integral Monte Carlo, are computationally intensive because the number of paths to be evaluated grows rapidly with time, while sampling often fails owing to the phase cancellation of wave-functions. 70,71 Much effort has been expended in attempting to overcome these problems, for example, by using the quasi-adiabatic path-integral (QUAPI) algorithm, [72][73][74][75][76][77][78][79] and to extend the applicability of the path-integralbased methods. Because these approaches can easily incorporate a semiclassical approximation for the bath 80 or introduce a modular scheme to effectively separate the system part and the environmental part, 81 they have advantages for the study of polyatomic systems treated in multidimensional coordinates.
Several equations of motion approaches that explicitly handle time-dependent random variables for the heat bath have also been developed. [82][83][84][85] These approaches are formally exact, but the realization of the random variables in the low-temperature regime is difficult, and the applicability of these approaches is still limited to simple systems.

PERSPECTIVE scitation.org/journal/jcp
The reduced hierarchical equations of motion (HEOM) theory is a method that can describe the dynamics of a system with a nonperturbative and non-Markovian system-bath interaction at finite temperature, even under strong time-dependent perturbations. [42][43][44][86][87][88][89][90][91][92] In this formalism, the effects of higher-order non-Markovian system-bath interactions are mapped into the hierarchical elements of the reduced density matrix. This formalism is valuable because it can be used to treat not only strong system-bath coupling but also quantum coherence or quantum entanglement between the system and bath ("bathentanglement"), in particular, for a system subject to a time-dependent external force and nonlinear response functions. [93][94][95] Although the HEOM approach is numerically very expensive in comparison with the other reduced equations of motion approaches, a variety of analytical and numerical techniques can be employed to integrate the HEOM. With these features, the HEOM approach offers wide applicability. For example, it is possible to study quantum heat-engine and quantum ratchet problems in which the nonequilibrium steady-state is described, respectively, by the long-time behavior of a system strongly coupled with two heat baths at different temperatures and by the long-time behavior of a system under periodic time-dependent external forces, as illustrated in Secs. II-VI. In this article, we present an overview of the HEOM theory, focusing on its theoretical background and applications, to help further the development of investigations in this field.
The organization of the remainder of this paper is as follows: In Sec. II, we present typical model systems for open quantum dynamics. In Sec. III, we explain the standard HEOM and their characteristic features. In Sec. IV, we demonstrate the accuracy of the HEOM by numerically "exact" tests. In Sec. V, we illustrate the variety of HEOM for various systems. In Sec. VI, we review various applications of the HEOM theory. We present future perspectives for HEOM theory in Sec. VII.

II. SYSTEM
We consider a situation in which a primary system interacts with heat baths that give rise to dissipation and fluctuation in the system. [25][26][27][28][29][30][31][32][33] The total Hamiltonian is expressed aŝ whereĤA is the Hamiltonian of the system andĤI+B is the bath Hamiltonian, which includes the system-bath interaction. The bath degrees of freedom (labeled a) are treated as Na harmonic oscillators, where the momentum, position, mass, and frequency of the jth oscillator in the ath bath are given byp a j ,x a j , m a j , and ω a j , respectively, V a is the system part of the interaction, and α a j is the coupling constant between the system and the jth oscillator in the ath bath. The above bath model is commonly applied to systems possessing discretized energy spaces for which it is expressed asĤA = ∑j ̵ hωj| j⟩⟨ j| + ∑j≠k ̵ hΔ jk (t)| j⟩⟨k| andV a = ∑ j,k V a jk | j⟩⟨k|, which are defined using the bra and ket vectors of the jth eigenenergy states of the system | j⟩ and ⟨ j|. A notable application of this kind is to the spin-boson Hamiltonian for j = 0 or 1. [30][31][32][33] The model is also commonly applied to systems possessing continuous configuration spaces for which it is expressed in general form asĤA =p 2 /2m + U(q; t) for a system with mass m and potentials U(q) andV a =V a (p,q) described in terms of the multidimensional momentump and coordinatê q. [25][26][27][28][29] A notable example of such an application is to the Brownian Hamiltonian. In the Brownian case, to avoid an unphysical energy divergence, a counterterm ∑ a ∑ j (α a j V a ) 2 /2m a j (ω a j ) 2 is introduced in the bath Hamiltonian such that the replacementx a j →x a j ≡x a j − α a jV a /m a j (ω a j ) 2 is made 25,26 to maintain translational symmetry in the case U(q) = 0. 42,89 In atomic spectroscopy, this divergence phenomenon is known as the Lamb shift. 4,5 In many physical processes, molecular motion and electronic excitation states are coupled and play important roles simultaneously. Using this bath model, we can further include the effects of the electronically excited states by extending the space of the system Hamiltonian asĤA = ∑ j | j⟩(p 2 /2m)⟨ j| + ∑ j,k | j⟩U jk (q; t)⟨k|. 13,21 With the features described above, the Brownian heat bath possesses wide applicability, despite its simplicity. This is because the influence of the environment can, in many cases, be approximated by a Gaussian process owing to the cumulative effect of a large number of weak environmental interactions in which case the ordinary central limit theorem is applicable, 46 while the distribution function of the harmonic oscillator bath itself also exhibits a Gaussian distribution. 42,46 Hereinafter, we assume that each bath is independent of the others. The ath heat bath can be characterized by the spectral distribution function (SDF), defined by and the inverse temperature β ≡ 1/k B T, where k B is the Boltzmann constant. Various environments, for example, those consisting of nanostructured materials, solvents, and/or protein molecules, can be modeled by adjusting the form of the SDF. For the heat bath to be an unlimited heat source possessing an infinite heat capacity, the number of heat bath oscillators Na is effectively made infinitely large by replacing Ja(ω) with a continuous distribution. When we reduce the bath degrees of freedom to have the reduced density matrixρA(t) = trB{ρA+B(t)}, the bath acts as a noise source for the system via the system-bath interaction. 42,86 The random variable defined asΩ(t) ≡ ∑ j αjxj describes the thermalization of the system through fluctuation and dissipation. In the case of a harmonic bath, the contributions of higher-order cumulants vanish, and the effects of dissipation are expressed in terms of the response function defined as L 1 (t) = i⟨Ω(t)Ω −ΩΩ(t)⟩ B / ̵ h, while those of thermal fluctuations are expressed as L 2 (t) = ⟨Ω(t)Ω +ΩΩ(t)⟩ B /2, where ⟨⋯⟩B represents the thermal average of the bath degrees of freedom. 31,[42][43][44] The thermal equilibrium state of the system is realized through the energy exchange arising from fluctuation and dissipation: the condition in which there is a thermal equilibrium state is completely described by the quantum version of the fluctuation-dissipation theorem, because of the Gaussian nature of the noise, which is equivalent to assuming a harmonic heat bath. 32,86 The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcp

III. REDUCED HIERARCHICAL EQUATIONS OF MOTION
To illustrate the characteristic features of the HEOM formalism in a simple manner, in this section, we consider a Drude SDF (an Ohmic distribution with a Lorentzian cutoff) for a single bath model, expressed as [42][43][44][86][87][88][89][90][91][92] where the constant γ represents the width of the spectral distribution of the collective bath modes and is the reciprocal of the relaxation time of the noise induced by the bath. The parameter η is the systembath coupling strength, which represents the magnitude of fluctuation and dissipation. By taking γ → ∞, the Drude form reduces to the Ohmic form J(ω) = ̵ hηω/π, which is often employed to obtain Markovian noise in the high-temperature case.
The relaxation function is temperature-independent and is expressed as L 1 (t) =c 0 e −γ|t| , withc 0 = ̵ hηγ 2 , while the noise correlation is given by L 2 (t) = c 0 e −γ|t| + ∑ k c k e −ν k |t| , with c 0 = ̵ hηγ 2 cot(β ̵ hγ/2)/2 and c k = −2ηγ 2 ν k /β(γ 2 − ν 2 k ) for the kth Matsubara frequency ν k ≡ 2πk/β ̵ h. 42,[86][87][88][89][90] A. The positivity condition The Markovian assumption is commonly employed in quantum open dynamics. Because the Markovian assumption is unrealistic in the physical problems owing to the ignorance of the non-Markovian features of quantum fluctuations, we encounter a breakdown of the positivity condition of the population states if we do not employ further approximations, most notably the rotating wave approximation (RWA) (or a secular or resonant approximation). 23,24,[37][38][39][40][41] To illustrate this point, we plot L 2 (t) in Drude cases for large γ in Fig. 1. 42,44 The profile of the response function L 1 (t) is similar to that of L 2 (t) in the high-temperature case and becomes a δ function [∝ δ(t)] as γ → ∞. At low temperatures, however, the noise correlation L 2 (t) becomes negative owing to the contribution from terms with c k e −ν k t for small k in the region of small t. This behavior, which we call low-temperature-induced non-Markovianity, is a characteristic feature of quantum thermal noise. 42,44 This indicates that the validity of the Markovian assumption in the quantum case is limited only in the high-temperature regime. Approaches that employ the Markovian master equation and the Redfield equation, which are usually applied to systems possessing discretized energy states, ignore or simplify such non-Markovian features of fluctuations, and this is the reason why the positivity condition is broken. [42][43][44] As a way to resolve this problem, the RWA is often employed, 39-41 but a system treated under this approximation will not satisfy the fluctuation-dissipation theorem. Thus, the use of such an approximation may introduce significant errors in the thermal equilibrium state and in the time evolution of the system toward equilibrium, as will be illustrated in Sec. IV. In the classical limit, with ̵ h tending to zero, L 2 (t) is always positive.
As can be seen from the form of L 1 (t) and L 2 (t), there are two types of non-Markovian components, as illustrated in the Drude case: one is of mechanical origin and is characterized by dissipation [i.e., L 1 (t)] and fluctuations [i.e., L 2 (t)] expressed in terms of e −γt , and the other is of quantum thermal origin and is characterized by fluctuations only [i.e., L 2 (t)], with the functions expressed in terms of c k e −ν k t for k ≥ 1, as the definition of ν k = 2πk/β ̵ h indicates. When γ is much larger than ν 1 , the mechanical contribution with e −γt vanishes for t > 1/ν 1 , and then, the effects from the quantum thermal noise arise. The quantum thermal fluctuations exhibit peculiar behavior compared with the mechanical fluctuations because the amplitude of the noise becomes negative for large γ, as illustrated in Fig. 1. While the mechanical noise term contributes to the relaxation and dephasing processes through fluctuation and dissipation, the quantum thermal noise term contributes to the dephasing process. This difference can be observed even in the free induction decay of a spin-boson system, if γ is large and the temperature is very low. 96

B. HEOM for density operators
The HEOM are derived by differentiating the reduced density matrixρA(t), defined by a path integral [42][43][44][86][87][88][89] or by time-ordered operators with the application of Wick's theorem. 97 The HEOM can describe the dynamics of the system under nonperturbative and non-Markovian system-bath interactions with any desired accuracy at finite temperature. [42][43][44] In their original formulation, these equations of motion were limited to the case in which the SDF takes the Drude form and the bath temperature is high. 86 Subsequently, these limitations have been removed, and the equations are applicable to arbitrary SDFs and temperatures. 87,89 While conventional approaches employing reduced equations of motion eliminate the bath degrees of freedom completely, the HEOM approach retains information with regard to bath dynamics, including system-bath coherence or system-bath entanglement (bathentanglement) in the hierarchy elementsρ (n) j 1 ,...,j K , where the set of indices n and {j k } arise from the hierarchical expansion of the noise correlation functions in terms of e −γt and e −ν k t , respectively. [42][43][44] The zeroth element is identical to the reduced density matrix,ρA(t) =ρ where ω 0 is the characteristic frequency of the system and K is a cutoff satisfying K ≫ ω 0 /ν 1 , the HEOM are expressed as 42,43,[86][87][88][89] ∂ ∂tρ is the renormalization operator to compensate for the effects of the cutoff K, and for an integer k ≥ 1. Here, the hyper-operatorsÔ ×f ≡Ôf −fÔ andÔ ○f ≡Ôf +fÔ for any operator Ô and operandf are those introduced by Kubo 98 and Tanimura, 86 respectively.
In this formalism, the effects of the higher-order non-Markovian system-bath interactions are mapped into the hierarchical elements of the reduced density matrix (see Fig. 2). The operator Φ demolishes the bath states fromρ involve the first-order higher (−) and lower (+) interactions thanρ (n) j 1 ,...,j K (t). Thus, one can regard Eq. (5) as a type of rate law among the density operators with three different bath-excited states,ρ (n±1) is determined by its time evolution operator (LA + nγ + ∑ K k=1 j k ν k +Ξ) and the incoming and outgoing contributions of the hierarchical elements connected throughΦ,Θ 0 , and depict the elements of the first and second members of the hierarchy, which involveρ (1) 0,0,...,0 (t) andρ (2) 0,0,...,0 (t). In each diagram, the left and right solid lines represent the time evolution of the system, and the red dashed lines represent the bath excitations. The white and black circles correspond to the system-bath interactions involved inΦ andΘj (j = 0, 1, . . .), respectively. Here, the quantum entanglement between the system and the bath (bathentanglement) is illustrated by the red dashed curves. A detailed explanation can be found in Refs. 42-44. Θ k . 42 The HEOM approach can treat the noise correlation even if it becomes negative (as shown in Fig. 1) in a nonperturbative and non-Markovian manner because this approach retains the information on the system-bath correlation in the hierarchy elements with the set of indices n and {j k }. Moreover, these hierarchical elements allow us to evaluate the bathentanglement and were utilized to calculate the heat flow between the system and the bath. [99][100][101][102] The HEOM approach differs conceptually from conventional perturbative expansion approaches, where the zeroth member does not include any system-bath interactions and where higher members take into account higher-order system-bath interactions: in the HEOM approach,ρ (0) 0,0,...,0 (t) includes all orders of the systembath interactions, and it is the exact solution of the Hamiltonian, Eq. (1). [42][43][44][86][87][88][89][90] As an example, a numerical solution of the HEOM for a three-level heat engine model is presented in Fig. 3. 101,102 C. Quantum hierarchical Fokker-Planck equations For a system described in coordinate and momentum space, it is often more convenient to use the phase space representation than the energy eigenstate representation, in particular to implement various boundary conditions, most notably periodic, open, and inflow-outflow boundary conditions, and to take the classical limit to identify pure quantum effects. For the density matrix element ρ (n) j 1 ,...,j K (q, q ′ ; t), we introduce the Wigner distribution function, which is the quantum analog of the classical distribution function in phase space and is defined as 103,104 The Journal of Chemical Physics The Wigner distribution function is a real function, in contrast to the complex density matrix. In terms of the Wigner distribution, the quantum Liouvillian for the system Hamiltonian takes the form 104 −LQMW where The quantum hierarchical Fokker-Planck equations (QHFPE) with the counterterm can then be expressed as 44,105-108 where the dashed operators are the Wigner representations of the operators in Eq. (5). For linear-linear system-bath coupling, V(q) = q, they are expressed aŝ As an example, a numerical solution of the QHFPE for a self-current oscillation of a resonant tunneling diode is presented in Fig. 4. [105][106][107] The QHFPE for linear-linear (LL) plus square-linear (SL), V(q) = q + cq 2 , and exponential-linear (EL), V(q) = e −cq , system-bath coupling cases for any constant c are given in Refs. [109][110][111][112][113][114], respectively. It should be noted that because the counterterm is incorporated into the equations of  In the LL high-temperature case, we can ignore the lowtemperature correction terms, and the above equations for where ζ ≡ η/m. These equations are a generalization of the Caldeira-Leggett quantum Fokker-Planck equations (QFPE) 116,117 for non-Markovian noise. 89,90 The classical limit of the equations of motion can be obtained by taking the limit ̵ h → 0, which leads to the replacement of the Liouvillian throughL QM →L cl (p, q) ≡ (p/m)(∂/∂q) +f (q)(∂/∂p), where f (q) = ∂U(q)/∂q. The classical limit of Eq. (9) is the classical hierarchical Fokker-Planck equations (CHFPE), which represent a generalization of the Kramers equation 118,119 for non-Markovian noise. 89,90 Because β ̵ hγ = 0 is always satisfied in the classical case, the validity of the CHFPE approach is not limited to the high-temperature regime in the classical case.

D. Continued fraction expression and truncation of the hierarchical equations
The hierarchy of equations introduced above continues to infinity and is "formally exact," but these equations are not easy to solve numerically. If the system-bath coupling η is small, the contribution of the elements of higher members of the hierarchy becomes smaller than that of the elements of the lower members; thus, we can safely disregard the deeper hierarchy. This is not the case, however, if the system-bath interaction is strong. In the high-temperature case, we can solve the HEOM analytically using the continued fraction expressions for the resolvent of a simple system. 86,87 Under the high-temperature condition of Eq. (5), the Laplace-transformed ele- where the nth resolvent is defined in recursive form as 42,120 Gn[s] = n + 1 Nonlinear response functions can also be expressed as products of Gn[s] (see Fig. 5), which allows us to evaluate the bathentanglement analytically. 120-123 For a complex system, most notably a system described in coordinate space, or a system driven by a time-dependent external force, it is easier to solve the equations of motion with a finite number of the elements K by adopting one of the following truncation schemes. The first of these is the time-nonlocal scheme utilizing the asymptotic relation among the higher-order hierarchy members: As can be seen from Eq. (10), we haveĜK[s] ≈ (K + 1)/[s + (K + 1)γ] for (K + 1)γ ≫ ω 0 , where ω 0 is the characteristic frequency of the system. 42,89,90,124 The second is the time-local scheme ignoring the higher-order equations of motion with higher Matsubara frequency FIG. 5. The ultrafast nonlinear laser spectrum is calculated from the nonlinearresponse functions. In a high-temperature case, the third-order optical response can be expressed analytically in terms of the HEOM resolvent in continued fraction form. 120,121 Here, a chemical reaction process is discussed for a proton transfer system (a) coupled to a bath. Then, the two-dimensional (2D) vibrational spectrum I(ω 1 , t 2 , ω 3 ) is calculated for different t 2 (see the Appendix). A snapshot of the 2D spectrum is presented in (b). The diagonal peaks arise, respectively, from the transitions A and B depicted in (a). The volume of the off-diagonal peaks has been shown to be proportional to the transition rate between the ground and the first excited state. By evaluating this volume as a function of T 2 , we can evaluate the chemical reaction rate directly from the measurement (c). Reproduced with permission from A. Ishizaki  terms, while a renormalization operator is introduced to take into account the effects of the cutoff. 88 The truncated HEOM can be evaluated with the desired accuracy by depicting the asymptotic behavior of the hierarchical 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 becomes negligibly small when K is sufficiently large.
Here, we illustrate their hybrid truncation scheme for the quantum Fokker-Planck equation (QFPE). 44 In this scheme, we set the number of Matsubara frequencies to be included in the calculation as K, which is chosen to satisfy K ≫ ω 0 /ν 1 . The upper limit for the number of hierarchy members for given γ is 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 hierarchical equations by replacing Eq. (8) with while, for the case n = Kγ, we employ The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcp is the renormalization operator in the Wigner representation. The time-local scheme corresponds to Eq. (11) and the time-nonlocal scheme to Eq. (12). We can then evaluate W (n) j 1 ,...,j K (p, q; t) through numerical integration of the above equations. In the Markovian limit γ → ∞, which is taken after the high-temperature limit, yielding the condition β ̵ hγ ≪ 1, we have the QFPE 116,117 ∂ ∂t which is identical to the quantum master equation without the RWA. 42 Because we assume that the relation β ̵ hγ ≪ 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. To overcome this limitation, we must include low-temperature correction terms to obtain the low-temperature-corrected QFPE (LT-QFPE). 125

E. Multistate quantum hierarchical Fokker-Planck equations
Up to this point, we have separately discussed two types of dynamics employing distinct approaches: one is a discrete energylevel system discussed from the viewpoint of the HEOM, and the other is a molecular or atomic coordinate system discussed from the viewpoint of the QHFPE. As mentioned in Sec. II, molecular motion and electronic excited states are coupled and play important roles simultaneously in many important physical processes. In the HEOM approach, this extension is straightforward and presented as the multistate QFPE, which can be applied to both optical and nonadiabatic transition problems, taking into account nonperturbative and non-Markovian system-bath interactions. In this case, the Wigner distribution function is expanded in the electronic basis set as [126][127][128][129] where W jk (p, q) is the jk element of the Wigner distribution function. We represent the Wigner functions for different electronic states in matrix form as W (p, q) ≡{W jk (p, q)}. Similarly, the transition operator is written in Wigner representation form as The quantum Liouvillian can also be expressed in matrix form as 126 By using the above transformation, the hierarchy operatorsΦ ′ ,Θ 0 , Θ ′ k , andΞ ′ in Eq. (8) are cast in matrix form asΦ ′ ,Θ 0 ,Θ ′ k , and Ξ ′ for the system-bath interaction defined as V (q) = {V jk (q)}. In this way, we can obtain the multistate quantum hierarchical Fokker-Planck equations (MS-QHFPE) for W (n) j 1 ,...,j K (p, q). [128][129][130][131][132] Comparisons among the transient absorption spectrum (TAS) calculated from the Ohmic MS-QHFPE with low-temperature correction terms (LT-MS-QFPE), the Zusman equation, fewest switch surface hopping, and the Ehrenfest approach were made to demonstrate the reliability of the numerically "exact" approach. 125 A numerical solution of the low-temperature multistate quantum Smoluchowski equation (LT-MS-QSE), which is the strong-damping limit of LT-MS-QFPE, is presented in Fig. 6.
A quantum electron transport problem described by the spinless Anderson-Holstein model was also investigated in a similar manner. 133

F. The bathentanglement states and the positivity of the HEOM
In the HEOM formalism, the information concerning the system-bath coherence (bathentanglement) is stored in the hierarchical elements, which allows us to simulate the quantum entangled dynamics between the system and bath (see Fig. 2). Thus, for example, the correlated initial equilibrium state can be set by running the HEOM program until all of the hierarchy elements reach a steady-state and then use these elements as the initial state: the steady-state solution of the first hierarchy elementρ (0) 0,0,...,0 (t = ∞) agrees with the correlated thermal equilibrium state defined bŷ ρ eq A = trB{exp(−βĤtot)}, while the other elements describe the bathentanglement states. Because the conventional Markovian and non-Markovian reduced equations of motion approaches eliminate the bath degrees of freedom completely, they cannot properly take into account the bathentanglement states, as illustrated in Fig. 7. 44,99,134,135 Although the effects of bathentanglement are not easy to observe as long as we are working on linear response measurements, it is possible to study them experimentally by measuring the nonlinear response of the system (see the Appendix).
Because of their structure, the HEOM can treat any noise correlation even it becomes negative (as shown in Fig. 1), and the HEOM continue to satisfy the positivity condition, as demonstrated by various numerical simulations (see also Sec. IV). Owing to the complex structure of the HEOM, however, analytical verification of their positivity has so far been limited to a simple case. 136 To obtain a numerically converged solution as the thermal equilibrium state, it is important to express the noise correlation function in terms of the decaying functions as cj exp(−ζjt) or cj exp(−ζj ± iωj), where cj, ζj, and ωj characterize the noise correlation for the jth hierarchical elements, to maintain the positivity of the HEOM formalism. [137][138][139][140][141][142][143][144][145][146][147][148] Although the HEOM can be constructed . Because of the factorization assumption, the TCL approach cannot take into account the bathentanglement at the time at which the external forces are applied to the system (the blue arrows with black dashed lines), whereas the HEOM can take these entangled states into account accurately by utilizing the hierarchy elements illustrated in Fig. 2. The simulation results corresponding to (a-i) and (a-ii) are presented in Figs. 8(d-i) and 8(d-ii), and those corresponding to (b-i) and (b-ii) are presented in Figs. 14(a) and 14(b), respectively. The blue dashed curves in (b-i) represent the bathentanglement states that contribute to the rephrasing echo signal in the 2D spectrum, as presented in Fig. 14(a).
hierarchical structure, in particular in the low-temperature case, numerical integration of the HEOM is computationally intensive in terms of both memory and central processing unit (CPU), although, in many cases, there is no other way to obtain reliable results. Therefore, great efforts have been made to reduce the computational costs by improving the algorithmic and numerical techniques employed. An efficient approach is to suppress the number of elements in the hierarchy by constructing the HEOM using Padé [157][158][159] or Fano 160,161 spectral decompositions instead of Matsubara frequency decomposition. By introducing rescaling factors, it is possible to reduce the number of elements in the hierarchy. 162 A numerical algorithm based on the optimization of hierarchical basis sets 163 and a tensor network has also been examined with the aim of reducing the number of calculations required. 164,165 Like other reduced equation approaches, the HEOM are formulated in terms of the reduced density matrix, which requires N × N memory space for a system with N eigenstates. This makes the scalability of the system size very low, in particular when the system is described in Wigner space. Therefore, wavefunction-based HEOM approaches whose scalabilities are similar to that of the Schrödinger equation have been developed [166][167][168][169][170][171][172][173][174][175][176][177] (see Sec. V C). In addition to the scalability of the memory, these approaches are advantageous because they allow us to employ the wide range of numerical techniques developed for the Schrödinger equation.
Advances in computer technology have also led to HEOM calculations becoming faster and larger. Thus, computer codes for a graphic processing unit (GPU) and parallel computing through the Message Passing Interface (MPI) utilizing distributed memory (DM) have been developed to treat large systems and various SDFs. Commonly used HEOM of this kind are the GPU-HEOM 178-180 and DM-HEOM. [181][182][183][184] To develop an integral routine for the HEOM for a time step Δt, an easy and efficient way is to employ the math libraries developed for targeting computer architectures. For example, to utilize the GPU architecture, one can construct the large matrices for the entire HEOM reduced density operators and the Liouvillians, and then manipulate them using the math library suitable for these matrices. 180 The most efficient time-integration routine currently available is the low-storage Runge-Kutta method, 185 which has been used to study multistate nonadiabatic dynamics described in two-dimensional configuration space. 129

IV. NUMERICALLY "EXACT" TESTS FOR NON-MARKOVIAN AND NONPERTURBATIVE DYNAMICS
The capabilities of the reduced dynamics theory can be verified through nonperturbative and non-Markovian tests on the basis of analytical solutions for the Brownian oscillator model with arbitrary SDF, J(ω). 44 These tests can be applied to the system described in the coordinate representation,ĤA =p 2 /2m + mω 2 0q 2 /2, or to the system in the energy eigenstate representation,ĤA = ̵ hω 0 (â + â − + 1/2), where â ± are the creation-annihilation operators for the eigenstates. The tests are based on the solutions for (a) the steady-state distribution, (b) the symmetric autocorrelation function C(t) ≡ ⟨q(t)q +qq(t)⟩/2, (c) the linear response function R (1) (t) ≡ i⟨[q(t),q] ⟩/ ̵ h, and (d) the nonlinear response function R  186,187 which is the observable of 2D terahertz Raman spectroscopy 188 (the simulation method is explained in the Appendix). These are tests of the ability of the theory to account for (a) the bathentanglement thermal equilibrium state, (b) fluctuations, (c) dissipation, and (d) dynamical bathentanglement. Test (d) is particularly important if we wish to study dynamics under time-dependent external forces. 130,131,137,189 Here, we illustrate these results for the Drude case [Eq. (4)] calculated from analytically exact expressions for the Brownian oscillator system, 28-30 the QHFPE, 44 TTR [ω1, ω2]. 186,187 The results obtained from the QHFPE approach (i) almost overlap with the analytical result, while the results from the TCL Redfield approach (ii) miss some peaks owing to the factorized nature of the TCL description, as illustrated in Fig. 7

PERSPECTIVE
scitation.org/journal/jcp the original factorized initial state: the difference between the two distributions arises from "static bathentanglement" and represents the nonfactorized effect of the thermal equilibrium state. Figure 8(b) depicts the autocorrelation (symmetric correlation) functions in the frequency domain. The red-and blue-dotted, and blue-dashed curves represent the results obtained from the analytical expression, the QHFPE, the TCL Redfield equation, and the TCL Redfield equation with the RWA, respectively. Figure 8(c) shows the linear response function R (1) [ω] for strong system-bath coupling, ζ = 3.0. This function is temperatureindependent in the harmonic case. While the QHFPE results (red curves) coincide with the exact results (black dots), the TCL Redfield results without the RWA (blue curves) and with the RWA (blue dashed curves) are close only near the maximum peak, regardless of temperature. The low-frequency parts of the spectra arise from the slow dynamics of the reduced system near the thermal equilibrium state, and the discrepancy between the TCL results and the exact results arises from the slow equilibration process. Figure 8(d) shows the second-order response function R (2) TTR [ω 1 , ω 2 ] corresponding to the intermediate coupling case considered in Fig. 8(b). The results here were obtained from (i) the QHFPE approach and (ii) the TCL Redfield approach without the RWA. The QHFPE results overlap with those from the analytical expression. Owing to "dynamical bathentanglement," we observe peaks near (ω 1 , ω 2 ) = (0, 1) and (1, 1), whereas the conventional reduced equation of motion approaches cannot take system-bath coherence into account owing to their use of a factorized description of the bath state, as shown in Fig. 7(a).
As illustrated in this section, we have been able to obtain very accurate results from the HEOM approach, as judged by tests (a)-(d). By contrast, we have found that the TCL Redfield equation has limited applicability when judged in a similar way.

A. HEOM for arbitrary spectral distribution functions
To extend the applicability of the HEOM in a numerical "exact" manner, various HEOM have been developed. A straightforward extension of the HEOM approach is to consider a more complicated SDF to describe the environmental effects of realistic chemical and biochemical systems. When the noise correlation function is expressed in terms of damped oscillators, which can be done for the Brownian, [137][138][139][140] Ohmic, 125 Lorentz, 141 super-Ohmic, 142 and Drude-Lorentz 143 SDFs (and combinations of these; see Refs. 144 and 145), we can obtain the reduced equations of motion in the hierarchical structure. Alternatively, we can extend the HEOM for arbitrary SDFs using Fourier, 89 Gauss-Legendre, or Chebyshevquadrature spectral decompositions (eHEOM). [149][150][151][152][153][154] This approach allowed the investigation of a spin-boson system coupled to a bath with sub-Ohmic SDF at zero temperature. 153,154 An extension of the HEOM for a combination of exponential/non-exponential bath correlation functions has also been presented. 190

B. Stochastic HEOM
The difficulty of solving the HEOM arises from the fluctuation terms described by L 2 (t), which give a negative contribution at low temperatures (Fig. 1). In the classical generalized Langevin formalism, dissipation is described by a damping kernel, while fluctuations are described by noise whose correlation function is related to the damping kernel through the fluctuation-dissipation theorem. Although the trajectories of a particle under negatively correlated noise cannot be evaluated using the Langevin formalism, we can evaluate the distribution function of the trajectories by extending the HEOM formalism. 42 Thus, various stochastic equations of motion have been derived in which a hierarchical structure is employed for the dissipation part, while the fluctuation part is treated as noise. [191][192][193][194] The efficiency of this approach relies on its treatment of the noise correlation. A methodology that treats the exponentially decaying part of the noise using a hierarchical formulation in the same manner as the dissipation has been developed (the stochastic HEOM). Because this approach does not depend on a hierarchy of Matsubara frequencies or a Padé decomposition, the memory requirement of computations is dramatically reduced in comparison with that of the regular HEOM, in particular for a system interacting with multiple heat baths. An extension of the stochastic approach to non-Gaussian noise that includes noise from a spin bath has also been suggested. 194 However, in the stochastic approach, although the memory requirement is reduced, sampling of the trajectories becomes difficult for the lower-temperature case or for longer-time simulations, in particular when the system is driven by a time-dependent external force. Thus, under such conditions, the accuracy of calculations becomes lower than that of the regular HEOM.

C. Wave-function-based HEOM
As shown in the derivation of the influence functional with a correlated initial condition, the Feynman-Vernon influence functional can also be defined on the basis of the wave-function using complex time counter integrals. 43 This indicates that we can derive HEOM for a wave-function: such approaches include the stochastic hierarchy of pure states (HOPS), [166][167][168][169] the stochastic Schrödinger equation (SSE), 170 the hierarchy of stochastic Schrödinger equations (HSSE), [171][172][173][174][175][176] and the hierarchical Schrödinger equations of motion (HSEOM), 177 all of whose scalabilities are similar to that of the Schrödinger equation. These equations are advantageous not only because the scale of the memory required becomes N for an N-level system but also because various numerical techniques developed for the Schrödinger equation can be employed for configuration space or energy eigenstate representations or a mixture of these. Technically, these wave-function-based approaches are not regarded as equation of motion approaches because they are not defined by the time t but by the complex time counter path 0 Fig. 9). Thus, whereas in the conventional HEOM approach, the HEOM density operatorsρ (n) j 1 ,...,j K (t+Δt) are computed fromρ (n) j 1 ,...,j K (t), in the wave-function-based HEOM approach, we must repeat the full integration for t along −iβ ̵ h for different Δt. Moreover, because the convergence of the trajectories becomes slow in the stochastic approaches and because the required Chebyshev function set becomes larger for longer simulations in the HSEOM approach, the wave-function-based approach is not suitable for studying a system with slow relaxation or a system subjected to a slowly varying time-dependent external force. This indicates that the regular HEOM are the method of choice to study long-term dynamics if the system is not too large, while the wave-function-based HEOM give access to short-and intermediatetime dynamics over a wider range of temperatures for a large system.

D. HEOM for a grand canonical electron bath
Although the HEOM were originally obtained by considering a harmonic heat bath in a canonical ensemble, Jin and his collaborators 195 showed that the HEOM could also be derived from Gaussian-fermionic and Gaussian-bosonic bath models in a grand canonical ensemble by introducing the chemical potential μ to serve as a second thermodynamic variable in addition to the temperature. The HEOM for an electron bath are appropriate for studying charge transport problems with electron-electron interactions, most notably a quantum dot and a molecule coupled to two electrodes and representing a continuum of noninteracting electronic states 196 (see Fig. 10).
We consider a HamiltonianĤA described in terms of the creation and annihilation operators of the jth electronic statedj;σ and d † j;σ , where σ represents the spin state. The electrode bath is expressed as whereĉ k;σ ′ andĉ † k;σ ′ are the creation and annihilation operators of an electron in an electrode with wavevector k and spin σ ′ . The systemelectrode interaction is expressed aŝ where f represents the site of the system that is coupled to an electrode and V α k;σ ′ describe the tunneling rates for the left and right electrodes, with α = L and R, respectively. The electrode system is then characterized by the tunneling efficiency function, which plays a similar role to the harmonic SDF.
For the sake of conciseness, here, we consider a spinless electronic bath where the system-electrode coupling and the bandwidths of the electrodes are assumed to be symmetric with respect to both sides of the lead. Then, the effective tunneling distribution is assumed to have a simple Drude form, where ζα is the coupling strength and μα is the chemical potential for α = L and R. This allows us to express the electronic bath relaxation function in the form with γ α,s for m ≥ 1. The index s ∈ { +, −} corresponds to creation and annihilation of electrons, i.e.,d + j;σ =d † j;σ andd − j;σ =dj;σ. Thus, the HEOM for a system coupled to the electrode baths can be expressed as [196][197][198] ∂ ∂tρ a 1 ,...,a n (t) γ α k ,s k m k )ρa 1 ,...,a n (t) − i ∑ a [d¯s j;σρa1,...,an,a (t) − (−) nρ a 1 ,...,a n , a(t)ds j;σ ] − i n ∑ k=1 (−) n−k [η α k ,s k m kd s k j k ;σ kρ a 1 ,...,a k−1 ,a k+1 ,..., a n (t) + (−) n (η α k ,s k m k ) * ρ a 1 ,...,a k−1 , a k+1 ,..., a n (t)d s k j k ;σ k ], The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcp where a k = {j k , α k , σ k , m k , s k } and σ k is the index of the system side of the spin in the j k th electronic state, m k is the index for the Matsubara frequencies, and s k can be + or − withs k = −s k . Because the higher-order cumulants vanish in the treatment of a Gaussian-electron source, it is not necessary to utilize the influence functional. 199 At first glance, the above HEOM appear to be similar to the conventional HEOM, but their structure is more complex than those of the regular HEOM because the number of elements in this kind of hierarchy quickly increases as the system size increases, analogously to the HEOM derived on the basis of Fourier decomposition. 87 Extensions to take account of a realistic band structure 200,201 and a truncation scheme specific for a fermionic bath have been investigated. 202 A computer code for a fermionic bath, HEOM-Quick, has been developed and distributed. 203 The electrode heat bath considered here is a noncorrelated electron environment based on a Gaussian single-particle picture, where fermionic fluctuationdissipation at second order can be applied. This condition is convenient when a density functional treatment of the system is employed. Thus, the real-time evolution of the reduced single-electron density matrix at the tight-binding level has been studied by combining time-dependent density functional theory with HEOM theory (TDDFT-HEOM). 204,205 An extension to include electronic-vibrational coupling has also been presented. 206 The expression for the HEOM in this case was derived using a small polaron transformation to take account of strong electronic-vibrational coupling. [207][208][209]

E. HEOM for different system-bath models
The HEOM approach has also been applied to the case of solid state materials described by a deformation potential 87 and the Holstein Hamiltonian. 210,211 In the Holstein case, because we can study only a small system and the number of phonon modes associated with the system site is finite, special treatment is necessary to maintain the stability of the equations. 211 The HEOM for a molecular rotor system in which the rotational symmetry of the system is maintained have also been derived. 212,213 This extension is important to account for the rotational bands of linear and nonlinear spectra.

F. Phenomenological and approximate approaches
Although it is not possible to construct a reduced equation of motion in the HEOM structure for various system-bath Hamiltonians through a procedure of self-induction, one can still adopt the HEOM structure in a phenomenological approach. 214,215 For example, if the heat bath is not harmonic or not Gaussian owing to the anharmonicity of the bath oscillators or to nonlinearity of the system-bath interactions, we have to include hierarchical elements with odd orders of noise correlation functions, most typically ⟨[x 2 j (t 1 ),xj(0)] ⟩ and ⟨[xj(t 2 ), [xj(t 1 ),xj(0)]] ⟩ for the bath coordinatexj for various time sequences ti. The contributions from even orders of noise correlation functions, for example, ⟨[xj(t 3 ), [xj(t 2 ), [xj(t 1 ),xj(0)]]]⟩, must be taken into account because the higher-order cumulants do not vanish for a non-Gaussian bath. Accordingly, we have to introduce an SDF for these higher-order noise correlation functions, as has been investigated in molecular dynamics (MD) simulations. 216 Complex profiles of these correlation functions have been investigated as the observables of multidimensional vibrational spectroscopy. 186,217 For electron transport in nanosystems, non-Gaussian dynamics were investigated using the HEOM approach in the framework of full counting statistics. 209 It should be noted that in a non-Gaussian case, either a fermionic or bosonic case, the regular fluctuation-dissipation theorem cannot be applied because the bath response is highly nonlinear. Nevertheless, assuming that the non-Gaussian effects are minor, we can still adopt the HEOM structure as the starting point to simulate the system dynamics. 218,219 Non-Gaussian effects can then be included as non-Gaussian corrections in the HEOM formalism. 220 As an approximate approach, the HEOM have also been utilized to construct a time-dependent kernel of an open quantum dynamics equation for a charge carrier transfer process with Holstein-Peierls interactions. 221

G. Imaginary-time HEOM
Although we can obtain the correct (entangled) thermal equilibrium state of the reduced system, as the steady-state solution of the HEOM, it is easier to solve the imaginary-time HEOM by integrating over the inverse temperature. 43,44 Because the structures of the imaginary-time HEOM and the wave-function-based HEOM are similar, the cost of numerical integration is much cheaper than that of the regular HEOM. Moreover, because the solution of the imaginary-time HEOM is the partition function rather than the distribution function, we can evaluate thermodynamic variables, most importantly the Helmholtz free energy and the Boltzmann entropy, as a function of temperature. 43,44 If we wish to obtain a steady-state solution under a periodic external force, however, we must use the real-time HEOM. Therefore, a method has been developed to obtain steady-state solutions of the HEOM very efficiently. 222,223 VI. APPLICATIONS A. Proton, electron, and exciton transfer problems Many chemical physics and biochemical physics problems involve environments that are complex and strongly coupled to a molecular system at finite temperature. Moreover, the non-Markovian effects of environments arising from intermolecular and intramolecular interactions play essential roles. Therefore, a great deal of effort has been dedicated to studying the problems of quantum as well as classical dynamics with nonperturbative and non-Markovian system-environment interactions from the HEOM approach. Commonly studied problems of this kind are chemical reactions, 89,90,122 proton transfer, 115,224 electron transfer, 138,139,144,225 electron-coupled proton transfer, 226 excitoncoupled electron transfer, 227 charge separation, 228 exciton-hole separation, 229 and nonadiabatic transitions in photochemical processes, 128,[130][131][132] including those involving conical intersection [230][231][232][233] and exciton transfer.  In particular, the HEOM approach has been used to simulate the exciton transfer process of the photosynthesis antenna system with the aim of investigating how natural systems can realize such highly efficient yields, presumably by manipulating quantum mechanical processes. This is because the  97 and because it allows the calculation of multidimensional electronic spectra, which are the observables of experiments of this kind. 179 The numerical costs of solving the HEOM for photosynthesis problems are high because each exciton or electron site is independently coupled with its own heat bath. For exciton transfer problems, on the other hand, because the systems are in the high-temperature regime, where the quantum nature of thermal noise plays a minor role, we can integrate the HEOM by ignoring the low-temperature correction terms. Thus, although the Fenna-Matthews-Olsen (FMO) system has often been used as a benchmark test, it is not suitable to verify the description of open quantum dynamics (however, it is convenient for testing the scalability of calculations, in particular for a realistic SDF). As a test of open quantum dynamics, the non-Markovian tests presented in Sec. IV are more suitable.
The ability of the HEOM to deal with time-dependent external forces is ideal for investigating the possibility of Floquet engineering for exciton transfer processes 261 and electron transfer 262 in which the system Hamiltonian involves a periodic external force. On the basis of the HEOM approach, the effects of nondissipative local phonon modes have been investigated. Examples of this kind are systems described by the Holstein model 210,211 and the Holstein-Tavis-Cummings (HTC) model. 263

B. Nonlinear and multidimensional spectroscopies
Although it is only relatively recently that the sensitivity of a nonlinear response function has been investigated as a nonlinear quantum measure, 93-95 the importance of non-Markovian effects as well as that of the bathentanglement effects had been realized since the 1980s, when the ultrafast dynamics of molecular motion were first measured by ultrafast nonlinear laser spectroscopy. 21 Thus, nonlinear spectra, in particular multidimensional spectra, can be used as a measure of bathentanglement, which characterizes the difference betweenρtot(t) andρA(t) ⊗ρ eq B . The theoretical foundation for the computation of nonlinear spectra is briefly explained in the Appendix. Because the HEOM approach was originally developed to simulate nonlinear spectra, 120,121,130,131 it is able to calculate nonlinear response functions, which the conventional reduced equation of motion approaches cannot do, as illustrated in Figs. 8(d) and 15.
At present, nonlinear spectroscopic measurements are the only experimental approach for observing the ultrafast dynamics of molecular motions, including exciton and electron transfer, and so it is important to calculate nonlinear spectra to confirm the descriptions provided by models. 241 Because the HEOM approach does not have to employ the eigenstate representation of the system, it can treat any profile of time-dependent external fields for which the eigenstate of the system cannot be defined. Thus, optical Stark spectroscopy, 130,131,137 doubleslit experiments, 265 and optimal control problems 266 have been investigated simply by integrating the HEOM with time-dependent perturbations.
As an extension of the stochastic approach, the HEOM has been used to calculate linear and nonlinear NMR spectra 42,180,267 and muon spin (μSR) spectra. 42,96 C. Quantum transport problems At present, there are two HEOM approaches for the treatment of quantum transport problems. One is based on the Wigner distribution function, where the particles can come in or come out through the boundary condition of the system. The heat bath then acts as the thermal energy source for the system. The other is based on the grand canonical treatment of the heat baths in which the baths act as infinite sources or absorbers of particles and the flow of the particles among the baths is controlled by the chemical potentials.
The Wigner-function-based approach has been applied to the study of chemical reactions, 90 quantum ratchets, 108 molecular motors, 132 (see Fig. 6), photodissociation, 130,131 and resonant tunneling processes. [105][106][107] Calculations have been performed using the periodic, open, and inflow-outflow conditions. The investigation of quantum ratchets showed that the current is suppressed as the tunneling rate increases. 108 In the investigation of the resonant tunneling diode, a self-current oscillation was discovered in the negativeresistance region, although the Hamiltonian is time-independent (see Fig. 4).

D. Quantum information and quantum thermodynamics
Because the HEOM can provide an exact numerical treatment of the dynamics defined by the Hamiltonian, it is possible to carry out desktop experiments to verify fundamental statements about quantum information and quantum thermodynamics in a practical manner under extreme quantum conditions. Commonly studied problems of quantum thermodynamics are quantum heat flow 100,156,[289][290][291] and quantum heat engines. 101,102 Because the quantum systems we consider are on the nanoscale, we cannot ignore the quantum contributions from the system-bath interaction. Thus, for example, it has been shown that the second law of thermodynamics is broken if we define the heat current using the system energy instead of the bath energy. 100 As a problem of quantum information, quantum entanglement is important not only for the system itself but also for the systembath interaction. A frequently used quantum measure for a twoqubit system is the concurrence. From the HEOM approach, we observe a rapid revival of the concurrence after the sudden death of entanglement because the quantum entanglement between the two qubits builds up through the bathentanglement between the spins and the baths (see Fig. 13). The sensitivity of a nonlinear response function has also been investigated as a nonlinear quantum measure. 93,94 The effects of geometrical phase 292 and quantum synchronization 293 have been investigated simply by integrating the HEOM with time-dependent driving fields.

VII. FUTURE PERSPECTIVES
Thirty years have passed since HEOM theory was developed as a bridge between the perturbative and Markovian quantum master equation theory and phenomenological stochastic theory. With this in mind, we summarize here the limitations and major challenges of the currently available HEOM theories.
The advantage of the HEOM approach lies in the structure of the equations of motion: a set of simultaneous differential equations for the elements of the hierarchy can utilize the bathentanglement even under a strong time-dependent external force. When the external perturbation is switched off, the system approaches the correct thermal equilibrium state at finite temperature, which is important for the study of relaxation dynamics. The HEOM approach allows us to calculate nonlinear response functions, which is significant for the detection of dynamical bathentanglement. Tremendous efforts have been devoted by many researchers to devise extensions of the HEOM approach that will aid in the development of theoretical background, numerical techniques, and applications and thereby further the study of open quantum dynamics.
Thus, we can now study reasonably large systems, including those consisting of 20-30 spins and those described by an anharmonic two-dimensional potential, in a numerically "exact" manner under quantum mechanically extreme conditions. While the capability of this approach is still limited, it should become possible to utilize such extended degrees of freedom to provide descriptions of realistic systems and of realistic anharmonic heat baths.
For example, on the system side, we should be able to employ a quantum chemistry approach or a first-principles MD approach to describe the system dynamics, while the effects of the molecular environment are modeled using a harmonic oscillator bath. Machine learning approaches will be helpful for the construction of realistic system-bath models. 294 On the bath side, to treat non-Gaussian heat baths and spin or correlated fermionic baths, a practical approach is to introduce a reasonably large subsystem between the system and the bath that consists of an ensemble of anharmonic oscillators, spins, or electrons that can interact with each other. If the subsystem degrees of freedom are sufficiently large and the interactions between the subsystem and bath are weak, we should be able to study the effects of an anharmonic or non-bosonic bath whose thermal states are not characterized by the fluctuation-dissipation theorem. Using a model of this kind, we can investigate, for example, the quantum dynamics of an impurity in a complex molecular matrix that is, in turn, coupled to a phonon bath. The numerical cost of facilitating the treatment of large subsystems is extremely high, and therefore, efficient algorithms have to be employed.
The key feature of the HEOM formalism arises from the definition of the hierarchical elements, where the zeroth hierarchical element includes all orders of system-bath interactions and is the exact solution of the total Hamiltonian, with the higher members then including the lower-order system-bath interactions, in contrast to conventional perturbative approaches. Although the structure of the HEOM becomes complex and may not be truncated in a simple manner, there is no inherent restriction on constructing the HEOM for any system that interacts with an environment. To study quantum coherence among spatially distributed sites comparable to the thermal de Broglie wavelength, it would be important to construct the HEOM for the Hamiltonian in momentum space, for example for the Fröhlich Hamiltonian, in the way derived for the deformation potential Hamiltonian. 87 In this way, we can explore the interplay of quantum coherence in time and space in a uniform manner.
Finally, it should be mentioned that the framework of the HEOM theory can also be applied to the relativistic quantum theory of the electromagnetic field and to other field theories in which the fundamental interactions arise from the exchange of gauge bosons. The nonperturbative nature of the HEOM approach may provide new insight into the problem of irreversibility not only in time but also in space that we experience in our universe.

APPENDIX: NONLINEAR RESPONSE FUNCTIONS
In quantum mechanics, any physical observable can be expressed as an expectation value of a physical operator. In an optical measurement, the observable at time t is expressed as P(t) = tr{μρtot(t)} or P(t) = tr{Πρtot(t)}, whereρtot(t) involves the interaction between the excitation fields and the system through the dipole operatorμ or the polarizability operatorΠ. 42 We expand ρtot(t) in terms of the nth-order response function, which involves n excitations. In laser spectroscopy, any order of the response function is thus expressed by 2 n elements with different configurations corresponding to different time evolutions of the density matrix element. 21 The first-order contribution is the linear response function, R (1) (t 1 ) = ⟨[μ(t 1 ),μ] ⟩/ ̵ h. The signal for the harmonic case is presented in Fig. 8(c). Using the Liouville notation, this is expressed as R (1) (t 1 ) = −i tr{μĜ(t 1 )μ ×ρeq }/ ̵ h, whereĜ(t) is the propagator, which does not involve the laser interactions. 120 Each time evolution for the second-and third-order cases is characterized by the quantum Liouville paths for n = 2 and 3, as illustrated in Fig. 14. The second-order contribution is the observable of 2D Raman, 2D terahertz Raman, and 2D IR Raman spectroscopy. For 2D terahertz Raman spectroscopy, the second-order response function is expressed as FIG. 14. Double-sided Feynman diagrams for the second-and third-order response functions: (a) elements of R (2) (t 2 , t 1 ) and (b) elements of R (3) (t 3 , t 2 , t 1 ). In each diagram, time runs from bottom to top, and t i represents the time intervals for the ith sequence between the successive laser-system interactions. The left line represents the time evolution of the ket, whereas the right line represents that of the bra. The excited states are represented by the red lines. The complex conjugate paths of these, which can be obtained by interchanging the ket and bra diagrams, are not shown here (see Ref. 21).
In the HEOM approach, the density matrix is replaced by a reduced one, and the Liouvillian in G(t) is replaced using Eqs. (5) or (8). We then evaluate the response function, for example, the RII (t 2 , t 1 ) by the following steps: (i) The system is initially in the equilibrium state, which is obtained by numerically integrating the HEOM until a steady-state is reached. (ii) The system is excited by the first Raman interaction at t = 0 by applyingΠ × to all of the hierarchical elements to take into account the system-bath entangled states in Figs. 2(a)-2(c) that overlap with the diagrams presented in Fig. 14(a), as illustrated in Fig. 7(a-i). (iii) The time evolution of the perturbed elements is then computed by integrating the HEOM for the time period t 1 . (iv) The second excitationμ × is then applied, and the perturbed elements are again computed by integrating the HEOM for the time period t 2 .
The response function R RII (t 2 , t 1 ) is calculated from the expectation value ofμ. If necessary, we calculate R The higher-order response functions can be calculated in a similar manner. If we wish to calculate the response function for a specific Liouville path presented in Fig. 14, we selectively applyμ from the right or left instead ofμ × by adjusting the method outlined in steps (i)-(iv) to the targeting Liouville path. 144 In Fig. 15, we present the 2DES, which is obtained from the double Fourier transform of R (3) (t 3 , t 2 , t 1 ) for t 1 and t 3 for the diagrams in Figs. 14(b-i) and 14(b-iv) calculated from the HEOM approach and the TCL Redfield approach. In the high-temperature case, we can also evaluate the nonlinear response function from the HEOM approach analytically as the products of the resolvent G[s], which is the Laplace transform of G(t). [120][121][122] The sensitivity of a nonlinear response function has been investigated as a nonlinear quantum measure. [93][94][95] DATA AVAILABILITY Data sharing is not applicable to this article as no new data were created or analyzed in this study.