Correlation Functions From Tensor Network Influence Functionals: The Case of the Spin-Boson Model

The forms of the equilibrium correlation functions given in Eq. (1) and Eq. (2) can be viewed as combinations of real-time and imaginary-time propagations and only differ in the position of the observable $B$. Taking advantage of this, we have developed a method based on propagations along the Kadanoff-Baym-like (KB) contour. The contour is schematically shown in Fig. 1a, where the points $a$, $b$, and $b^{\prime}$ correspond to locations where measurements of observables are made. Whereas $G_{AB}(t)$ requires insertions at points $a$ and $b$, the symmetrized correlator $C_{AB}(t)$ requires points $a$ and $b^{\prime}$.

Alternatively, the symmetrized correlator $C_{AB}(t)$ can be obtained via another approach explored by Bose [20], based on the complex-time (CT) contour shown in Fig. 1b, where the propagations occur over complex times. This contour differs from the KB contour in that it can only be used to compute the symmetrized correlation function for a single time $t$, though it was claimed [20] that this single calculation has low computational cost. We will revisit and expand upon this point in more detail.

With a specified choice of contour, $G_{AB}(t)$ and $C_{AB}(t)$ can be Trotterized and written in terms of influence functionals $F$. That is, the expressions Eq. (1) and Eq. (2) are recast into the general form,

Here, $n$ labels discrete positions along the contour, and $K$ encapsulates both the bare system propagations and any measurement of observables. Explicit expressions of $K$ and $F$ for thermal and symmetrized correlation functions can be found in section B of the Supplementary Materials.

In the spin-boson model, the harmonicity of the bath and its linear coupling to a diagonalizable operator of the system (with eigenvalues $s^{\pm}_{k}$) leads to the simple form,

where the coefficients $\eta_{k,k^{\prime}}$, which we loosely refer to as temporal correlations, are dependent on the choice of the contour.

We formulate our approach to calculating equilibrium correlation functions on the KB contour for its capability to directly compute both $G_{AB}(t)$ and $C_{AB}(t)$. Specifically, this is achieved via the process tensor formulation [35, 10] of the Time-Evolving Matrix Product Operator (TEMPO) algorithm, which constructs the influence functionals $F[\{s_{n}^{\pm}\}]$ as an MPS while maintaining the causal structure of evolution along the contour. This is achieved by making use of the commuting nature of all terms in $F$, which allows one to write the IF as the contraction of small tensors [9, 10] as illustrated in Fig. 2a. These small tensors are defined generically as

and similarly for $b$ and $c$. The different labels for the $b$, $c$, and $d$ tensors indicate correlations between real-real times, imaginary-imaginary times, and real-imaginary times respectively. After the IF has been constructed, correlation functions can be calculated by inserting measurements along the time contour, c.f. Eq. (4) and illustrated in Fig. 2b. The correlations are extracted over the entire temporal range $\{n\delta t|n=0,\ldots,N\}$ that the IF describes; in contrast, the CT contour limits the calculation of symmetrized correlation functions to a single time point per run, since the contour changes for every different time at which the operator $A(t)$ is measured.

Finally, to get an idea of how to reduce computational effort in the process tensor approach, which hinges upon the bond dimension of the MPS, we now examine the behavior of temporal correlations in the IF. These temporal correlations, occurring between different points along the CT contour, are given by the $\eta_{k,k^{\prime}}$, which are integrals of the bath correlation function, $L(t)$ (refer to Section. B of the Appendix for explicit expressions). One generally finds [20] that the magnitude of $L(t)$ is greatest along the imaginary time axis and decays with increasing distance from this axis, e.g., the points $b$ and $b^{\prime}$ in Fig. 1a should be more correlated than $b^{\prime}$ and $a$, as $t_{b}-t_{b^{\prime}}$ is purely imaginary while $t_{a}-t_{b^{\prime}}$ acquires a large real component. In practice we find that there is large entanglement around the boundary separating the indices $\{s^{\pm}_{i}\}$ and $\{s^{\beta}_{j}\}$ and between the indices $\{s^{\beta}_{j}\}$ themselves, in the notation of Fig. 2a. This issue of entanglement—and of large bond dimension by proxy—can be ameliorated by noting that the measurements required to calculate $G_{AB}$ or $C_{AB}$ will never be inserted at points on the contour between $b$ and $b^{\prime}$ in Fig. 1a. We can take advantage of this redundancy by using the causal nature of the construction of the influence functionals [10] to sum over the unneeded indices in $\{s^{\beta}_{i}\}$ once they no longer have influence on the construction of the remainder of the IF. This process is depicted in Fig. 2b.

