Multiply quantized vortices in fermionic superfluids: angular momentum, unpaired fermions, and spectral asymmetry

TL;DR

This paper investigates how multiply quantized vortices affect the orbital angular momentum in fermionic superfluids, revealing significant reductions due to unpaired fermions and spectral flow effects, supported by analytical and numerical solutions.

Contribution

It provides a detailed analysis of angular momentum reduction in multiply quantized vortices, highlighting the role of unpaired fermions and spectral asymmetry in the BCS regime.

Findings

Angular momentum is reduced from the BEC regime value in the BCS regime.

Unpaired fermions arise due to spectral flow along vortex states.

Results are supported by analytical and numerical solutions of Bogoliubov-de-Gennes equations.

Abstract

We compute the orbital angular momentum $L_z$ of an s-wave paired superfluid in the presence of an axisymmetric multiply quantized vortex. For vortices with winding number $|k| > 1$, we find that in the weak-pairing BCS regime $L_z$ is significantly reduced from its value $\hbar N k/2$ in the Bose-Einstein condensation (BEC) regime, where $N$ is the total number of fermions. This deviation results from the presence of unpaired fermions in the BCS ground state, which arise as a consequence of spectral flow along the vortex sub-gap states. We support our results analytically and numerically by solving the Bogoliubov-de-Gennes equations within the weak-pairing BCS regime.

Equations

$$\hat{H}=\sum_{\begin{subarray}{c}l\\ n,n^{\prime}\end{subarray}}\left(\begin{array}[]{c}a_{n,l+k\uparrow}^{\dagger}\\ a_{n,-l\downarrow}\end{array}\right)^{T}\left(\begin{array}[]{cc}\epsilon_{n,l+k}\delta_{n,n^{\prime}}&\Delta^{(l)}_{n,n^{\prime}}\\ \Delta^{(l)*}_{n,n^{\prime}}&-\epsilon_{n,-l}\delta_{n,n^{\prime}}\end{array}\right)\left(\begin{array}[]{c}a_{n^{\prime},l+k\uparrow}\\ a_{n^{\prime},-l\downarrow}^{\dagger}\end{array}\right)$$

$$\left(\begin{array}[]{c}b_{m}^{(l)}\\ d_{\bar{m}}^{(l)\dagger}\end{array}\right)=\sum_{n=1}^{M}\left(\begin{array}[]{cc}S_{1,mn}^{(l)}&S_{2,mn}^{(l)}\\ S_{3,\bar{m}n}^{(l)}&S_{4,\bar{m}n}^{(l)}\end{array}\right)\left(\begin{array}[]{c}a_{n,l+k\uparrow}\\ a_{n,-l\downarrow}^{\dagger}\end{array}\right),$$

$$|BCS⟩\sim {\otimes }_l\prod _{m=1}^M{b}_m^{(l)}\prod _{\bar{m}=1}^M{d}_{\bar{m}}^{(l)}|0⟩,$$

$$|BCS⟩\sim {\otimes }_l\text{sideset}{}^{\prime }{\prod }_m{\tilde{b}}_m^{(l)}\text{sideset}{}^{\prime }{\prod }_{\bar{m}}{\tilde{d}}_{\bar{m}}^{(l)}|0⟩.$$

$$\begin{array}{rl}{|BCS⟩}_l=& \left(\prod _{i=1}^{M_{\uparrow }^{(l)}}{\tilde{a}}_{i,l+k\uparrow }^{†}\right)\left(\prod _{i=1}^{M_{\downarrow }^{(l)}}{\tilde{a}}_{i,-l\downarrow }^{†}\right)\\ & \times \exp \left(\sum _{j>M_{\uparrow }^{(l)}}^M\sum _{j^{\prime }>M_{\downarrow }^{(l)}}^M{\tilde{a}}_{j,l+k\uparrow }^{†}{\mathcal{K}}_{j,j^{\prime }}^{(l)}{\tilde{a}}_{j^{\prime },-l\downarrow }^{†}\right)|0⟩,\end{array}$$

$$\mathcal{L}=-\frac{1}{2}\sum _l\left(l+\frac{k}{2}\right){\eta }_l.$$

$$H=\int d^{2}r\,\sum_{\sigma=\uparrow,\downarrow}\psi_{\sigma}^{\dagger}\left(-\frac{\nabla^{2}}{2}+V(r)-\mu\right)\psi_{\sigma}+\big{[}\int d^{2}r\,\psi_{\uparrow}^{\dagger}\Delta(r)e^{ik\varphi/2}(p_{x}+ip_{y})^{\nu}e^{ik\varphi/2}\psi_{\downarrow}^{\dagger}+\text{h.c.}\big{]}$$

$$\Psi=\left(\begin{array}[]{c}\psi_{\uparrow}\\ \psi_{\downarrow}^{\dagger}\end{array}\right),\quad\Psi^{\dagger}=\left(\psi_{\uparrow}^{\dagger}\,,\psi_{\downarrow}\right),$$

$$\{{\Psi }_i(\mathbf{r}),{\Psi }_j^{†}({\mathbf{r}}^{\prime })\}={\delta }_{ij}\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$$\begin{array}{rl}H& =\int d^{2}r{\Psi }^{†}\left(-\frac{{\nabla }^2}{2}+V(r)-\mu \right){\tau }_3\Psi \\ & +\int d^{2}r{\Psi }^{†}\Delta (r)\left(e^{ik\phi /2}(p_x+i{p}_y{)}^{\nu }{e}^{ik\phi /2}{\tau }_{+}+e^{-ik\phi /2}(p_x-i{p}_y{)}^{\nu }{e}^{-ik\phi /2}{\tau }_{-}\right)\Psi .\end{array}$$

$${\hat{L}}_z$$

$$\hat{N}$$

$$(p_x\pm i{p}_y)=-i{e}^{\pm i\phi }\left(\frac{\partial }{\partial r}\pm \frac{i}{r}\frac{\partial }{\partial \phi }\right).$$

$$[{\Psi }^{†}\hat{A}\Psi ,{\Psi }^{†}\hat{B}\Psi ]={\Psi }^{†}[\hat{A},\hat{B}]\Psi .$$

$$[\widehat{L_z},\hat{H}]$$

$$[\hat{N},\hat{H}]$$

$$\hat{\mathcal{L}}={\hat{L}}_z-\frac{k+\nu }{2}\hat{N}$$

$$\sum_{n^{\prime}=1}^{M}H^{(l)}_{n,n^{\prime}}\left(\begin{array}[]{c}u_{n^{\prime}m}^{(l)}\\ v_{n^{\prime}m}^{(l)}\end{array}\right)=E_{m}^{(l)}\left(\begin{array}[]{c}u_{nm}^{(l)}\\ v_{nm}^{(l)}\end{array}\right),$$

$$E_m^{(l)}>0,m=1,\dots ,M_{+}^{(l)},E_{m+M_{+}^{(l)}}^{(l)}<0,m=1,\dots ,M_{-}^{(l)}.$$

$$a_{n,l+k\uparrow }$$

$$a_{n,-l\downarrow }^{†}$$

$$d_{\bar{m}}^{(l)}\equiv b_{\bar{m}+M_{+}^{(l)}}^{(l)†},\bar{m}=1,\dots ,M_{-}^{(l)},$$

$$b_m^{(l)}|BCS⟩$$

$$d_{\bar{m}}^{(l)}|BCS⟩$$

$${|BCS⟩}_l=\prod _{m=1}^M{b}_m^{(l)}\prod _{\bar{m}=1}^M{d}_{\bar{m}}^{(l)}|0⟩,$$

