Simple master equations for describing driven systems subject to classical non-Markovian noise

Driven quantum systems subject to non-Markovian noise are typically difficult to model even if the noise is classical. We present a systematic method based on generalized cumulant expansions for deriving a time-local master equation for such systems. This master equation has an intuitive form that directly parallels a standard Lindblad equation, but contains several surprising features: the combination of driving and non-Markovianity results in effective time-dependent dephasing rates that can be negative, and the noise can generate Hamiltonian renormalizations even though it is classical. We analyze in detail the highly relevant case of a Rabi-driven qubit subject to various kinds of non-Markovian noise including $1/f$ fluctuations, finding an excellent agreement between our master equation and numerically-exact simulations over relevant timescales. The approach outlined here is more accurate than commonly employed phenomenological master equations which ignore the interplay between driving and noise.


Introduction
As quantum technology advances, it is becoming increasingly more important to have tools for modelling dissipation and noise in driven quantum systems that are accurate, physically transparent, and easily implemented numerically. The standard approach is to describe dissipation using a time-local Lindblad master equation (LME) [1,2]; for simple systems involving a few qubits, such equations can be numerically simulated using off-the-shelf packages (see e.g., [3]). Unfortunately, Lindblad master equations are inaccurate descriptions in many relevant situations, even if the system is just a single qubit. If the noise is non-Markovian, and if the qubit is driven, then using a Lindblad ap- proach generally requires an uncontrolled approximation. While more complicated formalisms exist for non-Markovian, driven problems (see e.g., [4]), they typically lack the transparency and simplicity of an LME. One approach is to attempt to model the effects of highly-correlated noise using an ad-hoc Lindblad master equation with either constant rates, or with rates that are explicitly time-dependent (even when the noise is stationary). These approaches are common in superconducting quantum information (see e.g. [5,6,7,8]), where the dominant dephasing noise typically has a 1/f character [9,10,11,12]. While these equations are often physically well-motivated and retain the simplicity of an LME, they rely on an uncontrolled approximation, so their ultimate validity is unclear. Other alternative attempts of capturing the effects of non-Markovian noise been also proposed (e.g., [13]), however, they have not yet been widely adopted.
In this paper, we discuss a systematic, controlled approach for describing the dynamics of driven quantum systems subject to non-Markovian noise with a time-dependent but time-local master equation that is reminiscent of a LME. We focus on the general case where the dissipation is due to classical non-Markovian noise, but is still non-trivial, as the system Hamiltonian (which includes driving) does not commute with the noise Hamiltonian. We use the machinery of cumulant expansions (well-known in other contexts, see e.g., [14,15,16,17]) to derive Lindbladstyle master equations in this general setting. Our approach transparently shows how the combination of non-Markovianity and driving leads to effects that are missed in more naive approaches: the noise gives rise to a renormalization of the system Hamiltonian (even though there is no Lamb shift physics as the noise is classical), and is also described by two time-dependent rates, one of which is generically negative. As we show, this parameterization of the master equation in terms of generalized time-dependent rates provides a wealth of physical intuition. We call this resulting equation a pseudo-Lindblad master equation (PLME) 1 .
To illustrate and test our approach, we present a detailed study of a Rabi-driven qubit subject to various kinds of classical non-Markovian noise. We find that even in cases where the noise has no well-defined correlation time (i.e., 1/f or quasistatic noise), the PLME provides an accurate description of the full dynamics over times that can be much longer than the Rabi period. We also show that the negative rates we obtain in the PLME are in excellent correspondence to the negative rates that exactly characterize the instantaneous generator (as obtained from numerical simulation).
While cumulant expansion techniques have been employed to model classical noise in specific chemical systems [19,20] and in specific qubit systems [21], our work differs in several respects. In particular, the use of these techniques to give a microscopic derivation of a pseudo-Lindblad structure (with the emergence of negative dissipation rates due to colored noise and driving) is new to our work. We note that Ref. [22] showed how a very general class of quantum maps can be expressed in pseudo-Lindblad form; however, unlike our work, the connection to explicit microscopic models of noise and driving was not discussed.
The remainder of this paper is organized as follows. 1 A simialr approach has been recently also presented in [18].
In Sec. 2, we introduce the generic classically-stochastic quantum system of interest. In Sec. 3, we derive the pseudo-Lindblad master equation. We specialize to the specific example of a Rabi-driven qubit subject to non-Markovian dephasing noise in Sec. 4. We then explore different kinds of non-Markovian noise in subsequent sections: quasistatic noise in Sec. 5, noise with a Lorentzian spectrum in Sec. 6, and 1/f noise in Sec. 7. We conclude in Sec. 8.
2 General setup: Driven quantum system subject to classical non-Markovian noise Consider a generic setup where a driven finitedimensional quantum system (with a time-dependent system Hamiltonian) is coupled to classical noise η(t) via a Hermitian system operator g(t)Â. Here, g(t) represents a possibly time-dependent coupling constant. Working in the interaction picture with respect to the noise-free HamiltonianĤ 0 (t), the dynamics is described by the stochastic Hamiltonian . Here, we allow for a general time-dependent Hamiltonian H 0 (t) that need not commute with itself at different times. Therefore, whereρ η (t) is the system density matrix for a particular realization of the noise η(t).
In what follows, we will focus on the case where η(t) is zero-mean Gaussian noise. It is fully characterized by its autocorrelation function Our goal is to obtain a closed, time-local evolution equation forρ(t), the system density matrix averaged over all realizations of the noise. This is easily done in either the case where (a) the noise is white-noise (see e.g., [4]), or (b) [Â(t),Â(t )] commutes with itself at all times (see [23,24] and discussion below). These correspond respectively to the Markovian noise limit, and to the limit of undriven, pure dephasing dynamics respectively. In all other cases, deriving an accurate time-local equation of motion for our system is a non-trivial task (despite only having classical Gaussian noise).

