General solution to inhomogeneous dephasing and smooth pulse dynamical decoupling

We focus on the case of a driven qubit described by the following Hamiltonian

$$\begin{eqnarray}&&{ \mathcal H }(t)=\displaystyle \frac{{\rm{\Omega }}(t)}{2}{\sigma }_{z}+\delta \beta {\sigma }_{x}.\end{eqnarray} \tag{ 1 }$$

Here, Ω(t) is the control field, while δβ is an unknown stochastic fluctuation. ${ \mathcal H }(t)$ is one of the most widely used models of decoherence due to slow noise [46–48]. We can view it as describing a qubit with a tunable energy splitting given by Ω(t), with δβ interpreted as a fluctuating transverse field. Alternatively, we can think of Ω(t) as the amplitude of a resonant, monochromatic driving field, in which case δβ represents a fluctuation in the detuning caused by fluctuations in the qubit energy splitting which give rise to inhomogeneous dephasing. In the second case, ${ \mathcal H }(t)$ is the Hamiltonian in a frame which rotates with the driving field (in the rotating wave approximation). This Hamiltonian could also describe two driven levels within a larger Hilbert space. For example, the dynamics of three-level lambda systems can often be mapped to an effective driven two-level system [49–51].

We parametrize the evolution operator according to

$$\begin{eqnarray}U(t)=\left(\begin{array}{cc}{u}_{1}(t) & -{u}_{2}^{* }(t)\\ {u}_{2}(t) & {u}_{1}^{* }(t)\end{array}\right).\end{eqnarray} \tag{ 2 }$$

It is generally not possible to solve the Schrödinger equation, ${\rm{i}}\dot{U}(t)={ \mathcal H }(t)U(t)$, to obtain u₁ and u₂ analytically aside from a few special choices of Ω(t). These special choices include the cases of δ-function pulses and square pulses, which is the primary reason why most dynamical decoupling sequences are based on these ideal waveforms. There exist methods for solving this equation more generally using a type of reverse-engineering [42, 52, 53], but these techniques do not apply for the form of ${ \mathcal H }(t)$ we are considering here. Although we cannot obtain closed-form solutions in general, we can obtain solutions as a power series in δβ [39]

$$\begin{eqnarray}{u}_{1}(t) & = & {{\rm{e}}}^{-{\rm{i}}\phi (t)/2}({g}_{0}(t)-{g}_{2}(t)\delta {\beta }^{2}+{g}_{4}(t)\delta {\beta }^{4}-...),\\ {u}_{2}(t) & = & -{{\rm{ie}}}^{{\rm{i}}\phi (t)/2}({g}_{1}^{* }(t)\delta \beta -{g}_{3}^{* }(t)\delta {\beta }^{3}+...),\end{eqnarray} \tag{ 3 }$$

where $\phi (t)={\int }_{0}^{t}{\rm{\Omega }}(\tau ){\rm{d}}\tau$ is the qubit rotation angle, and the coefficients of δβ obey a recurrence relation:

$$\begin{eqnarray}&&{g}_{n}(t)={\int }_{0}^{t}{{\rm{e}}}^{{\rm{i}}\phi (\tau )}{g}_{n-1}^{* }(\tau ){\rm{d}}\tau ,\end{eqnarray} \tag{ 4 }$$

with g₀(t) = 1. To obtain robust quantum operations with operation time T, we must find functions ϕ(t) (from which we obtain the driving field via ${\rm{\Omega }}(t)=\dot{\phi }$) which yield vanishing error coefficients, g_m(T) = 0, for all m ≤ n, where n is the desired degree of robustness. Since g_m is complex, this gives two real constraints on Ω(t) for each odd m. For even m, there is only one real constraint because unitarity requires $| {u}_{1}{| }^{2}+| {u}_{2}{| }^{2}=1$, which from equation (3) implies that at any time t,

$$\begin{eqnarray}&&\mathrm{Re}[{g}_{n}]=\displaystyle \sum _{k=1}^{n/2-1}{(-1)}^{k-1}\mathrm{Re}[{g}_{k}{g}_{n-k}^{* }]-\displaystyle \frac{{(-1)}^{n/2}}{2}| {g}_{n/2}{| }^{2}.\end{eqnarray} \tag{ 5 }$$

Thus for even $n,\mathrm{Re}[{g}_{n}(T)]=0$ is automatically satisfied if g_m(T) = 0 for all m ≤ n/2, and we are left with one real independent constraint at this order: $\mathrm{Im}[{g}_{n}(T)]=0$. For a sequence of δ-function π pulses, ${{\rm{e}}}^{{\rm{i}}\phi (\tau )}=\pm 1$ is the associated characteristic function [46], and it is easy to show that g_n(T) = 0 for spin echo, CPMG, and other well known sequences. In what follows, we will show that there exists a geometrical framework that allows us to find all possible driving fields that satisfy these constraints order by order.

We begin by considering the first-order constraint, g₁(T) = 0. The solution space of this constraint can be constructed by first parametrizing this function in terms of plane Cartesian coordinates:

