HOQST: Hamiltonian Open Quantum System Toolkit

We present an open-source software package called"Hamiltonian Open Quantum System Toolkit"(HOQST), a collection of tools for the investigation of open quantum system dynamics in Hamiltonian quantum computing, including both quantum annealing and the gate-model of quantum computing. It features the key master equations (MEs) used in the field, suitable for describing the reduced system dynamics of an arbitrary time-dependent Hamiltonian with either weak or strong coupling to infinite-dimensional quantum baths. This includes the Redfield ME, the polaron-transformed Redfield ME, the adiabatic ME, the coarse-grained ME, and the universal Lindblad ME. HOQST also includes the stochastic Schrodinger equation with spin-fluctuators. We present an overview of the theories behind the various MEs and provide examples to illustrate typical workflows in HOQST. We present an example that shows that HOQST can provide order of magnitude speedups compared to QuTiP, for problems with time-dependent Hamiltonians. The package is ready to be deployed on high-performance computing (HPC) clusters and is aimed at providing reliable open-system analysis tools for noisy intermediate-scale quantum (NISQ) devices. The HOQST Github repository (https://github.com/USCqserver/OpenQuantumTools.jl) provides the starting point for users. Detailed information can be found in the README file.


INTRODUCTION
The theory of open quantum system has been an important subfield of quantum physics during the past decades with a rich collection of well-established methods [1][2][3]. Since perfect isolation of quantum systems is impossible, any quantum mechanical system must be treated as an open system in practice. The theory of open quantum system thus plays a major role in various applications of quantum physics, e.g., quantum optics [4,5], quantum control [6], and quantum computing (QC) [7]. It becomes even more relevant in the context of Hamiltonian quantum computing (HQC), broadly defined as analog QC performed via continuously and smoothly driven Hamiltonians, as opposed to discrete gate-model QC, where Hamiltonians are driven discontinuously. Well-known example of HQC include adiabatic quantum computing (AQC) [8,9] and quantum annealing (QA) [10,11], as well as holonomic QC [12,13]. For example, in QA the Hamiltonian needs to move continuously from the initial driver Hamiltonian to the final problem Hamiltonian, and therefore, unlike idealized gate-model quantum computers whose description often involves effective noise channels (completely positive maps), practical quantum annealers are better described by noise models derived directly from first principles [14][15][16][17][18]. As the entire field of QC is now in the noisy intermediate-scale quantum (NISQ) era [19], an efficient and evolving framework of open quantum system simulation is essential for advancing our understanding of noise in quantum devices, as well for helping in the search for more effective error suppression and correction techniques [20]. Moreover, the distinction between analog and discrete models of QC is to some degree arbitrary, since in reality, even gate-model QC involves continuous driving due to the finite bandwidth of signal generators and controllers. We thus view the gate-model of QC as part of HQC for the purposes of this work.
At present, there is an increasing number of software tools being developed for open system simulations. An important example of open-source software in this area is "Quantum Toolbox in Python" (QuTiP) [21]. It is one of the first packages in the field to adopt the modern software engineering paradigm and is actively maintained on Github since its release, with new features and enhancements being continuously added. However, since QuTiP is designed to be as general as possible, it lacks several tools and the computational performance required to address the new challenges we are now facing in the field of HQC.
Inspired by this challenge, and by the success of QuTiP, we present here a complementary and alternative open system simulation framework, which we call "Hamiltonian Open Quantum System Toolkit" (HOQST). As the name suggests, the goal of HOQST is not to cover the entire field of open quantum system simulation but to focus on Hamiltonian QC, while retaining the flexibility to simulate systems subject to arbitrary time-dependent Hamiltonians. This focus gives us the ability to adopt domain-specific design choices and optimizations. The resulting implementation distinguishes itself from other available software by offering the following advantages: • HOQST is written in the Julia programming language [22], which is designed for high performance computing.
• HOQST is built upon the ordinary differential equations (ODE) package DifferentialEquations.jl [23]; thus HOQST also benefits from progress in the field of ODE solvers.
• Focusing solely on Hamiltonian QC, HOQST features several recently published master equations.
• HOQST includes tools that work beyond the weak coupling limit.
• HOQST provides a native interface for HPC clusters.
HOQST is developed following the Julia design philosophy: we intend it to be as user-friendly as possible without compromising performance. Although there is room for optimization, the first release of HOQST features reliable and efficient implementations of several key master equations (MEs) adopted in the HQC field, together with a highly modularized framework suitable for future development. Since the HOQST project started as an attempt to build a tool to simulate quantum annealing, it displays a certain bias towards QA. However, we reemphasize that it is broadly applicable to open quantum systems evolving subject to any time-dependent Hamiltonian. Besides the HOQST package itself, we provide error bounds on the MEs included. We also present examples to illustrate the typical workflows of HOQST. In particular, we focus on a three-qubit entanglement witness experiment performed using a quantum annealer [24]. Previous studies of open system models were unable to reproduce some of the key experimental features; we demonstrate that HOQST does now offer this capability.

RESULTS
Throughout this work we consider a quantum mechanical system S coupled to a bath B. The total Hamiltonian is assumed to have the following form where H S and H B denote, respectively, the free system and bath Hamiltonians. H I is the system-bath interaction, which is often written as where A α and B α are dimensionless Hermitian operators acting on the system and the bath, respectively, and exclude both I S and I B (the identity operators on the system and the bath, respectively). The parameters g α are sometimes absorbed into B α but are kept explicit in HOQST, are have units of energy. In addition, we assume for simplicity the factorized initial condition ρ(0) = ρ S (0) ⊗ ρ B for the joint system-bath state at the initial time t = 0, where ρ B is a Gibbs state at inverse temperature β though we note that factorization is not necessary for a valid description of open system dynamics [25,26]. We work in units ofh = 1 and k B = 1 so that β has units of inverse energy, or time. Before proceeding, we transform the original Hamiltonian Eq. (1) into a rotating frame defined by U (t), i.e.
The specific form of the unitary operator U (t) will lead to different master equations and will be discussed in detail in subsequent sections. The goals of the rotation are to remove the pure-bath Hamiltonian H B and to identify terms that remains small in different system bath coupling regimes. As long as U (t) acts non-trivially on the bath, we may assume without loss of generality we that the rotating frame HamiltonianH has the following form: whereH S acts only the system andH I acts jointly on the system and the bath. The Liouville von Neumann equation in this rotating frame is whereL(t) denotes the Liouvillian superoperator. Once again we can always writẽ whereÃ α andB α are, respectively, system and bath operators (excluding identity). However, it is important to note thatÃ α andB α do not necessarily correspond to U † (t)A α U (t) and U † (t)B α U (t) in Eq. (2) because U (t) may not preserve the tensor product structure (i.e., we allow for denote the expectation value of any rotating frame bath operatorX with respect toρ B . Then the two-point correlation function is If the correlation function is time-translation-invariant then the noise spectrum of the bath can be properly defined by taking the Fourier transform The widely used Ohmic bath case is

Timescales
We define the two timescales to measure the range of applicability [27] 1 Here t f is the total evolution time, used as a cutoff which can often be taken as ∞. The quantity τ SB is the fastest system decoherence timescale, or timescale over which the system density matrix ρ S changes due to the coupling to the bath, in the interaction picture. The quantity τ B is the characteristic timescale of the decay of C(τ ). Note that the expression for τ B becomes an identity if we choose |C(τ )| ∝ e −τ /τB and take the limit t f → ∞. τ SB and τ B are the only two parameters relevant for determining the range of applicability of the various master equations discussed here [27], with the Universal Lindblad Equation (ULE; see Methods) case being no exception [28]. For convenience we collect the corresponding error bounds here, before discussing the various MEs. Namely, the error bound of the Redfield master equation is where ρ true (t) denotes the true (approximation-free) state.
The error bound of the Davies-Lindblad master equation is where δE =min i =j |E i −E j | is the level spacing, with E i the eigenenergies of the system Hamiltonian H S . This original version of this ME [29] does not directly allow for time-dependent driving, and we shall consider an adiabatic variant that does, the adiabatic master equation (AME) [30] (see Methods). The same error bound should apply in this case, since the difference is only in that the Lindblad operators are rotated in the AME case with the (adiabatically) changing eigenstates of H S , and the norms used to arrive at Eq. (15) are invariant under this unitary transformation. The error bound of the coarse-grained master equation (GCME; see Methods) is The error of the polaron transform master equation (PTRE; see Methods) can be separated into two parts. The first part comes from truncating the expansion to 2nd order. It can be bounded using the same expression as in Eq. (14), with timescales defined by the polaron frame correlation K(τ ): A detailed discussion of the above quantities is presented in Supplementary Note 2 and 3 (the explicit form of K(τ ) is also described in the Polaron Transform subsection in the Methods), and as far as we know the error bounds we derive here for the PTRE are new. We mention here that if the system-bath coupling strength g α in Eq. (2) is sent to infinity, both 1/τ SB and τ B /τ SB go to 0. Thus the PTRE works in the strong coupling regime. The second part of the error is caused by ignoring the 1st and 2nd order inhomogeneous terms, which themselves are due to the polaron transformation breaking the factorized initial condition. We do not have a bound on this error yet, but numerical studies suggest it is small when ρ S (0) is diagonal [31]. These bounds assume a Gaussian bath. For a non-Gaussian bath extra timescales relating to higher-order correlation functions generally appear, and the error bounds will contain additional terms.

Capabilities of HOQST and comparison with other quantum simulators
Recent developments in the field of QC have led to an explosion of quantum software platforms [32], such as Qiskit [33,34], pyQuil [35] and ProjectQ [36]. Because some of these platforms include the capability to simulate noisy quantum circuits, we briefly compare their respective noise models and solver types in Table I. Furthermore, for packages that support arbitrary time-dependent Hamiltonian and rely on MEs as solvers, we list their compatible MEs in Table II.  advanced/02-3 qubit entanglement witness 3-qubit entanglement witness experiment meaningful. As a useful example, we choose the alternating-sectors-chain (ASC) [38] as the benchmark problem. The Hamiltonian of the N -qubit experiment is where τ = t/t f is the dimensionless time and a(s) and b(s) are the annealing schedules shown in Fig. 1(a). The alternating-sectors-chain Hamiltonian is where the coupling strength J i alternates between sectors of size n To keep the problem size manageable, we fix n = 1 and vary the system size N . The open-system model is given by In (b), we show the annealing parameters s(τ ) and sp(τ ) (see Eq. (22)) used in the entanglement witness experiment. There are three stages in the experiment: i. first evolve s(τ ) from 1 to s * for τ ∈ [0, τ1] and then evolve sp(τ ) from 1 to s * p for τ ∈ [τ1, 2τ1]; ii. pause for a time τ2; iii. reverse the first stage. In our simulation, we choose s * = 0.339, τ1 * t f = 10µs [37] and s * p = 0.612 such that 2A(s * p ) ≈ 1MHz [24]. The value of τ2 is varied to obtain the tunneling rate.
where each qubit couples to an independent bath via σ z with equal coupling strength g, and H B is the bath Hamiltonian. The bath is chosen to be Ohmic [Eq. (12)] with coupling strength ηg 2 S /h 2 = 1.2 × 10 −4 , cutoff frequency f c = 4GHz, and temperature T = 12mK [38]. In the benchmark simulation, we solve the AME (frequency form Redfield equation) using both HOQST and QuTiP. Because of the large computational cost, the full AME simulation can hardly scale beyond a few qubits. HOQST provides interfaces to solve the AME in a low energy subspace. It can greatly speed up the computation if the evolution is confined within a small, low energy subspace. We also include this version of the AME solver in our benchmark (see Supplementary Note 6 for details). Finally, we ignore the Lamb shift in all the simulations since QuTiP does not include it (a significant drawback since in general the Lamb shift can have a strong effect [30]).
The benchmark result is shown in Fig. 2, where a significant runtime improvement of HOQST over QuTiP is observed. In fact, the QuTiP runtime became excessive for N > 4, while HOQST remains nearly an order of magnitude faster even for N = 5. Using HOQST's subspace truncation capability allowed us to continue simulations up to N = 10 without exceeding a runtime of 10 3 s. However, we note that we deliberately chose a benchmark problem that would demonstrate HOQST's advantage. There are areas of overlap between QuTiP and HOQST (e.g., the time-independent Lindblad equation) where the packages would perform similarly.
As a final remark, we emphasize that HOQST simulates the full open-system dynamics in the sense that the solution of the master equation is obtained up to a precision allowed by the underlying ODE algorithm. The computational cost of such a simulation scales exponentially with respect to the system size, and, without further assumptions or approximations, no algorithm with better scaling has been discovered. On the other hand, quantum Monte Carlo or tensor-network based algorithms could achieve superior scaling under additional assumptions or approximations. For example, if we know that the state stays within the space of matrix product state during the evolution, and the Liouvillian superoperator of the master equation could be effectively expressed in terms of matrix product operators, a tensor-network version of the ODE algorithm can solve the open-system dynamics efficiently. However, generically such assumptions are hard to satisfy, even approximately. Thus, it is not meaningful to benchmark HOQST against tensor-network based algorithm since the latter involves much stronger assumptions.

Entanglement witness experiment modeling
HOQST has a large collection of tutorials located at a dedicated Github repo , which are summarized in Table III. As an illustrative yet non-trivial example, we next discuss the simulation of a three-qubit quantum annealing entanglement witness experiment. The entanglement witness experiment was proposed to provide evidence of entanglement in a D-Wave quantum annealing device [24]. An open system analysis of these experiments was performed using the AME [37], but failed to reproduce the observed width of the tunneling rate peaks. This is the impetus for us revisiting this experiment here. As we shall show, the new tools provided in HOQST allow us to much more closely match the experimental data than was possible before.
FIG. 2. Total runtime vs system size for the benchmark simulation of the alternating sectors chain problem. The x-axis (N ) is the system size. Each data point corresponds to the average runtime of 7-10 runs. Error bars represent 5 standard deviations. The "QuTiP" and "HOQST (Full)" points represent the runtime of full adiabatic master equation (AME) simulations using the corresponding packages. The "HOQST (Subspace)" points represent the runtime of AME simulation in the lowest 20-level subspace. The benchmark was done on a desktop computer with an Intel(R) Core(TM) i7-6700 @3.40GHz CPU and 16 GB memory. The software versions were QuTiP 4.6.2, HOQST 0.6.3, python 3.9.7 and Julia 1.6.3. The operating system was Ubuntu 20.04.3.
The crux of the experiment is actually a form of tunneling spectroscopy [39], where the goal is to find the energy gaps of the Hamiltonian −a(s) i σ i x + b(s)H Ising . This is done by observing the location of a peak in the tunneling rate as measured using a probe qubit. The Hamiltonian of the 3-qubit-version of the experiment is where a(s) and b(s) are the annealing schedules, and s(τ ) and s p (τ ) are functions of the dimensionless time τ = t/t f , known as annealing parameters. The Hamiltonian consists of two system qubits coupled to an ancilla system qubit, as shown in Fig. 3. The aforementioned location of the tunneling rate peak can be controlled by varying h P , and this information can be used to extract the energy gaps as a function of s. We refer interested readers to Refs. [24,37] for more information. 1 The Ising Hamiltonian of the 3-qubit entanglement witness experiment. The figure shows the graphical representation of HIsing in Eq. (22). Here each circle represents a qubit (qubit 1, 2 and an ancilla denoted as P). hi is the local field strength and Jij is the coupling strength between qubit i and j. The Hamiltonian can be written as i hiσ z i + ij Jijσ z i σ z j . The goal of the experiment is to demonstrate entanglement between qubits 1 and 2, by probing the ancilla qubit. In our simulations, h1, h2, J1P and JS are fixed at J1P = −h1 = −1.8, JS = −2.5, h2 = 0.
The annealing schedules and annealing parameters are illustrated in Fig. 1. To extract the tunneling rate, we first perform the simulation with the initial all-one state |ψ(0) = |1 ⊗3 for different h P and t 2 = τ 2 t f values. The population of the all-one state at the end of anneal is then obtained as a function of h P and τ 2 : Lastly, we fit P |1 (h P , t 2 ) to the function ae bt2 + ce dt2 , from which the rate Γ can be estimated For the open system model, we assume the qubits are coupled to independent baths [37] but the bath coupling to the probe qubit g P is much stronger than the coupling to the two system qubits g 1 = g 2 = g S . In addition, we assume the bath is Ohmic [Eq. (12)] with coupling strength ηg 2 S /h 2 = 1.2732 × 10 −4 , cutoff frequency f c = 4GHz, and temperature T = 12.5mK. We performed numerical simulations with different models of B P : 1. An Ohmic bath with interaction strength g P = 10g S , using different flavors of the AME (see Methods, Sec. Adiabatic master equation).
2. Hybrid Ohmic bath whose coupling strength to the Ohmic component is g P = 10g S and varying macroscopic resonant tunneling (MRT) width, using the PTRE (see the Polaron Transform subsection in the Methods).
The tunneling rates obtained via these different ME simulations are compared with the experimental results [24] in Fig. 4. The reported experimental parameters are T = 12.5mK and s * p = 0.612 [24][see Fig. 1(b)]. Simulations using these parameters and two different flavors of the AME, the one-sided AME [Eq. (42)] and the Lindblad form AME [Eq. (38)], are plotted. The results demonstrate that these two AME flavors only differ significantly near the small gap region, but neither one matches the experimental results. No further improvement is observed by varying the AME parameters: the linewidth remains too narrow to match the experiment.
In contrast with the AME, the PTRE with T = 12.5mK and s * p = 0.612 exhibits a larger Gaussian linewidth broadening and closer agreement with the experimental data (solid curves in Fig. 4). We note that if we increase W while fixing T , the PTRE curve is stretched to the right. This is the result of the fluctuation-dissipation theorem, where ε L scales quadratically with W . Such a shift can be compensated by increasing the temperature together with W [see the black dashed curve in Fig. (4)]. Despite the closer agreement using the PTRE at the reported experimental temperature and s * p values, there is still a mismatch between the theoretical and experimental amplitudes of the tunneling rate curves. A possible reason for this may be a discrepancy in the reported A(s * p ) and its true value, attributable to annealing schedule fluctuations and integrated control errors in the early model of D-Wave annealer used in the experiment. To account for this, we performed PTRE simulation with different T and s * p values from those reported [24]: T = 25mK and s * p = 0.59 (we also set W = 9mK). The result is plotted as the black dashed line in Fig. 4 and shows significantly better agreement with the experimental results. This highlights the fact that the slow bath coupled to the probe qubit may have a different temperature than the Ohmic one. Moreover, this results illustrates the power of HOQST's range of ME implementations.

CONCLUSION
In conclusion, we presented a software package called Hamiltonian Open Quantum System Toolkit (HOQST). It is user-friendly and written in Julia. It supports various master equations with a wide joint range of applicability, as well as stochastic Hamiltonians to model 1/f noise. We demonstrated that HOQST can achieve order of magnitude speedups over QuTiP for problems with a time-dependent Hamiltonian. We also illustrated the use of HOQST in simulating open quantum system dynamics in a 3-qubit entanglement witness experiment. Whereas previous modeling of this experiment was unable to capture the reported linewidth, HOQST's implementation of the polaron-transformed Redfield equation (PTRE) was able to do so. We also derived new error bounds for the PTRE.
We expect HOQST to be useful for researchers working in the field of open quantum systems, dealing with systems governed by time-dependent Hamiltonians. HOQST provides both basic and advanced numerical simulation tools in this area, which can be applied to simulate superconducting qubits of all types, trapped ions, NV centers, silicon quantum dot qubits, etc. Future releases of HOQST will expand both the suite of open system models and range of quantum control and computation models it supports.

METHODS
We present a brief review of the open system models and corresponding master equations supported by HOQST. We provide various additional technical details in the Supplementary Information.

CUMULANT EXPANSION
The cumulant expansion is a technique originally designed for the perturbation expansion of stochastic differential equations [40]. This technique can be generalized to the open quantum system setting and allows a systematic description of the reduced system dynamics [1]. By applying this technique in different rotating frames, master equations with different ranges of applicability can be derived [1,18,41]. Defining the projection operator P the formal cumulant expansion of the Liouville von Neumann equation is where the nth order generator K n (t) is and the quantities L (t)L(t 1 ) · · ·L(t n−1 ) oc are known as ordered cumulants [1]. In HOQST, we consider only the first and second order generators, which are given by: Incorporating higher order cumulants generally leads to more accurate results [1].

Redfield equation
The oldest and one of the most well-known MEs in this category is the Redfield equation [42] (also known as TCL2 [1], where TCLn stands for time-convolutionless at level n, arising from an expansion up to and including K n (t)), which directly follows from Eq. (26) after choosing the rotation U (t) to be where T + denotes the forward time-ordering operator. After rotating back to the Schrödinger picture, one of the most common forms of the Redfield equation iṡ where This is the form used in HOQST. The error bound for the Redfield ME is given in Eq. (14). We note that the current release of HOQST supports correlated baths for the Redfield and adiabatic master equation solvers. However, for simplicity, we henceforth focus on uncorrelated baths where C αβ (t) = δ αβ C α (t).
The most significant drawback of the Redfield equation is the fact that it does not generate a completely-positive evolution, and in particular can result in unphysical negative states (density matrices with negative eigenvalues). Though formal fixes for this problem have been proposed [43,44], to address the issue in HOQST an optional positivity check routine is implemented at the code level, which can stop the solver if the density matrix become negative. In addition, three variants of Redfield equation that guarantee positivity, namely the adiabatic master equation (AME) [30], the coarse-grained master equation (CGME) [27,45], and the universal Lindblad equation (ULE) [28], are included in HOQST. We detail these MEs next.

Coarse-grained master equation (CGME)
The CGME can be obtained from Eq. (30) by first time-averaging the Redfield part, i.e., shifting t → t + t 1 and applying 1 Ta Ta/2 −Ta/2 dt 1 . One then neglects a part of the integral to regain complete positivity. The result is [27]: where A α (t + t 1 ) = U † (t + t 1 , t)A α (t)U (t + t 1 , t) and the Lamb shift is given by The quantity T a is the coarse-graining time, a phenomenological parameter that can be manually specified or automatically chosen based on the bath correlation function [45]. HOQST uses a multidimensional h-adaptive algorithm [46] to perform the 2-dimensional integration. The error bound for the CGME is given in Eq. (16).

Universal Lindblad equation (ULE)
The ULE [28] is a Lindblad-form master equation that shares the same error bound as the Redfield equation, i.e., Eq. (14). A similar master equation with better accuracy, known as the geometric-arithmetic master equation (GAME), can also be derived by using a different formula for the Lamb shift [47]. The formal form of the ULE is identical to the Lindblad equation: where the time-dependent Lindblad operators are and the Lamb shift is In the above expression, g α (t) is called the jump correlation and is the inverse Fourier transform of the square root of the noise spectrum The integration limits of Eqs. (35) and (36) are problematic in practice because the unitary U S (t) does not go beyond [0, t f ]. In numerical implementation, we replace the integral limit with t f 0 . This is a good approximation when g(t) decays much faster than t f . This is the form of ULE used in HOQST.

Adiabatic master equation (AME)
To derive the AME, we replace U S (t − τ ) in Eq. (31) with the ideal adiabatic evolution and apply the standard Markov assumption and the rotating wave approximation (RWA). The resulting equation iṡ The AME is in Davies form [29] and the Lindblad operators are defined by where ε a is the instantaneous energy of the a'th level of the system Hamiltonian, i.e., H S (t) |ψ a (t) = ε a (t) |ψ a (t) . Finally, the Lamb shift term is where with P denoting the Cauchy principal value. The error bound is given in Eq. (15). If the RWA is not applied, the resulting equation is called the one-sided AME: where These two forms of the AME behave differently when the energy gaps are small because the RWA breaks down in such regions [27]. More importantly, like the Redfield equation, the one-sided AME does not generate a completely-positive evolution. The code-level positivity check routine works with this version of the AME as well.
Classical 1/f noise HOQST includes the ability to model 1/f noise, which is an important and dominant source of decoherence in most solid-state quantum NISQ platforms [48,49], in particular those based on superconducting qubits [50][51][52]. Fully quantum treatments of 1/f noise have been proposed [53,54], including for quantum annealing [18]. In HOQST we adopt the simpler approach of modeling 1/f noise as classical stochastic noise generated by a summation of telegraph processes, which has proved to be a good approximation to the fully quantum version [48]. Specifically, we provide a quantum-trajectory simulation of the following stochastic Schrödinger equation where each δ α (t) is a sum of telegraph processes where T i (t) switches randomly between ±b i with rate γ i . In the limit of N → ∞ and b i →b, if the γ i 's are loguniformly distributed in the interval [γ min , γ max ] (with γ max γ min ), the noise spectrum of δ α (t) approaches a 1/f spectrum within the same interval [48]. Empirically, we find that a good approximation can be achieved with relatively small N .

Hybrid model
The most significant drawback of a purely classical noise model is that if its steady state is unique then it is the maximally mixed state. To see this, we first realize that each trajectory of Eq. (44) generates a unitary acting on the space S(H S ) of density matrices Averaging over the trajectories over a distribution p(k) creates a unital (identity preserving) map from S(H S ) into itselfŪ If the steady state ρ ∞ is unique then we can define it as By unitality it would then follow that ρ ∞ = I, since we can choose ρ S (0) = I. However, this is not what is observed in real devices, e.g., in experiments with superconducting flux [16] or transmon [55] qubits. To account for this, HOQST includes a hybrid classical-quantum noise model: where δ α (t) is the same random process as in Eq. (45), and L is the superoperator generated by the cumulant expansion (27). At present, HOQST supports the combination of 1/f noise with both the Redfield and adiabatic master equations.

Polaron transform
If the bath operators in Eq. (2) are bosonic we can choose the joint system-bath unitary U (t) in Eq. (29) as [41] U p (t) = exp where A d α is the diagonal component of A α in the interaction Hamiltonian (2) and U B (t) is given in Eq. (29) (we use λ instead of g in B α to distinguish it from the expansion parameter in Eq. (7)). The corresponding second order ME (30) is known as the polaron-transformed Redfield equation (PTRE) [41] or the noninteracting-blip approximation (NIBA) [3]. The PTRE has a different range of applicability than the previous MEs we have discussed. Whereas the latter apply under weak-coupling conditions, the transformation defined in Eq. (51) leads to a complementary range of applicability under strong-coupling. This particular form of Eq. (51) does not preserve the factored initial state, so that inhomogeneous terms are present after the transformation. However, if ρ S (0) is diagonal then numerical studies of the effects of the inhomogeneous terms suggest that they can be ignored [31].
In addition, the PTRE can be extended beyond the spin-boson model by choosing a different form of the joint system-bath unitary (51), as [18,56] The two transformations in Eq. (51) and (52) lead to MEs with identical structure but slightly different expressions (see Supplementary Note 1 for details). Whether those differences make any physical significance is an interesting topic for further study. In this paper we choose to work with Eq. (51). Because the general form of the PTRE is unwieldy, we present its form for a standard quantum annealing model where a(t) and b(t) are the annealing schedules, and H driver and H prob are the standard driver and problem Hamiltonians, respectively: where the Pauli matrix σ x acting on qubit i is denoted by σ x i , etc. The transformed Hamiltonian is where The Redfield equation corresponding to Eq. (55) is where Here K αβ i (t, τ ) is the two-point correlation function in the polaron frame [akin to the correlation function defined in Eq. (9)], and κ i corresponds to the first order cumulant generator in Eq. (27) and is also known as the reorganization energy; it contributes a Lamb-shift-like term in Eq. (57). It is also worth mentioning that the polaron transformation (51) can be done partially, which means that in Eqs. (51) and (52), α can be summed over a subset of system-bath coupling terms.
To solve this form of the PTRE in HOQST, the user can define a new correlation function C αβ i (t, τ ) = a(t)a(τ )K αβ (t, τ ) and use the Redfield solver. An alternative approach is to make the Markov approximation in Eq. (58a) t 0 a(τ ) · · · dτ → a(t) ∞ 0 · · · dτ (60) and write Eq. (57) in Davies form [29] (see Supplementary Note 4 for more details). This leads to the same expression as the AME [Eq. (38)], but with different Lindblad operators: where now |ψ a is the energy eigenstate of the HamiltonianH S (t), and the noise spectrum is Then the AME solver can be used to solve this Lindblad-form PTRE. For example, the following ME can be derived for the entanglement witness problem following the aforementioned procedure : whereH and the Liouville operators L A and L P corresponds to the AME part and PTRE part of this equation, respectively: where the Lindblad operators are defined in Eq. (39) and Eq. (61) respectively. The function γ(ω) is the standard Ohmic spectrum and γ P (ω) is the polaron frame spectrum with a hybrid Ohmic form [18,56] discussed in Supplementary Note 5. We provide the explicit form of γ P (ω) here where and G L (ω) is the contribution of low frequency component, characterized by the MRT width W . Because W and ε L are connected through the fluctuation-dissipation theorem we have W 2 = 2ε L T ; thus, hybridizing low frequency noise with an Ohmic bath introduces one additional parameter.

Redfield backward integration
To solve the Redfield or Redfield-like master equation (31), one needs to integrate the unitary U S backward in time at each ODE step. Such integrations are computationally expensive for long evolution times and become the bottleneck of the solver. To improve the efficiency of the solver, we introduce an additional parameter T a as the lower integration limit: To justify this, note first that where · is any unitarily invariant norm. To obtain this inequality, we perform a change of variable t − τ → τ and make use of the fact that the operator A β (t) can always be normalized by absorbing a constant factor into the corresponding bath operator B β . Second, note that in most applications the bath correlation function C(τ ) decays fast compared with the total evolution time. As a result, the r.h.s. of Eq. (70) is small for sufficiently large t. The neglected part, i.e., the integral over [0, T a ], can thus be safely ignored so long as the r.h.s. of Eq. (70) is below the error tolerance of the numerical integration algorithm.
The same technique can also be applied to the ULE. The integration limits in Eqs. (35) and (36) can be localized around t, i.e. replaced by t+Ta t−Ta . However, choosing an appropriate T a is a process of trial and error. The user needs to determine its value in a case-by-case manner.

Precomputing the Lamb shift
Instead of evaluating the Lamb shift (41) at each ODE step, to speed up the computations all the ME solvers support precomputing the Lamb shift on a predefined grid and use interpolation to fill up the values between the grid points.

Adiabatic frame
For a typical annealing process, the total annealing time is usually much larger than the inverse energy scale of the problem t f 1 min s∈0,1 max(A(s), B(s)) .
Informally, the frequency of the oscillation between the real and imaginary part of the off-diagonal elements of ρ S in the neighborhood of s is positively proportional to both A(s) and B(s). As a consequence, directly solving the dynamics in the Schrödinger picture is challenging because the algorithm needs to deal with the fast oscillations induced by the Hamiltonian, thus impacting the step size. HOQST includes an optional pre-processing step to rotate the Hamiltonian into the adiabatic frame [57]. If the evolution is in the adiabatic limit, the off-diagonal elements of the density matrix in this frame should approximately vanish. The fast oscillation is absent and a large step size can be taken by the ODE solver. This technique provides advantages if the user wishes to repeatedly solve the same problem with different parameters. See Supplement Method 1: Adiabatic frame for a brief summary.

Quantum trajectories method
HOQST implements a quantum-trajectory solver for the AME [17]. Using the native distributed memory parallel computing interface of both Julia and DifferentialEquations.jl, the quantum-trajectory simulations can take advantage of HPC clusters with minimum changes in the code. In addition, classical 1/f noise can be infused into the AME trajectory solver to generate the hybrid dynamics described in Eq. (49).

Two different formulations
We illustrate two different formulations of the polaron transformation via the simplest Hamiltonian The standard polaron transformation rotates this Hamiltonian via the unitary [41] which leads to the interaction picture Hamiltoniañ where The only quantities that matter in the 2nd order TCL formalism are the average and two-point correlation function of ξ ± (t) given by where α, β ∈ {+, −} and U B (t) = exp{−iH B t}. To arrive at the second equality in Eq. (76a) we used [ρ B , U B ] = 0, which is true for ρ B in a Gibbs state. An alternative approach [18] is to first rotate the Hamiltonian (72) w.r.t. H B and then rotate it again by where It is worth noting that in this formalism the bath does not need to be bosonic. Let us define Then the effective Hamiltonian in this two-fold interaction picture is where Again, the 2nd order ME only depends on whereᾱ means the opposite symbol of α in {+, −}. Appendices D.4 and D.5 in [18] give detailed calculations for Eq. (82a) and (82b).

Correlation function in polaron frame
We now discuss Eqs. (76a, 76b, 82a, 82b) for the spin-boson model. We first introduce the bath spectral function [58] which is usually taken into the continuous limit. The standard result [58,59] for Eqs. (76a) and (76b) are and where For spectral functions which lead to κ → 0 and φ(t) → ∞ (e.g., an Ohmic bath), κ 2 e φ(t) converges to a finite number. Thus the standard result for an Ohmic bath is where Q 2 (t) and Q 1 (t) [58] are By substituting the Ohmic spectral function J(ω) = ηωe ω/ωc , the explicit form of Q 1 and Q 2 can be calculated. Here we take another approach by noticing that where C(t 1 , t 2 ) is the bath correlation function (defined in Eq. (9) of the main text) that can also be expressed in terms of the spectral function Thus Eq. (87b) can be rewritten as where γ(ω) is the spectral density defined in Eq. (11) of the the main text. For Eq. (82a) and (82b), we follow Ref. [18] and write down the the result here for comparison. The 1st order average is which is supposed to decrease exponentially fast for t > 0. The two-point correlation functions are The difference between Eqs. (87b) and (93b) are easily noticeable. The physical origins and significance of this difference is a subject for future studies.

SUPPLEMENTARY NOTE 2: ERROR ANALYSIS
We present an error analysis based on the standard form of the polaron transformation (73). Our strategy is to divide the error into separate parts and use Lemma 1 in Ref. [27] to estimate the total error of the master equation. Besides those already presented in [27], there are two additional error terms in the PTRE, the physical origin of which are the inhomogeneous terms resulting from the breakdown of the factorized initial condition by U p : In addition to the projector defined in Eq. (21) in the main text, we define its complementary Q ≡ I − P. The standard projection operator technique [1] leads to where the rightmost expression is the contribution of the inhomogeneous terms. Ref. [60] calculates I n (t) up to n = 4 and shows they are identical to those for K n (t) except that the final P in each term is replaced by a Q. Here we conjecture that this is also true for n > 4, and then use the same error bound as in Ref. [27] (called the Born approximation error there), to bound the error of truncating to the order of I 2 : To simplify the notation, we define ρ true (t) as the solution of Eq. (95) when n → ∞. The PTRE error bound can be written as the sum of three terms: where E BM is the error bound presented in Ref. [27]: 1 where the second term automatically satisfies the bound given in Eq. (98). To bound the first term, we first write out ρ(0) explicitly for the case of a single system qubit and where U p was defined in Eq. (47) in the main text, and Using our conjecture, the formal expansion of I n (t) is the same as K n (t) with the average operation in the rightmost term being replaced by After substituting Eq. (100c) into the above expression and expanding the commutators, one term in the summation is where we use +, − in the subscript instead of +1, −1 to simplify the notation. Thus, the next step is to calculate the average and two point correlation function under the new "average" operator (102): Eq. (105a) can be explicitly carried out using the following identities: and the Baker-Campbell-Hausdorff (BCH) formula After some algebraic manipulations, we get where f α n (t) = exp 4iαn k . For simplicity, we consider only the Ohmic bath case, for which ξ ± (t) = 0 and ξ ± (t)ξ l (0) decay exponentially for t > 0. As a result, we can ignore the first order term ξ ± (t) mn ≈ 0. The two-point correlation function is .
(109) where the second line equals 0 because ξ ± (t) = 0 and the noise is Gaussian. The above result can be generalized to the multi-point correlation function .
(110) For the m = n case, the correlation function under · mn is the same as the correlation function under · up to an additional phase factor. For the m = n case, the right hand side of Eq. (110) can be decomposed into two-point correlation functions by Wick's theorem. After the decomposition, each term would appear inside the integral as Again, because the particular pair ξ αs (t s )ξ l (0) decreases exponentially for t s > 0 and thus has approximately zero overlap with other two point correlation functions under the integral, we ignore all these terms. Finally, we can sum up every term over the indices m and n. For example, for terms that have the same form as Eq. (104), the summation is . Eq. (112) can be bounded by where we make use of the fact |õ| ≤ |ρ 00 S | + |ρ 11 S | = 1. Formally applying the same technique to every term after expanding I n (t)ρ(0), allows us to show that the first part of Eq. (99) also satisfies the error bound given in Eq. (98). 2 In conclusion, the error due to truncation of the inhomogeneous parts to second order has the same scaling as the error due to truncation of the homogeneous parts [Eq. (98)]. As a result, the total truncation error can be combined using a single big-O notation: We discuss the timescales in the polaron frame in the next section. Before proceeding, we remark that: 1. The same analysis also applies to the multi-qubit PTRE described in Eq. (53) of the main text by replacing ρ ij where k is the index for different qubits.
2. The analysis is also applicable for time-dependent ∆(t). In this case, the corresponding two-point correlation function is defined as C αβ (t 1 , t 2 ) = ∆(t 1 )∆(t 2 ) ξ α (t 1 )ξ β (t 2 ) , which can be upper bounded by where ∆ m = max t∈[0,t f ] ∆(t). This introduces a constant ∆ m in the definition of the timescales.
3. The error bound we analyzed in this section does not include the error due to ignoring the first and second order inhomogeneous terms E (2) I (which is usually the case in the standard PTRE treatment). Numerical studies show that they can be ignored if ρ S (0) is diagonal [31]. More rigorous error bounding is an interesting topic for future study.

SUPPLEMENTARY NOTE 3: TIME SCALES
The time scales given in Eq. (13) in the main text can be computed directly for the polaron frame two-point correlation function (87b). Because the error bound in Eq. (14) scales with 1/τ SB and τ B /τ SB , we explicitly write out those two quantities where ∆ m = max t∈[0,t f ] ∆(t). We can see that, when the system-bath coupling strength is sent to infinity, both of these quantities go to zero: 1/τ SB → 0 and τ B /τ SB → 0 as g → ∞. Thus, the PTRE works in the strong coupling regime.

SUPPLEMENTARY NOTE 4: LINDBLAD FORM OF THE PTRE VIA THE ADIABATIC APPROXIMATION
In this section, we derive an adiabatic form of the PTRE from Eq. (53) of the main text: where andŨ Three approximations need to be applied to obtain the final result:

Markov approximation
The first step is to move the function a(t) in Eq. (119a) outside the integral then perform a change of variable t − τ → τ . Eq. (119a) becomes: Finally, the integration limit is taken to infinity:

Adiabatic approximation
A detailed description of the adiabatic approximation is given in Ref. [30]. The core idea is to rewrite the unitary in Eq. (122) asŨ and then replaceŨ S with an appropriate adiabatic evolutioñ where U ad S (t, 0) is the ideal adiabatic evolution with a phase In the above expressions, |ε a (t) is the instantaneous eigenstate ofH S (t) and φ a (t) = i ε a (t) |ε a (t) is the Berry connection. After this procedure, Eq. (122) becomes: where ω ba = ε b − ε a , and L ab iβ (t) = ε a (t)|σ β i |ε b (t) |ε a (t) ε b (t)| .

Rotating wave approximation(RWA)
To perform the RWA, we first need to move the one-sided adiabatic PTRE into the interaction picture with respect toH S (t) via the approximated unitary in Eq. (124b). The non-Hamiltonian part of Eq. (118) is transformed into where Π ab (0) is a shorthand notation for |ε a (0) ε b (0)|. The RWA amounts to keeping only terms where either a = b, b = a or a = b, a = b . After the fast-oscillating terms are ignored, we move back to the polaron frame and make use of Eq. (39) of the main text to derive the Lindblad-form PTRE: where the new Lindblad operators are given by The notationᾱ again means the opposite symbol of α in {+, −}. Since σ ± i are Hermitian conjugates of each other: Finally, the Lamb shift term isH Because the techniques used here are the same as those in [27,30], the error bounds in these references are still valid.

SUPPLEMENTARY NOTE 5: HYBRID NOISE
The two point correlation function in the polaron frame [Eq. (91)] can be extended to a hybrid noise model [56]. The core assumption here is that the bath spectral density can be separated into low frequency and high frequency parts: The integral inside the exponent of Eq. (91) can also be separated into If γ L (ω) is concentrated on low frequencies, a 2nd order Taylor expansion of e −iωt is justified. As a result, the first term in Eq. (136) can be written as where and P stands for the Cauchy principal value. W and ε L , usually known as the MRT linewidth and reorganization energy, respectively, are experimentally measurable quantities that are connected through the fluctuation-dissipation theorem [56]: W 2 = 2ε L T . Finally, the two point correlation function under this hybrid noise model is The corresponding spectral density of K(t) is where G L (ω) and G H (ω) are the Fourier transforms of exp{4f L (t)} and exp{4f H (t)}, respectively. The low frequency component is Gaussian: while the high frequency part can be approximated as a single Lorentzian [18]: Finally, it should be mentioned that the fluctuation-dissipation theorem guarantees that the Kubo-Martin-Schwinger (KMS) condition be satisfied for G L (ω). So if γ H (ω) satisfies KMS condition, the spectral density in the polaron frame also satisfies the KMS condition γ K (|ω|) = e β|ω| γ K (−|ω|) . (143)

SUPPLEMENTARY NOTE 6: PSEUDOCODE FOR TRUNCATED-SUBSPACE AME
In this section we present the pseduocode for solving AME in the truncated-subspace. First, we define two sets of grid points {t i } N i=1 and {τ i } N +1 i=1 such that τ 1 < t i < τ 2 < t 2 · · · τ i < t 1 < τ i+1 < · · · < t N < τ N +1 , and τ 1 = 0, τ N +1 = t f . The values of {t i } N i=1 are chosen subject to the following constraint where |v k (t i ) is the kth eigenstate of the system Hamiltonian H S (t i ). Additionally, l is the number of energy levels to use in the simulation, n is the size of the energy subspace where the dynamics is confined and is the numerical error tolerance. This condition ensures the lowest n eigenstates of H S (τ i ) can almost be expressed by the lowest l eigenstates of H S (t i ). Then the simulation can be done in a piecewise fashion described in Algorithm 1. In practice, n is highly problem-dependent and needs to be estimated empirically based on the adiabatic theorem.

Algorithm 1 Truncated-subspace AME
Require: H(t) = k f k (t)M k ρ0(τ1) = ρ(0) ρ(0) is the initial state in the computaitonal basis i ← 1 while i ≤ N do Vi ← eigendecomposition of H(ti) truncated to lowest l levels Each column of Vi corresponds to one eigenvector Hr(t) ← k f k (t)V † i M k Vi An,r ← V † i AnVi An is the nth system-bath coupling operator ρr ← V † i ρi−1(τi)Vi ρi(t) ← solve AME defined by Hr and An,r from τi to τi+1 with initial state ρr i ← i + 1 end while