Majorana Braiding in Realistic Nanowire Y-Junctions and Tuning Forks

TL;DR

This study numerically investigates Majorana braiding in realistic nanowire Y-junctions and tuning forks, identifying optimal geometries and parameters that enhance braiding success probability for quantum computing applications.

Contribution

It provides the first detailed numerical analysis of braiding in realistic nanowire geometries, highlighting optimal angles and the advantages of tuning fork designs.

Findings

Optimal Y-junction half-angle around 20 degrees.

Tuning fork geometry reduces dynamical phase issues.

Braiding success is feasible with current experimental setups.

Abstract

Majorana fermions are predicted to arise at the ends of nanowire devices which combine superconductivity, strong spin-orbit coupling and an external magnetic field. By manipulating networks of these devices with suitable gating, it has been suggested that braiding operations may be performed which act as logic operations, suitable for quantum computation. However, the unavoidable misalignment of the magnetic field in any realistic device geometry has raised questions about the feasibility of such braiding. In this paper, we numerically simulate braiding operations in devices with Y-junction and tuning fork geometries using an experimentally motivated nanowire model. We study how the static and dynamical features vary with geometric parameters and identify parameter choices that optimise the probability of a successful braid. Notably, we find that there is an optimal Y-junction half-angle (about 20 degrees for our parameter values), which balances two competing mechanisms that reduce the energy gap to excitations. In addition, we find that a tuning fork geometry has significant advantages over a Y-junction geometry, as it substantially reduces the effect of dynamical phase oscillations that complicate the braiding process. Our results suggest that performing a successful braid is in principle possible with such devices, and lies within experimental reach.

Equations

$${\gamma }_{2j-1}=c_j^{†}+c_j^{},$$

$${\hat{H}}_{\mathrm{K}}=-\mu \sum _j\left[c_j^{†}{c}_j^{}\right]-\frac{1}{2}\sum _j\left[t{c}_j^{†}{c}_{j+1}^{}+\Delta e^{i\varphi }{c}_j^{}{c}_{j+1}^{}+H.c.\right].$$

$${\hat{H}}_{\mathrm{NW}}=\left(-\frac{{\hslash }^2}{2m^{\ast }}{\partial }_x^2-i{\alpha }_R{\partial }_x{\sigma }_y-\mu \right){\tau }_z+V_Z{\sigma }_x+{\Delta }_0{\tau }_x,$$

$$V_Z$$

$$|n_1=0,n_2=0⟩=c_1{c}_2|\Omega ⟩,$$

$$c_1$$

$$c_2$$

$$|n_1=1,n_2=1⟩=(-)c_1^{†}{c}_2^{†}|\Omega ⟩.$$

$${\hat{H}}_0$$

$$t(a)=\frac{{\hslash }^2}{2m^{\ast }{a}^2},$$

$${\partial }_x^2$$

$${\partial }_x{\sigma }_y$$

$${\hat{H}}_{\theta }$$

$$\mu (x,t)$$

$$r(x)$$

$${\theta }_c=\arcsin \left[{\Delta }_0/V_Z\right]\approx 32.2^{\circ },$$

$${\hat{H}}_{12}$$

$$U(t)$$

$$|\psi (0)⟩$$

$$|\psi (t)⟩$$

$$\displaystyle\left|\big{\langle}1_{-}\big{|}\psi(T)\big{\rangle}\right|^{2}$$

$$\displaystyle\left|\big{\langle}1_{+}\big{|}\psi(T)\big{\rangle}\right|^{2}$$

$$\left|{\psi(T/2)}\right\rangle=\frac{1}{2}\bigg{[}\left(\left|{1_{-}}\right\rangle+\left|{1_{+}}\right\rangle\right)+e^{i\phi}\left(\left|{2_{-}}\right\rangle+\left|{2_{+}}\right\rangle\right)\bigg{]},$$

$$|\psi (T/2+t)⟩$$

$${\Delta }_0(B)$$

$$V_{\mathrm{g}}(x,y,d;x_j,U_j)$$

$$g(u,v)$$

$$R(u,v)$$

$$V_{\max }$$

$$V_{\min }$$

$$\Delta (W,M)$$

Full text

I Introduction

Bound states of Majorana fermions are believed to exhibit non-abelian statistics Alicea (2012); Stanescu and Tewari (2013); Beenakker (2013); Sarma et al. (2015); Elliott and Franz (2015); Aasen et al. (2016); Aguado (2017); Lutchyn et al. (2018); Zhang et al. (2019a), and so offer the exciting possibility of realising a topological quantum computer Nayak et al. (2008); Kitaev (2003). This ambition, coupled with experimental advances, has focussed a great deal of recent effort towards realising emergent Majorana quasiparticles in condensed matter systems. Notably, suitable Majorana bound states are predicted to arise in certain fractional quantum Hall states Moore and Read (1991), in $p$-wave superconductors in one and two dimensions Read and Green (2000); Kitaev (2001), and in many other low-temperature solid state systems Sau et al. (2010). Beyond this, many varied theoretical proposals have been suggested which would allow the non-abelian properties of Majorana modes to be leveraged into a topological qubit Aasen et al. (2016). While compelling experimental evidence for the existence of Majorana modes has been detected in many of these settings (and in particular in nanowires Mourik et al. (2012); Das et al. (2012); Deng et al. (2016); Chen et al. (2017); Nichele et al. (2017); Deng et al. (2018); Gül et al. (2018); Shen et al. (2018); O’Farrell et al. (2018); Vaitiekėnas et al. (2018); Grivnin et al. (2018); Zhang et al. (2018); van Zanten et al. (2019)), a completely unambiguous signature of Majorana modes remains lacking. In this work, we study the braiding of a pair of Majorana fermions in a realistic nanowire device, which, if reproduced experimentally, could provide one such unambiguous signature.

The area of Majorana nanowires is now fairly mature Aguado (2017); Lutchyn et al. (2018); Zhang et al. (2019a), but much of our understanding of such systems stems from Kitaev’s toy model of a $p$-wave superconducting chain Kitaev (2001). In the topological regime of this model, the ends of the wire host a pair of localised Majorana zero modes (MZMs), which may be moved around by adjusting the local system parameters. By forming networks of such wires, non-abelian braiding operations may be performed. The Kitaev model, however, requires exotic superconducting pairing, and is therefore difficult to realise in a real material. Instead, most experiments are believed to approximate the more realistic nanowire models of Refs. Sau et al. (2010); Oreg et al. (2010); Lutchyn et al. (2010). In these systems, strong spin-orbit coupling, $s$-wave superconducting pairing, and a moderately strong magnetic field conspire to produce the conditions necessary for MZMs.

Nanowires based on this model have now been developed by several groups, both using direct epitaxial growth and by depleting regions of a two-dimensional electron gas to leave an effectively one-dimensional (1D) channel Mourik et al. (2012); Das et al. (2012); Deng et al. (2016); Chen et al. (2017); Nichele et al. (2017); Deng et al. (2018); Gül et al. (2018); Shen et al. (2018); O’Farrell et al. (2018); Vaitiekėnas et al. (2018); Grivnin et al. (2018); Zhang et al. (2018); van Zanten et al. (2019). Many of these experiments have reported transport signatures consistent with the existence of MZMs: notably, a MZM should lead to a robust zero-bias conductance peak quantised to $2e^{2}/h$ Flensberg (2010); Law et al. (2009) (although disorder and finite-temperature effects may disguise this Wimmer et al. (2011)). However, such a signature is not ‘smoking gun’ evidence of a MZM: similar peaks can also be caused by spurious Andreev bound states, non-topological states close to zero energy which form due to local variations in the device parameters Liu et al. (2017); Hell et al. (2018); Moore et al. (2018); Liu et al. (2018); Reeg et al. (2018); Aseev et al. (2019). It has been suggested that more convincing evidence for Majorana modes could be obtained through interference experiments Akhmerov et al. (2009); Fu and Kane (2009); Fu (2010); Sau et al. (2015), by performing simultaneous tunnelling into the device at each end of the wire Liu et al. (2017), or by studying spin-dependent transport signatures Ricco et al. (2018). Ultimately, however, the most direct way to probe the non-abelian nature of MZMs would be to perform a logic operation such as a braid. Such operations could form the basis of a topological quantum computer (albeit one that is not, on its own, universal Bravyi (2006)).

