Reexamining Ginzburg-Landau theory for neutron $^3P_2$ superfluidity in neutron stars

TL;DR

This paper extends the Ginzburg-Landau theory for neutron $^{3}P_{2}$ superfluidity in neutron stars by deriving an 8th order expansion, revealing a first-order phase transition that impacts neutron star interior models.

Contribution

It introduces an 8th order Ginzburg-Landau expansion for neutron $^{3}P_{2}$ superfluidity, showing the necessity of higher-order terms for a unique ground state and identifying a first-order phase transition.

Findings

8th order terms are necessary for a unique ground state.

Discovery of a first-order phase transition in the phase diagram.

Phase transition impacts neutron star interior structure.

Abstract

The Ginzburg-Landau (GL) effective theory is a useful tool to study a superconductivity or superfluidity near the critical temperature, and usually the expansion up to the 4th order in terms of order parameters is sufficient for the description of the second-order phase transition. In this paper, we discuss the GL equation for the neutron $^{3}P_{2}$ superfluidity relevant for interior of neutron stars. We derive the GL expansion up to the 8th order in the condensates and find that this order is necessary for the system to have the unique ground state, unlike the ordinary cases. Starting from the $LS$ potential, which provides the dominant attraction between two neutrons at the high density, we derive the GL equation in the path-integral formalism, where the auxiliary field method and the Nambu-Gor'kov representation are used. We present the detailed description for the trace calculation necessary in the derivation of the GL equation. As numerical results, we show the phase diagram of the neutron $^{3}P_{2}$ superfluidity on the plane spanned by the temperature and magnetic field, and find that the 8th order terms lead to a first-order phase transition, whose existence was predicted in the Bogoliubov-de Gennes equation but has not been found thus far within the framework of the GL expansion up to the 6th order.The first-order phase transition will affect the interior structures inside the neutron stars.

Equations

$$\displaystyle{\cal L}_{0}=\varphi^{{\dagger}}(t,\boldsymbol{x})\biggl{(}i\partial_{t}+\frac{\boldsymbol{\nabla}^{2}}{2m}+\mu+\boldsymbol{\mu}_{n}\!\cdot\!\boldsymbol{B}\biggr{)}\varphi(t,\boldsymbol{x}),$$

$$\displaystyle{\cal L}_{\mathrm{int}}=\sum_{i,j,k,l}\varphi_{\boldsymbol{q}i}^{\ast}\varphi_{\boldsymbol{p}j}\bigl{(}V_{LS}\bigr{)}_{ijkl}\varphi_{-\boldsymbol{q}k}^{\ast}\varphi_{-\boldsymbol{p}l},$$

$$\displaystyle\bigl{(}V_{LS}\bigr{)}_{ijkl}=i2G^{\prime}(\boldsymbol{S})_{ijkl}\!\cdot\!\bigl{(}\boldsymbol{p}\times\boldsymbol{q}\bigr{)},$$

$${\phi }_{\mathbf{p}}=\int {\mathrm{d}}^3\mathbf{x}{e}^{-i\mathbf{p}\cdot \mathbf{x}}\phi (t,\mathbf{x}),$$

$$\int \frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}\frac{{\mathrm{d}}^3\mathbf{q}}{(2\pi {)}^3}{\mathcal{L}}_{\mathrm{int}}$$

$$S(t,\mathbf{x})$$

$$V^a(t,\mathbf{x})$$

$$T^{ab}(t,\mathbf{x})$$

$$\displaystyle\phi^{ab}(t,\boldsymbol{x})\equiv-\varphi(t,\boldsymbol{x})^{t}\Sigma^{a{\dagger}}\bigl{(}\nabla^{b}_{x}\varphi(\tau,\boldsymbol{x})\bigr{)},$$

$${\Sigma }^0\equiv i{\sigma }^2,{\Sigma }^1\equiv i{\sigma }^1{\sigma }^2=-{\sigma }^3,{\Sigma }^2\equiv i{\sigma }^2{\sigma }^2=i1,{\Sigma }^3\equiv i{\sigma }^3{\sigma }^2={\sigma }^1,$$

$$G^{\prime }{X}^{†}(t,\mathbf{x})X(t,\mathbf{y})\to G{X}^{†}(t,\mathbf{x})X(t,\mathbf{x}){\delta }^{(3)}(\mathbf{x}-\mathbf{y}),$$

$$\displaystyle-G^{\prime}\int\mathrm{d}^{3}\boldsymbol{x}\,\mathrm{d}^{3}\boldsymbol{y}\,\biggl{(}\frac{2}{3}S^{{\dagger}}(t,\boldsymbol{x})S(t,\boldsymbol{y})+\frac{1}{2}\sum_{a}V^{a{\dagger}}(t,\boldsymbol{x})V^{a}(t,\boldsymbol{y})-\sum_{a,b}T^{ab{\dagger}}(t,\boldsymbol{x})T^{ab}(t,\boldsymbol{y})\biggr{)}$$

$$\displaystyle{\cal L}=\varphi^{{\dagger}}(t,\boldsymbol{x})\biggl{(}i\partial_{t}+\frac{\boldsymbol{\nabla}^{2}}{2m}+\mu+\boldsymbol{\mu}_{n}\!\cdot\!\boldsymbol{B}\biggr{)}\varphi(t,\boldsymbol{x})+G\sum_{a,b}T^{ab{\dagger}}(t,\boldsymbol{x})T^{ab}(t,\boldsymbol{x}),$$

$$\bar{Z}$$

$$\displaystyle\int{\cal D}A{\cal D}A^{{\dagger}}\exp\Biggl{(}-\frac{1}{G}\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,\sum_{a,b}A^{ab{\dagger}}(\tau,\boldsymbol{x})A^{ab}(\tau,\boldsymbol{x})\Biggr{)}$$

$$A^{ab}(\tau ,\mathbf{x})\to A^{ab}(\tau ,\mathbf{x})+G{T}^{ab}(\tau ,\mathbf{x}).$$

$$Z$$

$$\displaystyle\Psi(\tau,\boldsymbol{x})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\varphi(\tau,\boldsymbol{x})\\ \varphi^{\ast}(\tau,\boldsymbol{x})\end{array}\right),$$

$$\displaystyle Z=\int{\cal D}\Psi{\cal D}\Psi^{{\dagger}}{\cal D}A{\cal D}A^{{\dagger}}\exp\Biggl{[}-\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,\biggl{(}\Psi^{{\dagger}}(\tau,\boldsymbol{x})S(\tau,\boldsymbol{x})^{-1}\Psi(\tau,\boldsymbol{x})+\frac{1}{4G}\sum_{a,b}A^{ab{\dagger}}(\tau,\boldsymbol{x})A^{ab}(\tau,\boldsymbol{x})\biggr{)}\Biggr{]}.$$

