Arbitrary quantum control of qubits in the presence of universal noise

Broadly speaking, a qubit is a quantum system that exhibits two possible outcomes for the measurement of a particular physical observable: the z-component of the intrinsic angular momentum of a spin-1/2 particle being the archetypal example. Letting |0〉 and |1〉 denote the states in which one or other of these outcomes are returned with certainty, and assuming maximal knowledge of the system, an arbitrary qubit state may be written as a linear combination

$$\begin{equation} |\psi\rangle=\cos(\vartheta/2)|0\rangle+\mathrm{e}^{\mathrm{i}\varphi}\sin(\vartheta/2)|1\rangle. \end{equation} \tag{ 1 }$$

The angles (ϑ,φ) define a three-dimensional unit vector, the Bloch vector, that provides a convenient pictorial representation of the state [10]. Interpreting all possible states in this way, the two-dimensional qubit state space is mapped to a sphere of unit radius, called the Bloch sphere (figure 1).

Zoom In Zoom Out Reset image size

The temporal evolution of the qubit is described by the propagator U(t,0), a unitary operator that transforms some initial state |ψ(0)〉 to a final state |ψ(t)〉 = U(t,0)|ψ(0)〉, for t ⩾ 0. Within the Bloch sphere picture, this 'length-preserving' transformation corresponds to a simple rotation of the qubit Bloch vector. The propagator satisfies the Schrodinger equation (ℏ = 1)

$$\begin{equation} \mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}U(t,0)=H(t)U(t,0), \end{equation} \tag{ 2 }$$

where H(t) is the Hamiltonian operator, which we may assume to be traceless (i.e. ${\textrm {Tr}}{\left (H(t)\right )}=0$), for all t. Disregarding irrelevant global phase factors, the set of all single qubit propagators, in conjunction with the usual operator product, forms the Lie group SU(2) [36].

That a geometrical interpretation of the action of a propagator is possible in terms of a rotation of the Bloch vector is a consequence of the homomorphism (structure-preserving map) that may be constructed between SU(2) and the group of three-dimensional rotation operators, SO(3) [36]. In short, the action of each propagator U∈SU(2) (and its physically equivalent negative −U) on a single-qubit state |ψ〉 can be represented by the action of a rotation operator R∈SO(3) on the Bloch vector corresponding to |ψ〉. To make this relationship explicit, we note that an arbitrary solution to (2) may be written in the form $U(t,0)=\exp {[-\mathrm {i}\theta (t){\hat{\boldsymbol{n}}}(t){\boldsymbol {\sigma }}/2]}$, where ${\hat{\boldsymbol{n}}}(t)=(n_{x}(t),n_{y}(t),n_{z}(t))$ is a real row vector of unit length, σ = (σ_x,σ_y,σ_z)^T is a column vector comprised of the Pauli spin operators, and θ(t) is a real function of time. The rotation operator in SO(3) corresponding to U(t,0) is that which rotates the Bloch vector through an angle θ(t), about an axis defined by ${\hat{\boldsymbol{n}}}(t)$. Thus, any single-qubit propagator may be thought of as acting to cause a simple rotation of the Bloch vector.

To model the effect of noise on the capacity to control the evolution of a simple quantum system, we consider an ensemble of identically prepared noninteracting qubits evolving according to a Hamiltonian

$$\begin{equation} H(t)=H_{0}(t)+H_{\mathrm{c}}(t). \end{equation} \tag{ 3 }$$

Here, the effect of the environment is modeled by the generalized noise Hamiltonian

$$\begin{equation} H_{0}(t)=\boldsymbol{\beta}(t){\boldsymbol{\sigma}}, \end{equation} \tag{ 4 }$$

where $\boldsymbol {\beta }(t)=(\beta _{x}(t), \beta _{y}(t), \beta _{z}(t))$ is a three-element row vector, each component of which is a random process modeling classical noise in one of the three spatial dimensions. Control over the state of the qubit can be achieved through the application of an external field, represented by a control Hamiltonian

$$\begin{equation} H_{\mathrm{c}}(t)=\boldsymbol{h}(t){\boldsymbol{\sigma}}, \end{equation} \tag{ 5 }$$

where the components of the vector h(t) = (h_x(t),h_y(t),h_z(t)) describe the strength and direction of the control field as a function of time.

The presence of the noise term in (3) will, in general, affect adversely the ability to accurately 'steer' the qubit state via the control Hamiltonian H_c(t). As a simple illustration, consider the application of a constant control Hamiltonian H_c(t) = Ωσ_x/2 to an otherwise isolated qubit, where Ω ≡ π/τ for some positive real number τ. Under these circumstances, the total Hamiltonian is independent of time and equation (2) is easily integrated to give $U(\tau ,0)=\exp {[-\mathrm {i}\pi \sigma _{x}/2]}$. Thus, over a time interval [0,τ] the control field generates a rotation of the qubit Bloch vector through an angle of π radians, about the x-axis (a 'π_X-pulse'), as shown in figure 1(a). If we allow for nonzero but constant 'noise' by letting β(t) = (0,0,β_z), for example, then $U(\tau ,0)=\exp {[-\mathrm {i}\pi (\sigma _{x}+\nu \sigma _{z})/2]}$, where ν ≡ 2β_z/Ω. The propagator now describes a rotation through an angle θ' > π, about an axis that is tilted away from the xy-plane (figure 1(b)). The noise therefore causes both a change in the angle of rotation and a shift of the rotational axis, so that the final qubit state will not be the one intended.

In general, for noise that is time dependent, deriving a simple closed-form expression for the propagator U(t) is not possible, as the Hamiltonian will generally not commute with itself at different times (from here on we omit the first argument of a propagator when it is zero). Therefore, rather than attempt to solve (2) directly, we introduce the 'control propagator' U_c(t), defined as the solution to the noise-free Schrodinger equation idU_c(t)/dt = H_c(t)U_c(t). We may then write $U(t)=U_{\mathrm {c}}(t)\tilde {U}(t)$, with the 'error propagator' $\tilde {U}(t)$ capturing any deviation from the control evolution due to noise [20, 37]. Substituting $U(t)=U_{\mathrm {c}}(t)\tilde {U}(t)$ into (2), and defining the 'toggling frame' Hamiltonian

$$\begin{equation} \tilde{H}_{0}(t)\equiv U^{\dagger}_{\mathrm{c}}(t)H_{0}(t)U_{\mathrm{c}}(t), \end{equation} \tag{ 6 }$$

we find that $\tilde {U}(t)$ satisfies the modified Schrodinger equation

$$\begin{equation} \mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}\tilde{U}(t)=\tilde{H}_{0}(t)\tilde{U}(t). \end{equation} \tag{ 7 }$$

In the absence of noise, executing a target operation Q over a time interval [0,τ] simply requires one to engineer a control Hamiltonian such that U_c(τ) = Q. More generally, we retain the requirement that U_c(τ) = Q, so that $U(\tau )=Q\tilde {U}(\tau )$. The goal of DES is then to reduce $\tilde {U}(\tau )$ to the identity I.

The error propagator $\tilde {U}(\tau )$, being an element of SU(2), may be interpreted as a rotation of the qubit Bloch vector through some angle 2a(τ), about an axis $\hat{\mathbf{a}}(\tau )$. Hence

$$\begin{equation} \tilde{U}(\tau)=\exp{\left[-\mathrm{i}\mathbf{a}(\tau){\boldsymbol{\sigma}}\right]}, \end{equation} \tag{ 8 }$$

where $\mathbf {a}(\tau )\equiv a(\tau )\hat{\mathbf{a}}(\tau )$ is the real-valued 'error vector' with components a_i(τ), for i∈{x,y,z} and magnitude $|a(\tau )|=\left (\mathbf {a}(\tau )\mathbf {a}^{\mathrm {T}}(\tau )\right )^{1/2}$.

From equations (6)–(8), we see that the error vector is determined by the toggling frame Hamiltonian $\tilde {H}_{0}(t)$, which itself derives from the action of the control propagator U_c(t) on H₀(t). The aforementioned homomorphism between SU(2) and SO(3) enables us to identify U_c(t) with a three-dimensional 'control matrix' R(t)∈SO(3) [32]

$$\begin{equation} \tilde{H}_{0}(t)=\sum_{i=x,y,z}\beta_{i}(t)U^{\dagger}_{\mathrm{c}}(t)\sigma_{i}U_{\mathrm{c}}(t) \end{equation} \tag{ 9a }$$

$$\begin{equation} \hspace{2.5pc} =\sum_{i,j=x,y,z}\beta_{i}(t)R_{ij}(t)\sigma_{j} \end{equation} \tag{ 9b }$$

$$\begin{equation} \hspace{2.5pc}=\boldsymbol{\beta}(t)\boldsymbol{R}(t){\boldsymbol{\sigma}} . \end{equation} \tag{ 9c }$$

$$\begin{equation} R_{ij}(t)\equiv[\boldsymbol{R}(t)]_{ij}={\textrm{Tr}}{\left(U_{\mathrm{c}}^{\dagger}(t)\sigma_{i}U_{\mathrm{c}}(t)\sigma_{j}\right)}/2 \end{equation} \tag{ 10 }$$

