Electron-Tunneling-Assisted Non-Abelian Braiding of Rotating Majorana Bound States

TL;DR

This paper proposes a novel method to demonstrate non-Abelian statistics of Majorana bound states by using electron tunneling to actively manipulate fermion parity during their rotation, revealing new interference effects.

Contribution

It introduces a theory where tunneling-induced parity changes enable non-Abelian braiding detection, challenging the view that parity fluctuations are purely harmful.

Findings

Resonant tunneling current reveals interference from non-commuting braiding operations.

Tunneling events change fermion parity, enabling active control of non-Abelian states.

The approach demonstrates a new way to utilize parity fluctuations for quantum information processing.

Abstract

It has been argued that fluctuations of fermion parity are harmful for the demonstration of non-Abelian anyonic statistics. Here, we demonstrate a striking exception in which such fluctuations are actively used. We present a theory of coherent electron transport from a tunneling tip into a Corbino geometry Josephson junction where four Majorana bound states (MBSs) rotate. While the MBSs rotate, electron tunneling happens from the tip to one of the MBSs thereby changing the fermion parity of the MBSs. The tunneling events in combination with the rotation allow us to identify a novel braiding operator that does not commute with the braiding cycles in the absence of tunneling, revealing the non-Abelian nature of MBSs. The time-averaged tunneling current exhibits resonances as a function of the tip voltage with a period that is a direct consequence of the interference between the non-commuting braiding operations. Our work opens up a possibility for utilizing parity non-conserving processes to control non-Abelian states.

Equations

$$AB|\Psi ⟩\ne BA|\Psi ⟩\text{or}(AB-BA)|\Psi ⟩\ne 0.$$

$$H_C$$

$${\mathcal{H}}_C$$

$$\displaystyle\Delta({\mathbf{r}})=\left\{\begin{array}[]{cc}\Delta_{0}e^{i\phi_{1}}&0\leq r<R,\\ \Delta_{0}e^{-i4\theta+i\phi_{2}}&r>R,\end{array}\right.$$

$$\begin{array}{r}{\gamma }_1\to -s{\gamma }_2,{\gamma }_2\to -s{\gamma }_3,\\ {\gamma }_3\to s{\gamma }_4,{\gamma }_4\to -s{\gamma }_1,\end{array}$$

$$M(t_q,t_{q^{\prime }})=-i\text{Tr}\left[{\rho }_0{\hat{\gamma }}_1(t_q){\hat{\gamma }}_1(t_{q^{\prime }})\right],$$

$${\bar{U}}_c={\gamma }_1{U}_c{\gamma }_1,$$

$$M(t_q,t_{q^{\prime }})=-i\text{Tr}\left[{\rho }_0^{\prime }{\left({\bar{U}}_c\right)}^n{\left(U_c^{†}\right)}^n\right],$$

$$H_T(t)$$

$$\frac{d\bar{I}}{dV}$$

$$T(\epsilon )$$

$$S(\epsilon )$$

$$\frac{d\bar{I}}{dV}=\frac{e^2}{h}\frac{\pi \hslash }{8T_J{k}_{B}T}\sum _{l}T({\epsilon }_l){\text{sech}}^2\left(\frac{eV-{\epsilon }_l}{2k_{B}T}\right),$$

$$\hat{S}(\epsilon )=\frac{1}{2}+\sum _{n=1}^Q{e}^{in\tilde{\epsilon }}{\left({\bar{U}}_c\right)}^n{\left(U_c^{†}\right)}^n.$$

$$\displaystyle\hat{S}(\varepsilon)\big{|}_{\text{even}}$$

$$\displaystyle\hskip 35.0pt-e^{i4\tilde{\varepsilon}}\big{)}\times\sum^{\infty}_{m=0}e^{im(4\tilde{\varepsilon}+\pi)},$$

$$⟨I(t)⟩=$$

$$=$$

$${\hat{H}}_T(t^{\prime })$$

$$\hat{I}(t)$$

$$V_{1k}(t_0)=\int d^{2}r{t}_k(\mathbf{r}){\Psi }_{M1\downarrow }(\mathbf{r},t_0),$$

$$\begin{array}{rl}{\Gamma }_{kq{q}^{\prime }}(t,t^{\prime })& =|V_{1k}(t_0){|}^2{e}^{-\lambda |t-t_q|}{e}^{-\lambda |t^{\prime }-t_{q^{\prime }}|},\\ G_k(t,t^{\prime })& =-i⟨{\hat{c}}_k(t){\hat{c}}_k^{†}(t^{\prime })⟩=-i{e}^{-i({\epsilon }_k+eV)(t-t^{\prime })/\hslash }\left[1-n_F({\epsilon }_k)\right],\\ {\bar{G}}_k(t,t^{\prime })& =-i⟨{\hat{c}}_k^{†}(t){\hat{c}}_k(t^{\prime })⟩=-i{e}^{i({\epsilon }_k+eV)(t-t^{\prime })/\hslash }{n}_F({\epsilon }_k),\end{array}$$

$$M(t,t^{\prime })\approx M(t_q,t_{q^{\prime }})=-i\text{Tr}\left[{\rho }_0{\hat{\gamma }}_1(t_q){\hat{\gamma }}_1(t_{q^{\prime }})\right],$$

$$M(t_q,t_{q^{\prime }})=-i⟨{\Psi }_g(t_0)|{\hat{\gamma }}_1(t_q){\hat{\gamma }}_1(t_{q^{\prime }})|{\Psi }_g(t_0)⟩.$$

$$⟨I(t)⟩-⟨I(t-T_J)⟩\sim \sum _k{e}^{-i({\epsilon }_k+eV)(t-t_0)}{n}_F({\epsilon }_k)+\text{c.c.}\sim 0,$$

$$\bar{I}=\frac{1}{T_J}{\int }_{\tilde{t}-T_J}^{\tilde{t}}dt⟨I(t)⟩.$$

$$\bar{I}$$