The aforementioned ways of calculating $G_{AB}$ and $C_{AB}$ all construct the full equilibrium distribution $\exp(-\beta H)$, which is reflected by the need for imaginary time propagation. However, assuming that we are interested in the coupling of a finite sized system to a bath of infinitely many modes, one may presume the existence of a unique steady state given by the full equilibrium distribution. When such assumptions hold, one can convert the problem of calculating the equilibrium correlation function into the problem of calculating the steady state correlation function from an initially nonequilibrium state wherein the infinite bath is prepared at temperature $1/\beta$.

The process of propagation to infinite times is facilitated by the observation that the influence functionals are nearly translationally invariant, except near its temporal boundaries. Therefore, in the infinite time/steady state limit, the influence functionals can be represented as a uniform matrix product state (uMPS), with a repeating unit cell tensor $F_{\infty}$. While there are multiple ways of constructing this unit cell tensor, one approach noted by Link et al. [36] is to consider the repeating strip of tensors in the TEMPO network (see highlighted rectangle in Fig. 3a). Contracting these tensors using the infinite-time evolving block decimation (iTEBD) algorithm [37, 38, 39] directly yields $F_{\infty}$ (see Fig. 3a). Given the nature of this construction, $F_{\infty}$ will only include temporal correlations $\eta_{|k-k^{\prime}|,0}$ for $|k-k^{\prime}|=0,\ldots,N_{c}$, with $N_{c}\delta t$ defining the memory truncation time. In order to accurately capture the true long-time steady state, $N_{c}$ must be large enough such that the temporal correlations have almost vanished. Thus, cases in which temporal correlations decay slowly—such as with strongly sub-Ohmic baths—may require large $N_{c}$ and must be treated with care [36]. The overall cost of this approach thus depends on $N_{c}$ iterations of iTEBD, though we stress that this cost is not intrinsic to the steady state approach as there are different methods to obtain $F_{\infty}$.

Equipped with the uMPS $F_{\infty}$, we can now consider dynamics in the steady state. The $F_{\infty}$ is combined with free system propagations to construct an effective evolution tensor $T_{\delta t}$ (shown in the dashed orange box in Fig. 3b), describing the time evolution over a period of $\delta t$. The steady state is given by the right eigenvector $R$ of $T_{\delta t}$ with the eigenvalue of largest magnitude when there is a unique steady state of the dynamics. To $R$ there is an associated left eigenvector $L$ such that $|L\cdot R|>0$. From these definitions, the steady state correlation function is obtained via the tensor network ²²2Alternatively, one can directly obtain an analytical approximation for $G(t)$ as a sum of exponentials using the eigendecomposition of $T_{\delta t}$, which may be useful for analytic continuation. depicted in Fig. 3, corresponding to the expression,

We note, however, that this steady state correlation function may not capture the correct correlation function in cases where the long-time decay of the correlation function $|G_{AB}(t)-G_{AB}(\infty)|$ is algebraic, such as found at zero temperature in the spin-boson model [41]. The fact that MPS representations of the influence functionals with fixed bond dimension cannot reproduce algebraically decaying correlations should be expected, as this is a direct analogue of the situation found in real-space correlation functions in the ground state of gapless 1D systems [42]. Formally, to obtain a power law decay one would require that the bond dimension of $T$ grow as some polynomial of $k$, leading to an increased cost of the method. In practice, one can extract the asymptotic power-law behavior from the steady state approach, through a careful analysis of the spectrum of $T_{\delta t}$ [43, 44] as a function of increasing the bond dimension of the MPS.

Finally, we note that, unlike the previously stated approaches whereby the full thermal state is approximated in a systematically controllable way by the number of timesteps, the method we have outlined here for finding the steady state from real-time evolution is less straightforwardly controlled. In particular, the accuracy of this method will depend on the time discretization $\delta t$, the bond dimension $D$ of $F_{\infty}$, and the memory length $N_{c}$. While the first two parameters are all that control the accuracy of the methods involving complex-time propagations, the latter parameter of the memory length may play a significant role in the efficacy of the steady state approach.

For this low temperature ($\beta=5/\Delta$) and fast bath ($\omega_{c}=5\Delta$) regime, we expect that the equilibrium autocorrelation function of $\sigma^{z}{}$ should have qualitatively similar behavior to that found in of the well understood limit of zero temperature and $\Delta/\omega_{c}\to 0$. Particularly, there should be a crossover from underdamped to overdamped to incoherent decay as the coupling strength to the bath $\alpha$ increases. The results for three values of the couplings ($\alpha=0.1,0.3,1.0$) representative of these three regimes are shown in Fig. 4.

