Topological defects at the boundary of neutron $^{3}P_{2}$ superfluids in neutron stars

TL;DR

This paper investigates the surface effects and topological defects of neutron $^{3}P_{2}$ superfluids in neutron stars, revealing surface phase variations, anisotropic energy distributions, and the existence of vortices called boojums.

Contribution

It provides a detailed analysis of boundary phenomena and topological defects in neutron $^{3}P_{2}$ superfluids using Ginzburg-Landau theory, highlighting new surface properties and defect structures.

Findings

Surface phase can differ from bulk phase.

Surface energy density is anisotropic.

Existence of topological vortices (boojums) on the surface.

Abstract

We study surface effects of neutron $^{3}P_{2}$ superfluids in neutron stars. $^{3}P_{2}$ superfluids are in uniaxial nematic (UN), D$_{2}$ biaxial nematic (BN), or D$_{4}$ BN phase, depending on the strength of magnetic fields from small to large. We suppose a neutron $^{3}P_{2}$ superfluid in a ball with a spherical boundary. Adopting a suitable boundary condition for $^{3}P_{2}$ condensates, we solve the Ginzburg-Landau equation to find several surface properties for the neutron $^{3}P_{2}$ superfluid. First, the phase on the surface can be different from that of the bulk, and symmetry restoration or breaking occurs in general on the surface. Second, the distribution of the surface energy density has an anisotropy depending on the polar angle in the sphere, which may lead to the deformation of the geometrical shape of the surface. Third, the order parameter manifold induced on the surface, which is described by two-dimensional vector fields induced on the surface from the condensates, allows topological defects (vortices) on the surface, and there must exist such defects even in the ground state thanks to the Poincar\'{e}-Hopf theorem: although the numbers of the vortices and antivortices depend on the bulk phases, the difference between them is topologically invariant (the Euler number $\chi=2$) irrespective of the bulk phases. These vortices, which are not extended to the bulk, are called boojums in the context of liquid crystals and helium-3 superfluids. The surface properties of the neutron $^{3}P_{2}$ superfluid found in this paper may provide us useful information to study neutron stars.

Figures

Tables

Table 1. Table 1: The classification of nematic phases (a table taken from Ref.  Masuda and Nitta ( 2016 ) ). We show the range of r 𝑟 r , phases, unbroken symmetries H 𝐻 H , and order parameter manifolds M ≃ G / H similar-to-or-equals 𝑀 𝐺 𝐻 M\simeq G/H .

$r$	Phase	$H$	$M\simeq G/H$
$-1/2$	UN	$O\left(2\right)$	$U\left(1\right)\times \left[\mathrm{SO}\left(3\right)/O\left(2\right)\right]$
	$D_2$-BN	$D_2$	$U\left(1\right)\times \left[\mathrm{SO}\left(3\right)/D_2\right]$
$-1$	$D_4$-BN	$D_4$	$\left[U\left(1\right)\times \mathrm{SO}\left(3\right)\right]/D_4$

Table 2. Table 2: The averaged (dimensionless) surface energy density σ av subscript 𝜎 av \sigma_{\mathrm{av}} ( σ ~ av subscript ~ 𝜎 av \tilde{\sigma}_{\mathrm{av}} ) for the bulk phases (the UN phase, the D 2 -BN phase, and the D 4 -BN phase). E surf = 4 ​ π ​ R 2 ​ σ av subscript 𝐸 surf 4 𝜋 superscript 𝑅 2 subscript 𝜎 av E_{\mathrm{surf}}=4\pi R^{2}\sigma_{\mathrm{av}} is the surface energy in total for a neutron star with the radius R = 10 𝑅 10 R=10 km. The energy unit conversion 1 ​ keV ≈ 2.0 × 10 − 9 ​ erg 1 keV 2.0 superscript 10 9 erg 1\,\mathrm{keV}\approx 2.0\times 10^{-9}\,\mathrm{erg} is used. The last row indicates the geometrical symmetries for the spatial distribution of the energy density on the surface.

bulk phase	UN	D2-BN	D4-BN
${\tilde{\sigma }}_{\mathrm{av}}$	0.0095	0.0099	0.0085
${\sigma }_{\mathrm{av}}$ [keV/fm2]	0.26	0.27	0.23
$E_{\mathrm{surf}}$ [erg]	$6.5\times {10}^{29}$	$6.8\times {10}^{29}$	$5.8\times {10}^{29}$
geometrical sym.	D4	D2	D4

Table 3. Table 3: The positive and negative charges χ ± subscript 𝜒 plus-or-minus \chi_{\pm} for the vortices on the plane for the bulk phases (the UN phase, the D 2 -BN phase, and the D 4 -BN phase). We define χ = χ + + χ − 𝜒 subscript 𝜒 subscript 𝜒 \chi=\chi_{+}+\chi_{-} to indicate the sum of the positive charges and the negative charges in total. See also the bottom row in Fig.  4 .

bulk phase	UN	D2-BN	D4-BN
positive charge (${\chi }_{+}$)	10	12	8
negative charge (${\chi }_{-}$)	-8	-10	-6
total charge ($\chi$)	2	2	2

Equations

$$f[A]=f_0+f_{\mathrm{grad}}[A]+f_8^{(0)}[A]+f_2^{(\le 4)}[A]+f_4^{(\le 2)}[A]+\mathcal{O}(B^m{A}^n{)}_{m+n\ge 7},$$

$$\displaystyle f_{0}=-T\int\frac{\mathrm{d}^{3}\boldsymbol{p}}{(2\pi)^{3}}\ln\Bigl{(}\bigl{(}1+e^{-\xi_{\boldsymbol{p}}^{-}/T}\bigr{)}\bigl{(}1+e^{-\xi_{\boldsymbol{p}}^{+}/T}\bigr{)}\Bigr{)},$$

$$\displaystyle f_{\rm grad}[{A}]=K^{(0)}\Bigl{(}\nabla_{i}{A}^{ba\ast}\nabla_{i}{A}^{ab}+\nabla_{i}{A}^{ia\ast}\nabla_{j}{A}^{aj}+\nabla_{i}{A}^{ja\ast}\nabla_{j}{A}^{ai}\Bigr{)},$$

$$f_8^{(0)}[A]$$

$$f_2^{(\le 4)}[A]$$

$$f_4^{(\le 2)}[A]$$

$$K^{(0)}$$

$${\alpha }^{(0)}$$

$${\beta }^{(0)}$$

$${\gamma }^{(0)}$$

$${\delta }^{(0)}$$

$${\beta }^{(2)}$$

$${\beta }^{(4)}$$

$${\gamma }^{(2)}$$

$$\displaystyle A(d;\boldsymbol{n})=\left(\begin{array}[]{ccc}-F_{1}(d;\boldsymbol{n})&G_{3}(d;\boldsymbol{n})&G_{2}(d;\boldsymbol{n})\\ G_{3}(d;\boldsymbol{n})&-F_{2}(d;\boldsymbol{n})&G_{1}(d;\boldsymbol{n})\\ G_{2}(d;\boldsymbol{n})&G_{1}(d;\boldsymbol{n})&F_{1}(d;\boldsymbol{n})+F_{2}(d;\boldsymbol{n})\end{array}\right).$$

$$f_{\mathrm{grad}}[A]$$

$$-{\nabla }_d\frac{\delta f[A]}{\delta ({\nabla }_d{F}_{\alpha })}+\frac{\delta f[A]}{\delta F_{\alpha }}$$

$$-{\nabla }_d\frac{\delta f[A]}{\delta ({\nabla }_d{G}_{\beta })}+\frac{\delta f[A]}{\delta G_{\beta }}$$

$${\mathbf{n}}^{t}A(0)\mathbf{n}\equiv \sum _{i,j=1,2,3}{n}_i{A}_{ij}(0)n_j=0$$

$$\tilde{A}\equiv \frac{p_F}{T_{c0}}A,f_{\alpha }\equiv \frac{p_F}{T_{c0}}{F}_{\alpha },g_{\beta }\equiv \frac{p_F}{T_{c0}}{G}_{\beta },\mathbf{b}\equiv \frac{{\gamma }_n}{(1+F_0^a)T_{c0}}\mathbf{B},x\equiv \frac{m{T}_{c0}}{p_F}d,$$