$$=\frac{e}{h}{\int }_{-\infty }^{\infty }d\epsilon T(\epsilon )S(\epsilon )n_F(\epsilon -eV)+\frac{e}{h}{\int }_{-\infty }^{\infty }d\epsilon T(\epsilon )S(\epsilon )n_F(-\epsilon -eV)-\frac{e}{h}{\int }_{-\infty }^{\infty }d\epsilon T(\epsilon )S(\epsilon )$$

$$=\frac{e}{h}{\int }_{-\infty }^{\infty }d\epsilon T(\epsilon )[S(\epsilon )+S(-\epsilon )]n_F(\epsilon -eV)-\frac{e}{h}{\int }_{-\infty }^{\infty }d\epsilon T(\epsilon )S(\epsilon ),$$

$$T(\epsilon )$$

$$S(\epsilon )$$

$$S(\epsilon )+S(-\epsilon )$$

$$\displaystyle\left\{\hat{\gamma}_{1}(t_{q}),\hat{\gamma}_{1}(t_{q^{\prime}})\right\}=\left\{\begin{array}[]{cc}2e^{im\alpha}&\text{for}~{}q-q^{\prime}=4m,\\ 0&\text{otherwise}.\end{array}\right.$$

$$\frac{d\bar{I}}{dV}=\frac{e^2}{h}\frac{\pi \hslash }{8T_J{k}_{B}T}\sum _{l}T({\epsilon }_l){\text{sech}}^2\left(\frac{eV-{\epsilon }_l}{2k_{B}T}\right),$$

Full text

Introduction.— A braiding operation reveals the quantum statistics of identical particles Nayak2008 ; Wilczek2009 ; Stern2010 . Majorana zero-energy states bound to certain defects (e.g. vortices or edges) in topological superconductors are quasiparticles obeying non-Abelian statistics Hasan2010 ; Alicea2012 ; Beenakker2013 ; DasSarma2015 ; Aguado2017 . In an isolated system with $2N$ decoupled Majorana states, there is a $2^{N}$-fold degenerate ground state manifold $\{|\Psi\rangle\}$, and adiabatically moving one Majorana state around another acts as a unitary matrix on the manifold. Such unitary matrices of different braiding operations, $A$ and $B$, are in general non-commutative, so that the order of operations matter,

[TABLE]

Non-Abelian braiding is one of the hallmarks of topological quantum phases associated with non-Abelian statistics appearing in many contexts Stern2010 ; Zhu2011 ; Thomas2016 and also represents the basic resource for executing topologically protected gates for quantum computing Nayak2008 ; Aasen2016 .

The essence of the present work is to provide transport signatures of Majorana bound states (MBSs) induced by the non-commutativity shown in Eq. (1). The envisioned system is a Corbino geometry topological Josephson junction (JJ), formed by two $s$-wave superconductors on a topological insulator (TI) surface [see Fig. 1(a)]. Four vortices, each hosting a MBS, rotate along the junction, and the time-dependent tunneling conductance between the junction and a metallic tip is measured Park2015 . A ground state of the system evolves in the fourfold degenerate ground state manifold, governed by the rotation and the coherent electron tunneling processes. The evolution can be cast into two braiding operators (corresponding to $A$ and $B$ in Eq. (1)) which do not commute: one is a parity-conserving rotation and the other is a tunneling-assisted braiding. In the low bias voltage regime, the time-averaged conductance exhibits unusual peak positions, which we interpret as a direct signature of non-commutativity of the two braiding operators.

Tremendous amounts of proposals and experiments lead to great achievements in the realization Fu2008 ; Lutchyn2010 ; Oreg2010 ; Alicea2010 ; Choy2011 ; Perge2013 ; Hell2017 , manipulation Alicea2011 ; Flensberg2011 ; vanHeck2012 ; Grosfeld2011 ; Mi2013 ; Li2016-1 and detection Mourik2012 ; Das2012 ; Deng2012 ; Lee2014 ; Xu2015 ; Deng2016 ; Pawlak2016 ; Feldman2017 ; Gul2018 ; Rokhinson2012 ; Wiedenmann2016 ; Deacon2017 of MBSs in superconducting hybrid structures. In particular, a recent experiment exploiting a quantum anomalous Hall insulator-superconductor structure He2017 boosts interest in searches for transport signatures of non-Abelian braiding Lian2017 ; Beenakker2018 . Based on such hybrid structures, the authors of Refs. Lian2017 ; Beenakker2018 theoretically investigated transport properties of Mach-Zehnder-like interferometers of chiral Majorana modes. The overlap or fusion of two paths of Majorana modes whose relative dynamics is determined by braiding with the other Majoranas signals a unitary evolution (which is not a phase factor) of Majorana modes.

Different to these recent studies in Refs. Lian2017 ; Beenakker2018 , we demonstrate interference involving four rotating MBSs whose braiding operations are assisted by tunneling of electrons into or out of the MBSs and thus in which the fermion parity formed by the MBSs is not conserved. Such tunneling-assisted braiding has been to the best of our knowledge not considered before, on the contrary, electron tunneling was seen detrimental for topological quantum processing Budich2012 ; Leijnse2012 ; Woerkom2015 . We will show that, in our scheme, electron tunneling probes non-Abelian statistics via the tunneling conductance. Our scheme does not require control of fusions of Majorana states.

Theoretical model.— We consider a Corbino JJ deposited on the surface ($x$-$y$ plane) of a three dimensional TI [Fig. 1(a)]. The circular shaped junction with a radius $R$ is formed by thin films of inner ($S_{1}$) and outer ($S_{2}$) s-wave superconductors and contains four magnetic flux quanta, $4\Phi_{0}$ with $\Phi_{0}=h/(2e)$, inducing a phase difference across the junction (see Eq. (8)). The Bogoliubov-de Gennes (BdG) Hamiltonian for the TI surface proximity coupled to the Corbino JJ is given by app1

[TABLE]

and $\Phi({\mathbf{r}})=(\Phi_{\uparrow},\Phi_{\downarrow},\Phi^{\dagger}_{\downarrow},-\Phi^{\dagger}_{\uparrow})^{T}$ is the Nambu spinor and $\mathcal{H}_{0}=v_{F}(\sigma_{x}p_{x}+\sigma_{y}p_{y})$ with Pauli spin matrices $\sigma_{x,y}$ describes the surface states and $\mu$ is the chemical potential. The proximity-induced superconducting gap $\Delta({\mathbf{r}})$ is

[TABLE]

where $\phi_{1}$ and $\phi_{2}$ are spatially uniform phases in each superconducting region, and the polar-angle-dependent phase $-4\theta$ at $r>R$ is due to the presence of the four flux quanta Clem2010 . By solving the BdG equation $\mathcal{H}_{C}\Psi({\mathbf{r}})=E\Psi({\mathbf{r}})$, we find four Majorana wave functions $\Psi_{\text{M}j}({\mathbf{r}})$ with $j\in\{1,2,3,4\}$, at zero energy $E=0$. They are localized at ($r,\theta$) = ($R,\theta_{j}$) where $\theta_{j}=(3\pi-2\pi j)/4-(\phi_{1}-\phi_{2})/4$, at which the local phase difference across the junction is $\pi$. Detailed calculations of the Majorana wave functions for $\mu=0$ are given in Supplemental material Suppl .

If we change $\phi_{1}-\phi_{2}$ by $2\pi$, the four MBSs rotate by $\pi/2$ in a clockwise direction maintaining their relative distances, as plotted in Fig. 1(a), leading to a transformation $\gamma_{j}\rightarrow U_{c}\gamma_{j}U^{\dagger}_{c}$,

[TABLE]

where $\gamma_{j}=\int d^{2}r~{}\Psi^{\dagger}_{\text{M}j}({\mathbf{r}})\Phi({\mathbf{r}})$. $s=1(-1)$ corresponds to the change of $\phi_{1}(\phi_{2})$ by $2\pi(-2\pi)$. Graphical representation of the transformation is given in Fig. 1(b) for the $s=-1$ case. A rotation operator $U_{c}$ for the transformation can be constructed as a product of three pairwise braidings $U_{c}=U_{41}U_{12}U_{23}$ where $U_{ij}$ is the braiding exchange operator of $\gamma_{i}$ and $\gamma_{j}$ given by $U_{ij}=\text{exp}\left(s\pi\gamma_{i}\gamma_{j}/4\right)$ Ivanov2001 .

The adiabatic rotation can be achieved if a dc-bias voltage $V_{J}$ across the junction is much smaller than the excitation energy of the junction. For a finite $V_{J}$, $\phi_{1}-\phi_{2}$ varies in time $t$ as $\phi_{1}-\phi_{2}=\phi_{0}+2eV_{J}t/\hbar$ where $\phi_{0}$ is a spontaneously chosen constant. The states $\Psi_{\text{M}j}({\mathbf{r}},\phi_{1}(t),\phi_{2}(t))$ then become instantaneous eigenstates of $\mathcal{H}_{C}[\phi_{1}(t),\phi_{2}(t)]$ at zero energy, and $U_{c}$ can be considered as the time evolution operator of the MBSs from $t$ to $t+T_{J}$, where $T_{J}=\frac{\pi\hbar}{eV_{J}}$ is the time needed for the $\pi/2$-rotation.

Tunneling-assisted Majorana braiding.— To explore the effect of electron tunneling, we connect a metal tip to the Corbino JJ, as depicted in Fig. 2(a). The tip is located such that an electron can tunnel onto or off the Corbino JJ through $\gamma_{1}(t_{0})$ at $t=t_{0}$, and we assume that the tunnel coupling is switched on at $t=t_{0}$. A phase coherent time-dependent tunneling event between the tip and adiabatically rotating Majorana states can occur at discrete times $t_{q}=t_{0}+qT_{J}$, where $q=0,1,2,...$. Creation or annihilation of an electron via a Majorana state at $t=t_{q}$ is described by $\gamma_{1}(t_{0})|\Psi_{g}(t_{q})\rangle$ where $|\Psi_{g}(t_{q})\rangle=U_{c}^{q}|\Psi_{g}(t_{0})\rangle$ is the time-evolved initial state (being part of the ground-state manifold) of the MBSs from $t_{0}$ to $t_{q}$. Note that our proposal does not depend on the initial configuration of the ground state and other choices of Majorana states coupled to the tip at $t=t_{0}$. Hereafter, we will denote $\gamma_{1}(t_{0})$ by $\gamma_{1}$.

The time evolution of a Majorana state from $t=t_{q^{\prime}}$ to $t_{q}$ at which tunneling events occur is described by the Majorana Green’s function

[TABLE]

where $\hat{\gamma}_{1}(t_{q})=\left(U^{\dagger}_{c}\right)^{q}\gamma_{1}U_{c}^{q}$ and $\rho_{0}$ is a density matrix of the Majorana state at $t=t_{0}$. For a more comprehensive description of the tunneling effect, we introduce a tunneling-assisted braiding operator,

[TABLE]

consisting of three events: changing fermion-occupation-number parity due to the tunneling at $t=t_{q}$, followed by an evolution for a time $T_{J}$ with $U_{c}$, and then changing the parity again at $t=t_{q}+T_{J}$. The transformation governed by $\bar{U}_{c}$ is drawn in Fig. 2(b); comparing the cases without and with the tunneling in Figs. 1(b) and 2(b), respectively, notice that the tunneling effectively reverses the direction of the pairwise braiding when a braiding involves $\gamma_{1}$. Therefore, ${\bar{U}}_{c}=U_{14}U_{21}U_{23}$ can be considered – besides $U_{c}$ – as another genuine braiding operator. $M(t_{q},t_{q^{\prime}})$ then can be presented as

[TABLE]

where we used the cyclic property of the trace. $\rho^{\prime}_{0}=\left(U_{c}\right)^{q}\rho_{0}\left(U^{\dagger}_{c}\right)^{q}$ and $n=q-q^{\prime}$. We find that $U^{\dagger}_{c}$ and $\bar{U}_{c}$ do not commute, $\left[\bar{U}_{c},U^{\dagger}_{c}\right]\neq 0$. As a consequence, $M(t_{q},t_{q\prime})$ is not just a sum of phase factors but involves non-trivial state changes in the ground-state manifold. We show below that the non-commuting braidings result in observable interference signatures free of the necessity of physically fusing MBSs.

Transport signatures.— To obtain the tunneling current between the tip and the JJ in the weak coupling limit, we extend the formalism of Ref. Park2015 to four MBSs. The Hamiltonian of the tip is $H_{N}=\sum_{k\sigma}\varepsilon_{k}c^{\dagger}_{k\sigma}c_{k\sigma}$ where $c_{k\sigma}$ is the electron annihilation operator in the tip with momentum $k$ and spin $\sigma$. Since we are interested in the low-energy sector of the junction, tunneling between the tip and the MBSs is the only relevant process. Around $t=t_{q}$ where the coupling strength to $\gamma_{1}$ is maximal, we assume that the coupling increases and decreases exponentially as $\gamma_{1}$ approaches to and leaves from the tip, respectively, while its phase does not change significantly. Moreover, since the Majorana states are spin polarized, and couple only to electrons of the tip with their spin parallel to that of the Majorana states; electrons with opposite spin are reflected at the junction between the tip and the Corbino JJ and do not contribute to the tunneling current. Then the tunneling Hamiltonian becomes

[TABLE]

where $\lambda^{-1}$ is the tunneling duration and $V_{1k}$ is the coupling between the tip and $\gamma_{1}$. Here we have assumed $\lambda^{-1}\ll T_{J}$, implying that only nearest-neighbor coupling between the tip and the MBSs is taken into account.

Using the current expression $I(t)=-e\,dN_{T}/dt$ with the tip number operator $N_{T}=\sum_{k}c^{\dagger}_{k}c_{k}$ and lowest order perturbation theory in $H_{T}(t)$, the differential conductance of the time-averaged current measured after many rotation cycles of MBSs has the form,

[TABLE]

where $n_{F}$ is the Fermi-Dirac distribution and $eV$ is the bias voltage. The tunneling probability $T(\varepsilon)$ and the interference term $S(\varepsilon)$ are given by

[TABLE]

Here $\tilde{\varepsilon}=\varepsilon/(\hbar T^{-1}_{J})$ and the integer $Q\gg 1$, which will go to infinity later. $\Gamma=2\pi\rho|V_{1k}|^{2}$ where $\rho$ is the tip density of states. We assumed a wide-band approximation where $\rho$ and $V_{1k}$ are energy independent and we neglected the contributions proportional to $e^{-\lambda T_{J}/2}$; note that these small contributions do not change the positions of conductance peaks. The details for the calculation of $\bar{I}$ are given in Suppl . In the limit $Q\rightarrow\infty$, we obtain

[TABLE]

which shows peaks at $\varepsilon_{l}=\frac{\hbar}{4T_{J}}(2\pi l-\alpha)$ where $l$ is an integer and $\alpha=\pi$ arising from a $2\pi$-rotation of the four MBSs. This perturbative calculation is valid for $T(\varepsilon_{0})\hbar/(8T_{J})\ll k_{B}T\ll E_{g}$ where $E_{g}$ is the excitation energy of the junction.

The $d\bar{I}/dV$ in Eq. (17) is plotted in Fig. 3 for realistic parameters. It shows peaks at $eV=\varepsilon_{l}$. This is our main result. The peak positions are determined by $T_{J}$ and $\alpha$, but are independent of system details such as the initial Majorana state at $t=t_{0}$ and the tunneling strength $\Gamma$. Note that the periodicity $T_{J}$ of the system Hamiltonian in Eqs. (2) and (13) does not coincide with the periodicity of the ground state $4T_{J}$ from the fact that $U_{c}^{4}={\bar{U}}_{c}^{4}=\mathbbm{1}$. It is a consequence of the non-trivial state evolution within the ground-state manifold of 4MBSs requiring a matrix structure. As shown below, the $4T_{J}$-periodicity and the non-commutativity between $U_{c}$ and $\bar{U}_{c}$ result in peaks in $d{\bar{I}}/dV$ separated by $h/(4T_{J})$ and not by $h/T_{J}$ associated with the frequency of appearance of MBSs beneath the tip. The results are the same for the case of an anticlockwise rotation of four MBSs.

Non-Abelian statistics.— In order to clearly unveil such a link between the interference effect and the non-Abelian matrix structure, we analyse the term $S(\varepsilon)$ in the occupation number basis $\{|n_{1}n_{2}\rangle\}$, where $n_{1},n_{2}=0,1$ are occupation numbers for fermionic operators $f_{1}=\frac{\gamma_{1}+i\gamma_{2}}{2}$ and $f_{2}=\frac{\gamma_{3}+i\gamma_{4}}{2}$, see Suppl for more details on the occupation number representation. As Eq. (17) is independent of the initial condition, the specific form of the initial density matrix ($\rho_{0}$ or $\rho^{\prime}_{0}$) is unimportant. Substituting Eq. (12) into Eq. (16) leads to $S(\varepsilon)=\text{Re}\{\text{Tr}[\rho^{\prime}_{0}\hat{S}(\varepsilon)]\}$ where

[TABLE]

Note that the operations $U_{c}^{\dagger}$ and $\bar{U}_{c}$ do not commute, and thus the sum cannot be treated as a simple geometric series: $\sum^{\infty}_{n=1}e^{in\tilde{\varepsilon}}\left(\bar{U}_{c}\right)^{n}\left(U^{\dagger}_{c}\right)^{n}\nrightarrow\sum^{\infty}_{n=1}e^{in\tilde{\varepsilon}}e^{in\varphi}$. The operator $\hat{S}(\varepsilon)$ comes from the overlap between the following two processes of temporal length $QT_{J}$: In process I, an electron tunnels from the tip to $\gamma_{1}$ at $t_{0}+(Q-n)T_{J}$, and in process II, the tunneling happens at $t_{0}+QT_{J}$. Here $e^{in\tilde{\epsilon}}$ is the dynamical phase factor gained for the time interval $nT_{J}$. The interference between terms of different $n$ determines the peak positions of the conductance.

Let us assume that an even parity state, a mixture of $|00\rangle$ and $|11\rangle$, is prepared at $t=t_{0}$; the case of an odd parity state is obtained in a similar way. In the limit $Q\rightarrow\infty$, Eq. (18) for an even parity is given by

[TABLE]

where $\tau_{x,y,z}$ are Pauli matrices acting in the space of the even parity states, $|00\rangle$ and $|11\rangle$. In the second line, the summation is classified into four categories in each of which the Pauli matrix (including the identity matrix) is factored out, manifesting the interference with period $4T_{J}$. Using Eqs. (18) and (19) yields $S(\varepsilon)+S(-\varepsilon)=\sum_{m}\exp(im(4{\tilde{\varepsilon}}+\pi))\sim\sum_{l}\delta(4{\tilde{\varepsilon}}+\pi(2l+1))$, where $m,l$ are integers. Together with Eq. (14) we obtain our final result Eq. (17). We note that the period of $4T_{J}$ cannot be obtained by corresponding braiding operators that would commute, see Suppl . We also note that this non-Abelian interference effect cannot be envisaged in a system with two MBSs where non-commuting braiding operations do not occur Park2015 .

We remark that the suggested test of non-Abelian braiding statistics needs only a local measurement of MBSs that are at zero energy so that the way we fuse the 4 MBSs into the two fermions $f_{1}$ and $f_{2}$ is actually arbitrary. The period $4T_{J}$ also does not depend on a specific initial state (if the time-average is performed after times $t\gg T_{J}$) but is only a consequence of the non-commuting matrix structure of $U_{c}$ and ${\bar{U}}_{c}$. The extracted information of the state changes is due to interference that is generated because the MBSs rotate in the Corbino geometry JJ. This is fundamentally different compared to other braiding schemes which use the selective switching on and off of couplings between the Majorana bound states and the read-out of the non-Abelian state changes is done without physically moving the MBSs vanHeck2012 ; Bonderson2013 . In our scheme the rotation induces a dynamical coupling between the MBSs as we discuss in detail in the Supplemental material employing the Floquet picture. There we consider also the zero temperature case to all orders in the tunneling from the tip to the MBSs.

Discussion and conclusion.— We have demonstrated that a non-Abelian state evolution can be identified in tunneling conductance measurements between four rotating MBSs in a Corbino geometry topological Josephson junction and a metal tip. Unitary evolutions of the MBSs acting on even and odd parity subspaces, which are separable if the fermion parity is conserved, are intertwined by electron tunneling, inducing parity-conserving and tunneling-assisted braiding operators. Coherent interference between different orders of round trips of Majorana states governed by the parity-conserving and tunneling-assisted braiding operators yields a time-averaged conductance exhibiting peaks with a period of $h/(4T_{J})$ as a function of bias voltage between the metal tip and the Josephson junction, whereas the period of the Hamiltonian is $T_{J}$. This constitutes a clear signature of non-Abelian state evolution of four MBSs.

We explicitly showed that these results have their origin in the non-commutativity of the parity-conserving and tunneling-assisted braiding operators and are therefore independent on the way we fuse the MBSs into fermions which is fundamentally different from other recent proposals that use time-dependent couplings between the MBSs or Coulomb interaction to lift their degeneracies vanHeck2012 ; Bonderson2013 ; Aasen2016 . Here, an effective coupling between MBSs is induced dynamically by the rotation which only requires a dc-Josephson voltage applied between the two superconductors.

We expect that other kinds of exotic zero modes such as MBSs in time-reversal invariant topological superconductors Zhang2013 ; Keselman2013 ; Haim2014 ; Wolms2015 ; Wolms2016 ; Li2016-2 ; Schrade2018 and parafermions Fendley2012 ; Lindner2012 ; Cheng2012 ; Clarke2013 ; Vaezi2013 ; Barkeshli2014 ; Jelena2014 ; Maghrebi2015 ; Alicea2016 could be analyzed with our time-dependent tunneling scheme to manifest the quantum statistics of the corresponding modes.

The experimental realization may be challenging, but within reach of current experiments. Assuming the proximity-induced superconducting gap $\Delta_{0}=$1 meV that can be achieved, for example, in thin-films of Nb or NbN Lin2013 ; Du2017 , the excitation energy gap of Josephson vortices of the junction can be estimated by $E_{g}=\Delta_{0}\sqrt{4\xi/R}\sim 0.9$ meV for the radius of the junction $R=5\xi$ Park2015 ; Potter2013 , where $\xi$ is the superconducting coherence length. We require a coherent and adiabatic rotation of the MBSs so that $T_{J}$ (the time taken for the $\pi/2$ rotation) should satisfy $\hbar/E_{g}(=0.7~{}\text{ps})\ll T_{J}\ll t_{\text{qp}}(\gtrsim\mu\text{s})$ where $t_{\text{qp}}$ is the quasiparticle poisoning time Rainis2012 ; Higginbotham2015 . At the same time, the temperature should be much smaller than the separation between the conductance peaks $h/(4T_{J})$. MBSs can be spaced unequally apart in the presence of inhomogeneities in the junction. However, they do not affect the rotation time $T_{J}$ due to the periodicity of the system Hamiltonian and corresponding interference traces on the time scale of $4T_{J}$ due to non-Abelian evolution would remain. We believe that the Corbino geometry topological Josephson junction can also be realized in heterostructures of a thin-film topological insulator and a superconductor Hao2017 or Pb/Co/Si(111) two-dimensional topological superconductor Cren2018 .

Our findings provide a new way of looking at braiding experiments, by actively using parity switching events by tunneling, instead of avoiding them. This may define a new way to build non-Abelian operations for topological qubits utilizing coherent fluctuations of the fermion parity. Such a change of fermion parity could be achieved on demand during a definite time using charge pumps based on quantum dots in the single electron regime Fricke2014 coupled to the setup. Quantum dots could already be coupled to MBSs in experiment Deng2016 .

A. Time-averaged tunneling current

The time-dependent tunneling current between a metal tip and a Corbino geometry topological Josephson junction in the weak tunneling limit can be obtained using lowest order perturbation theory. To lowest order in $H_{T}(t)$, we find the tunneling current $\langle I(t)\rangle=-e\langle dN_{T}(t)/dt\rangle$,

[TABLE]

where $N_{T}(t)=\displaystyle\sum_{k}c^{\dagger}_{k}(t)c_{k}(t)$ is the metal tip number operator. $\hat{H}_{T}(t^{\prime})$ and $\hat{I}(t)$ which are expressed in the interaction picture are given by

[TABLE]

The tunneling Hamiltonian $\hat{H}_{T}(t^{\prime})$ switched on at time $t_{0}$ is valid in the low energy regime where MBSs are the only relevant states for the tunneling current and for $\lambda^{-1}\ll T_{J}$. The coupling coefficient $V_{1k}(t_{0})$ between the tip and $\hat{\gamma}_{1}(t_{0})$ is

[TABLE]

where $t_{k}({\mathbf{r}})$ is the tunneling coefficient between the tip and the junction and $\Psi_{M1\downarrow}({\mathbf{r}},t_{0})$ is the electron spin-down component of the Majorana wave function $\Psi_{M1}({\mathbf{r}})$ in Eq. (S79). In Eq. (S1), the time-dependent tunneling parameter $\Gamma_{kqq^{\prime}}(t,t^{\prime})$ and the tip-electron Green’s functions $G_{k}(t,t^{\prime})$ and $\bar{G}_{k}(t,t^{\prime})$ are given by

[TABLE]

where $\langle\cdot\rangle$ is the expectation value over a thermal ensemble of initial states at $t=t_{0}$, and $n_{F}(\varepsilon_{k})=1/[1+e^{\varepsilon_{k}/(k_{B}T)}]$ is the Fermi-Dirac distribution at $t=t_{0}$ with the temperature $T$. Since the tunneling current is exponentially small except for $t=t_{q}$ and $t^{\prime}=t_{q^{\prime}}$ due to the presence of the exponential factor of $\Gamma_{kqq^{\prime}}(t,t^{\prime})$, we can approximate the Majorana Green’s function,

[TABLE]

where $\rho_{0}$ is a density matrix of the Majorana state at $t=t_{0}$. If the Josephson junction is in one of the ground states $|\Psi_{g}(t_{0})\rangle$ at $t=t_{0}$, the density matrix has the form of $\rho_{0}=|\Psi_{g}(t_{0})\rangle\langle\Psi_{g}(t_{0})|$ and we get

[TABLE]

If $t$ is very far from $t_{0}$, we can find that the difference between $\langle I(t)\rangle$ and $\langle I(t-T_{J})\rangle$ is negligible,

[TABLE]

yielding a time-periodic behavior of the tunneling current $\langle I(t)\rangle=\langle I(t-T_{J})\rangle$. Without loss of generality, we assume that $t$ is in the interval $\left[(Q-1/2)T_{J},(Q+1/2)T_{J}\right]$ where $Q$ is a very large integer, $Q\gg 1$. Then the time-averaged tunneling current over an interval $\left[\tilde{t}-T_{J},\tilde{t}\right]$ is

[TABLE]

Let us change the variable in Eq. (S1) from $\varepsilon_{k}$ to $\varepsilon_{k}-eV$. After some algebra, we find $\bar{I}$ as

[TABLE]

where $T(\varepsilon)$ and $S(\varepsilon)$ are

[TABLE]

The second term in the third line in Eq.(S10) can be disregarded because it is independent of the bias voltage and does not contribute to the tunneling conductance. The term $S(\varepsilon)+S(-\varepsilon)$ in Eq. (S10) is written in terms of the Majorana Green’s function, and contains information of the non-commuting braiding operations. It yields

[TABLE]

The notation $[Q/4]$ denotes the integer part of the number $Q/4$ and we have used anti-commutation relations

[TABLE]

The phase factor $e^{im\alpha}$ with $\alpha=\pi$ comes from a $2\pi m$-rotation of the four MBSs, $\left(U^{\dagger}_{c}\right)^{4m}\gamma_{j}\left(U_{c}\right)^{4m}=(-1)^{m}\gamma_{j}$, and is physically due to crossing branch cuts emanating from the MBSs. In the limit $Q\rightarrow\infty$ (or $\tilde{t}-t_{0}\rightarrow\infty$), we obtain

[TABLE]

where $\varepsilon_{l}=\frac{\hbar}{4T_{J}}(2\pi l-\alpha)$ with integer $l$.

B. Occupation number representation

We describe the rotation of the four MBSs in occupation number space. We define two complex fermion operators,

[TABLE]

and four occupation number states which are degenerate at zero energy,

[TABLE]

Here, the state $|00\rangle$ is defined by $f_{1}|00\rangle=f_{2}|00\rangle=0$. The two states in each fermion-occupation-number parity subspace form a qubit, $|00\rangle$ and $|11\rangle$ for the even and $|10\rangle$ and $|01\rangle$ for the odd fermion parity subspace. In the basis $\{|00\rangle,|11\rangle,|10\rangle,|01\rangle\}$, $U_{c}$ in the main text is represented as

[TABLE]

where $U_{ce}(U_{co})$ is the evolution operator acting on the even (odd) parity space that rotates the qubit by $\pi/2(\pi)$ about the direction of $\hat{n}_{ce}(\hat{n}_{co})$ given by

[TABLE]

$\vec{\tau}=(\tau_{x},\tau_{y},\tau_{z})$ are Pauli matrices acting on the qubit, and $\bf{0}$ is $2\times 2$ null matrix. The qubit rotations induced by $U_{ce}$ and $U_{co}$ on the Bloch sphere are illustrated in Fig. S1.

The parity of the fermion occupation number is conserved in the transformation $U_{c}$. $\bar{U}_{c}=\gamma_{1}U_{c}\gamma_{1}$ in the same basis is represented by the interchange of $U_{ce}$ and $U_{co}$,

[TABLE]

We find that $U^{\dagger}_{c}$ and $\bar{U}_{c}$ (or $U_{c}$ and $\bar{U}_{c}$) do not commute,

[TABLE]

This is indicative of the different braiding evolutions of world lines corresponding to the two operator products $\bar{U}_{c}U^{\dagger}_{c}$ and $U^{\dagger}_{c}\bar{U}_{c}$.

C. Majorana wave functions

We provide the details of the calculation of Majorana wave functions $\Psi_{\text{M}j}({\mathbf{r}})$ with $j\in\{1,2,3,4\}$ in a Corbino geometry topological Josephson junction. We solve the BdG equation $\mathcal{H}_{C}\Psi({\mathbf{r}},\phi_{1},\phi_{2})=E\Psi({\mathbf{r}},\phi_{1},\phi_{2})$ for $E=0$ and $\mu=0$. Hereafter, we use the dimensionless length scale $r$ normalized by $\xi=\hbar v_{F}/\Delta_{0}$. For $r<R$, the wave function $\Psi_{r<R}(r,\theta)$ is given by

[TABLE]

where $I_{m}(r)$ is the modified Bessel function of the first kind, and $a_{m}$ and $b_{m}$ are coefficients. The wave function $\Psi_{r>R}(r,\theta)$ at $r>R$ is given by

[TABLE]

where $K_{n}(r)$ is the modified Bessel function of the second kind, and $c_{n}$ and $d_{n}$ are coefficients. We consider only the wave functions with spin down as those for spin up become non-normalizable solutions, and hence the coefficients $a_{m}$ and $c_{n}$ should be zero for all $m$ and $n$. In order to get the coefficients $b_{m}$ and $d_{n}$ we match the spin-down components at $r=R$,

[TABLE]

leading to

[TABLE]

and the following recurrence relations,

[TABLE]

where $l$ and $l^{\prime}$ are integers. From these recurrence relations we can construct four linearly independent solutions,

[TABLE]

where ${\eta}\in\{-2,-1,0,1\}$ and $\Psi^{\prime}_{\eta}(r,\theta)$ are

[TABLE]

Here the wave functions $\Psi^{<}_{4m+{\eta}}(r,\theta)$ at $r<R$ and $(\Psi^{>}_{4n+{\eta}}(r,\theta))$ at $r>R)$ are given by

