Benchmarking TCL4: Assessing the usability and reliability of fourth-order approximations Open Access

The study of open quantum systems is crucial for understanding the quantum mechanical interplay between a system and its surroundings. Exactly solvable models, such as the Jaynes–Cummings model¹ and the pure dephasing regime for most models, are rare because the interaction between the bath and the system introduces significant complexity, making analytical solutions challenging; thus, numerical methods become important. Numerical implementations of the Feynman–Vernon influence functional formalism, directly, such as the Time Evolution Matrix Product Operator (TEMPO),^2,3 or indirectly, such as Hierarchical Equations of Motion (HEOM),⁴ have gained popularity due to their accuracy and versatility, though each has its own compromises. In TEMPO, truncations of the memory length and data compressions with the singular value decomposition are necessary to manage the exponential growth in the dimensionality of the bath’s degrees of freedom, making computations challenging when the bath correlation time is long. HEOM, on the other hand, requires the bath correlation function to be expressible as a sum of exponentials to sustain the formalism of the auxiliary density operators (ADOs). To control memory usage, the hierarchy of ADOs must be truncated. HEOM becomes cumbersome when the bath correlation function cannot be well represented by exponentials, such as in environments with algebraically decaying bath correlation functions, which occur with spectral densities that feature an exponential cutoff,⁵ and the environment has a negative noise kernel, with the low temperature limit as an example.

While those methods can be practically applied to systems with limited size across all parameter regimes, the master equation methods using the Nakajima–Zwanzig projection technique can be applied to larger systems but usually within a certain parameter regime. Over the past several decades, various master equations have been proposed under different assumptions, each applicable to specific parameter regimes.^6–8 Interest in higher-order perturbative master equations has grown as researchers delve deeper into non-Markovian dynamics, where memory effects are significant as in realistic systems. The cumulant equation^9–11 systematically approximates the dynamics by expanding the dynamical map in terms of cumulants, describing the Liouvillian generator through a recursive relation of the second-order Bloch–Redfield generator, effectively capturing weak system–environment interactions and non-Markovian effects. Keldysh diagrammatic perturbation theory¹² recursively sums irreducible diagrams—with increasing interaction vertices and bath contractions—to transform the self-energy into effective Hamiltonian shifts and Lindblad dissipators. One of the most widely applied is the time-convolutionless (TCL) master equation,^13–15 which is theoretically exact and equivalent to perturbation to the Feynman–Vernon influence functional,^16,17 but practical implementation requires a truncation of perturbative expansion. The second-order TCL (TCL2), commonly known as the Bloch–Redfield equation, has been extensively explored and used in different areas of physics, such as quantum biology,¹⁸ quantum transport,^19,20 quantum computing, and recently the early universe in cosmology.²¹ However, due to the significant numerical computation required by the original triple integral of the TCL4 generator,^22,23 benchmarking work^10,23–25 for non-Markovian dynamics of TCL4 has been limited. Notably, we have recently developed a fast numerical implementation of TCL4,²⁶ reducing the triple integral to a single integration in terms of computational complexity, thus making parameter scans of TCL4 more accessible and feasible.

This work serves two purposes: comparing the two exact methods with TCL2 and TCL4 master equations using a spin-boson model and identifying the parameter space of temperature and bias where the accuracy of the TCL4 master equation is optimal, given the limited computational resources. To determine the domain of usability of the TCL4 perturbative master equation, the coupling strength to the bath will be chosen so that the master equation is near the edge of that domain.

