With the goal to study dissipative Landau–Zener (LZ) sweeps in realistic solid-state qubits, we utilize novel methods from non-Markovian open quantum system dynamics that enable reliable long-time simulations for sub-Ohmic environments. In particular, we combine a novel representation of the dynamical propagator, the uniform time evolving matrix product operator method, with a stochastic realization of finite temperature fluctuations. The latter greatly reduces the computational cost for the matrix product operator approach, enabling convergence in the experimentally relevant deeply sub-Ohmic regime. Our method allows the exact simulation of dynamical protocols with long operation times, such as the LZ sweep, in challenging parameter regimes that are realized in current experimental platforms.

Quantum dynamics are profoundly affected by environmental influences, leading to damping and coherence loss, which pose significant challenges for the design of scalable quantum hardware.1 Despite advancements, the limitations of current noisy intermediate-scale quantum devices have led to a shift in research focus from universal quantum computing toward more attainable goals, such as quantum simulation and annealing.2–4 These applications still require a precise control over individual quantum systems, necessitating a comprehensive theoretical understanding of the underlying noise processes. For instance, in quantum annealers, the dissipation poses a limit on the adiabaticity constraint, leading to compromises in the final state fidelity.2,5,6 However, in many currently available solid-state quantum devices, a thorough characterization of dissipation remains difficult, and the microscopic origins of the noise largely remain unclear.7 

An elementary protocol to test the influence of dissipation on a single qubit over different time scales is the famous Landau–Zener (LZ) problem.8,9 In this basic scheme, one considers the state transition probability of a qubit during a sweep across an avoided level crossing. In an ideal system, this can be described by the Hamiltonian
(1)
where σ̂i are Pauli matrices, for instance, σ̂z=|11||00|. The Landau–Zener transition probability P describes the probability for the qubit to be measured in the state |0⟩ at t = when it was prepared in the state |1⟩ at t = −. For pure unitary dynamics, an exact solution for this was found independently by Landau,10 Zener,11 Stueckelberg,12 and Majorana13 and is given by
(2)
In the adiabatic case (v ≪ Δ2), the probability becomes one, as the dynamics follows an avoided level crossing. However, in a real physical system, interaction with the environment has a strong impact on the LZ probability, especially for slow sweep velocities v.14–17 This makes the Landau–Zener protocol an interesting benchmark for evaluating the precision of theoretical models regarding the qubit dissipation. Due to a build-up of strong correlations between the system and the environment over time, conducting numerical simulations of dissipative dynamics can become highly nontrivial.

Previous exact numerical simulations of the dissipative Landau–Zener problem relied on the quasi-adiabatic path integral (QUAPI) approach and were restricted to environments with relatively short memory times. However, in recent experimental realizations of the Landau–Zener transition with flux qubits, strong low frequency noise contributions7,18,19 lead to highly non-Markovian dynamics for which the QUAPI method becomes unfeasible and master equation approaches become inaccurate. Here, we demonstrate that exact numerical simulations of such hard open system problems can be achieved by combining novel simulation techniques based on matrix product operator (MPO) methods20–31 with a stochastic realization of thermal fluctuations.32–34 This powerful numerical tool can be used in future studies to realize precise agreements of qubit dynamics in real world quantum devices with simulations and to deliver important insights into the characteristics of dissipation.

This paper is structured as follows: In Sec. II, we provide a brief review of open system dynamics and introduce the numerical tools that we use for our simulations. In Sec. III, we present our main results, exact simulations of the Landau–Zener transition in dissipative environments with very long memory times. In Sec. IV, we provide detailed introductions to the numerical methods and their practical implementation. Finally, we close with our conclusions and a review on future perspectives.

The standard approach to non-Markovian open system dynamics assumes a full total Hamiltonian of an extended system-and-bath model. We assume the following form:
(3)
where Ĥsys is the free system Hamiltonian and Ŝ is the system–bath coupling operator. For the dissipative Landau–Zener problem, the system Hamiltonian is given by (1) and the most dominant contribution to the dissipation is the so-called longitudinal noise, where Ŝ=σ̂z. The bath operator B̂(t) describes a continuum of bosonic modes,
(4)
with bosonic commutation [b̂(ω),b̂(ω)]=δ(ωω). This operator is related to the zero-temperature bath spectral density J(ω) = |g(ω)|2. For an initially uncorrelated Gaussian environment state, the bath correlation function α(ts) is sufficient to fully characterize the bath influence,
(5)
There is a manifold of methods available to solve the reduced dynamics according to Hamiltonian (3) for a given bath correlation function,35 i.e., to compute the dynamics of ρ(t) = trenvρtotal(t). The most recent advances utilize matrix product operator representations to efficiently capture temporal correlations.20,21,28,36–38 In this article, we use and promote the new simulation method introduced in Ref. 28, which is based on the time-evolving matrix product operator (TEMPO) framework.20 The starting point for TEMPO is the exact path integral representation of the reduced system density operator in terms of the Feynman–Vernon influence functional.32 TEMPO-based algorithms generate an optimized matrix product operator (MPO) representation of the influence functional with a low bond dimension such that the path integral can be evaluated. TEMPO has several crucial advantages over other exact approaches. In particular, only the bath correlation function needs to be provided to the algorithm and the accuracy of the method can be transparently controlled by changing the MPO compression threshold. The new scheme from Ref. 28 allows highly efficient simulations of long-time dynamics, as is required for describing the dissipative Landau–Zener transition with low sweep velocities.39 Further details on this method, which we call uniform TEMPO or “uniTEMPO,” are also provided in Sec. IV.

The crucial numerical bottleneck for the uniTEMPO computation is the final bond dimension of the MPO influence functional, as the computation time scales with the third power of this bond dimension. A long memory of the bath together with a strong system–bath coupling will require large bond dimensions, thus limiting the applicability of the approach. As demonstrated below, this problem arises in particular for strongly sub-Ohmic environments at high temperatures. Environments with such strong low-frequency contributions occur generically in solid state quantum devices. Thus, exact simulation techniques that can faithfully describe the open system evolution for large times are highly desired.