[TABLE]

and coefficients $B_{m{\eta}}$ and $D_{n{\eta}}$ are

[TABLE]

By superposing the solutions $\Psi^{\prime}_{\eta}(r,\theta)$ in Eq. (S56) and using particle-hole symmetry, we find four Majorana states $\Psi_{\text{M}j}$ satisfying $\Xi\Psi_{\text{M}j}=\Psi_{\text{M}j}$ where $\Xi=\sigma_{y}\tau_{y}\mathcal{C}$ is the particle-hole operator and $\mathcal{C}$ is the operator for complex conjugation, see Fig. S2. They are given by

[TABLE]

where the azimuthal angles $\theta=\theta_{j}$ at which $\Psi_{\text{M}j}$ are localized are given by $\theta_{j}=(3\pi-2\pi j)/4-(\phi_{1}-\phi_{2})/4$ and $N_{\eta}$ are normalization constants such that

[TABLE]

D. Floquet analysis

We confirm the conductance peak positions in the main text by using a Floquet analysis for four rotating MBSs including all orders in tunneling at zero temperature. This Floquet description is applicable as the time-dependent Hamiltonian discussed in the main text is periodic in time with periodicity $T_{J}$. See Ref. app-Cayssol2013 for a short review on the Floquet formalism applied to topological insulators and Ref. app-Park2015 in which two rotating MBSs are analyzed by using the Floquet formalism.