$$\displaystyle S(\tau,\boldsymbol{x})^{-1}\equiv\left(\begin{array}[]{cc}\partial_{\tau}-\dfrac{\boldsymbol{\nabla}^{2}}{2m}-\mu-\boldsymbol{\mu}_{n}\!\cdot\!\boldsymbol{B}&-{\displaystyle\sum_{a,b}}A^{ab}(\tau,\boldsymbol{x})t^{ab{\dagger}}(\boldsymbol{x})\\ {\displaystyle\sum_{a,b}}A^{ab{\dagger}}(\tau,\boldsymbol{x})t^{ab}(\boldsymbol{x})&\partial_{\tau}+\dfrac{\boldsymbol{\nabla}^{2}}{2m}+\mu+\boldsymbol{\mu}_{n}^{t}\!\cdot\!\boldsymbol{B}\end{array}\right),$$

$$\displaystyle t^{ab}(\boldsymbol{x})\equiv-\biggl{(}\frac{1}{2}\Sigma^{a{\dagger}}\nabla^{b}+\frac{1}{2}\Sigma^{b{\dagger}}\nabla^{a}-\frac{1}{3}\delta^{ab}\sum_{c}\Sigma^{c{\dagger}}\nabla^{c}\biggr{)},$$

$$\displaystyle\int{\cal D}\Psi{\cal D}\Psi^{{\dagger}}\exp\Biggl{[}-\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,\biggl{(}\Psi^{{\dagger}}(\tau,\boldsymbol{x})S(\tau,\boldsymbol{x})^{-1}\Psi(\tau,\boldsymbol{x})\biggr{)}\Biggr{]}=\mathrm{Det}\,\bigl{(}S(\tau,\boldsymbol{x})^{-1}\bigr{)},$$

$$\displaystyle Z=\int{\cal D}A{\cal D}A^{{\dagger}}\exp\Biggl{[}\mathrm{Tr}\ln\bigl{(}S(\tau,\boldsymbol{x})^{-1}\bigr{)}-\frac{1}{4G}\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,\sum_{a,b}A^{ab{\dagger}}(\tau,\boldsymbol{x})A^{ab}(\tau,\boldsymbol{x})\Biggr{]},$$

$$\displaystyle W=-\mathrm{Tr}\ln\bigl{(}S(\tau,\boldsymbol{x})^{-1}\bigr{)}+\frac{1}{4G}\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,\sum_{a,b}A^{ab{\dagger}}(\tau,\boldsymbol{x})A^{ab}(\tau,\boldsymbol{x}).$$

$$\displaystyle S_{m}(\tau,\boldsymbol{x})^{-1}=\left(\begin{array}[]{cc}-i\omega_{N}+\xi_{\boldsymbol{p}}-\boldsymbol{\mu}_{n}\!\cdot\!\boldsymbol{B}&0\\ 0&-i\omega_{N}-\xi_{\boldsymbol{p}}+\mu+\boldsymbol{\mu}_{n}^{t}\!\cdot\!\boldsymbol{B}\end{array}\right),$$

$$\hat{A}(\tau ,\mathbf{x})$$

$$\displaystyle W=-\mathrm{Tr}\ln\bigl{(}S_{m}(\tau,\boldsymbol{x})^{-1}\bigr{)}-\sum_{N\geq 1}\frac{(-1)^{N+1}}{N}\mathrm{Tr}\Bigl{(}\bigl{(}S_{m}(\tau,\boldsymbol{x})\hat{A}(\tau,\boldsymbol{x})\bigr{)}^{N}\Bigr{)}+\frac{1}{4G}\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,\sum_{a,b}A^{ab{\dagger}}(\tau,\boldsymbol{x})A^{ab}(\tau,\boldsymbol{x}).$$

$$W=\int \mathrm{d}\tau {\mathrm{d}}^3\mathbf{x}w[A(\tau ,\mathbf{x})].$$

$$f[A(\tau ,\mathbf{x})]=\frac{1}{2}w[A(\tau ,\mathbf{x})],$$

$$f[A]=f_0+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{)},$$

Full text

I Introduction

Neutron stars are one of the important astrophysical objects providing us laboratories of nuclear physics (see Refs. Graber et al. (2017); Baym et al. (2018) for recent reviews). The extreme environments in the neutron stars, such as high density state, rapid rotation, strong magnetic field, strong gravitational field, and so on, lead to interesting questions about the unconventional states of nuclear systems. For example, it is considered that there can exist a various kind of matter phases inside the neutron stars: neutron rich gas and crusts at the surface, and neutron superfluidity, hyperon matter, $\pi$, $K$ condensates, quark matter, etc., in the inside. Those states can be accessible from observations of the mass-radius relation, time-evolution of the surface temperature, neutrino emissions, and so on. The research of the equation of state is still an open question, and many efforts have been devoted to understand the massive neutron stars whose masses are almost twice as large as the solar mass Demorest et al. (2010); Antoniadis et al. (2013). One of the most recent observational developments were provided by the gravitational waves from the binary neutron star merger Abbott et al. (2017). This finding opens a new era of researches of neutron stars in terms of the multimessenger astronomy including the microscopic physics, such as nuclear and elementary particle physics.

Among the various phases in the neutron stars, the research of the neutron superfluidity has bridged the connection between the theoretical knowledge and the astrophysical observations (see Refs. Chamel (2017); Haskell and Sedrakian (2018); Sedrakian and Clark (2018) for recent reviews). It was discussed that the neutron superfluidity affects relaxation time after pulsar glitches (sudden speed-up events of neutron-star rotation) Reichley and Downs (1969). It was also discussed that the neutron superfluidity is related to rapid cooling by neutrino emissions, which is known as the modified Urca process Yakovlev et al. (1999). It was proposed that pulsar glitch phenomena can be explained by pinning of large amount of superfluid vortices Anderson and Itoh (1975). From a more microscopic point of view, it is worthwhile to remind us that nuclear forces can supply the attractions in different channels for the baryon number density from low density to high density (see Ref. Dean and Hjorth-Jensen (2003) for a recent review). Migdal considered the ${}^{1}S_{0}$ channel as the most attractive interaction at low density Migdal (1960). However, the ${}^{1}S_{0}$ channel becomes repulsive at higher densities due to the strong repulsion at short ranges of the nuclear force Wolf (1966). At such high density, the attractive channel is given by the ${}^{3}P_{2}$ channel from the $LS$ potential.111We notice that the ${}^{3}P_{0}$ and ${}^{3}P_{1}$ channels are repulsive.

