#  Majorana Braiding Dynamics on Nanowires

Cássio Sozinho Amorim  Department of Applied Physics, Nagoya University,
Nagoya 464-8603, Japan  Kazuto Ebihara  Department of Applied Physics, Nagoya
University, Nagoya 464-8603, Japan  Ai Yamakage  Department of Applied
Physics, Nagoya University, Nagoya 464-8603, Japan  Yukio Tanaka  Department
of Applied Physics, Nagoya University, Nagoya 464-8603, Japan  Masatoshi Sato
<EMAIL_ADDRESS>Department of Applied Physics, Nagoya University,
Nagoya 464-8603, Japan

###### Abstract

Superconductors hosting long-sought excitations called Majorana fermions may
be ultimately used as qubits of fault-tolerant topological quantum computers.
A crucial challenge toward the topological quantum computer is to implement
quantum operation of nearly degenerate quantum states as a dynamical process
of Majorana fermions. In this paper, we investigate the braiding dynamics of
Majorana fermions on superconducting nanowires. In a finite size system, a
non-adiabatic dynamical process dominates the non-Abelian braiding that
operates qubits of Majorana fermions. Our simulations clarify how qubits
behave in the real-time braiding process, and elucidate the optimum condition
of superconducting nanowires for efficient topological quantum operation.

###### pacs:

Introduction— Recent discovery of topological matters provides a novel
platform of quantum devices. In particular , topological superconductors
naturally realize yet-to-be discovered excitations called Majorana fermions as
a collective mode in condensed matter physics Tanaka et al. (2012); Qi and
Zhang (2011); Wilczek (2009). Because of the self-antiparticle nature, the
isolated Majorana zero modes display unusual physical properties such as non-
Abelian anyon statistics, which is of extreme interest in realization of
topological quantum computer in reality.

Topological superconductivity was originally recognized in $p$-wave spin-
triplet superconductors Read and Green (2000); Ivanov (2001); Kitaev (2001);
Sato (2010), however, advance on our understanding of topological matters
enables us to design it even in a conventional $s$-wave superconducting state
Sato (2003); Fu and Kane (2008); Sato et al. (2009); Sau et al. (2010a). A
recent proposed scheme to realize Majorana fermions by using the spin-orbit
interaction and Zeeman field Sato and Fujimoto (2009); Sato et al. (2009,
2010); Sau et al. (2010a); Alicea (2010) was eventually applied to a one-
dimensional nanowire with proximity induced $s$-wave pairing Lutchyn et al.
(2010); Oreg et al. (2010), which can be fabricated by the present
experimental technique Mourik et al. (2012); Deng et al. (2012); Das et al.
(2012); Williams et al. (2012); Veldhorst et al. (2012); Rokhinson et al.
(2012). Furthermore, varieties of proposals exist in order to improve the
experimental accessibility and controllability of Majorana modes Linder et al.
(2010); Sato and Fujimoto (2010); Alicea (2012); Stanescu and Tewari (2013);
Beenakker (2013); Hassler (2013); Sau et al. (2010b); Klinovaja and Loss
(2012); Sau et al. (2012); Chevallier et al. (2013); Sau and Sarma (2012);
Nakosai et al. (2013); Mi et al. (2013).

In topological quantum computation, quantum operations of qubits are
implemented as an exchange process of Majorana zero modes. Thus, a crucial
next step toward topological quantum computer is to understand such an
operation of collective excitations as a time-dependent dynamical process.

In this paper, we investigate the braiding dynamics of Majorana zero modes on
superconducting nanowires. Generalizing proposed methods of Majorana braiding
Alicea et al. (2011); Sau et al. (2011); Liang et al. (2012); Kotetes et al.
(2013); Halperin et al. (2012); Sau et al. (2010c); Liu et al. ; Zhang et al.
(2013); Li et al. ; Weithifer et al. ; Karzig et al. (2013); Hyart et al.
(2013); Kraus et al. (2013); Liu and Lobos (2013); Chiu et al. , we consider a
simpler cruciform junction of topologically non-trivial superconducting
nanowires. This simple system functions as a quantum NOT gate of a Majorana
qubit by switching gates connecting the wires to the cross point. Using this
model, we simulate the Majorana braiding by solving the time-dependent
Bogoliubov de Genne equation for the nanowires. A non-adiabatic dynamical
process dominates the non-Abelian braiding that operates qubits of Majorana
fermions. Our simulations clarify how qubits behave in the real-time braiding
process, and elucidate the optimum condition of superconducting nanowires for
efficient topological quantum operation.

Majorana Braiding— In the low energy limit, one-dimensional topological
superconductors reduce to a one-dimensional spinless $p_{x}$-wave
superconductor. We adopt the spinless $p_{x}$-wave superconductor as a model
to analyze universal aspects of Majorana dynamics on nanowires,

$\displaystyle{\cal
H}=-\mu\sum_{x=1}^{N}c^{\dagger}_{x}c_{x}-\sum_{x=1}^{N-1}\left(\lambda
c_{x}^{\dagger}c_{x+1}+\Delta e^{i\theta}c_{x}c_{x+1}+{\rm h.c.}\right),$ (1)