We demonstrate in the following that such difficult environments can still be handled numerically by combining the uniTEMPO method for zero temperature with an exact realization of finite temperature fluctuations in terms of stochastic unitary evolution. TEMPO approaches are particularly well suited for this stochastic approach because the main computational step, obtaining a MPO form for the influence functional, has to be performed only once. As has been observed in Ref. 21, when the MPO is available, the open system evolution can be computed very efficiently for different time-dependent system Hamiltonians.

In order to specify a bath for our numerical examples, we assume that the environment can be modeled with thermal harmonic modes. Note that, while it is usually safe to assume Gaussian statistics of the bath, this more restrictive assumption may not be justified. This holds in particular for solid state systems where the bath degrees of freedom may be non-harmonic. General baths may exhibit a different temperature dependence7,40 but can be dealt with analogously in simulations, assuming a temperature-dependent spectral density.41 For a harmonic bath at inverse temperature β = 1/(kBT), the bath correlation function reads
(6)
A generic model for the bath spectral density with tunable low-frequency behavior is the Ohmic-type spectral density. Here, we specifically choose (sub-)Ohmic spectral densities with an exponential high frequency cutoff,
(7)
where ωc is the cutoff frequency, the exponent s > 0 determines the low-frequency behavior, and α is a unit-less factor determining the coupling strength. For this model, the bath correlation function can be computed explicitly.42,43 From Eq. (6), we see that the bath correlation function splits into zero-temperature and finite-temperature components,
(8)
The zero temperature component causes dissipation and quantum fluctuations. Using the sub-Ohmic spectral density, one obtains
(9)
This function decays asymptotically as 1/ts+1. In contrast, the finite-temperature contribution is real valued,
(10)
where ζ(s, q) is the Hurwitz zeta function. It can be shown that the thermal bath correlation function decays only as 1/ts, that is one order slower than α0. For high temperatures, αβ dominates the bath correlation function because it has a much larger absolute value at t = 0 and a slower decay. This slow decay is especially drastic at low s, as can be seen in the example displayed in Fig. 1 (upper panel), where the different contributions to the bath correlation function are displayed for a sub-Ohmic bath. It is clearly visible that, while the zero-temperature correlation function decays quickly (weak memory), the finite temperature component has an extremely heavy tail. The slow decay effectively yields a huge memory time, hindering a straightforward numerical simulation of long-time dynamics.
FIG. 1.

Upper: Thermal and zero temperature components of the bath correlation function for a sub-Ohmic bath with s = 0.4. The thermal contribution to the correlation function strongly outweighs the zero-temperature components. Lower: Bond dimension of a MPO representation of the influence functional (via uniTEMPO) for a spin boson model with different sub-Ohmic environments. Depending on s, the bond dimension increases drastically with temperature. To ensure a reasonable comparison, the same MPO truncation threshold was used in all simulations.

FIG. 1.

Upper: Thermal and zero temperature components of the bath correlation function for a sub-Ohmic bath with s = 0.4. The thermal contribution to the correlation function strongly outweighs the zero-temperature components. Lower: Bond dimension of a MPO representation of the influence functional (via uniTEMPO) for a spin boson model with different sub-Ohmic environments. Depending on s, the bond dimension increases drastically with temperature. To ensure a reasonable comparison, the same MPO truncation threshold was used in all simulations.

Close modal

To exemplify this issue, we performed an example calculation for the spin boson model with different sub-Ohmic environments. We computed the MPO form for the influence functional obtained with uniTEMPO for a fixed compression tolerance. As can be seen in Fig. 1 (lower panel), the bond dimension of the influence functional MPO increases drastically as the bath temperature increases. This effect is particularly dramatic in the highly sub-Ohmic case (low s) because the tail of the thermal bath correlation function becomes extremely heavy. Note again that the computational demand scales with the third power of the bond dimension.

A central point that we like to highlight here is that the thermal part of the bath correlation function does not describe true dissipation, but instead, it can be seen as emerging from “classical” fluctuations. In fact, it is well established that the dynamics for a finite temperature bath can alternatively be described using zero temperature and an additional unitary stochastic driving included in the system Hamiltonian. Put differently, the thermal part of the bath correlation function αβ can be seen as describing stochastic unitary evolution. For this, one considers a real time-continuous Gaussian stochastic process y(t) characterized by the moments
(11)
where E[] denotes the process ensemble average. Then, we can define a corresponding stochastically driven Hamiltonian via
(12)
The average of stochastic dynamics according to this Hamiltonian including dissipation due to the zero temperature bath (open system dynamics with α0) will reproduce the exact thermal bath dynamics. A proof of this statement can be found in Refs. 32–34, as well as in Sec. IV. This stochastic treatment of temperature is very advantageous for the sub-Ohmic baths because we can perform open system evolution at arbitrary temperatures by utilizing uniTEMPO for a zero temperature bath and averaging over an ensemble of dynamics due to different stochastic system Hamiltonians. This way, we overcome the numerical issue that is displayed in Fig. 1 and perform exact simulations at arbitrary temperatures and very low exponents s. We demonstrate this in Sec. III, where we present exact simulations of the Landau–Zener transition in the challenging parameter regimes that are encountered in solid state quantum devices.