This attraction induces the neutron ${}^{3}P_{2}$ superfluidity as it was obtained by solving the gap equation Tabakin (1968); Hoffberg et al. (1970); Tamagaki (1970); Takatsuka and Tamagaki (1971); Takatsuka (1972); Amundsen and Ostgaard (1985); Takatsuka and Tamagaki (1993); 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). Thus, neutron ${}^{3}P_{2}$ superfluidity has been studied as one of the interesting states of the nuclear matter inside the neutron stars. Notice that, precisely to say, the neutron ${}^{3}P_{2}$ superfluidity exists with a small fraction of superconducting protons and normal electrons. As an astrophysical application, the rapid cooling time in neutron stars in Cassiopeia A was studied in terms of the the neutron ${}^{3}P_{2}$ superfluidity Heinke and Ho (2010); Shternin et al. (2011); Page et al. (2011), although the existence of the ${}^{3}P_{2}$ superfluidity is still elusive. It is also important to mention that, while the ${}^{1}S_{0}$ pairing is fragile in the presence of strong magnetic fields, the ${}^{3}P_{2}$ pairing is not. Therefore, one of possibilities to observe ${}^{3}P_{2}$ superfluidity is in magnetars, accompanied by large amount of magnetic fields. The magnetic field of magnetars reaches about $10^{15}$ G at the surface and may reach even $10^{18}$ G in the inside, and they are interesting astrophysical objects to study the neutron ${}^{3}P_{2}$ superfluidity under the superstrong magnetic field.222The origin of the strong magnetic fields was studied in terms of several mechanisms such as spin-dependent interactions between neutrons Brownell and Callaway (1969); Rice (1969); Silverstein (1969); Haensel and Bonazzola (1996), the pion domain wall Eto et al. (2013); Hashimoto (2015), the spin polarization in quark-matter core Tatsumi (2000); Nakano et al. (2003); Ohnishi et al. (2007), and so on. However, in recent years, the many-body calculation leads to a negative result for the realization of strong magnetic fields Bordbar and Bigdeli (2008).

The pairing in the neutron ${}^{3}P_{2}$ superfluidity is provided by the strong $LS$ potential between two neutrons as mentioned above. One of the particular features of the ${}^{3}P_{2}$ superfluidity is the richness of the symmetries of the condensate, which stems from the variety of the pairing patterns due to the various combinations between the relative angular momentum and the total spin. The ground state has the 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 baryon number, $S$ for spin rotation, $L$ for spatial rotation, $T$ for time-reversal symmetry, and $P$ for 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, ). It was discussed that, as a consequence of the symmetry breaking, low-energy excitations affect the cooling process by neutrino emission Bedaque et al. (2003); Bedaque and Nicholson (2013); Bedaque and Reddy (2014); Bedaque et al. (2015); Leinson (2010a, b, c, 2011a, 2011b, 2012, 2013, 2015). We comment that the cooling process is also related to quantum vortices Shahabasyan and Shahabasyan (2011). Recently the studies of the topological properties of the neutron ${}^{3}P_{2}$ superfluidity has been conducted: topological superfluidity and gapless Majorana fermions on the boundary of ${}^{3}P_{2}$ superfluids Mizushima et al. (2017), a quantized vortex Masuda and Nitta (2016), a soliton on it Chatterjee et al. (2017) and a half-quantized non-Abelian vortex Masuda and Nitta . Those studies were stimulated from the analogous states in condensed matter systems, such as $D$-wave superconductivity Mermin (1974), $P$-wave superfluidity in 3He liquid Vollhardt and Wölfle (2013); Volovik (2006), chiral $P$-wave superconductivity, e.g., in Sr2RuO4 Mackenzie and Maeno (2003), spin-2 Bose-Einstein condensates Kawaguchi and Ueda (2012), and so on. Thus, the research of the neutron ${}^{3}P_{2}$ superfluidity shares common interests in various condensed matter systems as an interdisciplinary subject.

The basic equation for the neutron ${}^{3}P_{2}$ superfluidity is the Bogoliubov–de-Gennes (BdG) equation  Tabakin (1968); Hoffberg et al. (1970); Tamagaki (1970); Takatsuka and Tamagaki (1971); Takatsuka (1972); Amundsen and Ostgaard (1985); Takatsuka and Tamagaki (1993); 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). Recently it was applied to investigate phase structures and topological properties in neutron ${}^{3}P_{2}$ superfluidity Mizushima et al. (2017). Around the critical temperature, the BdG equation can be reduced to the Ginzburg-Landau (GL) equation as the low-energy effective theory 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, ). The GL equation is obtained by integrating-out the fermionic degrees of freedom, and it is expressed in terms of the pairing gap behaving as a bosonic degree of freedom. It is known that, at the weak-coupling limit, the ground state of the neutron ${}^{3}P_{2}$ superfluidity is in the nematic phase, i.e., uniaxial nematic (UN) or biaxial nematic (BN) phase depending on the temperature and magnitude of the magnetic field Sauls and Serene (1978). The UN phase has an unbroken U(1) symmetry, while the BN phase is furthermore classified into two types according to symmetries, i.e., D2-BN and D4-BN phases, having an unbroken dihedral symmetries, D2 and D4, respectively. The UN phase is favored at zero magnetic field, while the BN phases are favored in the presence of finite magnetic fields Masuda and Nitta (2016); Mizushima et al. (2017), among which the D4-BN phase is favored above the critical magnetic field relevant for magnetars. Since the GL equation is much easier to deal with than the BdG equation, it has been applied to various problems although its use is limited to the vicinity of the critical temperature. For instance, it has been applied to study vortex structures of neutron ${}^{3}P_{2}$ superfluidity Richardson (1972); Muzikar et al. (1980); Sauls et al. (1982), inside whose core spontaneous magnetization was found Sauls et al. (1982); Masuda and Nitta (2016, ); Chatterjee et al. (2017). The GL equation is given in terms of the expansion series of the order parameter of the neutron ${}^{3}P_{2}$ superfluidity. Thus, it is an important question what order in the GL expansion is enough to capture physics of the neutron ${}^{3}P_{2}$ superfluidity. The 4th order terms are enough in usual cases but it is not so for the case of the neutron ${}^{3}P_{2}$ superfluidity as explained below. There is a continuous degeneracy in the ground states in the absence of the magnetic field as long as we consider the terms up to the 4th order for the order parameter. In other words, there is no energy difference among the UN, D2-BN, and D4-BN phases, and therefore we cannot determine the unique ground state.333At this order, there happens to exist an $\mathrm{SO}(5)$ symmetry as an extended symmetry of the potential term, which is absent in the Hamiltonian, and spontaneous breaking of such the symmetry leads a quasi-Nambu-Goldstone mode Uchino et al. (2010).

When we include the 6th order terms, only one phase is realized as the most stable state among them, i.e., the UN phase for zero magnetic field and the BN phases for finite magnetic field as mentioned above. The effect of the strong magnetic field beyond the leading-order term was analyzed in the phase diagram and the thermodynamic properties Yasui et al. (2019a, b). Such information will be useful to study the neutron ${}^{3}P_{2}$ superfluidity in magnetars. However, there is a serious drawback of the expansion up to the 6th order. The 6th order terms make the system unstable; There is the instability for large values of the order parameter and therefore the energy is not bounded from below.