Cumulant expansion and the pseudo-Lindblad master equation
To systematically derive an approximate time-local equation for the noise-averaged density matrixρ(t), we employ the generalized cumulant expansion first introduced by Kubo [14], and further developed by many others since (e.g., [15,16,17]). It is worth stressing that such technique is closely related to work on timeconvolutionless (TCL) master equations; in particular one can think of the cumulant-expansion method as a means to systematically obtain an n-th order contribution to the TCL equation generator [4,25,26]. Thus, our approach below, can be viewed as yielding the classical-noise limit of a TCL equation (truncated at a particular order in noise strength, and rewritten in a pseudo-Lindblad form). The technique can be viewed as yielding the classical-noise limit of the timeconvolutionless master equation [4]. Hence, while the basic formalism presented here is discussed in many places, our main goal is to make the general structure clear, in particular the interpretation of the final equation as a pseudo-Lindblad master equation with a modified system Hamiltonian and two dephasing rates (see Fig. 1). One of these rates is negative, which is a signature of the non-Markovianity of the noise process [22,25,4]. We point out that one common physical interpretation for these negative rates is related to the fact that (even a classical, as is the case here) non-Markovian environment can lead to a revival of coherence: as the system evolves, the information loss to the bath can be at times (often partially) undone and "flow back" into the system. In setups where the bath is quantum, this can be often heuristically understood by considering a simple model where the bath is thought of as an auxiliary mode that couples to its own (often Markovian) environment, and which also coherently couples to the original system in question. The revival of coherence can then be seen as resulting from the (partially) coherent system-auxiliary mode interaction, i.e., some of the information has a chance to traverse back the original system, before it gets permanently lost to the auxiliary mode's environment.
As is standard, we also assume that at t = 0, the initial system density matrixρ(0) is uncorrelated with the noise η(t) (i.e., it is not stochastic); this singles out t = 0, effectively breaking time-translation invariance.

2 nd -order pseudo-Lindblad master equation
We start by truncating the generalized cumulant expansion at 2 nd -order in the noise. We then obtain a timelocal evolution equation forρ(t) (the noise-averaged density matrix) that has a familiar double-commutator structure. We can write this equation in a suggestive manner 2 : whereÂ Both of these quantities have a simple physical interpretation: •Â avg (t) corresponds to the weighted average ofÂ(t) at times t < t. The weighting here is controlled by the noise autocorrelation function.
In that case,Â avg =Â, and Eq. (4) is a generalized dephasing master equation that causes coherences between eigenstates ofÂ to decay. For this special case, the cumulant expansion truncates at 2-order and the equation is exact.
To obtain further intuition, we will now rewrite Eq. (4) in a form that resembles a standard Lindblad master equation, i.e., express it in terms of an effective Hamiltonian and a set of jump operators. Given that Eq. (4) is trace and Hermiticity preserving, this is always possible [1], although the rates associated with dissipative processes will in general be time-dependent and possibly negative [4,22]. Our goal is to show how the parameters in this pseudo-Lindblad form reflect the properties of our system and provide physical intuition. The first step is to quantify both the size of the op-eratorÂ avg (t), and the degree to which it is orthogonal toÂ(t) (the instantaneous coupling operator). We achieve both using the Hilbert-Schmidt inner product, i.e., X ,Ŷ = Tr X † Y and its corresponding norm ||X|| 2 HS = X ,X . We first introduce the scalar quantities as well as the normalized operatorŝ Defining the angle φ(t) via we can now expressÂ avg in terms of operators parallel and perpendicular to the instantaneous system operator A(t): Thus even though our original stochastic Hamiltonian in Eq. (1) involved only a single system operatorÂ(t) at each instant, our dynamical equation involves two independent operators:Â(t) and the orthogonal operator b ⊥ (t).
Going forward, we will see that there are two timedependent parameters controlling our master equation: the angle φ(t) in Eq. (12) (which parameterizes how different A(t) was at earlier times compared to the present time), and an effective rateΓ(t): This is just a normalized version of the instantaneous dephasing rate Λ(t) introduced earlier.
We are now in a position to express the master equation in Eq. (4) in a form that is formally analogous to a standard Lindblad equation. We let D[X]ρ = XρX † − 1 2 {X †X ,ρ} denote the standard Lindblad dissipative superoperator. A few lines of algebra transform Eq. (4) into: Equation (15) is our general pseudo-Lindblad master equation (PLME), and is a central result of our work. The first term describes coherent dynamics generated byĤ ren (t), a noise-induced correction to the system Hamiltonian. It is given by the Hermitian operator: The second two terms of Eq. (15) describe generalized time-dependent dephasing processes, with Hermitian time-dependent operatorsê ± (t) that are rotated versions ofb(t),b ⊥ (t): The corresponding rates are given by Several comments are now in order: The Hamilto-nianĤ ren (t) is not a standard Lamb-shift renormalization [4], as there is no quantum bath here (only classical noise). Heuristically, 2 nd -order perturbation theory in the stochastic HamiltonianĤ(t) could yield new processes quadratic in η(t) that have a non-zero mean. Similar terms have been discussed in specific contexts [27,19]. Our derivation provides a more general description of such Hamiltonian renormalization terms, and shows that they are present whenever φ(t) = 0.
Furthermore, for any non-zero φ(t), we always have two effective dephasing processes involving two orthogonal operatorsê ± (t). This reflects the fact that dissipation involves both the system operatorÂ(t) at the current time, as well as at earlier times. For small φ(t), these operators are almost the same asb(t) andb ⊥ (t).
Finally, for any non-zero φ(t) the effective instantaneous dephasing rate Γ − (t) is always negative. While surprising, such negative dephasing rates are known to arise in problems where dephasing dynamics is nonmonotonic in time [28]. As we will show via comparison to full numerical simulations, the negative dephasing rate we predict is not an artefact of our approximation, but is instead in excellent agreement with features of the exact dissipative evolution.