$${\tilde{K}}^{(0)}$$

$${\tilde{\alpha }}^{(0)}$$

$${\tilde{\beta }}^{(0)}$$

$${\tilde{\gamma }}^{(0)}$$

$${\tilde{\delta }}^{(0)}$$

$${\tilde{\beta }}^{(2)}$$

$${\tilde{\beta }}^{(4)}$$

$${\tilde{\gamma }}^{(2)}$$

$$f[A]=N(0)T_{c0}^2\tilde{f}[\tilde{A}].$$

$$\displaystyle\tilde{A}^{\mathrm{bulk}}=\left(\begin{array}[]{ccc}-f_{1}^{\mathrm{bulk}}&0&0\\ 0&-f_{2}^{\mathrm{bulk}}&0\\ 0&0&f_{1}^{\mathrm{bulk}}+f_{2}^{\mathrm{bulk}}\end{array}\right),$$

$$\displaystyle\tilde{A}^{\mathrm{bulk}}=\tilde{A}_{0}\left(\begin{array}[]{ccc}r&0&0\\ 0&-1-r&0\\ 0&0&1\end{array}\right),$$

$$\tilde{A}(x)$$

$$\sigma (\mathbf{n})={\int }_0^{\infty }\Delta f(d;\mathbf{n})\mathrm{d}d=\frac{p_F^2{T}_{c0}}{2{\pi }^2}\tilde{\sigma }(\mathbf{n}),$$

$$\tilde{\sigma }(\mathbf{n})={\int }_0^{\infty }\Delta \tilde{f}(x;\mathbf{n})\mathrm{d}x,$$

$${\tilde{\sigma }}_{\mathrm{av}}\equiv \frac{1}{4\pi }{\int }_0^{\pi }\mathrm{d}\theta \sin \theta {\int }_0^{2\pi }\mathrm{d}\phi \tilde{\sigma }(\mathbf{n}).$$

$$M_{\mathrm{red}}\simeq \frac{S^1\times S^1}{{\mathbb{Z}}_2},$$

$${\tilde{\mathbf{A}}}_n(\theta ,\phi )={\tilde{A}}_{\theta }(\theta ,\phi ){\mathbf{n}}_{\theta }+{\tilde{A}}_{\phi }(\theta ,\phi ){\mathbf{n}}_{\phi },$$

$$\sum _i{\mathrm{index}}_{x_i}(\mathbf{v})=\chi (\mathcal{M}),$$

$$-{\nabla }_d\frac{\delta f}{\delta ({\nabla }_d{F}_1)}+\frac{\delta f}{\delta F_1}$$

Full text

I Introduction

Neutron stars provide us extreme environments for nuclear physics, such as high density state, rapid rotation, strong magnetic field, strong gravitational field, and so on, leading to the unveiling of new phases of nuclear systems (see Refs. Graber et al. (2017); Baym et al. (2018) for recent reviews). A variety of phases are considered to be inside neutron stars: there can be neutron rich gas and crusts at the surface, and there can be neutron superfluidity, hyperon matter, $\pi$ and/or $K$ condensates, and quark matter in the inner core. The most recent observational developments include the observations of massive neutron stars whose masses are almost twice as large as the solar mass Demorest et al. (2010); Antoniadis et al. (2013) and the observation of gravitational waves from a binary neutron star merger Abbott et al. (2017). One of the bulk phases inside neutron stars is neutron superfluidity, in which neutron pairs make a condensate in the ground state. The neutron superfluidity is directly related to astrophysical phenomena (see Refs. Chamel (2017); Haskell and Sedrakian (2018); Sedrakian and Clark (2019) for recent reviews). For example, it has been discussed that the neutron superfluidity affects relaxation time after pulsar glitches (i.e., sudden speed-up events of neutron-star rotation) Reichley and Downs (1969),111We notice that pulsar glitch phenomena can be explained by unpinning of a large amount of superfluid vortices Anderson and Itoh (1975). and the neutron superfluidity enhances rapid cooling by neutrino emissions from the inside of neutron stars (the modified Urca process) Yakovlev et al. (1999).

From the viewpoint of nuclear physics, it is an interesting property that nuclear forces have different attractive channels depending on the baryon number density, from low to high density (see Ref. Dean and Hjorth-Jensen (2003) for a recent review). In the early stage, Migdal considered the ${}^{1}S_{0}$ channel to be the most attractive one at low density Migdal (1960), though this channel turns out to be repulsive due to the strong core repulsion at higher densities Wolf (1966). At high density, instead, the attraction is supplied by the ${}^{3}P_{2}$ channel stemming from the $LS$ potential, leading to neutron ${}^{3}P_{2}$ superfluidity Tabakin (1968); Hoffberg et al. (1970); Tamagaki (1970); Takatsuka and Tamagaki (1971); Takatsuka (1972); Takatsuka and Tamagaki (1993); Amundsen and Ostgaard (1985); Baldo et al. (1992); Elgaroy 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).222We notice that the ${}^{3}P_{0}$ and ${}^{3}P_{1}$ channels are repulsive from the $LS$ potential, and hence that those channels do not contribute to the neutron pairing. More precisely, the neutron ${}^{3}P_{2}$ superfluidity exists with a small fraction of superconducting protons and normal electrons. However, such a mixture effect can be neglected in most cases.

Neutron ${}^{3}P_{2}$ superfluidity has important impacts on the observations of neutron stars. One example is the rapid neutrino cooling in neutron stars, as studied for Cassiopeia A Heinke and Ho (2010); Shternin et al. (2011); Page et al. (2011), although the existence of ${}^{3}P_{2}$ superfluidity is still elusive. Another example is the tolerance of the spin-triplet pairing in ${}^{3}P_{2}$ superfluidity against strong magnetic fields, while the spin-singlet pairing in the ${}^{1}S_{0}$ superfluidity is easily broken due to the Zeeman effect. The tolerance may give us a chance to observe ${}^{3}P_{2}$ superfluidity in 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.333In the literature, the origin of the strong magnetic fields has been studied in terms of several mechanisms such as spin-dependent interactions between two neutrons Brownell and Callaway (1969); Rice (1969); Silverstein (1969); Haensel and Bonazzola (1996), pion domain walls Eto et al. (2013); Hashimoto (2015), spin polarizations in the quark-matter core Tatsumi (2000); Nakano et al. (2003); Ohnishi et al. (2007) and so on. However, this problem remains still hard to deal with. It may be significant that the recent many-body calculation indicates a negative result for the generation of strong magnetic fields Bordbar and Bigdeli (2008).

In the pairing in neutron ${}^{3}P_{2}$ superfluidity, there are a variety of pairing combinations of relative angular momentum and total spin, leading to the existence of different phases which have symmetry breaking from the $\mathrm{U}(1)_{B}\times\mathrm{SO}(3)_{S}\times\mathrm{SO}(3)_{L}\times T\times P$ symmetry ($B$ for the baryon number, $S$ for spin rotation, $L$ for spatial rotation, $T$ for the time-reversal symmetry, and $P$ for the parity symmetry) Fujita and Tsuneto (1972); Richardson (1972); Sauls and Serene (1978); Muzikar et al. (1980); Sauls et al. (1982); Vulovic and Sauls (1984); Masuda and Nitta (2016, 2020). Examples of the symmetry breaking pattern will be presented shortly. As a consequence of the symmetry breaking, there appear low-energy excitations which affect the cooling process by neutrino emission Bedaque et al. (2003); Bedaque and Nicholson (2013); Bedaque and Reddy (2014); Bedaque and Sen (2014); Bedaque et al. (2015); Leinson (2010a, b, c, 2011a, 2011b, 2011c, 2012, 2013, 2015).444It should be noted that the cooling process is related not only to low-energy excitations but also to quantum vortices Shahabasyan and Shahabasyan (2011). Recently, studies of neutron ${}^{3}P_{2}$ superfluidity have been also devoted to understanding its topological properties: topological superfluidity and gapless Majorana fermions on the boundary of ${}^{3}P_{2}$ superfluids Mizushima et al. (2017), a quantized vortex Muzikar et al. (1980); Sauls et al. (1982); Masuda and Nitta (2016); Masaki et al. (2020), a soliton on it Chatterjee et al. (2017), and a half-quantized non-Abelian vortex Masuda and Nitta (2020). Those states have relevance to analogous states in condensed matter physics, such as $D$-wave superconductivity Mermin (1974), $P$-wave superfluidity in 3He liquid Vollhardt and Wölfle (2013), chiral $P$-wave superconductivity, e.g., in Sr2RuO4 Mackenzie and Maeno (2003), spin-2 Bose-Einstein condensates Kawaguchi and Ueda (2012), and so on.

