On stability issues of the HEOM method

The Hierarchical Equations of Motion (HEOM) method has become one of the cornerstones in the simulation of open quantum systems and their dynamics. It is commonly referred to as a non-perturbative method. Yet, there are certain instances, where the necessary truncation of the hierarchy of auxiliary density operators seems to introduce errors which are not fully controllable. We investigate the nature and causes of this type of critical error both in the case of pure decoherence, where exact results are available for comparison, and in the spin-boson system, a full system-reservoir model. We find that truncating the hierarchy to any finite size can be problematic for strong coupling to a dissipative reservoir, in particular when combined with an appreciable reservoir memory time.


Introduction
The dynamics of open quantum systems is a field in which has seen the development of mature methodologies and powerful abstractions such as completely positive maps [1], Lindblad generators [2], path integral methods [3] and stochastic approaches [4][5][6][7][8].All these approaches vary in their general applicability, in particular, to non-Markovian dynamics, as well as their potential as computationally feasible methods.Another computational method which has attracted increasing attention in recent years is the hierarchical equations of motion method (HEOM) [9][10][11][12][13], which extends the dynamics of the reduced density matrix by adding a (potentially large) number of auxiliary density operators (ADOs).Through these, one obtains a time-local version of the open-system dynamics without the usual construction of a dissipation superoperator.The latter approach typically requires application of Born, Markov and secular approximations, yielding a superoperator depending on both system and reservoir properties.The HEOM approach, on the other hand, offers a structure where the terms describing dissipation only depend on the system-reservoir coupling.This means that Hamiltonians with non-trivial level structure or complicated explicit time dependence can be adapted in the HEOM formalism far more easily than in master equations.This feature was leveraged in the application of HEOM to multi-dimensional spectroscopy [14].Recent developments greatly extending its utility in the regime of ultralow temperature [13] raise the expectation that this advantage will also be leveraged in the simulation of quantum devices.
All HEOM approaches rest on a multiexponential decomposition of the two-time correlation function of the quantum reservoir coupled to the system of interest.
In the case of a Gaussian reservoir, it provides all information needed to determine the reduced dynamics.In a strict analytic sense, all thermal (Matsubara) frequencies should be among the z j of such a decomposition.This suggests that the required number K of exponential terms will rise when lower temperatures are considered, eventually rendering the method infeasible in terms of computational cost.However, when requiring the sum on the right-hand side of eq. ( 1) only to be a close numerical approximation, the number of terms can be drastically lowered [13], allowing HEOM computations at arbitrary temperature.
In the resulting Free-Pole HEOM (FP-HEOM) approach, the coupled dynamics of the reduced density matrix ρ 00 (t) and ADOs ρ m,n (t) reads Here boldface indices m, n are multi-indices (m 1 , . . .m K ), and (n 1 , . . .n K ), and the superscripts "+" and "-" in conjunction with a subscript k mean that the k-th element of m or n is raised/lowered by 1.In FP-HEOM, both the coefficients d j and the decay rates z j can be complex numbers, however, the real part of each z j must be positive for a decaying C(t).
Theoretically, an infinite number of ADOs is required for any HEOM system to be exact, which is clearly impractical.Computationally, equations like (2) are truncated to systems of equations with a finite number of ADOs.Typically, it is tacitly assumed that increasing the number of ADOs will allow the result to converge.This seems to be the case in the majority of applications.However, cases are reported [15][16][17][18], where convergence seems to be very slow, and only pointwise with respect to time.In a few examples, the limit of large truncation level even appears to be divergent.The aim of the following sections is to develop a better understanding of the convergence of HEOM dynamics with respect to truncation depth and to point out potential remedies to convergence problems.
2 Pure dephasing in the hierarchy formalism

Analytic results
For pure dephasing, an analytic result in terms of quadratures of the reservoir correlation function is available, with real function (4) In the Markovian limit, φ(t) rises linearly with time, the slope being given by the spectral noise power at zero frequency.
Thus, the HEOM approach is not strictly needed in this context.However, this compact analytic result provides the opportunity to compare the dynamics resulting from the HEOM approach with analytic results.In the case of pure decoherence, a double multi-index is not needed; the resulting hierarchy equation is Initially, we consider the simple case of a real correlation function with exponential time dependence, C(t) = de −Γt .This can be realized as a limiting case of an ohmic quantum reservoir with very weak coupling and very high temperature, or by replacing the reservoir with classical noise (Ornstein-Ulhlenbeck process).Within this section, we will assume Γ > 0. Real C(τ ) and positive noise power require d to be positive as well, but we will relax this constraint here, since terms with negative or complex d typically appear for low-temperature quantum reservoirs, described by more than one exponential term.The corresponding HEOM system now has a single index, and we set ρ −1 ≡ 0. We consider the time evolution of off-diagonal matrix elements of the physical density matrix and the ADOs, setting y n (t) = q|ρ n (t)|q ′ , which results in with y −1 ≡ 0 and D = d(q − q ′ ) 2 .For a factorizing initial condition, y n (t = 0) = 0 for n > 0, this system has the closed-form solution with All of the y n (t) behave as ∝ e −Dt/Γ in the long-time limit.Their transient behavior, which becomes markedly visible for |D| ≫ Γ 2 , indicates the non-Markovian nature of dephasing through colored noise.In the opposite limit, the coupling between hierarchy levels is a small perturbation, and hierarchy elements with n > 1 are virtually negligible.
The right-hand side of the linear dynamical equation ( 7) defines a generator G equivalent to the Hamiltonian of an inverted quantum harmonic oscillator with a complex shift, Without any truncation, the eigenvalues of this generator can be obtained in complete analogy to the quantum mechanical oscillator; they are The analogy to the quantum oscillator is closest for the case D < 0, where G becomes hermitian with respect to the natural scalar product of sequence spaces.More generally, eq. ( 11) is valid for arbitrary complex parameters Γ and D, although the corresponding eigenvectors will generally be non-orthogonal.In particular, this is the case when both D and Γ are positive.