$$\begin{eqnarray}&&{g}_{1}(t)=x(t)+{\rm{i}}{y}(t).\end{eqnarray} \tag{ 6 }$$

Since the time derivative of g₁(t) is a pure phase (equation (4)), we find the following restriction between x(t) and y(t)

$$\begin{eqnarray}&&{\dot{x}}^{2}+{\dot{y}}^{2}=1.\end{eqnarray} \tag{ 7 }$$

The error coefficient, g₁(t), can thus be viewed as a curve in the plane spanned by x and y which starts at the origin at time t = 0: x(0) = y(0) = 0, and with slope ${\dot{x}}^{2}(0)+{\dot{y}}^{2}(0)=1$ (see figure 1). There is a one-to-one correspondence between these plane curves and possible driving fields which follows from plugging equation (6) into the left-hand-side of equation (4) and differentiating both sides twice:

$$\begin{eqnarray}&&\Omega(t)=\dot x\ddot y-\dot y\ddot x=\frac{\dot x\ddot y-\dot y\ddot x}{(\dot x^2+\dot y^2)^{3/2}}.\end{eqnarray} \tag{ 8 }$$

The above formula for Ω(t) has a surprisingly simple geometric interpretation: it is precisely equal to the signed curvature, κ, of the plane curve, which is defined to be the proportionality coefficient between the derivative of the curve's tangent vector and its normal vector: $\dot{{\boldsymbol{V}}}(t)=\kappa (t){\boldsymbol{N}}(t)$. Since the curvature is defined independently of parametrization, we may express it in terms of a general parameter λ along the curve

$$\begin{eqnarray}&&\kappa (\lambda )=\displaystyle \frac{x^{\prime} y^{\prime\prime} -y^{\prime} x^{\prime\prime} }{{(x{{\prime} }^{2}+y{{\prime} }^{2})}^{3/2}},\end{eqnarray} \tag{ 9 }$$

where the prime denotes differentiation with respect to λ. We are thus free to use any parametrization we choose to define the plane curve. Once we have made this choice, we can extract the driving field: Ω(t) = κ(λ(t)), where the mapping between λ and t is determined from equation (7):

$$\begin{eqnarray}&&t={\int }_{0}^{\lambda }{\rm{d}}\mu \sqrt{{[x^{\prime} (\mu )]}^{2}+{[y^{\prime} (\mu )]}^{2}}.\end{eqnarray} \tag{ 10 }$$

Notice that the integral on the right hand side measures distance along the curve, meaning that time t can be thought of as the arc-length parametrization of the curve. The driving field is thus the curvature of the plane curve expressed as a function of arc length (see figure 1). This also means that the evolution time T is equal to the total length of the curve. Also notice that the curvature, and hence the driving field, remain invariant under rigid rotations and translations of the curve. Thus, although ${\dot{g}}_{1}(0)={{\rm{e}}}^{{\rm{i}}\phi (0)}=1$ would seem to imply that the initial part of the curve must be tangential to the x axis, it is actually not necessary to enforce this property.

Zoom In Zoom Out Reset image size

Requiring the first-order error to vanish at t = T is equivalent to imposing the simple condition that the curve is closed: x(Λ) = y(Λ) = 0, where Λ ≡ λ(T). Also, the net rotation angle, ϕ(T), is determined by the angle, Δθ, subtended by the curve at the origin

$$\begin{eqnarray}&&\phi (T)={\int }_{0}^{T}{\rm{d}}{t}{\rm{\Omega }}(t)={\left.\arctan \left(\displaystyle \frac{y^{\prime} }{x^{\prime} }\right)\right|}_{0}^{{\rm{\Lambda }}}={\rm{\Delta }}\theta +\pi .\end{eqnarray} \tag{ 11 }$$

This is illustrated in figure 2(a) for an example curve with Δθ = π/2. The corresponding pulse, which implements a z-rotation by angle ϕ(T) = 3π/2, is shown in figure 2(b). The above formula for ϕ(T) is valid modulo 2π because of the multi-valued nature of the inverse tangent. Here, it is important to emphasize that any closed curve that starts and ends at the origin and which subtends the same angle, Δθ, at the origin yields a pulse that implements the same qubit rotation while canceling the leading-order error. Conversely, any pulse which accomplishes these tasks corresponds to a plane curve of this form. Note that in this framework, a single spin echo pulse is represented by the limiting case in which the closed curve collapses on itself to form a straight line segment extending away from the origin a finite distance r, as shown in figures 3(a), (b); the δ-function π pulse corresponds to the sharp turning point at r. More general δ-function pulse sequences also correspond to collapsed line segments, with one turning point per pulse (see figures 3(c), (d)). Square pulses correspond to circular arcs, while a full circle yields a square pulse that implements first-order dynamical decoupling (see figures 3(e), (f)). This geometrical interpretation not only makes it easy to visualize the complete solution space of the first-order constraint, it also makes it easy to systematically determine optimal pulses for a given rotation. This is because the driving field is proportional to the curvature, providing a simple way to visualize and to impose constraints on the smoothness of the pulse. In addition, the fact that the evolution time is equal to the length of the curve allows for optimization of the operation time given smoothness constraints.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This mapping between robust pulses and closed planar curves is reminiscent of a similar picture in which dynamical decoupling can be interpreted in terms of closed vector paths in the Lie algebra associated with the Hamiltonian [27]. Unlike the present work, however, there is no known method to find the smooth pulses associated with these closed paths, and one is typically restricted to considering piecewise-constant pulses that perform robust identity operations. Interestingly, closed-loop conditions are also known to arise for other types of Hamiltonians, for example ones describing coupled qubit-oscillator dynamics [54, 55].