Neutron ${}^{3}P_{2}$ superfluidity can be described by the Bogoliubov–de-Gennes (BdG) equation as the fundamental equation Tabakin (1968); Hoffberg et al. (1970); Tamagaki (1970); Takatsuka and Tamagaki (1971, 1993); Amundsen and Ostgaard (1985); Baldo et al. (1992); Elgaroy 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); Mizushima et al. (2020). In fact, the BdG equation was successfully applied to study phase structures and topological properties Mizushima et al. (2017). As a special case, the BdG equation can be reduced to the Ginzburg-Landau (GL) equation as the low-energy bosonic effective theory near the critical temperature Fujita and Tsuneto (1972); Richardson (1972); Sauls and Serene (1978); Muzikar et al. (1980); Sauls et al. (1982); Vulovic and Sauls (1984); Masuda and Nitta (2016, 2020); Yasui et al. (2019a, b); Mizushima et al. (2020). At the weak-coupling limit in the GL equation, the ground state is in the nematic phase, i.e., uniaxial nematic (UN) phase or biaxial nematic (BN) phase Sauls and Serene (1978). Each phase has different patterns of the symmetry breaking. The UN phase has an unbroken O(2) symmetry. The BN phase is separated furthermore into D2-BN and D4-BN phases, which have unbroken dihedral symmetries, D2 and D4, respectively. The phase realized in the ground state depends on the temperature and the magnetic field. The UN phase is favored at zero magnetic field, and the BN phases are favored at finite magnetic fields Masuda and Nitta (2016); Mizushima et al. (2017), More precisely, the D2-BN phase is favored at weak magnetic field, while the D4-BN phase is favored at strong magnetic field. Thus, the D2-BN and D4-BN phases would be relevant for magnetars. The GL equation has been also applied to study vortices in neutron ${}^{3}P_{2}$ superfluids: vortex structures of neutron ${}^{3}P_{2}$ superfluidity Richardson (1972); Muzikar et al. (1980); Sauls et al. (1982) and spontaneous magnetization in the core of the vortices Sauls et al. (1982); Masuda and Nitta (2016, 2020); Chatterjee et al. (2017).

It should be noted that the GL equation is an expansion series for small magnitude of the order parameter near the critical temperature. For terms up to the fourth order for the order parameter, the UN, D2-BN, and D4-BN phases are degenerate, and hence the ground state cannot be determined uniquely.555At the fourth order, an $\mathrm{SO}(5)$ symmetry happens to exist as an extended symmetry in the potential term, which is absent in the original Hamiltonian. In this case, the spontaneous breaking eventually leads to a quasi-Nambu-Goldstone mode Uchino et al. (2010).

Such degeneracy is resolved in the sixth order, where one of the BN, D2-BN, and D4-BN phases is realized in the ground state according to the temperature and the magnetic field. However, the sixth-order term cannot fully support the stability of the ground state. The next-to-leading order terms of the magnetic field were investigated for strong magnetic fields Yasui et al. (2019a). Recently, it was shown that the stability of the ground state is guaranteed at the eighth order of the condensate Yasui et al. (2019b). Interestingly, there exists a first-order phase transition at low temperature and weak magnetic field. The existence of first-order phase transitions is consistent with the result from the BdG equation Mizushima et al. (2017).

So far many studies of neutron ${}^{3}P_{2}$ superfluidity have been devoted to the bulk properties (except for Majorana fermions on the surface Mizushima et al. (2017)). In neutron stars, however, the neutron ${}^{3}P_{2}$ superfluid interfaces with the other phases, such as a neutron ${}^{1}S_{0}$ phase and nuclear crusts near the surface. Thus, it is an important question to ask how properties of the neutron ${}^{3}P_{2}$ superfluid change at the boundary. As an analogous situation, there have been several studies of surface effects of liquid crystals Urbanski et al. (2017) and 3He superfluids Ambegaokar et al. (1974); Buchholtz and Rainer (1979); Buchholtz and Zwicknagl (1981); Buchholtz (1986); Nagai and Hara (1988) (see also Refs. Vollhardt and Wölfle (2013); Volovik (2006) and the references therein), and applications of the geometrical structure of 3He liquid were developed: droplets Barton and Moore (1975), slabs Fetter (1976), pore geometry Takagi (1987), and so on. In neutron stars, the neutron ${}^{3}P_{2}$ superfluid will necessarily interface with neutron ${}^{1}S_{0}$ superfluidity at low density. As the most simple situation, however, we suppose that neutron ${}^{3}P_{2}$ superfluid interfaces with other phases at a sharp boundary at the surface of the neutron star. We notice that the phase with which the neutron ${}^{3}P_{2}$ superfluid interfaces is not necessarily the vacuum, but is possibly another phase such as the normal phase of neutron gas, the neutron ${}^{1}S_{0}$ superfluid phase, and so on. In spite of the simple situation, we will elucidate that neutron ${}^{3}P_{2}$ superfluidity exhibits nontrivial properties of symmetry breaking and topology at the surface. We first find that the phase structure on the surface can be different from that in the bulk, and therefore symmetry restoration or breaking can occur in the vicinity of or on the surface. We also find that the distribution of the surface energy density has an anisotropy depending on the polar angle in the sphere. Thus, this may lead to the deformation of the geometrical shape of the surface. Also, the order parameter manifold (OPM) induced on the surface, which is described by two-dimensional vector fields induced on the surface from the condensates, allows topological defects (vortices) on the surface. We show that there must exist such defects even in the ground state due to the Poincaré-Hopf theorem: The numbers of the vortices and antivortices depend on the bulk phases, but the difference between them is topologically invariant (the Euler number $\chi=2$) irrespective to the bulk phases. We point out that these vortices are not extended to the bulk since the first homotopy group is trivial for the OPM in the bulk while it is nontrivial for the OPM reduced on the boundary. Such defects are called boojums, which were named after Lewis Carroll’s poem about imaginary monsters, in the context of liquid crystals Urbanski et al. (2017) and helium-3 superfluids Vollhardt and Wölfle (2013); Volovik (2006). The surface properties of the neutron ${}^{3}P_{2}$ superfluids that we find in this paper will hopefully provide us useful information about neutron stars.

This paper is organized as follows. In Sec. II, we summarize the GL equation up to the eighth order. On the surface, we introduce the boundary condition for condensates of neutron ${}^{3}P_{2}$ superfluids. In Sec. III, we perform numerical analyses by solving the GL equations with a spherical boundary condition. After summarizing the symmetry breaking in the bulk space, we show results of new patterns of the symmetry breaking near the surface. Furthermore, we show anisotropic distributions of the surface energy density on the neutron star, and present the emergence of topological defects (vortices or boojums) on the surface. The final section is devoted to our conclusion and perspectives. In Appendix A, we show explicit forms of the Euler-Lagrange equation from the GL equation. In Appendix B, we briefly summarize the symmetries of neutron ${}^{3}P_{2}$ superfluidity. In Appendix C, we show the numerical results of the profile functions in the condensate.

II Formalism

II.1 Ginzburg-Landau equation

