Microscopic description of axisymmetric vortices in $^{3}P_{2}$ superfluids

TL;DR

This paper employs a microscopic approach to analyze quantized vortices in ${}^{3}P_{2}$ superfluids, revealing stable vortex configurations, Majorana bound states, and significant differences in core magnetization profiles.

Contribution

First microscopic study of vortices in ${}^{3}P_{2}$ superfluids using Eilenberger and BdG equations, identifying stable vortex types and Majorana states.

Findings

Multiple stable and metastable vortex configurations depending on magnetic field.

The $o$ vortex is most stable under strong magnetic fields.

Two zero-energy Majorana fermion bound states are found in the $o$ vortex core.

Abstract

We study quantized vortices in ${}^{3}P_{2}$ superfluids using a microscopic theory for the first time. The theory is based on the Eilenberger equation to determine the order parameters and the Bogoliubov-de Gennes (BdG) equation to obtain the eigenenergies and the core magnetization. Within axisymmetric vortex configurations, we find several stable and metastable vortex configurations which depend on the strength of a magnetic field, similar to a $v$ vortex and $o$ vortex in $^3$He superfluids. We demonstrate that the $o$ vortex is the most stable axisymmetric vortex in the presence of a strong magnetic field, and we find two zero-energy Majorana fermion bound states in the $o$-vortex core. We show that the profiles of the core magnetization calculated using the BdG equation are drastically different from those calculated using only the order parameter profiles known before.

Equations

$$H$$

$$H_1$$

$$H_2$$

$$T_{\alpha \beta }^{†}(\mathbf{r})$$

$$T_{\alpha \beta }(\mathbf{r})$$

$${\hat{t}}_{\alpha \beta }(-\mathrm{i}\bar{\nabla })$$

$$+{\hat{\sigma }}_{\beta }(-\mathrm{i}{\bar{\nabla }}_{\alpha })]-\frac{1}{3\sqrt{2}}{\delta }_{\alpha \beta }\hat{\mathbf{\sigma }}\cdot (-\mathrm{i}\bar{\nabla })\}.$$

$$\check{G}({\mathbf{r}}_1,{\mathbf{r}}_2;\mathrm{i}{\omega }_n)$$

$$={\int }_0^{\beta }\mathrm{d}\tau e^{\mathrm{i}{\omega }_n\tau }{\check{\tau }}_3⟨{\mathcal{T}}_{\tau }\vec{\Psi }({\mathbf{r}}_1,\tau ){\vec{\Psi }}^{†}({\mathbf{r}}_2)⟩.$$

$$-\check{1}\delta ({\mathbf{r}}_1-{\mathbf{r}}_2)=\left[\mathrm{i}{\omega }_n-{\check{H}}_{\mathrm{BdG}}({\mathbf{r}}_1)\right]{\check{\tau }}_3\check{G}({\mathbf{r}}_1,{\mathbf{r}}_2;\mathrm{i}{\omega }_n),$$

$${\check{H}}_{\mathrm{BdG}}(\mathbf{r})=\left(\begin{array}{cc}{h}_0(\mathbf{r})\hat{1}+\hat{U}(\mathbf{r})& \hat{\Delta }(\mathbf{r})\\ -{\hat{\Delta }}^{\ast }(\mathbf{r})& -h_0(\mathbf{r})\hat{1}-{\hat{U}}^{\mathrm{T}}(\mathbf{r})\end{array}\right).$$

$${\Delta }_{\alpha \beta }(\mathbf{R})$$

$$=gT\sum _n\int \frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^d}\mathrm{Tr}[{\hat{t}}_{\alpha \beta }(\mathbf{k}/k_{\mathrm{F}})\hat{F}(\mathbf{k},\mathbf{R};\mathrm{i}{\omega }_n)].$$

$$0$$

$$+[\mathrm{i}{\omega }_n{\check{\tau }}_3+\check{u}(\mathbf{R})+{\check{\sigma }}_{\Delta }({\mathbf{k}}_{\mathrm{F}},\mathbf{R}),\check{g}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)],$$

$${\check{\sigma }}_{\Delta }({\mathbf{k}}_{\mathrm{F}},\mathbf{R})=\left(\begin{array}{cc}\hat{0}& \hat{\Delta }({\mathbf{k}}_{\mathrm{F}},\mathbf{R})\\ -{\hat{\Delta }}^{†}({\mathbf{k}}_{\mathrm{F}},\mathbf{R})& \hat{0}\end{array}\right),$$

$$\check{g}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)={\oint }_{C_{\mathrm{qc}}}\frac{\mathrm{d}\xi }{\mathrm{i}\pi }\check{G}(\mathbf{k},\mathbf{R};\mathrm{i}{\omega }_n).$$

$$\hat{\Delta }({\mathbf{k}}_{\mathrm{F}},\mathbf{R})=\sum _{\alpha \beta }{\Delta }_{\alpha \beta }(\mathbf{R}){\hat{t}}_{\alpha \beta }^{\ast }({\bar{\mathbf{k}}}_{\mathrm{F}}),$$

$${\Delta }_{\alpha \beta }(\mathbf{R})=g{\nu }_{\mathrm{n}}\mathrm{i}\pi T\sum _n{⟨\mathrm{Tr}{\hat{t}}_{\alpha \beta }({\bar{\mathbf{k}}}_{\mathrm{F}})\hat{f}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)⟩}_{\mathrm{F}},$$

$$\check{g}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)=\left(\begin{array}{cc}\hat{g}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)& \hat{f}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)\\ -\hat{\bar{f}}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)& \hat{\bar{g}}({\mathbf{k}}_{\mathrm{F}},\mathbf{R};\mathrm{i}{\omega }_n)\end{array}\right).$$

$$A_{\alpha \beta }(\mathbf{R})=-\frac{{\Delta }_{\alpha \beta }(\mathbf{R})+{\Delta }_{\beta \alpha }(\mathbf{R})}{2\sqrt{2}}+\frac{{\sum }_{\gamma }{\Delta }_{\gamma \gamma }(\mathbf{R})}{3\sqrt{2}}{\delta }_{\alpha \beta }.$$

$$A(\mathbf{R})=\sum _{M=-2}^2{\gamma }_M(\rho )e^{\mathrm{i}(\kappa -M)\theta }{\Gamma }_M.$$

$$A^{\mathrm{Cyl}}(\mathbf{R})$$

$$k_{\rho }=\cos \theta k_1+\sin \theta k_2,$$

$$k_{\theta }=-\sin \theta k_1+\cos \theta k_2,$$

$${\hat{\sigma }}_{\rho }=\cos \theta {\hat{\sigma }}_1+\sin \theta {\hat{\sigma }}_2,$$

$${\hat{\sigma }}_{\theta }=-\sin \theta {\hat{\sigma }}_1+\cos \theta {\hat{\sigma }}_2.$$

$${\check{G}}_{\lambda }^{-1}({\mathbf{r}}_1,{\mathbf{r}}_2;\mathrm{i}{\omega }_n)={\check{G}}_0^{-1}({\mathbf{r}}_1,{\mathbf{r}}_2;\mathrm{i}{\omega }_n)-\lambda \check{\Sigma }({\mathbf{r}}_1,{\mathbf{r}}_2;\mathrm{i}{\omega }_n).$$

$${\mathcal{J}}_{\mathrm{sn}}=\frac{{\nu }_{\mathrm{n}}}{2}{\int }_0^1\mathrm{d}\lambda \mathrm{Sp}\left[{\check{\sigma }}_{\Delta }\left({\check{g}}_{\lambda }-\frac{1}{2}\check{g}\right)\right],$$

$$[\mathrm{i}{\omega }_n{\check{\tau }}_3+\check{u}+\lambda {\check{\sigma }}_{\Delta },{\check{g}}_{\lambda }]=0.$$

$$\frac{{\mathcal{J}}_{\mathrm{sn}}}{\Omega }$$

$${\bar{\mathfrak{J}}}_{\mathrm{sn}}(\rho )$$

$$E_{\nu }{\vec{u}}_{\nu }(\mathbf{r})$$

$${\vec{u}}_{\nu }(\mathbf{r})$$

$$\hat{\Delta }(\mathbf{r})$$

$$=\frac{1}{2}\sum _{\alpha \beta }\{A_{\alpha \beta }(\mathbf{r}),(-\mathrm{i}\bar{\nabla }{)}_{\beta }\}{\hat{\sigma }}_{\alpha }\mathrm{i}{\hat{\sigma }}_2$$

$$=\frac{1}{2}\sum _{\alpha \beta }\sum _M\{{\gamma }_M(\mathbf{r}),(-\mathrm{i}\bar{\nabla }{)}_{\beta }\}{\Gamma }_{M,\alpha \beta }{\hat{\sigma }}_{\alpha }\mathrm{i}{\hat{\sigma }}_2.$$

