Non-Abelian Majorana fermions in topological $s$-wave Fermi superfluids

TL;DR

This paper analytically derives Majorana zero-modes in vortices of a 2D s-wave Fermi superfluid with spin-orbit coupling, showing their topological nature and non-Abelian statistics, which are relevant for quantum computation.

Contribution

It provides an analytical solution for Majorana zero-modes in s-wave superfluids with spin-orbit coupling, highlighting their topological properties and non-Abelian statistics.

Findings

Majorana zero-modes are localized in vortices of the superfluid.

The zero-modes obey non-Abelian statistics.

Energy splitting due to Majorana hybridisation is calculated.

Abstract

By solving the Bogoliubov--De Gennes equations analytically, we derive the fermionic zero-modes satisfying the Majorana property that exist in vortices of a two-dimensional $s$-wave Fermi superfluid with spin-orbit coupling and Zeeman spin-splitting. The Majorana zero-mode becomes normalisable and exponentially localised to the vicinity of the vortex core when the superfluid is topologically non-trivial. We calculate the energy splitting due to Majorana hybridisation and identify that the $s$-wave Majorana vortices obey non-Abelian statistics.

Equations

$$\begin{array}{rl}& \hat{H}(t)-\sum _{\sigma }{\mu }_{\sigma }{\hat{N}}_{\sigma }\\ & =\sum _{\alpha \beta }\int \mathrm{d}\mathbf{r}{\hat{\psi }}_{\alpha }^{†}(\mathbf{r},t)\left[{\hat{K}}_{\alpha \alpha }(\mathbf{r}){\delta }_{\alpha \beta }+{\hat{K}}_{\alpha \beta }(\mathbf{r})\right]{\hat{\psi }}_{\beta }(\mathbf{r},t)\\ & +\frac{1}{2}\sum _{\delta \gamma \alpha \beta }\int \mathrm{d}\mathbf{r}\int \mathrm{d}{\text{r}}^{\prime }\\ & {\hat{\psi }}_{\delta }^{†}(\mathbf{r},t){\hat{\psi }}_{\gamma }^{†}({\text{r}}^{\prime },t)U_{\delta \gamma \alpha \beta }(\mathbf{r},{\text{r}}^{\prime }){\hat{\psi }}_{\alpha }({\text{r}}^{\prime },t){\hat{\psi }}_{\beta }(\mathbf{r},t),\end{array}$$

$$\begin{array}{rl}H(\mathbf{r},t)& =\left(\begin{array}{cccc}{\hat{K}}_{\uparrow \uparrow }& {\hat{K}}_{\uparrow \downarrow }& 0& \Delta (\mathbf{r},t)\\ {\hat{K}}_{\downarrow \uparrow }& {\hat{K}}_{\downarrow \downarrow }& -\Delta (\mathbf{r},t)& 0\\ 0& -{\Delta }^{\ast }(\mathbf{r},t)& -{\hat{K}}_{\uparrow \uparrow }^{\ast }& {\hat{K}}_{\downarrow \uparrow }\\ {\Delta }^{\ast }(\mathbf{r},t)& 0& {\hat{K}}_{\uparrow \downarrow }& -{\hat{K}}_{\downarrow \downarrow }^{\ast }\end{array}\right)\end{array}$$

$$\begin{array}{rl}\Delta (\mathbf{r})& =V⟨{\hat{\psi }}_{\uparrow }(\mathbf{r}){\hat{\psi }}_{\downarrow }(\mathbf{r})⟩\\ & =V\sum _{\nu }[u_{\nu ,\uparrow }(\mathbf{r})v_{\nu ,\downarrow }^{\ast }(\mathbf{r})f_{\nu }\\ & +u_{\nu ,\downarrow }(\mathbf{r})v_{\nu ,\uparrow }^{\ast }(\mathbf{r})(1-f_{\nu })],\end{array}$$

$$\begin{array}{rl}{n}_{\sigma }(\mathbf{r})& =\sum _{\nu }\left[{\left|u_{\nu ,\sigma }(\mathbf{r})\right|}^2{f}_{\nu }+{\left|v_{\nu ,\sigma }(\mathbf{r})\right|}^2(1-f_{\nu })\right],\end{array}$$

$$H(\mathbf{r}){\text{h}}_{\nu }=E_{\nu }{\text{h}}_{\nu }$$

$$u_{0,\sigma }={\mathrm{e}}^{-\mathrm{i}\theta }{v}_{0,\sigma }^{\ast },$$

$${\hat{\gamma }}_0^{†}=\int \mathrm{d}\mathbf{r}\sum _{\sigma }\left[u_{0,\sigma }(\mathbf{r}){\hat{\psi }}_{\sigma }^{†}(\mathbf{r})+v_{0,\sigma }(\mathbf{r}){\hat{\psi }}_{\sigma }(\mathbf{r})\right]$$

$$\Delta (\mathbf{r})=\Delta (r){\mathrm{e}}^{\mathrm{i}\ell \phi },$$

$$\begin{array}{rl}{H}^{({\mathcal{U}}_{\mathrm{s}})}(\mathbf{r})\equiv {\mathcal{U}}_{\mathrm{s}}H(\mathbf{r}){{\mathcal{U}}_{\mathrm{s}}}^{-1}& =\left(\begin{array}{cc}{D}_{\downarrow }& M\\ M^{†}& -D_{\uparrow }\end{array}\right),\end{array}$$

$$M=\left(\begin{array}{cc}-\Delta (\mathbf{r})& -\mathrm{i}\lambda ({\partial }_y-\mathrm{i}{\partial }_x)\\ -\mathrm{i}\lambda ({\partial }_y+\mathrm{i}{\partial }_x)& {\Delta }^{\ast }(\mathbf{r})\end{array}\right).$$

$${\hat{K}}_{\uparrow \uparrow }{u}_{\nu ,\uparrow }+{\hat{K}}_{\uparrow \downarrow }{u}_{\nu ,\downarrow }+\Delta (\mathbf{r})v_{\nu ,\downarrow }$$

$$-{\hat{K}}_{\uparrow \uparrow }{v}_{\nu ,\uparrow }+{\hat{K}}_{\downarrow \uparrow }{v}_{\nu ,\downarrow }-{\Delta }^{\ast }(\mathbf{r})u_{\nu ,\downarrow }$$

$${\hat{K}}_{\downarrow \uparrow }{u}_{\nu ,\uparrow }+{\hat{K}}_{\downarrow \downarrow }{u}_{\nu ,\downarrow }-\Delta (\mathbf{r})v_{\nu ,\uparrow }$$

$${\hat{K}}_{\uparrow \downarrow }{v}_{\nu ,\uparrow }-{\hat{K}}_{\downarrow \downarrow }{v}_{\nu ,\downarrow }+{\Delta }^{\ast }(\mathbf{r})u_{\nu ,\uparrow }$$

$${\hat{K}}_{\uparrow \uparrow }{u}_{0,\uparrow }+{\hat{K}}_{\uparrow \downarrow }{u}_{0,\downarrow }+\Delta (\mathbf{r})u_{0,\downarrow }^{\ast }$$

$${\hat{K}}_{\downarrow \uparrow }{u}_{0,\uparrow }+{\hat{K}}_{\downarrow \downarrow }{u}_{0,\downarrow }-\Delta (\mathbf{r})u_{0,\uparrow }^{\ast }$$