The condensate of the neutron ${}^{3}P_{2}$ superfluidity can be expressed by a symmetric and traceless three-by-three tensor $A$ as an order parameter of the symmetry breaking. The components of $A$ are denoted by $A^{ab}$ with the indices $a=1,2,3$ for the spin indices and $b=1,2,3$ for the space indices. The Ginzburg-Landau (GL) equation can be obtained by integrating out the neutron degrees of freedom and by adopting the loop expansion for the small coupling strength in the ${}^{3}P_{2}$ interaction for two neutrons Fujita and Tsuneto (1972); Richardson (1972); Sauls and Serene (1978); Muzikar et al. (1980); Sauls et al. (1982); Vulovic and Sauls (1984); Masuda and Nitta (2016, 2020); Yasui et al. (2019a, b). The GL equation is valid in the region in which the temperature $T$ is close to the critical temperature $T_{c0}$; $|1-T/T_{c0}|\ll 1$, where $T_{c0}$ is determined at zero magnetic field. The concrete form of the GL free energy reads

[TABLE]

as an expansion in terms of the condensate $A$ and the magnetic field $\boldsymbol{B}$. Each term is explained as follows. The first term $f_{0}$ is the sum of the free part and the spin-magnetic coupling term

[TABLE]

with $\displaystyle\xi_{\boldsymbol{p}}^{\pm}=\xi_{\boldsymbol{p}}\pm|\boldsymbol{\mu}_{n}||\boldsymbol{B}|$ and $\displaystyle\xi_{\boldsymbol{p}}=\boldsymbol{p}^{2}/(2m)-\mu$ for the neutron three-dimensional momentum $\boldsymbol{p}$, the neutron mass $m$, and the neutron chemical potential $\mu$. The bare magnetic moment of a neutron is $\boldsymbol{\mu}_{n}=-(\gamma_{n}/2)\boldsymbol{\sigma}$ with the gyromagnetic ratio $\gamma_{n}=1.2\times 10^{-13}$ MeV/T (in natural units, $\hbar=c=1$) and the Pauli matrices for the neutron spin $\boldsymbol{\sigma}$. The following terms include the condensate $A$: $f_{8}^{(0)}[{A}]$ consists of the terms including the field $A$ up to the eighth order with no magnetic field, $f_{2}^{(\leq 4)}[{A}]$ consists of the terms including the field $A$ up to the second order with the magnetic field up to $|\boldsymbol{B}|^{4}$, and $f_{4}^{(\leq 2)}[{A}]$ consists of the terms including the field $A$ up to the fourth order with the magnetic field up to $|\boldsymbol{B}|^{2}$. We show their explicit forms,

[TABLE]

for the gradient term (summed over the repeated indices) and

[TABLE]

for the potential and interaction terms.666In the present paper, we separate the derivative term and the nonderivative term for later convenience. This convention is different from the previous works by the present authors. The trace ($\mathrm{tr}$) is taken over the spin and space indices of $A$. The coefficients are given by

[TABLE]

We denote $N(0)=m\,p_{F}/(2\pi^{2})$ for the state-number density at the Fermi surface and $|\boldsymbol{\mu}_{n}^{\ast}|=(\gamma_{n}/2)/(1+F_{0}^{a})$ for the magnitude of the magnetic momentum of a neutron modified by the Landau parameter $F_{0}^{a}$. Notice that $\boldsymbol{\mu}_{n}^{\ast}$ is different from the bare magnetic moment of the neutron $\boldsymbol{\mu}_{n}$. It may be worthwhile to remember that the interaction between the neutron and the magnetic field ($\boldsymbol{B}$) supplies the energy splitting by the interaction Hamiltonian $-\boldsymbol{\mu}_{n}^{\ast}\cdot\boldsymbol{B}$. We notice that the Landau parameter stems from the Hartree-Fock approximation, which is not taken into account explicitly in the present mean-field approximation. In the expressions in Eq. (7), $\zeta(n)$ is the zeta function. In the present study, we will focus on the energy difference between the state with the superfluid state ($A\neq 0$) and the normal state ($A=0$). Thus, we do not consider explicitly the contribution from $f_{0}$ in Eq. (2), and hence we will neglect $f_{0}$ in the following discussion.

II.2 Boundary condition on the surface

We consider the boundary of a neutron ${}^{3}P_{2}$ superfluid and introduce the normal vector perpendicular to the surface: $\boldsymbol{n}=(n_{1},n_{2},n_{3})$ with $n_{1}=\sin\theta\cos\varphi$, $n_{2}=\sin\theta\sin\varphi$, and $n_{3}=\cos\theta$ as shown in Fig. 1. We assume that the geometrical shape of the surface can be locally approximated by a tangent plane, where the curvature can be neglected. This simplification should be justified when the curvature radius is sufficiently larger than the coherence length of the neutron ${}^{3}P_{2}$ superfluidity. We introduce the $d$ axis from the surface toward the center of the neutron star (Fig. 1). Then, the condensate $A$ is a $3\times 3$ matrix whose components are functions of $d(\geq 0)$:

[TABLE]

For precision, we should include $\boldsymbol{n}$ for the variables, because $A$ should depend on the choice of the surface with the normal vector $\boldsymbol{n}$. In the following equations, however, we will sometimes omit $\boldsymbol{n}$ for shorter notation: $A(d)=A(d;\boldsymbol{n})$, $F_{\alpha}(d)=F_{\alpha}(d;\boldsymbol{n})$, and $G_{\beta}(d)=G_{\beta}(d;\boldsymbol{n})$ with $\alpha=1,2$ and $\beta=1,2,3$. With the above setup, we rewrite the gradient terms in (3) as

[TABLE]

with $\nabla_{\!d}=\partial/\partial d$. We emphasize that $\boldsymbol{n}$ is assumed to be a constant vector and hence the derivatives with respect to $\theta$ and $\varphi$ are not included in Eq. (12). With the above-mentioned coordinate setting, we obtain the Euler-Lagrange (EL) equations for $A$,

[TABLE]

The concrete expressions of the left-hand sides are presented in detail in Appendix A. We impose the boundary conditions at $d=0$ and $d\rightarrow\infty$ for the EL equations. On the surface (at $d=0$), we adopt the condition that $A(0)$ satisfies

[TABLE]

as the boundary condition. This is consistent with the condition presented by Ambegaokar, de Gennes, and Rainer Ambegaokar et al. (1974).777In Ref. Ambegaokar et al. (1974), the authors considered the surface with $\boldsymbol{n}=(0,0,1)$ and concluded that, among the components in the matrix $A_{ij}(d)$ ($i,j=1,2,3$), $A_{11}(d)$ and $A_{22}(d)$ are symmetric functions for transforming $d$ to $-d$, hence they can continue to have finite values at the surface, while $A_{33}(d)$ is an antisymmetric function for transforming $d$ to $-d$, and hence it should vanish at the surface. This property is reflected in Eq. (15) for arbitrary $\boldsymbol{n}$.

At the center of the ${}^{3}P_{2}$ condensate (at $d\rightarrow\infty$), on the other hand, we require that the state approaches the bulk state, and hence that the condensate values should satisfy $F_{1}(d)\rightarrow F_{1}^{\mathrm{bulk}}$, $F_{2}(d)\rightarrow F_{2}^{\mathrm{bulk}}$, $G_{1}(d)\rightarrow G_{1}^{\mathrm{bulk}}=0$, $G_{2}(d)\rightarrow G_{2}^{\mathrm{bulk}}=0$, and $G_{3}(d)\rightarrow G_{3}^{\mathrm{bulk}}=0$, where $F_{\alpha}^{\mathrm{bulk}}$ ($\alpha=1,2$) and $G_{\beta}^{\mathrm{bulk}}$ ($\beta=1,2,3$) are the values in the ground state in the bulk.

II.3 Dimensionless form

For the convenience of the analysis, we introduce the dimensionless quantities $\tilde{A}$, $f_{\alpha}$ ($\alpha=1,2$), $g_{\beta}$ ($\beta=1,2,3$), $\boldsymbol{b}$, and $x$ defined by

[TABLE]

