Majorana Braiding Dynamics on Nanowires
Cássio Sozinho Amorim, Kazuto Ebihara, Ai Yamakage, Yukio Tanaka, Masatoshi Sato
MMajorana Braiding Dynamics on Nanowires
C´assio Sozinho Amorim, Kazuto Ebihara, Ai Yamakage, Yukio Tanaka, and Masatoshi Sato ∗ Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan (Dated: November 10, 2018)Superconductors hosting long-sought excitations called Majorana fermions may be ultimatelyused as qubits of fault-tolerant topological quantum computers. A crucial challenge toward thetopological quantum computer is to implement quantum operation of nearly degenerate quantumstates as a dynamical process of Majorana fermions. In this paper, we investigate the braidingdynamics 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 Majoranafermions. Our simulations clarify how qubits behave in the real-time braiding process, and elucidatethe optimum condition of superconducting nanowires for efficient topological quantum operation.
PACS numbers:
Introduction—
Recent discovery of topological mattersprovides a novel platform of quantum devices. In partic-ular , topological superconductors naturally realize yet-to-be discovered excitations called Majorana fermions asa collective mode in condensed matter physics [1–3]. Be-cause of the self-antiparticle nature, the isolated Majo-rana zero modes display unusual physical properties suchas non-Abelian anyon statistics, which is of extreme in-terest in realization of topological quantum computer inreality.Topological superconductivity was originally recog-nized in p -wave spin-triplet superconductors [4–7], how-ever, advance on our understanding of topological mat-ters enables us to design it even in a conventional s -wavesuperconducting state [8–11]. A recent proposed schemeto realize Majorana fermions by using the spin-orbit in-teraction and Zeeman field [10–14] was eventually appliedto a one-dimensional nanowire with proximity induced s -wave pairing [15, 16], which can be fabricated by thepresent experimental technique [17–22]. Furthermore,varieties of proposals exist in order to improve the ex-perimental accessibility and controllability of Majoranamodes [23–35].In topological quantum computation, quantum opera-tions of qubits are implemented as an exchange process ofMajorana zero modes. Thus, a crucial next step towardtopological quantum computer is to understand such anoperation of collective excitations as a time-dependentdynamical process.In this paper, we investigate the braiding dynamicsof Majorana zero modes on superconducting nanowires.Generalizing proposed methods of Majorana braiding[36–50], we consider a simpler cruciform junction of topo-logically non-trivial superconducting nanowires. Thissimple system functions as a quantum NOT gate of a Ma-jorana qubit by switching gates connecting the wires tothe cross point. Using this model, we simulate the Majo-rana braiding by solving the time-dependent Bogoliubovde Genne equation for the nanowires. A non-adiabaticdynamical process dominates the non-Abelian braiding that operates qubits of Majorana fermions. Our simula-tions clarify how qubits behave in the real-time braidingprocess, and elucidate the optimum condition of super-conducting nanowires for efficient topological quantumoperation. Majorana Braiding—
In the low energy limit, one-dimensional topological superconductors reduce to aone-dimensional spinless p x -wave superconductor. Weadopt the spinless p x -wave superconductor as a modelto analyze universal aspects of Majorana dynamics onnanowires, H = − µ N (cid:88) x =1 c † x c x − N − (cid:88) x =1 (cid:0) λc † x c x +1 + ∆ e iθ c x c x +1 + h . c . (cid:1) , (1)where c x is a spinless fermion operator and µ , λ and∆ e iθ are the chemical potential, the hopping integral,and the p -wave pairing potential, respectively. ( λ > > p x -wave superconductor [6]. When | µ | < λ ,the p x -wave superconductor realizes a topologically non-trivial superconducting state, and thus it supports a Ma-jorana fermion on each end. In contrast, when | µ | > λ ,it becomes a topologically trivial state without Majoranaend modes. Below, we consider p x -wave superconductingnanowires in the topologically non-trivial phase.To braid the Majorana end states, we consider a cru-ciform junction illustrated in Fig.1 (a), where four topo-logically non-trivial nanowires (wire 1, 2, 3 and 4) areconnected by four gates (gate 1, 2, 3 and 4). The hop-ping integral λ and the paring potential ∆ at the gatesare tunable, so one can connect (disconnect) the wires byturning on (off) these parameters at the gates.Now let us illustrate how one can exchange the Ma-jorana end modes by switching these gates of the cruci-form junction. Initially, we prepare the configuration ofFig.2(a), where wires 2 and 4 are connected by turningon gates 2 and 4, while wires 1 and 3 are disconnected.There are six Majorana end states in the initial configu-ration since two Majorna end states at the inner edge of a r X i v : . [ c ond - m a t . m e s - h a ll ] M a y FIG. 1: (color online). (a) Cross-shaped topological super-conducting nanowire. The orange junctions connecting to thecentral site are cut/linked through gate potentials, effectivelyleaving 4 independently controlled wires. When gates 2 and 4are connected, while 1 and 3 are shut, six Majorana end modes γ i ( i = 1 , , , , ,
6) are obtained. (b) A typical energy spec-trum of MBSs as a function of t , where the time-dependenceis given by the gating process in Fig.2. Here we take ∆ = λ , µ = 0 . λ and N = 20. The color of the lines match the colorof the modes in Fig.1(a). The finite coupling of MBS on thesame wire slightly lifts the zero-mode degeneracy. wires 2 and 4 are gapped by the coupling at gates 2 and4. Counterclockwise exchange of the Majorana modes γ and γ , which are localized at the inner edges of wires 1and 3, can be implemented as follows: First, by turningon the gate 1 [Fig. 2(b)] and then turning off gate 2[Fig. 2(c)], γ moves to the inner edge of wire 2. Next, γ moves to the inner edge of wire 1 by turning on gate3 [Fig. 2(d)] and then turning off gate 1 [Fig. 2(e)].Finally, γ moves to the inner edge of wire 3 by turning ongate 2 [Fig. 2(f)] and then turning off gate 3 [Fig. 2(g)].The final gate configuration is identical to the initial one,but γ and γ are exchanged.In the initial configuration ( t = 0), wire 1 and wire 3are isolated from others. While wires 2 and 4 are con-nected to each other, they are also disconnected fromthe rest. Due to the finite length of wires, the mixingof Majorana modes occurs between γ and γ , γ and γ , and γ and γ , respectively [6, 51–54]. It induces thefollowing effective coupling between zero modes, H eff = i(cid:15)γ γ + i(cid:15)γ γ + i(cid:15) (cid:48) γ γ , (2)where the constants (cid:15) and (cid:15) (cid:48) are real because H shouldbe hermitian, and (cid:15) is larger than (cid:15) (cid:48) since the couplingbetween γ and γ is weaker than the others. Assumingthe standard anti-commutation relation of Majorana zeromodes, i.e. { γ i , γ j } = 2 δ ij , H can be recast into H eff = (cid:15) (cid:48) (cid:16) c † c − (cid:17) + (cid:15) (cid:16) c † c + 2 c † c − (cid:17) (3)with the Dirac operators c = γ + iγ , c = γ + iγ , c = γ + iγ , (4) FIG. 2: (color online). The initial conditions of our systemis shown in (a), which is operated according to dimensionlessgate parameters shown in (h), being 0 a completely separatedwire, and 1 a fully connected condition. T is our gate op-eration time parameter. Majorana γ is moved towards thecenter in (b) by connecting gate 1 and then to wire 2 in (c)by disconnecting gate 2. Majorana γ is moved from wire 3to 1 by a similar process in (d) and (e), later taking γ fromwire 2 to 3, obtaining an interchange of γ and γ . obeying { c i , c † j } = δ ij . Therefore, the mixing results inthree negative energy states | − E i (cid:105) ( i = 1 , ,
3) that areannihilated by c i , c i | − E i (cid:105) = 0 (5)and three positve partners | E i (cid:105) that are annihilated by c † i , c † i | E i (cid:105) = 0 , (6)with E = (cid:15) (cid:48) and E = E = (cid:15) .By switching the gates as described above, we can ex-change the Majorana zero modes γ and γ . A propergating process for the non-Abelian braiding does notmerely interchange γ and γ , but also provide a non-trivial relative phase between them, γ → − γ , γ → γ , or γ → γ , γ → − γ . In both cases, if we exchange γ and γ twice, γ and γ do not go back to the original,but they acquire the minus sign, γ → − γ , γ → − γ . Therefore, after exchange γ and γ twice, the Dirac op-erators c and c transform into their conjugates − c † and − c † as c = γ + iγ → − γ + iγ − c † .c = γ + iγ → − γ + iγ − c † , (7)The nontrivial transformation of c and c implies thatthe negative energy states | − E (cid:105) and | − E (cid:105) , which areannihilated by c and c , respectively, end up as the pos-itive energy partners | E (cid:105) and | E (cid:105) , and vice versa, afterthe exchange process. In other words, if we choose thesenegative energy states as an initial state, the final state isorthogonal to the initial one. This complete interferenceis a direct signal of the non-Abelian anyon statistics: In-deed if the system obeys an ordinary Abelian statistics,any exchange process results in a phase factor for anyinitial state, so the final state cannot be orthogonal tothe initial one. The above exchange process defines aquantum NOT gate for Majorana qubits ( | E (cid:105) , | − E (cid:105) )and ( | E (cid:105) , | − E (cid:105) ).Whereas the above procedure eventually works well asis shown below, the actual implementation needs a care-ful consideration for the gating. In Fig.1 (b), we showlower energy eigenvalues of the system as a function of t . The eigen energies E ( t ), E ( t ), E ( t ) and their neg-ative energy partners correspond to six Majorana zeromodes of the system, where E ( t ) and E ( t ) are degen-erate within numerical accuracy as well as their negativeenergy partners are. At t = 0, these eigen energies co-incide with E i ( i = 1 , ,
3) in the above, E i (0) = E i , and we also have E i (6 T ) = E i (12 T ) = E i since thesystem goes back to the initial configuration at t = 6 T and 12 T . We note here that there is no level crossing inthe energy spectrum in Fig.1 (b), as expected from thevon Neumann-Wigner theorem [55]. Therefore, a non-adiabatic transition is needed to achieve the non-Abelianbraiding discussed in the above, since any state cannotbe different from the original under an adiabatic pro-cess. Namely, the gating process in Fig.2 should not betoo slow. The non-adiabatic transition is not a classi-cal Landau-Zener transition, because the level spacingrarely depends on t and there is no level approaching toeach other at a particular time. We can also argue thata proper gating process should not be too fast at thesame time. A fast gating process may create bulk exci-tations on nanowires, which may give rise to problematicdecoherence of Majorana qubits. Therefore, the gatingprocess for non-Abelian braiding should be performed ata proper range of speed.Below we operate the gates 1, 2, 3 in accordance withthe time-sequence diagram in Fig.2 (h). The gating speedcan be controlled by an adiabatic parameter T : Thegate operation becomes slower (faster) and more adia-batic (non-adiabatic) for larger (smaller) T . A moderate T is required to realize the non-Abelian braiding. Braiding Dynamics—
We now numerically simulatethe Majorana braiding process in Fig.2. To numericallyevaluate our system, we take each wire length to be thesame with one central site linking them. Each gate isrepresented as a factor g i ∈ [0 ,
1] ( i = 1 , , ,
4) multi-plying the link on the gates in real space [60]. The dy-namics of the system is described by the time-dependentBogoliubov-de Genne equation i ¯ h ∂∂t Ψ( t ) = H ( t )Ψ( t ) , (8) where Ψ( t ) is the quasiparticle wavefunction in theNambu representation. The evaluation of the wavefunc-tion during a time ∆ t is given by Ψ( t + ∆ t ) = U ( t + ∆ t ; t )Ψ( t ) with the time-evolution operator U ( t + ∆ t ; t ), U ( t + ∆ t, t ) = T exp (cid:34) − i (cid:90) t + ∆ tt dτ H ( τ ) (cid:35) , (9)which is well-approximated as U ( t + ∆ t, t ) ≈ exp [ − i H ( t ) ∆ t ], within numerical errors for a sufficientlyshort ∆ t . To achieve a correct wavefunction change intime, we further expand the time-evolution operator interms of Chebishev polynomials [38, 56], which can beretrieved recursively, that is U ( t + ∆ t ; t ) = exp (cid:20) − i H ( t ) E ∆ tE (cid:21) = exp (cid:104) − i ˜ H ( t ) ∆ τ (cid:105) = ∞ (cid:88) k =0 c k ( ∆ τ ) T k ( ˜ H ( t )) , (10)where E ≡ max |(cid:104) Ψ |H| Ψ (cid:105)| normalizes the Hamiltonianto avoid singularities of Chebishev polynomials and c k ( ∆ τ ) = (cid:26) J ( ∆ τ ) ( k = 0)2( − i ) k J k ( ∆ τ ) ( k ≥
1) (11) T = 1 , T ( ˜ H ) = ˜ H,T k +1 ( ˜ H ) = 2 ˜ HT k ( ˜ H ) − T k − ( ˜ H ) , (12)constitute our expansion terms. Here ˜ H = H /E , ∆ τ = ∆ tE , and J k are the Bessel functions of first kind.For small ∆ τ , the coefficients c k ( ∆ τ ) rapidly converge tozero as k increases. Thus keeping the first few expansionterms in the right hand side of Eq.(10) is enough to reachnumerically reliable results.Figure 3 is one of the main results in this paper. InFig.3, we illustrate how the wavefunction evolves in timein our numerical simulation of the non-Abelian braid-ing. We choose | E (cid:105) as the initial state at t = 0, andtake T = 100 / ∆. It demonstrates that only the innerpart of the wave function moves in time, which exactlycorresponds to the movement of Majorana mode γ . At t = 6 T , although the gate configuration goes back to theinitial one, the inner part of the wave function moves tothe inner edge of wire 3, which indicates that γ is suc-cessfully interchanged with γ . Then finally, the innerpart goes back to the initial position at t = t f ≡ T .We project the same wavefunction into the instanta-neous eigenstates | ± E i ( t ) (cid:105) ( i = 1 , ,
3) in Fig.3 (f).Initially, the wavefunction consists of only | E (cid:105) , but af-ter the gating process starts, the wave function quicklyspreads over the eigenstates | − E (cid:105) and | ± E (cid:105) . Never-theless, the final state ends up at | − E (cid:105) , as expected asthe non-Abelian braiding process mentioned above.If one operates the gates too slowly or too quickly, thenon-Abelian braiding fails. For example, in the slow gat-ing with T = 100000 / ∆, both inner and outer Majorana FIG. 3: (color online). γ braiding with T = 100 / ∆. We take∆ = λ , µ = 0 . λ and N = 20. We start with the squaredwavefunction in (a) localized in wire 1, which can be seen tomigrate towards the center only for γ in (b). The completetransfer of γ to wires 2 and 3 are respectively seen in (c) and(d). Final state at (e) indicates the return of the MBS to wire1. (f) Projection of the wave function on the instantaneouseigenstates, P ± i ( t ) = (cid:104) Ψ( t ) | ± E i ( t ) (cid:105) ( i = 1 , , | E (cid:105) and aftersuperposing on other eigenstates, completely transfers to | − E (cid:105) . modes of wire 1 move together in time, in which the statemostly stays at the instantaneous eigenstate | E ( t ) (cid:105) , asexpected by the adiabatic theorem. On the other hand,for the quick gating with T = 1 / ∆, the wave functionextends over all eigenstates | ± E i ( t ) (cid:105) ( i = 1 , , FIG. 4: (color online). Behavior outside braiding (good) con-ditions. (a) and (b) represent the evolution of the systemon adiabatic limit ( T = 100000 / ∆), where both Majoranaedge-states migrate to the next wire, as can be seen from (b),consequently remaining on the same state, as can be observedfrom the projection of the wavefunction on the instantaneouseigenstates in (c). (d) and (e) account for the opposite, fastregime ( T = 1 / ∆), with bulk excitations visible in (d) anda more erratic spread over the eigenstates in (f). We take∆ = λ , µ = 0 . λ , and N = 20. To quantify the non-Abelian braiding, we introduceits success rate P s as the probability that the final state P s T ∆ P s T ∆ FIG. 5: (color online). Success rate P s of the Majorana braid-ing. Each line accounts for a different length (in site num-ber) of each wire (1-4), with x -axis being our operation time T , and y -axis the success rate P s . We take ∆ = 0 . λ and µ = 0 . λ . The inset shows P s for ∆ = λ . The red lines(80 and 20 sites/wire) draw a clear separation between itsadiabatic (slow, large T ) and dissipative (fast, small T ) do-mains, suggesting an appreciable necessity for both size andtime scale adjustment. | Ψ( t f ) (cid:105) is found to be the desired state | − E (cid:105) , P s = |(cid:104) Ψ( t f ) | − E (cid:105)| (13)In Fig.5, we plot the success rate P s versus the adiabaticparameter T , with various wire lengths. The data indi-cates that a longer wire is desirable for the non-Abelianbraiding. For longer wires, the success rate can reach themaximum, i.e. P s = 1 for large T . In a shorter wire, onthe other hand, a Majorana end mode is fairly coupledwith the Majorana mode on the other end, so they tendto move together, resulting in an adiabatic process evenfor a moderate T . We also find that a quicker gatingfails to achieve the non-Abelian braiding for any lengthof wires, since it excites undesirable bulk modes.Finally, from our numerical results, we evaluate the op-timal condition for non-Abelian braiding. We note that T should be larger than the inverse of the bulk gap 1 / ∆,not to excite bulk modes. Our numerical data deter-mines how large it should be. Figure 5 indicates that thelower bound T is not a merely O (1 / ∆), but it is evalu-ated as T min ∼ O (10 / ∆). For a typical superconduct-ing state with ∆ = O (1) K, T min can be a few nanosec-onds. On the other hand, the upper bound of T can bedetermined as follows. As we illustrated in the above,the non-Abelian braiding is realized as an non-adiabaticprocess between Majorana modes. Thus, for Majoranamodes with energy (cid:15) , (cid:15)T should not be too large. Ournumerical results imply that the success rate of the non-Abelian braiding P s reaches almost the maximum when (cid:15)T is less than O (10). The latter condition can be easilymet for long wires, since (cid:15) scales as (cid:15) ∼ e − N/l . It is alsofound in Fig. 5 that the quantum limit of superconduct-ing state, i.e. ∆ /λ = 1, requires less numbers of sitesfor the non-Abelian braiding, which is preferable if onerealizes topological supercondicting wires as a chain ofquantum dots.[33]The authors are grateful to J. D. Sau and K. T.Law for fruitful discussions. This work was supportedby the JSPS (No.25287085) and KAKENHI Grants-in-Aid (No.22103005) from MEXT. C.S.A. is supported byMEXT Scholarship (kokuhi-gaikokujin-ryugakusei 2013) ∗ Electronic address: [email protected][1] Y. Tanaka, M. Sato, and N. Nagaosa, J. Phys. Soc. Jpn. , 011013 (2012).[2] X. L. Qi and S. C. Zhang, Rev. Mod. Phys. , 1057(2011).[3] F. Wilczek, Nat. Phys. , 614 (2009).[4] N. Read and D. Green, Phys. Rev. B , 10267 (2000).[5] D. A. Ivanov, Phys. Rev. Lett. , 268 (2001).[6] A. Y. Kitaev, Physics-Uspekhi , 131 (2001).[7] M. Sato, Phys. Rev. B , 220504(R) (2010).[8] M. Sato, Phys. Lett. B , 126 (2003).[9] L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407(2008).[10] M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. , 020401 (2009).[11] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. D. Sarma,Phys. Rev. Lett. , 040502 (2010).[12] M. Sato and S. Fujimoto, Phys. Rev. B , 094504(2009).[13] M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B , 134521 (2010).[14] J. Alicea, Phys. Rev. B , 125318 (2010).[15] R. M. Lutchyn, J. D. Sau, and S. D. Sarma, Phys. Rev.Lett. , 077001 (2010).[16] Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. , 177002 (2010).[17] V. Mourik, K. Zuo, S. Frolov, S. Plissard, E. Bakkers,and L. Kouwenhoven, Science , 1003 (2012).[18] M. Deng, C. Yu, G. Huang, M. Larsson, P. Caro, andH. Xu, Nano Lett. , 6414 (2012).[19] A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, andH. Shtrikman, Nat. Phys. , 887 (2012).[20] J. Williams, A. Bestwick, P. Gallagher, S. Hong, Y. Cui,A. Bleich, J. Analytis, I. Fisher, and D. Goldhaber-Gordon, Phys. Rev. Lett. , 056803 (2012).[21] M. Veldhorst, M. Snelder, M. hoek, T. Gang, V. Guduru,X. L. Wang, U. Zeitler, W. G. van der Wiel, A. A. Gol-ubov, H. Hilgenkamp, et al., Nat. Mater. , 417 (2012).[22] L. P. Rokhinson, X. Liu, and J. K. Furdyna, Nat. Phys. , 795 (2012).[23] J. Linder, Y. Tanaka, T. Yokoyama, A. Sudbo, andN. Nagaosa, Phys. Rev. Lett. , 067001 (2010).[24] M. Sato and S. Fujimoto, Phys. Rev. Lett. , 217001(2010).[25] J. Alicea, Rep. Prog. Phys. , 076501 (2012).[26] T. D. Stanescu and S. Tewari, J. Phys.: Condens. Matter , 233201 (2013).[27] C. W. J. Beenakker, Annu. Rev. Con. Mat. Phys. , 113(2013).[28] F. Hassler, Quantum Information Processing. Lecturenotes of the 44th IFF Spring School 2013 (2013).[29] J. D. Sau, S. Tewari, R. M. Lutchyn, T. D. Stanescu, andS. D. Sarma, Phys. Rev. B , 214509 (2010). [30] J. Klinovaja and D. Loss, Phys. Rev. B , 085408(2012).[31] J. D. Sau, C. H. Lin, H.-Y. Hui, and S. D. Sarma, Phys.Rev. Lett. , 067001 (2012).[32] D. Chevallier, D. Sticlet, P. Simon, and C. Bena, Phys.Rev. B , 165414 (2013).[33] J. D. Sau and S. D. Sarma, Nat. Comm. , 964 (2012).[34] S. Nakosai, J. C. Budich, Y. Tanaka, B. Trauzettel, andN. Nagaosa, Phys. Rev. Lett. , 117002 (2013).[35] S. Mi, D. I. Pikulin, M. Wimmer, and C. W. J.Beenakker, Phys. Rev. B , 241405 (2013).[36] J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A.Fisher, Nat. Phys. , 412 (2011).[37] J. Sau, D. Clarke, and S. Tewari, Phys. Rev. B , 094505(2011).[38] Q. Liang, Z. Wang, and X. Hu, Euro Physics Letters ,50004 (2012).[39] P. Kotetes, G. Sch¨on, and A. Shnirman, Journal of theKorean Physical Society , 1558 (2013).[40] B. I. Halperin, Y. Oreg, A. Stern, G. Refael, J. Alicea,and F. von Oppen, Phys. Rev. B , 144501 (2012).[41] J. D. Sau, S. Tewari, and S. D. Sarma, Phys. Rev. A ,052322 (2010).[42] X. J. Liu, C. L. M. Wong, and K. T. Law,arXiv:1304.3765.[43] F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. , 056403 (2013).[44] J. Li, T. Neupert, B. A. Bernevig, and A. Yazdani,arXiv:1404.4058.[45] L. Weithifer, P. Recher, and T. L. Schimidt,arXiv:1309.4126.[46] T. Karzig, G. Refael, and F. von Oppen, Phys. Rev. X , 041017 (2013).[47] T. Hyart, B. van Heck, I. C. Fulga, M. Burrello, A. R.Akhmerov, and C. W. J. Beenakker, Phys. Rev. B ,035121 (2013).[48] C. V. Kraus, P. Zoller, and M. A. Baranov, Phys. Rev.Lett. , 203001 (2013).[49] X.-J. Liu and A. M. Lobos, Phys. Rev. B , 060504(2013).[50] C.-K. Chiu, M. Vazifeh, and M. Franz,arXiv:1403.0033v1.[51] M. Cheng, R. M. Lutchyn, V. Galitski, and S. D. Sarma,Phys. Rev. Lett. , 107001 (2009).[52] T. Mizushima and K. Machida, Phys. Rev. A , 023624(2010).[53] M. Cheng, R. M. Lutchyn, V. Galitski, and S. D. Sarma,Phys. Rev. B , 094504 (2010).[54] M. Cheng, V. Galitski, and S. D. Sarma, Phys. Rev. B , 104529 (2011).[55] J. von Neumann and E. P. Wigner, Z. Phys. , 467(1929).[56] H. Tal-Ezer and R. Kosloff, Jour. Chem. Phys. (1984).[57] D. G¨oddeke, R. Strzodka, and S. Turek, InternationalJournal of Parallel, Emergent and Distributed Systems(IJPEDS), Special issue: Applied parallel computing ,221 (2007).[58] M. Baboulin, A. Buttari, J. Dongarra, J. Kurzak, J. Lan-gou, J. Langou, P. Luszczek, and S. Tomov, Comp. Phys.Comm. , 2526 (2009).[59] S. L. Grand, A. W. G¨otz, and R. C. Walker, Comp. Phys.Comm. , 374 (2013).[60] For details, see Supplementary Material. Supplementary Material
Calculation Method
The Hamiltonian in Eq. (1) can be recast underBogoliubov-de Gennes representation as H BdG = − (cid:20) µ λ k x + cos k y ) (cid:21) τ z − ∆2 e iθ (sin k x τ y + sin k y τ x ) , (1)which in turn can be rewritten in real space in matrixform, defining our Hamiltonian as H = (cid:80) ij c † i H i , j c j onNambu basis, H i , i = (cid:18) − µ µ (cid:19) , H i ± ˆ e x , i = (cid:18) − λ/ ∓ ∆ e iθ / ± ∆ e − iθ / λ/ (cid:19) , H i ± ˆ e y , i = (cid:18) − λ/ ∓ ∆ e iθ / ∓ ∆ e − iθ / λ/ (cid:19) , (2)where i represents site position, assuming lattice con-stant to be 1. To numerically evaluate our system, wetake each wire length to be n sites long with one centralsite linking them. Each gate is represented as a factor g i ∈ [0 ,
1] ( i = 1 , , ,
4) multiplying the Hamiltonian el-ements between the central site and its neighbors in realspace, H ˆ e y , = g (cid:18) − λ/ − ie iθ ∆ / − ie − iθ ∆ / λ/ (cid:19) , H − ˆ e x , = g (cid:18) − λ/ e iθ ∆ / − e − iθ ∆ / λ/ (cid:19) , H − ˆ e y , = g (cid:18) − λ/ ie − iθ ∆ / ie − iθ ∆ / λ/ (cid:19) , H ˆ e x , = g (cid:18) − λ/ − e iθ ∆ / e − iθ ∆ / λ/ (cid:19) . (3)For a linear gate operation in the interval of time T asdescribed in this work, a gate being turned on (off) istaken to evolve as g i = t/T ( g i = 1 − t/T ), counting thetime t from the beginning of the operation. More generalfunctions may be used for operating both gates smoothlyat the same time. Finally, to achieve the wavefunctionchange in time, we expand the time-evolution operatorin terms of Chebishev polynomials[56], which can be re-trieved recursively, that is U ( t + ∆ t ; t ) = exp (cid:20) − i H ( t ) E ∆ tE (cid:21) = exp (cid:104) − i ˜ H ( t ) ∆ τ (cid:105) = ∞ (cid:88) k =0 c k ( ∆ τ ) T k ( ˜ H ( t )) , (4) E ≡ max |(cid:104) Ψ |H| Ψ (cid:105)| normalizes the Hamiltonian to avoidsingularities and c k ( ∆ τ ) = (cid:26) J ( ∆ τ ) ( k = 0)2( − i ) k J k ( ∆ τ ) ( k ≥
1) (5) T = 1 , T ( ˜ H ) = ˜ H,T k +1 ( ˜ H ) = 2 ˜ HT k ( ˜ H ) − T k − ( ˜ H ) , (6)constitute our expansion terms. J k are the Bessel func-tions of first kind, and their value is used to choose trun-cation point under double precision. While one can ap-ply the above operator successively, obtaining the desiredwavefunction in time t , this process is needlessly slow ifdone completely in double precision. Mixed precision[57–59] can be efficiently applied for a boost in calculationspeed if instead of calculating the wavefunction, only itsvariation is evaluated. Explicitly, this can be illustratedon the following steps:Ψ( t + ∆ t ) dl = U ( t + ∆ t ; t )Ψ( t ) dl → ∆ Ψ dl = ( U ( t + ∆ t ; t ) − t ) dl (7) ∆ Ψ = dl k (cid:48) (cid:88) k =0 c k T k ( ˜ H ) − fl Ψ( t ) fl = dl ( c − fl + k (cid:48) (cid:88) k =1 c k, fl T k ( ˜ H ) fl Ψ( t ) fl Ψ( t + ∆ t ) = dl Ψ( t ) dl + ∆ Ψ dl (8)Here the subscript dl (fl) represents a double (single) pre-cision conversion/variable. For example, equation 7 rep-resents a common time-evolution process done with dou-ble precision variables. The following expression claimsa conversion of Ψ( t ) from double to single precision, aswell as a single precision expansion of the time-evolutionoperator, while = dl implies that these data should be con-verted to double precision when adding to build up ∆ Ψ.Similarly, on the following line we point that the zero-order term should be changed to single precision aftercalculating it in double precision, and each other expan-sion term is calculated with single precision T k and c k .Shortly, each term in our expansion is calculated as avector in single precision, being later added as a doubleprecision variation to the double precision wavefunction.Note that the wavefunction is converted to single preci-sion for time evolution, but computed as double precisionin the end of each step. In fact, we start our evaluationwith an eigenstate of the Hamiltonian taken in doubleprecision, and only its variation, which corresponds tomost of the computations, is found in single precisionsteps.It is important to note that this method relies on theconstraint of small dt to work, as well as small andsmooth wavefunction variation, otherwise single preci-sion computation of the expansion terms may cumulate alarge error. Nevertheless, this constraint is also requiredfor the very numeric expansion of the time-evolution op-erator, in order to ignore its intrinsic time-ordering op-erator to a good approximation. In other words, thepossibility to evaluate time-evolution of the wavefunc-tion in real space-time already gives us the possibility ofa mixed precision method. Concretely, in our case each ∆ Ψ ∼ − , which in single precision allows for good enough numeric results in the 10 − ∼ − range, whichaccounts for first order terms, as well as 10 − ∼ − for the next order, all of them well fit in the limit ofdouble precision, up to 10 −15