We can construct simple, explicit examples of smooth pulses, which cancel the first-order error by modifying the well known lemniscates (figure-eight curves) of Bernoulli:

$$\begin{eqnarray}&&x(\lambda )=\displaystyle \frac{\alpha \sin (2\lambda )}{3+\cos (2\lambda )},\qquad y(\lambda )=\displaystyle \frac{2\sin \lambda }{3+\cos (2\lambda )},\end{eqnarray} \tag{ 12 }$$

and of Gerono:

$$\begin{eqnarray}&&x(\lambda )=(\alpha /2)\sin (2\lambda ),\qquad y(\lambda )=\sin \lambda .\end{eqnarray} \tag{ 13 }$$

In both cases, we choose the maximal value of λ to be Λ = π in order to retain only half the lemniscate, giving a closed curve that subtends an angle Δθ < π at the origin. We have included the free parameter α in order to make this angle, and hence the rotation angle ϕ(T), adjustable: Δθ = π − 2arccot(α). The example shown in figures 2(a), (b) corresponds to a modified Bernoulli curve with α = 1. Figures 2(c)–(f) shows additional examples of half-lemniscates and the corresponding pulses (obtained from the curvature, equation (9)) taken from both the Bernoulli and Gerono families of curves. Together, these solutions give, for each choice of rotation angle ϕ(T), two different pulses which implement the same rotation while canceling the leading-order error. Of course, all possible smooth pulses which achieve the same task can be generated from curves qualitatively similar to those shown in figures 2(c), (e). Notice that if we keep the pulse time fixed, then implementing larger rotation angles requires sharper peaks in the pulse, a fact which is clear from the geometrical construction since the plane curve is forced to have sharper bends as Δθ increases.

Because of the recursive nature of the error constraints evident in equation (4), requiring higher-order errors to cancel can be interpreted as imposing additional constraints on the shape of the plane curve. Consider the second-order constraint, g₂(T) = 0. From equation (5) we know that $\mathrm{Re}[{g}_{2}(T)]=0$ if the first-order error vanishes, g₁(T) = 0, i.e., if the plane curve is closed. The remaining second-order error can be written as

$$\begin{eqnarray}&&\mathrm{Im}[{g}_{2}(T)]={\int }_{0}^{{\rm{\Lambda }}}{\rm{d}}\lambda [x^{\prime} (\lambda )y(\lambda )-y^{\prime} (\lambda )x(\lambda )].\end{eqnarray} \tag{ 14 }$$

Remarkably, this integral is exactly twice the area enclosed by the plane curve, meaning that the second-order error vanishes if and only if the net area enclosed by the curve is zero. Note that the sign of the area is determined by the direction of the winding of the curve, so that the clockwise and counter-clockwise parts of the curve enclose areas of opposite sign. Thus we find that pulses which cancel both first- and second-order errors are in one-to-one correspondence with closed plane curves that enclose a vanishing net area. Full leminiscates such as those of Gerono and Bernoulli, equations (12), (13), satisfy these requirements, provided we now take Λ = 2π. The resulting curves and corresponding pulses are shown in figures 4(a), (b). Notice that Ω(t) < 0 in the second half of the evolution as follows from the negative curvature of one of the lemniscate lobes. This behavior can be interpreted as a sudden π phase shift in the driving field. An important difference compared to the first-order case is that now, each of these pulses implements dynamical decoupling, i.e., an identity operation with ϕ(T) = 0. This is because the angle subtended by the initial and final legs of the plane curves is equal to π regardless of the value of α. We see from the figure that it is easy to obtain arbitrarily many smooth pulses that perform dynamical decoupling up to second order.

Zoom In Zoom Out Reset image size

In order to perform nontrivial rotations while canceling errors up to second order, we need to consider more general zero-area curves that exhibit a kink at the origin. We introduce one family of such curves, parametrized as

$$\begin{eqnarray}x(\lambda ) & = & [a+b\cos (\lambda )]\sin (\lambda ),\\ y(\lambda ) & = & \lambda (2\pi -\lambda )+(4+b/a)[\cos (\lambda )-1],\end{eqnarray} \tag{ 15 }$$

with Λ = 2π, and where the rotation angle is given by

$$\begin{eqnarray}&&\phi (T)=2\arctan \left(\displaystyle \frac{a+b}{2\pi }\right)+\pi .\end{eqnarray} \tag{ 16 }$$