where $c_{x}$ is a spinless fermion operator and $\mu$, $\lambda$ and $\Delta
e^{i\theta}$ are the chemical potential, the hopping integral, and the
$p$-wave pairing potential, respectively. ($\lambda>0$, $\Delta>0$.) There are
two different topological phases in the spinless $p_{x}$-wave superconductor
Kitaev (2001). When $|\mu|<\lambda$, the $p_{x}$-wave superconductor realizes
a topologically non-trivial superconducting state, and thus it supports a
Majorana fermion on each end. In contrast, when $|\mu|>\lambda$, it becomes a
topologically trivial state without Majorana end modes. Below, we consider
$p_{x}$-wave superconducting nanowires in the topologically non-trivial phase.

To braid the Majorana end states, we consider a cruciform junction illustrated
in Fig.1 (a), where four topologically non-trivial nanowires (wire 1, 2, 3 and
4) are connected by four gates (gate 1, 2, 3 and 4). The hopping integral
$\lambda$ and the paring potential $\Delta$ at the gates are tunable, so one
can connect (disconnect) the wires by turning on (off) these parameters at the
gates.

Figure 1:  (color online). (a) Cross-shaped topological superconducting
nanowire. The orange junctions connecting to the central site are cut/linked
through gate potentials, effectively leaving 4 independently controlled wires.
When gates 2 and 4 are connected, while 1 and 3 are shut, six Majorana end
modes $\gamma_{i}$ $(i=1,2,3,4,5,6)$ are obtained. (b) A typical energy
spectrum of MBSs as a function of $t$, where the time-dependence is given by
the gating process in Fig.2. Here we take $\Delta=\lambda$, $\mu=0.7\lambda$
and $N=20$. The color of the lines match the color of the modes in Fig.1(a).
The finite coupling of MBS on the same wire slightly lifts the zero-mode
degeneracy.

Now let us illustrate how one can exchange the Majorana end modes by switching
these gates of the cruciform junction. Initially, we prepare the configuration
of Fig.2(a), where wires 2 and 4 are connected by turning on gates 2 and 4,
while wires 1 and 3 are disconnected. There are six Majorana end states in the
initial configuration since two Majorna end states at the inner edge of wires
2 and 4 are gapped by the coupling at gates 2 and 4.

Counterclockwise exchange of the Majorana modes $\gamma_{1}$ and $\gamma_{3}$,
which are localized at the inner edges of wires 1 and 3, can be implemented as
follows: First, by turning on the gate 1 [Fig. 2(b)] and then turning off gate
2 [Fig. 2(c)], $\gamma_{1}$ moves to the inner edge of wire 2. Next,
$\gamma_{3}$ moves to the inner edge of wire 1 by turning on gate 3 [Fig.
2(d)] and then turning off gate 1 [Fig. 2(e)]. Finally, $\gamma_{1}$ moves to
the inner edge of wire 3 by turning on gate 2 [Fig. 2(f)] and then turning off
gate 3 [Fig. 2(g)]. The final gate configuration is identical to the initial
one, but $\gamma_{1}$ and $\gamma_{3}$ are exchanged.

Figure 2: (color online). The initial conditions of our system is shown in
(a), which is operated according to dimensionless gate parameters shown in
(h), being 0 a completely separated wire, and 1 a fully connected condition.
$T$ is our gate operation time parameter. Majorana $\gamma_{1}$ is moved
towards the center in (b) by connecting gate 1 and then to wire 2 in (c) by
disconnecting gate 2. Majorana $\gamma_{3}$ is moved from wire 3 to 1 by a
similar process in (d) and (e), later taking $\gamma_{1}$ from wire 2 to 3,
obtaining an interchange of $\gamma_{1}$ and $\gamma_{3}$.

In the initial configuration $(t=0)$, wire 1 and wire 3 are isolated from
others. While wires 2 and 4 are connected to each other, they are also
disconnected from the rest. Due to the finite length of wires, the mixing of
Majorana modes occurs between $\gamma_{1}$ and $\gamma_{2}$, $\gamma_{3}$ and
$\gamma_{4}$, and $\gamma_{5}$ and $\gamma_{6}$, respectively Kitaev (2001);
Cheng et al. (2009); Mizushima and Machida (2010); Cheng et al. (2010, 2011).
It induces the following effective coupling between zero modes,

$\displaystyle{\cal H}_{\rm
eff}=i\epsilon\gamma_{1}\gamma_{2}+i\epsilon\gamma_{3}\gamma_{4}+i\epsilon^{\prime}\gamma_{5}\gamma_{6},$
(2)

where the constants $\epsilon$ and $\epsilon^{\prime}$ are real because ${\cal
H}$ should be hermitian, and $\epsilon$ is larger than $\epsilon^{\prime}$
since the coupling between $\gamma_{5}$ and $\gamma_{6}$ is weaker than the
others. Assuming the standard anti-commutation relation of Majorana zero
modes, i.e. $\\{\gamma_{i},\gamma_{j}\\}=2\delta_{ij}$, ${\cal H}$ can be
recast into

$\displaystyle{\cal H}_{\rm
eff}=\epsilon^{\prime}\left(2c_{1}^{\dagger}c_{1}-1\right)+\epsilon\left(2c_{2}^{\dagger}c_{2}+2c_{3}^{\dagger}c_{3}-2\right)$
(3)

with the Dirac operators

