Optimal noise-canceling shortcuts to adiabaticity: application to noisy Majorana-based gates

Adiabatic braiding of Majorana zero modes can be used for topologically protected quantum information processing. While extremely robust to many environmental perturbations, these systems are vulnerable to noise with high-frequency components. Ironically, slower processes needed for adiabaticity allow more noise-induced excitations to accumulate, resulting in an antiadiabatic behavior that limits the precision of Majorana gates if some noise is present. In a recent publication [Phys. Rev. B 96, 075158 (2017)], fast optimal protocols were proposed as a shortcut for implementing the same unitary operation as the adiabatic braiding. These shortcuts sacrifice topological protection in the absence of noise but provide performance gains and remarkable robustness to noise due to the shorter evolution time. Nevertheless, gates optimized for vanishing noise are suboptimal in the presence of noise. If we know the noise strength beforehand, can we design protocols optimized for the existing unavoidable noise, which will effectively correct the noise-induced errors? We address this question in the present paper. We find such optimal protocols using simulated-annealing Monte Carlo simulations. The numerically found pulse shapes, which we fully characterize, are in agreement with Pontryagin's minimum principle, which allows us to arbitrarily improve the approximate numerical results (due to discretization and imperfect convergence) and obtain numerically exact optimal protocols. The protocols are \textit{bang-bang} (sequence of sudden quenches) for vanishing noise, but contain continuous segments in the presence of multiplicative noise due to the nonlinearity of the master equation governing the evolution. We find that such noise-optimized protocols completely eliminate the above-mentioned antiadiabatic behavior.