instead of $A$, $F_{\alpha}$, $G_{\beta}$, $\boldsymbol{B}$, and $d$. Notice that $\tilde{A}=\tilde{A}(x)$, $f_{\alpha}=f_{\alpha}(x)$, and $g_{\beta}=g_{\beta}(x)$ are regarded as functions of the dimensionless distance $x$. With these new variables, we obtain the dimensionless form of the GL free energy, $\tilde{f}[\tilde{A}]$, which is given from $f[A]$ by replacing $A$, $F_{\alpha}$, $G_{\beta}$, $\boldsymbol{B}$, and $d$ with $\tilde{A}$, $f_{\alpha}$, $g_{\beta}$, $\boldsymbol{b}$, and $x$ as well as by replacing the dimensionful coefficients $K^{0}$, $\alpha^{(0)}$, $\beta^{(0)}$, $\gamma^{(0)}$, $\delta^{(0)}$, $\beta^{(2)}$, $\beta^{(4)}$, and $\gamma^{(2)}$ with the dimensionless ones

[TABLE]

with the normalized temperature $t=T/T_{c0}$. $f[A]$ and $\tilde{f}[\tilde{A}]$ are related by

[TABLE]

The boundary condition at $x=0$ (i.e., $d=0$) is expressed by $\boldsymbol{n}^{t}\tilde{A}(0)\boldsymbol{n}=0$.

III Numerical results

In the numerical calculation, we use the following parameter settings: the critical temperature $T_{c0}=0.2$ MeV, the neutron number density $n=0.17$ fm*-3* (the Fermi momentum $p_{F}=338$ MeV), and the Landau parameter $F_{0}^{a}=-0.75$. The value of $F_{0}^{a}$ is that of 3He liquid at low temperature. We consider without loss of generality that the magnetic field is applied along the $x_{2}$ axis: $\boldsymbol{B}=(0,B,0)$ or $\boldsymbol{b}=(0,b,0)$. This direction is chosen to minimize the total energy when the matrix $A$ ($\tilde{A}$) is expressed in the conventional form, as will be shown in Eq. (26).

III.1 Phase diagram in bulk space

Before investigating the surface effect, in this subsection, we analyze the phase diagram in the ground (uniform) state in the bulk space ($d\rightarrow\infty$), by neglecting the gradient term in the GL free energy (1). In the bulk space, we can express conventionally the order parameter $\tilde{A}$ as a constant matrix in diagonal form:

[TABLE]

where $f_{\alpha}^{\mathrm{bulk}}$ ($\alpha=1,2$) are values in the bulk space and the off-diagonal components are zero: $g_{\beta}^{\mathrm{bulk}}=0$ ($\beta=1,2,3$). Without loss of generality, as a convention, we may restrict the range of the values of $f_{1}^{\mathrm{bulk}}$ and $f_{2}^{\mathrm{bulk}}$ to satisfy $|f_{1}^{\mathrm{bulk}}|\leq|f_{2}^{\mathrm{bulk}}|\leq|f_{1}^{\mathrm{bulk}}+f_{2}^{\mathrm{bulk}}|$, i.e., $f_{2}^{\mathrm{bulk}}\leq f_{1}^{\mathrm{bulk}}$ and $f_{2}^{\mathrm{bulk}}\geq 0$. In the literature, $\tilde{A}^{\mathrm{bulk}}$ is often expressed as

[TABLE]

with the magnitude $\tilde{A}_{0}\geq 0$ and the internal parameter $-1\leq r\leq-1/2$. Both expressions are related through $f_{1}^{\mathrm{bulk}}=\tilde{A}_{0}r$ and $f_{1}^{\mathrm{bulk}}+f_{2}^{\mathrm{bulk}}=\tilde{A}_{0}$.

According to the values of $f_{1}^{\mathrm{bulk}}$ and $f_{2}^{\mathrm{bulk}}$, the order parameter $\tilde{A}^{\mathrm{bulk}}$ expresses different phases with various symmetries, as summarized in Table 1 and in Appendix B. When $f_{1}^{\mathrm{bulk}}=f_{2}^{\mathrm{bulk}}$ ($r=-1/2$), there is an O(2) symmetry around the $x_{3}$ axis, and this phase is called the uniaxial-nematic (UN) phase. When $f_{1}^{\mathrm{bulk}}<f_{2}^{\mathrm{bulk}}$ ($-1\leq r<-1/2$), the continuous symmetry is lost, and there remains only the discrete (dihedral) symmetry, whose phase is called the biaxial-nematic (BN) phase. The BN phase is furthermore classified into the case with $-1<r<-1/2$ and the case with $r=-1$. The former ($-1<r<-1/2$) has D2 symmetry and is called the D2-BN phase, while the latter ($r=-1$) has D4 symmetry and is called the D4-BN phase. The ground state is determined to minimize the GL free energy density (1) with respect to the variations of $f_{1}^{\mathrm{bulk}}$ and $f_{2}^{\mathrm{bulk}}$ (or $\tilde{A}_{0}$ and $r$)888It may be significant that, for the terms up to and including ${\cal O}(A^{4})$ at zero magnetic field in the GL free energy density (1), there happens to appear degeneracy for the UN phase, the D2-BN phase, and the D4-BN phase. In this situation, there is the SO(5) symmetry which is absent in the original Lagrangian, and the spontaneous breaking of the symmetry eventually leads to the existence of a quasi-Nambu-Goldstone mode Uchino et al. (2010)..

We show the obtained phase diagram in Fig. 2. We observe that there exist the UN phase at zero magnetic field ($b=0$), the D2-BN phase at weak magnetic field, and the D4-BN phase at strong magnetic field. Their phase boundaries are the second-order phase transition, except for the first-order phase transition at low temperature ($t\approx 0.772$) indicated by the green line in Fig. 2. Such a first-order phase transition was found recently when the eighth-order term ($\delta_{0}$ term) was included in the GL equation Yasui et al. (2019b).999The eighth-order term is important also to give stability in the ground state. It is worthwhile to mention that the existence of the first-order phase transition in the GL equation qualitatively agrees with the result from the analysis using the BdG equation Mizushima et al. (2017).

In the following subsections, we will consider the surface effects on the neutron ${}^{3}P_{2}$ superfluidity. We suppose three cases as the typical phases for the bulk space at $x\rightarrow\infty$ (i.e., $d\rightarrow\infty$). We choose the temperature $t=0.9$ and change the values of the magnetic field for each phase. We show the numerical parameter used for the magnetic field strength ($b$) and the values of $f_{1}^{\mathrm{bulk}}$, $f_{2}^{\mathrm{bulk}}$, and $r$ for each bulk phase:

(i) UN phase

$b=0$; $(f_{1}^{\mathrm{bulk}},f_{2}^{\mathrm{bulk}})=(0.64,0.64)$ and $r=-1/2$.

(ii) D2-BN phase

$b=0.15$; $(f_{1}^{\mathrm{bulk}},f_{2}^{\mathrm{bulk}})=(0.92,0.32)$ and $r=-0.74$.

(iii) D4-BN phase

$b=0.2$; $(f_{1}^{\mathrm{bulk}},f_{2}^{\mathrm{bulk}})=(1.12,0)$ and $r=-1$.

We have considered that the strength of the magnetic field can reach maximally $B\approx 10^{15}\,\mathrm{G}$ ($10^{11}\,\mathrm{T}$) at the surface of magnetars.101010Notice the unit conversion $1\,\mathrm{T}=10^{4}\,\mathrm{G}$ for the strength of a magnetic field. From Eq. (16), we obtain $b=0.24$ for $B=10^{15}$ G. We denote these three cases (i), (ii), and (iii) by circle, square, and triangle, respectively, in the phase diagram in Fig. 2.

III.2 Symmetry near the surface

For the situations (i), (iii), and (iii) in the previous subsection, we solve the EL equations (13) and (14) with the boundary conditions, i.e., $\boldsymbol{n}^{t}\tilde{A}(0)\boldsymbol{n}=0$ at $x=0$ from Eq. (15) and $\tilde{A}(x)\rightarrow\tilde{A}^{\mathrm{bulk}}$ in the bulk space ($x\rightarrow\infty$). We plot the obtained profile functions $f_{1}(x)$, $f_{2}(x)$, $g_{1}(x)$, $g_{2}(x)$, and $g_{3}(x)$ for several choices of the surface direction $\boldsymbol{n}=(n_{1},n_{2},n_{3})$ ($\theta$ and $\varphi$) in Figs. 6, 7, and 8 in Appendix C. We observe that the profile functions approach constant values in the bulk space as the boundary condition at $x\rightarrow\infty$, while they change drastically in the region x\mathrel{\hbox to0.0pt{\raise 2.20013pt\hbox{<}\hss}{\lower 2.20013pt\hbox{\sim}}}1 near the surface. We notice that the typical value of the healing distance from the surface is expressed by $r\approx\xi$ (i.e., $x\approx 1$) with $\xi=p_{F}/(mT_{c0})\simeq 360$ fm for $p_{F}=338$ MeV and $T_{c0}=0.2$ MeV. The changes of the profile functions induce the symmetry breaking near the surface. This can be seen directly by transforming $\tilde{A}(x)$ into diagonal form with some appropriate unitary matrix $U(x)$ at each position $x$:

[TABLE]

with $\tilde{A}_{0}(x)\geq 0$ and $-1\leq r(x)\leq-1/2$. Notice that $\tilde{A}_{0}$ and $r$ are functions of the position $x$, as they have been constant values in the bulk space [cf. Eq. (26)]. We plot $r(x)$ in Figs. 9, 10, and 11 in Appendix C for the bulk conditions (i), (ii), and (iii). We confirm that $r(x)$ is dependent on $x$, and that there exist various phases, i.e., the UN phase [$r(x)=-1/2$], the D2-BN phase [$-1<r(x)<-1/2$], and the D4-BN phase [$r(x)=-1$], that are different from those in the bulk.

The symmetry restoration or breaking near the surface can be understood in a reasonable way. First of all, we notice that the boundary condition at the surface, $\boldsymbol{n}^{t}\tilde{A}(0)\boldsymbol{n}=0$ in Eq. (15), does not allow the general transformation of spin and space for $\tilde{A}(0)$, because $\boldsymbol{n}$ is a fixed vector in space. Thus, the symmetry cannot be generally maintained from the bulk to the surface, and it should exhibit variation near the surface. However, there are some exceptional cases in which the symmetry is maintained from the bulk to the surface: $\boldsymbol{n}=(0,0,1)$ ($\theta=0$ and $\varphi$=0) with the bulk UN phase and $\boldsymbol{n}=(0,\pm 1,0)$ ($\theta=\pi/2$ and $\varphi=\pi/2$, $3\pi/2$) with the bulk D4-BN phase. In both cases, the profile functions are constant at any $x$, and hence the symmetry is kept invariant, as shown in Figs. 6 and 8 as well as in Figs. 9 and 11. In order to understand this invariance, we remember that the UN phase possesses the O(2) symmetry around the $x_{3}$ axis, and the D4-BN phase has the D4 symmetry around the $x_{2}$ axis. In the former case, we find that the boundary condition $\boldsymbol{n}^{t}\tilde{A}(0)\boldsymbol{n}=0$ with $\boldsymbol{n}=(0,0,1)$ holds for any transformation of $\tilde{A}(0)$ under the U(1) symmetry. In the latter case, we also find that the boundary condition $\boldsymbol{n}^{t}\tilde{A}(0)\boldsymbol{n}=0$ with $\boldsymbol{n}=(0,1,0)$ holds for any transformation of $\tilde{A}(0)$ under the D4 symmetry.

From Figs. 6, 7, and 8 and Figs. 9, 10, and 11, we notice that the D2-BN phase [$-1<r(0)<-1/2$] is realized at most points on the surface. However, there are special points, $(\theta,\varphi)=(0,0)$, $(\pi/2,0)$, $(\pi/2,\pi/2)$, $(\pi/2,\pi)$, $(\pi/2,3\pi/2)$, and $(\pi,0)$, where either of the UN phase [$r(0)=-1/2$] or the D4-BN phase [$r(0)=-1$] is realized. We call those points the UN points or the D4-BN points, and they appear differently for each bulk phase:

(i) UN phase

the UN points at $(\theta,\varphi)=(0,0)$, $(\pi,0)$ and the D4-BN points at $(\theta,\varphi)=(\pi/2,0)$, $(\pi/2,\pi/2)$, $(\pi/2,\pi)$, $(\pi/2,3\pi/2)$,

(ii) D2-BN phase

the D4-BN points at $(\theta,\varphi)=(0,0)$, $(\pi/2,0)$, $(\pi/2,\pi/2)$, $(\pi/2,\pi)$, $(\pi/2,3\pi/2)$, $(\pi,0)$,

(iii) D4-BN phase

the D4-BN points at $(\theta,\varphi)=(0,0)$, $(\pi/2,0)$, $(\pi/2,\pi/2)$, $(\pi/2,\pi)$, $(\pi/2,3\pi/2)$, $(\pi,0)$,

as displayed graphically in Fig. 3. We notice that the UN points appear only for the bulk UN phase, and not for the bulk D2-BN phase or the bulk D4-BN phase.

III.3 Energy density and topological defects on the surface

We consider the surface energy density, which is defined by the energy density per unit area on the surface. The surface energy stems from the gradient term and the difference between the potential value near the surface and the potential value in the bulk space. With the solutions of $f_{\alpha}(x)$ ($\alpha=1,2$) and $g_{\beta}(x)$ ($\beta=1,2,3$) from the EL equations (13) and (14), we obtain the surface energy density expressed as

[TABLE]

for $\Delta f(d;\boldsymbol{n})=f(d;\boldsymbol{n})-f_{\mathrm{bulk}}$ with the GL free energy density $f_{\mathrm{bulk}}$ in the bulk space ($d\rightarrow\infty$). In the preceding equation, for convenience, we have defined the dimensionless surface energy density by

[TABLE]

for $\Delta\tilde{f}(x;\boldsymbol{n})=\tilde{f}(x;\boldsymbol{n})-\tilde{f}_{\mathrm{bulk}}$ with the dimensionless GL free energy density $\tilde{f}_{\mathrm{bulk}}$ in the bulk space ($x\rightarrow\infty$).111111Here we recover $\boldsymbol{n}$ to emphasize that the surface energy density depends on the direction of the normal vector. In Eq. (31), the value of the coefficient is given by ${p_{F}^{2}T_{c0}}/({2\pi^{2}})\simeq 27$ $\mathrm{keV}/\mathrm{fm}^{2}$ for $p_{F}=338$ MeV and $T_{c0}=0.2$ MeV. We show the numerical results of the distributions of the surface energy density $\tilde{\sigma}(\boldsymbol{n})$ on the plane spanned by $\theta$ and $\varphi$ in the top row in Fig. 4. For the bulk D4-BN phase, for example, the minimum surface energy density [$\sigma(\boldsymbol{n})=0$] is realized at the point $\boldsymbol{n}=(0,\pm 1,0)$ ($\theta=\pi/2$ and $\varphi=\pi/2$, $3\pi/2$). Notice that zero density is reasonable because there is no change in $f_{\alpha}(x)$ ($\alpha=1,2$) and $g_{\beta}(x)$ ($\beta=1,2,3$), as discussed in the previous section. In this situation, because the magnetic field is applied along the $x_{2}$ axis, it is expected that the shape of the surface of the neutron star should be deformed to be an oblate spheroid with the short axis being aligned along the $x_{2}$ axis when the surface shape can be changed to be in balance with the Fermi pressure and the gravity. In Fig. 4, we observe that the geometrical distribution of the surface energy density on the sphere obeys the dihedral symmetries as a subgroup of the spherical symmetry: the D4 symmetry for the bulk UN phase, the D2 symmetry for the bulk D2-BN phase, and the D4 symmetry for the D4-BN phase, as summarized in Table 2.

We estimate the averaged values of the surface energy density. For this purpose, we define the averaged values of $\sigma(\boldsymbol{n})$ over $\theta$ and $\varphi$ by $\sigma_{\mathrm{av}}\equiv({p_{F}^{2}T_{c0}}/({2\pi^{2}}))\tilde{\sigma}_{\mathrm{av}}$ and

[TABLE]

From the values shown in Table 2, they are almost close to $\sigma_{\mathrm{av}}\approx 0.2$ keV/fm2 for each bulk phase. Among them, the minimum surface energy is provided by the bulk D4-BN phase. With the values of $\sigma_{\mathrm{av}}$, we obtain the surface energy in total, $E_{\mathrm{surf}}=4\pi R^{2}\sigma_{\mathrm{av}}$, for the neutron star with the radius $R$. Numerically, we obtain $E_{\mathrm{surf}}\approx 6\times 10^{29}$ erg for $R=10$ km. The values of $E_{\mathrm{surf}}$ for each bulk phase are summarized in Table 2.

