Removing instabilities in the hierarchical equations of motion: exact and approximate projection approaches

TL;DR

This paper addresses temperature-dependent instabilities in the hierarchical equations of motion (HEOM) for quantum systems, proposing projection methods to stabilize and accurately approximate dynamics with reduced computational cost.

Contribution

The authors introduce projection techniques, including eigenmode projection and Prony filtration, to eliminate instabilities in HEOM and improve computational efficiency for quantum dynamics simulations.

Findings

Projection onto stable eigenmodes removes instabilities effectively.

Few-level projected dynamics closely match exact results at zero temperature.

Prony filtration offers a practical alternative when eigenmode computation is costly.

Abstract

The hierarchical equations of motion (HEOM) provide a numerically exact approach for computing the reduced dynamics of a quantum system linearly coupled to a bath. We have found that HEOM contains temperature-dependent instabilities that grow exponentially in time. In the case of continuous-bath models, these instabilities may be delayed to later times by increasing the hierarchy dimension; however, for systems coupled to discrete, non-dispersive modes, increasing the hierarchy dimension does little to alleviate the problem. We show that these instabilities can also be removed completely at a potentially much lower cost via projection onto the space of stable eigenmodes; furthermore, we find that for discrete-bath models at zero temperature, the remaining projected dynamics computed with few hierarchy levels are essentially identical to the exact dynamics that otherwise might require an intractably large number of hierarchy levels for convergence. Recognizing that computation of the eigenmodes might be prohibitive, e.g. for large or strongly-coupled models, we present a Prony filtration algorithm that may be useful as an alternative for accomplishing this projection when diagonalization is too costly. We present results demonstrating the efficacy of HEOM projected via diagonalization and Prony filtration. We also discuss issues associated with the nonnormality of HEOM.

Figures

Equations

$$\hat{H}$$

$${\hat{H}}_s$$

$${\hat{H}}_b$$

$${\hat{H}}_{sb}$$

$$\hat{V}$$

$${\hat{H}}_{sb}$$

$$\hat{V}$$

$$J(\omega )$$

$${\hat{H}}_s$$

$${\hat{H}}_b$$

$${\hat{H}}_{sb}$$

$${\hat{V}}_n$$

$${\hat{H}}_{sb}$$

$${\hat{V}}_n$$

$$\hat{\sigma }(t)$$

$$\frac{d}{dt}{\hat{\rho }}_{m_{1\pm },...,m_{N_b\pm }}(t)$$

$$-i\sum _{n=1}^{N_b}{\omega }_0\left(m_{n-}-m_{n+}\right){\hat{\rho }}_{m_{1\pm },...,m_{N_b\pm }}(t)$$

$$\displaystyle+\sum_{n=1}^{N_{b}}\biggr{[}\Phi_{n}\biggr{(}\hat{\rho}_{m_{1\pm},...,m_{n+}+1,...,m_{N_{b}\pm}}(t)$$

$$\displaystyle+\hat{\rho}_{m_{1\pm},...,m_{n-}+1,...,m_{N_{b}\pm}}(t)\biggr{)}$$

$$+m_{n+}{\Theta }_{n+}{\hat{\rho }}_{m_{1\pm },...,m_{n+}-1,...,m_{N_b\pm }}(t)$$

$$\displaystyle+m_{n-}\Theta_{n-}\hat{\rho}_{m_{1\pm},...,m_{n-}-1,...,m_{N_{b}\pm}}(t)\biggr{]},$$

$$\mathcal{L}$$

$${\Phi }_n$$

$${\Theta }_{n\pm }$$

$$\displaystyle\mp\{\hat{V}_{n}-\hat{V}_{n-1},...\}\biggr{)}.$$

$$\sum _{n+=1}^{N_b}{m}_{n+}+\sum _{n-=1}^{N_b}{m}_{n-}$$

$${\hat{\rho }}_{m_{1\pm },...,m_{N_b\pm }}(t)$$

$$\sum _{n+=1}^{N_b}{m}_{n+}+\sum _{n-=1}^{N_b}{m}_{n-}$$

$$\frac{d}{dt}{\hat{\rho }}_{m_0,...,m_K}(t)$$

$$\displaystyle+\sum_{k=0}^{K}\biggr{[}\Phi\hat{\rho}_{m_{0},...,m_{k}+1,...,m_{K}}(t)$$

$$\displaystyle+\sum_{k=0}^{K}m_{k}\Theta_{k}\hat{\rho}_{m_{0},...,m_{k}-1,...,m_{K}}(t)\biggr{]},$$

