About the performance of perturbative treatments of the spin-boson dynamics within the hierarchical equations of motion approach

The hierarchical equations of motion (HEOM) provide a numerically exact approach for simulating the dynamics of open quantum systems coupled to a harmonic bath. However, its applicability has traditionally been limited to specific spectral forms and relatively high temperatures. Recently, an extended version called Free-Pole HEOM (FP-HEOM) has been developed to overcome these limitations. In this study, we demonstrate that the FP-HEOM method can be systematically employed to investigate higher-order master equations by truncating the FP-HEOM hierarchy at a desired tier. We focus on the challenging scenario of the spin-boson problem with a sub-Ohmic spectral distribution at zero temperature and analyze the performance of the corresponding master equations. Furthermore, we compare the memory kernel for population dynamics obtained from the exact FP-HEOM dynamics with that of the approximate NIBA (Non-Interacting-Blip Approximation).


Introduction
The second-order quantum master equation (QME) is a fundamental tool for studying the dynamics of open quantum systems and finds wide-ranging applications in diverse fields, including quantum optics, condensed matter physics, condensed phase chemistry, and biology [1][2][3][4][5][6].However, employing the QME in regimes with strong non-Markovian behavior, characterized by slow environmental evolution and significant retardation effects on the subsystem, presents substantial challenges.These challenges become even more pronounced in low-temperature conditions and structured environments, where non-Markovian effects are amplified.A notable example is the unconventional impact of sub-Ohmic noise at zero temperature [7][8][9].
To address the limitations of the QME in handling pronounced non-Markovian retardation effects, incorporating higher-order corrections that account for systembath correlations becomes essential [10][11][12].However, direct calculations involving these terms require computing high-dimensional time integrals, which is a daunting task [13,14].To the best of our knowledge, analytical calculations of high-order perturbations or kernels have only been achieved up to the sixth order [15,16].Previous work has employed numerical techniques such as hierarchical equations of motion (HEOM) to compute higher-order effects [11,12].However, these approaches have been primarily limited to elevated-temperature regimes due to the exponential increase in the number of auxiliary density operators (ADOs) involved, a challenge commonly referred to as the curse of dimensionality.This challenge becomes even more severe as the temperature approaches zero, making the analysis of memory kernels' properties and high-order perturbative master equations increasingly demanding.
To investigate the effects of high-order kernels in unconventional environments, we address key challenges such as zero-temperature calculations employing the widely used HEOM method [9,11,.The HEOM approach utilizes a set of ADOs to unravel system-bath correlations within an extended state space [24,30,31,39,43].Conventionally, these ADOs are constructed as high-dimensional arrays based on a series of exponential functions derived from the decomposition of the bath correlation function.The extended space, expanded by ADOs, grows exponentially with the number of modes, posing significant challenges in scenarios involving low temperatures and structured bath spectra.This is particularly relevant as standard analytical Matsubara frequencies (i.e., ν k = 2kπ/β) shift towards a continuum strip as temperatures decrease (β = 1/T ).Therefore, developing an efficient method for obtaining a limited number of well-behaved modes becomes crucial for successfully implementing HEOM in the study of quantum master equations.
Several decomposition schemes [9,41,44] have been developed to address the limitations of HEOM in low-temperature and general spectral density scenarios, such as the frequency-domain barycentric spectrum decomposition (BSD) [9].The effectiveness of these decomposition schemes has been demonstrated by their ability to reproduce phenomena like the Kondo resonance peak in Fermi baths [44] and the Shiba relation in sub-Ohmic boson baths [9].Notably, by optimizing the properties of auxiliary modes under specific constraints, these schemes allow for the integration of a minimal number of modes (poles) into various versions of HEOM and other methods, such as the pseudomode approach [45], the hierarchy of pure states (HOPS) [46], and the nonequilibrium dynamical mean-field theory (DMFT) [47], even in unconventional environments.
In this work, we employ the Free-Poles HEOM (FP-HEOM) [9], which combines the optimized BSD and a hierarchical structure, to showcase the performance of higherorder master equations in addressing unconventional environments.We also compare it with the well-established Non-Interacting-Blip Approximation (NIBA), a powerful but approximate scheme for non-perturbative treatment of the system-bath coupling in spin systems.The remainder of this paper is organized as follows: In Sections II and III, we provide a brief overview of the hierarchical equations of motion and discuss their limitations in simulating dynamics of open quantum systems.In Section IV, we present numerical simulation results to demonstrate the efficiency and accuracy of our proposed method.Finally, we summarize our findings and discuss potential avenues for future research in the concluding section.

Model Hamiltonian
The spin-boson model, which describes a two-state system interacting bilinearly with a harmonic bath, has been extensively studied in open quantum system dynamics [7].In this work, we employ this model as a prototypical example to benchmark our method, while generalizations to more complex models are straightforward.
We consider a spin-boson model with the Hamiltonian where the system degrees of freedom (DoF) are represented by dimensionless Pauli matrices σ x and σ z .∆ and 2ǫ denote the tunneling energy and energy bias between the two eigenstates of σ z , i.e., |± .The ith harmonic modes of the reservoir is characterized by its mass, coordinate, momentum, and frequency, i.e. m i , x i , p i , and ω i , respectively.The coupling between the system and the ith bath mode is denoted by c i with [mω 2 x] dimension.
The effective impact of the bath on the system is fully described by the couplingweighted spectral density In the following, we assume a generic spectral density [48] of the form where α is the dimensionless Kondo parameter characterizing the system-bath dissipation strength and ω c is the characteristic frequency of the bath.For s = 1 this spectral distributions describes an ohmic reservoir and it is calledf sub-ohmic for s < 1.
In the latter situation, the interplay of the relatively large portion of low frequency modes with the internal two level dynamics makes this model extremely challenging to simulate, particularly in the long time limit and close to or at T = 0.

Hierarchical equations of motion
We briefly describe the essence of the HEOM approach and refer to the literature for further details.The derivation assumes a factorized initial states of the total density at time zero, i.e.W (0) = ρ s (0) ⊗ e −βH b /Tr e −βH b .The generalization to correlated initial states can be found in Refs.[49,50].
In path integral representation [51], a formally exact expression for the reduced density is obtained by integrating over the bath degrees of freedom.The effective impact of the bath onto the system dynamics is captured by the Feynman-Vernon influence functional [51] which reads where q ± (τ ) denote forward and backward paths with respect to the eigenstates |± .The influence functional describes arbitrary long-ranged self-interactions in time of the system paths determined by the bath auto-correlation C(t) = ξ(t)ξ(0) with the collective degree of freedom ξ = i c i x i .Particularly at low temperatures correlation functions are known to decay algebraically so that direct simulations of the path integral, e.g., via path integral Monte-Carlo (PIMC) [8,52] are plagued by a degrading signal to noise ratio.
The hierarchical equations of motion (HEOM) approach deals with this problem by converting the time-non-local path integral into a set of time-local differential equations through the introduction of auxiliary density operators (ADOs).The basic ingredient is a decomposition of C(t) in a series of K exponential modes, where K has to be sufficiently bounded for the nested hierarchy to be practically solvable.To formulate a minimal set of modes is thus the prerequisite to access low temperatures and long times.As a matter of fact, this problem has been solved by us only recently [9] by implementing the barycentric representation of rational functions for the spectral noise power S β (ω), i.e. the Fourier transform of C(t) with complex-valued amplitudes d k and frequencies Here, S β (−ω and S β (ω) are related by the fluctuation dissipation theorem which can be cast into the form The barycentric representation constructs a function Sβ (z) in the complex plane such that it is along the real axis an approximant of S β (ω) (AAA algorithm, see [53]).The poles and the residues of Sβ appear as {d k } and {z k } in (5).For details of this Free Pole-HEOM see [9].One can show that the FP-HEOM thus operates with a minimal set of K modes with the correction δC(t) below a chosen threshold.
The exponential function's self-derivative property in (5) then allows to 'unravel' the time non-locality of the original Feynman Vernon path integral by introducing ADOs ρ m,n (t) with multi-indices n = (n 1 , . . ., n K ) and m = (m 1 , . . ., m K ).This then leads to a nested hierarchy of time-local evolution equations for the ADOs, i.e., ρm,n = −(iL s The subscript m ± k and n ± k denote {m 1 , . . ., m k ± 1, . . .m K } and {n 1 , . . ., n k ± 1, . . ., n K }, respectively.The bare system evolution is propagated by Eventually, the physical reduced density matrix ρs corresponds to the multi-index m = n = 0. To boost the numerical efficiency of the FP-HEOM (6), matrix product state (MPS) representations can be conveniently implemented which allows to tackle also asymptotic times down to temperature zero [9].
Note that the truncation of the nested hierarchy of ADOs after tier L such that formally implies to include system-bath coupling strength up to order α 2L [54].In fact, as we have shown previously [55,56], a truncation at L = 1 leads to a generalized Redfield equation.

Numerical Results
In this section, we utilize the nonperturbative Free-Pole HEOM (FP-HEOM) method to investigate the dynamics of a two-level system coupled to a sub-Ohmic bosonic bath at zero temperature (T = 0), characterized by a spectral density of the form J(ω) ∝ ω s , 0 ≤ s ≤ 1.Our analysis focuses on two main aspects: (i) the convergence properties of higher-order master equations with respect to the system-bath coupling strength α and spectral exponent s, and (ii) the nonperturbative memory kernels, which serve as effective generators for the system dynamics.
To ensure the accuracy of our calculations, we employ a barycentric decomposition with a tolerance threshold of δC(t) ≤ 10 −3 .The FP-HEOM is propagated using the time-dependent variational principle (TDVP) [57] with a matrix product states (MPS) representation, utilizing a maximum bond dimension of χ = 30.The simulations are performed on a single Intel Xeon Gold 6252 CPU @ 2.1 GHz core, with the computation time ranging from a few minutes to several hours, depending on the specific scenario under investigation.

Performance of Truncated Time Evolution Equations within FP-HEOM
Various approximate treatments for open quantum dynamics have been developed, the most prominent ones based on considering the system-reservoir coupling up to second order.When a time scale separation exists between the rapidly decaying reservoir correlations and the relaxation dynamics of the reduced density operator, these methods lead to the Redfield master equation [58].As mentioned above, an extended Redfield equation, referred to as Redfield-plus, appears naturally if the FP-HEIOM is limited to ADOs with k (n k + m k ) ≤ 1. Namely, in the interaction picture, the ADOs of the zeroth tier (reduced density matrix) read with the first-tier ADOs determined via Here, the index I indicates the interaction picture.It is worth noting that the structure of Eqs. ( 8) and ( 9) resembles the one known from generalized Floquet theory for driven systems [59][60][61].The factorized initial state gives a boundary condition with ρ0,0 (0) = ρ(0) and all ADOs with L > 0 are set to zero.The inhomogeneous differential equations in Eqs. ( 9) can be formally solved ρI and then plugged in into Eq.( 8).Summation over the reservoir modes as in Eq. ( 5) then leads to the following time-evolution equation in Born approximation Note that no Markovian coarse-graining in time is done here, not even on the level of the reduced density so that the reservoir induced retardation is fully taken into account in this order of system-bath coupling.One regains the conventional Redfield equation by setting ρ I (τ ) → ρ I (t), thus leading to a time-local evolution equation with time-dependent rates.Hence, we name the above integro-differential equation ( 11) Redfield-plus.Thus, the FP-HEOM can be viewed as an infinite-order extension of the Redfieldplus/Redfield approximation [54,62], where the truncation at tier L yields an evolution equation non-local in time of order α 2L .This enables a systematic analysis of the role of higher-order system-reservoir correlations, which are particularly subtle within perturbative theory.
As demonstrated in Fig. 1 with a coupling strength of α = 0.1, numerically converged "exact" results can be achieved when the truncation tier reaches L = 12.However, the spin dynamics can be approximated already quite decently in the 6th order (L = 3), while the Redfield-plus (L = 1) approximation has the tendency to diverge over longer timescales.We remark, that while the FP-HEOM can treat also strong system-bath couplings, here, we have chosen a sufficiently weak parameter α to be in line with the concept of a perturbation series.

Redfield-plus: Influence of Low-Frequency Modes
For a sub-ohmic bosonic bath, the mode distribution characterized by the parameter s plays a significant role in the feasibility of a perturbative expansion according to Redfield-plus, as we show here.By fixing α = 0.05, in Fig. 2, we observe that the Redfield-plus can nearly reproduce the exact dynamics only for exponents 0.5 s.However, as the weight of low-frequency modes increases with lowering s, the second order approximation becomes inadequate already on gradually shorter time scales.

Time-depenmdemnt Memory Kernel for Populations:
NIBA versus FP-HEOM For a spin system with no bias, i.e. ǫ = 0, immersed in a bosonic bath, the expectation is that the steady-state population is uniformly distributed between the two spin states, i.e., P (t → ∞) = σ z (t → ∞) → 0. However, it is well-known that at zero temperature a symmetry breaking can take place corresponding to a quantum phase transition from a delocalized to localized asymptotic state [8,9].Consequently, depending on the initial condition and reservoir parameters, P (t) = 0 over long timescales, which is a purely quantum phenomenon due to the absence of thermal fluctuations at zero temperature.This behavior can be conveniently analyzed using the time-dependent memory kernel K(t) which determines the spin dynamics according to We mention here that a powerful perturbative treatment to derive the kernel is the so-called Non-interacting Blip Approximation (NIBA) [7,63] and its extensions.There, kernels in powers of the tunnel splitting ∆ are obtained with the NIBA kernel being of second order in the tunnel splitting ∆ 2 .Accordingly, the NIBA is not based on a series expansion in α and thus accounts also for strong spin-bath coupling.It neglects long-range quantum coherences though and does not predict the correct equilibrium state for ǫ = 0.In Fig. 3 results for the NIBA memory kernel are depicted.From these data one would conclude a change in the dynamical behavior (monotonous decay versus oscillatory decay) to occur for values around s = 0.5.For smaller exponents strong oscillations emerge with decreasing s.In contrast, exact FP-HEOM data extracted from the exact dynamics [12] are shown in Fig. 4. For the chosen coupling strength, it is also observed that for s > 0.5, the memory decays clearly monotonous over time and remains always positive as also predicted by NIBA.However, quantitatively the exact decay appears to be much faster than within NIBA.For values of s below this threshold oscillatory pattern are seen as well for the FP-HEOM results, but with much smoother and with less oscillations as the NIBA prediction.We can thus conclude that while the NIBA provides qualitatively the correct physics, it quantitatively in this regime of paramneter space not reliable.
Physically, the changeover in dynamical behavior in K(t) can be attributed to a freezing of the population since the total weight of the kernel tends towards zero when integrated over a time span where P σ (t) does not change considerably.More specifically, in this regime of sufficiently small s, one may write for long times with P σ (τ ) ≈ P σ (t) in Eq. ( 12) that with rates Estimations from Fig. 4 allow to see that rates k σ,σ ′ tend indeed to zero corresponding to a slowing down of the relaxation dynamics up to the regime, where k σ,σ ′ = 0 so that the spin requires an infinitely long time to reach a Gibbs equilibrium state and instead displays, localization.This property aligns with the fact that hybridization between spin and reservoir induced by slow modes occurs, thus freezing the spin dynamics.

Conclusion
Since conventional second-order master equations fail in accurately capturing non-Markovian dynamics, this limitation becomes markedly significant for unconventional baths, particularly in conditions such as zero-temperature structured environments.One approach to circumvent this is to incorporate high-order corrections.However, the numerical computation of high-dimensional integrals in time poses considerable challenge as the accuracy of calculations become extremely sensitive to numerical errors.
In this work, we have employed the Hierarchical Equations of Motion (HEOM) approach to effectively address these challenges.This methodology enables the systematic unraveling of higher-order master equations, balancing numerical efficiency with high precision.Consequently, our investigation has concentrated on the zerotemperature sub-Ohmic regimes, providing relevant computations to underscore our approach.
Furthermore, it is significant to note that the exact memory kernels in the quantum master equation can be extracted from the FP-HEOM approach and then be compared with those from perturbative schemes such as the NIBA.The extraction further offers a useful tool for the analysis of quantum phase transition dynamics, establishing a comprehensive and reliable platform for further explorations into open quantum system dynamics.