In the present paper, to resolve this problem, we study the 8th order term of the order parameter which has not been studied so far. We find that the GL expansion up to the 8th order is the minimum order to uniquely determine the stable ground states (including the case of zero magnetic fields). We then investigate the phase diagram of the neutron ${}^{3}P_{2}$ superfluidity at this order, and find the existences of a 1st order phase transition line between the D2-BN and D4-BN phases and a tricritical point separating the 1st and 2nd order phase transition lines. Although a first-order line was previously predicted in the BdG equation under the quasiclassical approximation Mizushima et al. (2017), it was not found in the GL expansion up to the 6th order. Therefore, the GL expansion up to the 8th order is the minimum order to capture the essential physics of phase transitions in the ${}^{3}P_{2}$ superfluidity. As such, the first-order phase transition may affect the interior structures inside the neutron stars. We also present the detailed description about the derivation of the GL equation in the path-integral formalism, in which the auxiliary field method and the Nambu-Gor’kov representation are introduced.

This paper is organized as follows. In Sec. II, we introduce the $LS$ interaction between two neutrons, and give the tensor-type interaction relevant to the neutron ${}^{3}P_{2}$ superfluidity. In Sec. III, we give the detailed description for deriving the GL equation with the 8th order term by adopting the path-integral formalism. In Sec. V, we show the numerical results for the phase diagram of the neutron ${}^{3}P_{2}$ superfluidity, and discuss the effect of the 8th order term on the change of the phase diagram. The final section is devoted to our conclusion. In Appendix A, we present some details of the trace calculations used in the derivation of the GL equation. In Appendix B, we summarize symmetries of the 3P2 order parameters.

II Interaction Lagrangian

We introduce a nonrelativistic two-spinor field $\varphi(t,\boldsymbol{x})=(\varphi_{\uparrow}(t,\boldsymbol{x}),\varphi_{\downarrow}(t,\boldsymbol{x}))^{t}$ for a neutron. We consider the Lagrangian as a sum of the kinetic term and the the term for the interaction with the magnetic field,

[TABLE]

with the neutron mass $m$ and the chemical potential $\mu$. Notice that the interaction potential between the spin and the magnetic field is given by $\displaystyle V_{m}=-\varphi^{{\dagger}}\boldsymbol{\mu}_{n}\!\cdot\!\boldsymbol{B}\varphi$. Here $\displaystyle\boldsymbol{\mu}_{n}=-\gamma_{n}\hbar\boldsymbol{s}$ with the spin operator $\boldsymbol{s}=\boldsymbol{\sigma}/2$ with the Pauli matrices $\boldsymbol{\sigma}=(\sigma^{1},\sigma^{2},\sigma^{3})$ is the magnetic moment and $\gamma_{n}=1.2\times 10^{-13}$ MeV/T ($1\,\mathrm{T}=10^{4}\,\mathrm{G}$).

We consider the $LS$ interaction between two neutrons as the dominant force  Tabakin (1968); Hoffberg et al. (1970); Tamagaki (1970); Takatsuka and Tamagaki (1971); Takatsuka (1972); Amundsen and Ostgaard (1985); Takatsuka and Tamagaki (1993); 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). In the center-of-mass frame, we consider the scattering process $(\boldsymbol{p},j)+(-\boldsymbol{p},l)\rightarrow(\boldsymbol{q},i)+(-\boldsymbol{q},k)$ for the initial momentum $\boldsymbol{p},-\boldsymbol{p}$ and spin $j,l=\uparrow,\downarrow$ and the final momentum $\boldsymbol{q},-\boldsymbol{q}$ and spin $i,k=\uparrow,\downarrow$. Then, the interaction Lagrangian of the $LS$ interaction is provided by

[TABLE]

where we define the $LS$ operator defined by

[TABLE]

with the coupling constant $G^{\prime}>0$ and the spin operator $(\boldsymbol{S})_{ijkl}=(\boldsymbol{s})_{ij}\delta_{kl}+\delta_{ij}(\boldsymbol{s})_{kl}$ with $i,j,k,l=\uparrow,\downarrow$. For the later calculation, we transform the $LS$ interaction (2) in the momentum space to another form in the real space. For this purpose, adopting the Fourier transformation

[TABLE]

then we obtain

[TABLE]

where we have introduced scalar, vector and tensor fields by

[TABLE]

with $a,b,c=1,2,3$. Notice that the tensor field is symmetric and traceless. For the short notation, we define for the bilinear field for two neutrons

[TABLE]

with the spin operators defined by

[TABLE]

and the spatial derivative $\nabla^{b}_{x}=\partial/\partial x_{b}$. Here $\varphi^{t}\Sigma^{0}\varphi$ is invariant as the scalar and $\varphi^{t}\boldsymbol{\Sigma}\varphi$ with $\boldsymbol{\Sigma}=(\Sigma^{1},\Sigma^{2},\Sigma^{3})$ transforms as a vector. This is checked by the spin rotation of $\varphi$. We notice that Eq. (5) is the nonlocal interaction. For simplicity, we regard this interaction as the local one by replacing the nonlocal operator $X(t,\boldsymbol{x})=S(t,\boldsymbol{x})$, $V^{a}(t,\boldsymbol{x})$, and $T^{ab}(t,\boldsymbol{x})$ in the following way:

[TABLE]

where we have introduced a new coupling constant $G(>0)$ instead of the original one $G^{\prime}$. The mass dimension of $G$ is different from that of $G^{\prime}$. Then, Eq. (5) is reduced to

[TABLE]

Because the scalar and vector channels are repulsive for $G>0$ and these two channels will be irrelevant to the neutron superfluidity in a direct manner, we will neglect those two terms in the following discussions.444We keep it in mind that the magnetic field may cause the mixing between the scalar, vector, and tensor channels, when the strength of the magnetic field is strong enough. If such mixing occurs, then the tensor interaction may be more enhanced through the virtual processes with the scalar and vector channels. From the above consideration, we obtain the full Lagrangian relevant to the tensor interaction

[TABLE]

as the basic Lagrangian for the neutron ${}^{3}P_{2}$ superfluidity. In the next section, we will derive the GL free energy density from this Lagrangian.

III Derivation of the Ginzburg-Landau equation

III.1 Auxiliary field method

Based on the Lagrangian (13), we consider the generating functional. For the convenience of the calculation, we introduce the Wick rotation for the real time $t$ on the complex plane, and consider the imaginary time $\tau$ by replacing $it$ by $\tau$ in the Euclidean space and treat $\tau$ as a real parameter. Because we are interested in the finite temperature, we restrict the range of $\tau$ to the periodic range, $\tau\in[0,\beta]$ with the inverse temperature $\beta=1/T$. In the Wick rotation, $\partial_{t}$ is replaced to $\partial_{\tau}=\partial/\partial\tau$ whose relation is $\partial_{t}=i\partial_{\tau}$. The generating functional is given by

[TABLE]

in the path-integral formalism, where ${\cal D}\varphi^{\ast}$ and ${\cal D}\varphi$ indicate the path integrals for $\varphi^{\ast}$ and $\varphi$, respectively. The calculation of the path-integral is difficult to be performed in an analytic form due to the interaction term, and hence an approximate solution should be searched. For this purpose, we introduce a complex tensor field $A^{ab}(\tau,\boldsymbol{x})$ ($a,b=1,2,3$) as an auxiliary field, where $A^{ab}(\tau,\boldsymbol{x})$ is a symmetric and traceless tensor. It will turn out that $A^{ab}(\tau,\boldsymbol{x})$ is the order parameter relevant to the neutron pairing induced by the tensor term in Eq. (8). The inclusion of the auxiliary field into the Lagrangian is accomplished as follows. First, we multiply