$$u_{0,\sigma }={\mathrm{e}}^{\mathrm{i}{m}_{\sigma }\phi }{F}_{\sigma }(r){\mathrm{e}}^{-\frac{\Delta }{\lambda }r},$$

$$-F_{\uparrow }^{\prime \prime }$$

$$+\frac{m_{\uparrow }^2{F}_{\uparrow }}{r^2}+\frac{2m\lambda }{{\hslash }^2}\left(F_{\downarrow }^{\prime }+\frac{m_{\downarrow }{F}_{\downarrow }}{r}\right)=\frac{2m}{{\hslash }^2}{\mu }_{\uparrow }{F}_{\uparrow },$$

$$-F_{\downarrow }^{\prime \prime }$$

$$+\frac{m_{\downarrow }^2{F}_{\downarrow }}{r^2}+\frac{2m\lambda }{{\hslash }^2}\left(-F_{\uparrow }^{\prime }+\frac{m_{\uparrow }{F}_{\uparrow }}{r}\right)=\frac{2m}{{\hslash }^2}{\mu }_{\downarrow }{F}_{\downarrow }.$$

$$F_{\uparrow }(r)$$

$$+\left\{\frac{\lambda m}{{\hslash }^2}{F}_{\downarrow }^{\prime }(0)+\left[\frac{{\Delta }^2}{2{\lambda }^2}-\frac{(\mu +h)m}{2{\hslash }^2}\right]F_{\uparrow }(0)\right\}{r}^2$$

$$+\{\left[\frac{{\Delta }^3}{6{\lambda }^3}-\frac{\Delta (\mu +h)m}{2\lambda {\hslash }^2}-\frac{4\Delta \lambda m^2}{9{\hslash }^4}\right]F_{\uparrow }(0)$$

$$+\frac{11\Delta m}{9{\hslash }^2}{F}_{\downarrow }^{\prime }(0)\}{r}^3+\mathcal{O}(r^4),$$

$$F_{\downarrow }(r)$$

$$+\{\frac{2{\Delta }^2{\hslash }^4+{\lambda }^{2}m{\hslash }^2(h-\mu )-2{\lambda }^4{m}^2}{4{\lambda }^2{\hslash }^4}{F}_{\downarrow }^{\prime }(0)$$

$$+\frac{-8{\Delta }^2{\hslash }^{2}m+3h{\lambda }^2{m}^2+3{\lambda }^2\mu m^2}{12\lambda {\hslash }^4}{F}_{\uparrow }(0)\}{r}^3$$

$$+\mathcal{O}(r^4).$$

$$\lambda F_{\downarrow }^{\prime }(r)+\Delta F_{\downarrow }(r)-F_{\uparrow }(r)(\bar{\mu }+h)$$

$$\lambda F_{\uparrow }^{\prime }(r)+\Delta F_{\uparrow }(r)+F_{\downarrow }(r)(\bar{\mu }-h)$$

$$F_{\uparrow }(r)=C\text{exp}\left[-\frac{\Delta +\sqrt{h^2-{\bar{\mu }}^2}}{\lambda }r\right]=\frac{\bar{\mu }-h}{\sqrt{h^2-{\bar{\mu }}^2}}{F}_{\downarrow }(r),$$

$$\begin{array}{rl}{t}_{12}& =4\cos \left(\frac{{\theta }_1-{\theta }_2}{2}\right)\int \mathrm{d}\mathbf{r}G(\mathbf{r}){\sin }^2\left(\frac{{\phi }_1-{\phi }_2}{2}\right)\{\\ & {\Theta }_{\downarrow \uparrow }(\mathbf{r})-{\Theta }_{\uparrow \downarrow }(\mathbf{r})\}.\end{array}$$

$${\mathbf{u}}^{\prime }(r)=\mathcal{M}(r)\mathbf{u}(r),$$

$$\mathbf{u}(r)$$

$$\mathcal{M}(r)$$

$$\mathbf{u}(r)=\mathbf{u}(r_0){\mathrm{e}}^{{\int }_{r_0}^r\mathcal{M}(s)\mathrm{d}s},$$

$$\mathbf{u}(r)=\mathbf{u}(r_0){\mathrm{e}}^{{\sum }_{k=1}^{\infty }{\Omega }_k(s)},$$

Full text

Introduction – In two dimensions, the quasi-particle excitations of topologically ordered systems are generally non-Abelian anyons Fujimoto (2008). Possible realization of non-Abelian statistics has been studied in connection with the $\nu=5/2$ fractional quantum Hall (FQH) state Grosfeld and Stern (2006) and the vortex state of chiral $p_{\mathrm{x}}+\mathrm{i}p_{\mathrm{y}}$ superconductors and superfluids Read and Green (2000); Volovik (1999). Majorana fermions in the cores of superfluid vortex excitations have been actively considered in the context of $p$-wave pairing Tewari et al. (2008); Silaev and Volovik (2014), but can also occur in an $s$-wave superconductor coupled by the proximity effect to a topological insulator Wilczek (2009). Majorana vortices in the topological $s$-wave superconductor are non-Abelian, and in the same topological class Fujimoto (2008) with the Moore-Read Pfaffian FQH state Moore and Read (1991); Read and Moore (1992), $p_{\mathrm{x}}+\mathrm{i}p_{\mathrm{y}}$ superconductors Ivanov (2001); Kitaev (2006), and the gapped non-Abelian spin liquid phase of the Kitaev model Lee et al. (2007). A key benefit of non-Abelian exchange statistics is that the Majorana vortices can be used to realize braiding operations. Non-commutative braiding, which amounts to exchanging two Majorana vortices adiabatically Leijnse and Flensberg (2012), is a key ingredient needed for performing quantum logic operations on a fault-tolerant topological quantum computer Nayak et al. (2008).

In ultra-cold atom experiments Mukherjee et al. (2017); Kinast et al. (2004); Yefsah et al. (2013), however, the densities and temperatures are low, and the scattering between fermions takes place typically in the $s$-wave channel Petrov et al. (2004). Experimentally, the ultra-cold atomic $s$-wave Fermi gas is a highly flexible quantum many-body system whose interactions, spin balance, and trapping geometries can be tuned nearly arbitrarily Zwierlein et al. (2005); Ku et al. (2016). In the mean-field picture, it is possible to achieve band inversion leading to a topological phase Tewari et al. (2011); Fujimoto (2008) by combining the effects of two-dimensional spin-orbit (SO) coupling Cheuk et al. (2012); Wang et al. (2012), and spin imbalance Zwierlein et al. (2006); Partridge et al. (2006). Recently, tunable 2D SO coupling was demonstrated experimentally Huang et al. (2016). These advances open up a promising perspective for creating a highly-controlled atomic topological superfluid Liu and Hu (2012); Wei and Mueller (2012) in the laboratory in the near future.

Despite the promising recent experimental progress, Majorana zero-modes in vortices of an $s$-wave Fermi superfluid with SO coupling and a Zeeman field remain relatively poorly understood. Here, we derive analytically the Majorana vortex zero-mode, and use it to calculate the energy splitting due to inter-vortex tunnelling. The tunnelling generally lifts the zero-mode degeneracy Bonderson (2009), and has been studied for the Moore-Read state Baraban et al. (2009), Kitaev’s honeycomb model Lahtinen et al. (2008), and $p$-wave superconductors Cheng et al. (2009, 2010); Mizushima and Machida (2010). Hybridisation of Majorana fermions in dense vortex lattices gives rise to a band structure Murray and Vafek (2015); Liu and Franz (2015); Silaev (2013). Of particular interest for experimentally controlled quantum simulation are topologically non-trivial bands of Majoranas Kraus and Stern (2011), and the possibility of flattening them Toikka and Andreanov (2018).