The HEOM method is derived from the Feynman–Vernon influence functional formalism through integration by parts, which only works if the noise and dissipation kernels are exponential functions.⁴ This method reformulates the system–environment interaction by expanding it into a hierarchy of auxiliary density operators (ADOs), where the bath correlation functions are expressed as a sum of exponentials decay modes. Within this assumption, the ADOs are associated with the nearest ADOs by the differential equations. To ensure convergence and accuracy in HEOM simulations, several key parameters must be tuned. As shown in Fig. 1, the hierarchy depth N_k controls the number of ADOs included, with larger N_k capturing higher-order system–bath correlations but increasing computational costs. The number of exponential terms N_d determines how accurately the bath correlation function is approximated, where larger N_d improves precision but demands more resources. Therefore, a bath correlation function that cannot be accurately expressed in exponential modes poses a challenge for HEOM. Examples include low-temperature limits, where the non-Markovian negative noise kernel is difficult to express. Another example is a spectral density with an exponential cutoff, where the bath correlation function exhibits algebraic decay at large times. In conventional HEOM, one typically uses the Ishizaki–Tanimura low-temperature correction²⁸ to incorporate the contributions of high-frequency Matsubara modes. However, this scheme requires a very large number of auxiliary density operators to converge and may lead to numerical instabilities or even unphysical exponential decay when the hierarchy is truncated. Remedies for HEOM’s low-temperature instability and exponential cutoff issues include improved corrections and optimized representations of the bath correlation function. Fay²⁹ developed a projection-based correction that restores consistency with weak-coupling master equations using minimal auxiliary density operators. Xu et al.³⁰ introduced the free-pole HEOM (FP-HEOM) method, which clusters Matsubara poles via rational decomposition for efficient low-temperature simulations. Complementary fitting-based techniques, such as Prony filtration³¹ and Padé approximant decompositions,^4,32 further enhance convergence. These strategies collectively improve HEOM performance in challenging regimes and simulations. In this work, we use the conventional HEOM solver with the Matsubara expansion in the Qutip package^33,34 to perform HEOM computation.

Our TEMPO simulations utilize the open-source software detailed in the paper by Strathearn et al.² TEMPO employs the matrix product state/operator representation^35,36 to simulate the dynamics of open quantum systems, relying on singular value decomposition and the truncation of small singular values to reduce tensor dimensionality. The precision parameter λ_c sets the threshold for singular values retained in each decomposition, while the time evolution is dependent on Trotterization and discretization with a time step Δ_t. To capture the non-Markovian influence of the bath on the system, TEMPO propagates an augmented density tensor—a tensor documenting the possible system propagation choices from one time step to all subsequent steps—and retains K tensor legs at each time step to prevent exponential memory growth. Despite the TEMPO method being theoretically exact, the memory cutoff (K), precision parameter (λ_c), and time step (Δ_t) introduce inaccuracies into the computation. Although it is not possible to evaluate the effects of these parameters a priori, we empirically determine the confidence level of different parameters used in the simulation. Generally, we decide if the TEMPO results are trustworthy by observing the convergence when varying these parameters. However, in some cases, higher accuracy requirements (increasing K or λ_c) can lead to unstable data points in TEMPO. The goal of this article is to identify the validity region of TCL4. To achieve this, TEMPO parameters are considered acceptable if difference of TEMPO results of two sets of parameters exhibit variations much smaller than the discrepancies observed between TCL4 and TCL2. More specifically, if the parameter controlling convergence is increased and the trace distance of the density matrix (computed between the updated TEMPO parameters and the previous ones) is significantly smaller than the trace distance observed with TCL4 and TCL2, the new TEMPO parameters are finally adopted.

Figures 2 and 3 demonstrate the population and coherence dynamics at different temperatures and biases. As the bias and temperature (T) increase, the relaxation time of the system also increases. Specifically, at T = Ω, θ = π/20, and π/4, the population exhibits weak oscillations, whereas the population of θ = 9π/20 does not show an oscillation, indicating a transition from coherent transport to incoherent transport along θ. As θ increases, the system’s Hamiltonian H_S becomes more aligned with the system coupling operator σ_z, resulting in an increase in the dephasing rate and a reduction in the overall coherence of the system. At intermediate and high temperatures, population transport from the excited state becomes predominantly incoherent irrespective of bias as higher temperatures suppress the coherence by increasing the dephasing rate. Meanwhile, the temperature drives the steady-state upward across all three biases due to the higher bath correlation function. However, the relaxation rate decreases. In weak coupling theory, an increasing relaxation rate is typically observed as the temperature increases. This phenomenon is discussed in  Appendix A.

At low temperatures, HEOM faces significant challenges and shows poor agreement with TEMPO, despite aligning well with TEMPO at higher temperatures. HEOM relies on a finite set of basis functions, typically decaying exponentials, but at low temperatures, it requires substantially more exponential bath terms, making it computationally impractical. This issue stems from the scale-free noise spectrum at low frequencies, which is associated with zero-point quantum fluctuations. Although advances have aimed to address the limitations of HEOM in the low temperature regime by either including the correction terms or using curve fitting to shrink the bath terms needed,^4,29,37,38 here we use the traditional approach.