$$\frac{\partial }{\partial r_1}\pm \mathrm{i}\frac{\partial }{\partial r_2}=e^{\pm \mathrm{i}\theta }\left(\frac{\partial }{\partial \rho }\pm \mathrm{i}\frac{1}{\rho }\frac{\partial }{\partial \theta }\right).$$

$${\vec{u}}_{\nu =(n,\ell ,k_3)}(\rho )=\left(\begin{array}{c}{u}_{\uparrow n,\ell ,k_3}(\rho )\\ u_{\downarrow n,\ell +1,k_3}(\rho )\\ v_{\uparrow n,\ell +1-\kappa ,k_3}(\rho )\\ v_{\downarrow n,\ell -\kappa ,k_3}(\rho )\end{array}\right),$$

Full text

††preprint: aps/

I Introduction

Superfluidity and superconductivity are two of the most extraordinary states of matter. They are realized in materials or gases at low temperatures, such as metals, liquid heliums, and Bose-Einstein condensates of ultracold atomic gases. Apart from such situations realized in laboratory experiments, neutron stars offer much larger, stellar scale candidates of superfluidity and superconductivity Migdal (1960); Reichley and Downs (1969); Baym et al. (1969). The neutron density in the interior of neutron stars ranges from $10^{4}\mathrm{g/cm^{3}}$ to $10^{16}\mathrm{g/cm^{3}}$ and forms a hierarchical structure consisting of crusts and cores. Neutrons in the inner crust and outer core drip from nuclei and become a neutron fluid. The superfluidity is important to such a high density region because the temperature $T\sim 10^{6}\mathrm{K}$ is much lower than the Fermi temperature $T_{\mathrm{F}}\sim 10^{10}\mathrm{K}$ and the critical temperature $T_{\mathrm{c}}\sim 10^{8}\mathrm{K}$. The presence of Cooper pairs successfully describes rapid coolings Yakovlev et al. (1999) of neutron stars and slow relaxations Reichley and Downs (1969) after pulser glitches, i.e. phenomena in which angular momentum of neutron stars increases suddenly (see Refs. [Chamel, 2017; Haskell and Sedrakian, 2018; Sedrakian and Clark, 2019] for a recent review).

A type of the Cooper instability depends on the density. The attractive interaction is governed by the ${}^{1}S_{0}$ channel in the inner crust at lower density Migdal (1960), while in the outer core at higher density it becomes repulsive and the ${}^{3}P_{2}$ channel is dominant Tabakin (1968); Hoffberg et al. (1970); Tamagaki (1970); Takatsuka and Tamagaki (1971); Takatsuka (1972); Fujita and Tsuneto (1972); Richardson (1972); Amundsen and Ostgaard (1985); Takatsuka and Tamagaki (1993); Baldo et al. (1992); Elgarøy et al. (1996); Khodel et al. (1998); Baldo et al. (1998); Khodel et al. (2001); Zverev et al. (2003); Maurizio et al. (2014); Bogner et al. (2010); Srinivas and Ramanan (2016). The ${}^{3}P_{2}$ superfluid is the state in which the angular momenta of orbital ($L=1$) and spin ($S=1$) are aligned. The Ginzburg–Landau (GL) equation was obtained for ${}^{3}P_{2}$ superfluids Richardson (1972); Fujita and Tsuneto (1972); Sauls and Serene (1978); Muzikar et al. (1980); Sauls et al. (1982) and the ground state was found to be a nematic phase in the weak coupling regime Sauls and Serene (1978). The nematic phase has a continuous degeneracy 111 In spin-2 Bose-Einstein condensates, this is known to be related to the existence of quasi-Nambu-Goldstone modes Uchino et al. (2010). which is lifted by either magnetic field or sixth order terms in the GL free energy, and the uniaxial nematic (UN) phase is favored at zero temperature while $D_{2}$-biaxial nematic (BN) and $D_{4}$-BN phases are favored for finite and strong magnetic fields, respectively Masuda and Nitta (2016), relevant for magnetars. Low-energy excitations in ${}^{3}P_{2}$ superfluids affect the cooling process by neutrino emission Bedaque et al. (2003); Bedaque and Nicholson (2013); Bedaque and Reddy (2014); Bedaque et al. (2015); Leinson (2010a, b, c, 2011a, 2011b, 2012, 2013, 2015). The rapid cooling due to ${}^{3}P_{2}$ superfluids was studied for Cassiopeia A Heinke and Ho (2010); Shternin et al. (2011); Page et al. (2011), but a direct proof of the existence of the ${}^{3}P_{2}$ superfluidity is yet elusive. The ${}^{3}P_{2}$ superfluidity is more relevant for magnetars, in which the strength of the magnetic field reaches about $10^{15}$ G at the surface and possibly about $10^{18}$ G in the inside, because the tolerance of the spin-triplet pairing in the ${}^{3}P_{2}$ superfluidity is robust against the strong magnetic fields, in contrast to the spin-singlet pairing in the ${}^{1}S_{0}$ superfluidity fragile due to the Zeeman effect. The impact of strong magnetic fields on ${}^{3}P_{2}$ superfluid phases was studied in the GL equation Yasui et al. (2019a). The GL equation was also used for finding new exotic topological structures such as surface topological defects (boojums) Yasui et al. (2019b) and domain walls Yasui and Nitta (2020).

Since neutron stars rotate very rapidly, there exist a huge amount of quantized vortices. These vortices may play an important role on the glitches Anderson and Itoh (1975); One of scenarios of the glitches is described by unpinning of a large amount of vortices which transfers the angular momentum from the superfluid to the nonsuperfluid crust. Vortices in ${}^{3}P_{2}$ superfluids have been studied using the GL equation Richardson (1972); Fujita and Tsuneto (1972); Sauls and Serene (1978); Muzikar et al. (1980); Sauls et al. (1982); Masuda and Nitta (2016, 2020); Chatterjee et al. (2017). Because of a discrete symmetry of the condensates, we expect rich structures of vortices, such as fractional vortices and non-Abelian vortices Masuda and Nitta (2020). In addition to such nontrivial topology, spontaneous magnetization Sauls et al. (1982); Masuda and Nitta (2016, 2020); Chatterjee et al. (2017) is a crucial issue to explain the above phenomena. For further study of vortices beyond the GL equation, it is more important to take account of the fermion degrees of freedom because fermions form vortex bound states. Therefore, in this paper, we formulate a microscopic theory and calculate single vortex states with or without a magnetic field.

The microscopic model of ${}^{3}P_{2}$ superfluids was constructed long back by Richardson Richardson (1972) and Tamagaki and Takatsuka Tamagaki (1970); Takatsuka and Tamagaki (1971); Takatsuka (1972), but the first direct calculation was done recently Mizushima et al. (2017), which clarifies that nematic states of ${}^{3}P_{2}$ superfluid is a topological superfluid with time reversal symmetry (a class DIII in the classification of topological insulators and superconductors), allowing gapless Majorana fermion on the boundary. The existence of such fermion bound states is a noticeable character of the topological states. The microscopic theory also offers the phase diagram for the magnetic field and temperature and elucidates that there is a tricritical point on the phase boundary between $D_{2}$-BN and $D_{4}$-BN states Mizushima et al. (2017, 2020), which was later confirmed in the GL free energy up to the eighth order Yasui et al. (2019c). Moreover, cyclic and ferromagnetic phases are possible for total spin-2 condensates Mermin (1974), and these states have been shown to be Weyl semimetals for ${}^{3}P_{2}$ superfluids Mizushima et al. (2017); Mizushima and Nitta (2018), having gapless Weyl fermions in the bulk which may be relevant for cooling of neutron stars.

In general, for topological superconductors and superfluids, Majorana fermions may exist in the vortex core as well as their boundary. The topological aspect of fermion degrees of freedom emergent in vortices of ${}^{3}P_{2}$ superfluids has not yet been uncovered. In the superfluid 3He, which is a prototype of spin-triplet $p$-wave superfluid, all possible vortices are classified in terms of discrete symmetries preserved in vortex states Salomaa and Volovik (1987). The $o$-vortex and $v$-vortex states are the local minima of the thermodynamic potential in the superfluid 3He-B under rotation, where the former (latter) preserves (breaks) the magnetic $\pi$ rotation symmetry called the $P_{3}$ symmetry Salomaa and Volovik (1983); Passvogel et al. (1984); Thuneberg (1986); Kasamatsu et al. (2019); Regan et al. (2020). ${}^{3}P_{2}$ superfluids with vortices are categorized into the class D of the topological periodic table, and remaining discrete symmetry plays a critical role on the topological protection of the zero energy vortex-bound states Teo and Kane (2010); Shiozaki and Sato (2014). The $o$ vortex is the most symmetric vortex with spin-degenerated zero modes, which are protected by the $P_{3}$ symmetry Tsutsumi et al. (2015); Mizushima et al. (2016). In contrast, the $v$ vortex which spontaneously breaks the $P_{3}$ symmetry has no topologically protected zero modes. In the superfluid 3He-B under rotation, the $o$-vortex state is not thermodynamically stable against axisymmetric and nonaxisymmetric $v$-vortex states with no zero modes. Therefore, it is an important issue to study if a vortex with zero energy Majorana fermions is possible in ${}^{3}P_{2}$ superfluids.