for i,j∈{x,y,z}.

The control matrix may be written compactly as R(t) ≡ (R_x(t),R_y(t),R_z(t))^T, where R_i(t) ≡ (R_ix(t),R_iy(t),R_iz(t)) is a row vector comprised of the ith row of the control matrix. Using this notation, (9c) becomes

$$\begin{equation} \tilde{H}_{0}(t)=\sum_{i=x,y,z}\beta_{i}(t)\boldsymbol{R}_{i}(t){\boldsymbol{\sigma}}. \end{equation} \tag{ 11 }$$

We will use this form in below, where we seek to write the error vector in terms of the noise components and elements of the control matrix.

An explicit expression for the error vector a(τ) can be calculated by finding a solution to (7) that has the exponential form (8). At the cost of expressing the exponent as an infinite series, the Magnus expansion gives the desired result [38–40]. Specifically, given equation (7), we may write $\tilde {U}(\tau )=\exp {[-\mathrm {i}\Phi (\tau )]}$, where

$$\begin{equation} \Phi(\tau)=\sum_{\mu}^{\infty}\Phi_{\mu}(\tau). \end{equation} \tag{ 12a }$$

The first three terms in the expansion are

$$\begin{equation} \fl \Phi_{1}(\tau)=\int_{0}^{\tau}\mathrm{d}t\tilde{H}_{0}(t), \end{equation} \tag{ 12b }$$

$$\begin{equation} \fl \Phi_{2}(\tau)=-\frac{\mathrm{i}}{2}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\left[\tilde{H}_{0}(t_{1}),\tilde{H}_{0}(t_{2})\right], \end{equation} \tag{ 12c }$$

$$\begin{eqnarray} &&\fl \Phi_{3}(\tau)=-\frac{1}{6}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3} \left\{\left[\tilde{H}_{0}(t_{1}),\left[\tilde{H}_{0}(t_{2}),\tilde{H}_{0}(t_{3})\right]\right]+\left[\tilde{H}_{0}(t_{3}),\left[\tilde{H}_{0}(t_{2}),\tilde{H}_{0}(t_{1})\right]\right]\right\}\nonumber\\ \end{eqnarray} \tag{ 12d }$$

Substituting the toggling frame Hamiltonian (11) into the Magnus expansion, and using the vector identity [uσ,vσ] = 2i(u × v)σ, for u, $\mathbf {v}\in \mathbb {{R}}^3$, we find that we are able to write the error vector in the form of an infinite series

$$\begin{equation} \mathbf{a}(\tau)=\sum_{\mu}^{\infty}\mathbf{a}_{\mu}(\tau), \end{equation} \tag{ 13a }$$

where

$$\begin{equation} \mathbf{a}_{1}(\tau)=\sum_{i=x,y,z}\int^{\tau}_{0}\mathrm{d}t\beta_{i}(t)\mathbf{R}_{i}(t), \end{equation} \tag{ 13b }$$

$$\begin{equation} \mathbf{a}_{2}(\tau)=\sum_{i,j=x,y,z}\int^{\tau}_{0}\mathrm{d}t_{1}\int^{t_{1}}_{0}\mathrm{d}t_{2} \beta_{i}(t_{1})\beta_{j}(t_{2})\tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2}), \end{equation} \tag{ 13c }$$

$$\begin{equation} \mathbf{a}_{3}(\tau)=\frac{2}{3}\sum_{i,j,k=x,y,z}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\beta_{i}(t_{1})\beta_{j}(t_{2})\beta_{k}(t_{3})\tilde{\boldsymbol{R}}_{ijk}(t_{1},t_{2},t_{3}). \end{equation} \tag{ 13d }$$

$\tilde {\boldsymbol {R}}_{ij}(t_{1},t_{2}))\equiv [\boldsymbol {R}_{i}(t_{1})\times \boldsymbol {R}_{j}(t_{2})]$

$\tilde {\boldsymbol {R}}_{ijk}(t_{1},t_{2},t_{3})\equiv [\boldsymbol {R}_{i}(t_{1})\times [\boldsymbol {R}_{j}(t_{2})\times \boldsymbol {R}_{k}(t_{3})]] +[\boldsymbol {R}_{k}(t_{3})\times [\boldsymbol {R}_{j}(t_{2})\times \boldsymbol {R}_{i}(t_{1})]]$

$\beta _{i_{1}}(t_{1})\beta _{i_{2}}(t_{2})\ldots \beta _{i_{n}}(t_{n})\tilde {\boldsymbol {R}}_{i_{1}i_{2}\ldots i_{n}}(t_{1},t_{2},\ldots ,t_{n})$

$\tilde {\boldsymbol {R}}_{i_{1}i_{2}\ldots i_{n}}(t_{1},t_{2},\ldots ,t_{n})$

Given (13a)–(13d), we have a series expansion for the error vector in which the components of the noise vector and the elements of the control matrix appear explicitly. By defining an appropriate metric for the fidelity of a control operation, in terms of the error vector, we can use these results to evaluate the impact of noise on an arbitrary control operation on an ensemble of identical noninteracting qubits.

The individual terms in the series expansion of the fidelity (18) rely on time-domain correlation and cross-correlation functions and convolution with the multidimensional control matrix. Using the properties of the Fourier transform, it is experimentally more convenient to rewrite these terms in the frequency domain. Generally, we can define the Fourier transform $\mathcal {S}_{i_{1}\ldots i_{n}}(\omega _{1},\ldots ,\omega _{n})$ of an n-point cross-correlation function via

$$\begin{eqnarray}&&\fl \langle\beta_{i_{1}}(t_{1})\beta_{i_{2}}(t_{2})\ldots \beta_{i_{n}}(t_{n})\rangle \equiv\frac{1}{(2\pi)^{n}}\int \mathrm{d}\omega_{1}\ldots \int \mathrm{d}\omega_{n} \mathcal{S}_{i_{1}\ldots i_{n}}(\omega_{1},\ldots ,\omega_{n})\,\mathrm{e}^{\mathrm{i}(\omega_{1}t_{1}+\cdots+\omega_{n}t_{n})}. \end{eqnarray} \tag{ 20 }$$

The fidelity (18) can then be rewritten as

$$\begin{eqnarray}&&\fl \mathcal{F}_{\rm av}=1-\sum_{n=2}^{\infty} \left\{\frac{1}{(2\pi)^{n}}\sum_{i_{1}\ldots i_{n}}\int \mathrm{d}\omega_{1}\ldots\int \mathrm{d}\omega_{n} \mathcal{S}_{i_{1}\ldots i_{n}}(\omega_{1},\ldots ,\omega_{n})\mathcal{R}{i_{1}\ldots i_{n}}(\omega_{1},\ldots ,\omega_{n}) \right\}, \end{eqnarray} \tag{ 21 }$$

where $\mathcal {R}_{i_{1}\ldots i_{n}}(\omega _{1},\ldots ,\omega _{n})$ is determined solely by the control matrix.

For noise that is sufficiently weak, the n = 2 term in (21) will be the dominant contributor to fidelity loss in the short term. If the noise is also wide sense stationary then the two-point cross-correlation functions depend only on the time difference t₂ − t₁, so that

$$\begin{equation} \langle\beta_{i}(t_{1})\beta_{j}(t_{2})\rangle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega S_{ij}(\omega)\,\mathrm{e}^{\mathrm{i}\omega(t_{2}-t_{1})}, \end{equation} \tag{ 22 }$$

where S_ij(ω) is the cross-power spectral density between the random variables β_i(t) and β_j(t). We then have

$$\begin{equation} \mathcal{F}_{\rm av}\simeq1-\frac{1}{2\pi}\sum_{i,j,k=x,y,z}\int_{-\infty}^{\infty}\frac{{\rm d}\omega}{\omega^2}S_{ij}(\omega)R_{jk}(\omega)R^{*}_{ik}(\omega), \end{equation} \tag{ 23 }$$

where

$$\begin{equation} R_{ij}(\omega)\equiv-\mathrm{i}\omega\int_{0}^{\tau}\mathrm{d}tR_{ij}(t)\,\mathrm{e}^{\mathrm{i}\omega t} \end{equation} \tag{ 24 }$$

are the elements of the control matrix in the frequency domain (again we suppress the explicit dependence on the total time τ). The form of (23) is reminiscent of the filter function introduced in previous literature on DES, to which we will return later. Thus, in the weak noise regime (to be defined in the appendix, for Gaussian noise), it is possible to understand the fidelity of a qubit's unitary evolution in terms of spectral properties of the noise and the control [11, 21].

In the next section, we will show how the control matrix may be evaluated for the paradigmatic case of a piecewise-constant control sequence. In doing so, we explicitly demonstrate the connection between our generalized approach and the well-known filter-function formalism.