Experimentally, we might envisage performing a braid following the proposal of Ref. Alicea et al. (2011), where a series of side gates are attached to a three-legged device (or ‘Y-junction’), as in Fig. 1. By adjusting the voltage of the side gates, the local chemical potential may be changed, and different regions of the wire can be made topological or trivial as desired. If there are initially two pairs of MZMs, one Majorana from each pair can be braided following the protocol indicated in Fig. 1. Although the non-abelian nature of this operation has been confirmed in the ideal case Alicea et al. (2011), there has been some doubt about its feasibility in a real material system. One concern is that the magnetic field (directed along the mother branch of a ‘Y-junction’) is not parallel to the wire direction in the prongs of the device, which is known to reduce the size of the bulk gap Mourik et al. (2012); Osca et al. (2014); Rex and Sudbø (2014); Nijholt and Akhmerov (2016). In addition, the effective superconducting pairing in the non-parallel legs can lead to the formation of a $\pi$-junction Alicea et al. (2011). Beyond this, decoherence and quasiparticle poisoning times may be too short to allow the successful completion of a braid Rainis and Loss (2012); Higginbotham et al. (2015); Albrecht et al. (2017); Zhang et al. (2019b); Clarke et al. (2017).

In this paper, we numerically study the feasibility of such a braid in a realistic nanowire system. In previous works, braiding of Majoranas has been performed successfully in numerical simulations of the Kitaev model Amorim et al. (2015); Sekania et al. (2017), while braiding of non-Majorana defects has also been carried out in simulations of the Su-Schrieffer-Heeger model Boross et al. (2019). In this work, we instead focus on the more realistic continuum nanowire model Sau et al. (2010); Oreg et al. (2010); Lutchyn et al. (2010), which allows us to take into account the misalignment of the external magnetic field and the device geometry directly. As our starting point, we use parameters from existing state-of-the-art simulations of nanowire devices, which have been shown to reproduce experimental conductance data extremely well Liu et al. (2017); Zhang et al. (2018). However, motivated by an ongoing experimental collaboration, we also incorporate material parameters and system geometries that are particularly relevant to devices grown using IBM’s template-assisted selective epitaxy (TASE) technique Borg et al. (2014); Czornomaz et al. (2015); Gooth et al. (2017a); Schmid et al. (2015); Gooth et al. (2017b). While we do not claim to draw categorical conclusions about braiding in a specific (existing or proposed) device, we believe our results provide qualitative and approximate quantitative statements about braiding in devices which are within current experimental reach, and which overall are encouraging.

The structure of this paper is as follows. In Sec. II we introduce the nanowire model we will be using throughout this paper and discuss the ideal braiding properties of Majorana fermions. In Sec. III we study the static (instantaneous) properties of a nanowire Y-junction as the braid is performed, and identify geometric parameters which optimise the important energy scales of the system. In Sec. IV we study the dynamical properties of the system by performing numerical braiding simulations. In particular, we consider how the braiding success varies as parameters and time scales are modified. In Sec. V we discuss how using a tuning fork geometry offers several advantages over a Y-junction geometry, and investigate the robustness of our conclusions to variations in the underlying parameters. Finally, in Sec. VI we summarise our results and provide some concluding remarks.

II Majorana Fermions in Nanowires

II.1 Nanowire Models

Majorana fermions arise in condensed matter systems as emergent fermionic quasiparticles which are their own ‘antiparticle’. Explicitly, Majorana fermions associated with operators $\gamma_{i}$ satisfy the relations $\gamma_{i}^{\dagger}=\gamma_{i}^{\phantom{\dagger}}$ and $\left\{\gamma_{i},\gamma_{j}\right\}=2\delta_{ij}$. They can always be defined (at least mathematically) by combining complex fermion creation and annihilation operators ($c^{\dagger}_{j}$ and $c^{\phantom{\dagger}}_{j}$) through the relations

[TABLE]

and in this way may be thought of as fractional excitations of the underlying electrons of a system. However, the transformation given above is not usually very useful—the resulting Majorana operators are often entangled in space and do not correspond to physically meaningful excitations of the system. Nevertheless, in the presence of superconductivity (where the particle number is only conserved modulo two and the quasiparticles are coherent superpositions of particles and holes) Majorana modes can arise as natural excitations. In particular, Majorana modes that are pinned to zero energy and which are spatially well-separated can arise as topological degrees of freedom.

The simplest 1D model which exhibits MZMs is the previously mentioned Kitaev chain Kitaev (2001), which has the Hamiltonian

[TABLE]

In this Hamiltonian, $c^{\dagger}_{j}$ creates a spinless fermion on site $j$, $\mu$ is the chemical potential, $t>0$ describes the intersite hopping, and $\Delta e^{i\phi}$ describes the $p$-wave superconducting pairing. When written in terms of Majorana fermions, this Hamiltonian describes pairing between Majorana fermions on the same site and between Majorana fermions from neighbouring sites. If $\left|\mu\right|<t$, the system is in the topological regime, and an unpaired MZM arises at each end of the chain. 111We note that Eq. (2) is also related to the XY spin chain through the Jordan-Wigner transformation—see, for example, Ref. Derzhko, 2008.

Unfortunately, $p$-wave superconductivity is difficult to realise in the laboratory, but more realistic models which support MZMs were introduced in Refs. Sau et al. (2010); Oreg et al. (2010); Lutchyn et al. (2010). In these cases, the underlying system is a 1D semiconducting nanowire with strong Rashba spin-orbit coupling, usually corresponding experimentally to InAs or InSb. The nanowire is then proximity coupled to an $s$-wave superconductor, achieved experimentally by depositing a material such as Al as a thin layer on top of the nanowire, as shown schematically in Fig. 1(c). In zero magnetic field, a system of this kind is a trivial quasi-1D superconductor. However, if a magnetic field aligned with the nanowire axis is increased in strength, the bulk gap may close and reopen in a topological phase transition, resulting in a topological superconductor (TSC) with MZMs at its ends. This process may be thought of as a competition between the spin-orbit coupling, the external field, and the superconducting pairing, which conspire to produce the effective $p$-wave pairing required for the Kitaev model.

The 1D Hamiltonian for the continuum nanowire model may be written in BdG form as

[TABLE]

which acts on the Nambu spinor $\hat{\psi}_{x}=\left(c^{\phantom{\dagger}}_{\uparrow x},c^{\phantom{\dagger}}_{\downarrow x},c^{\dagger}_{\uparrow x},-c^{\dagger}_{\downarrow x}\right)^{T}$. In this expression, the nanowire extends along the $x$-direction, $m^{*}$ is the effective mass of the nanowire material, $\alpha_{R}$ is the Rashba spin-orbit coupling parameter, $V_{Z}$ is the Zeeman term, and $\Delta_{0}$ is the $s$-wave pairing gap induced by proximity. The Pauli matrices $\{\tau_{x},\tau_{y},\tau_{z}\}$ act in particle-hole space, while $\{\sigma_{x},\sigma_{y},\sigma_{z}\}$ act in spin space. The model is in the topological regime when the topological criterion,

[TABLE]

is satisfied Sau et al. (2010); Oreg et al. (2010); Lutchyn et al. (2010). Deep in the topological phase, the Hamiltonian can be shown to reduce to that of Eq. (2) Alicea et al. (2011).

The continuum nanowire model has been used as a starting point for many theoretical and numerical studies of devices believed to host MZMs (see, for example, Refs. Liu et al. (2017, 2019)). By incorporating components such as leads, by adjusting the dimensionality, or by adding new terms to the Hamiltonian, a variety of different experimental setups can be simulated, and the effects on the stability and behaviour of the resulting Majorana modes studied. Notably, the continuum nanowire model has been used to produce numerical conductance simulations which agree extremely well with several state-of-the-art experimental measurements Liu et al. (2017); Zhang et al. (2018). In this way, numerical simulations have provided important supporting evidence for the existence of Majorana modes in real nanowire devices.

II.2 Majorana Braiding