$$\mathcal{L}$$

Full text

††preprint: AIP/123-QED

I Introduction

A grand challenge in the physical sciences lies in modeling quantum dynamics in the condensed phase(Berne et al., 1998; May and Kühn, 2011; Weiss, 2012; Nitzan, 2006). Progress has often been made by using efficient approximate approaches that invoke a “system-bath” separation and treat the bath classically(Billing, 1975; Chen and Reichman, 2016) or the system-bath interactions perturbatively(Redfield, 1965). However, in systems where quantum effects in the bath play a pronounced role and where no small coupling or energy parameter may be identified, exact solutions to a fully quantum system-bath model are desired. One of the most successful computational methods for calculating exact quantum dynamics is provided by the hierarchical equations of motion (HEOM)(Tanimura and Kubo, 1989; Ishizaki and Tanimura, 2005; Ishizaki and Fleming, 2009; Chen et al., 2009; Tanimura, 2014). First derived by Tanimura and Kubo(Tanimura and Kubo, 1989), HEOM is a reformulation of the Feynman-Vernon influence functional approach to quantum dissipative dynamics(Feynman and Vernon, 1963; Tanimura and Kubo, 1989; Makri and Makarov, 1995a, b; Ishizaki and Tanimura, 2005; Ishizaki and Fleming, 2009; Liu et al., 2014; Chen et al., 2015), and its solution yields the exact reduced dynamics of linearly coupled system-bath models. HEOM has successfully addressed a variety of applications; yet, its scope is limited since the original formulation of HEOM requires the bath to be represented by a continuous spectral density. Such a constraint is often inappropriate for describing phenomena captured in venerable models of quasi-particle dynamics in organic molecular crystals(Spano, 2010; Bakulin et al., 2016; Tempelaar and Reichman, 2018; Fujihashi et al., 2017; Morrison and Herbert, 2017), coupled excitonic and vibrational motion in light harvesting complexes(Womick and Moran, 2011; Christensson et al., 2012; Tiwari et al., 2013; Tempelaar et al., 2014), and transport in polar crystals with narrow phonon bandwidths(Fröhlich, 1954; Feynman, 1955; Feynman et al., 1962; Devreese and Alexandrov, 2009). In these cases where discrete bath modes play a significant role, efficient and exact methods for computing quantum dynamics are desirable.

Motivated by the aforementioned computational challenges, several groups have recently explored novel formulations of HEOM that treat a discrete spectral density(Liu et al., 2014; Chen et al., 2015). There are several potential advantages in developing a “discrete-bath HEOM”. First, even in the discrete-bath formulation, HEOM retains the benefit that the reduced equations automatically include all possible bath excitations. As a result, HEOM eliminates the issue of basis-set convergence present in techniques such as exact diagonalization and matrix product states, and instead relies on convergence with respect to hierarchy depth (number of hierarchy levels, see II.2). Second, including discrete bath modes in HEOM serves to expand the ability of this powerful methodology to tackle an important class of problems, with the added benefit that several popular HEOM software packages such as Parallel Hierarchy Integrator(Strümpfer and Schulten, 2012) (PHI), pyrho(Berkelbach, ), and potentially also GPU-HEOM(Kreisbeck et al., 2011) can be readily adapted for the discrete-bath case.

Nonetheless, existing formulations of discrete-bath HEOM are not without serious problems that we will explore and partially remedy in this work. Recently, Chen et al.(Chen et al., 2015) published discrete-bath HEOM simulations which converge with few hierarchy levels to the exact reduced dynamics for a 10-site Holstein model. We find that upon extending the time of these simulations, the dynamics are eventually plagued by the abrupt onset of an exponential instability. This instability is illustrated in Fig. 1. Unfortunately, increasing the number of hierarchy levels does little to delay this instability to later times and comes with a great computational cost. In the interest of honing HEOM to be a useful tool for modeling the dynamics of discrete-bath models, in this work we explore such instabilities and discuss approaches for removing them to facilitate longer-time simulations without the need for a large hierarchy depth. We also show that similar instabilities exist in the original “continuous-bath HEOM” at low temperatures. Our main finding, which is explicit in the discrete-bath case and plausible in the continuous-bath case, is that one can remove these instabilities without altering the exact long-time dynamics.