The existence of fermion degrees of freedom is also important at the macroscopic level, as pointed out by Jones Jones (2009). As is known in superfluid and superconducting systems, fermion bound states in the vortex core seriously affect vortex dynamics through the spectral flow force Kopnin and Salomaa (1991); Kopnin (2002); Volovik (2014). The vortex dynamics is a key role in interpreting a gradual decrease of angular momentum of a neutron star and its glitch. Understanding the self-consistent structure of vortices and the topological protection of vortex-bound states in ${}^{3}P_{2}$ superfluids may open a door to the issues in neutron stars.

In this paper, motivated by these earlier works, we investigate vortex states and fermion bound states in the vortex core in ${}^{3}P_{2}$ superfluids using the microscopic theory. The theory is based on the Eilenberger equation to determine the order parameters and the Bogoliubov-de Gennes (BdG) equation to obtain the eigenenergies and the core magnetization. We investigate several stable and metastable vortices and discuss their stability with respect to their free energies in the presence of a magnetic field. We also clarify if Majorana fermions exist or not in the vortex core using the BdG equation, and calculate spin densities around the vortex core. We find, in the presence of axisymmetry, an $o$ vortex is stable for strong magnetic field and allows spin-degenerate (two) zero-energy Majorana fermions in its core. This finding may be important to comprehension of the cooling rate and the changes of rotating rate of neutron stars. We also show that the profiles of the core magnetization calculated using the BdG equation are drastically different from those calculated using only the order parameter profiles in the GL theory Sauls et al. (1982); Masuda and Nitta (2016, 2020); Chatterjee et al. (2017).

The remaining part of this paper consists as follows. In Sect. II, we explain the microscopic equations of the ${}^{3}P_{2}$ superfluid and its axisymmetric condition. In Sect. III, we show the numerical results: We seek for order parameter profiles and their free energy densities on the basis of the quasiclassical scheme. We also obtain the eigenenergies of fermion excitations and core magnetizations which consist of order parameter modulations and fermion bound states using the BdG equation. In Sec. IV, we provide a summery and brief discussion about nonaxisymmetric vortices. In Appendix A, we summarize basis sets of the order parameter and their matrix representations. In Appendices B and C, we review the angular momentum operators of Cooper pairs and the rotation of a basis set, respectively. In Appendix D, we compare the fermion excitation spectra obtained by solving the BdG equation and those in the quasiclassical theory.

Throughout the paper, we specify the orthonormal spatial and spin directions by 1,2, and 3, and use the notations $\hat{\cdot}$, and $\check{\cdot}$ for a 2 by 2 and a 4 by 4 matrix, respectively. Particularly, $\hat{\sigma}_{\alpha=1,2,3}$ is the $\alpha$-th component of Pauli matrices. We also set $\hbar=k_{\rm B}=1$.

II Model and Method

II.1 Gor’kov equation

First we introduce the microscopic Hamiltonian with a zero range, and attractive ${}^{3}P_{2}$ force between neutronsRichardson (1972):

[TABLE]

The second line consists of the kinetic energy $h_{0}(\bm{r})=(-\frac{\nabla^{2}}{2m}-\mu)$ with a chemical potential $\mu$ and the Zeeman energy $\hat{U}(\bm{r})=-V_{\mathrm{Z}}\hat{\sigma}_{3}$ of the magnetic field parallel to the third direction. In the third line, the interaction strength $g$ is positive, and $T^{\dagger}(T)$ is a pair creation (annihilation) operator defined by

[TABLE]

where we use the dimensionless notation $\bar{\bm{\nabla}}\equiv k_{\mathrm{F}}^{-1}\bm{\nabla}$. A 2 by 2 matrix in spin space, $\hat{t}_{\alpha\beta}$, is defined for $\alpha,\beta=1,2,3$ as

[TABLE]

We see that $[\hat{t}_{\alpha\beta}(-\mathrm{i}\bar{\bm{\nabla}})]^{*}=\hat{t}_{\alpha\beta}^{*}(\mathrm{i}\bar{\bm{\nabla}})$ and that $\hat{t}_{\alpha\beta}$ is symmetric and traceless regarding the subscripts, namely, $T_{\alpha\beta}=T_{\beta\alpha}$. We also find that $t_{\alpha\beta,\sigma\sigma^{\prime}}=t_{\alpha\beta,\sigma^{\prime}\sigma}$. We derive the Gor’kov equation for Green’s functions defined by

[TABLE]

Here $\mathcal{T}_{\tau}$ is the time-ordering operator on the imaginary axis. We have defined $\psi_{\sigma}(\bm{r},\tau)=e^{H\tau}\psi_{\sigma}(\bm{r})e^{-H\tau}$ and $\psi_{\sigma}^{\dagger}(\bm{r},\tau)=e^{H\tau}\psi_{\sigma}^{\dagger}(\bm{r})e^{-H\tau}$ and the Nambu spinor $\vec{\Psi}(\bm{r},\tau)\!=\!{}^{t}(\psi_{\uparrow}(\bm{r},\tau),\psi_{\downarrow}(\bm{r},\tau),\psi_{\uparrow}^{\dagger}(\bm{r},\tau),\psi_{\downarrow}^{\dagger}(\bm{r},\tau))$. We perform the Hartree–Fock–Bogoliubov approximation for $H_{2}$ while neglecting its contribution to the one body potential: $\braket{T_{\tau}T_{\alpha\beta}(\bm{r}_{1},\tau)\psi_{\sigma^{\prime\prime}}^{\dagger}(\bm{r}_{1},\tau)\psi_{\sigma^{\prime}}^{\dagger}(\bm{r}_{2})}\approx\braket{T_{\alpha\beta}(\bm{r}_{1},\tau)}\braket{T_{\tau}\psi_{\sigma^{\prime\prime}}^{\dagger}(\bm{r}_{1},\tau)\psi_{\sigma^{\prime}}^{\dagger}(\bm{r}_{2})}$. Defining $\Delta_{\alpha\beta}(\bm{r})=g\braket{T_{\alpha\beta}(\bm{r})}=g\braket{T_{\alpha\beta}(\bm{r},\tau)}$, we introduce the mean field by $\hat{\Delta}(\bm{r})=-\sum_{\alpha\beta}\frac{1}{2}\{\Delta_{\alpha\beta}(\bm{r}),\hat{t}_{\alpha\beta}^{*}(\mathrm{i}\bar{\bm{\nabla}})\}$. The Gor’kov equation is given by

[TABLE]

The gap equation is represented using the Green function as

[TABLE]

II.2 Eilenberger equation

In a similar way, we have a right-Gor’kov equation. The quasiclassical transport equation can be obtained by (i) subtracting the right one from the left one, (ii) integrating the equation over the single particle energy $\xi_{\bm{k}}=\frac{k^{2}}{2m}-\mu$, and (iii) retaining the contribution from the pair potentials with Fermi momentum $\bm{k}_{\mathrm{F}}$. We obtain the Eilenberger equation as

[TABLE]

where $\check{u}=V_{\mathrm{Z}}\mathrm{diag}(\hat{\sigma}_{3},\hat{\sigma}_{3}^{*})$, and

[TABLE]

The contour $C_{\mathrm{qc}}$ stands for the mean of two-contour contributions: One is the counterclockwise contour in the half upper $\xi_{\bm{k}}$ plane, and the other is the clockwise contour in the half lower $\xi_{\bm{k}}$ planeEilenberger (1968); Masaki (2019). In the quasiclassical approximation, the gap equation is reduced to

[TABLE]

where $\bar{\bm{k}}_{\mathrm{F}}=\bm{k}_{\mathrm{F}}/k_{\mathrm{F}}$ and the quasiclassical anomalous propagator $\hat{f}$ is defined as

[TABLE]

The order parameter tensor $A_{\alpha\beta}(\bm{R})$, which is defined as $\hat{\Delta}(\bm{k}_{\mathrm{F}},\bm{R})=\sum_{\alpha\beta}A_{\alpha\beta}(\bm{R})\hat{\sigma}_{\alpha}\mathrm{i}\hat{\sigma}_{2}\bar{k}_{\beta}$, can be expressed as follows:

[TABLE]

II.3 Axisymmetric condition

This tensor $A$ has a representation using the third component of the angular momentum $M$ as $A(\bm{R})=\sum_{M=-2}^{2}\gamma_{M}(\bm{R})\Gamma_{M}$, where $\gamma_{M}(\bm{R})$ is a scalar function and $\Gamma_{M}$ is a 3 by 3 tensor. The representations of $\Gamma_{M}$ are explicitly shown in Appendix A. The axisymmetric condition, which is given by $(J_{3}-\kappa)A(\bm{R})=0$ for the total angular momentum $\kappa$, reads

[TABLE]

The definitions of angular momentum operators are given in Appendix B. It is instructive to see the tensor in the cylindrical representation, in which $A=R_{\theta}A^{\mathrm{Cyl}}R_{\theta}^{\mathrm{T}}$. Note that $R_{\theta}$ denotes the rotational matrix along the third axis by angle $\theta$. Since $R_{\theta}\Gamma_{M}R_{\theta}^{\mathrm{T}}=\Gamma_{M}e^{-\mathrm{i}M\theta}$, we have

[TABLE]

where $\Gamma_{M}^{\mathrm{Cyl}}$ has the same representation as $\Gamma_{M}$ but for basis $\bm{k}_{\mathrm{F}}=k_{\rho}\bm{e}_{\rho}+k_{\theta}\bm{e}_{\theta}+k_{3}\bm{e}_{3}$, and $\hat{\bm{\sigma}}=\hat{\sigma}_{\rho}\bm{e}_{\rho}+\hat{\sigma}_{\theta}\bm{e}_{\theta}+\hat{\sigma}_{3}\bm{e}_{3}$ and

[TABLE]

The rotation of the order parameter is also discussed in Appendix C. There are several choices of $\gamma_{M}(\rho\to\infty)$, and the schematic pictures of representative cases are shown in Fig. 1, using the $d$ vectors without vorticity $\kappa$: $d_{\mu}(\bm{k}_{\mathrm{F}},\theta)=\sum_{j}\sum_{M}\gamma_{M}(\infty)e^{-\mathrm{i}M\theta}\Gamma_{M,\mu j}\bar{k}_{j}$. Panels (a)–(c) are for the UN phase, while panel (d) is for the $D_{4}$-BN phase

II.4 Free energy

We calculate the free energy on the basis of the Luttinger–Ward formalism. By solving the Eilenberger equation combined with the gap equation, we have determined the self energy $\check{\sigma}$ self-consistently. Let us define an auxiliary Green function in the Gor’kov equation as

[TABLE]

Note that $\check{G}_{\lambda=0}=\check{G}_{0}$ and $\check{G}_{\lambda=1}=\check{G}$. Following Ref. [Vorontsov and Sauls, 2003], we obtain the difference of the thermodynamic potentials between the superfluid and normal states as

[TABLE]

where $\mathrm{Sp}[\cdots]=\mathrm{i}\pi T\sum_{n}\int\mathrm{d}\bm{R}\int\frac{\mathrm{d}\hat{\bm{k}}}{4\pi}\mathrm{Tr}[\cdots]$ and $\check{g}_{\lambda}$ is the solution to the equation

[TABLE]

When the system is axisymmetric, we obtain the free-energy density as a function of $\rho$ with $\Omega=\pi R^{2}\Omega_{3}$

[TABLE]

where $\mathrm{Tr}$ in the second line stands for the trace of the 2 by 2 matrix in the spin space. We show the dimensionless total free energy $\bar{\mathcal{J}}_{\mathrm{sn}}$ and free-energy density $\bar{\mathfrak{J}}_{\mathrm{sn}}(\rho)$ in Sect. III.

II.5 Bogoliubov–de Gennes equation

The BdG equation is derived from the equation $\mathrm{i}\partial_{t}\psi_{\sigma}^{(\dagger)}=[\psi_{\sigma}^{(\dagger)},H]$ using a mean field approximation and $\vec{\Psi}(\bm{r},t)=\sum_{\nu}e^{-\mathrm{i}E_{\nu}t}\vec{u}_{\nu}(\bm{r})\alpha_{\nu}$.

[TABLE]

Here we note that

[TABLE]

We use the cylindrical coordinate $(\rho,\theta,r_{3})$ and note the following:

[TABLE]

The eigenvector can be given by $\vec{u}_{\nu}(\rho,\theta,r_{3})=\frac{e^{\mathrm{i}k_{3}r_{3}}}{\sqrt{\Omega_{3}}}\check{U}_{\ell,\kappa}(\theta)\vec{u}_{n,\ell,k_{3}}(\rho)$ with a 4 by 4 diagonal matrix $\check{U}_{\ell,\kappa}(\theta)=\frac{1}{\sqrt{2\pi}}\mathrm{diag}(e^{\mathrm{i}\ell\theta},e^{\mathrm{i}(\ell+1)\theta},e^{\mathrm{i}(\ell+1-\kappa)\theta},e^{\mathrm{i}(\ell-\kappa)\theta})$ and the radial part to be determined:

[TABLE]

because of the equations $[\check{\mathcal{H}}_{\mathrm{BdG}},J_{3}\check{1}-\frac{\kappa}{2}\check{\tau}_{3}]=0$ and $[\check{\mathcal{H}}_{\mathrm{BdG}},(-\mathrm{i}\partial_{3})\check{1}]=0$. We may use the Bessel functions and their zeros for expansion of the radial part ($w=u$ or $\varv$):

[TABLE]

Here $\beta_{\ell,k}$ denotes the $k$ th zero of the Bessel function $J_{\ell}$. The orthonormalization condition is given by

[TABLE]

We obtain the Hamiltonian matrix to be diagonalized as $[\check{\mathcal{H}}_{\mathrm{BdG}}(\ell,k_{3})]_{k,k^{\prime}}=\int_{0}^{R_{0}}\mathrm{d}\rho\rho\check{\Phi}_{\ell,k}(\rho)\check{\mathcal{H}}_{\mathrm{BdG}}(\rho,\ell,k_{3})\check{\Phi}_{\ell,k^{\prime}}(\rho)$ with $\check{\Phi}_{\ell,k}(\rho)=\mathrm{diag}(\phi_{\ell,k},\phi_{\ell+1,k},\phi_{\ell+1-\kappa,k},\phi_{\ell-\kappa,k})(\rho)$, where the eigen equation takes the following form:

[TABLE]

We remark that there is a relation between the positive and negative eigenvalues. When we fix $\kappa=1$, the relation reads the one between the state with $\nu=(n,\ell,k_{3})$ and the state with $\bar{\nu}=(\bar{n},-\ell,-k_{3})$ for some $\bar{n}$:

[TABLE]

The Bogoliubov transformation is given by

[TABLE]

We calculate the expectation value of the third component of local spin $S_{3}(\bm{r})=\frac{1}{2}(\psi_{\uparrow}^{\dagger}(\bm{r})\psi_{\uparrow}(\bm{r})-\psi_{\downarrow}^{\dagger}(\bm{r})\psi_{\downarrow}(\bm{r}))$ as

[TABLE]

where $f(E_{\nu})=\braket{\alpha_{\nu}^{\dagger}\alpha_{\nu}}=1/(e^{E_{\nu}/T}+1)$ is the Fermi distribution function.

III Numerical results

In this section, by self-consistently solving the Eilenberger Eq. (11) and gap Eq. (15) with the boundary conditions as shown in Fig. 1, we show the core structure of stable $o$ and $v$ vortices in ${}^{3}P_{2}$ superfluids. Using the BdG Eq. (29) with the gap function obtained from the quasiclassical theory, we discuss excitation spectra and magnetizations around the vortices. We set $T=0.4T_{\mathrm{c}}$ and $\omega_{\mathrm{c}}=10T_{\mathrm{c}}$ for all the calculations shown here. The units of energy and length are, respectively, $T_{\mathrm{c}}$ and $\xi_{0}=\varv_{\mathrm{F}}/(2\pi T_{\mathrm{c}})$.