For a sequence of n consecutive unitary control operations (or control 'pulses') P₁,...,P_n, we divide the total sequence time τ into n corresponding time intervals, such that the lth operation P_l is executed over the interval [t_l−1,t_l], i.e. P_l = U_c(t_l,t_l−1), for l∈{1,2,...,n}, where t_n ≡ τ and t₀ ≡ 0 (figure 2). Letting P₀ ≡ I, and defining the cumulative operators Q_l ≡ P_lP_l−1...P₀ so that Q ≡ Q_n, the elements of the control matrix (10) take the form

$$\begin{equation} R_{ij}(t)=\frac{1}{2}\sum_{l=1}^{n}G^{(l)}(t) {\textrm{Tr}}{\left(Q_{l-1}^{\dagger}U_{\mathrm{c}}^{\dagger}(t,t_{l-1})\sigma_{i} U_{\mathrm{c}}(t,t_{l-1})Q_{l-1}\sigma_{j}\right)}, \end{equation} \tag{ 25 }$$

where the function G^(l)(t) ≡ Θ[t − t_l−1] − Θ[t − t_l] has unit value within the lth time interval and is zero elsewhere.

Zoom In Zoom Out Reset image size

Equation (25) can be rewritten in terms of the individual control matrices R^P_l(t − t_l−1) for each of the n operations P_l in the sequence, defined relative to the pulse start times t_l−1, by letting

$$\begin{equation} U_{\mathrm{c}}^{\dagger}(t,t_{l-1})\sigma_{i}U_{\mathrm{c}}(t,t_{l-1})=\sum_{j=x,y,z}R_{ij}^{P_{l}}(t-t_{l-1})\sigma_{j} \end{equation} \tag{ 26 }$$

for i∈{x,y,z}. Again, the components of R^P_l(t − t_l−1) are obtained from the Hilbert–Schmidt inner product

$$\begin{equation} R_{ij}^{P_{l}}(t-t_{l-1})\equiv[\boldsymbol{R}^{P_{l}}(t-t_{l-1})]_{ij}=\textstyle\frac{1}{2}{\textrm{Tr}}{\left(U_{\mathrm{c}}^{\dagger}(t,t_{l-1})\sigma_{i}U_{\mathrm{c}}(t,t_{l-1})\sigma_{j}\right)} \end{equation} \tag{ 27 }$$

for i,j∈{x,y,z}. Substituting (26) into (25), and using the linearity of the trace operation, we find that we can write the control matrix for the pulse sequence in the compact form

$$\begin{equation} \boldsymbol{R}(t)=\sum_{l=1}^{n}G^{(l)}(t)\boldsymbol{R}^{P_{l}}(t-t_{l-1})\boldsymbol{\Lambda}^{(l-1)}, \end{equation} \tag{ 28 }$$

where the matrix Λ^(l−1) has components

$$\begin{equation} \Lambda_{ij}^{(l-1)}\equiv[\boldsymbol{\Lambda}^{(l-1)}]_{ij}=\textstyle\frac{1}{2}{\textrm{Tr}}{\left(Q_{l-1}^{\dagger}\sigma_{i}Q_{l-1}\sigma_{j}\right)} \end{equation} \tag{ 29 }$$

for i,j∈{x,y,z}. Using (24) to define the control matrix in the frequency domain, we have

$$\begin{equation} \boldsymbol{R}(\omega)=\sum_{l=1}^{n}\,\mathrm{e}^{\mathrm{i}\omega t_{l-1}}\boldsymbol{R}^{P_{l}}(\omega)\boldsymbol{\Lambda}^{(l-1)}, \end{equation} \tag{ 30a }$$

where

$$\begin{equation} \boldsymbol{R}^{P_{l}}(\omega)\equiv-\mathrm{i}\omega\int_{0}^{t_{l}-t_{l-1}}\mathrm{d}t\,\mathrm{e}^{\mathrm{i}\omega t}\boldsymbol{R}^{P_{l}}(t) \end{equation} \tag{ 30b }$$

is the frequency-domain control matrix for the lth pulse, P_l.

With the above expressions and a knowledge of the statistical properties of the noise we can, in combination with the results of section 3, find an approximate value for the fidelity of an arbitrary piecewise-defined control sequence in a weak noise environment. This approach is straightforward to implement by hand or via numerics. Ultimately, the computational efficiency of calculating the control matrix and spectral overlap integrals far exceeds that of brute-force numerical simulations of the evolution of the Bloch vector [41].

In many of the physical systems with the potential for use in nascent quantum information technologies, the characteristic single-qubit dephasing time T₂ is considerably smaller than the relaxation time T₁, and governed by different physical processes. It is important, therefore, that we understand the effects of purely dephasing noise on single-qubit control. To this end, we consider a model system in which the noise vector β(t) has only a z-component and the noise Hamiltonian (4) reduces to

$$\begin{equation} H_{0}(t)=\beta_{z}(t)\sigma_{z}. \end{equation} \tag{ 31 }$$

In this context, only the third row of the control matrix is required to calculate the fidelity, and we refer to the corresponding row vector R_z(t) ≡ (R_zx(t),R_zy(t),R_zz(t)) as the 'control vector'. From (23), the lowest-order contribution of noise to the fidelity (18) is captured by the term

$$\begin{equation} \langle a_{1}^{2}\rangle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{z}(\omega)F^{(1)}_{z}(\omega), \end{equation} \tag{ 32 }$$

where the first-order dephasing 'filter function' F⁽¹⁾_z(ω) is simply the square modulus of the frequency-domain control vector

$$\begin{equation} F^{(1)}_{z}(\omega)\equiv \sum_{i=x,y,z}|R_{zi}(\omega)|^{2}. \end{equation} \tag{ 33 }$$

If we further assume that the noise is Gaussian, with a mean value of zero, then all correlation functions 〈β_z(t₁)...β_z(t_n)〉 for which n is odd are equal to zero. The Gaussian moment theorem also allows us to write all remaining higher-order correlations in terms of only the simplest two-point correlation function (see the appendix). In the frequency domain, the statistical properties of the noise are then captured entirely by the power spectral density S_z(ω), and the fidelity (18) can be written as the infinite sum [41]

$$\begin{eqnarray} &&\fl \mathcal{F}_{\rm av}=1-\left[\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{z}(\omega)F^{(1)}_{z}(\omega)+\frac{1}{(2\pi)^{2}}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{z}(\omega)\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega'}{\omega'^2} S_{z}(\omega')F^{(2)}_{zz}(\omega,\omega')+\cdots\right].\nonumber\\ \end{eqnarray} \tag{ 34 }$$

Here, the effect of the control on the qubit fidelity is described by the functions F^(p)_z..., for p = 1,2,..., which depend only on the control vector R_z. The complexity of the filter functions increases rapidly with increasing order. However, F⁽¹⁾_z(ω) can be quite simple in form and reasonably straightforward to calculate. In the weak noise regime we ignore the higher-order terms in (34) and write

$$\begin{equation} \mathcal{F}_{\rm av}(\tau)\simeq\textstyle\frac{1}{2}\left\{1+\exp{\left[-\chi(\tau)\right]}\right\}, \end{equation} \tag{ 35a }$$

where the rate of fidelity decay is given by

$$\begin{equation} \chi(\tau)\equiv\langle a_{1}^{2}\rangle=\frac{1}{\pi}\int_{0}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{z}(\omega)F^{(1)}_{z}(\omega). \end{equation} \tag{ 35b }$$

Thus, the decay rate is determined by a simple overlap integral between the power spectrum of the noise and a 'filter function' representing the control in the Fourier domain. This result is similar to the familiar expression for the coherence of a DD sequence in the case of unbounded controls [21] (see also section 5). The difference here being that the filter function F⁽¹⁾_z captures only the lowest-order nontrivial effect of the control on the fidelity in a dephasing environment, rather than providing an exact solution in the bang–bang limit.

For a piecewise-defined control sequence, subject to purely dephasing noise, equation (25) reduces to the following expression for the components of the control vector:

$$\begin{equation} R_{zi}(t)=\textstyle\frac{1}{2}\sum_{l=1}^{n}G^{(l)}(t) {\textrm{Tr}}{\left(Q_{l-1}^{\dagger}U_{\mathrm{c}}^{\dagger}(t,t_{l-1})\sigma_{z} U_{\mathrm{c}}(t,t_{l-1})Q_{l-1}\sigma_{i}\right)}. \end{equation} \tag{ 36 }$$

Similarly, (30a) and (30b) become

$$\begin{equation} \boldsymbol{R}_{z}(\omega)=\sum_{l=1}^{n}\mathrm{e}^{\mathrm{i}\omega t_{l-1}}\boldsymbol{R}_{z}^{P_{l}}(\omega)\boldsymbol{\Lambda}^{(l-1)} \end{equation} \tag{ 37a }$$

and

$$\begin{equation} \boldsymbol{R}_{z}^{P_{l}}(\omega)\equiv-\mathrm{i}\omega\int_{0}^{t_{l}-t_{l-1}}{\rm d}t\,\mathrm{e}^{{\rm i}\omega t}\boldsymbol{R}_{z}^{P_{l}}(t), \end{equation} \tag{ 37b }$$

respectively, where the components of the lth pulse control vector are