In this work we proceed as follows. In II we introduce the condensed phase models which we study, as well as the HEOM that describe their dynamics. In III we introduce our spectral approach for studying the stability of HEOM and illustrate several stable and unstable examples. In IV, we then apply our stability analysis to both discrete-bath and continuous-bath HEOM at nonzero temperatures, providing a brief discussion of instabilities in low-temperature continuous-bath HEOM that has only been alluded to(Moix and Cao, 2013) in the vast HEOM literature. In V we discuss a diagonalization approach for exactly projecting out the instabilities in discrete-bath HEOM. In VI we present an iterative method for accomplishing this projection that does not require the full diagonalization of the HEOM. Finally, in VII we comment on some difficulties that may arise due to the nonnormality of HEOM.

II Theory

II.1 Model Hamiltonians

To demonstrate instabilities in HEOM we consider several standard system-bath models for open quantum dynamics. Each of these models takes the form

[TABLE]

where $\hat{H}_{s}$ describes the degrees of freedom (DOFs) of the reduced system, $\hat{H}_{b}$ describes the DOFs of the phonon bath, and $\hat{H}_{sb}$ describes the coupling between the system and bath DOFs. In following sections, we will refer to the delta-function spin-boson model (DSB), the Holstein model(Holstein, 1959a, b), and the Su-Schrieffer-Heeger (SSH) model(Su et al., 1979) as discrete-bath models. Likewise, we will refer to the continuous spin-boson model (CSB) as a continuous-bath model. $\hat{a}$ and $\hat{b}$ are electron and phonon annihilation operators, respectively. We work in dimensionless units and set $\hbar=1$.

II.1.1 Spin-Boson Model

In this work, we consider a spin-boson model(Weiss, 2012; Chen and Reichman, 2016) with no energy bias between sites. The system has two energy levels,

[TABLE]

where $J$ is the inter-site coupling constant. This two-level system interacts with a harmonic oscillator bath

[TABLE]

where $\omega_{j}$ is the frequency of the $j$-th bath mode. We consider two forms of this model. In the DSB there is only one bath oscillator, such that $N_{b}=1$, $\omega_{j}\equiv\omega_{0}$, and $b_{j}\equiv b$. For the DSB, the coupling is defined as

[TABLE]

where

[TABLE]

In the CSB an infinite number of bath oscillators are included. For the CSB, the coupling is defined as

[TABLE]

where

[TABLE]

The coefficients $c_{j}$ are fixed via the spectral density

[TABLE]

which we choose in this work to be the Debye spectral density

[TABLE]

In practice, the spin-boson model may be used as a coarse-grained description for any system that is well approximated by a two-level system coupled linearly to a harmonic bath. Applications of the spin-boson model are enumerated by Weiss and include the description of qubits, tunneling phenomena, and electron transfer processes(Weiss, 2012).

II.1.2 Holstein model

We also consider a one-dimensional Holstein model(Holstein, 1959a, b; Mahan, 2000) with periodic boundary conditions. The system is described by a tight-binding Hamiltonian

[TABLE]

interacting with a harmonic oscillator bath

[TABLE]

with site-diagonal coupling

[TABLE]

where

[TABLE]

The Holstein model, termed a “molecular crystal model,” was introduced to extend the conventional continuum treatment of polarons(Pekar, 1946; Fröhlich, 1954) to account for the deformation of a discrete lattice(Holstein, 1959a, b). It reflects the decoupled nature of sites in a molecular crystal by including only local electron-phonon coupling under the assumption of Einstein phonons. The Holstein model has the advantage that it can be used for a range of coupling strengths to describe large and small polarons, alike(Holstein, 1959a, b). For an excellent review that discusses the relation between the Holstein model and the Fröhlich model as well as the DSB, see Devreese and Alexandrov(Devreese and Alexandrov, 2009). In addition to modeling electron dynamics, the Holstein model has also seen great success in modeling Frenkel exciton dynamics in organic molecular crystals(Merrifield, 1964; Spano, 2010; Hestand and Spano, 2018).

II.1.3 Su-Schrieffer-Heeger model (SSH)

As an alternative to the Holstein model, we also briefly consider the SSH model(Su et al., 1979) with periodic boundary conditions, which differs from the Holstein model in its off-diagonal system-bath coupling

[TABLE]

where

[TABLE]

The SSH model was originally proposed to describe solitons in polyacetylene(Su et al., 1979), and has also been employed to model charge transport in crystalline organic semiconductors(De Filippis et al., 2015). The SSH coupling (14) accounts for the modulation of electron-hopping rates based on the variable nuclear distance between sites. When the SSH coupling (14) is combined with the Holstein coupling (12), the resulting model which includes both local deformation and phonon-mediated hopping is known as the Holstein-Peierls model(Munn and Silbey, 1985).