The LZ probability P is typically considered as a function of the sweep velocity v. For the particular basis chosen in (2), it describes the probability for reaching the σ̂z eigenstate σ̂z|0=|0 at large positive times when starting in the eigenstate σ̂z|1=|1 at large negative times. In addition to this so-called diabatic basis, the instantaneous eigenstates of the system Hamiltonian H(t), |+(t)⟩, |−(t)⟩ (adiabatic basis), can be considered at any time during the sweep. Note that the lower energy adiabatic state |−(t)⟩ of the system approaches the diabatic state |1⟩ for large negative times, but the state |0⟩ for positive times. Due to the avoided level crossing, the eigenenergies corresponding to the adiabatic basis states are always separated by at least the gap Δ. Having introduced these two bases, the (ground-state) LZ-probability can be expressed either as the probability of the system to stay in the same adiabatic state it was prepared in, P=||Û(,)||2, or equivalently as the probability of a transition between the two diabatic basis states, P=|0|Û(,)|1|2. For a perfectly isolated qubit, this probability is close to one for slow sweeps, as the system is able to evolve adiabatically and always stays in a momentary eigenstate. As the sweep velocity is increased, the LZ probability decreases monotonically, reaching zero for very fast sweeps. Adding the influence of an environment to the system evolution, the probability may deviate strongly from the analytical solution. In the case of dissipative dynamics, the LZ probability is influenced in a highly nontrivial way and can display a nonmonotonic behavior.18,44

As a first example, we revisit the case of an Ohmic bath s = 1 for which previous results using the QUAPI method are available.44 For weak to intermediate coupling, the Ohmic case is computationally rather simple because the bath correlation function decays quickly. Our results are displayed in Fig. 2. One can see a nonmonotonic behavior of the LZ curves, which can be explained in terms of a competition between dissipation and thermal fluctuations. At fast sweep velocities, the effective interaction time with the environment is too short to cause deviations from the unitary behavior. At intermediate sweep velocities, thermal fluctuations are the main influence. These fluctuations allow excitation to the excited state and hence cause the LZ probability to decrease. If the sweep is slowed down even further, the dissipation causes a relaxation to the lower energy state, leading to a subsequent increase in the ground state probability. The dissipation due to the environment decreases with decreasing temperature. In fact, at zero temperature, it was shown analytically that the LZ probability is unaffected by the bath in the case of purely longitudinal coupling (Ŝ=σ̂z).8 

FIG. 2.

Landau–Zener probability P as a function of the sweep velocity v at different temperatures T. The bath is Ohmic with s = 1 with an exponential cutoff ωc = 10Δ and coupling strength α = 2 × 10−3. The results obtained using the uniTEMPO method (black lines) are shown in comparison with earlier QUAPI calculations (+ markers).44 

FIG. 2.

Landau–Zener probability P as a function of the sweep velocity v at different temperatures T. The bath is Ohmic with s = 1 with an exponential cutoff ωc = 10Δ and coupling strength α = 2 × 10−3. The results obtained using the uniTEMPO method (black lines) are shown in comparison with earlier QUAPI calculations (+ markers).44 

Close modal

In Fig. 2, we compare our simulations with QUAPI results from Ref. 44. We find slight deviations of our converged uniTEMPO calculations with the QUAPI data. We suspect that these deviations arise from the Trotterization used in the QUAPI simulations, where a fairly large time step δ = 0.1 was chosen to cope with the exponential scaling of the numerical effort in QUAPI. We reach close agreement with the QUAPI results using this large time step in uniTEMPO. Clearly, a smaller value for δ is needed for numerical convergence. For the Ohmic bath, it is still feasible to use the bare uniTEMPO method without a stochastic realization of temperature. In this case, the full thermal correlation function is decaying fast and the required uniTEMPO bond dimension remains small even at high temperatures (see Fig. 1).

We now turn to the challenging sub-Ohmic (s < 1) regime where the bath memory becomes large. In this case, we employ the stochastic approach for finite temperature fluctuations in order to still handle the extremely slow decaying thermal part of the bath correlation function. The additional zero-temperature “quantum” part of the bath can be captured efficiently by uniTEMPO. The drawback of this approach is that, at high temperatures, many trajectories are required to reduce the statistical error of the data points. In Figs. 35, these statistical errors are smaller than the symbol size and are therefore not displayed. Note that for long sweep times, computing the dynamics of single realizations takes long computation times, even though parallelization is trivial.

FIG. 3.

Landau–Zener probability P as a function of sweep velocity v at different temperatures T. The bath is sub-Ohmic with s = 0.5 with an exponential cutoff ωc = 10Δ. Subfigures (a) to (c) correspond to different coupling strength values: (a) α = 0.316 × 10−3, (b) α = 0.949 × 10−3, and (c) α = 0.316 × 10−2. The data points (crosses) are obtained from averaging over 1000 realizations of the stochastic Hamiltonian (12), and lines are error-weighted spine interpolations. The black dashed line shows the Landau–Zener probability without dissipation [Eq. (2)].

FIG. 3.

Landau–Zener probability P as a function of sweep velocity v at different temperatures T. The bath is sub-Ohmic with s = 0.5 with an exponential cutoff ωc = 10Δ. Subfigures (a) to (c) correspond to different coupling strength values: (a) α = 0.316 × 10−3, (b) α = 0.949 × 10−3, and (c) α = 0.316 × 10−2. The data points (crosses) are obtained from averaging over 1000 realizations of the stochastic Hamiltonian (12), and lines are error-weighted spine interpolations. The black dashed line shows the Landau–Zener probability without dissipation [Eq. (2)].

Close modal
FIG. 4.

Landau–Zener probability P as a function of sweep velocity v at different temperatures T. The bath is highly sub-Ohmic with s = 0.1 with an exponential cutoff ωc = 10Δ. Subfigures (a) to (c) correspond to different coupling strength values: (a) α = 0.126 × 10−3Δ, (b) α = 0.378 × 10−3Δ, and (c) α = 0.126 × 10−2Δ. The data points (crosses) are obtained from averaging over 1000 realizations of the stochastic Hamiltonian (12), and lines are error-weighted spline interpolations. The black dashed line shows the Landau–Zener probability without dissipation [Eq. (2)].

FIG. 4.