Physical system – The Hamiltonian of a two-dimensional Fermi gas with spin-orbit coupling and spin-imbalance reads Ketterson and Song (1999)

[TABLE]

where $\hat{N}_{\sigma}=\int\mathrm{d}\mathbf{r}\,\hat{\psi}_{\sigma}^{\dagger}(\mathbf{r},t)\hat{\psi}_{\sigma}(\mathbf{r},t)$ is the total number of atoms with spin $\sigma=\left\{\uparrow,\downarrow\right\}$, the operator-valued field $\hat{\psi}_{\sigma}^{\dagger}(\textbf{r},t)$ obeys Fermi anti-commutation relations and creates a spin-$\sigma$ fermion at location $\mathbf{r}$ at time $t$, $\hat{K}_{\sigma\sigma}(\mathbf{r})=-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu_{\sigma}$, $\bar{\mu}=\left(\mu_{\uparrow}+\mu_{\downarrow}\right)/2$, $h=\left(\mu_{\uparrow}-\mu_{\downarrow}\right)/2$ so that $\mu_{\uparrow}=\bar{\mu}+h$ and $\mu_{\downarrow}=\bar{\mu}-h$. Here $h$ is a Zeeman field creating an energy splitting between the spin components. The SO coupling is represented by $\hat{K}_{\uparrow\downarrow}=-\mathrm{i}\lambda(\partial_{y}+\mathrm{i}\partial_{x})$, $\hat{K}_{\downarrow\uparrow}=-\mathrm{i}\lambda(\partial_{y}-\mathrm{i}\partial_{x})$ with $\hat{K}_{\uparrow\downarrow}=\hat{K}_{\downarrow\uparrow}^{\dagger}$ and $\hat{K}_{\uparrow\downarrow}=-\hat{K}_{\downarrow\uparrow}^{*}$.

We now assume that the two-body interaction is spin-independent, $U_{\delta\gamma\alpha\beta}(\mathbf{r},{\textbf{r}^{\prime}})=U(\mathbf{r},{\textbf{r}^{\prime}})\delta_{\alpha\gamma}\delta_{\delta\beta}$, and represents contact interactions, $U(\mathbf{r},{\textbf{r}^{\prime}})=-V(\mathbf{r})\delta(\mathbf{r}-{\textbf{r}^{\prime}})$, where the sign convention is chosen such that $V(\mathbf{r})>0$ corresponds to attractive contact interactions.

Superfluidity in spin-balanced Fermi gases results from the formation of Cooper pairs, bound states of two fermions around the Fermi surface with opposite momenta $+{\textbf{k}}$ and $-{\textbf{k}}$. For spin-imbalanced pairing Wlazłowski et al. (2018), theoretical predictions such as FFLO Fulde and Ferrell (1964); Larkin and Ovchinnikov (1964) and Sarma phases Sarma (1963) together with deformed Fermi surface superfluidity Müther and Sedrakian (2002) have been presented. In the FFLO phase the Cooper pairs have a non-zero center-of-mass momentum which gives rise to a periodic spatial modulation of the order parameter and also possibly of density. While spin-imbalance can thus influence the pairing and subsequently the dynamics of topological defects by frustrating the Cooper pairing, in experiments the excess non-superfluid particles are typically spatially separated from the completely paired BCS superfluid Partridge et al. (2006); Liao et al. (2010); Shin et al. (2006).

In typical current experimental conditions in ultra-cold atoms, the pairs form in the $s$-wave state with zero angular momentum. The Pauli principle then necessitates that the spin state be a singlet. Since the Cooper pairs have no angular momentum and no net spin, the relative pair wavefunction can be fully characterised by a single complex amplitude $\Delta(\mathbf{r},t)$.

Ignoring now all the other quantum correlations apart from the pair correlations captured by the order parameter, the mean-field Hamiltonian describing quasi-particle excitations in the superfluid Fermi gas without spin-rotational invariance reads $\mathcal{H}(t)=\int\mathrm{d}\mathbf{r}\,\begin{pmatrix}\bm{\hat{\Psi}}^{\dagger}(\mathbf{r},t)&\bm{\hat{\Psi}}(\mathbf{r},t)\end{pmatrix}H(\mathbf{r},t)\begin{pmatrix}\bm{\hat{\Psi}}(\mathbf{r},t)&\bm{\hat{\Psi}}^{\dagger}(\mathbf{r},t)\end{pmatrix}^{\mathrm{T}}$, where $\bm{\hat{\Psi}}^{\dagger}(\mathbf{r},t)=\begin{pmatrix}\hat{\psi}_{\sigma}^{\dagger}(\mathbf{r},t)&\hat{\psi}_{\sigma^{\prime}}^{\dagger}(\mathbf{r},t)\end{pmatrix}$ is the 4-component Nambu spinor, and

[TABLE]

is the Bogoliubov–De Gennes (BdG) matrix representation. In what follows, we are not interested in explicit time dependence and therefore drop the labels.

In the mean-field picture, we replace the many-body problem (1) with the effective single-particle Hamiltonian (2) parametrised by the pair potential (order parameter) $\Delta(\mathbf{r})$. The pair potential together with the chemical potential is determined by optimising the effective Hamiltonian (2) such that it minimises the total free energy. The result for the pair potential is

[TABLE]

and for the atomic densities

[TABLE]

where $f_{\nu}=1/\left[\mathrm{e}^{E_{\nu}/(k_{\mathrm{B}}T)}+1\right]$ is the Fermi distribution for occupations and $T$ is the temperature. We set $T=0$. We have expressed the order parameter in terms of the amplitudes $u$ and $v$, defined by the Bogoliubov-Valatin transformations $\hat{\psi}_{\sigma}(\mathbf{r})=\sum_{\nu}\left[u_{\nu,\sigma}(\mathbf{r})\hat{\gamma}_{\nu}+v^{*}_{\nu,\sigma}(\mathbf{r})\hat{\gamma}_{\nu}^{\dagger}\right]$, $\hat{\gamma}_{\nu}^{\dagger}=\int\mathrm{d}\mathbf{r}\,\sum_{\sigma}\left[u_{\nu,\sigma}(\mathbf{r})\hat{\psi}_{\sigma}^{\dagger}(\mathbf{r})+v_{\nu,\sigma}(\mathbf{r})\hat{\psi}_{\sigma}(\mathbf{r})\right]$, which diagonalise the effective single-particle mean-field Hamiltonian.

Symmetries of the Hamiltonian – To obtain the static properties of Majorana fermions at the vortex cores in a topological $s$-wave Fermi superfluid, we need to solve the fundamental eigenvalue equation

[TABLE]

for the zero-mode $\nu=0$ with $E_{0}=0$. When considering equations of the form (5), it is instructive to first understand the symmetries of the underlying Hamiltonian.

