Boosting Majorana zero modes

One-dimensional topological superconductors are known to host Majorana zero modes at domain walls terminating the topological phase. Their nonabelian nature allows for processing quantum information by braiding operations which are insensitive to local perturbations, making Majorana zero modes a promising platform for topological quantum computation. Motivated by the ultimate goal of executing quantum information processing on a finite timescale, we study domain walls moving at a constant velocity. We exploit an effective Lorentz invariance of the Hamiltonian to obtain an exact solution of the associated quasiparticle spectrum and wave functions for arbitrary velocities. Essential features of the solution have a natural interpretation in terms of the familiar relativistic effects of Lorentz contraction and time dilation. We find that the Majorana zero modes remain stable as long as the domain wall moves at subluminal velocities with respect to the effective speed of light of the system. However, the Majorana bound state dissolves into a continuous quasiparticle spectrum once the domain wall propagates at luminal or even superluminal velocities. This relativistic catastrophe implies that there is an upper limit for possible braiding frequencies even in a perfectly clean system with an arbitrarily large topological gap. We also exploit our exact solution to consider domain walls moving past static impurities present in the system.

One-dimensional topological superconductors are known to host Majorana zero modes at domain walls terminating the topological phase. Their nonabelian nature allows for processing quantum information by braiding operations which are insensitive to local perturbations, making Majorana zero modes a promising platform for topological quantum computation. Motivated by the ultimate goal of executing quantum information processing on a finite timescale, we study domain walls moving at a constant velocity. We exploit an effective Lorentz invariance of the Hamiltonian to obtain an exact solution of the associated quasiparticle spectrum and wave functions for arbitrary velocities. Essential features of the solution have a natural interpretation in terms of the familiar relativistic effects of Lorentz contraction and time dilation. We find that the Majorana zero modes remain stable as long as the domain wall moves at subluminal velocities with respect to the effective speed of light of the system. However, the Majorana bound state dissolves into a continuous quasiparticle spectrum once the domain wall propagates at luminal or even superluminal velocities. This relativistic catastrophe implies that there is an upper limit for possible braiding frequencies even in a perfectly clean system with an arbitrarily large topological gap. We also exploit our exact solution to consider domain walls moving past static impurities present in the system.