$\displaystyle c_{1}=\frac{\gamma_{5}+i\gamma_{6}}{2},\quad
c_{2}=\frac{\gamma_{3}+i\gamma_{4}}{2},\quad
c_{3}=\frac{\gamma_{1}+i\gamma_{2}}{2},$ (4)

obeying $\\{c_{i},c_{j}^{\dagger}\\}=\delta_{ij}$. Therefore, the mixing
results in three negative energy states $|-E_{i}\rangle$ $(i=1,2,3)$ that are
annihilated by $c_{i}$,

$\displaystyle c_{i}|-E_{i}\rangle=0$ (5)

and three positve partners $|E_{i}\rangle$ that are annihilated by
$c_{i}^{\dagger}$,

$\displaystyle c_{i}^{\dagger}|E_{i}\rangle=0,$ (6)

with $E_{1}=\epsilon^{\prime}$ and $E_{2}=E_{3}=\epsilon$.

By switching the gates as described above, we can exchange the Majorana zero
modes $\gamma_{1}$ and $\gamma_{3}$. A proper gating process for the non-
Abelian braiding does not merely interchange $\gamma_{1}$ and $\gamma_{3}$,
but also provide a non-trivial relative phase between them,
$\gamma_{1}\rightarrow-\gamma_{3}$, $\gamma_{3}\rightarrow\gamma_{1},$ or
$\gamma_{1}\rightarrow\gamma_{3}$, $\gamma_{3}\rightarrow-\gamma_{1}$. In both
cases, if we exchange $\gamma_{1}$ and $\gamma_{3}$ twice, $\gamma_{1}$ and
$\gamma_{3}$ do not go back to the original, but they acquire the minus sign,
$\gamma_{1}\rightarrow-\gamma_{1},$ $\gamma_{3}\rightarrow-\gamma_{3}.$
Therefore, after exchange $\gamma_{1}$ and $\gamma_{3}$ twice, the Dirac
operators $c_{2}$ and $c_{3}$ transform into their conjugates
$-c_{2}^{\dagger}$ and $-c_{3}^{\dagger}$ as

$\displaystyle
c_{2}=\frac{\gamma_{3}+i\gamma_{4}}{2}\rightarrow\frac{-\gamma_{3}+i\gamma_{4}}{2}=-c_{2}^{\dagger}.$
$\displaystyle
c_{3}=\frac{\gamma_{1}+i\gamma_{2}}{2}\rightarrow\frac{-\gamma_{1}+i\gamma_{2}}{2}=-c_{3}^{\dagger},$
(7)

The nontrivial transformation of $c_{2}$ and $c_{3}$ implies that the negative
energy states $|-E_{2}\rangle$ and $|-E_{3}\rangle$, which are annihilated by
$c_{2}$ and $c_{3}$, respectively, end up as the positive energy partners
$|E_{2}\rangle$ and $|E_{3}\rangle$, and vice versa, after the exchange
process. In other words, if we choose these negative energy states as an
initial state, the final state is orthogonal to the initial one. This complete
interference is a direct signal of the non-Abelian anyon statistics: Indeed if
the system obeys an ordinary Abelian statistics, any exchange process results
in a phase factor for any initial state, so the final state cannot be
orthogonal to the initial one. The above exchange process defines a quantum
NOT gate for Majorana qubits $(|E_{2}\rangle,|-E_{2}\rangle)$ and
$(|E_{3}\rangle,|-E_{3}\rangle)$.

Whereas the above procedure eventually works well as is shown below, the
actual implementation needs a careful consideration for the gating. In Fig.1
(b), we show lower energy eigenvalues of the system as a function of $t$. The
eigen energies $E_{1}(t)$, $E_{2}(t)$, $E_{3}(t)$ and their negative energy
partners correspond to six Majorana zero modes of the system, where $E_{2}(t)$
and $E_{3}(t)$ are degenerate within numerical accuracy as well as their
negative energy partners are. At $t=0$, these eigen energies coincide with
$E_{i}$ $(i=1,2,3)$ in the above, $E_{i}(0)=E_{i},$ and we also have
$E_{i}(6T)=E_{i}(12T)=E_{i}$ since the system goes back to the initial
configuration at $t=6T$ and $12T$. We note here that there is no level
crossing in the energy spectrum in Fig.1 (b), as expected from the von
Neumann-Wigner theorem von Neumann and Wigner (1929). Therefore, a non-
adiabatic transition is needed to achieve the non-Abelian braiding discussed
in the above, since any state cannot be different from the original under an
adiabatic process. Namely, the gating process in Fig.2 should not be too slow.
The non-adiabatic transition is not a classical Landau-Zener transition,
because the level spacing rarely depends on $t$ and there is no level
approaching to each other at a particular time. We can also argue that a
proper gating process should not be too fast at the same time. A fast gating
process may create bulk excitations on nanowires, which may give rise to
problematic decoherence of Majorana qubits. Therefore, the gating process for
non-Abelian braiding should be performed at a proper range of speed.

Below we operate the gates 1, 2, 3 in accordance with the time-sequence
diagram in Fig.2 (h). The gating speed can be controlled by an adiabatic
parameter $T$: The gate operation becomes slower (faster) and more adiabatic
(non-adiabatic) for larger (smaller) $T$. A moderate $T$ is required to
realize the non-Abelian braiding.

