Finite temperature effects on Majorana bound states in chiral $p$-wave superconductors

TL;DR

This paper investigates how finite temperatures affect Majorana bound states in chiral p-wave superconductors, revealing that thermal effects cause algebraic decay in fermionic parity differences, but zero-temperature properties are recoverable below the minigap.

Contribution

It provides a detailed analysis of thermal effects on Majorana states, showing how temperature influences their stability and potential for quantum computing applications.

Findings

Thermal occupation causes algebraic decay of fermionic parity differences.

Zero-temperature properties are recovered below the minigap.

Temperatures above the minigap may still be compatible with topological quantum computation.

Abstract

We study Majorana fermions bound to vortex cores in a chiral $p$-wave superconductor at temperatures non-negligible compared to the superconducting gap. Thermal occupation of Caroli de Gennes-Matricon states, below the full gap, causes the free energy difference between the two fermionic parity sectors to decay algebraically with increasing temperature. The power law acquires an additional factor of $T^{-1}$ for each bound state thermally excited. The zero-temperature result is exponentially recovered well below the minigap (lowest-lying CdGM level). Our results suggest that temperatures larger than the minigap may not be disastrous for topological quantum computation. We discuss the prospect of precision measurements of pinning forces on vortices as a readout scheme for Majorana qubits.

Figures

Equations

$$\mathcal{H}\left(\begin{array}[]{c}u_{n}(\boldsymbol{r})\\ v_{n}(\boldsymbol{r})\end{array}\right)=\varepsilon_{n}/2\left(\begin{array}[]{c}u_{n}(\boldsymbol{r})\\ v_{n}(\boldsymbol{r})\end{array}\right),$$

$$\pazocal{H}=\begin{pmatrix}-\frac{1}{2m}\nabla^{2}-E_{F}&\frac{1}{2k_{F}}\big{\{}\Delta(\boldsymbol{r}),\partial_{x}+i\partial_{y}\big{\}}\\ -\frac{1}{2k_{F}}\{\Delta^{\ast}(\boldsymbol{r}),\partial_{x}-i\partial_{y}\}&\frac{1}{2m}\nabla^{2}+E_{F}\end{pmatrix}.$$

$$\hat{\gamma}_{n}=\int\mathrm{d}^{2}r\big{[}u_{n}^{\ast}(\boldsymbol{r})\hat{\psi}(\boldsymbol{r})+v_{n}^{\ast}(\boldsymbol{r})\hat{\psi}^{\dagger}(\boldsymbol{r})\big{]}.$$

$$\left(\begin{array}{c}{u}_0(\mathbf{r})\\ v_0(\mathbf{r})\end{array}\right)=\text{pazocal}N\frac{J_1\left(r\sqrt{2m{E}_F-1/{\xi }^2}\right)}{\cosh (r/\xi )}\left(\begin{array}{c}i{e}^{i\phi }\\ -i{e}^{-i\phi }\end{array}\right),$$

$$\Delta (\mathbf{r})={\Delta }_0\prod _{j=1}^{N}f(|\mathbf{r}-{\mathbf{R}}_j|)e^{i\arg (\mathbf{r}-{\mathbf{R}}_j)},$$

$${\epsilon }_0\approx \frac{4{\Delta }_0}{{\pi }^{3/2}}\frac{\cos \left(k_{F}R+\frac{\pi }{4}\right)}{\sqrt{k_{F}R}}{e}^{-R/\xi }.$$

$$\text{pazocal}{Z}_{\pm }$$

$$\displaystyle=\frac{1}{2}\pazocal{Z}_{0}\Big{[}1\pm\prod_{m}\tanh{(\beta\varepsilon_{m}/2)}\Big{]},$$

$$F_{\pm }=-{\beta }^{-1}\log \text{pazocal}{Z}_{\pm }.$$

$$\text{pazocal}Z$$

$$\times$$

$$\Delta F\equiv F_{-}-F_{+},$$

$$\Delta F={\epsilon }_0$$

$$+\frac{1}{\beta }\log \frac{e^{-\beta {\epsilon }_0}\cosh (\beta \frac{w_1+w_2}{2})+\cosh (\beta [{\delta }_{\epsilon }+\frac{w_2-w_1}{2}])}{e^{+\beta {\epsilon }_0}\cosh (\beta \frac{w_1+w_2}{2})+\cosh (\beta [{\delta }_{\epsilon }+\frac{w_2-w_1}{2}])}.$$

$$\Delta F\approx\varepsilon_{0}-\frac{4}{\beta}\cosh\big(\beta\frac{w_{1}+w_{2}}{2}\big{missing})\sinh(\beta\varepsilon_{0})e^{-\beta(\delta_{\varepsilon}+[w_{2}-w_{1}]/2)}.$$

$$\Delta F=\frac{1}{4}{\epsilon }_0{\beta }^2\left({\delta }_{\epsilon }^2+{\delta }_{\epsilon }(w_2-w_1)-w_1{w}_2\right)+\text{pazocal}O(\beta {\delta }_{\epsilon }{)}^4.$$

$$\Delta F$$