$$\left(\begin{array}[]{c}b_{m}^{(l)}\\ d_{\bar{m}}^{(l)\dagger}\end{array}\right)=\sum_{n=1}^{M}\left(\begin{array}[]{cc}S_{1,mn}^{(l)}&S_{2,mn}^{(l)}\\ S_{3,\bar{m}n}^{(l)}&S_{4,\bar{m}n}^{(l)}\end{array}\right)\left(\begin{array}[]{c}a_{n,l+k\uparrow}\\ a_{n,-l\downarrow}^{\dagger}\end{array}\right),$$

$${\tilde{a}}_{j,l+k\uparrow }$$

$${\tilde{a}}_{j,-l\downarrow }^{†}$$

$${\tilde{b}}_j^{(l)}$$

$${\tilde{d}}_j^{(l)†}$$

$$\left(\begin{array}[]{c}\tilde{b}_{m}^{(l)}\\ \tilde{d}_{\bar{m}}^{(l)\dagger}\end{array}\right)=\sum_{m=1}^{M}\left(\begin{array}[]{cc}\Lambda_{1,mn}^{(l)}&\Sigma_{2,mn}^{(l)}\\ \Sigma_{3,\bar{m}n}^{(l)}&\Lambda_{4,\bar{m}n}^{(l)}\end{array}\right)\left(\begin{array}[]{c}\tilde{a}_{n,l+k\uparrow}\\ \tilde{a}_{n,-l\downarrow}^{\dagger}\end{array}\right),$$

$$S_1^{(l)}{S}_1^{(l)†}+S_2^{(l)}{S}_2^{(l)†}$$

$$S_1^{(l)†}{S}_1^{(l)}+S_3^{(l)†}{S}_3^{(l)}$$

$$\left(S_2^{(l)}{S}_2^{(l)†}\right)y_{2j}^{(l)}$$

Full text

Quantized vortices are a hallmark of superfluids (SFs) and superconductors. These topological defects form in response to external rotation or magnetic field and play a key role in understanding a broad spectrum of phenomena, such as the Berezinskii-Kosterlitz-Thouless transition in two-dimensional (2D) SFs Berezinskii (1970); Kosterlitz and Thouless (1973), superconductor/insulator transitions Fisher (1990); Peskin (1978); Dasgupta and Halperin (1981), turbulence Feynman (1955), and dissipation Anderson (1966); Bardeen and Stephen (1965). In fermionic $s$-wave paired states, the structure of the ground state and low lying excitations of an axisymmetric singly quantized vortex has been established through analytical and numerical studies in both the strong-pairing regime (where the SF phase is understood as a Bose-Einstein condensate (BEC) of bosonic molecules) and in the weak-pairing Bardeen Cooper Schrieffer (BCS) regime. In the BEC regime, the microscopic Gross-Pitaevskii equation provides a reliable framework Gross (1961); Pitaevskii (1961), while in the BCS regime the (self-consistent) Bogoliubov-deGennes (BdG) theory is key in identifying the structure of the ground state Nygaard et al. (2003); Sensarma et al. (2006) and the spectrum of sub-gap fermionic excitations Caroli et al. (1964).

Multiply quantized vortices (MQVs) have however not received much attention. Generically in a homogeneous bulk system, the logarithmic repulsion between vortices, which scales as the square of the vortex winding number $k$, energetically favors an instability of a multiply quantized vortex into separated elementary unit vortices Pethick and Smith (2008). However, MQVs are of interest since under certain circumstances, the interaction between vortices is not purely repulsive and can support multi-vortex bound states, at least as metastable defects. This can happen, for instance, in type-II mesoscopic superconductors, where MQVs have been predicted Schweigert et al. (1998) and experimentally observed Geim et al. (2000); Kanda et al. (2004); Grigorieva et al. (2007); Cren et al. (2011). In addition, it has been argued that MQVs are expected to be energetically stable in multicomponent superconductors Dao et al. (2011); Babaev et al. (2017) and in chiral $p$-wave superconductors Garaud and Babaev (2015); Sauls and Eschrig (2009). In fermionic SFs, a doubly quantized vortex was predicted Volovik and Kopnin (1977) and observed in 3He-A Blaauwgeers et al. (2000). It has further been argued that fast rotating Fermi gases trapped in an anharmonic potential will support an MQV state Lundh (2006); Lundh and Cetoli (2009); Howe et al. (2009). Similar vortex states have been created in rotating BEC experiments Engels et al. (2003); Bretin et al. (2004); Leanhardt et al. (2002); Andersen et al. (2006).

Surprisingly, as we demonstrate in this Letter, there is a fundamental difference between a singly quantized vortex ($|k|=1$) and an MQV ($|k|>1$) in a weakly-paired fermionic $s$-wave SF. This difference is manifested most clearly in the orbital angular momentum (OAM) $L_{z}$, as illustrated in Fig. 1. At zero temperature in the BEC regime, a microscopic Gross-Pitaevskii calculation predicts $L_{z}=\hbar Nk/2$, where $N$ is the total number of fermions. Intuitively, this corresponds to a simple picture where an MQV induces a quantized OAM $k$ per molecule. For an elementary vortex, this result also holds in the BCS regime, as confirmed within the self-consistent BdG framework Nygaard et al. (2003); Sensarma et al. (2006). As we show in this Letter, for vortices with $|k|>1$ however, the BCS ground state contains unpaired fermions which carry OAM opposite to that carried by the Cooper pairs, thereby significantly reducing the total $L_{z}$ from its BEC value by an amount $\sim(k_{F}\xi)^{2}$, where $k_{F}$ is the Fermi momentum and $\xi$ the coherence length. While the proportionality constant is non-universal and depends on the vortex core structure, the scaling with $k_{F}$ and $\xi$ is robust, being independent of any boundary effects.

To derive our main result we consider a 2D111Due to the axial symmetry of the vortex line, it is sufficient to consider a two-dimensional BdG problem with a point vortex. $s$-wave paired SF in the weak-pairing BCS regime at zero temperature within the BdG framework. The mean-field Hamiltonian in the presence of an axisymmetric MQV with winding number $k$ is $\hat{H}=\int d^{2}r\,\Psi^{\dagger}\left[-\nabla^{2}/2+V(r)-\mu\right]\tau_{3}\Psi+\int d^{2}r\,\Psi^{\dagger}\Delta(r)\left(e^{ik\varphi}\tau_{+}+e^{-ik\varphi}\tau_{-}\right)\Psi$, where the Nambu spinor $\Psi=(\psi_{\uparrow},\psi_{\downarrow}^{\dagger})^{T}$ satisfies $\{\Psi_{i}(\mathbf{r}),\Psi_{j}^{\dagger}(\mathbf{r}^{\prime})\}=\delta_{ij}{\delta(\mathbf{r}-\mathbf{r}^{\prime})}$. Here, $\tau_{i}$ are Pauli matrices, $\tau_{\pm}=(\tau_{1}\pm i\tau_{2})/2$, $\hbar$ and the elementary fermion mass are set to unity, and $\mu$ is the chemical potential. In principle, $\Delta(r)$ should be determined self-consistently but since our results depend only weakly on its form, we use a fixed pairing term that for our numerical analysis is taken to be $\Delta(r)$ = $\Delta_{0}\tanh\left(r/\xi\right)$, where $\xi=k_{F}/\Delta_{0}$ and $\Delta_{0}$ is the BCS gap.