$$\begin{equation} R_{zj}^{P_{l}}(t-t_{l-1})=\textstyle\frac{1}{2}{\textrm{Tr}}{\left(U_{\mathrm{c}}^{\dagger}(t,t_{l-1})\sigma_{z}U_{\mathrm{c}}(t,t_{l-1})\sigma_{j}\right)} \end{equation} \tag{ 38 }$$

for j∈{x,y,z}.

If we assume that during each of the n time intervals comprising the sequence the applied control affects a steady rotation of the qubit Bloch vector through an angle θ_l, about an axis ${\hat{\boldsymbol{n}}}_{l}=\cos {(\phi _{l})}\hat {x}+\sin {(\phi _{l})}\hat {y}$ in the xy-plane, then for t∈[t_l−1,t_l]

$$\begin{equation} U_{\mathrm{c}}(t,t_{l-1})=\exp{\left[\frac{-\mathrm{i}}{2}\Omega_{l}\left(t-t_{l-1}\right)\sigma_{\phi_{l}}\right]}, \end{equation} \tag{ 39 }$$

where Ω_l ≡ θ_l/(t_l − t_l−1) is the rotation rate and σ_ϕl ≡ cos(ϕ_l)σ_x + sin(ϕ_l)σ_y. From (38) and (37b) we can show that the frequency-domain control vector for the lth operation has components

$$\begin{equation}\fl R_{zj}^{P_{l}}(\omega)=\frac{\omega}{\omega^2-\Omega_{l}^{2}}\left\{\delta_{zj}\left[\mathrm{i}\Omega_{l}g_{l}(\omega)-\omega f_{l}(\omega)\right] +\left[\Omega_{l}f_{l}(\omega)-\mathrm{i}\omega g_{l}(\omega)\right]{\textrm{Tr}}{(\sigma_{\phi_{l}}\sigma_{z}\sigma_{j})}\right\}, \end{equation} \tag{ 40 }$$

where f_l(ω) ≡ e^{iω(t_l−t_l−1)}cos(θ_l) − 1 and g_l(ω) ≡ e^{iω(t_l−t_l−1)}sin(θ_l). Although (40) has been derived assuming rotations take place about axes in the xy-plane only, periods of free evolution (no rotation) and the z-axis rotations are easily accounted for by setting θ_l = 0 and Ω_l = 0 in (40), whenever one of these operations occurs in the sequence. This is possible simply because, for dephasing noise, both of these operations commute with the noise operator.

With these analytical tools, the process of finding the operational fidelity of an arbitrary piecewise-constant control sequence in the presence of dephasing noise may be summarized as follows:

• (i)  
  calculate the matrix defined in equation (29) for each step in the chain of control operations;
• (ii)  
  substitute the above result and equation (40) into (37a);
• (iii)  
  repeat for all time steps l∈{1,...,n} to find the control vector in the Fourier domain;
• (iv)  
  find the first-order filter function from (33);
• (v)  
  use (35a) and (35b) to calculate approximate expression for the fidelity valid for weak dephasing noise.

In the next section, we employ this approach to study the impact of pulse imperfections on the performance of DD sequences used in quantum memory.

Beyond the bang–bang limit, the effectiveness of a DD sequence will depend on both the temporal spacing and form of the pulses used. To find optimal error-suppressing sequences for a given environment or physical system, either of these variables may be changed independently of the other, and these differences should be easily separated in an analytical framework. We are therefore motivated to express the sequence control vector Rz(ω) in terms that can be separated into those dependent only on pulse location and those dependent on pulse type (expressed in terms of the individual pulse control vectors RPlz(ω)).

As in section 4, we start with a piecewise-constant sequence of n pulses, however, to model a DD sequence we allow for a period of free evolution (i.e. no control operation) both before and after each pulse, so that immediately following the lth pulse, we have Q_l = P_lIP_l−1I...IP₀, for l∈{1,2,...,n}. Within the total sequence time τ (including both control pulses and free evolution periods), the lth pulse has a width τ_pl and is centered at a time δ_lτ, where 0 ⩽ δ_l ⩽ 1 and δ_n+1 = 1 (see figure 3). Using this notation the sequence control vector is

$$\begin{eqnarray} &&\fl \mathbf{R}_{z}(\omega)=\hat{\mathbf{z}}\left[1-\mathrm{e}^{\mathrm{i}\omega\tau}+\sum_{l=1}^{n}\mathrm{e}^{\mathrm{i}\omega\delta_{l}\tau} \left(\mathrm{e}^{\mathrm{i}\omega\tau_{p_{l}}/2}\boldsymbol{\Lambda}^{(l)}-\mathrm{e}^{-\mathrm{i}\omega\tau_{p_{l-1}}/2}\boldsymbol{\Lambda}^{(l-1)}\right)\right]+\sum_{l=1}^{n}\mathrm{e}^{\mathrm{i}\omega\delta_{l}\tau-\tau_{p_{l}}/2}\mathbf{R}_{z}^{P_{l}}(\omega) \boldsymbol{\Lambda}^{(l-1)},\nonumber\\ \end{eqnarray} \tag{ 42 }$$

the form of the pulse appearing explicitly only in the last term.

While DD sequences have traditionally been composed of identical π-pulses executed about a common axis, equation (42) is quite general. The only restriction being that, in the absence of noise, the net effect is to execute the identity operation. Hence, it may be applied to the analysis of unusual DD schemes such as the recently proposed KDD (Knill dynamical decoupling), based on the Knill pulse [44]. However, since our aim here is simply to provide an example application of our method, we will limit our attention to more conventional sequences. Specifically, we assume identical pulses, each of width τ_p and executing a net π-rotation of the qubit Bloch vector about the x-axis, i.e. P_l = σ_x, for l∈{1,2,...,n}. In this case, the pulse control vector R^P _z(ω) is independent of l and has only z and y components. The matrix Λ^(l−1) is diagonal with Λ^(l−1)_yy and Λ^(l−1)_zz alternating between 1 and −1 with each successive pulse. The resulting sequence control vector has the two components

$$\begin{equation} R_{zz}(\omega)=1-\mathrm{e}^{\mathrm{i}\omega\tau}+\left[2\cos{(\omega\tau_{p}/2)}-\mathrm{e}^{-\mathrm{i}\omega\tau_{p}/2}R_{zz}^{P}(\omega)\right] \sum_{l=1}^{n}(-1)^{l}\,\mathrm{e}^{\mathrm{i}\omega\delta_{l}\tau} \end{equation} \tag{ 43a }$$

and

$$\begin{equation} R_{zy}(\omega)=-\mathrm{e}^{-\mathrm{i}\omega\tau_{p}/2}R_{zy}^{P}(\omega)\sum_{l=1}^{n}(-1)^{l}\,\mathrm{e}^{\mathrm{i}\omega\delta_{l}\tau}. \end{equation} \tag{ 43b }$$

It is worth noting here that, when the form of the pulse is ignored altogether (i.e. when both R^P _zy(ω) and R^P _zz(ω) go to zero), the control vector reduces to a single component

$$\begin{equation} R_{zz}(\omega)=1-\mathrm{e}^{\mathrm{i}\omega\tau}+2\cos{(\omega\tau_{p}/2)} \sum_{l=1}^{n}(-1)^{l}\,\mathrm{e}^{\mathrm{i}\omega\delta_{l}\tau}, \end{equation} \tag{ 44 }$$

a result derived previously by assuming finite-width pulses during which noise was ignored [24, 25]. And, in the bang–bang limit of infinitely narrow pulses (τ_p → 0) we have the standard result (for even n)

$$\begin{equation} R_{zz}(\omega)=1-\mathrm{e}^{\mathrm{i}\omega\tau}+2\sum_{l=1}^{n}(-1)^{l}\,\mathrm{e}^{\mathrm{i}\omega\delta_{l}\tau} \end{equation} \tag{ 45 }$$

with F_z(ω) = |R_zz(ω)|².

Equations (43a) and (43b) capture modification of dephasing dynamics due to the duration and form of the pulses as well as depolarization effects occurring due to pulse errors derived from a fluctuating detuning from resonance (i.e. the impact of dephasing noise). We will move forward using these expressions in order to treat experimentally relevant cases of DD incorporating various pulse forms including error-suppressing dynamically corrected gates.

For a DD sequence we are interested in physical control operations that implement a logical NOT gate. The simplest approach is to employ 'primitive' π_X-pulses, implemented by applying the ideal control Hamiltonian H_c = Ωσ_x/2 over a time interval [0,τ_π], where Ω ≡ π/τ_π. This induces a simple rotation of the Bloch vector through an angle of π radians about the x-axis. In the presence of time-varying dephasing noise, the bare gate is modified by an error term that, to first order, results in a small additional rotation about an axis that has both y- and z-components. From equations (37a) and (40) we find that the two corresponding components of control vector are

$$\begin{equation} R^{P}_{zz}\equiv R^{({\rm Prim})}_{zz}(\omega)=\frac{\omega^{2}}{\omega^{2}-\Omega^{2}}\left(\mathrm{e}^{\mathrm{i}\omega\tau_{\pi}}+1\right) \end{equation} \tag{ 46a }$$