$$\approx {\epsilon }_0(\beta /2{)}^n\prod _{m=1}^n{\epsilon }_m\prod _{l=n+1}\tanh (\beta {\epsilon }_l/2)$$

$$\Delta F$$

$$\approx {\epsilon }_0\sqrt{2{\rho }_0{k}_{B}T}{e}^{-\frac{{\pi }^2}{4}{\rho }_0{k}_{B}T}{e}^{1+\text{pazocal}O(1/({\rho }_0{k}_{B}T{)}^2)}.$$

$${\tilde{F}}_{\pm }(k)\equiv \frac{1}{\sqrt{2\pi }}{\int }_0^{\infty }\mathrm{d}R{e}^{-ikR}{F}_{\pm }(R)$$

$$\Gamma \simeq k_{B}T{e}^{-{\delta }_{\epsilon }/(k_{B}T)}.$$

$${\epsilon }_{n-(\ell +1)/2}=-\left(n-\frac{\ell +1}{2}\right){\delta }_{\epsilon },$$

$${\delta }_{\epsilon }\approx \frac{2m{\Delta }_0^2}{k_F^2}\frac{{\int }_0^{\infty }\text{d}\rho f(\rho )/\rho \exp (-2{\int }_0^{\rho }\text{d}{\rho }^{\prime }f({\rho }^{\prime }))}{{\int }_0^{\infty }\text{d}\rho \exp (-2{\int }_0^{\rho }\text{d}{\rho }^{\prime }f({\rho }^{\prime }))},$$

$${\delta }_{\epsilon }\approx \frac{7\zeta (3)}{{\pi }^2}\frac{{\Delta }_0^2}{E_F},$$

$${\epsilon }_{\mathrm{top}}(R)$$

$${\epsilon }_{\mathrm{cl}}(R)$$

$${\tilde{\epsilon }}_{\mathrm{top}}(k)$$

$$\displaystyle+(i-\xi k)\sqrt{1-\left[\frac{k_{F}\xi}{\xi k-i}\right]^{2}}\sinh(\frac{1}{2}\mathrm{arctanh}\Big{(}\frac{k_{F}\xi}{\xi k-i}\Big{)})\Bigg{]}\Bigg{[}(i\xi k+1)^{3/2}\Big{(}1-\left[\frac{k_{F}\xi}{\xi k-i}\right]^{2}\Big{)}^{3/4}\Bigg{]}^{-1},$$

$$H=$$

$$=$$

Full text

I Introduction

In topological superconductors Majorana fermions appear as neutral zero energy excitations associated with defects in the superfluid-like domain walls and vortex cores Read and Green (2000); Gurarie and Radzihovsky (2007); Fu and Kane (2008); Akhmerov et al. (2009); Fu and Kane (2009). Majorana modes are candidate non-Abelian anyons Leinaas and Myrheim (1977); Wilczek (1982); Wen (1991), with direct applications to topological quantum computation O’Brien et al. (2018), and so these collective modes have stimulated intense fundamental and applied research Alicea (2012); Nayak et al. (2008).

Candidates for realizing Majorana fermions have included nanowire systems Kitaev (2001); Alicea et al. (2011); Mourik et al. (2012); Nadj-Perge et al. (2014); Higginbotham et al. (2015); Pawlak et al. (2016); Zhang et al. (2017); proximity-induced superconductivity on topological insulators or quantum anomalous Hall insulators Fu and Kane (2008); Koren et al. (2011); Xu et al. (2015); Chang et al. (2013); Kou et al. (2014); Bestwick et al. (2015); Feng et al. (2015); unconventional iron based superconductors Machida et al. (2018); Zhao et al. (2016); Wang et al. (2018) like (Li1-xFex)OHFeSe Liu et al. (2018); Niu et al. (2015). Interferometry protocols have been proposed to verify the existence of Majorana modes Lian et al. (2018); Beenakker et al. (2018); Fu and Kane (2009); Akhmerov et al. (2009), often subject to restrictions imposed by the (generally low) bulk quality of topological insulators Beidenkopf et al. (2011); Røising and Simon (2018).

In a $p$-wave superconductor two vortex bound states of finite separation $R$ hybridize and split in energy Cheng et al. (2009). The splitting effectively produces a finite-energy two-level system of the Majorana qubit, with the levels distinguished by the fermion parity: two Majoranas can ‘fuse’ into a state with either zero or one fermions (two states of different parity), like two spin-1/2 particles can combine into either a spin-0 or spin-1 state. In some cases, a barrier to exploiting the vortex bound states in topological superconductors is the mixing with the tower of excited states in the vortex cores, known as the Caroli de Gennes-Matricon (CdGM) states Caroli et al. (1964); Mizushima and Machida (2010); Möller et al. (2011). The CdGM states are characterized by a level spacing, the ‘minigap’, of size $\delta_{\varepsilon}\approx\Delta_{0}^{2}/E_{F}$, with $\Delta_{0}$ the full superconducting gap and $E_{F}$ the Fermi energy. Exciting these CdGM states does not lead to loss of quantum information, but it can hide the information so that it is very difficult to manipulate or to measure Akhmerov (2010) .