The BdG Hamiltonian $H(\mathbf{r},t)$ possesses a particle-hole symmetry (PHS) defined by $\{\Xi_{\mathrm{s-wave}},H(\mathbf{r},t)\}=0$ with $\Xi_{\mathrm{s-wave}}=\mathrm{e}^{\mathrm{i}\theta}\tau_{\mathrm{x}}\otimes 1_{\sigma}K$, where $K$ is the complex conjugation operator (in momentum space, $K$ additionally flips the sign of the momentum), $\tau_{\mathrm{y}}$ ($\sigma_{\mathrm{z}}$) is the Pauli matrix in particle-hole (spin) space, and $\theta\in\mathbb{R}$. As a result of the particle-hole symmetry, if $\textbf{h}_{\nu}=\begin{pmatrix}u_{\nu,\uparrow}&u_{\nu,\downarrow}&v_{\nu,\uparrow}&v_{\nu,\downarrow}\end{pmatrix}^{\mathrm{T}}$ is a solution of Eq. (5) with energy $E_{\nu}$, then $\Xi_{\mathrm{s-wave}}\textbf{h}_{\nu}=\mathrm{i}\,\mathrm{e}^{\mathrm{i}\theta}\begin{pmatrix}-v_{\nu,\uparrow}^{*}&v_{\nu,\downarrow}^{*}&u_{\nu,\uparrow}^{*}&-u_{\nu,\downarrow}^{*}\end{pmatrix}^{\mathrm{T}}$ is a solution with energy $-E_{\nu}$.

Additionally, if and only if $\Delta(\mathbf{r},t)\in\mathbb{R}$ and $h=0$, the BdG Hamiltonian possesses a time-reversal symmetry defined by $\left[T,H(\mathbf{r},t)\right]=0$, where $T=\mathrm{e}^{\mathrm{i}\theta}1_{\tau}\otimes\sigma_{\mathrm{y}}K$ ($T^{2}=-1$, for half-integer spin), and $\theta$ is arbitrary. We have defined $1_{\tau}$ to be the identity matrix in the particle-hole space, while $\sigma_{\mathrm{y}}$ is the Pauli matrix in spin space.

A straightforward calculation shows that non-degenerate zero-modes (that is, the modes $H(\mathbf{r})\textbf{h}_{0}=E_{0}\textbf{h}_{0}$ with $E_{0}=0$), must always satisfy the important symmetry

[TABLE]

which is the so-called Majorana property making the zero-mode special. The Majorana property ensures that the quasi-particle operators

[TABLE]

satisfy $\hat{\gamma}_{0}^{\dagger}=\hat{\gamma}_{0}$; they represent Majorana fermions.

Let us temporarily assume that time-reversal symmetry is respected, i.e. $\Delta(\mathbf{r})\in\mathbb{R}$ and $h=0$. If $\textbf{h}_{\nu}=\begin{pmatrix}u_{\nu,\uparrow}&u_{\nu,\downarrow}&v_{\nu,\uparrow}&v_{\nu,\downarrow}\end{pmatrix}^{\mathrm{T}}$ is a solution of Eq. (5) with energy $E_{\nu}$, then $T\textbf{h}_{\nu}$ is a solution with the same energy $E_{\nu}$. This is the pair of the zero-energy Majorana mode. With time-reversal symmetry, the Majorana probability densities are identically overlapping in space: $|\textbf{h}_{0}|^{2}=|T\textbf{h}_{0}|^{2}=2|u_{0,\uparrow}|^{2}+2|u_{0,\downarrow}|^{2}$. We need to break the time-reversal symmetry to separate them spatially, for example, by having a vortex in the order parameter.

A broken time-reversal symmetry $T$, a broken spin-rotational symmetry (due to spin-orbit coupling), and an unbroken PHS $\Xi_{\mathrm{s-wave}}$ result in the BdG Hamiltonian $H(\mathbf{r})$ belonging to the symmetry class D in the Altland-Zirnbauer classification Altland and Zirnbauer (1997). According to Ref. Teo and Kane (2010), point defects (such as the vortices in the order parameter considered here) with a Hamiltonian that belongs to the class D are associated with a $\mathbb{Z}_{2}$-valued topological invariant. This means that in the topologically non-trivial regime the number of exact Majorana zero-modes is given by $W\text{ mod }2$, where $W$ is the sum over all vortex winding numbers, with the rest hybridising in pairs forming a band structure. Generally, when multiple vortices are brought close together, the Majorana zero-modes hybridize into a band structure, leading to complex fermion states at positive and negative energy (Murray and Vafek, 2015; Liu and Franz, 2015; Silaev, 2013). We calculate the energy splitting in Eq. (18). For periodic vortex lattices, such as the square and the triangular lattice, the resulting low-energy theory is typically gapped and topologically non-trivial Liu and Franz (2015).

Eigenvalue equation for the zero-mode – We now include a vortex in the order parameter,

[TABLE]

where $r$ and $\varphi$ are polar coordinates centred on the vortex. Energy considerations require that $\Delta(\textbf{0})=0$ at the vortex core. The integer $\ell$ denotes the vorticity, and $\Delta(r)$ is a real function of $r$ that vanishes at $r=0$. In what follows, we will solve Eq. (5) analytically for the zero-mode $\textbf{h}_{0}$ that exists in the core of the vortex (8).

For convenience, let us first apply the unitary transformation $\mathcal{U_{\mathrm{s}}}$, which transforms the basis as $\mathcal{U_{\mathrm{s}}}\begin{pmatrix}u_{0,\uparrow}&u_{0,\downarrow}&v_{0,\uparrow}&v_{0,\downarrow}\end{pmatrix}^{\mathrm{T}}=\begin{pmatrix}u_{0,\downarrow}&v_{0,\downarrow}&v_{0,\uparrow}&u_{0,\uparrow}\end{pmatrix}^{\mathrm{T}}$, and shuffles the zeros in the rows and columns of the Hamiltonian $H(\textbf{r})$ [Eq. (2)] to give

[TABLE]

where $D_{\sigma}=\text{diag}(\hat{K}_{\sigma\sigma},-\hat{K}_{\sigma\sigma})$ and

[TABLE]

Setting $D_{\sigma}=0$ decouples the spins, and the Hamiltonian reduces to a $2\times 2$ BdG matrix corresponding to the Fu-Kane model at neutrality, whose Majorana modes have been considered Chiu et al. (2015). This is the low-energy regime near $\mu_{\sigma}=0$, which is the topological transition point in typical $p$-wave systems. The Majorana zero-modes have also been considered for the Fu-Kane model at $\mu\neq 0$ Cheng et al. (2010).

However, here we consider the full $4\times 4$ structure of Eq. (5):

[TABLE]

When looking for zero-modes, the system (11) can be simplified by using the symmetry property (6), that is, we seek simultaneous eigenstates of the PHS. There are two possibilities for $\theta$, namely $\theta=\theta_{0},\theta_{0}+\pi$, which both satisfy the Majorana relation. This means that with the final solution we must allow for both possibilities. In both cases, we obtain

[TABLE]

a coupled $2\times 2$ system, which must be solved for $u_{0,\uparrow}$ and $u_{0,\downarrow}$.

In polar coordinates, the spin-orbit coupling terms read $\hat{K}_{\uparrow\downarrow}(\mathbf{r})=\lambda\mathrm{e}^{-\mathrm{i}\varphi}\left[\partial_{r}-(\mathrm{i}/r)\partial_{\varphi}\right]$. Observing the azimuthal symmetry, let us try a separable ansatz with angular momentum eigenstates,

[TABLE]