and

$$\begin{equation} R^{P}_{zy}\equiv R^{({\rm Prim})}_{zy}(\omega)=\frac{\mathrm{i}\omega\Omega}{\omega^{2}-\Omega^{2}}\left(\mathrm{e}^{\mathrm{i}\omega\tau_{\pi}}+1\right). \end{equation} \tag{ 46b }$$

It is interesting to note the correspondence between the prefactors here and those arising in a master-equation treatment of a driven quantum system in the presence of dissipation [45].

Moving beyond the standard π_X-pulse of finite duration, we consider a gate that has been designed to provide robustness in the presence of random gate errors—the dynamically corrected gate. In previous work we demonstrated how such gates provide enhanced resistance to error in the presence of time-varying dephasing noise [41]. Here we include expressions for a dynamically corrected logical NOT gate consisting of three successive π_X pulses, the second of which is executed at a rotation rate Ω/2, half that of the other two [46]. The control Hamiltonian for this gate is

$$\begin{eqnarray} H_{\mathrm{c}}(t)=\left\{\begin{array}{@{}ll@{}} \Omega\sigma_{x}/2, & 0\leqslant t<\tau_{\pi},\\ \Omega\sigma_{x}/4, & \tau_{\pi}\leqslant t<3\tau_{\pi},\\ \Omega\sigma_{x}/2, & 3\tau_{\pi}\leqslant t\leqslant 4\tau_{\pi}, \end{array}\right. \end{eqnarray} \tag{ 47 }$$

so that τ_p = 4τ_π. In the case of dephasing noise, the control vector again has two components:

$$\begin{equation} R^{P}_{zz}\equiv R^{({\rm DCG})}_{zz}(\omega)=\omega^{2}\left[\frac{p_{1}(\omega)}{\omega^{2}-\Omega^{2}}-\frac{p_{2}(\omega)}{\omega^{2}-(\Omega/2)^{2}}\right], \end{equation} \tag{ 48a }$$

$$\begin{equation} R^{P}_{zy}\equiv R^{({\rm DCG})}_{zy}(\omega) =\mathrm{i}\omega\Omega\left[\frac{p_{1}(\omega)}{\omega^{2}-\Omega^{2}}-\frac{p_{2}(\omega)/2}{\omega^{2}-\left(\Omega/2\right)^{2}} \right], \end{equation} \tag{ 48b }$$

The inset to figure 4(a) shows the filter functions, F_z(ω) ≡ |R^P _zy(ω)|² + |R^P _zz(ω)|², for the primitive and dynamically corrected logical NOT operations described above, as a function of the angular frequency in units of τ⁻¹_π. The rate of growth of the filter function with ω captures the low-frequency 'filtering' performance of DD sequences. Beyond a certain frequency F_z(ω) → 1, above which noise is passed unimpeded. This corresponds to the physical observation that fluctuations fast relative to the shortest interpulse period cannot be effectively suppressed by DD.

Zoom In Zoom Out Reset image size

The filter function's behavior near zero frequency captures the leading-order error susceptibility. We may express F_z(ω)∝(ωτ_π)^2(α+1) as ωτ_π → 0, then α is the order of error suppression for the sequence [47]. Using these expressions, an examination of figure 4 reveals that the relative performance of the primitive and dynamically corrected gates is captured in the slope of the filter function as plotted on a log–log scale [23]; higher slope corresponds to a higher order of error suppression (for noise near zero frequency). In a filter-design framework, we may refer to this slope as a filter rolloff in the stopband, using the relation that the rolloff is 6(α + 1) dB per octave. At low frequencies the primitive π_X-pulse filter function is approximately quadratic, i.e. F_z(ω)∝(ωτ_π)², resulting in α = 0—the same as for free evolution. This is not a surprising result, since there is no error suppression built into the primitive gate. By contrast, the dynamically corrected gate (DCG) filter function scales as F_z(ω)∝(ωτ_π)⁴, giving α = 1, as expected for a construction designed to suppress error to first order.

We now combine the control vectors derived for the pulses themselves with those obtained for the generic DD sequences, considering three paradigmatic cases: (i) bang–bang decoupling, (ii) standard pulses with nonzero τ_p and (iii) DCGs. This is accomplished by calculating R^P _z and inserting the relevant Cartesian components into (43a) and (43b) for the full DD sequence. This study reveals how it is possible to evaluate complex perturbations to quantum control protocols in a straightforward manner using our method.

We compare the error-suppressing performance of CP (Carr and Purcell) and UDD (Uhrig dynamical decoupling) sequences as illustrative examples, accounting for each of the cases outlined above. The CP sequence, devised in the context of NMR [48], is a straightforward extension of the original Hahn spin echo [49]. In an n-pulse CP sequence, the lth pulse has the fractional location δ_l = (l − 1/2)/n. The bang–bang limit CP filter function for n = 6 is shown as a black line in figure 4(a) for n = 6. It exhibits an order of error suppression of α = 2, as do all CP sequences with n ⩾ 2, giving a roll off of 18 dB per octave. This is a limiting case in which pulse effects are entirely ignored.

If we now take into account the width of the pulses by assuming that they take the form of primitive NOT gates of width τ_p, as described above, then the components of the pulse control vector R^P _z(ω) are given by equations (46a) and (46b). The resulting filter functions, indicated by the blue lines in figure 4(a), show a reduced order of error suppression α = 2 → 1, with a corresponding reduction of the filter rolloff to 12 dB per octave. As the pulse width grows as a fraction of the entire sequence, the relative importance of the pulses grows, as expected. This is manifested as an expansion of the range over which the modification to the filter function due to pulse effects dominates the bang–bang filter function.

Despite the reduction in the order of error suppression, there is still noise cancelation, regardless of the pulse width, since the sequence is in essence an Eulerian dynamical decoupling (EDD) sequence [50]. Replacing the primitive pulses with dynamically corrected composite gates and the components of the pulse control vector with (48a) and (48b), we almost completely restore the original 18 dB per octave rolloff for the CP sequences (figure 4(a)). These filter functions lie approximately beneath the bang–bang filter function in this figure. This observation confirms that the DCG provides one order of error cancelation during the applied pulses. In essence, the deleterious effect of error accumulation during nonzero duration control pulses can be mitigated by choosing a compensating pulse design.

Studying UDD is more instructive because the n-pulse UDD sequence, with relative pulse locations given by δ_l = sin²[πl/(2n + 2)], ensures that (in the bang–bang limit) the first n-derivatives of R_zz(ω) vanish at ω = 0 [21, 43], giving an order of error suppression that increases linearly with n. The effects of pulse width and shape on a UDD sequence are therefore potentially more dramatic. For sufficiently long pulses we see, from figure 4(b), that the general benefits derived from the use of a UDD sequence can largely be destroyed as the error susceptibility of the pulses dominates the suppression provided by the pulse timing. Again, we observe that the reduction in the order of error suppression arising from addition of nonzero duration π_X-pulses can be partially offset by use of a DCG. However, because the DCG only suppresses error to first order, our ability to recover the original order of error suppression provided by a UDD sequence with large n is reduced.

Nonetheless, these results indicate that judicious choice of DCG or other compensating pulse protocols provides a path to mitigate the effects of pulse errors in DD sequences. This is especially important in long-storage settings where large pulse numbers may be employed in order to effectively suppress dephasing errors [51]. These observations represent an important validation of our method, and show that it provides a straightforward approach to account for a variety of pulse duration, modulation and shape effects.

A sufficient condition for the convergence of the Magnus expansion is [52]

$$\begin{equation} \int_{0}^{\tau}\mathrm{d}t\|\tilde{H}_{0}(t)\|_{\rm op}<\pi. \end{equation} \tag{ 49 }$$

The operator norm $\|\tilde {H}_{0}(t)\|_{\mathrm { op}}$ is the smallest number K ⩾ 0 such that $\|\tilde {H}_{0}(t)|\psi \rangle \|_{\mathcal {H}_{S}}\leqslant K\||\psi \rangle \|_{\mathcal {H}_{S}}$ for all $|\psi \rangle \in \mathcal {H_{S}}$, where $\||\psi \rangle \|_{\mathcal {H}_{S}}\equiv \sqrt {\langle \psi |\psi \rangle }$ is the usual vector norm defined on the system Hilbert space $\mathcal {H_{S}}$ [53]. It is a simple matter to derive the intuitively obvious result that

$$\begin{equation} \|\tilde{H}_{0}(t)\|_{\rm op}=\|\boldsymbol{\beta}(t)\|=\left(\beta_{x}(t)^{2}+\beta_{y}^{2}(t)+\beta_{y}^{2}(t)\right)^{1/2}, \end{equation} \tag{ 50 }$$

i.e. the 'size' of the toggling frame Hamiltonian is measured by the magnitude of the noise vector.

In an experiment, for each physical realization of the noise process, the ith noise component β_i(t) will have some maximum absolute value β^(m)_i, over the interval [0,τ]. We can then write ∥β(t)∥ ⩽ ∥β^(m)∥, for all t, where β^(m) = (β^(m)_x,β^(m)_y,β^(m)_z). It immediately follows that

$$\begin{equation} \int_{0}^{\tau}\mathrm{d}t\|\boldsymbol{\beta}(t)\|\leqslant \|\boldsymbol{\beta}^{(m)}\|\tau \end{equation} \tag{ 51 }$$

and convergence is assured if ∥β^(m)∥τ < π. The problem is, of course, that the maxima β^(m)_i are unknown quantities. Nonetheless, if the noise is stationary and has well-defined root mean square values δβ_i ≡ (〈β_i(0)²〉)^1/2, for i = x,y,z, then for small total operation times τ, letting β^(m)_i = C_mδβ_i, for a sufficiently large value of C_m > 0, will mean that contributions from elements of the ensemble not satisfying the convergence condition

$$\begin{equation} \xi<\pi/C_{m}, \end{equation} \tag{ 52 }$$

where ξ ≡ 〈β(0)β(0)^T〉^1/2τ can be ignored without significant error.

It can be shown (see the appendix) that this same parameter ξ represents a bound on the strength of high-order terms in the expansion for the fidelity (18) and serves as 'smallness parameter' providing insight into when the first-order filter function is sufficient to obtain a good approximation for the overall fidelity. Essentially, we require ξ² ≪ 1 for the first-order approximation to hold. In general, this is sufficient to guarantee convergence of the Magnus expansion for all but a statistically insignificant proportion of qubits in the ensemble.

Finally, we present an instructive example to illustrate that our techniques are amenable to detailed understanding of the limits of applicability and scaling behaviors of the calculated operational fidelity. To this end, we compare the performance of our analytical expressions with numerical simulations in which we explicitly time evolve the Schrodinger equation. As a test case, we simulate a primitive NOT-gate (π_X-pulse), in the presence of pure dephasing noise, $H(t)=\beta (t)\hat \sigma _z/2+\Omega _R\hat \sigma _x/2$, and study the behavior of the fidelity (18) as a function of the control field strength.

For concreteness we pick a Gaussian noise power spectrum

$$\begin{equation} S_{z}(\omega)=\sqrt{2\pi}\frac{\delta\beta^2}{4\sigma}\,\mathrm{e}^{-\omega^2/(2\sigma^2)}. \end{equation} \tag{ 53 }$$

Here δβ ≡ δβ_z, σ is the bandwidth of the noise field, and the spectrum is properly normalized so that it is connected to the auto-correlation C(t) of the noise through a Fourier transform

$$\begin{equation} C(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty}S_{z}(\omega)\,\mathrm{e}^{\mathrm{i}\omega t}\,\mathrm{d}\omega=\frac{\delta\beta^2}{4}\,\mathrm{e}^{-t^2/(2\tau_c^2)}, \end{equation} \tag{ 54 }$$

where τ_c = 1/σ. For our simulations we vary the spectral bandwidth of the Gaussian spectrum. The noise strength is parametrized by its root-mean-square deviation which we choose as δβ = 0.5 in each case.

To simulate the noise, we generate random classical noise trajectories by summing noise Fourier frequency components, each weighed by a random, normally distributed amplitude multiplied by the square root of the noise power at that frequency (as given by the power spectrum), and with a random phase. The gate fidelity is found by averaging over 100 realizations of the noise, and compared to the analytical results including only first-order terms in the Magnus expansion.

Analytic and numeric results are plotted in figure 5. We observe excellent agreement between the first-order filter function analytical prediction and detailed numerical simulation. We also find that there appears to be two different regimes of behavior in relation to dependence on noise bandwidth (or identically the correlation time). The calculated gate fidelity appears to be independent of noise bandwidth in the short time regime τ_π ≪ 1 and deviates from this behavior only for the large bandwidth case σ = 10. We can understand this behavior by analyzing in different temporal regimes, the behavior of the first-order term in the fidelity expansion (18) which is simply

$$\begin{equation} \langle a^2_1\rangle= \frac{1}{4}\int^{\tau_\pi}_0\mathrm{d}t_2\int^{\tau_\pi}_0\mathrm{d}t_1\delta\beta^2\,\mathrm{e}^{-(t_1-t_2)^2/(2\tau_c^2)}\cos{(\Omega_Rt_1)}\cos{(\Omega_Rt_2)} \end{equation} \tag{ 55 }$$

$$\begin{equation} \quad\qquad+ \frac{1}{4}\int^{\tau_\pi}_0\mathrm{d}t_2\int^{\tau_\pi}_0\mathrm{d}t_1\delta\beta^2\,\mathrm{e}^{-(t_1-t_2)^2/(2\tau_c^2)}\sin{(\Omega_Rt_1)}\sin{(\Omega_Rt_2)} \end{equation} \tag{ 56 }$$

$$\begin{equation} \qquad = \frac{1}{4}\int^{\tau_\pi}_0\mathrm{d}t_2\int^{\tau_\pi}_0\mathrm{d}t_1\delta\beta^2\,\mathrm{e}^{-(t_1-t_2)^2/(2\tau_c^2)}\cos{[\Omega_R(t_1-t_2)]}. \end{equation} \tag{ 57 }$$

For short control times, i.e. τ_π/τ_c ≪ 1, this integral reduces to

$$\begin{equation} \langle a^2_1\rangle= \delta\beta^2\tau_\pi^2\left[\frac{1}{\pi^2}-\left(\frac{12-\pi^2}{4\pi^2}\right)\left(\frac{\tau_\pi}{\tau_{c}}\right)^2+\cdots\right]. \end{equation} \tag{ 58 }$$

The lowest-order term in equation (58) is plotted as the green dash-dotted line in figure 5 (left). As demonstrated by both the analytical and the numerical results, one finds that the short-time behavior is independent of bandwidth to lowest order.

Zoom In Zoom Out Reset image size

Bandwidth-dependent effects should only appear when the control time becomes of the same order as the correlation time, reached for the case when σ = 10. To elucidate the behavior in that regime, we evaluate equation (57) approximately when τ_c/τ_π ≪ 1. To lowest order this gives

$$\begin{equation} \langle a^2_1\rangle\approx \sqrt{\pi}\frac{\delta\beta^2}{4}\tau_c\tau_\pi. \end{equation} \tag{ 59 }$$

This approximate solution is plotted as the green dashed line in the left panel of figure 5, illustrating that the full numerics and first-order analytics obey the expected scaling with control time, and that the error magnitude is set by the noise correlation time. For a fixed root-mean-square noise strength, and long control times, a primitive NOT gate is more robust under wide-bandwidth noise conditions.

Since ξ² = δβ²τ²_π ≪ 1 over most of the regime we investigated, higher-order terms in the Magnus expansion are expected to be negligible since they scale as ξ⁴ to lowest order (see the appendix). This expectation is confirmed by the numerics. The numerical results therefore support the scaling argument that ξ² ≪ 1 determines the bounds of validity of the first-order filter-function approximation. The inset to figure 5 (left) shows that the analytical result becomes a poor approximation to the numerics when the condition ξ² ≪ 1 is violated in the vicinity of τ_π ≈ 10. Note that growth in higher-order terms does not necessarily coincide with divergence of the Magnus expansion, but rather a breakdown of the first-order approximation; higher-order filter functions may be incorporated to account for residual error.

Finally, we note that the same arguments we have applied to the primitive NOT gate may also be applied to more general control operations, including those designed to compensate for noise, such as DCGs. First-order error correcting gates like the DCG will reduce the effect of the first-order term in the Taylor expansion of the noise correlation function relative to the higher-order bandwidth-dependent terms. As a consequence, bandwidth effects are expected to be more pronounced for corrected gates, which we have confirmed via numerics. Despite this issue we still obtain good agreement between numerics and analytic filter functions to within factors of order unity, as demonstrated in previous work [41], and in the right panel of figure 5.

In order to determine their relative contributions to the fidelity, we examine the terms in (18), beginning with

$$\begin{equation} \langle a_{1}^{2}\rangle = \sum_{i,j=x,y,z}\int^{\tau}_{0}\mathrm{d}t_{2}\int^{\tau}_{0}\mathrm{d}t_{1} \langle\beta_{i}(t_{1})\beta_{j}(t_{2})\rangle \boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{2}). \end{equation} \tag{ A.2 }$$