Consider a qubit made from two vortices each having one Majorana zero mode. To measure the state of this qubit, one might, in principle, measure the force between the vortices as we move them close together. At zero temperature there would be a difference in forces for the two different qubit states. At finite temperatures, excitation of the CdGM states reduces the force difference between the two different qubit states, although the quantum information is not lost until temperatures high enough that a bulk quasiparticle is excited that can carry away the fermionic parity Akhmerov (2010). It is thus generally assumed that temperature should be minimized as far as possible in realistic Majorana setups.

In this paper we examine the impact of thermal occupation of the CdGM states on the Majorana energy splitting in a two-vortex system. In a simple analytical model we find that the free energy difference between fermionic parity sectors lies exponentially close to the zero-temperature result for temperatures below the minigap, and decreases only as $T^{-2}$ with increasing temperature $T$ just above this temperature (and well below the superconducting gap). More generally, the contrast in free energy between parity sectors as two Majoranas are fused decreases with an additional factor of $T^{-1}$ for each thermally-occupied bound state. The polynomial dependence suggests that temperature need not be an immediately limiting factor for topological quantum computation. This result relies, however, on having a ‘sizable’ minigap. If the vortex core is swamped with a continuum of in-gap states, we find the parity contrast to be exponentially suppressed in temperature. We discuss experimental aspects, including limiting time scales from quasiparticle poisoning and thermal vortex motion, towards the end of the paper.

II Background: Majorana bound states in the $p+ip$ model

We consider an effective spinless $p_{x}+ip_{y}$ superconductor in two dimensions, described in the Bogoliubov-de Gennes (BdG) formalism in terms of a coupled eigensystem in the particle-hole basis Cheng et al. (2010),

[TABLE]

where $u_{n}$ ($v_{n}$) is the particle (hole) component of the eigenstate, and

[TABLE]

Here $\Delta(\boldsymbol{r})$ is the pairing function at position $\boldsymbol{r}=\left(x,y\right)$, $k_{F}$ and $E_{F}$ are the Fermi wavevector and Fermi energy, $m$ is the effective electron mass, and $\varepsilon_{n}/2$ is the energy of level $n$ of the BdG spectrum. The BdG Hamiltonian in Eq. (2) features a particle-hole-symmetry which ensures that the resulting spectrum is symmetric around zero, such that $\varepsilon_{n}>0$ denote the energy differences between the particle and the hole states. Formally, the model belongs to class D in the classification of non-interacting topological superconductors and insulators Schnyder et al. (2008); Altland and Zirnbauer (1997).

The quasiparticle annihilation (creation) operator $\hat{\gamma}_{n}^{(\dagger)}$ associated with level $\varepsilon_{n}$ of the BdG Hamiltonian is a superposition of the spinless electron annihilation (creation) operators $\hat{\psi}^{(\dagger)}(\boldsymbol{r})$,

[TABLE]

The operators associated with localized Majorana (zero) modes obey the defining criterion $\hat{\gamma}_{n}=\hat{\gamma}_{n}^{\dagger}$ and the anticommutation relations $\{\hat{\gamma}_{n},\hat{\gamma}_{m}\}=2\delta_{nm}$.

Vortices in the superconductor are regions around which the complex phase of the order parameter $\Delta$ winds through $2\pi$, with the gap magnitude going to zero in the vortex center. At position $\boldsymbol{r}$, defining $r=\absolutevalue{\boldsymbol{r}}$ and $\varphi=\arg\left(\boldsymbol{r}\right)$, this is described by the order parameter $\Delta(\boldsymbol{r})=\Delta_{0}e^{i\ell\varphi}f(r)$. Here $\Delta_{0}$ is the full (asymptotic) gap, $\ell$ the vorticity which is fixed to $\ell=1$ in the following, and $f(r)$ is the radial profile, where $f(r)\sim r^{\lvert\ell\rvert}$ close to the vortex core. Assuming a profile $f(r)=\tanh\left(r/\xi\right)$, where $\xi=v_{F}/\Delta_{0}$ is the coherence length and $v_{F}$ is the Fermi velocity, the Majorana zero mode solution of Eq. (1) for a vortex in an infinite system with $2mE_{F}>\Delta_{0}^{2}/v_{F}^{2}$ is explicitly given by Gurarie and Radzihovsky (2007); Cheng et al. (2010)

[TABLE]

manifestly satisfying the Majorana condition $u_{0}=v_{0}^{\ast}$. Here, $\pazocal{N}$ is a normalization constant, and $J_{1}$ is a Bessel function of the first kind. With $N$ vortices centred at positions $\boldsymbol{R}_{j}$, the order parameter may be approximated by

[TABLE]

when assuming no spatial or phase fluctuations of the vortices, limiting the validity to the (zero temperature) mean field result for the gap profile Cheng et al. (2009).

Within the ground state manifold defined by a collection of Majorana zero modes, the pairwise exchanges of excitations constitute a higher-dimensional representation of the (non-Abelian) braid group generators, which is encoded by the unitary operators $\hat{U}_{n,n+1}=\exp(-\frac{\pi}{4}\hat{\gamma}_{n}\hat{\gamma}_{n+1})$ Ivanov (2001); Alicea (2012). The braiding operators do not span all unitary gates needed for universal quantum computation, but the remaining set of gates can be implemented in a non-topological way with arbitrarily small errors Nayak et al. (2008); Bravyi (2006).