Simple limits
It is reassuring to confirm that Eq. (15) reduces to known results in some simple limits. First, consider the limiting case of Markovian noise, where the noise autocorrelation function approaches a δ-function: S(t) → σδ(t). In this case, the bath has no memory, and there-foreÂ avg (t) =Â(t). As a result,Â avg (t) has no component orthogonal toÂ(t), hence Z ⊥ (t) and φ(t) vanish. Our pseudo-Lindblad master equation thus reduces to: This is a Lindblad master equation describing generalized dephasing with an instantaneous rate σg(t) 2 . In this white-noise limit, this equation is exact, as there are no higher-order cumulants (as the lack of any memory means the dynamics is insensitive to whetherÂ(t) commutes with itself at different times).
A second simple case is that of a time-independent A(t). Again, in this case φ(t) = 0, and we recover the expected dephasing Lindblad equation: which is also exact in this case, as all higher cumulants vanish.
4 Example of the formalism: Rabi driven qubit subject to non-Markovian noise 4.1 Form of the PLME Next, we apply the general results derived above to an ubiquitous but surprisingly rich problem: a Rabi-driven qubit that is subject to classical non-Markovian dephasing noise. Treating the driving in the rotating-wave approximation, the rotating frame Hamiltonian of the system iŝ where Ω(t) is the Rabi drive amplitude. This system the form of our general problem, with g(t) = 1.
The corresponding interaction-frame system operator with Note that our treatment below goes beyond previous work on this problem that focused exclusively on the slow-noise limit [27,29].
We can now directly apply the results of Sec. 3. Using theÂ(t) given by Eq. (22), the PLME (based on keeping contributions up to 2 nd -order in the noise strength) follows from Eq. (15). It takes a transparent form in the rotating frame used to write H 0 (t) in Eq. (21): The instantaneous dephasing Γ ± (t) are given by Eq. (18), while from Eq. (16), the noise-induced Hamiltonian correction iŝ Here φ(t) is given by Eq. (12). The corresponding jump operatorsτ u/v (t) are rotated Pauli matrices, and read Note thatτ u/v (t) are the same as √ 2ê ± (t) in the general PLME Eq. (15), but expressed in the rotating frame used to write Eq. (21).
In the next sections we will explore the accuracy of the above approximate master equation, by considering its behavior given specific non-Markovian noise correlation functions.

Corrections to the PLME from 4 th -order cumulants
We can further improve the accuracy of our PLME by including higher-order noise terms that arise at 4 thorder in the cumulant expansion. The 2 nd -order equation, Eq. (24), has two dephasing dissipators involving Pauli matrices in the z − y plane. As one might expect, 4 th -order contributions modify the rates and forms of these dissipators. There is, however, a more interesting effect that emerges at 4 th -order: a new, third dissipator that corresponds to dephasing along the x axis (again, working in the rotating frame used to write Eq. (21)). This process corresponds to the term: In App. A we give the general expression for the 4 thorder cumulant contribution to the equation of motion, which in turn, lets us calculate the explicit form of Γ x (t). We stress that this new term has a simple heuristic explanation. We saw that at 2 nd order noise generates a renormalization of the coherent system Hamil-tonianĤ ren (t) ∝σ x , i.e., a correction to the effective Rabi drive amplitude. Not surprisingly, there are fluctuations in this noise-induced change of the drive amplitude. The dissipator in Eq. (28) simply describes the dephasing associated with these fluctuations. noise, with an autocorrelation function that is independent of time: While highly non-Markovian (i.e., noise correlations have an infinite time correlation), this noise is especially simple to work with, and results in compact analytical expressions. It also makes benchmarking our approach against exact numerics easy (as the propagator for any single noise realization can be calculated exactly). Note that in the absence of a finite correlation time (i.e., exponential decay of S(t) at long times), it is well known that the generalized cumulant expansion does not converge if one considers extremely long evolution times [15,15]. This deficiency is however not an issue for most quantum information applications, as one wants to understand finite-time evolution (i.e., during the execution of a gate or computation), and not some long-time dissipative steady state. For such short-to-intermediate timescales, the cumulant expansion can still provide an accurate description of the dynamics, whenever noise is not especially strong (something we show explicitly in what follows by comparing against numerically-exact results).