In contrast, the TCL2n master equation maintains a constant computational cost regardless of temperature. TCL4 is a result-reliable and computationally effective option at low temperatures, showing excellent consistency with TEMPO, especially for both population and coherence dynamics, as illustrated in Figs. 2(a)––2(c) and 3(a)–3(c). For example, the population elements between TCL4 and TEMPO differ by ∼2% at T = Ω, θ = π/20, whereas for TCL2, this difference is 15%. Regarding the coherences at T = Ω, TCL4 always aligns closely with TEMPO and TCL2 also aligns with TEMPO at low bias but starts to deviate from TEMPO at θ = π/4 and 9π/20. TCL4 effectively remedies the shortcomings of TCL2, increasing accuracy and correcting the equilibrium state at low temperatures.

Although TEMPO remains exact at low temperatures, its computational complexity grows exponentially with the simulation time, unlike the linear scaling seen with the HEOM and TCL methods. Furthermore, long-time simulations with TEMPO at low temperatures can introduce numerical instabilities, complicating its use for extended simulations.

As the temperature increases to T = 10Ω, the HEOM shows better agreement with TEMPO, as expected for this regime. For all the three biases, TEMPO and HEOM have around 0.01% difference of populations. Compared to the case of the low temperature limit, at T = 10Ω, the behavior of TCL4 against TCL2 diverges at different biases. In the case of TCL2, the population dynamics begins to match more closely the exact method at 10Ωt compared to the low temperature limit. TCL4 remains more accurate than TCL2 at low and medium bias levels but begins to lose reliability at higher biases, where the population shows notable deviations from the correct results. This trend becomes even more evident at higher temperatures. At T = 20Ω, TCL4 is only better than TCL2 in population dynamics at θ = π/20. For other populations and coherence dynamics at different biases, TCL4 is as erroneous as TCL2, though they deviate on opposite sides of the TEMPO curve due to the perturbative nature of the TCL equations. This suggests that TCL4 becomes unusable at excessively high temperatures.

Figure 4 presents the trace distance for coherence and population elements between TCL2, TCL4, and TEMPO at t = 10. As shown in Figs. 4(a) and 4(b), TCL2 exhibits significant errors in population dynamics at low temperatures, θ = 0 to θ = 9π/20, while TCL4 maintains much smaller errors. In contrast, TCL2 shows smaller errors in coherence dynamics, while the error in TCL4’s coherence dynamics remains below 1 × 10⁻⁴. In general, both TCL2 and TCL4 perform better in coherence dynamics than in population dynamics at low temperatures. As the temperature increases, the error in population dynamics for both TCL2 and TCL4 starts to increase at intermediate bias values although they continue to converge with TEMPO at zero and extreme bias. In contrast, the coherence dynamics follow a different trend, with the error increasing from zero bias to θ = 3π/8, as in Figs. 4(d) and 4(e). In Figs. 4(c) and 4(f), TCL4 does not consistently outperform TCL2. As the temperature exceeds 15Ω, TCL4 begins to perform worse than TCL2 in both coherence and population dynamics.

Figures 5(a) and 5(b) show the average time-averaged trace distance between TCLs and TEMPO. As in Sec. III B, TCL4 is significantly more accurate at low temperatures. The time-averaged trace distance between TCL2 and TEMPO is greater than 0.02 from π = 0 to 3π/8, while the error of TCL4 is below 0.005. The error of TCL2 approaches TCL4 when increasing θ, and at θ = π/2, the errors of TCL2 and TCL4 become equivalent, as shown in Fig. 5(c). This is a known result that in the pure dephasing regime, where H_s is parallel to σ_z, the second-order TCL equation can provide an exact solution. From low temperatures to T = 5Ω, the averaged trace distance of TCL2 generally decreases to values below 0.02. In this temperature region, TCL4 consistently outperforms TCL2. However, beyond T = 5Ω, the error in TCL2 increases as the temperature continues to rise, while TCL4 shows a similar trend. Notably, in Fig. 5(c), as the temperature is beyond T = 19Ω, a clear boundary emerges where TCL4 begins to perform worse than TCL2. In contrast to HEOM, high temperatures are the critical limitation for TCL4, whereas at low temperatures, TCL4 remains a reliable and useable method for systems coupled to an Ohmic bath.