With two vortices of finite separation $R=\lvert\boldsymbol{R}_{1}-\boldsymbol{R}_{2}\rvert$ in a topological superconductor, each containing a Majorana mode, the energy splitting between the Majorana modes is given by Cheng et al. (2009):

[TABLE]

This expression holds in the regime $R\gg\xi\gg k_{F}^{-1}$, where $k_{F}$ is the Fermi wavevector. The hybridization energy above, $\varepsilon_{0}=\varepsilon_{+}-\varepsilon_{-}$, is the energy difference between the two fermionic parity sectors, which, when the vortices are well-separated, are associated with the two-vortex wavefunctions $\boldsymbol{\Psi}_{\pm}=(\boldsymbol{\Psi}_{1}\pm i\boldsymbol{\Psi}_{2})/\sqrt{2}$, where \boldsymbol{\Psi}_{n}=\big{(}u_{n}(\boldsymbol{r}),\hskip 2.84526ptv_{n}(\boldsymbol{r})\big{)}^{\intercal} (with ⊺ denoting the transpose) is the wavefunction associated with vortex $n$.

III Effect of heating on CdGM states

The lowest-energy levels, the CdGM levels, in the vortex core occur symmetrically about the Fermi energy, and their level spacing is given by the minigap, $\delta_{\varepsilon}$. The presence of such in-gap states, which derive from the finite density of states at the Fermi level in the normal state of the core, can lead to reduced distinguishability of the Majorana parity sectors Caroli et al. (1964); Cheng et al. (2009). In this section we consider the quantitative effect of occupying these states at finite temperature. More details on the CdGM levels can be found in Appendix A.

III.1 The partition function

In the superconducting condensate, where the number of particles is not conserved, we employ the grand canonical ensemble in evaluations of thermal expectation values. Since the (total) fermionic parity, i.e. the number of fermions modulo $2$, is conserved the proper partition function is projected onto the even ($+$) and odd ($-$) parity sectors Tuominen et al. (1992); Jankó et al. (1994); Beenakker et al. (2013); Law and Lee (2011):

[TABLE]

where $\pazocal{Z}_{0}=\prod_{n}2\cosh{(\beta\varepsilon_{n}/2)}$ is the partition function without parity restrictions, $\beta=(k_{B}T)^{-1}$, and all the products run over the positive energy levels in the BdG spectrum. We ignore states with energies far above $k_{B}T$ as they are suppressed by Boltzmann factors at low temperatures. The free energy in the $\pm$ parity sectors can be found in the usual way:

[TABLE]

III.2 Low-temperature model

We consider two vortices at finite separation such that the Majorana modes in the vortex cores split according to Eq. (6). In the same fashion the excited CdGM states of the cores will also split (see Sec. III.4). As a general and effective description at low temperatures, where only the lowest levels in the tower of CdGM states are thermally activated, we assume that the system has six levels: $\pm\varepsilon_{n}/2$, with $n\in\{0,1,2\}$, $\varepsilon_{1}\equiv\delta_{\varepsilon}-w_{1}$ and $\varepsilon_{2}\equiv\delta_{\varepsilon}+w_{2}$, where $\delta_{\varepsilon}$ is the minigap. Based on the splitting of the zero modes (Eq. (6)) and numerics the deviation of the excited states from the minigap, $w_{1}$ and $w_{2}$ as defined above, are expected to decay with the vortex-vortex separation as $w_{i}\sim\exp(-R/\xi)$ for $i=1,2$.

The two parity sectors $P=\pm 1$ are treated separately, with the associated particle configurations being denoted by $n_{\pm}$. In the BdG-spectrum these configurations are depicted in Fig. 1. Notice that the number of particle excitations (disks above zero energy in the figure) is therefore even or odd in the respective parity channels.

We consider the hierarchy of energy scales $\varepsilon_{0}\simeq\lvert w_{1}\rvert,\hskip 2.84526pt\lvert w_{2}\rvert<\delta_{\varepsilon}$, which we expect to apply when the vortices are further apart than about $3\xi$ (see Sec. III.4). The partition function is given by

[TABLE]

The observable consequence of switching fermionic parity sector is that the free energy involved in bringing two vortices together changes. We term this difference in free energy between the two parity sectors the parity disparity,

[TABLE]

which in the current case evaluates to:

[TABLE]

As expected, the parity disparity decreases monotonically with increasing temperature. It reflects the feasibility of directly probing state of a Majorana qubit. In the limit of low temperature, $k_{B}T\lesssim\delta_{\varepsilon}$, the leading correction to the zero-temperature result is suppressed as $\sim\exp(-\beta\delta_{\varepsilon})$:

[TABLE]

At temperatures above the minigap we find $F_{\pm}=-\frac{1}{\beta}\log{4}+\pazocal{O}(\beta\delta_{\varepsilon})$ to leading order in the high-temperature expansion, simply reflecting thermal occupation of the four available configurations in each parity sector shown in Fig. 1. The linear corrections in $\beta\delta_{\varepsilon}$ are independent of parity, such that the parity disparity acquires a leading order correction at second order in $\beta$,

