# Hierarchical equations of motion

The Hierarchical equations of motion (HEOM) technique derived by Yoshitaka Tanimura and Ryogo Kubo in 1989,[1] and generalized to arbitrary temperature and arbitrary spectral densities by N. Dattani in 2012,[2][3] is a non-perturbative approach developed to study the evolution of a density matrix ${\displaystyle \rho (t)}$ of quantum dissipative systems. The method can treat system-bath interaction non-perturbatively as well as non-Markovian noise correlation times without the hindrance of the typical assumptions that conventional Redfield (master) equations suffer from such as the Born, Markovian and rotating-wave approximations. HEOM is applicable even at low temperatures where quantum effects are not negligible.

The hierarchical equation of motion for a system in a harmonic Markovian bath is[4]

${\displaystyle {\frac {\partial }{\partial t}}{\hat {\rho }}_{n}=-i({\hat {H}}_{A}+n\gamma ){\hat {\rho }}_{n}-{1 \over \hbar }{\hat {V}}^{\times }{\hat {\rho }}_{n+1}+{in \over \hbar }{\hat {\Theta }}{\hat {\rho }}_{n-1}}$

## Hierarchical Equations of Motion

HEOMs are developed to describe the time evolution of the density matrix ${\displaystyle \rho (t)}$ for an open quantum system. It is a non-perturbative, non-Markovian approach to propagating in time a quantum state. Motivated by the path integral formalism presented by Feynman and Vernon, Tanimura derive the HEOM from a combination of statistical and quantum dynamical techniques.[4][5][6] Using a two level spin-boson system Hamiltonian

${\displaystyle {\hat {H}}={\hat {H}}_{A}({\hat {a}}^{+},{\hat {a}}^{-})+V({\hat {a}}^{+},{\hat {a}}^{-})\sum _{j}c_{j}{\hat {x}}_{j}+\sum _{j}{\big [}{\ {\hat {p}}^{2} \over {2}}+{\frac {1}{2}}{\hat {x}}_{j}^{2}{\big ]}}$

Characterising the bath phonons by the spectral density ${\displaystyle J(\omega )=\sum \nolimits _{j}c_{j}^{2}\delta (\omega -\omega _{j})}$

By writing the density matrix in path integral notation and making use of Feynman-Vernon influence functional, all the bath coordinates in the interaction terms can be grouped into this influence functional which in some specific cases can be calculated in closed form. Assuming a high temperature heat bath with the Drude spectral distribution ${\displaystyle J(\omega )=\hbar \eta \gamma ^{2}\omega /\pi (\gamma ^{2}+\omega ^{2})}$ and taking the time derivative of the path integral form density matrix the equation and writing it in hierarchal form yields

${\displaystyle {\frac {\partial }{\partial t}}{\hat {\rho }}_{N}=-i({\hat {H}}_{A}+N\gamma ){\hat {\rho }}_{N}-{1 \over \gamma \hbar ^{2}}{\hat {V}}^{\times }{\hat {\Theta }}{\hat {\rho }}_{N}+{iN \over \hbar }{\hat {\Theta }}{\hat {\rho }}_{N-1}}$

where ${\displaystyle \Theta }$ destroys system excitation and hence can be referred to as the relaxation operator.

${\displaystyle {\hat {\Theta }}=-{\frac {n\gamma }{\beta }}{\big (}{\hat {V}}^{\times }-i{\frac {\beta \hbar \gamma }{2}}{\hat {V}}^{\circ }{\big )}}$

The second term in ${\displaystyle {\hat {\Theta }}}$ is the temperature correction term with the inverse temperature ${\displaystyle \beta =1/k_{B}T}$ and the "Hyper-operator" notation is introduced.

${\displaystyle {\hat {A}}^{\times }{\hat {\rho }}={\hat {A}}{\hat {\rho }}-{\hat {\rho }}{\hat {A}}}$

${\displaystyle {\hat {A}}^{\circ }{\hat {\rho }}={\hat {A}}{\hat {\rho }}+{\hat {\rho }}{\hat {A}}}$

As with the Kubo's Stochastic Liouville Equation in hierarchal form, the counter ${\displaystyle n}$ can go up to infinity which is a problem numerically, however Tanimura and Kubo provide a method by which the infinite hierarchy can be truncated to a finite set of ${\displaystyle N}$ differential equations where ${\displaystyle N}$ is determined by some constraint sensitive to the characteristics of the system i.e. frequency, amplitude of fluctuations, bath coupling etc. The "Terminator" defines the depth of the hierarchy. A simple relation to eliminate the ${\displaystyle {\hat {\rho }}_{n+1}}$ term is found. ${\displaystyle \ {\hat {\rho }}_{N+1}=-{\hat {\Theta }}{\hat {\rho }}_{N}/\hbar \gamma }$.[7]

