Dynamical decoupling of a qubit with always-on control fields

Quantum information processing seeks to exploit quantum mechanical systems to store and manipulate information in novel ways. Typically, however, this quantum behavior exists at the microscopic level, such as a single atom or electron, making control by the experimenter challenging. Moreover, such a system is never perfectly isolated. Stray electromagnetic fields, thermal fluctuations, etc can perturb the quantum information carrier, imparting noise into the system. This noise can disrupt the quantum information processor and negate any benefits to operating in a quantum system.

Quantum states can be protected from noise with active control fields. Before the field of quantum information emerged, a wide variety of techniques for dynamical decoupling (DD) of a spin ensemble were developed in the nuclear magnetic resonance (NMR) community. The first discovery was the Hahn spin echo sequence [1]. By repeating the spin echo sequence several times in succession, one finds the Carr–Purcell sequence [2]. This sequence was soon improved by selectively choosing the axes around which one rotates the spin vector, yielding the workhorse Carr–Purcell–Meiboom–Gill sequence [3]. Ever more complicated sequences were developed in ensuing decades, allowing detailed characterization or manipulation of the environments of nuclear spin ensembles, or the compensation against inhomogeneous control fields.

DD found a renewed purpose in quantum information science. For example, quantum error correction systems [4–10] demand an error-per-operation of less than about 1% [7, 11]. Such low error rates can be difficult to achieve in practice since many candidate physical systems will lose their qubit coherence very quickly owing to environment-induced decoherence; as a result, quantum error correction may fail due to short qubit lifetimes. DD schemes offer a relatively simple solution to this problem: by applying a determined sequence of control pulses, one can significantly enhance the lifetime of a qubit.

Early work in DD sought to suppress decoherence by applying a periodic sequence of instantaneous 'bang–bang' pulses, which periodically flip the state of a qubit and undo the coupling to the bath [12–15]. This work was extended to Hamiltonian engineering by using decoupling sequences to selectively enable coupling Hamiltonians [16, 17]. Subsequently, attention turned to using Eulerian graphs to design sequences using bounded control operations that were robust to many types of systematic control errors [18]. This led to the notion of dynamically-corrected gates [19], which combine techniques from DD and composite pulse sequences from NMR to produce error-suppressing quantum gates.

The idealized 'bang–bang' control pulse is instantaneous, like the Dirac delta function. However, real physical pulses can only approximate 'bang–bang' control, because such idealized control would require infinite power. The implications of a real, continuous-time function can be significant. Most DD sequences are designed to correct noise errors between pulses, but offer no intrinsic protection to errors during pulses. Early studies of always-on DD arose in the problem of dipolar decoupling in solid-state NMR [20], culminating from substantial work in extending pulsed decoupling techniques to 'soft', or low-power, pulses. A modern consideration of this problem in the quantum information context is found in [21]. More recently, optimal control approaches to customizing continuous, bounded controls to decouple a specific but arbitrary noise bath were considered in [22–24]. Of particular relevance to the present work are approaches that extend optimized bang–bang sequences such as the Uhrig dynamical decoupling (UDD) sequence [25, 26], which is derived by canceling general noise Hamiltonians order-by-order in perturbation theory. Extensions of UDD to continuous and bounded control were studied by Uhrig and Pasini [27–29]. This work introduced pulses that correct errors during their own evolution to finite order in a time-dependent perturbation expansion. These sequences are termed 'realistic UDD', or RUDD. Subsequent work showed how to produce arbitrary gates decoupled from environment noise using bounded-strength controls [30, 31], as well as the application of modifications based on optimal control given known system drifts [32]. DD with such 'finite-width' pulses represents efforts to improve existing sequences by accounting for the real continuous-time nature of control fields.

DD sequences where the control field is always or almost always on pose some unique challenges for sequence engineering and characterization. Methods for characterization that are exact in the 'bang–bang' limit of DD and approximately correct with short but finite pulse widths require continuous analogues. In this work, we examine methods to both create and characterize always-on sequences, by considering continuous extensions of the UDD sequence. We introduce a new DD sequence called Chebyshev and Fourier expansion (CAFE), and compare it to similar RUDD sequences.

The paper is organized as follows. Section 2 gives a derivation of the theory behind CAFE, presenting it as a continuous extension of the bang–bang UDD sequence. Ultimately, CAFE is derived from a particular heuristic rather than from a fundamental determination of optimality, but this is typical of continuous sequences. Critical for evaluating the heuristic, then, are methods of characterizing the sequence, and careful choices of sequences for comparison, as discussed in section 3. Section 4 tests the CAFE sequence as well as RUDD sequences for performance in decoupling a qubit from a simulated quantum bath, represented by a small dipole-coupled spin lattice. Section 5 studies the frequency response of CAFE and compares analytical and simulated filter functions for this sequence and RUDD sequences. Section 6 simulates the scenario where a qubit is coupled to semiclassical noise and where control fields for decoupling are noisy, again for both CAFE and RUDD sequences.

We consider a qubit that is exposed to an environment that causes phase errors, or a dephasing bath. For simplicity, we parameterize our system in terms of Pauli operators σ^j, so that the phase errors correspond to σ^z operators and the control of the qubit is proportional to the σ^x operator. The total system Hamiltonian is

$$\begin{equation} \mathcal{H}(t) = \mathbbm{1} \otimes B_0 + \sigma^z \otimes B_Z + \alpha(t) \sigma^x \otimes \mathbbm{1}, \end{equation} \noindent \tag{ 1 }$$

where $\mathbbm {1}$ refers to identity in its subspace. In this expression, $\mathbbm {1}\otimes B_0$ is a pure bath evolution term, σ^z⊗B_Z is the qubit–bath coupling term, and $\alpha (t) \sigma ^x \otimes \mathbbm {1}$ is the time-dependent term representing the experimenter's control. In the following analysis, we find it useful to define the time-integrated control function

$$\begin{equation} \beta(t) = 2 \int_0^t \alpha(s)\,\mathrm{d}s. \end{equation} \noindent \tag{ 2 }$$

The quantity β(t) has an intuitive meaning—it represents the total angle (commonly called 'pulse area') that the qubit state is rotated on the Bloch sphere, in this case around the X-axis. We now transform our system to the interaction picture, so that the interaction Hamiltonian becomes

$$\begin{equation} \mathcal{H}_{\mathrm{int}}(t) = \cos\left[\beta(t)\right] \sigma^z \otimes B(t) + \sin\left[\beta(t)\right] \sigma^y \otimes B(t), \end{equation} \noindent \tag{ 3 }$$

where

$$\begin{equation} B(t) = \mathrm{e}^{\mathrm{i}B_0 t} B_Z \,\mathrm{e}^{-\mathrm{i}B_0 t}. \end{equation} \noindent \tag{ 4 }$$

