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.

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–10 electron and charge transfer rates,11–13 nonadiabatic transitions,14,15 quantum device systems,16 ratchet rectification,17,18 and superconducting quantum interference device (SQUID) rings19,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–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 spectroscopy35 and NMR spectroscopy36 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 equation2 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–41 the thermal equilibrium state and the dynamics of the original total Hamiltonian are altered.42–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 Redfield47,48 and time-convolution-less (TCL) Redfield equations49–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 approach52,53 and the reaction coordinate mapping approach54–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) approach58–60 and the time-dependent Davydov ansatz (TDDA),61–63 and asymptotic approaches, such as the effective-mode (EM) approach,64,65 the density matrix renormalization group (DMRG),66,67 and the time-evolving 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–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 wave-function-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–79 and to extend the applicability of the path-integral-based methods. Because these approaches can easily incorporate a semiclassical approximation for the bath80 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–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.

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–44,86–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–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.

We consider a situation in which a primary system interacts with heat baths that give rise to dissipation and fluctuation in the system.25–33 The total Hamiltonian is expressed as

Ĥtot=ĤA+ĤI+B,
(1)

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,

ĤI+B=aj=1Na(p^ja)22mja+12mja(ωja)2(x^ja)2αjaV^ax^ja,
(2)

where the momentum, position, mass, and frequency of the jth oscillator in the ath bath are given by p^ja, x^ja, mja, and ωja, respectively, V^a is the system part of the interaction, and αja 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ℏωj| j⟩⟨ j| + jkΔjk(t)| j⟩⟨k| and V^a=j,kVjka|jk|, 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–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^) and V^a=V^a(p^,q^) described in terms of the multidimensional momentum p^ and coordinate q^.25–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 aj(αjaVa)2/2mja(ωja)2 is introduced in the bath Hamiltonian such that the replacement x^jax¯^jax^jaαjaV^a/mja(ωja)2 is made25,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|jUjk(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

Ja(ω)j=1Na(αja)22mjaωjaδ(ωωja),
(3)

and the inverse temperature β ≡ 1/kBT, where kB 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αjx^j 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 L1(t)=iΩ^(t)Ω^Ω^Ω^(t)B/, while those of thermal fluctuations are expressed as L2(t)=Ω^(t)Ω^+Ω^Ω^(t)B/2, where ⟨⋯⟩B represents the thermal average of the bath degrees of freedom.31,42–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

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 as42–44,86–92

J(ω)=ηπγ2ωγ2+ω2,
(4)

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 system–bath coupling strength, which represents the magnitude of fluctuation and dissipation. By taking γ, the Drude form reduces to the Ohmic form J(ω) = ℏηω/π, which is often employed to obtain Markovian noise in the high-temperature case.

The relaxation function is temperature-independent and is expressed as L1(t)=c¯0eγ|t|, with c¯0=ηγ2, while the noise correlation is given by L2(t)=c0eγ|t|+kckeνk|t|, with c0 = ℏηγ2 cot(βℏγ/2)/2 and ck=2ηγ2νk/β(γ2νk2) for the kth Matsubara frequency νk ≡ 2πk/βℏ.42,86–90

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–41 To illustrate this point, we plot L2(t) in Drude cases for large γ in Fig. 1.42,44 The profile of the response function L1(t) is similar to that of L2(t) in the high-temperature case and becomes a δ function [∝ δ(t)] as γ. At low temperatures, however, the noise correlation L2(t) becomes negative owing to the contribution from terms with ckeνkt 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–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 tending to zero, L2(t) is always positive.

FIG. 1.

Noise correlation (fluctuation) function L2(t) for the Drude spectral distribution, Eq. (4), depicted as a function of the dimensionless time t in the fast-modulation case (γ = 10) for the different inverse temperatures βℏ. The profile of L2(t) in the high-temperature case becomes similar to that of the response function L1(t). While L2(t) exhibits a Markovian nature in the Ohmic case (γ) at high temperatures, it becomes negative at low temperatures (large βℏ) owing to the contribution of the Matsubara frequency terms. Because of this feature, the quantum noise cannot be Markovian. Reproduced with permission from Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). Copyright 2006 Author(s), licensed under a Creative Commons Attribution (CC BY) 4.0 Unported license.

FIG. 1.

Noise correlation (fluctuation) function L2(t) for the Drude spectral distribution, Eq. (4), depicted as a function of the dimensionless time t in the fast-modulation case (γ = 10) for the different inverse temperatures βℏ. The profile of L2(t) in the high-temperature case becomes similar to that of the response function L1(t). While L2(t) exhibits a Markovian nature in the Ohmic case (γ) at high temperatures, it becomes negative at low temperatures (large βℏ) owing to the contribution of the Matsubara frequency terms. Because of this feature, the quantum noise cannot be Markovian. Reproduced with permission from Y. Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). Copyright 2006 Author(s), licensed under a Creative Commons Attribution (CC BY) 4.0 Unported license.

Close modal

As can be seen from the form of L1(t) and L2(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., L1(t)] and fluctuations [i.e., L2(t)] expressed in terms of eγt, and the other is of quantum thermal origin and is characterized by fluctuations only [i.e., L2(t)], with the functions expressed in terms of ckeνkt for k ≥ 1, as the definition of νk = 2πk/βℏ 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 

The HEOM are derived by differentiating the reduced density matrix ρ^A(t), defined by a path integral42–44,86–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–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 ρ^j1,,jK(n), where the set of indices n and {jk} arise from the hierarchical expansion of the noise correlation functions in terms of eγt and eνkt, respectively.42–44 The zeroth element is identical to the reduced density matrix, ρ^A(t)=ρ^0,0,,0(0)(t), which includes all orders of the system–bath interactions, and it is the exact solution of the Hamiltonian, Eq. (1). For the set of n and {jk} with N=ω0/min(γ,ν1)>(n+Σk=1Kjk), where ω0 is the characteristic frequency of the system and K is a cutoff satisfying Kω0/ν1, the HEOM are expressed as42,43,86–89

tρ^j1,,jK(n)(t)=iL^A+nγ+k=1Kjkνk+Ξ^ρ^j1,,jK(n)(t)+Φ^ρ^j1,,jK(n+1)(t)+k=1Kρ^j1,,jk+1,,jK(n)(t)+nΘ^0ρ^j1,,jK(n1)(t)+k=1KjkΘ^kρ^j1,,jk1,,jK(n)(t),
(5)

where L^A is the Liouvillian of ĤA, defined as L^AĤA×/,

Ξ^=k=K+1ckνkV^×V^×2

is the renormalization operator to compensate for the effects of the cutoff K, and

Φ^=iV^×,  Θ^0=i(c0V^×ic¯0V^),  Θ^k=ickV^×

for an integer k ≥ 1. Here, the hyper-operators O^×f^O^f^f^O^ and O^f^O^f^+f^O^ for any operator Ô and operand f^ are those introduced by Kubo98 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 ρ^j1,,jK(n+1)(t) and ρ^j1,,jk+1,,jK(n)(t) to ρ^j1,,jK(n)(t), whereas Θ^0 and Θ^k create the bath states from ρ^j1,,jK(n1)(t) and ρ^j1,,jk1,,jK(n)(t) to ρ^j1,,jK(n)(t), respectively; ρ^j1,,jK(n1)(t) and ρ^j1,,jk1,,jK(n)(t) involve the first-order higher (−) and lower (+) interactions than ρ^j1,,jK(n)(t). Thus, one can regard Eq. (5) as a type of rate law among the density operators with three different bath-excited states, ρ^j1,,jK(n±1)(t) and ρ^j1,,jk±,,jK(n)(t); the time evolution of ρ^j1,,jK(n)(t) is determined by its time evolution operator (L^A+nγ+k=1Kjkνk+Ξ^) and the incoming and outgoing contributions of the hierarchical elements connected through Φ^, Θ^0, and Θ^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 {jk}. 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–102 

FIG. 2.