Since we show several kinds of single vortices, we summarize our classification rule here. Their differences appear owing to the symmetry of order-parameter components around the core and the boundary conditions. Regarding the symmetry around the core, we construct $o$ and $v$ vortices; They are distinguished whether a $P_{3}$ symmetry exists or not. This classification is based on the context of the superfluid 3He-B phaseSalomaa and Volovik (1987), and details are discussed in the next subsection. As for the boundary conditions, we label $3$, $\rho$, $\theta$ for vortices in the UN or $D_{2}$-BN phase, which depend on the boundary conditions, while use no label for the $D_{4}$-BN phase. Let $V_{\mathrm{Zc}}$ be a transition magnetic field between the $D_{2}$-BN and $D_{4}$-BN states. For magnetic field $V_{\mathrm{Z}}<V_{\mathrm{Zc}}$, where either UN or $D_{2}$-BN state realizes, the name is determined by the direction of the maximum eigenvalue (max. EV) of the order parameter tensor. We consider three representative directions: the third, radial and angular directions. At zero magnetic field, the schematic images of the boundary conditions for these three cases are shown in Figs. 1(a)–1(c), which stand for $3$, $\rho$, and $\theta$ vortices, respectively. In the presence of magnetic field parallel to the third axis, the boundary conditions are obtained by transforming Figs. 1(b) and 1(c) continuously, though we do not explicitly show their schematic images. Note that when the magnetic field is parallel to the direction of the max. EV, such an order parameter is unstable and changes its direction of the max. EV. Therefore we do not study magnetic field effects on the $3$ vortex. For $V_{\mathrm{Z}}\geq V_{\mathrm{Zc}}$, the $D_{4}$-BN state is realized and the boundary conditions of $\rho$ and $\theta$ become the same. We name vortices in that region “$D_{4}$-BN vortex”. When we have several configurations with the same boundary condition, we just add further labels 1 and 2 to the above rules to distinguish them, e.g., $\theta 1$ and $\theta 2$ and so on.

III.1 Zero magnetic field

First we consider the zero-field case, where the equilibrium state is the UN state. The state can be characterized by the direction of the max. EV of the order parameter matrix $A$. A simple form of single vortex states is given when the max. EV points to the third direction, i.e., $A=A_{0}\mathrm{diag}(-1/2,-1/2,1)$ [Fig. 1(a)], which we call the $3$ vortex. The boundary condition is imposed at $\rho=R_{\mathrm{c}}$ as $A(R_{\mathrm{c}},\theta)=\gamma_{0}e^{\mathrm{i}\kappa\theta}\Gamma_{0}$ with the vorticity $\kappa$, where $\gamma_{0}$ is the value of the order parameter in the uniform state. We fix $\kappa=1$ and first show the possible axisymmetric vortex solutions and compare the results with those for 3He-B. In the earlier work using GL theory up to the sixth order termMasuda and Nitta (2016), the axisymmetric form of $A(\bm{R})=\gamma_{0}(\rho)e^{\mathrm{i}\kappa\theta}\Gamma_{0}$ is only considered. We investigate other solutions in the form of Eq. (18) in an analogy to the case of 3He-B.

Figures 2(a) and 2(b) show the spatial profiles of the order parameter. The axisymmetric form of $A(\bm{R})=\gamma_{0}(\rho)e^{\mathrm{i}\kappa\theta}\Gamma_{0}$ in Ref. [Masuda and Nitta, 2016] corresponds to $\gamma_{M=0}(\rho)$ in the panel (a) with $\gamma_{M=\pm 2}=0$. We remark that it does not satisfy the gap Eq. (15) since the right-hand side in Eq. (15) has nonzero components with $M=\pm 2$. Instead, we obtain the profile shown in panel (a) with nonzero induced components $\gamma_{M=\pm 2}(\rho)$, which is the so-called “$o$ vortex” in the context of 3He B-phase. The $o$ vortex is the most symmetric vortex in terms of the three discrete symmetries: They are called $P_{1}$, $P_{2}$, and $P_{3}$ symmetries with $P_{1}P_{2}P_{3}=1$. The physical meaning of these symmetries are, respectively, inversion, magnetic mirror, and magnetic $\pi$ rotation symmetries and the details are referred to as Refs. [Salomaa and Volovik, 1985, 1987; Tsutsumi et al., 2015; Mizushima et al., 2016]. The presence of these symmetries is equivalent for ${}^{3}P_{2}$ superfluids to the case in which the components of $M=\pm 2$ and [math] are real, and those of $M=\pm 1$ are zero. Perturbations with nonzero $\gamma_{\pm 1}$ change it to “$v$ vortex” shown in panel (b). The $v$ vortex is characterized by the presence of $P_{2}$ symmetry and broken $P_{1}$ ($P_{3}$) symmetry. It is represented by five real components $\gamma_{M=-2,\cdots,2}$. Since the $v$ vortex has the unwinding component, $M=1$, the vortex core is filled with $\gamma_{1}$. The power law of $\gamma_{\pm 1}$ in the asymptotic region is $1/\rho$, while $\gamma_{0,\pm 2}$ approaches to the bulk values with $1/\rho^{2}$. This asymptotic behavior is the same as that of the $v$ vortex in the 3He B-phase. The difference may appear owing to the restriction of the total angular momentum to the $J=2$ sector. The $v$ vortex in the 3He B-phase is considered to have the A-phase core since the A-phase component at $\rho=0$ is more dominant than the $\beta$-phase component. Only these two components are allowed to be nonzero at the origin, and they have opposite signs. On the other hand, in the case of ${}^{3}P_{2}$ superfluids, these two components always have the same magnitude with the same sign because the order parameter tensor is symmetric. We also calculate free energy densities of these vortices. Figure 3(a) shows the free energy densities of these vortices measured from the uniform solution, i.e., $\delta\mathfrak{J}_{\mathrm{sn}}(\rho)\to 0(\rho\to\infty)$ with $1/\rho^{2}$. The integrated values shown in its legend imply that the $v$ vortex is the most stable. It seems to gain the condensation energy at the vortex center by filling its core with $\gamma_{1}$. The $o$ vortex has less free energy owing to the condensation energy of $\gamma_{\pm 2}$ than the vortex constructed by $\gamma_{0}$ without any induced components, i.e., $\gamma_{M\neq 0}=0$, which is labeled “w/o ind.” in Fig. 3(a).

We show other possible boundary conditions of axisymmetric vortices in the absence of the magnetic field. They are more relevant for finite magnetic field because the above vortex states are unstable against the magnetic field parallel to the third axis. We consider two representative directions of the max. EVs which are the radial and angular directions [Figs. 1(b) and 1(c), respectively]. At zero field, we construct solutions for the angular direction only and we will discuss the case of the radial direction later in the presence of the magnetic field. As far as we investigated, we have found two kinds of solutions $\theta 1$ (Figs. 2(c) and 2(d)) and $\theta 2$ (Figs. 2(e) and 2(f)), each of which has vortices with a core ($o$ vortex) and without a core ($v$ vortex). We label panels (c)–(f) “$\theta 1$-o”, “$\theta 2$-o”, “$\theta 1$-v”, and “$\theta 2$-v”, respectively. In all the cases, there are structures similar to a domain wall ring in $10\lesssim\rho/\xi_{0}\lesssim 30$, and structures in the domain including the vortex center are well-described by those in Figs. 2(a) and 2(b). In the middle and right panels in Fig. 2, the order parameters in the domains including their vortex cores have opposite signs. The order parameters around the domain wall rings are also different. In the middle panels in Fig. 2, $\gamma_{0}(\rho)$ does not cross the zero values. We discuss which vortices are the most stable on the basis of the free energy. The free-energy density for each profile is shown in Fig. 3(b). We see that the free energies of the $\theta 1$ vortices (Figs. 2(c) and 2(d)) are lower than those of the $\theta 2$ vortices (Figs. 2(e) and 2(f)). For the $\theta 1$ vortices, the core region and the bulk region are gradually connected, accompanied by loss of the free energy with a gentle tail, as indicated by the red-solid, and pink-dashed curves in Fig. 3. The $\theta 2$-$o$ and $\theta 2$-$v$ vortices have bump structures around the connecting regions in their free energy densities, as indicated by the blue-solid, and light-blue-dashed curves. The difference in these structures determine the free-energy difference between the $\theta 1$ and $\theta 2$ vortices. The difference between $o$ and $v$ vortices can be seen in the vortex core, where a $v$ vortex has the condensation energy. In Fig. 3(c), the power law in the asymptotic region is $1/\rho^{2}$ and its coefficient depends on the boundary condition, i.e., it does not depend on induced components or the difference in profiles around the cores. This hydrodynamic part proportional to $1/\rho^{2}$ gives a logarithmic divergence $\log R_{\mathrm{c}}$ in the free energy $\delta\mathcal{J}_{\mathrm{sn}}$. Therefore, the 3 vortex has much lower energy than any of $\theta$ vortices in the absence of magnetic fields, although it may be unstable in the presence of a magnetic field.