The directions of the condensate $A$ at the surface have unique topological properties. In order to see this, we define the vector $\boldsymbol{A}_{n}(\theta,\varphi)\equiv A(0)\boldsymbol{n}$ with the normal vector $\boldsymbol{n}$ at the surface [$A(0)$ is a $3\times 3$ matrix].

There are two remarks: First, we remember that the $d$ vector is defined by $A\hat{\boldsymbol{p}}$ with $\hat{\boldsymbol{p}}=\boldsymbol{p}/|\boldsymbol{p}|$ for the three-dimensional momentum $\boldsymbol{p}$. Because the axis direction perpendicular to the surface is considered to be a one-dimensional system, we reasonably regard $\hat{\boldsymbol{p}}=\boldsymbol{n}$, and hence conclude that $\boldsymbol{A}_{n}(\theta,\varphi)$ [$\tilde{\boldsymbol{A}}_{n}(\theta,\varphi)$] is the same as the $d$ vector at the surface.

Second, we also comment that, more precisely, the vector field $\boldsymbol{A}_{n}$ parametrizes an order parameter space reduced at the boundary. Namely, at the boundary, the order parameter space is reduced due to the boundary condition to its submanifold,

[TABLE]

where two $S^{1}$’s denote the U(1) phase and the spatial rotation around $\boldsymbol{n}$, and $\mathbb{Z}_{2}=\{1,-1\}$ is introduced to remove identical two points; a simultaneous transformation of a $\pi$ phase rotation and a $\pi$ spatial rotation around $\boldsymbol{n}$ does not change the order parameters and should be removed. Since we restrict $A$ to be real-valued, $\boldsymbol{A}_{n}(\theta,\varphi)$ is a real vector parametrizing one $S^{1}$, and $\boldsymbol{A}_{n}(\theta,\varphi)$ is a two-dimensional vector on the plane orthogonal to the normal vector $\boldsymbol{n}$ because $\boldsymbol{n}\!\cdot\!\boldsymbol{A}_{n}(\theta,\varphi)=0$ is induced from the boundary condition (15). This two-dimensional vector can be expressed by $\boldsymbol{A}_{n}(\theta,\varphi)=(T_{c0}/p_{F})\tilde{\boldsymbol{A}}_{n}(\theta,\varphi)$ with

[TABLE]

for $\boldsymbol{n}_{\theta}$ and $\boldsymbol{n}_{\varphi}$ defined by $\boldsymbol{n}_{\theta}\equiv\partial\boldsymbol{n}/\partial\theta$ and $\boldsymbol{n}_{\varphi}\equiv\partial\boldsymbol{n}/\partial\varphi$, respectively. Notice that $\boldsymbol{n}_{\theta}$ and $\boldsymbol{n}_{\varphi}$ are the unit vectors perpendicular to each other and also to $\boldsymbol{n}$, i.e., $\boldsymbol{n}\!\cdot\!\boldsymbol{n}_{\theta}=\boldsymbol{n}\!\cdot\!\boldsymbol{n}_{\varphi}=\boldsymbol{n}_{\theta}\!\cdot\!\boldsymbol{n}_{\varphi}=0$. We show a plot of $\tilde{\boldsymbol{A}}_{n}=(\tilde{A}_{\theta},\tilde{A}_{\varphi})$ on the plane spanned by $\theta$ and $\varphi$ ($0\leq\theta<\pi$ and $0\leq\varphi<2\pi$) in the bottom row in Fig. 4. Interestingly, we observe that there exist topological defects (vortices) on the plane. Each defect has a positive or negative charge $\pm 1$. The charges of the defects are defined according to the circulating directions of the vector fields as shown in Fig. 5. We denote the total positive (negative) charge by $\chi_{+}$ ($\chi_{-}$). In Table 3, we summarize the values of $\chi_{+}$ and $\chi_{-}$ for each bulk phase: $\chi_{+}=10$ and $\chi_{-}=-8$ for the bulk UN phase, $\chi_{+}=12$ and $\chi_{-}=-10$ for the bulk D2-BN phase, and $\chi_{+}=8$ and $\chi_{-}=-6$ for the bulk D4-BN phase. We define $\chi=\chi_{+}+\chi_{-}$ to denote the sum of $\chi_{+}$ and $\chi_{-}$, which is the difference between the numbers of the defects with a charge $+1$ and charge $-1$. Importantly, we find that $\chi=2$ is realized uniquely for all the different bulk phases. This is an inevitable consequence of a topological property of two-dimensional vector fields on a sphere, i.e., the Poincaré-Hopf (hairy ball, no-wind) theorem. We consider a vector field $\boldsymbol{v}$ with isolated zeros on a two-dimensional manifold ${\cal M}$. Then, the Poincaré-Hopf theorem requires the relation

[TABLE]

where $\mathrm{index}_{x_{i}}(\boldsymbol{v})$ is the index ($\pm 1$ charges) of $\boldsymbol{v}$ at the position $x_{i}$ of isolated zeros on ${\cal M}$, $\Sigma_{i}$ indicates a sum over all the isolated zeros ($i$), and $\chi({\cal M})$ is the Euler characteristic of ${\cal M}$. $\chi({\cal M})=2$ if ${\cal M}$ is a sphere $S^{2}$ relevant for us. For example, for the bulk UN phase, we have $\sum_{i}\mathrm{index}_{x_{i}}(\boldsymbol{A}_{n})=8+(-6)=2$, which is indeed identical to $\chi=2$. This relation is checked for the bulk D2-BN phase and for the bulk D4-BN phase. Therefore, we understand that $\chi=2$ should hold for any bulk phases. On the other hand, $\chi_{+}$ and $\chi_{-}$ can vary according to the symmetries in the bulk phases. In order for $\chi=2$ to hold, creations or annihilations should happen for a pair of vortices with charges $\pm 1$. From Table 3, we see that two pairs of vortices are created from the bulk UN phase to the bulk D2-BN phase, and that four pairs of vortices are annihilated from the bulk D2-BN phase to the bulk D4-BN phase.

The final comment is that vortices induced on a boundary are called boojums in the context of liquid crystals Urbanski et al. (2017) and the 3He superfluids Vollhardt and Wölfle (2013); Volovik (2006). It is interesting that these boojums do not extend to the bulk since closed paths winding around $M_{\rm red}$ become trivial loops in the full order parameter space $M$, where $M_{\rm red}\subset M$. We notice that there recently appeared a study on the surface vortices for spinor Bose-Einstein condensates (BECs) with spin-orbital-angular momentum coupling Cheng et al. . The presented results indicate that there are interesting common properties, such as the symmetry on surface, the distribution of topological defects, and so on, as we have discussed in our study. Thus, it would be valuable to investigate further similarities as well as differences between the neutron ${}^{3}P_{2}$ superfluids and spinor BECs.

IV Conclusion and perspectives

We have studied surface effects of neutron ${}^{3}P_{2}$ superfluidity in the UN, D2-BN, and D4-BN phases in neutron stars. We have supposed the situation in which the neutron ${}^{3}P_{2}$ superfluid exists in a large ball with a spherical boundary, and introduced a boundary condition suitable for the condensate at the sphere’s surface. Solving the GL equation with the boundary condition, we have found several interesting properties of the surface effects of the neutron ${}^{3}P_{2}$ superfluid. First, we have shown that the symmetry in the bulk space can be restored or broken to other symmetries at the surface. Second, the distribution of the surface energy density on the sphere has an anisotropy depending on the polar angle. This will lead to the geometrical deformation of the surface of the neutron ${}^{3}P_{2}$ superfluidity from the spherical shape to an oblate spheroid. Third, we have investigated the two-dimensional vector field defined from the condensate at the surface, and have shown that there must exist topological defects (vortices) with the charges $\pm 1$ on the sphere. It should be emphasized that those defects appear in the ground state but not as an excited state. While the number of the defects (vortices) is dependent on the symmetries in the bulk phases, i.e., the UN, D2-BN, and D4-BN phases, the difference between the numbers of the defects with a charge $+1$ and charge $-1$ remains topologically invariant ($\chi=2$) thanks to the Poincaré-Hopf theorem. These vortices are called boojums in the context of liquid crystals and 3He superfluids.

