Simple master equations for describing driven systems subject
to classical non-Markovian noisePeter Groszkowski1,2 , Alireza Seif1 , Jens Koch3 , and A. A. Clerk1
1 Pritzker School of Molecular Engineering, University of Chicago, Chicago, IL, USA
2 National Center for Computational Sciences, Oak Ridge National Laboratory, TN 37831, USA
3 Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA

Driven quantum systems subject to nonMarkovian noise are typically difficult to modeleven if the noise is classical. We present asystematic method based on generalized cumulant expansions for deriving a time-local masterequation for such systems. This master equation has an intuitive form that directly parallels a standard Lindblad equation, but contains
several surprising features: the combination ofdriving 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 ofa 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 employedphenomenological master equations which ignore the interplay between driving and noise.

1 Introduction
As quantum technology advances, it is becoming increasingly more important to have tools for modellingdissipation and noise in driven quantum systems thatare accurate, physically transparent, and easily implemented numerically. The standard approach is to describe dissipation using a time-local Lindblad masterequation (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 apPeter Groszkowski: piotrekg@gmail.com(a) Stochastic evolution(b) Oversimplified master equation

Drive

Drive

SystemNoise

SystemNoise(c) More accurate treatment
Modified Drive
System

Noise

(d)
Rates

arXiv:2207.03980v2 [quant-ph] 2 Apr 2023

April 2, 2023

Time

Figure 1: (a) Schematic of a driven qubit subject to classical,
non-Markovian noise η(t) via noise-coupling operator Â(t) withstrength g(t). (b) Standard approach to treating the noiseaveraged dynamics: one uses a master equation with a single
Lindblad dissipator D[Â(t)] with corresponding rate Γ(t) (derived in the absence of driving). The neglect of the interplay
between noise and driving can lead to large errors. (c) We
derive a substantially more accurate method for modeling dynamics via a pseudo-Lindblad master equation. In this case,the rates Γ± (t), the corresponding jump operators ê± (t), andthe system Hamiltonian are all modified by the drive. (d) Our
method results in instantaneous dephasing rates that can be
negative [see Eq. (15)].

proach generally requires an uncontrolled approximation.
While more complicated formalisms exist for nonMarkovian, 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 ofhighly-correlated noise using an ad-hoc Lindblad masterequation with either constant rates, or with rates thatare explicitly time-dependent (even when the noise isstationary). These approaches are common in superconducting quantum information (see e.g. [5, 6, 7, 8]),where the dominant dephasing noise typically has a 1/fcharacter [9, 10, 11, 12]. While these equations are often
physically well-motivated and retain the simplicity of an

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

1

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 notyet 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 atime-dependent but time-local master equation thatis reminiscent of a LME. We focus on the generalcase where the dissipation is due to classical nonMarkovian 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 othercontexts, 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 (eventhough there is no Lamb shift physics as the noise isclassical), 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 termsof generalized time-dependent rates provides a wealthof physical intuition. We call this resulting equation apseudo-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 thateven 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], ourwork differs in several respects. In particular, the useof these techniques to give a microscopic derivation of apseudo-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 expressedin pseudo-Lindblad form; however, unlike our work, theconnection to explicit microscopic models of noise anddriving 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 thepseudo-Lindblad master equation. We specialize to thespecific example of a Rabi-driven qubit subject to nonMarkovian 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
noiseConsider 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)Â. Here, g(t) represents a possibly time-dependent coupling constant.
Working in the interaction picture with respect to thenoise-free Hamiltonian Ĥ0 (t), the dynamics is describedby the stochastic HamiltonianĤ(t) = η(t)g(t)Â(t),(1)

where Â(t) = Û0 (t, t0 )† ÂÛ0 (t, t0 ). Here, we allow fora general time-dependent Hamiltonian H0 (t) that neednot commute with itself at different
 times. Therefore,RtÛ0 (t, t0 ) = T exp −i t0 Ĥ0 (s)ds , where T indicates
the time-ordering operation and t0 is an arbitrary reference time. This yields the evolution equation (setting~ = 1)
hi∂t ρ̂η (t) = −i Ĥ(t), ρ̂η (t) ,(2)where ρ̂η (t) is the system density matrix for a particularrealization 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 functionhη(t)η(t0 )i = S(t − t0 ).

(3)

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 (seee.g., [4]), or (b) [Â(t), Â(t0 )] 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).

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

2

3 Cumulant expansion and the pseudoLindblad 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 manyothers since (e.g., [15, 16, 17]). It is worth stressingthat such technique is closely related to work on timeconvolutionless (TCL) master equations; in particular
one can think of the cumulant-expansion method asa 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 (truncatedat a particular order in noise strength, and rewritten
in a pseudo-Lindblad form). The technique can beviewed as yielding the classical-noise limit of the timeconvolutionless master equation [4]. Hence, while thebasic 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 amodified system Hamiltonian and two dephasing rates(see Fig. 1). One of these rates is negative, which isa signature of the non-Markovianity of the noise process [22, 25, 4]. We point out that one common physicalinterpretation for these negative rates is related to thefact that (even a classical, as is the case here) nonMarkovian environment can lead to a revival of coherence: as the system evolves, the information loss to thebath can be at times (often partially) undone and “flowback” into the system. In setups where the bath isquantum, this can be often heuristically understood byconsidering a simple model where the bath is thoughtof as an auxiliary mode that couples to its own (often Markovian) environment, and which also coherently
couples to the original system in question. The revivalof 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 backthe original system, before it gets permanently lost to
the auxiliary mode’s environment.
As is standard, we also assume that at t = 0, theinitial 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.

3.1 2nd -order pseudo-Lindblad master equation
We start by truncating the generalized cumulant expansion at 2nd -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 :hhii∂t ρ̂(t) = −Λ(t) Â(t), Âavg (t), ρ̂(t) ,
(4)
where
Rt
Âavg (t) =0

dt1 S(t − t1 )g(t1 )Â(t1 )Rt
dt1 S(t − t1 )g(t1 )
0

(5)

Z t
dt1 S(t − t1 )g(t1 ),Λ(t) = g(t)(6)

0Both of these quantities have a simple physical interpretation:• Âavg (t) corresponds to the weighted average of Â(t)at times t0 < t. The weighting here is controlled by
the noise autocorrelation function.
• Λ(t) is the instantaneous dephasing rate we would
obtain if Â(t) = Â (i.e. if Â commutes with Ĥ0 (t)).
In that case, Âavg = Â, and Eq. (4) is a generalizeddephasing master equation that causes coherences
between eigenstates of Â to decay. For this specialcase, the cumulant expansion truncates at 2-orderand 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 isalways possible [1], although the rates associated with
dissipative processes will in general be time-dependentand possibly negative [4, 22]. Our goal is to show how
the parameters in this pseudo-Lindblad form reflect theproperties of our system and provide physical intuition.
The first step is to quantify both the size of the operator Âavg (t), and the degree to which it is orthogonal to Â(t) (the instantaneous coupling operator). We
achieve both using the Hilbert-Schmidt
inner product,