Double-sided Feynman diagrams of density matrices with the three lowest system–bath interactions. (a-i)–(a-iii) depict the elements of ρ^0,0,,0(0)(t) in terms of the second-, fourth-, and sixth-order system–bath interactions, whereas (b-i′)–(b-iii′) and (c-ii′′)–(c-iii′′) depict the elements of the first and second members of the hierarchy, which involve ρ^0,0,,0(1)(t) and ρ^0,0,,0(2)(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.

FIG. 2.

Double-sided Feynman diagrams of density matrices with the three lowest system–bath interactions. (a-i)–(a-iii) depict the elements of ρ^0,0,,0(0)(t) in terms of the second-, fourth-, and sixth-order system–bath interactions, whereas (b-i′)–(b-iii′) and (c-ii′′)–(c-iii′′) depict the elements of the first and second members of the hierarchy, which involve ρ^0,0,,0(1)(t) and ρ^0,0,,0(2)(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.

Close modal

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 system–bath interactions, and it is the exact solution of the Hamiltonian, Eq. (1).42–44,86–90 As an example, a numerical solution of the HEOM for a three-level heat engine model is presented in Fig. 3.101,102

FIG. 3.

The three-level heat engine model considered here consists of three states, denoted by |0⟩, |1⟩, and |2⟩, coupled to two bosonic baths. The system is driven by a periodic external field with frequency Ω. This model acts as a quantum heat engine when population inversion between the two excited states |1⟩ and |2⟩ occurs. We illustrate the time dependences of the heat current for η1 = (a) 0.01, (b) 0.1, and (c) 1 with η2 = 0.001 defined from the system energy, Q̇S1 and Q̇S2, and the power, Ẇ, that arise from transitions under one cycle of the external force. The time delays observed for the heat currents imply that the transition |0⟩ → |1⟩ → |2⟩ is cyclic. This behavior can be regarded as a microscopic manifestation of a quantum heat engine. Details are presented in Refs. 101 and 102. Reproduced with permission from A. Kato and Y. Tanimura, J. Chem. Phys. 145, 224105 (2016). Copyright 2016 AIP Publishing LLC.

FIG. 3.

The three-level heat engine model considered here consists of three states, denoted by |0⟩, |1⟩, and |2⟩, coupled to two bosonic baths. The system is driven by a periodic external field with frequency Ω. This model acts as a quantum heat engine when population inversion between the two excited states |1⟩ and |2⟩ occurs. We illustrate the time dependences of the heat current for η1 = (a) 0.01, (b) 0.1, and (c) 1 with η2 = 0.001 defined from the system energy, Q̇S1 and Q̇S2, and the power, Ẇ, that arise from transitions under one cycle of the external force. The time delays observed for the heat currents imply that the transition |0⟩ → |1⟩ → |2⟩ is cyclic. This behavior can be regarded as a microscopic manifestation of a quantum heat engine. Details are presented in Refs. 101 and 102. Reproduced with permission from A. Kato and Y. Tanimura, J. Chem. Phys. 145, 224105 (2016). Copyright 2016 AIP Publishing LLC.

Close modal

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 ρj1,,jK(n)(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 as103,104

Wj1,,jK(n)(p,q;t)12πdxeipx/ρj1,,jK(n)qx2,q+x2;t.
(6)

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 form104 

L^QMWj1,,jK(n)(p,q)pmqWj1,,jK(n)(p,q)1dp2πUW(pp,q)Wj1,,jK(n)(p,q),
(7)

where

UW(p,q)=20dxsinpxUq+x2Uqx2.

The quantum hierarchical Fokker–Planck equations (QHFPE) with the counterterm can then be expressed as44,105–108

tWj1,,jK(n)(p,q;t)     =L^QM+nγ+k=1Kjkνk+Ξ^Wj1,,jK(n)(p,q;t)      +Φ^Wj1,,jK(n+1)(p,q;t)+k=1KWj1,,jk+1,,jK(n)(p,q;t)      +nΘ¯^0Wj1,,jK(n1)(p,q;t)      +k=1KjkΘ^kWj1,,jk1,,jK(n)(p,q;t),
(8)

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 as

Φ^=p,Θ¯^0=ηpm+c0cotβγ2p,Θ^k=ckp,Ξ^ηβk=K+1ck2p2.

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–107 The QHFPE for linear–linear (LL) plus square–linear (SL), V(q) = q + cq2, and exponential–linear (EL), V(q) = ecq, system–bath coupling cases for any constant c are given in Refs. 109–114 and Ref. 115, respectively. It should be noted that because the counterterm is incorporated into the equations of motion to maintain the translational symmetry of the reduced equations of motion, the QHFPE are not merely a Wigner transformation of the regular HEOM, Eq. (5).

FIG. 4.

A snapshot of the Wigner distribution function for a self-current oscillation of a resonant tunneling diode. Current flows into the system from the emitter side of the boundary (q = 0 nm). Then, part of the current is scattered by the emitter side of the barrier. The scattered electrons flow in a tornado-like manner to the central peak in the emitter basin owing to dissipation. The shaking motion of the effective potential periodically accelerates the component at the central peak to the tunneling state, and the current thus exhibits steady oscillations, although the Hamiltonian is time-independent. Details are presented in Refs. 105–107. Reproduced with permission from A. Sakurai and Y. Tanimura, J. Phys. Soc. Jpn. 82, 033707 (2013). Copyright 2013 Author(s), licensed under a Creative Commons Attribution 4.0 Unported license.

FIG. 4.

A snapshot of the Wigner distribution function for a self-current oscillation of a resonant tunneling diode. Current flows into the system from the emitter side of the boundary (q = 0 nm). Then, part of the current is scattered by the emitter side of the barrier. The scattered electrons flow in a tornado-like manner to the central peak in the emitter basin owing to dissipation. The shaking motion of the effective potential periodically accelerates the component at the central peak to the tunneling state, and the current thus exhibits steady oscillations, although the Hamiltonian is time-independent. Details are presented in Refs. 105–107. Reproduced with permission from A. Sakurai and Y. Tanimura, J. Phys. Soc. Jpn. 82, 033707 (2013). Copyright 2013 Author(s), licensed under a Creative Commons Attribution 4.0 Unported license.

Close modal

In the LL high-temperature case, we can ignore the low-temperature correction terms, and the above equations for W(n)(p,q)=W0,,0(n)(p,q) further reduce to89,90

tW(n)(p,q)=(L^QM+nγ)W(n)(p,q)+pW(n+1)(p,q)+nγζp+mβpW(n1)(p,q),
(9)

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 → 0, which leads to the replacement of the Liouvillian through L^QML^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 equation118,119 for non-Markovian noise.89,90 Because βℏγ = 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.

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 element ρ^(0)(t) with a factorized initial condition ρ^(0)(0) is expressed as ρ^(0)[s]=Ĝ0[s]ρ^(0)(0), where the nth resolvent is defined in recursive form as42,120

Ĝn[s]=n+1s+(n+1)γ+iL^A+Φ^Ĝn+1[s]Θ^0.
(10)

Nonlinear response functions can also be expressed as products of Ĝn[s] (see Fig. 5), which allows us to evaluate the bathentanglement analytically.120–123 

FIG. 5.

The ultrafast nonlinear laser spectrum is calculated from the nonlinear-response 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, t2, ω3) is calculated for different t2 (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 T2, we can evaluate the chemical reaction rate directly from the measurement (c). Reproduced with permission from A. Ishizaki and Y. Tanimura, J. Chem. Phys. 123, 014503 (2005). Copyright 2005 AIP Publishing LLC.

FIG. 5.

The ultrafast nonlinear laser spectrum is calculated from the nonlinear-response 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, t2, ω3) is calculated for different t2 (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 T2, we can evaluate the chemical reaction rate directly from the measurement (c). Reproduced with permission from A. Ishizaki and Y. Tanimura, J. Chem. Phys. 123, 014503 (2005). Copyright 2005 AIP Publishing LLC.

Close modal

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 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(1/γ) for ν1 > γ and KγK for ν1γ. Then, for the case k=1Kjk>K, we truncate the hierarchical equations by replacing Eq. (8) with

tWj1,,jK(n)(p,q;t)=(L^QM+Ξ^)Wj1,,jK(n)(p,q;t),
(11)

while, for the case n = Kγ, we employ

tW0,,0(Kγ)(p,q;t)=(L^QM+KγγΦ^Θ¯^0+Ξ^)W0,,0(Kγ)(p,q;t)KγγΘ¯^0W0,,0(Kγ1)(p,q;t),
(12)

where

Ξ^mζβk=K+1ck2γ2γ2νk22p2

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 Wj1,,jK(n)(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 βℏγ ≪ 1, we have the QFPE116,117

tW(0)(p,q;t)=L^QMW(0)(p,q;t)+ζpp+mβpW(0)(p,q;t),
(13)

which is identical to the quantum master equation without the RWA.42 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. To overcome this limitation, we must include low-temperature correction terms to obtain the low-temperature-corrected QFPE (LT-QFPE).125 

Up to this point, we have separately discussed two types of dynamics employing distinct approaches: one is a discrete energy-level 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 as126–129 

W(p,q)=j,k|jWjk(p,q)k|,
(14)

where Wjk(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) ≡{Wjk(p, q)}. Similarly, the transition operator is written in Wigner representation form as

Xij±(p,q;t)=±idxeipx/Uijqx2;t.
(15)

The quantum Liouvillian can also be expressed in matrix form as126 

L^QMW(p,q;t)jk=pmqWjk(p,q;t)+1dp2π×mXjm+(pp,q;t)Wmk(p,q;t)+Xmk(pp,q;t)Wjm(p,q;t).
(16)

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) = {Vjk(q)}. In this way, we can obtain the multistate quantum hierarchical Fokker–Planck equations (MS-QHFPE) for Wj1,,jK(n)(p,q).128–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.

FIG. 6.

The dynamics of a molecular motor system and its nonlinear optical response (see the  Appendix) were investigated using the LT-MS-QSE approach. (a) The adiabatic Born–Oppenheimer potential energy surfaces (BOPESs) are plotted as a function of θ. The black and red curves represent the ground (S0) and excited (S1) BOPESs, respectively. The thermal stable state and photoproduct state are labeled “A” and “B.” Molecular structures of a typical two-step photodriven molecular rotary motor are also depicted. (b) Pump–probe spectrum (PPS) for the photoisomerization process from A to B. (c) Transient absorption spectrum (TAS) for the thermalization process from B to A. In both spectra, the red and blue areas represent emission and absorption, respectively. PPS is sensitive for the detection of the photoisomerization process, while the TAS is sensitive for the detection of the thermalization process. Reproduced with permission from T. Ikeda and Y. Tanimura, J. Chem. Phys. 147, 014102 (2017). Copyright 2017 AIP Publishing LLC.

FIG. 6.

The dynamics of a molecular motor system and its nonlinear optical response (see the  Appendix) were investigated using the LT-MS-QSE approach. (a) The adiabatic Born–Oppenheimer potential energy surfaces (BOPESs) are plotted as a function of θ. The black and red curves represent the ground (S0) and excited (S1) BOPESs, respectively. The thermal stable state and photoproduct state are labeled “A” and “B.” Molecular structures of a typical two-step photodriven molecular rotary motor are also depicted. (b) Pump–probe spectrum (PPS) for the photoisomerization process from A to B. (c) Transient absorption spectrum (TAS) for the thermalization process from B to A. In both spectra, the red and blue areas represent emission and absorption, respectively. PPS is sensitive for the detection of the photoisomerization process, while the TAS is sensitive for the detection of the thermalization process. Reproduced with permission from T. Ikeda and Y. Tanimura, J. Chem. Phys. 147, 014102 (2017). Copyright 2017 AIP Publishing LLC.

Close modal

A quantum electron transport problem described by the spinless Anderson–Holstein model was also investigated in a similar manner.133 

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 by ρ^Aeq=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).

FIG. 7.

Schematic diagrams involved in (a) the second-order and (b) the third-order response functions, which are defined as the expectation values of the system under the two and three impulsive excitations illustrated by the blue arrows, respectively (see the  Appendix). The accurate and HEOM descriptions of the bathentanglement (system–bath entangled) states are illustrated in (a-i) and (b-i), and the time-convolution-less (TCL) Redfield descriptions of these states are illustrated in (a-ii) and (b-ii). 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).

FIG. 7.

Schematic diagrams involved in (a) the second-order and (b) the third-order response functions, which are defined as the expectation values of the system under the two and three impulsive excitations illustrated by the blue arrows, respectively (see the  Appendix). The accurate and HEOM descriptions of the bathentanglement (system–bath entangled) states are illustrated in (a-i) and (b-i), and the time-convolution-less (TCL) Redfield descriptions of these states are illustrated in (a-ii) and (b-ii). 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).

