Robust quantum control using smooth pulses and topological winding

The greatest challenge in achieving the high level of control needed for future technologies based on coherent quantum systems is the decoherence induced by the environment. Here, we present an analytical approach that yields explicit constraints on the driving field which are necessary and sufficient to ensure that the leading-order noise-induced errors in a qubit's evolution cancel exactly. We derive constraints for two of the most common types of noise that arise in qubits: slow fluctuations of the qubit energy splitting and fluctuations in the driving field itself. By theoretically recasting a phase in the qubit's wavefunction as a topological winding number, we can satisfy the noise-cancelation conditions by adjusting driving field parameters without altering the target state or quantum evolution. We demonstrate our method by constructing robust quantum gates for two types of spin qubit: phosphorous donors in silicon and nitrogen-vacancy centers in diamond.

In this work, we are only interested in the case of single-axis driving, for which we have ϕ = 0 and β is a constant. In this case, the Hamiltonian reduces to the form where we have defined b z (t) = Ω(t). This Hamiltonian is in a form that is relevant for several types of experimentally relevant two-level quantum systems, for example singlet-triplet spin qubits in semiconductor quantum dots. [2] In many types of qubits, particularly those used in electron spin resonance or nuclear magnetic resonance experiments, the effective qubit Hamiltonian can be expressed as where E is the qubit energy splitting and Ω(t)e iωt represents a monochromatic pulse with envelope Ω and frequency ω.
This rotating-frame Hamiltonian has the form of Eq. (5) with β = (ω − E)/2. In the context of NV centers or electron spin qubits in phosphorous donors in silicon, fluctuations in β result from fluctuations in E, which in turn are caused by Overhauser noise arising from nuclear spins in the environment. In both these contexts, as well as in the case of nuclear spin qubits in phosphorous donors, noise in Ω(t) would result from e.g., power fluctuations in the pulse generator. Such fluctuations were identified as a possible cause of nuclear spin control infidelities. [3] Although H rot has the same form as Eq. (5), it is important to note that we want the evolution operator in the lab frame (i.e., the evolution operator associated with H qubit above) to be the identity matrix at t = 0. The evolution operators corresponding to H qubit and H rot are related by so that our initial condition for U rot should then be U rot (0) = T † (0) = e i π 4 σy = 1 √ 2 One way to incorporate this new initial condition would be to simply multiply U rot on the right by T † (0). This will re-arrange the components of the evolution operator shown in Eq. (3) but does not modify the error cancelation procedure we will develop using the Hamiltonian given in Eq. (5).