The statistical nature of the HEOM approach allows information about the bath noise and system response to be encoded into the equation of motion doctoring the infinite energy problem of Kubo's SLE by introducing the relaxation operator ensuring a return to equilibrium.

### Low Temperature

HEOM can be derived for a variety of spectral distributions i.e. Drude,[8] Brownian,[9] Lorentzian,[10], and Sub-Ohmic [11], as well as combinations of such distributions at any temperature (see Table 1 of [2]), or even arbitrary bath response functions [3] at any temperature.

In the Drude case, by modifying the correlation function that describes the noise correlation function strongly non-Markovian and non-perturbative system-bath interactions can be dealt with[4][8]

${\displaystyle {\frac {\partial }{\partial t}}{\hat {\rho }}_{n,j_{1},..,j_{K}}=-{\bigg [}{i \over \hbar }{\hat {H}}_{A}^{\times }+n\gamma +\sum _{k=1}^{K}(j_{k}\nu _{k}-{1 \over {\nu _{k}\hbar ^{2}}}{\hat {V}}^{\times }{\hat {\Theta }}_{k})+{\hat {\Gamma }}_{0}{\bigg ]}{\hat {\rho }}_{n,j_{1},..,j_{K}}-{i \over \hbar }{\hat {V}}^{\times }{\bigg [}{\hat {\rho }}_{(n+1),j_{1},..,j_{K}}+\sum _{k=1}^{K}{\hat {\rho }}_{n,j_{1},..,(j_{k}+1),..,j_{K}}{\bigg ]}-{in \over \hbar }{\hat {\Theta }}_{0}{\hat {\rho }}_{(n-1),j_{1}..j_{K}}-\sum _{k=1}^{K}{ij_{k} \over \hbar }{\hat {\Theta }}_{k}{\hat {\rho }}_{n,j_{1},..,(j_{k}-1),..,j_{K}}}$ In this equation, only ${\displaystyle {\hat {\rho }}_{0,0,...,0}}$ contains all order of system bath interactions with other elements ${\displaystyle {\hat {\rho }}_{n,j_{1}...j_{K}}}$ being auxiliary terms, moving deeper into the hierarchy, the order of interactions decreases, which is contrary to usual perturbative treatments of such systems. ${\displaystyle {\hat {\Theta }}_{k}=c_{k}{\hat {V}}^{\times }}$ where ${\displaystyle c_{k}}$ is a constant determined in the correlation function.

${\displaystyle {\hat {\Gamma }}_{0}\equiv {\eta \over {\beta \hbar ^{2}}}{\big (}1-{\beta \over \gamma n}c_{0}{\big )}{\hat {V}}^{\times }{\hat {V}}^{\times }}$

This ${\displaystyle {\hat {\Gamma }}_{0}}$ term arises from the Matsubara cut-off term introduced to the correlation function and thus holds information about the memory of the noise.

Below is the terminator for the HEOM

${\displaystyle {\frac {\partial }{\partial t}}{\hat {\rho }}_{n,j_{1},...,j_{K}}\simeq -{\big (}{i \over \hbar }{\hat {H}}_{A}^{\times }-\sum _{k=1}^{K}{1 \over {\nu _{k}\hbar ^{2}}}{\hat {V}}^{\times }{\hat {\Theta }}_{k}+{\hat {\Gamma }}_{0}{\big )}{\hat {\rho }}_{n,j_{1},...,j_{K}}}$

Performing a Wigner transformation on this HEOM, the quantum Fokker-Planck equation with low temperature correction terms emerges.[7][12]

### Arbitrary Spectral Density and Arbitrary Temperature

It was pointed out by Dattani et al. in 2012 that the HEOM method can be employed as long as the bath correlation function is written as a sum of exponentials,[2] and therefore an arbitrarily complicated spectral distribution function ${\displaystyle J(\omega )}$ can be fitted to any of the forms listed in Table 1 of [2] whose bath response function can analytically be written as a sum of exponentials, and then the HEOM can be applied for that spectral density at arbitrary temperature.

In a subsequent paper,[3] it was suggested that the bath response function be fitted directly to a sum of exponentials rather than fitting the spectral density to one of the forms in Table 1 of [2] and then calculating the bath response function as a sum of exponentials analytically. This direct approach requires non-linear fitting to a complex function that has both real and imaginary parts and therefore the choice of starting parameters for the non-linear fitting is non-trivial, but a recipe for choosing starting parameters was given in,[3] and the non-linear complex-valued fitting procedure can be carried out completely by the program BATHFIT.[13]

### Algorithm and Computational Cost

When the open quantum system is represented by ${\displaystyle M}$ levels and ${\displaystyle M}$ baths with each bath response function represented by ${\displaystyle K}$ exponentials, a hierarchy with ${\displaystyle {\mathcal {N}}}$ layers will contain:

${\displaystyle {\frac {\left(MK+{\mathcal {N}}\right)!}{\left(MK\right)!{\mathcal {N}}!}}}$

matrices, each with ${\displaystyle M^{2}}$ complex-valued (containing both real and imaginary parts) elements. Therefore, the limiting factor in HEOM calculations is the amount of RAM required, since if one copy of each matrix is stored, the total RAM required would be:

${\displaystyle 16M^{2}{\frac {\left(MK+{\mathcal {N}}\right)!}{\left(MK\right)!{\mathcal {N}}!}}}$

bytes if using double-precision, and half of that if using single-precision (however, numerically exact open quantum system dynamics has been shown to require double-precision in Fig. 2 of [14]).

Further RAM may be necessary depending on how the system of differential equations is numerically solved. The algorithm which requires the least amount of memory to store the auxiliary density operators (ADOs), yet is still stable numerically, was presented in 2015 by Wilkins and Dattani,[15] and has been implemented for GPUs by Tsuchimoto and Tanimura.[16] An implementation for CPUs is given in the original paper.[15]

### Applications

The largest number of levels coupled to the bath, that have been successfully treated with the HEOM using one bath, with the bath response function represented by one exponential is 4096 (equivalent to 12 qubits) in an application to the Ising model.[16] The largest number of levels treated using two exponentials and ${\displaystyle M}$ baths is 24 in an application to the Fenna-Matthews-Olson complex.[15]

## References

1. ^ Tanimura, Yoshitaka; Kubo, Ryogo (1989), "Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath", J. Phys. Soc. Jpn., 58: 101–114
2. Dattani, Nike (2012), "Analytic influence functionals for numerical Feynman integrals in most open quantum systems" (PDF), Quantum Physics Letters, 1: 35–45
3. ^ a b c d
4. ^ a b c Tanimura, Yoshitaka (1990), "Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath", Phys. Rev. A, 41: 6676–6687
5. ^ Tanimura, Yoshitaka (2006), "Stochastic Liouville, Langevin, Fokker-Planck, and Master Equation Approaches to Quantum Dissipative Systems", J. Phys. Soc. Jpn., 75: 082001
6. ^ Tanimura, Yoshitaka (2014), "Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities", J. Chem. Phys., 141: 044114
7. ^ a b Tanimura, Yoshitaka; Wolynes, Peter (1991), "Quantum and classical Fokker-Planck equations for a Gaussian-Markovian noise bath", Phys. Rev. A, 43: 4131–4142
8. ^ a b Ishizaki, Akihito; Tanimura, Yoshitaka (2005), "Quantum Dynamics of System Strongly Coupled to Low-Temperature Colored Noise Bath: Reduced Hierarchy Equations Approach", J. Phys. Soc. Jpn., 74: 3131–3134
9. ^ Tanaka, Midori; Tanimura, Yoshitaka (2009), "Quantum Dissipative Dynamics of Electron Transfer Reaction System: Nonperturbative Hierarchy Equations Approach", J. Phys. Soc. Jpn., 78: 073802 (2009)
10. ^ Ma, Jian; Sun, Zhe; Wang, Xiaoguanag; Nori, Franco (2012), "Entanglement dynamics of two qubits in a common bath", Phys. Rev. A, 85: 062323 (2012)
11. ^ Duan, Chenru; Zhoufei, Tang; Jianshu, Cao; Jianlan, Wu (2017), "Zero-temperature localization in a sub-Ohmic spin-boson model investigated by an extended hierarchy equation of motion", Phys. Rev. B, 95: 214308
12. ^ Tanimura, Yoshitaka (2015), "Real-time and imaginary-time quantum hierarchical Fokker-Planck equations", J. Chem. Phys., 141: 044114
13. ^ Dattani, Nike (2012), "BATHFIT", GitHub
14. ^ Dattani, Nike (2013), "FeynDyn: A MATLAB program for fast numerical Feynman integral calculations for open quantum system dynamics on GPUs", Computer Physics Communications, 184 (12): 2828–2833, arXiv:, Bibcode:2013CoPhC.184.2828D, doi:10.1016/j.cpc.2013.07.001
15. ^ a b c Wilkins, David; Dattani, Nike (2015), "Why Quantum Coherence Is Not Important in the Fenna–Matthews–Olsen Complex", Journal of Chemical Theory and Computation, 11 (7): 3411–3419, arXiv:, doi:10.1021/ct501066k
16. ^ a b Tsuchimoto, Masashi; Tanimura, Yoshitaka (2015), "Spins Dynamics in a Dissipative Environment: Hierarchal Equations of Motion Approach Using a Graphics Processing Unit (GPU)", Journal of Chemical Theory and Computation, 11 (7): 3859–3865