Holonomic implementation of CNOT gate on topological Majorana qubits

The CNOT gate is a two-qubit gate which is essential for universal quantum computation. A well-established approach to implement it within Majorana-based qubits relies on subsequent measurement of (joint) Majorana parities. We propose an alternative scheme which operates a protected CNOT gate via the holonomic control of a handful of system parameters, without requiring any measurement. We show how the adiabatic tuning of pair-wise couplings between Majoranas can robustly lead to the full entanglement of two qubits, insensitive with respect to small variations in the control of the parameters.


Introduction
Topological quantum computation (TQC) is an approach to quantum computing that aims at minimizing decoherence at the hardware level, by exploiting topological properties of non-local degrees of freedom composed of non-Abelian anyons [1][2][3]. The latter are exotic quasiparticle excitations, which feature non-trivial exchange statistics, described by multidimensional representations of the braid group. A collection of non-Abelian anyons is embedded in a degenerate ground state manifold, which allows to non-locally store quantum information and to process it by implementing unitary transformations via braiding.
The building block of Majorana-based TQC is the Majorana qubit, consisting of four MZMs. By physically braiding those MZMs it is possible to implement all single-qubit Clifford gates [21][22][23]. Those gates are topologically protected, as their outcome exclusively depends on the topology of the trajectories adiabatically followed by the anyons in a 2 + 1 dimensional space. Importantly, the braiding of a single pair of MZMs can be realized in several ways, which are all equivalent to a physical exchange of the two non-Abelian anyons [24][25][26][27][28][29][30]. Indeed, by considering the presence of additional (hybridized) ancilla Majoranas, we can perform braiding by properly tuning pair-wise couplings between different MZMs [31,32], or by performing sequential projective parity measurements [8,[33][34][35][36][37][38]. Non-Clifford operations such as the T gate cannot be realized via Majorana braiding and necessarily rely on implementations that are not topologically protected and eventually require error correction schemes such as magic state distillation [23,39].
To achieve universal quantum computation, single-qubit gates must be supplemented with an entangling gate, such as the CNOT gate. Unfortunately, this two-qubit Clifford gate cannot be realized within a scalable architecture by exclusively using Majorana braiding operations [22,40]. The measurement-based approach allows us to overcome this issue by implementing the CNOT gate by performing high-fidelity projective measurements of the (joint) Majorana parities [8,35,[41][42][43][44]. However, while measurement-based TQC has proven to be extremely valuable for the future development of a fully scalable topological quantum computer, the required measurement protocols still represent a formidable challenge [35,45,46]. For the time being, it is therefore desirable to devise and characterize alternative schemes, which do not rely on high-fidelity measurements but still allow to robustly entangle distinct topological qubits.
In this work, we propose a measurement-free realization of the CNOT gate based on a holonomic approach. The key idea of holonomic quantum computation is to exploit non-Abelian geometrical phases to implement unitary operations on a degenerate eigenspace of the underlying Hamiltonian [47]. Those gauge-invariant phases emerge when the parameters of the system are tuned along degeneracy-preserving closed loops in parameter space. This approach is rather general and has been successfully exploited in non-topological quantum computation schemes [47][48][49]. Therefore, it is interesting to utilize holonomic techniques also in TQC. Indeed, the braiding of Majoranas itself can be interpreted as a holonomic process, where the system follows specific, topologically-protected loops in the three-dimensional parameter space of pair-wise Majorana couplings [8,31]. The advantage of the holonomic description of the braiding is that it can be easily generalized, both by considering loops with a different structure, within the same parameter space, and/or by considering a different parameter space altogether. In the first case, a careful modification of the loops can effectively implement non-Clifford (and non-topological) gates, such as the T gate [50]. We consider the second case and demonstrate that it is possible to implement entangling gates, such as the CNOT gate, by working in specific parameter spaces of a two-qubit system. In presence of physical constraints on the fermion parity of individual qubits, provided, for instance, by a finite charging energy [30,32,36], our holonomic entangling scheme is robust with respect to the presence of otherwise detrimental couplings between MZMs and/or limited control of system parameters.
The article is organized as follows. In Sec. 2, we introduce the structure of the specific Majorana qubits under consideration. We also briefly review the concepts of holonomic quantum computation. In Sec. 3, we propose two different holonomic implementations of two-qubit entangling gates, whose robustness is analyzed in Sec. 4. Finally, we summarize and discuss our findings in Sec. 5.

Model
For the sake of generality, we do not focus on a specific physical implementation of Majorana qubits. Instead, we discuss a generic low-energy effective model consisting of several couples of Majorana modes (or Majoranas). The latter are described by the self-adjoint operators γ j = γ † j which obey If they commute with the system Hamiltonian [H, γ j ] = 0, we refer to them as MZMs. Two Majorana operators define a (possibly non-local) fermion whose associated parity operator reads Two degenerate states with opposite parity P jk |0 jk = |0 jk and P jk |1 jk = −|1 jk can be therefore associated with every pair of MZMs. As the global fermion parity of an isolated system is fixed, a working Majorana qubit requires the presence of four different MZMs γ i (with i = 1, 2, 3, 4) [22]. Without loss of generality, we consider the global parity of the Majorana qubit to be even and choose the following basis for the computational space of the qubit |0 = |0 12 0 34 , Our goal is to operate on the Majorana qubit by following an holonomic approach, i.e. by changing in time some of its parameters. However, a system consisting only of four MZMs does not allow for such a manipulation, as it does not feature tunable parameters. To overcome this issue, we add two additional Majoranas, γ 0 and γ 5 , to the system and consider the general Hamiltonian which describes pair-wise couplings within a "five-point star" scheme, as sketched in Fig. 1(a). When no operations are performed on the qubit, it is in the idle configuration which features a single non-vanishing coupling strength c 5 = Θ > 0 and leaves us with the four MZMs γ i (i = 1, 2, 3, 4). Without loss of generality, we consider the total fermion parity of the system to be even P = P 12 P 34 P 05 = 1. In this case, the computational space is spanned by the two degenerate ground states The excited states, with the same total fermion parity, are separated from the ground states by an energy gap ∆E = 2Θ. The advantage of this five-point star architecture is that it is tunable. Indeed, it is possible to tune the coupling strengths c j away from the idle configuration without destroying the qubit, keeping the fixed-parity computational space degenerate and separated from the excited states by a finite energy gap. A sufficient condition for the stability of the qubit is that, at each time t, either one or two of the five different couplings strengths c j must be nonvanishing [31,50,51].

Holonomic description of Majorana braiding
We now briefly review a protocol which allows us to braid a couple of MZMs, say γ 1 and γ 2 , thus implementing a topological quantum gate. This protocol has been extensively discussed in Ref. [8,31,50] To this end, we focus on the Hamitonian Starting from the idle configuration at t = τ 0 = 0, we tune the three coupling strengths χ j (t) (with j = x, y, z) along the closed loop Γ in the three-dimensional parameter space shown Fig. 1(b). The loop consists of six straight lines connecting six key configurations reached at times t = τ k (with k = 0, . . . , 6), see Tab. 1 for more details. At the end of the protocol t = T = τ 6 , the system goes back to the original idle configuration, i.e. χ j (τ 6 ) = χ j (τ 0 ). In between two subsequent configurations, two parameters are kept fixed while the third one smoothly interpolates between 0 and its maximum value (or viceversa). For the sake of concreteness, in what follows, we consider the interpolating function The closed loop Γ(t) followed by the three coupling strengths χ x , χ y and χ z specified in Eq. (7) during the clockwise braiding of MZMs γ 1 and γ 2 . The "anchor points" Γ(τ j ), which are listed in Tab. 1, are highlighted with small bullets.
with g(0) = 0 and g(1) = 1. This protocol clearly satisfies the (above-mentioned) sufficient condition which ensures the preservation of the qubit computational space. We can therefore study the adiabatic time evolution of a generic initial state |i = α|0 + β|1 . Once the system is back in the idle configuration, at time t = T , the final state |f must be related to the initial one by a U (2) transformation |f (T ) = U Γ |i . Because of the degeneracy of the computational space, the dynamical phase picked by |0 during the time evolution is the same as the one picked by |1 . Therefore, a non-trivial U Γ can only emerge as a non-Abelian Berry phase, which depends on the geometrical properties of the loop Γ. Notably, it is independent of the specific values of τ j and the interpolating function g. Moreover, the conservation of [P 34 , H 1↔2 (t)] = 0 and [P, H 1↔2 (t)] = 0 implies that U Γ is diagonal in the basis {|0 , |1 }. The Berry phase picked by the state |0 (|1 ) equals (minus) the solid angle Ω Γ enclosed by the loop Γ [8,31,50]. For the loop depicted in Fig. 1(b), up to an overall phase, this corresponds to the unitary matrix Therefore, the holonomic scheme we briefly reviewed corresponds to a quantum phase gate (also known as π/4 gate). This holonomic protocol is completely equivalent to the physical clockwise braiding of γ 1 and γ 2 . This can be easily understood by tracking the motion of the zero-energy modes, initially associated with γ 1 and γ 2 , throughout the evolution of the system along Γ [31,51]. It is therefore not surprising that this holonomic protocol inherits the same topological protection Table 1: Values of the coupling strengths χ i (t) at the key configurations t = τ α (with α = 0, . . . , 6) along the closed loop shown in Fig. 1 featured by the physical braiding of MZMs. A single topological MZM can be completely decoupled from other Majoranas with exponential accuracy, for instance, by varying its distance from the other Majoranas. This has two important consequences: (i) It guarantees the exponential protection of the degeneracy of the computational space, thus preventing the onset of unwanted dynamical phase differences. (ii) It exponentially confines the loop Γ on the three coordinate planes χ j = 0. Therefore, the enclosed solid angle is always Ω Γ = π/2, regardless of deviations from the ideal "cubic" shape of the loop Γ in Fig. 1(b) which may arise due to a limited control on the non-vanishing couplings χ k [8,31,50]. Hence, the unitary operation implemented by the holonomic protocol can approach U Γ in Eq. (9) to exponential accuracy.
With respect to the physical braiding, the advantage of the holonomic approach is that it can be easily generalized to different classes of loops and/or to different parameter spaces. In the first case, for example, by halving the solid angle enclosed by the loop Γ, it is possible to implement a T gate (also known as π/8 gate). While the full topological protection is lost, this implementation of the T gate can still take advantage of geometrical protection. Indeed, perturbations on the loop which do not change the enclosed solid angle Ω do not affect the gate outcome [50]. In the next section, we explain how two Majorana qubits can be entangled by holonomic schemes.

Entangling gates
We consider a two-qubit system, consisting of two copies of the five-point star architecture as shown in Fig. 2(a). It features twelve different Majorana modes γ α j , with j = 0, . . . 5 and the qubit index α = A, B. In complete analogy with the single-qubit case, we define inter-coupling pairings P α jk = iγ α j γ α k and introduce the idle configuration, which features two non-vanishing Fig. 2(a)]. Moreover, without loss of generality, we assume each qubit to obey a global even fermion parity, i.e.
The computational space of the whole system in the idle configuration is thus spanned by four degenerate ground states where the single-qubit even-parity states |0 α and |1 α are defined as in Eq. (6).
While holonomic single-qubit gates can be easily implemented by manipulating the intraqubit couplings P A ij , two-qubit gates require the manipulation of terms in the Hamiltonian which act on both qubits, such as the inter-qubit couplings I jk = iγ A i γ B j . Importantly, these terms do not commute with the fermion parity of each qubit, i.e. [P α , I jk ] = 0. Therefore, they can potentially push the two-qubit system out of its four-dimensional computational space, spanned by the even-even states in Eq. (10). This is precisely what happens if we braid Majoranas γ A 1 and γ B 2 , for example by holonomically tuning parameters P B 05 , P B 02 , and I 10 . Such an operation would entangle the two subsystems A and B but it would also induce leakages out of the two-qubit computational space. This specific example is a consequence of the more general non-entangling rule, introduced by Bravyi in 2006 [22] stating that no entangling gates between distinct qubits can be realized just by braiding Ising anyons. For the sake of completeness, we acknowledge that it is possible to entangle a couple of Majorana qubits, defined in a peculiar way such that they share a pair of Ising anyons, just by using Majorana braiding [52]. However, this approach is not scalable since the image of the braiding group is a nontrivial subgroup of the Clifford group for n ≥ 3 qubits, meaning that there exist Clifford gates which cannot be realized only by braiding [40].
To overcome the limitations posed by the non-entangling rule, it is necessary to devise novel holonomic protocols that go beyond the simple braiding of MZMs and preserve the fermion parity P α of each individual qubit. To this end, it is possible to follow two different approaches. One possibility is to manipulate more complicated operators that act on both qubits while still commuting with each P α . Alternatively, it is possible to ensure the conservation of fermion parity by external means, for instance by making unfavorable single-electron tunnelings between the qubits. In what follows, we propose two protocols which are based on the first and second approach, respectively.

The 4γ protocol
To devise an holonomic procedure which fully entangles the two qubits, we introduce the operator It represents an interacting term which involves the two qubits but preserves their individual parities [O 4γ , P α ] = 0. We argue that it allows us to realize a holonomic entangling gate by considering the time-dependent Hamiltonian where the coupling strengths χ j (t) adiabatically follow the loop Γ depicted in Fig. 1(b) and described in Tab. 1. Before presenting a detailed and rigorous analysis of this protocol, it is useful to develop some intuition about its entangling capabilities. To this end, we notice that the Hamiltonian H 4γ closely resembles the one which would braid Majoranas γ B 1 and γ B 2 within the same qubit [see Eq. (7)]. The key difference is the presence of the operator O 4γ , which replaces the simple intra-qubit coupling P B 01 with the product P A 12 P B 01 . This allows to effectively implement a controlled braiding. Indeed, the parity of Majoranas γ A 1 and γ A 2 on the control qubit A, controls the direction of the braiding of γ B 0 and γ B 1 on the target qubit B: for P A 12 = ±1, the effective braiding on qubit B is either clockwise (+1) or anti-clockwise (−1), as shown in Fig. 2(b) with solid purple and dotted green lines, respectively. This justifies the entangling capability of the procedure. In order to fully characterize the holonomic protocol, it is necessary to study the adiabatic time evolution of the four even-even ground states in Eq. (10). Initially, at time t = 0, we assume the system to be in one of the states |ψ nm (t = 0) = |n A m B , with n, m ∈ {0, 1}. Since the parity operators P A 34 and P B 34 are conserved throughout the whole time evolution, they allow us to conveniently label the states at every time t, that is Importantly, the four states |ψ nm (t) are always degenerate. This can be proven by identifying zero energy operators which allow us to transform between these four states. For 0 ≤ t ≤ τ 2 (see Tab. 1 and Fig. 1), we find that represent these operators as they commute with the Hamiltonian [H(t), R A ] = [H(t), R B ] = 0 and anticommute with the respective conserved quantity {P α 34 , R α } = 0. Analogously, zero energy operators can also be found for the following steps, i.e. for τ 2 ≤ t ≤ τ 6 = T . The energy spectrum of the system, obtained by diagonalizing the full Hamiltonian H 4γ (t), confirms the degeneracy of states |ψ nm and the presence of a finite energy gap, which separates them from excited states [see Fig. 2(c)].
The adiabatic theorem, together with the conservation of the parities P α 34 , allows us to relate the final states |ψ nm (T ) to the initial ones via a U (4) diagonal matrix U ent Γ . Up to a global phase, this transformation only depends on the geometric properties of the loop Γ in the parameter space. In order to find U ent Γ , we compute the Berry curvature associated with each state (see App. A.1 for more details) Therefore, analogously to the simple holonomic braiding, the Berry phases picked up by each element of the basis in Eq. (10) only depends on the solid angle enclosed by the closed loop Γ in parameter space. In particular, we get which, up to a global phase, corresponds to the unitary transformation The entangling power of this gate is EP(U ent Γ ) = 2/9, which means that it is able to fully entangle the two qubits [53]. For instance, starting from the product state |i = (|0 A + |1 A )⊗ (|0 B + |1 B ) /2, the holonomic procedure generates the maximally entangled state These results have also been confirmed numerically, by simulating the adiabatic time evolution of the system with the QuTip package [54,55] [see, for example, Fig. 4(a)]. By applying additional single-qubit Hadamard (H) and phase (S) gates, which can be implemented using, for example, holonomic braiding as discussed in Sec. 2.1, it is possible to obtain the control-Z gate as well as the CNOT gate CNOT = (I ⊗ H)CZ(I ⊗ H).
We mention in passing that the unitary operator U ent Γ can also be obtained by using a measurement-only approach [34]. The latter would require four subsequent projective forced measurements of the same parity operators, which we tune in H 4γ . That is where Π + P = (1 + P )/2 is the projector on the eigenstate of the parity operator P with eigenvalues +1. While our holonomic protocol does not rely on the realization of projective (joint) parity measurements, it poses (at least) two important experimental challenges: (i) The implementation of a tunable four-Majorana interaction term O 4γ and (ii) the detrimental effect of parasitic couplings. The latter corresponds to unwanted pair-wise couplings, which can appear during the manipulation of the Majoranas. For example, at time t = τ 1 , Majorana Fig. 2(a)] and it simultaneously interacts with γ A 2 , γ A 1 and γ B 1 [blue area in Fig. 2(a)]. In this configuration, it might be difficult to prevent the onset of additional small parasitic couplings between the involved Majoranas, such as Unfortunately, these terms do not conserve the parity of individual qubits and/or split the ground state degeneracy, spoiling the holonomic entangling procedure. In the next section, we propose an alternative holonomic approach that overcomes those two aforementioned drawbacks.