Restricting attention to pulses that start and end at zero, Ω(0) = Ω(T) = 0, requires the parameters a and b to be related through either b = −a or b = −6a. The former condition produces 3π rotations while the latter allows for a wide range of rotation angles. Examples of curves with b = −6a and their corresponding pulses are shown in figures 4(c), (d). Figure 5(a) shows examples of curves with b = −a, with the associated pulses given in figure 5(b). Since our cancellation constraints, equation (4), are derived assuming quasistatic noise in which δβ is an unknown but constant variation in the Hamiltonian, it is important to consider the performance of the resulting pulses in the presence of time-dependent noise fluctuations. Perhaps the most ubiquitous type of classical time-dependent noise in the context of quantum technologies is 1/f noise [12, 15, 56, 57]. Figure 5(c) shows the operation fidelity [58] for the pulse corresponding to a = 2/3 = −b (blue curve in figure 5(b)) as a function of noise strength. For comparison, the figure also shows the performance of a square pulse with the same duration. It is clear that even for time-dependent noise, pulses derived from equation (4) outperform naive pulses over a wide range of noise strengths if the noise is sufficiently slow. Further details regarding the numerical simulation can be found in the Methods section.

Zoom In Zoom Out Reset image size

To further illustrate the simplicity of our geometrical approach to dynamical decoupling and dynamically corrected gates, we have written a program called DDdraw [59] that allows the user to draw plane curves by hand using a computer mouse. DDdraw displays the corresponding pulse, curve area (second-order cancellation constraint), and rotation angle in real time, making it easier to tailor the pulse as desired and to ensure that the errors are canceled up to second order. Figure 6 shows an example of a second-order dynamical decoupling pulse generated from a hand-drawn plane curve (figure 6(a)) using DDdraw.

Zoom In Zoom Out Reset image size

The advantage of smooth pulses over square or δ-function pulses comes from experimental limitations in pulse generation. Microwave pulses produced from arbitrary waveform generators cannot vary arbitrarily quickly; typical rise time limits are on the scale of 100 ps–1 ns. This leads to unavoidable errors when square or δ-function pulses are attempted in experiments, and these errors can be appreciable for pulse times on the order of nanoseconds. Limitations on pulse amplitude can come from physical constraints on the system Hamiltonian (e.g., voltage-controlled exchange pulses are constrained by the upper bound on the strength of the exchange interaction set by the confinement energy [60]) or from the need to avoid overheating the system (as is the case for e.g., magnetically driven electron spins in quantum dots [3]). This leads to further errors in the case of δ-function dynamical decoupling since the pulses may be leveled off at relatively low amplitudes because of these restrictions.

To demonstrate the significance of these effects and the advantage afforded by smooth pulses, we compare our smooth pulses with square and δ-function decoupling pulses subject to both amplitude and rise time restrictions. We first consider the case of pure dynamical decoupling where the goal is to perform an identity operation while canceling noise errors to second order. Ideally, this can be achieved with many different δ-function sequences, but for concreteness we focus on 2-pulse CPMG (which cancels errors to arbitrary order in the case of quasistatic noise). We account for amplitude and rise time restrictions by considering the following 'pseudo' δ-function sequence made from hyperbolic secant pulses:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{CPMG}}}(t)=\displaystyle \frac{\gamma }{T}[{\rm{sech}} (\gamma (t-T/4))-{\rm{sech}} (\gamma (t-3T/4))],\end{eqnarray} \tag{ 17 }$$

where we choose γ = 47.2 so that the area of each pulse is π, and T is the total evolution time. This 2-pulse sequence is shown in figure 7(a). Here, we have chosen the rise time to be 2% of the total evolution time. In this figure, we also show one of our smooth pulses, the Bernoulli lemniscate with α = 1 shown in figure 4(b), which is also designed to perform an identity operation while canceling second-order errors. In figure 7(b) we compare the infidelities for these two types of dynamical decoupling as a function of transverse noise strength. Here, the noise is taken to be quasistatic Gaussian noise, and the strength is characterized by the width σ_β of the Gaussian distribution. We find that the infidelity is reduced by several orders of magnitude over a broad range of noise strengths when the smooth pulse is used.

Zoom In Zoom Out Reset image size

Next, we perform a similar comparison, but now for a non-trivial rotation. Ideally, such a rotation can be performed while canceling noise errors using a composite square pulse such as the one shown in figure 8(a). This composite pulse is designed to cancel errors to second order while implementing a ϕ = −π/2 rotation about z. We impose rise time restrictions by deforming this composite pulse into a smoothened, pseudo-square version (red curve in the figure). We again take the rise time to be 2% of the total evolution time for concreteness. The precise form of this composite pulse is

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{comp}}}(t)=f(t)\displaystyle \sum _{i=1}^{5}{a}_{i}[1-\tanh (b(-{l}_{i}+t-{t}_{0i}))][\tanh (b({l}_{i}+t-{t}_{0i}))+1],\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&f(t)=\displaystyle \frac{\tanh (4{bt})[\tanh (4b(t-T))-1]+\tanh (4b(t-T))+\tanh (4{bT})}{\tanh (4{bT})-\tanh (2{bT})[2+\tanh (2{bT})]},\end{eqnarray} \tag{ 19 }$$