I. INTRODUCTION
Originally, Majorana fermions describe fermionic excitations in relativistic quantum field theories which are their own antiparticles [1]. More recently, Majorana bound states or zero-energy Majorana fermions (often also loosely referred to as Majorana fermions) have become a popular and rapidly developing research field in quantum condensed matter physics [2,3]. An important impetus for this field is provided by topological quantum information processing which, in its simplest incarnations, might be based on Majorana bound states [4,5]. While, mutatis mutandis, Majorana bound states retain the field-theory concept of a fermion which is its own antiparticle, their original relativistic nature typically plays no role. In fact, Majorana bound states are predominantly considered as static and localized excitations and their motion, if considered at all, is usually treated as adiabatic. (For notable exceptions, see Refs. 6-9.) Indeed, adiabatic motion of Majorana bound states is underlying their nonabelian braiding statistics [5,[10][11][12][13][14] which is a cornerstone for their use in topological quantum computation [5]. While information storage would rely on nonlocal qubits exploiting the 2 N -fold ground state degeneracy in the presence of 2N Majorana bound states, it is envisioned that information processing proceeds via adiabatic braiding operations of the Majorana bound states. Unlike the more familiar abelian cases (bosons, fermions, or anyons) in which the wavefunction is multiplied by a phase factor upon particle interchange, braiding of Majorana fermions causes a unitary rotation of the initial wavefunction in the space of degenerate ground states.
An essential building block of braiding operations consists of moving a Majorana bound state through the sys-tem at a constant velocity v. While adiabatic motion corresponds to the limit of small v, it may be desirable to perform quantum information processing on finite time scales and thus at finite values of the velocity v. This motivates us in this paper to consider the motion of Majorana bound states at arbitrary velocities. The centerpiece of our work is an exact solution of this problem for a particular model which exploits an emergent Lorentz invariance of the underlying equations and thus reintroduces aspects of relativity into the dynamics of Majorana bound states (albeit again in a somewhat different manner than in the original relativistic-field-theory context).
Several key features of our results can indeed be interpreted in terms of familiar effects of special relativity such as Lorentz contraction and time dilation.
Majorana bound states were initially discovered in the condensed matter context as excitations of correlated electron phases, most notably in certain fractional quantum Hall states [15]. Much of the recent excitement in the field stems from a series of proposals that predict Majorana bound states in hybrid systems made of more conventional materials [16][17][18][19][20][21]. Significant theoretical and experimental efforts have been expended on effectively one-dimensional systems with strong spin-orbit interactions proximity coupled to conventional superconductors. This can be realized following a seminal proposal of Fu and Kane [16,22] which is based on topological insulator edge states or following an alternative route which utilizes semiconductor quantum wires [19,20]. Indeed, several experiments on quantum wires may have already provided evidence for Majorana bound states [23][24][25][26]. While braiding operations are ill-defined in strictly onedimensional systems, the nonabelian statistics of Majorana bound states does survive in networks of quantum wires [27]. In this paper, we consider Majorana bound states moving along such a one-dimensional system. In onedimensional systems, Majorana bound states (as well as a discrete set of further Andreev bound states) are localized at domain walls between superconducting phases of different topology. These domain walls can be induced in the system by varying parameters -such as magnetic field, proximity-induced pairing strength, supercurrent, or chemical potential (as controlled by a keyboard of gate electrodes) -along the wire [19,20,[27][28][29]. By appropriately changing these parameters in time, one can make the domain wall move along the wire at arbitrary velocities v (see Fig. 1). We find that the nature of the solution is very different for subluminal (v < u) and superluminal (v > u) velocities of the domain wall. Here, u denotes the edge-mode velocity of the topological-insulator edge which takes the role of an effective speed of light. In particular, the Majorana bound state exists only for subluminal velocities and dissolves into a continuum of states for luminal or superluminal velocities. While our exact solution is for a particular model of a domain wall, we argue that essential features apply much more generally. Specifically, important aspects carry over not only to more general models of the domain wall but also from topological-insulator edges to quantum-wire-based systems since for parameters in the vicinity of the topological phase transition, the low-energy theory takes on a universal form.
The paper is organized as follows. In Sec. II, we review the model and the phase diagram of a topological superconductor based on a topological-insulator edge proximity coupled to an s-wave superconductor and discuss the spectrum and the wavefunctions of Andreev bound states (including the Majorana bound state) for a smooth domain wall. In Sec. III, we establish the effective Lorentz invariance of the Hamiltonian and exploit this invariance to find an exact solution for the problem of a moving domain wall. For subluminal velocities, this solution relies on Lorentz boosting the system to a renormalized static problem, while for superluminal velocities, the solution proceeds via a reference frame in which the problem becomes translationally invariant. We also show that the exact solution has an interesting counterpart in graphene in crossed electric and magnetic fields. Sec. IV is devoted to a discussion of various physical consequences of our exact solution. We discuss the solution for both, subluminal and superluminal velocities as well as the relativistic collapse of the spectrum when reaching the effective speed of light. We also consider the effect of static impuri-ties which, in the comoving frame, act as time-dependent perturbations to the Andreev bound states. In Sec. V, we discuss various generalizations of our results, including the generalization from topological-insulator edges to systems based on semiconductor quantum wires. Sec. VI collects our conclusions.

A. Topological-insulator edge
We begin by reviewing the Hamiltonian of a topological-insulator edge proximity coupled to a superconductor and subject to an applied magnetic field [22]. Our starting point is a pairing Hamiltonian in secondquantized form,Ĥ Here, we indicate operators acting on the many-body Hilbert space with a hat. The resulting Bogoliubov-de Gennes (BdG) Hamiltonian H(x) takes the form where σ i (τ i ) are Pauli matrices acting on spin (particlehole) space, while u denotes the edge-mode velocity and µ, ∆, and B are the chemical potential, the superconducting pairing strength, and the magnetic field, respectively. Throughout this paper we set = 1.
A competition between ∆ and B leads to two topologically distinct phases, depending on whether the gap is dominated by B or ∆. For spatially constant µ, B and ∆, the two topological phases can be described by the sign of the gap function ∆ top = ∆ 2 + µ 2 − |B| such that ∆ top < 0 (∆ top > 0) for the B-(∆-)dominated phase. A transition between these phases necessarily requires a closing of the gap ∆ top . In inhomogeneous systems, it is possible to realize regions of distinct topology. The closing of the gap at the domain wall between these regions leads to localized Majorana bound states.

B. Solution for a static domain wall
Choosing µ = 0 allows for a straightforward decoupling of the four-component Bogoliubov-de Gennes equation [20] (see Sec. V B for a more general discussion). By observing that H commutes with σ x τ x , the BdG equation can be separated into two two-component subspaces described by the Hamiltonians where the (∓) spaces are spanned by the σ x τ x eigenstates {(∓1, 0, 0, 1), (0, ∓1, 1, 0)}, respectively. (We denote the Pauli matrices in the resulting subspaces also by σ i .) For definiteness, we assume that B, ∆ > 0 in the vicinity of the domain wall which we take to be located at x = 0. By comparing with Eq. (3), we see that with these choices, only the (−) subspace has a vanishing gap, while the (+) subspace remains essentially unaffected by the domain wall, retaining a finite gap of B(x) + ∆(x) ∼ 2B(0) in its vicinity. Then, the relevant gap function which vanishes at x = 0 takes the form ∆ top (x) = ∆(x) − B(x). We can therefore focus on the low-energy Hamiltonian which has the form of a two-component Dirac Hamiltonian whose mass term changes sign at x = 0. It is well known [30] Eq. (5) then takes the form of a harmonic oscillator Hamiltonian. The constant shift of ubσ y can cancel the zero-point energy and one thus finds the energy spectrum of H − as with a single zero-energy (Majorana) eigenstate n = 0.
Here, n is a (possibly negative) integer, ω = √ 2ub, and we assumed u, b > 0. The finite-energy eigenstates are given by where the g n (x) denote harmonic oscillator eigenfunctions with oscillator length ξ = u/b. The Majorana wavefunction is localized at the domain wall, has a Gaussian form, and can be read off from Eq. (7) by setting n = 0, g −1 = 0, and including an additional factor of √ 2 for correct normalization. The result takes the form The There, the (±) subspaces switch roles relative to x = 0 such that the (+) subspace describesγ 0 , while the (−) subspace remains gapped.

A. Lorentz invariance
We now consider a domain wall moving along the topological-insulator edge at a constant velocity v. This can be described by the Hamiltonian (4) with a timedependent gap function To solve this time-dependent problem, we search for a transformation into an appropriate moving frame. The linear dispersion of the BdG Hamiltonian (4) suggests that this can be achieved by a Lorentz transformation.
To explicitly establish the behavior of the timedependent BdG equation under Lorentz transformations, it is useful to multiply by σ x . This brings Eq. (10) into a "covariant" form which treats space and time coordinates on an equal footing, Here, the matrices α 0 = σ x and α 1 = −iσ y indeed fulfill the Dirac algebra α i , α j = 2g ij with the metric tensor g ij = diag(1, −1). Eq. (11) therefore takes the form of a two-component Dirac equation with an effective speed of light u and a space-time dependent mass term −∆ top (x− vt).
We first set the mass constant and review the known invariance [31] of the Dirac equation under a Lorentz boost with the usual relativistic notations β Λ = v Λ /u and γ Λ = 1/ 1 − β 2 Λ for a boost velocity |v Λ | < u. In the four-vector notation x µ = {ut, x} and ∂ µ = {∂ ut , ∂ x }, the Lorentz boost in Eqs. (12) and (13) implies that the derivatives transform as ∂ µ = Λ ν µ ∂ ν , where Λ ν µ = ∂ µ x ν . To make the Dirac equation (11) invariant under the Lorentz boost Λ, we must also transform the spinor as ψ = Sψ. Multiplying Eq. (11) by S leads to where we used that the mass term is proportional to the unit matrix. The boosted Dirac equation (14) is identical to the original one for Sα µ S −1 Λ ν µ = α ν . This condition determines the transformation S of the Dirac spinor. One can motivate its form by the following analogy. For a Dirac particle at rest, the Pauli spinor transforms with S R = exp (σ x σ y θ/2) under rotations in the x − y plane. A Lorentz boost, on the other hand, can be viewed as a hyperbolic rotation in Minkowski spacetime (i.e., the x − t plane) by an angle θ Λ = artanh (β Λ ). It is therefore natural to assume that S = exp α 1 α 0 θ Λ /2 which can indeed be checked explicitly. In the following, we will be primarily interested in where we assumed β Λ > 0 for the second equality. A crucial difference between Eq. (11) and the entirely Lorentz-invariant Dirac equation is the space-time dependence of the mass term. Even though the mass is no longer constant, it is still proportional to the unit matrix such that the above arguments remain untouched. The Dirac equation in the boosted frame can then be obtained by applying the Lorentz boost Λ to the argument x − vt of the mass term, yielding the boosted Hamiltonian We can now significantly simplify the problem by appropriate choices of β Λ . We note in passing that Tsvelik employed Lorentz invariance in a related manner to discuss Majorana fermions interacting with fast bosonic fields [7].

B. Subluminal motion (β < 1)
When the domain wall moves at a subluminal velocity v < u, we can choose a boosted frame which is moving at the same velocity, v Λ = v. In this comoving frame, the boosted Hamiltonian (16) becomes time independent and takes the form of Eq. (4) with a renormalized slope b = b/γ of the topological gap. The renormalization of b can be understood as a consequence of the familiar length contraction of special relativity. In the comoving frame, the domain wall is at rest and one measures its proper length. Consequently, the size of the domain wall in the lab frame is contracted by a factor of 1/γ. Since, however, the form of ∆ top is defined by the lab-frame Hamiltonian, the size of the domain wall increases in the comoving frame, thus reducing its slope b by a factor of 1/γ.
With the knowledge of b , the solutions of the BdG equation in the comoving frame can be read off from Eqs. (6) and (7). The wavefunctions take the form with the corresponding energy spectrum The factors of 1/ √ γ in these equations follow from the b -dependence of the oscillator length ξ and the frequency ω.
To complete the solution of the original timedependent Hamiltonian (11), we still need to transform from the comoving frame back into the lab frame. This is achieved through the inverse Lorentz boost We thus obtain the lab-frame solution for the wavefunctions, where we defined n βγ 2 /u. The renormalized spectrum in the lab frame is given by This provides an exact solution for the problem of a domain wall moving at an arbitrary subluminal velocity v < u.
The renormalization of the energy spectrum in Eq. (22) can also be understood by analogy with special relativity. We already argued that the length contraction leads to a factor of γ −1/2 in the energy spectrum of the comoving frame [see Eq. (18)]. The additional factor of γ −1 originates from the time dilation (i.e., suppression of frequencies) when transforming from the comoving frame back into the lab frame.
The wavefunctions (20) can now be employed to define (time-independent) Bogoliubov operators (cp. App. A) in the lab frame. Note thatγ n ] † . This relation follows by explicit calculation and reflects that ψ where CT denotes the product of charge conjugation C and time reversal T (see App. A). Specifically, this implies that even for the moving domain wall, the n = 0 solution defines a Majorana operatorγ It is interesting to express the original many-body HamiltonianĤ (projected to the (−) subspace) in terms of the operatorsγ n . Expanding the Nambu field opera-torsΨ(x, t) in theγ Here and in the following, we use the explicit time de-pendenceĤ − (t) to denote many-body operators in the Heisenberg picture. In view of the structure of the explicit wavefunctions ψ ±n (x, t) in Eq. (20), this can be written asĤ which involves the total momentum operator We will now investigate the case of a domain wall moving at a superluminal velocity v > u. Although we can no longer boost to a comoving reference frame, superluminal motion is a perfectly physical scenario in the present context. From Eq. (16), it is clear that we cannot find a transformation to a time-independent ∆ top with β Λ < 1. In the language of special relativity, the space-time dependence of ∆ top switched from space-like to time-like. This, however, allows for a different simplification by choosing β Λ = 1/β. In this reference frame, the Hamiltonian is spatially independent and the topological gap becomes purely time dependent, taking the form ∆ top (−ut /γ) withγ = 1/ β 2 − 1. The problem becomes translationally invariant and thus, momentum is a good quantum number. We therefore make the planewave ansatz We emphasize that in contrast to the case of subluminal motion, these solutions are no longer localized. In fact, this also leads to extended states ψ k (x, t) in the lab frame which are labeled by the continuous set of momentum quantum numbers k and thus constitute a continuous spectrum. For a linear ∆ top and fixed momentum k, the timedependent partφ k (t ) of the wavefunction in the boosted frame is described by the Hamiltonian which equals that of a (rotated) Landau-Zener problem. The solutionsφ k (t ) are known to be composed of complex-argument parabolic cylinder functions [32]. To develop intuition for these solutions and to understand the connection to the subluminal case, it is helpful to multiply Eq. (28) by σ z . Up to an exchange of time and space variables, the resulting time-dependent BdG equation has the same form as for the static problem Eq. (4). The crucial difference is that the mass term is now imaginary.
As a consequence, squaring the Hamiltonian results in a deconfining harmonic potential (see Fig. 2) and one obtains temporally extended states as expected for unitary time evolution. Specifically, we find that where the D ν (z) denote parabolic cylinder functions and ω = ω/ √γ . Moreover, e ± = {1, ±1} are the eigenvectors of σ x and α ± are free constants to accommodate boundary conditions. The parabolic cylinder functions can be expressed as D ν (z) = 2 −ν/2 H ν (z) exp −z 2 /4 in terms of Hermite functions H ν (z) and show oscillatory behavior ∼ exp[∓i(ωt ) 2 /4γ] as a result of the deconfining harmonic potential (Fig. 2).
The solution ψ k (x, t) in the lab frame is again obtained by an inverse Lorentz transformation, Of particular interest is the k = 0 case, which follows by setting H 0 (z) = 1 or by directly integrating Eq. (29). The corresponding solution reads Eqs. (20) and (8)]. By lifting the restriction β < 1, we can analytically continue the subluminal solution to the superluminal case via γ → −iγ and S −1 → diag(1, i)S −1 . Note that this relation between the suband superluminal case is reminiscent of the physics of Cherenkov radiation (see, e.g., Ref. [33]) that is emitted by electrons moving faster than the speed of lightc of a dielectric medium. There, finite-frequency electromagnetic fields originating from a moving electron change their character from evanescent to propagating waves because the relativistic factor 1 − (v/c) 2 , in the same sense as here, turns imaginary when v >c.
Applying the above analytic continuation, the subluminal Majorana solution continues exactly to one of the k = 0 superluminal solutions in Eq. (32), namely the (−)-term. However, in contrast to the subluminal case, Eq. (32) contains a second independent solution. The latter would emerge from an analytic continuation of the unphysical solution ∝ exp +γ (x − vt) 2 /2ξ 2 . In addition to the continuous nature of the spectrum and the absence of bound states, this is another manifestation of the fact that the topological character is lost in the superluminal case. In fact, this is also confirmed by looking at the Bogoliubov operators associated with the (±)-terms of Eq. (32), These are no longer Majorana operators but instead fulfill γ † 0+ =γ 0− , suggesting a connection to ordinary Bogoliubov quasiparticles with opposite energies. This can be made explicit for systems where the linear dependence of ∆ top saturates at some energy scale ∆ ∞ . In this case, these solutions do indeed become plane waves asymptotically far from the domain wall with a finite energy ±∆ ∞ .