The double tunneling protocol
In this section, we take advantage of the framework developed for the study of the 4γ protocol and devise a new holonomic scheme that features high resilience against parasitic couplings and only requires to tune pair-wise Majorana couplings. As discussed before, such an approach necessarily relies on external means to preserve the parity of each individual qubit. To this end, we consider the additional time-independent Hamiltonian which makes single-electron tunneling between the two qubits energetically unfavorable (in the regime of a significant magnitude of E). Importantly, such a Hamiltonian naturally emerges in systems where each qubit consists of a floating superconducting island with a finite charging energy. This feature has previously been exploited to protect qubits from quasiparticle poisoning and to allow joint-parity measurements of more than two Majoranas [31,35,36,46]. The presence of a finite H E allows us to trade the interacting term O 4γ = I 20 I 11 [see Eq. (11)] with a simpler and more feasible one even though the latter does not commute with the qubit parities P α . In order to develop some intuition, we focus at first on the limit of large E so that the individual inter-qubit tunnelings I 20 and I 11 are highly suppressed. In this limit, only virtual cotunneling processes can happen such as whose effect is analogous to the interacting term O 4γ . Hence, the time-dependent Hamiltonian is likely to function in a similar way as previously discussed for H 4γ . This is confirmed by the careful analysis of the protocol H dt (t) which we carry out below. Importantly, we argue that it is not necessary to work in the limit of large E, as the holonomic entangling scheme can be successfully implemented even for a finite E ∼ Θ (see Fig. 4).
We characterize the protocol based on H dt (t) along the lines of the previous section. In particular, we observe that the parity operators P A 34 and P B 34 are still conserved and that operators like R A and R B , transforming between the states of the computational space, are still at . This guarantees the degeneracy of the four groundstates |ψ nm (t) throughout the whole protocol, as nicely confirmed by the spectrum plotted in Fig. 3(b) and obtained by the exact diagonalization of the full Hamiltonian H dt (t). Hence, after the adiabatic evolution of the system, the final states |ψ nm (T ) are related to the initial ones |ψ nm (0) via a U (4) diagonal matrix which, up to a global phase, depends only on the geometrical properties of the closed loop Γ. In order to determine this unitary transformation, we compute the Berry phase picked up by each state by integrating the Berry connection along Γ (see App. A.2). Even though the Berry curvature differs from the simple form in Eq. (16), we find that, up to a global phase, the implementation of the protocol based on Hamiltonian H dt (t) results in the same unitary transformation U ent Γ we obtained within the 4γ protocol. Importantly, this result does not depend on E as long as the latter is finite. The protocol is, therefore, able to maximally entangle the two qubits and it can be straightforwardly turned into a CNOT gate by adding single-qubit Clifford gates according to Eq. (21).
All these results have been confirmed numerically by simulating the time evolution of the system and testing the validity of Eq. (19) [see Fig. 4(a)]. Importantly, numerical simulations represent also a valuable tool to fully characterize the protocol by inspecting its robustness against non-adiabatic effects, parasitic couplings, and poor control of the tuning parameters. These aspects, which are of fundamental importance when it comes to practical implementations of the proposed entangling scheme, are carefully discussed in the next section.