Among all possible vortices, the $o$-vortex states have topologically-protected zero-energy modes bound to vortices regardless of the boundary condition at $\rho\rightarrow\infty$. Following Refs. [Tsutsumi et al., 2015; Mizushima et al., 2016; Shiozaki and Sato, 2014], we start with the semiclassical approximation where the Hamiltonian varies slowly in the real-space coordinate. The spatial modulation due to a vortex line is considered as adiabatic changes in the Hamiltonian as a function of the real-space coordinate surrounding the defect with an angle $\theta$. Then, the Hamiltonian is obtained in the base space, $(\bm{k},\theta)$, as

[TABLE]

where $\varepsilon({\bm{k}})$ is composed of the kinetic energy $h_{0}\hat{1}$ and the Zeeman energy $\hat{U}$. $\hat{\Delta}({\bm{k}},\theta)$ is the asymptotic form of the vortex order parameter at $\rho\rightarrow\infty$ that satisfies the boundary condition in Fig. 1. The $o$ vortex preserves the $P_{3}$ symmetry, that is, the magnetic $\pi$ rotation symmetry. Then, one can construct the chiral operator $\check{\Gamma}$ as a combination of the particle-hole exchange operator and $P_{3}$ operator, and the BdG Hamiltonian $\check{\mathcal{H}}({\bm{k}},\theta)$ preserves the chiral symmetry, $\{\check{\Gamma},\check{\mathcal{H}}({\bm{k}},\theta)\}=0$ for $k_{2}=k_{3}=0$. As long as the chiral symmetry is preserved, one can define the one-dimensional winding number for $\theta$ as

[TABLE]

where $w_{\rm 1d}(\theta=0)=2$ and $w_{\rm 1d}(\theta=\pi)=-2$ for the $o$-vortex state regardless of the boundary condition. Then, the topological invariant is defined as the difference of $w_{\rm 1d}(\theta)$

[TABLE]

which ensures the presence of the two zero energy states at $k_{3}=0$. Hence, a pair of zero energy states is guaranteed by the $P_{3}$ symmetry in the $o$-vortex state.

Solving the BdG equation (29), we here investigate the excitation spectra and core magnetizations for the $o$ and $v$ vortices in the cases of max. EV 3, $\theta 1$, and $\theta 2$. The effects of magnetic field are discussed in the next subsection. Here we set $k_{\mathrm{F}}\xi_{0}=4$ and $R_{\mathrm{0}}/\xi_{0}=80$. Figure 4 shows the energy spectra of quasiparticles with $k_{3}=0$. The local density of states (LDOS), which is comparable with the eigen spectra, is shown in Appendix D for $\theta 2$ vortices. We do not show the helical edge modes which appear in all the cases. Another common feature is that there are two spin degenerated Majorana zero modes with $\ell=0$ in $o$-vortex cores, while they mix up and split in $v$-vortex cores, as is known in 3He-B, and the anomalous branches cross at $\ell=\pm\ell_{\mathrm{c}}$. The presence of the zero-energy modes is consistent with topological consideration shown in Eq. (42). In addition, the cases of $o$ and $v$ vortices with max. EV along the third direction (in Figs. 4(a) and 4(b), respectively) are almost the same as those of 3He-B. Here we show the spin polarization rate using the color bar, which is defined by

[TABLE]

The Hamiltonian is block diagonalized by spin sectors $P_{\mathrm{s}}=\pm 1$ for the $o$ vortex at $k_{3}=0$. In the 3He-B phase, the spin-splitting of the vortex core bound states has been found in Refs. Silaev, 2009; Khaymovich and Silaev, 2010. On the other hand, for the $v$ vortex $P_{\mathrm{s}}$ takes neither $-1$ nor $1$ around small $\ell$ because the mixing of the spin sectors is caused by $\gamma_{1}$ and $\gamma_{-1}$ which are induced in the $v$-vortex core.

Then we see the case of max. EV-$\theta 1$, shown in the middle panels. In Fig. 4(c), there are three chiral anomalous branches with spin down crossing at $\ell=0$ and $\ell\simeq\pm 6$. The only one branch appears for the spin up component with the opposite angular velocity, which is the slope at $\ell=0$. The induced components with $\gamma_{M=\pm 1}$ gap out a pair of chiral branches with spin-up and spin-down crossing at $\ell=0$, as seen in Fig. 4(d). The other two chiral anomalous branches with mainly spin down component are present, but they do not possess topological zero modes.

In the case of max. EV-$\theta 2$, excitation spectra with small angular momentum, $|\ell|\lesssim 15$, are similar to those of max. EV-$3$. The sign of the angular velocities of two spin sectors near $\ell=0$ are the same. Topological structures of the anomalous branches are the same as those of max. EV-$\theta 1$; as $\ell$ decreases (increases), three branches (a single branch) carry(carries) spin-down (spin-up) quasiparticles from negative energy to positive energy (see Figs. 4(c) and 4(e)). In the presence of $\gamma_{M=\pm 1}$ components, two of the anomalous branches turn to two branches which cross zero of energy an even number of times. This topological structure is the same as that in Fig. 4(d).

In summary of the excitation spectra, the adiabatic pumping under a virtual process with increasing an angular momentum, two quasiholes are carried into the negative energy region as a net change in any case. In terms of spins, two spin down holes are created for max. EV-$\theta 1$ and $\theta 2$, while nothing changes for max. EV-$3$.

Next we see the polarization of spins around the vortex core. We calculate it using Eq. (39) in addition to the calculation using the formula given by

[TABLE]

The magnetization is basically zero when the system is particle-hole symmetric. Equation (44) incorporates the contribution from the slope of the density of states at the Fermi energy. The formula based on the microscopic theory Eq. (39) includes other particle-hole asymmetric contributions such as the one in the gradient expansions in the mixed representation of the order parameter. In order to calculate the spin density and the excitation spectra beyond the first order quantum correction correctly, we need to take account of the feedback effect of the magnetization profile to the order parameter profile.

Figure 5 shows the polarization of spins $\braket{S_{3}(\rho)}$ for several cases with comparison between Eqs. (39) and (44). We set an energy cut-off of the summation in Eq. (39) to $15T_{\mathrm{c}}$. In the case of $o$ vortices shown in Figs. 5(a), 5(c), and 5(e), the magnetization at the origin is always zero in the GL picture because the vortex has a normal core [Fig. 5(a)]. On the other hand, in the microscopic point of view, the quasiparticles with the angular momentum $\ell=0$ and $-1$ contribute to the spin density at the origin and thus there are finite spontaneous magnetizations as shown by red circles and in the inset of Fig. 5(a). In Figs. 5(c) and 5(e), we see finite magnetizations far from the vortex core even for the GL picture. They are attributed to the large difference between $M=2$ and $-2$ components owing to the domain structure. In addition, the magnitudes of these components are large because they are finite in the bulk region. Sign reversals for $\rho\gg\xi_{0}$ are related to the change of the magnitude between $M=2$ and $-2$. The number of sign reversals in Figs. 5(a), 5(b), and 5(c) are, respectively, 1, 0, and 2. This character is consistent with the order parameter profiles.

In the case of the $v$ vortices shown in Figs. 5(b), 5(d), and 5(f), finite magnetization occurs at the origin even in the GL picture because there is the finite order parameter as well. However the peak position is shifted in the microscopic point of view by $\rho=b_{\mathrm{c}}$, because an anomalous branch crosses at the finite angular momentum $\ell_{\mathrm{c}}\approx k_{\mathrm{F}}b_{\mathrm{c}}$.

III.2 Finite magnetic field

In this subsection, we work out effects of the magnetic field. In the 3He B-phase, there are several reports on the effects on $v$- and double core vortices on the basis of the GL theoryKasamatsu et al. (2019); Regan et al. (2020). In this paper, we restrict ourselves in the axisymmetric cases, and study the field effects on $o$ and $v$ vortices. We show the free energies as functions of the magnetic field in Figs. 6(a) and 6(b). For the magnetic field lower than the critical field between the $D_{2}$-BN and $D_{4}$-BN phases, we distinguish the species of vortices by the direction of the max. EVs: $\rho$, $\theta$, and $3$. Note that we have seen that two solutions $\theta 1$ and $\theta 2$ are possible at zero field. The free energies of $\theta 2$ vortices become lower than those of $\theta 1$ vortices for $V_{\mathrm{Z}}/T_{\mathrm{c}}\gtrsim 0.15$. This may be understood as follows: The magnetic field destabilizes the 3 vortex. An area of the core consisting of the 3 vortex reduces, namely the position of the bump structure goes inside. The loss of the free energy at the bump region decreases because the circumstance along the bump structure becomes small. On the other hand, the difference in the free-energy density decays slowly for the $\theta 1$ vortex, and it is not affected very much by changing the domain position.