Landau–Zener probability P as a function of sweep velocity v at different temperatures T. The bath is highly sub-Ohmic with s = 0.1 with an exponential cutoff ωc = 10Δ. Subfigures (a) to (c) correspond to different coupling strength values: (a) α = 0.126 × 10−3Δ, (b) α = 0.378 × 10−3Δ, and (c) α = 0.126 × 10−2Δ. The data points (crosses) are obtained from averaging over 1000 realizations of the stochastic Hamiltonian (12), and lines are error-weighted spline interpolations. The black dashed line shows the Landau–Zener probability without dissipation [Eq. (2)].

Close modal
FIG. 5.

Landau–Zener probability P as a function of sweep velocity v at different mixing angles θ. The bath is highly sub-Ohmic with s = 0.1 at a temperature T = 30Δ with a cutoff ωc = 100Δ and a coupling strength α0 = 0.317 × 10−5Δ. The results were obtained by averaging over 500 realizations of the stochastic Hamiltonian (12). The black line shows the Landau–Zener probability without dissipation [Eq. (2)].

FIG. 5.

Landau–Zener probability P as a function of sweep velocity v at different mixing angles θ. The bath is highly sub-Ohmic with s = 0.1 at a temperature T = 30Δ with a cutoff ωc = 100Δ and a coupling strength α0 = 0.317 × 10−5Δ. The results were obtained by averaging over 500 realizations of the stochastic Hamiltonian (12). The black line shows the Landau–Zener probability without dissipation [Eq. (2)].

Close modal

As concrete examples, we calculated the dissipative Landau–Zener problem for s = 0.5 in Fig. 3 and s = 0.1 in Fig. 4 at different coupling strengths and temperatures. Qualitatively, we find a similar behavior as in the Ohmic case in Fig. 2. Moreover, we find that the Landau–Zener probability is not very sensitive to the exact value of s, i.e., the low frequency behavior, if the values for the coupling strength α used lead to a comparable influence of the dissipation [as chosen in Figs. 3 and 4, panels (a)–(c)]. For s = 0.1, thermal fluctuations are generally more dominant at low frequencies, leading to a lower probability for reaching the ground state at low sweep velocities compared with s = 0.5.

Finally, we consider as an extension the case of a mixed coupling to the bath, where the bath coupling operator is not given by σ̂z. This so-called transverse noise has previously been studied in Refs. 8, 9, and 44. To interpolate between longitudinal noise and transverse noise, one can set the coupling operator to
(13)
where θ is called the “mixing angle.” Pure longitudinal coupling is recovered for θ = 0, and pure transverse coupling is recovered for θ = π/2.

In the example calculation displayed in Fig. 5, we consider a strongly sub-Ohmic bath at finite temperatures coupled via different mixing angles. Similar to the Ohmic case considered in Ref. 44, we find that transverse coupling leads to a stronger effect on the LZ probability compared with longitudinal coupling. This, in turn, requires a larger bond dimension for uniTEMPO in order to achieve convergence.

Overall, we did not find significant qualitative differences in the behavior of the LZ probability for sub-Ohmic baths compared with Ohmic baths. It should, however, be noted that environments of real quantum devices, such as superconducting flux qubits, are known to have strong low-frequency contributions within a deeply sub-Ohmic regime (the so-called 1/f noise).7,19,45,46 With the combination of stochastic fluctuations and uniTEMPO, we are able to perform exact calculations in this challenging regime over long evolution times with a relatively low computational effort.

In Sec. IV, we provide detailed explanations of the theoretical framework and the numerical methods that were used in the simulations.

In order to describe the evolution of the system state, we utilize a path integral theory, where the bath is captured in an exact way via the Feynman–Vernon influence functional. In the case of a Gaussian bath, the influence functional is a Gaussian with the bath correlation function as covariance. The basis of the influence functional theory is a Trotter splitting of the full unitary evolution operator generated by Ĥ(t)=Ĥsys(t)+Ĥint(t) [see Eq. (3)]. We assume a small but finite time step δ and second order Trotter splitting,
(14)
In the following, we consider the general eigenbasis |n⟩ of the coupling operator Ŝ|n=Sn|n, n = 1…d. It is also convenient to work in a Liouville space picture, where operators are denoted as vectors. For this vectorization, we use a single greek index μ = (μlμr) as a label for two (left and right) eigenstates such that matrix elements of operators are denoted as
(15)
In this notation, super-operators are denoted as matrices. For instance, the unitary channel generated by Ĥsys is written as
(16)
To describe the system evolution within the second order Trotter splitting, we further define the three-index tensor,
(17)
where one may notice that a sum over the center index gives the unitary channel for a single time step,
(18)
Within the Trotterization, we can compute the evolution of the system state for N time steps δ using the influence functional FNμ1μN, an N-index tensor.20,28 The open system evolution is given by the discrete path integral,
(19)
For a more convenient graphical notation, we can introduce pictorial representations for the objects in this equation. The density operator is a single-leg tensor,
(20)
A unitary channel is a linear map for density operators, so it should have two tensor legs. However, we introduced an additional “center” index in Eq. (17) in order to contract the channel with the influence functional such that it becomes a three-leg object,
(21)
The influence functional is a large N-leg tensor,
(22)
With this notation, we can write Eq. (19) in a pictorial tensor network notation,
(23)
where connected tensor legs imply summation.

For our model (3), the influence functional is well known analytically. Since the operator B̂(t) obeys Gaussian statistics with respect to the Gaussian bath initial state, a cumulant expansion truncates after the second order, which yields a Gaussian form of F.32 