and the parameters ${t}_{0i}=({x}_{i}+{x}_{i-1})/2$ and $2{l}_{i}={x}_{i}-{x}_{i-1}$ give respectively the center and the duration of the ith square pulse, which starts and ends at t = x_i−1 and t = x_i, respectively. Here, we choose the parameters x/T = (0, 0.268, 0.299, 0.701, 0.732, 1), aT = (1.422, 3.357, −3.357, 3.357, 1.422), and bT = 110.581, which yields the pulse shown in figure 8(a). This figure also shows one of our smooth error-suppressing pulses which implements the same rotation and cancels noise to second order while obeying the same rise time restriction. This pulse is obtained from the family of curves introduced in equation (15) by choosing a = 2π/5 and b = −6a. The advantage of this pulse over the composite one is evident in figure 8(b), which shows that the infidelity for the smooth pulse is orders of magnitude smaller than that of the composite pulse for a wide range of noise strengths. In quantum computing and related fields, where infidelities typically need to be suppressed to the 10⁻³ level or below, this advantage can become very important.

Zoom In Zoom Out Reset image size

The relationship between the evolution time and the distance along the curve has an additional geometrical interpretation: the plane curve can be viewed as a lightlike (or null) curve in 3D Minkowski space. This fact can be seen by rewriting equation (7) as −dt² + dx² + dy² = 0, which is the infinitesimal proper distance along a lightlike curve. In this interpretation, any Ω(t) which cancels the first-order error maps onto a lightlike curve which starts and ends at the origin of 3D Minkowski space. Such a curve is necessarily confined within a lightcone with apex at x = y = t = 0, as shown in figure 9. This lightcone is a reflection of the quantum speed limit, which refers to the minimum time it takes to evolve from one quantum state to another [53, 61–65]. For a two-level system, the quantum speed limit is saturated for zero driving, Ω(t) = 0, which corresponds to a straight line in the plane, or equivalently a lightlike curve lying on the 3D lightcone. This is consistent with our earlier observation that δ-function pulse sequences correspond to line segments in the plane: in the 3D Minkowski picture, the evolution between δ-pulses is represented by lightlike curves (straight lines parallel to the lightcone), with a sharp bend for each pulse. This is illustrated in figure 9 for the case of spin echo.

Zoom In Zoom Out Reset image size

The higher-order error-cancellation constraints give further restrictions on the shape of the plane curves. In the Methods section, we show that the constraints at arbitrary order n are built from two basic classes of integrals:

$$\begin{eqnarray}&&{p}_{k}(\lambda )={\int }_{0}^{\lambda }{\rm{d}}\mu {r}^{k}\theta ^{\prime} ,\quad {q}_{k}(\lambda )={\int }_{0}^{\lambda }{\rm{d}}\mu {r}^{k}\theta ^{\prime} {{\rm{e}}}^{{\rm{i}}\theta },\end{eqnarray} \tag{ 20 }$$

where r and θ are plane polar coordinates: re^iθ = x + iy, and k is an integer. For example, third-order cancellation introduces the new constraint, q₃(Λ) = 0, while for fourth-order cancellation, we need p₄(Λ) = 0. These can be understood as requiring that certain signed volumes bounded by the plane curves must vanish. Consider for instance the constraint p_k(Λ) = 0. Since the area enclosed by the plane curve can be written as $\tfrac{1}{2}{\int }_{0}^{{\rm{\Lambda }}}{\rm{d}}\mu {r}^{2}(\mu )\theta ^{\prime} (\mu )$, we recognize that p_k(Λ) is proportional to the volume lying between the plane curve area and the function F(r) = r^k−2 over the plane. The constraints q_k(Λ) = 0 admit a similar interpretation. In the Methods section, we show that the following family of deformed Gerono lemniscates parametrized by even integer m,

$$\begin{eqnarray}&&x(\lambda )=\sin \lambda ,\qquad y(\lambda )=\displaystyle \frac{1}{2}\sin (m\lambda )\sin \lambda \end{eqnarray} \tag{ 21 }$$

satisfy the error constraints up to order m − 1. These curves and their associated pulses are shown in figure 10, where it is evident that higher-order errors are canceled by inserting additional self-intersections in the curve. For example, the curve shown in figure 10(a) yields a pulse (figure 10(e)) that cancels only the first-order error error, while the curve in figure 10(b), which contains two additional self-intersections, gives a pulse (figure 10(f)) that cancels up through third order. The self-intersections of the curves allow integral constraints of the form p_2n(T) = 0 to be satisfied since the lobes of the curve on either side of the intersection then have areas of opposite sign. Each additional lobe that appears in figures 10(a)–(d) has precisely the area need to cancel p_2n(T) for all n < m/2. Integrals of the type ${q}_{2n+1}(T)$ for n < m/2 cancel due to the inversion symmetry of the curves (see Methods). It is apparent from figure 10 that for fixed pulse time, the pulses become sharply peaked as higher-order errors are canceled. This is a consequence of needing to insert more turns and self-intersections in the plane curves in order to cancel higher-order errors, and it suggests that an all-orders error cancellation can only be achieved with δ-function pulses. This is consistent with a no-go theorem [39] stating that all-orders error cancellation is not possible with finite pulses implemented within a finite time.