Robustness of the entangling protocol
So far, we have studied the entangling protocol based on the Hamiltonian H dt (t) under ideal conditions. Indeed, we have considered the time evolution of the system while the three control parameters χ j (t) are adiabatically and precisely tuned along the loop Γ. This raises the question to which degree the resulting unitary operation U ent Γ is robust with respect to non-ideal effects, which arise under more realistic conditions.

Adiabaticity
Strictly speaking, the existence of a finite energy gap ∆E between the ground state manifold and the excited states guarantees the applicability of the adiabatic theorem only in the limit T → ∞. Protocols with a shorter time duration T might indeed feature diabatic transitions (e.g. of Landau-Zener type) which push the system out of computational space and spoil the holonomic quantum gate.
In order to analyze and quantify this effect, we numerically simulate the time evolution of the initial state |i = (|0 A + |1 A ) ⊗ (|0 B + |1 B ) /2 for several different durations T of the whole protocol. The final states |f (T ) can then be compared with the expected result |f (∞) = U ent Γ |i [see Eq. (19)] by computing the overlaps | f (∞)|f (T ) | 2 . The latter are plotted in Fig. 4(a) for E = Θ and show that the time evolution of the system behaves adiabatically for T 100∆ −1 . Faster implementations of the protocol would result in |f (T ) = |f (∞) . The oscillations featured by the overlap | f (∞)|f (T ) | 2 for short T can be understood as interference patterns resulting from subsequent transitions between ground and excited states, in analogy with the Landau-Zener-Stückelberg effect [51,56].
Let us define the "adiabatic threshold" T ad as the minimal duration of the holonomic process for which diabatic effects become negligible, say | f (∞)|f (T ) | 2 < 10 −3 for all T ≥ T ad .
This time scale clearly depends on the spectrum of the system during the time evolution which, in turn, depends on the energy E and on the maximum coupling strengths Θ i featured by the Γ loop specified in Tab. 1. In Fig. 4(b), we numerically compute T ad as a function of E for three different values of Θ y while keeping Θ x = Θ fixed. Notably, both large and small values of E are detrimental for the adiabaticity of the protocol: For a given Θ y , the smallest values of T ad are actually reached for E of the order of E ∼ Θ y . This can be understood by observing that, for E, Θ y < (1 + √ 2)Θ, the smallest energy gap ∆E throughout the loop Γ is reached for t = τ 2 and it reads This gap is maximized for Θ y = E, which is qualitatively consistent with the fact that T ad is minimized for E ∼ Θ y . The lack of a quantitative agreement stems from the fact that the energy spectrum of the system [see Fig. 3(b)] features a more complicated structure with respect to simple Landau-Zener transitions.
In any case, Fig. 4(b) shows that, for a wide range of parameters, T ad lays between 100Θ −1 and 200Θ −1 . A reasonable estimation for the coupling strengths is Θ ∼ 10 GHz [32,35,57] which corresponds to a reasonable timescale T ad ∼ 10 ns. Longer durations result in even smaller diabatic errors, which are known to be polynomially suppressed in T [57]. In particular, we find a suppression approximately proportional to T −4 , as shown in the inset of Fig. 4(a). In complete analogy with standard braiding schemes [57], depending on the required gate fidelity, the assessment and mitigation of diabatic errors might play an important role in realistic implementations of our entangling scheme.