with the assumption that $\Delta$ has no $r$-dependence. Physically, this approximation treats the vortices as points with only a phase profile. Substitution into the system (12) shows that we can eliminate the angular dependence by taking $m_{+}=\zeta$ and $m_{-}=-\ell$, where $m_{\pm}\equiv\pm m_{\sigma^{\prime}}-m_{\sigma}$, $\zeta=\pm 1$ with $\zeta=+1$ for $\sigma=\uparrow$ and $\zeta=-1$ for $\sigma=\downarrow$, and $\sigma^{\prime}=\uparrow,\downarrow$ if $\sigma=\downarrow,\uparrow$ respectively. This implies $m_{\sigma^{\prime}}=(\ell+\zeta)/2$ and $m_{\sigma}=(\ell-\zeta)/2$. The only requirement here is that $\ell\in\mathbb{Z}$. Explicitly, the coupled system (12) then reads

[TABLE]

Below, we solve analytically the coupled system (14) for $F_{\uparrow}(r)$ and $F_{\downarrow}(r)$.

Series solution for $F_{\sigma}(r)$ – The vortex core ($r=0$) is an irregular singular point of Eq. (14) where numerical methods are inherently unstable. An obvious solution is to start the numerical integration at a point $r_{0}>0$ thus avoiding the singularity, but in this case we need to know accurately the initial condition at the rather arbitrary point $r_{0}$. It is therefore highly desirable to have analytic solutions, at least around the vortex core. In what follows we solve the system (14) using the Wasow method Wasow (2002), a generalisation of the Fröbenius series method for coupled systems of differential equations with irregular singular points.

The details for solving the system (14) are shown in Appendix A. For definitess, we set $\ell=+1$. The result is a series solution for $F_{\sigma}(r)$ around the origin that can be evaluated analytically to arbitrary order, and reads

[TABLE]

While $F_{\uparrow}(0)$ can be finite, the boundary condition $F_{\downarrow}(0)=0$ is forced to keep $F_{\downarrow}$ finite at the origin, which results from the series solution for $F_{\downarrow}$ starting with the power $r^{-1}$ whereas that for $F_{\uparrow}$ starts with $r^{0}$. That one spin component is zero and the other finite at the vortex core has been observed in all numerical studies of the Majorana zero-mode in $s$-wave Fermi gases Liu et al. (2012); Iskin (2012). In principle the system (14) needs four initial conditions, but requiring that $F_{\downarrow}$ be finite at the origin consumes two of them leaving only the two initial conditions $F_{\uparrow}(0)$ and $F^{\prime}_{\downarrow}(0)$ in the solution (15).

Existence of the Majorana zero-mode – The existence of the analytical series solutions (15) alone does not guarantee the existence of a topologically protected Majorana zero-mode. The condition for the physical existence of the zero-mode is that it remain normalisable as $r\to\infty$.

In the limit $r\to\infty$, the kinetic energy and terms in $1/r$ can be dropped in system (14) to obtain

[TABLE]

Physically, we are interested in only the solutions for system (16) that are normalisable. For example, if $\lambda>0$ and $\Delta<0$, the only normalisable asymptotic solution is

[TABLE]

where $C$ is a constant. Even then, the solution (17) is normalisable only when $h>\sqrt{\bar{\mu}^{2}+|\Delta|^{2}}$. The value at equality corresponds to the critical Zeeman field marking the phase transition to the topological regime with a uniform order parameter, $\Delta(\mathbf{r})=\Delta$ Brand et al. (2018), which here is associated with the boundedness of the Majorana zero-mode at infinity. However, it is not straightforward to find initial conditions $F_{\uparrow}(0)$ and $F^{\prime}_{\downarrow}(0)$ such that the zero-mode has the desired large-$r$ asymptotics.

In principle, $\bar{\mu}$ and $\Delta$ must be obtained self-consistently in conjunction with solving the BdG eigenvalue equation. Considering a uniform vortex-free system, we can solve Eqs. (3), (4), and (11) self-consistently to find that the superfluid is topologically non-trivial with the parameters $\lambda=5$, $h=10$, $\Delta=-1.587$, $\bar{\mu}=-5.987$. Here all units are measured with respect to $\hbar=m=1$. We have used the pair binding energy of $E_{\mathrm{b}}=0.05$ in the renormalisation of the gap equation, which otherwise diverges logarithmically. The pair binding energy determines the state of the Cooper pairs across the BCS-BEC crossover Brand et al. (2018). From Eq. (17), therefore, we know that with these parameters there is an exponentially decaying asymptotic solution, and we have performed a direct search for $F_{\uparrow}(0)$ and $F^{\prime}_{\downarrow}(0)$ that match with this asymptotic behaviour.

We retain terms upto $\mathcal{O}(r^{10})$, and use the analytical solution (15) to evaluate accurate initial conditions for a numerical integration of the system (14). The truncated analytical series solution agrees with the numerical integration for $r\lesssim 0.5$ (Fig. 1). The analytic series expansion can be developed to arbitrary precision. As expected, the numerical solution becomes unstable as the vortex core is approached, but more importantly it can be used to study the boundedness of the mode as $r\to\infty$. We find a bounded zero-mode in vortices with unit positive circulation with the specific parameter values above for $3.488885<F_{\uparrow}(0)<3.488886$ and $F^{\prime}_{\downarrow}(0)=-12.518$. The full solution including azimuthal dependence is then given by Eq. (13).

Majorana energy splitting – We now use the Majorana zero-mode to calculate the energy splitting in an $s$-wave Fermi gas that results from hybridisation between the Majorana modes of two vortices located at points $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$. Then, $\Delta(\mathbf{r})\approx\Delta_{i}(\mathbf{r})$ near vortex $i$ ($i=1,2$), where $\Delta_{i}(\mathbf{r})$ is the single-vortex pair function (8) for vortex $i$. Introduction of a single-vortex background phase $\theta_{i}$, $\Delta_{i}(\mathbf{r})=\Delta\mathrm{e}^{\mathrm{i}(\ell_{i}\varphi_{i}+\theta_{i})}$ with $\varphi_{i}$ the azimuthal angle measured with respect to vortex $i$, amounts to the Majorana mode changing by $u_{0,\sigma}(\mathbf{r}-\mathbf{r}_{i})\to\mathrm{e}^{\mathrm{i}\frac{\theta_{i}}{2}}u_{0,\sigma}(\mathbf{r}-\mathbf{r}_{i})$, leaving invariant the system (12).

In the tight-binding approximation, the Hamiltonian $H_{\Delta}$ describing hopping between the two vortices, with a straightforward generalisation to lattices of Majorana vortices, is given by $H_{\Delta}=\mathrm{i}\,t_{12}\,\hat{\gamma}_{0,1}\hat{\gamma}_{0,2}$, where $\hat{\gamma}_{0,i}$ is the Majorana mode (7) at vortex $i$, and $t_{12}=\langle\textbf{g}_{0}^{(2)}|H^{(\mathcal{U_{\mathrm{s}}})}|\textbf{g}_{0}^{(1)}\rangle=\langle\textbf{g}_{0}^{(2)}\left|\begin{pmatrix}D_{\downarrow}&M\\ M^{\dagger}&-D_{\uparrow}\end{pmatrix}\right|\textbf{g}_{0}^{(1)}\rangle$ is the overlap integral that gives the energy splitting. Here $\textbf{g}_{0}^{(i)}(\mathbf{r})=\begin{pmatrix}u_{0,\downarrow}(\mathbf{r}-\mathbf{r}_{i})&u^{*}_{0,\downarrow}(\mathbf{r}-\mathbf{r}_{i})&u^{*}_{0,\uparrow}(\mathbf{r}-\mathbf{r}_{i})&u_{0,\uparrow}(\mathbf{r}-\mathbf{r}_{i})\end{pmatrix}^{\mathrm{T}}$ is the Majorana zero-mode centered at vortex $i$.