Majorana fermions, being examples of non-abelian Ising anyons, are useful for quantum computing due to their ability to store information nonlocally. In the ideal case, when a pair of localised Majorana fermions are brought together, one of two things may happen: the pair may annihilate, or the pair may fuse to form an ordinary complex fermion. However, when the Majorana fermions are well separated, there is no local measurement which can be made to distinguish which outcome (or fusion channel) will result from this process. Instead, the outcome depends on the topology of the trajectories that the Majorana fermions have followed up to this point. A benefit of this nonlocal encoding of the state is that local perturbations, which are prevalent in any realistic experimental system, cannot easily affect this topology, which endows the state with an inherent robustness. On the other hand, to be useful for quantum computing, the state will need to remain coherent over timescales long enough to perform a topological operation.

In this paper, we focus on one such operation: a braid between two Majorana fermions (i.e. two sequential exchanges) which alternates between the two fusion channels. Although the braid only directly involves two Majorana modes, we require a set of four MZMs in total in order to alter the fusion channel, as the particle-hole symmetry requires the total fermion parity to be conserved. Abstractly, we can imagine creating two pairs of MZMs from the vacuum ($\{\gamma_{1},\gamma_{2}\}$ and $\{\gamma_{3},\gamma_{4}\}$), as shown at the top of Fig. 1. At this stage, if we were to bring the two MZMs from each pair back together they would annihilate. We can therefore label this state as

[TABLE]

where $n_{i}=c^{\dagger}_{i}c^{\phantom{\dagger}}_{i}$ is an occupation number, $\left|{\Omega}\right\rangle$ is the initial vacuum, and we have used the complex fermion operators $c_{1}=\frac{1}{2}\left(\gamma_{1}+i\gamma_{2}\right)$ and $c_{2}=\frac{1}{2}\left(\gamma_{3}+i\gamma_{4}\right)$. [Note, however, that there is a gauge freedom in this choice]. For later use, we define this initial state as the $\left|{0}\right\rangle$ state of a topological qubit.

To perform the braid, we adiabatically exchange the Majorana fermions $\gamma_{2}$ and $\gamma_{3}$ twice with the same orientation, so that each one returns to its initial position but with their trajectories intertwined, as indicated in Fig. 1. It may be shown that a braid of this form results in the transformation $\gamma_{2}\to-\gamma_{2}$, $\gamma_{3}\to-\gamma_{3}$ Alicea et al. (2011), and so

[TABLE]

in terms of complex fermion operators. In this way, after the braid the pairs $\{\gamma_{1},\gamma_{2}\}$ and $\{\gamma_{3},\gamma_{4}\}$ now each fuse to form a complex fermion, and the starting state has been transformed to

[TABLE]

If we identify this state as the $\left|{1}\right\rangle$ state of a topological qubit, then the braid may be interpreted as a Pauli $\sigma_{x}$ operation or a NOT gate. We observe that total parity is conserved throughout this process, as required by particle-hole symmetry.

The ideal case discussed above and illustrated in Fig. 1(a) considers an adiabatic braid involving only four degenerate Majorana fermions. In a real system, any MZMs will necessarily be part of a much larger spectrum of states, the MZMs themselves will not be at exactly zero energy, and any braiding process will take place over a finite amount of time. All of these factors can negatively affect a braid’s success.

In this paper, we will reproduce this braid in nanowire Y-junctions and tuning forks described by the Hamiltonian in Eq. (3), following the sequence of moves shown in Fig. 1(b). In this setup, we set the locations of the topological regions of the wire by altering the local chemical potential $\mu$ so that the topological criterion (Eq. (4)) is locally satisfied, forming MZMs at the boundary. By changing the local chemical potential as a function of time, we move the boundaries of the topological regions and, consequently, the positions of the MZMs. In contrast to the ideal case, a real nanowire has a continuum of bulk excited states: the braid must be performed slowly enough that the initial state does not mix with these excitations. On the other hand, the Majorana modes themselves will not be exactly at zero energy. Instead, the overlap between Majorana modes at different ends of the device will lead a small energy splitting, while the ideal braid assumes that the Majorana modes exist in a degenerate zero-energy subspace. Our braid must be performed fast enough that this assumption still effectively holds, and the dynamical phases that are introduced by this energy splitting must be carefully taken into account. We will study the effects of these considerations on the braiding process in the next few sections.

We note that while in our case we can study the success of a braid using numerical measures (such as wavefunction overlaps), evaluating a braid in a real device is a more challenging endeavour. Several techniques have been suggested to measure the fusion channel of a pair of Majorana fermions experimentally: In one proposal, the pair of Majoranas can be fused across a Josephson junction and the resulting Josephson current measured Alicea et al. (2011). Alternatively, a Josephson junction with a valve may be used as a ‘parity-to-charge converter’, which can either be used to detect the charge of the fusion channel directly or the channel can be measured indirectly through cyclic current measurements Aasen et al. (2016); Clarke et al. (2017).

III Static Properties of Majorana Y-Junctions

III.1 System Setup

To perform a braid, we require a device with at least three legs so that the Majorana fermions can remain well separated and can move around one another. For our purposes, we focus on two simple such designs: a ‘Y-junction’ and a ‘tuning fork’, as shown in Fig. 2. Nanowires can be grown with either of these geometries using TASE Gooth et al. (2017b). We assume that the three legs of each device have the same length $L$, and that the prongs of the $Y$-junction each make an angle of $\theta$ with the the direction of the magnetic field, which we take to lie in the positive $x$ direction. The tuning fork design is additionally labelled by the parameter $\eta$, which gives the proportion of each right-hand leg that is angled (before becoming horizontal again). We are particularly interested in studying how the feasibility of a braid varies with these geometric parameters.

We simulate these devices numerically by discretising the continuum nanowire model of Eq. (3) on a quasi-1D lattice with lattice spacing $a$. The restriction to 1D enables us to simulate braiding processes with reasonable computational efficiently, but means that subband and orbital magnetic effects are neglected Nijholt and Akhmerov (2016). We discuss the consequences of this approximation in Sec. V. After discretisation, the Hamiltonian for the left leg of the Y-junction is written

[TABLE]

with

[TABLE]

and where $n$ labels each site in the lattice.

In the right-hand legs of the Y-junction, however, the magnetic field and nanowire axis are no longer parallel. The continuum Hamiltonian for a nanowire making an angle $\theta$ to the $x$-axis may be obtained from Eq. 3 by making the substitutions

[TABLE]

which is equivalent to replacing $x\to x^{\prime}$ and rotating the spin-orbit coupling direction so that it remains perpendicular to the direction of propagation along the wire. In our setup, the magnetic field and spin-orbit coupling direction both lie in the $xy$-plane; however, see Refs. Osca et al., 2014; Rex and Sudbø, 2014 for a nanowire Hamiltonian in which the magnetic field may vary in three dimensional space.

After discretisation, sections of the nanowire which make an angle of $\theta$ with the horizontal are described by the Hamiltonian

[TABLE]

We will find that the misalignment of the field and the nanowire leg has a significant effect on the eigenstates of the system. In addition to the above transformations, we must also treat the central site of the Y-junction with care: since it connects to all three legs, we replace the diagonal discretisation offset term $2t(a)\tau_{z}$ with $3t(a)\tau_{z}$, so that the model has a well-defined continuum limit as $a\to 0$.

We build the discretised system and perform exact diagonalisation using the KWANT python package Groth et al. (2014), taking the lattice spacing $a$ small enough to avoid discrete artefacts. The parameters we choose are inspired by existing numerical works, based on InSb nanowires and aluminium superconductors, which have been shown to agree extremely well with experimental data Liu et al. (2017); Zhang et al. (2018). We also take additional input from the device geometries and materials used in the TASE process Gooth et al. (2017b). Our starting parameters (which we take as a ‘best case’ for such devices) are $m^{*}=0.015\leavevmode\nobreak\ m_{e}$, $\alpha_{R}=0.5\leavevmode\nobreak\ {\rm eV\leavevmode\nobreak\ \AA}$, $L=2\leavevmode\nobreak\ \mu m$, $\Delta_{0}=0.8\leavevmode\nobreak\ {\rm meV}$, and $V_{Z}=1.5\leavevmode\nobreak\ {\rm meV}$. This value of $V_{Z}$ corresponds to a magnetic field strength of $B\approx 1.3$ T, assuming $V_{Z}=\frac{1}{2}g_{\rm eff}\mu_{B}B$ with $g_{\rm eff}\approx 40$ Liu et al. (2017); Zhang et al. (2018). We choose a lattice spacing of $a=10^{-8}\leavevmode\nobreak\ {\rm m}$ (so that the device consists of 600 sites), which we find is small enough that the low energy band structure and dynamical properties of the system are stable and immune from discreteness effects. These best-case parameters will be used throughout the simulations in Secs. III and IV.