II. FLUCTUATIONS IN β
In this section, we want to consider the situation where the qubit splitting (or external field detuning) β exhibits stochastic fluctuations that are constant throughout the duration of the pulse, i.e., β = β 0 + δβ, where β 0 is a known constant and δβ is an unknown constant. We would like to engineer pulses that execute a target evolution while minimizing the errors caused by the unknown variation δβ. In order to set up this problem, it turns out to be beneficial to use a new parametrization of the driving field and the corresponding evolution operator. We first derive this new parametrization and then return to the question of how to incorporate stochastic fluctuations into our formalism for generating analytical solutions of two-level quantum dynamics.
The first step in deriving the new parametrization is to transform to a rotating frame in the x-basis: The functions D ± then solve the following set of equations which follow from the Schrödinger equation for the evolution operator U :Ḋ These two equations can be solved easily for Ω(t): which implies that we may writeḊ where u = u(t, β) is an unknown complex function. These two equations are easily solved: where we have imposed the initial conditions D ± (0) = 1/ √ 2. Equations (12) also implyḊ Plugging Eqs. (15) into this equation yields Differentiating both sides of this equation with respect to time and performing some algebraic manipulations results in Writing we have It is straightforward to express D ± in terms of w: Any choice of w(t) such that the Ω(t) computed from Eq. (20) is real immediately yields an analytical solution to the Schrödinger equation. It is straightforward to find the necessary restriction on w. Writing w = w r + iw i in Eq. (20), the condition Im[Ω] = 0 impliesẇ This equation is easily integrated to give In this expression, c is an integration constant determined by the boundary condition w r (0) = 1 2 log tan c, which we will leave arbitrary, at least for now. We must impose to ensure that w r is real. Plugging Eq. (23) into Eq. (20), we find an expression for Ω in terms of w i : This result can in fact be obtained from the method given in Ref. [4] by choosing So far, we have derived an alternate algorithm for generating analytical solutions to the Schrödinger equation: An analytical solution can be produced by choosing any w i (t) such that 0 ≤ 2β t 0 dt sin(2w i (t )) + 2c ≤ π. We now consider the case where β = β 0 + δβ, where β 0 is a known constant and δβ is a small stochastic variation. Since our formula for Ω(t), Eq. (20), contains an explicit dependence on β, we must take care to choose w(t) in such a way that the resulting Ω(t) obeys Ω(t, β) = Ω(t, β 0 ) + O(δβ 2 ). This condition is necessary since Ω(t) cannot itself depend on the stochastic variable δβ. As long as we are only concerned with first-order variations in δβ, it is sufficient to require the first-order variation of Ω(t) to vanish. Note that the present analysis could be extended to derive constraints which ensure that higher-order errors in the evolution operator vanish by requiring the corresponding higher-order variations of Ω(t) to vanish.
We can arrange for the first-order variation of Ω(t) to vanish by writing w = w 0 + δβw 1 and varying Eq. (20) to obtain The first-order variation will vanish if w 1 is given by where k is an integration constant. What this result means is that we may still use our formalism for constructing analytical solutions even in the presence of the stochastic fluctuation δβ. However, unless we arrange for higher-order variations in Ω(t) to vanish as well, we can only obtain the evolution up to first order in δβ. In most physical situations, δβ β 0 , and determining the evolution up to first order in δβ is sufficient. The procedure works by first choosing w 0,i to fix the desired driving field, Ω(t). Choosing w 0,i also fixes w 0,r and w 0 , and the latter then determines w 1 . The corresponding evolution operator to first order in δβ is obtained by using w = w 0 + δβw 1 in Eq. (21).
We have just seen that the w-parametrization is useful for removing the δβ dependence from Ω(t). However, now that we have laid out the necessary steps, it turns out that further simplifications occur if we return to the χ-parametrization. In particular, we may simplify the expression for w 1 . To see this, begin with the general relation between the two parameterizations: Since we wish to impose χ(0) = 0, we will set c = 0 from now on. Using this expression for χ(t), it is straightforward to show so that where Comparing this result for cosh(2w 0 ) with the definition of the function ξ(t) defined earlier, we see that where here ξ 0 (t) is defined in terms of χ 0 (t) (and not the full χ(t)): Substituting this result into Eq. (28) and setting the integration constant k therein to zero then yields The first-order variation, δw = δβw 1 , leads to the following variation in χ(t): This completes the derivation of Eq. (3a) of the main text. Below, we will show how to use this final expression for δ β χ to construct dynamically corrected qubit operations with smooth analytical pulses.

III. FLUCTUATIONS IN Ω(t)
To see how we can include first-order fluctuations of the form we again start from the w-parametrization from Eq. (20): Writing w(t) = w 0 (t) + δ w 1 (t), expanding the right hand side to first order in δ and equating all the first-order terms, we find This equation is easily solved, with the result where k is an integration constant. Using Eq. (33), we can rewrite this as In order to translate this result into a variation δ χ(t), we use Eq. (29) which yields Using that and integrating by parts, we have This completes the derivation of Eq. (3b) of the main text.