Still, a direct computation of the path sum is practically impossible because there are exponentially many terms to sum over, each with a different complex phase. Tensor network approaches make it possible to find an approximate representation of the influence functional as a matrix product operator such that the sums can be performed iteratively and, hence, with a linear effort in N. The so-called PT-TEMPO method uses an exact representation of the influence functional as a two-dimensional tensor network to obtain such a form. The new approach introduced in Ref. 28, which we call “uniform TEMPO” (uniTEMPO), achieves a uniform MPO form for the influence functional by using infinite tensor network methods. With uniTEMPO, the influence functional can be represented in the form
(24)
where, for a given index μ, fμ is a square matrix (bond dimension χ) and vl/r are vectors realizing finite-time boundary conditions. The tensors f are all identical and independent of N. Such a representation can be highly efficient in capturing the time-nonlocal interaction with the bath. It results from a tensor network contraction algorithm, which performs automated compression to achieve an optimized bond dimension χ. With the MPO form (24), the open system evolution is given by
(25)
which can be computed straightforwardly in an iterative way from left to right (linear effort in N).
The crucial factor for the speed of the calculation is the bond dimension χ and the effort for generating the MPO form of F. Note that this last step requires matrix factorizations with computational complexity O(χ3). Usually, one can assume that the bath correlation function decays to zero and, hence, the bath has a finite memory time. We denote the number of time steps δ after which the bath correlation function has decayed sufficiently as Nc, i.e.,
(26)
The uniTEMPO contraction algorithm from Ref. 28 requires a total of Nc matrix factorizations to obtain fμ and vl/r. This makes the algorithm numerically superior to previous TEMPO-based methods that rely on finite MPO contraction algorithms.20,21,47 For the simulations in Figs. 25, obtaining the MPO form of the influence functional takes only a few seconds. The computational advantage is partly due to the number of required matrix factorizations needed to perform the contraction [finite network algorithms required O(Nc2) or O(NclogNc) effort22]. However, it is important to note that focusing only on such a formal scaling is misleading. This is because the individual matrix factorizations are performed on tensors with largely varying bond dimensions. In fact, choosing a very large Nc typically has little effect on the speed of the uniTEMPO algorithm because the contraction is performed in a reverse memory time direction such that bond dimensions at the beginning of the contraction are very small (see the supplementary material of Ref. 28). In contrast, in standard TEMPO, the contraction of the tensor network was performed in the evolution time direction, leading to large bond dimensions during substantial parts of the network contraction. We also like to highlight the low memory requirements of uniTEMPO: Only a single tensor must be stored at any given time in the computation.

An important issue that occurs in the case of sub-Ohmic spectral densities is the algebraic decay of the bath correlation function. True algebraic correlations cannot be realized in a MPO with finite bond-dimensions. In order to ensure the stability of the uniTEMPO algorithm, we include a smooth exponential cutoff to the bath correlation function after a large cutoff time. Exponential memory-time cutoffs are standard in common numerical methods for open system evolutions, such as the hierarchical equations of motion (HEOM). QUAPI and the finite PT-TEMPO algorithms use a hard memory cutoff, which, in contrast to the exponential cutoff, can lead to numerical instabilities when used in uniTEMPO. In general, the dynamics of the open system is unaffected by the cutoff if the cutoff time is larger than the considered evolution time. For longer evolution times, convergence of the results with respect to the cutoff should always be checked.

Starting from the influence functional theory, we can easily establish a proof for the stochastic realization of the finite-temperature contribution in the bath correlation function via Eq. (12). Because this contribution [αβ; see (8)] is real-valued, the corresponding Gaussian influence functional can be written as32 
(27)
where η is the discretized bath correlation function and is given as
(28)
The thermal influence functional (27) only takes real and positive values. Thus, the quadratic form in the exponent can be factorized by using a Hubbard–Stratonovich transformation.32,33 For this, one defines a real mean-zero Gaussian process yn with
(29)
Then, basic Gaussian integration yields the identity
(30)
Note that the transformed influence functional without the average is now time-local. In fact, it represents a time-dependent unitary evolution including the process yn as a classical “drive.” Upon defining an analogous time-continuous stochastic process y(t), the influence functional can be associated with a stochastic unitary evolution generated by the stochastic Hamiltonian (12).

The derivation above is applicable to any real-valued bath correlation function whenever the coupling to the bath is Hermitian. Moreover, for a non-Hermitian coupling to the bath, an alternative derivation of the stochastic decoupling based on the thermal P-function of the environment can still be applied in the case of both Gaussian34 and spin environments.48,49

Numerically generating stationary Gaussian processes with heavy-tailed temporal correlations is a nontrivial task. In particular, standard Fourier or convolutional methods become computationally expensive for correlations that decay as 1/ts with s ≪ 1. To circumvent a large numerical effort, we developed a novel scheme to efficiently sample a special class of heavy-tailed stationary Gaussian processes using a finite number of Ornstein–Uhlenbeck processes.50 Consider the special correlation function
(31)
which decays algebraically as 1/ty+1 (for the thermal noise, we need y = s − 1). The parameters θmax and c0 can be chosen freely to match the desired asymptotic correlation. When we apply a numerical quadrature scheme to the integral (31), it reduces to a finite sum over exponentials (in t),
(32)
with weights wk and abscissas θk. Such exponentially decaying correlation functions can be realized by stationary Ornstein–Uhlenbeck processes for which the transition probability is known analytically. Thus, the numerical quadrature scheme for the integral representation of the correlation function cf(t) delivers an approximation of the corresponding stochastic process as a sum over a finite number of stationary Ornstein–Uhlenbeck processes. This strategy allows an efficient and highly accurate generation of realizations of stationary Gaussian processes with a very long memory.
We use such a special process to realize the algebraic tail of the full thermal correlation αβ(t) and to simulate the remaining fast-decaying correlation with the Fourier method implemented in Ref. 51. In particular, we write the thermal bath correlation function as
(33)
where cf(t) is the correlation function of a special long memory process according to Eq. (31), realizing the algebraic tail of αβ(t). The remaining kernel α̃(t) then has a short memory such that samples can be generated with the Fourier approach. Once we have generated independent realizations of the short-time and heavy-tail processes, their sum will be a realization of the process with the desired correlation αβ(t).