We note that our choices for $L$ and $\Delta_{0}$ are slightly larger than what has been assumed in previous simulations, but are values we believe are achievable using the TASE technique and using nitride-based superconductors Lutchyn et al. (2018); we study the effect of reducing these parameters in Sec. V. We change between the topological and trivial superconducting regimes by adjusting the chemical potential, choosing $\mu_{\rm top}=0$ and $\mu_{\rm triv}=6\leavevmode\nobreak\ {\rm meV}$, where the topological phase transition occurs at $\mu_{c}=1.279\leavevmode\nobreak\ {\rm meV}$ for the horizontal wire. The parameters $\theta$ and $\eta$, which define the geometry of the devices, will be varied.

As described in Sec. II.2, a braid can be performed by moving the topological regions of the wire (and consequently, the Majorana fermions at their boundaries) through the sequence shown in Fig. 1(b). Numerically, this is achieved by changing the local chemical potential $\mu(x,t)$ as a function of position and time. However, there is considerable freedom in deciding the steepness and functional form of the chemical potential at the boundary between different phases, where it varies between $\mu_{\rm top}$ and $\mu_{\rm triv}$. We use a sine-squared ramping potential based on that of Ref. Sekania et al. (2017), which takes the form

[TABLE]

where $r(x)$ is the linear ramp function

[TABLE]

In this expression, $\beta=5\leavevmode\nobreak\ (\mu{\rm m})^{-1}$ sets the steepness of the ramp (chosen to minimise the numerical Majorana wavefunction size), and $x_{c}(t)$ identifies the position of the base of the ramp as a function of time. In Ref. Sekania et al. (2017) the sine-squared ramp was found to yield better braiding results than the bare linear ramp, and it should also give a better approximation to the smooth chemical potential profiles that would be realised using side gates in an experiment. In Appendix A, we demonstrate that a ramp potential similar to this may be generated using a realistic arrangement of side gates. There, we estimate that approximately twenty gates would be required on each leg of the nanowire device to reproduce a ramping potential to within about 5%.

An illustration of two phase boundaries described by the sine-squared ramping profile is shown in Fig. 3. During the braiding process, these phase boundaries are moved throughout the device so that the two Majoranas attached to them form a braid. [The other two Majorana fermions remain at the ends of the nanowire legs throughout the process].

III.2 Low-lying Bulk States

For a successful braid, we require (at the very least) that the quasi-degenerate Majorana subspace remains separated from other excited states by an energy gap $\Delta E$. The size of this energy gap will impose limits on the speed at which we can perform the braid without causing diabatic excitations. For this reason, we first study the instantaneous low-lying energy levels near the Majorana subspace as the braid progresses. An example of this is shown in Fig. 4 for the Y-junction with parameters as given above, and with $\theta=15^{\circ}$. In this plot, the progress of the exchange operation is labelled by the parameter $0\leq s\leq 1$, which may loosely be thought of as a dimensionless time parameter. [Note, however, that the energy levels plotted are of the instantaneous Hamiltonian, and so $s$ is not associated with any physical time scale].

We see from Fig. 4 that, as hoped, the instantaneous Majorana energy levels remain well-separated from the bulk states throughout the braid. However, there appear to be two distinct types of low-lying excitation: first, there is a broad spectrum of ‘general’ bulk states, and secondly there is a single pair of low-lying states which approach the Majorana levels for $0.5\lesssim s\lesssim 0.9$.

The broad spectrum corresponds to bulk states that are also present in a system with periodic boundary conditions, and which numerically are found to be delocalised across the device. These define the bulk energy gap, which is a (complicated) function of the parameters of Hamiltonian (12). Importantly, this bulk gap decreases when $\theta$ increases, as the external magnetic field and the effective spin-orbit coupling field (which is perpendicular to ${\mathbf{B}}$ when $\theta=0^{\circ}$) start to become aligned. This effect is apparent in the band structure of a single straight nanowire, shown for a continuous periodic system in Fig. 5(a): at $\theta=0^{\circ}$ the electron and hole bands are gapped and symmetric about $k=0$, but when $\theta=20^{\circ}$, the bands become tilted and the bulk gap reduces. This is associated with the breaking of the symmetry $\sigma_{x}\hat{H}_{\theta}({\mathbf{k}})\sigma_{x}=\hat{H}_{\theta}(-{\mathbf{k}})$ of Hamiltonian (12) when $\theta\neq 0$ [see Ref. Nijholt and Akhmerov (2016) for further details of symmetries in TSC nanowire models]. We find numerically (for both the periodic system and the Y-junction device) that the bulk band gap collapses approximately linearly with $\theta$, from its bare value near $\Delta_{0}$ at $\theta=0^{\circ}$, to zero at $\theta\approx 35^{\circ}$, as shown in Fig. 5 (b). Indeed, Refs. Osca et al., 2014; Rex and Sudbø, 2014 consider this band tilting for the infinite continuum nanowire model, and find a bulk gap closure that is consistent with our numerical results for the Y-junction device. In particular, Ref. Rex and Sudbø, 2014 demonstrates that the critical angle should satisfy (in our notation)

[TABLE]

which agrees very well with our numerical results.

The second type of low-energy modes, the single pair of states which emerges from the bulk during the evolution, is specific to the three-legged geometry: numerically, these states have most of their density concentrated at the junction between the three legs. In the period where these modes are most significant (between $0.5\lesssim s\lesssim 0.9$ in Fig. 4), the two right-hand legs are in the topological regime and may be interpreted as two separate TSC nanowires forming a junction. When two such nanowires (each hosting a pair of MZMs) are brought together, the two MZMs that meet at the junction couple through the (approximate) Hamiltonian term

[TABLE]

where $\gamma_{1,2}$ are the relevant Majorana operators, $\Gamma$ is a coupling strength, and $\phi_{A}$ and $\phi_{B}$ are the effective superconducting pairing phases in each nanowire Kitaev (2001); Alicea et al. (2011). In general, this means that when two TSC nanowires are brought together, the two MZMs at the junction are gapped out, leaving a single, longer TSC nanowire with a single pair of MZMs. However, in the fine-tuned case where $\phi_{A}-\phi_{B}=\pi$ (a ‘$\pi$-junction’ Alicea et al. (2011)), the coupling term vanishes and the resulting system has four Majorana zero modes.

For the realistic nanowire model we are using, the induced superconducting pairing depends on the spin-orbit coupling direction Oreg et al. (2010); Lutchyn et al. (2010); Alicea et al. (2011), which is different in each leg. With a half-angle of $\theta$, the effective superconducting pairing phases are proportional to $e^{i\phi_{\rm UR}}\propto e^{-i\theta}$ and $e^{i\phi_{\rm LR}}\propto e^{-i\left(\pi-\theta\right)}$ (where $\phi_{\rm UR}$ and $\phi_{\rm LR}$ refer to the upper-right and lower-right legs, respectively). This means that as $\theta$ decreases, the two legs get closer and closer to forming a $\pi$-junction: the single pair of low-lying states correspond to these incipient Majorana modes that would arise exactly at $\theta=0^{\circ}$. The resulting gap is therefore greater at larger values of $\theta$. Expanding Eq. (19), we expect the energy gap caused by incipient Majorana modes to increase linearly with $\theta$ (or as $\sin\theta$ for large enough values of $\theta$).

These two sources of low-lying modes, the bulk gap closure and the incipient Majoranas, become dominant at different values of $\theta$. The optimum (largest) energy gap above the MZM subspace will be obtained for a half-angle which is greater than zero (to avoid a $\pi$-junction), but not so great that the bulk gap closure starts to become limiting. For our base parameter values, we find that this is achieved for a critical value $\theta_{c}\approx 20^{\circ}$, as demonstrated in Fig. 5(b). Note, however, that the energy gap at this angle is lower than its value for a single straight nanowire at $\theta=0^{\circ}$. In Sec. V, we study how this angle dependence changes for the tuning fork geometry and with different parameter values.