Braiding Dynamics— We now numerically simulate the Majorana braiding process
in Fig.2. To numerically evaluate our system, we take each wire length to be
the same with one central site linking them. Each gate is represented as a
factor $g_{i}\in[0,1]$ ($i=1,2,3,4$) multiplying the link on the gates in real
space 111For details, see Supplementary Material.. The dynamics of the system
is described by the time-dependent Bogoliubov-de Genne equation

$\displaystyle i\hbar\frac{\partial}{\partial t}\Psi(t)={\cal H}(t)\Psi(t),$
(8)

where $\Psi(t)$ is the quasiparticle wavefunction in the Nambu representation.
The evaluation of the wavefunction during a time ${\mathit{\Delta}}t$ is given
by $\Psi(t+{\mathit{\Delta}}t)={\cal U}(t+\mathit{\Delta}t;t)\Psi(t)$ with the
time-evolution operator ${\cal U}(t+{\mathit{\Delta}}t;t)$,

$\displaystyle{\cal
U}(t+\mathit{\Delta}t,t)=T\exp\left[-i\int^{t+\mathit{\Delta}t}_{t}d\tau{\cal
H}(\tau)\right],$ (9)

which is well-approximated as ${\cal
U}(t+\mathit{\Delta}t,t)\approx\exp\left[-i{\cal
H}(t)\mathit{\Delta}t\right]$, within numerical errors for a sufficiently
short $\mathit{\Delta}t$. To achieve a correct wavefunction change in time, we
further expand the time-evolution operator in terms of Chebishev polynomials
Tal-Ezer and Kosloff (1984); Liang et al. (2012), which can be retrieved
recursively, that is

$\displaystyle\mathcal{U}(t+\mathit{\Delta}t;t)$ $\displaystyle=$
$\displaystyle\exp\left[-i\frac{\mathcal{H}(t)}{E_{0}}\mathit{\Delta}tE_{0}\right]=\exp\left[-i\tilde{H}(t)\mathit{\Delta}\tau\right]$
(10) $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}c_{k}(\mathit{\Delta}\tau)T_{k}(\tilde{H}(t)),$

where $E_{0}\equiv\max|\langle\Psi|\mathcal{H}|\Psi\rangle|$ normalizes the
Hamiltonian to avoid singularities of Chebishev polynomials and