II.2 The hierarchical equations of motion

We now formally define HEOM, the exact quantum dynamics method of interest in this work. HEOM consists of a set of coupled linear differential equations that govern the time evolution of a hierarchy of indexed matrices. At the root of the hierarchy lies the reduced density matrix of the system,

[TABLE]

The dynamics of the DSB, Holstein model, and SSH model are described by the following HEOM(Liu et al., 2014; Chen et al., 2015).

[TABLE]

where

[TABLE]

and

[TABLE]

For the Holstein model and the DSB,

[TABLE]

while for the SSH model

[TABLE]

For the Holstein and SSH models, $N_{b}=N$. For the DSB, $N_{b}=1$ and $\hat{V}_{n}=\hat{V}$. Throughout this paper we define the inverse temperature $\beta=\left(k_{B}T\right)^{-1}$ and work in units where $k_{B}=1$.

The $l$-th hierarchy level consists of all matrices $\hat{\rho}_{m_{1\pm},...,m_{N_{b}\pm}}(t)$ in Eq. (17) for which

[TABLE]

In this study we truncate this infinite hierarchy of coupled differential equations after $L$ hierarchy levels with a “time-nonlocal” closure(Chen et al., 2009; Liu et al., 2014; Chen et al., 2015), where we set

[TABLE]

for

[TABLE]

In practice, solutions to (17) are to be converged with respect to the hierarchy depth $L$.

The CSB is described by the following HEOM(Ishizaki and Tanimura, 2005), which are similar in structure to Eq. (17) but produce markedly different dynamics due to the incorporation of an infinite bath.

[TABLE]

where

[TABLE]

and

[TABLE]

Eq. (25) incorporates an infinite Matsubara series, resulting from a high-temperature expansion, that is closed via truncation after $K$ Matsubara terms. Another popular closure for the Matsubara series in HEOM was derived by Ishizaki and Tanimura(Ishizaki and Tanimura, 2005). The Ishizaki-Tanimura closure approximately accounts for Matsubara terms with $k>K$ (for sufficiently large $K$) by replacing rapidly decaying factors of $\nu_{k}e^{-\nu_{k}t}$ with $\delta\left(t\right)$; this closure also has the added benefit of improved stability, although instabilities are still present. However, in the interest of using a continuous-bath HEOM that closely resembles the discrete-bath formulation in Eq. (17) we will not employ the Ishizaki-Tanimura closure in this work.

Here, the $l$-th hierarchy level consists of all matrices $\hat{\rho}_{m_{0},...,m_{K}}(t)$ in Eq. (17) for which

[TABLE]

Again, we use a time-nonlocal closure after $L$ hierarchy levels(Ishizaki and Tanimura, 2005; Chen et al., 2009) such that

[TABLE]

for

[TABLE]

In practice, solutions to (25) are to be converged with respect to the hierarchy depth $L$, as well as the number of Matsubara terms $K$.

Alternate closures exist for Eqs. (17) and (25) besides those shown in Eqs. (24) and (34). Continuous-bath HEOM studies regularly employ a closure that relies on the exponential suppression of deeper hierarchy levels(Tanimura and Wolynes, 1991; Ishizaki and Tanimura, 2005); we do not investigate this closure here since it is not applicable to discrete-bath models. Furthermore, the time-local closure(Xu et al., 2005; Chen et al., 2009) is applicable for both discrete-bath HEOM and continuous-bath HEOM; however, since it does not appear to suppress instabilities in discrete-bath HEOM and also is not amenable to the spectral analysis in III, we do not explore the time-local closure here.

For the initial condition employed in this work, similar to many other studies, we set the population of the first site $\rho_{0,...,0}^{00}(t=0)=1$, and all the other hierarchical matrix elements are set to zero.

III Spectral analysis

The HEOM presented in Eqs. (17) and (25) may be represented as linear systems of the form

[TABLE]

where $\vec{\rho}(t)$ is a vector containing all the hierarchy elements $\rho_{m_{1\pm},...,m_{N_{b}\pm}}^{ij}(t)$ or $\rho_{m_{0},...,m_{K}}^{ij}(t)$ and $A$ is a nonnormal matrix containing the coupling between all of the hierarchy elements. With this flattened representation of HEOM, spectral analysis can be used to study the stability of solutions. The solution to Eq. (35) may be written as