Due to the pairing term, neither the total particle number $\hat{N}=\int d^{2}r\,\Psi^{\dagger}\tau_{3}\Psi$ nor the OAM $\hat{L}_{z}=\int d^{2}r\,\Psi^{\dagger}\left(-i\partial_{\varphi}\right)\Psi$ commutes with $\hat{H}$, and so neither are separately conserved. Instead, as pointed out in Salomaa and Volovik (1987); Volovik (1995), the generalized OAM operator $\mathscr{\hat{L}}=\hat{L}_{z}-k\hat{N}/2$ generates a symmetry and thus, the BdG ground state and all quasi-particle excitations carry a sharp $\mathscr{\hat{L}}$ quantum number. More generally, in a chiral SF with pairing symmetry $\sim(p_{x}+ip_{y})^{\nu}$ and with an MQV, the conserved operator is $\hat{L}_{z}-(k+\nu)\hat{N}/2$ (see Supplemental Material sup ). While the OAM of vortex-free chiral paired SFs ($k=0$) was analysed in Tada et al. (2015); Ojanen (2016); Volovik (2015), here we focus on $s$-wave SFs ($\nu=0$) with MQVs, noting that our results readily generalize to chiral states with MQVs.

Physically, $\mathscr{\hat{L}}$ measures the deviation of OAM in the BCS ground state from its expectation value $L_{z}^{\text{BEC}}=kN/2$ in the BEC regime (with $N=<\hat{N}>$). The suppression of $L_{z}$ in the BCS regime will hence be reflected in the eigenvalue $\mathscr{L}$ of $\mathscr{\hat{L}}$, evaluated in the ground state of the BdG Hamiltonian. We consider a disc geometry with Dirichlet boundary conditions, i.e., $V(r<R)=0$ and $V(r>R)=\infty$. Expanding the fermionic operators in a single particle basis as $\psi_{\sigma}(\bm{r})=\sum_{n,l}a_{nl\sigma}\Phi_{nl}(\bm{r})$ where $\Phi_{nl}$ satisfies $\left[-\nabla^{2}/2+V(r)-\mu\right]\Phi_{nl}(\bm{r})=\epsilon_{nl}\Phi_{nl}(\bm{r})$, the Hamiltonian becomes

[TABLE]

with $\Delta^{(l)}_{n,n^{\prime}}=\int d^{2}r\,\Phi^{*}_{n,l+k}\Delta(r)e^{ik\varphi}\Phi^{*}_{n^{\prime},-l}$ and where $n,l$ are the radial and angular momentum quantum numbers respectively. Denoting the single-particle Hamiltonian matrix as $H^{(l)}$, particle-hole (PH) symmetry connects the different $l$-sectors through $H^{(l)*}=-CH^{(-l-k)}C^{-1}$ and the spectrum is hence PH symmetric about $l=-k/2$.

The ground state of the BdG Hamiltonian is constructed using a generalized Bogoliubov transformation Labonté (1974); Ring and Schuck (1980) whose main steps we present here (see Supplemental Material sup for details). First, we regularize the BdG Hamiltonian $H^{(l)}$ by introducing a cutoff $M\gg 1$ on $n,n^{\prime}$. Generically, $H^{(l)}$ will have a different number of positive and negative eigenvalues, $M_{+}^{(l)}$ and $M_{-}^{(l)}$ respectively. The (unitary) Bogoliubov transformation is then written as

[TABLE]

where $m=1,\dots M_{+}^{(l)}$, $\bar{m}=1,\dots M_{-}^{(l)}$, and $M_{+}^{(l)}+M_{-}^{(l)}=2M$. The Bogoliubov operator $b_{m}^{(l)}$ annihilates a quasi-particle with positive energy $E_{m}^{(l)}$, $\mathscr{L}$-charge 222The quasi-particle operators carry a sharp $\mathscr{L}$-charge $l-k/2$, rather than an $l$ quantum number. Nevertheless, since the former differs from $l$ by a constant shift, it is convenient to continue labelling the states by $l$. $l+k/2$, and spin $\uparrow$. Alternatively, by PH symmetry we can interpret it as the creation operator for a spin $\downarrow$ state with negative energy $-E_{m}^{(l)}$ and $\mathscr{L}$-charge $-l-k/2$. In addition, we introduce the operator $d_{\bar{m}}^{(l)}$ that creates a spin $\uparrow$ state with negative energy $E_{M_{+}^{(l)}+\bar{m}}^{(l)}$ and $\mathscr{L}$-charge $l+k/2$.

In terms of these operators, the ground state $\ket{BCS}\sim\otimes_{l}\ket{BCS}_{l}$ is defined as the vacuum for all positive energy quasi-particles and thus satisfies $b_{m}^{(l)}\ket{BCS}=0$ and $d_{\bar{m}}^{(l)}\ket{BCS}=0$. For systems with $M_{+}^{(l)}=M_{-}^{(l)}$, the ground state $\ket{BCS}$ closely resembles a Fermi sea with all negative energy states occupied

[TABLE]

where $\ket{0}$ is the Fock vacuum for $a_{n,l\sigma}$. This ground state can be understood in terms of Cooper pairs, where spin $\uparrow$ quasi-particles with $\mathscr{L}$-charge $v=l+k/2$ (created by $d^{(l)}$) are paired with quasi-particles of the opposite spin $\downarrow$ and with the opposite $\mathscr{L}$-charge $-v$ (created by $b^{(l)}$). Re-expressing the quasi-particle operators in terms of elementary fermions, we find a familiar exponential form, $\ket{BCS}_{l}=\exp\left(a_{n,l+k\uparrow}^{\dagger}K^{(l)}_{n,n^{\prime}}a_{n^{\prime},-l\downarrow}^{\dagger}\right)\ket{0}$, where $K^{(l)}$ is an $M\times M$ matrix (derived in Supplemental Material sup ), and the sum over $n,n^{\prime}$ is implicit. Since $b^{(l)}$ and $d^{(l)}$ carry opposite $\mathscr{L}$-charge, the ground state Eq. (3) has a vanishing $\mathscr{L}$ eigenvalue.

When $M_{+}^{(l)}\neq M_{-}^{(l)}$ however, the ground state is no longer given by Eq. (3) since there will exist an imbalance between the number of quasi-particles with $\mathscr{L}$-charge $l+k/2$ and with $\mathscr{L}$-charge $-l-k/2$. This mismatch is quantified by the spectral asymmetry of the energy spectrum $\eta_{l}=\sum_{m}\text{sgn}(E_{m}^{(l)})=M_{+}^{(l)}-M_{-}^{(l)}$, where $\{E_{m}^{(l)}\}_{m\in\mathbb{N}}$ are the eigenvalues of $H^{(l)}$. In order to demonstrate that the presence of a non-trivial $\eta_{l}$ leads to unpaired fermions in the ground state, we perform a judiciously chosen unitary rotation on $a_{n,l\sigma}$ to a new basis of fermions $\tilde{a}_{j,l\sigma}$ via a conventional (non-Bogoliubov) rotation which does not mix creation and annihilation operators (see Supplemental Material sup ). Through a separate unitary rotation, we simultaneously transform the Bogoliubov operators $b^{(l)},d^{(l)}$ into a new basis $\tilde{b}^{(l)},\tilde{d}^{(l)}$. The new fermions $\tilde{a}$ and Bogoliubov quasi-particles $\tilde{b},\tilde{d}$ are related through a Bogoliubov transformation which, as always, takes the schematic form $\tilde{b}=U\tilde{a}+V\tilde{a}^{\dagger}$, where the matrix-valued coefficients $U,V$ satisfy $|U|^{2}+|V|^{2}=1$. Following Ring and Schuck (1980), we find that the preceding transformations naturally distinguish between operators for which either $U$ vanishes exactly: $U=0,V=1$ (occupied levels), or $V$ vanishes exactly: $V=0,U\neq 0$ (empty levels), with the remaining operators, for which both $U,V\neq 0$, describing paired levels. In the new basis, the ground state is superficially similar to Eq. (3) since it can be expressed as