We see that, in the interaction picture, the system–bath coupling terms are functions of β(t). In the case of 'bang–bang' π-pulse control, β(t) is a piecewise-constant function that is everywhere a multiple of π, so that $\sin \left [\beta (t)\right ] = 0$ everywhere and we can replace $\cos \left [\beta (t)\right ]$ with a 'switching function' y(t) that takes on the values ±1. For example, this is the formalism considered in [25]. We will later refer to the theory behind the bang–bang version of UDD to derive always-on CAFE sequences.

2.1. Dynamical decoupling by order in T

Our task is to develop a class of control functions β(t) so that the system evolution of the qubit over time interval t∈[0,T], given by

$$\begin{equation} \mathcal{U}(0,t) = \mathcal{T}\exp\left(-\mathrm{i}\int_0^t \mathcal{H}(\tau) \,\mathrm{d}\tau\right), \end{equation} \noindent \tag{ 5 }$$

closely approximates the identity operation on the qubit (to within a complex phase factor) at the final time t = T, where $\mathcal {T}$ is the time-ordering operator. We approach this problem by expanding the interaction picture propagator in the Magnus expansion, as follows:

$$\begin{equation} \mathcal{U}_{\mathrm{int}}(0,T) = \exp\left(-\mathrm{i} \sum_{k=0}^\infty \Omega_k\right), \end{equation} \tag{ 6 }$$

where

$$\begin{equation} \Omega_1 = \int_0^T \mathcal{H}_{\mathrm{int}}(s)\,\mathrm{d}s, \end{equation} \tag{ 7 }$$

$$\begin{equation} \Omega_2 = \frac{1}{2} \int_0^T \int_0^{s_1} [\mathcal{H}_{\mathrm{int}}(s_1),\mathcal{H}_{\mathrm{int}}(s_2)] \,\mathrm{d}s_1\, \mathrm{d}s_2, \end{equation} \tag{ 8 }$$

and so forth. To simplify analysis, we transform the integration variables to dimensionless form as

$$\begin{equation} u_j = \frac{2}{T}s_j - 1, \end{equation} \tag{ 9 }$$

so that the domain of integration is now [ − 1,1]. As a result, Ω_n∝Tⁿ, so the expansion in (6) is a power series in T.

We say that a sequence which satisfies Ω_j ≃ 0 for j = 1,...,n 'decouples to nth-order' because the first nonvanishing error is of order $T^{n+1} \left \| \mathcal {H}_{\mathrm {int}} \right \|^{n+1}$. Approximating equality with zero in each term is necessary because the Ω_j terms each depend on integrals over successively higher-order commutators of the interaction-picture noise term B(t) from (4), a stochastic time-dependent operator. The random degrees of freedom in B(t) require us to summarily state that the Ω_j up to some order n are independently nearly zero in order to achieve decoupling for all possible instances of B(t). We use approximate equality because the CAFE sequence we will derive is ultimately limited by numerical precision and physically realizable control fields, which may violate strict equality of the 'decoupling constraints'. However, we defend the introduction of any approximation errors with our simulation results. This approach to the DD problem is akin to time-dependent perturbation theory, and is useful when the noisy environment is not strongly perturbing the system and hence the Magnus expansion (in the interaction picture defined above) converges rapidly.

The 'zeroth-order' constraint which must be satisfied is that the control field itself (in absence of any bath effects) is identity. This could be modified, if desired, to give some single-qubit rotation, but we do not pursue this possibility here. Hence,

$$\begin{equation} \beta(T) = 0 \; (\mathrm{mod} \; 2 \pi). \end{equation} \tag{ 10 }$$

This additional constraint must be appended because the initial set of constraints was derived in the interaction picture. For simplicity, we will use β(T) = 0 strictly in this work.

The first- and second-order decoupling constraints can be derived from (3) and (6) by nullifying any terms which are not identity operations on the qubit. For example, the first-order constraints (∝T) are

$$\begin{equation} \int_{-1}^{1} \sin [\beta (u)]\,\mathrm{ d}u =0, \end{equation} \tag{ 11 }$$

$$\begin{equation} \int_{-1}^{1} \cos [\beta (u)]\,\mathrm{ d}u =0. \end{equation} \tag{ 12 }$$

The second-order constraints (∝T²) can be derived from the Magnus expansion as

$$\begin{equation} \int_{-1}^{1} \int_{-1}^{u_{1}} \left\{\sin [\beta (u_{1})] - \sin [\beta (u_{2})] \right\}\,\mathrm{d}u_{1}\,\mathrm{d}u_{2} =0, \end{equation} \tag{ 13 }$$

$$\begin{equation} \int_{-1}^{1} \int_{-1}^{u_{1}} \left\{\cos [\beta (u_{1})] - \cos [\beta (u_{2})] \right\} \,\mathrm{d}u_{1}\, \mathrm{d}u_{2} =0, \end{equation} \tag{ 14 }$$

$$\begin{equation} \int_{-1}^{1} \int_{-1}^{u_{1}} \sin \left[\beta (u_{1}) - \beta (u_{2})\right] \,\mathrm{d}u_{1} \,\mathrm{d}u_{2}=0. \end{equation} \tag{ 15 }$$

These continuous DD constraints are equivalent to those studied elsewhere, for example in [27].

2.2. A continuous control function that mimics Uhrig dynamical decoupling

Finding arbitrary time-dependent control functions that satisfy the continuous DD constraints, even at first-order, is not trivial. We begin our search by studying the cumulative pulse area β_UDD(t) for the UDD sequence because this sequence is known to decouple by perturbative order using bang–bang pulses [33]. By transforming to the domain u∈[ − 1,1], the time-dependent rotation in a UDD sequence of N pulses is given by

$$\begin{eqnarray} \beta_{\mathrm{UDD}} (u) & = \pi \int_{-1}^{u} \sum_{j = 1}^{N} \delta\left[u_1 + \cos\left(\frac{j \pi}{N+1}\right)\right] \mathrm{d}u_1 \nonumber \\ &= \pi \sum_{j = 1}^{N} \Theta\left[ u + \cos\left(\frac{j\pi}{N+1}\right)\right], \end{eqnarray} \tag{ 16 }$$

where δ(u) is the Dirac delta function and Θ(u) is the Heaviside step function.

The integrated control function for UDD, β_UDD(u), can be well approximated by a smooth function of the form

$$\begin{equation} \beta_{\mathrm{s}} (u) = (N+1) \cos ^{-1} (-u), \end{equation} \tag{ 17 }$$

where N corresponds to the number of pulses in the UDD sequence, which is plotted in figure 1 for N = 6. Note that Yang and Liu [33] make a convenient variable substitution, which is to operate in a transformed coordinate basis given by θ = cos⁻¹(−u). In this basis, β_s is proportional to θ. We utilize this transformation later. By differentiation of (17), one can determine that the corresponding control function is given by

$$\begin{equation} \alpha_{\mathrm{s}} (u) = \frac{N+1}{2 \sqrt{1 - u^{2}}}. \end{equation} \tag{ 18 }$$

Moreover, the interaction picture Hamiltonian can now be expressed in terms of Chebyshev polynomials when using this continuous control function since