D. Mapping to graphene
Eqs. (6) and (7) for the spectrum and wavefunctions of the static domain wall bear a strong resemblance to Landau levels in graphene [34]. This is not accidental as it is indeed possible to map the static domain-wall Hamiltonian H − (x) to that of graphene in a magnetic field. In the vicinity of the K point, the latter takes the form where π i = p i + eA i denotes the kinetic momentum in terms of the vector potential A. The mapping to this Hamiltonian exploits the fact that up to rescaling, the components of the kinetic momentum operator are canonically conjugate variables, [π x , π y ] = −i/l 2 B . We start to demonstrate this mapping explicitly by rewriting the domain-wall Hamiltonian of Eq. (4) with a static ∆ top (x) = bx as Performing a spin rotation about the x-axis to rotate σ z into σ y and renaming x → y then brings Eq. (35) into the same form as the graphene Hamiltonian of Eq. (34), The static domain-wall problem therefore maps into the p x = 0 solutions of graphene subject to the vector potential eA = {−y/ξ 2 , 0, 0} and Fermi velocity v F = u. This Landau-gauge vector potential describes a constant magnetic field eB G = 1/l 2 B pointing in z-direction such that the oscillator length ξ = u/b of the localized domainwall states maps into the magnetic length l B for the graphene Hamiltonian. To obtain the wavefunctions of the original static domain-wall problem from the corresponding solution for graphene [34], one has to undo the spin rotation. This is the origin of the factor (1−iσ x )/ √ 2 in Eq. (7).
Interestingly, also the moving domain-wall problem has an analog in graphene. Repeating the same steps as above, the time-dependent domain-wall Hamiltonian can be mapped to graphene subject to a time-dependent vector potential eA = − (y − vt) /ξ 2ê x . The time dependence of A implies that in addition to the magnetic field, there is an in-plane electric field E = −∂ t A = −vB Gêx pointing in the x-direction. The magnitude of the electric field reflects by the velocity v of the domain-wall motion.
The moving domain wall therefore maps to graphene in crossed electric and magnetic fields. The latter problem has been solved in Ref. [35] and the results map exactly to our solutions in Eqs. (20) and (22). It is interesting to observe that also the transition from the discrete subluminal to the continuous superluminal spectrum has an analog in graphene. There, the subluminal case corresponds to a magnetic-field-dominated regime with a Landau level spectrum. On the other hand, the superluminal case corresponds to dominating electric fields and metallic transport. In geometric terms, the momentum-space trajectories are closed and ellipsoidal for v F |B| > |E| and turn into open hyperbolas once |E| > v F |B| [36] (see Fig. 3).