[TABLE]

This inverse square law decay in the parity disparity with increasing temperature is weak enough that it will likely not preclude direct measurement at temperatures above the minigap. However, this toy model can only describe the regime $k_{B}T\lesssim 2\delta_{\varepsilon}$: at higher temperatures further CdGM states will be thermally excited.

Differences between the two parity channels are expected to be washed out at high temperatures. At low temperatures, the fact that the configurations $2_{+}$ and $3_{+}$ have excitation energies on the order of $\varepsilon_{0}$ greater than $2_{-}$ and $3_{-}$ (Fig. 1), respectively, distinguishes the two sectors.

III.3 Arbitrary numbers of core states

Now consider an arbitrary number of CdGM states thermally active in the two vortex cores. If the vortices are well-separated it is reasonable to assume that $\varepsilon_{0}<\varepsilon_{1}\approx\varepsilon_{2}<\varepsilon_{3}\approx\varepsilon_{4}<\dots$, where $\varepsilon_{0}$ is the Majorana energy, and $\varepsilon_{i>0}$ are the excited CdGM levels. The approximate equalities become exact as $R/\xi\to\infty$. If the temperature is much greater than at least one of the levels, $\varepsilon_{n}\ll k_{B}T\ll\Delta_{0}$ for $n\geqslant 0$, the parity disparity can be approximated by

[TABLE]

since $\tanh{(\beta\varepsilon_{m}/2)}\approx\beta\varepsilon_{m}/2\ll 1$ $\forall m\leqslant n$. Thus, the temperature decay above level $\varepsilon_{n}$ is algebraic: $\Delta F\sim T^{-n}$. The above decay suggests that temperatures on the order of the minigap can be tolerated but that exciting additional bound states causes the parity disparity to decrease rapidly.

The derivation holds in the limit of a well-defined and sizable minigap. If instead the vortex core contains a continuum of in-gap levels, the parity disparity of Eq. (14) is recast Higginbotham et al. (2015) and approximated 111Using that $F(x)=\int_{x}^{\infty}\mathrm{d}y\log\coth{y}=\frac{\pi^{2}}{24}$ $+$ \frac{1}{2}\big{(}\log\tanh{x}\log\tanh{(1+\tanh{x})} $+$ $\mathrm{Li}_{2}(1-\tanh{x})$ $+$ \mathrm{Li}_{2}(-\tanh{x})\big{)}, where $\mathrm{Li}_{2}$ is the dilogarithm function. For small $x$ this function behaves as $F(x)=\frac{\pi^{2}}{8}+x(\log{x}-1)+\pazocal{O}(x^{3})$ in the high-temperature limit as

[TABLE]

Here $\rho_{0}\equiv 1/\delta_{\varepsilon}$ is the density of states in the vortex core. The first term in the middle row of Eq. (15) appears when replacing the sum of equally spaced levels with an integral using the trapezoidal rule. The exponential suppression in $k_{B}T$ implies that vortex cores with densely packed in-gap states, as on topological insulator surfaces Xu et al. (2015) or in superconducting Pb monolayers Ménard et al. (2018), would likely remain unworkable in the context of readout schemes and topological quantum computation.

III.4 Numerical results

Consider two vortices with separation $R$ in a spinless $p+ip$ superconductor. We first assume that $T=0$ and that the two vortices have the same (positive) phase winding, as in Eq. (5). It should be noted that this is formally different from the vortex-antivortex and the antivortex-antivortex configuration in a $p+ip$ superconductor.

We approach the problem numerically by solving the BdG equations with a finite element method for a large range of inter-vortex distances on a slab no smaller than $[-10\xi,10\xi]\times[-9\xi,9\xi]$ with Dirichlet boundary conditions. In Appendix B the calculation is repeated with one vortex replaced by a hole, with a flux quantum penetrating the hole. We fix the (dimensionless) parameters of the $p+ip$ model to $E_{F}=3$, $\Delta_{0}=1$, and $mk_{F}^{2}=54$. The lowest-lying vortex states at zero temperature are shown in Fig. 2. Real-space colour maps of $\lvert\boldsymbol{\Psi}_{n}(\boldsymbol{r})\rvert=(\lvert u_{n}(\boldsymbol{r})\rvert^{2}+\lvert v_{n}(\boldsymbol{r})\rvert^{2})^{1/2}$ are shown in the figure insets. By following the trick of Ref. Caroli et al., 1964 the first predicted CdGM levels when using the approximation discussed in Appendix A are $\varepsilon_{1}/2\approx 0.142\Delta_{0}$, $\varepsilon_{2}/2\approx 0.284\Delta_{0}$, $\varepsilon_{3}/2\approx 0.426\Delta_{0}$, which are indicated with gray dotted lines in Fig. 2 and agree reasonably well with the numerical results, despite not strictly being in the BCS regime ($\Delta_{0}\ll E_{F}$) where the approximation is valid.

