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

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

$$\begin{eqnarray}&&H({\bf{r}}){{\bf{h}}}_{\nu }={E}_{\nu }{{\bf{h}}}_{\nu }\end{eqnarray} \tag{ 5 }$$

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({\bf{r}},t)$ possesses a particle-hole symmetry (PHS) defined by $\{{{\rm{\Xi }}}_{{\rm{s}}-\mathrm{wave}},H({\bf{r}},t)\}\,=0$ with ${{\rm{\Xi }}}_{{\rm{s}}-\mathrm{wave}}={{\rm{e}}}^{{\rm{i}}\theta }{\tau }_{{\rm{x}}}\otimes {1}_{\sigma }K$, where K is the complex conjugation operator (in momentum space, K additionally flips the sign of the momentum), ${\tau }_{{\rm{y}}}$ (${\sigma }_{{\rm{z}}}$) is the Pauli matrix in particle-hole (spin) space, and $\theta \in {\mathbb{R}}$. As a result of the particle-hole symmetry, if ${{\bf{h}}}_{\nu }={\left(\begin{array}{cccc}{u}_{\nu ,\uparrow } & {u}_{\nu ,\downarrow } & {v}_{\nu ,\uparrow } & {v}_{\nu ,\downarrow }\end{array}\right)}^{{T}}$ is a solution of equation (5) with energy E_ν, then ${{\rm{\Xi }}}_{{\rm{s}}-\mathrm{wave}}{{\bf{h}}}_{\nu }={\mathrm{ie}}^{{\rm{i}}\theta }{\left(\begin{array}{cccc}-{v}_{\nu ,\uparrow }^{* } & {v}_{\nu ,\downarrow }^{* } & {u}_{\nu ,\uparrow }^{* } & -{u}_{\nu ,\downarrow }^{* }\end{array}\right)}^{{\rm{T}}}$ is a solution with energy $-{E}_{\nu }$.

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

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

$$\begin{eqnarray}&&{u}_{0,\sigma }={{\rm{e}}}^{-{\rm{i}}\theta }{v}_{0,\sigma }^{* },\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&{\hat{\gamma }}_{0}^{\dagger }=\int {\rm{d}}{\bf{r}}\displaystyle \sum _{\sigma }\left[{u}_{0,\sigma }({\bf{r}}){\hat{\psi }}_{\sigma }^{\dagger }({\bf{r}})+{v}_{0,\sigma }({\bf{r}}){\hat{\psi }}_{\sigma }({\bf{r}})\right]\end{eqnarray} \tag{ 7 }$$

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