A. Renormalization effects
In Sec. III B, we found that the most immediate effect of subluminal domain-wall motion is a renormalization of the bound states, changing both the spectrum and the spatial extent of the quasiparticle wavefunctions. For a fixed quantum number n, the extent of the lab-frame wavefunctions reduces by a relativistic factor of 1/ √ γ relative to the static domain wall, cf. Eq. (21). It is interesting to note that this contrasts with the wavefunctions in the comoving frame which actually become more extended by a factor of √ γ, cf. Eq. (17). This dichotomy is explained by the fact that this increase is overcompensated by the Lorentz contraction by a factor γ when transforming back into the lab frame. As a visualization of these renormalization effects, Fig. 4 shows the form of the Majorana wavefunction at different velocities, in the comoving and lab frame, respectively. Fig. 4 also emphasizes the effect of the spinor rotation S −1 when transforming back from the comoving to the lab frame. With the first (second) component of the low energy subspace corresponding to the operatorsψ † ↑(↓) andψ ↑(↓) , this rotation shifts the weight between spin-up and spin-down particles contributing to the Majorana mode. Since S −1 depends linearly on β [cf. Eq (15)], this shift is the leading consequence of the finite domain wall movement for small velocities.
To estimate the extent of the wavefunction in the lab frame, we note that a state with quantum number n consists of harmonic oscillator wavefunctions with quantum numbers |n| and |n|−1. For a static domain wall, the size of the state is thus of order x 0 = 2|n| + 1ξ, where ξ denotes the oscillator length. From Eq. (21), we then find for the moving domain wall that the extent of the wavefunctions is given by Interestingly, this implies that the spatial extent of the Majorana bound state is given by ξ/ √ γ which tends towards perfect localization as the domain-wall velocity approaches the effective speed of light, β → 1.
While the shrinking of the Majorana wavefunction might be viewed as beneficial for the protection of a topological qubit, a finite domain-wall velocity also entails disadvantages which arise from the renormalization of the spectrum. To start with, the level spacings shrink with increasing domain-wall velocity, effectively squeezing the bound-state spectrum (see Fig. 5) and making the system more susceptible to perturbations. In addition, we should not only compare the spatial extent of wavefunctions with the same quantum numbers, but also the extent of states with the same energies. In fact, when considering states at a fixed energy E, the squeezing of the spectrum implies that we should compare a state with quantum number n for the static domain wall to a state with quantum number γ 3 n of the moving domain wall. When doing so, the spatial extent of the states is larger for the moving domain wall by a factor γ.

