Equations of motion for the density matrix
HEOMs are developed to describe the time evolution of the density matrix
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.[2][3][4]
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}\left[{\ {\hat {p}}^{2} \over {2}}+{\frac {1}{2}}{\hat {x}}_{j}^{2}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dcf75e38f40c4f1b8c7ce4635124b9bae46c4aba)
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 Drude spectral function

and a high temperature heat bath, taking the time derivative of the system density matrix, and writing it in hierarchical form yields (
)

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

with the inverse temperature
and the following "super-operator" notation:

The counter
provides for
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
differential equations. This "terminator"
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
term is[5]

The closing line of the hierarchy is thus:
.
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.
Implementations
The HEOM method is implemented in a number of freely available codes. A number of these are at the website of Yoshitaka Tanimura[6] including a version for GPUs [7] which used improvements introduced by David Wilkins and Nike Dattani.[8] The nanoHUB version provides a very flexible implementation.[9] An open source parallel CPU implementation is available from the Schulten group.[10]