A geometric protocol for a robust Majorana magic gate

A universal quantum computer requires a full set of basic quantum gates. With Majorana bound states one can form all necessary quantum gates in a topologically protected way, bar one. In this manuscript we present a protocol that achieves the missing, so called, $\pi/8$ 'magic' phase gate. The protocol is based on the manipulation of geometric phases in a universal manner, and does not require fine tuning for distinct physical realizations. The protocol converges exponentially with the number of steps in the geometric path. Furthermore, the magic gate protocol relies on the most basic hardware previously suggested for topologically protected gates, and can be extended to any-phase-gate, where $\pi/8$ is substituted by any $\alpha$.


I. INTRODUCTION AND MAIN RESULTS
In two landmark papers, Kitaev suggested that nonabelian anyons could be used to store and process quantum information in a topologically protected way [1,2]. Furthermore, he outlined how one would try to realize the simplest of these non-abelian states, Majorana zero energy bound states (Majoranas for short), in a solid state system. Since, much activity has been dedicated to realizing Majoranas in quantum Hall states as well as quantum wells in proximity to superconductors, both theoretically [3][4][5][6][7][8][9] and experimentally [10][11][12][13][14][15][16][17], with significant recent success. Moreover, the experimental efforts recently shifted from a mere detection of Majorana signatures to concrete steps towards the realization of platforms that reveal their non-Abelian statistics and allow for quantum information processing via braiding [18,19].
Nevertheless, a stubborn roadblock still prevents us from proposing a topologically protected Majorana-based platform that is capable of universal quantum computation. Kitaev and Bravyi demonstrated that all gates could be realized in a platform that could carry out the Clifford gates (Hadamard, π/4, and controlled-not (CNOT) gates), and, crucially, possesses a magic state e −iπ/8 |0 + e iπ/8 |1 [20,21]. (Here |0 and |1 present the two quantum states of the qubit.) A four-Majorana network can realize a not-operation (σ x ) by braiding, and a Hadamard and π/4 gate through exchange. CNOT can also be implemented employing projective measurements [22]. Despite much inspirational effort [20,[23][24][25][26], there is still no protected or precise practical way to produce Majorana magic states, or the equivalent π/8 gate.
The quest for a magic gate is hampered by a pervasive challenge of quantum computing. Decoherence, and even more so, the lack of precise control of quantum information processing systems, necessitates the development of elaborate error-correction strategies, and quantum state distillation techniques. Topological quantum computing was developed as the ultimate fault-tolerant scheme, where environment noise is unable to decohere the quantum state of a qubit, since it is encoded nonlocally and spread over the entire platform. Also, gates that can be realized using topological manipulations such as braiding or exchange, are completely insensitive to the imprecision in the control of the system's parameters. When it comes to Majorana platforms, however, the absence of a topologically protected scheme for a magic gate requires us to revert to non-topological procedures [18,[27][28][29][30][31][32][33], and, therefore to rely on conventional error-correction schemes [20] which come at the cost of a significant overhead [34].
Indeed, procedures proposed so far for the realization of the Majorana magic gate require precise control of the coupling constants in the system. For instance, a relative phase between the two states of a two-Majorana parityqubit could be produced by bringing the two Majoranas close to each other; the tunneling between them produces a relative phase winding, as in, e.g., Refs. 27, 28, and 30. The integrated dynamical phase winding is dictated by the strength of tunneling, and the time of proximity, which need to be precisely controlled to achieve the coveted π/4 difference. A particularly clever way to produce a Majorana phase gate is through interference. A Majorana state bound to a quantum-mechanical vortex could be made to split between a path that carries out an exchange gate, and one that does not [26]. If the splitting is exactly equal, a π/8 gate will result. While experiments are progressing at a precipitous rate, such a level of control is unlikely to be reached soon. Furthermore, any improvements in the fault-tolerance of Majorana magic-gate realizations with respect to systematic machine errors, will dramatically reduce the amount of hardware necessary for state distillation.
In this work we present a robust scheme for obtaining a Majorana π/8 gate, which is insensitive to such machine control imprecision. The protection against such errors arises from universal geometric properties of the Majorana Hilbert space, alongside with the topological properties of the system. Our starting point is a 4-Majorana system arranged in a Y-junction. The Y-junction is the archetype model for a general Majorana exchange [35] and has been the simplest proposed platform for carrying out braiding. Moreover, as explained in Ref. [18,36], it could be realized with the most accessible Majorana supporting building blocks so far, which are spin-orbit coupled nanowires in proximity with superconductors. Our  FIG. 1. The Y-junction system. Lines and label indicate the model Hamiltonian of four coupled Majoranas γi, while the background shows a possible realization using wires (gray) proximity coupled to s-wave superconductors (green). We assume that couplings at the arms are determined through external imprecise controls. The true couplings are ∆ = δ+ f ( δ). It is indeed beneficial to think of the coupling to the x, y and z arms as geometric objects, namely, vectors in a three dimensional space. In addition to the couplings along the arms, any physical system will also exhibit couplings between the tips (dashed blue lines). These unavoidable couplings introduce a parity-dependent dynamical phase, and, along with the control uncertainty, are leading sources of error. multi-step scheme guarantees, under very broad assumptions, that the gate converges to the desired π/8 gate with an error δα such that ln δα = −O(N ), where N is the number of steps. The crucial assumption is that in the translation from ideal to actual control of the system, the spectral weight of the error function, i.e., uncertainty in the actual parameters of the system, is small at frequencies above the control clock rate.
The magic gate scheme we outline can be realized in any system, where the coupling between Majorana states can be controlled, even if imprecisely. Notwithstanding, we require the ability to decouple the Majorana states from each other with high precision, which is also required for topological protected Clifford gates, and therefore should be the case for any topological quantum information processing platform. The main setup, shown in Fig. 1, is described by three Majorana coupling parameters, ∆ x , ∆ y , and ∆ z . Exchange of Majoranas' positions is performed by changing the coupling between them in a specific time sequence using gates or fluxes in a system of finite size superconductors [36,37]. Our scheme, does not require any modification of this hardware; rather, it shows that certain time sequences of the coupling constants can result in exponentially high accuracy even for calculations that do not enjoy topological protection.
The main hindrance to a precise π/8 gate is the uncertainty in the values of ∆ a , a = x, y, z. It is this obstacle that our scheme completely eliminates. We assume the following: 1. The system controls are temporally constant for the duration of the gate.
Both assumptions can be relaxed. Achieving the same precision would, however, impose more stringent constraints on the rate with which the gate could be performed.
Our realization of the π/8 gate requires methods that are reminiscent of universal dynamical decoupling [38], and NMR echoes. (Though, here we deal with geometric phases rather then dynamical ones.) To eliminate the error due to the unknown device functions, f a ( δ), we will describe a trajectory in the δ space with 2N turning points. These will eliminate the first N − 1 coefficients in a Chebyshev-Fourier expansion of the errors. Since under rather broad assumptions (see App. B) Chebyshev expansion coefficients decay exponentially, the error due to the machine uncertainty can be made to vanish exponentially in the number of steps N . Note, however, that the underlying topological protection of Majoranas provides boundary conditions that are crucial to unlock this exponential behavior (see Sec. III and App. D).
A second grave problem arises from unavoidable dynamical phases due to couplings not included in the ideal, three-couplings, Y-junction setup. These dynamical phases, however, can be eliminated by repeating the gate protocol after applying a NOT gate to the Majoranas. Just as in a π spin-echo in NMR [39], this would cancel the error due to the extra coupling, as long as the system control functions are constant in time. (Additional steps to eliminate error when this is not the case can be applied.) In what follows, we first describe the ideal Y-junction system, and how to use it to get a non-protected π/8 gate (Sec. II). Next, we describe the Chebyshev universal geometric decoupling trajectory (Sec. III). We then present the echo method to eliminate the dynamical phase error (Sec. IV). Another possible source of error is the retardation effects in the system, which we discuss in Sec. V. Before concluding, we demonstrate our π/8 scheme through numerical simulations for particular generic error functions, and discuss the necessary scales of the coupling and motion rates (Sec. VI). The simulations are done by solving the full time dependent Schrödinger equation, without any assumptions on the adiabaticity of the process.

II. π/8 GATE IN AN IDEALIZED SYSTEM
The platform at the root of our scheme is the Yjunction system (see Fig. 1). It contains four Majoranas. Three of them, γ x , γ y and γ z , are located at the tips of the Y-junction and interact only with the fourth Majorana, γ 0 , which is at the center of the junction. The Hamiltonian for this system is: where we conveniently defined the Majorana vector γ = (γ x , γ y , γ z ) and the coupling vector ∆ = (∆ x , ∆ y , ∆ z ).
The topological nature of the Majoranas enters through the exponential dependency of the Y-junction couplings on physical parameters. For example, it depends exponentially on the distance between the Majoranas, and in the flux controlled qubit it depends exponentially on the flux applied to SQUIDs placed near the Y's arms [36]. Therefore, it is easy to essentially turn off one of the couplings. The Majorana state at the tip of this coupling is then an exact zero-mode of the Hamiltonian (i.e., it commutes with the Hamiltonian). This robustness (i.e. ∆ a = 0 if δ a = 0) lies at the heart of the protected Majorana braiding process and is also crucial for the π/8 gate discussed in this paper.

A. Exchange process
The presence of this zero energy mode allows for exponentially (topologically) protected exchange operations.
For simplicity let us assume that ∆ = ∆ 2 x + ∆ 2 y + ∆ 2 z ≡ ∆ remains constant through out the process. Consider the following trajectory: we start with ∆ z ≈ ∆ ∆ x , ∆ y , then move to ∆ y ≈ ∆ ∆ z , ∆ x in a continuous fashion while keeping ∆ x ∆. This is followed by similar moves to ∆ x ≈ ∆ ∆ y , ∆ z (while keeping ∆ z ∆) and finally returning the system to its original state ∆ z ∆ x , ∆ y (while keeping ∆ y ∆). This sequence is easily visualized as the arm of a clock indicating which coupling is strong; the sequence describes the arm making a full clockwise rotation. By doing so we carried out an exchange of the Majoranas γ x and γ y (see the upper panel of Fig. 2). [40] How do we mathematically see that the full turn of the 'clock arm' corresponds to performing an exchange? We can elegantly show this by taking advantage of the geometric analogy of ∆ to a vector in a three dimensional space described by spherical coordinates [41]. While ∆ is analogous to the radius-vector, we can additionally make use of the polar angle (θ) and the azimuthal angle (φ) of the spherical coordinates and their unit vectors. Noticing that ∆ = | ∆|(sin θ cos φ, sin θ sin φ, cos θ), and denotinĝ e θ andê φ as the unit vectors in the θ and φ directions, we now define: Clearly these are zero-modes: The exchange process consists of the unit vector ∆/| ∆| marking an octant on the unit sphere. The octant is bounded between the φ = 0, θ = π/2 and φ = π/2 planes. See the lower panel of Fig. 2.
These two zero-modes combine into a single Fermi annihilation operator: FIG. 2. A visualization of the exchange process as the turning arm of a clock. First ∆z ∆x, ∆y, to indicate that the line presenting the coupling between γ0 and γz is bold. Then ∆y ∆x, ∆z, then ∆x ∆y, ∆z and finally ∆z ∆x, ∆y again, so that the arm of the clock completes a full turn. This process can also be visualized as a line covering an octant on a unit sphere. The Berry phase difference of the two parity sectors accumulated in this process is equal to the covered solid angle, −π/2. We show in the text (see App. A) that this gives rise to a −π/4 phase gate, meaning a phase ∓π/4 for each fusion channel. (The − sign appears due to the clockwise orientation of the trajectory, and the convention we chose.) Since we can make one of the coupling constants exponentially smaller than the other two the trajectory in the parameter space is glued to the edges of the octant making the accumulated Berry phase difference equal to −π/2 with exponential accuracy. and we are interested in the difference of the accumulated Berry phase, 2α, between the process where the system is in its ground state, |0 , defined by a |0 = 0, and its partner, which is a † |0 = |1 . As we show in App. A, the phase difference accumulated during the process of exchange, or any process (described by a contour c) for that matter is Just like the Berry phase of a spin system, the phase gate angle α is equal to half of the solid angle Ω c enclosed by a contour c. Indeed, for the exchange process of Fig. 2 we obtain: (The − sign appears because the contour is counter clockwise). This corresponds to the operator U exchange = e − π 4 γ φ γ θ -a π/4 gate.

B. A naive π/8 gate
When considering how to perform a π/8 gate, the calculation of the exchange gate is very suggestive. The exchange entails a π/4 gate; all we need is half the angle. For half the angle, we simply need to cover half the area of the octant.
While the geometric construction seems to be taking advantage of simple area consideration, the result contains deeper roots. Conceptually, we would obtain a π/4 phase difference between the |0 and |1 states. Instead of carrying out an exchange between the two Majoranas, we would have managed to do the following feat: split one of the Majoranas into an equal superposition, where one part carries out the exchange, and the other doesn't. Next, we reunite the two parts. The interference between the two processes will yield to the relative phase 1 √ 2 (1 + i) = e iπ/4 . This process, which combines the weirdness of quantum mechanics with that of Majoranas is precisely what the Y-junction sequence presented above is performing.
Unlike the exchange process, however, there is no protection for the φ = π/4 plane. Equivalently, we can keep δ x = δ y within our control module, but clearly ∆ x − ∆ y = f x − f y = 0. The error in the device control may introduce an arbitrary error in our computation. An additional complication arises due to the need to go through the center region of the octant. In the π/8 trajectory, we can not avoid a region where all three Majorana couplings have similar strengths. Invariably, they give rise to a direct next nearest neighbour coupling between the Majoranas at the Y-junction tip [42]. In this case the ground state degeneracy is split, and the relative phase between the |0 and |1 states receives a time-dependent dynamical phase on top of the pathdependent Berry phase contribution. In the following sections we demonstrate how these errors could be universally corrected.

III. SYSTEMATIC ERROR ELIMINATION USING UNIVERSAL GEOMETRIC DECOUPLING
In this section we will analyze an universal scheme which allows to dramatically reduce errors of the 'magic' π/8 phase gate as compared to the naive implementation.
Intuitively, it seems that smooth errors due to the imprecise coupling constants tend to be canceled in contours that have the snake like shape as in Fig. 4. The unwanted geometric phase accumulated on the way from θ = 0 to θ = π/2 is subtracted by a similar perturbation on the way back from θ = π/2 to θ = 0. We will treat the snakelike trajectory of Fig. 4 as a variation trajectory and optimize the turning points φ N 1 , φ N 2 , . . . , φ N n , n = 1, . . . , 2N in order to minimize the error in the accumulated phase.
We have to choose the turning points φ N n such that α = π/8. It is useful to perform the transformation x = 2 π φ and y = 1 − cos θ then the topologically protected boundaries are transformed to the boundary of a square x ∈ [0, 1] and y ∈ [0, 1]. We find The last equation is correct for any closed contour c. The contour c in the x − y plane is the image of the contour c in the φ − θ surface. For example, the topologically protected contour of Fig. 2 simply follows the boundaries of a unit square in the x − y plane. We now want to find a contour c in the x − y square that gives a c = 1/2 with an exponential accuracy. For the perfect snake contour of Fig. 4, this leads to a single condition for the 2N turning points x n = 2 π φ n . The idea of our geometric decoupling scheme is to use the remaining 2N − 1 degrees of freedom to systematically reduce the effect of errors. The latter will change the contour from c to C with a parametric representation (X(t), Y (t)) different from the ideal desired contour (x(t), y(t)). When assuming that the ideal contour doesn't stop in regions with finite errors there is a one to one correspondence between t and (x, y) which allows to parameterize the error functions as δx(x, y) = X(x, y)−x and δy(x, y) = Y (x, y)−y. Importantly, due to the topological protection, the functions δx and δy must vanish on the square boundaries (δx on x = 0, 1 and δy on y = 0, 1).
Using the undisturbed coordinates (x, y), the area A C = − C XdY encircled by the disturbed contour can be written as a sum over the 2N vertical sweeps n where y changes from 0 to 1 (for odd n and from 1 to 0 for even n) while x is fixed at x n : The alternating sum over x n yields the desired enclosed area of the perfect contour a c and is reproduced in Eq. (9) due to the vanishing boundary conditions of δy . Equation (9) then further suggests that the remaining effect of the errors δx and δy can be mapped to that of a snake contour with straight vertical trajectories (∂ y δx = 0) where the turning points are shifted by δx eff (x n ) relative to the perfect implementation. Note that the vanishing boundary conditions of δx at x = 0, 1 carry over to δx eff . We now want to systematically cancel the alternating sum over δx eff (x n ). To this end we use the expansion where at the square boundaries the set of orthonormal functions P m (x), m = 1, . . . , ∞ vanishes, explicitly P m (0) = P m (1) = 0. Inserting this expansion into Eq. (10) yields By assumption, the error function f and its mapping to δx eff all have physical origins, we can therefore assume that they are smooth and analytic. In addition, they are bounded due to the topologically protected gluing to the square boundaries so we may conclude that lim m→∞ A m = 0, for any orthonormal set of basis functions P m (x), m = 1, . . . ∞. Choosing x N n = 2 π φ N n properly we can eliminate the first M = 2N −1 components of the expansion which protects the phase gate by reducing the error to A protected π/8-phase gate can therefore be implemented when aiming for turning points that fulfill the (13) with a c = 1/2. These are 2N non-linear equations for 2N unknowns x N n , n = 1, 2, . . . 2N . Remarkably, if solutions of Eqs. (13) exist, the errors can be canceled up to order 2N − 1 independent of the expansion coefficients A m . The scheme is therefore universal and independent of the details of the errors as long as they are smooth. In fact, Eqs. (13) resemble similarities to the concept of universal dynamical decoupling [38].
For a good choice of expanding functions P m we rely on the common knowledge in numerical analysis that expansions of analytic and bounded functions in terms of Chebyshev polynomials [45] converge very quickly with expansion coefficients decaying exponentially A m ∼ e −m . (See also App. B, Legendre or Laguerre polynomials are as good as the Chebyshev ones, however for the Fourier expansion discussed in App. C we expect that the convergence is slower.) Expanding in terms of Chebyshev polynomials therefore leads to an exponential suppression of the gate errors in the number of turns δa c ∼ e −2N .
Interestingly, the topological protection is crucial for enabling such an exponential decay of errors in the number of turns. Solutions of Eqs. (13) for non-vanishing a c only exist if x is linearly independent from the set of functions {P m (x)}. Fortunately, the function x violates the topological boundary conditions P m (0) = P m (1) = 0 and is therefore orthogonal to the basis functions P m . If the topological protection is relaxed and the errors are expanded in more general P m , an exponential decay of the expansion coefficients only allows for solutions when a c is exponentially small in N (see Appendix D for details).
As a particular choice of P m that uses the power of Chebyshev polynomials and is compatible with the topological boundary conditions we use with T m (x) being Chebyshev's polynomials of the first kind. Note that the order m = 0 that involves only linear terms in x vanishes identically which again reflects the orthogonality of the function x with respect to {P m }.
The solutions of Eqs. (13) for a c = 1/2 (π/8 gate) can be expressed analytically and are given by The general solution for other a c 's can be found numerically and will be discussed in Sec. VI below.

IV. THE DYNAMICAL-PHASE ERROR AND ITS ELIMINATION WITH A PARITY ECHO
The Chebyshev protocol above efficiently eliminates the systematic machine error, but at the same time introduces an equally potent source of error: an uncontrollable dynamical phase. The geometric decoupling method of Sec. III requires that at certain times all the couplings ∆ i are comparable. In these time spans, it is unavoidable that substantial direct couplings emerge between the tip-Majoranas, γ x , γ y , and γ z . These couplings induce a finite energy splitting between the two, otherwise degenerate, parity sectors of the system. This splitting integrated over the gate's duration, will distort the relative phase we seek to control. This distortion can be completely eliminated by carrying out a parity-echo: canceling the dynamical phase accrued in the gate by the opposite dynamical phase accrued from the same gate when reversing the parity for the low-lying parity sector. In order to add instead of subtract the wanted geometric phase the second gate is applied with reversed trajectory.
Let us first consider the strength and origin of the parasitic couplings. Majoranas i and j will in general be coupled by a term of order ∆ i ∆ j /∆, where∆ is the energy scale of high energy modes that were integrated out to obtain the four-Majorana low-energy Hamiltonian, typically of the order of the superconducting gap. In the spherical polar coordinates, the parasitic couplings induce the term: where γ θ and γ φ are defined in Eq. This strategy for mitigating the dynamical phase error, however, results in strong constraints on the speed with which the π/8 gate could be carried out.
A superior strategy employs spin-echo-like schemes [39]. If we switch the sign of the Hamiltonian (here δH) for half the duration of the protocol, the dynamical phases from the two halves of the procedure exactly cancel, regardless of how strong they are. Indeed, the Majorana structure of the energy splitting, Eq. (16), allows to switch the sign of δH without fine tuning of any parameters by applying a parity flip (NOT gate of the qubit): γ θ γ φ → −γ θ γ φ . One possible implementation of this "parity echo" includes (1) carrying out a geometrically robust π/16 gate (by solving Eqs. (13) for a c = 1/4 instead of a c = 1/2), (2) performing a parity flip, and (3) carrying out the same π/16 gate with a reversed direction of the contour. The contribution to the geometric phase switches sign twice (parity flip and reverse of direction) thus adding up to an overall phase of π/8. The dynamical phase, however, is direction independent, and thus cancels out after the echo is completed. Alternative parity echoes, and echoes based on a sign change of F (θ, φ) by manipulating θ and φ ("angular echo") are discussed in App. F. Needless to say, these simple echo procedures rely on the system not changing during the gate execution. More complicated echoes (similar to the original idea of dynamical decoupling) could also allow us to suppress finite frequency noise effects in the dynamical phase.

V. RETARDATION EFFECTS ELIMINATION USING ECHO
In addition to the geometric variations of X(x, y) and Y (x, y), inductive and capacitive effects might introduce velocity dependent changes. For the vertical snake protocol (Fig. 4) this could, in particular, introduce terms of the form X(x, y,ẏ) = n X (n) (x)ẏ n . All the terms with an even power are cancelled by our protocol as they simply lead to an overall (lateral) shift of X for the trajectories of constant longitude. The odd terms are problematic as they tend to influence the northwards (at x N 2n ) and southwards (at x N 2n+1 ) trajectories in an opposite and correlated manner so that the total area enclosed by the trajectory will deterministically change.
Fortunately, these terms with an odd power are eliminated by the same parity echo protocols of Sec. IV that eliminate possible dynamical phases. The idea of the parity echo is to cover half of the solid angle for the geometric phase on the forward run and the other half on the backward run. If odd order terms lead to a velocity dependent deformation of the contour that changes the Berry phase of the first run, there will be an opposite deformation on the way back where velocities are reversed. In the total geometric phase of the echo, the velocity dependent effects, therefore, cancel out. There might be higher order errors from an imperfect cancellation of the dynamical phase since the forward and backward contour are now slightly different. Note, however, that when implementing the same velocity profileẏ(t) (up to the switching sign) for each trajectory of constant longitude, the error for the imperfect cancellation due to the odd order velocity dependent terms will be again an alternating sum of the 2N trajectories which cancels in the limit of large N [46].

A. Verification of the Chebyshev protocol
First, let us demonstrate the robustness of the above protocols in the absence of a dynamical phase. For this purpose, we conducted a simulation of the full time evolution of the system where the change of the Hamiltonian H(t) [cf. Eq. (1)] in parameter space is described by the vertical snake contour of Fig. 4. We extract the gate angle α from the time evolution operator obtained by numerically solving the time dependent Schrödinger equation corresponding to H(t).
Adiabaticity is certainly a concern in our protocol [47][48][49]. In order to reach the adiabatic regime more easily we slow down the speed of parameter change close to the sharp turning points of the protocol. This allows us to stay well inside the adiabatic regime for moderate time spans τ between the turns. Throughout this paper we use τ = 25/∆ which yields non-adiabatic phase errors < 10 −10 .
As expected, we find that a perfect implementation of the protocol with turning points φ N n = (π/2)x N n , as in Eq. (15) gives a π/8 gate up to (the simulating) machine accuracy. Systematic errors in the control of the Majorana system would give rise to deviations of the implemented Φ N n (φ N n ) from the perfect turning points of Eq. (15). Fig. 5(a) shows an error model that leads to substantial deviations, corresponding to an analytic control function Φ(φ) [50]. Despite this control error, and the strong deviations of the turning points (of up to 20%), the relative phase error of the gate δα = (α − π/8)/(π/8), vanishes exponentially fast with increasing N [see Fig. 5(c)]. More complicated error functions Φ(φ, θ) that include cross correlations between θ and Φ yield (as expected by the general argument in Sec. III) a similar exponential decay and are discussed in App. G. Notice that we chose here an error function that is analytic on the real axis with complex poles. As explained in App. B this leads to a simple exponential de-cay rather than an exp[−N log(N )] decay when the poles are absent.
As noted above, our procedure is not limited to phasegates with α = π/8. Although simple analytic solutions of φ n will in general not exist for α = π/8, Eqs. (13) can be solved numerically for any α = a c π/4 . We checked that these numerical solutions indeed possess the same stability and protection as the π/8 case, reproducing essentially the same behavior as in Fig. 5.

B. Simulations of the parity-echo procedure
Let us next include in our simulations the errors due to the dynamical-phase discussed in Sec. IV. Finite nextneighbour Majorana couplings give rise to dynamicallyinduced phase errors that have to be addressed independently of the control-function errors. To study their effect we added a term to the Hamiltonian. In the absence of steps to mitigate the dynamical phase error, as expected, the performance of the Chebyshev protocol is limited, and we find that δα cannot be reduced beyond 2N rε∆τ [see Fig. 5 where r ≈ 0.06 is a numerical constant given by averaging the energy splitting over time (it is much smaller than 1 because our protocols slow down close to the turning points where the energy splitting is small). The echo sequences discussed in Sec. IV, however, can cancel the dynamical phase effects. In Fig. 5(e), we carry out the parity-echo protocol: first, applying a π/16 gate, then performing a parity switch (NOT gate), and finally applying the π/16 gate in reverse. This, as the numerical results of Fig. 5(e) show, fully restores the ε = 0 behavior.
Interestingly, finite next neighbour couplings as in Eq. (17) also alter the geometric aspects of the problem, since they also modify the Berry phase acquired by the system. This is due to the change of the fermionic low energy mode a (zero-mode at ε = 0) which leads to extra contributions to the Berry phase and thus geometric errors. Since these errors, however, vanish at the edges of the octant in parameter space where the zero-mode is recovered, they are automatically taken into account and corrected by our snake contours [51]. When applying the parity echo, we therefore observe an expected e −N log(N ) decay of δα even for 'perfect' implementations Φ = φ.

C. Consequences of time dependence in the system
Finally, we note that there are some remaining errors from the temporal change of the system between the parity-echo protocol steps. We model this numerically by changing ε to ε + δε on the second part of the echo. We find that the performance of the gate is then limited by 2N rδε∆τ [see Fig. 5(f)]. Note that it is in principle possible to reduce ε and the corresponding δε by choosing a coupling strength ∆ much smaller than the superconducting gap. Also note that δε itself can be very small. A likely source of a time-dependent change of ε is charge noise, which leads to typical δε ∼ 10 −3 ε [52] and even can be reduced to δε ∼ 10 −6 ε when tuning the system to the charge insensitive sweet spot.

VII. SUMMARY AND DISCUSSION
In this manuscript we have suggested a geometric multi-step protocol realizing a π/8 phase gate with an exponential accuracy in the number of steps N , and generalized it to any α phase gate, which may be useful as practical shortcut protocols for producing desired phase gates. We have demonstrated the protocol on a setup [36], where 4 Majoranas are situated at the 3 tips and the center of a Y-shaped junction. This makes the scheme particularly appealing, as one may use any hardware realizing Majoranas that does not support universal quantum computation and make it universal with the multi-step protocol we suggest.
Manipulating the system through a sequence of coupling constants between the 4 Majoranas induces various phase gates. We mapped the coupling sequences to contours on the Bloch sphere, and, in particular, we showed that the topologically-protected exchange process (the π/4 phase gate) corresponds to a contour encircling an octant of the Bloch sphere. A contour covering a solid angle Ω produces an α = Ω/2 phase gate. Due to the topological nature of the exchange process deviations from that contour are exponentially small in the laboratory physical parameters.
Contours that cover parts of the octant can be interpreted as a split of one Majorana where another is exchanged with only one of the split parts, leading (in case when the two parts have an equal weights) to a π/8 phase gate. The exact splitting portions depend on details and are not protected by any symmetry or topology. The algorithm we suggest in the form of a contour with N switchbacks does not enjoy the full topological protection such as the π/4 phase gate. However, because parts of the trajectory are at the boundary of the octant we were able to show that the geometric phase accumulated in the process is π/8 with an exponential accuracy. We demonstrated that realizations of different phase gates are also possible. Although the topological protection at the boundary is crucial in our scheme, it might still be interesting to study whether extensions exist that also improve the performance of geometric quantum computation schemes in non-topological systems [53][54][55].
Because parts of the contour are in the vicinity of the octant's middle, next nearest neighbour coupling effects may not be negligible. We have analyzed parity and angular echo protocols eliminating these dynamical effects, as well as parasitic retardation effects.
Assuming that the induced superconductor gap energy is ∆ SC = 3K ≈ (60GHz), the couplings between the Majoranas is at max ∆ = (5mK) ≈ 100MHz (leading to next neighbour couplings ε∆ ≈ 0.2MHz), a π/8 gate (with 2N = 10 turns) can be operational at turn frequencies 1/τ ≈ 10MHz. The latter fulfills the adiabatic condition τ ∆ 1 and leads to a small dynamical phase error which we estimate to be δα = 0.01. Applying a parity echo with temporal fluctuations δε < 0.01ε would then allow to reach relative errors < 10 −4 (Note that this is a conservative estimate. When the temporal fluctuations are only caused by charge noise the overall errors would be 10 −5 and even 10 −8 at the charge insensitive sweet spot, see Sec. VI C).
The scheme we suggest relies on a variational protocol and can still be improved upon. In particular, we did not try to optimize the contours used in the calculation, and the snake-like contours may not be the most optimal. Also, we employed only a one step protocol for the cancellation of the dynamical effects. We are currently considering possible improvements along these lines.
Finally, we would like to emphasize that the scheme proposed here is universal. It should be operational for all realizations of Majoranas and all models of environmental noise.
[42] Notice that if the couplings between γ0 and γa is exponentially small then the next nearest neighbour couplings are exponentially smaller then them. So in principle these dynamical corrections are exponentially small. It however requires exponential slow rate of exchange which is undesirable.
[43] Notice that the term c X(x, y)∂xY (x, y)dx vanishes since the horizontal parts of the contour (cf. Fig. 4)  To define the α phase gate we assume that there are two Majoranas γ θ and γ φ forming an annihilation fermion operator a = (γ θ + iγ φ )/2.

α -phase gate and geometrical phases
To calculate the phase one accumulates in physical models with a particular trajectory of the model's parameters, we note that iγ θ γ φ = 2a † a − 1, so that to find α associated with a trajectory we need to calculate the difference in the geometrical phases for empty a-state 0 ∆(t) , defined by c 0 ∆(t) = 0 and an occupied state 1 ∆(t) = c † 0 ∆(t) .
for these instantaneous states the geometric phase accumulated in a contour c is given by, The last term in the equation vanishes since 0| a † = (a |0 ) † = 0. To find the geometric phase in terms of the trajectory in the θ, φ plane we use again the definition a = (γ θ + iγ φ )/2 so that we have using now the mathematical identities ∂ θ γ θ = −γ r , ∂ θ γ φ = 0, and ∂ φ γ θ = cos θγ φ , ∂ φ γ φ = − cos φγ x − sin φγ y = − cos θγ θ − sin θγ r and using the expression for the Berry phase in terms of the operator in Eq. (A1) we find that the variation with respect to the polar angle θ vanishes a, ∂ θ a † = −1/4 {(γ θ + iγ φ ) , γ r } = 0 and the variation with respect to the azimuthal angle φ gives: We therefore have: With Ω c being the solid angle enclosed by the contour c.

Evolution of operators
The operators γ θ and γ φ defined in Eq. (2) evolve explicitly due to the instantaneous change of the coupling constant and due to the Berry phase. To find their evolution it is convenient to use the evolution operator [49, 56], U(t) = n=0,1 e iβn(t) |n(t) n(0)| , with |n(t) the instantaneous eigenstates and β n (t) = E n t + t 0 n(t)| ∂ t |n(t) dt the acquired phase. The dynamical phase E n t may be omitted here since the eigen-energy of both the occupied and the empty states are zero, E n = 0, n = 0, 1. Evolving the operator a † = (γ θ − iγ φ )/2 we find with Ω c (t) = i t 0 ( 1(t)| ∂ t |1(t) − 0(t)| ∂ t |0(t) ) dt = − t 0 cos θ dφ dt dt and the bar aboveā indicates that it is a physical operator with the additional evolution due to Berry's phase. Using the relation of Eq. (4), we establish the following relation between the instantaneous Majoranas,γ, and the physical Majoranas γ that evolve with the Berry phase: We find that the relation between the instantaneous and the physical operators is given by an additional rotation with an angle equal to the Berry phase. Notice that for a closed contour with Berry phase Ω c = π/2 the instantaneous γ's return to their original positions while the physicalγ's exchange their positions and γ θ acquires an additional sign. For Ω c = π/4 the physical operators are a superposition of the instantaneous ones.

Appendix B: Properties of Chebyshev polynomials
The Chebyshev polynomials of the first kind are defined as: A direct substitute of this definition in Eq. (14) proves that: To show that x N n in Eq. (15) solve the set of 2N non linear equations [Eq. (13)] we use the fact that T m satisfy the discrete orthogonality condition: where x k = cos π 2k+1 2N are the N Chebyshev nodes. For completeness we quote here the following theorems on the Chebyshev expansion, their proof can be found in many text books, for example, you may consult Ref.
Consider the following expansion of the function f (x) with the partial sum of the form f (x) ≈ S n (x) =

Appendix C: Other protocols
Here we discuss other possible protocols that reduce the systematic control errors. A straight forward variation of the vertical snake contour of Fig. 4 is to perform horizontal sweeps with turning points θ N n (see Fig. 6). Since XdY = − Y dX on any closed path, exchanging the role of X and Y has no effect on the arguments in Sec. III. Then we will have a set of y N n instead of x N n . One should, however, note that due to the nontrivial relation y = cos(θ), Y (y) = cos(Θ(arccos(y))) is not necessarily analytic when Θ(θ) is analytic. It might then be useful to instead implement modified turning points y N n → ξ(y N n ) (with some function ξ(0) = 0 and ξ(1) = 1) such that Y (ξ(y)) is analytic.
It is also possible to choose different basis functions for the expansion of the errors (while keeping in the vertical snake contour). In particular, a choice that agrees with the boundary conditions P m (0) = 1 and P m (1) = 1 is P m (x) = sin(mπx). Using the symmetry properties of the sine function we conclude that with we satisfy Eqs. (13) with a c = 1/2. Interestingly, in this case, we can also find analytic solutions x N n = (n − 1/2 + (−1) n (a c − 1/2)/(2N ) for arbitrary a c . Unfortunately, the above Fourier expansion does not enjoy the same exponential convergence as the Chebyshev polynomials and errors in the gate would decay polynomially in the number of turns N .

Appendix D: Topological protection and exponential convergence
In this section we detail the connection between the exponential convergence of our scheme with the turn number N and the underlying topological protection of the Majoranas. The topological boundary conditions enter at two crucial points in the derivation of Sec. III. First, an over or undershooting δy(x n , 1) = 0 would add an extra error to Eq. (9), that goes beyond an effective shift of turning points. In the absence of the topological boundary conditions the expansion δy(x n , 1) = ∞ m=0 B m T m (2x n − 1) will in general have a constant contribution B 0 = 0. Canceling this error then would require 2N n=1 (−1) n x n = 0, which only allows trivial a c = 0.
Even when assuming that δy fulfills the topological boundary conditions, we arrive at similar limitations, when considering the effect of unrestricted δx. Expanding the latter requires a general Chebyshev expansion δx = ∞ m=0 A m T m (2x−1) (as opposed to the constrained version in the main text). The problem is that we can also expand x = ∞ m=0 C m T m (2x−1). We then find that the first Eq. in (13) takes the form where the orders m = 0 . . . 2N − 2, drop out because of the error canceling Eqs. in (13). With T m being Chebyshev polynomials of first kind, the only non-vanishing expansion coefficient is C m=1 = 1 such that we again find a c = 0. In principle we could use different basis functions for the expansion but if their coefficients decay (as required) exponentially with the order of the expansion, Eq. (D2) constraints a c to be at most exponentially small in 2N − 1.
where the functions c θ/φ vanish when θ/φ are 0 or π/2 due to the topological protection at the edge of the octant in parameter space. In particular we choose two polynomials for c θ/φ obtained from interpolating 5 random numbers in between the vanishing boundary conditions. Figure 8(a) shows the resulting contour for variations |c θ c φ | 0.01. As expected, we still find an expo-nential decay of the phase error with increasing number of turns (for the same reason as in the main text the strength of the errors is not crucial for the exponential behavior). Fig. 8(b) shows the corresponding behavior for ε = 0 with an decay that is only slightly weaker as in the absence of the cross terms. Finite values of ε can be corrected by echo protocols similar to the main text.