In the present study, we have considered the simplest situation for the boundary condition for neutron ${}^{3}P_{2}$ superfluidity. In a more realistic situation, however, the neutron ${}^{3}P_{2}$ superfluid can interface with a neutron ${}^{1}S_{0}$ superfluid as well as with nuclear crusts composed of a lattice of neutron-rich nuclei. Deeper inside in neutron stars, a connection of the neutron ${}^{3}P_{2}$ phase to other exotic phases such as hyperon matter, quark matter, and so on should be interesting. One of the important questions to ask is how the topological objects are connected between those phases. When a neutron star rotates, Abelian quantum vortices in the hadron matter and non-Abelian quantum vortices (color magnetic flux tubes) Balachandran et al. (2006); Nakano et al. (2008); Eto and Nitta (2009) in the quark matter can be connected through colorful boojums (endpoints of vortices) Cipriani et al. (2012); Alford et al. (2019); Chatterjee et al. (2019); Cherman et al. (2019) (see Ref. Eto et al. (2014) for a review of non-Abelian quantum vortices). These boojums are different from those studied in this paper. Although we have limited ourselves to the BN, D2-BN, and D4-BN phases for the neutron ${}^{3}P_{2}$ superfluid, it will be also important to study the possibilities of cyclic and ferromagnetic phases Mermin (1974); Kawaguchi and Ueda (2012). The former phase can lead to one-third quantized non-Abelian vortices Semenoff and Zhou (2007) forming a network in collision Kobayashi et al. (2009), while the latter phase could be relevant to inner structures of magnetars. In terms of the topological matter, the cyclic and ferromagnetic phases have gapless Weyl fermions in the bulk (Weyl semimetals) Mizushima et al. (2017); Mizushima and Nitta (2018). The above mentioned subjects should be studied carefully in the future.

Acknowledgment

We would like to thank Michikazu Kobayashi and Takeshi Mizushima for discussions. This work is supported by the Ministry of Education, Culture, Sports, Science (MEXT) Supported Program for the Strategic Research Foundation at Private Universities “Topological Science” (Grant No. S1511006). This work is also supported in part by JSPS Grants-in-Aid for Scientific Research [KAKENHI Grant No. 17K05435 (S.Y.), No. 19K14713 (C.C.), No. 16H03984 (M.N.), and No. 18H01217 (M.N.)], and also by MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas “Topological Materials Science” No. 15H05855 (M.N.).

Appendix A Euler-Lagrange equations

We present the concrete expressions of the left-hand sides in the EL equations (13) and (14):

[TABLE]

for $F_{1}$,

[TABLE]

for $F_{2}$,

[TABLE]

for $G_{1}$,

[TABLE]

for $G_{2}$, and

[TABLE]

for $G_{3}$.

Appendix B Symmetries

We summarize the properties of the symmetries of the order parameter $A$ in Eq. (26). First of all, we remember that $A$ possesses the symmetry

[TABLE]

in the Lagrangian, where $e^{i\alpha}\in\mathrm{U}(1)$ and $O(\theta,\boldsymbol{n})$, $O(\theta^{\prime},\boldsymbol{n}^{\prime})\in\mathrm{SO}(3)$ with $\boldsymbol{n}$ ($\boldsymbol{n}^{\prime}$) the rotation axis and $\theta$ ($\theta^{\prime}$) the rotation angle around $\boldsymbol{n}$ ($\boldsymbol{n}^{\prime}$). Here, $O(\theta,\boldsymbol{n})$ is the rotation in the spin space, and $O(\theta^{\prime},\boldsymbol{n}^{\prime})$ is the rotation in the real space. $\tilde{\boldsymbol{x}}$ is the vector rotated by $O(\theta^{\prime},\boldsymbol{n}^{\prime})$ from $\boldsymbol{x}$. The symmetries (42) are spontaneously broken to subgroups when the state is in the nematic phase as presented in Eq. (26). Instead, there exist the O(2), D2, and D4 symmetries, and the corresponding phases are called the UN ($r=-1/2$), D2-BN ($-1<r<-1/2$), and D4-BN ($r=-1$) phases, respectively. We will give the concrete forms of those symmetries in the following.

B.1 UN phase ($r=-1/2$)

For $r=-1/2$, the order parameter $A(\boldsymbol{x})$ is written as

[TABLE]

This is invariant under the rotation around the $x_{3}$-axis $A_{\mathrm{UN}}(\boldsymbol{x})\rightarrow O(\theta)A_{\mathrm{UN}}(\boldsymbol{x}^{\prime})O^{t}(\theta)$ with $e^{i\alpha}=1$ ($\alpha=0$), where $O(\theta)\in\mathrm{O}(2)$ is the rotation operator

[TABLE]

with the rotation angle $\theta$ ($0\leq\theta<2\pi$). This is the uniaxial nematic (UN) phase. The rotation in the spin space and the rotation in the real space are locked to each other, $O(\theta,\boldsymbol{n})=O(\theta^{\prime},\boldsymbol{n}^{\prime})$ with $\boldsymbol{n}=\boldsymbol{n}^{\prime}$ and $\theta=\theta^{\prime}$. The locking also occurs for the BN phase.

B.2 D2-BN phase ($-1<r<-1/2$)

For $-1<r<-1/2$, the order parameter $A(\boldsymbol{x})$ is written as

[TABLE]

This is invariant under the D2 symmetry. The generators of the D2 group are given by

[TABLE]

with $e^{i\alpha}=1$ ($\alpha=0$), $\boldsymbol{1}_{3}$ is a unit matrix, and $I_{1}$, $I_{2}$, and $I_{3}$ are defined by

[TABLE]

$I_{i}$ ($i=1,2,3$) indicates a $\pi$ rotation around the $i$th axis. We confirm easily that $A_{\mathrm{D}_{2}\mathrm{BN}}(\boldsymbol{x})$ is invariant under the transformation $A_{\mathrm{D}_{2}\mathrm{BN}}(\boldsymbol{x})\rightarrow OA_{\mathrm{D}_{2}\mathrm{BN}}(\boldsymbol{x}^{\prime})O^{t}$ ($O\in\mathrm{D}_{2}$). This phase is called the D2-biaxial nematic (D2-BN) phase.

B.3 D4-BN phase ($r=-1$)

For $r=-1/2$, the order parameter $A$ is written as

[TABLE]

This is invariant under the D4 symmetry. The generators of the D4 group are given by

[TABLE]

where $\boldsymbol{1}_{3}$ and $I_{i}$ ($i=1,2,3$) have been defined in the D2 group, and $R_{2}$ is newly defined by

[TABLE]

$R_{2}$ indicates the $\pi/2$ rotation around the second axis ($x_{2}$ axis). It should be noted that the phase \bigl{\{}e^{i\alpha}\bigr{\}}=\bigl{\{}1,-1\bigr{\}}\in\mathbb{Z}_{2} ($\alpha=0,\pi$) is locked with the spin rotation and the spatial rotation. We confirm easily that $A_{\mathrm{D}_{4}\mathrm{BN}}$ is invariant under the transformation $A_{\mathrm{D}_{4}\mathrm{BN}}(\boldsymbol{x})\rightarrow e^{i\alpha}OA_{\mathrm{D}_{4}\mathrm{BN}}(\boldsymbol{x}^{\prime})O^{t}$ [$(e^{i\alpha},O)\in\mathrm{D}_{4}$]. This is the D4-biaxial nematic (D4-BN) phase.

Appendix C Numerical results for the profile functions

From Sec. III.2, we summarize the numerical results for $f_{\alpha}(x)$ ($\alpha=1,2,3$), $g_{\beta}(x)$ ($\beta=1,2,3$) in Figs. 6, 7, and 8 and $r(x)$ in Figs. 9, 10, and 11 for the UN, D2-BN, and D4-BN phases in the bulk.