Truncated HEOM dynamics
The procedure of first obtaining exact solutions y n (t), then discussing which ones can be neglected is not equivalent to the approach which is taken when applying HEOM as a numerical algorithm: There, eq. ( 6) must be truncated, typically by setting ρ n = 0 for n larger than some n max .For a theoretical analysis, the truncation could also be imagined as a rank-one modification of the generator G which makes it block diagonal, one block with n ∈ [0, n max ] and one (semi-infinite) block with n > n max .It is to be noted that there are cases in which the rank-one modification mentioned above is not a small perturbation: The level spacing of G is Γ, and the modification eliminates matrix elements proportional to √ n max D. The potential discrepancies between exact and truncated dynamics are illustrated in Fig. 1, where the quantities |y n (t)| exp(+Dt/Γ) are shown for D = ±1.5Γ 2 and D = ±2.5Γ 2 .Significant deviations are observed for large positive D, where decoherence effects are underestimated in HEOM.The different level of error depending on the sign of D probably reflects the different nature of the eigensystem of G for different signs of D -orthogonal vs. non-orthogonal.
When using a multi-exponential representation of a low-temperature quantum correlation function C(t), one typically obtains one or more coefficients d with negative real part, leading to rising exponentials e |D|t/Γ in the resulting dynamics, which are compensated by other terms in the multi-exponential representation.Inaccuracies due to the hierarchy truncation can jeopardize this cancellation, as is occasionally observed in HEOM simulations [15][16][17].In order to reproduce this difficulty in a clean example, we consider the degenerate case of a formal two-exponential representation The exact solution for y 0 (t) is obviously constant, y 0 (t) = y 0 (0).Formally applying eq. ( 5) with K = 2, d 1,2 = ±d and z k = Γ results in a HEOM approximation of the exact result, the quality of which depends on the truncation parameter n max .The data in Fig. 2 shows that a high truncation parameter can ensure this with a small realtive error of approximately 8 • 10 −5 for a relatively large parameter D = 10Γ 2 (a), while further raising D to the value 20Γ 2 (b) breaks the cancellation, leading to an exponential divergence even for extremely deep hierarchies.
3 Stability issues in full HEOM dynamics

Time-domain analysis
In practical HEOM computations, two slightly different truncation schemes are used.In the socalled N truncation scheme [19], all ADOs where any index exceeds a given n max are set to zero.In the L truncation scheme [19], ADOs where the sum of all indices exceeds a given n max are neglected.For the dynamics governed by eq. ( 2), this yields n 2K max ADOs for N truncation and nmax−1 r=0

2K−1+r r
ADOs for L truncation, a significantly smaller number than the former.Converging the dynamics in a numerical simulation is simply a matter of increasing the truncation parameter n max in either of the truncation methods until the dynamics no longer changes.In the following we investigate how this procedure succeeds or fails in the context of a spin-boson system with strong coupling and a sluggish, roughly resonant reservoir with Hamiltonian Ĥs = ∆σ x and coupling operator q = σz .For the influence of the reservoir a ohmic spectral density J(ω) = αω/(1+(ω/ω c ) 2 ) with Drude form cutoff was used.For a given coupling parameter α, inverse temperature β and cutoff frequency ω c a barycentric fit [13,20] with tolerance 10 −2 was performed, yielding a correlation function with a single (K = 1) complex exponential mode C(t) = de −zt , which is more than adequate for a high-temperature reservoir.The exploration of large n max is thus facilitated.In Fig. 3 the dynamics of the population of the exited state can be observed for the strong dissipation case α = 5.The truncation parameter n max is steadily increased in an attempt to achieve physical dynamics for both the N-and L-Truncation.Within the numerically accessible parameter range, there appears to be pointwise convergence in the N-truncation up to a threshold near n max = 100, with convergence attainable only on a finite time interval.
The population dynamics in the L-truncation displays somewhat different convergence properties.Low truncation orders do not seem to be unstable, but due to the strong dissipation also not yet converged.An increase of the truncation parameter eventually causes divergences instead of reaching a physical result; these divergences appear to persist for arbitrarily deep hierarchies.