Rates and Hamiltonian renormalization
For simplicity, we assume a constant Rabi drive amplitude, Ω(t) = Ω. We then find that the two instantaneous dephasing rates in the 2 nd -order PLME of Eq. (24) are given by Eq. (18): while the angle φ(t) that controls the form of the two dephasing dissipators is simply given by φ(t) = −Ωt/2. We see that both dephasing rates have a simple oscillatory form, with Γ + (t) ≥ 0 and Γ − (t) ≤ 0. Further, both rates vanish periodically at times t n = 2πn/Ω. This reflects the basic physics of continuous dynamical decoupling [30,31,32]: the Rabi drive tends to periodically average out the dominant dephasing effects of the noise. The remaining noise-induced term in the 2 nd -order PLME is the Hamiltonian renormalization Similar to the dephasing effects, we see that this noiseinduced Hamiltonian also vanishes periodically in time.
We can also compute the dominant symmetry breaking term that emerges from the 4 th -order cumulant corrections to the PLME, i.e., the induced x dephasing term in Eq. (28). We find that the rate of this process is always positive, but not periodic:

Limiting forms for weak driving and/or short times
It is instructive to examine the form of the PLME parameters in the short-time or weak-drive limit, Ωt 1. In this case, we find (dropping terms O (Ωt) 3 ) whileĤ ren (t), reduces tô To 0 th -order in the drive amplitude, both the negative dephasing rate and the Hamiltonian correction vanish, and we are left with a single, linear-in-time positive dephasing rate. This corresponds to the expected generalized dephasing rate for the undriven problem (a limit where, as discussed, our equation becomes exact). At 2 nd order in the drive, things become more interesting: the magnitude of Γ + (t) is suppressed, and we generate a non-zero negative dephasing rate Γ − (t).

Comparison of instantaneous dissipation rates: PLME vs numerically-exact results
While negative dephasing rates may be alarming at first, they are a common feature of non-Markovian dynamics [22]. To instill confidence, we numerically calculate the instantaneous rates that characterize the full dissipative evolution of our system (without any approximation). The definition of such rates is discussed in many places (e.g., Refs. [22,28]). We numerically calculate the full quantum map V(t) that propagatesρ(0) toρ(t), and then obtain the instantaneous generator of the dynamics via [28] As was done in Sec. 3, L inst (t) can be brought into a Lindblad-like form. This lets us extract the generalized instantaneous dephasing rates that characterize the full dynamical map. Dephasing rate (10  15), as well as the rate Γx(t) (green), calculated by keeping the 4 th -order noise contributions in the cumulant expansion. For comparison, the black curves are the rates obtained from the exact system evolution, showing good agreement with our analytically derived approximations. We especially stress that that the large-in-magnitude negative rate Γ−(t) is present in the true dynamics of the system. In the simulations, we take σ/Ω = 0.05. Figure 2 shows the results for our analytical predictions for Γ ± (t) (orange and blue curves) which closely reproduce the two large-in-magnitude rates obtained from the exact evolution (black curves). We stress the central point of these results: the highly negative Γ − (t) that is present in our PLME, is also there in the true evolution. As it is as large in amplitude as Γ + (t), its presence can be important in accurately modeling the qubit's evolution. Furthermore, is worth stressing that the exact dynamics also predicts a third rate (black dotted curve), which is accurately captured by Γ x (t), obtained from keeping the 4 th -order terms in the cumulant expansion.