i.e., hX̂, Ŷ i = Tr X † Y and its corresponding norm||X̂||2HS = hX̂, X̂i. We first introduce the scalar quantities
Z0 (t) = ||Â(t)||HS ,(7)

Z(t) = ||Âavg (t)||HS ,(8)Z⊥ (t) = ||Âavg (t) − hb̂(t), Âavg (t)ib̂(t)||HS ,(9)

2 This interpretation is naturally only sensible if the denominator in Eq. (5) is non-zero for all times of interest. For moregeneral cases, one can simply set Λ(t) = 1, and define Âavg (t)with the denominator in Eq. (5) set to 1. All of our results below
then continue to apply.

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

3

as well as the normalized operators
1Â(t),Z0 (t)

1 
Âavg (t) − hb̂(t), Âavg (t)ib̂(t) .
b̂⊥ (t) =Z⊥ (t)b̂(t) =

(10)
(11)

The corresponding rates are given by

Defining the angle φ(t) viaZ⊥ (t) = Z(t) sin φ(t),e (cos φ(t) ± 1) .
Γ± (t) = Γ(t)(12)

we can now express Âavg in terms of operators paralleland perpendicular to the instantaneous system operatorÂ(t):
hiÂavg (t) = Z(t) cos φ(t)b̂(t) + sin φ(t)b̂⊥ (t) . (13)Thus even though our original stochastic Hamiltonianin Eq. (1) involved only a single system operator Â(t)at each instant, our dynamical equation involves two independent operators: Â(t) and the orthogonal operatorb̂⊥ (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 presenttime), and an effective rate Γ̃(t):e = 1 Λ(t)Z(t)Z0 (t)Γ(t)
2time-dependent operators ê± (t) that are rotated versions of b̂(t), b̂⊥ (t):
!



φ(t)
sin
ê+ (t)cos φ(t)
b̂(t)
2
2
.
=
ê− (t)b̂⊥ (t)− sin φ(t)cos φ(t)
2
2
(17)

(14)

This is just a normalized version of the instantaneousdephasing 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̂ † − 21 {X̂ † X̂, ρ̂} denote the standard Lindblad dissipative superoperator. A few lines of algebra transformEq. (4) into:hi∂t ρ̂(t) = − i Ĥren (t), ρ̂(t)X+2Γs (t)D [ês (t)] ρ̂(t).
(15)
s=+,−

Equation (15) is our general pseudo-Lindblad masterequation (PLME), and is a central result of our work.
The first term describes coherent dynamics generatedby Ĥren (t), a noise-induced correction to the system
Hamiltonian. It is given by the Hermitian operator:hie sin φ(t) b̂(t), b̂⊥ (t) .
Ĥren (t) = −iΓ(t)(16)
The second two terms of Eq. (15) describe generalized
time-dependent dephasing processes, with Hermitian(18)

Several comments are now in order: The Hamiltonian Ĥren (t) is not a standard Lamb-shift renormalization [4], as there is no quantum bath here (only classical noise). Heuristically, 2nd -order perturbation theory in the stochastic Hamiltonian Ĥ(t) could yield newprocesses 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, andshows that they are present whenever φ(t) 6= 0.
Furthermore, for any non-zero φ(t), we always have
two effective dephasing processes involving two orthogonal operators ê± (t). This reflects the fact that dissipation involves both the system operator Â(t) at thecurrent time, as well as at earlier times. For small φ(t),these operators are almost the same as b̂(t) and b̂⊥ (t).
Finally, for any non-zero φ(t) the effective instantaneous dephasing rate Γ− (t) is always negative. Whilesurprising, such negative dephasing rates are known toarise in problems where dephasing dynamics is nonmonotonic in time [28]. As we will show via comparison
to full numerical simulations, the negative dephasingrate we predict is not an artefact of our approximation,
but is instead in excellent agreement with features ofthe exact dissipative evolution.

3.2 Simple limits
It is reassuring to confirm that Eq. (15) reduces toknown results in some simple limits. First, consider thelimiting case of Markovian noise, where the noise autocorrelation function approaches a δ-function: S(t) →σδ(t). In this case, the bath has no memory, and therefore Âavg (t) = Â(t). As a result, Âavg (t) has no component orthogonal to Â(t), hence Z⊥ (t) and φ(t) vanish.
Our pseudo-Lindblad master equation thus reduces to:hi∂t ρ̂(t) = σD g(t)Â(t) ρ̂(t).
(19)
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 Â(t)commutes with itself at different times).

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

4

A second simple case is that of a time-independentÂ(t). Again, in this case φ(t) = 0, and we recover the
expected dephasing Lindblad equation:∂t ρ̂(t) = 4Γ̃(t)D[b̂]ρ̂(t) = 2Λ(t)D[Â]ρ̂(t),(20)

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 thesystem isĤrf (t) =Ω(t) xσ̂ + η(t)σ̂ z = Ĥ0 (t) + η(t)σ̂ z ,2(21)

where Ω(t) is the Rabi drive amplitude. This system
the form of our general problem, with g(t) = 1.
The corresponding interaction-frame system operatorÂ(t) is thenÂ(t) = cos(θ(t))σ̂ z + sin(θ(t))σ̂ y ,(22)Here φ(t) is given by Eq. (12). The corresponding jump
operators τ̂ u/v (t) are rotated Pauli matrices, and read
!
 u