[TABLE]

to the generating functional $\bar{Z}$, and shift $A^{ab}(\tau,\boldsymbol{x})$ by

[TABLE]

Notice that the quantity in Eq. (15) is a constant number, because the path-integral is performed for $A^{ab}$ and $A^{ab{\dagger}}$ as it is represented by ${\cal D}A$ and ${\cal D}A^{{\dagger}}$, respectively. Thus, the multiplication of Eq. (15) to the generating functional does not change the contents of physics. Notice also that the above shift of $A^{ab}(\tau,\boldsymbol{x})$ is always possible due to the arbitrariness of parametrizations of the auxiliary field. In the end, we obtain a new form of the generating functional

[TABLE]

where the neutron field $\varphi$ and the auxiliary field $A^{ab}$ are coupled to each other. $Z$ is physically the same as the original generating functional $\bar{Z}$. In the next subsection, we will integrate out the neutron field and construct the effective Lagrangian containing only the auxiliary field at the low-energy scale.

III.2 Nambu-Gor’kov representation

In performing the path-integral for the neutron field $\varphi$ in Eq. (17), it is convenient to represent $\varphi$ in the Nambu-Gor’kov representation. For this purpose, we define a four-component spinor field

[TABLE]

which contains the particle component $\varphi$ in the first and second columns and the hole component $\varphi^{\ast}$ in the third and fourth columns. With this field, the generating function can be expressed by

[TABLE]

Here, we have defined the inverse of the propagator in the Nambu-Gor’kov space by

[TABLE]

where we have introduced

[TABLE]

for the coupling of the spin and the momentum. It is convenient to notice $\boldsymbol{\sigma}^{t}=-\sigma^{2}\boldsymbol{\sigma}\sigma^{2}$ in the following calculation. The first term in the exponential in Eq. (21) is a bilinear form for $\Psi$ and $\Psi^{{\dagger}}$, and hence it is possible to perform the integral analytically. The result is

[TABLE]

where $\mathrm{Det}$ indicates the determinant of the matrix whose indices are given not only by the spin of the neutron but also by the imaginary time $\tau$ and the three-dimensional space $\boldsymbol{x}$. Furthermore, adopting the formula $\mathrm{Det}\,M=e^{\mathrm{Tr}\ln M}$ which is applicable to a matrix $M$, we can express the generating functional (21) as

[TABLE]

where only the path-integral for $A$ and $A^{{\dagger}}$ is left. However, the above form of the generating functional does not provide directly an answer to the question about the ground state of the system. This is because the first term in the exponential \mathrm{Tr}\ln\bigl{(}S(\tau,\boldsymbol{x})^{-1}\bigr{)} contains the field $A(\tau,\boldsymbol{x})$ and $A^{{\dagger}}(\tau,\boldsymbol{x})$, and hence it is still difficult to perform the path-integral for $A$ and $A^{{\dagger}}$. Instead of the complete calculation, we will expand \mathrm{Tr}\ln\bigl{(}S(\tau,\boldsymbol{x})^{-1}\bigr{)} for a small $A(\tau,\boldsymbol{x})$ and $A^{{\dagger}}(\tau,\boldsymbol{x})$, and will try to obtain approximately an analytical equation applicable for the small $A(\tau,\boldsymbol{x})$ around the critical temperature. We also expand \mathrm{Tr}\ln\bigl{(}S(\tau,\boldsymbol{x})^{-1}\bigr{)} in terms of a small magnetic field $\boldsymbol{B}$. In the previous work by ones of the authors, the higher-order terms of the magnetic field was investigated, and the influence of the strong magnetic field on the phase diagram was studied Yasui et al. (2019a, b). The present calculation is the extension to the 8th order of $A(\tau,\boldsymbol{x})$.

III.3 Effective potential

For the generating function $Z$ in Eq. (27), we define the effective potential $W$ by $Z=e^{-W}$ where only the one-particle irreducible diagrams are contained and the other diagrams are removed. By regarding $A(\tau,\boldsymbol{x})$ as a classical field (mean-field) and comparing $W$ with the functions in the exponential in Eq. (27), we obtain the effective potential

[TABLE]

Notice that the mean-field approximation was adopted, where the fluctuations about $A(\tau,\boldsymbol{x})$ are not considered and hence the strong-coupling effect is not included. In order to expand the trace in terms of $A(\tau,\boldsymbol{x})$ and the magnetic field $\boldsymbol{B}$, we decompose the inverse of the propagator as $\displaystyle S(\tau,\boldsymbol{x})^{-1}=S_{m}(\tau,\boldsymbol{x})^{-1}+\hat{A}(\tau,\boldsymbol{x})$. Here $S_{m}(\tau,\boldsymbol{x})^{-1}$ is the part including the free term and the spin-magnetic interaction term

[TABLE]

and $\hat{A}(\tau,\boldsymbol{x})$ is the part which contains the auxiliary field $A^{ab}(\tau,\boldsymbol{x})$

[TABLE]

By adopting the expansion for the logarithmic function, we finally obtain an expanded form of the effective potential

[TABLE]

We notice that trace, $\mathrm{Tr}$, takes the sum over all the variables such as spin, imaginary time and space. The concrete procedure of the trace calculation is presented in Appendix A. For the later convenience, we define the density of the effective potential $w[A(\tau,\boldsymbol{x})]$ by

[TABLE]

This quantity will be used in the next subsection.

III.4 Ginzburg-Landau free energy

By using the density of the effective potential $w[A(\tau,\boldsymbol{x})]$ in Eq. (36), we define the GL free energy $f[A(\tau,\boldsymbol{x})]$ by

[TABLE]

where the coefficient $1/2$ is multiplied in order to avoid the double-counting of the fermionic degrees of freedom in the Nambu-Gor’kov representation (cf. Eq. (20)). By adopting the expansion for the boson field $A^{ab}$ and the magnetic field $\boldsymbol{B}$ in Eq. (35), we obtain the free energy, given by

[TABLE]

as a series of $A$ and $\boldsymbol{B}$. Each function is defined as follows. $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$. The other terms include the condensate $A^{ab}$. $f_{8}^{(0)}[{A}]$ consists of the terms including the field $A^{ab}$ up to the 8th order with no magnetic field, $f_{2}^{(\leq 4)}[{A}]$ consists of the terms including the field $A^{ab}$ up to the 2nd order with the magnetic field up to $|\boldsymbol{B}|^{4}$, and $f_{4}^{(\leq 2)}[{A}]$ consists of the terms including the field $A^{ab}$ up to the 4th order with the magnetic field up to $|\boldsymbol{B}|^{2}$. Their explicit forms are given by

[TABLE]

Here, the GL coefficients can be obtained as

[TABLE]