In addition to a large energy gap to excitations, we also desire the modes corresponding to Majorana fermions to be as close to zero energy as possible: in the ideal case, all four MZM states would be exactly degenerate. In a real system, however, the the Majorana zero modes are extended over some localisation length $\xi$, and modes from different regions of the device overlap in space, leading to an energy splitting. For two modes separated in space by a distance $\ell$, the splitting is proportional to $\Delta E_{\rm MZM}\propto e^{-\ell/\xi}$. In contrast to the ideal case, this energy splitting means that our braid operation will need to be performed fast enough that the different MZMs can be considered as a single subspace (i.e. the braid should be able to mix states freely within this subspace). Even then, the splitting will generally introduce different dynamical phases for different states in the subspace, which will need to be taken into account.

We can see the energy splitting between different Majorana states in the lower panel of Fig. 4: for example, there are two modes in red and blue with approximate energy levels $E\approx\pm 0.65\leavevmode\nobreak\ \mu{\rm eV}$ that remain constant for $0\lesssim s\lesssim 0.4$. These correspond to the two MZMs initially localised on the lower-right leg of the Y-junction, which remain stationary throughout this part of the braid. They have a larger energy splitting than the other MZM modes because they are closer together and because they exist in the region of the wire with a smaller bulk gap (since $\theta=15^{\circ}$ here). We also observe oscillations whenever two Majorana modes get closer together or move further apart.

Since the Majorana energy splitting will be important in interpreting the braiding results of the next section, we plot the energy of the largest Majorana mode (relative to zero) at $s=0.125$ as a function of $\theta$ in Fig. 5(c). We choose this specific gap in the low-energy subspace, as it is the largest constant gap and persists long enough to have noticeable effects. We expect the oscillatory splitting, on the other hand, to generate dynamical phases that on average cancel out.

We note that the energy splitting shown in Fig. 5(c) increases approximately exponentially with $\theta$ and does not show any particular behaviour at the critical angle of $\theta_{c}\approx 20$. This is consistent with the fact that the localisation length of the MZM is inversely proportional to the (local) bulk pairing gap Alicea (2012). In this case, we found numerically [as in Fig. 5(b)] that the bulk gap decays approximately linearly as $\Delta E_{\rm bulk}\propto\Delta_{0}-c\theta$ for some constant $c$. This implies a MZM splitting proportional to $\Delta E_{\rm MZM}\propto e^{-\kappa\left[\Delta_{0}-c\theta\right]}$ (for some approximately constant $\kappa$), which reproduces the behaviour of Fig. 5(c).

IV Braiding Simulations in Majorana Y-Junctions

IV.1 Braiding Numerics

We now see how these instantaneous energy considerations play out when we simulate a finite-time braiding operation. As discussed above, we perform a braid numerically by sweeping the chemical potential as a function of space and time, following the ramp function defined in Eq. (13) and the protocol illustrated in Fig. 1(b). We will assume that the braid takes place over a time scale $0\leq t\leq T$, where the total braiding time $T$ will be varied. In our units (where the Hamiltonian parameters are given in meV), one unit of time corresponds to $4.14$ ps.

In the ideal case, we would calculate the time-evolution operator at a given time $t$ through the relation

[TABLE]

where $\mathcal{T}$ indicates time ordering and where $H(t)$ is the time-dependent Hamiltonian that incorporates the changes to the chemical potential $\mu(x,t)$. Numerically, however, we must split the time-evolution operator into a finite number of steps ($N$) and use the relation

[TABLE]

where $\Delta t=t/N$, to obtain a discrete approximation to $U(t)$. In practice, we increase $N$ until the expression for $U(t)$ converges for a given time $t$. To calculate the factors $e^{-iH(j\Delta t)}$ efficiently, we use a Chebychev expansion which converges rapidly with just a small number of terms Tal-Ezer and Kosloff (1984). We vary the chemical potential $\mu(x,t)$ by changing the offset of the ramp function $x_{c}(t)$ linearly with time during the braiding process [see Eq. (13) and Fig. 1(b)].

To perform a braid in an interacting system, we would take the many-body ground state (consisting of a Slater determinant of all occupied single-particle states including MZMs) and act on it with the many-body unitary evolution $U^{\rm mb}_{\rm braid}(T)$, which is a braid that takes place over a time period $T$. For our noninteracting system, we will instead consider the simpler action of the single-particle unitary evolution $U_{\rm braid}(T)$ on a single-particle state that acts as a proxy for the many-body ground state, following the method of Ref. Amorim et al., 2015. In the absence of any interacting terms in the Hamiltonian, this procedure provides an equivalent probe of the non-abelian statistics of the MZMs.

Specifically, we recall that at $t=0$ there are four states close to zero energy that correspond to the four MZMs we expect to form under the starting arrangement of the chemical potential. The two states closest to zero energy (corresponding to orange and green in the lower panel of Fig. 4) are equal superpositions of the two Majorana fermions on the left leg of the Y-junction (labelled $\gamma_{1}$ and $\gamma_{2}$ in Fig. 1(b)). These states are equal superpositions of $\gamma_{1}$ and $\gamma_{2}$ due to the splitting that arises from their nonzero overlap. We write the negative/positive energy states as $\left|{1_{-}}\right\rangle$ / $\left|{1_{+}}\right\rangle$ and associate them with the (complex fermion) operators $c_{1}^{\phantom{\dagger}}=\frac{1}{2}\left(\gamma_{1}+i\gamma_{2}\right)$ and $c_{1}^{\dagger}=\frac{1}{2}\left(\gamma_{1}-i\gamma_{2}\right)$ respectively, as in the discussion in Sec. II.2. Similarly, we associate the other two states close to zero energy (blue and red in the lower panel of Fig. 4) with operators $c_{2}^{\phantom{\dagger}}=\frac{1}{2}\left(\gamma_{3}+i\gamma_{4}\right)$ and $c_{2}^{\dagger}=\frac{1}{2}\left(\gamma_{3}-i\gamma_{4}\right)$ and write the corresponding states as $\left|{2_{-}}\right\rangle$ and $\left|{2_{+}}\right\rangle$. These states are associated with Majorana fermions $\gamma_{3}$ and $\gamma_{4}$, which are initially localised on the lower-right leg of the Y-junction. Note that these states are defined only for the chemical potential conformation at $t=0$.

As discussed in Sec. II.2, the ideal braiding operation should map operators $c_{1}\to c_{1}^{\dagger}$ and $c_{2}\to-c_{2}^{\dagger}$ (up to a choice of signs). In terms of states, this would map $\left|{1_{-}}\right\rangle\to\left|{1_{+}}\right\rangle$ and $\left|{2_{-}}\right\rangle\to-\left|{2_{+}}\right\rangle$ (up to an overall phase). We can therefore test the success of the braid by starting at $t=0$ with the state

[TABLE]

and acting on this with the unitary evolution operator to obtain

[TABLE]

At the end of the evolution we compare the overlap of the state $\left|{\psi(T)}\right\rangle$ with the initial states $\left|{1_{-}}\right\rangle$ and $\left|{1_{+}}\right\rangle$. If the braid is successful, we should find

[TABLE]

indicating that all the weight from state $\left|{1_{-}}\right\rangle$ has been transferred to state $\left|{1_{+}}\right\rangle$. In general, however, the overlaps will not take these ideal values, and there may be additional nonzero overlaps between $\left|{\psi(T)}\right\rangle$ and other states (either excited states or other states from the MZM subspace).

IV.2 Braiding Results

If we perform the braid adiabatically (by numerically projecting onto the subspace of states close to zero energy at each instant of the evolution), then we obtain wavefunction overlaps that are exactly as in Eq. (24): in this ideal case, we have assumed that the braiding time $T\to\infty$ and that the Majorana states are degenerate, and the braiding operation reproduces the ideal theoretical prediction. In the more realistic case, however, where we perform the evolution unitarily, the gap to excited states and energy splitting within the Majorana subspace leads to deviations from this ideal behaviour. In Fig. 6, we plot the square overlaps \left|\big{\langle}\phi\big{|}\psi(T)\big{\rangle}\right|^{2} of the final state $\left|{\psi(T)}\right\rangle$ with a variety of initial states $\left|{\phi}\right\rangle$, for $\theta=15^{\circ}$, $20^{\circ}$ and $25^{\circ}$. We recall that from considerations of the instantaneous energy gaps, $\theta=20^{\circ}$ is close to the optimal case (with maximum energy gap to excitations), while smaller or larger angles lead to a suppression of the gap due to incipient Majorana modes or bulk gap closure, respectively.