Benchmarking the PLME against numerically-exact results
In order to asses the usefulness of the PLME, we need a means to quantify its accuracy against the exact system evolution. For our purposes, we will compare against the numerically-exact evolution, obtained by numerically averaging unitary evolution over realizations of the stochastic process η(t).
There are various metrics that can be used to assess the accuracy of an approximate open quantum systems evolution. One approach is to consider the performance of the approximation only for a small, specific subset of initial states [33]. A more complete approach is to try to compare the full structure of the exact and approximate evolution maps (e.g., Ref. [34]).
Here, we will focus on a measure that characterizes the entire evolution map, and is not tied to a particular choice of initial states (although in App. B, for completeness, we also show the evolution of specific Pauli operator expectation values given a particular initial state). For some evolution time t, let J exact (t) denote the Choi matrix of the exact evolution map, and J approx (t) the Choi matrix of the evolution map obtained via the approximation of interest (i.e., the PLME). We then define the approximation error at time t by: Here ... is the diamond norm [35,36]. We choose the diamond norm to measure the size of the deviation between the approximate and true maps, as this gives (t) a direct operational meaning: it is directly tied to the single-shot channel distinguishability between these two maps [36]. is 0 when the map associated with the approximate dynamics cannot be distinguished from the exact map after a single evolution (no matter how the initial state is optimized), and equals 2 when the two maps are maximally distinguishable [36]. At a heuristic level, it characterizes the worst case performance of our approximation, i.e., the largest error maximized over all possible initial states. With this definition in hand, we can now quantify the approximation error of the PLME for different choices of system parameters and noise spectra. To put these results in context, we will also consider the approximation error associated with a much simpler kind of time-local master equation. In this approach (which we will term the 0 th -order master equation), one obtains a Lindblad-form master equation in two steps. First, one ignores the Rabi driving, and derives a dephasing dissipator (with an Ω-independent rate). Then, one simply puts back the Rabi driving in the Hamiltonian term of the master equation.
Working in the interaction frame (and focusing on quasistatic noise), the 0 th -order master equation takes the form: We stress that this kind of approach is commonly used because of its simplicity [5,6,7,8]. Unlike the PLME, it results in only a single, Ω-independent dephasing rate. The approximation error (t) of the PLME as well as the naive 0 th -order master equation of Eq. (38) are shown in Fig. 3, for quasistatic noise that is weak compared to the Rabi frequency, σ/Ω = 0.05. For comparison, besides only looking at the 2 th -order PLME explicitly shown in Eq. (24), here, for quasistatic noise, we also include the performance of a PLME obtained by keeping all the 4 th order contributions, which result in: corrections to the rates Γ ± (t), the corresponding jump operators (that rotate in the z − y plane), the Hamiltonian renormalizationĤ ren (t), and importantly, the x-dephasing rate Γ x (t), and its jump operator.
We see that both forms of the PLME (obtained through either 2 nd -or 4 th -order cumulant expansions) have a significantly smaller approximation error than the 0 th -order equation. As expected, (t) for both schemes grows with increasing time. Nonetheless, for weak noise, the PLME approximation error remains low even if t 1/Ω. To be concrete, we see that for the chosen parameters, the 0 th -order approach has an approximation error (t) < 10 −3 only if tΩ 1. In contrast, the 2 th -order PLME has (t) < 10 −3 at times where tΩ 5, while the 4 th -order PLME, over even longer longer times, where tΩ 17. We also point out that in App. B, we present an example of the expectation values evolution for a particular initial state.
For further insights, we can also consider the behavior of (t) for the different approaches in the t → 0 limit. Unfortunately, the asymptotics of the diamond norm are difficult to compute. As a proxy, we can instead consider a different error measure where the diamond norm is approximated by the 1-norm 3 . We find from numerics that this quantity behaves analogously to (t) in the t → 0 limit. Further, the 1-norm metric is amenable to analytic estimates. Using this equivalence,  : Approximation errors (t) of the 2 nd -order PLME (blue) and the 0 th -order approximation (orange) for a Rabi driven qubit subject to quasistatic noise, in the short time limit Ωt 1. Black dotted and dashed curves show (t) associated with the 0 th -order approximation scales as t 3 . In contrast, for the PLME is scales as t 5 . One also finds a different asymptotic scaling with the noise strength σ (σ 2 for the 0 th -order equation, σ 4 for the PLME. These results demonstrate the gain in accuracy associated with using the PLME). All curves correspond to σ/Ω = 0.2.
we find that for short times, the approximation error for the 2 nd -order PLME scales as: whereas for the naive 0 th -order equation we obtain From Fig. 4 we see that Eqs. (39) and (40) accurately described the asymptotic behavior of (t) in the small t limit. The more favorable approximation error scaling of the PLME is not surprising, as in contrast to the 0 th -order approximation of Eq. (38), it fully captures 2 nd -order noise contributions to the qubit's dynamics. Finally, while the PLME yields a highly accurate approximation, its form is not guaranteed to generate a completely positive map. In particular, for certain initial conditions, the evolved density matrix may acquire small, negative eigenvalues. The size of this unphysical negativity is necessarily small for short times and weak noise, as it is bounded by the error metrics we have considered above. We discuss the positivity of our approximation in more detail in App. C.

Rabi driven qubit subject to noise with a finite correlation time
In this section we study the Rabi-driven qubit (see Eq. (21)), now subject to classical dephasing noise  that has a finite correlation time τ c . We take η(t) to be an Ornstein-Uhlenbeck process, implying: This yields a Lorentzian noise spectral density as encountered generically in various settings (i.e., classical telegraph noise produced by a two-state fluctuator, in the Gaussian approximation). It also provides a convenient test case, as we can interpolate between the quasistatic noise of the previous section, and noise that is Markovian, simply by tuning τ c .