The two parameter regions where TCL2 exhibits significantly higher error compared to other regions are the low-temperature limit and the high-temperature regime with increased bias. This discrepancy is reflected in the norm ratio, as shown in Fig. 6(b), which presents the norm ratio between the fourth-order and second-order generators as an indicator of the failure of perturbative expansion. The regions where the norm ratio peaks correspond closely to those where the time-averaged trace distance between TCL2 and TEMPO becomes prominent, as seen by comparing Fig. 5(a) with Fig. 6(b). The breakdown of TCL2 at low temperatures is attributed to positivity violations due to the negative noise kernel at early times; TCL4 substantially mitigates these positivity violations. Furthermore, positivity violation primarily impacts the population in the density matrix, as illustrated in Fig. 4(a). At high temperatures, invalidity arises from coherence elements, with the fourth-order generator improving accuracy up to T < 20Ω; however, TCL4 itself becomes unreliable beyond this point.

The only term that introduces temperature into the computation is the hyperbolic function $cothω2kBT$ in the integrand of the noise kernel. As shown in Fig. 6(a), when T increases, the noise kernel grows larger, leading to greater fluctuations. Consequently, the perturbation to the influence functional, which accounts for both fluctuation and dissipation, requires higher-order corrections due to the larger fluctuation terms. In the high temperature limit (T → ∞), the time-convolutionless (TCL2n) approach becomes impractical as it requires an infinite number of perturbation terms. It should also be noted that our coupling strength is relatively strong compared to the norm ratio between the fourth-order generator and the second-order generator. A lower coupling strength would result in a higher temperature range for the validity of TCL4.

We have seen an effect in Sec. III that the relaxation rate decreases with temperature, which is usually not the case in the weak coupling regime. A higher temperature means a higher spectral density and a stronger relaxation rate is expected, but here the parameters are near the edge of the weak coupling regime, which can change the temperature dependence of the relaxation rate.

This can be demonstrated with the relaxation and dephasing rates at second order, as detailed in  Appendix A. The three eigenvalues $λ=−J−Δ−JΔ,12−J−Δ−JΔ±J−Δ+JΔ2−4Δ2−ΔS−Δ+ΔSΔ$⁠. In the case $(J−Δ+JΔ)2<4Δ2−ΔS−Δ+ΔSΔ$⁠, the relaxation rate solely depends on the real part of spectral density (J_Δ + J_−Δ), which increases with rising temperature. However, if $(J−Δ+JΔ)2>4Δ2−ΔS−Δ+ΔSΔ$⁠, then the dynamics is overdamped,¹⁵ and the relaxation rate can decrease with increasing temperature. In the case of θ > 0, the term (S_−Δ − S_Δ) is automatically incorporated into the real part of dephasing, which does not increase monotonically with temperature. Figure 7 measures the relaxation rate for TEMPO and TCLs through curve fitting. We observe that the relaxation rate is at maximum at T = 2Ω and then quickly drops across the higher temperatures, hindering relaxation.

In many derivations of time-local master equations for open quantum systems, the assumption of invertibility for certain dynamical maps is central. For instance, Van Kampen’s approach, which expresses the reduced density matrix in terms of an exponential function of cumulants, implicitly requires the quantum dynamical map to be invertible.^39,40 Additionally, Chruściński and Kossakowski introduced a time-local generator that explicitly depends on the invertibility of the quantum dynamical map,⁴¹ with the generator defined as $L=Λ̇Λ−1$⁠. The recent cumulant equation,⁹ as well, requires an inverse of the dynamical map to reach a time-local equation. The TCL master equation, on the other hand, relies on the invertibility of an associated superoperator as expressed in Eq. (15). It remains a question whether the TCL equation maintains validity when the dynamical map becomes non-invertible.

We use quantum process tomography to reconstruct the dynamical map, starting from four initial states: $1000$⁠, $0100$⁠, $0010$⁠, and $0001$⁠. Dynamics are obtained using TCL4 at T = 10Ω. Subsequently, the dynamical map is constructed from these density matrices using the Kraus representation.⁴² At T = 10Ω, the dynamical map becomes non-invertible at t = 5.82/Ω, 5.73/Ω, and 5.62/Ω for all three biases investigated (θ = π/20, π/4, 9π/20). We have evaluated the dynamical maps also using the HEOM and found that the loss in invertibility of the maps occurs similarly at t = 5.85/Ω, 5.76/Ω, and 5.67/Ω for the same biases. TCL4 dynamics presented in Sec. III A remain robust compared to TEMPO and HEOM simulation. This indicates that invertibility is not a necessary condition for the validity of the TCL2n master equation. On the contrary, even when the dynamical map is non-invertible, TCL4 continues to align very closely with nonperturbative exact methods. We do not observe any direct relation between the loss of invertibility of the dynamical map and the validity of TCL4.