B. Tunneling spectroscopy of the moving bound states
To further clarify the physical significance of the renormalized spectrum, we note that the E (v) n can, in principle, be probed by tunneling spectroscopy. Unlike for previous proposals for probing Majorana bound states in tunneling experiments [37][38][39][40], we consider a tunneling contact that extends over a long distance L along the topological-insulator edge. This can, for instance, be realized by tunneling from a parallel wire. The motivation for this choice is that the bound states of the moving domain wall can be resolved in energy only when they are probed over a sufficiently long time interval.
Specifically, we consider a tunneling source described by the Hamiltonian where ε 0 accounts for the offset of the band bottom of the source wire relative to the chemical potential µ of the topological insulator edge which we choose to be µ = 0 for simplicity. Tunneling between the source and the domain-wall bound states is described by the tunneling Hamiltonian with η 0 measuring the tunneling strength. As we are interested in the low-lying bound states, we project H T onto the low energy (−) subspace by using which follows directly from the definition of the (∓) subspaces. Here, ψ (in terms of ε k = k 2 /2m−ε 0 ), the tunneling Hamiltonian takes the form Here, we used that E (v) n is measured relative to the chemical potential µ and defined ξ k = ε k − µ S as well as eV = µ S − µ in terms of the chemical potential µ S of the source. Moreover, are time-independent tunneling matrix elements. The tunneling Hamiltonian in Eq. (41) is effectively that of a more conventional static tunneling problem. If we assume that the initial state at t = 0 obeys [γ (v) n ] † γ (v) n ∝ δ nn , the tunneling current can be obtained by a standard calculation. As a result, the energy spectrum E (v) n manifests itself as (zero-temperature) peaks in the differential tunneling conductance, where we introduced the source Fermi momentum k F = 2m(eV + ε 0 ). Essential features of this result can be understood by physical considerations. First, the strength of the tunneling peaks directly reflects the fact that for the model under consideration, tunneling is momentum conserving. Thus, only the Fourier components ±k F of the bound state wavefunctions determine the strengths of the tunneling peaks. Second, the energy shifts ±vk F can be thought of as Doppler shifts associated with the relative motion of domain wall and source wire.