with $N(0)=m\,p_{F}/(2\pi^{2})$ the state-number-density at the Fermi surface and $|\boldsymbol{\mu}_{n}^{\ast}|=(\gamma_{n}/2)/(1+F_{0}^{a})$ the magnitude of the magnetic momentum of a neutron modified by the Landau parameter $F_{0}^{a}$. The Landau parameter stems from the Hartree-Fock approximation which are not considered explicitly in the present formalism. $\zeta(n)$ is the zeta function. We use the condition that the temperature $T$ is close to the critical temperature at zero magnetic field $T_{c0}$; $|1-T/T_{c0}|\ll 1$. The critical temperature is related to the strength of the coupling constants through

[TABLE]

with the Euler’s constant $\gamma$, where the energy scale $D$ relevant for the typical scales in the dynamics is introduced as the cutoff parameter in the momentum space below and above the Fermi surface. We notice that the parameters of the coupling constant $G$ and the momentum cutoff $D$ are controlled by the critical temperature $T_{c0}$ through the above relation. Therefore, $T_{c0}$ is the only parameter to control the energy scale of the neutron ${}^{3}P_{2}$ superfluidity. We remind us that this simplification is just due to the one-loop approximation for the neutron part in the calculation of the GL free energy density. In the derivation of the above equations, we have used the quasi-classical approximation for the momentum integrals. The concrete calculations are shown in Appendix A. We emphasize that the 8th order terms, i.e., the terms with the coefficient $\delta^{(0)}$, provide a new ingredient which has not been calculated so far in the literature.

IV Stability of the ground state

In this section, as one of the most important applications of the 8th order terms that we have obtained in the last section, we show that the GL expansion up to the 8th order admits the stable and unique ground state once the temperature and magnetic field are fixed, while the expansion up to 6th order has the instability. We plot the GL free energy density, $\Delta f=f-f_{0}$ (the difference between the full equation of the GL free energy density and the part from no condensate terms), as a function of $A_{0}$ for fixed $t$, $b$, and $r$ in Fig. 1. We exemplify the two cases: (i) $t=0.9$, $b=0$, and $r=-1/2$, and (ii) $t=0.9$, $b=1$, and $r=-1$. The former denotes the UN phase, and the latter the D4-BN phase. The GL free energy function with $\delta_{0}=0$ is presented by the dashed lines in the figure. The (locally) minimum points of the GL free energy density are presented by the magenta blobs. As for the location of the minimum points, there is no much difference between the $\delta_{0}=0$ and $\delta_{0}\neq 0$ cases, as the behaviors for the smaller $A_{0}$ are observed to be almost the same for the both cases. However, we notice that the behaviors at larger $A_{0}$ are different to each other: $\Delta f$ has the instability in the $\delta_{0}=0$ case, and $\Delta f$ has no such instability in the $\delta_{0}\neq 0$ case. The former leads to the stability of the ground state as the local minimum, while the latter to the global minimum. The existence of the global minimum is important to ensure the stability of the ground state. It is also practically useful to search the minimal points in the variational calculations in the wide range of the parameters. Thus, the GL expansion up the 8th order is necessary for the system having the stable and unique ground state, and it helps us to investigate easily various properties of the neutron ${}^{3}P_{2}$ superfluidity.

V Phase diagram

With the GL free energy density (38), we show the phase diagram of the neutron ${}^{3}P_{2}$ superfluidity on the plane spanned by the temperature and the magnetic field. We consider the nematic phases in which the order parameter $A^{ab}$ is parametrized by

[TABLE]

with $A_{0}\geq 0$ being the magnitude of the condensate and $r$ is a real parameter which can be restricted to $-1\leq r\leq-1/2$ by using the $SO(3)\times U(1)$ symmetry without the loss of generality. We notice that it is always possible to transform any symmetric and traceless matrix into the diagonal form in which the diagonal elements are aligned in the order of the magnitudes of the absolute values. For the different values of $r$, the order parameter possesses different symmetries: the uniaxial-nematic (UN) phase ($r=-1/2$) and the biaxial-nematic (BN) phase ($-1\leq r<-1/2$). The UN phase has the U(1) symmetry. The BN phase is further classified to the D2-BN phase ($-1<r<-1/2$) and the D4-BN phase ($r=-1$), having the D2 and D4 dihedral symmetries, respectively. The explicit forms of group transformations of those symmetries are summarized in Appendix B. We substitute Eq. (48) to the GL free energy density (38), and perform the variational calculation with respect to $A_{0}$ and $r$. In the present study, we choose the direction of the magnetic field along the $y$-axis, $\boldsymbol{B}=(0,B,0)$, to minimize the total energy.

We plot the phase diagrams of $A_{0}$ and $r$ on the plane of the temperature ($t=T/T_{c0}$) and the magnetic field ($b=\gamma_{n}B/((1+F_{0}^{a})T_{c0})$) in Fig. 2. Notice that $t$ and $b$ are dimensionless quantities. In Figs. 2(a) and 2(c), we show the previous results for the metastable states without the 8th order terms by setting $\delta_{0}=0$ in Eq. (38). In Figs. 2(b) and 2(d), we show the new results with the 8th order terms by setting $\delta_{0}\neq 0$ as in Eq. (43). The upper and lower panels are for $A_{0}$ and $r$, respectively. In the case of $\delta_{0}=0$, the UN phase appears at the temperature lower than the critical temperature ($t<1$) and the zero magnetic field ($b=0$), and the D2-BN phase appears at the finite magnetic field, and it turns to be the D4-BN phase at the stronger magnetic field. The $\delta_{0}\neq 0$ case has almost no difference to the $\delta_{0}=0$ case around the critical temperature ($t\approx 1$). However, the particular difference appears at the lower temperature and the smaller magnetic field. In the $\delta_{0}=0$ case, the UN phase and the D2-BN phase are almost unchanged from the higher temperature to the lower temperature for the fixed magnetic field. In contrast, in the $\delta_{0}\neq 0$ case, the UN and D2-BN phases disappear at the low temperature, and the D4-BN phase appears. This change is the first-order phase transition. The transition lines are shown by the cyan lines in Fig. 2. The critical temperature at zero magnetic field is $t_{c}=0.772$. The existence of the first-order phase transition is also shown in Fig. 3. The discontinuous changes in $A_{0}$ and $r$ at $t_{c}=0.772$ for the zero magnetic field ($b=0$) can be seen clearly. The first-order phase transition has not been known in the analysis of the GL theory up to the 6th order for the neutron ${}^{3}P_{2}$ superfluidity, while it was found in the analysis based on the BdG equation with quasi-classical approximation Mizushima et al. (2017). Thus, the 8th order terms for $\delta_{0}\neq 0$ in the GL equation capture all phase transitions in the neutron ${}^{3}P_{2}$ superfluidity at least qualitatively.

VI Summary and discussion