The CdGM wavefunctions are exponentially localized (with decay length $\xi$) around the two vortex cores, with spatial oscillations set by $2\pi k_{F}^{-1}$. Close to the vortex cores the expected small argument behaviour of the appropriate Bessel functions is recovered (e.g. Eq. (4)) Mizushima and Machida (2010). In the well-separated limit the two vortex cores each host a single Majorana zero mode, and the CdGM levels become doubly degenerate as reflected in Fig. 2. Using the Majorana basis Cheng et al. (2010) for intermediate separation reveals that the wavefunctions have started to separate into two localized states at separations as small as about $R/\xi\approx 3$. At larger separations the states $\boldsymbol{\Psi}_{\pm}$, corresponding to opposite fermionic parity, are reasonable approximations to the true ground state.

The near alignment of the crossing of the energy levels seen in Fig. 2 is linked to the radial profiles of the vortex states which are given by Bessel functions of the first kind Mizushima and Machida (2010), with an argument that increases slightly as a function of the excitation number. The Bessel-type wavefunctions enforce the energy splittings, which derive from the overlaps between the respective states Cheng et al. (2010), to acquire similar oscillations.

We note that on a finite slab with Dirichlet boundary conditions, the net angular momentum induced by the two vortices enforces boundary states peaked around the edges of the sample, with an energy spacing set by the boundary length, $v_{F}/L$. When the vortices are located far from the edge compared to $\xi$, these states are not affected by the presence of the vortices or their separation (see Appendix C). Considering instead periodic boundary conditions in both directions (solving the system on a torus) Read and Green (2000) and replacing one vortex by an antivortex Kraus et al. (2009), removes the boundary states completely.

In Figs. 3 and 4 we display the influence of temperature in the presence of in-gap states for the two-vortex system, when Eq. (8) is applied to the numerically found energy levels of Fig. 2. As the temperature surpasses the minigap, the oscillations begin to smear. Fig. 4 displays the numerically exact temperature dependence of the parity disparity at $R/\xi=3.0$ in the main figure. The black straight lines represent the simple power law derived below Eq. (14), when replacing $\tanh(\beta\varepsilon_{n}/2)$ by $\beta\varepsilon_{n}/2$ for $k_{B}T>\varepsilon_{n}$ and $1$ for $k_{B}T<\varepsilon_{n}$. Measuring the fusion channel of the corresponding qubit as defined by the two Majoranas, and hence the effect of braiding, thus becomes correspondingly difficult at temperatures well above the minigap.

Appendix B contains the same finite-temperature calculation with one vortex replaced by a hole. Qualitatively, the results are very similar to the vortex-vortex case. However, when one vortex is replaced by a hole, and reflection symmetry about $(\boldsymbol{R}_{1}+\boldsymbol{R}_{2})/2$ in the order parameter is lost, the previous double degeneracy of the levels in the $R/\xi\to\infty$ limit is lifted. The level spacing of the hole states is reduced with the hole circumference (cf. Ref. Akzyanov et al., 2014), generally making the issue of a small minigap worse when the hole radius is larger than $\xi$. We note also that the frequency of the CdGM level oscillations with the vortex-hole separation roughly doubles when comparing to the vortex-vortex case.

IV Measurement considerations

In practice, the topological splitting of the Majoranas constitutes a small contribution compared to the vortex-vortex repulsion (see Appendix C for a listing of the expressions). The topological contribution will in principle remain clear under the Fourier transform of the free energy. The (Friedel-like) oscillations of frequency $k_{F}$ from the topological splitting manifest as a bump, with distinguishable features in the two parity channels, revealed for instance in $\mathrm{Re}\{\tilde{F}_{-}-\tilde{F}_{+}\}$, with

[TABLE]

being the spatial Fourier transform of $F_{\pm}$. We show this in Fig. 5.

By subtracting the free energy before and after a parity flip one is left with a finite result only if the vortices contain non-Abelian anyons. Although the $k_{F}$ bump in the energy difference would constitute strong evidence of successful braiding, the subtraction of the two curves will require a sensitivity set by $\varepsilon_{0}$ in the measured free energy. The measurement series of $F_{-}-F_{+}$ can in principle be derived from the corresponding pinning force when one vortex is brought close to another, as measured before ($F_{\pm}$) and after ($F_{\mp}$) braiding with a third vortex.

In the following we consider the putative topological superconductor (Li1-xFex)OHFeSe to evaluate the scope of applicability for force measurements to readout of Majorana qubits Liu et al. (2018). The material has a small Fermi energy $E_{F}\simeq 50-60$ meV and a superconducting gap on the order of $2\Delta_{0}\simeq 20$ meV. Tunneling measurements on superconducting vortices identify a zero-bias peak, attributed to a Majorana zero mode, clearly separated from excited CdGM levels, with the unusually large minigap $\delta_{\varepsilon}\simeq 1$ meV $\approx 11$ K.

Numerically, we find the parity disparity at temperatures below the minigap to be roughly $\varepsilon_{0}\approx 0.04\Delta_{0}$ at $R/\xi=3$. To measure the corresponding pinning force difference between the parity sectors would require a force sensitivity of $\delta F\lesssim\varepsilon_{0}/\xi\approx 0.05$ pN, with $\xi\simeq 1.4$ nm Liu et al. (2018). This precision is at least an order of magnitude below previously reported force measurement thresholds Straver et al. (2008); Kremen et al. (2016); Veshchunov et al. (2016), although optimizing these measurements was not the primary goal of those studies.