For convenience, we define $M^{(1)}\equiv M-\tilde{M}$ such that $M^{(1)}$ coincides with Eq. (10) when the order parameter is given by just a single vortex at position $\mathbf{r}_{1}$. It follows that near vortex 1 $\tilde{M}\approx 0$, and the dominant contribution to the overlap integral comes from the vicinity of vortex 2. Taking $\ell_{1}=\ell_{2}=1$, using Eq. (13), and introducing the definitions $G(\mathbf{r})\equiv-\Delta\,\mathrm{exp}\left[-\frac{f^{2}}{\lambda}\left(|\mathbf{r}-\mathbf{r}_{2}|+|\mathbf{r}-\mathbf{r}_{1}|\right)\right]$, $\Theta_{\sigma{\sigma^{\prime}}}(\mathbf{r})\equiv F_{\sigma}(|\mathbf{r}-\mathbf{r}_{2}|)F_{\sigma^{\prime}}(|\mathbf{r}-\mathbf{r}_{1}|)$, by definition of the zero-mode at vortex 1 we obtain

[TABLE]

Compared with $p$-wave pairing Biswas (2013), under Majorana exchange $t_{12}$ is now symmetric with respect to $\theta_{1}-\theta_{2}$, but still anti-symmetric due to the spin degree of freedom.

It was pointed out by Fujimoto Fujimoto (2008) that the existence of the Majorana zero-mode with SO coupling does not automatically guarantee their non-Abelian statistics. Majorana interchange $1\leftrightarrow 2$ gives $t_{12}=-t_{21}$, and the U(1) gauge transformation $\theta_{i}$ of the Majorana fermion has the important property that when $\theta_{i}$ changes from $0\to 2\pi$, the Majorana changes sign, $\hat{\gamma}_{0,i}\to-\hat{\gamma}_{0,i}$. Braiding of vortices $i$ and $j$ changes the superfluid phase at one vortex by $2\pi$ amounting to $\hat{\gamma}_{0,i}\to\hat{\gamma}_{0,j}$, $\hat{\gamma}_{0,j}\to-\hat{\gamma}_{0,i}$, and it was shown by Ivanov Ivanov (2001) that this property together with quantisation of $\ell_{i}$ gives rise to non-Abelian statistics.

Conclusions – We have derived analytically the Majorana zero-mode in a topological Fermi superfluid with $s$-wave pairing, two-dimensional spin-orbit coupling, and a Zeeman field. We find an exponentially localised Majorana zero-mode at the vortex core only in the topologically non-trivial regime of the superfluid. We find that in the $s$-wave Fermi superfluid the Majorana fermions obey non-Abelian exchange statistics, and the energy splitting due to Majorana hybridisation is determined by both spin sectors. Knowing the Majorana zero-mode analytically paves the way for studies of quantum many-body correlations and simulation of topological quantum matter in lattices of Majorana vortices in a new experimentally well-controlled setup.

Appendix A Solution for the Majorana zero-mode

A.1 Formulation as a first-order system

Introducing the first derivatives $x=F_{\uparrow}^{\prime}$ and $y=F_{\downarrow}^{\prime}$ as auxiliary variables, we can trivially regroup the system (14) for $x^{\prime}$ and $y^{\prime}$ in terms of $x,y,F_{\uparrow},F_{\downarrow}$. This gives the equivalent matrix equation

[TABLE]

where

[TABLE]

where $\Gamma_{\sigma}\equiv-\frac{f^{4}}{\lambda^{2}}+\frac{f^{2}}{\lambda r}+\frac{m_{\sigma}^{2}}{r^{2}}-\frac{2m}{\hbar^{2}}\mu_{\sigma}$. This is a set of linear equations because $\mathcal{M}(r)$ does not depend on the components of $\mathbf{u}(r)$, but the coefficient matrix $\mathcal{M}(r)$ is non-constant depending on $r$.

The point $r=0$ is an irregular singular point of Eq. (19) because the coefficient matrix $\mathcal{M}(r)$ has a pole of order 2 at the origin.

Generally, if $\mathcal{M}(r)\mathcal{M}(r^{\prime})=\mathcal{M}(r^{\prime})\mathcal{M}(r)\,\,\forall r,r^{\prime}$, then Eq. (19) can be easly solved in terms of the matrix exponential

[TABLE]

where $\mathbf{u}(r_{0})$ is a 4-component constant vector. However, the commutation property does not hold here. There is no general closed-form solution for differential equations of the form of Eq. (19) where $\mathcal{M}(r)$ does not satisfy the commutation property. The Magnus series provides systematically an exact solution in terms of an infinite series of nested commutators:

[TABLE]

where

[TABLE]

but this approach suffers from the irregular singular point at the origin as well.

Near $r=0$, we can always write

[TABLE]

where $\mathcal{M}_{\nu}$ are constant $4\times 4$ matrices holomorphic at $r=0$ for which the series converges component-wise in a neighbourhood of the origin. Here $g=2$, and $\mathcal{M}_{\nu}\neq 0$ for only $\nu=0,1,2$.

A.2 Solution for the eigenvalue equation (14)

Generalising the vector $\mathbf{u}$ into the matrix $X$, and transforming $x=1/r$ changes Eq. (24) into

[TABLE]

where $q=g-2$ and

[TABLE]

with the only non-zero matrices being

[TABLE]

The matrix $\mathcal{M}(x)$ is holomorphic at $x=\infty$ meaning that there exists a convergent expansion of the form (26) for sufficiently large $x_{0}$ such that $|x|>x_{0}$. If the integer $q+1>0$, the singular point is irregular, and if $q+1=0$, the singular point is regular. For us $q=0$. Our goal is to reduce Eq. (25) (where $q=0$) through formal transformations into a system with $q=-1$, which corresponds to a regular singular point at the origin $x=\infty$, and can be solved in terms of more standard Fröbenius series ansatz methods.

To form our starting point, we transform $\mathcal{M}_{0}$ into the Jordan canonical form (JCF) (we define the JCF as having the 1’s on the superdiagonal). The matrix $\mathcal{M}_{0}$ is brought to the JCF by the non-singular constant similarity transformation $\mathcal{P}^{(0)}$,

[TABLE]

where for definiteness we have fixed $\ell=+1$. $\mathcal{A}_{0}$ is now in the canonical Jordan block diagonal form

[TABLE]

where $H_{j}$ ($j=1,2,3$) are shifting matrices where $H_{1}$ is of dimension two and $H_{2}$, $H_{3}$ are of dimension 1. The same transformation is applied to all the matrices, defining the new starting point $\mathcal{A}={\mathcal{P}^{(0)}}^{-1}(-\mathcal{M})\mathcal{P}^{(0)}$, $X\equiv\mathcal{P}^{(0)}Y$ such that

[TABLE]

where $q=0$, and

[TABLE]