In all regimes, the various methods are broadly able to reproduce the reference calculations from ML-MCTDH (Fig. 4). All calculations with the process tensor approach were performed with the same number real- and imaginary-time discretizations, $\delta t=0.05/\Delta$ and $\delta\tau=0.125/\Delta$ respectively. Where available, the steady state correlation function is computed with $\delta t=0.04/\Delta$, with a memory cutoff $N_{c}=10^{5}$. The symmetrized correlation function computed on the CT contour are all calculated with a variable discretization $N$ such that the size of the timestep along the contour $|t-i\beta/2|/N\approx 0.05$. The lower portions of each panel show how the results in the upper portions differ from ML-MCTDH, with statistical uncertainties depicted by the shaded green regions. Here, the statistical uncertainties are taken as $2\times$ the standard error of the mean of the $24576$ trajectories assuming independent samples. These errors are likely underestimated at longer times such that the apparent deviations from ML-MCTDH in Fig. 4(b,c) for $t\Delta\gtrsim 7$ are exaggerated. Additionally, we find that the errors in the process tensor approach in Fig. 4a are governed by the timestep discretization $\delta t$ rather than by MPS compression. ⁴⁴4We show the additional convergence behaviors in Sec. A of the Supplementary Materials: the errors decrease as $O(dt^{2})$ for both the correlation functions on the KB contour and for the steady state correlation function.

While we have established the accuracy of the various methods presented in Sec. II for calculating the correlation functions $G(t)$ or $C(t)$, we note that the computational costs to achieve a desired level of accuracy differs greatly between the approaches. We quantify these differences in terms of computational complexities, memory requirements, and additional costs of the calculations, such as the computation of the $\eta_{k,k^{\prime}}$ coefficients entering into the influence functionals.

We first assess the effort required to converge each method to a fixed absolute error of 0.05 for $t\Delta\in[0,10]$. We show in Fig. 5 the effort to reach this level of convergence through the maximum number $N_{\text{tot}}$ of complex-valued elements in the IF-MPS, which quantifies the memory requirement across coupling strengths $\alpha$, as well as the corresponding bond dimension $D_{\text{max}}$. For the CT contour, since each calculation yields $C(t)$ at a specific time point, we perform convergence tests for each time $t$ considered. This convergence is evaluated with respect to both the discreteness of the Trotterization and the level of compression of the MPS. The data for the CT contour reports the largest number of elements observed within the specified time range.

In terms of computational storage, Fig. 5 shows that the process tensor approach is significantly more demanding than either the steady state and CT contour approaches. The process tensor also carries large bond dimensions which, when considered with the fact that the associated MPS consists of $N$ tensors as opposed to the lone unit cell tensor of the steady state’s uMPS, results in overall longer runtimes for tensor operations relative to the other two approaches. It can be seen that the costs associated with the process tensor grow monotonically with increasing coupling strength $\alpha$. This trend holds similarly for the steady state approach.

By contrast, the overall storage cost and bond dimensions for the CT contour approach are roughly constant across coupling strengths, though it is important to note that this is the cost associated with the calculation of $C(t)$ for a single value of $t$. At first glance, Fig. 5 may suggest that the steady state approach is a more expensive method compared to the CT contour approach. However, as noted by Link et al. [36], the bond dimension of the steady state uMPS remains relatively low during the majority of the iTEBD iterations, and only reaches the maximal bond dimension reported in Fig. 5 towards the end. In practice, we have found that all of the calculations in Fig. 5 using the steady state method took at most one minute.

In addition to the cost of constructing and storing the influence functionals in MPS form, it is important to note that there is a cost associated with the calculation of the $\eta_{k,k^{\prime}}$. These $\eta_{k,k^{\prime}}$ are computed for every pair of points along the chosen contour. Hence for the process tensor approach along the KB contour (Fig. 1a), there will be $O(MN)+O(N)+O(M)$ values of $\eta_{k,k^{\prime}}$ needed if there are $2M$ steps in imaginary time and $2N$ steps in real time. For the calculations shown in Fig. 4(a,b,c), $M\sim 10$ while $N\sim 10^{2}$, so that about $10^{3}$ values of $\eta_{k,k^{\prime}}$ are sufficient to compute the correlator over all $N$ time points. In the case of the steady state approach, for which the relevant temporal correlations enjoy a time-translationally invariant form $\eta_{k,k^{\prime}}=\eta_{|k-k^{\prime}|}$, only $N_{c}$ values of $\eta$’s are needed. This $N_{c}$ is chosen such that $\eta_{N_{c}}$ has largely decayed to zero; throughout this work $N_{c}\delta t=400/\Delta$, from which $G(t)$ can be calculated for all $t\Delta\lesssim 400$. By contrast, IFs on a CT contour (Fig. 1b) composed of $2N$ steps will require $O(N^{2})$ evaluations of $\eta_{k,k^{\prime}}$, one for each time point that $C(t)$ is to be evaluated over. We find that, in practice, this can pose significant costs for the CT contour approach when large values of $N$ are needed to converge the calculations. This is examined in closer detail in the following section.