The Floquet Hamiltonian $H_{F}$ is defined by the time-independent Hamiltonian that would yield the same evolution as with $U_{c}$ in the main text after one period $T_{J}$,

[TABLE]

where $\hat{n}_{ce}$ and $\hat{n}_{co}$ are given in Eq. (S21), and $l$ is an integer. The last term $-2\pi l\hbar/T_{J}\mathbb{I}$ in $H_{F}$ only shifts the energy levels and can be ignored for the moment. At the end of the calculation, we will restore this term. The representation of $H_{F}$ in terms of Majorana operators is

[TABLE]

where $E_{0}=\pi\hbar/(4T_{J})$ and the subindices $i$ and $j$ range from $1$ to $4$. $t_{ij}$ describing the effective coupling between MBSs caused by the rotation is the $(i,j)$ component of the matrix ${\bf t}$ in the basis $(\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4})$,

[TABLE]

We calculate the differential conductance of the metal tip coupled to the Majorana network described by $H_{F}$ following the Keldysh technique calculation used in Ref. app-Flensberg2010 . We consider the case where the tip is coupled only to $\gamma_{1}$. The differential conductance then is given by the formula

[TABLE]

where $G^{R}_{11}(\omega)$ is the $(1,1)$ component of the $4\times 4$ matrix ${\bf G}^{R}(\omega)$ given by