Spectral features of extended-state-propagation
Since eq. ( 2) is linear, HEOM instabilities can also be analyzed through the eigenvalues of the generator of the extended-state propagation.For this purpose, it is convenient to vectorize the set of ADOs ρm,n → |ρ m,n and to define a global state |ρ = (. . ., |ρ m,n , . . . ) t m,n∈T , where T is the set of all multi-indices (m, n) allowed by the truncation method.The FP-HEOM is now equivalently expressed by a linear map characterized by a Matrix M (t), d|ρ /dt = M (t)|ρ .M (t) is sparse, since it changes ADO indices at most by one.When it comes to investigating the long time stability, we will only consider the case of a time independent Hamiltonian, which allows the map to be expressed through a time-independent matrix M .From it is obvious that an exponential divergence occurs when M has an eigenvalue with positive real part.Without simulating the system in the time domain we can simply inspect for such eigenvalues.
Determining the eigenvalues can be done either analytically or numerically [16].We use The HEOM indices 0 . . .7 are color coded as dark blue, red, yellow, purple, green, light blue, cardinal and blue.Drawn lines represent HEOM result for finite truncation (n max = 10); dots represent the exact result (8).
the ARPACK library [21] for its efficient iterative algorithm in the search of eigenvalues with the largest real part for sparse matrices.For the numerical analysis of the instability we consider the same spin-boson system from the time-domain analysis.With stability thus defined, Fig. 4 shows a phase diagram with stable and unstable regions in the parameter plane spanned by n max and α, with different phase boundaries for L truncation and N truncation.In both cases, instabilities arise for α above the boundary.All α below the line will result in stable, but not necessarily converged dynamics.In the L-truncation the allowed maximum coupling initially decreases with the depth of the hierarchy n max and tends towards some constant.For the N-truncation it seems to be the other way around, where the maximum allowed damping increases linearly with n max until some threshold depth and becoming constant afterwards.From Fig. 4 it appears that for high hierarchy depths there is a critical coupling constant α c ≈ 3 which becomes independent of n max .Thus, merely increasing the hierarchy depth n max does not appear to be a viable strategy for curing divergence problems.
In the case of L truncation, comparably high values of α are allowed for smaller n max , which however, may be too small for numerical convergence of the dynamics.So in the worst case scenario, it might be impossible to converge the result without running into instabilities.This behaviour is consistent with the dynamics in Fig. 3. Fig. 5 shows the full eigenvalue spectrum of two specific points for n max = 50 in the N-truncation, namely α = 1 and α = 5.In the moderate damping case α = 1 both the N-and L-truncation are stable and converged, whereas the strong damping case α = 5 does not allow for stable dynamics.The spectrum for the stable case α = 1 shows the absence of eigenvalues with positive real part and a eigenvalue with Re(λ) = 0 corresponding to the  shows many eigenvalues with positive real part much greater than zero corresponding to unstable dynamics.For increasing n max , table 1 shows the largest real part of any eigenvalue of M .It appears that these tend to a constant, again suggesting that a convergent and correct result cannot be achieved.Different remedies have been suggested against the problem of divergent cases.Dunn et al. [16] proposed a filtering algorithm which eliminates from M the eigenvectors corresponding to unstable eigenvalues, with some success in the case  where the modification of M is of low rank.Considering the large number of eigenvalues with positive real part shown in Fig. 5d), with most real parts large compared to system parameters, the applicability of this procedure seems limited.Further development of truncation strategies approximating the omitted ADOs by a function of the retained ADOs [10,22] instead of setting them to zero might lead to progress in solving this problem, possibly even allowing the use of smaller n max than with a "hard" truncation.The reformulation of HEOM dynamics through a bath coordinate mapping [23] also seems a promising strategy in the case of a reservoir which is sluggish or near-resonant to the system.Finally, the combination of HEOM for fast reservoir components and an optimized stochastic approach for slow reservoir components [24] can provide a solution to the problem.

Conclusions
The appearance of divergent errors in some HEOM computations is succinctly related to reservoir parameters and truncation depth.The simple, analytically solvable case of pure dephasing reveals that a hard truncation at some maximum index is not necessarily a perturbative change of the dynamics, and it does not become perturbative when the maximum index is raised.A numerical study of the spin-boson system clearly confirms this picture.Problems of this nature seem mostly confined to the combined regime of strong coupling and slow reservoir dynamics.Apart from this regime, the HEOM method remains a valid and efficient method, as evidenced by the vast majority of applications where divergences are not observed.For the problematic parameter regime, our work should stimulate the development of alternative strategies for the finite closure of HEOM equations or a HEOM mode structure less susceptible to divergence problems.

Fig. 5 :
Fig. 5: Full spectrum (left) and magnification (right) of the eigenvalues of the linear FP-HEOM map.Figures a)and b) correspond to the stable example points in Fig.4, and c) and d) correspond to the unstable example points respectively.

Table 1 :
Most unstable eigenvalue as a function of truncation depth.See text for details.