PLME parameters
The 2 th -order PLME of Eq. (24) has the effective dephasing rates given by: The jump operators corresponding to these rateŝ τ u (t),τ v (t) are determined by Eqs. (26) and (27) Note that φ(t) exhibits initial oscillatory behavior, but for times longer than the correlation time (i.e., t τ c ), it has the simple form φ(t) ≈ −Ωτ c . Hence, for noise that is only very weakly correlated, we see that θ(t) ≈ Ω(t − τ c /2); the jump operator associated with the rate Γ + (t) is equivalent to the system operatorÂ(t), but acting at a slightly "earlier time".
We also obtain a Hamiltonian renormalization In the τ c → ∞ limit, one recovers the results for quasistatic noise presented in Sec. 5. In contrast, one can obtain the Markovian white-noise limit by taking σ 2 = D/τ c and letting τ c → 0. In this limit Γ − (t) and H ren (t) vanish as expected. Equation (42) implies that a 2 nd -order PLME, for any non-zero correlation time, we still have two instantaneous dephasing rates Γ ± (t), with Γ − (t) < 0. Similarly to the case of quasistatic noise, at 4 th -order, one also observes a third rate Γ x (t), that corresponds to dephasing along the qubit's x-axis. We explicitly show the analytical form of Γ x (t) in App. D. Figure 5(a) compares the PLME instantaneous dephasing rates against those obtained from the numerically-exact evolution (see Sec. 5.2), for different values of the correlation time τ c . For the weak noise strength considered here (σ/Ω = 0.05), there is an excellent agreement even for times that are large compared to 1/Ω. Similar to our analysis of quasistatic noise, we can study the approximation error (see Eq. (37)) of the 2 ndorder PLME, and compare it to the simpler 0 th -order approximate master equation. This latter approximation equation is analogous to Eq. (38): one first derives a dephasing dissipator ignoring the Rabi driving in the Hamiltonian, and only adds back the Rabi drive as a Hamiltonian term in the final master equation (see Fig. 1(b)). Explicitly, this equation in the interaction frame reads: (47) Figure 5(b) compares the approximation errors for these two approaches as a function of time, for three different choices of the noise correlation time τ c . We see that the PLME has significantly smaller approximation error than the naive 0 th -order equation, even for relatively short correlation times.