We can construct $\rho$ vortices for finite $V_{\mathrm{Z}}$, which have domain structures as well as $\theta$ vortices. There are $o$- and $v$ vortices under the axisymmetry. We show the profile of the $v$ vortex at $V_{\mathrm{Z}}=0.2T_{\mathrm{c}}$ in Fig. 7(a). It has a domain structure similar to a $\theta 2$-$v$ vortex. It should be noted that the gradient energy of a $\rho$ vortex is higher than that of a $\theta$ vortex in contrast to the GL theory. This difference may be attributed to the higher order corrections in the GL expansion. In the $D_{4}$-BN phase, the boundary conditions of $\rho$- and $\theta$ vortices are equivalent, and the branches of $\rho$ and $\theta$ vortices in the state space merge into one branch. Actually Fig. 6 shows that the free energy of the $\theta 2$ vortex becomes the same as that of the $\rho$ vortex at $V_{\mathrm{Z}}=V_{\mathrm{Zc}}$. On the other hand, the $\theta 1$-$o$ and $\theta 1$-$v$ vortices are destabilized for a decay into $\theta 2$-$o$ and $\theta 2$-$v$ vortices, respectively, by applying the field. In particular, the $\theta 1$-$o$ vortex becomes unstable for $V_{\mathrm{Z}}\gtrsim 0.5T_{\mathrm{c}}$ and it is different from the $\rho$-$o$ vortex even for $V_{\mathrm{Z}}\geq V_{\mathrm{Zc}}$. We remark that structures of the excitation spectra for $\rho$ vortices are similar to those for $\theta 2$ vortices as expected from the similarity in their profiles.

We discuss the field effects in the $D_{4}$-BN phase. The axisymmetry imposes the unique boundary condition except for the global gauge transformation (see Fig. 1(d)) because the amplitude of the order parameter in the uniform state is isotropic in the momentum space. For $V_{\mathrm{Z}}\lesssim 0.9T_{\mathrm{c}}$, the solutions for $o$- and $v$ vortices can be obtained as ones connected smoothly to the $\rho$ and $\theta 2$ vortices. Their free energies are displayed in Figs. 6(a) and 6(b) by inverted triangles. We find that there are different solutions only for the $o$ vortices at $V_{\mathrm{Z}}=0.85T_{\mathrm{c}}$ and $0.9T_{\mathrm{c}}$. We call these two $o$ vortices in the $D_{4}$-BN phase “$D_{4}$-BN-$o$1” and “$D_{4}$-BN-$o$2” vortices. The former smoothly connects to $\rho$-$o$ and $\theta 2$-$o$ vortices, as shown by the red-inverted triangles in Figs. 6(a) and 6(b), while the latter is shown by blue rhombuses. For $V_{\mathrm{Z}}\gtrsim 0.9T_{\mathrm{c}}$, the $D_{4}$-BN-$o$1 vortex is destabilized and the $D_{4}$-BN-$o$2 vortex becomes the unique $o$-vortex solution. In Figs. 7(b)–7(d), we show the solutions of $D_{4}$-BN-$v$, $D_{4}$-BN-$o$1, and $D_{4}$-BN-$o$2 vortices, respectively, at $V_{\mathrm{Z}}=0.9T_{\mathrm{c}}$. We see slight differences between Figs. 7(c) and 7(d), e.g. the initial slopes of $M=2$ components, but we do not discuss it in details. Note that the global $\pi$ gauge transformation is necessary to see a smooth connection between Figs. 2(e) and 7(c) and that between Figs. 2(f) and 7(b).

From Fig. 6(b), we observe that the free-energy difference between the $o2$- and $v$ vortices decreases after discontinuous transitions from the $D_{4}$-BN-$o1$ to $D_{4}$-BN-$o2$ vortices. Finally the difference vanishes continuously at around $V_{\mathrm{Z}}\sim 1.4T_{\mathrm{c}}$; the $M=\pm 1$ components of the $v$ vortices are no longer finite, and the axisymmetric vortex recovers the $P_{1(3)}$ symmetry. The spatial profile of the order parameter at $V_{\mathrm{Z}}=1.5T_{\mathrm{c}}$ is shown in Fig. 7(e). As the Majorana fermions in the $o$-vortex core of the superfluid 3He B-phase is protected by the $P_{3}$ symmetry, we have confirmed the existence of the Majorana fermion in the $D_{4}$-BN phase within the axisymmetry as shown in Fig. 7(f). It should be emphasized that the microscopic calculation of a single vortex in multicomponent superfluids had not been done for finite Zeeman field so far even in the context of the superfluid 3He-B. We have microscopically demonstrated for the first time that the strong magnetic field actually eliminates $M=\pm 1$ components, which break the $P_{3}$ symmetry.

We give a few comment on axisymmetric vortices in the presence of the magnetic field. In the case of the superfluid 3He-B phase, the GL theory suggests that the nonaxisymmetric double-core vortex is still stable even though the $v$ vortex becomes unstable in the presence of a strong magnetic field. We also emphasize that situations may be different in the ${}^{3}P_{2}$ superfluids. In this case, the magnetic field affects components which are unaffected in the 3He-B phase, through the property of symmetric tensor, and thereby we expect that the strong magnetic field excludes the possibility of the double core vortex which is stable in the 3He-B phase. On the other hand, a single vortex in $D_{4}$-BN phase is split into two half-quantized vortices in the GL theoryMasuda and Nitta (2020) using different boundary conditions. It remains important to study a possibility of nonaxisymmetric vortices in ${}^{3}P_{2}$ superfluids and the presence of topological zero modes on the basis of microscopic theory.

IV Summary and discussion

We have studied axisymmetric vortices in ${}^{3}P_{2}$ superfluids, using microscopic theory: the Eilenberger equation to determine the order parameters and the BdG equation to study the eigenenergies and the core magnetization. We have found that several features as a multicomponent superfluid are common to the superfluid 3He-B phase, though they are overlooked in the GL theory for ${}^{3}P_{2}$ superfluids, e.g. the existence of the $v$ vortex. We have shown that the profiles of the core magnetization calculated using the BdG equation are drastically different from those calculated using only the order parameter profiles calculated in the GL theory. We have demonstrated that the $o$ vortex is the most stable axisymmetric vortex in the presence of a strong magnetic field, and have found two zero energy Majorana fermions in the $o$-vortex core. This observation is based on the microscopic calculation in the presence of the magnetic field for multicomponent superfluids and the $P_{3}$ symmetry argument.

One of open questions is whether two Majorana zero modes in the $o$ vortex give a nontrivial non-Abelian statistics, in contrast to the conventional case of one Majorana fermion zero mode in a vortex resulting in non-Abelian statistics among vortices Ivanov (2001), which can be generalized to odd numbers of Majorana fermionsYasui et al. (2011); Hirono et al. (2012); Eto et al. (2014). In the presence of the mirror symmetry, a possibility of the non-Abelian statistics of spin-degenerated Majorana zero modes are consideredSato et al. (2014).

The case of a nonaxisymetric case is one of the most important extensions of the present work. Although the quadrupole deformation leads to the nonaxisymmetric double-core vortex in the 3He-B phase, the order parameter tensor of ${}^{3}P_{2}$ superfluids is restricted to symmetric traceless one and it is nontrivial whether the same deformation may occur or not. 222 The effects of the magnetic field on a double-core vortex is not trivial. In the case of 3He-B phase, the $\beta$ phase component does not conflict with the magnetic fieldKasamatsu et al. (2019). However, the components of the $\beta$ phase and A phase are always the same in ${}^{3}P_{2}$ superfluids, and thus the induced components of the double core vortex might be expected to reduce as well as $v$ vortex.

Thus we should take account of a general symmetric deformation of an $o$ vortex to a nonaxisymmetric vortex also in terms of the protection of Majorana fermionsTokuyasu et al. (1990); Kobayashi et al. (2009a); Klein et al. (2014); Silaev et al. (2015). It is also important whether 1/2 quantized non-Abelian vortices in the $D_{4}$-BN phase Masuda and Nitta (2020) admit zero-energy fermions in their cores, and if so whether it may give a novel non-Abelian statistics.

In this paper, we have focused on vortices in the nematic phases, which are the most stable ${}^{3}P_{2}$ superfluid phases in the weak coupling regime under no rotation. However, neutron stars are dense neutron matter under extreme conditions, such as rapid rotation. Strong coupling effect and rapid rotation might essentially change the superfluid phase diagram, and make the cyclic phase and the ferromagnetic phase competitive to the nematic phasesMizushima et al. (2017); Mizushima and Nitta (2018). In particular, it has been found that the cyclic phase admits the 1/3 quantized non-Abelian vorticesSemenoff and Zhou (2007); Kobayashi et al. (2009b). It is important to explore the microscopic structures of vortices and existence of Majorana zero modes in the cyclic and the ferromagnetic phase in ${}^{3}P_{2}$ superfluids.