In recent years, significant progress has been made in the emerging area of shortcuts to adiabaticity [29][30][31], namely, performing dynamical quantum operations which yield the same outcome as the adiabatic evolution but in a shorter time. One of the approaches to finding such shortcuts is through optimal control [32][33][34][35][36][37][38][39][40][41][42][43]. An optimally fast scheme was derived in the absence of noise in Ref. [44], which produces the same unitary operation as adiabatic (infinite-time) braiding of MZMs in a short time. While sacrificing strict topological protection, this scheme offers practical advantages in balancing performance and robustness. The advantages are especially important when the control fields are noisy.
Noise, which may contain high-frequency components, is generically unavoidable and not protected against by topology. Furthermore, noise-induced errors accumulate over time and become more harmful in the adiabatic limit. Hence, the fast protocols of Ref. [44] suffer from much fewer noise-induced errors than slow adiabatic protocols that take a long time. In fact, noise results in an antiadiabatic behavior, where slowing down the dynamics can actually increase the diabatic excitations [45].
Nevertheless, the protocols of Ref. [44] are designed assuming vanishing noise, and suffer from a similar noise-induced heating rate characteristic of adiabatic protocols. This heating is limited since the fast protocols expose the system to noise for a shorter time than adiabatic protocols. In this paper, we obtain protocols that are optimal for an a priori known and unavoidable noise. They are optimized to cancel out the effects of the existing noise. Finding these noise-canceling optimal shortcuts helps us determine the precision limit of a Majorana-based gate in the presence of noise, as the protocols represent the best we can do given a certain total time, an existing noise strength, and an experimentally relevant range for the tunable inter-Majorana hybridization energies.
Throughout the paper, we use multiplicative noise, which is dominant for topological qubits. We focus on white noise for simplicity. A ubiquitous source of noise in superconducting devices, namely, the pink 1/ f noise, has a slowly decaying noise spectrum. It was shown in Ref. [44] that pink and white noise have a qualitatively similar antiadiabatic effect on topological gates. The hybridization energies ∆ j are tunable. (b) Adiabatic exchange of initially decoupled γ 1 and γ 2 , i.e., ∆ 1 (0) = ∆ 2 (0) = 0, can be achieved by, e.g., the linear protocol above in the limit of τ → ∞.
We address several open questions: What is the minimum error we can get if the protocols are optimized for the existing strength of noise? Is exact preparation of the unitary corresponding to the adiabatic evolution possible in the presence of noise? Will the antiadiabatic effect still occur when using noise-canceling optimal protocols? To answer these questions, we focus on a widely used minimal effective model of MZM braiding [46][47][48], whose evolution is governed by a Lindblad master equation [45,[49][50][51]. We find optimal protocols by numerical minimization (through simulated-annealing Monte Carlo simulations) of a cost function based on the trace distance of the noise-averaged final density matrix from the perfectly adiabatic target density matrix.
We also perform an in-depth analysis of the mathematical optimal-control theory of our system. In agreement with Pontryagin's minimum principle [52], we find that our protocols are bang-bang for vanishing noise (and close to bang-bang for weak noise). Upon increasing the noise strength, protocols with non-bang-bang segments occur due to the nonlinearities of multiplicative noise, still in full agreement with Pontryagin's principle. We fully characterize these optimal protocols. The connection to Pontryagin's principle is an important check that the Monte Carlo simulations have indeed found the optimal protocol. It can also help eliminate the discretization errors and imperfect convergence in the simulations, yielding numerically exact optimal protocols. Furthermore, in the regime of weak noise and short times, where Monte Carlo has divergence issues, a much more efficient bang-bang ansatz is used, which can yield very high-performance protocols.
Our main results are as follows. As a function of total time, we find three regimes in the behavior of the minimum cost function (obtained from optimal protocols customized for each total time and noise strength). In the first regime, the cost function does not change from its initial value, implying there is no permissible evolution with a total time smaller than a critical time τ c that can reduce the cost function. In the second regime, the cost function rapidly decreases as we increase the total time (in the absence of noise, this rapid decent continues all the way to zero), and in the third regime, the cost function exhibits a slow decay for finite noise (for vanishing noise we have a plateau with a vanishing cost function). Extrapolating the cost function of the third regime to the infinite-time limit, we find the cost function monotonically decreases with the total time τ and eventually vanishes as τ → ∞. Antiadiabatic behavior is completely eliminated by using these noise-canceling optimal protocols. Other significant contributions of the paper are the full characterization of the optimal pulse shapes. The patterns that emerge in the finite-time studies can inform numerical simulations for longer total times, where the computations become more expensive.
The outline of the paper is as follows. In Sec. II, we describe the effective model of the system and the unitary operator corresponding to perfect adiabatic braiding. We also present an error function, suitable for noisy systems, for the deviations from adiabaticity. In Sec. III, we discuss the Lindblad master equation used for the evolution of the noise-averaged density matrix in the presence of multiplicative (relevant to Majorana-based gates) white noise. We also discuss the antiadiabatic behavior of the system when standard protocols are used [see Fig. 1(b)]. In Sec. IV, we discuss the Monte Carlo and bang-bang optimization methods, and present our numerical results for the minimum error function as a function of noise strength and total time, as well as its extrapolation to the infinite-time limit. We find three distinct regimes in the behavior of the cost function. Sec. V is focused on the application of the mathematical theory of optimal control and Pontryagin's minimum principle to our system. In terms of the time-dependent density matrix and certain conjugate momenta, we derive the analytical form of the optimal protocols. In Sec. VI, we verify our numerically found optimal protocols against the analytical results of Sec. V, and characterize the patterns of the optimal protocols in different regimes. We close the paper in Sec. VII with brief conclusions. Some details of the application of Pontryagin's minimum principle are presented in two appendices.

II. MODEL OF A MAJORANA-BASED GATE AND ITS ERROR FUNCTION
We use an effective model of Majorana braiding relevant to devices such as the top-transmon [46][47][48]. The Hamiltonian can be written in terms of four Majorana operators as where γ j = γ † j and {γ i , γ j } = 2δ i j . Three Majorana operators can hybridize with one central operator γ 0 as in Fig. 1(a). Each ∆ j represents the hybridization energy between γ 0 and γ j for j = 1, 2, 3. We assume we have time-dependent control over these three Hamiltonian parameters and can tune each of them to an arbitrary value in the range 0 ≤ ∆ j (t) ≤ D. Hereafter, we set the maximum hybridization energy scale D = 1.
By defining two Dirac fermions c = (γ 1 + iγ 2 )/2 and d = (γ 0 + iγ 3 )/2, we can write the Hamiltonian as a 4 × 4 matrix in the basis |0 , d † |1 , |1 , d † |0 , where |1 ≡ c † |0 and |0 is the vacuum. The first (last) two basis states have even (odd) fermion parity, which is a symmetry of the Hamiltonian. We then represent H(t) by the block-diagonal matrix with where σ x,y,z are the Pauli matrices and 1 1 is the 2 × 2 identity matrix. The top (bottom) block corresponds to even (odd) fermion parity.
As illustrated in Fig. 1(b), for the case of turning the coupling constants on and off as linear functions of time, each step (e.g., turning ∆ 3 off and ∆ 1 on in the first step) is performed adiabatically in the standard scheme (e.g., by increasing τ). To generate the same gate in an optimal manner, we minimize an appropriate cost function (which vanishes for the above adiabatic process in the limit τ → ∞) for an arbitrary cyclic process of fixed total time τ with All three hybridization parameters can be arbitrary functions of time with only the following constraint on the range: The ideal gate transforms an initial superposition c 0 |0 + c 1 |1 into e iφ (c 0 |0 + ic 1 |1 ) at time t = τ, where φ is an arbitrary overall phase. As we are interested in a noisy system, a cost function defined in terms of noise-averaged quantities would be most helpful. Furthermore, we would like our cost function to be independent of φ. Density matrices are insensitive to the overall phase of the wavefunction and convenient for averaging over noise, making them a useful tool for characterizing gate precision. As in Ref. [44], we pick a particular equal-weight superposition |ψ 0 = 1 √ 2 (|0 + |1 ) as the initial state. The initial density matrix is then given by The target state has a φ-independent density matrix σ = 1 2 (|0 0| − i|0 1| + i|1 0| + |1 1|). We choose the trace distance between σ and the noise-averaged density matrix ρ(τ) as the error function of the gate: Notice that, due to noise averaging, ρ(τ) generically describes a mixed state. Even though we have a pure state with ρ(0) = ρ 0 at time t = 0 and the system evolves coherently for each realization of noise. A comment is in order regarding the choice of ρ 0 . The Hamiltonian is block-diagonal for one realization of noise so |0 and |1 are not coupled for a single realization. If the operation transforms ρ 0 to σ, we do not have leakage out of the ground-state manifold and the correct relative phase is imparted to these ground states. In the limit of small C[σ, ρ(τ)], we expect our protocols to correctly transform arbitrary superpositions of |0 and |1 , not just |ψ 0 defined above Eq. (8). It is important, however, that ρ 0 describes a superposition to capture the relative phase.

III. LINDBLAD MASTER EQUATION FOR MULTIPLICATIVE WHITE NOISE
Our tunable control parameters are ∆ j (t). However, the physical system may experience different parameters, such as ∆ S j (t) instead of ∆ j (t), which we try to impart to the system: We have used a multiplicative random fluctuation to account for the topological nature of the qubits. When ∆ j (t) is suppressed due to the separation distance between two Majorana modes, the wavefunction overlap between these Majorana modes is exponentially suppressed and the noise in the hybridization parameter should be similarly suppressed. Therefore, we do not expect any additive noise for a Majorana-based topological qubit. For simplicity, we focus on white noise with where W is the strength of noise. Cross-talk between different control parameters is also neglected, which is a reasonable assumption for a topological qubit. While 1/ f noise is more realistic than white noise, it has a slow power-law decay and considerable support over high frequencies. Previous studies suggest qualitative similarities between the two noise spectra [44]. White noise allows us to determine the noise-averaged density matrix by efficiently solving a single deterministic master equation and is therefore very convenient for performing optimization over the shape of the protocol, which requires solving the dynamics a large number of times. The master equation for the multiplicative white noise above has the following Lindblad form [44]: where H(t) is given in Eq. (1). The matrix representation of the master equation in our computational basis [see Eq. (2)] is given by with O i matrices defined in Eq. (3). Solving the above differential equation with the initial condition ρ(0) = ρ 0 [see Eq. (8)] gives the final noise-averaged density matrix at time τ.
Before proceeding, we calculate the above cost function [see Eq. (9)] using the linear protocol of Fig. 1(b) in the adiabatic limit of long times 10 < τ < 100. The results are presented in Fig. 2. In the absence of noise, the cost function approaches zero in the limit τ → ∞. Topological protection does not save us from the detrimental effects of noise, which, ironically, accumulate over the long times required for adiabaticity. As seen in Fig. 2, this gives rise to an antiadiabatic behavior, where longer times lead to greater errors. The antiadiabatic behavior provides an important motivation for this work. Finding protocols that minimize the cost function requires a parameterization of the protocol in terms of a finite number of variables. The cost function C can be uniquely determined by three bounded functions ∆ j (t) for j = 1, 2, 3 and 0 < t < τ. In other words, C is a functional of these ∆ j (t) functions. By discretizing time, as shown in Fig. 3, we can approximate an arbitrary protocol as a piecewise-constant function, treating the height of these pieces ( f n in Fig. 3) as the optimization parameters. Increasing the number of intervals allows us to scan the entire function space. Such direct minimization of the cost function is in the spirit of a machine-learning approach [53].
We then perform simulated-annealing Monte Carlo simulations where, in each step, a random time interval t n is chosen and the three optimization parameters corresponding to the interval are varied by a small amount. If all the new values are within the permissible bounds, the move is accepted according to Metropolis rules; It it accepted with probability one if it reduces C, (i.e., ∆C < 0) and with probability e −β∆C if it increases C (∆C > 0). The fictitious inverse temperature β is slowly increased as the simulations run.

B. Solving the master equation
The simulations rely on calculating C for arbitrary piecewise-constant protocols in each step. By arranging the elements of the 4 × 4 density matrix into a vectorρ of length 16, we can write the master equation (12) as where K is a 16×16 matrix. The formal solution for constant ∆ j , corresponding to each piece of the piecewise-constant protocol, is given byρ where we have divided the total time τ into N intervals of length ∆t = τ/N. The above expression leads tō ρ(τ) = e K N ∆t . . . e K 2 ∆t e K 1 ∆tρ 0 .
C. The optimal cost function from simulated annealing We start our simulations by fixing ∆t = 0.02. The number of variational parameters is then proportional to the total time τ. This direct optimization method does not make any assumptions about the shape of the protocols a priori. As it turns out, the optimal protocols for zero noise are bang-bang while the protocols for the noisy cases contain continuously changing segments. The results for the minimum cost function obtained from Monte Carlo simulations (with protocols individually optimized for each τ and W) are shown in Fig. 4. We find three distinct regimes: I: The minimum cost function for all W initially remains constant upon increasing the total time. While τ is less than a critical time τ c ∼ 1.3, we cannot change the minimum cost from its initial value using quantum evolution.
II: The second regime is characterized by the rapid decrease of the cost function as a roughly linear function of total time. For vanishing noise (W = 0), the cost function reaches zero at the end of the second regime at τ ∼ 2.6. For larger strengths of noise, C min is finite at the end of the second regime and increases with noise. III: Finally, a third regime occurs for longer times. For zero noise, the minimum cost function remains identically zero in this regime. For finite W, we observe a slow decay in C min as total time increases. The antiadiabatic limit is completely eliminated by these optimal protocols.
An important question concerns the fate of the minimum cost function as we increase the total time. In particular, will it saturate to a finite value or decay to zero? We extrapolate the data in the third regime to τ → ∞ using a polynomial fit to 1/τ. The extrapolation (using a cubic polynomial of 1/τ) suggests that C min likely decays to zero in the limit of infinite time, as shown in Fig. 5. Upon extrapolation, we find very small negative or positive values at 1/τ = 0, with 0 always contained within the error bars. Note that C is by construction nonnegative. Extrapolation with a fit to a quadratic function of 1/τ yields very small positive values (compared with the error bar). The linear fit is not as high quality, and a quartic fit produces very large error bars.
Each data point in Fig. 4 has a corresponding time-dependent optimal protocol ∆ j (t). We will present our results on pulse shapes after discussing the properties of the optimal protocols, using the Pontryagin's minimum principle.

D. Optimization with the bang-bang variational ansatz
The Monte Carlo simulations show excellent convergence except at the boundary between the first and the second regimes. We believe this is due to the presence of an infinite number of protocols in the first regime and proximity to this degenerate space in the beginning of the second regime. We can improve the results in this region using the bang-bang ansatz. We found that the protocols in the beginning of the second regime have a bang-bang character. This allows us to to use an efficient optimization with the bang-bang ansatz [54][55][56] to find the minimum cost function in this regime.
The bang-bang ansatz, motivated by Monte Carlo simulations, is shown in Fig. 6(a) and contains six variational parameters. It is easy to write the cost function C as a function of these t j parameters. We minimized C(t 1 , . . . t 6 ) using the interior-point minimization algorithm and found the minimum values of Fig. 6(b). We have better convergence and a lower minimum cost function (than Monte Carlo) at the interface of the first and the second regimes, and comparable performance in the first and the second regime. In the third regime, however, the bang-bang algorithm significantly underperforms the direct simulatedannealing, as the actual optimal protocols become very different from these bang-bang protocols.

V. PROPERTIES OF THE OPTIMAL PROTOCOLS FROM PONTRYAGIN'S MINIMUM PRINCIPLE
We begin by reviewing the Pontryagin's theorem for a generic dynamical system. For dynamical variables x(t) = (x 1 , x 2 , · · · , x M ) driven by control parameters ∆(t) = (∆ 1 , ∆ 2 , · · · , ∆ N ) such that we can define conjugate momenta p(t) = (p 1 , p 2 , · · · , p M ), evolving with where f = ( f 1 , f 2 , · · · , f M ). The optimal controls ∆ * that minimize a cost function C[x(τ)] at the final time t = τ then satisfy for all 0 < t < τ, where the optimal-control Hamiltonian is defined as The boundary conditions for the conjugate momenta are given by In the present problem, we treat the real and imaginary parts of the elements of the density matrix as our dynamical variables. As the density matrix is Hermitian and has trace unity, these dynamical variables are not independent. The independence of the dynamical variables is not a requirement for the theory, however. Let us define ρ αβ R ≡ Re(ρ αβ ), ρ αβ I ≡ Im(ρ αβ ).
We then represent the conjugate momentum to ρ αβ R (ρ αβ I ) by Π αβ R (Π αβ I ). As shown in Appendix A, the conjugate momenta satisfy a similar master equation to that of the density matrix, except the Lindblad double commutator appears with the opposite sign: We have combined the conjugate momenta for the real and imaginary parts of ρ αβ into one complex momentum Π αβ ≡ Π αβ R +iΠ αβ I . These terms can be collected in a matrix The optimal-control Hamiltonian can then be written as where F j and G j are known functions defined in Appendix. A.
In the absence of noise, i.e., W = 0, H is a linear function of ∆ j . Minimizing the optimal-control Hamiltonian (25) at a given time [see Eq. (19)] gives If F j identically vanishes over a finite interval, we have a singular segment and the optimal protocol can take intermediate values.
Otherwise, we have a bang-bang protocol, with the control parameter equal to either its minimum or its maximum allowed value. For W > 0, there is a third possibility for the minimum, which can occur at If 0 ∆ d j 1, we need to see which possible value of ∆ j , (0, 1, ∆ d j ) minimizes H, i.e., correspond to min(0, F j + W 2 G j , −F 2 j /4W 2 G j ). If ∆ d j is outside this range, we just need to compare the boundary values at ∆ j = 0, 1. For non-zero values of noise, the nonlinearity above can lead to protocols that are not bang-bang, even in the absence of singular segments.
The protocols ∆ j (t) uniquely determines both ρ and Π (and consequently F j and G j ) as a function of time. The initial condition for ρ is known and the deterministic master equation yields ρ(t) for all later times. In particular, for a piecewise-constant protocol, using Eq. (14), we can write For the conjugate momenta Π, we similarly know the equations of motion (23). However, we still need the boundary conditions to determine Π for all times. We can get these boundary conditions from Eq. (21), which, in the present problem, gives Explicit expressions for the above gradients are given in Eqs. (B5) and (B6) of Appendix B. As discussed in Appendix B, it is important to slightly generalize the cost function so that the above gradients are real. Again, solving the deterministic master equation (23) backward in time yields Π(t). For the piecewise-constant case, we have where the negative sign for W 2 in K(−W 2 ) follows from Eq. (23).

A. First regime
In the first regime, the simulations converged to many different protocols, all of which have the same cost function as the identically vanishing protocol independent of noise. To check this result, we created a large number random permissible protocols, without optimizing the variational parameters. We then plotted the resulting cost values in a two-dimensional histogram, shown in Fig. 7. We found that for all times in the first regime, the lower limit to the cost function coincided with the cost for τ = 0. While the reachable set grows with increasing τ, the growth is away from the target density matrix.

B. Second regime
Let us first focus on the noiseless case W = 0. In this case, we expect bang-bang protocols from the application of Pontryagin's minimum principle [see Eq. (26)]. Indeed, we find this behavior in out Monte Carlo results. In Fig. 8, we show two representative protocols in the second regime for W = 0. The optimal ∆ j is determined by sgn(F j ). For shorter times we see that ∆ 3 (t) = 1 while ∆ 2 and ∆ 1 have two switchings at intermediate times. For longer times, all ∆ j control parameters have two switchings. The times of these jumps are in agreement with Eq.(26). FIG. 7: Numerically generated approximate histogram of C for randomly generated protocols. Below a critical total time, no random protocol seems to decrease C, as indicated by the red lines at the initial C. Upon turning on the noise, we find similar protocols with segments where ∆ j is either 0 or 1. However, the sudden jumps get replaced by continuous segments, which become smoother upon increasing either τ or W. These continuous pieces are in full agreement with Pontryagin's minimum principle. They minimize H and are given by Eq. (27).
This remarkable agreement indicates that the Monte Carlo simulations have indeed converged to the optimal protocol despite their approximate discrete nature. The almost overlapping continuous pieces from Eq. (27) are closer to the true optimal protocol than the piecewise-constant Mote Carlo result. As a given protocol determines F j and G j , and, consequently, a new protocol containing 0, 1 or −F j /2W G j , this improvement can be iterated to arbitrary precision, leading to numerically exact optimal protocols. We note that although our Monte Carlo simulations are very accurate, the Pontryagin theorem is a powerful practical tool in finding the optimal protocol from a much less accurate simulation (or other minimization routine), which produces an approximate optimal protocol. If the approximate protocol is not too far from optimal, we have found that the above-mentioned iteration of Pontryagin's theorem can yield the optimal protocol.

C. Third regime
In the third regime, the protocols become more and more continuous. As shown in Fig. 10, the general pattern is as follows. For ∆ 1 , we have an S-shaped protocol (ascending, descending and ascending again). Upon increasing the total time, the pulse shape is squished in the vertical direction. When this S-shaped pulse is negative, the optimal protocol has a constant segment with ∆ 1 = 0. For ∆ 2 , there is a similar S-shaped protocol (now descending, ascending, and descending again). When this S-shaped pulse is negative (larger than 1), the optimal protocol has a constant segment with ∆ 2 = 0(1). ∆ 3 has a fixed pattern shown in Fig. 10.

VII. CONCLUSIONS
Motivated by an antiadiabatic behavior in the presence of high-frequency noise, which may be detrimental for all adiabatic quantum information processing schemes, we investigated optimal shortcuts to adiabatic evolution in the presence of white noise for a particular Majorana-based gate. We formulated the noise-averaged dynamics in terms of a master equation and performed Monte Carlo simulations to minimize a cost function of the final density matrix, which vanishes for perfectly adiabatic evolution in the absence of noise.
We found three distinct regimes as a function of the total time. As we are interested in small cost functions and precise gates, the third regime is of particular interest. We found that this optimization completely eliminates the antiadiabatic behavior.
We discovered very rich behavior in the pulse shapes of our optimal protocols. Starting with bang-bang protocols for vanishing noise, we found that the optimal protocols developed continuous segments as noise was increased. We performed an in-depth analysis of the application of the Pontryagin's minimum principle, which not only provides a mathematical understanding of the pulse shapes, but also allows us to arbitrarily improve approximate numerically obtained optimal protocols. This method may be of value to various quantum optimization schemes and variational quantum algorithms.
With the rapid development of quantum technologies, striking a balance between performance and robustness is crucially important. Topological quantum computing relies on adiabaticity and offers remarkable robustness. However, topological protection can be lost in the limit of long times due to the presence of noise or other factors like quasiparticle poisoning. Finding fast dynamical protocols that cancel the effects of noise suggest that even without strict topological protection, we may be able to perform high-performance and reliable information processing with Majorana-based platforms.
The equations of motion for the conjugate momenta can be written as which leads to Noting that all O j matrices are Hermitian, O 1 is imaginary and O 2,3 are real, we have