Hierarchical equations of motion

The hierarchical equations of motion (HEOM) technique derived by Yoshitaka Tanimura and Ryogo Kubo in 1989, is a non-perturbative approach developed to study the evolution of a density matrix $\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

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

where the superscript $^{\times }$ denoting a commutator and the temperature-dependent super-operator ${\hat {\Theta }}$ are defined below. The parameter $\gamma$ is the frequency width of the Drude spectral function $J(\omega )$ (see below).

Equations of motion for the density matrix

HEOMs are developed to describe the time evolution of the density matrix $\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. Using a two level spin-boson system Hamiltonian

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

By writing the density matrix in path integral notation and making use of Feynman–Vernon influence functional, all the bath coordinates $x_{j}$ in the interaction terms can be grouped into this influence functional which in some specific cases can be calculated in closed form. Assuming a Drude spectral function

$J(\omega )=\sum \nolimits _{j}c_{j}^{2}\delta (\omega -\omega _{j})={\frac {\hbar \eta \gamma ^{2}\omega }{\pi (\gamma ^{2}+\omega ^{2})}}$

and a high temperature heat bath, taking the time derivative of the system density matrix, and writing it in hierarchical form yields ($n=0,1,\ldots$)

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

Here $\Theta$ reduces the system excitation and hence is referred to as the relaxation operator:

${\hat {\Theta }}=-{\frac {\eta \gamma }{\beta }}\left({\hat {V}}^{\times }-i{\frac {\beta \hbar \gamma }{2}}{\hat {V}}^{\circ }\right)$

with the inverse temperature $\beta =1/k_{B}T$ and the following "super-operator" notation:

${\begin{aligned}{\hat {A}}^{\times }{\hat {\rho }}&={\hat {A}}{\hat {\rho }}-{\hat {\rho }}{\hat {A}}\\{\hat {A}}^{\circ }{\hat {\rho }}&={\hat {A}}{\hat {\rho }}+{\hat {\rho }}{\hat {A}}\end{aligned}}$

The counter $n$ provides for $n=0$ the system density matrix. As with Kubo's stochastic Liouville equation in hierarchical form, it goes up to infinity in the hierarchy which is a problem numerically. Tanimura and Kubo, however, provide a method by which the hierarchy can be truncated to a finite set of $N$ differential equations. This "terminator" $N$ defines the depth of the hierarchy and is determined by some constraint sensitive to the characteristics of the system, i.e. frequency, amplitude of fluctuations, bath coupling etc. A simple relation to eliminate the ${\hat {\rho }}_{n+1}$ term is

${\hat {\rho }}_{N+1}=-{\hat {\Theta }}{\hat {\rho }}_{N}/\hbar \gamma .$

The closing line of the hierarchy is thus:

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

The HEOM approach allows information about the bath noise and system response to be encoded into the equations of motion. It cures the infinite energy problem of Kubo's stochastic Liouville equation by introducing the relaxation operator that ensures a return to equilibrium.

Computational cost

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

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

matrices, each with $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:

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

bytes (assuming double-precision).

Implementations

The HEOM method is implemented in a number of freely available codes. A number of these are at the website of Yoshitaka Tanimura including a version for GPUs which used improvements introduced by David Wilkins and Nike Dattani. The nanoHUB version provides a very flexible implementation. An open source parallel CPU implementation is available from the Schulten group.

See also

• Quantum master equation
• Open quantum system
• Fokker–Planck equation
• Quantum dynamical semigroup
• Quantum dissipation