[TABLE]

Assuming $A$ to be a diagonalizable matrix and employing the eigen-decomposition of $A$, we can write Eq. (36) as

[TABLE]

where $\Lambda$ is a diagonal matrix containing the eigenvalues $\{\lambda_{i}\}$ of $A$ and $V$ is a matrix whose columns $\vec{v}_{i}$ are the normalized eigenvectors of $A$. We can write Eq. (37) in the eigenbasis of $A$ as

[TABLE]

where the components of $\vec{d}=V^{-1}\vec{\rho}(0)$ are the expansion coefficients of $\vec{\rho}(0)$ in the eigenbasis of $A$. From this representation, it is clear that terms in the sum with $\Re[\lambda_{i}]>0$ are asymptotically unstable, and will be referred to here as the unstable modes.

III.1 Unstable HEOM

In Fig. 2 and Fig. 3 we calculate the eigen-decomposition for the flattened representation of Eq. (17) and plot $\{\lambda_{i}\}$ for two discrete-bath models: the 2-site Holstein model and the 3-site SSH model. Likewise, in Fig. 4 we do the same for the CSB described by Eq. (25). Notice that eigenvalues are present in the right half-plane in all three cases: these eigenvalues correspond to unstable modes. Furthermore, computing $\vec{d}$, the decomposition of $\vec{\rho}(0)$ in the eigenbasis of $A$, reveals that some of these unstable modes have nonzero weights in the initial condition, leading to asymptotic instability in the dynamics. In the following sections we will interpret these unstable modes and discuss computational strategies for removing them.

While we have shown in Fig. 4 a particularly unstable example of continuous-bath HEOM, in many practical continuous-bath cases one can suppress any instabilities by converging with respect to $L$ and $K$ and employing the Ishizaki-Tanimura closure(Ishizaki and Tanimura, 2005) for the Matsubara series. To the contrary, instabilities are much harder to suppress via convergence with respect to $L$ in discrete-bath HEOM.

III.2 Asymptotically Stable HEOM

To show an example of asymptotically stable HEOM, we will now contrast the former unstable examples with the HEOM for the CSB (25) at a temperature where no unstable modes are present. In Fig. 5 we again show the spectrum of the flattened representation of Eq. (25)

(the HEOM for the CSB), only this time at a higher temperature. Since the high-temperature $K=0$ approximation for continuous-bath HEOM is known to be equivalent to the Zusman equation(Shi et al., 2009), we note the resemblance here to the eigentree structure from the Zusman equation spectral analysis reported by Jung et al.(Jung et al., 1999). It is clear that all eigenvalues are confined to the left half-plane of Fig. 5. As a result, the corresponding dynamics are asymptotically stable.

At this point, we would be remiss to not acknowledge the structure and symmetry present in the spectral plots of Figs. 2, 3, 4, and 5. While we omit a discussion of the spectral dependence on hierarchy depth, number of sites, and choice of model, we exemplify these dependences in the animations shown in the supplementary material.

IV Temperature Dependence of HEOM Spectra

The spectra of both discrete-bath HEOM and continuous-bath HEOM admit a rich temperature dependence. In Figs. 6, 7, and 8

we plot for the DSB and the CSB the real part of the most unstable eigenvalue, $\max_{i}\Re\left[\lambda_{i}\right]$, as a function of temperature. Note the similarity between the qualitative behavior in Fig. 6 and Figs. 7 and 8. Both discrete-bath and continuous-bath HEOM reveal unstable regions, i.e. temperature ranges where $\max_{i}\Re\left[\lambda_{i}\right]>0$, intermingled with stable regions. Furthermore, both have unstable regions concentrated at lower temperatures.

While there are gross similarities between these cases, there are several important differences of note. First, the behavior shown in Fig. 6 seems to be piecewise continuous while that of Figs. 7 and 8 contains many asymptotes. A simple analytically solvable example that helps us rationalize the appearance of these asymptotes for the CSB will be discussed in the Appendix. Second, while the discrete-bath HEOM instabilities in Fig. 6 are governed by a single convergence parameter $L$, the behavior appears quite complicated: note how the unstable regions merge at lower temperatures and tend to grow more unstable as $L$ is increased. Since the unstable regions shift as a function of hierarchy depth, the behavior of the HEOM solution may be erratic as the hierarchy depth is varied for a fixed temperature, since instabilities may appear and disappear and vary in severity at any fixed temperature. In contrast, for continuous-bath HEOM we see two distinct types of behavior governed by the parameters $K$ and $L$. In Fig. 7 we see that as $K$ is incremented, unstable asymptotes are annihilated one at a time without changing the temperatures of the remaining asymptotes, effectively lowering the upper bound on temperatures at which instabilities become problematic. In Fig. 8 we observe that as $L$ increases, the unstable regions simply change in shape and the temperatures at which asymptotes occur are invariant. Thus, we find that for continuous-bath HEOM** $K$** is predominantly responsible for controlling the most severe instabilities; an increase in** $L$** alone cannot remove the instability.