Let us temporarily assume that time-reversal symmetry is respected, i.e. ${\rm{\Delta }}({\bf{r}})\in {\mathbb{R}}$ and $h=0$. If ${{\bf{h}}}_{\nu }={\left(\begin{array}{cccc}{u}_{\nu ,\uparrow } & {u}_{\nu ,\downarrow } & {v}_{\nu ,\uparrow } & {v}_{\nu ,\downarrow }\end{array}\right)}^{{\rm{T}}}$ is a solution of equation (5) with energy E_ν, then $T{{\bf{h}}}_{\nu }$ is a solution with the same energy E_ν. This is the pair of the zero-energy Majorana mode. With time-reversal symmetry, the Majorana probability densities are identically overlapping in space: $| {{\bf{h}}}_{0}{| }^{2}=| T{{\bf{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. We note that the focus here is on an infinite system where we ignore the boundary physics. Strictly speaking, in the case considered here, there can only be one or zero Majorana zero-modes depending on the topological properties of the superfluid. In experiments the Fermi superfluid is always trapped, where in the topological regime there is a localised state at the edge and at the vortex core which can be identified as Majorana fermions in the thermodynamic limit. This case has been studied self-consistently numerically for harmonic trapping in [60].

A broken time-reversal symmetry T, a broken spin-rotational symmetry (due to spin–orbit coupling), and an unbroken PHS ${{\rm{\Xi }}}_{{\rm{s}}-\mathrm{wave}}$ result in the BdG Hamiltonian $H({\bf{r}})$ belonging to the symmetry class D in the Altland-Zirnbauer classification [61]. According to [62], 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\mathrm{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 [47–49]. We calculate the energy splitting in equation (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 [48].

We now include a vortex in the order parameter

$$\begin{eqnarray}&&{\rm{\Delta }}({\bf{r}})={\rm{\Delta }}(r){{\rm{e}}}^{{\rm{i}}{\ell }\varphi },\end{eqnarray} \tag{ 8 }$$

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

For convenience, let us first apply the unitary transformation ${{ \mathcal U }}_{{\rm{s}}}$, which transforms the basis as ${{ \mathcal U }}_{{\rm{s}}}{\left(\begin{array}{cccc}{u}_{0,\uparrow } & {u}_{0,\downarrow } & {v}_{0,\uparrow } & {v}_{0,\downarrow }\end{array}\right)}^{{T}}={\left(\begin{array}{cccc}{u}_{0,\downarrow } & {v}_{0,\downarrow } & {v}_{0,\uparrow } & {u}_{0,\uparrow }\end{array}\right)}^{{T}}$, and shuffles the zeros in the rows and columns of the Hamiltonian $H({\bf{r}})$ (equation (2)) to give

$$\begin{eqnarray}{H}^{({{ \mathcal U }}_{{\rm{s}}})}({\bf{r}})\equiv {{ \mathcal U }}_{{\rm{s}}}H({\bf{r}}){{ \mathcal U }}_{{\rm{s}}}^{-1}=\left(\begin{array}{cc}{D}_{\downarrow } & M\\ {M}^{\dagger } & -{D}_{\uparrow }\end{array}\right),\end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray}M=\left(\begin{array}{cc}-{\rm{\Delta }}({\bf{r}}) & -{\rm{i}}\lambda ({\partial }_{y}-{\rm{i}}{\partial }_{x})\\ -{\rm{i}}\lambda ({\partial }_{y}+{\rm{i}}{\partial }_{x}) & {{\rm{\Delta }}}^{* }({\bf{r}})\end{array}\right).\end{eqnarray} \tag{ 10 }$$

Setting ${D}_{\sigma }=0$ decouples the spins, and the Hamiltonian reduces to a 2 × 2 BdG matrix corresponding to the Fu-Kane model at neutrality, whose Majorana modes have been considered [63]. 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 \ne 0$ [45].

However, here we consider the full 4 × 4 structure of equation (5):

$$\begin{eqnarray}&&{\hat{K}}_{\uparrow \uparrow }{u}_{\nu ,\uparrow }+{\hat{K}}_{\uparrow \downarrow }{u}_{\nu ,\downarrow }+{\rm{\Delta }}({\bf{r}}){v}_{\nu ,\downarrow }={E}_{\nu }{u}_{\nu ,\uparrow },\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&-{\hat{K}}_{\uparrow \uparrow }{v}_{\nu ,\uparrow }+{\hat{K}}_{\downarrow \uparrow }{v}_{\nu ,\downarrow }-{{\rm{\Delta }}}^{* }({\bf{r}}){u}_{\nu ,\downarrow }={E}_{\nu }{v}_{\nu ,\uparrow },\end{eqnarray} \tag{ 11b }$$

$$\begin{eqnarray}&&{\hat{K}}_{\downarrow \uparrow }{u}_{\nu ,\uparrow }+{\hat{K}}_{\downarrow \downarrow }{u}_{\nu ,\downarrow }-{\rm{\Delta }}({\bf{r}}){v}_{\nu ,\uparrow }={E}_{\nu }{u}_{\nu ,\downarrow },\end{eqnarray} \tag{ 11c }$$

$$\begin{eqnarray}&&{\hat{K}}_{\uparrow \downarrow }{v}_{\nu ,\uparrow }-{\hat{K}}_{\downarrow \downarrow }{v}_{\nu ,\downarrow }+{{\rm{\Delta }}}^{* }({\bf{r}}){u}_{\nu ,\uparrow }={E}_{\nu }{v}_{\nu ,\downarrow }.\end{eqnarray} \tag{ 11d }$$

When looking for zero-modes, 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 θ, 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

$$\begin{eqnarray}&&{\hat{K}}_{\uparrow \uparrow }{u}_{0,\uparrow }+{\hat{K}}_{\uparrow \downarrow }{u}_{0,\downarrow }+{\rm{\Delta }}({\bf{r}}){u}_{0,\downarrow }^{* }=0,\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{\hat{K}}_{\downarrow \uparrow }{u}_{0,\uparrow }+{\hat{K}}_{\downarrow \downarrow }{u}_{0,\downarrow }-{\rm{\Delta }}({\bf{r}}){u}_{0,\uparrow }^{* }=0,\end{eqnarray} \tag{ 12b }$$

a coupled 2 × 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 }({\bf{r}})=\lambda {{\rm{e}}}^{-{\rm{i}}\varphi }\left[{\partial }_{r}-({\rm{i}}/r){\partial }_{\varphi }\right]$. Observing the azimuthal symmetry, let us try a separable ansatz with angular momentum eigenstates

$$\begin{eqnarray}&&{u}_{0,\sigma }={{\rm{e}}}^{{\rm{i}}{m}_{\sigma }\varphi }{F}_{\sigma }(r){{\rm{e}}}^{-\tfrac{{\rm{\Delta }}}{\lambda }r},\end{eqnarray} \tag{ 13 }$$

with the assumption that ${\rm{\Delta }}$ has no r-dependence. Physically, this approximation treats the vortices as points with only a phase profile. Substitution into 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

$$\begin{eqnarray}&&-{F}_{\uparrow }^{{\prime\prime} }-\left(\displaystyle \frac{1}{r}-2\displaystyle \frac{{\rm{\Delta }}}{\lambda }\right){F}_{\uparrow }^{{\prime} }+\left(-\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}+\displaystyle \frac{{\rm{\Delta }}}{\lambda r}\right){F}_{\uparrow }+\displaystyle \frac{{m}_{\uparrow }^{2}{F}_{\uparrow }}{{r}^{2}}+\displaystyle \frac{2m\lambda }{{{\rm{\hslash }}}^{2}}\left({F}_{\downarrow }^{{\prime} }+\displaystyle \frac{{m}_{\downarrow }{F}_{\downarrow }}{r}\right)=\displaystyle \frac{2m}{{{\rm{\hslash }}}^{2}}{\mu }_{\uparrow }{F}_{\uparrow },\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&-{F}_{\downarrow }^{{\prime\prime} }-\left(\displaystyle \frac{1}{r}-2\displaystyle \frac{{\rm{\Delta }}}{\lambda }\right){F}_{\downarrow }^{{\prime} }+\left(-\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}+\displaystyle \frac{{\rm{\Delta }}}{\lambda r}\right){F}_{\downarrow }+\displaystyle \frac{{m}_{\downarrow }^{2}{F}_{\downarrow }}{{r}^{2}}+\displaystyle \frac{2m\lambda }{{{\rm{\hslash }}}^{2}}\left(-{F}_{\uparrow }^{{\prime} }+\displaystyle \frac{{m}_{\uparrow }{F}_{\uparrow }}{r}\right)=\displaystyle \frac{2m}{{{\rm{\hslash }}}^{2}}{\mu }_{\downarrow }{F}_{\downarrow }.\end{eqnarray} \tag{ 14b }$$

The prime denotes differentiation with respect to r. Below, we solve analytically the coupled system (14) for ${F}_{\uparrow }(r)$ and ${F}_{\downarrow }(r)$.

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

The details for solving system (14) are shown in appendix. 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

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\uparrow }(r) & = & {F}_{\uparrow }(0)+\displaystyle \frac{{\rm{\Delta }}{F}_{\uparrow }(0)}{\lambda }r+\left\{\displaystyle \frac{\lambda m}{{{\rm{\hslash }}}^{2}}{F}_{\downarrow }^{{\prime} }(0)+\left[\displaystyle \frac{{{\rm{\Delta }}}^{2}}{2{\lambda }^{2}}-\displaystyle \frac{(\mu +h)m}{2{{\rm{\hslash }}}^{2}}\right]{F}_{\uparrow }(0)\right\}{r}^{2}\\ & & +\left\{\left[\displaystyle \frac{{{\rm{\Delta }}}^{3}}{6{\lambda }^{3}}-\displaystyle \frac{{\rm{\Delta }}(\mu +h)m}{2\lambda {{\rm{\hslash }}}^{2}}-\displaystyle \frac{4{\rm{\Delta }}\lambda {m}^{2}}{9{{\rm{\hslash }}}^{4}}\right]{F}_{\uparrow }(0)+\displaystyle \frac{11{\rm{\Delta }}m}{9{{\rm{\hslash }}}^{2}}{F}_{\downarrow }^{{\prime} }(0)\right\}{r}^{3}+{ \mathcal O }({r}^{4}),\end{array}\end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\downarrow }(r) & = & {F}_{\downarrow }^{{\prime} }(0)r+{\rm{\Delta }}\left(\displaystyle \frac{{F}_{\downarrow }^{{\prime} }(0)}{\lambda }-\displaystyle \frac{2{F}_{\uparrow }(0)m}{3{{\rm{\hslash }}}^{2}}\right){r}^{2}+\left\{\displaystyle \frac{2{{\rm{\Delta }}}^{2}{{\rm{\hslash }}}^{4}+{\lambda }^{2}m{{\rm{\hslash }}}^{2}(h-\mu )-2{\lambda }^{4}{m}^{2}}{4{\lambda }^{2}{{\rm{\hslash }}}^{4}}{F}_{\downarrow }^{{\prime} }(0)\right.\\ & & \left.+\displaystyle \frac{-8{{\rm{\Delta }}}^{2}{{\rm{\hslash }}}^{2}m+3h{\lambda }^{2}{m}^{2}\,+\,3{\lambda }^{2}\mu {m}^{2}}{12\lambda {{\rm{\hslash }}}^{4}}{F}_{\uparrow }(0)\right\}{r}^{3}+{ \mathcal O }({r}^{4}).\end{array}\end{eqnarray} \tag{ 15b }$$

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⁰. 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 [60, 65]. In principle 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}_{\downarrow }^{{\prime} }(0)$ in the solution (15).

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