Close modal

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 ± j), where cj, ζj, and ωj characterize the noise correlation for the jth hierarchical elements, to maintain the positivity of the HEOM formalism.137–148 Although the HEOM can be constructed for an arbitrary SDF in the form of a Fourier expansion87 or in terms of special functions,149–156 the HEOM that are derived from a finite set of non-time-decaying functions do not converge in time and thus may not describe phenomena toward the thermal equilibrium state.

The HEOM are simultaneous differential equations expressed in terms of the density matrix elements. Owing to the complex 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–159 or Fano160,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 sets163 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, wave-function-based HEOM approaches whose scalabilities are similar to that of the Schrödinger equation have been developed166–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-HEOM178–180 and DM-HEOM.181–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 

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ω02q^2/2, or to the system in the energy eigenstate representation, ĤA = ℏω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^+q^q^(t)/2, (c) the linear response function R(1)(t)i[q^(t),q^]  /, and (d) the nonlinear response function RTTR(2)(t2,t1)=[[q^2(t1+t2),q^(t1)],q^]  /2,186,187 which is the observable of 2D terahertz Raman spectroscopy188 (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 and the TCL Redfield equation.49–51 For cases (b)–(d), we plot the Fourier-transformed results denoted by C[ω], R(1)[ω], and RTTR(2)[ω1,ω2], respectively.

Figure 8(a) illustrates the factorized initial conditions (blue curves) and steady-state solutions (red curves) calculated from (i) the QHFPE and (ii) the TCL Redfield equation. After integrating the QHFPE and the TCL Redfield equation for a sufficiently long time, the distribution reaches a steady-state. In the case of the QHFPE, the obtained steady-state is identical within numerical error to the thermal equilibrium state obtained from the analytical expression, whereas those from the TCL Redfield equation are similar to 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.

FIG. 8.

Numerically “exact” tests for the description of (a) the bathentanglement thermal equilibrium state, (b) fluctuations, (c) dissipation, and (d) dynamical bathentanglement calculated for a Brownian oscillator system with frequency ω0 = 1.0.44 The values of the bath parameters for the Drude SDF were chosen as γ = 1.0, ζ = 1, and βℏ = 3.0, unless otherwise indicated. (a) Factorized initial conditions (blue curves) and steady-state solutions (red curves), calculated from (i) the QHFPE and (ii) the TCL Redfield equation. (b) Autocorrelation function of the Brownian oscillator system for βℏ = 1.0 and (c) the linear response functions for ζ = 3.0. Here, 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. (d) Second-order response function of the Brownian oscillator system, RTTR(2)[ω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(a). Reproduced with permission from Y. Tanimura, J. Chem. Phys. 142, 144110 (2015). Copyright 2015 Author(s), licensed under a Creative Commons Attribution 4.0 Unported license.

FIG. 8.

Numerically “exact” tests for the description of (a) the bathentanglement thermal equilibrium state, (b) fluctuations, (c) dissipation, and (d) dynamical bathentanglement calculated for a Brownian oscillator system with frequency ω0 = 1.0.44 The values of the bath parameters for the Drude SDF were chosen as γ = 1.0, ζ = 1, and βℏ = 3.0, unless otherwise indicated. (a) Factorized initial conditions (blue curves) and steady-state solutions (red curves), calculated from (i) the QHFPE and (ii) the TCL Redfield equation. (b) Autocorrelation function of the Brownian oscillator system for βℏ = 1.0 and (c) the linear response functions for ζ = 3.0. Here, 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. (d) Second-order response function of the Brownian oscillator system, RTTR(2)[ω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(a). Reproduced with permission from Y. Tanimura, J. Chem. Phys. 142, 144110 (2015). Copyright 2015 Author(s), licensed under a Creative Commons Attribution 4.0 Unported license.

Close modal

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 temperature-independent 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 RTTR(2)[ω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.

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–140 Ohmic,125 Lorentz,141 super-Ohmic,142 and Drude–Lorentz143 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 Chebyshev–quadrature spectral decompositions (eHEOM).149–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 

The difficulty of solving the HEOM arises from the fluctuation terms described by L2(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–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.

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–169 the stochastic Schrödinger equation (SSE),170 the hierarchy of stochastic Schrödinger equations (HSSE),171–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 → −iβℏtiβℏ → − iβℏ (see Fig. 9). Thus, whereas in the conventional HEOM approach, the HEOM density operators ρ^j1,,jK(n)(t+Δt) are computed from ρ^j1,,jK(n)(t), in the wave-function-based HEOM approach, we must repeat the full integration for t along −iβℏt + Δtiβℏ → − iβℏ 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 intermediate-time dynamics over a wider range of temperatures for a large system.

FIG. 9.

The complex counter path for the hierarchical Schrödinger equations of motion (HSEOM)177 depicted with the system–bath interaction for the Liouville path of the second-order response function presented in Fig. 7(a–i).

FIG. 9.

The complex counter path for the hierarchical Schrödinger equations of motion (HSEOM)177 depicted with the system–bath interaction for the Liouville path of the second-order response function presented in Fig. 7(a–i).

Close modal

Although the HEOM were originally obtained by considering a harmonic heat bath in a canonical ensemble, Jin and his collaborators195 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 states196 (see Fig. 10).

FIG. 10.

Atomistic (a) and schematic (b) models of a molecular junction, showing the left and right electrodes (leads) and the molecular bridge. Reproduced with permission from M. Thoss and F. Evers, J. Chem. Phys. 148, 030901 (2018). Copyright 2018 AIP Publishing LLC.

FIG. 10.

Atomistic (a) and schematic (b) models of a molecular junction, showing the left and right electrodes (leads) and the molecular bridge. Reproduced with permission from M. Thoss and F. Evers, J. Chem. Phys. 148, 030901 (2018). Copyright 2018 AIP Publishing LLC.

Close modal

We consider a Hamiltonian ĤA described in terms of the creation and annihilation operators of the jth electronic state d^j;σ and d^j;σ, where σ represents the spin state. The electrode bath is expressed as

ĤB=k;σĉk;σĉk;σ,
(17)

where ĉk;σ and ĉk;σ are the creation and annihilation operators of an electron in an electrode with wavevector k and spin σ′. The system–electrode interaction is expressed as

ĤI=f;σk;σVk;σĉk;σd^f;σ+ĉk;σd^f;σ,
(18)

where f represents the site of the system that is coupled to an electrode and Vk;σα 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,

Γα[ϵ]=2πj;σ(Vj;σα)2δ(ϵϵj;σ),
(19)

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,

Γα[ω]=2πζαγ(ωμα)2+γ2,
(20)

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

Γαs(t)=m=0nηmα,seγmα,st,
(21)

with

γ0α,s=γsiμα,η0α,s=πζαγ1+eiβγ,
γmα,s=β1(2m1)πsiμα,

and

ηmα,s=i2πζαγ2β[(sμαiγmα,s)2+γ2]

for m ≥ 1. The index s ∈ { +, −} corresponds to creation and annihilation of electrons, i.e., d^j;σ+=d^j;σ and d^j;σ=d^j;σ. Thus, the HEOM for a system coupled to the electrode baths can be expressed as196–198 

tρ^a1,,an(t)     =iL^A+k=1nγmkαk,skρ^a1,,an(t)      iad^j;σs¯ρ^a1,,an,a(t)()nρ^a1,,an,a(t)d^j;σs¯      ik=1n()nkηmkαk,skd^jk;σkskρ^a1,,ak1,ak+1,,an(t)      +()n(ηmkαk,s¯k)*ρ^a1,,ak1,ak+1,,an(t)d^jk;σksk,
(22)

where ak = {jk, αk, σk, mk, sk} and σk is the index of the system side of the spin in the jkth electronic state, mk is the index for the Matsubara frequencies, and sk can be + or − with s¯k=sk. 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 structure200,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 fluctuation–dissipation 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–209 

The HEOM approach has also been applied to the case of solid state materials described by a deformation potential87 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.

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^j2(t1),x^j(0)]   and [x^j(t2),[x^j(t1),x^j(0)]]   for the bath coordinate x^j for various time sequences ti. The contributions from even orders of noise correlation functions, for example, [x^j(t3),[x^j(t2),[x^j(t1),x^j(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 

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

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 exciton-coupled electron transfer,227 charge separation,228 exciton–hole separation,229 and nonadiabatic transitions in photochemical processes,128,130–132 including those involving conical intersection230–233 and exciton transfer.234–260 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 HEOM approach can be applied seamlessly across a wide coupling range from the weak Redfield regime to strong Förster exciton–environment coupling under various non-Markovian conditions,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 processes261 and electron transfer262 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 model210,211 and the Holstein–Tavis–Cummings (HTC) model.263 

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)ρ^Beq.

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.

Since the 1990s, 2D Raman and 2D terahertz Raman vibrational spectra for intermolecular modes109–113,187 and 2D IR spectra for high-frequency intramolecular modes114 have been calculated using the QHFPE with LL + SL system–bath interactions. As shown in Fig. 11, the MD result of 2D IR–Raman was consistently described using the classical hierarchical Fokker–Planck equations (CHFPE) approach and the BO model with anharmonic and nonlinear system–bath coupling.264 An example of a calculated 2D vibrational spectrum (2DVS) is also presented in Fig. 5. Because the HEOM formalism treats the quantum and classical systems with any form of potential from the same point of view, it allows the identification of purely quantum mechanical effects through comparison of classical and quantum results in the Wigner distribution.114 

FIG. 11.

2D Raman–IR–IR spectra obtained from molecular dynamics (MD) simulation and from the CHFPE on the basis of the multimode nonlinear Brownian oscillator (BO) model whose parameters were chosen to reproduce 2D IR–Raman spectra obtained from the MD simulation. (a) MD result for the contributions of intermolecular hydrogen-bond (HB) and intramolecular OH stretching mode–mode coupling. (b) CHFPE result. From the BO model analysis, we found that the existence of the negative peak A indicates strong coupling between the HB and OH stretching modes, while the peak B arises from the nonlinearity of the molecular dipole and polarizability among these two modes. Reproduced from H. Ito and Y. Tanimura, J. Chem. Phys. 144, 074201 (2016). Copyright 2016 AIP Publishing LLC.

FIG. 11.

2D Raman–IR–IR spectra obtained from molecular dynamics (MD) simulation and from the CHFPE on the basis of the multimode nonlinear Brownian oscillator (BO) model whose parameters were chosen to reproduce 2D IR–Raman spectra obtained from the MD simulation. (a) MD result for the contributions of intermolecular hydrogen-bond (HB) and intramolecular OH stretching mode–mode coupling. (b) CHFPE result. From the BO model analysis, we found that the existence of the negative peak A indicates strong coupling between the HB and OH stretching modes, while the peak B arises from the nonlinearity of the molecular dipole and polarizability among these two modes. Reproduced from H. Ito and Y. Tanimura, J. Chem. Phys. 144, 074201 (2016). Copyright 2016 AIP Publishing LLC.

Close modal

The HEOM have also been used to calculate nonlinear electronic spectra. By choosing the Liouville paths presented in the  Appendix for targeting measurements, the signals of 2D electronic spectra (2DES)143–145,179,241–248 were calculated using the regular HEOM, while 2D electronic and vibrational (2DEVS) spectra were calculated using the MS-QHFPE128,130,131 and LT-MS-QFPE.125,132 Examples of calculated nonlinear spectra for pump–probe spectroscopy and transient absorption spectroscopy are presented in Figs. 6(b) and 6(c).

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 double-slit experiments,265 and optimal control problems266 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 spectra42,180,267 and muon spin (μSR) spectra.42,96

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–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 negative-resistance region, although the Hamiltonian is time-independent (see Fig. 4).

The approach based on the grand canonical ensemble has been applied to inelastic quantum transport problems modeled by a molecular system coupled to electron baths.196,206,268–271 More specifically, Coulomb blockade,272 the Kondo effect,222,273–275 quantum dots,197,276 Aharonov–Bohm interferometers,163,277 Anderson localization,278,279 molecular junctions,198,280 and vibronic reaction dynamics at metal surfaces280 have been studied. By combining density functional theory with the HEOM (DFT-HEOM), the electrical structure of graphene nanoribbons281 and the local spin states of adsorbed and embedded organometallic molecules that exhibit Kondo phenomena [including d-CoPc/Au(111),282 FePc/Au(111),283 few-layer CoPc/Au(111),284 Au–Co(tpy-SH)2–Au composites,285,286 CoPc and FePc on Pb(111),287 and FeOEP on Pb(111)288] have been investigated (see Fig. 12). In this treatment, the system can dissipate or obtain energy only through the exchange of particles with the baths.

FIG. 12.

(a) Side view of the optimized structure of the triple-layer CoPc/Pb(111) composite calculated from DFT-HEOM. The spin density distribution is also shown, with the isosurface for a spin density of 0.01 (−0.01) Å−1 shaded in red (blue). The lower panels display the PDOS of Co d orbitals in the (b) first, (d) second, and (f) third layers. [(c) and (e)] Schematic diagrams of the Co 3d orbitals in the second and third layers, respectively. Reproduced with permission from Y. Wang, X. Zheng, and J. Yang, J. Chem. Phys. 145, 154301 (2016). Copyright 2016 AIP Publishing LLC.

FIG. 12.

(a) Side view of the optimized structure of the triple-layer CoPc/Pb(111) composite calculated from DFT-HEOM. The spin density distribution is also shown, with the isosurface for a spin density of 0.01 (−0.01) Å−1 shaded in red (blue). The lower panels display the PDOS of Co d orbitals in the (b) first, (d) second, and (f) third layers. [(c) and (e)] Schematic diagrams of the Co 3d orbitals in the second and third layers, respectively. Reproduced with permission from Y. Wang, X. Zheng, and J. Yang, J. Chem. Phys. 145, 154301 (2016). Copyright 2016 AIP Publishing LLC.

Close modal

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 flow100,156,289–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 system–bath interaction. A frequently used quantum measure for a two-qubit 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 phase292 and quantum synchronization293 have been investigated simply by integrating the HEOM with time-dependent driving fields.

FIG. 13.

Decay of quantum entanglement (concurrence) for two qubits coupled to independent baths as a function of time from a maximum entangled state. We present the results calculated from the HEOM with the bathentanglement initial condition (solid curve), the HEOM from factorized initial condition (blue dashed curve), and the Markovian Redfield equation with the RWA (red dashed-dotted curve). Reproduced with permission from A. G. Dijkstra and Y. Tanimura, Phys. Rev. Lett. 104, 250401 (2010). Copyright 2010 The American Physical Society.

FIG. 13.

Decay of quantum entanglement (concurrence) for two qubits coupled to independent baths as a function of time from a maximum entangled state. We present the results calculated from the HEOM with the bathentanglement initial condition (solid curve), the HEOM from factorized initial condition (blue dashed curve), and the Markovian Redfield equation with the RWA (red dashed-dotted curve). Reproduced with permission from A. G. Dijkstra and Y. Tanimura, Phys. Rev. Lett. 104, 250401 (2010). Copyright 2010 The American Physical Society.

Close modal

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.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

The author is grateful to all researchers who have contributed to the development of HEOM theory. Financial support from the Kyoto University Foundation is acknowledged.

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 2n 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)(t1)=[μ^(t1),μ^]  /. The signal for the harmonic case is presented in Fig. 8(c). Using the Liouville notation, this is expressed as R(1)(t1)=itr{μ^Ĝ(t1)μ^×ρ^eq}/, 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 RTTR(2)(t2,t1)=tr{Π^Ĝ(t2)μ^×Ĝ(t1)μ^×ρ^eq}/2.295 The diagrams for the second-order response function are presented in Fig. 14(a). The complex conjugate diagrams of these are not depicted. The signal for the harmonic case is presented in Fig. 8(d–i).

FIG. 14.

Double-sided Feynman diagrams for the second- and third-order response functions: (a) elements of R(2)(t2, t1) and (b) elements of R(3)(t3, t2, t1). In each diagram, time runs from bottom to top, and ti 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).

FIG. 14.

Double-sided Feynman diagrams for the second- and third-order response functions: (a) elements of R(2)(t2, t1) and (b) elements of R(3)(t3, t2, t1). In each diagram, time runs from bottom to top, and ti 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).

Close modal

The third-order response function describes 2D infrared (2DIR) spectroscopy as R(3)(t3,t2,t1)=itr{μ^Ĝ(t3)μ^×Ĝ(t2)μ^×Ĝ(t1)μ^×ρ^eq}/3, where μ^ represents the dipole operator.110–113 The same expression can describe third-order electronic measurements, which include two-dimensional electronic spectroscopy (2DES), by regarding μ^ as the electronic transition dipole operator.144 By using the third-order diagrams, the pump–probe spectrum (PPS) and transient absorption spectrum (TAS) are, for example, calculated from the diagrams presented in Figs. 14(b-ii) and 14(b-iii).

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 2D Raman–IR–IR response function RRII(2)(t2,t1) by the following steps:

  • The system is initially in the equilibrium state, which is obtained by numerically integrating the HEOM until a steady-state is reached.

  • 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).

  • The time evolution of the perturbed elements is then computed by integrating the HEOM for the time period t1.

  • The second excitation μ^× is then applied, and the perturbed elements are again computed by integrating the HEOM for the time period t2.

The response function RRII(2)(t2,t1) is calculated from the expectation value of μ^. If necessary, we calculate RRII(2)[ω1,ω2] by performing a fast Fourier transform for t1 and t2 [Figs. 8(iv) and 11]. 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) (t3, t2, t1) for t1 and t3 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–122 

FIG. 15.

Two-dimensional electronic spectrum of a two-level system calculated from (a) the HEOM approach and (b) the TCL Redfield approach for ω0 = 2100 cm−1, ζ = 0.5ω0, γ = 0.005ω0, and βℏω0 = 10. Because the TCL approach cannot properly take account of the system–bath coherence (entanglement) over the t2 period, as illustrated by the blue dashed lines in Fig. 7(b–i), it does not account for the dephasing of the echo signal along the ω1 = ω2 direction. Reproduced with permission from A. Ishizaki and Y. Tanimura, Chem. Phys. 347, 185 (2008). Copyright 2008 Elsevier.

FIG. 15.

Two-dimensional electronic spectrum of a two-level system calculated from (a) the HEOM approach and (b) the TCL Redfield approach for ω0 = 2100 cm−1, ζ = 0.5ω0, γ = 0.005ω0, and βℏω0 = 10. Because the TCL approach cannot properly take account of the system–bath coherence (entanglement) over the t2 period, as illustrated by the blue dashed lines in Fig. 7(b–i), it does not account for the dephasing of the echo signal along the ω1 = ω2 direction. Reproduced with permission from A. Ishizaki and Y. Tanimura, Chem. Phys. 347, 185 (2008). Copyright 2008 Elsevier.

Close modal

The sensitivity of a nonlinear response function has been investigated as a nonlinear quantum measure.93–95 

1.
R. K.
Wangsness
and
F.
Bloch
,
Phys. Rev.
89
,
728
(
1953
).
3.
4.
M. O.
Scully
and
W. E.
Lamb
,
Phys. Rev.
159
,
208
(
1967
).
5.
B. R.
Mollow
and
M. M.
Miller
,
Ann. Phys.
52
,
464
(
1969
).
7.
F.
Haake
,
Quantum Statistics in Optics and Solid-State Physics
, Springer Tracts in Modern Physics Vol. 66 (
Springer
,
Berlin
,
1973
), p.
98
.
8.
P. G.
Wolynes
,
Phys. Rev. Lett.
47
,
968
(
1981
).
9.
G. A.
Voth
,
D.
Chandler
, and
W. H.
Miller
,
J. Chem. Phys.
91
,
7749
(
1989
).
10.
P.
Hänggi
,
P.
Talkner
, and
M.
Borkovec
,
Rev. Mod. Phys.
62
,
251
(
1990
).
11.
R. A.
Marcus
,
Rev. Mod. Phys.
65
,
599
(
1993
).
12.
V.
May
and
O.
Kühn
,
Charge and Energy Transfer Dynamics in Molecular Systems
(
Wiley VCH
,
Berlin
,
2003
).
13.
M.
Schroeter
,
S. D.
Ivanov
,
J.
Schulze
,
S. P.
Polyutov
,
Y. J.
Yan
,
T.
Pullerits
, and
O.
Kühn
,
Phys. Rep.
567
,
1
(
2015
).
14.
A.
Garg
,
J. N.
Onuchic
, and
V.
Ambegaokar
,
J. Chem. Phys.
83
,
4491
(
1985
).
15.
Y. J.
Yan
,
M.
Sparpaglione
, and
S.
Mukamel
,
J. Phys. Chem.
92
,
4842
(
1988
).
16.
B. A.
Mason
and
K.
Hess
,
Phys. Rev. B
39
,
5051
(
1989
).
17.
P.
Reimann
,
M.
Grifoni
, and
P.
Hänggi
,
Phys. Rev. Lett.
79
,
10
(
1997
).
18.
P.
Hänggi
and
F.
Marchesoni
,
Rev. Mod. Phys.
81
,
387
(
2009
).
19.
Y.-C.
Chen
,
J. Low Temp. Phys.
65
,
133
(
1986
).
20.
K.
Mitra
,
F. W.
Strauch
,
C.
Lobb
,
J. R.
Anderson
, and
F. C.
Wellstood
,
Phys. Rev. B
77
,
214512
(
2008
).
21.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
New York
,
1995
).
22.
P.
Hamm
and
M. T.
Zanni
, in
Concepts and Methods of 2D Infrared Spectroscopy
(
Cambridge University Press
,
1995
).
23.
E. B.
Davies
,
Quantum Theory of Open Systems
(
Academic Press
,
1976
).
24.
25.
A. O.
Caldeira
and
A. J.
Leggett
,
Phys. Rev. Lett.
46
,
211
(
1981
).
26.
A. O.
Caldeira
and
A. J.
Leggett
,
Ann. Phys.
149
,
374
(
1983
).
27.
D.
Waxman
and
A. J.
Leggett
,
Phys. Rev. B
32
,
4450
(
1985
).
28.
H.
Grabert
,
U.
Weiss
, and
P.
Talkner
,
Z. Phys. B: Condens. Matter
55
,
87
(
1984
).
29.
H.
Grabert
,
P.
Schramm
, and
G.-L.
Ingold
,
Phys. Rep.
168
,
115
(
1988
).
30.
U.
Weiss
,
Quantum Dissipative Systems
, 4th ed. (
World Scientific
,
Singapore
,
2012
).
31.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
,
Rev. Mod. Phys.
59
,
1
(
1987
).
32.
R.
Kubo
,
M.
Toda
, and
N.
Hashitsume
,
Statistical Physics
(
Springer-Verlag
,
1985
), Vol. 2.
33.
H. P.
Breuer
and
F.
Petruccione
,
The Theory of Open Quantum Systems
(
Oxford University Press
,
New York
,
2002
).
34.
Thermodynamics in the Quantum Regime
, edited by
F.
Binder
,
L. A.
Correa
,
C.
Gogolin
,
J.
Anders
, and
G.
Adesso
(
Springer
,
2018
).
35.
P.
Meystre
and
M.
Sargent
 III
,
Elements of Quantum Optics
, 2nd ed. (
Springer-Verlag
,
1991
).
36.
H.
Günther
,
NMR Spectroscopy: Basic Principles, Concepts and Applications in Chemistry
, 3rd ed. (
Wiley VCH
,
2013
).
37.
38.
K. F. F.
Romero
,
P.
Talkner
, and
P.
Hanggi
,
Phys. Rev.
69
,
052109
(
2004
).
39.
A.
Frigerio
,
J. T.
Lewis
, and
J. V.
Pulè
,
Adv. Appl. Math.
2
,
456
(
1981
).
40.
A.
Frigerio
,
J. T.
Lewis
, and
J. V.
Pulè
,
J. Approx. Theory
45
,
310
(
1985
).
41.
V.
Reimer
,
M. R.
Wegewijs
,
K.
Nestmann
, and
M.
Pletyukhov
,
J. Chem. Phys.
151
,
044101
(
2019
).
42.
Y.
Tanimura
,
J. Phys. Soc. Jpn.
75
,
082001
(
2006
).
43.
Y.
Tanimura
,
J. Chem. Phys.
141
,
044114
(
2014
).
44.
Y.
Tanimura
,
J. Chem. Phys.
142
,
144110
(
2015
).
45.
E. N.
Zimanyi
and
R. J.
Silbey
,
Philos. Trans. R. Soc., A
370
,
3620
(
2012
).
46.
N. G.
Van Kampen
,
Stochastic Processes in Physics and Chemistry
(
North-Holland
,
Amsterdam
,
1981
).
47.
W. M.
Zhang
,
T.
Meier
,
V.
Chernyak
, and
S.
Mukamel
,
J. Chem. Phys.
108
,
7763
(
1998
).
48.
M.
Yang
and
G. R.
Fleming
,
Chem. Phys.
275
,
355
(
2002
).
49.
F.
Shibata
,
Y.
Takahashi
, and
N.
Hashitsume
,
J. Stat. Phys.
17
,
171
(
1977
).
50.
S.
Chaturvedi
and
F.
Shibata
,
Z. Phys. B: Condens. Matter Quanta
35
,
297
(
1979
).
51.
M.
Ban
,
S.
Kitajima
, and
F.
Shibata
,
Phys. Lett. A
374
,
2324
(
2010
).
52.
B. M.
Garraway
,
Phys. Rev. A
55
,
4636
(
1997
).
53.
L.
Mazzola
,
S.
Maniscalco
,
J.
Piilo
,
K.-A.
Suominen
, and
B. M.
Garraway
,
Phys. Rev. A
80
,
012104
(
2009
).
54.
J.
Iles-Smith
,
N.
Lambert
, and
A.
Nazir
,
Phys. Rev. A
90
,
032114
(
2014
).
55.
G.
Schaller
,
J.
Cerrillo
,
G.
Engelhardt
, and
P.
Strasberg
,
Phys. Rev. B
97
,
195104
(
2018
).
56.
P.
Strasberg
,
G.
Schaller
,
T. L.
Schmidt
, and
M.
Esposito
,
Phys. Rev. B
97
,
205405
(
2018
).
57.
N.
Lambert
,
S.
Ahmed
,
M.
Cirio
, and
F.
Nori
,
Nat. Commun.
10
,
3721
(
2019
).
58.
H.-D.
Meyer
,
U.
Manthe
, and
L. S.
Cederbaum
,
Chem. Phys. Lett.
165
,
73
(
1990
).
59.
U.
Manthe
,
H. D.
Meyer
, and
L. S.
Cederbaum
,
J. Chem. Phys.
97
,
3199
(
1992
).
60.
H.
Wang
and
M.
Thoss
,
J. Phys. Chem. A
111
,
10369
(
2007
).
61.
B.
Luo
,
J.
Ye
,
C.
Guan
, and
Y.
Zhao
,
Phys. Chem. Chem. Phys.
12
,
15073
(
2010
).
62.
N.
Zhou
,
L.
Chen
,
Z.
Huang
,
K.
Sun
,
Y.
Tanimura
, and
Y.
Zhao
,
J.Phys. Chem. A
120
,
1562
(
2016
).
63.
L.
Chen
,
M. F.
Gelin
, and
W.
Domcke
,
J. Chem. Phys.
150
,
024101
(
2019
).
64.
L. S.
Cederbaum
,
E.
Gindensperger
, and
I.
Burghardt
,
Phys. Rev. Lett.
94
,
113003
(
2005
).
65.
R.
Martinazzo
,
B.
Vacchini
,
K. H.
Hughes
, and
I.
Burghardt
,
J. Chem. Phys.
134
,
011101
(
2011
).
66.
67.
S. R.
White
and
A. E.
Feiguin
,
Phys. Rev. Lett.
93
,
076401
(
2004
).
68.
A. W.
Chin
,
Á.
Rivas
,
S. F.
Huelga
, and
M. B.
Plenio
,
J. Math. Phys.
51
,
092109
(
2010
).
69.
J.
Prior
,
A. W.
Chin
,
S. F.
Huelga
, and
M. B.
Plenio
,
Phys. Rev. Lett.
105
,
050404
(
2010
).
70.
R.
Egger
and
C. H.
Mak
,
Phys. Rev. B
50
,
15210
(
1994
).
71.
J.
Cao
,
L. W.
Ungar
, and
G. A.
Voth
,
J. Chem. Phys.
104
,
4189
(
1996
).
72.
N.
Makri
,
J. Math. Phys.
36
,
2430
(
1995
).
73.
N.
Makri
and
D. E.
Makarov
,
J. Chem. Phys.
102
,
4600
(
1995
).
74.
N.
Makri
and
D. E.
Makarov
,
J. Chem. Phys.
102
,
4611
(
1995
).
75.
M.
Thorwart
,
P.
Reimann
, and
P.
Hänggi
,
Phys. Rev. E
62
,
5808
(
2000
).
76.
E.
Bukhman
and
N.
Makri
,
J. Phys.Chem. A
111
,
11320
(
2007
).
77.
V.
Jadhao
and
N.
Makri
,
J. Chem. Phys.
129
,
161102
(
2008
).
78.
N.
Makri
,
J. Chem. Phys.
141
,
134117
(
2014
).
79.
D.
Segal
,
A. J.
Millis
, and
D. R.
Reichman
,
Phys. Rev. B
82
,
205323
(
2010
).
80.
N.
Makri
,
Int. J. Quan. Chem.
115
,
1209
(
2015
).
81.
N.
Makri
,
J. Chem. Phys.
148
,
101101
(
2018
).
82.
L.
Diósi
and
W. T.
Strunz
,
Phys. Lett. A
235
,
569
(
1997
).
83.
J. T.
Stockburger
and
C. H.
Mak
,
J. Chem. Phys.
110
,
4983
(
1999
).
84.
J. T.
Stockburger
and
H.
Grabert
,
Phys. Rev. Lett.
88
,
170407
(
2002
).
85.
J.
Shao
,
J. Chem. Phys.
120
,
5053
(
2004
).
86.
Y.
Tanimura
and
R.
Kubo
,
J. Phys. Soc. Jpn.
58
,
101
(
1989
).
88.
A.
Ishizaki
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
74
,
3131
(
2005
).
89.
Y.
Tanimura
and
P. G.
Wolynes
,
Phys. Rev. A
43
,
4131
(
1991
).
90.
Y.
Tanimura
and
P. G.
Wolynes
,
J. Chem. Phys.
96
,
8485
(
1992
).
91.
R.-X.
Xu
,
P.
Cui
,
X.-Q.
Li
,
Y.
Mo
, and
Y.
Yan
,
J. Chem. Phys.
122
,
041103
(
2005
).
92.
P.
Han
,
R.-X.
Xu
,
B.
Li
,
J.
Xu
,
P.
Cui
,
Y.
Mo
, and
Y. J.
Yan
,
J. Phys. Chem. B
110
,
11438
(
2006
).
93.
A. G.
Dijkstra
and
Y.
Tanimura
,
Phys. Rev. Lett.
104
,
250401
(
2010
).
94.
A. G.
Dijkstra
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
81
,
063301
(
2012
).
95.
A. G.
Dijkstra
and
Y.
Tanimura
,
Philos. Trans. R. Soc., A
370
,
3658
(
2012
).
96.
H.
Takahashi
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
89
,
064710
(
2020
).
97.
A.
Ishizaki
and
G. R.
Fleming
,
J. Chem. Phys.
130
,
234110
(
2009
).
98.
T.
Takagahara
,
E.
Hanamura
, and
R.
Kubo
,
J. Phys. Soc. Jpn.
43
,
802
(
1977
).
99.
L.
Zhu
,
H.
Liu
,
W.
Xie
, and
Q.
Shi
,
J. Chem. Phys.
137
,
194106
(
2012
).
100.
A.
Kato
and
Y.
Tanimura
,
J. Chem. Phys.
143
,
064107
(
2015
).
101.
A.
Kato
and
Y.
Tanimura
,
J. Chem. Phys.
145
,
224105
(
2016
).
102.
A.
Kato
and
Y.
Tanimura
, “
Hierarchical equations of motion approach to quantum thermodynamics
,” in
Thermodynamics in the Quantum Regime
, Fundamental Theories of Physics Vol. 195, edited by
F.
Binder
, et al
(
Springer
,
2018
), p.
579
.
103.
104.
W. R.
Frensley
,
Rev. Mod. Phys.
62
,
745
(
1990
).
105.
A.
Sakurai
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
82
,
033707
(
2013
).
106.
A.
Sakurai
and
Y.
Tanimura
,
New J. Phys.
16
,
015002
(
2014
).
107.
F.
Grossmann
,
J. Chem. Phys.
141
,
144305
(
2014
).
108.
A.
Kato
and
Y.
Tanimura
,
J. Phys. Chem. B
117
,
13132
(
2013
).
109.
Y.
Tanimura
and
A.
Ishizaki
,
Acc. Chem. Res.
42
,
1270
(
2009
).
110.
T.
Steffen
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
69
,
3115
(
2000
).
111.
Y.
Tanimura
and
T.
Steffen
,
J. Phys. Soc. Jpn.
69
,
4095
(
2000
).
112.
T.
Kato
and
Y.
Tanimura
,
J. Chem. Phys.
117
,
6221
(
2002
).
113.
T.
Kato
and
Y.
Tanimura
,
J. Chem. Phys.
120
,
260
(
2004
).
114.
A.
Sakurai
and
Y.
Tanimura
,
J. Phys. Chem. A
115
,
4009
(
2011
).
115.
J.
Zhang
,
R.
Borrelli
, and
Y.
Tanimura
,
J. Chem. Phys.
152
,
214114
(
2020
).
116.
A. O.
Caldeira
and
A. J.
Leggett
,
Physica A
121
,
587
(
1983
).
117.
L.-D.
Chang
and
D.
Waxman
,
J. Phys. C: Soild State Phys.
18
,
5873
(
1985
).
119.
H.
Risken
,
The Fokker-Planck Equation
, 2nd ed. (
Springer
,
Berlin
,
1989
).
120.
Y.
Tanimura
and
R.
Kubo
,
J. Phys. Soc. Jpn.
58
,
1199
(
1989
).
121.
Y.
Tanimura
,
T.
Suzuki
, and
R.
Kubo
,
J. Phys. Soc. Jpn.
58
,
1850
(
1989
).
122.
A.
Ishizaki
and
Y.
Tanimura
,
J. Chem. Phys.
123
,
014503
(
2005
).
123.
Y.
Tanimura
and
R.
Kubo
,
J. Phys. Soc. Jpn.
58
,
3001
(
1989
).
124.
H.-D.
Zhang
and
Y.
Yan
,
J. Chem. Phys.
143
,
214112
(
2015
).
125.
T.
Ikeda
and
Y.
Tanimura
,
J. Chem. Theory Comput.
15
,
2517
(
2019
).
126.
Y.
Tanimura
and
S.
Mukamel
,
J. Chem. Phys.
101
,
3049
(
1994
).
127.
V.
Chernyak
and
S.
Mukamel
,
J. Chem. Phys.
105
,
4565
(
1996
).
128.
T.
Ikeda
and
Y.
Tanimura
,
J. Chem. Phys.
147
,
014102
(
2017
).
129.
T.
Ikeda
and
Y.
Tanimura
,
Chem. Phys.
515
,
203
(
2018
).
130.
Y.
Tanimura
and
Y.
Maruyama
,
J. Chem. Phys.
107
,
1779
(
1997
).
131.
Y.
Maruyama
and
Y.
Tanimura
,
Chem. Phys. Lett.
292
,
28
(
1998
).
132.
T.
Ikeda
,
A. G.
Dijkstra
, and
Y.
Tanimura
,
J. Chem. Phys.
150
,
114103
(
2019
).
133.
W.
Dou
,
C.
Schinabeck
,
M.
Thoss
, and
J. E.
Subotnik
,
J. Chem. Phys.
148
,
102317
(
2018
).
134.
A.
Ishizaki
and
Y.
Tanimura
,
Chem. Phys.
347
,
185
(
2008
).
135.
L.
Song
and
Q.
Shi
,
J. Chem. Phys.
143
,
194106
(
2015
).
136.
B.
Witt
,
Ł.
Rudnicki
,
Y.
Tanimura
, and
F.
Mintert
,
New J. Phys.
19
,
013007
(
2017
).
137.
Y.
Tanimura
and
S.
Mukamel
,
J. Phys. Soc. Jpn.
63
,
66
(
1994
).
138.
M.
Tanaka
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
78
,
073802
(
2009
).
139.
M.
Tanaka
and
Y.
Tanimura
,
J. Chem. Phys.
132
,
214502
(
2010
).
140.
J.-J.
Ding
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
136
,
224103
(
2012
).
141.
J.
Ma
,
Z.
Sun
,
X.
Wang
, and
F.
Nori
,
Phys. Rev. A
85
,
062323
(
2012
).
142.
H.
Liu
,
L.
Zhu
,
S.
Bai
, and
Q.
Shi
,
J. Chem. Phys.
140
,
134106
(
2014
).
143.
C.
Kreisbeck
and
T.
Kramer
,
J. Phys. Chem. Lett.
3
,
2828
(
2012
).
144.
Y.
Tanimura
,
J. Chem. Phys.
137
,
22A550
(
2012
).
145.
C.
Kreisbeck
,
T.
Kramer
, and
A.
Aspuru-Guzik
,
J. Phys. Chem. B
117
,
9380
(
2013
).
146.
J.-J.
Ding
,
J.
Xu
,
J.
Hu
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
135
,
164107
(
2011
).
147.
L.
Ye
,
H. D.
Zhang
,
Y.
Wang
,
X.
Zheng
, and
Y. J.
Yan
,
J. Chem. Phys.
147
,
074111
(
2017
).
148.
A.
Ishizaki
,
J. Phys. Soc. Jpn.
89
,
015001
(
2020
).
149.
H.
Tian
and
G.
Chen
,
J. Chem. Phys.
137
,
204114
(
2012
).
150.
Z.
Tan
,
Z.
Ouyang
,
Z.
Gong
,
H.
Wang
, and
J.
Wu
,
J. Chem. Phys.
143
,
224112
(
2015
).
151.
B.
Popescu
,
H.
Rahman
, and
U.
Kleinekathöfer
,
J. Phys. Chem. A
120
,
3270
(
2016
).
152.
H.
Rahman
and
U.
Kleinekathöfer
,
J. Chem. Phys.
150
,
244104
(
2019
).
153.
C.
Duan
,
Z.
Tang
,
J.
Cao
, and
J.
Wu
,
Phys. Rev. B
95
,
214308
(
2017
).
154.
C.
Duan
,
Q.
Wang
,
Z.
Tang
, and
J.
Wu
,
J. Chem. Phys.
147
,
164112
(
2017
).
155.
Q.
Wang
,
Z.
Gong
,
C.
Duan
,
Z.
Tang
, and
J.
Wu
,
J. Chem Phys.
150
,
084114
(
2019
).
156.
C.
Duan
,
C.-Y.
Hsieh
,
J.
Liu
, and
J.
Cao
,
J. Phys. Chem. Lett.
11
,
4080
(
2020
).
157.
J.
Hu
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
133
,
101106
(
2010
).
158.
B.-L.
Tian
,
J.-J.
Ding
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
133
,
114112
(
2010
).
159.
J.
Hu
,
M.
Luo
,
F.
Jiang
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
134
,
244106
(
2011
).
160.
L.
Cui
,
H.-D.
Zhang
,
X.
Zheng
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
151
,
024110
(
2019
).
161.
H.-D.
Zhang
,
L.
Cui
,
H.
Gong
,
R.-X.
Xu
,
X.
Zheng
, and
Y.
Yan
,
J. Chem. Phys.
152
,
064107
(
2020
).
162.
Q.
Shi
,
L.
Chen
,
G.
Nan
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
130
,
084105
(
2009
).
163.
D.
Hou
,
S.
Wang
,
R.
Wang
,
L.
Ye
,
R.
Xu
,
X.
Zheng
, and
Y.
Yan
,
J. Chem. Phys.
142
,
104112
(
2015
).
164.
Q.
Shi
,
Y.
Xu
,
Y.
Yan
, and
M.
Xu
,
J. Chem. Phys.
148
,
174102
(
2018
).
165.
R.
Borrelli
,
J. Chem. Phys.
150
,
234102
(
2019
).
166.
S.
Suess
,
A.
Eisfeld
, and
W. T.
Strunz
,
Phys. Rev. Lett.
113
,
150403
(
2014
).
167.
D.
Suess
,
W. T.
Strunz
, and
A.
Eisfeld
,
J. Stat. Phys.
159
,
1408
(
2015
).
168.
V.
Link
and
W. T.
Strunz
,
Phys. Rev. Lett.
119
,
180401
(
2017
).
169.
R.
Hartmann
and
W. T.
Strunz
,
J. Chem. Theory Comput.
13
,
5834
(
2017
).
170.
K.
Song
,
L.
Song
, and
Q.
Shi
,
J. Chem. Phys.
144
,
224105
(
2016
).
171.
Y.
Ke
and
Y.
Zhao
,
J. Chem. Phys.
145
,
024101
(
2016
).
172.
Y.
Ke
and
Y.
Zhao
,
J. Chem. Phys.
146
,
174105
(
2017
).
173.
Y.
Ke
and
Y.
Zhao
,
J. Chem. Phys.
146
,
214105
(
2017
).
174.
Y.
Ke
and
Y.
Zhao
,
J. Chem. Phys.
149
,
014104
(
2018
).
175.
M.
Lian
,
Y.-C.
Wang
,
Y.
Ke
, and
Y.
Zhao
,
J. Chem. Phys.
151
,
044115
(
2019
).
176.
Y.-C.
Wang
,
Y.
Ke
, and
Y.
Zhao
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
9
,
e1375
(
2019
).
177.
K.
Nakamura
and
Y.
Tanimura
,
Phys. Rev. A
98
,
012109
(
2018
).
178.
C.
Kreisbeck
,
T.
Kramer
,
M.
Rodríguez
, and
B.
Hein
,
J. Chem. Theory Comput.
7
,
2166
(
2011
).
179.
B.
Hein
,
C.
Kreisbeck
,
T.
Kramer
, and
M.
Rodríguez
,
New J. Phys.
14
,
023018
(
2012
).
180.
M.
Tsuchimoto
and
Y.
Tanimura
,
J. Chem. Theory Comput.
11
,
3859
(
2015
).
181.
J.
Strümpfer
and
K.
Schulten
,
J. Chem. Theory Comput.
8
,
2808
(
2012
).
182.
C.
Kreisbeck
,
T.
Kramer
, and
A.
Aspuru-Guzik
,
J. Chem. Theory Comput.
10
,
4045
(
2014
).
183.
T.
Kramer
,
M.
Noack
,
A.
Reinefeld
,
M.
Rodríguez
, and
Y.
Zelinskyy
,
J. Comput. Chem.
39
,
1779
(
2018
).
184.
T.
Kramer
,
M.
Noack
,
J. R.
Reimers
,
A.
Reinefeld
,
M.
Rodríguez
, and
S.
Yin
,
Chem. Phys.
515
,
262
(
2018
).
185.
Y.-a.
Yan
,
Chin. J. Chem. Phys.
30
,
277
(
2017
).
186.
Y.
Tanimura
and
S.
Mukamel
,
J. Chem. Phys.
99
,
9496
(
1993
).
187.
T.
Ikeda
,
H.
Ito
, and
Y.
Tanimura
,
J. Chem. Phys.
142
,
212421
(
2015
).
188.
P.
Hamm
and
A.
Shalit
,
J. Chem. Phys.
146
,
130901
(
2017
).
189.
J.
Xu
,
R.-X.
Xu
, and
Y. J.
Yan
,
New J. Phys.
11
,
105037
(
2009
).
190.
T.
Ikeda
and
G. D.
Scholes
,
J. Chem. Phys.
152
,
204101
(
2020
).
191.
Y.-A.
Yan
,
F.
Yang
,
Y.
Liu
, and
J.
Shao
,
Chem. Phys. Lett.
395
,
216
(
2004
).
192.
J. M.
Moix
and
J.
Cao
,
J. Chem. Phys.
139
,
134106
(
2013
).
193.
C.-Y.
Hsieh
and
J.
Cao
,
J. Chem. Phys.
148
,
014103
(
2018
).
194.
C.-Y.
Hsieh
and
J.
Cao
,
J. Chem. Phys.
148
,
014104
(
2018
).
195.
J.
Jin
,
S.
Welack
,
J.
Luo
,
X.-Q.
Li
,
P.
Cui
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
126
,
134113
(
2007
).
196.
J.
Jin
,
X.
Zheng
, and
Y.
Yan
,
J. Chem. Phys.
128
,
234703
(
2008
).
197.
R.
Härtle
,
G.
Cohen
,
D. R.
Reichman
, and
A. J.
Millis
,
Phys. Rev. B
88
,
235426
(
2013
).
198.
M.
Thoss
and
F.
Evers
,
J. Chem. Phys.
148
,
030901
(
2018
).
199.
L.
Han
,
V.
Chernyak
,
Y.-A.
Yan
,
X.
Zheng
, and
Y. J.
Yan
,
Phys. Rev. Lett.
123
,
050601
(
2019
).
200.
A.
Erpenbeck
,
L.
Götzendörfer
,
C.
Schinabeck
, and
M.
Thoss
,
Eur. Phys. J.
227
,
1981
(
2019
).
201.
A.
Erpenbeck
,
C.
Hertlein
,
C.
Schinabeck
, and
M.
Thoss
,
J. Chem. Phys.
149
,
064106
(
2018
).
202.
L.
Han
,
H.-D.
Zhang
,
X.
Zheng
, and
Y.
Yan
,
J. Chem. Phys.
148
,
234108
(
2018
).
203.
L.
Ye
,
X.
Wang
,
D.
Hou
,
R.-X.
Xu
,
X.
Zheng
, and
Y.
Yan
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
6
,
608
(
2016
).
204.
X.
Zheng
,
G.
Chen
,
Y.
Mo
,
S.
Koo
,
H.
Tian
,
C.
Yam
, and
Y.
Yan
,
J. Chem. Phys.
133
,
114101
(
2010
).
205.
H.
Xie
,
F.
Jiang
,
H.
Tian
,
X.
Zheng
,
Y.
Kwok
,
S.
Chen
,
C.
Yam
,
Y.
Yan
, and
G.
Chen
,
J. Chem. Phys.
137
,
044113
(
2012
).
206.
F.
Jiang
,
J.
Jin
,
S.
Wang
, and
Y. J.
Yan
,
Phys. Rev. B
85
,
245427
(
2012
).
207.
C.
Schinabeck
,
A.
Erpenbeck
,
R.
Härtle
, and
M.
Thoss
,
Phys. Rev. B
94
,
201407
(
2016
).
208.
C.
Schinabeck
,
R.
Härtle
, and
M.
Thoss
,
Phys. Rev. B
97
,
235429
(
2018
).
209.
C.
Schinabeck
and
M.
Thoss
,
Phys. Rev. B
101
,
075422
(
2020
).
210.
L.
Chen
,
Y.
Zhao
, and
Y.
Tanimura
,
J. Phys. Chem. Lett.
6
,
3110
(
2015
).
211.
I. S.
Dunn
,
R.
Tempelaar
, and
D. R.
Reichman
,
J. Chem. Phys.
150
,
184109
(
2019
).
212.
L.
Chen
,
M. F.
Gelin
, and
W.
Domcke
,
J. Chem. Phys.
151
,
034101
(
2019
).
213.
Y.
Iwamoto
and
Y.
Tanimura
,
J. Chem. Phys.
151
,
044105
(
2019
).
214.
H.-D.
Zhang
and
Y.
Yan
,
J. Phys. Chem. A
120
(
19
),
3241
(
2016
).
215.
R.-x.
Xu
,
X.-c.
Tao
,
Y.
Wang
,
Y.
Liu
,
H.-d.
Zhang
, and
Y.
Yan
,
Chin. J. Chem. Phys.
31
,
608
(
2018
).
216.
M.
Shiga
and
S.
Okazaki
,
J. Chem. Phys.
111
,
5390
(
1999
).
217.
K.
Okumura
and
Y.
Tanimura
,
J. Chem. Phys.
106
,
1687
(
1997
).
218.
Y.
Yan
,
J. Chem. Phys.
141
,
054105
(
2014
).
219.
Y.
Wang
,
Z. J.
Pan
,
H.-D.
Zhang
, and
Y. J.
Yan
,
Chem. Phys.
515
,
94
(
2018
).
220.
M.
Sarovar
and
M. D.
Grace
,
Phys. Rev. Lett.
109
,
130401
(
2012
).
221.
Y.
Yan
,
M.
Xu
,
Y.
Liu
, and
Q.
Shi
,
J. Chem. Phys.
150
,
234101
(
2019
).
222.
X.
Zheng
,
J.
Luo
,
J.
Jin
, and
Y.
Yan
,
J. Chem. Phys.
130
,
124508
(
2009
).
223.
H.-D.
Zhang
,
Q.
Qiao
,
R.-X.
Xu
,
X.
Zheng
, and
Y.
Yan
,
J. Chem. Phys.
147
,
044105
(
2017
).
224.
Q.
Shi
,
L.
Zhu
, and
L.
Chen
,
J. Chem. Phys.
135
,
044505
(
2011
).
225.
Q.
Shi
,
L.
Chen
,
G.
Nan
,
R.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
130
,
164518
(
2009
).
226.
K.
Song
and
Q.
Shi
,
J. Chem. Phys.
146
,
184108
(
2017
).
227.
S.
Sakamoto
and
Y.
Tanimura
,
J. Phys. Chem. Lett.
8
,
5390
(
2017
).
228.
Y.
Yan
,
L.
Song
, and
Q.
Shi
,
J. Chem. Phys.
148
,
084109
(
2018
).
229.
A.
Kato
and
A.
Ishizaki
,
Phys. Rev. Lett.
121
,
026001
(
2018
).
230.
L.
Chen
,
M. F.
Gelin
,
V. Y.
Chernyak
,
W.
Domcke
, and
Y.
Zhao
,
Faraday Discuss.
194
,
61
(
2016
).
231.
H.-G.
Duan
and
M.
Thorwart
,
J. Phys. Chem. Lett.
7
,
382
(
2016
).
232.
D.-L.
Qi
,
H.-G.
Duan
,
Z.-R.
Sun
,
R. J. D.
Miller
, and
M.
Thorwart
,
J. Chem. Phys.
147
,
074101
(
2017
).
233.
E.
Mangaud
,
B.
Lasorne
,
O.
Atabek
, and
M.
Desouter-Lecomte
,
J. Chem. Phys.
151
,
244102
(
2019
).
234.
A.
Ishizaki
and
G. R.
Fleming
,
Proc. Natl. Acad. Sci. U. S. A.
106
,
17255
(
2009
).
235.
M.
Sarovar
,
A.
Ishizaki
,
G. R.
Fleming
, and
K. B.
Whaley
,
Nat. Phys.
6
,
462
(
2010
).
236.
A.
Ishizaki
,
T. R.
Calhoun
,
G. S.
Schlau-Cohen
, and
G. R.
Fleming
,
Phys. Chem. Chem. Phys.
12
,
7319
(
2010
).
237.
P.
Nalbach
,
A.
Ishizaki
,
G. R.
Fleming
, and
M.
Thorwart
,
New J. Phys.
13
,
063040
(
2011
).
238.
J.
Strümpfer
and
K.
Schulten
,
J. Chem. Phys.
131
,
225101
(
2009
).
239.
J.
Strümpfer
and
K.
Schulten
,
J. Chem. Phys.
134
,
095102
(
2011
).
240.
A. G.
Dijkstra
and
Y.
Tanimura
,
New J. Phys.
14
,
073027
(
2012
).
241.
A. G.
Dijkstra
and
Y.
Tanimura
,
J. Chem. Phys.
142
,
212423
(
2015
).
242.
A. G.
Dijkstra
and
Y.
Tanimura
,
New J. Phys.
12
,
055005
(
2010
).
243.
D.
Green
,
B. S.
Humphries
,
A. G.
Dijkstra
, and
G. A.
Jones
,
J. Chem. Phys.
151
,
174112
(
2019
).
244.
J.
Xu
,
R.-x.
Xu
,
D.
Abramavicius
,
H.-d.
Zhang
, and
Y.-j.
Yan
,
Chin. J. Chem. Phys.
24
,
497
(
2011
).
245.
L.
Chen
,
R.
Zheng
,
Q.
Shi
, and
Y.
Yan
,
J. Chem. Phys.
132
,
024505
(
2010
).
246.
J.
Xu
,
H.-D.
Zhang
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
138
,
024106
(
2013
).
247.
A. G.
Dijkstra
and
V. I.
Prokhorenko
,
J. Chem. Phys.
147
,
064102
(
2017
).
248.
L.
Chen
,
R.
Zheng
,
Y.
Jing
, and
Q.
Shi
,
J. Chem. Phys.
134
,
194508
(
2011
).
249.
J.
Zhu
,
S.
Kais
,
P.
Rebentrost
, and
A.
Aspuru-Guzik
,
J. Phys. Chem. B
115
,
1531
(
2011
).
250.
L.
Banchi
,
G.
Costagliola
,
A.
Ishizaki
, and
P.
Giorda
,
J. Chem. Phys.
138
,
184107
(
2013
).
251.
L.
Zhu
,
H.
Liu
, and
Q.
Shi
,
New J. Phys.
15
,
095020
(
2013
).
252.
M.
Schröter
,
T.
Pullerits
, and
O.
Kühn
,
Ann. Phys.
527
,
536
(
2015
).
253.
Y.
Fujihashi
,
G. R.
Fleming
, and
A.
Ishizaki
,
J. Chem. Phys.
142
,
212403
(
2015
).
254.
Y.
Fujihashi
,
L.
Chen
,
A.
Ishizaki
,
J.
Wang
, and
Y.
Zhao
,
J. Chem. Phys.
146
,
044101
(
2017
).
255.
Y.
Fujihashi
,
M.
Higashi
, and
A.
Ishizaki
,
J. Phys. Chem. Lett.
9
(
17
),
4921
(
2018
).
256.
J.
Seibt
and
T.
Mančal
,
Chem. Phys.
515
,
129
(
2018
).
257.
S.
Saito
,
M.
Higashi
, and
G. R.
Fleming
,
J. Phys. Chem. B
123
,
9762
(
2019
).
258.
Y.
Yao
,
W.
Yang
, and
Y.
Zhao
,
J. Chem. Phys.
140
,
104113
(
2014
).
259.
J. M.
Moix
,
Y.
Zhao
, and
J.
Cao
,
Phys. Rev. B
85
,
115412
(
2012
).
260.
Z.
Tong
,
Z.
Huai
,
Y.
Mei
, and
Y.
Mo
,
J. Chem. Phys.
152
,
135101
(
2020
).
261.
N. T.
Phuc
and
A.
Ishizaki
,
J. Phys. Chem. Lett.
9
,
1243
(
2018
).
262.
N. T.
Phuc
and
A.
Ishizaki
,
Phys. Rev. B
99
,
064301
(
2019
).
263.
N. T.
Phuc
and
A.
Ishizaki
,
Phys. Rev. Res.
1
,
033019
(
2019
).
264.
H.
Ito
and
Y.
Tanimura
,
J. Chem. Phys.
144
,
074201
(
2016
).
265.
M. F.
Gelin
,
Y.
Tanimura
, and
W.
Domcke
,
J. Chem. Phys.
139
,
214302
(
2013
).
266.
E.
Mangaud
,
R.
Puthumpally-Joseph
,
D.
Sugny
,
C.
Meier
,
O.
Atabek
, and
M.
Desouter-Lecomte
,
New J. Phys.
20
,
043050
(
2018
).
267.
T.
Joutsuka
and
Y.
Tanimura
,
Chem. Phys. Lett.
457
,
237
(
2008
).
268.
Y.
Mo
,
X.
Zheng
,
G.
Chen
, and
Y.
Yan
,
J. Phys.: Condens Matter
21
,
355301
(
2009
).
269.
Z. H.
Li
,
N.
Tong
,
Z.
Zheng
,
D.
Hou
,
J.
Wei
,
J.
Hu
, and
Y.
Yan
,
Phys. Rev. Lett.
109
,
266403
(
2012
).
270.
L.
Han
,
A.
Ullah
,
Y.-A.
Yan
,
X.
Zheng
,
Y.
Yan
, and
V.
Chernyak
,
J. Chem. Phys.
152
,
204105
(
2020
).
271.
A.
Ullah
,
L.
Han
,
Y.-A.
Yan
,
X.
Zheng
,
Y.
Yan
, and
V.
Chernyak
,
J. Chem. Phys.
152
,
204106
(
2020
).
272.
X.
Zheng
,
J.
Jin
, and
Y.
Yan
,
New J. Phys.
10
,
093016
(
2008
).
273.
X.
Zheng
,
J.
Jin
,
S.
Welack
,
M.
Luo
, and
Y.
Yan
,
J. Chem. Phys.
130
,
164708
(
2009
).
274.
Y.
Cheng
,
W.
Hou
,
Y.
Wang
,
Z.
Li
,
J.
Wei
, and
Y.
Yan
,
New J. Phys.
17
,
033009
(
2015
).
275.
Y.
Cheng
,
J.
Wei
, and
Y.
Yan
,
Europhys. Lett.
112
,
57001
(
2015
).
276.
R.
Härtle
and
A. J.
Millis
,
Phys. Rev. B
90
,
245426
(
2014
).
277.
L.-V.
Ye
,
D.
Hou
,
X.
Zheng
,
Y. J.
Yan
, and
M. D.
Ventra
,
Phys. Rev. B
91
,
205106
(
2015
).
278.
R.
Härtle
,
G.
Cohen
,
D. R.
Reichman
, and
A. J.
Millis
,
Phys. Rev. B
92
,
085430
(
2015
).
279.
J.
Okamoto
,
L.
Mathey
, and
R.
Härtle
,
Phys. Rev. B
94
,
235411
(
2016
).
280.
A.
Erpenbeck
,
C.
Schinabeck
,
U.
Peskin
, and
M.
Thoss
,
Phys. Rev. B
97
,
235452
(
2018
).
281.
A.
Erpenbeck
and
M.
Thoss
,
J. Chem. Phys.
151
,
191101
(
2019
).
282.
H.
Xie
,
Y.
Kwok
,
Y.
Zhang
,
F.
Jiang
,
X.
Zheng
,
Y.
Yan
, and
G.
Chen
,
Phys. Status Solidi B
250
,
2481
(
2013
).
283.
Y.
Wang
,
X.
Zheng
,
B.
Li
, and
J.
Yang
,
J. Chem. Phys.
141
,
084713
(
2014
).