$$\begin{equation} \cos [(N+1) \cos ^{-1} (-u)] = T_{N+1} (u), \end{equation} \tag{ 19 }$$

$$\begin{equation} \sin [(N+1) \cos ^{-1} (-u)] = \sqrt{1-u^{2}} U_{N} (u), \end{equation} \tag{ 20 }$$

where T_N+1 is the (N + 1)th-order Chebyshev polynomial of the first kind and U_N is the Nth-order Chebyshev polynomial of the second kind. The latter is particularly noteworthy because $\sqrt {1-u^{2}}$ is the weight function for U_N, so that $\sqrt {1-u^{2}} U_{N} (u)$ is orthogonal to any polynomial in u of order N − 1 or lower on the interval [ − 1,1]. This property allows us to assert that the DD constraints containing just a 'sine' function, as in (11) and (13), are made zero by this choice of function [33]. We can summarily state this as

$$\begin{equation} \underbrace{\int_{-1}^{1} \mathrm{d}u_{1} \cdots \int_{-1}^{u_{m-1}} \mathrm{d}u_{m}}_{\times m} \sin \left[\beta_{\mathrm{s}} (u_{k})\right] = 0 \end{equation}\noindent \tag{ 21 }$$

for m ⩽ N − 1 and 1 ⩽ k ⩽ m [33]. Furthermore, half of the corresponding 'cosine' constraints in the u-domain are zero by parity rules, since T_N+1(u) has definite parity. The remaining constraints are not zero, so this control function has a nonzero constraint at first- or second-order. Figure 1 compares $\sin \left [\beta _{\mathrm {s}} (u)\right ]$ to $\cos \left [\beta _{\mathrm {UDD}} (u)\right ]$, also known as the UDD 'switching function' y(t). This indicates that β_s preserves many of the decoupling properties of UDD, as supported by numerical simulations to be discussed later.

Zoom In Zoom Out Reset image size

Before proceeding to derive new sequences based on smooth approximations of UDD, we note that the connection of UDD's switching function to Chebyshev polynomials allows an alternative interpretation of its decoupling abilities. To arrive at a switching function, we treat the time-dependent operator B(t) as a random classical function, following standard procedures surrounding the Born–Markov approximation. In a particular trajectory, a qubit in state c₀|0〉 + c₁|1〉 evolves to $c_0|0\rangle +c_1\exp (\mathrm {i}\phi )|1\rangle$ with the accrued phase during the UDD sequence

$$\begin{equation} \phi = \int_0^T B(t) y(t)\,\mathrm{d}t. \end{equation} \tag{ 22 }$$

This integral may be rewritten

$$\begin{equation} \phi = \frac{2(N+1)}{\pi}\sum_{j=0}^{N+1} (-1)^j w_j \int_0^{\tau_j} B(\tau) \,\mathrm{d}\tau, \end{equation} \tag{ 23 }$$