$$\begin{eqnarray}&&\lambda {F}_{\downarrow }^{{\prime} }(r)+{\rm{\Delta }}{F}_{\downarrow }(r)-{F}_{\uparrow }(r)(\bar{\mu }+h)=0,\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&\lambda {F}_{\uparrow }^{{\prime} }(r)+{\rm{\Delta }}{F}_{\uparrow }(r)+{F}_{\downarrow }(r)(\bar{\mu }-h)=0.\end{eqnarray} \tag{ 16b }$$

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

$$\begin{eqnarray}&&{F}_{\uparrow }(r)=C\,{\rm{\exp }}\left[-\displaystyle \frac{{\rm{\Delta }}+\sqrt{{h}^{2}-{\bar{\mu }}^{2}}}{\lambda }r\right]=\displaystyle \frac{\bar{\mu }-h}{\sqrt{{h}^{2}-{\bar{\mu }}^{2}}}{F}_{\downarrow }(r),\end{eqnarray} \tag{ 17 }$$

where C is a constant. Even then, the solution (17) is normalisable only when $h\gt \sqrt{{\bar{\mu }}^{2}+| {\rm{\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, ${\rm{\Delta }}({\bf{r}})={\rm{\Delta }}$ [30], 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}_{\downarrow }^{{\prime} }(0)$ such that the zero-mode has the desired large-r asymptotics.

In principle, $\bar{\mu }$ and Δ must be obtained self-consistently in conjunction with solving the BdG eigenvalue equation. Considering a uniform vortex-free system, we can solve equations (3), (4), and (11) self-consistently to find that the superfluid is topologically non-trivial with the parameters $\lambda =5,h=10,{\rm{\Delta }}=-1.587,\bar{\mu }=-5.987$. Here all units are measured with respect to ${\rm{\hslash }}=m=1$. We have used the pair binding energy of ${E}_{{\rm{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 [30]. From equation (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}_{\downarrow }^{{\prime} }(0)$ that match with this asymptotic behaviour.

We retain terms upto ${ \mathcal O }({r}^{15})$, and use the analytical solution (15) to evaluate accurate initial conditions for a numerical integration of system (14). The truncated analytical series solution agrees with the numerical integration for $r\lesssim 0.5$ (figure 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.488\,885\lt {F}_{\uparrow }(0)\lt 3.488\,886$ and ${F}_{\downarrow }^{{\prime} }(0)=-12.518$. The full solution including azimuthal dependence is then given by equation (13).

Zoom In Zoom Out Reset image size

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 ${{\bf{r}}}_{1}$ and ${{\bf{r}}}_{2}$. Then, ${\rm{\Delta }}({\bf{r}})\approx {{\rm{\Delta }}}_{i}({\bf{r}})$ near vortex i ($i=1,2$), where ${{\rm{\Delta }}}_{i}({\bf{r}})$ is the single-vortex pair function (8) for vortex i. Introduction of a single-vortex background phase ${\theta }_{i},{{\rm{\Delta }}}_{i}({\bf{r}})={\rm{\Delta }}{{\rm{e}}}^{{\rm{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 }({\bf{r}}-{{\bf{r}}}_{i})\to {{\rm{e}}}^{{\rm{i}}\tfrac{{\theta }_{i}}{2}}{u}_{0,\sigma }({\bf{r}}-{{\bf{r}}}_{i})$, leaving invariant system (12).

In the tight-binding approximation, the Hamiltonian H_Δ describing hopping between the two vortices, whose generalisation to lattices of Majorana vortices requires incorporating a ${{\mathbb{Z}}}_{2}$ symmetry (section 4.2), is given by ${H}_{{\rm{\Delta }}}\,={\rm{i}}{t}_{12}\,{\hat{\gamma }}_{\mathrm{0,1}}{\hat{\gamma }}_{\mathrm{0,2}}$, where ${\hat{\gamma }}_{0,i}$ is the Majorana mode (7) at vortex i, and ${t}_{12}=\langle {{\bf{g}}}_{0}^{(2)}| {H}^{({{ \mathcal U }}_{{\rm{s}}})}| {{\bf{g}}}_{0}^{(1)}\rangle =\langle {{\bf{g}}}_{0}^{(2)}\left|\left(\begin{array}{cc}{D}_{\downarrow } & M\\ {M}^{\dagger } & -{D}_{\uparrow }\end{array}\right)\right|{{\bf{g}}}_{0}^{(1)}\rangle$ is the overlap integral that gives the energy splitting. Here ${{\bf{g}}}_{0}^{(i)}({\bf{r}})={\left(\begin{array}{cccc}{u}_{0,\downarrow }({\bf{r}}-{{\bf{r}}}_{i}) & {u}_{0,\downarrow }^{* }({\bf{r}}-{{\bf{r}}}_{i}) & {u}_{0,\uparrow }^{* }({\bf{r}}-{{\bf{r}}}_{i}) & {u}_{0,\uparrow }({\bf{r}}-{{\bf{r}}}_{i})\end{array}\right)}^{{\rm{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 equation (10) when the order parameter is given by just a single vortex at position ${{\bf{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 equation (13), and introducing the definitions $G({\bf{r}})\equiv -{\rm{\Delta }}\exp \left[-\tfrac{{\rm{\Delta }}}{\lambda }\left(| {\bf{r}}-{{\bf{r}}}_{2}| +| {\bf{r}}-{{\bf{r}}}_{1}| \right)\right],{{\rm{\Theta }}}_{\sigma {\sigma }^{{\prime} }}({\bf{r}})\equiv {F}_{\sigma }(| {\bf{r}}-{{\bf{r}}}_{2}| ){F}_{{\sigma }^{{\prime} }}(| {\bf{r}}-{{\bf{r}}}_{1}| )$, by definition of the zero-mode at vortex 1 we obtain

$$\begin{eqnarray}&&{t}_{12}=4\cos \left(\displaystyle \frac{{\theta }_{1}-{\theta }_{2}}{2}\right)\int {\rm{d}}{\bf{r}}G({\bf{r}}){\sin }^{2}\left(\displaystyle \frac{{\varphi }_{1}-{\varphi }_{2}}{2}\right)\left\{{{\rm{\Theta }}}_{\downarrow \uparrow }({\bf{r}})-\,{{\rm{\Theta }}}_{\uparrow \downarrow }({\bf{r}})\right\}.\end{eqnarray} \tag{ 18 }$$

Equation (18) has both oscillatory behaviour and exponential decay in system size (figure 2). The latter is captured by $G({\bf{r}})$ and the former is captured by ${{\rm{\Theta }}}_{\sigma {\sigma }^{{\prime} }}$, which has strong and non-trivial oscillatory dependence on the system parameters as indicated explicitly in system (15) for F_σ.

Zoom In Zoom Out Reset image size

Compared to systems with p-wave pairing [66], or the Fu–Kane model [63], in the s-wave Fermi superfluid considered here under Majorana exchange t₁₂ is symmetric with respect to ${\theta }_{1}-{\theta }_{2}$, but still anti-symmetric overall due to the spin degree of freedom: swopping the vortices amounts to ${{\rm{\Theta }}}_{\downarrow \uparrow }({\bf{r}})-{{\rm{\Theta }}}_{\uparrow \downarrow }({\bf{r}})\to {{\rm{\Theta }}}_{\uparrow \downarrow }({\bf{r}})-\,{{\rm{\Theta }}}_{\downarrow \uparrow }({\bf{r}})$. Anti-symmetry of t₁₂ reflects the non-Abelian statistics, and for Majorana vortices in p-wave systems or the Fu–Kane model, where ${t}_{12}\sim \sin \left(\tfrac{{\theta }_{1}-{\theta }_{2}}{2}\right)$, it arises from anti-symmetry with respect to ${\theta }_{1}-{\theta }_{2}$.

It was pointed out by Fujimoto [1] that the existence of the Majorana zero-mode with SO coupling does not automatically guarantee their non-Abelian statistics. Here, 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π while leaving the phase at the other vortex invariant, 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 [13] that this property together with quantisation of ℓ_i gives rise to non-Abelian statistics.

In other words, we can perform a local ${{\mathbb{Z}}}_{2}$ gauge transformation ${\hat{\gamma }}_{0,j}\to -{\hat{\gamma }}_{0,j}$ because the Majorana commutation algebra $\{{\hat{\gamma }}_{0,i},{\hat{\gamma }}_{0,j}\}=2{\delta }_{{ij}},{\hat{\gamma }}_{0,j}^{\dagger }={\hat{\gamma }}_{0,j}$ remains invariant. To account for the ${{\mathbb{Z}}}_{2}$ ambiguity when generalising to multiple Majorana vortices we introduce the ${{\mathbb{Z}}}_{2}$ gauge factors ${s}_{{ij}}={{\rm{e}}}^{{\rm{i}}{\phi }_{{ij}}}=\pm {\rm{i}}$ to the tight-binding Hamiltonian H_Δ

$$\begin{eqnarray}&&{H}_{{\rm{\Delta }}}=\displaystyle \sum _{\langle i,j\rangle }{s}_{{ij}}{t}_{{ij}}\,{\hat{\gamma }}_{0,i}{\hat{\gamma }}_{0,j},\end{eqnarray} \tag{ 19 }$$

where ${t}_{{ij}}=-{t}_{{ji}}$ and given in equation (18). We note that it is Hermitian, ${H}_{{\rm{\Delta }}}^{\dagger }={H}_{{\rm{\Delta }}}$. Due to the exponential decay of the Majorana zero-mode and the hopping amplitude t_ij (figure 2), we focus on nearest-neighbour hopping $\langle i,j\rangle$ only. Since the Majorana fermion corresponds to a zero-energy solution of BdG equations (11), the on-site energy in the tight-binding model is zero leaving only hopping terms. Hamiltonian (19) has been considered in the literature for vortex lattices in the Kitaev spin model on the honeycomb lattice with non-Abelian vortex excitations [14]. Lahtinen [67] extended [43], and mapped the Kitaev spin model to an effective lattice model of Majorana fermions for the vortex excitations described by Hamiltonian (19) by mapping the spin-$\tfrac{1}{2}$ particles to Majorana fermions.

The existence of the Majorana zero-mode (15) has important consequences for the collective dynamics of networks of Majorana vortices. In a vortex lattice, owing to the factors s_ij in equation (19), a closed loop can enclose a plaquette of the ${{\mathbb{Z}}}_{2}$ gauge flux, which is defined by the product of the factors s_ij along the trajectory. Using the same Hamiltonian H_Δ (but different t₁₂), Liu and Franz [48] showed that in the context of Majorana zero-modes in a vortex lattice in a two-dimensional p-wave superconductor and the Fu–Kane model characterising an s-wave superconductor coupled to a Dirac material through the proximity effect, the ${{\mathbb{Z}}}_{2}$ gauge flux through a general polygon formed by n vortices is given by

$$\begin{eqnarray}&&\displaystyle \sum _{\mathrm{polygon}}{\phi }_{{ij}}=\displaystyle \frac{\pi }{2}(n-2).\end{eqnarray} \tag{ 20 }$$

The same result can therefore be expected to hold also in the s-wave Fermi superfluid despite the different symmetry with respect to the U(1) gauge factors ${\theta }_{i}$. Originally, Grosfeld and Stern [2] derived rule (20) for Majorana zero-modes in the context of the Moore–Read $\nu =5/2$ fractional quantum Hall state [11]. The effective theory for the Moore–Read Pfaffian $\nu =5/2$ fractional quantum Hall state is analogous to a spinless ${p}_{{\rm{x}}}+{\rm{i}}{p}_{{\rm{y}}}$-wave superconductor of composite fermions [2]. It follows directly from equation (20) that there is a non-zero ${{\mathbb{Z}}}_{2}$ gauge flux of $\pi /2$ and π through an elementary triangular and square plaquette, respectively, of the Majorana lattice. A finite gauge flux that alternates between neighbouring plaquettes is typical of Chern systems [68] with the occupied bands in the gapped phases of Hamiltonian (19) being characterised by a non-zero Chern number.

It should be noted that for strict thermodynamic stability of the vortex lattice we need a gauge field ${\bf{A}}$, which in the superconductor context is a magnetic field [48], and in the superfluid context represents rotation [69–71]. The phase difference $({\theta }_{1}-{\theta }_{2})/2$ in equation (18) is then replaced by a gauge-invariant version ${\rm{\Delta }}\theta ={\int }_{{{\bf{r}}}_{1}}^{{{\bf{r}}}_{2}}\left(\tfrac{1}{2}{\rm{\nabla }}\theta -{\bf{A}}\right)\cdot {\rm{d}}{\bf{l}}$, where the integral is along the line segment joining the two vortices [63]. While these terms must be included in any calculation that aims to accurately study the band structure of many-Majorana physics in a thermodynamically stable vortex lattice, the conclusions drawn here are not sensitive to the thermodynamic stability.

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

$$\begin{eqnarray}&&{{\bf{u}}}^{{\prime} }(r)={ \mathcal M }(r){\bf{u}}(r),\end{eqnarray} \tag{ A1 }$$

where

$$\begin{eqnarray}&&{\bf{u}}(r)=\left(\begin{array}{c}{F}_{\uparrow }\\ x\\ {F}_{\downarrow }\\ y\end{array}\right),\end{eqnarray} \tag{ A2a }$$

$$\begin{eqnarray}{ \mathcal M }(r)=\left(\begin{array}{cccc}0 & 1 & 0 & 0\\ {{\rm{\Gamma }}}_{\uparrow } & -\left(\tfrac{1}{r}-2\tfrac{{\rm{\Delta }}}{\lambda }\right) & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}}\tfrac{{m}_{\downarrow }}{r} & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}}\\ 0 & 0 & 0 & 1\\ \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}}\tfrac{{m}_{\uparrow }}{r} & -\tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & {{\rm{\Gamma }}}_{\downarrow } & -\left(\tfrac{1}{r}-2\tfrac{{\rm{\Delta }}}{\lambda }\right)\end{array}\right),\end{eqnarray} \tag{ A2b }$$

where ${{\rm{\Gamma }}}_{\sigma }\equiv -\tfrac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}+\tfrac{{\rm{\Delta }}}{\lambda r}+\tfrac{{m}_{\sigma }^{2}}{{r}^{2}}-\tfrac{2m}{{{\rm{\hslash }}}^{2}}{\mu }_{\sigma }$. This is a set of linear equations because ${ \mathcal M }(r)$ does not depend on the components of ${\bf{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 equation (A1) 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 equation (A1) can be easily solved in terms of the matrix exponential

$$\begin{eqnarray}&&{\bf{u}}(r)={\bf{u}}({r}_{0}){{\rm{e}}}^{{\int }_{{r}_{0}}^{r}{ \mathcal M }(s){\rm{d}}s},\end{eqnarray} \tag{ A3 }$$

where ${\bf{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 equation (A1) 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:

$$\begin{eqnarray}&&{\bf{u}}(r)={\bf{u}}({r}_{0}){{\rm{e}}}^{{\displaystyle \sum }_{k=1}^{\infty }{{\rm{\Omega }}}_{k}(s)},\end{eqnarray} \tag{ A4 }$$

where

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{1}(s)={\int }_{{r}_{0}}^{r}{ \mathcal M }({s}_{1}){\rm{d}}{s}_{1},\end{eqnarray} \tag{ A5a }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{2}(s)=\displaystyle \frac{1}{2}{\int }_{{r}_{0}}^{r}\,{\rm{d}}{s}_{1}{\int }_{{r}_{0}}^{{s}_{1}}\,{\rm{d}}{s}_{2}\ \left[{ \mathcal M }({s}_{1}),{ \mathcal M }({s}_{2})\right],\end{eqnarray} \tag{ A5b }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{3}(s)=\ldots ,\end{eqnarray} \tag{ A5c }$$

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

Near r = 0, we can always write

$$\begin{eqnarray}&&{{\bf{u}}}^{{\prime} }(r)=\left(\displaystyle \frac{1}{{r}^{g}}\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal M }}_{\nu }{r}^{\nu }\right){\bf{u}}(r),\end{eqnarray} \tag{ A6 }$$

where ${{ \mathcal M }}_{\nu }$ are constant 4 × 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 }\ne 0$ for only $\nu =0,1,2$.

Generalising the vector ${\bf{u}}$ into the matrix X, and transforming $x=1/r$ changes equation (A6) into

$$\begin{eqnarray}&&{x}^{-q}{X}^{{\prime} }(x)=-{ \mathcal M }(x)X(x),\end{eqnarray} \tag{ A7 }$$

where $q=g-2$ and

$$\begin{eqnarray}&&{ \mathcal M }(x)=\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal M }}_{\nu }{x}^{-\nu }\end{eqnarray} \tag{ A8 }$$

with the only non-zero matrices being

$$\begin{eqnarray}{{ \mathcal M }}_{0}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ {m}_{\uparrow }^{2} & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & {m}_{\downarrow }^{2} & 0\end{array}\right),\end{eqnarray} \tag{ A9a }$$

$$\begin{eqnarray}{{ \mathcal M }}_{1}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ \tfrac{{\rm{\Delta }}}{\lambda } & -1 & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}}{m}_{\downarrow } & 0\\ 0 & 0 & 0 & 0\\ \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}}{m}_{\uparrow } & 0 & \tfrac{{\rm{\Delta }}}{\lambda } & -1,\end{array}\right),\end{eqnarray} \tag{ A9b }$$

