2. Scientific Background
2.1. Key definitions of HEOM
The dynamics of a quantum system embedded into an environment (“bath”)
can often be described by the Hamiltonian of the form:
\(H=\sum_{n,m=0}^{N-1}{\left.\ |n\right\rangle H_{\text{nm}}\left\langle m|\right.\ }+\sum_{n=0}^{N-1}{\left.\ |n\right\rangle\sum_{b=0}^{N_{n}-1}\left(\frac{p_{b,n}^{2}}{2}+\frac{1}{2}\omega_{b,n}^{2}x_{b,n}^{2}\right)\left\langle n|\right.\ }+\sum_{n=0}^{N-1}{\left.\ |n\right\rangle F_{n}\left\langle n|\right.\ }\).
(1)
Here, the first sum represents the electronic Hamiltonian of the quantum
system, the second term represents the phonon Hamiltonian describing the
dynamics of the bath modes on every electronic state, and the last term
represents the linear electron-phonon (vibronic) coupling Hamiltonian
describing the perturbation of the quantum system by the environmental
degrees of freedom. The matrix elements \(H_{\text{nn}}\) correspond to
state energies, and \(H_{\text{nm}}=H_{\text{mn}},\ n\neq m\)correspond to electronic (or excitonic, depending on the interpretation)
couplings. The numbers are assumed to be independent of the bath degrees
of freedom (DOF). The number N represents the total number of quantum
states, \(N_{n}\)is the number of the bath modes coupled to the n-th
quantum state, \(p_{b,n}\) and \(x_{b,n}\) are the mass-weighted normal
modes’ momenta and coordinates for a mode b in a quantum staten , and \(\omega_{b,n}\) are the normal mode frequencies for a
mode b in a quantum state n . Finally, the variable\(F_{n}\) describes the way a quantum state n is coupled to
various modes and can be expressed as the following:
\(F_{n}=\sum_{b=0}^{N_{b}-1}{f_{b,n}x_{b,n}}\). (2)
Often, the quantum states \(\left\{\left.\ |n\right\rangle\right\}\)are understood as site (molecular) states, but one can also associate
them with a set of quantum states of a single system. The bath
operators, \(F_{n}\), can be more general and can be considered as a
user-defined control parameter describing system-bath interactions.
The parameters, \(\left\{f_{b,n}\right\}\), reflect the type of the
electron-phonon interactions the bath induces, which can be quantified
by the bath spectral density:
\(J_{n}\left(\omega\right)=\frac{\pi}{2}\sum_{b=0}^{N_{n}-1}\frac{f_{b,n}^{2}}{\omega_{b,n}}\delta\left(\omega-\omega_{b,n}\right)\).
(3)
In many chemical applications, the baths are well described by the Debye
spectral density, which describes an overdamped Brownian motion of the
energy gaps in a high-temperature regime:
\(J_{n}\left(\omega\right)=\frac{\text{ηγ}\omega^{2}}{\omega^{2}+\gamma^{2}}\).
(4)
Here, \(\eta\) is the bath reorganization energy, and\(\gamma=\frac{\Gamma}{\hslash}\) is the frequency that corresponds to
the system-bath coupling energy, \(\Gamma\). Throughout this account and
in the Libra software, we utilize the atomic units, in which\(\hslash=1\), so \(\hslash\) may be omitted in some equations, and,
numerically, \(\gamma=\Gamma\) (although the two variables have
different units). Within the HEOM framework, the dynamics of the
system’s reduced density matrix (RDM), \(\rho_{\mathbf{0}}\), can be
described by the equations:25,31,32
\({\dot{\rho}}_{\mathbf{n}}=-i\left[H,\rho_{\mathbf{n}}\right]-\sum_{m=0}^{M-1}{\left(\sum_{k=0}^{K}{n_{\text{mk}}\gamma_{\text{mk}}}\right)\rho_{\mathbf{n}}}+\rho_{\mathbf{n}}^{\left(+\right)}+\rho_{\mathbf{n}}^{\left(-\right)}+T_{\mathbf{n}}\).
(5a)
\(\rho_{\mathbf{n}}^{(+)}=-i\sum_{m=0}^{M-1}\left[Q_{m},\sum_{k=0}^{K}\rho_{\mathbf{n}_{\text{mk}}^{\mathbf{+}}}\right]\).
(5b)
\(\rho_{\mathbf{n}}^{(-)}=-i\sum_{m=0}^{M-1}{\sum_{k=0}^{K}{n_{\text{mk}}(F_{\text{mk}}c_{\text{mk}}\rho_{\mathbf{n}_{\text{mk}}^{\mathbf{-}}}-c_{\text{mk}}^{*}\rho_{\mathbf{n}_{\text{mk}}^{\mathbf{-}}}F_{\text{mk}})}}\).
(5c)
\(T_{\mathbf{n}}=\sum_{m=0}^{M-1}{\Delta_{K}\left[Q_{m},\left[Q_{m},\rho_{\mathbf{n}}\right]\right]}\).
(5d)
The HEOM method provides an exact, non-perturbative way of describing
the dynamics of a quantum system (the reduced density matrix evolution)
in the presence of complex environment, whose degrees of freedom are
removed through integration. Although we report the implementation of
the HEOM for the spectral density given in Eq. 4, the formulations for
more general spectral densities are possible.4,33–36
The first term on the right-hand side of Eq. 5a describes the quantum
dynamics of an isolated quantum system, the second term – the “quantum
friction” - is due to the presence of the environment. The terms Eq. 5b
and 5c describe coupling between ADMs of various tiers, and the term Eq.
5d describes the correction to the HEOM truncation.
Here, \(\left\{\gamma_{k}\right\}\) and \(\left\{c_{k}\right\}\) are
the Matsubara frequencies and expansion coefficients, respectively,
which describe the decay of the autocorrelation function of the
collective coordinate of the bath:
\(C\left(t>0\right)=\sum_{k=0}^{K}{c_{k}exp(-\gamma_{k}t)}\).
(6)
Here, K defines the number of Matsubara frequencies. The Matsubara
frequencies are defined by the system-bath “interaction” time,\(1/\gamma\) (so that \(\gamma\) is a characteristic frequency of the
bath) and by temperature \(\beta=\frac{1}{k_{B}T}\):
\(\gamma_{0}=\gamma\). (7a)
\(\gamma_{n}=\frac{2\text{πn}}{\beta}=2\text{πn}k_{B}T,n\geq 1\),
(7b)
With the definition of the spectral density, Eq. 4, the Matsubara
expansion coefficients are given by:31
\(c_{0}=\frac{1}{2}\text{γη}\ \left(\left[\tan\left(\frac{\gamma}{2k_{B}T}\right)\right]^{-1}-i\right)\ \),
(8a)
\(c_{k}=\frac{4\text{nπηγ}}{\left(2\text{kπ}\right)^{2}-\left(\text{βγ}\right)^{2}}=\frac{4\text{nπηγ}}{\beta^{2}\left[\left(\frac{2\text{nπ}}{\beta}\right)^{2}-\left(\gamma\right)^{2}\right]}=2\eta k_{B}T\frac{\gamma_{0}\gamma_{n}}{\gamma_{n}^{2}-\gamma_{0}^{2}},\ k\geq 1\).
(8b)
The parameter \(\Delta_{K}\) entering Eq. 5d is a residual sum that is
used to truncate the hierarchy:
\(\Delta_{K}=\sum_{n=0}^{\infty}\frac{c_{K+n}}{\gamma_{K+n}}\).
(9)
Finally, the matrices\(Q_{n}=\sum_{a,b=0}^{M-1}{\left.\ |a\right\rangle Q_{ab,n}\left\langle b|\right.\ }\)are the matrices describing how a phonon n is coupled to all
other electronic states,\(\{\left.\ |a\right\rangle,\ a=0,\ldots,M-1\)}. Here, M is the
number of phonon modes. The off-diagonal terms \(Q_{ab,n}\) correspond
to the coupling of a phonon \(n\) to electronic couplings between pairs
of states \(a\) and \(b\). In the simplest case, when each electronic
state is coupled to a single (distinct) phonon, the operator \(Q\) takes
the form: \(Q_{n}=\left.\ |n\right\rangle\left\langle n|\right.\ \)and \(M=N\).
Eq. 5d is a good approximation when the hierarchy of equations is
truncated at a finite number of equations. It is not present in the
original formulation, with the infinite hierarchy of equations.
2.2. Indexing and hierarchy
tracking
Within HEOM, a set of auxiliary density matrices (ADMs),\(\{\rho_{\mathbf{n}}\}\), is introduced, with the \(\rho_{\mathbf{0}}\)being the RDM of the quantum system, and is the main property of
interest. The ADMs are indexed by the multi-dimensional index (bolded):
\(\mathbf{n}=\left(n_{00},n_{01},\ldots,\ n_{0K},\ n_{10},n_{11},\ldots,\ n_{1K},\ldots,n_{M-1,0},n_{M-1,1},\ldots,\ n_{M-1,K}\right)\).
(10a)
Here, \(M\) is the number of phonon modes in the system and K+1 is the
number of Matsubara frequencies. For simplicity, we take the indexing
convention that quantum states’ indexing starts with 0 and ends at\(M-1\), whereas the Matsubara frequencies’ indexing starts with 0 and
ends at K. Thus, the length of the
multi-dimensional index vector is \(d=M*(K+1)\).
\(\mathbf{n}_{\mathbf{\text{mk}}}^{\mathbf{+}}=\left(n_{00},\ldots,\ n_{0K},\ldots,n_{m0},\ldots,\ n_{\text{mk}}+1,\ldots,n_{\text{mK}},\ldots,n_{M-1,0},n_{M-1,1},\ldots,n_{M-1,K}\right)\).
(10b)
\(\mathbf{n}_{\mathbf{\text{mk}}}^{\mathbf{-}}=\left(n_{00},\ldots,\ n_{0K},\ldots,n_{m0},\ldots,\ n_{\text{mk}}-1,\ldots,n_{\text{mK}},\ldots,n_{M-1,0},n_{N-1,1},\ldots,n_{M-1,K}\right).\)(10c)
The tier of the ADM, \(n\), is defined as the sum of all elements of the
multi-component index vector \(\mathbf{n}\):
\(n=\text{tier}\left(\mathbf{n}\right)=\sum_{m=0}^{M-1}{\sum_{k=0}^{K}n_{\text{mk}}}\).
(11)
2.3. Scaled HEOM
In the current version of Libra, we have also implemented the “scaled”
HEOM of Shi et al.25 According to the method, the
scaled ADMs are constructed as:
\({\tilde{\rho}}_{\mathbf{n}}=\left(\prod_{m=0}^{M-1}{\prod_{k=0}^{K}{n_{\text{mk}}!\left|c_{\text{mk}}\right|^{n_{\text{mk}}}}}\right)^{-1/2}\rho_{\mathbf{n}}\).
(12)
Note that \({\tilde{\rho}}_{\mathbf{0}}=\rho_{\mathbf{0}}\), that is
one may propagate the dynamics in terms of the scaled ADMs, but the
evolution of the first (zero tier) member of the hierarchy would
correspond to that of the original reduced density matrix of the quantum
system. In terms of the scaled ADMs, the HEOM takes the form:
\(\frac{{d\tilde{\rho}}_{\mathbf{n}}}{\text{dt}}=-i\left[H,{\tilde{\rho}}_{\mathbf{n}}\right]-\sum_{m=0}^{M-1}{\left(\sum_{k=0}^{K}{n_{\text{mk}}\gamma_{\text{mk}}}\right){\tilde{\rho}}_{\mathbf{n}}}+{\tilde{\rho}}_{\mathbf{n}}^{\left(+\right)}+{\tilde{\rho}}_{\mathbf{n}}^{\left(-\right)}+{\tilde{T}}_{\mathbf{n}}\).
(13a)
\({\tilde{\rho}}_{\mathbf{n}}^{(+)}=-i\sum_{m=0}^{M-1}\left[Q_{m},\sum_{k=0}^{K}{\sqrt{\left(n_{\text{mk}}+1\right)\left|c_{\text{mk}}\right|}{\tilde{\rho}}_{\mathbf{n}_{\text{mk}}^{\mathbf{+}}}}\right]\).
(13b)
\({\tilde{\rho}}_{\mathbf{n}}^{(-)}=-i\sum_{m=0}^{M-1}{\sum_{k=0}^{K}{\sqrt{n_{\text{mk}}/\left|c_{\text{mk}}\right|}(F_{\text{mk}}c_{\text{mk}}{\tilde{\rho}}_{\mathbf{n}_{\text{mk}}^{\mathbf{-}}}-c_{\text{mk}}^{*}{\tilde{\rho}}_{\mathbf{n}_{\text{mk}}^{\mathbf{-}}}F_{\text{mk}})}}\).
(13c)
\({\tilde{T}}_{\mathbf{n}}=\sum_{m=0}^{M-1}{\Delta_{K}\left[Q_{m},\left[Q_{m},{\tilde{\rho}}_{\mathbf{n}}\right]\right]}\).
(13d)
These equations are isomorphic to Eqs. 5 and reduce to them when\(\sqrt{\left(n_{\text{mk}}+1\right)\left|c_{\text{mk}}\right|}\rightarrow 1\)in Eq. 13b and\(\sqrt{n_{\text{mk}}/\left|c_{\text{mk}}\right|}\rightarrow n_{\text{mk}}\)in Eq. 13c. The main advantage of the scaled HEOM Eqs. 13 is that the
scaled densities become negligible for high tiers much more quickly than
the unscaled ones. Based on this property, the truncation and filtering
can become much more efficient, allowing the converged results to be
achieved using smaller \(K\) and \(L\) values, where \(L\) is the
maximal tier of ADM.
2.4. Filtering
The main complexity of the HEOM is the large (factorially-growing)
number of ADMs and the corresponding equations to solve. To combat this,
the hierarchy is truncated at a certain tier level, \(L\), and a minimal
number of Matsubara frequencies should be used. However, the results of
the calculations should be converged with respect to both parameters.
The use of a finite number of ADMs is compensated by the truncation
terms, Eqs. 5d and 13d.
In addition to the use of scaled version of HEOM and truncation
correction, the current implementation allows filtering equations based
on the \(\left|\rho_{\mathbf{n}}\right|<tol_{1}\) or\(\left|{\dot{\rho}}_{\mathbf{n}}\right|<tol_{2}\) criteria (for
scaled or unscaled ADMs and their derivatives). The norm of the matrix
(whether ADM or its time-derivative) is computed as the maximal
magnitude of all matrix elements.25 The first
criterion is used to turn off the computation of any terms involving
very small ADMs. The second criterion is used to turn off the
computation of the time-derivatives of the ADMs when such derivatives
are expected to be very small. We construct two lists to keep track of
the information on whether to flag out (with 0) or flag in (with 1) the
use of the corresponding ADMs or recalculation of the corresponding
time-derivatives. These lists can be updated with a user-defined
frequency to minimize the overhead in their updates. This approach is
similar to the Verlet-list technique often used in classical molecular
dynamics simulations.
2.5. Spectra
The system’s RDM, propagated via HEOM,\(\rho_{\mathbf{0}}\left(t\right)\), can be used to compute the line
shapes of the absorption spectrum as a Fourier transform of the dipole
autocorrelation function (ACF):
\(I\left(\omega\right)=\frac{1}{\pi}\text{\ Re}\left[\int_{0}^{\infty}{\text{dt}e^{i\omega t}}\left\langle u\left(t\right)u\left(0\right)\right\rangle_{g}\right]\),
(14)
The ACF is computed as a trace of the product of the time-dependent
density matrix, \(\rho_{\mathbf{0}}\left(t\right)\), and the
(time-independent) transition dipole moment matrix, \(\mu\), assuming
one starts in the ground electronic state:
\(\left\langle u\left(t\right)u\left(0\right)\right\rangle_{g}=Tr\left[\rho_{\mathbf{0}}\left(t\right)\mu\right]\).
(15)
The initial density matrix is chosen as the ground state,\(\rho_{0}\left(0\right)=\left.\ |g\right\rangle\left\langle g|\right.\ \).