Finite tuning accuracy and parasitic couplings
A finite tuning accuracy of parameters χ j (t) results in deviations from the loop Γ specified in Tab. 1. Importantly, in complete analogy with standard holonomic braiding schemes [8,31], both our 4γ and double tunneling protocols are robust with respect to errors on the control of χ j (t). This is due to two important reasons: (i) The topological nature of the Majorana qubits guarantees the existence of parameters (e.g. the separation lengths) which exponentially suppress the coupling/interactions between the MZMs. As a result, the loops in parameter space can be confined to the three coordinate planes χ j = 0 (j = x, y, z) with exponential accuracy. (ii) On each of these coordinate planes, the Berry curvature associated with the holonomic protocols has no components perpendicular to the plane itself. This statement is proven in App. A.2 and it can be readily seen from Eq. (16) for the 4γ protocol. Poor control over deviations from Γ within each coordinate plane has, therefore, no influence on the resulting unitary operation U ent Γ . Finally, let us discuss the robustness with respect to possible parasitic couplings. In contrast to the 4γ protocol, the double tunneling scheme only features three of them. In particular, for τ 0 < t < τ 2 and τ 2 < t < τ 4 , the only unwanted couplings between the involved Majoranas are the inter-qubit terms I 25 and I 22 , respectively. As for the last stage of the protocol τ 4 < t < τ 6 , the only possible parasitic coupling is the intra-qubit term P B 25 . For the sake of clarity, these three parasitic couplings are shown with dashed-dotted lines in Fig.  3(a).
Importantly, the presence of those three parasitic coupling is not detrimental for the implementation of the holonomic entangling gate. Let us focus, for simplicity, on the first stage of the protocol τ 0 < t < τ 2 . In this case, the only possible parasitic coupling is I 25 .
As it commutes with both the zero energy operators R A and R B , the degeneracy of the four ground states is preserved. The same holds for the other two stages of the protocol, i.e. τ 2 < t < τ 4 and τ 4 < t < τ 6 . Moreover, we numerically verify that the presence of any of the three parasitic couplings does not modify the geometrical phases acquired by the qubits state throughout the holonomic protocol. To this end, we simulate the time evolution of the initial state |i according to the Hamiltonian H pc (t) = H dt (t) + v 1 (t)I 25 + v 2 (t)I 22 + v 3 (t)P B 25 (28) where the parameters v j (t) evolve according to v j (t) = and control the strength of the three parasitic couplings. In the adiabatic regime, we verify that the final state still satisfies |f (T ) = |f (∞) , independently of the maximum value V acquired by v j (t). As an example, in Fig. 4(a), we plot the overlap | f (∞)|f (T ) | 2 computed with V = Θ/2 (orange thin line).