e
e
τ̂ (t)cos θ(t)sin θ(t)σ̂ z=, (26)eeτ̂ v (t)σ̂ y− sin θ(t)cos θ(t)withe = φ(t) + θ(t).
θ(t)
(27)
2
√
Note that τ̂ u/v (t) are the same as 2ê± (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.

4.2 Corrections to the PLME from 4th -order cumulants
We can further improve the accuracy of our PLME by
including higher-order noise terms that arise at 4th order in the cumulant expansion. The 2nd -order equation, Eq. (24), has two dephasing dissipators involving
Pauli matrices in the z − y plane. As one might expect,
4th -order contributions modify the rates and forms ofthese dissipators. There is, however, a more interesting
effect that emerges at 4th -order: a new, third dissipatorthat corresponds to dephasing along the x axis (again,
working in the rotating frame used to write Eq. (21)).
This process corresponds to the term:withΓx (t)D[σ̂ x ]ρ̂(t).
Z t
θ(t) =

dt0 Ω(t0 ).

(23)

0

Note that our treatment below goes beyond previous
work on this problem that focused exclusively on theslow-noise limit [27, 29].
We can now directly apply the results of Sec. 3. Using
the Â(t) given by Eq. (22), the PLME (based on keeping contributions up to 2nd -order in the noise strength)
follows from Eq. (15). It takes a transparent form in
the rotating frame used to write H0 (t) in Eq. (21):hi∂t ρ̂(t) = − i Ĥren (t), ρ̂(t) + Γ+ (t)D [τ̂ u (t)] ρ̂(t)+ Γ− (t)D [τ̂ v (t)] ρ̂(t).

(24)

The instantaneous dephasing Γ± (t) are given by
Eq. (18), while from Eq. (16), the noise-induced Hamiltonian correction ise sin φ(t)σ̂ x .
Ĥren (t) = −Γ(t)

(25)

(28)

In App. A we give the general expression for the 4th order 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 2nd order noise generates a renormalization of the coherent system Hamiltonian Ĥ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 thedephasing associated with these fluctuations.

5 Rabi driven qubit subject to quasistatic noise
In this section, we apply the PLME description to a
Rabi-driven qubit subject to longitudinal classical noise(see Eq. (21)), and study what is perhaps the most nonMarkovian of all Gaussian noise processes: quasistatic

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

5noise, with an autocorrelation function that is independent of time:S(t) = σ 2 .

(29)

While highly non-Markovian (i.e., noise correlations
have an infinite time correlation), this noise is especiallysimple to work with, and results in compact analyticalexpressions. It also makes benchmarking our approachagainst 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 wellknown that the generalized cumulant expansion does
not converge if one considers extremely long evolution
times [15, 15]. This deficiency is however not an issue formost quantum information applications, as one wants tounderstand 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-intermediatetimescales, 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).

5.1 PLME parameters
5.1.1

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 2nd -order PLME of
Eq. (24) are given by Eq. (18):Γ± (t) =σ2(sin(Ωt) ± 2| sin(Ωt/2)|) ,Ω(30)

while the angle φ(t) that controls the form of the twodephasing 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 tn = 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 thenoise.
The remaining noise-induced term in the 2nd -order
PLME is the Hamiltonian renormalizationĤren (t) =σ2(1 − cos(Ωt)) σ̂ x .
Ωterm in Eq. (28). We find that the rate of this process
is always positive, but not periodic:(31)

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 4th -order cumulant corrections to the PLME, i.e., the induced x dephasingΓx (t) =

5.1.2

2σ 4
2Ωt − 2Ωt cos(Ωt)Ω3+ sin(2Ωt) − 2 sin(Ωt) .

(32)

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 )
52
Γ+ (t) ≈ 2σ t 1 − (Ωt) ,
48
2



(33)

(Ωt)2
,
16

(34)

Ωσ 2 t2 x
σ̂ .
2

(35)

Γ− (t) ≈ −2σ 2 twhile Ĥren (t), reduces toĤren (t) ≈To 0th -order in the drive amplitude, both the negativedephasing 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
2nd order in the drive, things become more interesting:the magnitude of Γ+ (t) is suppressed, and we generatea non-zero negative dephasing rate Γ− (t).

5.2 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 fulldissipative 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 ofthe dynamics via [28]Linst (t) ≈ lim

1

dt→0 dtV(t + dt)V −1 (t) − 1l .

(36)

As was done in Sec. 3, Linst (t) can be brought into a
Lindblad-like form. This lets us extract the generalized
instantaneous dephasing rates that characterize the fulldynamical map.

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

6Γ− (t)Γx (t)

5.0

100
Approximation error

Dephasing rate (10−3 Ω)

Γ+ (t)

2.5
0.0
−2.5
−5.0

10−2
10−4
PLME
10−6

0th -order Approx.
PLME (4th -order)

0

10

20
30
Time t (1/Ω)

1040

Figure 2: Time evolution of Γ± (t) (blue and orange) from the
2nd -order PLME of Eq. (15), as well as the rate Γx (t) (green),
calculated by keeping the 4th -order noise contributions in the
cumulant expansion. For comparison, the black curves are therates 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 thesimulations, 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 thecentral 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), itspresence can be important in accurately modeling thequbit’s evolution. Furthermore, is worth stressing thatthe exact dynamics also predicts a third rate (black dotted curve), which is accurately captured by Γx (t), obtained from keeping the 4th -order terms in the cumulant
expansion.

5.3 Benchmarking
thenumerically-exact results

PLME

againstIn order to asses the usefulness of the PLME, we need ameans to quantify its accuracy against the exact systemevolution. For our purposes, we will compare against
the numerically-exact evolution, obtained by numerically averaging unitary evolution over realizations of thestochastic 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 performanceof the approximation only for a small, specific subset of
initial states [33]. A more complete approach is to try tocompare the full structure of the exact and approximate
evolution maps (e.g., Ref. [34]).
Here, we will focus on a measure that characterizes

−8

0

10

20
30
Time t (1/Ω)

40

Figure 3: The accuracy of different approximations for calculating the qubit’s evolution under quasistatic noise. The approximation error of Eq. (37) is shown as a function of the evolution
time for the PLME when 2nd -order (blue) or 4th -order (green)terms are kept in the cumulant expansion, as well as the for the
0th -order approximation of Eq. (38) (orange). We see that bothversions of the PLME yield large accuracy improvements overthe naive 0th equation. All curves correspond to σ/Ω = 0.05.

the entire evolution map, and is not tied to a particular choice of initial states (although in App. B, forcompleteness, we also show the evolution of specificPauli operator expectation values given a particular initial state). For some evolution time t, let Jexact (t)denote the Choi matrix of the exact evolution map,
and Japprox (t) the Choi matrix of the evolution map
obtained via the approximation of interest (i.e., thePLME). We then define the approximation error at timet by:(t) = kJexact (t) − Japprox (t)k .

(37)