$$\begin{eqnarray}{{ \mathcal M }}_{2}=\left(\begin{array}{cccc}0 & 1 & 0 & 0\\ -\tfrac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}-\tfrac{2m}{{{\rm{\hslash }}}^{2}}{\mu }_{\uparrow } & 2\tfrac{{\rm{\Delta }}}{\lambda } & 0 & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}}\\ 0 & 0 & 0 & 1\\ 0 & -\tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & -\tfrac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}-\tfrac{2m}{{{\rm{\hslash }}}^{2}}{\mu }_{\downarrow } & 2\tfrac{{\rm{\Delta }}}{\lambda },\end{array}\right).\end{eqnarray} \tag{ A9c }$$

The matrix ${ \mathcal M }(x)$ is holomorphic at $x=\infty$ meaning that there exists a convergent expansion of the form (A8) for sufficiently large x₀ such that $| x| \gt {x}_{0}$. If the integer $q+1\gt 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 equation (A7) (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)}$

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal A }}_{0} & = & {{{ \mathcal P }}^{(0)}}^{-1}(-{{ \mathcal M }}_{0}){{ \mathcal P }}^{(0)}=\left(\begin{array}{cccc}0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\\ {{ \mathcal P }}^{(0)} & = & \left(\begin{array}{cccc}0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0\\ 0 & -1 & 0 & 0\\ 1 & 0 & 0 & 0\end{array}\right),\end{array}\end{eqnarray} \tag{ A10 }$$

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

$$\begin{eqnarray}{{ \mathcal A }}_{0}=\left(\begin{array}{ccc}{H}_{1} & 0 & 0\\ 0 & {H}_{2} & 0\\ 0 & 0 & {H}_{3}\end{array}\right)={H}_{1}\oplus {H}_{2}\oplus {H}_{3},\end{eqnarray} \tag{ A11 }$$

where H_j ($j=1,2,3$) are shifting matrices where H₁ 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

$$\begin{eqnarray}&&{x}^{-q}{Y}^{{\prime} }(x)={ \mathcal A }(x)Y(x),\end{eqnarray} \tag{ A12 }$$

where q = 0, and

$$\begin{eqnarray}{{ \mathcal A }}_{0}=\left(\begin{array}{cccc}0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ A13a }$$

$$\begin{eqnarray}{{ \mathcal A }}_{1}=\left(\begin{array}{cccc}1 & \tfrac{{\rm{\Delta }}}{\lambda } & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & 1 & -\tfrac{{\rm{\Delta }}}{\lambda }\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ A13b }$$

$$\begin{eqnarray}{{ \mathcal A }}_{2}=\left(\begin{array}{cccc}-\tfrac{2{\rm{\Delta }}}{\lambda } & \tfrac{2m(h-\mu )}{{{\rm{\hslash }}}^{2}}-\tfrac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}} & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & 0\\ 1 & 0 & 0 & 0\\ -\tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & 0 & -\tfrac{2{\rm{\Delta }}}{\lambda } & \tfrac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}+\tfrac{2m(h+\mu )}{{{\rm{\hslash }}}^{2}}\\ 0 & 0 & -1 & 0\end{array}\right).\end{eqnarray} \tag{ A13c }$$

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 equation (A12) into