where the τ_j are the locations of UDD pulses, i.e. τ_j = Tsin²(jπ/[2(N + 1)]), or, on the u axis, u_j = −cos(jπ/(N + 1). The weights are given by w_j = π/2(N + 1) for j = 0,N + 1 and w_j = π/(N + 1) otherwise. These weight functions are those of Gauss–Lobatto–Chebyshev quadrature, which says that for order M,

$$\begin{equation} \sum_{j=0}^{M} w_j f(u_j) = \int_{-1}^1 \frac{f(u)}{\sqrt{1-u^2}}\,\mathrm{d}u + R_M, \end{equation} \tag{ 24 }$$

where R_M is a remainder term which vanishes quickly with M and u_j = −cos(jπ/M). Hence we associate M, the order of Gauss–Lobatto–Chebyshev quadrature, with N + 1, where N is the number of π pulses in a UDD sequence. The other term we must consider is (−1)^j. This may be written as (−1)^j = T_N+1(u_j), evident from the formula T_n(u) = cos(ncos⁻¹(u)). Therefore, if we define

$$\begin{equation} f(u) = \frac{(N+1)T}{\pi}T_{N+1}(u)\int_{-1}^{u}B[T(v+1)/2]\,\mathrm{d}v, \end{equation} \tag{ 25 }$$

then (23) is exactly in the same form as the left-hand side of (24). We may therefore replace our sum with an integral with only higher-order corrections. Then, integration by parts and a return to the t domain yields

$$\begin{equation} \phi = \frac{4}{\pi}\int_{0}^T \mathrm{d}t B(t)U_N\left(\frac{2t}{T}-1\right)\sqrt{\frac{t}{T}\left(1-\frac{t}{T}\right)} + R_{N+1}. \end{equation} \tag{ 26 }$$

Recalling that $\sqrt {(t/T)(1-t/T)}=\sqrt {1-u^2}$ is the weight function for the orthogonal polynomials U_N(−u), we see that UDD effectively extracts the Nth term of an expansion of B(t) in Chebyshev polynomials of the second kind. Lower-order components of this expansion are decoupled. Higher-order components still appear in the remainder term R_N+1.

This result provides two major conclusions. First, it provides an alternative viewpoint as to how UDD achieves decoupling: UDD eliminates lower order Chebyshev polynomial components of the bath field. Second, direct comparison of (26) and (22) shows that, up to a remainder R_N+1 which vanishes for high N, the discontinuous switching function y(t) is well approximated by the continuous function $(4/\pi )U_N(2t/T-1)\sqrt {(t/T)(1-t/T)}=(4/\pi )\sin [\beta _{\mathrm {s}}(t)],$ as shown in figure 1(b). This connection motivates our effort to find smooth versions of the UDD sequence.

2.3. Variational parameters for improved continuous DD control functions

We will now develop a method for improving the DD properties of the continuous-control sequence given by β_s. We consider a technique for expanding the difference β_s − β_UDD as a Fourier series, where we subsequently truncate the series and tune the remaining terms for optimal performance. Recent work by Uhrig has adapted the UDD sequence to replace instantaneous pulses with time-dependent control functions, while preserving many of the error-suppression properties [28]. This followed a mathematical structure derived by Yang and Liu [33] used to prove that the original UDD sequence works as intended. Recall the variable substitution θ = cos⁻¹(−u). In the θ-domain, β_UDD(θ) is a 'staircase' function while β_s(θ) = Nθ is a linear approximation to β_UDD. The difference between the two sequences is a periodic sawtooth wave, which has a Fourier series given by

$$\begin{equation} \beta_{\mathrm{s}}(\theta) - \beta_{\mathrm{UDD}}(\theta) = \sum_{j = 1}^{\infty} \frac{\sin (2Nj \theta)}{j}. \end{equation} \tag{ 27 }$$

The first attempt at improving the control function would be to truncate the Fourier series with the first several terms (as shown in figure 2), but this results in only a modest improvement in decoupling. However, there is an insight to be gained here, because it turns out that any 'sine' Fourier components given by λ_psin(Npθ) do not affect the 'sine' DD constraints, in the sense that

$$\begin{equation} \underbrace{\int_{0}^{\pi} \mathrm{d}\theta_{1} \ldots \int_{0}^{\theta_{m-1}} \mathrm{d}\theta_{m}}_{\times m} \sin \left[\beta_{\mathrm{s}} (\theta_{k}) + \sum_{p} \lambda_{p}\sin (Np \theta_{k})\right] = 0 \end{equation} \tag{ 28 }$$

for 1 ⩽ k ⩽ m and any real coefficients λ_p (see appendix A). We can exploit this property to nullify some of the other DD constraints by varying the λ_p parameters without violating the above 'sine' DD constraints. In particular, the 'cosine' DD constraints do depend on the λ_p parameters, so we set up a system of nonlinear equations involving the first m constraints that are not already zero by identity and the first m variational parameters (λ_p). These equations involve nested integrals over complicated functions, and so finding genuine roots to the equations can be numerically intensive.

Zoom In Zoom Out Reset image size

The optimized set of parameters we examine more closely is given by five DD constraints that are not identically zero for the control function

$$\begin{equation} \beta_{\lambda}^{(5)} (\theta) = 4\theta + \sum_{k=1}^{5} \lambda_{k}\sin (4k \theta), \end{equation}\noindent \tag{ 29 }$$

where λ₁,...,λ₅ are variational parameters we can change to satisfy the five DD constraints. As discussed previously, this functional form automatically nullifies 'sine' DD constraints as in (11) and (13). Furthermore, this form has definite parity (even), so half of the 'cosine' constraints can be eliminated by parity rules. The system of equations we solve for is given by:

$$\begin{equation} \int_{-1}^{1} \cos [\beta (u_{1})] \,\mathrm{d}u_{1} = 0, \end{equation} \tag{ 30 }$$

$$\begin{equation} \int_{-1}^{1} \mathrm{d}u_{1} \int_{-1}^{u_{1}} \mathrm{d}u_{2} \sin [\beta (u_{1}) - \beta (u_{2})] = 0, \end{equation} \tag{ 31 }$$

$$\begin{equation} \int_{-1}^{1} \left(u_{1}\right)^2 \cos [\beta (u_{1})]\, \mathrm{d}u_{1} = 0, \end{equation} \tag{ 32 }$$

$$\begin{equation} \int_{-1}^{1} \mathrm{d}u_{1} \int_{-1}^{u_{1}} \mathrm{d}u_{2} \left(u_{2}\right)^{2} \sin [\beta (u_{1}) - \beta (u_{2})] = 0, \end{equation} \tag{ 33 }$$

$$\begin{equation} \int_{-1}^{1} \left(u_{1}\right)^4 \cos [\beta (u_{1})] \,\mathrm{d}u_{1} = 0. \end{equation}\noindent \tag{ 34 }$$

Importantly, the parity rules noted in section 2.2 only apply to equations in the u-domain, not in real time. Still, this set of equations nullifies all of the first- and second-order constraints in (11)–(15). Moreover, several error terms at third-, fourth- and fifth-order are also canceled. These constraints are chosen for being the error terms up to fifth order in T which: (i) are not automatically zero, and (ii) require single- or double-integration. Higher-order integrals posed problems for the nonlinear numerical solver we used; in principle, further work in developing robust programs to solve nonlinear systems of equations particular to this decoupling problem can yield control sequences which decouple to higher orders than third. Alternative formalisms based on time-dependent perturbation theory, as considered for example in [34], might also be employed for higher-order sequences.

The λ_p parameters were treated as variables to solve the system of equations in (30)–(34). The initial approximate solution was taken to be the truncated Fourier series, and table 1 lists the initial parameters and their values after optimization. Figure 2 shows the optimized continuous control function with five λ_p parameters. We call such decoupling sequences with optimized parameters CAFE. To specify precisely which CAFE sequence is being considered, we denote each as CAFE^(N,m) where N is the order of the continuous approximation to UDD, as appears in (17) and subsequent references to β_s like (29). The quantity m is the number of variational λ_p parameters. For example, figure 2 shows the control function α⁽⁵⁾(t) for CAFE^(3,5).

Table 1. The first five λ_p parameters for a truncated Fourier series and the optimized sequence found by numerical solution of five decoupling constraints.

	Fourier	Optimized search
λ1	0	0.0017
λ2	1	0.9121
λ3	0	−0.2869
λ4	1/2	1.3520
λ5	0	0.4920

2.4. Piecing together sequences: splice and invert

The previous section described how to improve sequences by inserting m variational parameters and searching for a set that minimizes m previously nonzero DD constraints. The resulting control function α^(m)(t) has the undesirable property of having poles at the start and endpoints, which are impractical to implement in an experiment. However, the function is very steep here as well, so it seems reasonable to 'chop off' this section of the sequence with minimal impact on the decoupling properties. This process is illustrated in figure 3(a). The problem here is that one must also ensure that the net action of the control function is the identity gate. Otherwise, the sequence will decouple the qubit well, but the state of the qubit will be rotated deterministically around the Bloch sphere. This implies that the pulse area in β(t) that is removed must be equal to 2π. The pulse area can be calculated readily, but ensuring this area is precisely 2π in a real experiment seems unreasonable, given how large the first derivative of α(t) is at the endpoints of the sequence.

Zoom In Zoom Out Reset image size

A rather simple and effective solution is to chop at zeros of α(t) near the start/end of the sequence, as shown in figure 3(a). For α(t) being an even function of time, the zeros lie symmetrically about the sequence, so we can splice together a window of time with control α(t) followed by an equal window of time with control −α(t) which has been chopped in the same manner; for α(t) being an odd function of time, one must time-reverse control in one of the windows. Figure 3(b) shows how to use this method to make a simple bounded sequence by splicing and inverting two segments of the CAFE sequence from figure 2. The resulting sequence has several advantages: (i) it is finite and continuous, and the first derivative is continuous, at all points; (ii) the total pulse area is automatically zero since the second half of the sequence reverses any action in the first half; (iii) this sequence preserves much of its DD properties and performs quite well in simulation. We expand the notation for this decoupling function to CAFE^(N,m,r) × L where r is the root-number at which the sequence is truncated (counting from the end) and L is the number of times the CAFE base sequence is repeated. For example, figure 3(b) shows a CAFE^(3,5,2) × 2 sequence. The advantages of symmetrizing or antisymmetrizing sequences using similar methods are well established in decoupling; for an early treatment from NMR see [35] and for a more recent discussion, see [21].

We evaluate the performance of the CAFE sequences in comparison to a very similar set of sequences, namely the RUDD sequences of Uhrig and Pasini [27–29]. However, we make two alterations to the original RUDD sequence to make the comparison as appropriate as possible. First, the RUDD sequence as originally derived [27] begins and ends with very sharp pulses of area 2π. Our numeric studies as well as those of Uhrig and Pasini [28, 29] have shown that these pulses have very little effect on the decoupling abilities of the RUDD sequence. We therefore omit these pulses in all comparisons shown in this paper. Another modification we make is to invert the sign of each subsequent pulse in the RUDD sequence. We denote a RUDD sequence modified in this way by RUDD*-N-d, where N is the number of pulses and d is the duty cycle.

The inversion of subsequent copies of CAFE was introduced to ensure continuity, but this process also introduced a further advantage of adding robustness to low-frequency error in the amplitude of the control field. Such errors will be discussed in more detail in section 6. In the case of RUDD with perfect pulses, the modification of inverting each pulse has no effect on the performance of the sequence. However, it introduces the same robustness to low-frequency control noise, and so therefore provides a closer comparison to the CAFE sequences.

In this paper, we will limit our discussion to CAFE^(3,5,r) × L sequences, varying only the root number r and the number of repetitions L. We will compare these sequences to RUDD*-N-d sequences of different pulse numbers N and duty cycles d. In general, increasing L or reducing r improves the decoupling ability of CAFE, while increasing N or reducing d improves the decoupling ability of RUDD*. Any of these actions, however, increases the maximum slew rate, $\max _t |\mathrm {d}\alpha /\mathrm {d}t|$, of the sequence. For RUDD* a high duty cycle d means that pulses are stretched to have longer duration, which reduces slew rate since integrated pulse area is constant (π-pulses). Limits on slew rate are likely to be the primary reason for using continuous sequences such as CAFE and RUDD*, and therefore we use this criterion to choose which CAFE^(3,5,r) × L sequence to compare to which RUDD*-N-d sequence. In particular, for numeric studies we will look at two groups of sequences, each of roughly equal maximum slew rate. One group has a lower slew rate and is summarized in figure 4; a second with higher slew rate is summarized in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our first evaluation of the CAFE sequence considers a fully quantum mechanical description of a bath. For this we employ the canonical central-spin problem, in which our qubit is coupled to six spins which are themselves coupled via dipole–dipole interactions, i.e.

$$\begin{equation} B_0 = 4J \sum_{j \ne k} \frac{r_{jk}^2\textbf{I}_j\cdot\textbf{I}_k-3(\textbf{r}_{jk}\cdot \textbf{I}_j)(\textbf{r}_{jk}\cdot\textbf{I}_k)}{r_{jk}^5}, \end{equation} \tag{ 35 }$$

$$\begin{equation} B_Z = \sum_j C_j I^z_j, \end{equation} \tag{ 36 }$$

where $\textbf{I}_j$ is the spin vector (using I = 1/2) for the bath spins and $\textbf{r}_{jk}$ is the spatial vector connecting them. The cubic geometry we choose is shown in figure 6(a). This type of simulation tests the DD sequences against a system with more degrees of freedom (see (11)–(15)) than the Born–Markov classical bath simulations we present later. However, these additional degrees of freedom also complicate analysis, so we present no analytic expectation for the result. We demonstrate only that CAFE successfully decouples a qubit from a fully-quantum environment in a manner highly comparable to RUDD*.

Zoom In Zoom Out Reset image size

The plots in figure 6(b) show the projections of the qubit Bloch vector in the interaction picture after averaging over the bath spins. Specifically, this simulation initializes into state

$$\begin{equation} \rho_j(0)=\frac{\sigma^I+\sigma^j}{2}\otimes \frac{\mathbbm{1}}{2^6}, \end{equation} \tag{ 37 }$$

where σ^j is either a Pauli matrix or single-qubit identity (j = I,x,y,z) and $\mathbbm {1}$ is the bath identity matrix. Numeric integration then solves dρ_j/dt = −i[H,ρ_j] for each j = x,y,z, and the three qubit directions k = x,y,z are plotted as Tr[σ^kρ_j(t)] at each time t for varying initial conditions j.

To evaluate infidelity, we consider the procedure of quantum process tomography (QPT), in which the final state of the qubit, after the sequence is complete at time T, may be written

$$\begin{equation} \rho(T) = \sum_{j,k=\{I,x,y,z\}} \chi_{jk}\sigma^j \rho(0) \sigma^k. \end{equation} \tag{ 38 }$$

In single-qubit QPT, measurements of the three components of the Bloch vector with four initial conditions—ρ_x, ρ_y, ρ_z, and ρ_−z—are sufficient to construct the matrix χ_jk in general [36], but since this system features unitary evolution prior to the final trace, projections for initial condition ρ_−z are exactly the negative of projections for initial condition ρ_z. Further, for our particular interaction and bath, the Ising-like symmetry assures Tr[σ^yρ_x(t)] = Tr[σ^zρ_x(t)] = 0 and Tr[σ^xρ_j(t)] = 0 for all j ≠ x. Hence the four choices of projection and initial condition shown in figure 6(b) are sufficient to construct χ_jk. The infidelity for DD is defined as 1 − χ_II; which may be interpreted as the probability that something other than identity happened to the qubit. In general, for a single qubit in which the process is unitary evolution involving a bath followed by tracing that bath, one may readily show

$$\begin{equation} \mathcal{I} = 1-\chi_{\mathrm{II}} = \frac{1}{4}\left\{3-\sum_{j=\{x,y,z\}}\mathrm{Tr}[\sigma^j\rho_j(T)]\right\}. \end{equation}\noindent \tag{ 39 }$$

Figure 7 shows the infidelity as a function of J, the overall strength of the dipole–dipole coupling within the bath, for both the low-slew-rate sequences and high-slew-rate CAFE and RUDD* sequences shown in figures 4 and 5. Varying J results in varying the correlation time of the bath. At low J, the bath is very slow, and here dephasing without DD is maximal and determined entirely by the number of possible configurations of the bath spins. As J increases, even without DD the infidelity decreases due to motional narrowing, in which spin-flips within the bath reduce its coupling to the qubit. With CAFE and RUDD* sequences, it is clearly seen that DD strongly reduces the infidelity at low J, corresponding to long correlation times. The efficacy of decoupling reduces, at varying rates depending on the sequence, as J is increased; improved fidelity at even higher J is again due to motional narrowing effects.

Zoom In Zoom Out Reset image size

Figure 6 shows that RUDD* sequences and CAFE sequences are qualitatively different in how they decouple in time; in particular RUDD remains pulsed, while CAFE is fully continuous. However, figure 7 shows that they have roughly the same decoupling power per slew rate as JT is varied. To analyze their response to varying bath correlation times more quantitatively, we shift to a classical, Born–Markov type bath and a filter function analysis.

When studying the effects of noise on information, often it is instructive to analyze the system in the frequency domain. DD schemes are effectively high-pass filters, meaning they suppress low-frequency components of the system–environment coupling while not reducing high-frequency noise. The 'cut-off frequency' between these regions is often closely related to the rate at which one can apply control pulses, but the precise relationship is less clear in examples like UDD, where pulses are not equidistant. To accurately describe how a decoupling sequence protects against noise with a given spectral density, we characterize each sequence by a 'filter function' [25, 26, 37, 38], which characterizes the degree to which noise at a given frequency is suppressed. Filter functions are a meaningful way to compare different DD sequences because they provide insight into which method is best suited to a given noise spectrum.

We require that the filter function be a linear response which is independent of the magnitude of the environment noise, but this is only appropriate under certain conditions. First, we invoke the Born–Markov approximation and represent the dephasing bath with a scalar function B(t), so that the Hamiltonian is now

$$\begin{equation} \mathcal{H}(t) = B(t) \sigma^z + \alpha(t) \sigma^x, \end{equation} \noindent \tag{ 40 }$$

where B(t) is a stationary random noise field with zero mean and spectral density S(ω) given by

$$\begin{equation} S(\omega) = \int_{-\infty}^{\infty} \mathrm{d}t \langle B(t) B(\tau) \rangle \cos(\omega (t-\tau)). \end{equation} \noindent \tag{ 41 }$$

Second, we assume that the magnitude of the bath noise is sufficiently weak that we can calculate the qubit infidelity as a linear function of S(ω). This assumption is not valid in general; it requires either a weak-noise limit made explicit in appendix B, or more restrictive assumptions such as the instantaneous π-pulses used in [25, 26, 37]. Relaxing this assumption requires higher-order terms in the filter function, which have been considered in [31]. Under our assumption, we may define the filter function F(z) by stating that the qubit infidelity $\mathcal {I} = 1 - \chi _{\mathrm {II}}$, or total error, is given by

$$\begin{equation} \mathcal{I} = \int_0^{\infty} \frac{\mathrm{d}\omega}{2 \pi \omega^2} S(\omega) F(\omega T). \end{equation} \noindent \tag{ 42 }$$

The filter function for a given decoupling sequence can be calculated numerically in simulation, by modeling the noisy environment with a sinusoid. By measuring the degradation in qubit fidelity averaged over an ensemble of random phases for the sinusoid at frequency ω, one obtains the linear response F(ωT). For control fields without a σ^z-component, the filter function may be theoretically expressed as

$$\begin{equation} F(\omega T) = \left|\omega T\int_0^1 \mathrm{e}^{\mathrm{i}[\omega Tu+\beta(Tu)]}\,\mathrm{d}u \right|^2. \end{equation} \tag{ 43 }$$

For a general derivation of this expression in the continuous case, see appendix B. This expression clearly reduces to the canonical expression for the filter function derived for sequences composed of ideal π-pulses such as UDD [25, 26, 37, 38]. For such sequences, e^iβ(t) = cos[β(t)] = y(t) is the switching function. Note that if the smooth approximation for y(t) corresponding to UDD, i.e. $(4/\pi ) U_N(2t/T-1)\sqrt {(t/T)(1-t/T)}$ as discussed in section 2.2, is substituted into (43), the integral may be analytically performed to yield F(ωT) ≈ 16(N + 1)²J²_N+1(ωT/2) for UDD-N, where J_N(z) is the Bessel function of order N. This expression yields a very close approximation to the exact UDD filter function, as stated without derivation in [25].

Figure 8 shows filter functions comparing different versions of CAFE and RUDD, calculated both with explicit integration and via Monte Carlo simulation. Examining figure 8, we see that despite dissimilar time-domain construction, CAFE provides an intermediate limit between low-N, low-duty-cycle RUDD sequences and high-N, high-duty-cycle RUDD sequences. The sharp dips in the filter function are quite different between the two cycles, but these play little role in realistic noise bath models which typically integrate out such features. We turn to such models next.

Zoom In Zoom Out Reset image size

To further analyze the effectiveness of CAFE sequences in comparison to RUDD* we simulate the total evolution of arbitrary qubit states coupled to classical dephasing baths. As in the filter function analysis, we invoke the Born–Markov approximation leading to the scalar noise field B(t) in (40) with noise spectral density S(ω) given by (41). We now consider bath noise with both Gaussian and Lorentzian spectral density. The reason we consider a classical bath is that by doing so, we can readily define parameters such as absolute magnitude and correlation time for a scalar noise spectral density, as opposed to the quantum bath which has far too many degrees of freedom to permit simplified analysis.

Our simulations incorporate two sources of noise: environmental dephasing noise and noisy control fields. Control noise is inevitable in a practical experimental setting. Even a sophisticated decoupling sequence could yield mediocre results if such control noise is not taken into account; moreover, a DD sequence could possibly introduce more noise than it corrects. Due to these possibilities, we investigate the effect of control noise on CAFE and RUDD* by simulating their decoupling effectiveness with and without noise added to the control functions.

Figure 9 shows the behavior of our example CAFE and RUDD* sequences under Born–Markov noise with a spectral density of Gaussian form, i.e. $S(\omega )\propto \exp [-(\omega \tau_{\rm c})^2]$, both with perfect control and with 1% control noise. The control noise is modeled as a time-dependent fluctuation on the amplitude of the control field α(t), taken with root-mean-square magnitude 1% of α(t) and a Gaussian spectral density with correlation time T/2 for the plots shown. Theoretical curves, found by numeric integration, are shown for comparison as solid lines. No theoretical curves are calculated for the case of nonzero control noise; the zero-control noise curves are repeated for reference. The approximations inherent in these curves fail most prominently at very low correlation times, where higher-order noise terms neglected in the filter-function analysis become important, and at very long correlation times, where a noise floor results from numeric inaccuracy. That numeric noise floor may be explicitly raised by deliberately adding control noise.

Zoom In Zoom Out Reset image size

Figure 10 shows a similar simulation with the spectral density of Lorentzian form, i.e. S(ω)∝1/[1 + (ωτ_c)²]. The soft roll-off of this spectral density causes stronger variation from the filter-function analysis, since higher-frequency components poorly modeled by the lowest-order filter function analysis are more important. This discrepancy highlights the conclusion that for baths with soft frequency cut-offs and continuous sequences, low-order filter functions are of limited use, and numeric simulation is an important tool for analysis.

Zoom In Zoom Out Reset image size

To compare sequences for this Born–Markov bath case, in the case of no control noise, CAFE's performance is roughly intermediate between the RUDD* sequences, as already evident from the filter functions analysis of figure 8. One should note the strong similarity between the results of figure 9 with no control noise and those of figure 7; the nearly identical behavior gives a clear example of a fully quantum bath behaving identically to a classical one, supporting the frequent use of the Born–Markov approximation to describe spin baths.

In the cases where control noise is added, we see that all sequences roughly follow their expected behavior at short correlation times, where decoupling is already least effective. In the regime of long correlation times, where DD performs best, a noise-floor is reached which depends on the sequences; this noise-floor is consistently better for the RUDD* sequences. It is tempting to associate this lower noise floor with a lower duty cycle for the RUDD* sequences with respect to CAFE's 100%, but this conclusion is questioned by the highly similar behavior of RUDD* sequences of very different duty cycle. As discussed in section 3, the noise floor of the RUDD sequences without the modification of alternating the sign of pulses, not shown, is much worse. This simple method of alternating pulse signs has certainly added some robustness, but many improvements in optimization against control noise are certainly possible in future work.

The CAFE sequence originated from the analysis of UDD as if it were continuous. Using principles from numerical quadrature, we showed that UDD's essential behavior could be captured by continuous approximations to its normally discontinuous description. The continuous approximation to UDD provided a heuristic for continuous sequence design. The RUDD sequences follow a similar path, but are based on a different heuristic [27–29]. Both sequence families can be extended to higher order. A promising feature of these types of sequences appears to be that, for these two families at least, sequences can be locally altered to eliminate short-time spiking behavior without appreciable detriment to their decoupling ability. Although different heuristics were employed, the two sequence families behave remarkably similarly, as we have demonstrated using a variety of analyses.

The reason analysis of continuous-control decoupling is more complicated than the analysis of 'bang–bang' sequences is because noise can disturb the qubit during control operations. More than a complication of analysis, however, this situation is a reality of experiment. We believe, based on our results, that filter functions remain relevant for continuous sequences, and they provide an intuitive picture for what components of the environment noise are suppressed. Furthermore, Monte Carlo simulations of semi-classical noise bath spectra and fully-quantum spin lattices indicate how well a sequence performs at protecting quantum information. By adding noise to the control field, we can test another source of error, the experimenter and the experimenter's equipment, and how well the sequence copes. We hope the heuristics underlying the construction of CAFE and the characterization methods we have used may be useful to future research in continuous DD.

Sponsored by United States Department of Defense. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressly or implied, of the United States Department of Defense or the US Government. Approved for public release, distribution unlimited.

We begin with the definitions from before that

$$\begin{equation} \beta_{\mathrm{s}} (t) = (N+1) \cos ^{-1} (-u), \end{equation}\noindent \tag{ A.1 }$$

which in θ-domain (θ = cos⁻¹(−u)) becomes

$$\begin{equation} \beta_{\mathrm{s}} (\theta) = (N+1) \theta. \end{equation} \tag{ A.2 }$$

We can now define a control function with m λ-parameters given as

$$\begin{equation} \beta_{\lambda}^{(m)} (\theta) = \beta_{\mathrm{s}} (\theta) + \sum_{k} \lambda_{k}\sin [(N+1)k \theta]. \end{equation} \noindent \tag{ A.3 }$$

β^(m)_λ is a control sequence with m variational parameters designed to cancel m DD constraints not zero for β_s alone. Here we show that this construction will not affect the 'sine'-like DD constraints given by

$$\begin{equation} \underbrace{\int_{0}^{T} \mathrm{d}t_{1} \cdots \int_{0}^{t_{m-1}} \mathrm{d}t_{m}}_{\times m} \sin [\beta (t_{k})] = 0 \end{equation} \noindent \tag{ A.4 }$$

for 1 ⩽ k ⩽ m. Using integration by parts, the above set of constraints can be reduced to the system of equations

$$\begin{equation} \int_{0}^{T} t^{n} \sin (\beta) \mathrm{d}t = 0 \end{equation} \noindent \tag{ A.5 }$$

for n = 0,1,...,(m − 1). Working in the θ-domain, this set of constraints is equivalent to

$$\begin{equation} \int_{0}^{\pi} [1-\cos (\theta)]^{n} \sin (\beta) \sin (\theta)\,\mathrm{d}\theta = 0, \end{equation} \noindent \tag{ A.6 }$$

which is also equivalent to

$$\begin{equation} \int_{0}^{\pi} \cos (p\theta) \sin (\beta) \sin (\theta) \,\mathrm{d}\theta = 0 \end{equation} \noindent \tag{ A.7 }$$

for p = 0,1,...,(m − 1). By substituting the expression in (A.3),

$$\begin{eqnarray} &&\fl \int_{0}^{\pi} \cos (p\theta) \sin [\beta_{\lambda}^{(m)}] \sin (\theta)\, \mathrm{d}\theta\nonumber\\&&\ms =\int_{0}^{\pi} \cos (p\theta) \sin (\theta) \sin [(N+1)\theta] \cos \biggl[\sum_{k} \lambda_{k}\sin ((N+1)k \theta)\biggr] \mathrm{d}\theta\nonumber\\&&\ms \quad+\int_{0}^{\pi} \cos (p\theta) \sin (\theta) \cos [(N+1)\theta] \sin \biggl[\sum_{k} \lambda_{k}\sin ((N+1)k \theta)\biggr] \mathrm{d}\theta. \end{eqnarray} \noindent \tag{ A.8 }$$

Note that $\cos [\sum _{k} \lambda _{k}\sin ((N+1)k \theta )]$ can be written as the infinite sum $\sum _{j = 0}^{\infty } a_{j} \cos ((N+1)j\theta )$ for some real coefficients a_j since cos(ϕ) is an even function in ϕ. Similarly, $\sin [\sum _{k} \lambda _{k}\sin ((N+1)k \theta )]$ can be written as $\sum _{j = 1}^{\infty } b_{j} \sin ((N+1)j\theta )$ for some real coefficients b_j since sin(ϕ) is an odd function in ϕ. The DD constraints are therefore satisfied if

$$\begin{equation} \int_{0}^{\pi} \cos(p\theta) \sin(\theta) \sin[A(N+1)\theta] \cos[B(N+1)\theta]\mathrm{ d}\theta = 0 \end{equation} \noindent \tag{ A.9 }$$

for integers A, B satisfying A > 0, B ⩾ 0, and p < N. This equation can be directly evaluated to always be true under these conditions.

In this appendix we derive a general expression for the filter function of a continuous DD sequence. We generalize slightly from the main text, allowing a Hamiltonian of the form

$$\begin{equation} H(t) = 2\epsilon B(t)S^z + \textbf{A}(t)\cdot\textbf{S}. \end{equation} \tag{ B.1 }$$

Here $\textbf{A}$ is a vector of control fields, one for each component, and $\textbf{S}$ is the spin vector, i.e. S^j = σ^j/2. The unitless parameter is in place to keep track of orders of bath interactions, since the filter function is defined only at the lowest nonvanishing order. We use to define that the filter function F(z) in terms of the infidelity $\mathcal {I}$ at vanishing . Again, $\mathcal {I}$ is defined as 1 − χ_II for the QPT-derived matrix χ, due to evolution at finite . Then F(z) is defined so that

$$\begin{equation} \lim_{\epsilon\rightarrow 0} \epsilon^{-2} \mathcal{I} = \int \frac{\mathrm{d}\omega}{2\pi\omega^2} S(\omega) F(\omega T). \end{equation} \tag{ B.2 }$$

This type of -based definition is not required for pulses composed of δ-function pulses, but in the continuous case where noise occurs during the control, the limit as  → 0 is critical for the linear-response theory upon which the filter function is based.

We begin by transforming to an interaction picture, for which we need the control unitary $U_0(t)=\mathcal {T}\exp [-\mathrm {i}\int _0^t \textbf{A}(\tau )\cdot \textbf{S}]\mathrm {d}\tau$. In terms of this unitary, our interaction picture Hamiltonian is $\tilde {H}(t)=U_0^\dagger (t) H(t) U_0^{\phantom {\dagger }}(t)$. If $\textbf{A}(t)$ changes in direction over time, U₀(t) in general may describe complicated trajectories over the Bloch sphere. In this most general case, it is not clear how to arrive at an analytic form for a filter function, and one generally must rely on numeric solvers. Here, we consider the simple case that $\textbf{A}(t) = 2\alpha (t)\hat {\textbf{n}}$, for constant direction described by unit vector $\hat {\textbf{n}}$. The CAFE and RUDD sequences considered in this paper have all used $\hat {\textbf{n}}$ along the x direction; our definition of infidelity and assumption about the bath interaction indicate that any direction along the Bloch sphere equator is equivalent. In this appendix, we generalize slightly by allowing a constant z-component to $\hat {\textbf{n}}$ as well. Then,

$$\begin{equation} \fl \tilde{H}(t) = 2\epsilon B(t) \{S^z\cos[\beta(t)]+\hat{\textbf{n}}^z\hat{\textbf{n}}\cdot\textbf{S}\{1-\cos[\beta(t)]\} -\textbf{S}\cdot \hat{\textbf{z}}\times\hat{\textbf{n}}\sin[\beta(t)]\}, \end{equation} \tag{ B.3 }$$

using β(t) as defined by (2).

We consider a single-qubit density matrix in the interaction picture

$$\begin{equation} \rho(t) = \frac{1}{2} + \sum_{j}\rho^j(t) S^j, \end{equation} \tag{ B.4 }$$

which evolves via $\tilde {H}(t)$, to lowest nonvanishing order in time-dependent perturbation theory with ensemble averaging and standard Born–Markov separation of timescales, as

$$\begin{equation} \frac{\mathrm{d}\rho}{\mathrm{d}t} = -\int_0^t \mathrm{d}\tau \langle [\tilde{H}(t),[\tilde{H}(\tau),\rho(\tau)]]\rangle. \end{equation} \tag{ B.5 }$$

Correspondingly each Bloch-vector component ρ^j evolves according to the following, where each repeated superscript is summed:

$$\begin{eqnarray} \frac{\mathrm{d}\rho^j}{\mathrm{d}t} &=& -2\int_0^t \mathrm{d}\tau\; \mathrm{Tr}\{S^j \langle [\tilde{H}(t),[\tilde{H}(\tau),S^k]]\rangle \}\rho^k(\tau)\nonumber\\ &=& \epsilon^2\int_0^t \mathrm{d}\tau \int\frac{\mathrm{d}\omega}{2\pi} \cos[\omega(t-\tau)] S(\omega) M^{k\ell}(\tau)M^{\ell j}(t)\rho^k(\tau)\\ &=& -\epsilon^2\int_0^t \mathrm{d}\tau \int\frac{\mathrm{d}\omega}{2\pi} \cos[\omega(t-\tau)] S(\omega) \Gamma^{jk}(t,\tau)\rho^k(\tau),\nonumber\ \end{eqnarray} \tag{ B.6 }$$

where

$$\begin{eqnarray} &&\fl M^{k \ell}(t)=2\{\varepsilon^{zk\ell}\cos[\beta(t)] +\hat{\textbf{n}}^z\hat{\textbf{n}}^m\epsilon^{mk\ell}\{1-\cos[\beta(t)]\} +[\hat{\textbf{n}}^k\delta^{z\ell}-\delta^{kz}\hat{\textbf{n}}^{\ell}]\sin[\beta(t)]\}. \end{eqnarray}\noindent \tag{ B.7 }$$

Here ε^jkℓ is the fully antisymmetric Levi–Cevita tensor density and δ^jk is the identity matrix.

Now, we time-integrate ρ^j to lowest order in , to find

$$\begin{equation} \fl \rho^j(T) =\left[\mathbbm{1}-\epsilon^2\int_0^T \mathrm{d}t \int_0^t \mathrm{d}\tau \int \frac{\mathrm{d}\omega}{2\pi} \cos[\omega(t-\tau)]S(\omega)\Gamma^{jk}(t,\tau)\right]\rho^k(0)+O(\epsilon^4). \end{equation}\noindent \tag{ B.8 }$$

We use this expression in QPT, which in general gives

$$\begin{eqnarray} &&\fl \chi_{\mathrm{II}}(T) = \textstyle\frac{1}{8}\left[2+2\rho_x^x(T)+2\rho_y^y(T)+\rho_z^z(T) -\rho_{-{z}}^z(T) -\rho_z^x(T)-\rho_z^y(T)-\rho_{-{z}}^x(T)-\rho_{-{z}}^y(T)\right],\nonumber\\ \end{eqnarray} \tag{ B.9 }$$

where ρ^k_j = Tr{σ^kρ_j(T)} for initial condition ρ_±j(0) = (σ^I ± σ^j)/2. As discussed in section 4, the unitary evolution of this model assures ρ^j_z = −ρ^j_−z, simplifying the fidelity expression to

$$\begin{equation} \chi_{\mathrm{II}}(T) = \frac{1+\rho_x^x(T)+\rho_y^y(T)+\rho_z^z(T)}{4}. \end{equation} \tag{ B.10 }$$

As a result, the infidelity is simply written

$$\begin{equation} \mathcal{I} = \frac{\epsilon^2}{4}\int_0^T \mathrm{d}t \int_0^t \mathrm{d}\tau \int \frac{\mathrm{d}\omega}{2\pi} \cos[\omega(t-\tau)]S(\omega)\mathrm{Tr}\{{\bf{\Gamma}}(t,\tau)\}+O(\epsilon^4), \end{equation} \tag{ B.11 }$$

and our generalized filter function is therefore simply

$$\begin{equation} F(z) = \frac{z^2}{4}\int_0^1 \mathrm{d}u \int_0^u \mathrm{d}v \cos[z(u-v)]\mathrm{Tr}\{{\bf{\Gamma}}(Tu,Tv)\}. \end{equation} \tag{ B.12 }$$

Upon taking the trace, the matrix Γ^jk(t,τ) simplifies substantially to

$$\begin{equation} \mathrm{Tr}[{\bf{\Gamma}}(t,\tau)] =8\{\cos[\beta(t)-\beta(\tau)] +[\hat{\textbf{n}}^z]^2[1-\cos[\beta(t)-\beta(\tau)]]\}. \end{equation} \tag{ B.13 }$$

In this paper, we only consider control fields for which $\hat {\textbf{n}}^z(t)=0$ at all times. In this case our general expression for F(z) is

$$\begin{equation} F(z) = 2z^2\int_0^1 \mathrm{d}u \int_0^u \mathrm{d}v \cos[z(u-v)]\cos[\beta(Tu)-\beta(Tv)]. \end{equation} \tag{ B.14 }$$

This easily simplifies to (43) in section 5.