Zoom In Zoom Out Reset image size

The above findings also imply that if we keep the maximal pulse amplitude Ω_max fixed (and hence the maximal allowed curvature of the plane curve), then the pulse duration T must increase as the order of the error cancellation is increased. This in turn implies that the pulse bandwidth would decrease since it behaves like 1/T. For the particular class of error-canceling pulses shown in figures 10(e)–(h), we have calculated Ω_maxT for the first 6 pulses (corresponding to m = 2, 4, 6, 8, 10, 12) in the series (the last of which cancels up to 11th order errors) and found that the maximal bandwidths are $1/T={{\rm{\Omega }}}_{\max }\times [0.02,0.002,3.6\times {10}^{-4},9.6\times {10}^{-5},3.3\times {10}^{-5},1.4\times {10}^{-5}]$. This data can be fit to the function $\exp (-{{am}}^{b})$ where a = 2.41 and b = 0.68, which shows that the maximal bandwidth decreases subexponentially with increasing error-cancellation order.

Here, we provide details regarding our time-dependent 1/f noise simulations. We simulate 1/f noise as a sum over many random telegraph noise (RTN) sources:

$$\begin{eqnarray}&&\delta \beta (t)=\displaystyle \sum _{i=0}^{N}{w}_{i}{\eta }_{i}(t),\end{eqnarray} \tag{ 22 }$$

where η_i(t) represents a stochastic noise fluctuation for a single RTN source, N + 1 is the total number of RTN sources, and w_i is a weight factor that we must determine. The above expansion implies that

$$\begin{eqnarray}\langle \delta \beta (t)\delta \beta (0)\rangle & = & \displaystyle \sum _{i=0}^{N}\displaystyle \sum _{j=0}^{N}{w}_{i}{w}_{j}\langle {\eta }_{i}(t){\eta }_{j}(0)\rangle \\ & = & \displaystyle \sum _{i=0}^{N}\langle {\eta }_{i}(t){\eta }_{i}(0)\rangle {w}_{i}^{2}=\displaystyle \sum _{i=0}^{N}{{\rm{e}}}^{-2| t| /{\tau }_{i}}{w}_{i}^{2},\end{eqnarray} \tag{ 23 }$$

where in the final step, we used that different RTN sources are uncorrelated and that the two-point correlation function of the ith RTN source decays exponentially with characteristic timescale defined as τ_i. Taking the Fourier transform of both sides, we find that the noise power spectrum is given by

$$\begin{eqnarray}&&S(\omega )=\displaystyle \sum _{i=0}^{N}\displaystyle \frac{4{\tau }_{i}}{4+{\omega }^{2}{\tau }_{i}^{2}}{w}_{i}^{2}.\end{eqnarray} \tag{ 24 }$$

We want to determine the RTN parameters such that this spectrum approximates that of 1/f noise:

$$\begin{eqnarray}&&S(\omega )\approx \displaystyle \frac{{A}^{2}}{\omega }.\end{eqnarray} \tag{ 25 }$$

More precisely, this relation should hold in the continuum limit, $N\to \infty$. This will be the case if we choose the τ_i to be equally spaced:

$$\begin{eqnarray}&&{\tau }_{i}=\displaystyle \frac{{\rm{i}}}{N}({\tau }_{\max }-{\tau }_{\min })+{\tau }_{\min },\end{eqnarray} \tag{ 26 }$$

and if we choose the weight factor to be

$$\begin{eqnarray}&&{w}_{i}^{2}=\displaystyle \frac{{A}^{2}}{N\pi }\displaystyle \frac{{\tau }_{\max }-{\tau }_{\min }}{{\tau }_{i}}.\end{eqnarray} \tag{ 27 }$$

Equations (22), (26) and (27) are what we use to approximate 1/f noise using an ensemble of RTN sources in numerical simulations. For the numerical results shown in the main text, we take $N={10}^{5},{\tau }_{\min }=44T,{\tau }_{\max }=441T$, and vary A so that the time-averaged noise strength, $\langle | \delta \beta (t)| \rangle$, varies from 0 to 9.2/T, averaging over 50 instances of the noise for each value. The square pulse used in the comparison is taken to have the same duration, T, as our smooth error-correcting pulse, as well as the same area of 3π. We use the definition of fidelity given in [58].

Starting from the recursion relation, equation (4), and performing integrations by parts multiple times yields one of four expressions depending on the value of n mod 4:

$$\begin{eqnarray}{g}_{n} & = & {g}_{n-1}^{* }{g}_{1}-{g}_{n-2}{g}_{2}^{* }+\ldots -{g}_{(n+2)/2}{g}_{(n-2)/2}^{* }\\ & & +\displaystyle \int {{\rm{d}}{t}{g}}_{n/2}^{* }{\dot{g}}_{n/2},\qquad n=4m+2,\\ {g}_{n} & = & {g}_{n-1}^{* }{g}_{1}-{g}_{n-2}{g}_{2}^{* }+\ldots +{g}_{(n+2)/2}^{* }{g}_{(n-2)/2}\\ & & -\displaystyle \int {{\rm{d}}{t}{g}}_{n/2}{\dot{g}}_{n/2}^{* },\qquad n=4m+4,\\ {g}_{n} & = & {g}_{n-1}^{* }{g}_{1}-{g}_{n-2}{g}_{2}^{* }+\ldots +{g}_{(n+3)/2}^{* }{g}_{(n-1)/2}\\ & & -\displaystyle \int {{\rm{d}}{t}{g}}_{(n+1)/2}{\dot{g}}_{(n-1)/2}^{* },\qquad n=4m+1,\\ {g}_{n} & = & {g}_{n-1}^{* }{g}_{1}-{g}_{n-2}{g}_{2}^{* }+\ldots -{g}_{(n+3)/2}{g}_{(n-1)/2}^{* }\\ & & +\displaystyle \int {{\rm{d}}{t}{g}}_{(n+1)/2}^{* }{\dot{g}}_{(n-1)/2},\qquad n=4m+3,\end{eqnarray} \tag{ 28 }$$

where m is a non-negative integer. We see that if we assume g_k(T) = 0 for k < n, then the constraint g_n(T) = 0 reduces to a condition involving two lower-order error coefficients:

$$\begin{eqnarray}{\displaystyle \int }_{0}^{T}{{\rm{d}}{t}{g}}_{n/2}^{* }{\dot{g}}_{n/2} & = & 0,\qquad n\ \mathrm{even},\\ {\displaystyle \int }_{0}^{T}{{\rm{d}}{t}{g}}_{(n+1)/2}^{* }{\dot{g}}_{(n-1)/2} & = & 0,\qquad n\ \mathrm{odd}.\end{eqnarray} \tag{ 29 }$$

Taking n = 2 and working with plane polar coordinates, g₁ = re^iθ, we obtain

$$\begin{eqnarray}&&{g}_{2}(t)={\int }_{0}^{t}{\rm{d}}\tau {g}_{1}^{* }{\dot{g}}_{1}=\displaystyle \frac{1}{2}r{(t)}^{2}+{\rm{i}}{\int }_{0}^{t}{\rm{d}}\tau {r}^{2}\dot{\theta }.\end{eqnarray} \tag{ 30 }$$

We recognize the integral on the right as the area enclosed by the plane curve when t = T. Using this result for g₂ in equation (29) and assuming r(T) = 0, we find that the third-order constraint simplifies to

$$\begin{eqnarray}&&{g}_{3}(T)=\displaystyle \frac{4{\rm{i}}}{3}{\int }_{0}^{T}{{\rm{d}}{t}{r}}^{3}{{\rm{e}}}^{{\rm{i}}\theta }\dot{\theta }=0.\end{eqnarray} \tag{ 31 }$$

Similarly, if we also assume g₂(T) = 0, then cancellation of the fourth-order constraint requires

$$\begin{eqnarray}&&{g}_{4}(T)={\rm{i}}{\int }_{0}^{T}{{\rm{d}}{t}{r}}^{4}\dot{\theta }=0,\end{eqnarray} \tag{ 32 }$$

while for fifth order we have

$$\begin{eqnarray}&&{g}_{5}(T)=\displaystyle \frac{8{\rm{i}}}{15}{\int }_{0}^{T}{{\rm{d}}{t}{r}}^{5}{{\rm{e}}}^{{\rm{i}}\theta }\dot{\theta }-\displaystyle \frac{8}{3}{\int }_{0}^{T}{{\rm{d}}{t}{r}}^{2}\dot{\theta }{\int }_{0}^{t}{\rm{d}}{t}^{\prime} {r}^{3}{{\rm{e}}}^{{\rm{i}}\theta }\dot{\theta },\end{eqnarray} \tag{ 33 }$$

and for sixth order,

$$\begin{eqnarray}&&{g}_{6}(T)=\displaystyle \frac{2{\rm{i}}}{9}{\int }_{0}^{T}{{\rm{d}}{t}{r}}^{6}\dot{\theta }+\displaystyle \frac{16}{9}{\int }_{0}^{T}{{\rm{d}}{t}{r}}^{3}{{\rm{e}}}^{{\rm{i}}\theta }\dot{\theta }{\int }_{0}^{t}{\rm{d}}{t}^{\prime} {r}^{3}{{\rm{e}}}^{-{\rm{i}}\theta }\dot{\theta }.\end{eqnarray} \tag{ 34 }$$