[TABLE]

where ${\bf\Gamma}$ describes the tip-MBS tunneling which is the $4\times 4$ matrix whose components are given by $({\bf\Gamma})_{ij}=\Gamma\delta_{1i}\delta_{1j}$. $G^{R}_{11}(\omega)$ is then computed as

[TABLE]

Substituting this into Eq. (S86) gives the differential conductance at zero temperature,

[TABLE]

which exhibits peaks at $eV=\pm E_{0}$ and $\pm 3E_{0}$. If we restore the term $-2\pi l\hbar/T_{J}$ in Eq. (S83), the peaks would be at

[TABLE]

Therefore, we conclude that the Floquet theory gives a consistent result with the time-averaged differential conductance shown in Fig. 3 in the main text.

E. Case of commuting braiding operations

For an unambiguous demonstration of the relation of the conductance peak positions to the presence of non-Abelian operations, we consider similar tunneling experiments, that is, a tip is coupled to $\gamma_{1}$ at $t=t_{0}$ and the system Hamiltonian is periodic in time $T_{J}$, but with commutating operations of four MBSs. We explicitly show that the resulting conductance peak positions are different from those in Eq. (12) in the main text.

Let us introduce two different evolution operators, $W$ and $\bar{W}=\gamma_{1}W\gamma_{1}$, similar to the parity-conserving braiding operator $U_{c}$ and the tunneling-assisted braiding operator $\bar{U}_{c}$, respectively. The only difference compared to $U_{c}$ and $\bar{U}_{c}$ is that $W$ and $\bar{W}$ commute such that $[W,\bar{W}]=0$. This commutativity condition allows us to find the generic form of $W$ (and thus of $\bar{W}$) to be,