$\displaystyle
c_{k}(\mathit{\Delta}\tau)=\left\\{\begin{array}[]{ll}J_{0}(\mathit{\Delta}\tau)&(k=0)\\\
2(-i)^{k}J_{k}(\mathit{\Delta}\tau)&(k\geq 1)\end{array}\right.$ (13)
$\displaystyle T_{0}=1,\quad T_{1}(\tilde{H})=\tilde{H},$ $\displaystyle
T_{k+1}(\tilde{H})=2\tilde{H}T_{k}(\tilde{H})-T_{k-1}(\tilde{H}),$ (14)

constitute our expansion terms. Here $\tilde{H}={\cal H}/E_{0}$,
$\mathit{\Delta}\tau=\mathit{\Delta}tE_{0}$, and $J_{k}$ are the Bessel
functions of first kind. For small $\mathit{\Delta}\tau$, the coefficients
$c_{k}(\mathit{\Delta}\tau)$ rapidly converge to zero as $k$ increases. Thus
keeping the first few expansion terms in the right hand side of Eq.(10) is
enough to reach numerically reliable results.

Figure 3: (color online). $\gamma_{2}$ braiding with $T=100/\Delta$. We take
$\Delta=\lambda$, $\mu=0.7\lambda$ and $N=20$. We start with the squared
wavefunction in (a) localized in wire 1, which can be seen to migrate towards
the center only for $\gamma_{1}$ in (b). The complete transfer of $\gamma_{1}$
to wires 2 and 3 are respectively seen in (c) and (d). Final state at (e)
indicates the return of the MBS to wire 1. (f) Projection of the wave function
on the instantaneous eigenstates, $P_{\pm i}(t)=\langle\Psi(t)|\pm
E_{i}(t)\rangle$ ($i=1,2,3$). We can observe our wavefunction begins
completely on $|E_{3}\rangle$ and after superposing on other eigenstates,
completely transfers to $|-E_{3}\rangle$.

Figure 3 is one of the main results in this paper. In Fig.3, we illustrate how
the wavefunction evolves in time in our numerical simulation of the non-
Abelian braiding. We choose $|E_{3}\rangle$ as the initial state at $t=0$, and
take $T=100/\Delta$. It demonstrates that only the inner part of the wave
function moves in time, which exactly corresponds to the movement of Majorana
mode $\gamma_{1}$. At $t=6T$, although the gate configuration goes back to the
initial one, the inner part of the wave function moves to the inner edge of
wire 3, which indicates that $\gamma_{1}$ is successfully interchanged with
$\gamma_{3}$. Then finally, the inner part goes back to the initial position
at $t=t_{\rm f}\equiv 12T$.

We project the same wavefunction into the instantaneous eigenstates $|\pm
E_{i}(t)\rangle$ $(i=1,2,3)$ in Fig.3 (f). Initially, the wavefunction
consists of only $|E_{3}\rangle$, but after the gating process starts, the
wave function quickly spreads over the eigenstates $|-E_{3}\rangle$ and $|\pm
E_{3}\rangle$. Nevertheless, the final state ends up at $|-E_{3}\rangle$, as
expected as the non-Abelian braiding process mentioned above.

If one operates the gates too slowly or too quickly, the non-Abelian braiding
fails. For example, in the slow gating with $T=100000/\Delta$, both inner and
outer Majorana modes of wire 1 move together in time, in which the state
mostly stays at the instantaneous eigenstate $|E_{3}(t)\rangle$, as expected
by the adiabatic theorem. On the other hand, for the quick gating with
$T=1/\Delta$, the wave function extends over all eigenstates $|\pm
E_{i}(t)\rangle$ $(i=1,2,3)$, and it never goes back to the initial state. In
the latter case, the wave function in the final state also spreads over
nanowires in space, which suggests that bulk modes are excited during the
gating. We exemplify these unsuccessful braiding in Fig. 4.

Figure 4: (color online). Behavior outside braiding (good) conditions. (a) and
(b) represent the evolution of the system on adiabatic limit
($T=100000/\Delta$), where both Majorana edge-states migrate to the next wire,
as can be seen from (b), consequently remaining on the same state, as can be
observed from the projection of the wavefunction on the instantaneous
eigenstates in (c). (d) and (e) account for the opposite, fast regime
($T=1/\Delta$), with bulk excitations visible in (d) and a more erratic spread
over the eigenstates in (f). We take $\Delta=\lambda$, $\mu=0.7\lambda$, and
$N=20$. Figure 5: (color online). Success rate $P_{\rm s}$ of the Majorana
braiding. Each line accounts for a different length (in site number) of each
wire (1-4), with $x$-axis being our operation time $T$, and $y$-axis the
success rate $P_{s}$. We take $\Delta=0.1\lambda$ and $\mu=0.7\lambda$. The
inset shows $P_{\rm s}$ for $\Delta=\lambda$. The red lines (80 and 20
sites/wire) draw a clear separation between its adiabatic (slow, large $T$)
and dissipative (fast, small $T$) domains, suggesting an appreciable necessity
for both size and time scale adjustment.

To quantify the non-Abelian braiding, we introduce its success rate $P_{\rm
s}$ as the probability that the final state $|\Psi(t_{\rm f})\rangle$ is found
to be the desired state $|-E_{3}\rangle$,

$\displaystyle P_{\rm s}=|\langle\Psi(t_{\rm f})|-E_{3}\rangle|^{2}$ (15)

In Fig.5, we plot the success rate $P_{\rm s}$ versus the adiabatic parameter
$T$, with various wire lengths. The data indicates that a longer wire is
desirable for the non-Abelian braiding. For longer wires, the success rate can
reach the maximum, i.e. $P_{\rm s}=1$ for large $T$. In a shorter wire, on the
other hand, a Majorana end mode is fairly coupled with the Majorana mode on
the other end, so they tend to move together, resulting in an adiabatic
process even for a moderate $T$. We also find that a quicker gating fails to
achieve the non-Abelian braiding for any length of wires, since it excites
undesirable bulk modes.

Finally, from our numerical results, we evaluate the optimal condition for
non-Abelian braiding. We note that $T$ should be larger than the inverse of
the bulk gap $1/\Delta$, not to excite bulk modes. Our numerical data
determines how large it should be. Figure 5 indicates that the lower bound $T$
is not a merely $O(1/\Delta)$, but it is evaluated as $T_{\rm min}\sim
O(10^{2}/\Delta)$. For a typical superconducting state with $\Delta=O(1)$ K,
$T_{\rm min}$ can be a few nanoseconds. On the other hand, the upper bound of
$T$ can be determined as follows. As we illustrated in the above, the non-
Abelian braiding is realized as an non-adiabatic process between Majorana
modes. Thus, for Majorana modes with energy $\epsilon$, $\epsilon T$ should
not be too large. Our numerical results imply that the success rate of the
non-Abelian braiding $P_{\rm s}$ reaches almost the maximum when $\epsilon T$
is less than $O(10)$. The latter condition can be easily met for long wires,
since $\epsilon$ scales as $\epsilon\sim e^{-N/l_{0}}$. It is also found in
Fig. 5 that the quantum limit of superconducting state, i.e.
$\Delta/\lambda=1$, requires less numbers of sites for the non-Abelian
braiding, which is preferable if one realizes topological supercondicting
wires as a chain of quantum dots.Sau and Sarma (2012)

The authors are grateful to J. D. Sau and K. T. Law for fruitful discussions.
This work was supported by the JSPS (No.25287085) and KAKENHI Grants-in-Aid
(No.22103005) from MEXT. C.S.A. is supported by MEXT Scholarship (kokuhi-
gaikokujin-ryugakusei 2013)

## References

  * Tanaka et al. (2012) Y. Tanaka, M. Sato, and N. Nagaosa, J. Phys. Soc. Jpn. 81, 011013 (2012).
  * Qi and Zhang (2011) X. L. Qi and S. C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Wilczek (2009) F. Wilczek, Nat. Phys. 5, 614 (2009).
  * Read and Green (2000) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * Ivanov (2001) D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
  * Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Sato (2010) M. Sato, Phys. Rev. B 81, 220504(R) (2010).
  * Sato (2003) M. Sato, Phys. Lett. B 575, 126 (2003).
  * Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Sato et al. (2009) M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. 103, 020401 (2009).
  * Sau et al. (2010a) J. D. Sau, R. M. Lutchyn, S. Tewari, and S. D. Sarma, Phys. Rev. Lett. 104, 040502 (2010a).
  * Sato and Fujimoto (2009) M. Sato and S. Fujimoto, Phys. Rev. B 79, 094504 (2009).
  * Sato et al. (2010) M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B 82, 134521 (2010).
  * Alicea (2010) J. Alicea, Phys. Rev. B 81, 125318 (2010).
  * Lutchyn et al. (2010) R. M. Lutchyn, J. D. Sau, and S. D. Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg et al. (2010) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Mourik et al. (2012) V. Mourik, K. Zuo, S. Frolov, S. Plissard, E. Bakkers, and L. Kouwenhoven, Science 336, 1003 (2012).
  * Deng et al. (2012) M. Deng, C. Yu, G. Huang, M. Larsson, P. Caro, and H. Xu, Nano Lett. 12, 6414 (2012).
  * Das et al. (2012) A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Nat. Phys. 8, 887 (2012).
  * Williams et al. (2012) J. Williams, A. Bestwick, P. Gallagher, S. Hong, Y. Cui, A. Bleich, J. Analytis, I. Fisher, and D. Goldhaber-Gordon, Phys. Rev. Lett. 109, 056803 (2012).
  * Veldhorst et al. (2012) M. Veldhorst, M. Snelder, M. hoek, T. Gang, V. Guduru, X. L. Wang, U. Zeitler, W. G. van der Wiel, A. A. Golubov, H. Hilgenkamp, et al., Nat. Mater. 11, 417 (2012).
  * Rokhinson et al. (2012) L. P. Rokhinson, X. Liu, and J. K. Furdyna, Nat. Phys. 8, 795 (2012).
  * Linder et al. (2010) J. Linder, Y. Tanaka, T. Yokoyama, A. Sudbo, and N. Nagaosa, Phys. Rev. Lett. 104, 067001 (2010).
  * Sato and Fujimoto (2010) M. Sato and S. Fujimoto, Phys. Rev. Lett. 105, 217001 (2010).
  * Alicea (2012) J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
  * Stanescu and Tewari (2013) T. D. Stanescu and S. Tewari, J. Phys.: Condens. Matter 25, 233201 (2013).
  * Beenakker (2013) C. W. J. Beenakker, Annu. Rev. Con. Mat. Phys. 4, 113 (2013).
  * Hassler (2013) F. Hassler, Quantum Information Processing. Lecture notes of the 44th IFF Spring School 2013 (2013).
  * Sau et al. (2010b) J. D. Sau, S. Tewari, R. M. Lutchyn, T. D. Stanescu, and S. D. Sarma, Phys. Rev. B 82, 214509 (2010b).
  * Klinovaja and Loss (2012) J. Klinovaja and D. Loss, Phys. Rev. B 86, 085408 (2012).
  * Sau et al. (2012) J. D. Sau, C. H. Lin, H.-Y. Hui, and S. D. Sarma, Phys. Rev. Lett. 108, 067001 (2012).
  * Chevallier et al. (2013) D. Chevallier, D. Sticlet, P. Simon, and C. Bena, Phys. Rev. B 87, 165414 (2013).
  * Sau and Sarma (2012) J. D. Sau and S. D. Sarma, Nat. Comm. 3, 964 (2012).
  * Nakosai et al. (2013) S. Nakosai, J. C. Budich, Y. Tanaka, B. Trauzettel, and N. Nagaosa, Phys. Rev. Lett. 110, 117002 (2013).
  * Mi et al. (2013) S. Mi, D. I. Pikulin, M. Wimmer, and C. W. J. Beenakker, Phys. Rev. B 87, 241405 (2013).
  * Alicea et al. (2011) J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Nat. Phys. 7, 412 (2011).
  * Sau et al. (2011) J. Sau, D. Clarke, and S. Tewari, Phys. Rev. B 84, 094505 (2011).
  * Liang et al. (2012) Q. Liang, Z. Wang, and X. Hu, Euro Physics Letters 99, 50004 (2012).
  * Kotetes et al. (2013) P. Kotetes, G. Schön, and A. Shnirman, Journal of the Korean Physical Society 62, 1558 (2013).
  * Halperin et al. (2012) B. I. Halperin, Y. Oreg, A. Stern, G. Refael, J. Alicea, and F. von Oppen, Phys. Rev. B 85, 144501 (2012).
  * Sau et al. (2010c) J. D. Sau, S. Tewari, and S. D. Sarma, Phys. Rev. A 82, 052322 (2010c).
  * (42) X. J. Liu, C. L. M. Wong, and K. T. Law, eprint arXiv:1304.3765.
  * Zhang et al. (2013) F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 111, 056403 (2013).
  * (44) J. Li, T. Neupert, B. A. Bernevig, and A. Yazdani, eprint arXiv:1404.4058.
  * (45) L. Weithifer, P. Recher, and T. L. Schimidt, eprint arXiv:1309.4126.
  * Karzig et al. (2013) T. Karzig, G. Refael, and F. von Oppen, Phys. Rev. X 3, 041017 (2013).
  * Hyart et al. (2013) T. Hyart, B. van Heck, I. C. Fulga, M. Burrello, A. R. Akhmerov, and C. W. J. Beenakker, Phys. Rev. B 88, 035121 (2013).
  * Kraus et al. (2013) C. V. Kraus, P. Zoller, and M. A. Baranov, Phys. Rev. Lett. 111, 203001 (2013).
  * Liu and Lobos (2013) X.-J. Liu and A. M. Lobos, Phys. Rev. B 87, 060504 (2013).
  * (50) C.-K. Chiu, M. Vazifeh, and M. Franz, eprint arXiv:1403.0033v1.
  * Cheng et al. (2009) M. Cheng, R. M. Lutchyn, V. Galitski, and S. D. Sarma, Phys. Rev. Lett. 103, 107001 (2009).
  * Mizushima and Machida (2010) T. Mizushima and K. Machida, Phys. Rev. A 82, 023624 (2010).
  * Cheng et al. (2010) M. Cheng, R. M. Lutchyn, V. Galitski, and S. D. Sarma, Phys. Rev. B 82, 094504 (2010).
  * Cheng et al. (2011) M. Cheng, V. Galitski, and S. D. Sarma, Phys. Rev. B 84, 104529 (2011).
  * von Neumann and Wigner (1929) J. von Neumann and E. P. Wigner, Z. Phys. 30, 467 (1929).
  * Tal-Ezer and Kosloff (1984) H. Tal-Ezer and R. Kosloff, Jour. Chem. Phys. 81 (1984).
  * Göddeke et al. (2007) D. Göddeke, R. Strzodka, and S. Turek, International Journal of Parallel, Emergent and Distributed Systems (IJPEDS), Special issue: Applied parallel computing 22, 221 (2007).
  * Baboulin et al. (2009) M. Baboulin, A. Buttari, J. Dongarra, J. Kurzak, J. Langou, J. Langou, P. Luszczek, and S. Tomov, Comp. Phys. Comm. 180, 2526 (2009).
  * Grand et al. (2013) S. L. Grand, A. W. Götz, and R. C. Walker, Comp. Phys. Comm. 184, 374 (2013).

Supplementary Material

##  Appendix A Calculation Method

The Hamiltonian in Eq. (1) can be recast under Bogoliubov-de Gennes
representation as

$\displaystyle\mathcal{H}^{\mathrm{BdG}}$ $\displaystyle=$
$\displaystyle-\left[\frac{\mu}{2}+\frac{\lambda}{2}\left(\cos k_{x}+\cos
k_{y}\right)\right]\tau_{z}$ (1)
$\displaystyle-\frac{\Delta}{2}e^{i\theta}\left(\sin k_{x}\tau_{y}+\sin
k_{y}\tau_{x}\right),$

which in turn can be rewritten in real space in matrix form, defining our
Hamiltonian as
$H=\sum_{ij}c^{\dagger}_{i}\mathcal{H}_{{\bm{i}},{\bm{j}}}c_{j}$ on Nambu
basis,

$\displaystyle\mathcal{H}_{{\bm{i}},{\bm{i}}}=\left(\begin{array}[]{cc}-{\mu}&0\\\
0&{\mu}\end{array}\right),$ (4)
$\displaystyle\mathcal{H}_{{\bm{i}}\pm\hat{\bm{e}}_{x},{\bm{i}}}=\left(\begin{array}[]{cc}-\lambda/2&\mp\Delta
e^{i\theta}/2\\\ \pm\Delta e^{-i\theta}/2&\lambda/2\end{array}\right),$ (7)
$\displaystyle\mathcal{H}_{{\bm{i}}\pm\hat{\bm{e}}_{y},{\bm{i}}}=\left(\begin{array}[]{cc}-\lambda/2&\mp\Delta
e^{i\theta}/2\\\ \mp\Delta e^{-i\theta}/2&\lambda/2\end{array}\right),$ (10)

where ${\bm{i}}$ represents site position, assuming lattice constant to be 1.
To numerically evaluate our system, we take each wire length to be $n$ sites
long with one central site linking them. Each gate is represented as a factor
$g_{i}\in[0,1]$ ($i=1,2,3,4$) multiplying the Hamiltonian elements between the
central site and its neighbors in real space,

$\displaystyle{\cal
H}_{\hat{\bm{e}}_{y},0}=g_{1}\left(\begin{array}[]{cc}-\lambda/2&-ie^{i\theta}\Delta/2\\\
-ie^{-i\theta}\Delta/2&\lambda/2\end{array}\right),$ (13) $\displaystyle{\cal
H}_{-\hat{\bm{e}}_{x},0}=g_{2}\left(\begin{array}[]{cc}-\lambda/2&e^{i\theta}\Delta/2\\\
-e^{-i\theta}\Delta/2&\lambda/2\end{array}\right),$ (16) $\displaystyle{\cal
H}_{-\hat{\bm{e}}_{y},0}=g_{3}\left(\begin{array}[]{cc}-\lambda/2&ie^{-i\theta}\Delta/2\\\
ie^{-i\theta}\Delta/2&\lambda/2\end{array}\right),$ (19) $\displaystyle{\cal
H}_{\hat{\bm{e}}_{x},0}=g_{4}\left(\begin{array}[]{cc}-\lambda/2&-e^{i\theta}\Delta/2\\\
e^{-i\theta}\Delta/2&\lambda/2\end{array}\right).$ (22)

For a linear gate operation in the interval of time $T$ as described in this
work, a gate being turned on (off) is taken to evolve as $g_{i}=t/T$
($g_{i}=1-t/T$), counting the time $t$ from the beginning of the operation.
More general functions may be used for operating both gates smoothly at the
same time. Finally, to achieve the wavefunction change in time, we expand the
time-evolution operator in terms of Chebishev polynomialsTal-Ezer and Kosloff
(1984), which can be retrieved recursively, that is

$\displaystyle\mathcal{U}(t+\mathit{\Delta}t;t)$ $\displaystyle=$
$\displaystyle\exp\left[-i\frac{\mathcal{H}(t)}{E_{0}}\mathit{\Delta}tE_{0}\right]=\exp\left[-i\tilde{H}(t)\mathit{\Delta}\tau\right]$
(23) $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}c_{k}(\mathit{\Delta}\tau)T_{k}(\tilde{H}(t)),$