We consider only the case in which the three noise components are independent Gaussian random processes, so that 〈β_i(t₁)β_j(t₂)〉 = δ_ij〈β_i(t₁)β_i(t₂)〉 and

$$\begin{equation} \langle a_{1}^{2}\rangle = \sum_{i=x,y,z}\int^{\tau}_{0}\mathrm{d}t_{2}\int^{\tau}_{0}\mathrm{d}t_{1} \langle\beta_{i}(t_{1})\beta_{i}(t_{2})\rangle \boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{2}). \end{equation} \tag{ A.3 }$$

Letting $\langle \beta _{i}(t_{1})\beta _{i}(t_{2})\rangle =\langle \beta _{i}^2(0)\rangle \overline {\langle \beta _{i}(t_{1})\beta _{i}(t_{2})\rangle }$, where 〈β²_i(0)〉 is the mean square value of the ith component of the noise and $|\overline {\langle \beta _{i}(t_{1})\beta _{i}(t_{2})\rangle }|\leqslant 1$, we have

$$\begin{eqnarray} |\langle a_{1}^{2}\rangle|&=&\left|\sum_{i=x,y,z}\langle\beta_{i}^2(0)\rangle\int^{\tau}_{0}\mathrm{d}t_{2}\int^{\tau}_{0}\mathrm{d}t_{1} \overline{\langle\beta_{i}(t_{1})\beta_{i}(t_{2})\rangle} \boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{2})\right|\nonumber\\ &\leqslant&\sum_{i=x,y,z}\langle\beta_{i}^2(0)\rangle\left|\int^{\tau}_{0}\mathrm{d}t_{2}\int^{\tau}_{0}\mathrm{d}t_{1} \overline{\langle\beta_{i}(t_{1})\beta_{i}(t_{2})\rangle} \boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{2})\right|\nonumber\\ &\leqslant&\xi^2, \end{eqnarray} \tag{ A.4 }$$