C. Stability of Majorana bound states
Subluminal motion at a constant velocity can be mapped to a static process and hence does not create excitations or destroy the stability of the Majorana bound state. In particular, consider a system that at t = 0 starts in a Fock state (e.g., the ground state) |Φ with respect to the quasiparticlesγ (v) n as determined by a quasiparticle distribution f n defined through [γ n |Φ = f n |Φ . This quasiparticle distribution will then stay unchanged at any later time t with respect to the quasiparticleŝ γ n,0 translated along the system by a distance vt. The absence of a motion-induced change of the quasiparticle distribution can be checked explicitly by observing that n t) which follows directly from Eq. (23).
Even in the absence of acceleration of the domain wall, the Majorana bound states (and hence the associated topological qubits) become unstable when the domain wall moves at superluminal velocities. As shown in Sec. III C, this situation is necessarily described by a time-dependent Hamiltonian and lacks the notion of a localized Majorana mode. This provides a "speed limit" for Majorana bound states and thus imposes an upper bound f max on the braiding frequency f b . If the braiding operation requires the domain walls to move along a pathlength l b , the maximal braiding frequency is The length of the braiding path must at least exceed the typical size of the Majorana bound states, l b ξ, to ensure spatially separated Majorana bound states. The highest possible f max is therefore reached in the limit l b ∼ ξ and of the order of ω, which controls the energy of the first excited bound state level. This coincides with a naive estimate that the braiding frequency should be smaller than the minigap. We emphasize, however, that conventionally, this emerges from an argument about acceleration-induced excitations which is distinctly different from the origin of Eq. (45). Moreover, the length of the braiding paths may in general differ in magnitude from the spatial extent of the Majorana bound states due to other design requirements or the need to braid Majorana bound states which are not nearest neighbors. Then, the condition in Eq. (45) is the more stringent one.

D. Impurities
So far, we considered a moving domain wall in a clean system. Experimental systems will also contain localized impurities which do not move along with the domain wall. Viewed in the comoving frame of reference, these impurities act as time-dependent perturbations, as illustrated schematically in Fig. 6. Thus, these impurities may cause transitions between domain-wall bound states which destroy Majorana-based topological qubits over time. We now apply our exact solution of the moving domain wall for the case of subluminal motion to discuss the transition probabilities for such impurity-induced excitations.
Specifically, we consider a general short-range impurity due to local variations of the chemical potential, the magnetic field, or the proximity-induced superconducting pairing. Within the low-energy subspace, all these impurity types locally modify the magnitude of the topological gap and thus take the same form (Strictly speaking, the term involving δµ anticipates our discussion of the finite-µ case in Sec. V below.) We assume that the spatial extent of the impurity is small compared to the size of the domain-wall bound states. Then, we can approximate the impurity Hamiltonian as where ν measures the impurity strength. Note that this Hamiltonian is written in the lab frame and that its matrix structure would be modified when transformed into the comoving reference frame. We now consider a sufficiently weak impurity for which we can compute the transition amplitude T f i = f |U (∞, −∞) |i in Born approximation. Here, U denotes the single-particle (BdG) time-evolution operator. The transition amplitude can be expanded in the Born series where T (n) f i describes the n-th order in the Born approximation. Here, we focus on the first-order term in terms of the transition matrix elements Note that the transition amplitudes T (1) f i are Lorentz invariant due to the Lorentz invariance of the volume element dtdx of spacetime. Thus, the transition amplitudes are independent of the reference frame and we could equally well consider the transition amplitudes in the comoving reference frame where T (1) f i describes transitions between stationary eigenstates.
Evaluating the transition matrix element in the lab frame yields where we used the exact wavefunctions (20) and defined ω f i = E f − E i as well as s n = sign(n). We set s 0 = 0 to capture the absence of a g −1 term in the Majorana solution. (We leave a factor of √ 2 implicit which would be needed for correct normalization in the n = 0 case.) Reflecting the finite interaction time between domain wall and impurity, the transition matrix element V f i (t) is appreciable only during a finite time interval and decays as a Gaussian exp(−γv 2 t 2 /ξ 2 ) for large times, i.e., when the domain wall is far from the impurity. The time integral in Eq. (49) can be performed analytically and is controlled by Franck-Condon-like matrix elements where M = max(n, n ), m = min(n, n ), and L M −m m denotes the associated Laguerre polynomials. Since the harmonic-oscillator functions are eigenfunction of the Fourier transform, M nn (q) is essentially an overlap of harmonic oscillator functions which are spatially shifted relative to each other by √ 2qξ. In terms of these overlaps, the first-order Born approximation is given by (54) The transition amplitudes |T (1) f i | are plotted for various domain-wall velocities in Figs. 7 and 8. The principal characteristics are a strong suppression of transitions between positive and negative energy states as well as a suppression for large energy differences between the initial and final states. These observations can be understood based on the correspondence between the M f i (q) and the overlaps of harmonic-oscillator functions shifted by √ 2qξ. This shift creates an overall distance between the wavefunction centers given by ∆x f i = √ 2 s f |f | − s i |i| ξ/β. On the other hand, the sum of the characteristic spatial extents of the wavefunctions is r f i = √ 2 |f | + 1/2 + |i| + 1/2 ξ. Thus, we expect a Gaussian suppression of the overlaps M f i (q) once |∆x f i | > r f i . Since β < 1, this condition is always fulfilled for s f s i = −1, i.e., for transitions between positive and negative energy states. This explains the suppression of these transitions as seen in Fig. 7. For s f s i = 1, the condition ∆x f i = r f i corresponds to f ≈ i(1 ± β) 2 /(1 ∓ β) 2 which defines the limiting straight lines in Fig. 7 (shown as dashed lines in the figure).
We are particularly interested in the stability of the Majorana bound state, as described by the transition amplitudes with either i or f equal to zero. Note that these amplitudes also describe the dominant excitations from the stationary state in which all negative energy states are occupied and all positive energy states are empty. This state most closely resembles the ground state of the time-independent system and transitions within the sets of positive or negative energy states are forbidden by the  (v) f (instead of the quantum numbers i and f as in Fig. 7). This emphasizes the √ n dependence of the energies E (v) n as well as the squeezing of the energy spectrum for β → 1. For β = 0.9, the squeezing implies that there are γ 3 ≈ 12 times as many states in a fixed energy range compared to the static case.
As the difference ∆x f 0 − r f 0 increases with increasing f , the transition amplitudes fall off exponentially with final-state quantum number f . The leading amplitude is thus given by which shows a strong Gaussian suppression for small velocities, β 1. In fact, in this limit, this suppression is a generic feature of all transition amplitudes since they involve an exponential of the form We conclude that for slow domain-wall velocities, impurities are very inefficient in creating excitations and disturbing Majorana-based topological qubits. Interestingly, the first-order Born approximation in Eq. (56) remains well behaved in the adiabatic limit β → 0, where one would expect the Born approximation to break down. By a saddle point approximation valid for β 1, one can show that this property persists to higher orders of the Born approximation. Indeed, these higher order amplitudes are of the order of