$$\begin{eqnarray}&&{x}^{-q}{Z}^{{\prime} }(x)={ \mathcal B }(x)Z(x)\end{eqnarray} \tag{ A14 }$$

with $q\geqslant 0$ and

$$\begin{eqnarray}&&{x}^{-q}{{ \mathcal K }}^{{\prime} }(x)={ \mathcal A }(x){ \mathcal K }(x)-{ \mathcal K }(x){ \mathcal B }(x).\end{eqnarray} \tag{ A15 }$$

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

$$\begin{eqnarray}&&{{ \mathcal B }}_{0}={{ \mathcal A }}_{0},\end{eqnarray} \tag{ A16a }$$

$$\begin{eqnarray}&&{{ \mathcal K }}_{0}=I.\end{eqnarray} \tag{ A16b }$$

Using equation (A16), substitution of the series ansätze into equation (A15) and comparison of like coefficients results in the recursion relation

$$\begin{eqnarray}&&{{ \mathcal A }}_{0}{{ \mathcal K }}_{0}-{{ \mathcal K }}_{0}{{ \mathcal A }}_{0}=0,\end{eqnarray} \tag{ A17a }$$

$$\begin{eqnarray}&&{{ \mathcal A }}_{0}{{ \mathcal K }}_{\nu }-{{ \mathcal K }}_{\nu }{{ \mathcal A }}_{0}=\displaystyle \sum _{s=0}^{\nu -1}({{ \mathcal K }}_{s}{{ \mathcal B }}_{\nu -s}-{{ \mathcal A }}_{\nu -s}{{ \mathcal K }}_{s})-(\nu -q-1){{ \mathcal K }}_{\nu -q-1}\qquad (\nu \gt 0),\end{eqnarray} \tag{ A17b }$$