Finally, applications of fermion zero modes to neutron star physics such as contribution to cooling rate and vortex dynamics remain as an important future workJones (2009).

Acknowledgements.

The authors thank Michikazu Kobayashi and Shigehiro Yasui for useful discussions. Y.M. thanks Robert Regan for useful comments. This work was supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT)-Supported Program for the Strategic Research Foundation at Private Universities “Topological Science” (Grant No. S1511006). This work was also supported in part by the Japan Society for the Promotion of Science (JSPS) Early-Career Scientists Grant No. JP19K14662 (Y.M.), JSPS Grant-in-Aid for Scientific Research (KAKENHI) [Grants No. JP16K05448 (T.M.), No. JP16H03984 (M.N.), and No. JP18H01217 (M.N.)], and JSPS KAKENHI on Innovative Areas “Topological Materials Science” (Grant No. JP15H05855).

Appendix A Representation of Cooper pairs

In this Appendix, we summarize the basis for order parameter and matrix representations. We consider the two basis sets for $A_{\mu j}$. In one basis, $\mu$ and $j$ take $1,\cdots,3$, and in the other basis they take $+,0,-$. We distinguish them by calling Cartesian (Car) and Angular Momentum (AM), respectively. We denote them by $A^{\mathrm{Car}}$ and $A^{\mathrm{AM}}$. Here the relations between two basis sets are given by

[TABLE]

with $\hat{\tau}_{\mu=1,2,3}\equiv\mathrm{i}\hat{\sigma}_{\mu=1,2,3}\hat{\sigma}_{2}$. The matrix for the transformation is given by

[TABLE]

with

[TABLE]

We review the decomposition of triplet $p$-wave pairing with respect to the total angular momentum. Cooper pair has $L=1$ and $S=1$, and thus the total angular momentum $\bm{J}=\bm{L}+\bm{S}$ takes $J=0,1,2$: $(L=1)\otimes(S=1)=(J=0)\oplus(J=1)\oplus(J=2)$. The order parameter tensor $A_{\mu\nu}^{\mathrm{Car}}$, which is defined by $d_{\mu}=A_{\mu\nu}^{\mathrm{Car}}\bar{k}_{\nu}$, satisfies the following properties:

[TABLE]

Next we see the decomposition of traceless and symmetric tensor, $A^{J=2}$, into the irreducible representations: $\Gamma_{M=-2,\cdots,2}$: $A=\sum_{M}\gamma_{M}\Gamma_{M}$. The basis $\Gamma_{M}$ in Cartesian representation is defined so that

[TABLE]

Therefore we obtain in Cartesian representation

[TABLE]

The basis set in terms of the order parameter is explicitly written down, through $\hat{\Delta}_{M}(\bm{k})=\Gamma_{M,\mu i}\bar{k}_{i}\hat{\tau}_{\mu}$, as

[TABLE]

The orthonormal basis sets in momentum and spin spaces are given, respectively, by

[TABLE]

where $\braket{\cdot}{\cdot}_{k}=\frac{3}{4\pi}\int\mathrm{d}\cos\theta\mathrm{d}\phi\cdots$, and $\braket{\cdot}{\cdot}_{s}=\mathrm{Tr}{\cdots}$. For these notations, $\hat{\Delta}_{M}(\bm{k})$ are represented as

[TABLE]

In angular momentum representation, $\Gamma_{M}$ are described as

[TABLE]

Appendix B Angular momentum

Here we summarize the definitions of angular momentum operators for Cooper pairs, which are given by

[TABLE]

The antisymmetric tensor in the last term can be regarded as the matrix acting on the right subscript of $A$. The spin angular momentum operator is given by

[TABLE]

The total angular momentum is given by $J=L+S=L^{\mathrm{ext}}+J^{\mathrm{int}}$. For the irreducible representations, the $3$rd component of the angular momentum $J_{3}^{\mathrm{int}}$ can be calculated as

[TABLE]

Appendix C Rotation of triad

In this Appendix, we discuss rotations of a basis set of the order parameter for intuitive understanding of our boundary conditions. Let $\widehat{u}$, $\widehat{v}$, $\widehat{w}$ be the triad. Here $\widehat{w}=\widehat{3}$, and $\widehat{u}$, $\widehat{v}$ are obtained by rotation $\widehat{1}$ and $\widehat{2}$ about $\widehat{3}$. Tensors discussed above are described using $\widehat{1}$, $\widehat{2}$, $\widehat{3}$. For example, in the Cartesian representation, $\Gamma_{M=2,\mu\nu}^{\mathrm{Car}}=(\widehat{1}_{\mu}\widehat{1}_{\nu}-\widehat{2}_{\mu}\widehat{2}_{\nu}+\mathrm{i}\widehat{1}_{\mu}\widehat{2}_{\nu}+\mathrm{i}\widehat{2}_{\mu}\widehat{1}_{\nu})/2$. We consider the following triad:

[TABLE]

Instead, we can write $(\widehat{u},\widehat{v},\widehat{w})=(\widehat{1},\widehat{2},\widehat{3})R_{\varphi}$, where $R_{\varphi}$ is a rotation matrix around $\widehat{3}$ by angle $\varphi$. Let us summarize this kind of rotation at first. Note that we take the third component of $L$ and $S$ in the same direction $\widehat{3}$. The rotation operator in real vector (e.g. $\bm{k}$, $\bm{R}$) is performed using $U_{L}(\varphi)=\exp[-\mathrm{i}\varphi L_{3}]$. [$U_{L}(\varphi):(\rho,\theta,r_{3})\to(\rho,\theta-\varphi,r_{3})$]. The operation on a function in the Wigner representation $\psi(\bm{k},\bm{R})$ is given by $U_{L}(\varphi)\psi(\bm{k},\bm{R})=\psi(R_{\varphi}^{-1}\bm{k},R_{\varphi}^{-1}\bm{R})$. Therefore, the order parameter is transformed as

[TABLE]

The simultaneous rotation in real and spin spaces is performed using $U_{J}(\varphi)=\exp(-\mathrm{i}\varphi J_{3})$ as

[TABLE]

Then we can check $R_{\varphi}\Gamma_{M}R_{\varphi}^{T}=e^{-\mathrm{i}M\varphi}\Gamma_{M}$.

Noting $\hat{1}_{i}=\delta_{1i}$, we see that

[TABLE]

This reads

[TABLE]

Therefore the rotation of the triad is described by that of $\Gamma$.

Appendix D Comparison between local density of states and eigen spectum

The angular momentum of the quasiparticle state is related to the impact parameter in the quasiclassical theory: $|\ell|=k_{\mathrm{F}}b$, where the Bohr–Sommerfeld quantization is necessary for the discreteness of $\ell$. Therefore, the relation of the bound state spectra with the angular momentum is also approximated by the peak positions of the LDOS. The LDOS can be calculated through the retarded quasiclassical Green functions as

[TABLE]

where $\bm{k}_{\mathrm{F}}=k_{\mathrm{F}}(\cos\alpha\hat{e}_{1}+\sin\alpha\hat{e}_{2})$, and $\eta$ is an infinitesimal positive value set to $0.02T_{\mathrm{c}}$. We only consider $k_{3}=0$, and compare the LDOS with the eigen spectra obtained by solving the BdG equation for $\theta 2$-vortices. The eigen spectra are again shown in Figs. 8(a) and 8(b) for the $o$ vortex and the $v$ vortex, respectively by changing the horizontal axis to the aboslute value of the angular momentum $|\ell|$ divided by the quasiclassical parameter $k_{\mathrm{F}}\xi_{0}=4$ for comparison. In Fig. 8, the top panels are for the $o$ vortex and the bottom panels are for the $v$ vortex. We show the total LDOS given by Eq. (77) in Figs. 8(c) and 8(d). The contribution from each spin component is also given in panels (e)–(f). In the case of the $o$ vortex, the spin-up and the spin-down sector are not coupled. Hence we see the correspondence between the red (blue) symbols in panel (a) and panel (e) (panel (g)). On the other hand, in the case of the $v$ vortex, the spin sectors are coupled near the vortex core owing to the induced components $\gamma_{\pm 1}(\rho)$. The effect appears in $P_{\mathrm{s}}$ in Fig. 8(b). Correspondingly, the spin dependent LDOS for the $v$ vortex $\nu_{\uparrow}$ and $\nu_{\downarrow}$ have the intensity on the same energy near the vortex core $b\lesssim 5\xi_{0}$, while they have different branches for $b\gtrsim 10\xi_{0}$ since $P_{\mathrm{s}}\simeq\pm 1$ in panel (b). After all, for both the $o$ vortex and the $v$ vortex, the consistency between the spectra and the LDOS are confirmed.