Critical endpoint and universality class of neutron 3 P 2 superfluids in neutron stars

X iv :1 90 8. 07 94 4v 1 [ nu cl -t h] 2 1 A ug 2 01 9 Critical endpoint and universality class of neutron P2 superfluids in neutron stars Takeshi Mizushima, ∗ Shigehiro Yasui, † and Muneto Nitta ‡ Department of Materials Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan Department of Physics & Research and Education Center for Natural Sciences, Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa 223-8521, Japan (Dated: August 22, 2019)


I. INTRODUCTION
A neutron star is a compact star which is composed almost entirely of neutrons under extreme conditions such as high density, rapid rotation, and a strong magnetic field (see Refs. [1,2] for recent reviews). The most recent discoveries include the observations of massive neutron stars whose masses are almost twice as large as the solar mass [3,4] and the observation of gravitational waves from a binary neutron star merger [5]. In the inner structure, neutron superfluidity and proton superconductivity are key ingredients for understanding the evolution of neutron stars (see Refs. [6][7][8] for recent reviews). As the superfluid and superconducting components reorganize low-lying elementary excitations, their presence profoundly affects neutrino emissivities and specific heats and can explain the long relaxation time observed in the sudden speed-up events of neutron stars [9][10][11] and the enhancement of neutrino emission at the onset of superfluid transition [12][13][14][15][16][17]. Sudden changes of spin periods observed in pulsars (pulsar glitches) may also be explained by the existence of superfluid components with quantized vortices [18,19].
perfluids are expected to be realized in the inner cores of neutron stars. Furthermore, the neutron 3 P 2 superfluids have tolerance against the strong magnetic field, such as 10 15 − 10 18 G in magnetars, because the spin-triplet pairing is not broken through the spin-magnetic field interaction by Zeeman effects. 3 It has recently been proposed that the observation of the rapid cooling of the neutron star in Cassiopeia A may be explained by enhanced neutrino emissivities due to the formation and dissociation of neutron Cooper pairs in the 3 P 2 channel which is a short-ranged attraction in the total angular momentum J = 2 [16,17,52] (see also Refs. [53][54][55]). Theoretically, neutron 3 P 2 superfluids provide a fertile ground for exploring exotic superfluidity. The superfluid states with J = 2 are classified into several phases: Nematic, cyclic, and ferromagnetic phases [27,28,[56][57][58][59][60][61]. The nematic phase is further divided into the uniaxial nematic (UN) phase and the dihedral-two and dihedral-four biaxial nematic (D 2 -BN and D 4 -BN) phases. All these phases are accompanied by topologically protected Bogoliubov quasiparticles. The nematic phase is a prototype of class-DIII topological superconductors and a harbor of Majorana fermions [62]. The other phases are non-unitary states with broken time-reversal symmetry and promising platforms to host Weyl superfluidity [62,63]. In addition to such exotic fermions, the 3 P 2 order parameters also bring about rich bosonic excitations [64][65][66][67][68][69][70][71][72][73][74][75][76], which might be relevant to the cooling process by neutrino emission 4 , as well as exotic topological defects, including spontaneously magnetized vortices [27,57,58,60] and vortices with Majorana fermions [78], solitonic excitations on a vortex [79], and half- 3 The origin of the strong magnetic fields in neutrons stars or in magnetars has been studied in several types of mechanisms such as spin-dependent interactions [42][43][44][45], pion domain walls [46,47], spin polarizations in the quark-matter in the neutron star core [48][49][50] and so on. However, this problem is not settled yet. Recently, a negative result for the generation of strong magnetic fields was reported in the study based on the nuclear many-body calculations [51]. 4 The cooling process is related not only to low-energy excitations but also to quantum vortices [77]. quantized non-Abelian vortices [61], domain walls [80], and surface topological defects (boojums) on the boundary of 3 P 2 superfluids [81]. Those states share common properties in the condensed matter systems, such as D-wave superconductors [82], P-wave superfluidity in 3 He liquid [83][84][85], chiral P-wave superconductivity e.g. in Sr 2 RuO 4 [86,87] and U-based ferromagnetic superconductors [88], spin-2 Bose-Einstein condensates [89], and so on.
The neutron 3 P 2 superfluidity can be described by the Fermi liquid theory which is composed of the set of self-consistent equations based on the Luttinger-Ward thermodynamic functional. Microscopically, the most fundamental equation of the neutron 3 P 2 superfluidity is provided by the Bogoliubov-de Gennes (BdG) equation where the order parameter, i.e., gap function, should be solved self-consistently with the wavefunctions of the gapped neutrons [22-25, 30, 32-40]. The BdG equation was successfully applied to study the topological properties of the neutron 3 P 2 superfluidity [62]. The phase digram with respect to the magnetic field and temperature was obtained in Ref. [62], where the first and second-order phase transitions between the D 2 -BN and D 4 -BN phases are present and these transitions meet at a critical endpoint (CEP), as shown in Figs. 1(a) and 1(b). 5 The existence of the CEP is cerntainly important for transport coefficients and equations of state of neutron matter when neutron stars are cooled down.
Around the transition temperature from the normal phase to the superfluid phase, the Ginzburg-Landau (GL) theory can be induced by the systematic expansion of the functional with respect to the order parameter field and the magnetic field [27, 28, 56-61, 81, 90-92]. Unlike the ordinary cases, the GL expansion up to the 4th order in terms of the order parameter cannot determine the unique ground state, but there exists a continuous degeneracy among the UN, D 2 -BN and D 4 -BN phases. 6 The GL expansion up to the 6th order deter- 5 See Eqs. (30) and (34) for the definitions of G (n) . 6 We notice that, at the 4th order, there happens to exist an SO(5) symme-mines the unique ground state [60] but it is stable only locally and there exists the instability for a large value of the order parameter. Recently, in order to solve this problem, the GL equation up to the 8th order term in the condensates was obtained [92], in which it was shown that the 8th order term ensures the global stability with respect to the variation of the order parameter in the ground state. As a byproduct, it was also found that the phase diagram in the expansion up to the 8th oder possesses the CEP as shown in Fig. 1(c), in contrast to the GL equation up to the 6th order in which no CEP exists, although the positions of the CEPs in the BdG and GL formalism are rather different as shown in Fig. 1(d).
In this paper, we study the critical exponents at the CEP in the BdG equation and in the GL equation. Under the scaling hypothesis, a set of critical exponents, (α, β, γ, δ), 7 at the CEP should satisfy the the universal relations, i.e., the Rushbrooke, Griffiths, and Widom equalities In the both cases of the BdG equation and of the GL equation, we extract the critical behavior of neutron 3 P 2 superfluids by directly computing all the critical exponents at the CEP. We demonstrate that the CEP in the GL approach properly captures the critical phenomena in the BdG equation, and the exponents satisfy all three equalities reasonably in both the BdG and GL equations within a numerical error. We find that the try as an extended symmetry in the potential term, which is absent in the original Hamiltonian. In this case, the spontaneous breaking eventually generate a quasi-Nambu-Goldstone mode which should be irrelevant to the excitations in the true ground state [93]. This is nothing but the origin of the continuous degeneracy. 7 See Eqs. (57)- (59) for the definitions of (α, β, γ, δ). 3 P 2 superfluid at the CEP exhibits critical behavior with nontrivial critical exponents in such a manner that the exponents associated with the critical behaviors of the specific heat and magnetization exhibits α ∼ 0.6 and γ ∼ 0.5. In particular, the exponent γ < 1 is unique and essentially different from γ ≥ 1 in ordinary universality classes [94,95], except for a few models, e.g., O(n) models with n < 0 [96] and the tricritical Ising model coupled to massless Dirac fermions [97]. This indicates the CEP in neutron 3 P 2 superfluids belongs to a new universality class.
The organization of this paper is as follows. In Sec. II, we present the superfluid Fermi liquid theory, where the selfconsistent equations for the gap functions and Fermi liquid corrections are described in detail. This theory is based on the quasiclassical approximation which is relevant to 3 P 2 superfluids of neutrons. Based on the theory, we show that the phase diagram of 3 P 2 superfluids under strong magnetic fields has the CEP and compute the critical exponents, indicating a new universality class. Furthermore, in Sec. III, we present the GL theory up to the 8th-order expansion to examine the critical phenomena at the CEP, showing that the critical exponents in the GL theory coincide with those in the BdG theory within a certain accuracy. Sec. IV is devoted to a summary and discussion.

A. General formalism
Here we start with the Hamiltonian for neutrons interacting through the potential V c,d a,b , where r 12 ≡ r 1 − r 2 denotes the relative coordinate and ψ a and ψ † a (a =↑, ↓ for spins) denote the fermionic field operators. The single-particle energy for a neutron under a magnetic field B is given by with ξ 0 (k) = k 2 /(2m) − µ for the neutron mass m and the chemical potential µ. Here γ n = 1.2 × 10 −13 MeV/T is the gyromagnetic ratio for a neutron, 8 and σ = (σ 1 , σ 2 , σ 3 ) denotes the Pauli matrices in the spin space. In Eq. (4), V c,d a,b (r 12 ) contains microscopic informations on neutron-neutron interaction potentials. The repeated Roman and Greek indices imply the sum over the spin degrees of freedom and the threedimensional spatial component (x, y, z), respectively. In this paper, we set = k B = 1. 8 Notice the unit conversion 1 T = 10 4 G for the strength of a magnetic field. Let us define the Nambu-Gor'kov (NG) Green's function in terms of a grand ensemble average of the fermion-field operators in the Nambu space, Ψ ≡ (ψ ↑ , ψ ↓ , ψ † ↑ , ψ † ↓ ) tr , as , where x i ≡ (r i , τ i ) with the three dimensional space position r i and the imaginary time τ i for the neutron i = 1, 2. a tr denotes the transpose of the matrix a. In this paper, we consider homogeneous neutron matter and transform the space-time position x to momentum p and Matsubara frequency ε n = (2n + 1)πT (n = 0, ±1, ±2, . . . ): x → ( p, ε n ). The self-consistent formalism is derived from the Luttinger-Ward thermodynamic functional which is given in terms of the full NG Green's function G and the self-energy Σ as where with the trace (Tr) taken over the spin space and the NG (particle-hole) space. The inverse propagator for free fermions is given by , and V ext is an external field including a magnetic Zeeman term in Eq. (5). Here we use τ = (τ 1 , τ 2 , τ 3 ) to denote the matrices in the NG space. The Green's function and the self-energy are related to the functional Φ[G] by the stationary conditions with respect to the Green's function, δΩ/δG tr = 0, and the self-energy, δΩ/δΣ tr = 0. The former is recast into the definition of the self-energy in terms of the functional derivative The Dyson's equation for the full Green's function is obtained from the latter stationary condition as The above set of equations from Eq. (6) to Eq. (9) provide a starting point for deriving the quasiclassical Fermi liquid theory for 3 P 2 superfluids.

B. Quasiclassical approximation
In general, the quasiclassical approximation provides a powerful tool for describing phenomena when the characteristic lengths are much greater than the Fermi wavelength, λ F ∼ 2π/p F (p F the Fermi momentum), and characteristic frequencies are much smaller than the Fermi energy, ω ≪ ε F / (ε F the Fermi energy) [98,99]. The typical scales in the superfluid state of 3 He and superconducting states are the coherence length, ξ c ≡ v F /2πk B T c and the excitation gap ∆ 0 ∼ k B T c . The quasiclassical theory uses the fact that all relevant parameters, such as temperature T and external potentials V, are very small relative to the atomic scales which are given by Fermi temperature T F , Fermi energy ε F and Fermi momentum p F . This difference in scales allows one to perform an asymptotic expansion of full many-body propagators in small parameters T/T F and |V|/ε F , and it leads eventually to integrate out all quantities that vary on the atomic scales.
A key feature of the quasiclassical approximation is that G is sharply peaked at the Fermi surface, and depends weakly on energies far away from it. We use this assumption to split the propagator into low and high energy parts, G = G low + G high , where G low ( p, ε n ) = G( p, ε n ) for |ε| < ε c and G low ( p, ε n ) = 0 for |ε| > ε c . The cutoff energy ε c is taken to be ε c ≪ ε F . As shown in Fig. 2, we introduce the renormalized vertices (filled circles) that sum an infinite set of diagrams composed of the high energy part of the propagator and the bare vertices (open circles). The low energy part of the propagator obeys the Dyson equation, where G low 0 denotes the low energy part for the free propagator. This equation will play an important role in the following discussion.
For the convenience of the analysis, we define the quasiclassical propagator for the low energy part, g( p F , ε n ), as an integral over a shell, |v F · ( p − p F )|, in momentum space near the Fermi surface: where The propagator is normalized by dividing by the weight of the quasiparticle pole in the spectral function, a. The quasiclassical propagator matrix is parameterized as where f 0 and f represent the spin-singlet and spin-triplet components of anomalous propagators, and g 0 and g represent the spin-singlet and spin-triplet components of normal propagators. The Matsubara propagators maintain the following sets of symmetry relations in the NG space, To convert the Dyson equation (9) to a transport-like equation, we first perform the "left-right subtraction trick" for quasiclassical propagators, i.e., The inverse Green's function for free fermions, G −1 0 , is replaced with a −1 (ε − ξ( p)τ 3 ) if we include renormalization of the normal propagator by the zeroth order self-energy in the small parameter, i.e., T c /T F ≪ 1 or ε c /ε F ≪ 1. We remember that a is the weight of the quasiparticle pole in the spectral function. The kinetic equation is then reduced to (16) where R is the center-of-mass coordinate for the two neutrons. An important property of the self-energies is their weak dependence on momentum. We suppose that their characteristic momentum scale is set by the Fermi momentum, p F . For the quasiclassical renormalized perturbation, we can introduce v ext and σ MF , which are related to an external potential V ext and self-energy Σ taken at the Fermi level by respectively, with p ≈ p F . The factors, a and τ 3 , are included in v ext and σ MF for convenience. After the ξ p -integration, Eq. (16) reduces to which is the equation to determine the quasiclassical propagator g. Notice that Eq. (18) holds for homogeneous systems.
The mean-field self-energies σ MF are composed of the Fermiliquid corrections (diagonal parts) and the pair potentials (offdiagonal parts) as The spin-triplet pair potentials are parametrized by and∆ in terms of the three-dimensional vector which is called the d-vector for the spin-triplet superfluidity.
In the above equation, the sum is taken over µ. The explicit form of the d-vector will be expressed by the Green's function in Eq. (37) in the next subsection. The solution of Eq. (18) is not uniquely determined per se, because a + bg satisfies the same equation as g (a and b are arbitrary constants). To determine uniquely a solution for g, Eq. (18) must be supplemented by the normalization condition on the quasiclassical propagator: g 2 = −π 2 I with the unit matrix I in the NG space.

C. Mean-field self-energies and self-consistent equations
The interaction between two neutrons is modified by the "polarization effects". In the vicinity of the Fermi surface, a perturbation that couples to the quasiparticle states generates a polarization of the fermionic vacuum. Such polarization leads to a correction to the self-energy with respect to the energy of a fermionic quasiparticle. The leading-order correction is given by mean-field interaction energy associated with a particle-hole excitation. As we mentioned above, the two-body interaction between fermionic quasiparticles is represented by a four-point renormalized vertex (Fig. 2), (22) which is composed of the amplitudes for spin-independent scattering (Γ (s) ) and spin-dependent exchange scattering (Γ (a) ), where p ≡ ( p, ε). It sums the bare two-body interactions to FIG. 2. Leading order contributions to quasiclassical self-energies. Filled vertices, which couple to low-energy propagators (solid lines), show the particle-hole and particle-particle vertices, Γ ph and Γ pp , that sum all orders of the bare interaction (open circle) and highenergy intermediate states, G high (dashed lines). The particle-hole and particle-particle vertices determine the leading-order quasiparticle self-energy and pair potential, respectively. all orders involving all possible intermediate states of highenergy fermions. We use the notation ( p, −p) and ( p ′ , −p ′ ) to stand for the in-coming and out-going momenta, respectively, for the fermion 1 and 2. As the quasiclassical approximation takes account of quasiparticles confined to a low-energy shell near the Fermi surface, the vertex function can be evaluated with p = p F and ε → 0. The resulting vertex function reduces to Then, the scalar (Σ 0 ) and vector (Σ) components in the diagonal parts in the mean-field self-energy σ MF , Eq. (19), are determined as respectively, where n denotes the Matsubara sum with the cutoff energy ε c , and is the average over the Fermi surface with and φ ′ p for p ′ . Let f ( p 1 , p 2 ) be the particle-hole interaction between two nucleons which is generalized in spin (σ) and isospin (τ) spaces as where p 1 ( p 2 ) are the three-dimensional momentum for the in-coming (out-going) particle, σ 1 (σ 2 ) and τ 1 (τ 2 ) stand for SU(2) spin and isospin interactions for the particle 1 (2). The first two terms with F( p 1 , p 2 ) and F ′ ( p 1 , p 2 ) represent symmetric (spin-independent) quasiparticle scattering processes, while the latter two terms with G( p 1 , p 2 ) and G ′ ( p 1 , p 2 ) represent the antisymmetric (spin-dependent) quasiparticle scattering processes. The single-particle momenta are taken at the Fermi surface, | p i | ≈ | p F | for i = 1, 2. We introduce the factor N F = mp F /(2π 2 ) for the density-of-state of the fermion on the Fermi surface, so that the parameters F, F ′ , G, and G ′ are dimensionless quantities. Near the Fermi surface, F, F ′ , G, and G ′ (= G) are approximately regarded as functions only of the angle between p 1 and p 2 , and thus they can be expanded in terms of the Legendre polynomials P ℓ (x) (ℓ = 0, 1, 2, · · · ): for the function G( p F · p ′ F ) depending on p F · p ′ F . Here G ℓ is the coefficient in the channel ℓ. For neutron matter, one has τ 1 · τ 2 = 1 and the interaction potentials can be reduced to a more compact form by defining F (n) = F + F ′ and G (n) = G + G ′ . The superscript (n) stands for the neutron matter. The self-energies describe the Fermi liquid corrections due to symmetric (A (s) ) and antisymmetric (A (a) ) quasiparticle scattering processes. The symmetric and antisymmetric quasiparticle scattering amplitudes are parametrized with the Landau's Fermi-liquid parameters F (n) ℓ and G (n) ℓ as for the spin-symmetric case, and for the spin-asymmetric case. We notice that, among several coefficients, F (n) ℓ=1 and G (n) ℓ=0 give the Fermi liquid corrections to mass and spin susceptibility of a free neutron.
By taking into account the high-energy vertex corrections, the Zeeman energy in Eq. (5), −(1/2)γ n σ· B, in the NG space is recast into with the factor 1/(1 + G (n) 0 ). Here we introduce the magnetization density as a sum of the magnetization in the normal state and the correction by the superfluid state. The first term is explicitly given by M N = χ N B, where χ N = (1/2)γ 2 n N F /(1 + G (n) 0 ) is the spin susceptibility renormalized by the Fermi-liquid correction (G (n) 0 ) in the normal state. The nonvanishing magnetization density M(R) is fed back to the effective magnetic field B eff through the Fermi-liquid correction (G (n) 0 ), by referring Eqs. (31) and (32). Thus, the effective magnetic field including the corrections of spin-polarization density is given by This gives rise to a nonlinear effect of the Zeeman magnetic field.
In the same manner, the polarization effects exist for the four-fermion vertex which is denoted by Γ pp ab;cd ( p, ε; p ′ , ε ′ ) with the three dimensional momentum p ( p ′ ), the energy ε (ε ′ ), and spin indices a, c =↑, ↓ (b, d =↑, ↓) for the in-coming (out-going) particles. This vertex is irreducible in the particleparticle channel that sums bare two-body interactions to all orders involving all possible intermediate states of high-energy fermions (Fig. 2). As fermion pairs with binding energy |∆| ≪ ε c are confined to a low-energy band near the Fermi surface |ε| ≤ ε c ≪ ε F , the particle-particle vertex varies slowly on p in the neighborhood of the Fermi surface. Thus, the vertex reduces to functions only of the relative momenta, with the Fermi momenta p F and p ′ F . The particle-particle vertex function is decomposed into the spin-singlet (e: even parity) and spin-triplet (o: odd parity) functions for the particleparticle channels: µν (iσ y σ ν ) cd . (36) Using the effective interaction potential, one obtains the gap equation with µ = 1, 2, 3, which determines the equilibrium d-vector, d = (d 1 , d 2 , d 3 ). Here we have taken only the negative part of the 3 P 2 channel (odd parity) as the effective pairing interaction for dense neutrons, and have discarded the even-parity channel, i.e., V (e) = 0. The interaction in the 3 P 2 channel is supposed to be the short-range one so that the momentum dependence can be safely neglected.
For the representation of the interaction potential, let us introduce the spherical tensors, {t (m) µi } m=−J,··· ,+J with µ, i = 1, 2, 3, that form bases for representations of the rotational symmetry SO(3). Here m = −J, · · · , +J are the eigenvalues of J z . For the 3 P 2 channel, i.e., J = 2, the interaction potential can be expressed as the separable form of the symmetric and traceless tensors as where v > 0 is the coupling constant of the zero-range attractive 3 P 2 interaction. Here p F,i (p ′ F, j ) is the ith ( jth) component of the three-dimensional momentum p F ( p ′ F ) for the in-coming (out-going) states of the scattering neutron. It is important that the momentum dependence appears because the P-wave interaction potential is adopted. Eq. (38) is recast into with the traceless and symmetric tensor which obeys T µν ( p F ) = T νµ ( p F ) and tr(T ( p F )) ≡ µ T µµ ( p F ) = 0 [100]. The order parameter of spin-triplet superfluids, d( p F , R), is parameterized as with the rank-2 tensor A µi , where the index µ (i) denotes the spin (orbital) degrees of freedom of the Cooper pair and we have introducedp = (p 1 ,p 2 ,p 3 ) withp ≡p F /p F . In general, a rank-2 tensor can be expanded as a sum of the terms of the total angular momentum J = 0, 1, and 2 as where the scalar function A (0) ≡ tr(A)/3, the antisymmetric matrix A (1) µi ≡ (A µi − A iµ )/2, and the symmetric traceless matrix A (2) µi ≡ (A µi + A iµ )/2 − 1 3 δ µi tr(A) are the eigenstates of J = 0, 1, and 2, respectively. Thus, the number of independent components in A µi is then given as 3⊗3 = 1⊕3⊕5, where the numbers in the right-hand side represent the multiplicities of eigenstates of the total angular momentum J = 0, 1, 2, respectively. In the followings, we consider the neutron 3 P 2 superfluidity. Thus, we neglect the J = 0 and 1 components (A (0) = A (1) µi = 0), and express the tensor of the 3 P 2 order parameter by A µi = A (2) µi . A µi is then determined by solving the gap equation with where the average calculation for the momentum has been adopted by Eq. (26). In calculating Eqs. (43) and (44), we utilize the fact that the cutoff energy ε c and the coupling constant v are related to measurable quantity, i.e., the bulk transition temperature T c , through linearized gap equation Eliminating ε c and v from the above gap equation, Eq. (44) reduces to This is free from the ultraviolet divergence. Thus, the regularization of the gap equation leads to rapidly convergent series defined in terms of T c .

D. Thermodynamic potential
The thermodynamic potential in Eq. (6), which is the Φfunctional, generates the diagonal components and the offdiagonal components in the self-energy (8). To derive the thermodynamic potential within the quasiclassical approximation, we subtract the normal-state contributions from the Luttinger-Ward functional as ∆Ω ≡ Ω[G, Σ] − Ω[G N , Σ N ], where G N and Σ N are the Green's function and self-energy in the normal state, respectively. In this approximation, the Luttinger-Ward thermodynamic potential is then given by [101,102] where ∆Φ[g] is the Φ-functional confined to the low-energy region of the phase space. In the diagrammatic representation, ∆Φ is formally constructed by a number of low-energy propagators (G low ), and the higher-energy propagator (G high ) is renormalized into vertices as in Fig. 2. In Eq. (47) we set The quasiclassical auxiliary function g λ is given by replacing σ MF → λσ MF in Eq. (18). Here we determine the quasiclassical Φ-functional so as to consistently generate the self-energy through the functional derivative, σ MF = 2δ∆Φ[g]/δg T . It is found that the Φ-functional is constructed as which generates the self-consistent equations (24), (25), and (37).

E. Critical exponents and new universality class
The superfluid states subject to the total angular momentum J = 2 are classified into several phases: Uniaxial/biaxial nematic (UN/BN) phases, the ferromagnetic phase, and the cyclic phase [62,103,104]. The nematic phases preserve the time-reversal symmetry (TRS) and occupy the almost region of the phase diagram under a uniform magnetic field, while the latter two are nonunitary states with spontaneously broken time-reversal-symmetry. The ground state at the weak coupling limit is the uniaxial/biaxial nematic phase in which the rank-2 tensor A µi is represented by  0 . The consequence of the CEP is captured by thermodynamic quantities. First, in Fig. 3(b), we plot the heat capacity in the superfluid state per volume, C(T, B), which is obtained from the Luttinger-Ward thermodynamic potential, ∆Ω[g], as where the heat capacity of the normal gas of neutrons is given by C N (T ) = (2π 2 /3)N F k 2 B T . The heat capacity contains critical information on the thermal evolution of neutron stars [105]. The heat capacity shows the jump at the lower T . Another quantity which captures the consequence of the CEP is the magnetization M. This is defined as the first derivative of ∆Ω, which coincides with Eq. (32). It is seen from Fig. 3(c) that the T -dependence of M has the jump in the higher B region, indicating the first-order phase transition from the D 2 -BN phase to the D 4 -BN phase. The jump in M decreases as the magnetic field approaches the CEP. The discontinuity of M implies the divergence of the spin susceptibility, To extract the critical behaviors of the 3 P 2 superfluids, we compute the critical exponents around the CEP at (T cep , B cep ). We note that the contributions of the normal gas of neutrons to Eqs. (51) and (52), C N (T ) and M N (B), are negligible in the vicinity of the CEP, and the critical behaviors of the heat capacity C V , the magnetization M, and the spin susceptibility χ, are governed by the superfluid contributions, Then, we consider the set of the critical exponents (α, β, γ, δ) from the scaling behavior, which are parametrized by for T < T cep and B < B cep . Under the scaling hypothesis, the set of the critical exponents, (α, β, γ, δ), satisfies three equalities in Eqs. (1)-(3), i.e., Rushbrooke, Griffiths, and Widom equalities, that should hold at the CEP for any systems irrespective to the different interactions and dimensions. These three relations relate the critical exponents of magnetic systems, which have the endpoint of a first-order phase transition in non-zero temperatures. Figure 4 shows the scaling behavior of the specific heat C V (T, B), the magnetization M(T, B) and the spin susceptibility χ(T, B) around the CEP, (T cep , B cep ), which are directly computed with self-consistent solutions of the superfluid Fermi liquid theory. From these data, we read the values of the critical exponents as for G (n) 0 = −0.7. We find that the values of Eq. (61) satisfy the three equalities Eqs. (1)-(3) within the error range of 10% at most: indicating that the superfluid Fermi liquid theory properly captures the critical behavior of the CEP in neutron 3 P 2 superfluids. For G (n) 0 = −0.4, we read (α, β, γ, δ) = (0.60, 0.45, 0.59, 2.3) from the data in Fig. 4, which satisfy the above equalities within the error range of 10% at most, such that α + 2β + γ = 2.09, α + β(1 + δ) = 2.09, and − γ β + δ = 0.99. In Table I, we summarize the values of the critical exponents, {α, β, γ, δ}, for the Landau parameters G (n) 0 = −0.7 and −0.4. It turns out that the resulting exponents are insensitive to the screening effect of the external magnetic field due to the spinpolarized molecular field.
We propose that the set of the critical exponents, (α, β, γ, δ), in Eq. (61) belongs to a new type of university class. First of all, one may notice the large value of α (α ∼ 0.6). It is known that the value of α is usually much smaller than one in the phase transitions accompanied with a continuous symmetry breaking at least in known models thus far. However, the value of α can be larger in the cases accompanied with discrete symmetry breaking, such as the Potts model in two dimensions [94,95]. In our case, the CEP appears in the phase transition with a discrete symmetry breaking, i.e., from D 2 to D 4 . Thus, it may be natural to have the large value of α in the neutron 3 P 2 superfluid. Phenomenologically, the large α indicates that the heat capacity is much enhanced at the CEP (cf. Eq. (57)), which may affect the cooling process in the evolution of neutron stars. Naively to say, the large heat capacity will lead to a slow cooling in the evolution of neutron stars.
Another feature of the critical exponents in Eq. (61) is that the value of γ is smaller than one (γ ∼ 0.5). In the literature, there are only a limited number of examples which indicate γ < 1. One example is the O(n) model with n < 0 [96]. The O(n) model induces the Ising model at n = 1 and the percolation model at n = 0. If the value of γ is expressed in the asymptotic series up to the second-order terms in the vicinity of four dimensions, it is found that γ can be smaller than one if n is extrapolated to the negative region (n < 0). Another example for γ < 1 is the tricritical Ising model coupled to massless Dirac fermions [97]. In conclusion, the large α and the small γ are the unique feature of the critical exponents in the neutron 3 P 2 superfluid, implying a new universality class.

III. GINZBURG-LANDAU THEORY FOR THE CRITICAL ENDPOINT
A. Ginzburg-Landau free energy Let us turn to a discussion based on the Ginzburg-Landau (GL) theory [27, 28, 56-61, 81, 90-92]. In the weak coupling limit for the neutron-neutron interaction, we obtain the GL free energy density as an expansion series in terms of the condensate A µi and the magnetic field B [90,92]. We have adopted the quasiclassical approximation for the momentum integrals for the neutron loops. Notice that the free energy part for the noninteracting neutron is not included, because they are irrelevant to the the condensate. Each term in Eq. (65) with the derivative ∇ i for the spatial direction i = 1, 2, 3 and the GL coefficients defined by ζ(n) is the zeta function. In the above expression, µ n has been replaced to µ * n = (γ n /2)σ/(1+G (n) 0 ), i.e., the magnetic momentum of a neutron modified by the Landau parameter G (n) 0 . We notice that the Landau parameter stems from the Hartree-Fock approximation which are not taken into account explicitly in the present procedure for the fermion-loop expansion. The choice of the value of G (n) 0 does not affect the values of the critical exponents, because the magnetic field is scaled uniformly. This is different from the analysis in the BdG equation in Sec. II, where the spin-polarization leads to the non-linear effect for the Zeeman magnetic field (cf. Eq. (34)).
We notice that the β (4) and γ (2) terms were derived beyond the leading-order term for the magnetic field [90], and the δ (0) term was calculated to recover the global stability of the ground state which was absent at the 6th order [92]. We emphasize that, as discussed in detail in Ref. [92], the δ (0) term (the 8th order term) produces the first-order transition in the GL equation, which was absent in the analysis up to the 6th order term, and hence it leads to the existence of the CEP in the GL equation. In the derivation of the GL equation, we have supposed that the temperature T is close to the critical temperature at zero magnetic field T c , and hence the applicable region of the GL equation is limited in |1 − T/T c | ≪ 1. Notice that the critical temperature is the unique energy scale in the GL theory in the above.
B. Critical endpoint of 3 P 2 superfluid phase diagrams Let us consider the phase diagram drawn by the variational calculation for the GL free energy. For the magnetic field as B = (0, 0, B), we consider to minimize the GL free energy with respect to ∆ and r in Eq. (50). We show the phase diagram on the plane spanned by the temperature (T/T c ) and the magnetic field (γ n B/(πT c )) in Fig. 1(d). The CEP is (T cep /T c , γ n B cep /(πT c )) = (0.774597, 0.004465). The phase boundary at T < T cep and B < B cep is the first order transition, as indicated by the cyan lines in the figure. It should be noted that the existence of the CEP is due to the 8th order term as it was discovered in the previous work [92]. Thus, the GL theory shares common properties with the BdG theory about the existence of the CEP, though the positions of the CEP are different.
Those values agree with the values in the right-hand-sides in Eqs. (1), (2), and (3) by the exact values within the 10% numerical error.
In Table I, we summarize the values of the critical exponents from the BdG equation with different values of G (n) 0 and the ones from the GL equation. Interestingly, we observe that they are close to each other, even though the positions of the CEP on the T -B plane are different (cf. Fig. 1). The coincidence between the two suggests that the GL equation up to the 8th order term (δ 0 term) captures the essence of the CEP of the neutron 3 P 2 superfluid. Thus, the GL equation also supports a new universality class discussed in Sec. II E.

IV. SUMMARY AND DISCUSSION
We have discussed the critical exponents at the CEP in the phase diagram of the neutron 3 P 2 superfluidity, which can exist inside of neutron stars. Adopting the BdG equation with TABLE I. Critical exponents (α, β, γ, δ) computed by the superfluid Fermi liquid theory with G (n) 0 = −0.7 and −0.4 in the BdG theory and the GL theory. The Fermi liquid correction with G (n) < 0 leads to the screening effect of a magnetic field due to spin-polarized molecular field. Notice that the values of the critical exponents in the GL theory are independent of G (n) 0 . the spin-polarization effect, we have obtained the critical exponents and have confirmed that they satisfy the universal relations, i.e. the Rushbrooke, Griffiths, and Widom equalities, which hold for the spin systems. We have argued that the set of the critical exponents with large α and small γ belongs to a new universality class. One of the interesting features of the obtained critical exponents is the large α and small γ, in which the former indicates the slow cooling in the evolution of the neutron stars. We have checked that the spin-polarization effect induces the unique values of the critical exponents within the 10% numerical errors. We also have investigated the critical exponents in the GL equation up to the 8th order, and have confirmed that they satisfy the same universality relations, again within 10% errors. In spite of the different locations of the CEP in the phase diagram for the BdG equation and for the GL equation, we have found that the values of the critical exponents from the GL equation are properly regarded to be the same to the ones from the BdG equation, although there are still some discrepancies between the two due to the limited number of terms in the GL equation. Thus, we reach the conclusion that the GL equation up to the 8th order captures correctly the behaviors at the CEP in the neutron 3 P 2 superfluids. For more advanced study in future, we leave comments on the Fermi liquid parameter G (n) 0 in dense neutron matter. Bäckman et al. [106] computed the Fermi liquid parameters including the spin channel and isospin channel of quasiparticle scattering processes and find that G (n) 0 in Eq. (29) may be negligible. This result stems from the short-range property of the ρ-meson exchange interaction. It will be important to more carefully study the short-range behavior of the nucleonnucleon interaction, which will be attributed by the quarkexchange contributions at the nucleon core, the core polarization in the nuclear medium, and so on. As for another question, we may consider how the evolution of the neutron stars are influenced by the enhancement of the heat capacity, the magnetization, and the spin susceptibility at the CEP. Those information will be useful to research the internal structures of the neutron stars through the astrophysical observations.