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

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.


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 [1][2][3] . The neutron density in the interior of neutron stars ranges from 10 4 g/cm 3 to 10 16 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 ∼ 10 6 K is much lower than the Fermi temperature T F ∼ 10 10 K and the critical temperature T c ∼ 10 8 K. The presence of Cooper pairs successfully describes rapid coolings 4 of neutron stars and slow relaxations 2 after pulser glitches, i.e. phenomena in which angular momentum of neutron stars increases suddenly (see Refs. [5][6][7] 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 1 , while in the outer core at higher density it becomes repulsive and the 3 P 2 channel is dominant [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] . 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 13,14,[26][27][28] and the ground state was found to be a nematic phase in the weak coupling regime 26 . The nematic phase has a continuous degeneracy 29 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 30 , relevant for magnetars. Low-energy excitations in 3 P 2 superfluids affect the cooling process by neutrino emission [31][32][33][34][35][36][37][38][39][40][41][42] . The rapid cooling due to 3 P 2 superfluids was studied for Cassiopeia A [43][44][45] , 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 46 . The GL equation was also used for finding new exotic topological structures such as surface topological defects (boojums) 47 and domain walls 48 .
Since neutron stars rotate very rapidly, there exist a huge amount of quantized vortices. These vortices may play an important role on the glitches 49 ; 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 13,14,[26][27][28]30,50,51 . Because of a discrete symmetry of the condensates, we expect rich structures of vortices, such as fractional vortices and non-Abelian vortices 50 . In addition to such nontrivial topology, spontaneous magnetization 28,30,50,51 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 14 and Tamagaki and Takatsuka [10][11][12] , but the first direct calculation was done recently 52 , 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 52,53 , which was later confirmed in the GL free energy up to the eighth order 54 . Moreover, cyclic and ferromagnetic phases are possible for total spin-2 condensates 55 , and these states have been shown to be Weyl semimetals for 3 P 2 superfluids 52,56 , 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 3 He, which is a prototype of spin-triplet p-wave superfluid, all possible vortices are classified in terms of discrete symmetries preserved in vortex states 57 . The o-vortex and v-vortex states are the local minima of the thermodynamic potential in the superfluid 3 He-B under rotation, where the former (latter) preserves (breaks) the magnetic π rotation symmetry called the P 3 symmetry 58-62 . 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 63,64 . The o vortex is the most symmetric vortex with spin-degenerated zero modes, which are protected by the P 3 symmetry 65,66 . In contrast, the v vortex which spontaneously breaks the P 3 symmetry has no topologically protected zero modes. In the superfluid 3 He-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 67 . As is known in superfluid and superconducting systems, fermion bound states in the vortex core seriously affect vortex dynamics through the spectral flow force [68][69][70] . 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 pres-ence 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 28,30,50,51 .
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·, anď · for a 2 by 2 and a 4 by 4 matrix, respectively. Particularly, σ α=1,2,3 is the α-th component of Pauli matrices. We also set = k B = 1.

B. 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 ξ k = k 2 2m − µ, and (iii) retaining the contribution from the pair potentials with Fermi momentum k F . We obtain the Eilenberger equation as The contour C qc stands for the mean of two-contour contributions: One is the counterclockwise contour in the half upper ξ k plane, and the other is the clockwise contour in the half lower ξ k plane 71,72 . In the quasiclassical approximation, the gap equation is reduced tô wherek F = k F /k F and the quasiclassical anomalous propagatorf is defined aš The order parameter tensor A αβ (R), which is defined aŝ ∆(k F , R) = αβ A αβ (R)σ α iσ 2kβ , can be expressed as follows:

C. Axisymmetric condition
This tensor A has a representation using the third component of the angular momentum M as A(R) = 2 M=−2 γ M (R)Γ M , where γ M (R) is a scalar function and Γ M is a 3 by 3 tensor. The representations of Γ M are explicitly shown in Appendix A. The axisymmetric condition, which is given by (J 3 − κ)A(R) = 0 for the total angular momentum κ, reads 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 θ A Cyl R T θ . Note that R θ denotes the rotational matrix along the third axis by angle θ. Since where Γ Cyl M has the same representation as Γ M but for basis k F = k ρ e ρ + k θ e θ + k 3 e 3 , andσ =σ ρ e ρ +σ θ e θ +σ 3 e 3 and The rotation of the order parameter is also discussed in Appendix C. There are several choices of γ M (ρ → ∞), and the schematic pictures of representative cases are shown in 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σ self-consistently. Let us define an auxiliary Green function in the Gor'kov equation aš Note thatǦ λ=0 =Ǧ 0 andǦ λ=1 =Ǧ. Following Ref. [73], we obtain the difference of the thermodynamic potentials between the superfluid and normal states as where Sp[· · · ] = iπT n ∫ dR ∫ dk 4π Tr[· · · ] andǧ λ is the solution to the equation When the system is axisymmetric, we obtain the free-energy density as a function of ρ with Ω = πR 2 Ω 3 where 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 energyJ sn and free-energy densityJ sn (ρ) in Sect. III.

E. Bogoliubov-de Gennes equation
The BdG equation is derived from the equation i∂ t ψ ( †) σ = [ψ ( †) σ , H] using a mean field approximation and ì Ψ( Here we note that We use the cylindrical coordinate (ρ, θ, r 3 ) and note the following: The eigenvector can be given by ì diag(e iℓθ , e i(ℓ+1)θ , e i(ℓ+1−κ)θ , e i(ℓ−κ)θ ) and the radial part to be determined: because of the equations We may use the Bessel functions and their zeros for expansion of the radial part (w = u or ): Here β ℓ,k denotes the k th zero of the Bessel function J ℓ . The orthonormalization condition is given by We obtain the Hamiltonian matrix to be diagonalized where the eigen equation takes the following form: We remark that there is a relation between the positive and negative eigenvalues. When we fix κ = 1, the relation reads the one between the state with ν = (n, ℓ, k 3 ) and the state with ν = (n, −ℓ, −k 3 ) for somen: The Bogoliubov transformation is given by We calculate the expectation value of the third component of local spin where f (E ν ) = α † ν α ν = 1/(e E ν /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 c and ω c = 10T c for all the calculations shown here. The units of energy and length are, respectively, T c and ξ 0 = F /(2πT 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 3 He-B phase 57 , and details are discussed in the next subsection. As for the boundary conditions, we label 3, ρ, θ 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 Zc be a transition magnetic field between the D 2 -BN and D 4 -BN states. For magnetic field V Z < V 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, ρ, and θ 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 Z ≥ V Zc , the D 4 -BN state is realized and the boundary conditions of ρ and θ 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., θ1 and θ2 and so on.

A. 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 diag(−1/2, −1/2, 1) [ Fig. 1(a)], which we call the 3 vortex. The boundary condition is imposed at ρ = R c as A(R c , θ) = γ 0 e iκθ Γ 0 with the vorticity κ, where γ 0 is the value of the order parameter in the uniform state. We fix κ = 1 and first show the possible axisymmetric vortex solutions and compare the results with those for 3 He-B. In the earlier work using GL theory up to the sixth order show the spatial profiles of the order parameter. The axisymmetric form of A(R) = γ 0 (ρ)e iκθ Γ 0 in Ref. [30] corresponds to γ M=0 (ρ) in the panel (a) with γ M=±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 = ±2. Instead, we obtain the profile shown in panel (a) with nonzero induced components γ M=±2 (ρ), which is the so-called "o vortex" in the context of 3 He 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 π rotation symmetries and the details are referred to as Refs. [57, 65, 66, and 74]. The presence of these symmetries is equivalent for 3 P 2 superfluids to the case in which the components of M = ±2 and 0 are real, and those of M = ±1 are zero. Perturbations with nonzero γ ±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 γ M=−2, ···,2 . Since the v vortex has the unwinding component, M = 1, the vortex core is filled with γ 1 . The power law of γ ±1 in the asymptotic region is 1/ρ, while γ 0,±2 approaches to the bulk values with 1/ρ 2 . This asymptotic behavior is the same as that of the v vortex in the 3 He 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 3 He B-phase is considered to have the A-phase core since the A-phase component at ρ = 0 is more dominant than the β-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., δJ sn (ρ) → 0(ρ → ∞) with 1/ρ 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 γ 1 . The o vortex has less free energy owing to the condensation energy of γ ±2 than the vortex constructed by γ 0 without any induced components, i.e., γ M 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 θ1 (Figs. 2(c) and 2(d)) and θ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) "θ1-o", "θ2-o", "θ1-v", and "θ2-v", respectively. In all the cases, there are structures similar to a domain wall ring in 10 ρ/ξ 0 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, γ 0 (ρ) 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 θ1 vortices (Figs. 2(c) and 2(d)) are lower than those of the θ2 vortices (Figs. 2(e) and 2(f)). For the θ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 θ2-o and θ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 θ1 and θ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/ρ 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/ρ 2 gives a logarithmic divergence log R c in the free energy δJ sn . Therefore, the 3 vortex has much lower energy than any of θ 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 ρ → ∞. Following Refs. [64][65][66], 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 θ. Then, the Hamiltonian is obtained in the base space, (k, θ), aš where ε(k) is composed of the kinetic energy h 01 and the Zeeman energyÛ.∆(k, θ) is the asymptotic form of the vortex order parameter at ρ → ∞ that satisfies the boundary condition in Fig. 1. The o vortex preserves the P 3 symmetry, that is, the magnetic π rotation symmetry. Then, one can construct the chiral operatorΓ as a combination of the particle-hole exchange operator and P 3 operator, and the BdG Hamiltoniaň H (k, θ) preserves the chiral symmetry, {Γ,Ȟ (k, θ)} = 0 for k 2 = k 3 = 0. As long as the chiral symmetry is preserved, one can define the one-dimensional winding number for θ as where w 1d (θ = 0) = 2 and w 1d (θ = π) = −2 for the ovortex state regardless of the boundary condition. Then, the topological invariant is defined as the difference of w 1d (θ) 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, θ1, and θ2. The effects of magnetic field are discussed in the next subsection. Here we set k F ξ 0 = 4 and R 0 /ξ 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 θ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 ℓ = 0 in o-vortex cores, while they mix up and split in v-vortex cores, as is known in 3 He-B, and the anomalous branches cross at ℓ = ±ℓ 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 3 He-B. Here we show the spin polarization rate using the color bar, which is defined by The Hamiltonian is block diagonalized by spin sectors P s = ±1 for the o vortex at k 3 = 0. In the 3 He-B phase, the spin-splitting of the vortex core bound states has been found in Refs. 75 and 76. On the other hand, for the v vortex P s takes neither −1 nor 1 around small ℓ because the mixing of the spin sectors is caused by γ 1 and γ −1 which are induced in the v-vortex core.
Then we see the case of max. EV-θ1, shown in the middle panels. In Fig. 4(c), there are three chiral anomalous branches with spin down crossing at ℓ = 0 and ℓ ≃ ±6. The only one branch appears for the spin up component with the opposite angular velocity, which is the slope at ℓ = 0. The induced components with γ M=±1 gap out a pair of chiral branches with spin-up and spin-down crossing at ℓ = 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-θ2, excitation spectra with small angular momentum, |ℓ| 15, are similar to those of max. EV-3. The sign of the angular velocities of two spin sectors near ℓ = 0 are the same. Topological structures of the anomalous branches are the same as those of max. EV-θ1; as ℓ 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 γ M=±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-θ1 and θ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 13 13 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 S 3 (ρ) for several cases with comparison between Eqs. (39) and (44). We set an energy cut-off of the summation in Eq. (39) to 15T 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 ℓ = 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 ρ ≫ ξ 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 ρ = b c , because an anomalous branch crosses at the finite angular momentum ℓ c ≈ k F b c .

B. Finite magnetic field
In this subsection, we work out effects of the magnetic field. In the 3 He B-phase, there are several reports on the effects on v-and double core vortices on the basis of the GL theory 61,62 . 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: ρ, θ, and 3. Note that we have seen that two solutions θ1 and θ2 are possible at zero field. The free energies of θ2 vortices become lower than those of θ1 vortices for V Z /T c 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 θ1 vortex, and it is not affected very much by changing the domain position. We can construct ρ vortices for finite V Z , which have domain structures as well as θ vortices. There are o-and v vortices under the axisymmetry. We show the profile of the v vortex at V Z = 0.2T c in Fig. 7(a). It has a domain structure similar to a θ2-v vortex. It should be noted that the gradient energy of a ρ vortex is higher than that of a θ 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 ρand θ vortices are equivalent, and the branches of ρ and θ vortices in the state space merge into one branch. Actually Fig. 6 shows that the free energy of the θ2 vortex becomes the same as that of the ρ vortex at V Z = V Zc . On the other hand, the θ1-o and θ1-v vortices are destabilized for a decay into θ2-o and θ2-v vortices, respectively, by applying the field. In particular, the θ1-o vortex becomes unstable for V Z 0.5T c and it is different from the ρ-o vortex even for V Z ≥ V Zc . We remark that structures of the excitation spectra for ρ vortices are similar to those for θ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. 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 Z ∼ 1.4T c ; the M = ±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 Z = 1.5T c is shown in Fig. 7(e). As the Majorana fermions in the o-vortex core of the superfluid 3 He 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 3 He-B. We have microscopically demonstrated for the first time that the strong magnetic field actually eliminates M = ±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 3 He-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 3 He-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 3 He-B phase. On the other hand, a single vortex in D 4 -BN phase is split into two halfquantized vortices in the GL theory 50 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.  Fig. 1(b), and (b) -(d) the cases of v-, o1-and o2 vortices in D 4 -BN phase at V Z /T c = 0.9, respectively. We also show the order parameter profile and excitation spectra at V Z = 1.5T c in panels (e) and (f), respectively. The profile in the vicinity of the vortex center is shown in the inset of each panel.

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 3 He-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 77 , which can be generalized to odd numbers of Majorana fermions [78][79][80] . In the presence of the mirror symmetry, a possibility of the non-Abelian statistics of spin-degenerated Majorana zero modes are considered 81 .
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 3 He-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. 82 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 fermions [83][84][85][86] . It is also important whether 1/2 quantized non-Abelian vortices in the D 4 -BN phase 50 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 phases 52,56 . In particular, it has been found that the cyclic phase admits the 1/3 quantized non-Abelian vortices 87,88 . 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 work 67 .
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 J = L + S takes J = 0, 1, 2: (L = 1) ⊗ (S = 1) = (J = 0) ⊕ (J = 1) ⊕ (J = 2). The order parameter tensor A Car µν , which is defined by d µ = A Car µνkν , satisfies the following properties: A T = A and TrA = 0 for J = 2.
Instead, we can write ( u, v, w) = ( 1, 2, 3)R ϕ , where R ϕ is a rotation matrix around 3 by angle ϕ. Let us summarize this kind of rotation at first. Note that we take the third component of L and S in the same direction 3. The rotation operator in real vector (e.g. k, R) is performed using U L (ϕ) = exp[−iϕL 3 ].

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: |ℓ| = k F b, where the Bohr-Sommerfeld quantization is necessary for the discreteness of ℓ. 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 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 γ ±1 (ρ). The effect appears in P s in Fig. 8(b). Correspondingly, the spin dependent LDOS for the v vortex ν ↑ and ν ↓ have the intensity on the same energy near the vortex core b 5ξ 0 , while they have different branches for b 10ξ 0 since P s ≃ ±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.