Here k...k 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 tothe single-shot channel distinguishability between these
two maps [36].  is 0 when the map associated with theapproximate dynamics cannot be distinguished from theexact map after a single evolution (no matter how theinitial state is optimized), and equals 2 when the twomaps are maximally distinguishable [36]. At a heuristic
level, it characterizes the worst case performance of ourapproximation, i.e., the largest error maximized over all
possible initial states.
With this definition in hand, we can now quantify theapproximation error of the PLME for different choicesof system parameters and noise spectra. To put theseresults in context, we will also consider the approximation error associated with a much simpler kind oftime-local master equation. In this approach (which we

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

7˙ = 2σ 2 tD [cos(Ωt)σ̂ z + sin(Ωt)σ̂ y ] ρ̂(t),
ρ̂(t)

(38)

We stress that this kind of approach is commonly usedbecause 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 wellas the naive 0th -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 2th -order PLME
explicitly shown in Eq. (24), here, for quasistatic noise,
we also include the performance of a PLME obtainedby keeping all the 4th order contributions, which resultin: corrections to the rates Γ± (t), the corresponding
jump operators (that rotate in the z − y plane), theHamiltonian renormalization Ĥ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 2nd - or 4th -order cumulant expansions)have a significantly smaller approximation error thanthe 0th -order equation. As expected, (t) for bothschemes grows with increasing time. Nonetheless, forweak noise, the PLME approximation error remains loweven if t  1/Ω. To be concrete, we see that for the chosen parameters, the 0th -order approach has an approximation error (t) < 10−3 only if tΩ . 1. In contrast,
the 2th -order PLME has (t) < 10−3 at times wheretΩ . 5, while the 4th -order PLME, over even longer
longer times, where tΩ . 17. We also point out thatin 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 diamondnorm 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 findfrom numerics that this quantity behaves analogouslyto (t) in the t → 0 limit. Further, the 1-norm metric isamenable to analytic estimates. Using this equivalence,
3 We consider a 1-norm of a superoperator associated with thedifference between the average evolution propagators of the maps
representing the two processes being compared. Here, Pa 1-norm of
a superoperator V(t) is defined as: ||V(t)||1 = maxj|v (t)|,
i ij
with vij (t) representing the ij-th element of V(t).

101
Approximation error

will term the 0th -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 0th -order master equation takes
the form:

10−1

PLME
0th -order Approx.
4
2
15 (σt) Ωt
2
2
3 (σt) Ωt

10−3
10−5
10−7 −2
10

10−1
100
Time t (1/Ω)

101

Figure 4: Approximation errors (t) of the 2nd -order PLME
(blue) and the 0th -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) associatedwith the 0th -order approximation scales as t3 . In contrast, for
the PLME is scales as t5 . One also finds a different asymptotic
scaling with the noise strength σ (σ 2 for the 0th -order equation,
σ 4 for the PLME. These results demonstrate the gain in accuracy associated with using the PLME). All curves correspondto σ/Ω = 0.2.

we find that for short times, the approximation error
for the 2nd -order PLME scales as:
2
PLME (t) ≈(σt)4 Ωt,(39)
15
whereas for the naive 0th -order equation we obtain2(σt)2 Ωt.
(40)
3
From Fig. 4 we see that Eqs. (39) and (40) accurately
described the asymptotic behavior of (t) in the small tlimit. The more favorable approximation error scalingof the PLME is not surprising, as in contrast to the0th -order approximation of Eq. (38), it fully captures
2nd -order noise contributions to the qubit’s dynamics.
Finally, while the PLME yields a highly accurate approximation, its form is not guaranteed to generate acompletely 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 andweak noise, as it is bounded by the error metrics wehave considered above. We discuss the positivity of ourapproximation in more detail in App. C.
0th (t) ≈6 Rabi driven qubit subject to noise witha finite correlation time
In this section we study the Rabi-driven qubit(see Eq. (21)), now subject to classical dephasing noise

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

8

0

4

− π2

3

−π

5

0

20

40

2

−5

0

10

20
30
Time t (1/Ω)

0.5
0

20

−1

40

0

10

20
30
Time t (1/Ω)

40

0

100

100

10−1

10−1

10−1

10−2

10−2

10−2

10−3

10−3

10−3

10−5

0

10

20
30
Time t (1/Ω)

40

10−4
10−5

Ωτc = 1
0

10

20

40

10

20
30
Time t (1/Ω)

40

0th -order Approx.

100Ωτc = 100Ωτc = 1/5Ωτc = 1

PLME

10−4

− π2

−π

40

0.0

(b)

Approximation error

0

− π2

−π

0
Ωτc = 10Γx (t)

1.0

10

Γ− (t)

0
φ(t)

φ(t)

Dephasing rate (10−3 Ω)

Γ+ (t)

10

φ(t)

(a)

20
30
Time t (1/Ω)

40

10−4
10−5

Ωτc = 1/5
0

10

20
30
Time t (1/Ω)

40

Figure 5: Accuracy and performance of the PLME for the case of Lorentzian noise. (a) 2nd -order PLME rates Γ± (t) (blue andorange) and Γx (t) (green curve), obtained by keeping 4th -order terms in the cumulant expansion (see App. D). The correlation
time τc is largest in the left panel and decreases towards the right. In particular, we have Ωτc = 10, (left column), Ωτc = 1 (middle
column) and Ωτc = 1/5 (right column). Black curves in each plot present rates obtained from the true system evolution (hereestimated by averaging over 20,000 noise realizations). The time evolution of the angle φ(t) (red) is shown in the insets. As thenoise becomes less correlated (i.e., Ωτc decreases), all rates and φ(t) show fewer oscillations and quickly approach constant values.
Furthermore, the negative rate Γ− (t) becomes smaller and less relevant, thus making Γ+ (t) a dominant source of incoherent
evolution. (b) The approximation error of Eq. (37) as a function of time quantifying the accuracy of the evolution obtained fromthe 2nd -order PLME (blue) or the 0th -order approximation of Eq. (47) (orange). In all plots the noise strength is taken to beσ = 0.05Ω.

that has a finite correlation time τc . We take η(t) to be
an Ornstein-Uhlenbeck process, implying:S(t) = σ 2 exp (−|t|/τc ) .

(41)

This yields a Lorentzian noise spectral density as encountered generically in various settings (i.e., classicaltelegraph noise produced by a two-state fluctuator, in
the Gaussian approximation). It also provides a convenient test case, as we can interpolate between thequasistatic noise of the previous section, and noise thatis Markovian, simply by tuning τc .

6.1 PLME parameters
The 2th -order PLME of Eq. (24) has the effective dephasing rates given by:
σ 2 τcΓ± (t) =e−t/τc (Ωτc sin(Ωt) − cos(Ωt))(Ωτc )2 + 1p+ 1 ± (Ωτc )2 + 11
× (e−2t/τc − 2e−t/τc cos(Ωt) + 1) 2 . (42)The jump operators corresponding to these ratesτ̂ u (t), τ̂ v (t) are determined by Eqs. (26) and (27), with
the non-Markovian angle φ(t) satisfyingcos φ(t) =Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