A.2.1 Prepare for a shearing transformation that simplifies the problem

The transformation $Y(x)=\mathcal{K}(x)Z(x)$, where the matrix $\mathcal{K}(x)$ is holomorphic and has a non-vanishing determinant at $x=\infty$ changes Eq. (30) into

[TABLE]

with $q\geq 0$ and

[TABLE]

We seek series solutions $\mathcal{K}(x)=\sum_{\nu=0}^{\infty}\mathcal{K}_{\nu}x^{-\nu}$, $\mathcal{B}(x)=\sum_{\nu=0}^{\infty}\mathcal{B}_{\nu}x^{-\nu}$, where we recall $\mathcal{A}(x)=\sum_{\nu=0}^{\infty}\mathcal{A}_{\nu}x^{-\nu}$. We set

[TABLE]

Using Eq. (34), substitution of the series ansätze into Eq. (33) and comparison of like coefficients results in the recursion relation

[TABLE]

where the last term in Eq. (35b) is absent for $\nu-q-1<0$, that is, for $\nu<1$. The recursion relation is of the form

[TABLE]

where $K_{\nu}=\sum_{s=1}^{\nu-1}\mathcal{K}_{s}\mathcal{B}_{\nu-s}-\sum_{s=0}^{\nu-1}\mathcal{A}_{\nu-s}\mathcal{K}_{s}-(\nu-q-1)\mathcal{K}_{\nu-q-1}$ depends only on the $\mathcal{K}_{j}$, $\mathcal{B}_{j}$ with $j<\nu$. If all the eigenvalues of $\mathcal{A}_{0}$ are distinct, then $\mathcal{B}(x)$ will be diagonal and the problem is uncoupled and easily solved. If at least two eigenvalues are distinct, we can take all $\mathcal{B}_{\nu}$ ($\nu>0$) zero or block-diagonal. However, this is not possible here because $\mathcal{A}_{0}$ has only one distinct eigenvalue, zero.

Instead, we partition each Eq. (36) into blocks of the same order as the Jordan blocks $H_{j}$ for $\mathcal{A}_{0}$, and call these blocks $\mathcal{K}_{\nu}^{jk}$ with $j,k=1,2,\ldots,s$. For us $s=3$. Then each relation (36) corresponds to $s^{2}=9$ relations

[TABLE]

It can be proven that the equation $AX-XB=0$ possesses solutions other than $X=0$ if and only if $A$ and $B$ have at least one common eigenvalue. Since this result guarantees the existence of non-trivial solutions to the corresponding homogeneous equations of Eq. (37), each equation of Eq. (37) can be soluble only if the matrices $\mathcal{B}_{\nu}^{jk}$ satisfy some restrictive condition, explained below. Consider the auxiliary equation

[TABLE]

where $H$ and $K$ are shifting matrices of orders $h$ and $k$ respectively, and $M$ is a $h\times k$ matrix whose first $h-1$ rows are given constant vectors while the entries $\alpha_{1},\alpha_{2},\ldots,\alpha_{k}$ of the last row are variables. It can be proven that the numbers $\alpha_{1},\alpha_{2},\ldots,\alpha_{k}$ can be determined uniquely in such a way that Eq. (38) is solvable for the $h\times k$ matrix $X$. We now apply this result to Eq. (37). Without loss of generality we take all rows but the last of $\mathcal{B}_{\nu}^{jk}$ to be zero. Then, the series

[TABLE]

is determined by solving the recursion relations (37) successively for $\nu=1,2,\ldots$. While it will in general be divergent, it can be proven that it is the asymptotic expansion of some holomorphic matrix function $\mathcal{K}(x)$ for sufficiently large $x_{0}$ such that $|x|>x_{0}$.

The point of the transformation $\mathcal{K}$ is to obtain a differential equation (32) such that the matrices $\mathcal{B}$, by construction, satisfy the following three properties: (i) $\mathcal{B}(x)=\sum_{\nu=0}^{\infty}\mathcal{B}_{\nu}x^{-\nu}$; (ii) $\mathcal{B}_{0}=H_{1}\oplus H_{2}\oplus\cdots\oplus H_{s}$ (the $H_{j}$ are shifting matrices); and, most importantly, (iii) that the only non-zero entries in $\mathcal{B}_{\nu}$ for $\nu>0$ occur in the rows corresponding to the last rows of the Jordan blocks $H_{k}$ ($k=1,2,\ldots,s$) in the representation of $\mathcal{A}_{0}$. The transformation $\mathcal{K}$ induces zeros into the matrices preparing them for a shearing transformation. The result is an asymptotic series solution valid at $x\to\infty$ such that $\mathcal{B}(x)=\mathcal{B}_{0}+\mathcal{B}_{1}x^{-1}+\mathcal{B}_{2}x^{-2}+\mathcal{B}_{3}x^{-3}+\ldots$, where

[TABLE]

We work until $\nu=\nu_{\mathrm{max}}$. The final series solutions will then be exact upto $\nu=\nu_{\mathrm{max}}-1$ for $F_{\uparrow}(r)$ and $\nu=\nu_{\mathrm{max}}-2$ for $F_{\downarrow}(r)$. For spin-$\uparrow$ the derivative (component 2) will agree with a direct derivative of component 1 upto $\nu=\nu_{\mathrm{max}}-2$ and for spin-$\downarrow$ the derivative (component 4) will agree with a direct derivative of component 3 upto $\nu=\nu_{\mathrm{max}}-3$.

A.2.2 Simplify the problem by a shearing transformation

Having obtained the matrix $\mathcal{B}$, Eq. (32) is ready for the shearing transformation

[TABLE]

with a temporarily unknown positive parameter $\xi$, which takes Eq. (32) into

[TABLE]

where the matrix $\mathcal{C}(x)=S^{-1}(x)\mathcal{B}(x)S(x)-x^{-q}S^{-1}(x)S^{\prime}(x)$. The matrix elements $c_{jk}$ ($j,k=1,\ldots,4$, $n=4$) of $\mathcal{C}$ read

[TABLE]

Here $\delta_{jk}$ is the Kronecker delta, and $b_{jk}$ are the matrix elements of $\mathcal{B}$. The parameter $\xi$ must be chosen appropriately to induce, where possible, non-zero elements below the main diagonal into the leading order matrix $\mathcal{C}_{0}$. Above the main diagonal it is equal to $\mathcal{B}_{0}$.

The parameter $\xi$ must be chosen judiciously as follows. Any $b_{jk}\neq 0$ is of the form

[TABLE]

where $b_{jk\nu}\neq 0$ and each positive integer $\alpha_{jk}\geq 1$ except for the special elements with a 1 from a shifting matrix in the Jordan decomposition for $\mathcal{B}_{0}$ for which $\alpha_{jk}=0$. Since we could not diagonalise $\mathcal{B}_{0}$ and instead obtained a Jordan matrix, at least one such special element is present. Before the shearing transformation ($\xi=0$) the special elements have the lowest $\alpha_{jk}$, and the purpose of the transformation is to add non-zero elements on or below the diagonal to $\mathcal{B}_{0}$ by a suitable choice of $\xi$. After the shearing transformation the expansions of non-zero off-diagonal elements $c_{jk}$ ($\delta_{jk}=0$) begin with the power $-\alpha_{jk}+\xi(j-k)$. The special elements on the superdiagonal begin with the power $x^{-\xi}$. There exists a smallest rational $\xi_{0}=q^{(1)}/p^{(1)}>0$ with $q^{(1)},p^{(1)}$ coprime, for which the special elements have the same leading power as an element below the main diagonal, $\xi_{0}=\alpha_{jk}-\xi_{0}(j-k)$ for some $j,k<j$. We take $\xi=\xi_{0}$ in the shearing transformation; the result of the above process is $\xi=q^{(1)}/p^{(1)}$ with $q^{(1)}=1$, $p^{(1)}=2$ coprime. Fractional powers thus unavoidably appear in the shearing transformation (41).