Notice that g₄(T) and g₆(T) are both purely imaginary (the latter assumes g₃(T) = 0) as expected from equation (5). The fact that all these constraints consist of integrals of ${\dot{p}}_{k}(t),{\dot{q}}_{k}(t)$ (which are defined in equation (20)) and their complex conjugates arises because the functions g_n(t) are themselves comprised of such terms This in turn can be proven using induction in conjunction with equation (4).

In this section, we show that the deformed Gerono lemniscate of order m given in equation (21) solves the error-cancellation constraints up to order $m-1$. First, we show that integrals of the type p_2n(T/2) vanish for m > 2n, where p_k(t) is defined in equation (20). To do this, we begin by rewriting this integral:

$$\begin{eqnarray}{p}_{2n}(T/2) & \propto & {\displaystyle \int }_{0}^{\pi }{\sin }^{2n}(\lambda ){({\sin }^{2}(m\lambda )+4)}^{n-1}\,{\rm{d}}\sin (m\lambda )\\ & = & \displaystyle \sum _{i=0}^{n-1}{a}_{i}{\displaystyle \int }_{0}^{\pi }{\sin }^{2n}(\lambda ){\sin }^{2i}(m\lambda )\,{\rm{d}}\sin (m\lambda )\\ & = & \displaystyle \sum _{i=0}^{n-1}\displaystyle \frac{{a}_{i}}{2i+1}{\displaystyle \int }_{0}^{\pi }{\sin }^{2n}(\lambda ){\rm{d}}{\sin }^{2i+1}(m\lambda )\\ & = & \displaystyle \sum _{i=0}^{n-1}{\tilde{a}}_{i}{\displaystyle \int }_{0}^{\pi }{\sin }^{2i+1}(m\lambda ){\sin }^{2n-1}(\lambda )\cos (\lambda ){\rm{d}}\lambda \\ & = & \displaystyle \sum _{i=0}^{n-1}\displaystyle \sum _{j=0}^{i}{\tilde{a}}_{i}{b}_{j}{\displaystyle \int }_{0}^{\pi }\sin ((2j+1)m\lambda ){\sin }^{2n-1}(\lambda )\cos (\lambda ){\rm{d}}\lambda .\end{eqnarray} \tag{ 35 }$$

In this formula, each term has the form

$$\begin{eqnarray}{I}_{2n+1} & = & {\displaystyle \int }_{0}^{\pi }\sin (2{\ell }\lambda ){\sin }^{2n+1}(\lambda )\cos (\lambda ){\rm{d}}\lambda \\ & = & -{\ell }{\displaystyle \int }_{0}^{\pi }\cos (2{\ell }\lambda ){\sin }^{2n+2}(\lambda ){\rm{d}}\lambda -{{nI}}_{2n+1}.\end{eqnarray} \tag{ 36 }$$

Thus we get ${I}_{2n+1}=-\tfrac{{\ell }}{n+1}{J}_{2n+2}$, where we have defined

$$\begin{eqnarray}&&{J}_{2n}={\int }_{0}^{\pi }\cos (2{\ell }\lambda ){\sin }^{2n}(\lambda ){\rm{d}}\lambda .\end{eqnarray} \tag{ 37 }$$

Using integration by parts, we can also obtain

$$\begin{eqnarray}&&{J}_{2n+2}=\displaystyle \frac{-{\ell }}{n+1}{I}_{2n+1}+\displaystyle \frac{2n+1}{2n+2}{J}_{2n},\end{eqnarray} \tag{ 38 }$$

so we have in the end

$$\begin{eqnarray}&&{I}_{2n+1}=\displaystyle \frac{-(2n+1){\ell }}{2(n+1+{\ell })(n+1-{\ell })}{J}_{2n}.\end{eqnarray} \tag{ 39 }$$

We see from equation (37) that J₀ = 0, so from the relation above it follows that for n + 1 < ℓ we have ${I}_{2n+1}=0$. If we combine this conclusion with the formula in equation (35), we find that as long as $m\gt 2n$ the deformed Gerono lemniscate, equation (21), is able to cancel ${p}_{2n}(T/2)$. The cancellation of these integrals implies that ${p}_{2n}(T)$ also vanishes.

The deformed Gerono lemniscates also cancel the integrals ${q}_{2n+1}(T)$, where q_k(t) is defined in equation (20). This follows from the curves' inversion symmetry. Every point, λ, along the curve has an inversion-related partner, $\tilde{\lambda }$, such that $\theta (\tilde{\lambda })=\pi +\theta (\lambda ),r(\tilde{\lambda })=r(\lambda )$, and $\theta ^{\prime} (\tilde{\lambda })=\theta ^{\prime} (\lambda )$. This implies that the contribution to ${q}_{2n+1}(T)$ from the x < 0 half of the curve exactly cancels that from the x > 0 half.

In addition to terms like p_n(T) or q_n(T), the nth order error coefficient also contains terms with nested integrals involving p_k(t) and q_k(t) (see e.g., equations (33) and (34)). We do not have a proof that such terms also vanish identically for the deformed Gerono lemniscates, but we have verified this numerically up to seventh order and conjecture that this trend continues to arbitrary order.