As previously mentioned, there are distinct drawbacks to the calculation of symmetrized correlation functions that are independent of the chosen contour or MPS representation. The first is that calculations along the CT contour become increasingly difficult as the temperature is lowered or as the time increases. This is due to the need to maintain a certain level of discretization of the contour in order to control the Trotter error, which may be compounded by the need to calculate $C(t)$ to longer times as $C$ decays more slowly with increasing $\alpha$ or lower temperature. Second, in view of Eq. (3), features in the spectra obtained from $C(t)$ may become exponentially suppressed with inverse temperature $\beta$ relative to what may be found from $G(t)$. In such cases, errors may be exponentially amplified when attempting to recover $G(t)$ from $C(t)$. It would be crucial therefore to ensure that the calculation of $C(t)$ is properly converged, particularly with respect to Trotter errors. In this work, we leave aside questions of the numerical stability of converting $C(t)$ to $\widetilde{G}(\omega)$, and hence will not dwell on this second point.

We will now focus on the first point, examining the low temperature regime and comparing the use of the steady-state approach to calculate $G(t)$, versus the calculation of $C(t)$ on the CT contour. For numerical illustration, we examine more closely the long time behavior of $G(t)$ and $C(t)$ in Fig. 4(c,f). To resolve the decay behavior at intermediate times, we look at the correlation functions at times $t\Delta=10,20,40$ to examine the effort required to converge calculations with respect to the discretization in real ($\delta t$) or complex time ($\delta\tau=|t-i\beta/2|/N$ for number of timesteps $N$ along one leg of the contour in Fig. 1b). This is shown in Fig. 6. Even though the convergence behaviors differ between the steady state and CT contour approaches, they are converged at roughly similar sizes of Trotter steps along their respective contours. This highlights a potential pitfall of the CT contour approach, that the accurate resolution of $C(t)$ may become more difficult with increasing $t$. If the correlation function decays as $f(t)$ and the overall Trotter error scales as $O(t^{2}/N^{2})$, then in order to have fixed level of relative error at time $t$ would require $N\sim O(tf(t)^{-1/2})$. In the intermediate time regime $t\lesssim\beta$, it is expected [41] that $f\sim t^{-2}$. Thus the total number of $\eta_{k,k^{\prime}}$ computations per time $t$ would scale as $O(t^{4})$, which would make the CT contour approach much more computationally expensive than Fig. 5 might suggest.

Finally, we examine the utility of the steady state approach for calculating equilibrium correlation functions by focusing on the regime that cannot be exactly calculated by methods relying on imaginary time propagation. We will focus on the zero temperature limit of the spin-boson model, at the regime of couplings below the putative localization transition ($\alpha_{c}\simeq 1.25$ [9] for $\omega_{c}=5\Delta$), comparing the results against direct ground state simulations using ML-MCTDH.

In Fig. 7a, we show the zero temperature results for $\operatorname{Re}G(t)$ across coupling strengths from top to bottom, $\alpha=0.1,0.2,0.3,0.4,0.5,0.6$, which have been shifted for clarity. The results largely agree well with ML-MCTDH calculations, though closer inspection at larger couplings $\alpha\sim 0.6$ shows slight deviations from the reference results (Fig. 7b). It can be seen that the steady state correlation functions are converging towards the ML-MCTDH calculations as the bond dimension increases. Additionally, given that the zero temperature correlation functions in the unbiased spin-boson model are expected to decay algebraically [41], it is reasonable to anticipate the departure of the steady state calculations—which approximate $G(t)$ as the sum of exponentials—from the true result due to finite bond dimensions [49, 50]. This behavior can be explicitly seen in Fig. 7c, where the long-time behavior of $G(t)$ is exponential for small bond dimensions, e.g. $D=266$, and converges to the expected $\sim t^{-2}$ decay as $D$ increases. Lastly, a previous investigation by Rams et al. [43] suggests some utility of extrapolating to the infinite bond dimension limit of infinite MPS’s using a power-law ansatz, $x_{D}\sim x_{\infty}+cD^{p}$, even away from criticality. Applying this ansatz to our approximations to the correlation function, we recover reasonably good agreement with ML-MCTDH (Fig. 7b).