where we have introduced the smallness parameter

$$\begin{equation} \xi^2\equiv\sum_{i=x,y,z}\langle\beta_{i}^2(0)\rangle\tau^2=\langle\boldsymbol{\beta}(0)\boldsymbol{\beta}^{\mathrm{T}}(0)\rangle\tau^2. \end{equation} \tag{ A.5 }$$

For Gaussian noise, correlation functions evaluated at an odd numbers of time points vanish, so 〈a₁a^T₂〉 = 0.

Moving on to the first of the three terms that involve four-point noise correlation functions

$$\begin{eqnarray} &&\fl \langle a_{2}^{2}\rangle=\sum_{i,j,i',j'} \int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2} \int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4} \langle \beta_{i}(t_{1})\beta_{j}(t_{2})\beta_{i'}(t_{3})\beta_{j'}(t_{4})\rangle\tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{i'j'}(t_{3},t_{4}).\nonumber\\ \end{eqnarray} \tag{ A.6 }$$

Here we can apply the Gaussian moment theorem to write

$$\begin{eqnarray} &&\fl \langle \beta_{i}(t_{1})\beta_{j}(t_{2})\beta_{i'}(t_{3})\beta_{j'}(t_{4})\rangle =\langle \beta_{i}(t_{1})\beta_{j}(t_{2})\rangle\langle\beta_{i'}(t_{3})\beta_{j'}(t_{4})\rangle+\langle \beta_{i}(t_{1})\beta_{i'}(t_{3})\rangle\langle\beta_{j}(t_{2})\beta_{j'}(t_{4})\rangle\nonumber\\ &&+\langle \beta_{i}(t_{1})\beta_{j'}(t_{4})\rangle\langle\beta_{i'}(t_{3})\beta_{j}(t_{2})\rangle. \end{eqnarray} \tag{ A.7 }$$

Using the independence of the noise components, this becomes

$$\begin{eqnarray} &&\fl \langle \beta_{i}(t_{1})\beta_{j}(t_{2})\beta_{i'}(t_{3})\beta_{j'}(t_{4})\rangle =\delta_{ij}\delta_{i'j'}\langle \beta_{i}(t_{1})\beta_{i}(t_{2})\rangle\langle\beta_{i'}(t_{3})\beta_{i'}(t_{4})\rangle\nonumber\\ &&+\delta_{ii'}\delta_{jj'}\langle \beta_{i}(t_{1})\beta_{i}(t_{3})\rangle\langle\beta_{j}(t_{2})\beta_{j}(t_{4})\rangle +\delta_{ij'}\delta_{i'j}\langle \beta_{i}(t_{1})\beta_{i}(t_{4})\rangle\langle\beta_{i'}(t_{3})\beta_{i'}(t_{2})\rangle, \end{eqnarray} \tag{ A.8 }$$

so that

$$\begin{eqnarray}&&\fl \langle a_{2}^{2}\rangle=\sum_{ij}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\left\{\langle \beta_{i}(t_{1})\beta_{i}(t_{2})\rangle \langle\beta_{j}(t_{3})\beta_{j}(t_{4})\rangle \tilde{\boldsymbol{R}}_{ii}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jj}(t_{3},t_{4})\right.\nonumber\\ &&\left.+\langle \beta_{i}(t_{1})\beta_{i}(t_{3})\rangle \langle\beta_{j}(t_{2})\beta_{j}(t_{4})\rangle \tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ij}(t_{3},t_{4})\right.\nonumber\\ &&\left.+\langle \beta_{i}(t_{1})\beta_{i}(t_{4})\rangle \langle\beta_{j}(t_{3})\beta_{j}(t_{2})\rangle \tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ji}(t_{3},t_{4})\right\} \end{eqnarray} \tag{ A.9 }$$

and

$$\begin{eqnarray}&&\fl |\langle a_{2}^{2}\rangle|\leqslant \sum_{ij}\langle\beta_{i}^2(0)\rangle\langle\beta_{j}^2(0)\rangle \left|\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\nonumber\right.\\ &&\left.\times\left\{\overline{\langle \beta_{i}(t_{1})\beta_{i}(t_{2})\rangle} \overline{\langle\beta_{j}(t_{3})\beta_{j}(t_{4})\rangle} \tilde{\boldsymbol{R}}_{ii}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jj}(t_{3},t_{4})\right.\right.\nonumber\\ &&\left.\left.+\overline{\langle \beta_{i}(t_{1})\beta_{i}(t_{3})\rangle} \overline{\langle\beta_{j}(t_{2})\beta_{j}(t_{4})\rangle} \tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ij}(t_{3},t_{4})\right.\right.\nonumber\\ &&\left.\left.+\overline{\langle \beta_{i}(t_{1})\beta_{i}(t_{4})\rangle} \overline{\langle\beta_{j}(t_{3})\beta_{j}(t_{2})\rangle} \tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ji}(t_{3},t_{4})\right\}\right|. \end{eqnarray} \tag{ A.10 }$$

Each of the three integrands in (A.10) has a magnitude less than or equal to unity, so

$$\begin{equation} |\langle a_{2}^{2}\rangle|\leqslant3\sum_{ij}\langle\beta_{i}^2(0)\rangle\langle\beta_{j}^2(0)\rangle\tau^4/4 =3\xi^4/4. \end{equation} \tag{ A.11 }$$

Now

$$\begin{equation} \langle \mathbf{a}_{1}\mathbf{a}^{\mathrm{T}}_{3}\rangle=\frac{2}{3}\sum_{i'ijk}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{\tau}\mathrm{d}t_{2} \int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\langle \beta_{i'}(t_1)\beta_{i}(t_2)\beta_{j}(t_3)\beta_{k}(t_4)\rangle, \end{equation} \tag{ A.12 }$$

which, using the Gaussian moment theorem and the independence of the noise components, we have

$$\begin{eqnarray} &&\fl \langle\mathbf{a}_{1}\mathbf{a}^{\mathrm{T}}_{3} \rangle=\frac{2}{3}\sum_{ij}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{\tau}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\left\{\langle \beta_{i}(t_{1})\beta_{i}(t_{2})\rangle \langle\beta_{j}(t_{3})\beta_{j}(t_{4})\rangle \boldsymbol{R}_{i}(t_{1})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ijj}(t_{2},t_{3},t_{4})\right.\nonumber\\ &&\left.+\langle \beta_{i}(t_{1})\beta_{i}(t_{3})\rangle \langle\beta_{j}(t_{2})\beta_{j}(t_{4})\rangle \boldsymbol{R}_{i}(t_{1})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jij}(t_{2},t_{3},t_{4})\right.\nonumber\\ &&\left.+\langle \beta_{i}(t_{1})\beta_{i}(t_{4})\rangle \langle\beta_{j}(t_{3})\beta_{j}(t_{2})\rangle \boldsymbol{R}_{i}(t_{1})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jji}(t_{2},t_{3},t_{4})\right\}. \end{eqnarray} \tag{ A.13 }$$

Following the same procedure used to derive (A.11), we find that

$$\begin{equation} |\langle\mathbf{a}_{1}\mathbf{a}^{\mathrm{T}}_{3}\rangle|\leqslant\xi^{4}/3. \end{equation} \tag{ A.14 }$$

Finally,

$$\begin{eqnarray}&&\fl \langle a_{1}^{4}\rangle=\sum_{i,j,i',j'} \int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{\tau}\mathrm{d}t_{2} \int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{\tau}\mathrm{d}t_{4} \langle \beta_{i}(t_{1})\beta_{j}(t_{2})\beta_{i'}(t_{3})\beta_{j'}(t_{4})\rangle\nonumber\\ &&\times\left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{2})\right)\left(\boldsymbol{R}_{i'}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{j'}(t_{4})\right), \end{eqnarray} \tag{ A.15 }$$

which becomes