Next, we turn to the low-temperature behavior in continuous-bath and discrete-bath HEOM as it relates to the aforementioned instabilities. In continuous-bath HEOM, low temperatures are manifestly problematic since the HEOM are derived using a high-temperature Matsubara expansion. At $T=0$ continuous-bath HEOM as expressed in Eq. (25) is ill-defined due to the explicit factor of $\cot\left(\beta\gamma/2\right)$ that appears in Eq. 28. For small non-zero temperatures the situation is still problematic; while the HEOM are defined at temperatures between the asymptotes of $\cot\left(\beta\gamma/2\right)$, due to the high density of asymptotes (per unit temperature) at low temperature it is necessary to use a large $K$ to annihilate the offending asymptotes and obtain converged dynamics. The situation is quite different in discrete-bath HEOM, where there is no high-temperature expansion and no notion of Matsubara convergence. It would seem therefore that discrete-bath HEOM should be amenable to facile low temperatures simulations as claimed by Chen et al.(Chen et al., 2015). While the discrete-bath HEOM are indeed well-defined at all temperatures, they do not eliminate the issue of instabilities at low temperature. Specifically, the low-temperature instability that in continuous-bath HEOM is controlled by the number of Matsubara terms reappears in discrete-bath HEOM as an instability that is controlled by the hierarchy depth. Now that we have investigated the temperature-dependence of the instabilities, we focus our attention on methods for overcoming instabilities in discrete-bath HEOM at $T=0$.

V Projecting away instabilities exactly

Let us momentarily abandon the goal of extracting a physical $\vec{\rho}(t)$ from the unstable HEOM solution and instead only focus on obtaining a stable $\vec{\rho}(t)$. This task can be accomplished by projecting out the unstable modes from the dynamics as follows. Consider the diagonal matrix

[TABLE]

We construct the projected time evolution matrix $\bar{T}(t)$ by zeroing out any elements $e^{\lambda_{i}t}$ of $T(t)$ that are larger than unity in modulus. Then the following projected HEOM solution will be asymptotically stable,

[TABLE]

Although there is no guarantee that this projection $\vec{\rho}(t)\rightarrow\bar{\vec{\rho}}(t)$ will not unrecognizably alter the dynamics, it turns out that in cases we have studied for small and intermediate system-bath coupling in discrete-bath models, $\bar{\vec{\rho}}(t)$ rapidly converges to the exact $\vec{\rho}(t)$ as the hierarchy depth is increased. This is unmistakably evident in Fig. 9 , where for the 3-site Holstein model the projection transforms a severely unstable $L=8$ trajectory into what quantitatively resembles the exact dynamics as computed by diagonalization of the Hamiltonian. Similar success is observed for the 2-site Holstein model and the DSB. Since the net effect of unstable modes disappears from the dynamics as $L$ is increased, we interpret these unstable modes as a spurious, unphysical consequence of hierarchy truncation. Meanwhile, since the stable modes with $\Re\left[\lambda_{i}\right]\leq 0$ approach the exact dynamics as $L$ is increased, we interpret these stable modes as corresponding to the physical dynamics, which justifies our projection scheme as a method for obtaining stable, exact dynamics from discrete-bath HEOM.

VI Projecting away instabilities iteratively with Prony filtration

The matrix $A$ in Eq. (35) contains $\tilde{N}^{2}=\left(N^{2}\sum_{n=0}^{L}\frac{(n+2N_{b})!}{n!\left(2N_{b}\right)!}\right)^{2}$ elements(Shi et al., 2009) for the discrete-bath case. Therefore, the $O(\tilde{N}^{3})$ diagonalization-based projection algorithm proposed in V is completely intractable for all but the smallest systems, and only then with sufficiently weak system-bath coupling due to the increased hierarchy depth required to treat stronger system-bath coupling. As such there is a need for approximate or iterative computational techniques for removing unstable modes without requiring an explicit computation of all eigenmodes. Here we discuss one such algorithm. Since (35) is a first order linear differential equation, each hierarchy element may be written as a sum,