In this work, we presented exact numerical simulations of the dissipative Landau–Zener transition probability for strongly non-Markovian sub-Ohmic environments. This setup directly corresponds to solid-state quantum devices, where 1/f noise is the primary source of dissipation.19 Previous simulations of this dissipation over extended periods have been limited by strong memory effects, which preclude the short-time truncation necessary for many numerical open system methods. We overcome these limitations by combining novel MPO-based simulation techniques28 with a long-known but rarely utilized stochastic sampling of thermal fluctuations.32,34 In particular, the new uniTEMPO approach enables the simulation of dynamics with strong dissipation over arbitrarily long times and with very low memory requirements. In addition, it can be easily integrated with the stochastic finite-temperature method, thereby mitigating convergence issues associated with the heavy-tail statistics of thermal noise.

Our framework provides a robust tool for verifying the quantitative agreement between experimental results and theoretical predictions without relying on additional approximations commonly used in deriving master equations. This capability allows a more accurate assessment of the noise characterization in qubit devices over extended periods and may thereby offer deeper insights into the underlying noise processes affecting the qubit dynamics.

The simulation method is, of course, completely independent of the specific driving protocol and can be used to simulate arbitrary dynamics of single or multiple qubits. For instance, one can straightforwardly study transitions in more intricate sweep protocols.52 For many-qubit systems, the numerical effort will scale exponentially with the number of qubits. However, for local dissipation, this exponential scaling could be lifted in favorable cases by using a MPO representation for the qubit state as well, similar to the approaches used in Refs. 25, 53, and 54.

It is a pleasure to thank Xi Dai, Adrian Lupascu, and Johannes Kassel for insightful discussions in connection with this work. All authors are grateful to the Center for Information Services and High Performance Computing [Zentrum für Informationsdienste und Hochleistungsrechnen (ZIH)] at TU Dresden for providing its facilities for high throughput calculations.

The authors have no conflicts to disclose.