where the last term in equation (A17b) is absent for $\nu -q-1\lt 0$, that is, for ν < 1. The recursion relation is of the form

$$\begin{eqnarray}&&{{ \mathcal A }}_{0}{{ \mathcal K }}_{\nu }-{{ \mathcal K }}_{\nu }{{ \mathcal A }}_{0}={{ \mathcal B }}_{\nu }+{K}_{\nu },\qquad \nu \gt 0\end{eqnarray} \tag{ A18 }$$

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 < ν. 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 }$ (ν > 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 equation (A18) 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 (A18) corresponds to s² = 9 relations

$$\begin{eqnarray}&&{H}_{j}{{ \mathcal K }}_{\nu }^{{jk}}-{{ \mathcal K }}_{\nu }^{{jk}}{H}_{k}={{ \mathcal B }}_{\nu }^{{jk}}+{K}_{\nu }^{{jk}},\qquad \nu \gt 0.\end{eqnarray} \tag{ A19 }$$

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 equation (A19), each equation of equation (A19) can be soluble only if the matrices ${{ \mathcal B }}_{\nu }^{{jk}}$ satisfy some restrictive condition, explained below. Consider the auxiliary equation

$$\begin{eqnarray}&&{HX}-{XK}=M,\end{eqnarray} \tag{ A20 }$$

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

$$\begin{eqnarray}&&\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal K }}_{\nu }{x}^{-\nu }\end{eqnarray} \tag{ A21 }$$

is determined by solving the recursion relations (A19) successively for ν = 1, 2, .... 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₀ such that $| x| \gt {x}_{0}$.