We have discussed the GL free energy density for the neutron ${}^{3}P_{2}$ superfluidity in terms of the $LS$ potential between two neutrons. Adopting the path-integral formalism, the auxiliary field, and the Nambu-Gor’kov representation, we have derived the GL free energy density up to the 8th order terms of the condensate. We also have presented the detailed description of the calculation. We have found that the 8th order terms support the ground state as the global minimum, and it is practically useful for the variational calculation for the GL free energy density. This resolves the problem of the GL expansions up to the 6th order in which the system is unstable for large values of the order parameters. We have shown the phase diagram of the nematic phases, i.e., the UN phase, the D2-BN phase, and the D4-BN phase on the plane spanned by the temperature and the magnetic field. It has tuned out that the D4-BN phase appears at the low temperature and the smaller magnetic field, and that there is a first-order phase transition between the D4-BN and D2-BN phases. Those properties were not known in the previous studies in the GL free energy density, in which the maximum order in the condensate was limited to the 6th order. Therefore, the 8th order terms capture a qualitatively correct phase diagram for the neutron ${}^{3}P_{2}$ superfluidity.

We comment that it was known already in the BdG equation that the first-order phase transition exists at the low temperature and the finite magnetic field Mizushima et al. (2017). When we compare our result with the one from the BdG equation, however, we notice that the values of the temperature and the magnetic field in which the first-order phase transition appears are quite different: the transition line attaches to the temperature axis in our result, while it attaches to the axis of the magnetic field in the analysis by the BdG equation. We will report a further detailed analysis soon.

The existence of the first-order transition can affect the properties of the neutron stars. For example, it is one of the interesting questions to ask how the quantum vortices are affected by the first-order phase transition. For future work, a connection of the neutron ${}^{3}P_{2}$ phase to the other phases such as hyperon matter, quark matter and so on should be studied. When a neutron star rotates, a connection between 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 was discussed in the presence or absence of boojums which are endpoints of vortices Cipriani et al. (2012); Alford et al. (2019); Chatterjee et al. (2019a, b); Cherman et al. . The boojums, if they exist, may affect dynamical properties of the neutron stars. See Ref. Eto et al. (2014) for a review of non-Abelian quantum vortices in quark matter. It is also important to study cyclic and ferromagnetic phases, which are possible for general GL theory with total angular momentum two Mermin (1974); Kawaguchi and Ueda (2012). The former makes it possible to realize one-third quantized non-Abelian vortices Semenoff and Zhou (2007) which form a network in collision Kobayashi et al. (2009), while the latter could be relevant to the properties of the magnetic field in magnetars. The cyclic and ferromagnetic phases are topological matter (Weyl semimetals) which can have gapless fermions in the bulk Mizushima et al. (2017); Mizushima and Nitta (2018). It is also an interesting to apply our formulation to superfluidity to higher spin systems Venderbos et al. (2018).

It is also an interesting question to investigate the higher-partial wave component, i.e., ${}^{3}F_{2}$, which couples to the ${}^{3}P_{2}$ component through the tensor interaction Amundsen and Ostgaard (1985); Takatsuka and Tamagaki (1993); Zverev et al. (2003). The mixing between the ${}^{3}P_{2}$ component and the ${}^{3}F_{2}$ component induces the enhancement of the amplitude of the gap value and enhances the stability of the neutron superfluidity, and it can provide richer phases than the pure ${}^{3}P_{2}$ component Zverev et al. (2003). Those new phases should be studied in detail in future studies.

Acknowledgment

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). C. C. acknowledges support as an International Research Fellow of the Japan Society for the Promotion of Science (JSPS) (Grant No. 16F16322). This work is also supported in part by JSPS Grant-in-Aid for Scientific Research (KAKENHI Grant No. 16H03984 (M. N.), No. 18H01217 (M. N.), No. 17K05435 (S. Y.), and No. 19K14713 (C. C.)), and also by MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas “Topological Materials Science” (Grant No. 15H05855 (M. N.)).

Appendix A Trace calculation

We present some examples of the trace calculation of Eq. (35). We demonstrate the calculations of the log term in Sec. A.1 and the calculation of the power series in Sec. A.2 to derive the GL free energy density (38).

A.1 Log term

We follow the several steps in the calculation. First, we express the trace, ${\mathrm{Tr}}$, by the state $|\tau,\boldsymbol{x}\rangle$ in the real space (imaginary time and three-dimensional space). Second, we insert the complete set in the momentum space whose basis is denoted by $|i\omega_{n},\boldsymbol{p}\rangle$ with the Matsubara frequency $\omega_{n}$ and the three-dimensional momentum $\boldsymbol{p}$. Third, we utilize the plane-wave wave function $\displaystyle\langle i\omega_{n},\boldsymbol{p}|\tau,\boldsymbol{x}\rangle=e^{-i\omega_{n}\tau+i\boldsymbol{p}\cdot\boldsymbol{x}}$. With the above steps, we perform the calculation of the logarithmic term as

[TABLE]

where we have introduced the Matsubara frequency $\omega_{n}=(2n+1)\pi/\beta$ ($n\in\mathbb{Z}$) with the inverse temperature $\beta=1/T$, and the particle energy $\xi_{\boldsymbol{p}}^{\pm}\equiv\xi_{\boldsymbol{p}}\pm|\boldsymbol{\mu}_{n}||\boldsymbol{B}|$ with $\xi_{\boldsymbol{p}}\equiv{\boldsymbol{p}^{2}}/({2m})-\mu$. The summation of the Matsubara frequencies are performed as

[TABLE]

where we have used the trick to introduce the integral about $\xi$ and have used the relation

[TABLE]

for an analytic function $f(p_{0})$. The path $C$ turns around $i\omega_{n}$ counterclockwise. In the end of the above calculation, we obtain

[TABLE]

where the constant terms are neglected. We also have neglected the term \frac{1}{\beta}\ln\bigl{(}1+e^{\beta\xi_{\boldsymbol{p}}^{\pm}}\bigr{)} which can be regarded as a shift by treating the particle and hole on the equal footing in the Nambu-Gor’kov representation. By dividing the above quantity by the total volume $\displaystyle\int\mathrm{d}\tau\,\mathrm{d}^{3}\boldsymbol{x}\,1=\beta V$ ($V$ being the volume of the three-dimensional space) and by multiplying the factor $1/2$ for removing the double-counting in the Nambu-Gor’kov representation, we obtain $f_{0}$ which appears in Eq. (39).

A.2 Power series

We apply the essentially same procedure of the calculation of the power series ($N$th order) in Eq. (35). The different point is the necessity of the multiple number of states for the completeness relations, such as $|\tau_{1},\boldsymbol{x}_{1}\rangle$, $|\tau_{2},\boldsymbol{x}_{2}\rangle$, $\dots$, $|\tau_{N},\boldsymbol{x}_{N}\rangle$ in the real space and $|i\omega_{1},\boldsymbol{p}_{1}\rangle$, $|i\omega_{2},\boldsymbol{p}_{2}\rangle$, $\dots$, $|i\omega_{N},\boldsymbol{p}_{N}\rangle$ in the momentum space. The $N$th order term is calculated as follows:

[TABLE]

where we have used the completeness relations

[TABLE]

in the real space and

[TABLE]