IV.1 Timescales and vortex motion

One of the leading threats to topological quantum computation comes from free electrons mixing with the Majorana mode, causing an uncontrolled qubit parity flip Nayak et al. (2008); Akhmerov (2010). The timescale associated with this process is termed the ‘poisoning time’ Rainis and Loss (2012). The error rate $\Gamma$ is set by the energy scale Das Sarma et al. (2005)

[TABLE]

Using $\delta_{\varepsilon}=0.1\,$meV as a conservative estimate for the minigap of (Li1-xFex)OHFeSe, and $T=100$ mK, the expected poisoning time is $t_{p}\simeq\hbar/\Gamma\approx 9$ $\mu$s. Traversing a loop of radius $10\xi$, say, would require a vortex speed of $v\gtrsim 1$ cm s*-1*, which is in principle within reach of existing optical techniques Veshchunov et al. (2016). Another source of qubit decoherence at finite temperature, potentially relevant for this setup, is phonon interactions Aseev et al. (2019).

Thermal fluctuations change the vortex-vortex distance dynamically, causing smearing over an oscillating function (Eq. (6)). This can greatly reduce the contrast between the parity sectors Cheng et al. (2009) unless the vortex pinning potential changes significantly on the order of the oscillation lengthscale $2\pi k_{F}^{-1}\simeq 21$ nm. This approach would therefore require the vortices to be artificially pinned in tight potentials Embon et al. (2015) on this scale, if the temperature is on the order of the minigap.

Requiring adiabatic motion, to avoid the excitation of quasiparticles, introduces a lower time limit on the braiding operations. However, this timescale is small for the compound under consideration, $t_{a}\simeq\hbar/\delta_{\varepsilon}\approx 0.7\,$ps. The braiding operations are therefore restricted to timescales $t_{a}\ll t\ll t_{p}$. This could be achievable in the near future if improvements continue to be made to individual vortex manipulation Veshchunov et al. (2016); Straver et al. (2008); Kremen et al. (2016); Embon et al. (2015) for topological superconductors with sizable minigaps Liu et al. (2018); Wang et al. (2018).

V Conclusions

In this paper we have demonstrated the effect, on the Majorana energy splitting of two bound states, of thermally exciting the CdGM states in vortex cores. The parity disparity (the difference in free energies between the two Majorana parity channels when the Majoranas are brought close together), which reflects the state of a Majorana qubit, lies exponentially close to the zero-temperature result for temperatures well below the minigap $\sim\Delta_{0}^{2}/E_{F}$. Below the full gap, thermal excitation of higher CdGM states causes the amplitude to decay algebraically in temperature, with an additional factor of $T^{-1}$ for each CdGM state thermally occupied on average. If the in-gap states are densely packed the suppression in temperature becomes exponential.

The relatively weak decay of the parity disparity with increasing temperature means that temperature does not necessarily need to be minimized in experimental scenarios. In fact, the Majorana modes can in principle coexist with excited states without loss of quantum information in any readout scheme based on the total fermion parity Akhmerov (2010). Local heating can potentially be used as another degree of freedom. This could prove useful, for example, if the magnitude of the external magnetic field is limited by other factors (working at a higher temperature can decrease the critical field $H_{c1}$). Another possible application is to intentionally excite the CdGM modes in order to decouple the vortex motion from the state of the Majorana qubit, thus helping to increase decoherence times of the qubit.

As a promising candidate topological superconductor we draw specific attention to the intrinsic type-II superconductor (Li1-xFex)OHFeSe, recently suggested as a possible platform for topological quantum computation Liu et al. (2018). With the results presented in the above sections, this material should have an experimentally accessible temperature range in which one can aim at probing the parity disparity. In the near-term, we hope the knowledge that finite temperature effects need not be disastrous for the measurement and control of Majoranas may open other avenues of enquiry.

Acknowledgments

We thank Jay D. Sau for suggesting the outcome of a small minigap in Eq. (15). H.S.R. acknowledges discussions with J.M. Leinaas, and is grateful to T. Zhang for useful comments. H.S.R.  is supported by the Aker Scholarship. T.M. is supported by the Deutsche Forschungsgemeinschaft via the Emmy Noether Programme ME 4844/1-1 and through SFB 1143. S.H.S. is supported by EPSRC grant numbers EP/I031014/1 and EP/N01930X/1. F.F. acknowledges helpful discussions with P.J.W. Moll, S. Speller, and A. Akhmerov, and support from the Astor Junior Research Fellowship of New College, Oxford.

Appendix A The CdGM states

The low-energy spectrum of the vortex is given by Mizushima and Machida (2010)

[TABLE]

with $n\in\mathbb{Z}$ the angular momentum of the state, $\delta_{\varepsilon}$ the minigap, and $\ell\in\mathbb{Z}$ the vorticity. The excited core-localized CdGM states thus disperse linearly in momentum in $p$-wave superconductors. For odd vorticity the vortex hosts a zero energy mode of zero angular momentum relative to the condensate. The minigap for a $p$-wave superconductor can be estimated by approximating the vortex as a hard step in the gap function. By continuity of the wavefunction at the step radius one can deduce Caroli et al. (1964); Mizushima and Machida (2010)