[TABLE]

up to an overall phase factor which does not affect the tunneling current. The rigorous derivation of this form of $W$ is given in the next section. $\bar{W}$ is obtained by interchanging $W_{e}$ and $W_{o}$ in the $W$ matrix. Here the general commuting braiding operations are characterized by the unit vector $\hat{n}_{w}$ and angles $\beta$ and $\beta^{\prime}$. Two related situations are drawn in Fig. S3.

The interfering terms in this case commute and thus can be expressed as

[TABLE]

indicating that the relative dynamics between $n$ and $n+1$ cycles adds a phase factor to the eigenstate of $\hat{n}_{w}\cdot\vec{\tau}$. In order to find its consequence, we calculate the anti-commutation relation,

[TABLE]

where $\hat{\gamma}_{1}(t_{q})=\left(W^{\dagger}\right)^{q}\gamma_{1}W^{q}$. By substituting this into $S(\varepsilon)+S(-\varepsilon)$ of Eq. (S10), the peak positions of the conductance in the low-bias voltage regime are found as

[TABLE]

where $j$ is an integer. If $0\leq|\beta-\beta^{\prime}|<\pi$, the peak separations $2|\beta-\beta^{\prime}|\hbar/T_{J}$ and $(2\pi-2|\beta-\beta^{\prime}|)\hbar/T_{J}$ appear alternately. If $\pi\leq|\beta-\beta^{\prime}|<2\pi$, then the separations of $(2|\beta-\beta^{\prime}|-2\pi)\hbar/T_{J}$ and $(4\pi-2|\beta-\beta^{\prime}|)\hbar/T_{J}$ are seen alternately. Note that for any value of $|\beta-\beta^{\prime}|$, these peak configurations cannot give rise to the results shown in Eq. (12) in the main text, manifesting the noncommutative structure of non-Abelian statistics.

