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