SchWARMA: A model-based approach for time-correlated noise in quantum circuits

Temporal noise correlations are ubiquitous in quantum systems, yet often neglected in the analysis of quantum circuits due to the complexity required to accurately characterize and model them. Autoregressive moving average (ARMA) models are a well-known technique from time series analysis that model time correlations in data. By identifying the space of completely positive trace reserving (CPTP) quantum operations with a particular matrix manifold, we generalize ARMA models to the space of CPTP maps to parameterize and simulate temporally correlated noise in quantum circuits. This approach, denoted Schr\"odinger Wave ARMA (SchWARMA), provides a natural path for generalization of classic techniques from signal processing, control theory, and system identification for which ARMA models and linear systems are essential. This enables the broad theory of classical signal processing to be applied to quantum system simulation, characterization, and noise mitigation.


I. INTRODUCTION
The circuit model of quantum computing [1][2][3] has been used extensively to prove algorithmic speedups [4][5][6] as well as the theory of fault-tolerant quantum computing [7][8][9].This model makes an implicit assumption that a quantum computation can be broken up into a series of (possibly imperfect) incremental gates.This assumption generally implies that noise processes are uncorrelated in space and time.However, the presence of time-correlated noise, which is the focus of this Article, is well established in actual physical systems, e.g., 1/f α noise in superconducting qubits [10][11][12][13][14].
Characterization of time-correlated noise processes typically involves non-parametric reconstructions of the entire noise spectra [11,[15][16][17][18][19][20] requiring a large number of experiments.Furthermore, outside of a few closed-form cases, there is no more efficient way of forward simulating the dynamics of time-correlated noise short of brute-force numerical integration of the stochastic Liouville equation [21] for the entire circuit.The time step for such a simulation is typically orders of magnitude smaller than the nominal "gate duration," resulting in a large number of matrix operations per gate.This approach is computationally very expensive, which has prevented much in the way of analysis of time-correlations in quantum circuits beyond a few qubits.
Understanding the impact of time-correlated errors on quantum computation is important as they are known to have an impact on quantum error correction.These errors have been considered in the context of proving fault-tolerant threshold theorems [22][23][24], although exactly how well quantum error correction will suppress noise in the presence of these errors is not well understood, but we do know of some examples where their presence can lead to fairly detrimental affects [25].
In this Article, we introduce a model of temporally correlated noise that allows the wide-ranging field of classical time-series analysis to be applied to noisy circuit simulation and quantum noise spectroscopy.Our model draws on the technique of autoregressive moving average (ARMA) models [26,27], linking open quantum systems to discrete time linear systems theory and model-based spectrum estimation [28].With simulation and estimation of noisy quantum dynamics in mind, we generalize ARMA to the space of completely positive trace preserving (CPTP) maps, which we denote Schrödinger Wave ARMA (SchWARMA).We convey the efficacy and versatility of this method through a set of examples of temporally correlated noise-affected dynamics from quantum control and error correction.
SchWARMA has a natural connection to semi-classical stochastic Hamiltonian noise in quantum systems, but is more general as it can also model correlated errors in dissipative system-baths as in a stochastic form of Lindblad master equation [29].In this sense, the errors considered here are fundamentally Markovian with respect to any quantum bath, as the dynamics produced are correlated yet remain completely positive (CP) divisible [30].Techniques such as [31] that seek to model non-Markovian quantum dynamics are fundamentally different from this noisy-circuit approach, as non-CP error channels interspersed with perfect gates can produce nonphysical states.That said, our technique can be generalized to simulate system-bath models of non-Markovian quantum dynamics at an obvious reduction of the number of qubits simulated, but we leave this to future work.

II. THE SCHWARMA MODEL
An ARMA model relates the output y k , at discrete steps k, to a series of inputs x k and prior outputs via [26,27], The set {a i } defines the autoregressive portion of the model, and {b j } the moving average portion with p and q+1 elements of each set respectively.The power spectrum S y (ω) of the output y of an ARMA model driven by i.i.d.Gaussian noise is arXiv:2010.04580v2 [quant-ph] 12 Aug 2021 and ARMA models can approximate any discrete-time power spectrum to arbitrary accuracy [32].ARMA models (and their many generalizations [32]) are commonly used in classical time series analysis to parametrically define, estimate, and simulate time-correlated processes [27], and have wideranging applications all the way from actuarial science [33] to zoology [34].We generalize ARMA models to CPTP maps by using Stiefel manifolds [35], spaces of matrices with orthonormal columns.The tangent space of these manifolds enable a natural connection to the mathematics of ARMA on Lie algebras (i.e., the tangent space to a Lie group) formally developed in [36] and to the Lindblad master equation as shown in App. A. Using these relationships, we build a theory of ARMA relevant for quantum information processing that models temporally-correlated noise in a circuit-model formalism.
We begin by showing how CPTP maps are related to Stiefel manifolds.A CPTP map Φ operating on N × N density operators ρ can be expressed by Φ(ρ) = k M k ρM † k where M k are N × N complex matrices called Kraus operators [37] with k M † k M k = I N .This form is not unique, but any CPTP map can be expressed using no more than N 2 Kraus operators [37].Given K Kraus operators {M k }, we define , a KN × N matrix with orthonormal columns.
Complex n × m matrices (n ≥ m) with orthonormal columns define a Stiefel manifold [35], denoted V m (C n ).Complex Stiefel manifolds are homogeneous spaces of unitary matrix groups, , where U (n) denotes the n × n unitary matrix (Lie) group.Elements of V m (C n ) are associated with equivalence classes of unitary matrices with the first m columns identical.We note there is a known correspondence between Kraus operators and columns of unitary matrices [38][39][40][41].Thus, the KN × N matrix S ∈ V N (C KN ) relates Kraus operators to points on the corresponding Stiefel manifold.In order to describe perturbations around this point on the manifold, we introduce the tangent space of V N (C KN ).The tangent vectors are given by matrices X = SA + S ⊥ B with with A N × N skew-Hermitian, B (KN −N )×N arbitrary, and S ⊥ denotes an orthogonal complement to S such that [S S ⊥ ] is unitary [42].An exponential map at S maps X back to the manifold via where I n,m denotes the first m columns of the n × n identity matrix I n , and 0 is a zero-matrix of appropriate dimension [42,Eq. (2.42)].Therefore, this map provides a mechanism to generate new elements of V N (C KN ) for specific choices of matrices A, B and S ⊥ .For quantum information applications, we typically describe a noisy quantum channel as a perturbation around some desired unitary operation U .This corresponds to choosing S = [U † , 0] † , and provides a natural choice for S ⊥ = [0, I KN −N ] † .The exponential map in Eq. ( 3) generates the perturbed quantum channel.This form of the map highlights a natural geometric interpretation: B is a block matrix of K − 1 individual N × N matrices B k , so r nonzero blocks result in CPTP maps with Kraus rank K = 1+r, except periodically at the initial point U .This connects to the Linblad master equation through the infinitesimal limit of Eq. ( 3) that yields the time-evolution of the system density matrix (see App.A for more details), where H = −A/i.We now generalize ARMA techniques to the evolution of CPTP maps.Consider a combination of L independent ARMA models with an upper index where y ( ) k denotes the output of the th model at time step k, and so on for a k are independent zeromean Gaussian random variables with unit variance (other distributions may be used, as was done in [25]).Let X ( ) denote fixed elements of the tangent space of U k in V N (C KN ) corresponding to the target unitary U k .Then, we define the CPTP map output S k of a SchWARMA model at time k by direct generalization of Eq. (2.20) in [36] from Lie Groups to Stiefel manifolds where for ideal unitary U , its corresponding Stiefel representation U k , and tangent space element with A an N × N skew-symmetric complex matrix, and B an KN × N arbitrary complex matrix.The S k can then be composed as CPTP maps to implement correlated noisy quantum dynamics.
The SchWARMA model provides a coarse-grained representation of time-correlations in a quantum channel.Specifically, it describes noisy quantum dynamics at discrete time steps (indexed by k in Eq. (5) and subsequent examples), such as after circuit gates.To our knowledge, the only other way to do this is via Suzuki-Trotter methods [43][44][45][46] to numerically solve the Liouiville equation.This requires step sizes much smaller than the fluctuation time-scale of the Hamiltonian such that ∆t ||∂H(t)/∂t|| −1 [47].For high-bandwidth noise or control pulses, this can require many time steps to accurately simulate the trajectory.SchWARMA, on the other hand, estimates the impact of time-dependent noise as an error channel after a perfect gate, with simulation error proportional to the infidelity squared when the noise does not commute with the desired gate.
By separating control and noise, even in the presence of time-correlations, SchWARMA can provide significant advantage for simulating noisy quantum systems with minimal impact to simulation accuracy.In the following, we discuss examples using SchWARMA models for dephasing and amplitude damping channels, depolarizing channel, as well as continuously driven quantum systems.Details of the error analysis are provided in App. C. Furthermore, we note that this model-based formulation also offers theoretical paths to mitigate time-correlated noise.

III. SPECIFIC SCHWARMA MODELS
Here we provide details of how we implement the SchWARMA models for the various examples we provide in the next section.First, consider an example where we wish to simulate correlated noise within a quantum circuit such as the example shown in Fig. 1 that contains arbitrary one and two qubit gates specified as G i .For the noise model, we create two SchWARMA models (denoted by an additional upper index on S k ) given by Eq. (5b) with U = I 2 chosen to be the identity as we want to sample our error channel near the identity.This allows one to simply insert random noise channels into the circuit with the amplitude of the noise determined by the ARMA coefficients given in Eq. (5a).Specific choices of noise channel are considered next.

G3 S
(1) 2 G4 S (1) 3 FIG. 1. Circuit schematic demonstrating how one might implement a SchWARMA model within a circuit of arbitrary one and two qubit gates labelled with Gi.It allows one to consider perfect gates followed by noise that retains the time correlation.The S (q) k terms are the error channel defined in Eq. (5b) with additional upper index denoting the qubit.

A. Z Dephasing
We define a Z dephasing SchWARMA model with: so N = 2 and K = 1, corresponding to unitary dynamics on a qubit.Since K = 1, we have that exp U k is the usual matrix exponential and this gives the map to apply at circuit location k by the simple expression where we have suppressed the index, since there is only a single SchWARMA model needed.This is a Z rotation by a random angle given by a classical ARMA model.

B. Multi-axis Dephasing
For the slightly more complicated multi-axis noise that we consider for our surface code simulations in the next section, the SchWARMA model is given by: Again, N = 2 and K = 1, corresponding to unitary dynamics, but this time with L = 3 corresponding to multiple underlying ARMA models.Again, we have that exp U k is the usual matrix exponential for qubit unitary dynamics, and this gives the map to apply at circuit location k by the expression = cos(g k ) − iy where , N = 1, and sinc(x) = sin(x)/x.

C. Amplitude Damping
Finally, we consider the SchWARMA model for the nonunitary amplitude damping channel.The tangent space element for that SchWARMA model is : where the notation 0 m×n indicates a matrix of zeros with m rows and n columns, N = 2, and K = 2. Thus S k is given by: yielding an 4 × 2 sized matrix containing the two 2 × 2 Kraus operators.Since the tangent space element X in ( 13) has all zeros in the skew-symmetric A component, in this case the driving ARMA model y k can actually be complex-valued.These two Kraus operators at SchWARMA index k are: We see that this SchWARMA model is nothing more than the standard amplitude damping model, but with a decay coefficient given by the square of the sin of the absolute value of the classical ARMA model coefficient.This SchWARMA model simulates an amplitude damping channel with a correlated, fluctuating decay rate, allowing for the analysis of temporal fluctuations in the decay rate of the qubit [48].However, it does not consider time-dynamics of the system-bath interaction which is still assumed to happen on infinitely fast time-scales.
We note that this process allows for some quantum processes that appear to be non-physical with, for example, negative amplitude damping coefficients.All channels created using this method are CPTP and physical.Since manifolds are only locally linear, when one takes a large step in the tangent space, the curvature of the manifold is not respected.Therefore, a direction that appears to be amplitude damping at the start of the curve will appear like something else relative to some point further along the curve defined by the arc.
Therefore, for |y k | < π/2 we have classic amplitude damping, which includes both the non-unital shift towards |0 as well as anisotropic contractions in the X, Y , and Z dimensions of the Bloch sphere to maintain physicality.For y k = ±π/2 all states get mapped to the |0 state.For y k between π/2 and π, the non-unital shift shrinks back to the origin and the contraction process reverses, except that the X and Y axes are both flipped.At y k = π we obtain the Z gate.This sequence is then reversed for π < y k < 2π.The primary point here is that for non-unitary channels such as the amplitude damping channel, the step size should be small enough such that 0 ≤ |y k | < π/2.For high-fidelity quantum operations, this will never be a concern.

D. Lindbladian Dephasing and Depolarizing
When averaged over random instances, the stochastic unitary SchWARMA models in sections III A and III B above will produce, at each time step k, dephasing channels and depolarizing channels respectively, where the p j are determined by the statistics of the underlying ARMA models at step k.At the SchWARMA level, this is fundamentally different from Lindblad style decay into an infinitely fast bath using Lindblad operators as described in App.A (i.e., B k = 0).This is evident, for example, in the fact that SchWARMA errors corresponding to stochastic unitary evolution offer the potential of control (for timecorrelated errors), whereas SchWARMA errors motivated by Lindblad evolutions will be purely dissapative regardless of the spectrum of the driving SchWARMA model.Furthermore, the correspondence between SchWARMA tangent vectors and the Lindblad master equation tells us that these Lindbladianstyle depolarizing channels can be simulated in SchWARMA by U = I 2 and with corresponding (possibly complex-valued) ARMA outputs y k for j = x, y, z.This model results in SchWARMA steps where (c.f., Eq.A4) and the resulting error can be verified numerically to be four Kraus operators corresponding to a depolarizing channel.Clearly, the x and y terms can be dropped off and the z element "compacted" to a 4 × 2 matrix to restrict the errors to dephasing.

IV. EXAMPLES
We now demonstrate the utility of SchWARMA for analyzing correlated noise within a variety of domains relevant to quantum information processing.We first build from singlequbit examples in quantum noise spectroscopy and control, showing strong agreement with established theory.Then we show the utility of SchWARMA in a multi-qubit circuit simulation, as it results in substantial computational savings over a stochastic Liouville approach.Code to reproduce these examples and others can be found in our mezze package [49].
A. Example 1: Quantum Noise Spectroscopy In this first example, we verify the SchWARMA noise spectrum is consistent with quantum noise spectroscopy.We model the dephasing dynamics of a qubit with the stochastic Hamiltonian H(t) = η(t)σ z + Ω(t)σ x using SchWARMA, where η is a stationary and Gaussian semi-classical noise and Ω infinitely fast control (as in a quantum circuit model).Stochastic dephasing noise leads to decay in the survival probability p of a qubit both prepared and measured in the |+ state (assuming no state preparation and measurement error), given by p = 1/2[1 + exp(− dω|F (ω)| 2 S y (ω))] where F (ω) is the filter function (Fourier transform of the control's modulation function, for single axis X control) [50,51].Under these assumptions, the noise spectrum can be inferred by applying different control inputs.Here, we use different gate sequences to define a regression problem in the fashion of [52][53][54].
To perform effective quantum noise spectroscopy, it is necessary to have spectrally diverse (in terms of the filter function) control sequences.To generate these, we used W/2 different sequences of common length W gates (128 for the bandlimited and multipole noise, 2048 for the f −α noise) whose discrete-time modulation functions [50,51] were defined by f k (m) = sign(cos(π(k − 1)m/W )) for k, m = 1, . . ., W/2, meaning that an X gate is applied when the sign of the modulation function changed, or an I gate otherwise.Assuming ideal state preparation and measurement in |+ , the Fourier transform of f k (m), F k (ω), can be used to define a system of equations to estimate the power spectrum S Y (ω) by inverting the system of equations S Y (2kπ/W ) = m |F m (2kπ/W )| 2 χ m , where χ m −log(2p m −1), with each control sequence repeated sufficiently many times to achieve an accurate estimate of the survival probability p m .Note that nonnegative least squares should be used to ensure nonnegativity of the estimated power spectrum.Additionally, in order to reduce sampling error, we did not perform repeated Bernoulli trials to construct pm , instead we took the average of the exact survival probabilities over 1000 independent SchWARMA trajectories.
Exploiting the universality of ARMA to produce any power spectra, SchWARMA simulations of quantum noise spectroscopy for several spectra are shown in Fig. 2. The noise generation was handled using Z-dephasing SchWARMA models.The bandlimited signals were MA models defined by a 128-tap Hamming window of desired bandwidth and cosine modulated to the appropriate center frequency (as needed).The multipole noise was generated by a purely AR SchWARMA model with coefficients set by placing three complex poles of the transfer function on the unit circle at locations corresponding to the desired peaks in the normalized frequency space and their complex conjugates.SchWARMA coefficients for the f −α noise were generated using the methods in [55].The resulting estimates strongly agree with the true spectra, with the primary limiting factors being Monte Carlo accuracy and conditioning of the regression equations.Next, consider a single qubit in the presence of multi-axis noise and control, (21) where η i (t) are zero-mean wide-sense stationary Gaussian processes with correlation functions η i (τ )η j (t) = f i (|τ − t|)δ ij , φ is a user-controlled parameter that specifies the axis of rotation, and Ω(t) is the time-dependent amplitude of the control pulse.Furthermore, let the η i have correlation functions where τ c,i is the correlation-time of the noise for axis i and γ i is the corresponding noise-amplitude.Decoupling control protocols can be an effective mechanism for reducing decoherence.The key to their error suppression is the fact that no error process happens instantaneously, despite this being a common approximation in open quantum systems theory.If one can control the system on a time-scale fast relative to the error dynamics, then error suppression is possible with these protocols [56][57][58][59].Furthermore, some protocols are more effective at suppressing different classes of noise compared to others.In particular, the XY 4 protocol [60] suppresses multi-axis semi-classical errors to second order, whereas repeated XX gates only suppress Y and Z errors to second order.Fig. 3 (left) shows SchWARMA reproducing these phenomena in terms of process fidelity, using the Hamiltonian in Eq.21 with separate ARMA models driving the η i , and instantaneous π control pulses.The process fidelity F p (U, Φ) at each time-step between the ideal unitary operator U , and the average noisy map, is defined as where * and denote complex conjugation and transpose, respectively.
For amplitude damping, decoupling protocols will not reduce decoherence here as we assume an infinitely fast bath, but consider the non-unital action of the map.All CPTP maps have an affine form [61], visualized in Bloch space as a combination of a unital contraction and rotation operations along with a non-unital shift from the origin.This shift is defined by the vector β, For amplitude damping β shifts the origin towards |0 .However, decoupling sequences will flip the direction of this shift, resulting in slower rate of unitality (1-||β||) decay, as depicted in Fig. 3 (right).This demonstrates SchWARMA simulating the impact of temporally-correlated classical noise processes that affect the rate of decoherence even for non-unital channels such as those observed in [48].

C. Example 3: Quantum Circuit Simulation
Next, to demonstrate the utility of using SchWARMA to simulate noisy multi-qubit quantum circuits, we model Z and X check circuits for the surface code [62].We compare an approach using both a full density-matrix simulation of the noisy pulse sequence with a first-order Trotter-Suzuki technique [43][44][45][46], and a SchWARMA simulation.The SchWARMA simulation assumes the gates to be instantaneous and follows each gate with a random unitary rotation from the SchWARMA model: for each qubit j at time-step k.
For the full Trotter-based simulation, each single-qubit gate is generated by a Hamiltonian: (23) where η i (t) are zero-mean wide-sense stationary Gaussian processes with correlation functions η i (τ )η j (t) = f i (|τ − t|)δ ij , φ is a user-controlled parameter that specifies the axis of rotation, Ω(t) is the time-dependent amplitude of the control pulse, and η i is defined in Eq. (22).Entanglement in the system is generated via a two qubit ZZ entangling operation The X stabilizer is given on the left and the Z stabilizer is given on the right.Measurements of the ancilla are omitted from these circuits.
FIG. 5.The CNOT decomposed into elementary gates.The Z rotations are taken to be virtual within all of our simulations and take no time and occur without error.
where Ω (k,j) ZZ (t) is a controllable coupling term and k and j denote the qubit index.To create equivalent noise models between the Trotter simulation and SchWARMA simulation, we used the methods outlined in App.B.
The standard X and Z check circuits for the surface code are given in Fig. 4: We compile the CNOT into the circuit shown in Fig. 5.The gates are defined as: and the controlled ZZ gate between qubits i and j is For the Z check operator, all the single qubit X and Y rotations between the CNOT gates cancel each other.This allows for a further simplification of the compiled circuit of just single-qubit rotations on the ancilla to start, a sequence of ZZ 90 gates, followed by single qubit rotations on the ancilla.There are leftover virtual Z 1/2 gates to be done on the data qubits, but we leave these as they would be virtual and taken care of during the next round of error correction in any system.We chose this gate-set as an exemplar and with the knowledge that it is a commonly seen Hamiltonian in experimental systems.
For the full density matrix simulation, we simulate the entire circuit using a first-order Trotter-Suzuki [43][44][45][46] approximation to the dynamics.We generate rectangular pulses for each of the gates in our gate-set.We average over different noise realizations and compute for average process map for the entire circuit.An illustrative example is provided in Fig. 6 where we plot the single and two-qubit control pulses for qubit 6.These two plots provide an illustrative example of how the full simulation of the surface code circuits are done.In (a), we show the single-qubit control amplitudes, Ω(t), for the Hamiltonian in Eq. ( 23) and the two-qubit control amplitudes, A(t), for the Hamiltonian in Eq. ( 24) for the qubits labelled in the legend.We see singlequbit X and Y controls being applied at the beginning and end of the circuit.These correspond to the application of the Hadamard gate for the X stabilizer.In the middle, the entangling operations are performed.In (b), we plot a single realization of a noise trajectory, ηi(t), for the multi-axis noise term in Eq. ( 23) for the same qubit used in (a).We choose noise values γx = γy = γz = 10 −4 and the correlation time τc,x = τc,y = τc,z = 32.Those are in units of inverse gate length and gate length respectively.The noise is always on across the entire duration of the circuit.These values yield a highly time-correlated trajectory.
0 (the top qubit in the circuit diagram) and a single realization of a noise trajectory for the terms in Eqs. ( 23) and (24).In contrast, the SchWARMA simulations are much simpler.For these, we do standard circuit simulation where we simply apply the gates as perfect unitary operations and follow them with SchWARMA generated noise after each gate.We average over the same number of noise realizations as in the Trotter-Suzuki simulation method.This does not require generating the full pulse wave forms, as we had to for the Trotter method, yet provides the same circuit level noise model.
We compare the process infidelity 1 − F p of Monte Carlo averaged CPTP maps for a five qubit circuit and plot the absolute error (difference in infidelity) between the SchWARMA and Trotter simulation methods in Fig. 7 (1000 Monte Carlo trials per point) for different noise amplitudes and correlation times for the Z check circuit (The X check circuit results are nearly identical and not displayed).SchWARMA shows very good agreement with the full Trotter simulation.The error observed between the two methods is consistent Correlation Time (Gate Length) with Monte Carlo sampling error as the slope of the error fit is 10 −1.5 ≈ 1/1000.Error proportional to the infidelity squared due to the SchWARMA model is negligible compared to sampling error when the noise is low, see C. The SchWARMA simulations take 1/∆t less computational time, where ∆t is the Trotter time step, which in our case is a factor of 1000 computational speedup.This example highlights the utility of using SchWARMA to simulate correlated noise within quantum circuits, with a complexity no worse than a state-vector simulator such as the ones used in Refs.[63,64].

D. Example 4: Continuously Driven Systems
SchWARMA can also be utilized to simulate continuously driven system dynamics in the presence of correlated noise, such as faulty adiabatic dynamics.We demonstrate this by simulating where H c (t) generates the pure continuously driven (controlled) dynamics and H e (t) is the error-generating Hamiltonian.Here, H e (t) describes a temporally correlated noise process.Generically, we can represent the error Hamiltonian as H E (t) = k η k (t)S k with η i (t)η j (t ) denoting the twopoint correlation function for the classical noise process.Note that S k are system operators (not necessarily Pauli operators) designating the manner in which the noise affects the system.Simulations are performed by Trotterizing H(t) as where ∆t sets the time resolution of the simulation and N T = T /∆t denotes the number of time steps.In order to examine the effect of correlated noise on the system, one must average over many trajectories of the noise and therefore, perform this full dynamics simulation numerous times.SchWARMA simulations are performed by partitioning the controlled dynamics and the correlated noise evolution.The total time evolution is composed of products of ideal evolution followed by SchWARMA-generated noise that effectively captures the effect of H e (t).More concretely, the evolution is explicitly approximated by where N sch = T /∆t is the number of SchWARMA time steps, ∆t = κ∆t is the effective SchWARMA sampling time and Correlated noise is generated via U E (t j ) = exp(−i γ y (γ) j S γ ), where the variables y (γ) j are defined by the underlying SchWARMA model and thus, they encapsulate the temporal properties of the noise.Simulating the dynamics in this manner allows one to leverage the noiseless approximation of adiabatic dynamics as a base for the faulty simulations.The noisy dynamics are then captured by averaging a desired metric over realizations of Eq. ( 27), where the noise is applied at times j∆t that may be much larger than the Trotter time resolution ∆t.
We consider three different continuously driven systems to demonstrate the efficacy of SchWARMA: (a) a 2-level Landau-Zener system (LZS) [65,66], (b) a 3-level LZS [67], and (c) N -qubit adiabatic Grover's Search [68].SchWARMA is shown to be in agreement with full dynamics simulation of each system, while only requiring a fraction of the simulation cost.In the subsequent simulations, ∆t/∆t = 100; thus, while the ideal simulation cost remains the same, the cost of a noisy simulation is approximately reduced by a factor of 100.
We model the faulty LZ system as and H e (t) = η x (t)S x , where S x and S z represent spin operators for either the spin-1/2 or spin-1 system.We simulate the dynamics from t 0 to t 0 to effectively capture the LZ behavior from t → −∞ to t → ∞ and compare our results to the analytical expressions obtained in Ref. [69] for the transition probability from |0 to |1 in the presence of slow correlated noise.Fig. 8 conveys that the analytical expression for transition probability as a function of total time is in good agreement with SchWARMA-driven simulations for both systems (a) and (b).The total time, noise correlation time, and noise power are all normalized with respect to the drive strength α for distinct values of λ = ∆ 2 /2α.The total drive time and noise correlation time are given by T = T 0 / √ α and τ c = τ 0 / √ α, respectively.The normalized times T 0 = 600 and τ 0 = 100 are chosen for the simulations shown in Fig. 8.We express the noise power f (0) = 0.003α in terms of the driving rate.Parameter definitions are chosen in accordance with Ref. [69].
In addition to the LZS, we consider the adiabatic implementation of Grover's search algorithm: where α(t) describes the time-optimal Grover control [68], |+ denotes the equal superposition state, and |m represents the target, marked state.The noise is modeled as H E (t) = η(t) i σ z i ; thus, we consider collective dephasing with respect to the computational basis.Simulations shown in Fig. 8 include comparisons between fully Trotterized dynamics and SchWARMA simulations as a function of correlation time with respect to ∆ min , the minimum spectral gap between the ground state and first excited state.The fidelity metric is , where |Φ(t) denotes the instantaneous ground state of H ad (t) and ψ(T ) is the time-evolved state at total time T .The noise power and correlation times are also normalized with respect to ∆ min .The comparisons indicate that SchWARMA is in good agreement with analytical predictions and full dynamics simulations for a variety of N -qubit system sizes.

V. DISCUSSION AND FUTURE WORK
We presented an approach, motivated by classic ARMA modeling, that parameterizes and simulates time-correlated dynamics in quantum systems.To highlight the versatility of this technique, we demonstrated its applicability in simulating dephasing noise (single and multi-axis) and non-unital dissipative noise, applied to single qubits as well as multiqubit circuits.The single-qubit (and qutrit) simulations using SchwARMA produce strong agreement with established theory, and the multi-qubit simulations tracked accurately to stochastic Liouville simulations, at substantial computational savings.
We emphasise that SchWARMA is not limited to the examples considered above, or even to channels near the identity as we have done.In particular, semi-classical noisy Hamiltonian dynamics are converted to SchWARMA models by multiplying the Hamiltonian elements by i, c.f., the single-and multi-axis dephasing examples above.The multiplication by i converts the Hermitian Hamiltonians to the skew-symmetric constraint needed for that portion of the Stiefel manifolds tangent space.This applies to arbitrary N -level systems including multi-qubit entangling noise.To extend the approach to arbitrary stochastic noise on Kraus operators, the log map for the Stiefel manifold [70] can be used.Here, this corresponds to the standard matrix logarithm on the Kraus-stacked Stinespring form to compute its corresponding Lie algebra element (multiplied by i), where the first N columns correspond to the tangent space of the corresponding Stiefel manifold.
Beyond the application of SchWARMA to general CPTP maps and arbitrary stationary Gaussian noise, the underlying ARMA formalism itself can be further generalized to model additional phenomena.At its core, SchWARMA is itself a non-linear vector ARMA (VARMA) model [32,36], and generalizing the underlying VARMA model can be done in a similar fashion to standard ARMA models, which would then scale the the tangent space elements and be matrix exponentiated.There is a rich field of generalizations of ARMA models in the literature [32], which may be appropriate for different physical phenomenon than standard ARMA models.For example, ARIMA models nonstationary noise, SARMA and PARMA models are often used to model stochastic variations around some otherwise periodic behaviors, and GARCH models model volatility in the overall noise power [71].Additional nonlinearities could be used to model non-Gaussian semiclassical noise, such as that found in [72].We also expect that time-varying model coefficients and cross-correlated driving terms x ( ) k will find use in the modeling of control noise and spatial correlations, respectively.
The results presented here open up many avenues for future work.We have shown it to have utility in experimental noise injection [73] and analysis of quantum error correction [25].The ability to generate temporally-correlated errors at a low computational cost can be immediately applied to other simulations and experiments relevant to quantum error correction, such as for threshold analysis in the presence of absolutely continuous errors [64].Similarly, SchWARMA could be used to analyze and extend techniques for quantum system characterization [74][75][76][77] and quantum noise spectroscopy in the presence of correlated noise.Given the connections between ARMA and the fields of spectrum estimation, classical signal processing, and control, we believe there are a number of straightforward generalizations from the classical realm to that of open quantum systems.In particular, this approach could be a "discrete time" analogue to the filter function formalism [51,78], allowing for noise filtering and mitigation.
This indicates several relevant factors for SchWARMA.First, to identify tangent space elements that generate target non-unitary Kraus operators (with some parametric form), we need to first find the Kraus operator that can be set closest to I N , then set B k to the Kraus operators that are not that near-I N operator.Equivalently, if we have a known master equation in the above canonical form that we wish to emulate with SchWARMA, the Stiefel manifold tangent elements B k should be set to L α .Again, while the Lindblad and SchWARMA approaches are equivalent to first order, SchWARMA temporarily propagates the bath and so the trajectories will diverge.That said, it is clear that small SchWARMA steps will approximate the Lindblad master equation.Furthermore, as the SchWARMA approach as formulated always produces CPTP maps (that may arise from non-CP divisible trajectories), a CP-divisible trajectory could be derived from the Lindblad master equation that produces equivalent maps at the gate time-scale.
This suggests an alternative (and essentially equivalent) usage of SchWARMA for stochastic Lindblad master equations where individual time-correlated steps are defined by |ρ k+1 = L k |ρ k , with where H k is a "unitary" SchWARMA Hamiltonian, and y are the outputs of (possibly complex) ARMA models as in the original SchWARMA formulation.This formulation is clearly dissipative and CP-divisible, as the coefficients Lindblad coefficient ||y k || 2 are nonnegative.It is clear, however, that maintaining the SchWARMA "bath" offers a potential path forward for non-CP divisible dynamics.
Consider the case of the noisy Hamiltonian H(t) = η(t)σ, for some bounded Hermitian operator σ.Let η(t) = µ(t) + η(t) with µ(t) a deterministic integrable function and η(t) a stochastic process satisfying E[η(t)] = 0 and E[η(t)η(s)] ≤ C, i.e., uniformly bounded for all s, t ∈ R. Examples of such processes include single axis dephasing evolutions in the absence of control, or single axis multiplicative control noise.Then, for any time T > 0 and a specific stochastic trajectory of η, U (T ) = exp(−i T 0 dt η(t)σ).Thus, the expected superoperator L(T ) = E [U (T ) * ⊗ U (T )].Since σ is Hermitian, it can be diagonalized σ = ΣΛΣ † , so that U (T ) = exp(−i (η(t)σ)) = Σ exp(−i η(t)Λ)Σ † , and furthermore, the Liouvillian can be diagonalized by and so the dynamics of L are ultimately determined by exp(i η(t)(λ j − λ k )) where λ j , λ k are eigenvalues of σ.Thus, by averaging over random trajectories, we have that L(T ) is determined by the expectations E[exp(i T 0 dt η(t)(λ j − λ k )].Now since η(t) is a stochastic process, T 0 dt η(t) S(T ) is a random variable and the E[exp(iS(T )(λ j − λ k ))] is its characteristic function evaluated at (λ j − λ k ), which exists.Now, if we assume that η(t) is a member of a closed (under addition) functional family, such as a Levy α-stable distribution (an example of such a distribution is the Gaussian distribution), then S(T ) are all the same form of random variable (e.g., Cauchy, Gaussian), and we can exploit known results of characteristic functions.
In the case that η (and thus S) are Gaussian, we have that By construction, we have that E[S(T )] = T 0 dtµ(t) which exists by assumption, and Var(S(T is the continuous time Fourier transform of the autocorrelation of η, i.e., the power spectrum of η. The above expression is useful for determining SchWARMA models for dephasing noise, as it indicates that if we are interested in discrete times kT (e.g., as in a circuit with fixed gate time), then we need to find an ARMA model with autocovariance r s [k] that satisfies the continuous time version of Eq.B2-B4, i.e., replacing the summations on the right hand side of Eq.B2-B4 with integrals as in Eq.C4.Again, we reiterate that since the possible power spectra of AR, MA, and ARMA models are each dense in the space of valid power spectra (and thus valid autocovariances), finding an AR, MA, or ARMA model that is arbitrarily close to given autocovariance is always possible in principle, and many practical methods for finding solutions exist such as [85].Again, we emphasize that this is a matching on the discrete kT timescale of interest, and that the continuous time power spectrum of a derived ARMA model will always be periodic in 1/(2T ), but result in equivalent sampled dynamics at time steps kT , meaning this difference is irrelevant on the gate time scale.
The above discussion highlights a subtle distinction about conversion between continuous and discrete time-scales, or between different discrete-time scales.Clearly, a given functional form for an autocovariance or power spectrum, such as squared exponential, is not preserved by the conversion due to the time-integration.This agrees with the intuition that a noise source with correlation much shorter than a gate time should be well modeled by white noise at the gate time scale.That is not to say, however, that the short-time correlations are completely ignored by a SchWARMA model, they are instead aggregated at the coarser time scale (and to the accuracy of the ARMA fit).This also means that parameter sweeps at the gate time scale do not map perfectly to parameter sweeps at the continuous time scale, but coarser features such as correlation time, will remain accurate to the order of the discretization.

Average error operator after a single gate
In this section, we consider a single gate, U c that is applied over a time T in the presence of a native noise process.Using the Magnus expansion, we derive the expression for an average error operator applied after the perfect gate to model the dynamics.This average error operator can be then reduced to a SchWARMA model to characterize the error induced by the noise at this gate time-scale.
Consider the gate U c (0, T ) = T exp T 0 dtH c (t) being applied by the Hamiltonian H c for time T , with an unwanted noise process given by the Hamiltonian H η (t).We can define the control and noise Hamiltonians by a given set of basis operators {P p }, where, Ũ (t) = U † c (0, t) Ũ (0, t)U c (0, t).The operator, Ũ (0, T ) characterizes the error in the application of this gate.
In the following, we use the Magnus expansion to obtain an approximate expression for the error applied.
Going to an interaction frame, using idU/dt = (H c + H η )U , and −idU † c /dt = U † c H c , we have the following expression for Ũ , In the following analysis, we truncate the Magnus expansion, which is a reasonable assumption to make in the limit of weak noise, T max t { H η (t) } 1. Truncating at the first term in the Magnus expansion, we have the following expression for the error operator, Ũ (0, T ) ≈ Ũ1 (0, T ) dt s q (t) P q (C18) = exp −i q S q (T )P q , where, we introduce the notation, S q (T ) = T 0 dt s q (t) = T 0 dt η p (t)c pq (t).Note that in principle, the above expression can be generalized to take into account higher order terms in the Magnus expansion, but we leave that for future work.Utilizing (C19) in the above expression, one can obtain expressions analogous to (C2) to derive the effective noise model for this quantum operation.
The error in performing the truncation to first order in the Magnus expansion can be estimated from the norm of the second term in the expansion (this is appropriate again in the limit of weak noise), {s q1 (t 1 )s q2 (t 2 )} [P q1 , P q2 ] Furthermore, the expression in (C20) suggests that this approximation is accurate to the cross-covariance of the noise sources.For noise sources that act independently, E [s q1 (t 1 )s q2 (t 2 )] ∼ δ q1q2 one needs to consider higher order contributions and thus the error scales atleast as the covariance squared.Since the infidelity of the map is proportional to the variance in this case, the SchWARMA model error here will be proportional to the infidelity squared.This is consistent with other results derived from filter functions, e.g., [87].

Generalization to discrete steps
Let us consider the case where we now have an evolution from t = 0 to t = N g T , and we are interested in the discrete time-steps kT .In this case, one can generalize the derivation, by examining the evolution from time kT to (k + 1)T , U (kT, (k + 1)T ) = T exp
FIG.6.These two plots provide an illustrative example of how the full simulation of the surface code circuits are done.In (a), we show the single-qubit control amplitudes, Ω(t), for the Hamiltonian in Eq. (23) and the two-qubit control amplitudes, A(t), for the Hamiltonian in Eq. (24) for the qubits labelled in the legend.We see singlequbit X and Y controls being applied at the beginning and end of the circuit.These correspond to the application of the Hadamard gate for the X stabilizer.In the middle, the entangling operations are performed.In (b), we plot a single realization of a noise trajectory, ηi(t), for the multi-axis noise term in Eq. (23) for the same qubit used in (a).We choose noise values γx = γy = γz = 10 −4 and the correlation time τc,x = τc,y = τc,z = 32.Those are in units of inverse gate length and gate length respectively.The noise is always on across the entire duration of the circuit.These values yield a highly time-correlated trajectory.

FIG. 7 .
FIG. 7. Process infidelities for surface code Z check circuits (a) and the absolute error between the SchWARMA and Trotter simulations (b) as a function of noise parameters τc,i ((a) x-axis), γi (color) defined below Eq. (21), and process infidelity ((b) x-axis).The noise strengths are γx,y,z ∈ {10 −12 , 10 −10 , • • • , 10 −2 } (bottom to top on (a) plot) in units of inverse gate length.In (a), the Trotter simulations are the line plot while the SchWARMA simulations are the marks with Monte Carlo (1000 samples) error bars.In (b), the absolute difference between the SchWARMA and Trotter simulations are plotted with a fit line that is consistent with Monte Carlo sampling error.

5 FIG. 8 .
FIG.8.SchWARMA (symbols) vs. analytical results and full dynamics simulation of continuously driven systems (both shown as solid lines).The comparison indicates good agreement between SchWARMA and alternative approaches.LZ results are shown for an (a) 2-level and (b) 3-level system with noise power f (0) = 0.003α for different driving rates α.Grover's search algorithm comparisons [panel (c)] focus on varying correlation time τc using a fixed noise power of f (0) = 0.001∆ 2 min .