F. Derivation of a generic form of $W$

We argue that the generic form of the matrix $W$ satisfying $\left[W,\bar{W}\right]=0$, shown in Eq. (S91), is

[TABLE]

where

[TABLE]

Here $\sigma_{0}$ is the $2\times 2$ identity matrix and ${\bf 0}$ is the null matrix. In this appendix, we prove this argument. Let us define $W$ as

[TABLE]

where the unit vectors $\hat{n}_{e}$ and $\hat{n}_{o}$ are

[TABLE]

$\bar{W}$ is then given by

[TABLE]

Case 1. $\beta$ or $\beta^{\prime}$ is equal to a multiple of $\pi$.

We show that the form of $W$ in Eq. (S95) holds for the following cases

• •

$\beta=m\pi$ and $\beta^{\prime}\neq n\pi$

• •

$\beta\neq m\pi$ and $\beta^{\prime}=n\pi$

• •

$\beta=m\pi$ and $\beta^{\prime}=n\pi$,

where $m$ and $n$ are integers. It is enough to consider the first case where $\beta=m\pi$ and $\beta^{\prime}\neq n\pi$ as the proofs for the other two cases are similar. The matrix $W_{e}$ in this case is

[TABLE]

independent of $\hat{n}_{e}$ and the commutation relation between $W$ and $\bar{W}$ is zero,

[TABLE]

We can rewrite $W_{e}$ as

[TABLE]

which is a valid expression if $\beta=m\pi$. Therefore, $W$ is written as

[TABLE]

which completes the proof in this case by changing the notation $\hat{n}_{o}$ by $\hat{n}_{w}$.

Case 2. $\beta\neq m\pi$ and $\beta^{\prime}\neq n\pi$.

In this case, we solve the problem

[TABLE]

Specifically, we need to solve

[TABLE]

Because $\text{sin}\beta\text{sin}\beta^{\prime}\neq 0$ in this case, $\left[\hat{n}_{e}\cdot\vec{\tau},\hat{n}_{0}\cdot\vec{\tau}\right]$ should be zero, which leads to

[TABLE]

As the Pauli matrices $\tau_{x,y,z}$ form an orthogonal basis, we have

[TABLE]

and

[TABLE]

which allows us to find

[TABLE]

where $s=1$ or $-1$. Therefore, $\hat{n}_{o}=s\hat{n_{e}}$ and $W$ in this case is obtained by

[TABLE]

By redefining notations $s\beta^{\prime}$ and $\hat{n}_{e}$ by $\beta^{\prime}$ and $\hat{n}_{w}$, we obtain Eq. (S95).