e−t/τc (Ωτc sin(Ωt) − cos(Ωt)) + 1
,
K

(43)

9

and
sin φ(t) =

e

−t/τc

(Ωτc cos(Ωt) + sin(Ωt)) − Ωτc, (44)
K

with
K=

q

(Ω2 τc2 + 1)(1 − 2 cos(Ωt)e−t/τc + e−2t/τc ).
(45)

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 noisee ≈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 Â(t), but
acting at a slightly “earlier time”.
We also obtain a Hamiltonian renormalization(στc )2 ΩĤren (t) =(Ωτc )2 + 1sin(Ωt)−t/τc× 1−ecos(Ωt) −σ̂ x .
Ωτc(46)
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) andĤren (t) vanish as expected.
Equation (42) implies that a 2nd -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 4th -order, one
also observes a third rate Γx (t), that corresponds to dephasing along the qubit’s x-axis. We explicitly show theanalytical form of Γx (t) in App. D.
Figure 5(a) compares the PLME instantaneousdephasing rates against those obtained from thenumerically-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 canstudy the approximation error (see Eq. (37)) of the 2nd order PLME, and compare it to the simpler 0th -orderapproximate master equation. This latter approximation equation is analogous to Eq. (38): one first derives a dephasing dissipator ignoring the Rabi drivingin the Hamiltonian, and only adds back the Rabi drive
as a Hamiltonian term in the final master equation (seeFig. 1(b)). Explicitly, this equation in the interaction
frame reads:∂t ρ̂(t) =2σ 2 τc 1 − e−t/τc× D [cos(Ωt)σ̂ z + sin(Ωt)σ̂ y ] ρ̂(t).