V. GENERALIZATIONS
So far we focused on the Hamiltonian (2) for a topological-insulator edge with a domain wall described by a linearly varying gap function ∆ top and zero chemical potential. As we show above, this system is particularly attractive since it allows for an exact analytical solution. At the same time, it is important to understand to which degree this solution also describes the physics of moving domain walls more generally. We will show in the subsequent sections that indeed, essential features of our exact solution carry over to much more general situations.

A. General domain-wall structures
While a domain wall with a linearly varying gap allows for an exact solution, it is more realistic to consider domain walls for which the gap varies linearly in the vicinity of the domain wall, but eventually saturates to a constant value ∆ ∞ far from the domain wall. In fact, both Figs. 1 and 6 sketch domain walls with such a structure.
There are characteristic consequences of the saturation of the topological gap at some distance away from the domain wall, even in the static case. The saturation of the topological gap implies that the domain wall binds only a finite number of discrete states and that the spectrum becomes continuous at energies larger than the maximum topological gap ∆ ∞ . Moreover, the boundstate wavefunctions exhibit an exponential decay at large distances from the domain wall and the Gaussian decay of the linear-domain-wall model is limited to intermediate distances. At the same time, unless the domain wall is abrupt on the scale of the oscillator length (as defined by the slope of the topological gap at the domain-wall position), the low-energy spectrum of the static domain wall remains accurately described by the linear-domainwall model.
The Lorentz transformation of Eq. (16) can be applied to any domain wall which moves at a uniform velocity and is thus described by a topological gap of the form ∆ top (x − vt). The mappings to a static problem for β < 1 and to a spatially homogeneous and timedependent Hamiltonian for β > 1 are possible for any domain-wall structure ∆ top (x − vt). This implies that the transition from a low-energy discrete spectrum with a Majorana zero mode for subluminal domain-wall motion to a continuous spectrum for superluminal motion is a generic feature, independent of the domain-wall structure. Thus, the maximal braiding frequency f max = u/l b in Eq. (45) carries over to this more general case. Moreover, due to the contraction of the bound-state wavefunctions of a moving domain wall, the number of bound states increases with the velocity of the domain wall and diverges when the domain-wall velocity approaches the effective speed of light u.
The asymptotically exponential decay of the boundstate wavefunctions far from the saturating domain wall also has interesting consequences for the impurityinduced transition rates. Consider a domain wall with a linearly varying topological gap near its center -with slope b and corresponding frequency ω = √ 2ub -and a saturated gap of ∆ ∞ at large distances. Then, the Gaussian time dependence V f i (t) ∼ exp − √ γvt/ξ 2 of the transition matrix element becomes modified into a simple exponential exp −∆ ∞ √ γvt/u for sufficiently large times √ γvt > ∆ ∞ /b. In the first-order Born approximation, the transition amplitude is essentially the Fourier transform of V f i (t) and thus, we find [42] T (1) when the linear gap saturates abruptly at ∆ ∞ . In this case, the exponential suppression in 1/β crosses over to an algebraic one, albeit with a prefactor which is exponentially small in the number of subgap states (∆ ∞ /ω) 2 . This crossover occurs at a new characteristic velocity β ∞ = ω f i /∆ ∞ . Interestingly, this implies that an increasing number of subgap states has a beneficial effect on the protection of braiding operations against disorder. Note, however, that the behavior for β < β ∞ is nonuniversal as it depends on how ∆ top (x) saturates. When one replaces the abrupt cutoff by a Fermi function profile of the topological gap, an exponential suppression in 1/β even persists for β < β ∞ (see Fig. 9).