Limiting long-time behavior
For evolution times t much longer than the noise correlation time, our instantaneous effective dephasing rates become time-independent, with the form: Note that in this limit, the second negative rate Γ − (t) remains non-zero. We can get further insight by considering the case of a weak drive, i.e., Ωτ c 1, in addition to taking the long time limit. Dropping terms of order O (Ωτ c ) 5 , we have andĤ In this weak Rabi-drive limit, the negative rate Γ − (t) vanishes as (Ωτ c ) 2 , whereas the positive rate tends to the expected dephasing rate for simple undriven dephasing (i.e., it is determined by the noise spectral density at zero frequency). In the opposite strong driving limit, where Ωτ c 1, at leading order, both Γ ± (t) are needed to capture the qubit's dynamics properly. In particular, we see that in this limit (and at long times t τ c ): whileĤ (53) In this case, |Γ + (t)| ≈ |Γ − (t)| and the dephasing jump operator associated with Γ − (t), can partially "undo" the loss generated by the jump operator of Γ + (t).
Heuristically, this is also a manifestation of the same kind of continuous dynamical decoupling we discussed in our analysis of quasistatic noise: the Rabi frequency is large enough to efficiently average away the comparatively slow dephasing noise.
7 Rabi driven qubit subject to 1/f noise We now investigate the PLME for a Rabi-driven qubit subject to classical noise with a 1/f spectral density, a kind of noise ubiquitous in many systems (e.g., superconducting qubits [9]), with substantial ongoing efforts to mitigate its effects (see e.g., [5,37,38]). The noise spectral density of interest has the form with ω l and ω h representing the low and high frequency cutoffs respectively. In the time domain, the corresponding autocorrelation function reads with Ci(t) = − ∞ t dt cos(t )/t . In the second line we have taken the limit of a large high frequency cutoff (i.e., ω h t → ∞), something we will do throughout this section.
Once again, using the discussion in Secs. 3 and 4, we can calculate the rates Γ ± (t) as well as the noise-induced Hamiltonian termĤ ren (t) that enter the PLME. We will also consider the case Ω ± ω l ≈ Ω to arrive at simpler analytic expressions (though this simplification is not a requirement of our treatment). We find: where, for brevity, we have introduced In Fig. 6(a), besides plotting an explicit form of φ(t), we also show that the analytical approximations for the rates Γ ± (t) exhibit good agreement with rates obtained from the numerically-exact evolution. As discussed in Dephasing rate (10  : Accuracy and performance of the 2 nd -order PLME for the case of 1/f noise. (a) PLME rates Γ±(t) (blue and orange) along with rates obtained from the true qubit evolution (black), here estimated by averaging over 20,000 noise realizations. The red dotted curve (range given on the right axis) shows the oscillatory nature of the angle φ(t). (b) Approximation error [Eq. (37)] of the 2 nd -order PLME (blue) and the 0 th -order approximation (orange) [Eq. (62)]. We stress that when numerically solving the PLME, we do not make the assumption that Ω ± ω l ≈ Ω, made to simplify the expressions in main text. Parameters are σ = 0.01Ω and ω l = 10 −3 Ω.
Sec. 4.2, at 4 th -order in the cumulant expansion one also obtain a time-dependent x-dephasing rate Γ x (t); its expression is unwieldy so we do not include it here. For the parameters used in Fig. 6, Γ x (t) is much smaller than the rates Γ ± (t) (hence we do not show it on this plot).
where in the last line, for completeness, we show the commonly used approximation obtained by only keeping the leading-order contribution in ω l t (note, that in simulations we use the form shown in the first line of Eq. (63)). From the plot, it is clear that the PLME is a more accurate description of the dynamics than the 0 th -order approximation. As was the case with the other highlycorrelated noise types, the 0 th -order approximation, completely misses the coherent renormalization Hamil-tonianĤ ren (t), the jump operator associated with the negative rate Γ − (t), as well as, overestimates the decoherence resulting from the positive rate Γ + (t).
It is tempting to assume that the effects of 1/f noise are similar to that of quasistatic noise, as both are correlated over long time-scales. Further, in the freeinduction decay limit (where Ω → 0), they both give rise to dephasing that depends exponentially on t 2 . Our analysis shows that the correspondence between these two kinds of noise is far from exact when we include a Rabi drive. For quasistatic noise (at the level of 2 ndorder cumulant expansion contributions), the two dephasing rates Γ ± (t) are non-decaying oscillatory functions, and have equal maximum amplitudes. Further, the angle φ(t) simply advances linearly with time.
Our results for 1/f noise show that the symmetry in amplitudes of the rates Γ ± (t) is broken. In particular, from Fig. 6(a), we see that in contrast to quasistatic noise, for 1/f noise the instantaneous dephasing rate Γ + (t) is more dominant (in amplitude) than the negative rate Γ − (t). Further, the oscillations in both Γ ± (t) as well as the angle φ(t) get smaller as time evolves (on relevant times scales, where ω l t 1). Accounting for these types of differences correctly may prove to be important in applications where accurate modeling of the effects of these different environment types on the qubit dynamics is required.

Conclusions
This work introduces and derives a new time-local pseudo-Lindblad master equation (PLME) that accurately models the effects of Gaussian non-Markovian classical noise on driven quantum systems. The approach is based on a generalized cumulant expansion, and yields a master equation that can be written in terms of time-dependent Hamiltonian and Lindblad dissipators. We find that the combination of non-Markovianity and driving leads generically to instantaneous dephasing rates that can be negative. Through a detailed comparison against numerics for a Rabi-driven qubit, we find that the PLME can give an extremely accurate description of dynamics over timescales typically required to perform qubit control operations. We considered three forms of non-Markovian noise: quasistatic, Lorentzian as well as 1/f , which is often dominant in various physical scenarios. As our PLME is time-local, and requires no stochastic averaging, it may help simplify and speed up accurate modeling of noisy quantum systems, that may otherwise be notoriously difficult to simulate.

A Generalized cumulant expansion
In this Appendix, we briefly outline the derivation of the cumulant expansion that we use to obtain the PLME. While many different versions and derivations of such expansions exist, we follow a particularly physically transparent one that is similar to Budimir and Skinner [19].
Our starting point is the general stochastic master equation shown in Eq. (1). We can formally solve it (for any given noise realization) as: where T represents the time-ordering operator, and where we have defined the superoperator L(t) via for some operatorQ. Next, one can formally perform an exact average over the Gaussian noise η(t ) yieldinĝ with where S(t) represents the spectral density of the bath. We next assume that U(t) is invertible, and use the above equation to derive a formally exact time-local equation for ∂ tρ (t) The generalized cumulant expansion now arises from a perturbative expansion of the superoperator R(t) in powers of the noise, i.e., in powers of η(t )L(t ), thus letting us write The leading non-zero term is at 2 nd -order in the noise. At this order, U −1 is approximated as unity, and one obtains Eq. (4) from the main text. As discussed, this is exactly the form of a Redfield master equation. This equation is also exact in the special cases where either A(t) commutes with itself at different times, or when the noise is delta correlated. In these cases, there is a cancellation at all higher orders between the expansion of U(t) and that of U −1 (t).
For the general case, the next non-vanishing contribution comes at 4 th -order in the noise, i.e., the fourth cumulant. It takes the general form: where, for notational convenience, we have introduced S ij to mean S(t i − t j ), with t 0 = t. It is easy to see that using Eq. (65), the above superoperator commutators, can be readily written in terms of commutators of the system Hamiltonian, and hence ofÂ(t). The non-commutativity ofÂ(t) at different times means that U −1 (t) does not cancel all contributions at this order; we see explicitly that all of the terms above involve commutators of the superoperator L(t ).

B Expectation value evolution
In the main text, in order to compare the accuracy of the PLME and the 0 th -order approximation to the  numerically-exact evolution, we studied the approximation error (t) (see Eq. (37)). Such a measure is conservative, but gives a useful lower bound to how well a given modeling approach performs.
To present the reader with a little more insight, in this Appendix we also show the expectation values σ y and σ z (noting that σ x has more trivial dynamics and thus for brevity, is not presented explicitly) as a function of time, of a particular initial state, here taken as |0 (positive energy eigenstate ofσ z ). Figures 7, 8  and 9 show the results for quasistatic, Lorentizan as well as 1/f noise, respectively. In each case, the evolution is obtained in the interaction frame (see discussed around Eq. (22)), and thus in the limit of vanishing noise strength (i.e., σ → 0) one should expect both σ y and σ z to stay constant as time evolves. We once again highlight the fact that the very model of a driven qubit we consider here, can be conveniently interpreted as a case of a continues dynamical-decoupling protocol: the transverse drive perpetually flips the qubit, and with each half-oscillation the effects of longitudinal noise can be partly "mitigated" (to what extent that happens largely depends on the noise type).
The shown plots complement our main text discussion: namely, when the noise is highly non-Markovian, the 0 th -order approximation provides a poor description of the evolution, while the PLME, in contrast, closely captures it both qualitatively and quantitatively. This is directly observed for the quasistatic, 1/f , as well as long-correlated (i.e., Ωτ c = 10) Lorentzian noise. As the noise correlation time decreases (see Fig. 8 panels (b), (e), and (c), (f)), the 0 th -order approximation starts to perform better, i.e., as τ c → 0, the system tends towards the white-noise limit, one where both approximations (PLME and the 0 th -order) reach parity (see discussion in Sec. 3 and Sec. 6). Finally we point out that, naturally, the exact form of the curves in Figs. 7 to 9 will change as one chooses a different initial state, the general trend discussed above, however, should hold even in those different cases.

C Detailed study of PLME positivity for quasistatic noise
In this Appendix we address the question of whether the evolution that our PLME leads to, is actually physical, i.e., can be described by a completely positive and trace-preserving (CPTP) map [39,28]? While in the case of a standard Lindblad master equation (even with time-dependent, but non-negative rates), this is guaranteed [1,28], we find that the PLME presented here, at times, can lead to dynamics that is not always completely positive. Even though this may seem as especially problematic, in this Appendix, we argue that this non-physicality of the 2 nd -and 4 th -order PLME evolution may be only relevant for some limited subset of initial states, and final evolution times, and hence in practice may not be especially problematic.
To explore this aspect of our theory further, we consider a Rabi-drive qubit, subject to quasistatic noise (see Sec. 5 for a detailed discussion) and describe the map associated with the 2 th -order PLME evolution in terms of a (numerically obtained) process matrix (often whereP k ∈ {1l,σ x ,σ y ,σ z } form an operator-basis for a single qubit. Any map is guaranteed to be completely positive as long as the corresponding χ-matrix is positive-semidefinite [28]. Hence, in Fig. 10 we show a plot of the smallest eigenvalue of χ as a function of the evolution time and noise strength. It is clear that there are certain times, at which the evolution from the PLME dynamics is not physical. As expected, at those times, the non-complete-positivity is more of a problem at higher noise strengths. This is because the leadingorder error that our approach makes in modeling the true evolution, scales as the quartic power of the noisestrength σ -see discussion in Sec. 5. We also observe Figure 10: The smallest eigenvalue of the process matrix χ(t) associated with the PLME evolution under quasistatic noise, plotted as a function of evolution time and the noise strength σ.
Negative regions correspond to evolution that is non-physical. For two particular configurations (black and orange dots on the density plot), we show a subset of initial pure states ("patches" on Bloch spheres) that result in non-physical final states. As one might expect, this subset of states is smaller at weaker noise strengths.
a periodic behavior that follows the Rabi drive period ∼ 1/Ω. In particular the PLME map starts to be noncompletely-positive approximately at times ∼ 2πm/Ω, with m ∈ Z + , and stays non-physical for roughly a half of a Rabi period. This behavior is mainly due the interplay between the dissipators associated with rates Γ ± (t). In particular, for some set of initial states, the negative rate may be a more dominant of the two, dictating much of the evolution, which can in fact at times "stretch" the qubit's state vector to be larger than 1, resulting in a potentially non-physical final state. Naturally this should not happen (and does not happen in the true dynamics of the system), however, our analytic form of Γ − (t) is only an approximation to the true rates of the instantaneous Liouvillian (see Fig. 2 of the main text). In Fig. 10 we also choose two particular maps (marked by black an orange dots on the density plot) and show on corresponding Bloch spheres a set of pure initial states (colored patches) that would lead to non-physical final states (i.e., density matrices that are not positivesemidefinite) when the given map is applied. We see that the subset of these problematic states is typically expected to be small, as long as the noise strength is not too large.
We stress that this type of artifact of our PLME is not unique when it comes to modeling the effects of highly non-Markovian noise. Other methods (for example, see [34]), can also result in dynamics that leads to final states which may not be fully physical. Ultimately, as long as the error made relative to the true system evolution is small enough, an approximation obtained from our PLME may still be a preferred over the much more widely used approaches, that can result physical final states, but at a cost of much less accurate dynamics. Furthermore, the evolution obtained with a PLME can be readily made physical, by projecting the final evolved state onto a "closest" physical state. Alternatively, a similar projection procedure can be also performed at the process map level. Such approaches can be readily implemented, and in practice can have little effect on the dominant error measures as long as the noise is not too strong.