Fig. 6 plots the square overlap with the ideal final state $\left|{1_{+}}\right\rangle$ (in solid blue), the square overlap with the starting state $\left|{1_{-}}\right\rangle$ (in dashed red), and the total square overlap with the initial Majorana subspace (in solid black), for each value of $\theta$. These overlaps are plotted as a function of the total braiding time $T$ and are the overlaps that would be measured immediately after the braiding protocol has been completed. Unlike in the adiabatic case, the probability of the braid being successful (given by \left|\big{\langle}1_{+}\big{|}\psi(T)\big{\rangle}\right|^{2}) oscillates as a function of the total braiding time $T$, and is suppressed significantly if the braid is performed too quickly. The dark grey shading indicates ranges of $T$ for which the braiding may approximately be viewed as ‘successful’, in which the overlap with ideal final state $\left|{1+}\right\rangle$ is larger than the overlap with the initial state $\left|{1-}\right\rangle$. The lighter grey shading (including the dark grey areas) indicates values of $T$ for which the process may approximately be viewed as adiabatic with respect to the excited states, for which the overlap with the total Majorana subspace is greater than 95%.

The low success rate at small $T$ is in fact a reduction of the square overlap with the entire MZM subspace, indicated by the black lines in Fig. 6. This small-$T$ feature shows that the final braided state lies outside of the initial low-energy subspace, and is therefore due to diabatic excitations to higher energy states that arise during the braid if it is performed too quickly. Indeed, the overlap with the whole MZM subspace is highest for $\theta=20^{\circ}$, which has the largest energy gap to excitations. Interestingly, the equivalent overlap for $\theta=25^{\circ}$ is higher than that for $\theta=15^{\circ}$, even though the excitation gap is smaller (see Fig. 5(b)): this perhaps suggests that the braided state is more strongly excited into incipient Majorana modes than into generic bulk states. We also note that at very small $T$, the overlap with the MZM subspace initially decreases for each value of $\theta$. This is due to the braid being so fast that most of the weight of the braided state $\left|{\psi(T)}\right\rangle$ remains in the initial state $\left|{1_{-}}\right\rangle$ (which is itself in the MZM subspace), as the system does not have enough time to adjust.

We now turn to the oscillatory behaviour of the overlaps with $\left|{1_{-}}\right\rangle$ and $\left|{1_{+}}\right\rangle$, shown in red and blue, respectively, in Fig. 6 (and emphasised by the dark grey shading). These indicate oscillations of amplitude between the different states within the (quasi)-zero-energy subspace (and includes smaller amplitude oscillations with states $\left|{2_{+}}\right\rangle$ and $\left|{2_{-}}\right\rangle$, not shown). These arise due to the fact that the four MZM states are not exactly degenerate at zero energy throughout the braid. Since the primary oscillation period decreases as $\theta$ increases, they must be associated with an energy splitting that increases as a function of $\theta$.

In Fig. 5(c) we identified the largest such (constant) energy splitting, arising for $0\leq s\leq 0.4$ in Fig. 4, which indeed increases exponentially with $\theta$. We can use this to obtain a heuristic explanation for the oscillatory behaviour. First, we recall that the braid we have simulated begins in state $\left|{1_{-}}\right\rangle$, which is an equal superposition of Majorana fermions $\gamma_{1}$ and $\gamma_{2}$. Halfway through the braid (i.e. at $s=1$ in Fig. 4), we have exchanged Majoranas $\gamma_{2}$ and $\gamma_{3}$, and so in the ideal case, the instantaneous state should be an equal superposition of the original $\gamma_{1}$ and a MZM located where $\gamma_{3}$ was initially. We can write this instantaneous state in terms of the initial states as

[TABLE]

where $\phi$ is the (unspecified) phase difference between the two pairs of terms. The combination of states $\left|{1_{-}}\right\rangle$ and $\left|{1_{+}}\right\rangle$ corresponds to $\gamma_{1}$, while the combination of states $\left|{2_{-}}\right\rangle$ and $\left|{2_{+}}\right\rangle$ corresponds to $\gamma_{3}$.

For the next section of the drive (corresponding approximately to $1\leq s\leq 1.4$), the instantaneous state is not an eigenstate, and so different terms in the wavefunction expansion will pick up different dynamical phases. Assuming the instantaneous energy levels are constant during this time (i.e. ignoring the oscillatory pieces in Fig. 4), we can write the evolved state as

[TABLE]

where have defined $E_{1}=E_{1_{+}}=-E_{1_{-}}$ and $E_{2}=E_{2_{+}}=-E_{2_{-}}$ as the energies of the corresponding eigenstates. From Fig. 4, the energy $E_{2}$ is the dominant energy scale, and will cause a shift in the relative phase of $\left|{2_{-}}\right\rangle$ and $\left|{2_{+}}\right\rangle$. This will cause oscillations between Majorana modes $\gamma_{3}$ and $\gamma_{4}$, which will in turn alter the final state of the braid.

This picture suggests that the oscillations will have a dominant characteristic time period of approximately $T^{\prime}=\pi/E_{2}$, which is comparable with the numerical oscillations observed in the braiding overlaps. Specifically, extracting the MZM splitting from the data in Fig. 5, we obtain predicted time periods of $T_{1}\approx 4800,1800,650$ for $\theta=15^{\circ},20^{\circ},25^{\circ}$, respectively. Extracting the time periods from Fig. 6 directly instead gives $T_{1}\approx 3000,1250,440$. The differences likely arise from the oscillatory regions of the MZM energy splitting (where the energy levels are changing with time), and from other energy splittings that arise during the braid. In addition to this behaviour with period $T_{1}$, other oscillations with approximate time period $T_{2}=2\pi/(E_{2}-E_{1})\approx 2T_{1}$ are also visible in the overlap with states $\left|{2_{-}}\right\rangle$ and $\left|{2_{+}}\right\rangle$ (not shown).

Overall, these numerical simulations suggest that a topological braiding operation is feasible in a realistic nanowire device, but that dynamical phases may need to be accounted for carefully. Notably, the final braided state remains within the low-energy subspace for a total braid duration of approximately $T\gtrsim 1000$ in our units. This corresponds to $T\gtrsim 4$ ns, which compares favourably to measured quasiparticle poisoning times which are typically on the order of microseconds Higginbotham et al. (2015); Albrecht et al. (2017).

The oscillatory behaviour of the final state overlaps demonstrates the difficulties that may arise from dynamical phases. To use the simulated devices directly, it would likely be necessary to sweep the braiding time $T$ to find a maximum braiding success probability. This may require repeated measurements to determine the complete profile of the overlap curves, and may be complemented by performing repeated braiding operations (as suggested in Ref. Clarke et al., 2017). More generally, these dynamical phase effects can be minimised by reducing the energy splitting between Majorana modes. In particular, the most significant splitting arises from the two Majorana modes localised in the angled prong of the Y-junction at the beginning of the braid. This can be reduced by increasing the length of the device legs (so that the Majorana modes are further apart) or by increasing the effective superconducting gap (to reduce the extent of the Majorana modes). However, it may be more experimentally favourable to make just the lower-right leg of the device longer (possibly at the expense of the others), so that the most significant energy splitting is reduced. It may also be advantageous to use a narrower prong angle $\theta$, which reduces the Majorana energy splitting (albeit at the expense of a smaller gap to excitations). Finally, the critical angle $\theta_{c}$ is known to increase as the ratio $\Delta_{0}/V_{Z}$ approaches one Rex and Sudbø (2014) (at the expense of a smaller band gap at $k=0$); there may therefore be some advantages to reducing the external magnetic field.

V Tuning Forks and Other Parameter Regimes

In the discussions above, we considered a Y-junction device with experimentally motivated ‘best case’ parameters introduced in Sec. III.1. In this section, we discuss the benefits that may arise from instead using a tuning fork geometry [as in Fig. 2(b)]. We also consider how robust our conclusions are to changes in the underlying parameters.