[TABLE]

where $c_{ij}$ and $\lambda_{j}$ are complex numbers. The projection in V is equivalent to subtracting off all terms in the sum for which $\Re\left[\lambda_{j}\right]>0$. Consider the following algorithm for subtracting these terms approximately:

Numerically integrate the HEOM using an explicit time-stepping algorithm such as fourth-order Runge-Kutta until a time $t_{2}$ when $\max_{i}|\rho_{i}(t_{2})|\gg 1$. It is necessary that by $t_{2}$ one or a small number of modes $\{\vec{F}_{j}\}\subset\{\vec{f}_{j}\}$ have grown many orders of magnitude larger than the other modes. 2. 2.

For each hierarchy element, approximate the weights $c_{ij}=\left(\vec{c}_{j}\right)_{i}$ and complex frequencies $\lambda_{j}$ of the dominant modes $\{\vec{F}_{j}\}$ using Beylkin and Monzón’s approximate Prony method(Beylkin and Monzon, 2005). This algorithm for fitting a function to a sum of complex exponential functions is described in detail in section 4 of Beylkin and Monzón(Beylkin and Monzon, 2005). 3. 3.

Choose a time $t_{1}<t_{2}$ such that $\max_{i,j}|c_{ij}e^{\lambda_{j}t_{1}}|\ll 1$. Subtract the unstable modes from the hierarchy at time $t_{1}$:

[TABLE] 4. 4.

Resume numerical integration from time $t_{1}$ until the next instability occurs.

In spirit, this algorithm is similar to excited state methods in quantum mechanics that project out low energy states by imaginary time propagation(Blume et al., 1997). In Fig. 10 we depict a single iteration of the algorithm above. Iterating this “Prony-filtering” algorithm, one can piece together the same projected HEOM solution that would have been provided via a single diagonalization using the approach in V. Such a trajectory is illustrated in Fig. 11 for the DSB; in Fig. 12 we show the corresponding spectrum of $A$, and in Fig. 13 we show the essentially perfect agreement between the diagonalization and Prony-filtering approaches for calculating stable projected dynamics.

The Prony filtering approach certainly has a significant prefactor due to the cost of Runge-Kutta time-stepping and the approximate Prony analysis of step 2. However, in a similar spirit to the power iteration method or Krylov methods for calculating dominant eigenvectors, this filtration approach has the advantage of not requiring an $O(\tilde{N}^{3})$ diagonalization to compute all of the eigenvectors of $A$. Therefore, the filtering algorithm holds the promise of better scalability compared to the diagonalization-based projection algorithm.

VII Nonnormality of HEOM

Many of the matrices typically encountered in quantum mechanics are Hermitian. Hermitian matrices are diagonalizable and have the properties that all eigenvalues are real, and that the eigenvectors form an orthonormal set, i.e. the matrix is normal. We have already seen that the realness of the eigenvalues is violated for the generator of hierarchy evolution $A$. It turns out that the latter property is also violated.

The nonorthogonality of the eigenmodes in HEOM suggests that removal of instabilities may be ill-conditioned, especially when using deeper hierarchies as is necessary to treat strong coupling. In Fig. 14

we show how the numerical error in computing the eigen-decomposition $A-VDV^{-1}$ grows as a function of hierarchy depth. This numerical error is not surprising given the nonnormality of $A$, and it explicitly indicates that the diagonalization-based projection in V may be not only expensive, but also ill-conditioned. We also show in Fig. 14 how the condition number of $V$ grows with hierarchy depth. The condition number is defined as

[TABLE]

where the norm $||V||_{2}$ is the largest singular value of $V$. This metric exposes how the linear dependence of the eigenmodes of $A$ grows with hierarchy depth, and it illustrates how their removal via either diagonalization-based projection or Prony filtration may be difficult. For calculating the metrics in Fig. 14 we have employed hierarchy scaling(Shi et al., 2009) with the aim of reducing the diagonalization error and condition number. Although the error is smaller than if we were to use unscaled HEOM, clearly this scaling does not eliminate the growth with hierarchy depth demonstrated in Fig. 14.