[TABLE]

Importantly however, the restricted products here run only over paired and occupied levels. Bogoliubov operators $\tilde{b},\tilde{d}$ for empty states, which are linear super-positions of $\tilde{a}$’s, annihilate the bare vacuum $\ket{0}$ and are thus disallowed in Eq. (4). Conversely, occupied states contribute to Eq. (4) but since these states create unitarily rotated fermions with certainty, $\tilde{b},\tilde{d}\sim\tilde{a}^{\dagger}$, they do not participate in pairing. The expression (4) is in turn equivalent to (see Supplemental Material supmatfordetails

[TABLE]

where $M^{(l)}_{\downarrow}$ and $M^{(l)}_{\uparrow}$ are the number of occupied (and also empty) $\tilde{b}_{m}^{(l)}$ and $\tilde{d}_{\bar{m}}^{(l)}$ levels respectively. In terms of these parameters, the spectral asymmetry $\eta_{l}=2(M^{(l)}_{\downarrow}-M^{(l)}_{\uparrow})$, with $M^{(l)}_{\uparrow,\downarrow}=\max(0,M-M^{(l)}_{+,-})$.

The exponential part of $\ket{BCS}$ explicitly illustrates the singlet pairing while $M_{\sigma}^{(l)}\neq 0$ signals the presence of unpaired fermions in the ground state. The eigenvalue of $\mathscr{\hat{L}}$ can now be obtained directly from Eq. (4) by summing the individual contributions of the filled quasi-particle states and noting that $\tilde{b}^{(l)},\tilde{d}^{(l)}$ carry the same $\mathscr{L}$-charges as $b^{(l)},d^{(l)}$. While contributions from the paired levels cancel out, the occupied levels lead to

[TABLE]

Alternatively, this equation can be derived directly from Eq. (5), and has previously appeared in the literature in the context of chiral SFs Tada et al. (2015); Ojanen (2016), where $k$ is replaced by the chirality $\nu$. Physically, Eq. (6) quantifies the contribution of unpaired fermions to the OAM.

The physics originating from unpaired fermions in the ground state of a paired state was previously identified and studied in nuclear physics Ring and Schuck (1980), FFLO superfluids Sheehy and Radzihovsky (2007), and chiral superfluids paired in higher partial waves Tada et al. (2015); Volovik (2015); Huang et al. (2014, 2015); Ojanen (2016). We now demonstrate that for a weakly-paired $s$-wave SF with an MQV, a nontrivial $\eta_{l}$ and the associated unpaired fermions arise as a consequence of vortex core states.

In the BCS regime, the spectrum of the vortex core (vc) states for a singly quantized vortex $|k|=1$ was calculated analytically by Caroli-deGennes-Matricon (CdGM) Caroli et al. (1964) who found a single branch $E_{\text{vc}}^{(l)}$ (per spin projection) that crosses the Fermi level. This branch is PH symmetric with respect to itself, $E_{\text{vc}}^{(l)}=-E_{\text{vc}}^{(-l-1)}$ and at low energies ($E_{\text{vc}}\ll\Delta_{0}$) behaves linearly $E^{(l)}_{\text{vc}}=-\omega_{0}(l+1/2)$, where the mini-gap $\omega_{0}\sim\Delta_{0}/(k_{F}\xi)$. By numerically diagonalizing $H^{(l)}$ for $k=1$, we find that $\eta_{l}=0$ for all $l$ and hence there are no unpaired fermions in the BCS ground state of an $s$-wave paired SF with an elementary vortex. Eq. (6) then predicts $\mathscr{L}=0$ and thus the ground state expectation value $L_{z}=N/2$, which agrees with self-consistent BdG calculations Nygaard et al. (2003). The physics here is analogous to that of weakly paired $p+ip$ SFs, where there is a single PH symmetric edge mode that carries no OAM Stone and Roy (2004); Stone and Anduaga (2008); Tada et al. (2015).

For an MQV with winding number $k$, the CdGM method can be generalized and the vortex core spectrum analytically calculated within the BdG framework (see Supplemental Material sup ). In agreement with an argument relating the number of vortex core branches to a topological invariant Volovik (1993), we find that $|k|$ branches (per spin projection) cross the Fermi level. At low energies these branches disperse linearly, $E_{j}(l)=-\omega_{0}(l-l_{j})$, where $j=1,\dots,k$ indexes the branches and the $l_{j}$’s are the angular momenta at which the branches cross the Fermi level. This is consistent with results obtained by numerically diagonalizing the BdG Hamiltonian $H^{(l)}$ (for $k=2$, see Fig. 2) and with previous results on MQVs in superconductors, obtained through quasi-classical approximations Volovik (1993); Mel’nikov and Silaev (2006); Mel’nikov et al. (2008) and numerical simulations Tanaka et al. (1993, 2002); Rainer et al. (1996); Virtanen and Salomaa (1999).

Since in the BEC regime the spectrum is completely gapped for any $k$, we find $\eta_{l}=0$ for all $l$ and thus the ground state OAM is exactly $L_{z}=kN/2$. On the other hand, in the weakly-paired regime the energy spectrum of an MQV exhibits a nontrivial spectral asymmetry. We consider the case $k=2$ first (Fig. 2), where there exist two vortex core branches with linear dispersions at low energies, $E_{\text{vc},\pm}^{(l)}\sim-\omega_{0}(l-l_{\pm})$ with $l_{+}>l_{-}$. Under PH symmetry, these branches are exchanged as $E_{\text{vc},+}^{(l)}=-E_{\text{vc},-}^{(-l-2)}$ which fixes $l_{-}=-(l_{+}+2)$. As shown in Fig. 2, we find that at these crossing points $\eta_{l}$ acquires a non-zero value: $\eta_{l}=-2$ for $l_{-}<l<-1$ and $\eta=+2$ for $-1<l<l_{+}$, with $\eta_{l}=0$ at $l=-1$. Intuitively, this can be understood as follows—at large negative $l$, the branches are merged into the bulk and since there are no sub-gap states, $\eta_{l}=0$. On increasing $l$, the branches begin separating from the bulk but since both have positive energy, $\eta_{l}$ still vanishes. At $l_{-}$ however, one of the branches crosses the Fermi energy, creating a difference of precisely two between the number of negative and positive energy eigenvalues of $H^{(l)}$. At $l=-1$, $\eta_{l}$ necessarily vanishes due to PH symmetry, which also fixes $\eta_{l}$ for $l>-1$. In contrast with $|k|=1$, the branches are not PH symmetric with respect to themselves, allowing the spectral asymmetry to acquire a non-zero value in the BCS regime. The fact that $\eta_{l}$ changes from the BEC to the BCS regime can also be understood as a consequence of spectral flow along the vortex core states, since $\eta_{l}$ (and hence $\mathscr{L}$) cannot change its value in any other way.

A non-zero spectral asymmetry $\eta_{l}$ appears generally for any $|k|\geq 2$ within the BCS regime: for even $k$ (see Fig. 2), there are $|k|/2$ pairs of branches such that the branches within each pair are PH symmetric with each other. $\eta_{l}$ then changes by $\pm 2$ whenever one of these branches crosses the Fermi level; for odd $k$ (see Fig. 2) there are $(|k|-1)/2$ pairs that contribute to a non-trivial $\eta_{l}$, since the branches within each pair go into each other under a PH transformation, while the remaining branch is PH symmetric with respect to itself and therefore does not contribute to $\eta_{l}$.

Having established the existence of a non-vanishing $\eta_{l}$, we see that there must exist unpaired fermions in the BCS ground state for $|k|\geq 2$, and as a consequence of Eq. (6), $\mathscr{L}$ acquires a non-trivial ground state eigenvalue. For $k=2$, this is $\mathscr{L}=-l_{+}^{2}-l_{+}$, where we used PH symmetry to relate $l_{-}$ to $l_{+}$. Importantly, the analytic calculation of the vortex core states (performed in Supplemental Material sup ) demonstrates that the positions of the crossing points are located at $l_{\pm}\sim k_{F}\xi$ with the pre-factor fixed by the form of $\Delta(r)$. This scaling persists in self-consistent numerical calculations Tanaka et al. (2002); Rainer et al. (1996); Virtanen and Salomaa (1999). Eq. (6) along with this scaling thus establishes the reduction of the OAM of the $k=2$ MQV in the weakly paired regime. To leading order in $k_{F}\xi$, $\mathscr{L}=L_{z}-N\sim-\left(k_{F}\xi\right)^{2}$. As a result, the OAM is significantly suppressed from $L_{z}^{BEC}=N$ since $k_{F}\xi\gg 1$ in the BCS regime ($\Delta_{0}\ll E_{F}$). This analysis confirms that the unpaired fermions carry angular momentum opposite to that carried by the Cooper pairs. On a disc, $N\approx(k_{F}R)^{2}/2$, leading to $L_{z}/N\approx 1-\alpha\left(\xi/R\right)^{2}$, where $\alpha$ is an $O(1)$ constant fixed by $\Delta(r)$. As an independent check, we have verified this behavior by numerically calculating $L_{z}/N$ using the full BdG solution (see Supplemental Material sup for details). In Fig. 3, the quadratic scaling is shown to be in good agreement with the numerical data. We thus expect a substantial reduction of the OAM in the BCS regime, where $\xi$ can be comparable to $R$ Geim et al. (2000). We also expect that when two elementary vortices merge into a $k=2$ MQV Mel’nikov and Silaev (2006); Mel’nikov et al. (2008), the ground state OAM decreases from $L_{z}=N$ by an amount $\sim(k_{F}\xi)^{2}$.

A central feature of our result is that the suppression of $L_{z}$ for $|k|\geq 2$ is independent of any boundary effects and is solely determined by the splitting between the vortex core branches. Given this insensitivity to boundary details, we expect our results to hold for more general sample geometries, which may lack axial symmetry. Unlike the ground state energy, which might depend strongly on the gap profile, the OAM thus exhibits universal scaling behavior in the weak pairing BCS regime. The lack of dependence of the OAM on the system boundary is in stark contrast with weakly-paired chiral (e.g., $d+id$) SFs, where it was shown Tada et al. (2015); Volovik (2015) that the OAM is suppressed due to the topological edge modes, but that this effect is strongly dependent on the edge details Huang et al. (2014, 2015); Ojanen (2016); Tada (2015). Our analysis hence suggests that $s$-wave SFs with MQVs may prove to be a more robust platform for investigating the intriguing suppression of the OAM in paired SFs. While the OAM has been measured in SFs Chevy et al. (2000); Hodby et al. (2003); Riedl et al. (2011), we also expect signatures of unpaired fermions—which create a current localized around the vortex core that flows counter to the superflow—in local supercurrent density measurements in MQV states Kanda et al. (2004).

Acknowledgements:

We acknowledge useful discussions with Egor Babaev, Masaki Oshikawa, Michael Stone, Yasuhiro Tada, and Grigory Volovik. A.P. thanks W. Cairncross for helpful comments on the draft. A.P. and V.G. acknowledge support by NSF grants DMR-1205303 and PHY-1211914. The work of S.M. is supported by the Emmy Noether Programme of German Research Foundation (DFG) under grant No. MO 3013/1-1. This research was supported in part by the National Science Foundation under Grant No. DMR-1001240 (L.R.), through the KITP under Grant No. NSF PHY-1125915 (S.M. and L.R.) and by the Simons Investigator award from the Simons Foundation (L.R.). We thank the KITP for its hospitality during our stay as part of the “Universality in Few-Body Systems” (SM), ”Synthetic Quantum Matter” (L.R.), and sabbatical (L.R.) programs, when part of this work was completed.

I Supplemental Material

I.1 A. Proof of $[\mathscr{\hat{L}},\hat{H}]=0$

Consider the BdG Hamiltonian for a two-dimensional superfluid (SF) with an MQV carrying the winding number $k$

[TABLE]

Here the case $\nu=0$ corresponds to the $s$-wave spin-singlet paring, while spin-singlet chiral SFs are obtained by taking non-zero even values of $\nu$. The fermionic operators satisfy canonical anti-commutation relations: $\{\psi_{\sigma}(\mathbf{r}),\psi_{\sigma^{\prime}}^{\dagger}(\mathbf{r}^{\prime})\}=\delta_{\sigma\sigma^{\prime}}\delta(\mathbf{r}-\mathbf{r^{\prime}})$. In terms of Nambu spinors

[TABLE]

that satisfy

[TABLE]

the Hamiltonian becomes

[TABLE]

The angular momentum and particle number operators are

[TABLE]

respectively, where $\tau_{i}$ are the Pauli matrices and $\tau_{\pm}=\frac{1}{2}(\tau_{1}\pm i\tau_{2})$. We work in polar coordinates, where

[TABLE]

In addition, the anti-commutation relations Eq. (S3) lead to the following identity

[TABLE]

Putting everything together, it is straightforward now to show that

[TABLE]

From this it follows that

[TABLE]

commutes with $\hat{H}$. The above calculation easily generalizes to the case of spin-triplet chiral SFs (where $\nu$ is odd) with MQVs, where the same operator Eq. (S9) is conserved.

I.2 B. Construction of ground state wave function

A generalized framework for deriving the ground state of a paired Hamiltonian through a Bogoliubov transformation was constructed in Labonté (1974). Here, we present a self-contained discussion and obtain the ground state wave functions for SFs with MQVs. While we focus on $s$-wave SFs here, this construction can be readily generalized to chiral SFs with vortices as well.

The eigenstates $(u,v)^{T}$ of the BdG Hamiltonian satisfy

[TABLE]

where we have introduced a cut-off $M\gg 1$ on the radial quantum numbers $n,n^{\prime}$. Suppose the number of positive and negative eigenvalues of $H^{(l)}$ are $M_{+}^{(l)}$ and $M_{-}^{(l)}$ respectively, with $M_{+}^{(l)}+M_{-}^{(l)}=2M$. In the absence of spectral asymmetry, $M_{+}^{(l)}=M_{-}^{(l)}$, but in general, $M_{+}^{(l)}\neq M_{-}^{(l)}$. Let us now order the energies such that $E_{1}^{(l)}\geq\dots\geq E_{2M}^{(l)}$, with

[TABLE]

Next, we introduce the (inverse) Bogoliubov transformation

[TABLE]

Here, $b_{m}^{(l)}$ are Bogoliubov operators that annihilate a spin $\uparrow$ state with energy $E^{(l)}_{m}$ and $\mathscr{L}$-charge $l+k/2$. We can exploit the PH symmetry of the system to alternatively interpret $b_{m}^{(l)}$ as the creation operator for a spin $\downarrow$ state with energy $-E_{m}^{(l)}$ and $\mathscr{L}$-charge $-l-k/2$. In addition, we have introduced another set of Bogoliubov operators

[TABLE]

such that the operator $d_{\bar{m}}^{(l)}$ creates a spin $\uparrow$ state with energy $E^{(l)}_{\bar{m}+M_{+}^{(l)}}<0$ and $\mathscr{L}$-charge $l+k/2$. As a matter of principle, we note that since the pairing Hamiltonian does not commute with the angular momentum operator $\hat{L}_{z}$, the Bogoliubov quasi-particles $b$ and $d$ carry $\mathscr{L}$-charge rather than an $l$ quantum number and thus the energy eigenvalues $E_{k}^{(l)}$ should be labelled instead by their $\mathscr{L}$ quantum number, $v=l-k/2$. However, since $v$ is simply $l$ shifted by a constant, it is more convenient to continue labelling the eigenvalues and quasi-particles by $l$.

The BCS ground state is defined as the vacuum with respect to all positive energy quasi-particles, $\ket{BCS}\sim\otimes_{l}\ket{BCS}_{l}$, and must thus satisfy

[TABLE]

When $M_{+}^{(l)}=M_{-}^{(l)}$, the ground state may be expressed as the state with all negative energy excitations occupied, with

[TABLE]

where $\ket{0}$ is the Fock vacuum with respect to elementary fermions $a_{n,l\sigma}$. The paired nature of this state is evident in the wave function Eq. (S15) as expressed in terms of the Bogoliubov operators, since a spin $\downarrow$ quasi-particle with $\mathscr{L}$-charge $-l-k/2$ (created by $b_{m}^{(l)}$) is paired with a spin $\uparrow$ quasi-particle with $\mathscr{L}$-charge $l+k/2$ (created by $d_{\bar{m}}^{(l)}$).

However, when $M_{+}^{(l)}\neq M_{-}^{(l)}$, there will be some states left unpaired as a consequence of the asymmetry in the spectrum. In order to elucidate the nature of the ground state in the presence of unpaired fermions, it is instructive to transform to a particular bases of elementary fermions and of Bogoliubov quasi-particles in which the structure of the ground state becomes especially transparent. Before presenting technical details of this procedure, we briefly describe the steps involved.

Fig. S1 illustrates the series of transformations that we perform in order to express $\ket{BCS}$ in a transparent form. First, we unitarily rotate the elementary fermions operators $a_{n,l\sigma}$ into a new basis of fermionic operators $\tilde{a}_{n,l\sigma}$ (where the $\tilde{a}$’s are linear combinations of only $a$’s but not $a^{\dagger}$’s). While the states destroyed by operators $\tilde{a}_{n,l\sigma}$ carry well defined spin and angular momentum quantum numbers, $\sigma$ and $l$ respectively, they are not energy eigenstates of the non-interacting Hamiltonian. In a similar spirit, we will unitarily rotate the quasi-particle operators $b,d$ into a new basis of quasi-particle operators $\tilde{b},\tilde{d}$. Importantly, since this operation does not mix quasi-particle creation and annihilation operators, the ground state is still defined as the vacuum with respect to positive energy excitations and hence satisfies $\tilde{b_{m}^{(l)}}\ket{BCS}=0$ and $\tilde{d}_{\bar{m}}^{(l)}\ket{BCS}=0$.

Having established these new bases, we will then relate the rotated fermions $\tilde{a}$ to the rotated quasi-particles $\tilde{b},\tilde{d}$ through a Bogoliubov transformation akin to Eq. (I.2). This relation then allows us to express $\ket{BCS}$ in terms of the rotated fermions $\tilde{a}$ in a manner that makes explicit the nature of pairing in the ground state since it naturally distinguishes between paired and unpaired states. We note that the purpose of these transformations is not to diagonalize the BdG Hamiltonian $H^{(l)}$ (which is diagonalized by the original Bogoliubov transformation Eq. (I.2)) but rather to explicitly construct the ground state wave function for an $s$-wave SF with an MQV.

We now discuss the above procedure in detail. To begin, we first invert Eq. (I.2) and write the unitary Bogoliubov transformation as

[TABLE]

where $S_{1}^{(l)},S_{2}^{(l)}$ are $M_{+}^{(l)}\times M$ and $S_{3}^{(l)},S_{4}^{(l)}$ are $M_{-}^{(l)}\times M$ dimensional matrices, respectively. Next, we perform a singular value decomposition on the matrices $S_{2}$ and $S_{3}$, $S_{i}^{(l)}=U_{i}^{(l)}\Sigma_{i}^{(l)}V_{i}^{(l)\dagger}$ where $U_{i}^{(l)},V_{i}^{(l)}$ are unitary matrices and $\Sigma_{i}^{(l)}$ is a rectangular diagonal matrix with non-negative real entries. We then perform unitary rotations on the elementary fermions $a$ and the Bogoliubov quasi-particles, $b$ and $d$,

[TABLE]

It is straightforward to check that these transformed operators satisfy the canonical anti-commutation relations. In this new basis, the Bogoliubov transformation is expressed as

[TABLE]

where $\Lambda_{1}^{(l)}=U_{2}^{(l)\dagger}S_{1}^{(l)}V_{3}^{(l)}$ and $\Lambda_{4}^{(l)}=U_{3}^{(l)\dagger}S_{4}^{(l)}V_{2}^{(l)}$. Before proceeding with the construction of the ground state, it is necessary to establish some properties of $\Sigma_{2}^{(l)}$ and $\Sigma_{3}^{(l)}$. In particular, we will show now that in the presence of a non-trivial spectral asymmetry $\eta_{l}=M_{+}^{(l)}-M_{-}^{(l)}$, either one of $\Sigma_{2}^{(l)}$ or $\Sigma_{3}^{(l)}$ have $|\eta_{l}|/2$ entries on the diagonal which are equal to one. To prove this, let us assume without loss of generality that $M_{+}^{(l)}<M<M_{-}^{(l)}$. Since the Bogoliubov transformation is unitary, the transformation matrix S^{(l)}\equiv\left(\begin{array}[]{cc}S^{(l)}_{1}&S^{(l)}_{2}\\ S^{(l)}_{3}&S^{(l)}_{4}\end{array}\right) satisfies $S^{(l)\dagger}S^{(l)}=S^{(l)}S^{(l)\dagger}=\mathbb{I}_{2M\times 2M}$. This leads to the conditions

[TABLE]

Since $S_{3}^{(l)\dagger}S_{3}^{(l)}$ and $S_{2}^{(l)}S_{2}^{(l)\dagger}$ are Hermitian matrices, they have real eigenvalues and eigenvectors,

[TABLE]

Eq. (I.2) then implies

[TABLE]

The singular values of $S_{1}^{(l)}$, however, are the square roots of the non-zero eigenvalues of both $S_{1}^{(l)\dagger}S_{1}^{(l)}$ and $S_{1}^{(l)}S_{1}^{(l)\dagger}$. Since $S_{1}^{(l)\dagger}S_{1}^{(l)}$ has exactly $M-M_{+}^{(l)}$ more eigenvalues than $S_{1}^{(l)}S_{1}^{(l)\dagger}$, those extra eigenvalues must necessarily be zero. This in turn implies that $S_{3}^{(l)\dagger}S_{3}^{(l)}$ has $M_{\uparrow}^{(l)}=M-M_{+}^{(l)}$ unity eigenvalues, or equivalently, that $\Sigma^{(l)}_{3}$ has precisely $M_{\uparrow}^{(l)}$ unity entries on the diagonal. In addition, since $\Sigma^{(l)}_{3}$ is a rectangular diagonal $M_{-}^{(l)}\times M$ matrix, its last $M^{(l)}_{-}-M=M-M_{+}^{(l)}=M_{\uparrow}^{(l)}$ rows contain only zeros.

Following this discussion, in general for any $M_{+}^{(l)}$ and $M_{-}^{(l)}$ we can define the quantities $M_{\uparrow}^{(l)}=\text{max}\left(M-M_{+}^{(l)},0\right)$ and $M_{\downarrow}^{(l)}=\text{max}\left(M-M_{-}^{(l)},0\right)$. Using the fact that $M_{+}^{(l)}+M_{-}^{(l)}=2M$, we can also show that $M_{+}^{(l)}=M+M_{\downarrow}^{(l)}-M_{\uparrow}^{(l)}$ and $M_{-}^{(l)}=M-M_{\downarrow}^{(l)}+M_{\uparrow}^{(l)}$. From this, it follows that

[TABLE]

Combining the above results with the unitarity of the Bogoliubov transformation Eq. (S18), we conclude that in the presence of spectral asymmetry we will have cases where $\tilde{b}^{(l)}=\tilde{a}_{-l\downarrow}^{\dagger}$ or $\tilde{d}^{(l)}=\tilde{a}_{l+k\uparrow}^{\dagger}$, which will give rise to unpaired fermions in the ground state.

Specifically, we find that the Bogoliubov transformation Eq. (S18) splits into three classes

[TABLE]

[TABLE]

and

[TABLE]

where we have omitted the $(l)$ superscript to simplify notation. Physically, these three classes correspond to three different kinds of quasi-particle operators Ring and Schuck (1980):

• •

Occupied levels Eq. (I.2) are those where $\tilde{b}=\tilde{a}_{\downarrow}^{\dagger}$ or $\tilde{d}=\tilde{a}_{\uparrow}^{\dagger}$. These operators create a unitarily rotated fermion with unit probability.

• •

Paired levels Eq. (I.2) are those for which $\tilde{b}$ and $\tilde{d}$ are non-trivial superpositions of $\tilde{a}_{\uparrow}$ and $\tilde{a}_{\downarrow}^{\dagger}$.

• •

Empty levels Eq. (I.2) $\tilde{b}$ and $\tilde{d}$ are linear superpositions of $\tilde{a}_{\uparrow}$’s and $\tilde{a}_{\downarrow}$’s, respectively. These operators annihilate the Fock vacuum $\ket{0}$.

The ground state is still the vacuum with respect to all positive energy quasi-particles and, in terms of the unitarily rotated Bogoliubov quasi-particles, it is given by $\ket{BCS}\sim\otimes_{l}\ket{BCS}_{l}$ with

[TABLE]

where the restricted product runs over all paired and occupied levels. Since empty levels annihilate the Fock vacuum, they are not included in Eq. (S26). By construction this state satisfies

[TABLE]

for all $m=1,\dots,M_{+}^{(l)}$ and $\bar{m}=1,\dots,M_{-}^{(l)}$. We can further simplify the ground state since it factorizes into unpaired and paired terms,

[TABLE]

The first two terms in the product clearly indicate that for $M_{\sigma}^{(l)}\neq 0$, there are unpaired fermions in the ground state. Following Labonté (1974), we can now write the paired part of the ground state in terms of creation operators $\tilde{a}^{\dagger}$ of elementary (unitarily rotated) fermions

[TABLE]

where the kernel $\mathcal{K}$ satisfies

[TABLE]

with $j,j^{\prime\prime}>M_{\downarrow}^{(l)}$.

In order to derive Eq. (6) in the main text, we note that we can express the conserved operator $\mathscr{\hat{L}}=\hat{L}-k\hat{N}/2$ in terms of the elementary fermions as

[TABLE]

where the last equality follows since the $\tilde{a}$’s are unitarily related to the $a$’s. Since the set of operators that appear in the exponential part of $\ket{BCS}_{l}$ anti-commute with those that appear in the products (unpaired fermions), we can consider the action of $\mathscr{\hat{L}}$ on these separately. First, we evaluate the commutator

[TABLE]

and similarly,

[TABLE]

which follows from the usual anti-commutation relations satisfied by the rotated fermions $\tilde{a}$. Next, we consider the action of $\mathscr{\hat{L}}$ on the exponential (paired) sector of the wave function. This is done by first calculating the contribution from the spin $\uparrow$ sector,

[TABLE]

and then observing that the contribution from the spin $\downarrow$ sector

[TABLE]

exactly compensates for that coming from the spin $\uparrow$ sector. Hence, we see that $\mathscr{\hat{L}}$ commutes with the exponential part of the wave function and so the paired levels do not contribute to $\mathscr{L}$,

[TABLE]

Since the $\tilde{a}$’s are unitarily related to the $a$’s, $\ket{0}$ is also the Fock vacuum with respect to $\tilde{a}_{j,l\sigma}$. Putting the above together, we hence find that the eigenvalue $\mathscr{L}$ of the operator $\mathscr{\hat{L}}$ when evaluated in the ground state $\ket{BCS}\sim\otimes_{l}\ket{BCS}_{l}$, with $\ket{BCS}_{l}$ given by Eq. (S29), is

[TABLE]

where $\eta_{l}=2(M_{\downarrow}^{(l)}-M_{\uparrow}^{(l)})$.

I.3 C. Analytic solution for vortex core states

We generalize the CdGM method Caroli et al. (1964) for obtaining the vortex core bound states in an $s$-wave superconductor to the case of multiply quantized vortices. For a step-like pair-potential, an explicit solution was obtained previously in Bardeen et al. (1969); Berthod (2005). The procedure presented here applies more generally to any pairing term $\Delta(r)$ and agrees with that calculation where the regimes of validity overlap.

We start with an $s$-wave state with a vortex of vorticity $k$. The (symmetric) gap function is thus $\Delta(r)e^{ik\varphi}$ where $\lim_{r\to 0}\Delta(r)=0$ and $\lim_{r\to\infty}\Delta(r)=\Delta_{0}$. The BdG equations are hence

[TABLE]

Separating the angular and radial dependence of the BdG solutions, we let

[TABLE]

where $l\in\mathbb{Z}\left(\mathbb{Z}+\frac{1}{2}\right)$ if $k\in\text{Even(Odd)}$, such that Eq. (S38) becomes

[TABLE]

We now rewrite these equations in the form

[TABLE]

where we have defined $\alpha=\sqrt{l^{2}+\frac{k^{2}}{4}}$ and $\beta=lk$.

In order to derive an analytic solution for these coupled equations, we introduce a radius $r=r^{*}$ such that $\frac{1}{k_{F}}\ll r^{*}\ll\xi$, where $k_{F}=\sqrt{2\mu}$ is the Fermi momentum and $\xi=k_{F}/\Delta_{0}$ is the coherence length. We then consider the BdG equations (S38) separately in the limits where $r\ll r^{*}$ and where $r\gg r^{*}$. Demanding that the wave function be continuous, we then match the solutions from these two regimes at $r=r^{*}$ to arrive at a solution that holds over the entire range.

We first consider the limit where $r\ll r^{*}$. For physically relevant pairing terms, $\Delta(r)\to 0$ in this limit and we are hence justified in ignoring the pairing term. The BdG equations (S38) thus decouple in this limit,

[TABLE]

The solutions to these equations are Bessel functions parametrized by $\sqrt{2(\mu\pm E)}$. Since we are interested in understanding the nature of the vortex core states close to zero energy, we make an additional approximation and consider energies such that $E\ll\mu$. Reinstating all the proper units, this corresponds to

[TABLE]

where we have defined $p=\frac{Em_{e}}{\hbar^{2}k_{F}}$. Thus, the solutions for $r\ll r^{*}$ are

[TABLE]

where $C_{1}$ and $C_{2}$ are arbitrary constants and where the Bessel functions of the second kind $Y$ are chosen to have vanishing amplitudes since we require a well behaved solution in the limit $r\to 0$.

Next, we consider the case where $r\gg r^{*}$. In this limit, we expect that the pairing term is approximately constant $\Delta(r)\to\Delta_{0}$. Hence, we write the BdG solutions as rapidly oscillating Hankel functions enveloped by functions that vary slowly i.e.,

[TABLE]

where $H^{(1,2)}(r)$ are Hankel functions of the first and second kind. Substituting this ansatz into the BdG equations (S38) and considering only the $H^{(1)}$ component of the solution (since the other follows from this immediately), we find

[TABLE]

Here, $H_{\alpha}\equiv H^{(1)}_{\alpha}$. Since we are in the regime $k_{F}r\gg 1$, we use the asymptotic expansion for the Hankel functions

[TABLE]

and drop the terms $f^{\prime\prime},g^{\prime\prime},\frac{f^{\prime}}{r},\frac{g^{\prime}}{r}$, which is justified since we are considering the limit $r\gg\frac{1}{k_{F}}$. With these further approximations, we find that the slowly varying envelope functions $f,g$ satisfy

[TABLE]

To solve these coupled equations, we treat the right hand side as a perturbation. To zeroth order,

[TABLE]

Imposing the condition that the solutions remain well behaved as $r\to\infty$, we find that the solutions are

[TABLE]

In order to find the solution to first order in $E$, we make the ansatz

[TABLE]

which when substituted into Eq. (I.3) leads to

[TABLE]

Approximating $\sin(\psi)\approx\psi$, this is equivalent to

[TABLE]

the solution to which is

[TABLE]

With this, we find that the leading order solution of Eq. (S38) in the limit $r\gg r^{*}$ is

[TABLE]

with $\psi(r)$ given by Eq. (S54).

We have thus constructed the general solution of the BdG equations (S38), with the solution in the limits $r\ll r^{*}$ and $r\gg r^{*}$ given by Eq. (I.3) and Eq. (S55) respectively. In order to completely fix the undetermined coefficients and to get the quantization condition on the energy, we demand that the solution be continuous and match the solution in the regime $r<r^{*}$ with that in the regime $r>r^{*}$.

Since we wish to extend the solution for $r\ll r^{*}$ Eq. (I.3) to the vicinity of $r^{*}$, we use the asymptotic expansion

[TABLE]

such that Eq. (I.3) becomes

[TABLE]

Similarly, we want to extend the solution from the opposite regime ($r\gg r^{*}$) towards the vicinity of $r^{*}$ and we use the asymptotic expansion

[TABLE]

such that Eq. (S55) becomes

[TABLE]

where $\gamma(r)=k_{F}r+\frac{\alpha^{2}-\frac{1}{4}}{2k_{F}r}-\frac{2\alpha+1}{4}\pi$. We must now match the solutions from Eqs. (I.3) and (S59) at $r\sim r^{*}$ in order to find a continuous solution.

We first match $u(r)$. Making the ansatz $B_{1}=\frac{C_{1}}{2}e^{i\kappa}$ and $B_{2}=\frac{C_{1}}{2}e^{-i\kappa}$ and matching the solutions at $r\sim r^{*}$ leads to a condition on $\psi(r)$

[TABLE]

Next, we match $v(r)$, which leads to

[TABLE]

Here, $n\in\mathbb{Z}$. Comparing Eqs. (S60) and (S61), we find that the parameter $\kappa$ is

[TABLE]

Thus, we find that at $r\sim r^{*}$ the function $\psi(r)$ is approximately

[TABLE]

However, recall that we earlier found that $\psi(r)$ is given by Eq. (S54) while constructing the BdG solution in the regime $r\gg r^{*}$. Hence, in order to have a consistent solution, we must compare the approximate solution Eq. (S63) with Eq. (S54). Since we want to understand the behavior of $\psi(r)$ in the regime where $r\sim r^{*}\gg\frac{1}{k_{F}}$, as a first approximation we can drop the term $e^{\frac{2}{k_{F}}\int_{0}^{r}dr^{\prime}\,\Delta(r^{\prime})}$. This is justified since the pairing term $\Delta(r)$ approaches 0 at least linearly and, as we are working in the regime where $r\ll\xi$, we can approximate this term as $e^{\frac{\Delta_{0}r}{k_{F}}}=e^{\frac{r}{\xi}}\sim 1$. Thus, in order to understand the behavior of $\psi(r)$ in the vicinity of $r^{*}$ we need only consider the integral

[TABLE]

We write $I=I_{1}+I_{2}$ where

[TABLE]

We can now approximate $I_{1}$ as

[TABLE]

and can further write

[TABLE]

Integrating $I_{2}$ by parts, we find

[TABLE]

We can then approximate $I_{2}$ as

[TABLE]

since we only make an exponentially small error in extending the integral (which is further suppressed by a factor of $k_{F}$) over the entire range. Thus, we approximate Eq. (S64) as

[TABLE]

In the vicinity of $r^{*}$, we hence find that the function $\psi(r)$ (Eq. (S54)) is approximately

[TABLE]

Comparing this asymptotic behavior with Eq. (S63), we find the vortex core energies

[TABLE]

where

[TABLE]

For an elementary vortex $k=1$, $l$ is a half integer, and thus with $n=0$, we reproduce the CdGM solution Caroli et al. (1964). For an MQV with $k=2$, we see that we cannot get a zero-energy solution since $n\in\mathbb{Z}$. Moreover, for any $k$, we find $|k|$ branches of vortex core states by taking the appropriate values of $n$ such that $E\ll\Delta_{0}$ since our calculation is only valid in this regime. While our method does not reproduce the detailed structure of the sub-gap states away from zero energy, it allows us to analytically estimate the spectral asymmetry, since we can extract the separation between the branches from the spectrum found above.

In order to better understand the nature of the vortex core states, we now consider the specific pairing profile used in our numerical analysis,

[TABLE]

Since the coherence length $\xi=\frac{k_{F}}{\Delta_{0}}$, we find that

[TABLE]

where $a=7\zeta(3)/\pi^{2}$ and $b=\pi/2$. From this, we then see that the mini-gap $\omega_{0}\sim\frac{\Delta_{0}}{k_{F}\xi}$ and furthermore, we find (pseudo) zero-energy states when

[TABLE]

where $\gamma=b/a\approx 1.8$. Imposing the condition $|l|\lesssim k_{F}\xi$, we see that we should restrict to $n=0,-1,-2,\dots,1-k$ which gives us exactly $k$ branches of vortex core states. Our calculation thus demonstrates that the angular momenta where the branches cross zero energy, called crossing points in the main text, are separated by an amount $\sim k_{F}\xi$, in agreement with our numerical results. Furthermore, we observe that taking a different form of the pair profile simply changes the constant $\gamma$ (e.g., for $\Delta(r)=\Delta_{0}\theta(r-\xi)$, we find $\gamma=4.5$) but does not affect the scaling of the crossing points $l_{n}$ with $k_{F}$ and $\xi$.

We note that while the spectrum for chiral states with singly quantized vortices has previously been calculated Möller et al. (2011), the method presented here easily generalizes to chiral states with MQVs.

I.4 D. Observables from BdG solutions

The particle number $\hat{N}$ and OAM $\hat{L}_{z}$ in the BCS ground state can be found by numerically diagonalizing the BdG Hamiltonian $H^{(l)}$. To obtain a finite spectrum we introduce a cutoff $M\gg 1$ on the radial quantum numbers such that $H^{(l)}$ is a $2M\times 2M$ Hermitian matrix. The eigenstates $(u,v)^{T}$ of the BdG Hamiltonian satisfy

[TABLE]

and are normalized as $\sum_{n=1}^{M}\left(|u^{(l)}_{nm}|^{2}+|v^{(l)}_{nm}|^{2}\right)=1$. Given these, we obtain

[TABLE]

where the sum over $m$ is restricted to the positive part of the energy spectrum.