In Fig. 7(a), we show the minimum gap to excited energy states (during the braid) as a function of half-angle $\theta$, for a variety of choices for $\eta$, the geometric parameter that determines when the tuning fork straightens out (see Fig. 2(b)). We see that as $\eta$ decreases, and more of the device becomes aligned with the magnetic field, the larger the maximum gap becomes (and the larger the value of $\theta$ is for which this occurs). This arises because the tuning fork geometry has higher energy bulk states as compared to the Y-junction: more of the device is aligned with the magnetic field, and so the band tilting effect shown in Fig. 5(a) is less significant. In this way, the maximum excitation gap arises for larger values of $\theta$, when the effects of the bulk gap closure and incipient Majorana modes are equal. For $\eta=0.25$, the maximum gap arises at $\theta\approx 23^{\circ}$.

The tuning fork geometry also leads to smaller energy splitting between the Majorana modes, as shown (for $s=0.125$ during the braid) in Fig. 7(b). This is because the splitting is exponentially small in the Majorana localisation length, which in turn is inversely related to the bulk energy gap. For $\eta=0.25$, the energy splitting is almost two orders of magnitude smaller than for the Y-junction ($\eta=1$).

These two features suggest that tuning fork devices are more promising geometries for realising and braiding Majorana modes than Y-junctions. The larger bulk gap means that the braid should not need to be performed as slowly to be in the adiabatic regime, and the reduced Majorana energy splitting should reduce the unhelpful effects of dynamical phases. Indeed, in Fig. 7(c) we plot the final braiding overlaps for the tuning fork geometry and see that both of these advantages are apparent: The overlap with the Majorana subspace is higher for shorter braiding times $T$ and the dominant period of oscillation is significantly longer, and is barely visible in the figure. Overall, the final braided state is close to the ideal state for almost all of the simulated braiding times. This suggests that the tuning fork geometry has significant advantages over the Y-junction geometry.

Unfortunately, we cannot improve these energy scales indefinitely by taking $\eta\to 0$. Not only would this pathological limit be impossible to engineer, but it would also lead exactly to a $\pi$-junction forming at the node, which would interfere with the braiding process. The optimum value of $\eta$ is likely to depend on the spatial extent of the Majorana modes: when two Majorana modes meet at the junction, they should have most of their weight on an angled part of the device, to avoid forming a $\pi$-junction.

In passing, we note that the tuning fork geometry is likely to have an additional advantage over the Y-junction geometry in terms of orbital effects, which lie beyond the scope of the one-dimensional model. Specifically, Peierls phases are expected to be gained by electrons hopping perpendicularly to the external magnetic field, which is most significant in the angled legs of the device. According to Ref. Nijholt and Akhmerov, 2016, this is expected to reduce the gap to excitations in a similar manner to the interplay between the Zeeman term and spin-orbit coupling. This effect should be smaller for the tuning fork geometry, which has more of the device aligned with the external field. A detailed study of the orbital effect of the magnetic field in three-legged devices (using 3D simulations) remains an interesting avenue for future research.

In order to study the robustness of our conclusions, we now relax our best-case parameter choices and observe how the band gap and Majorana energy splitting vary. Specifically, we alter the induced pairing gap, the length of the device legs, the spin-orbit coupling strength, the external magnetic field, and aspects of the braiding protocol. The results are plotted in Fig. 8.

We see that in general, modifying these parameters away from their ‘best case’ values decreases the size of the bulk gap (upper panel of Fig. 8) and increases the value of the Majorana energy splitting (lower panel of Fig. 8). The most significant parameters affecting the Majorana splitting are the induced superconducting pairing gap $\Delta_{0}$, the length of the device legs $L$, and the external magnetic field strength $B$. Changes in $\Delta_{0}$ or $L$ of about 25% increase the Majorana splitting energy by about one order of magnitude. Modifying the other parameters (notably changing the steepness and strength of the chemical potential profile) has a comparatively smaller effect on the splitting.

A change in the magnetic field strength $B$ affects the Hamiltonian in two important ways: first, it directly changes the Zeeman energy through the relation $V_{Z}=\frac{1}{2}g_{\rm eff}\mu_{B}B$, and secondly, it indirectly reduces the superconducting pairing energy as the external magnetic field approaches the critical field of the superconductor, $B_{c}$. We incorporate this suppression through the relation

[TABLE]

which reproduces the pairing collapse observed in experiments well Zhang et al. (2018). For the nitride-based superconductors we have assumed in this text, the critical field strength is about $B_{c}=5$ T Pushp et al. (2018), which we have used in Fig. 8. Specifically, for the pink curve in Fig. 8, we have increased the field from $B\approx 1.3$ T to $B\approx 2.6$ T. This changes the Zeeman energy from $V_{Z}=1.5$ meV to $V_{Z}=3.0$ meV and reduces the pairing strength from $\Delta_{0}=0.8$ to $\Delta_{0}=0.68$. As a result, the system has a much larger Majorana energy splitting and a substantially reduced bulk gap. From these energetic considerations, it seems sensible to choose an external magnetic field which is large enough for the system to be in the topological regime, but not so large that the pairing gap starts to collapse significantly.

We note from Fig. 8 that the bulk gap is more robust to changes in the underlying parameters than the Majorana energy splitting, although the precise value of the optimum half-angle $\theta_{c}$ is slightly different in each case. The most dramatic changes arise when the superconducting pairing gap $\Delta_{0}$ is reduced or the magnetic field strength $B$ is increased, with a smaller bulk energy gap having a significant effect on the overall gap profile as a function of $\theta$. This would be particularly important for devices with aluminium superconductors, whose pairing gap is usually taken to be about $0.2\leavevmode\nobreak\ {\rm meV}$, and whose critical field $B_{c}$ is close to $1.1$ T Liu et al. (2017); Zhang et al. (2018), both significantly lower than for nitride superconductors.

Following the plots in Fig. 8, we would expect suitable aluminium-based devices to require a much narrower half-angle (in the range $5^{\circ}$–$10^{\circ}$). Perhaps surprisingly, however, the magnitude of the maximum gap does not decrease linearly with $\Delta_{0}$, and so aluminium-based devices should still be suitable for braiding operations. In addition, the smaller required values of $\theta$ should also suppress the Majorana splitting, to partially compensate for the increase in localisation length due to the smaller bulk gap. These features could likely also be further improved by using a tuning fork geometry.

Overall, these results suggest that our conclusions on the feasibility of braiding are reasonably robust to changes in the parameters of the underlying model. From an experimental perspective, it would seem most important to achieve a large induced superconducting gap and long device legs (although this latter condition is likely to amplify environmental effects). Variations in the model parameters affect the Majorana energy splitting more strongly than the overall gap to excitations, suggesting that the effects of dynamical phases will differ from device to device and may need to be characterised carefully in each case. However, these issues are likely to be improved if a tuning fork geometry is used instead of a pure Y-junction, following the discussion above.

VI Conclusion

In this work, we have studied the feasibility of performing braiding operations in nanowire devices in the shape of Y-junctions and tuning forks. As a stepping stone towards understanding a real braiding experiment, we have used a continuum model to describe these devices, with parameters extracted from existing transport experiments. This model has previously been shown to give good agreement with experimental measurements.

We began by computing the energetic properties of a Y-junction as a putative braid is performed, noting that the gap to excitations, which affects the speed at which a braid can be performed, depends on two factors. First, the bulk gap may collapse due to the misalignment of the external magnetic field and the effective spin-orbit coupling field. This effect is most prominent when the Y-junction makes a wide angle to the external field. However, if the angle between the legs is too narrow, then incipient Majorana modes may form at the junction, which also reduces the gap to excitations. The optimum angle which maximises this gap depends on the device parameters, but for our choice of values was close to $\theta_{c}\approx 20^{\circ}$.

A realistic device also exhibits energy splitting between low-lying Majorana modes. This is determined by the localisation length of the Majorana modes, which in turn is inversely related to the bulk gap. The energy splitting was found to increase exponentially with the Y-junction angle, in approximate correspondence with the closure of the bulk gap.

Both of these static considerations have important implications for the feasibility of a real braiding operation. Notably, a smaller gap to excitations means that the braid must be performed more slowly, while nonzero Majorana splitting generates dynamical phases that produce oscillations in the final braided state. We studied these effects in detail by simulating a braiding operation for a Y-junction with three different values of $\theta$ and calculating the overlap of the final state with the theoretical prediction. Overall, our results suggest that such devices are capable of performing nonabelian braiding operations, but that dynamical phase effects may need to be carefully taken into account. In order to be successful, the braiding operation should take place over a time scale of at least 4 ns, which compares favourably to existing estimates of quasiparticle poisoning times Higginbotham et al. (2015); Albrecht et al. (2017).