Conclusions
Two-qubit entangling gates, such as the CNOT gate, are essential building blocks for quantum computations. In the realm of topological Majorana qubits, these gates are almost exclusively considered within the framework of measurement-based topological quantum computation. Instead, we provide a complementary approach and propose two holonomic entangling protocols. They represent non-trivial extensions to the two-qubit case of a known technique to implement (topological) single-qubit gates. With respect to the corresponding measurement-based approaches, our holonomic protocols do not require the implementations of forced and projective measurements of (joint) Majorana parities. We believe that this might be an advantage, especially since efficient and high-fidelity measurement schemes are still lacking. We fully characterize the protocols and prove the high degree of robustness of the double tunneling scheme. Importantly, the latter only requires the ability to tune pair-wise Majorana couplings within qubits characterized by a finite charging energy. Those ingredients are actually featured by several proposed setups for TQC, such as Josephson junction arrays [31,32], Majorana box qubits [36] and the tetron/hexon schemes [35], thus increasing the potential experimental relevance of our proposal. belong to. They read The four bases we have chosen are n = 0; m = 0 → (40) Figure 5: Functions F(R, θ) (in blue) and ξ(R, θ) (in red) for R = .
the planes themselves. This guarantees that the double tunneling protocol, despite a more involved structure of the Berry connection, is still topologically protected. To find the unitary transformation U ent Γ associated with the holonomic protocol, we compute the non-Abelian Berry phase associated with the Γ loop by integrating the Berry curvature on a surface enclosed by Γ (or, equivalently, by integrating the Berry connection A nm = i ψ nm ( χ)|∇ χ |ψ nm ( χ) on the loop Γ).