Felix Kahlert: Conceptualization (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Valentin Link: Conceptualization (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Richard Hartmann: Conceptualization (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Walter T. Strunz: Conceptualization (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

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

1.
Y.
Kim
,
A.
Eddins
,
S.
Anand
,
K. X.
Wei
,
E.
van den Berg
,
S.
Rosenblatt
,
H.
Nayfeh
,
Y.
Wu
,
M.
Zaletel
,
K.
Temme
, and
A.
Kandala
, “
Evidence for the utility of quantum computing before fault tolerance
,”
Nature
618
,
500
505
(
2023
).
2.
J.
Preskill
, “
Quantum computing in the NISQ era and beyond
,”
Quantum
2
,
79
(
2018
); arXiv:1801.00862v3.
3.
A. J.
Daley
,
I.
Bloch
,
C.
Kokail
,
S.
Flannigan
,
N.
Pearson
,
M.
Troyer
, and
P.
Zoller
, “
Practical quantum advantage in quantum simulation
,”
Nature
607
,
667
676
(
2022
).
4.
K.
Bharti
,
A.
Cervera-Lierta
,
T. H.
Kyaw
,
T.
Haug
,
S.
Alperin-Lea
,
A.
Anand
,
M.
Degroote
,
H.
Heimonen
,
J. S.
Kottmann
,
T.
Menke
,
W.-K.
Mok
,
S.
Sim
,
L.-C.
Kwek
, and
A.
Aspuru-Guzik
, “
Noisy intermediate-scale quantum algorithms
,”
Rev. Mod. Phys.
94
,
015004
(
2022
).
5.
R.
Menu
,
J.
Langbehn
,
C. P.
Koch
, and
G.
Morigi
, “
Reservoir-engineering shortcuts to adiabaticity
,”
Phys. Rev. Res.
4
,
033005
(
2022
).
6.
E. C.
King
,
L.
Giannelli
,
R.
Menu
,
J. N.
Kriel
, and
G.
Morigi
, “
Adiabatic quantum trajectories in engineered reservoirs
,”
Quantum
8
,
1428
(
2024
); arXiv:2311.11937v2.
7.
C. M.
Quintana
,
Y.
Chen
,
D.
Sank
,
A. G.
Petukhov
,
T. C.
White
,
D.
Kafri
,
B.
Chiaro
,
A.
Megrant
,
R.
Barends
,
B.
Campbell
,
Z.
Chen
,
A.
Dunsworth
,
A. G.
Fowler
,
R.
Graff
,
E.
Jeffrey
,
J.
Kelly
,
E.
Lucero
,
J. Y.
Mutus
,
M.
Neeley
,
C.
Neill
,
P. J. J.
O’Malley
,
P.
Roushan
,
A.
Shabani
,
V. N.
Smelyanskiy
,
A.
Vainsencher
,
J.
Wenner
,
H.
Neven
, and
J. M.
Martinis
, “
Observation of classical-quantum crossover of 1/f flux noise and its paramagnetic temperature dependence
,”
Phys. Rev. Lett.
118
,
057702
(
2017
).
8.
M.
Wubs
,
K.
Saito
,
S.
Kohler
,
P.
Hänggi
, and
Y.
Kayanuma
, “
Gauging a quantum heat bath with dissipative Landau–Zener transitions
,”
Phys. Rev. Lett.
97
,
200404
(
2006
).
9.
K.
Saito
,
M.
Wubs
,
S.
Kohler
,
Y.
Kayanuma
, and
P.
Hänggi
, “
Dissipative Landau–Zener transitions of a qubit: Bath-specific and universal behavior
,”
Phys. Rev. B
75
,
214308
(
2007
).
10.
L. D.
Landau
, “
A theory of energy transfer. II
,” in
Collected Papers of L.D. Landau
(Pergamon, Oxford, England,
1965
), pp.
63
66
, ISBN: 978-0-08-010586-4.
11.
C.
Zener
, “
Nonadiabatic crossing of energy levels
,”
Proc. R. Soc. A
137
,
696
702
(
1932
).
12.
E. C. G.
Stueckelberg
, “
Theorie der unelastischen Stösse zwischen Atomen
,” in
E.C.G. Stueckelberg, An Unconventional Figure of Twentieth Century Physics
(Birkhäuser, Basel, Switzerland,
2009
), pp.
117
171
, ISBN: 978-3-7643-8878-2.
13.
E.
Majorana
, “
Oriented atoms in a variable magnetic field
,”
Nuovo Cimento
9
,
43
50
(
1932
).
14.
P.
Nalbach
and
M.
Thorwart
, “
Landau–Zener transitions in a dissipative environment: Numerically exact results
,”
Phys. Rev. Lett.
103
,
220401
(
2009
).
15.
P.
Nalbach
and
M.
Thorwart
, “
Competition between relaxation and external driving in the dissipative Landau–Zener problem
,”
Chem. Phys.
375
,
234
242
(
2010
).
16.
L.
Arceci
,
S.
Barbarino
,
R.
Fazio
, and
G. E.
Santoro
, “
Dissipative Landau–Zener problem and thermally assisted quantum annealing
,”
Phys. Rev. B
96
,
054301
(
2017
).
17.
D.
Wang
and
D.
Xu
, “
Nonadiabatic evolution and thermodynamics of a time-dependent open quantum system
,”
Phys. Rev. A
104
,
032201
(
2021
).
18.
X.
Dai
,
R.
Trappen
,
H.
Chen
,
D.
Melanson
,
M. A.
Yurtalan
,
D. M.
Tennant
,
A. J.
Martinez
,
Y.
Tang
,
E.
Mozgunov
,
J.
Gibson
,
J. A.
Grover
,
S. M.
Disseler
,
J. I.
Basham
,
S.
Novikov
,
R.
Das
,
A. J.
Melville
,
B. M.
Niedzielski
,
C. F.
Hirjibehedin
,
K.
Serniak
,
S. J.
Weber
,
J. L.
Yoder
,
W. D.
Oliver
,
K. M.
Zick
,
D. A.
Lidar
, and
A.
Lupascu
, “
Dissipative Landau–Zener tunneling: Crossover from weak to strong environment coupling
,” arXiv.2207.02017 (
2022
).
19.
R.
Trappen
,
X.
Dai
,
M. A.
Yurtalan
,
D.
Melanson
,
D. M.
Tennant
,
A. J.
Martinez
,
Y.
Tang
,
J.
Gibson
,
J. A.
Grover
,
S. M.
Disseler
,
J. I.
Basham
,
R.
Das
,
D. K.
Kim
,
A. J.
Melville
,
B. M.
Niedzielski
,
C. F.
Hirjibehedin
,
K.
Serniak
,
S. J.
Weber
,
J. L.
Yoder
,
W. D.
Oliver
,
D. A.
Lidar
, and
A.
Lupascu
, “
Decoherence of a tunable capacitively shunted flux qubit
,” arXiv.2307.13961 (
2023
).
20.
A.
Strathearn
,
P.
Kirton
,
D.
Kilda
,
J.
Keeling
, and
B. W.
Lovett
, “
Efficient non-Markovian quantum dynamics using time-evolving matrix product operators
,”
Nat. Commun.
9
,
3322
(
2018
).
21.
G. E.
Fux
,
E. P.
Butler
,
P. R.
Eastham
,
B. W.
Lovett
, and
J.
Keeling
, “
Efficient exploration of Hamiltonian parameter space for optimal control of non-Markovian open quantum systems
,”
Phys. Rev. Lett.
126
,
200401
(
2021
).
22.
M.
Cygorek
,
J.
Keeling
,
B. W.
Lovett
, and
E. M.
Gauger
, “
Sublinear scaling in non-Markovian open quantum systems simulations
,”
Phys. Rev. X
14
,
011010
(
2024
).
23.
A.
Lerose
,
M.
Sonner
, and
D. A.
Abanin
, “
Influence matrix approach to many-body Floquet dynamics
,”
Phys. Rev. X
11
,
021040
(
2021
).
24.
M.
Cygorek
,
M.
Cosacchi
,
A.
Vagov
,
V. M.
Axt
,
B. W.
Lovett
,
J.
Keeling
, and
E. M.
Gauger
, “
Simulation of open quantum systems by automated compression of arbitrary environments
,”
Nat. Phys.
18
,
662
668
(
2022
).
25.
G. E.
Fux
,
D.
Kilda
,
B. W.
Lovett
, and
J.
Keeling
, “
Tensor network simulation of chains of non-Markovian open quantum systems
,”
Phys. Rev. Res.
5
,
033078
(
2023
).
26.
R.
Chen
,
X.
Xu
, and
C.
Guo
, “
Grassmann time-evolving matrix product operators for equilibrium quantum impurity problems
,”
New J. Phys.
26
,
013019
(
2024
).
27.
R.
Chen
,
X.
Xu
, and
C.
Guo
, “
Grassmann time-evolving matrix product operators for quantum impurity models
,”
Phys. Rev. B
109
,
045140
(
2024
).
28.
V.
Link
,
H.-H.
Tu
, and
W. T.
Strunz
, “
Open quantum system dynamics from infinite tensor network contraction
,”
Phys. Rev. Lett.
132
,
200403
(
2024
).
29.
M.
Cygorek
and
E. M.
Gauger
, “
ACE: A general-purpose non-Markovian open quantum systems simulation toolkit based on process tensors
,”
J. Chem. Phys.
161
(7), 074111 (
2024
).
30.
G. E.
Fux
,
P.
Fowler-Wright
,
J.
Beckles
,
E. P.
Butler
,
P. R.
Eastham
,
D.
Gribben
,
J.
Keeling
,
D.
Kilda
,
P.
Kirton
,
E. D. C.
Lawrence
,
B. W.
Lovett
,
E.
O’Neill
,
A.
Strathearn
, and
R.
de Wit
, “
OQuPy: A Python package to efficiently simulate non-Markovian open quantum systems with process tensors
,”
J. Chem. Phys.
161
(12),
124108
(
2024
).
31.
C.
Guo
and
R.
Chen
, “
Infinite Grassmann time-evolving matrix product operator method in the steady state
,”
Phys. Rev. B
110
,
045106
(
2024
).
32.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
1965
).
33.
W. T.
Strunz
, “
Linear quantum state diffusion for non-Markovian open quantum systems
,”
Phys. Lett. A
224
,
25
30
(
1996
).
34.
R.
Hartmann
and
W. T.
Strunz
, “
Exact open quantum system dynamics using the hierarchy of pure states (HOPS)
,”
J. Chem. Theory Comput.
13
,
5834
5845
(
2017
).
35.
I.
de Vega
and
D.
Alonso
, “
Dynamics of non-Markovian open quantum systems
,”
Rev. Mod. Phys.
89
,
015001
(
2017
).
36.
J.
Prior
,
A. W.
Chin
,
S. F.
Huelga
, and
M. B.
Plenio
, “
Efficient simulation of strong system–environment interactions
,”
Phys. Rev. Lett.
105
,
050404
(
2010
).
37.
I.
Vilkoviskiy
and
D. A.
Abanin
, “
Bound on approximating non-Markovian dynamics by tensor networks in the time domain
,”
Phys. Rev. B
109
,
205126
(
2024
).
38.
M.
Cygorek
and
E. M.
Gauger
, “
Understanding and utilizing the inner bonds of process tensors
,” arXiv.2404.01287 (
2024
).
39.
J.
Iles-Smith
,
O.
Diba
, and
A.
Nazir
, “
Capturing non-Markovian polaron dressing with the master equation formalism
,”
J. Chem. Phys.
161
(13),
134111
(
2024
).
40.
S.
Sendelbach
,
D.
Hover
,
A.
Kittel
,
M.
Mück
,
J. M.
Martinis
, and
R.
McDermott
, “
Magnetism in SQUIDs at millikelvin temperatures
,”
Phys. Rev. Lett.
100
,
227006
(
2008
).
41.
N.
Makri
, “
The linear response approximation and its lowest order corrections: An influence functional approach
,”
J. Phys. Chem. B
103
,
2823
2829
(
1999
).
42.
E.
Crowder
,
L.
Lampert
,
G.
Manchanda
,
B.
Shoffeitt
,
S.
Gadamsetty
,
Y.
Pei
,
S.
Chaudhary
, and
D.
Davidović
, “
Invalidation of the Bloch-Redfield equation in the sub-Ohmic regime via a practical time-convolutionless fourth-order master equation
,”
Phys. Rev. A
109
(5),
052205
(
2024
).
43.
R.
Hartmann
, “
Exact open quantum system dynamics—Investigating environmentally induced entanglement
,” Ph.D. thesis,
Dresden, Tech. U.
,
2021
.
44.
S.
Javanbakht
,
P.
Nalbach
, and
M.
Thorwart
, “
Dissipative Landau–Zener quantum dynamics with transversal and longitudinal noise
,”
Phys. Rev. A
91
,
052103
(
2015
).
45.
T.
Lanting
,
M. H.
Amin
,
C.
Baron
,
M.
Babcock
,
J.
Boschee
,
S.
Boixo
,
V. N.
Smelyanskiy
,
M.
Foygel
, and
A. G.
Petukhov
, “
Probing environmental spin polarization with superconducting flux qubits
,” arXiv.2003.14244 (
2020
).
46.
D. A.
Rower
,
L.
Ateshian
,
L. H.
Li
,
M.
Hays
,
D.
Bluvstein
,
L.
Ding
,
B.
Kannan
,
A.
Almanakly
,
J.
Braumüller
,
D. K.
Kim
,
A.
Melville
,
B. M.
Niedzielski
,
M. E.
Schwartz
,
J. L.
Yoder
,
T. P.
Orlando
,
J. I.-J.
Wang
,
S.
Gustavsson
,
J. A.
Grover
,
K.
Serniak
,
R.
Comin
, and
W. D.
Oliver
, “
Evolution of 1/f flux noise in superconducting qubits with weak magnetic fields
,”
Phys. Rev. Lett.
130
,
220602
(
2023
).
47.
M. R.
Jørgensen
and
F. A.
Pollock
, “
Exploiting the causal tensor network structure of quantum processes to efficiently simulate non-Markovian path integrals
,”
Phys. Rev. Lett.
123
,
240602
(
2019
).
48.
V.
Link
and
W. T.
Strunz
, “
Stochastic Feshbach projection for the dynamics of open quantum systems
,”
Phys. Rev. Lett.
119
,
180401
(
2017
).
49.
V.
Link
,
K.
Luoma
, and
W. T.
Strunz
, “
Non-Markovian quantum state diffusion for spin environments
,”
New J. Phys.
25
,
093006
(
2023
).
50.
V.
Link
, Longmemoryprocess, github.com/val-link,
2023
.
51.
OpQuSyD
, stocproc, github.com/OpQuSyD,
2024
.
52.
F.
Zhuang
,
J.
Zeng
,
S. E.
Economou
, and
E.
Barnes
, “
Noise-resistant Landau–Zener sweeps from geometrical curves
,”
Quantum
6
,
639
(
2022
); arXiv:2103.07586v3.
53.
S.
Flannigan
,
F.
Damanet
, and
A. J.
Daley
, “
Many-body quantum state diffusion for non-Markovian dynamics in strongly interacting systems
,”
Phys. Rev. Lett.
128
,
063601
(
2022
).
54.
X.
Gao
,
J.
Ren
,
A.
Eisfeld
, and
Z.
Shuai
, “
Non-Markovian stochastic Schrödinger equation: Matrix-product-state approach to the hierarchy of pure states
,”
Phys. Rev. A
105
,
L030202
(
2022
).