In descending powers of $x^{-1/p^{(1)}}$, the matrix $\mathcal{C}$ will begin with the power $x^{-q^{(1)}/p^{(1)}}$. The matrix $x^{\xi}\mathcal{C}$ as $x\to\infty$ has at least one non-zero entry on or below the main diagonal and above the main diagonal it is equal to $\mathcal{B}_{0}$:

[TABLE]

We remove the fractional powers by multiplying both sides by $x^{\xi}$ and introducing the new independent variable $x_{1}$ such that

[TABLE]

to obtain from Eq. (42) the differential equation

[TABLE]

where

[TABLE]

The branch of the multi-valued function $x^{p^{(1)}}$ can be chosen freely.

While new elements have appeared on and below the main diagonal in $\mathcal{D}_{0}$, the upper triangular part above the main diagonal, in particular the superdiagonal, of $\mathcal{D}_{0}$ is equal to that of $\mathcal{B}_{0}$ by construction. This completes the purpose of the shearing transformation.

If $h<0$, the problem would be solved because then either the singular point would be regular ($h=-1$), or there would not be any singular point. If $h\geq 0$ and $\mathcal{D}_{0}$ has at least two distinct eigenvalues, then the problem reduces to a set of similar problems of lower order. However, the eigenvalues of $\mathcal{D}_{0}$ are all equal to zero, and $\mathcal{D}_{0}$ is nilpotent. It can be proven Wasow (2002) that if the process outlined above is repeatedly carried out, eventually we obtain a nilpotent $\mathcal{D}_{0}$ that is itself a shifting matrix i.e. $s=1$.

A.2.3 Obtain an equation with a regular singular point

The entire process must be reapplied, starting from the JCF for $\mathcal{D}_{0}$. Another shearing transformation with the parameter $1/2$ is required. Repeating the process once more, the third shearing transformation comes with an integer-valued shearing parameter of $1$, which means that via such a chain of transformations we have reduced the coupled system (30),

[TABLE]

to the system

[TABLE]

where

[TABLE]

and $x_{1}=x_{2}^{2}/4$. We have reached our goal that we set below Eq. (25): Eq. (50) has only a regular singular point at infinity (or the vortex core in terms of $r$). This equation can be solved with more standard methods.

Through a sequence of standard double transformations that raise the eigenvalues of a Jordan block of $\mathcal{A}^{(3)}_{0}$ by an integer followed by a non-singular constant similarity transformation of $\mathcal{A}^{(3)}_{0}$ to JCF, we obtain

[TABLE]

such that no eigenvalues of the leading matrix $\mathcal{A}^{(6)}_{0}$ differ by a positive integer. In fact, it consists of two identical Jordan blocks: $\mathcal{A}^{(6)}_{0}=\begin{pmatrix}13&1\\ 0&13\\ \end{pmatrix}\oplus\begin{pmatrix}13&1\\ 0&13\\ \end{pmatrix}$.

Let us now transform back to radial coordinates with $r_{2}=1/x_{2}$ (here $x_{1}=x_{2}^{2}/4$) to obtain from Eq. (52) the equation

[TABLE]

where $\mathcal{A}^{(7)}(r_{2})=-\mathcal{A}^{(6)}(1/r_{2})$ and $\mathcal{A}^{(7)}(0)$ is holomorphic.

Since $\mathcal{A}^{(7)}(r_{2})$ is holomorphic at $r_{2}=0$ and since no two eigenvalues of $\mathcal{A}^{(7)}_{0}$ differ by a positive integer, Eq. (53) has a fundamental matrix solution of the form

[TABLE]

where $\mathcal{K}^{(7)}(0)$ is holomorphic. Its power series representation $\mathcal{K}^{(7)}=\sum_{\nu=0}^{\infty}\mathcal{K}^{(7)}_{\nu}r_{2}^{\nu}$ can be calculated by rational operations from the coefficients $\mathcal{A}^{(7)}_{\nu}$ in the series $\mathcal{A}^{(7)}=\sum_{\nu=0}^{\infty}\mathcal{A}^{(7)}_{\nu}r_{2}^{\nu}$.

The matrix $\mathcal{K}^{(7)}(r_{2})$ in the transformation $Y^{(7)}=\mathcal{K}^{(7)}(r_{2})\Upsilon$, which results in the equation

[TABLE]

is calculated using the recursion relation (35) with $q=-1$ and $r=1/x$ (also changing the sign of the last term in Eq. (56b)), viz.

[TABLE]

where

[TABLE]

Given that the matrix $\mathcal{A}^{(7)}_{0}$ in the convergent expansion $\mathcal{A}^{(7)}=\sum_{\nu=0}^{\infty}\mathcal{A}^{(7)}_{\nu}r_{2}^{\nu}$ has no eigenvalues that differ from each other by positive integers, there exists a formal convergent series $\mathcal{K}^{(7)}(r_{2})$ with $\mathcal{K}^{(7)}_{0}=I$ such that the formal transformation $Y^{(7)}=\mathcal{K}^{(7)}(r_{2})\Upsilon$ reduces the differential equation (53) to the form of Eq. (55) such that $\mathcal{B}^{(7)}_{\nu}=0$ for all $\nu>0$. Therefore, the recursion relations (56) can be solved for $\mathcal{K}^{(7)}(r_{2})$ such that $\mathcal{B}^{(7)}_{\nu}=0$ for $\nu>0$. The solution for $Y^{(7)}(r_{2})$ is then given by Eq. (54).

Since $\left[\frac{\mathcal{B}^{(7)}_{0}}{r_{1}},\frac{\mathcal{B}^{(7)}_{0}}{r_{2}}\right]=0$, the general solution to the equation

[TABLE]

is given by the matrix exponential (21),

[TABLE]

Here a set of four linearly independent vectors forms the fundamental matrix solution $\Upsilon(r_{2})$, and $\Upsilon(r_{0})$ is the fundamental matrix solution at $r_{2}=r_{0}$. We can set the auxiliary parameter $r_{0}=1$ and choose $\Upsilon(r_{0})$ in such a way as to choose the boundary conditions for $F_{\sigma}(r)$ at $r=0$. The general (vector) solution is then given by linear combinations of the columns.

The logarithms of the fundamental solution matrix end up in columns 2 and 4 making them either divergent or zero at the origin $r=0$. The columns 1 and 3, in contrast, can be chosen to remain well-behaved and finite at the vortex core $r=0$. We do not include the columns 2 and 4 in what follows, that is, we take a linear combination only of columns 1 and 3.

Performing all the transformations carried out in an inverse order, we can compute the matrix solution $X$ for the original problem (14) from knowing $\Upsilon(r_{2})$. A judicious choice for $\Upsilon(1)$ (with $r_{0}=1$) then gives the explicit solution shown in Eq. (15) in the main text.