The point of the transformation ${ \mathcal K }$ is to obtain a differential equation (A14) such that the matrices ${\bf{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 ν > 0 occur in the rows corresponding to the last rows of the Jordan blocks H_k (k = 1, 2, ..., 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

$$\begin{eqnarray}{{ \mathcal B }}_{0}=\left(\begin{array}{cccc}0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\qquad {{ \mathcal B }}_{1}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & 1 & -\tfrac{{\rm{\Delta }}}{\lambda }\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ A22a }$$

$$\begin{eqnarray}{{ \mathcal B }}_{2}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & -\tfrac{3{\rm{\Delta }}}{\lambda } & 0 & 0\\ -\tfrac{4m\lambda }{{{\rm{\hslash }}}^{2}} & -\tfrac{2{\rm{\Delta }}m}{{{\rm{\hslash }}}^{2}} & -\tfrac{2{\rm{\Delta }}}{\lambda } & \tfrac{{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}+\tfrac{2m(h+\mu )}{{{\rm{\hslash }}}^{2}}\\ 0 & 0 & -1 & 0\end{array}\right),\end{eqnarray} \tag{ A22b }$$

$$\begin{eqnarray}{{ \mathcal B }}_{3}=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ \tfrac{6{\rm{\Delta }}}{\lambda } & \tfrac{3{{\rm{\Delta }}}^{2}}{{\lambda }^{2}}+\tfrac{4m\left(m{\lambda }^{2}+(\mu -h){{\rm{\hslash }}}^{2}\right)}{{{\rm{\hslash }}}^{4}} & -\tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & -\tfrac{2{\rm{\Delta }}m}{{{\rm{\hslash }}}^{2}}\\ \tfrac{6{\rm{\Delta }}m}{{{\rm{\hslash }}}^{2}} & \tfrac{4m\left({{\rm{\hslash }}}^{2}{{\rm{\Delta }}}^{2}-{hm}{\lambda }^{2}+m{\lambda }^{2}\mu \right)}{\lambda {{\rm{\hslash }}}^{4}} & -\tfrac{4{m}^{2}{\lambda }^{2}}{{{\rm{\hslash }}}^{4}} & 0\\ 0 & 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ A22c }$$

$$\begin{eqnarray}&&{{ \mathcal B }}_{4}=\cdots .\end{eqnarray} \tag{ A22d }$$

We work until $\nu ={\nu }_{\max }$. The final series solutions will then be exact upto $\nu ={\nu }_{\max }-1$ for ${F}_{\uparrow }(r)$ and $\nu ={\nu }_{\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 }_{\max }-2$ and for spin-$\downarrow$ the derivative (component 4) will agree with a direct derivative of component 3 upto $\nu ={\nu }_{\max }-3$.

A.2.2. Simplify the problem by a shearing transformation

Having obtained the matrix ${\bf{B}}$, equation (A14) is ready for the shearing transformation

$$\begin{eqnarray}&&Z=S(x)V\equiv \mathrm{diag}(1,{x}^{-\xi },{x}^{-2\xi },{x}^{-3\xi })V\end{eqnarray} \tag{ A23 }$$

with a temporarily unknown positive parameter ξ, which takes equation (A14) into

$$\begin{eqnarray}&&{x}^{-q}{V}^{{\prime} }(x)={ \mathcal C }V=\left(\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal C }}_{\nu }{x}^{-\nu }\right)V,\end{eqnarray} \tag{ A24 }$$

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

$$\begin{eqnarray}&&{c}_{{jk}}={b}_{{jk}}{x}^{\xi (j-k)}+(j-1)\xi {\delta }_{{jk}}{x}^{-q-1}.\end{eqnarray} \tag{ A25 }$$

Here δ_jk is the Kronecker delta, and b_jk are the matrix elements of ${\bf{B}}$. The parameter ξ 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 ξ must be chosen judiciously as follows. Any ${b}_{{jk}}\ne 0$ is of the form

$$\begin{eqnarray}&&{b}_{{jk}}={x}^{-{\alpha }_{{jk}}}\displaystyle \sum _{\nu =0}^{\infty }{b}_{{jk}\nu }{x}^{-\nu },\end{eqnarray} \tag{ A26 }$$

where ${b}_{{jk}\nu }\ne 0$ and each positive integer ${\alpha }_{{jk}}\geqslant 1$ except for the special elements with a 1 from a shifting matrix in the Jordan decomposition for ${{ \mathcal B }}_{0}$ for which α_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 α_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 ξ. After the shearing transformation the expansions of non-zero off-diagonal elements c_jk (δ_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)}\gt 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 ξ = ξ₀ 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 (A23).

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}$:

$$\begin{eqnarray}\mathop{\mathrm{lim}}\limits_{x\to \infty }{x}^{\xi }{ \mathcal C }(x)=\left(\begin{array}{cccc}0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & \tfrac{2m\lambda }{{{\rm{\hslash }}}^{2}} & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right).\end{eqnarray} \tag{ A27 }$$

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

$$\begin{eqnarray}&&x={{p}^{(1)}}^{1/(\xi -q-1)}{x}_{1}^{{p}^{(1)}}={x}_{1}^{2}/4,\end{eqnarray} \tag{ A28 }$$

to obtain from equation (A24) the differential equation

$$\begin{eqnarray}&&{x}_{1}^{-h}{V}^{{\prime} }({x}_{1})={ \mathcal D }({x}_{1})V=\left(\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal D }}_{\nu }{x}_{1}^{-\nu }\right)V,\end{eqnarray} \tag{ A29 }$$

where

$$\begin{eqnarray}&&h={p}^{(1)}q+{p}^{(1)}-{p}^{(1)}\xi -1=0,\end{eqnarray} \tag{ A30a }$$

$$\begin{eqnarray}\begin{array}{rcl}{d}_{{ij}} & = & {\left({{p}^{(1)}}^{1/(\xi -q-1)}{t}^{{p}^{(1)}}\right)}^{\xi }\left[{b}_{{jk}}{{p}^{(1)}}^{\xi (j-k)/(\xi -q-1)}{x}_{1}^{{q}^{(1)}(j-k)}+(j-1)\displaystyle \frac{{q}^{(1)}}{{p}^{(1)}}{\delta }_{{jk}}{{p}^{(1)}}^{(-q-1)/(\xi -q-1)}{x}_{1}^{{p}^{(1)}(-q-1)}\right]\\ & = & {b}_{{jk}}{\left(\displaystyle \frac{{x}_{1}}{2}\right)}^{j-k+1}+{\delta }_{{jk}}(j-1){x}_{1}^{-1}.\end{array}\end{eqnarray} \tag{ A30b }$$

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\lt 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\geqslant 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 [64] 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 (A12)