[TABLE]

where $\rho=r/\xi$ is dimensionless length. Using $f(\rho)=\tanh(\rho)$ for the vortex profile of unit vorticity Tinkham (1975), and assuming $k_{F}^{2}\approx 2mE_{F}$ (which holds in the BCS regime $\Delta_{0}\ll E_{F}$), the above formula yields the level spacing

[TABLE]

with $\zeta$ the Riemann zeta function. The above formulas yield the gray dotted levels indicated in Fig. 2.

Appendix B The vortex-hole system

In this appendix we list numerical results, similar to those shown in Figs. 2 and 3, for a vortex-hole system. We use the order parameter magnitude $\lvert\Delta(\boldsymbol{r})\rvert=\Delta_{0}f_{1}(\lvert\boldsymbol{r}-\boldsymbol{R}_{1}\rvert)f_{2}(\lvert\boldsymbol{r}-\boldsymbol{R}_{2}\rvert)$, where $\Delta_{0}$ is the full gap, and $\boldsymbol{R}_{1}$ ($\boldsymbol{R}_{2}$) is the position of the vortex (hole). For the gap profiles we take $f_{1}(r)=\tanh(r/\xi)$ and $f_{2}(r)=\frac{1}{2}\left(1+\tanh(\alpha[r^{2}/\xi^{2}-\eta^{2}])\right)$ (which is chosen for numerical convenience and approximates a Heaviside step function when $\alpha$ is large), with $\eta=0.6$ and $\alpha=10$. The vorticity is $+1$ for both the vortex and the hole, and we solve the BdG equations on a finite slab with Dirichlet boundary conditions as described in Sec. III.4. The parameters of the $p+ip$ model are taken to be the same as in the aforementioned section.

In Fig. 6 the vortex-hole energy levels are shown, along with colour maps of the three lowest-lying wavefunctions. We note that unlike the vortex-vortex case where the CdGM levels are doubly degenerate in the $R/\xi\to\infty$ limit, the levels corresponding to the vortex and the hole are distinguishable. Accordingly, the wavefunctions of the excited states localize around either the hole (right in the figure) or the vortex (left in the figure). In Fig. 7 we show the impact of finite temperature on the parity disparity in the case of a vortex-hole system. Qualitatively, the temperature smearing here is similar to that of the vortex-vortex system in Fig. 3.

Appendix C Energy contributions

Here we list the full expression for the spatial Fourier transform (Eq. (16)) of the topological energy contribution (Eq. (6)), which in Sec. IV was suggested as a key signature for experimentally detecting the result of braiding. Second, we address the a priori potential problem of the bulk levels depending on the inter-vortex separation to conspire against the contribution of the Majoranas (Eq. (6)).

C.1 Fourier transforms of the two-vortex energy contributions

The contributions to the vortex-vortex energy are

[TABLE]

where $\varepsilon_{0,\mathrm{top}}=4\Delta_{0}\pi^{-3/2}$ from Eq. (6), $K_{0}$ is a hyperbolic Bessel function, $\varepsilon_{0,\mathrm{cl}}=\Phi_{0}^{2}/(4\pi\lambda)^{2}$ where $\Phi_{0}=h/(2e)$ is the flux quantum, and $\lambda$ is the London penetration depth Tinkham (1975). Defining the Fourier transform as in Eq. (16) yields for the topological contribution

[TABLE]

of which the real part is shown in Fig. 5. From the oscillating part of $\varepsilon_{\mathrm{top}}(R)$ there is a characteristic peak in $\tilde{\varepsilon}_{\mathrm{top}}(k)$ at $k=k_{F}$ that is enhanced with increasing $\xi$.

C.2 The background energy

The Hamiltonian associated with the mean field description of Eq. (2) can be written in diagonal form:

[TABLE]

where $\hat{\psi}^{(\dagger)}(\boldsymbol{r})$ is the annihilation (creation) operator of a spinless electron at position $\boldsymbol{r}$, and the quasiparticle annihilation (creation) operator associated with level $\varepsilon_{n}$ is $\hat{\gamma}^{(\dagger)}_{n}$.

The oscillations in the energy of the Majorana modes, as given in Eq. (6), could in principle be threatened by a conspiracy of the bulk energy levels, even at zero temperature (deriving from the last term in Eq. (27)). The sum of bulk energies could a priori have an oscillatory dependence on the inter-vortex separation $R$, like the individual Majorana and the CdGM levels, and thereby drown out the topological contribution. On a finite slab with Dirichlet boundary conditions, Fig. 8 shows the result of the energy sum for the $150$ lowest energy states, including all the low-lying edge states.

The oscillations of the Majorana levels at zero temperature $\pm\varepsilon_{0}$ are clearly visible on top of a $\sim\log(R/\xi)$ background that emerges. This background energy we assign to the ‘classical’ vortex-vortex repulsion that arises from the mutual Lorentz force of the flux lines Tinkham (1975). We checked that in the case of a vortex-antivortex system on a torus, where there are no edge states, the same background energy emerges with the opposite sign, i.e. an attractive contribution.