B. Nonzero chemical potential
So far, our considerations focused on the particularly simple case of zero chemical potential, µ = 0. For this choice, the Hamiltonian (2) decouples exactly into two subspaces [see Eq. (3)]. In the vicinity of the domain wall, the low-energy physics is determined entirely by one of these subspaces while the other one remains gapped. This decoupling is no longer exact when µ is nonzero. However, for a saturating domain wall with ∆ ∞ 2B, we can still project the problem in a controlled manner to a low-energy subspace by expanding the Hamiltonian about the critical point where the topological gap vanishes. To do so, we start by diagonalizing the original Hamiltonian (2) for p = 0 and ∆ top = ∆ 2 + µ 2 −B = 0. Since this puts the system on the phase boundary, we find a low-energy subspace with This low-energy Hamiltonian is of the same form as for µ = 0 but with renormalized parameters which implies that our considerations carry over essentially unchanged to nonzero chemical potentials. Specifically, the topological gap takes the form and there is a downward renormalization of the mode velocityũ This implies that the effective speed of light of the system is reduced, and with it the critical velocity at which the domain-wall spectrum collapses into a continuum.
C. Quantum wires with spin-orbit coupling Spin-orbit coupling can also induce topological superconductivity in semiconductor quantum wires proximitycoupled to s-wave superconductors [19,20]. The corresponding BdG Hamiltonian takes the form Our results essentially carry over to this system as well and the reasoning follows the arguments of Sec. V B for nonzero chemical potential in the topological-insulator model. Indeed, the gap closing at the topological phase transition line of Eq. (61) always occurs at p = 0. For this reason, we can expand Eq. (61) around the critical point by the same sequence of steps as in Sec. V B. We first identify the low-energy subspace from considering p = 0 and ∆ top = ∆ 2 + µ 2 − B = 0. Thus, the low-and high-energy subspaces are again spanned by the same basis vectors. If the domain wall saturates at a sufficiently small ∆ ∞ , we can drop the terms originating from the quadratic dispersion relative to the linear spin-orbit terms. This yields the same low-energy model as in Eq. (58). Quantitatively, a straight-forward estimate yields the condition ∆ ∞ ε SO ∆ 2 /µB, where ε SO = mu 2 measures the strength of the spin-orbit coupling in the wire. Note that this condition is always fulfilled for zero chemical potential.
Finally, we note that nonlinearities in the spectrum can lead to high energy and large momentum states with energies ε k < uk. For quantum wires in the limit ε SO B, this condition applies to the states near k F ≈ 2mu where there is a gap of order ∆. These states introduce a new velocity scale v h = ε k /k and are always irrelevant for the low-energy physics for velocities v < v h . For v > v h , similar to the Landau criterion [41], the moving domain wall can in principle mix these large momentum states with the low energy states at k ≈ 0 that contribute to the bound states. For smooth domain walls, this mixing will however be exponentially suppressed.

VI. CONCLUSIONS
Every realization of topological quantum computation on a finite time scale has to deal with the dynamic evolution of the corresponding nonabelian quasiparticles. Here, we discussed a crucial ingredient of this dynamics, a moving Majorana-carrying domain wall. Understanding the effects of motion on the Majorana bound state as well as on other components of the spectrum is an essential prerequisite for determining the rate and the space needed to perform Majorana manipulations with high fidelity.
Alongside these practical motivations of our work, we uncover a set of intriguing connections between the physics of a Majorana bound state moving at a constant velocity, and various seemingly unrelated phenom-ena. First and foremost, we recapture the 'relativistic' nature of the Majorana state. We show that once the Schrödinger equation is transformed to the canonical Dirac form, a boost of Majorana states is carried out by means of a standard representation of the Lorentz group, with the mode velocity playing the role of the speed of light. The characteristics of the Majorana bound state exhibit both Lorentz contraction and time dilation, leading to renormalizations of its size and its energy separation from finite-energy states. The ability to apply Lorentz transformations enables us to obtain an exact analytical solution for Majorana-carrying domain walls moving at arbitrary velocities. Second, we find a direct mapping between the problems of a moving domain wall and graphene in crossed electric and magnetic fields. This mapping may help to gain intuition for both systems.
Our exact solution implies practical bounds on the speed with which Majorana states can be manipulated. Most importantly, Majorana states should not be accelerated to more than the effective speed of lightũ of the underlying Dirac Hamiltonian. A domain wall moving faster thanũ leads to an instability similar to superluminal particles in a dielectric that emit Cherenkov radiation. As a result, the discrete domain-wall bound states become continuous and delocalize over the entire system, which leads to a loss of the notion of Majorana bound states. Interestingly, the ensuing upper limit on the braiding frequency is in general more stringent than a naive estimate based on the magnitude of the (mini)gap (cf. Sec. IV C). To push the instability of the spectrum to as large a velocity as possible, one should work close to µ = 0. There,ũ is given by the bare mode velocity, while a finite chemical potential leads to a downward renormalization ofũ by a factor of 1 − (µ/B) 2 .
As an application of our exact solution, we study the effect of static short-range impurities which will in general lead to excitations of moving domain-wall bound states and further restrict the applicable velocities. While large velocities v ∼ũ are clearly detrimental, the Majorana states remain stable for moderate values of β = v/ũ. For an ideal linear domain wall, for which the gap increases indefinitely away from its center, excitations of the Majorana bound state are in fact strongly suppressed as a power of e −1/β 2 (as shown in Sec. IV D). For domain walls of a finite size x D , this strong suppression still holds down to velocities β ∼ ξ/x D (in terms of the size of the Majorana bound state ξ). For smaller velocities β < ξ/x D , the suppression becomes nonuniversal, depending on specifics of the domain-wall shape.
As a consequence of the universal low-energy physics underlying the formation of Majorana bound states, the applicability of our theory is rather broad. Specifically, essential aspects of the analysis carry over to various forms of smooth domain walls and (with some restrictions, see Sec. V C) to semiconductor wires as well. (The effective "speed of light" is then given by the spin-orbit velocity of the quantum wire).
The detailed understanding of boosted Majorana states in general one-dimensional systems is central for the quest to optimize Majorana manipulations. Our exact analytic solution for constant velocities provides an important stepping stone for the study of more complex braiding protocols. This opens the search for the best strategies to translate Majorana states while preserving the encoded quantum information.