$$\begin{eqnarray}&&\fl \langle a_{1}^{4}\rangle=\sum_{ij}\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{\tau}\mathrm{d}t_{2}\int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{\tau}\mathrm{d}t_{4}\nonumber\\ &&\times\left\{\langle \beta_{i}(t_{1})\beta_{i}(t_{2})\rangle \langle\beta_{j}(t_{3})\beta_{j}(t_{4})\rangle \left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{2})\right)\left(\boldsymbol{R}_{j}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{4})\right)\right.\nonumber\\ &&\left.+\langle \beta_{i}(t_{1})\beta_{i}(t_{3})\rangle \langle\beta_{j}(t_{2})\beta_{j}(t_{4})\rangle \left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{2})\right)\left(\boldsymbol{R}_{i}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{4})\right)\right.\nonumber\\ &&\left.+\langle \beta_{i}(t_{1})\beta_{i}(t_{4})\rangle \langle\beta_{j}(t_{3})\beta_{j}(t_{2})\rangle \left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{2})\right)\left(\boldsymbol{R}_{j}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{4})\right)\right\}. \end{eqnarray} \tag{ A.16 }$$

We can then show that

$$\begin{equation} |\langle a_{1}^{4}\rangle|\leqslant3\xi^{4}. \end{equation} \tag{ A.17 }$$

Obviously, for Gaussian noise, the overall trend here is that those terms containing n-point correlation functions make a maximum contribution of the order of ξⁿ to the fidelity. If we limit our attention to a regime in which ξ² ≪ 1, then higher-order terms will have little effect and the fidelity may be expressed in terms of the first-order error vector only.

The noise power spectral density S_i(ω) of the ith noise component may be defined by

$$\begin{equation} \langle\beta_{i}(t_{1})\beta_{i}(t_{2})\rangle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega S_{i}(\omega)\,\mathrm{e}^{\mathrm{i}\omega(t_{2}-t_{1})}. \end{equation} \tag{ A.18 }$$

Noting that

$$\begin{equation} \langle\beta_{i}^2(0)\rangle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega S_{i}(\omega)=\frac{1}{\pi}\int_{0}^{\infty}\mathrm{d}\omega S_{i}(\omega), \end{equation} \tag{ A.19 }$$

the smallness parameter ξ may be expressed in terms of the noise spectrum. In particular, if the ith noise component has a cutoff frequency ω_ci, then

$$\begin{equation} \langle\beta_{i}^2(0)\rangle=\frac{1}{\pi}\int_{0}^{\omega_{ci}}\mathrm{d}\omega S_{i}(\omega) \end{equation} \tag{ A.20 }$$

and

$$\begin{equation} \xi^2\equiv\frac{\tau^{2}}{\pi}\sum_{i=x,y,z}\int_{0}^{\omega_{ci}}\mathrm{d}\omega S_{i}(\omega). \end{equation} \tag{ A.21 }$$

The dependence of the fidelity on the spectral properties of the noise and of the control can also be made explicit using (A.18). Substituting (A.18) into (A.3), we have

$$\begin{equation} \langle a_{1}^{2}\rangle=\frac{1}{2\pi}\sum_{i=x,y,z}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{i}(\omega)F^{1}_{i}(\omega), \end{equation} \tag{ A.22 }$$

where

$$\begin{equation} F^{(1)}_{i}(\omega)\equiv\boldsymbol{R}_{i}(\omega)\boldsymbol{R}_{i}^{\mathrm{T}*}(\omega) \end{equation} \tag{ A.23 }$$

is the first-order generalized filter function for the ith noise component, defined in terms of the ith row of the frequency-domain control matrix

$$\begin{equation} \boldsymbol{R}_{i}(\omega)\equiv-\mathrm{i}\omega\int_{0}^{\tau}\mathrm{d}t\,\mathrm{e}^{\mathrm{i}\omega t}\boldsymbol{R}_{i}(t). \end{equation} \tag{ A.24 }$$

Similarly, from (A.9)

$$\begin{eqnarray} &&\fl \langle a_{2}^{2}\rangle=\frac{1}{(2\pi)^2}\sum_{ij}\int_{-\infty}^{\infty}\mathrm{d}\omega S_{i}(\omega) \int_{-\infty}^{\infty}\mathrm{d}\omega' S_{j}(\omega')\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\nonumber\\ &&\times\left\{\mathrm{e}^{\mathrm{i}\omega(t_{2}-t_{1})}\,\mathrm{e}^{\mathrm{i}\omega'(t_{4}-t_{3})}\tilde{\boldsymbol{R}}_{ii}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jj}(t_{3},t_{4}) +\mathrm{e}^{\mathrm{i}\omega(t_{3}-t_{1})}\,\mathrm{e}^{\mathrm{i}\omega'(t_{4}-t_{2})}\tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ij}(t_{3},t_{4})\right.\nonumber\\ &&\left.+\mathrm{e}^{\mathrm{i}\omega(t_{4}-t_{1})}\,\mathrm{e}^{\mathrm{i}\omega'(t_{3}-t_{2})}\tilde{\boldsymbol{R}}_{ij}(t_{1},t_{2})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ji}(t_{3},t_{4})\right\} \end{eqnarray} \tag{ A.25 }$$

while, from (A.13)

$$\begin{eqnarray} &&\fl \langle\mathbf{a}_{1}\mathbf{a}^{\mathrm{T}}_{3} \rangle=\frac{1}{(2\pi)^2}\sum_{ij}\int_{-\infty}^{\infty}\mathrm{d}\omega S_{i}(\omega) \int_{-\infty}^{\infty}\mathrm{d}\mathrm{}\omega' S_{j}(\omega')\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{\tau}\mathrm{d}t_{2}\int_{0}^{t_{2}}\mathrm{d}t_{3}\int_{0}^{t_{3}}\mathrm{d}t_{4}\nonumber\\ &&\times\left\{{\rm e}^{{\rm i}\omega(t_{2}-t_{1})}\,{\rm e}^{{\rm i}\omega'(t_{4}-t_{3})}\boldsymbol{R}_{i}(t_{1})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{ijj}(t_{2},t_{3},t_{4}) +{\rm e}^{{\rm i}\omega(t_{3}-t_{1})}\,{\rm e}^{{\rm i}\omega'(t_{4}-t_{2})}\boldsymbol{R}_{i}(t_{1})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jij}(t_{2},t_{3},t_{4})\right.\nonumber\\ &&\left.+\,{\rm e}^{{\rm i}\omega(t_{4}-t_{1})}\,{\rm e}^{{\rm i}\omega'(t_{3}-t_{2})}\boldsymbol{R}_{i}(t_{1})\tilde{\boldsymbol{R}}^{\mathrm{T}}_{jji}(t_{2},t_{3},t_{4})\right\} \end{eqnarray} \tag{ A.26 }$$

and from (A.16)

$$\begin{eqnarray}&&\fl \langle a_{1}^{4}\rangle=\frac{1}{(2\pi)^2}\sum_{ij}\int_{-\infty}^{\infty}\mathrm{d}\omega S_{i}(\omega) \int_{-\infty}^{\infty}\mathrm{d}\omega' S_{j}(\omega')\int_{0}^{\tau}\mathrm{d}t_{1}\int_{0}^{\tau}\mathrm{d}t_{2}\int_{0}^{\tau}\mathrm{d}t_{3}\int_{0}^{\tau}\mathrm{d}t_{4}\nonumber\\ &&\times \left\{{\rm e}^{{\rm i}\omega(t_{2}-t_{1})}\,{\rm e}^{{\rm i}\omega'(t_{4}-t_{3})}\left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{2})\right)\left(\boldsymbol{R}_{j}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{4})\right)\right. \end{eqnarray} \tag{ A.27 }$$

$$\begin{eqnarray} & &\left.+{\rm e}^{{\rm i}\omega(t_{3}-t_{1})}\,{\rm e}^{{\rm i}\omega'(t_{4}-t_{2})}\left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{2})\right)\left(\boldsymbol{R}_{i}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{4})\right)\right.\nonumber\\ & &\left.+{\rm e}^{{\rm i}\omega(t_{4}-t_{1})}\,{\rm e}^{{\rm i}\omega'(t_{3}-t_{2})}\left(\boldsymbol{R}_{i}(t_{1})\boldsymbol{R}^{\mathrm{T}}_{j}(t_{2})\right)\left(\boldsymbol{R}_{j}(t_{3})\boldsymbol{R}^{\mathrm{T}}_{i}(t_{4})\right)\right\}. \end{eqnarray} \tag{ A.28 }$$

Extrapolating these results, we find the following alternative form for the series expansion of the fidelity

$$\begin{eqnarray}&&\fl \mathcal{F}_{\rm av}=1-\left[\frac{1}{2\pi}\sum_{i}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{i}(\omega)F^{(1)}_{i}(\omega)\nonumber\right.\\ &&\left.+\frac{1}{(2\pi)^{2}}\sum_{ij}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{\omega^2}S_{i}(\omega)\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega'}{\omega'^2} S_{j}(\omega')F^{(2)}_{ij}(\omega,\omega')+\cdots\right], \end{eqnarray} \tag{ A.29 }$$

where the generalized filter functions F^(p)_i1i2...(ω,ω',...) are derived solely from the control matrix.