IV. CONTROLLED ROTATIONS VIA TOPOLOGY
The χ formalism reviewed above allows us to construct arbitrary rotations by choosing χ(t) appropriately. However, fixing the target evolution precisely remains challenging because of the fact that the phase ξ(t) = (ξ + +ξ − )/2 appearing in the evolution operator (3) is given as an integral of a complicated nonlinear expression involving χ(t): Even though we get to choose the form of χ(t) with only the relatively weak constraint |χ| ≤ |β| to worry about, it is difficult to predict what values of ξ f = ξ(t f ) at the final time t f can be achieved with a given choice, and ultimately this generally requires numerical evaluation of the integral above. Moreover, if we wish to include tunable parameters in χ(t) in order to cancel errors and/or improve rotational fidelities, this task is hindered by the complicated nature of the above expression for ξ(t) since it will be hard to keep ξ f fixed as such parameters are varied. In this section, we will show how these problems can be circumvented by incorporating the concept of topological winding into our formalism.
The key observation is to notice that the invariance of ξ f under parameter variations in χ(t) is equivalent to the statement that ξ f is a quantized topological winding number. This quantization can be imposed by choosing χ(t) in such a way that the integrand in ξ is proportional to the derivative of the argument of a complex function W (t): where λ is a real constant that we are free to choose. If W (t) winds around the origin of the complex W plane an integral number (n) of times as time evolves from t = 0 to t = t f , then The value of ξ f will then be quantized (in units of 2πλ) according to This formula shows that we can fix ξ f to be any value we like by choosing λ appropriately and by picking a W (t) which exhibits nontrivial winding. Including tunable parameters in χ(t) is tantamount to including them in W (t); adjusting these parameters will deform the contour traced by W (t) but will leave the winding number n the same so long as the contour does not cross the origin in the process. Once W (t) is chosen, we then solve Eq. (46) to obtain χ(t), from which the corresponding driving field, Ω(t), follows from Eq. (4). The procedure outlined in the previous paragraph allows us to hold fixed ξ f (and hence the rotation axis and angle) while we adjust parameters in χ(t). However, it introduces a new difficulty in that the step of solving Eq. (46) can be challenging due to its strongly nonlinear nature. We can simplify this task by expressing W (t) as a function of χ(t): This allows us to rewrite Eq. (46) asχ Since the term in curly brackets is now purely a function of χ, we may solve this equation by simple integration: Instead of solving a nonlinear differential equation, the task has now been reduced to performing an ordinary integral and inverting the result. To ensure that ξ f is quantized, we now need to choose Y (χ) to be such that it winds an integral number of times around the origin of the complex Y plane as χ evolves from 0 to its final value, χ f . We can further streamline the determination of b z (t) by expressing this quantity in terms of χ(t) using Eq. (50). Explicitly, we find where We can also express ξ(t) as a function of χ: In practice, it turns out to be easier to obtain simpler-looking functions b z (t) if we work directly with the function rather than with Y (χ). The driving field is fully determined by Φ(χ): To obtain the driving field as a function of time, we need to solve the following equation for χ (see Eq. (51)): The integrand on the right-hand-side can be identified as the line element obtained from the following metric for a 2-sphere: where we view χ as the azimuthal angle and Φ as the polar angle. This in turn implies that the duration of the pulse, is just the length of the curve defined by Φ(χ) on the surface of the 2-sphere (times 2/β).

V. DERIVATION OF GENERAL ERROR-CANCELATION CONSTRAINTS
We are now ready to combine all of the ingredients developed in the previous sections to systematically construct rotations that dynamically correct for errors induced by fluctuations in either β or Ω(t). In order to obtain nontrivial results, it is important that we use the unperturbed quantities χ 0 (t), β 0 , and ξ 0 (t) in the topological winding formalism developed above.