In addition to the existence of asymptotic instabilities governed by the spectrum of $A$, it is also possible that numerical instabilities occur at shorter times due to large transient behavior which is characteristic of nonnormal dynamics. A classic example of such numerically unstable transient behavior, as well as its pseudospectral analysis, can be found in a control theory study of Boeing 767 aircraft(Burke et al., 2003; Trefethen and Embree, 2005). We will leave for future work an investigation of whether the HEOM instabilities witnessed for larger Holstein models are caused by such numerically unstable transients or whether they are simply due to the asymptotic growth of the unstable eigenmodes. It is likely that hierarchy scaling(Shi et al., 2009) substantially reduces such transients by scaling down the magnitude of elements deep in the hierarchy.

VIII Concluding remarks

While HEOM is a powerful method that often converges quickly to the numerically exact dynamics over a significant time range, we have shown evidence that HEOM trajectories for both continuous-bath and discrete-bath models at sufficiently low temperature will eventually hit an exponential wall of instability that completely corrupts the description of the time evolution. While this instability can typically be converged away in continuous-bath HEOM, we find that in discrete-bath HEOM deepening the hierarchy does little to delay instabilities, such that novel projection schemes are desired. Two methods, direct and iterative, have been presented to project out the instabilities, and for discrete-bath HEOM it has been shown that the remaining projected solution converges to the exact dynamics without requiring many hierarchy levels. We have discussed challenges that may arise associated with the computational cost and increasing nonnormality of larger and more complex HEOM simulations. As of now, we still lack a complete analytical understanding of the properties of HEOM that lead to the instabilities. We also fall short of a generic scalable solution for removing these instabilities. Perhaps with the advent of new HEOM algorithms such as distributed memory HEOM(Kramer et al., 2018) and matrix product state compressed HEOM(Shi et al., 2018), one may be able to accelerate numerical integration of the HEOM sufficiently to facilitate projection approaches such as the Prony filtration approach introduced here. The challenge of obtaining efficient, stable HEOM solutions will surely benefit from future work that explores alternative closures to the HEOM which reduce the instabilities without corrupting the remaining dynamics, the relation between the breaking of positivity in HEOM(Witt et al., 2017) – as evidenced by the negative populations in this work – and the instabilities, the nature of instabilities in other novel HEOM formulations(Tang et al., 2015; Duan et al., 2017; Nakamura and Tanimura, 2018; Erpenbeck et al., 2018), and computational techniques for removing the unstable modes from nonnormal, unstable linear systems.

Supplementary Material

See supplementary material for animations depicting the spectral dependence of HEOM on hierarchy depth, number of sites, and choice of model.

Acknowledgements.

The authors thank Gregory Beylkin, Seogjoo Jang, Ramin Khajeh, Benedikt Kloss, Matthew Reuter, and Qiang Shi for helpful and enlightening discussions. I.S.D. acknowledges support from the United States Department of Energy through the Computational Sciences Graduate Fellowship (DOE CSGF) under grant number: DE-FG02-97ER25308. D.R.R. acknowledges funding from NSF Grant No. CHE-1839464.

Appendix - Analytical Treatment of CSB

In this appendix we will use a small analytical example to demonstrate how unstable modes arise in HEOM for the CSB. Consider the HEOM time evolution operator for the CSB with $K=0,L=2$,

[TABLE]

If we consider the case $J=0$, we can get closed form expressions for the eigenvalues of $A$:

[TABLE]

Consider further the case where $\lambda>0,\gamma>0$. The last of these eigenvalues gives rise to an unstable mode whenever

[TABLE]

This condition is equivalent to

[TABLE]

This analytical example reveals alternating temperature regions of stability and instability for the CSB, akin to those demonstrated in Figs. 7 and 8, with boundaries located at

[TABLE]

The asymptotes in the cotangent function at $\beta=2n\pi/\gamma$ correspond to an infinite eigenvalue; at such temperatures the HEOM specified by Eq. (53) are completely undefined. Furthermore, for $\beta$ just slightly less than $2n\pi/\gamma$, these HEOM will be exceptionally unstable due to the large magnitude of the last eigenvalue in Eq. (54) near the asymptotes of the cotangent function. On the contrary, for $n\pi<\beta\gamma<\left(n+\frac{1}{2}\right)\pi$ we see that these HEOM are asymptotically stable. These findings agree with the spectral data shown in Figs. 7 and 8.

We should also note that these asymptotes are not present in Fig. 6, since the HEOM in Eq. (17) contain $\coth\left(\beta\omega_{0}/2\right)$ rather than $\cot\left(\beta\gamma/2\right)$ and the hyperbolic cotangent does not contain asymptotes at finite temperatures. Thus, unlike the behavior exhibited in Figs. 7 and 8, that of Fig. 6 is piecewise continuous.