$$\begin{eqnarray}&&{Y}^{{\prime} }(x)={ \mathcal A }(x)Y(x)=\left(\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal A }}_{\nu }{x}^{-\nu }\right)Y,\end{eqnarray} \tag{ A31 }$$

to the system

$$\begin{eqnarray}&&{x}_{2}{Y}^{(3)}{\prime} ({x}_{2})={{ \mathcal A }}^{(3)}{Y}^{(3)}({x}_{2})=\left(\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal A }}_{\nu }^{(3)}{x}_{2}^{-\nu }\right){Y}^{(3)},\end{eqnarray} \tag{ A32 }$$

where

$$\begin{eqnarray}{{ \mathcal A }}_{0}^{(3)}=\left(\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & 9 & 1 & 0\\ 0 & 0 & 9 & 0\\ 0 & 0 & 0 & 13\end{array}\right),\qquad {{ \mathcal A }}_{1}^{(3)}=\cdots ,\end{eqnarray} \tag{ A33 }$$

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

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

$$\begin{eqnarray}&&{x}_{2}{Y}^{(6)}{\prime} ({x}_{2})={{ \mathcal A }}^{(6)}{Y}^{(6)}\end{eqnarray} \tag{ A34 }$$

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

Let us now transform back to radial coordinates with ${r}_{2}=1/{x}_{2}$ (here ${x}_{1}={x}_{2}^{2}/4$) to obtain from equation (A34) the equation

$$\begin{eqnarray}&&{r}_{2}{{Y}^{(7)}}^{{\prime} }({r}_{2})={{ \mathcal A }}^{(7)}({r}_{2}){Y}^{(7)}=\left(\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal A }}_{\nu }^{(7)}{r}_{2}^{\nu }\right){Y}^{(7)},\end{eqnarray} \tag{ A35 }$$

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 }}_{0}^{(7)}$ differ by a positive integer, equation (A35) has a fundamental matrix solution of the form

$$\begin{eqnarray}&&{Y}^{(7)}({r}_{2})={{ \mathcal K }}^{(7)}({r}_{2}){\rm{\Upsilon }}({r}_{2}),\qquad {{ \mathcal K }}_{0}^{(7)}\,=\,I,\end{eqnarray} \tag{ A36 }$$

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

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

$$\begin{eqnarray}&&{r}_{2}{{\rm{\Upsilon }}}^{{\prime} }({r}_{2})={{ \mathcal B }}^{(7)}({r}_{2}){\rm{\Upsilon }}=\left(\displaystyle \sum _{\nu =0}^{\infty }{{ \mathcal B }}_{\nu }^{(7)}{r}_{2}^{\nu }\right){\rm{\Upsilon }},\end{eqnarray} \tag{ A37 }$$

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

$$\begin{eqnarray}&&{{ \mathcal A }}_{0}^{(7)}{{ \mathcal K }}_{0}^{(7)}-{{ \mathcal K }}_{0}^{(7)}{{ \mathcal A }}_{0}^{(7)}=0,\end{eqnarray} \tag{ A38a }$$

$$\begin{eqnarray}&&{{ \mathcal A }}_{0}^{(7)}{{ \mathcal K }}_{\nu }^{(7)}-{{ \mathcal K }}_{\nu }^{(7)}{{ \mathcal A }}_{0}^{(7)}=\displaystyle \sum _{s=0}^{\nu -1}({{ \mathcal K }}_{s}^{(7)}{{ \mathcal B }}_{\nu -s}^{(7)}-{{ \mathcal A }}_{\nu -s}^{(7)}{{ \mathcal K }}_{s}^{(7)})+\nu {{ \mathcal K }}_{\nu }^{(7)}\qquad (\nu \gt 0),\end{eqnarray} \tag{ A38b }$$

where

$$\begin{eqnarray}&&{{ \mathcal B }}_{0}^{(7)}={{ \mathcal A }}_{0}^{(7)},\end{eqnarray} \tag{ A39a }$$

$$\begin{eqnarray}&&{{ \mathcal K }}_{0}^{(7)}=I.\end{eqnarray} \tag{ A39b }$$

Given that the matrix ${{ \mathcal A }}_{0}^{(7)}$ in the convergent expansion ${{ \mathcal A }}^{(7)}={\sum }_{\nu =0}^{\infty }{{ \mathcal A }}_{\nu }^{(7)}{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 }}_{0}^{(7)}=I$ such that the formal transformation ${Y}^{(7)}={{ \mathcal K }}^{(7)}({r}_{2}){\rm{\Upsilon }}$ reduces the differential equation (A35) to the form of equation (A37) such that ${{ \mathcal B }}_{\nu }^{(7)}=0$ for all ν > 0. Therefore, the recursion relations (A38) can be solved for ${{ \mathcal K }}^{(7)}({r}_{2})$ such that ${{ \mathcal B }}_{\nu }^{(7)}=0$ for ν > 0. The solution for ${Y}^{(7)}({r}_{2})$ is then given by equation (A36).

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

$$\begin{eqnarray}&&{r}_{2}{{\rm{\Upsilon }}}^{{\prime} }({r}_{2})={{ \mathcal B }}_{0}^{(7)}{\rm{\Upsilon }}\end{eqnarray} \tag{ A40 }$$

is given by the matrix exponential (A3)

$$\begin{eqnarray}&&{\rm{\Upsilon }}({r}_{2})={\rm{\Upsilon }}({r}_{0})\exp \left\{{{ \mathcal B }}_{0}^{(7)}\mathrm{ln}\left(\displaystyle \frac{{r}_{2}}{{r}_{0}}\right)\right\}.\end{eqnarray} \tag{ A41 }$$

Here a set of four linearly independent vectors forms the fundamental matrix solution ${\rm{\Upsilon }}({r}_{2})$, and ϒ(r₀) is the fundamental matrix solution at r₂ = r₀. We can set the auxiliary parameter r₀ = 1 and choose ϒ(r₀) in such a way as to choose the boundary conditions for F_σ(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 ϒ(r₂). A judicious choice for ϒ(1) (with r₀ = 1) then gives the explicit solution shown in equation (15) in the main text.