(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 0th -order equation, even for relatively
short correlation times.

6.2 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:
!
1
1
2
±p
. (48)
Γ± (t) ≈ σ τc(Ωτc )2 + 1(Ωτc )2 + 1
Note that in this limit, the second negative rate Γ− (t)remains non-zero.
We can get further insight by considering the case ofa weak drive, i.e., Ωτc  1, in addition to taking thelong time limit. Dropping terms of order O (Ωτc )5 ,we have


11
3
2
4
2
(49)
Γ+ (t) ≈ 2σ τc 1 − (Ωτc ) + (Ωτc ) ,416
Γ− (t) ≈ −2σ 2 τc




5
1
(Ωτc )2 − (Ωτc )4 ,
4
16

(50)

and

Ĥren (t) ≈ σ 2 τc Ωτc − (Ωτc )3 σ̂ x .

(51)

In this weak Rabi-drive limit, the negative rate Γ− (t)vanishes as (Ωτc )2 , whereas the positive rate tends tothe 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 thequbit’s dynamics properly. In particular, we see that inthis limit (and at long times t  τc ):


σ2
1
1
Γ± (t) ≈ ±
1±
−
,
(52)
Ω
Ωτc
2(Ωτc )2while
Ĥren (t) =

σ2
Ω


1−

1
1
+
2
(Ωτc )(Ωτc )4σ̂ x .

(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 samekind 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.

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

10

(a)

Γ+ (t)

4

2

S(t) =2σ (Ci (ωh |t|) − Ci (ωl |t|))

(55)

2

≈ −2σ Ci (ωl |t|)
R∞

0

0

0

with Ci(t) = − t dt cos(t )/t . In the second line wehave taken the limit of a large high frequency cutoff(i.e., ωh t → ∞), something we will do throughout thissection.
Once again, using the discussion in Secs. 3 and 4,
we can calculate the rates Γ± (t) as well as thenoise-induced Hamiltonian term Ĥren (t) that enter the
PLME. We will also consider the case Ω ± ωl ≈ Ω toarrive at simpler analytic expressions (though this simplification is not a requirement of our treatment). We
find:

p
2σ 2 
−Gu (t) ∓ sign(G(t)) Gu (t)2 + Gv (t)2 ,Γ± (t) =Ω(56)
andĤren =2σ 2
Gv (t)σ̂ x ,
Ω

(57)

φ(t)

0

3
2

−1

1
0

−2

−1
−2
−3

0

10

20
30
Time t (1/Ω)

40

−3

(b)

100
Approximation error

with ωl and ωh representing the low and high frequency
cutoffs respectively. In the time domain, the corresponding autocorrelation function reads

Dephasing rate (10−3 Ω)We now investigate the PLME for a Rabi-driven qubit
subject to classical noise with a 1/f spectral density, akind 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(
2π
, if ωl ≤ ω ≤ ωhσ 2 |ω|(54)
S[ω] =
0,otherwiseΓ− (t)

Angle φ(t)7 Rabi driven qubit subject to 1/f noise

10−1
10−2
10−3
PLME

10−4
10

−5

0th -order Approx.
0

10

20
30
Time t (1/Ω)

40

Figure 6: Accuracy and performance of the 2nd -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 2nd -order PLME (blue) and the0th -order approximation (orange) [Eq. (62)]. We stress thatwhen 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 Ω.

where, for brevity, we have introducedG(t) = Ci (ωl t) ωl t − sin (ωl t) ,(58)

Gu (t) = sin (Ωt) Ci (ωl t) − Si (Ωt) , Ω
Gv (t) = cos (Ωt) Ci (ωl t) + ln
− Ci (Ωt) ,
ωl

(59)

with Si(t) =Rt0cos φ(t) = p(60)

dt0 sin(t0 )/t0 . The angle φ(t) satisfies

Gu,
G2u + G2v

sin φ(t) = p

Gv. (61)
G2u + G2v

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 inSec. 4.2, at 4th -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 thisplot).
Figure 6(b) exhibits a performance comparison of the
2nd -order PLME and a 0th -order approximation masterequation, which (in the interaction frame) reads∂ ρ̂(t) =Γ0th (t)D [cos(Ωt)σ̂ z + sin(Ωt)σ̂ y ] ρ̂(t).

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

(62)
11

with
Γ0th (t) =4σ 2 tsin(ωl t)− Ci (tωl )ωl t≈4σ 2 t ln(ωl t),(63)

where in the last line, for completeness, we show thecommonly used approximation obtained by only keeping the leading-order contribution in ωl t (note, that insimulations we use the form shown in the first line of
Eq. (63)).
From the plot, it is clear that the PLME is a moreaccurate description of the dynamics than the 0th -orderapproximation. As was the case with the other highlycorrelated noise types, the 0th -order approximation,completely misses the coherent renormalization Hamiltonian Ĥren (t), the jump operator associated with thenegative 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 arecorrelated over long time-scales. Further, in the freeinduction decay limit (where Ω → 0), they both giverise to dephasing that depends exponentially on t2 . Our
analysis shows that the correspondence between thesetwo kinds of noise is far from exact when we include a
Rabi drive. For quasistatic noise (at the level of 2nd order 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 inamplitudes 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 (onrelevant times scales, where ωl t  1). Accounting forthese types of differences correctly may prove to be important in applications where accurate modeling of theeffects of these different environment types on the qubit
dynamics is required.

8 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 interms of time-dependent Hamiltonian and Lindblad

dissipators. We find that the combination of nonMarkovianity and driving leads generically to instantaneous dephasing rates that can be negative. Through adetailed comparison against numerics for a Rabi-driven
qubit, we find that the PLME can give an extremely accurate description of dynamics over timescales typicallyrequired to perform qubit control operations. We considered three forms of non-Markovian noise: quasistatic,
Lorentzian as well as 1/f , which is often dominant invarious 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 tosimulate.

Acknowledgments
This work was primary supported by the Departmentof Energy BES Quantum Information Science Program
under Award No. DE-SC002015. A.S. is supported bya Chicago Prize Postdoctoral Fellowship in Theoretical Quantum Science. J.K. was supported by the ARO
under Grant No. W911NF-1910016. A. A. C. acknowledges support from the Simons Foundation through aSimons Investigator award (Grant No. 669487, A. A.
C.).

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 physicallytransparent one that is similar to Budimir and Skinner[19].
Our starting point is the general stochastic masterequation shown in Eq. (1). We can formally solve it(for any given noise realization) as:Z tρ̂η (t) = T exp

0

0

0



dt η(t )L(t ) ρ̂η (0)

(64)

0

where T represents the time-ordering operator, andwhere we have defined the superoperator L(t) via
hi
L(t)Q̂ = −i g(t)Â(t), Q̂ ,(65)

for some operator Q̂. Next, one can formally perform
an exact average over the Gaussian noise η(t0 ) yielding

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

ρ̂(t) ≡ U(t)ρ̂(0).

(66)
12

(a)

with

 Z tZ t1dt1 dt2 S(t1 − t2 )L(t1 )L(t2 ) ,U(t) = T exp
2 0 0
(67)

(69)

The leading non-zero term is at 2nd -order in the noise.
At this order, U −1 is approximated as unity, and one
obtains Eq. (4) from the main text. As discussed, thisis exactly the form of a Redfield master equation. Thisequation is also exact in the special cases where eitherÂ(t) commutes with itself at different times, or when
the noise is delta correlated. In these cases, there is acancellation at all higher orders between the expansionof U(t) and that of U −1 (t).
For the general case, the next non-vanishing contribution comes at 4th -order in the noise, i.e., the fourth
cumulant. It takes the general form:
Z t
Z t1
Z t2
R(4) (t) =L(t)
dt1
dt2
dt3
(70)
0
0
0

× S02 S13 [L(t1 ), L(t2 )]L(t3 )+ S03 S12 L(t1 )[L(t2 ), L(t3 )]+ S03 S12 [L(t1 ), L(t3 )]L(t2 ) ,where, for notational convenience, we have introduced
Sij to mean S(ti − tj ), with t0 = t. It is easy to seethat using Eq. (65), the above superoperator commutators, can be readily written in terms of commutatorsof the system Hamiltonian, and hence of Â(t). The
non-commutativity of Â(t) at different times means thatU −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(t0 ).

B Expectation value evolution
In the main text, in order to compare the accuracyof the PLME and the 0th -order approximation to thePLME
0th -order Approx.

−0.15

PLME (4th -order)
exact

−0.20

(68)

The generalized cumulant expansion now arises froma perturbative expansion of the superoperator R(t) in
powers of the noise, i.e., in powers of η(t0 )L(t0 ), thus
letting us writeR(t) = R(2) (t) + R(4) (t) + . . .

−0.10

0

10

20
30
Time t (1/Ω)

40

(b)

1.0
0.8

PLME

0.6

PLME (4th -order)exact

0th -order Approx.

hσ̂z i∂t ρ̂(t) = [∂t U(t)] U −1 (t)ρ̂(t) ≡ R(t)ρ̂(t)−0.05

hσ̂y iwhere S(t) represents the spectral density of the bath.
We next assume that U(t) is invertible, and use theabove equation to derive a formally exact time-local
equation for ∂t ρ̂(t)

0.00

0.4
0.2
0.0
0

10

20
30
Time t (1/Ω)

40

Figure 7: Evolution of the qubit affected by quasistatic noise,
shown in terms of expectation values: panel (a) hσ̂ y i and (b)hσ̂ z i, as a function of time. Different curves are obtained using the PLME when 2nd -order (blue, solid) or 4th -order (green,dash-dotted) terms are kept in the cumulant expansion, the 0th order approximation (orange, dashed), as well as numericallyexact evolution (black, dotted). The qubit starts in an initialstate |0i (positive energy eigenstate of σ̂ z ), and the simulation is performed in the interaction frame (see Eq. (22)), withσ/Ω = 0.05.

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 agiven modeling approach performs.
To present the reader with a little more insight, in
this Appendix we also show the expectation values hσ̂ y i
and hσ̂ z i (noting that hσ̂ x i has more trivial dynamics
and thus for brevity, is not presented explicitly) as afunction of time, of a particular initial state, here taken
as |0i (positive energy eigenstate of σ̂ z ). Figures 7, 8
and 9 show the results for quasistatic, Lorentizan aswell 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 vanishingnoise strength (i.e., σ → 0) one should expect bothhσ̂ y i and hσ̂ z i to stay constant as time evolves. We once

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

13

(a)

(b)

(c)
0.00

0.00

Ωτc = 10−0.05

Ωτc = 1−0.02

hσ̂y i−0.15

0th -order Approx.
exact

−0.20
0

10

20
30
Time t (1/Ω)−0.004

−0.06

PLME

−0.006

−0.08
40

(d)

−0.10

0

10

20
30
Time t (1/Ω)

40

(e)
1.0

hσ̂z i

0.8
0.6
0.4

PLME

0.2

0th -order Approx.
exact
0

10

20
30
Time t (1/Ω)Ωτc = 1

0

10

20
30
Time t (1/Ω)

1.00

0.95

0.99

0.90

0.98

40

Ωτc = 1/5

0.97

0.85

0.96

0.80
40

−0.008(f)
1.00Ωτc = 10Ωτc = 1/5
−0.002

−0.04

−0.10

0.000

0

10

20
30
Time t (1/Ω)

40

0

10

20
30
Time t (1/Ω)

40

Figure 8: Evolution of the qubit affected by Lorentzian noise, shown in terms of expectation values: panels (a-c) show hσ̂ y i,while (d-f) hσ̂ z i, as a function of time. Different curves are obtained using the PLME (blue, solid), the 0th -order approximation(orange, dashed), as well as numerically-exact evolution (black, dotted), obtained from averaging over 20000 noise realizations.
The correlation time τc is largest in the left panels and decreases towards the right. In particular, we have in (a), (d): Ωτc = 10,(b), (e): Ωτc = 1 and (c), (f): Ωτc = 1/5. The qubit starts in an initial state |0i (positive energy eigenstate of σ̂ z ), and thesimulation is performed in the interaction frame (see Eq. (22)), with σ = 0.05Ω.

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 longitudinalnoise 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 0th -order approximation provides a poor descriptionof the evolution, while the PLME, in contrast, closely
captures it both qualitatively and quantitatively. Thisis directly observed for the quasistatic, 1/f , as well as
long-correlated (i.e., Ωτc = 10) Lorentzian noise. As thenoise correlation time decreases (see Fig. 8 panels (b),(e), and (c), (f)), the 0th -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 0th -order) reach parity (seediscussion 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 holdeven 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 andtrace-preserving (CPTP) map [39, 28]? While in thecase 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, attimes, can lead to dynamics that is not always completely positive. Even though this may seem as especially problematic, in this Appendix, we argue that thisnon-physicality of the 2nd - and 4th -order PLME evolution may be only relevant for some limited subset of
initial states, and final evolution times, and hence inpractice 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 themap associated with the 2th -order PLME evolution interms of a (numerically obtained) process matrix (often

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

14

(a)

0.000

hσ̂y i−0.025−0.050PLME−0.0750th -order Approx.
exact−0.1000102030Time t (1/Ω)

40

(b)

1.0

hσ̂z i

0.8
0.6
PLME
0.4

0th -order Approx.
exact

0.2
0

10

20
30
Time t (1/Ω)

40

Figure 9: Evolution of the qubit affected by 1/f noise, shown interms of expectation values: panel (a) hσ̂ y i and (b) hσ̂ z i, as afunction of time. Different curves are obtained using the PLME
(blue, solid), the 0th -order approximation (orange, dashed), as
well as numerically-exact evolution (black, dotted), obtained
from averaging over 20000 noise realizations. The qubit starts
in an initial state |0i (positive energy eigenstate of σ̂ z ), and thesimulation is performed in the interaction frame (see Eq. (22)),with σ = 0.01Ω and ωl = 10−3 Ω.

referred to as the χ-matrix) which readsX
Et (ρ) →χnm (t)P̂n ρ̂P̂m ,

(71)

n,m

where P̂k ∈ {1l, σ̂ x , σ̂ y , σ̂ z } form an operator-basis fora 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 thosetimes, 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 thetrue 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 thedensity 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 weakernoise 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 halfof 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, thenegative 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 ratesof the instantaneous Liouvillian (see Fig. 2 of the maintext).
In Fig. 10 we also choose two particular maps (marked
by black an orange dots on the density plot) and showon corresponding Bloch spheres a set of pure initialstates (colored patches) that would lead to non-physicalfinal states (i.e., density matrices that are not positivesemidefinite) when the given map is applied. We seethat the subset of these problematic states is typicallyAccepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

15

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 ofhighly 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 overthe 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 thefinal 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 havelittle effect on the dominant error measures as long as
the noise is not too strong.

D Dephasing rate Γx (t) for Lorentziannoise
In this brief Appendix, we write down the complete
expression for the Γx (t) x-dephasing rate (calculated
by keeping 4th -order terms in the cumulant expansion)for the case of Lorentzian noise. It reads 5 2 2 22tΩ τc + 5
τc3 e− τc
4 τc Ω−
Γx (t) =2σ33 β1(Ω2 τc2 + 1)(Ω2 τc2 + 1)
t2τc2 e− τc
−
β2
(72)
2
(Ω2 τc2 + 1) Ωwith
β1 =Ω4 τc4 − sin (2tΩ) Ω3 τc3 + Ω2 (3 cos (2tΩ) + 2) τc2+ 3 sin (2tΩ) Ωτc − cos (2tΩ) + 1,

(73)

and
β2 =Ω2 τc2 (cos (tΩ) Ωt + sin (tΩ)) + 2Ω2 tτc sin (tΩ)− Ωt cos (tΩ) + sin (tΩ) .

(74)

As we expect, it scales with the fourth power of thenoise strength σ.

References[1] Vittorio Gorini, Andrzej Kossakowski, and E. C. G.
Sudarshan. “Completely positive dynamical semigroups of n-level systems”. Journal of Mathematical Physics 17, 821–825 (1976).

[2] G. Lindblad. “On the generators of quantum dynamical semigroups”. Communications in Mathematical Physics 48, 119–130 (1976).
[3] J. R. Johansson, P. D. Nation, and Franco Nori.
“Qutip: An open-source python framework for the
dynamics of open quantum systems”. Computer
Physics Communications 183, 1760–1772 (2012).
[4] Heinz-Peter Breuer and Francesco Petruccione.
“The theory of open quantum systems”. Oxford
University Press, USA. (2002).
[5] Nicolas Didier, Eyob A. Sete, Joshua Combes, andMarcus P. da Silva. “ac flux sweet spots in parametrically modulated superconducting qubits”.
Phys. Rev. Applied 12, 054015 (2019).
[6] Agustin Di Paolo, Thomas E. Baker, Alexandre
Foley, David Sénéchal, and Alexandre Blais. “Efficient modeling of superconducting quantum circuits with tensor networks”. npj Quantum Information 7, 1–11 (2021).
[7] M. A. Yurtalan, J. Shi, G. J. K. Flatt, and A. Lupascu. “Characterization of multilevel dynamics
and decoherence in a high-anharmonicity capacitively shunted flux circuit”. Physical Review Applied 16, 054051 (2021).
[8] F. Setiawan, Peter Groszkowski, Hugo Ribeiro, andAashish A. Clerk. “Analytic design of accelerated
adiabatic gates in realistic qubits: General theory and applications to superconducting circuits”.
PRX Quantum 2, 030306 (2021).
[9] G. Ithier, E. Collin, P. Joyez, P. J. Meeson,
D. Vion, D. Esteve, F. Chiarello, A. Shnirman,
Y. Makhlin, J. Schriefl, and G. Schon. “Decoherence in a superconducting quantum bit circuit”.
Phys. Rev. B 72, 134519 (2005).
[10] Alexandre Blais, Arne L. Grimsmo, Steven M.
Girvin, and Andreas Wallraff. “Circuit quantum
electrodynamics”. Reviews of Modern Physics 93,
025005 (2021).
[11] E. Paladino, Y. M. Galperin, G. Falci, and B. L.
Altshuler. “1/f noise: Implications for solid-state
quantum information”. Reviews of Modern Physics
86, 361 (2014).
[12] E. Paladino, L. Faoro, G. Falci, and Rosario Fazio.
“Decoherence and 1/f noise in josephson qubits”.
Physical review letters 88, 228304 (2002).
[13] Tobias Becker, Ling-Na Wu, and André Eckardt.
“Lindbladian approximation beyond ultraweak
coupling”. Physical Review E 104, 014110 (2021).
[14] Ryogo Kubo. “Generalized cumulant expansion
method”. Journal of the Physical Society of Japan
17, 1100–1120 (1962).

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

16

[15] N. G. Van Kampen. “A cumulant expansion forstochastic linear differential equations. i”. Physica
74, 215–238 (1974).
[16] N. G. Van Kampen. “A cumulant expansion forstochastic linear differential equations. ii”. Physica
74, 239–247 (1974).
[17] Ronald Forrest Fox. “Critique of the generalized
cumulant expansion method”. J. Math. Phys. 17,
7 (1976).
[18] Gulácsi, Balázs and Burkard, Guido. “Smokinggun signatures of non-Markovianity of a superconducting qubit”. arXiv:2302.09092 (2023).
[19] J. Budimir and J. L. Skinner. “On the relationship between T1 and T2 for stochastic relaxation models”. Journal of Statistical Physics 49,
1029–1042 (1987).
[20] M. Aihara, H. M. Sevian, and J. L. Skinner. “Nonmarkovian relaxation of a spin-1/2 particle in afluctuating transverse field: Cumulant expansion
and stochastic simulation results”. Physical Review A 41, 6596–6601 (1990).
[21] Yuan-Chi Yang, S. N. Coppersmith, and MarkFriesen. “Achieving high-fidelity single-qubit gatesin a strongly driven charge qubit with 1/f charge
noise”. npj Quantum Information 5, 12 (2019).
[22] Michael J. W. Hall, James D. Cresser, Li Li, andErika Andersson. “Canonical form of master equations and characterization of non-markovianity”.
Phys. Rev. A 89, 042120 (2014).
[23] John H. Reina, Luis Quiroga, and Neil F. Johnson.
“Decoherence of quantum registers”. Phys. Rev. A
65, 032326 (2002).
[24] G. Massimo Palma, Kalle-Antti Suominen, and Artur Ekert. “Quantum computers and dissipation”.
Proceedings of the Royal Society of London. Series
A: Mathematical, Physical and Engineering Sciences 452, 567–584 (1996).
[25] Heinz-Peter Breuer, Bernd Kappler, and Francesco
Petruccione. “Stochastic wave-function method fornon-markovian quantum master equations”. Physical Review A 59, 1633 (1999).
[26] S. Chaturvedi and F. Shibata.
“Timeconvolutionless projection operator formalism
for elimination of fast variables. Applications
to Brownian motion”. Zeitschrift für Physik B
Condensed Matter 35, 297–308 (1979).
[27] V. V. Dobrovitski, A. E. Feiguin, R. Hanson, andD. D. Awschalom. “Decay of rabi oscillations
by dipolar-coupled dynamical spin environments”.
Phys. Rev. Lett. 102, 237601 (2009).

[28] Heinz-Peter Breuer, Elsi-Mari Laine, Jyrki Piilo, and Bassano Vacchini. “Colloquium: Nonmarkovian dynamics in open quantum systems”.
Rev. Mod. Phys. 88, 021002 (2016).
[29] K. Rabenstein, V. A. Sverdlov, and D. V. Averin.
“Qubit decoherence by gaussian low-frequency
noise”. Journal of Experimental and Theoretical
Physics Letters 79, 646–649 (2004).
[30] P. Facchi, D. A. Lidar, and S. Pascazio. “Unification of dynamical decoupling and the quantum
zeno effect”. Phys. Rev. A 69, 032314 (2004).
[31] J-M Cai, Boris Naydenov, Rainer Pfeiffer, Liam P.
McGuinness, Kay D. Jahnke, Fedor Jelezko, Martin B. Plenio, and Alex Retzker. “Robust dynamical decoupling with concatenated continuous driving”. New Journal of Physics 14, 113023 (2012).
[32] Kevin C. Miao, Joseph P. Blanton, Christopher P. Anderson, Alexandre Bourassa, Alexander L Crook, Gary Wolfowicz, Hiroshi Abe,
Takeshi Ohshima, and David D. Awschalom. “Universal coherence protection in a solid-state spin
qubit”. Science 369, 1493–1497 (2020).
[33] Haimeng Zhang, Bibek Pokharel, E. M. LevensonFalk, and Daniel Lidar.
“Predicting nonmarkovian superconducting-qubit dynamics fromtomographic reconstruction”. Physical Review Applied 17, 054018 (2022).
[34] Richard Hartmann and Walter T. Strunz. “Accuracy assessment of perturbative master equations:Embracing nonpositivity”. Physical Review A 101,
012103 (2020).
[35] Dorit Aharonov, Alexei Kitaev, and Noam Nisan.
“Quantum circuits with mixed states”. In Proceedings of the thirtieth annual ACM symposium on
Theory of computing. Pages 20–30. (1998).
[36] John Watrous. “The theory of quantum information”. Cambridge University Press. USA (2018).
1st edition.
[37] Ziwen Huang, Pranav S Mundada, András Gyenis, David I Schuster, Andrew A Houck, and Jens
Koch. “Engineering dynamical sweet spots to protect qubits from 1/f noise”. Physical Review Applied 15, 034065 (2021).
[38] Ziwen Huang, Xinyuan You, Ugur Alyanak,
Alexander Romanenko, Anna Grassellino, andShaojiang Zhu. “High-order qubit dephasing at
sweet spots by non-gaussian fluctuators: Symmetry breaking and floquet protection”. Physical Review Applied 18, L061001 (2022).

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

17

[39] M. A. Nielsen and I. L. Chuang. “Quantum computation and quantum information”. Cambridge
University Press. Cambridge, UK (2000).

Accepted in Quantum 2023-03-02, click title to verify. Published under CC-BY 4.0.

18