in the momentum space with the Matsubara frequency $\omega_{n_{i}}=(2n_{i}+1)\pi/\beta$ for each $i=1,2,\dots,N$. The trace, $\mathrm{tr}$, is taken over for the spin of the neutron and the up and down-components in the Nambu-Gor’kov representation. In the above equations, we have introduced the propagator of a nucleon under the finite magnetic field

[TABLE]

and the matrix of the auxiliary field

[TABLE]

which are (partially) represented in the momentum space and the Matsubara frequency. The imaginary time and the positions $(\tau_{i},\boldsymbol{x}_{i})$ are different for $i=1,2,\dots,N$. As for the kinetic term with $N=2$, we use the expansion

[TABLE]

by assuming that the difference between $(\tau_{2},\boldsymbol{x}_{2})$ and $(\tau_{1},\boldsymbol{x}_{1})$ is small. The higher order terms in the derivatives are neglected. In the case of $N\geq 3$, we regard that all the different points $(\tau_{i},\boldsymbol{x}_{i})$ are identical to $(\tau,\boldsymbol{x})$, and apply the approximation

[TABLE]

for each $i=1,2,\dots,N$. We also adopt the expansion for the magnetic field $\boldsymbol{B}$ up to the appropriate order, such as

[TABLE]

whenever it is necessary.

We leave comments on the tips for the integrals for the three-dimensional momentum. Parametrizing $\boldsymbol{p}=(p\sin\theta\cos\varphi,p\sin\theta\sin\varphi,p\cos\theta)$ in the polar coordinate and assuming the isotropy in the momentum space, we calculate the polar integrals as

[TABLE]

with $\mathrm{d}\Omega_{\boldsymbol{p}}=\sin\theta\mathrm{d}\theta\mathrm{d}\varphi$, where $f(p)$ is a function of $p=|\boldsymbol{p}|$. We use the notation that $S^{a_{1}a_{2}\dots a_{n}}$ indicates the totally symmetric tensor for $a_{1},a_{2},\dots,a_{n}$ ($a_{i}=1,2,3$ for $i=1,2,\dots,n$), such as

[TABLE]

As for the integrals for $p$, we adopt the quasiclassical approximation

[TABLE]

with $N(0)={m\,p_{F}}/({2\pi^{2}})$ being the number of the state density at the Fermi surface, where we measure the energy from the chemical potential as $\xi_{\boldsymbol{p}}=\boldsymbol{p}^{2}/(2m)-\mu$. In the first approximate equality, we have replaced the measure part of the momentum integral from $p^{2}/(2\pi^{2})$ to $p_{F}^{2}/(2\pi^{2})$ with the Fermi momentum $p_{F}$, and in the second approximate equality, we have introduced $\xi_{\boldsymbol{p}}$ instead of $p$ and have determined the integral range by $-\infty<\xi_{\boldsymbol{p}}<\infty$. Such an extension of the integral range will be reasonable at low temperature. Now we take the summation over the Matsubara frequencies $\omega_{n}$ and perform the momentum integrals for $p$. For the term linear to $\mathrm{tr}\,A^{\ast}A$, we consider

[TABLE]

and obtain

[TABLE]

Here, we have introduced the momentum cutoff $D$ above and below the Fermi surface and have used the integral

[TABLE]

with the Euler’s constant $\gamma$. For the other terms, we adopt the following procedure

[TABLE]

for a positive integer $k$. Notice that the momentum integral is convergent without introducing the finite momentum cutoff $D$. We have used $\omega_{n}=(2n+1)\pi\,T$ for the temperature $T$, the Gamma function $\Gamma(x)$, and the zeta function $\zeta(x)=\sum_{n=1}^{\infty}1/n^{x}$.

Finally, we comment that the equations represented by trace, $\mathrm{tr}$, for the symmetric, traceless, and complex matrix $A$ with $3\times 3$ dimension have some identity relations. For $N=4$, we have the identity Sauls and Serene (1978)

[TABLE]

and for $N=6$ we have the identities

[TABLE]

which are sometimes useful to reduce the number of the terms in the final result. In the end of the above calculations, we obtain the GL free energy density (42).

Appendix B Symmetries

In this Appendix, we summarize the symmetries of the order parameter $A$ parametrized in Eq. (48). We remember that the symmetry of $A$ is given

[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}$). We notice that $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 above symmetries are spontaneously broken to subgroups in general in the nematic phase in Eq. (48). Instead, there are the U(1) symmetry, the D2 symmetry, and the D4 symmetry, which are the symmetries in the UN phase ($r=-1/2$), the D2-BN phase ($-1<r<-1/2$), and the D4-BN phase ($r=-1$). We will explain the concrete forms of those symmetries in the following.

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

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

[TABLE]

This is invariant under the rotation around the $z$-axis

[TABLE]

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 phase is called the uniaxial nematic (UN) phase, because the remaining symmetry is $\mathrm{O}(2)\approx\mathrm{U}(1)$. We notice that the rotation in the spin space and the rotation in the real space is 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 as we will explain in the next subsection.

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

For $-1<r<-1/2$, the order parameter $A(\tau,\boldsymbol{x})$ is expressed 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$). Each element is defined by

[TABLE]

We notice that $\boldsymbol{1}_{3}$ is a unit matrix, and $I_{i}$ ($i=1,2,3$) indicate the $\pi$-rotation around the $i$th axis. It is easily checked that $A_{\mathrm{D}_{2}\mathrm{BN}}(\tau,\boldsymbol{x})$ is invariant under the transformation $A_{\mathrm{D}_{2}\mathrm{BN}}(\tau,\boldsymbol{x})\rightarrow OA_{\mathrm{D}_{2}\mathrm{BN}}(\tau,\tilde{\boldsymbol{x}})O^{t}$ ($O\in\mathrm{D}_{2}$). This phase is called the D2-biaxial nematic (D2-BN) phase, because the remaining symmetry is the $\mathrm{D}_{2}$ symmetry.

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

For $r=-1$, the order parameter $A(\tau,\boldsymbol{x})$ is expressed 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 defined by

[TABLE]

which indicates the $\pi/2$-rotation around the 2nd axis ($y$-axis). We notice the phase \bigl{\{}e^{i\alpha}\bigr{\}}=\bigl{\{}1,-1\bigr{\}}\in\mathbb{Z}_{2} ($\alpha=0,\pi$) is locked with the rotation in the spin space and in the real space. It is easily checked that $A_{\mathrm{D}_{4}\mathrm{BN}}(\tau,\boldsymbol{x})$ is invariant under the transformation $A_{\mathrm{D}_{4}\mathrm{BN}}(\tau,\boldsymbol{x})\rightarrow e^{i\alpha}OA_{\mathrm{D}_{4}\mathrm{BN}}(\tau,\tilde{\boldsymbol{x}})O^{t}$ ($(e^{i\alpha},O)\in\mathrm{D}_{4}$). This phase is called the D4-biaxial nematic (D4-BN) phase, because the remaining symmetry is the $\mathrm{D}_{4}$ symmetry.