A. β noise
From the general form of the evolution operator, Eq. (3), we see that in the case of a general single-axis driving field, the first-order error due to δβ-noise will vanish only if the following three conditions on the first-order variations hold: Combining the topological winding formalism developed in the previous section with our expression for δ β χ(t) in the case of fluctuations in β, Eq. (36), we have We will make this expression slightly more compact by defining the function yielding Introducing the shorthand notation δχ f = δχ(t f ), etc., the end-point fluctuations are then given by and imposing the first two error-cancelation conditions from Eq. (60) amounts to requiring which is equivalent to the single complex condition Notice that in these conditions, we have omitted the possibility that Φ (χ f ) = 0, which would have automatically guaranteed that δχ f = (χ f /β 0 )δβ, meaning that we could ignore the second condition in Eq. (65). The reason for this is that when Φ (χ f ) = 0 (which impliesχ f = β 0 ), we do not completely cancel the first-order error in the final evolution operator because of a square root appearing in its first-order variation stemming from the arcsine part of ξ ± : Ifχ f = β 0 , then the square root in the denominator vanishes, implying that this part of the first-order variation in the evolution will not vanish in general. The only way to guarantee that the first-order error in the evolution operator is completely removed is to not imposeχ f = β 0 . It remains to analyze the third condition in Eq. (60). To calculate δ β ξ f , we begin with the general expression for ξ(t): The square root inside the integrand can be expanded as follows: while the remaining factor becomes We therefore have In the case of β-noise, we substitute the expressions for δ β χ(t) and δ βχ (t) from Eq. (63) into Eq. (71) to obtain Using the relation β 2 0 −χ 2 0 =χ 0 Φ (χ 0 ) sin(2χ 0 ), we can eliminate all of the square roots: With the help of the relationχ we can write each term as an integral over χ: This result can be simplified dramatically if we perform an integration by parts on the last term: The second term here precisely cancels a similar-looking term in Eq. (75), leaving us with where we have also used that The requirement that δ β ξ f = 0 then yields the following condition: However, the second term in this condition is already required to vanish according to the additional constraints we have derived earlier (see Eq. (65)), so the first term must also vanish independently: In summary, the following two constraints are necessary and sufficient for ensuring that the leading-order error in the evolution operator due to β-noise vanishes identically: From the general form of the evolution operator, Eq. (3), we see that in the case of a general single-axis driving field, the first-order error due to Ω-noise will vanish only if the following three conditions on the first-order variations hold: In Sec. III, we found the fluctations in χ(t) caused by fluctuations in Ω(t), Eq. (44). Using the results of Sec. IV, this can be rewritten as where we have definedg such that g(t) =g(χ 0 (t)). Defining the function and setting k = 0 as is appropriate for fluctuations that originate from the control field Ω(t), this becomes As in the case of β-noise, we again want to avoid setting Φ (χ f ) = 0 since this will lead to an incomplete cancelation of the first-order error. We are then left with the following condition necessary for canceling the error: We now turn to the third condition from Eq. (82). For driving field noise, we may again use the first-order variation of ξ(t), Eq. (71), but now drop the term involving δβ: Using the expressions for δ χ(t) and δ χ(t) given earlier in Eq. (85), this becomes (88) We next remove the square roots and turn the integrations over t into integrations over χ using Integrating the second term by parts, we find that the third term here cancels the first term in Eq. (89), leaving The end-point variation is then When we impose the earlier condition from Eq. (86), the second term vanishes, in turn requiring that the first term also vanishes independently. In summary, the following two constraints are necessary and sufficient for ensuring that the leading-order error in the evolution operator due to Ω-noise vanishes identically: Here, we collect all the preceding results and summarize the general approach to constructing robust control fields. These control fields are completely determined by a function Φ(χ) which we are free to choose. We can make the initial evolution operator the identity matrix by imposing Φ(0) = 0 and Φ (0) = 0. Choosing η = −1 in Eq. (3) and focusing on the case of single-axis driving, we find that the target evolution operator is determined by the value of Φ and its first derivative at a chosen final value of χ, which we denote by χ f : A driving field Ω(t) which generates this target evolution via the Hamiltonian H = Ω(t)σ z + βσ x can be found by performing the following integration inverting the result to obtain χ 0 (t), and plugging this into The most general way to ensure that the resulting pulse has finite duration is to require that the numerator of the expression for Ω(χ), Eq. (96), vanishes at χ = χ f . We should think of this condition as a condition on the final value of the second derivative of Φ(χ): Viewing the condition this way is appropriate since we have already fixed χ f and Φ (χ f ) in accordance with our desired U target . We can further ensure that the resulting pulse amplitude has zero slope at the beginning and end of the pulse by imposing that the third derivative of Φ(χ) vanishes at χ = 0 and χ = χ f . This will force the pulse to have flattened tails, which in turn can reduce sensitivity to timing errors. If we want the leading-order errors in the evolution to vanish, then we must constrain not only the boundary conditions of Φ(χ) but its full behavior from χ = 0 to χ = χ f as well. In the case of β-noise, Φ(χ) must satisfy whereas in the case of Ω(t)-noise, it must satisfy where the functiong(χ) is determined by the precise nature of the noise one wishes to consider: δΩ =g(χ 0 (t))δ . For example, pulse amplitude fluctuations are described by choosingg(χ) = Ω(χ). If one wishes to cancel errors due to both types of noise simultaneously, then Φ(χ) must satisfy both sets of constraints. The fact that the arbitrary time-dependence in the Ω-noise fluctuation carries all the way through the calculation and leads to the appearance of the arbitrary functiong(χ) in Eq. (99) suggests that the analysis could be extended to the case of non-static noise. In particular, one could consider performing a Gaussian path integral average ofg, weighted by a nontrivial noise power spectrum. Sinceg appears linearly inside a single integration in Eq. (99), it is likely that this average could be performed analytically using standard Gaussian integration techniques. This analysis will be carried out elsewhere.
A straightforward way to solve the above constraints on Φ(χ), either Eqs. (98) or Eqs. (99), is to start with an ansatz for Φ(χ) of the form where a 1 and a 2 are constants and f is such that f (0) = f (0) = f (χ f ) = f (χ f ) = 0. We can set any desired U target by choosing χ f , a 1 , and a 2 appropriately. We may then include additional parameters in f (χ) and tune these until the error-cancelation constraints are satisfied. For example, if we choose then we may freely tune a 3 and a 4 without changing U target . This ansatz is used to construct the robust driving fields shown in Fig. 1 of the main text for the case of a driving field subject to noise in its amplitude. Interestingly, it is also possible to solve the β-noise cancelation constraints analytically for a certain set of target rotations. In particular, we can solve the first constraint in Eq. (98) when χ f = nπ/4 for some integer n by choosing Φ(χ) to have the general form where ζ(θ + nπ/2) = ζ(θ) is any periodic function with period nπ/2, and λ is an arbitrary real parameter. We can automatically solve the second constraint in Eq. (98) for any non-constant ζ(θ) by setting λ to This solution gives rise to a class of target evolutions determined by n and λζ (nπ); in particular, the phases in the evolution operator are If n is an even integer but not a multiple of 4, then U target reduces to a π-rotation about the x axis. If n is a multiple of 4, then U target is proportional to the identity. In this case, the control field implements a smooth-pulse version of dynamical decoupling. The fact that the β-noise constraints can be solved analytically greatly simplifies the problem in the case where one wishes to cancel both types of noise simultaneously. To do so, one could start with a more general ansatz for Φ: where all the ζ k (θ) are periodic functions of period nπ/2. This ansatz automatically solves the β-noise constraints provided the λ k satisfy one linear constraint coming from the second condition in Eq. (98). The remaining λ k can then be tuned until the Ω-noise constraints are also satisfied.