In Sec. V, we showed that a tuning fork geometry has significant advantages over the Y-junction geometry. By increasing the proportion of the nanowire that is aligned with the external magnetic field, the gap to excitations can be increased and the Majorana energy splitting reduced. In turn, this reduces the prominence of dynamical phase oscillations, and means that dynamical excitations to higher energy states are suppressed, increasing the speed with which the braid can be successfully performed. Tuning forks are also expected to be more robust to orbital field effects. Beyond this, we studied the robustness of our simulations to small variations in the parameters of our underlying model, finding that the induced pairing gap, and to a lesser extent the length of the device, are the most important factors in determining the success of a braid.

Our numerical results provide qualitative and quantitative statements about the feasibility of braiding in realistic nanowire devices, which we believe are overall encouraging. In particular, our simulations suggest that such a braid lies within experimental capabilities, although many steps undoubtedly remain before this becomes a reality. To further bridge this divide between theory and experiment, a number of extensions to this numerical work could be made. First, while we have based our calculations on a model which has previously shown good agreement with experiments, it remains a one-dimensional model that is incapable of capturing certain effects (notably the effect of subbands and the orbital effect of the field). It would be interesting to extend this type of simulation to a full three-dimensional model which includes these extra features, and perhaps more realistic gating setups. In addition, we introduced superconducting pairing to the device through a bare pairing term in the Hamiltonian. There are now a number of works which have introduced methods to carefully model the superconductor-nanowire interface (see, for example, Ref. Winkler et al., 2018), which could be incorporated into a simulation of this kind. Finally, there are many different mechanisms which affect the lifetime of a Majorana qubit (see, for example, Refs. Rainis and Loss, 2012; Aseev et al., 2018). A direct calculation of such lifetimes for a specific Y-junction geometry would allow the feasibility of a braiding operation to be more sharply evaluated. Ultimately, of course, the feasibility of braiding as a tool for quantum computation can only be confirmed by a successful experiment.

Acknowledgements.

We gratefully acknowledge discussions and ongoing collaborations with B. Madon, M. A. Mueed, F. Nichele, H. Riel and M. Ritter. We also thank S. Rex and J. Perk for a number of very useful comments on the manuscript. This work was supported by DARPA Topological Excitations in Electronics (TEE) project No. 140D6318C0028. The authors also acknowledge computational resources supported by the NSF under CAREER DMR-1455368.

Appendix A Side Gate Arrangement

In this appendix, we show that a chemical potential ramp profile similar to the one used in our simulations could be achieved using a realistic arrangement of side gates. In particular, a rough calculation suggests that a gate density of approximately one gate every 100 nm should be sufficient to realise a ramp profile to within fluctuations of approximately 5%. For the devices we consider in the main text, this corresponds to about twenty gates per leg.

To simplify the calculation, we assume that the gates are all identical rectangular plates with width $2W$, and are arranged symmetrically about the nanowire device as shown in Fig. 9a. Each pair of gates has a separation distance of $2S$, while neighbouring pairs are separated by a distance $2M$. Ignoring any screening effects, a pair of gates centred at $(x=x_{j},y=0,z=0)$ and held at a potential $U_{j}$ generates an electrostatic potential energy at $(x,y,d)$ given by Kjaergaard et al. (2016); Davies et al. (1995)

[TABLE]

with

[TABLE]

and

[TABLE]

We will calculate the total potential energy at the centre of the nanowire device, ignoring screening effects, by summing the potentials from a series of gate pairs held at different potentials $U_{j}$ and centred at different locations $x_{j}$. We will assume that the centre of the nanowire is at a height of $d=20\leavevmode\nobreak\ {\rm nm}$ (as is the case for nanowires grown using TASE), and that each pair of gates has a separation of $2S=80\leavevmode\nobreak\ {\rm nm}$, which is two nanowire widths. We will allow $W$ and $M$ to be varied. Since the overall energy scale can be altered using the finger gates and the back gate, we will try to find a gate arrangement that reproduces a ramp function between the arbitrary units of zero and one. To simplify the calculation we will use a linear ramp function (rather than the sine-squared function used in Eq. 13), but the results should be similar for any ramp function that varies over the same length scale.

We first consider oscillations in the electrostatic potential away from the ramp itself, in the region where the total potential should be constant (and, in our units, equal to one) [see Fig. 9b for an illustration of these oscillations]. In this region, all gates will be held at the same potential $U_{j}=U$, and the total potential in the centre of the nanowire can be found by summing a series of potentials of the form in Eq. (28), each with a different offset $x_{j}=2j(W+M)$. The potential will vary as a function of $x$ along the nanowire device, taking its largest value at $x=x_{j}$ (aligned with the centre of a pair of gates) and a minimum value at $x=x_{j}+W+M$ (aligned with the midpoint between two neighbouring pairs of gates). Approximating the number of gates as infinite, we can write

[TABLE]

where we have suppressed the variables $y=0$, $d=20\leavevmode\nobreak\ {\rm nm}$ and $U_{j}=U$ which are the same for all gates, and where pairs of gates are centred at $x_{j}=2j(W+M)$.

We quantify the fluctuations in the potential with the ratio

[TABLE]

which, if the oscillations were sinusoidal, would give the maximum absolute deviation relative to the mean. Fig. 9c shows the regions in the $W$-$M$ plane for which $\Delta<5\%$ and $\Delta<2.5\%$. The bottom left of the figure, corresponding to more, narrower gates, reproduces the constant potential more accurately. In the main text, each leg of the nanowire device was taken to have a length of $2\leavevmode\nobreak\ \mu{\rm m}$. With gate dimensions as in Fig. 9a, each leg can therefore support approximately $N$ gates, where

[TABLE]

We superimpose contours of constant $N$ on Fig. 9c to indicate how varying $W$ and $M$ alters the number of gates per leg. In particular, we can achieve a fluctuation ratio of $\Delta<5\%$ using $W=M=25\leavevmode\nobreak\ {\rm nm}$, or approximately 20 gates per leg.

However, we must also check that this gate arrangement can reproduce the varying part of the ideal potential profile. To do this, we take $W=M=25\leavevmode\nobreak\ {\rm nm}$ and try varying the individual gate potentials $U_{j}$ to reproduce a linear ramp function. Explicitly, we write the total potential as

[TABLE]

and try to reproduce the linear ramp function

[TABLE]

(see Eq. (17)) by minimising the (scaled) mean square error (MSE),

[TABLE]

We carry out the minimisation by taking a series of eleven gates and performing gradient descent to find a minimum of the MSE as a function of $\{U_{j}\}$. As our starting point, we take each $U_{j}$ to be equal to the value of the ideal ramp function at the corresponding $x_{j}$.

Fig. 10 shows the reproduced ramp function for a ramp starting at the midpoint of a gate (i.e. with $x_{c}=0$) and a ramp starting at the midpoint between neighbouring pairs of gates (i.e. $x_{c}=50$). Also shown is the ideal linear ramp function $\mu(x)$ in grey, along with a shaded region corresponding to $\mu(x)\pm 0.05$ (which is the fluctuation tolerance we expect in the region of constant potential). In the first case, the gate arrangement is able to reproduce the linear ramp function within the tolerance band. However, in the second case, the reproduced profile lies slightly outside the tolerance band at the top and bottom of the ramp. Despite this, in both cases, the realistic gate arrangement is able to reproduce a potential profile that varies between zero and one in a roughly linear fashion. By continuously tuning the gate potentials $U_{j}$ between these optimal values (shown in the central panel of Fig. 10), the ramp function will move smoothly across the device as required to perform a braid.

Overall, this rough calculation suggests that approximately 20 gates would be needed per leg of the nanowire device to reproduce our ideal ramp function within a reasonable tolerance. Equivalently, the gates should have a width of approximately 50 nm and be separated by a gap of approximately 50 nm. An interesting avenue for further work would be to study the effects of the unwanted oscillations on the braiding process, and to see whether the 5% tolerance used here is sufficient or even if it can be relaxed. Of course, a complete discussion of this issue, beyond the scope of this work, would also require a detailed study of the screening effects that take place within the nanowire device, as well as the variation of the chemical potential across the nanowire cross section.