Benchmarking of the fourth-order time-convolutionless (TCL4) master equation highlights its ability to improve accuracy over the second-order TCL (TCL2) master equation, particularly for population and coherence dynamics in the low-temperature regime. TCL4 effectively reduces positivity violations seen in TCL2, which arise from the negative noise kernel at early times. In this regime, TCL4 demonstrates strong agreement with exact numerical methods, offering a computationally efficient alternative for capturing non-Markovian dynamics in open quantum systems. At higher temperatures and biases, the reliability of TCL4 decreases as fluctuations in the noise kernel grow and higher-order corrections become necessary. While TCL4 maintains better accuracy than TCL2 under intermediate conditions, performance degrades under high-temperature and high-bias conditions, often converging with or underperforming TCL2. These limitations align with the perturbative nature of TCL4, where the norm ratio of higher-order generators indicates the breakdown of the perturbative expansion in these parameter spaces.

Comparisons with exact numerical methods further emphasize the advantages and limitations of TCL4. TEMPO provides exact results but suffers from exponential scaling in computational complexity, making it impractical for long simulation times or large systems. HEOM offers robust performance at moderate to high temperatures but struggles at low temperatures due to the extensive basis required for representing the bath correlation function. In contrast, TCL4 delivers consistent computational efficiency across all regimes, offering a viable alternative in parameter spaces where exact methods are computationally prohibitive.

The authors express their gratitude to the School of Physics at the Georgia Institute of Technology for their invaluable support, which included access to computational resources and seed funding that significantly contributed to the success of this project.

The authors have no conflicts to disclose.

Jiahao Chen: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Elyana Crowder: Software (equal); Writing – review & editing (equal). Lian Xiang: Conceptualization (supporting); Data curation (supporting); Investigation (supporting); Software (supporting); Writing – review & editing (equal). Dragomir Davidovic: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The time-dependent spectral density is calculated from the correlation function through a direct integration from 0 to t_end, as shown in Eq. (7). This requires that t_N should be larger than t_end.

Unlike the second-order TCL generator, which only requires the spectral density at the local time, TCL4 involves a convolution of the spectral density over all past times. For the Markovian master equation, one only needs to calculate the generator at a large time where the generator saturates. The computation of a single generator requires a time complexity of O(M³n), where M is the dimension of the Bohr frequency and n is the number of timesteps up to the end time. If dt is fixed, the computational time scales linearly with the simulation time.

For the non-Markovian case, the TCL4 generators must be calculated at each timestep, resulting in a time complexity of O(M³n(n − 1)). In comparison, the TEMPO method requires a time complexity of O(nχ³M⁶K),⁴³ where χ is the bond dimension.

Here, we study the unbiased spin-boson model in the sub-Ohmic bath at zero temperature and compare our results with the example provided by FP-HEOM.³⁰ In Fig. 8, TCL4 yields reliable results at s = 0.5 and 0.7 but significantly differs from FP-HEOM in the deep sub-Ohmic regime at s = 0.1 and 0.3. Indeed, TCL4 fails to provide a localized steady state in these cases. However, a higher order TCL equation maybe required to reach the localization regime.⁴⁴ 

In complement of the Drude cutoff, here we present the comparison of the TCL2, TCL4, and TEMPO methods with exponential cutoff bath in Figs. 9 and 10. T = 0.1Ω, 5Ω, and 15Ω are used. A similar pattern of error is found in most regions. TCL4 manifests its advantage at low temperatures and its inability to operate at high temperatures of both coherence and population. However, the difference is at high temperatures and high dephasing area, where T = 15Ω and θ = 9π/20, where TCL4 shows a signature of bouncing coherence.

Figures 11(a) and 11(b) show Δ between TCL2 and TEMPO and between TCL4 and TEMPO. The diagram is more complicated than the one with Drude spectral density as one more region shows that TCL4 is worse than TCL2. The low-temperature dominance of TCL4 and the high-temperature, low-bias region are trending similarly to a bath with a Drude cutoff. In the case of pure dephasing (θ = π/2), both TCL4 and TCL2 show no deviation from TEMPO, confirming the exactness of the second-order TCL equation for pure dephasing. However, as soon as the bias shifts downward from θ = π/2, the equations of TCL2 and TCL4 cease to be exact and enter an error regime where TCL4 gains more deviation from the result of TEMPO. Figure 11(c) clearly manifests this area at temperatures above 4Ω, where TCL4 introduces more error in addition to the second-order generator due to coherence reoccurrence, meaning the coherence becomes bouncing along the exact result.