VII. ANTISYMMETRIC PULSES
If χ is an odd function, χ(−t) = −χ(t), which implies that Ω(−t) = −Ω(t), then the total evolution operator which evolves the qubit from t = −t f to t = t f takes a particularly simple form which leads to fewer error-cancelation constraints that need to be satisfied. The antisymmetry of χ(t) immediately implies that the phases in the evolution operator enjoy the following property under time reversal: Using this relation, it can be shown that the full evolution is described by the operator where U (t f ) describes the evolution from t = 0 to t = t f . Eq. (107) leads to the expressions U tot,11 = cos(2χ f ), This evolution operator corresponds to rotations about axes in the xy plane. The axis and angle of rotation can be expressed in terms of χ f and Φ (χ f ) according to where θ is the angle between the axis of rotation (which lies in the xy plane) and the x-axis, and φ is the angle of rotation. Thus finding a Φ(χ) that implements a desired evolution amounts to imposing the following condition on the derivative of Φ(χ) at χ f = φ/4: Because U tot does not depend on ξ(t f ) in the case of an antisymmetric pulse, we do not need to impose the errorcancelation constraint associated with fluctuations in ξ(t f ). This means that we may ignore the second constraints in Eqs. (98) and (99), and we need only solve the first constraints in each set to obtain robust antisymmetric control fields.

VIII. UNIVERSAL SET OF ROTATIONS
To perform an arbitrary single qubit rotation, it suffices to be able to rotate by an arbitrary angle around one chosen axis, plus being able to do a π rotation around an axis π/4 (45 degrees) apart from that chosen axis. The latter gate serves a similar role as the Hadamard gate which turns any rotation around one axis into a rotation with the same angle around an axis perpendicular to it.
In this section we present the parameters required to perform a set of rotations around an axis in the x-y plane lying at an angle θ from the x-axis, namely R(θ, φ) = exp −i(sin θσ y + cos θσ x ) φ 2 .
For each of the two cases discussed in the main text (pulses correcting driving field errors and δβ errors), we give pulses that implement rotations R(5π/12, π) and R(π/6, φ) for a range of φ. The two rotation axes are π/4 apart as required. Due to the strong nonlinearity of the problem, it is not practical to scan over the entire [0, 2π) range of φ continuously; we therefore demonstrate a set of rotations with angles between π and 3π with step 0.2π. Since our method is completely general and can take as many parameters as needed, it is straightforward to find the pulses for any single-qubit rotation.
Parameters for the set of rotations discussed above are provided in Table I.

IX. CONSTRUCTION OF ORDINARY PULSES AND DEFINITION OF FIDELITY
In Figs. 2 and 3 of the main text, we verify that the first-order error in the evolution operator has been canceled by comparing the infidelity of the designed pulse with that of a square pulse sequence that implements the same target evolution but which has not been designed to combat errors. In the case of single-axis driving, we can systematically construct a square pulse sequence that generates any target evolution operator using the following general form involving two identical square pulses: where R β (Ω; t) ≡ e −it(Ωσz+βσx) and τ = π/(2 √ 2β 0 ). Given U target , we can solve this equation numerically to obtain the parameters t a , t b , t c . In Figs. 2c and 3c of the main text, we use the definition of infidelity as in Ref. [5]: Here, U (t f ) is the actual evolution operator including the errors to all orders. In the case of β-noise for example, the evolution corresponding to the uncorrected pulses is given simply by while U (t f ) for the designed pulses is computed by solving the Schrödinger equation numerically with Hamiltonian H = Ω(t)σ z + (β 0 + δβ)σ x for each value of the error strength, δβ.