$E_{0}\equiv\max|\langle\Psi|\mathcal{H}|\Psi\rangle|$ normalizes the
Hamiltonian to avoid singularities and

$\displaystyle
c_{k}(\mathit{\Delta}\tau)=\left\\{\begin{array}[]{ll}J_{0}(\mathit{\Delta}\tau)&(k=0)\\\
2(-i)^{k}J_{k}(\mathit{\Delta}\tau)&(k\geq 1)\end{array}\right.$ (26)
$\displaystyle T_{0}=1,\quad T_{1}(\tilde{H})=\tilde{H},$ $\displaystyle
T_{k+1}(\tilde{H})=2\tilde{H}T_{k}(\tilde{H})-T_{k-1}(\tilde{H}),$ (27)

constitute our expansion terms. $J_{k}$ are the Bessel functions of first
kind, and their value is used to choose truncation point under double
precision. While one can apply the above operator successively, obtaining the
desired wavefunction in time $t$, this process is needlessly slow if done
completely in double precision. Mixed precisionGöddeke et al. (2007); Baboulin
et al. (2009); Grand et al. (2013) can be efficiently applied for a boost in
calculation speed if instead of calculating the wavefunction, only its
variation is evaluated. Explicitly, this can be illustrated on the following
steps:

$\displaystyle\Psi(t+\mathit{\Delta}t)_{\rm
dl}=\mathcal{U}(t+\mathit{\Delta}t;t)\Psi(t)_{\rm dl}$
$\displaystyle\rightarrow\mathit{\Delta}\Psi_{\rm
dl}=(\mathcal{U}(t+\mathit{\Delta}t;t)-1)\Psi(t)_{\rm dl}$ (28)
$\displaystyle\mathit{\Delta}\Psi=_{\rm
dl}\left(\sum_{k=0}^{k^{\prime}}c_{k}T_{k}(\tilde{H})-1\right)_{\rm
fl}\Psi(t)_{\rm fl}$ $\displaystyle=_{\rm dl}\left[(c_{0}-1)_{\rm
fl}+\sum_{k=1}^{k^{\prime}}c_{k,{\rm fl}}T_{k}(\tilde{H})_{\rm
fl}\right]\Psi(t)_{\rm fl}$ $\displaystyle\Psi(t+\mathit{\Delta}t)=_{\rm
dl}\Psi(t)_{\rm dl}+\mathit{\Delta}\Psi_{\rm dl}$ (29)

Here the subscript ${\rm dl}$ (${\rm fl}$) represents a double (single)
precision conversion/variable. For example, equation 28 represents a common
time-evolution process done with double precision variables. The following
expression claims a conversion of $\Psi(t)$ from double to single precision,
as well as a single precision expansion of the time-evolution operator, while
$=_{\rm dl}$ implies that these data should be converted to double precision
when adding to build up $\mathit{\Delta}\Psi$. Similarly, on the following
line we point that the zero-order term should be changed to single precision
after calculating it in double precision, and each other expansion term is
calculated with single precision $T_{k}$ and $c_{k}$. Shortly, each term in
our expansion is calculated as a vector in single precision, being later added
as a double precision variation to the double precision wavefunction. Note
that the wavefunction is converted to single precision for time evolution, but
computed as double precision in the end of each step. In fact, we start our
evaluation with an eigenstate of the Hamiltonian taken in double precision,
and only its variation, which corresponds to most of the computations, is
found in single precision steps.

It is important to note that this method relies on the constraint of small
$dt$ to work, as well as small and smooth wavefunction variation, otherwise
single precision computation of the expansion terms may cumulate a large
error. Nevertheless, this constraint is also required for the very numeric
expansion of the time-evolution operator, in order to ignore its intrinsic
time-ordering operator to a good approximation. In other words, the
possibility to evaluate time-evolution of the wavefunction in real space-time
already gives us the possibility of a mixed precision method. Concretely, in
our case each $\mathit{\Delta}\Psi\sim 10^{-7}$, which in single precision
allows for good enough numeric results in the $10^{-7}\sim 10^{-14}$ range,
which accounts for first order terms, as well as $10^{-9}\sim 10^{-16}$ for
the next order, all of them well fit in the limit of double precision, up to
$10^{-15}$. Therefore, summing these terms up on double precision avoids
greater errors from ignoring their contribution, which forcedly would happen
if they were added completely in single precision.