Taming Quantum Noise for Efficient Low Temperature Simulations of Open Quantum Systems

TL;DR

This work extends the hierarchical equations of motion (HEOM) to arbitrary temperatures and general reservoirs by introducing a rigorously optimized pole decomposition of the bath correlation function $C(t)$ via a rational barycentric representation implemented with the AAA algorithm. The resulting Free-Pole HEOM (FP-HEOM) uses a compact set of poles $z_k=i\omega_k+\gamma_k$ to represent $C(t)$ as $C(t)=\sum_{k=1}^K d_k e^{-z_k t}$, enabling accurate, long-time open-system dynamics with linear scaling in the effective mode number when combined with matrix product state techniques. The approach yields substantial performance gains over established methods, demonstrated on sub-ohmic spin-boson models at $T=0$ (including verification of the Shiba relation) and on structured reservoirs such as bandgap environments, with detailed convergence and cost analyses. By enabling stable, high-precision simulations across temperatures and reservoir types, FP-HEOM provides a practical benchmark tool and broad applicability to open quantum systems. Data supporting the results are available on request.

Abstract

The hierarchical equations of motion (HEOM), derived from the exact Feynman-Vernon path integral, is one of the most powerful numerical methods to simulate the dynamics of open quantum systems that are embedded in thermal environments. However, its applicability is restricted to specific forms of spectral reservoir distributions and relatively elevated temperatures. Here we solve this problem and introduce an effective treatment of quantum noise in frequency space by systematically clustering higher order Matsubara poles equivalent to an optimized rational decomposition. This leads to an elegant extension of the HEOM to arbitrary temperatures and very general reservoirs in combination with efficiency, high accuracy and long-time stability. Moreover, the technique can directly be implemented in alternative approaches such as Green's function, stochastic, and pseudo-mode formulations. As one highly non-trivial application, for the sub-ohmic spin-boson model at vanishing temperature the Shiba relation is quantitatively verified which predicts the long-time decay of correlation functions.

Figures (8)

• Figure 1: Flowchart for the optimized pole distribution and decomposition of the correlation $C(t)$ in the FP-HEOM.
• Figure 2: Dependence of the number of free poles $K$ to represent the noise power on the optimization tolerance $\delta$ in sub-ohmic case. The red dot corresponds the situation used in Figs. 1 and 2 in the main text. The parameters are: $s=1/2$, $\alpha = 0.05$, $\omega_c = 20$, $\mathcal{A} = [-10^3, -10^{-4}]\cup [10^{-4},10^3]$, and $T = 0$.
• Figure 3: Coherence to localization transition dynamics affected by $\alpha\in[0.05, 1.00]$. Solid lines denote FP-HEOM simulations, while black dot lines denote TD-DMRG simulations ren2022time. The simulation parameters are: $s=1/2$, $\epsilon = 0$, $\Delta = 1$, $\omega_c = 20$, and $T=0$.
• Figure 4: Shiba relations in sub-Ohmic case corresponds to $s = 0.7, 0.9, 1.0$. Black lines denote Shiba predictions, while red lines denote FP-HEOM simulations. The optimization parameters are: $\mathcal{A} = [-10^3, -10^{-4}]\cup [10^{-4},10^3]$, and accuracy $\delta = 10^{-8}$. Other simulation parameters are: $\epsilon = 0$, $\Delta = 1$, $\omega_c = 20$, $\alpha = 0.05$, ${T = 0}$.
• Figure 5: Relaxation dynamics of a TLS in bandgap environment. Simulation parameters are: $\kappa_1 = \kappa_2 = 2$, $\xi_1 = \xi_2 = 1$, $\omega_1 = 2$, $\omega_2 = 4$, $\epsilon = 6$, and $\Delta = 1$ so that the bare TLS transition frequency is $\delta\omega=2 \sqrt{\Delta^2+\epsilon^2/4}> \omega_2>\omega_1$.