¹S₀ Pairing Gaps, Chemical Potentials and Entrainment Matrix in Superfluid Neutron-Star Cores for the Brussels–Montreal Functionals

Abstract

Temperature and velocity-dependent ${}^1$S${}_0$ pairing gaps, chemical potentials and entrainment matrix in dense homogeneous neutron–proton superfluid mixtures constituting the outer core of neutron stars, are determined fully self-consistently by solving numerically the time-dependent Hartree–Fock–Bogoliubov equations over the whole range of temperatures and flow velocities for which superfluidity can exist. Calculations have been made for $npe\mu$ in beta-equilibrium using the Brussels–Montreal functional BSk24. The accuracy of various approximations is assessed and the physical meaning of the different velocities and momentum densities appearing in the theory is clarified. Together with the unified equation of state published earlier, the present results provide consistent microscopic inputs for modeling superfluid neutron-star cores.

1. Introduction

Different superfluid and superconducting phases are predicted to exist in neutron stars (see, e.g., [1] for a review). In particular, the (electrically charge neutral) outer core is expected to be made of a neutron–proton superfluid mixture in beta-equilibrium with a normal gas of leptons. More speculative superfluid phases involving other particles such as hyperons or quarks might occur in the innermost region of the core of a neutron star but will not be considered here. Although the interior of the star is highly degenerate, thermal effects may still play an important role in the rotational evolution of the star [2,3,4]. This stems from the fact that the critical temperatures $T_{cq}$ above which superfluidity is destroyed ($q=n,p$ for neutrons and protons respectively) are much smaller than the Fermi temperatures $T_{Fq}$, and may thus be comparable to the actual temperature T of the star. Superfluidity leads to a very complicated dynamics, characterized by the coexistence of different flows (see, e.g., [5] for a recent review). The core of a neutron star involves at least three distinct fluids: the neutron and proton superfluids, as well as a normal fluid made of leptons and excitations. Due to nuclear interactions, the neutron and proton superfluids in the core do not flow freely. They are mutually coupled by entrainment effects of the same kind as the ones discussed by Andreev and Bashkin [6] in the context of superfluid ${}^4$He-${}^3$He mixtures: the mass currents ${\mathit{\rho }}_{\mathit{q}}$ are expressible as linear combinations of the velocity ${\mathit{v}}_{\mathit{N}}$ of the normal fluid and of the so-called “superfluid velocities” ${\mathit{V}}_{\mathit{q}}$ as

$$\begin{array}{c}\hfill {\mathit{\rho }}_{\mathit{n}}=\left({\rho }_n-{\rho }_{nn}-{\rho }_{np}\right){\mathit{v}}_{\mathit{N}}+{\rho }_{nn}{\mathit{V}}_{\mathit{n}}+{\rho }_{np}{\mathit{V}}_{\mathit{p}},\end{array}$$

(1)

$$\begin{array}{c}\hfill {\mathit{\rho }}_{\mathit{p}}=\left({\rho }_p-{\rho }_{pp}-{\rho }_{pn}\right){\mathit{v}}_{\mathit{N}}+{\rho }_{pp}{\mathit{V}}_{\mathit{p}}+{\rho }_{pn}{\mathit{V}}_{\mathit{n}},\end{array}$$

(2)

where ${\rho }_q=m{n}_q$ is the mass density of nucleon of charge type q with associated nucleon number density $n_q$ (ignoring the small difference between neutron and proton masses, which is denoted simply by m). As shown by Carter and Khalatnikov in the context of Landau’s canonical two-fluid model of superfluid ${}^4$He [7] (see also [8]), ${\mathit{V}}_{\mathit{q}}$ are not true velocities but physically represent average nucleon momenta per unit mass. Indeed, it can be easily seen that the true velocities should be defined in terms of the mass currents as ${\mathit{v}}_{\mathit{q}}={\mathit{\rho }}_{\mathit{q}}/{\rho }_q$ and in general these velocities do not coincide with ${\mathit{V}}_{\mathit{q}}$. To avoid any misinterpretation, the superfluid “velocities” are thus written with a capital letter. Please note that ${\mathit{v}}_{\mathit{N}}$ is a true velocity. The importance of entrainment effects is measured by the (symmetric) matrix ${\rho }_{q{q}^{\prime }}$, which is a key microscopic input for modeling the dynamics of neutron stars, see, e.g., [5] and references therein. The entrainment matrix is expected to depend not only on the composition and the baryon density $n=n_n+n_p$, but also on the superfluid velocities as well as the temperature, and this may have an impact on neutron-star oscillations [9,10,11,12]. The influence of temperature and velocity on entrainment effects has been previously studied within Landau’s theory [13,14,15]. However, the Landau parameters, the critical temperatures, and the composition had to be given. Recently, we have derived the entrainment matrix self-consistently for arbitrary superfluid velocities and temperatures within the nuclear-energy-density functional theory by solving exactly the time-dependent Hartree–Fock–Bogoliubov (TDHFB) equations [16,17]. The expressions we have obtained are quite general and applicable to a wide variety of functionals.

In this paper, we have calculated various properties of homogeneous neutron–proton superfluid mixtures in the outer core of neutron stars by solving numerically the self-consistent TDHFB equations using the Brussels–Montreal functional BSk24 [18] for which unified equations of state are already available [19,20,21,22] as well as gravitoelectric and gravitomagnetic tidal Love numbers up to $\ell =5$ [23,24]. More importantly, unlike most available functionals, BSk24 has been accurately calibrated to realistic microscopic calculations of ${}^1$S${}_0$ pairing gaps in infinite homogeneous nuclear matter (at zero temperature and in the absence of currents). We have determined within the same microscopic framework not only the entrainment matrix but also the ${}^1$S${}_0$ pairing gaps and chemical potentials of the superfluid mixture without any approximation, varying the superfluid velocities, the temperature and the baryon density. In this way, we have also been able to assess the accuracy of various approximations. We have focused on ${}^1$S${}_0$ nuclear superfluid phases, which are reliably predicted to exist in the outer core of neutron stars (see, e.g., [25] and references therein). We have not considered ${}^3$PF${}_2$ neutron superfluidity. As a matter of fact, it has been recently shown that in regions where both types of pairing can potentially exist, the ${}^3$PF${}_2$ superfluid phase is completely excluded by the ${}^1$S${}_0$ phase unless strong magnetic fields are present [26].

The formalism to describe nuclear superfluidity is presented in Section 2. After briefly recapitulating the general principles of the TDHFB theory in Section 2.1 and the functionals in Section 2.2, the exact solution in homogeneous nuclear matter is given in Section 2.3, where explicit expressions for various quantities entering the calculations of superfluid properties are derived. The physical interpretation of the different velocities and momentum densities are clarified in Section 2.4. In Section 2.5, it is shown how the TDHFB theory can be recast into Landau’s theory after introducing a series of approximations. Applications to neutron stars are presented in Section 3. The main features of the Brussels–Montreal functionals are recapitulated in Section 3.1. After describing our numerical implementation of the TDHFB equations in Section 3.2, detailed numerical results for various properties of the superfluid mixture are presented and analyzed in Section 3.3, Section 3.4, Section 3.5, Section 3.6, Section 3.7 and Section 3.8. The accuracy of various approximations and interpolations are also discussed. Our conclusions are given in Section 4.

2. Nuclear Superfluidity within the Time-Dependent Hartree–Fock–Bogoliubov Theory

2.1. General Principles

The TDHFB theory [27] provides a unified microscopic framework for studying the dynamics of various nuclear systems, ranging from atomic nuclei to the dense nuclear matter present in neutron stars—the main focus of this work.

Introducing the one-body density matrix $n_q^{ij}=\langle c_q^{j†}{c}_q^i\rangle =n_q^{ji\ast }$ and pairing tensor ${\kappa }_q^{ij}=\langle c_q^j{c}_q^i\rangle =-{\kappa }_q^{ji}$ defined in terms of thermal averages of products of creation ($c_q^{i†}$) and destruction ($c_q^i$) operators for nucleons of charge type q in a quantum state i (using the symbol † for Hermitian conjugation and ∗ for complex conjugation), and assuming that the energy E of a nucleon-matter element of volume V is a function of $n_q^{ij}$, ${\kappa }_q^{ij}$ and ${\kappa }_q^{ij\ast }$ lead to the following equations of motion [27]:

$$\begin{array}{c}\hfill i\hslash \frac{\partial n_q^{ij}}{\partial t}=\sum _k\left(h_q^{ik}{n}_q^{kj}-n_q^{ik}{h}_q^{kj}+{\kappa }_q^{ik}{\Delta }_q^{kj\ast }-{\Delta }_q^{ik}{\kappa }_q^{kj\ast }\right),\end{array}$$

(3)

$$\begin{array}{c}\hfill i\hslash \frac{\partial {\kappa }_q^{ij}}{\partial t}=\sum _k\left[(h_q^{ik}-{\lambda }_q{\delta }^{ik}){\kappa }_q^{kj}+{\kappa }_q^{ik}(h_q^{kj\ast }-{\lambda }_q{\delta }^{kj})-{\Delta }_q^{ik}{n}_q^{kj\ast }-n_q^{ik}{\Delta }_q^{kj}\right]+{\Delta }_q^{ij},\end{array}$$

(4)

where ${\lambda }_q$ denote the chemical potentials. The matrices $h_q^{ij}$ and ${\Delta }_q^{ij}$, defined by

$$\begin{array}{c}\hfill h_q^{ij}=\frac{\partial E}{\partial n_q^{ji}}=h_q^{ji\ast },\end{array}$$

(5)

$$\begin{array}{c}\hfill {\Delta }_q^{ij}=\frac{\partial E}{\partial {\kappa }_q^{ij\ast }}=-{\Delta }_q^{ji},\end{array}$$

(6)

generally depend themselves on $n_q^{ij}$, ${\kappa }_q^{ij}$ and ${\kappa }_q^{ij\ast }$ and must therefore be determined self-consistently.

2.2. Functionals of Local Densities and Currents

The class of energy functionals $E(n_q^{ij},{\kappa }_q^{ij},{\kappa }_q^{ij\ast })$ that we consider here depend on the density matrices and pair tensors only through the following local densities and currents:

(i)

the nucleon densities at position $\mathit{r}$ at time t ($\sigma =\pm 1$ distinguishing the two spin states),

$$\begin{array}{c}\hfill n_q(\mathit{r},t)=\sum _{\sigma =\pm 1}{n}_q(\mathit{r},\sigma ;\mathit{r},\sigma ;t),\end{array}$$

(7)

(ii)

the kinetic-energy densities (in units of ${\hslash }^2/2m)$ at position $\mathit{r}$ and time t,

$$\begin{array}{c}\hfill {\tau }_q(\mathit{r},t)=\sum _{\sigma =\pm 1}\int {\mathrm{d}}^3{\mathit{r}}^{\mathrm{\prime }}\delta (\mathit{r}-{\mathit{r}}^{\mathbf{\prime }})\mathsf{\nabla }·{\mathsf{\nabla }}^{\mathbf{\prime }}{n}_q(\mathit{r},\sigma ;{\mathit{r}}^{\mathbf{\prime }},\sigma ;t),\end{array}$$

(8)

(iii)

the momentum densities (in units of ℏ) at position $\mathit{r}$ and time t,

$$\begin{array}{c}\hfill {\mathit{j}}_{\mathit{q}}(\mathit{r},t)=-\frac{i}{2}\sum _{\sigma =\pm 1}\int {\mathrm{d}}^3{\mathit{r}}^{\mathbf{\prime }}\delta (\mathit{r}-{\mathit{r}}^{\mathbf{\prime }})(\mathsf{\nabla }-{\mathsf{\nabla }}^{\mathbf{\prime }})n_q(\mathit{r},\sigma ;{\mathit{r}}^{\mathbf{\prime }},\sigma ;t),\end{array}$$

(9)

(iv)

the abnormal densities at position $\mathit{r}$ and time t

$$\begin{array}{c}\hfill {\tilde{n}}_q(\mathit{r},t)=\sum _{\sigma =\pm 1}{\tilde{n}}_q(\mathit{r},\sigma ;\mathit{r},\sigma ;t),\end{array}$$

(10)

where

$$n_q(\mathit{r},\sigma ;{\mathit{r}}^{\prime },{\sigma }^{\prime };t)=\sum _{i,j}{n}_q^{ij}{\phi }_i^{\left(q\right)}(\mathit{r},\sigma ){\phi }_j^{\left(q\right)}{({\mathit{r}}^{\prime },{\sigma }^{\prime })}^{\ast },$$

(11)

$${\tilde{n}}_q(\mathit{r},\sigma ;{\mathit{r}}^{\prime },{\sigma }^{\prime };t)=-{\sigma }^{\prime }\sum _{i,j}{\kappa }_q^{ij}{\phi }_i^{\left(q\right)}(\mathit{r},\sigma ){\phi }_j^{\left(q\right)}({\mathit{r}}^{\prime },-{\sigma }^{\prime }),$$

(12)

are the so-called normal and abnormal density matrices respectively [28], and ${\phi }_i^{\left(q\right)}(\mathit{r},\sigma )$ is the single-particle wavefunction associated with the state i. The abnormal densities are the local order parameters of the neutron and proton superfluid phases [17]. These densities are complex and the gradient of their respective phase ${\varphi }_q(\mathit{r},t)$ defines the superfluid velocities as follows:

$$\begin{array}{c}\hfill {\mathit{V}}_{\mathit{q}}(\mathit{r},t)=\frac{\hslash }{2m}\mathsf{\nabla }{\varphi }_q(\mathit{r},t).\end{array}$$

(13)

The matrices (5) and (6) can be alternatively expressed as [17]

$$\begin{array}{c}\hfill h_q^{ij}\left(t\right)=\sum _{\sigma }\int {\mathrm{d}}^3\mathit{r}{\phi }_i^{\left(q\right)}{(\mathit{r},\sigma )}^{\ast }{h}_q(\mathit{r},t){\phi }_j^{\left(q\right)}(\mathit{r},\sigma ),\end{array}$$

(14)

$$\begin{array}{c}\hfill {\Delta }_q^{ij}\left(t\right)=-\sum _{\sigma }\sigma \int {\mathrm{d}}^3\mathit{r}{\phi }_i^{\left(q\right)}{(\mathit{r},\sigma )}^{\ast }{\Delta }_q(\mathit{r},t){\phi }_j^{\left(q\right)}{(\mathit{r},-\sigma )}^{\ast },\end{array}$$

(15)

where

$$\begin{array}{ccc}\hfill h_q(\mathit{r},t)& =& -\mathsf{\nabla }·\frac{{\hslash }^2}{2m_q^{\oplus }(\mathit{r},t)}\mathsf{\nabla }+U_q(\mathit{r},t)-\frac{i}{2}\left[{\mathit{I}}_{\mathit{q}}(\mathit{r},t)·\mathsf{\nabla }+\mathsf{\nabla }·{\mathit{I}}_{\mathit{q}}(\mathit{r},t)\right],\hfill \end{array}$$

(16)

$$\begin{array}{c}\hfill \frac{{\hslash }^2}{2m_q^{\oplus }(\mathit{r},t)}=\frac{\delta E}{\delta {\tau }_q(\mathit{r},t)},U_q(\mathit{r},t)=\frac{\delta E}{\delta n_q(\mathit{r},t)},{\mathit{I}}_{\mathit{q}}(\mathit{r},t)=\frac{\delta E}{\delta {\mathit{j}}_{\mathit{q}}(\mathit{r},t)},\end{array}$$

(17)

$$\begin{array}{c}\hfill {\Delta }_q(\mathit{r},t)=2\frac{\delta E}{\delta {\tilde{n}}_q{(\mathit{r},t)}^{\ast }}=2\frac{\delta E}{\delta |{\tilde{n}}_q{(\mathit{r},t)|}^2}{\tilde{n}}_q(\mathit{r},t).\end{array}$$

(18)

The last equality in Equation (18) arises from the requirement that the energy must be real. Please note that the pair potential ${\Delta }_q(\mathit{r},t)$ is a complex field sharing the same phase ${\varphi }_q(\mathit{r},t)$ as the abnormal density ${\tilde{n}}_q(\mathit{r},t)$. Let us recall that this phase enters through the definition of the superfluid velocities (13).

2.3. Application to Homogeneous Systems

Considering a homogeneous neutron–proton superfluid mixture with stationary flows in the normal-fluid rest frame where ${\mathit{v}}_{\mathit{N}}=\mathit{0}$, the TDHFB equations can be solved exactly [17]. In particular, the entrainment matrix reads (${\delta }_{q{q}^{\prime }}$ denotes the Kronecker symbol and $\rho ={\rho }_n+{\rho }_p$ is the total mass density)

$$\begin{array}{c}\hfill {\rho }_{q{q}^{\prime }}={\rho }_q\left(1-{\mathcal{Y}}_q\right)\left(\frac{m}{m_q^{\oplus }}{\delta }_{q{q}^{\prime }}+\frac{{\mathcal{I}}_{q{q}^{\prime }}}{\hslash }\right),\end{array}$$

(19)

with

$$\begin{array}{c}\hfill \frac{m}{m_q^{\oplus }}=1+\frac{2\rho }{{\hslash }^2}\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\right)+\frac{4{\rho }_q}{{\hslash }^2}\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1},\end{array}$$

(20)

$$\begin{array}{cc}\hfill {\mathcal{I}}_{nn}& =\frac{2}{\hslash }{\rho }_n\left(1-{\mathcal{Y}}_n\right)\mathsf{\Theta }\left[\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\left(\frac{8}{{\hslash }^2}\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}{m}_p^{\oplus }{n}_p{\mathcal{Y}}_p-1\right)-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}\right],\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\mathcal{I}}_{pp}& =\frac{2}{\hslash }{\rho }_p\left(1-{\mathcal{Y}}_p\right)\mathsf{\Theta }\left[\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\left(\frac{8}{{\hslash }^2}\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}{m}_n^{\oplus }{n}_n{\mathcal{Y}}_n-1\right)-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}\right],\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\mathcal{I}}_{np}& =\frac{2}{\hslash }{\rho }_p\left(1-{\mathcal{Y}}_p\right)\mathsf{\Theta }\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\mathcal{I}}_{pn}& =\frac{2}{\hslash }{\rho }_n\left(1-{\mathcal{Y}}_n\right)\mathsf{\Theta }\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \mathsf{\Theta }& \equiv \left[1-\frac{2}{{\hslash }^2}\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}+\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\right)\left(m_n^{\oplus }{n}_n{\mathcal{Y}}_n+m_p^{\oplus }{n}_p{\mathcal{Y}}_p\right)\right.\hfill \\ \hfill & {\left.+{\left(\frac{4}{{\hslash }^2}\right)}^2\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}{m}_n^{\oplus }{n}_n{m}_p^{\oplus }{n}_p{\mathcal{Y}}_n{\mathcal{Y}}_p\right]}^{-1}.\hfill \end{array}$$

(25)

Here $E_{\mathrm{nuc}}^j$ represents the nuclear-energy terms contributing to the mass currents. Galilean invariance requires that these terms depend on the following combinations:

$$X_0(\mathit{r},t)=n_0(\mathit{r},t){\tau }_0(\mathit{r},t)-{\mathit{j}}_{\mathit{0}}{(\mathit{r},t)}^2,$$

(26)

$$X_1(\mathit{r},t)=n_1(\mathit{r},t){\tau }_1(\mathit{r},t)-{\mathit{j}}_{\mathit{1}}{(\mathit{r},t)}^2.$$

(27)

The subscripts 0 and 1 denote isoscalar and isovector quantities, respectively, namely sums over neutrons and protons for the former (e.g., $n_0\equiv n=n_n+n_p$) and differences between neutrons and protons for the latter (e.g., $n_1=n_n-n_p$). The temperature and velocity-dependent functions ${\mathcal{Y}}_q$ are defined by ($k_{\mathrm{B}}$ being the Boltzmann constant)

$$\begin{array}{c}\hfill {\mathcal{Y}}_q(T,{\mathsf{𝕍}}_{\mathit{q}})\equiv \frac{\hslash }{m_q^{\oplus }{n}_q{\mathsf{𝕍}}_q^2}\frac{1}{V}\sum _{\mathit{k}}\mathit{k}·{\mathsf{𝕍}}_{\mathit{q}}tanh\left(\frac{{\mathfrak{E}}_{\mathit{k}}^{\left(q\right)}}{2k_{\mathrm{B}}T}\right),\end{array}$$

(28)

where we have introduced the effective superfluid velocities

$$\begin{array}{c}\hfill {\mathsf{𝕍}}_{\mathit{q}}\equiv \frac{m}{m_q^{\oplus }}{\mathit{V}}_{\mathit{q}}+\frac{{\mathit{I}}_{\mathit{q}}}{\hslash },\end{array}$$

(29)

and ${\mathfrak{E}}_{\mathit{k}}^{\left(q\right)}$ represent the energies of quasiparticle excitations, given by

$$\begin{array}{c}\hfill {\mathfrak{E}}_{\mathit{k}}^{\left(q\right)}=\hslash \mathit{k}·{\mathsf{𝕍}}_{\mathit{q}}+\sqrt{{\epsilon }_{\mathit{k}}^{\left(q\right)2}+{\Delta }_q^2},\end{array}$$

(30)

with

$$\begin{array}{c}\hfill {\epsilon }_{\mathit{k}}^{\left(q\right)}=\frac{{\hslash }^2{\mathit{k}}^2}{2m_q^{\oplus }}+\frac{1}{2}{m}_q^{\oplus }\left({\mathsf{𝕍}}_{\mathit{q}}+\frac{{\mathit{I}}_{\mathit{q}}}{\hslash }\right)·\left({\mathsf{𝕍}}_{\mathit{q}}-\frac{{\mathit{I}}_{\mathit{q}}}{\hslash }\right)+U_q-{\lambda }_q.\end{array}$$

(31)

In turn, the vector ${\mathit{I}}_{\mathit{q}}$ is expressible in terms of the superfluid velocities as follows

$$\begin{array}{c}\hfill {\mathit{I}}_{\mathit{q}}=\sum _{q^{\prime }}{\mathcal{I}}_{q{q}^{\prime }}{\mathit{V}}_{{\mathit{q}}^{\prime }}.\end{array}$$

(32)

The pairing gaps (as defined as the nonvanishing matrix elements of the pair potential, see [17] ) are obtained from the self-consistent equations

$$\begin{array}{c}\hfill {\Delta }_q(T,{\mathsf{𝕍}}_{\mathit{q}})=-\frac{2}{V}\frac{\delta E}{\delta |{\tilde{n}}_q{|}^2}\sum _{\mathit{k}}\frac{{\Delta }_q(T,{\mathsf{𝕍}}_{\mathit{q}})}{\sqrt{{\epsilon }_{\mathit{k}}^{\left(q\right)2}+{\Delta }_q{(T,{\mathsf{𝕍}}_{\mathit{q}})}^2}}tanh\left(\frac{{\mathfrak{E}}_{\mathit{k}}^{\left(q\right)}}{2k_{\mathrm{B}}T}\right),\end{array}$$

(33)

where it is understood that the summation must be regularized to remove ultraviolet divergences, as will be discussed below. The gap equations must be solved together with the particle number conservation conditions

$$\begin{array}{c}\hfill n_q=\frac{1}{V}\sum _{\mathit{k}}\left[1-\frac{{\epsilon }_{\mathit{k}}^{\left(q\right)}}{\sqrt{{\epsilon }_{\mathit{k}}^{\left(q\right)2}+{\Delta }_q^2}}tanh\left(\frac{{\mathfrak{E}}_{\mathit{k}}^{\left(q\right)}}{2k_{\mathrm{B}}T}\right)\right].\end{array}$$

(34)

As can be seen from Equation (31), Equations (28), (33) and (34) all depend on the reduced chemical potentials defined by

$$\begin{array}{c}\hfill {\mu }_q={\lambda }_q-U_q-\frac{1}{2}{m}_q^{\oplus }\left({\mathsf{𝕍}}_{\mathit{q}}+\frac{{\mathit{I}}_{\mathit{q}}}{\hslash }\right)·\left({\mathsf{𝕍}}_{\mathit{q}}-\frac{{\mathit{I}}_{\mathit{q}}}{\hslash }\right)\end{array}$$

(35)

so that neither the pairing gaps nor the entrainment matrix require the explicit form of the potentials $U_q$.

From now on, we will take the continuum limit, i.e., we will replace discrete summations over wave vectors $\mathit{k}$ by integrations as follows:

$$\begin{array}{c}\hfill \frac{1}{V}\sum _{\mathit{k}}\cdots \to \int \frac{{\mathrm{d}}^3\mathit{k}}{{\left(2\pi \right)}^3}\cdots =\int \frac{\mathrm{d}{\mathsf{\Omega }}_{\mathit{k}}}{4\pi }{\int }_{-{\mu }_q}^{+\infty }\mathrm{d}\epsilon \mathcal{D}\left(\epsilon \right)\cdots \end{array}$$

(36)

with ${\mathsf{\Omega }}_{\mathit{k}}$ the solid angle in $\mathit{k}$-space and $\mathcal{D}\left(\epsilon \right)$ the density of single-particle states per one spin state given by

$$\begin{array}{c}\hfill \mathcal{D}\left(\epsilon \right)=\frac{m_q^{\oplus }}{2{\pi }^2{\hslash }^3}\sqrt{2m_q^{\oplus }(\epsilon +{\mu }_q)}.\end{array}$$

(37)

Integrating over solid angle and changing variables, Equation (28) can thus be expressed as

$$\begin{array}{ccc}\hfill {\mathcal{Y}}_q& =\hfill & \hfill \frac{3}{8}\frac{{\overline{T}}_q}{{\overline{\mathbb{V}}}_q^2}{\int }_0^{+\infty }\mathrm{d}x\sqrt{x}log\left\{\left[1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}\right]\left[1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}\right]\right\}\\ \hfill & \hfill & \hfill +\frac{3}{16}\frac{{\overline{T}}_q^2}{{\overline{\mathbb{V}}}_q^3}{\int }_0^{+\infty }\mathrm{d}x\left\{{\mathrm{Li}}_2\left[-{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}\right]-{\mathrm{Li}}_2\left[-{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}\right]\right\}\end{array}$$

(38)

where

$$\begin{array}{c}\hfill {\mathbb{E}}_x^{\left(q\right)}=\sqrt{{\left(x-{\overline{\mu }}_q\right)}^2+{\overline{\Delta }}_q^2},\end{array}$$

(39)

${\displaystyle {\mathrm{Li}}_2\left(x\right)={\int }_1^{1-x}\frac{logu}{1-u}\mathrm{d}u}$ is the dilogarithm function, and we have introduced the dimensionless ratios

$$\begin{array}{c}\hfill {\overline{T}}_q\equiv \frac{T}{T_{Fq}},{\overline{\mathbb{V}}}_q\equiv \frac{{\mathbb{V}}_q}{V_{Fq}},{\overline{\mu }}_q\equiv \frac{{\mu }_q}{{\epsilon }_{Fq}},{\overline{\Delta }}_q\equiv \frac{{\Delta }_q}{{\epsilon }_{Fq}}.\end{array}$$

(40)

The Fermi temperature is defined by $T_{Fq}={\epsilon }_{Fq}/k_{\mathrm{B}}$ with the Fermi energy

$${\epsilon }_{Fq}=\frac{{\hslash }^2{k}_{Fq}^2}{2m_q^{\oplus }}$$

(41)

and Fermi wave number $k_{Fq}={\left(3{\pi }^2{n}_q\right)}^{1/3}$; the Fermi velocity is given by

$$V_{Fq}=\frac{\hslash k_{Fq}}{m_q^{\oplus }}.$$

(42)

Similarly, the gap Equation (33) and the particle number conservation Equation (34) become, respectively

$$\begin{array}{ccc}\hfill {\Delta }_q& =& {\displaystyle -\frac{m_q^{\oplus }{k}_{Fq}}{2{\pi }^2{\hslash }^2}\frac{{\overline{T}}_q}{{\overline{\mathbb{V}}}_q}\frac{\delta E}{\delta |{\tilde{n}}_q{|}^2}{\Delta }_q{\int }_0^{({\mu }_q+{\epsilon }_{\mathsf{\Lambda }})/{\epsilon }_{Fq}}\frac{\mathrm{d}x}{{\mathbb{E}}_x^{\left(q\right)}}}\hfill \\ & & \times log\left[cosh\left(\frac{{\mathbb{E}}_x^{\left(q\right)}}{2{\overline{T}}_q}+\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\sqrt{x}\right)\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left(\frac{{\mathbb{E}}_x^{\left(q\right)}}{2{\overline{T}}_q}-\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\sqrt{x}\right)\right],\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill \frac{4}{3}& =& {\int }_0^{+\infty }\mathrm{d}x\left\{\sqrt{x}-\frac{{\overline{T}}_q}{{\overline{\mathbb{V}}}_q}\frac{x-{\overline{\mu }}_q}{2{\mathbb{E}}_x^{\left(q\right)}}\right.\hfill \\ & & \left.\times log\left[cosh\left(\frac{{\mathbb{E}}_x^{\left(q\right)}}{2{\overline{T}}_q}+\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\sqrt{x}\right)\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left(\frac{{\mathbb{E}}_x^{\left(q\right)}}{2{\overline{T}}_q}-\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\sqrt{x}\right)\right]\right\},\hfill \end{array}$$

(44)

and ${\epsilon }_{\mathsf{\Lambda }}$ is a cutoff above the Fermi level to regularize the ultraviolet divergences (see Section 3.1). Expressing the hyperbolic functions in terms of the exponential function, we can alternatively rewrite (43) and (44) as

$$\begin{array}{ccc}\hfill {\Delta }_q& =& {\displaystyle -\frac{m_q^{\oplus }{k}_{Fq}}{{\pi }^2{\hslash }^2}\frac{\delta E}{\delta |{\tilde{n}}_q{|}^2}{\Delta }_q{\int }_0^{({\mu }_q+{\epsilon }_{\mathsf{\Lambda }})/{\epsilon }_{Fq}}\mathrm{d}x\frac{\sqrt{x}}{{\mathbb{E}}_x^{\left(q\right)}}}\hfill \\ & & \times \left\{1+\frac{{\overline{T}}_q}{2{\overline{\mathbb{V}}}_q\sqrt{x}}log\left[\frac{1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}}{1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}}\right]\right\},\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill \frac{4}{3}& =& {\int }_0^{+\infty }\mathrm{d}x\sqrt{x}\left\{1-\frac{x-{\overline{\mu }}_q}{{\mathbb{E}}_x^{\left(q\right)}}\right.\hfill \\ & & \left.\times \left[1+\frac{{\overline{T}}_q}{2{\overline{\mathbb{V}}}_q\sqrt{x}}log\left(\frac{1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}}{1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}}\right)\right]\right\}.\hfill \end{array}$$

(46)

We have made use of the identity $log\left[cosh\left(a+b\right)\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left(a-b\right)\right]=2b+log\left[1+e^{-2\left(a+b\right)}\right]-log\left[1+e^{-2(a-b)}\right]$.

It is worth remarking that although the pairing gaps and the entrainment matrix depend in general on the directions of the superfluid velocities ${\mathit{V}}_{\mathit{q}}$, this dependence is entirely contained in the norm of the effective superfluid velocities ${\mathsf{𝕍}}_{\mathit{q}}$. Using Equations (29) and (32), it can be seen that the two kinds of velocities are related by

$$\begin{array}{c}\hfill {\mathsf{𝕍}}_{\mathit{q}}=\sum _{q^{\prime }=n,p}\left(\frac{m}{m_q^{\oplus }}{\delta }_{q{q}^{\prime }}+\frac{{\mathcal{I}}_{q{q}^{\prime }}}{\hslash }\right){\mathit{V}}_{{\mathit{q}}^{\prime }}.\end{array}$$

(47)

It is important to realize that this relation is highly non-linear because the matrix elements ${\mathcal{I}}_{q{q}^{\prime }}$, defined by Equations (21)–(25), depend themselves on the effective superfluid velocities through the functions ${\mathcal{Y}}_q$. For this reason, the mapping between ${\mathsf{𝕍}}_{\mathit{n}}$, ${\mathsf{𝕍}}_{\mathit{p}}$ and ${\mathit{V}}_{\mathit{n}},{\mathit{V}}_{\mathit{p}}$ is quite complicated. It is, therefore, much more convenient to express the results in terms of ${\mathsf{𝕍}}_{\mathit{q}}$ instead of ${\mathit{V}}_{\mathit{q}}$. In particular, it can be seen that the neutron (proton) pairing gaps depend only the norms of neutron (proton) effective superfluid velocity. It is only when the chemical potentials ${\lambda }_q$ are needed rather than the reduced ones ${\mu }_q$ that the directions of the superflows become important since ${\lambda }_q$ are obtained from Equation (35) using Equation (32), namely

$$\begin{array}{c}\hfill {\lambda }_n={\mu }_n+\left(\frac{1}{2}\frac{m}{m_n^{\oplus }}+\frac{{\mathcal{I}}_{nn}}{\hslash }\right)m{\mathit{V}}_{\mathit{n}}^2+\frac{{\mathcal{I}}_{np}}{\hslash }m{\mathit{V}}_{\mathit{n}}·{\mathit{V}}_{\mathit{p}}+U_n,\end{array}$$

(48)

$$\begin{array}{c}\hfill {\lambda }_p={\mu }_p+\left(\frac{1}{2}\frac{m}{m_p^{\oplus }}+\frac{{\mathcal{I}}_{pp}}{\hslash }\right)m{\mathit{V}}_{\mathit{p}}^2+\frac{{\mathcal{I}}_{pn}}{\hslash }m{\mathit{V}}_{\mathit{p}}·{\mathit{V}}_{\mathit{n}}+U_p.\end{array}$$

(49)

The potentials $U_q$ are functions of the nucleon densities $n_q$, the momentum densities ${\mathit{j}}_{\mathit{q}}$ and the kinetic densities ${\tau }_q$. The momentum density ${\mathit{j}}_{\mathit{q}}$ can be expressed as [17]

$$\begin{array}{c}\hfill {\mathit{j}}_{\mathit{q}}=\frac{{\rho }_q}{\hslash }\left[(1-{\mathcal{Y}}_q){\mathit{V}}_{\mathit{q}}-\frac{m_q^{\oplus }}{m}{\mathcal{Y}}_q\sum _{q^{\prime }=n,p}\frac{{\mathcal{I}}_{q{q}^{\prime }}}{\hslash }{\mathit{V}}_{{\mathit{q}}^{\prime }}\right].\end{array}$$

(50)

Using Equation (8), we find for the kinetic-energy density:

$$\begin{array}{cc}\hfill {\tau }_q=& \frac{3}{4}{\left(3{\pi }^2\right)}^{2/3}{n}_q^{5/3}{\int }_0^{({\mu }_q+{\epsilon }_{\mathsf{\Lambda }})/{\epsilon }_{Fq}}\mathrm{d}xx\left\{\sqrt{x}-\frac{{\overline{T}}_q}{{\overline{\mathbb{V}}}_q}\frac{x-{\overline{\mu }}_q}{2{\mathbb{E}}_x^{\left(q\right)}}\right.\hfill \\ \hfill & \left.\times log\left[cosh\left(\frac{{\mathbb{E}}_x^{\left(q\right)}}{2{\overline{T}}_q}+\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\sqrt{x}\right)\left(\frac{{\mathbb{E}}_x^{\left(q\right)}}{2{\overline{T}}_q}-\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\sqrt{x}\right)\right]\right\}\hfill \\ \hfill & +\frac{1}{2}\left(\frac{2m_q^{\oplus }}{{\hslash }^2}\right){\rho }_q{\mathit{V}}_{\mathit{q}}·\left({\mathit{V}}_{\mathit{q}}-2{\mathcal{Y}}_q{\mathsf{𝕍}}_{\mathit{q}}\right).\hfill \end{array}$$

(51)

In the regime ${\overline{T}}_q\ll 1$, ${\overline{\mathbb{V}}}_q\ll 1$, ${\overline{\Delta }}_q\ll 1$ and ${\overline{\mu }}_q\approx 1$, the second term in the right-hand side of Equation (51) becomes negligible and the integral reduces to the Thomas-Fermi expression ${\tau }_q\approx \frac{3}{5}{\left(3{\pi }^2\right)}^{2/3}{n}_q^{5/3}$. The kinetic-energy density can be equivalently expressed in terms of the exponential function as

$$\begin{array}{cc}\hfill {\tau }_q=& \frac{3}{4}{\left(3{\pi }^2\right)}^{2/3}{n}_q^{5/3}{\int }_0^{({\mu }_q+{\epsilon }_{\mathsf{\Lambda }})/{\epsilon }_{Fq}}\mathrm{d}x{x}^{3/2}\left\{1-\frac{x-{\overline{\mu }}_q}{{\mathbb{E}}_x^{\left(q\right)}}\right.\hfill \\ \hfill & \left.\times \left[1+\frac{{\overline{T}}_q}{2{\overline{\mathbb{V}}}_q\sqrt{x}}log\left(\frac{1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}}{1+{\mathrm{e}}^{-\left({\mathbb{E}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\sqrt{x}\right)/{\overline{T}}_q}}\right)\right]\right\}\hfill \\ \hfill & +\frac{1}{2}\left(\frac{2m_q^{\oplus }}{{\hslash }^2}\right){\rho }_q{\mathit{V}}_{\mathit{q}}·\left({\mathit{V}}_{\mathit{q}}-2{\mathcal{Y}}_q{\mathsf{𝕍}}_{\mathit{q}}\right).\hfill \end{array}$$

(52)

2.4. Physical Interpretation of the Different Velocities and Momentum Densities

Using Equation (66) of [17], it can be immediately seen that the true velocities associated with the transport of nucleons (mass) are related to the effective superfluid velocities (29) through the relation

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathit{q}}\equiv \frac{{\mathsf{\rho }}_{\mathit{q}}}{{\rho }_q}=(1-{\mathcal{Y}}_q){\mathsf{𝕍}}_{\mathit{q}}.\end{array}$$

(53)

Let us recall that these velocities are measured relative to the normal-fluid rest frame. At zero temperature and subcritical superflow of nucleons of type q, the functions ${\mathcal{Y}}_q$ will be shown to vanish in Section 3.5: in this case, the effective superfluid velocity thus actually represents the true velocity ${\mathit{v}}_{\mathit{q}}={\mathsf{𝕍}}_{\mathit{q}}$. At finite temperatures, the excitation of quasiparticles entails a finite fraction ${\mathcal{Y}}_q>0$: nucleons thus move with a lower speed at $T>0$ than at $T=0$. If nucleons of type q are nonsuperfluid, ${\mathcal{Y}}_q=1$ as we will see in Section 3.5, therefore their true velocity vanishes $v_q=0$: nucleons move with the normal fluid (however, the other nucleon species can flow with a different velocity if it is superfluid). The function ${\mathcal{Y}}_q$ thus measures the relative importance of quasiparticle excitations for the transport of nucleons of type q.

As already mentioned earlier, the superfluid “velocity” ${\mathit{V}}_{\mathit{q}}$ defined by the gradient of the phase of the condensate through Equation (13) represents the momentum per unit mass of the superfluid. The superfluid momentum density of the nucleon species q, given by ${\rho }_q{\mathit{V}}_{\mathit{q}}$, does not coincide with the momentum density $\hslash {\mathit{j}}_{\mathit{q}}$ introduced in Equation (54). This stems from the fact that the latter not only accounts for the superfluid momentum density but also includes the contribution from quasiparticles. This can be directly seen from Equation (50), which can be equivalently written as

$$\begin{array}{c}\hfill \hslash {\mathit{j}}_{\mathit{q}}={\rho }_q{\mathit{V}}_{\mathit{q}}-{\mathcal{Y}}_q{\rho }_q\frac{m_q^{\oplus }}{m}{\mathsf{𝕍}}_{\mathit{q}}.\end{array}$$

(54)

The second term can be interpreted as the momentum density of quasiparticles. Indeed, this contribution vanishes if ${\mathcal{Y}}_q=0$, i.e., in the absence of quasiparticle excitations. It is only in this limiting case that the total momentum density $\hslash {\mathit{j}}_{\mathit{q}}$ coincides with the superfluid momentum density ${\rho }_q{\mathit{V}}_{\mathit{q}}$. In general, it can be shown using the self-consistent solutions of the TDHFB equations presented in the previous section that the total mass current is equal to the total momentum density

$$\begin{array}{c}\hfill {\mathit{\rho }}_{\mathit{n}}+{\mathit{\rho }}_{\mathit{p}}=\hslash ({\mathit{j}}_{\mathit{n}}+{\mathit{j}}_{\mathit{p}}),\end{array}$$

(55)

as required by Galilean invariance (this identity can be more easily demonstrated using the general expression of the mass currents [16]).

The distinction between the different velocities and momentum densities becomes irrelevant if both nucleon species are nonsuperfluid since ${\mathcal{Y}}_n={\mathcal{Y}}_p=1$ implies that ${\mathit{v}}_{\mathit{n}}$, ${\mathit{v}}_{\mathit{p}}$, $\hslash {\mathit{j}}_{\mathit{n}}/{\rho }_n$ and $\hslash {\mathit{j}}_{\mathit{p}}/{\rho }_p$ all vanish in the fluid rest frame, i.e., all nucleons move with the normal fluid, as expected. Likewise, in the limiting case of a single superfluid constituent at zero temperature and subcritical superflow, we have ${\mathit{v}}_{\mathit{q}}={\mathit{V}}_{\mathit{q}}={\mathsf{𝕍}}_{\mathit{q}}=\hslash {\mathit{j}}_{\mathit{q}}/{\rho }_q$.

2.5. Landau’s Approximations

The neutron–proton superfluid mixture can be alternatively described using Landau’s theory [13,14,15]. The TDHFB theory can be reduced to a similar form after introducing a series of approximations. Specifically, assuming that the critical temperatures and the critical superfluid velocities are small compared to their Fermi counterpart,

• instead of solving Equation (44), the reduced chemical potentials (35) are approximated by their associated Fermi energies (${\mu }_q\approx {\epsilon }_{Fq}$), thus ignoring any dependence on temperature, currents, and pairing gaps;
• the single-particle energies (31) are calculated at zero temperature, in the absence of currents ignoring any dependence on the pairing gaps, and expanding linearly around the Fermi surface (denoting by $\breve{Q}$ the approximate expression for a quantity Q)

  $$\begin{array}{c}\hfill {\epsilon }_{\mathit{k}}^{\left(q\right)}\approx {\breve{\epsilon }}_{\mathit{k}}^{\left(q\right)}\equiv \hslash V_{Fq}(k-k_{Fq});\end{array}$$

  (56)
• the quasiparticle energies (30) are similarly expanded as

  $$\begin{array}{c}\hfill {\mathfrak{E}}_{\mathit{k}}^{\left(q\right)}\approx {\breve{\mathfrak{E}}}_{\mathit{k}}^{\left(q\right)}+\hslash k_{Fq}{\mathsf{𝕍}}_{\mathit{q}}cos{\theta }_{\mathit{k}},{\breve{\mathfrak{E}}}_{\mathit{k}}^{\left(q\right)}=\sqrt{{\breve{\epsilon }}_{\mathit{k}}^{\left(q\right)2}+{\breve{\Delta }}_q^2};\end{array}$$

  (57)
• the density of single-particle states $\mathcal{D}\left(\epsilon \right)$ in $\mathit{k}$-space integrations (36) is approximated by its value on the Fermi surface, $\mathcal{D}\left(\epsilon \right)\approx \breve{\mathcal{D}}\left(0\right)$ with

  $$\begin{array}{c}\hfill \breve{\mathcal{D}}\left(0\right)=\frac{k_{Fq}{\breve{m}}_q^{\oplus }}{2{\pi }^2{\hslash }^2};\end{array}$$

  (58)
• the derivatives of the nuclear-energy terms $E_{\mathrm{nuc}}^j$ entering Equations (20)–(25), are evaluated in the absence of currents;

In previous studies [13,14,15], the pairing gaps ${\breve{\Delta }}_q$ were obtained in the weak-coupling approximation ${\breve{\Delta }}_q\ll {\epsilon }_{Fq},{\epsilon }_{\mathsf{\Lambda }}$ at zeroth order from the following approximate equation (see Appendix A)

$$\begin{array}{c}\hfill log\left(\frac{{\breve{\Delta }}_q}{{\breve{\Delta }}_q^{\left(0\right)}}\right)\approx {\int }_0^{+\infty }\frac{\mathrm{d}x}{{\breve{\mathbb{E}}}_x^{\left(q\right)}}\left\{\frac{{\overline{T}}_q}{2{\overline{\mathbb{V}}}_q}log\left[cosh\left(\frac{{\breve{\mathbb{E}}}_x^{\left(q\right)}}{2{\overline{T}}_q}+\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\right)\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left(\frac{{\breve{\mathbb{E}}}_x^{\left(q\right)}}{2{\overline{T}}_q}-\frac{{\overline{\mathbb{V}}}_q}{{\overline{T}}_q}\right)\right]-1\right\},\end{array}$$

(59)

or in terms of the exponential function

$$\begin{array}{c}\hfill log\left(\frac{{\breve{\Delta }}_q}{{\breve{\Delta }}_q^{\left(0\right)}}\right)\approx \frac{{\overline{T}}_q}{2{\overline{\mathbb{V}}}_q}{\int }_0^{+\infty }\frac{\mathrm{d}x}{{\breve{\mathbb{E}}}_x^{\left(q\right)}}log\left[\frac{1+{\mathrm{e}}^{-\left({\breve{\mathbb{E}}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\right)/{\overline{T}}_q}}{1+{\mathrm{e}}^{-\left({\breve{\mathbb{E}}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\right)/{\overline{T}}_q}}\right],\end{array}$$

(60)

where ${\breve{\Delta }}_q^{\left(0\right)}$ denotes the pairing gaps at $T=0$ in the absence of currents. The latter were determined using the BCS relation [29] (introducing the Euler–Mascheroni constant $\gamma \simeq 0.577216$)

$$\begin{array}{c}\hfill {\breve{\Delta }}_q^{\left(0\right)}=\frac{k_{\mathrm{B}}\pi }{{\mathrm{e}}^{\gamma }}{\breve{T}}_{cq}^{\left(0\right)},\end{array}$$

(61)

by fixing arbitrarily the associated critical temperatures ${\breve{T}}_{cq}^{\left(0\right)}$.

Moreover, the functions ${\mathcal{Y}}_q$ were replaced by the functions ${\mathsf{\Phi }}_q$ of [15], which can be expressed as

$$\begin{array}{ccc}\hfill {\mathsf{\Phi }}_q& =& \frac{3}{4}\frac{{\overline{T}}_q}{{\overline{\mathbb{V}}}_q^2}{\int }_0^{+\infty }\mathrm{d}xlog\left\{\left[1+{\mathrm{e}}^{-\left({\breve{\mathbb{E}}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\right)/{\overline{T}}_q}\right]\left[1+{\mathrm{e}}^{-\left({\breve{\mathbb{E}}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\right)/{\overline{T}}_q}\right]\right\}\hfill \\ & +& \frac{3}{8}\frac{{\overline{T}}_q^2}{{\overline{\mathbb{V}}}_q^3}{\int }_0^{+\infty }\mathrm{d}x\left\{{\mathrm{Li}}_2\left[-{\mathrm{e}}^{-\left({\breve{\mathbb{E}}}_x^{\left(q\right)}-2{\overline{\mathbb{V}}}_q\right)/{\overline{T}}_q}\right]-{\mathrm{Li}}_2\left[-{\mathrm{e}}^{-\left({\breve{\mathbb{E}}}_x^{\left(q\right)}+2{\overline{\mathbb{V}}}_q\right)/{\overline{T}}_q}\right]\right\},\hfill \end{array}$$

(62)

where

$$\begin{array}{c}\hfill {\breve{\mathbb{E}}}_x^{\left(q\right)}=\sqrt{x^2+{\left(\frac{{\breve{\Delta }}_q}{{\epsilon }_{Fq}}\right)}^2}.\end{array}$$

(63)

Introducing the critical effective superfluid velocities [30]

$$\begin{array}{c}\hfill {\breve{\mathbb{V}}}_{cq}^{\left(0\right)}=\frac{\mathrm{e}}{2}\frac{{\breve{\Delta }}_q^{\left(0\right)}}{\hslash k_{Fq}},\end{array}$$

(64)

the approximate pairing gap Equation (59) and the functions (62) can be equivalently expressed in terms of the reduced temperature ${\tilde{T}}_q\equiv T/{\breve{T}}_{cq}^{\left(0\right)}$ and the reduced effective superfluid velocity ${\tilde{\mathbb{V}}}_q\equiv {\mathsf{𝕍}}_{\mathit{q}}/{\breve{\mathbb{V}}}_{cq}^{\left(0\right)}$ as follows:

$$\begin{array}{cc}\hfill log\left(\frac{{\breve{\Delta }}_q}{{\breve{\Delta }}_q^{\left(0\right)}}\right)\approx & {\int }_0^{+\infty }\frac{\mathrm{d}y}{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}\left\{\frac{2}{\pi }{\mathrm{e}}^{\gamma -1}\frac{{\tilde{T}}_q}{{\tilde{\mathbb{V}}}_q}log\left[cosh\left(\frac{\pi }{2}{\mathrm{e}}^{-\gamma }\left(\frac{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}{{\tilde{T}}_q}+\frac{\mathrm{e}}{2}\frac{{\tilde{\mathbb{V}}}_q}{{\tilde{T}}_q}\right)\right)\right.\right.\hfill \\ \hfill & \times \mathsf{s}\mathsf{e}\mathsf{c}\mathsf{h}\left.\left.\left(\frac{\pi }{2}{\mathrm{e}}^{-\gamma }\left(\frac{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}{{\tilde{T}}_q}-\frac{\mathrm{e}}{2}\frac{{\tilde{\mathbb{V}}}_q}{{\tilde{T}}_q}\right)\right)\right]-1\right\},\hfill \end{array}$$

(65)

or using Equation (60)

$$\begin{array}{cc}\hfill log\left(\frac{{\breve{\Delta }}_q}{{\breve{\Delta }}_q^{\left(0\right)}}\right)\approx & \frac{2}{\pi }{\mathrm{e}}^{\gamma -1}\frac{{\tilde{T}}_q}{{\tilde{\mathbb{V}}}_q}{\int }_0^{+\infty }\frac{\mathrm{d}y}{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}log\left[\frac{1+{\mathrm{e}}^{-\pi {\mathrm{e}}^{-\gamma }\left({\tilde{\mathbb{E}}}_y^{\left(q\right)}+\frac{\mathrm{e}}{2}{\tilde{\mathbb{V}}}_q\right)/{\tilde{T}}_q}}{1+{\mathrm{e}}^{-\pi {\mathrm{e}}^{-\gamma }\left({\tilde{\mathbb{E}}}_y^{\left(q\right)}-\frac{\mathrm{e}}{2}{\tilde{\mathbb{V}}}_q\right)/{\tilde{T}}_q}}\right],\hfill \end{array}$$

(66)

and

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_q& =\frac{12}{\pi }{\mathrm{e}}^{\gamma -2}\frac{{\tilde{T}}_q}{{\tilde{\mathbb{V}}}_q^2}{\int }_0^{+\infty }\mathrm{d}ylog\left\{\left[1+exp\left(-\pi {\mathrm{e}}^{-\gamma }\left(\frac{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}{{\tilde{T}}_q}-\frac{\mathrm{e}}{2}\frac{{\tilde{\mathbb{V}}}_q}{{\tilde{T}}_q}\right)\right)\right]\right.\hfill \\ \hfill & \times \left.\left[1+exp\left(-\pi {\mathrm{e}}^{-\gamma }\left(\frac{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}{{\tilde{T}}_q}+\frac{\mathrm{e}}{2}\frac{{\tilde{\mathbb{V}}}_q}{{\tilde{T}}_q}\right)\right)\right]\right\}\hfill \\ \hfill & +\frac{24}{{\pi }^2}{\mathrm{e}}^{2\gamma -3}\frac{{\tilde{T}}_q^2}{{\tilde{\mathbb{V}}}_q^3}{\int }_0^{+\infty }\mathrm{d}y\left\{{\mathrm{Li}}_2\left[-exp\left(-\pi {\mathrm{e}}^{-\gamma }\left(\frac{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}{{\tilde{T}}_q}-\frac{\mathrm{e}}{2}\frac{{\tilde{\mathbb{V}}}_q}{{\tilde{T}}_q}\right)\right)\right]\right.\hfill \\ & -\left.{\mathrm{Li}}_2\left[-exp\left(-\pi {\mathrm{e}}^{-\gamma }\left(\frac{{\tilde{\mathbb{E}}}_y^{\left(q\right)}}{{\tilde{T}}_q}+\frac{\mathrm{e}}{2}\frac{{\tilde{\mathbb{V}}}_q}{{\tilde{T}}_q}\right)\right)\right]\right\},\hfill \end{array}$$

(67)

with

$$\begin{array}{c}\hfill {\tilde{\mathbb{E}}}_y^{\left(q\right)}=\sqrt{y^2+{\left(\frac{{\breve{\Delta }}_q}{{\breve{\Delta }}_q^{\left(0\right)}}\right)}^2}.\end{array}$$

(68)

These alternative formulations show that ${\breve{\Delta }}_q/{\breve{\Delta }}_q^{\left(0\right)}$ and ${\mathsf{\Phi }}_q$ are universal functions of suitably rescaled temperature and effective superfluid velocity, independently of the nucleon species under consideration, the composition, and the details of the adopted nuclear-energy-density functional.

3. Application to Neutron Stars

Although the entrainment matrix can be written in the deceptively simple analytical form (19), its dependencies on the temperature and on the superfluid velocities remain implicit and highly nontrivial. To obtain actual values, numerical solutions of Equations (43) and (44) are needed. In this work, we have considered the Brussels–Montreal functionals, whose main features are described in Section 3.1. Results are presented in the subsequent sections.

3.1. Brussels–Montreal Functionals

The Brussels–Montreal functionals from BSk16 and beyond (see [31,32] for a brief overview) were constructed from extended Skyrme effective nucleon-nucleon interactions, whose parameters were precision-fitted to essentially all experimental nuclear data on atomic masses and charge radii while ensuring realistic properties of homogeneous nuclear matter (neutron-matter equation of state, effective masses, symmetry energy, incompressibility coefficient, pairing gaps).

The functional derivatives of the energy $E_{\mathrm{nuc}}^j$ with respect to $X_0$ and $X_1$ appearing in the effective masses, the matrix ${\mathcal{I}}_{q{q}^{\prime }}$ and the entrainment matrix are expressible in terms of the parameters of the effective interaction as [16]

$$\begin{array}{c}\hfill \frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}=\frac{3}{16}{t}_1+\frac{1}{4}{t}_2\left(\frac{5}{4}+x_2\right)+\frac{3}{16}{t}_4{n}^{\beta }+\frac{1}{4}{t}_5\left(\frac{5}{4}+x_5\right)n^{\gamma }\end{array}$$

(69)

$$\begin{array}{c}\hfill \frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}=-\frac{1}{8}{t}_1\left(\frac{1}{2}+x_1\right)+\frac{1}{8}{t}_2\left(\frac{1}{2}+x_2\right)-\frac{1}{8}{t}_4\left(\frac{1}{2}+x_4\right)n^{\beta }+\frac{1}{8}{t}_5\left(\frac{1}{2}+x_5\right)n^{\gamma }.\end{array}$$

(70)

The potentials in homogeneous matter read (recalling the shorthand notations $n\equiv n_0$, $\mathit{j}\equiv {\mathit{j}}_{\mathit{0}}$ and $\tau \equiv {\tau }_0$)

$$\begin{array}{ccc}\hfill U_q& =& t_0\left[\left(1+\frac{1}{2}{x}_0\right)n-\left(\frac{1}{2}+x_0\right)n_q\right]+\frac{1}{4}{t}_1\left[\left(1+\frac{1}{2}{x}_1\right)\tau -\left(\frac{1}{2}+x_1\right){\tau }_q\right]\hfill \\ & +& \frac{1}{4}{t}_2\left[\left(1+\frac{1}{2}{x}_2\right)\tau +\left(\frac{1}{2}+x_2\right){\tau }_q\right]\hfill \\ & +& \frac{1}{12}{t}_3{n}^{\alpha -1}\left[\left(1+\frac{1}{2}{x}_3\right)(2+\alpha )n^2-\left(\frac{1}{2}+x_3\right)\left(2n{n}_q+\alpha \sum _{q^{\prime }=n,p}{n}_{q^{\prime }}^2\right)\right]\hfill \\ & +& \frac{1}{4}{t}_4{n}^{\beta -1}\left[\left(1+\frac{1}{2}{x}_4\right)(1+\beta )n\tau -\left(\frac{1}{2}+x_4\right)\left(n{\tau }_q+\beta \sum _{q^{\prime }=n,p}{n}_{q^{\prime }}{\tau }_{q^{\prime }}\right)\right]\hfill \\ & +& \frac{1}{4}{t}_5{n}^{\gamma -1}\left[\left(1+\frac{1}{2}{x}_5\right)(1+\gamma )n\tau +\left(\frac{1}{2}+x_5\right)\left(n{\tau }_q+\gamma \sum _{q^{\prime }=n,p}{n}_{q^{\prime }}{\tau }_{q^{\prime }}\right)\right]\hfill \\ & +& \frac{1}{8}{t}_4\beta n^{\beta -1}\left[(x_4-1){\mathit{j}}^2+4\left(\frac{1}{2}+x_4\right){\mathit{j}}_{\mathit{q}}·({\mathit{j}}_{\mathit{q}}-\mathit{j})\right]\hfill \\ & -& \frac{1}{8}{t}_5\gamma n^{\gamma -1}\left[3\left(x_5+1\right){\mathit{j}}^2+4\left(\frac{1}{2}+x_5\right){\mathit{j}}_{\mathit{q}}·\left({\mathit{j}}_{\mathit{q}}-\mathit{j}\right)\right]\hfill \\ & +& \frac{1}{4}\sum _{q^{\prime }=n,p}\frac{\partial v^{\pi q^{\prime }}}{\partial n_q}{\tilde{n}}_{q^{\prime }}^2.\hfill \end{array}$$

(71)

The functional derivative of the energy E with respect to the square modulus of the abnormal density ${\tilde{n}}_q$ is related to the strength $v^{\pi q}$ of the effective pairing interaction as

$$\begin{array}{c}\hfill \frac{\delta E}{\delta |{\tilde{n}}_q{|}^2}=\frac{1}{4}{v}^{\pi q}.\end{array}$$

(72)

In most existing functionals, $v^{\pi q}$ is expressed [33] as the sum of a constant “volume” term and a “surface term” proportional to the density n to some power with parameters adjusted empirically to reproduce the average pairing gaps in some finite nuclei [34]. Such functionals may thus lead to unreliable predictions when applied to homogeneous nuclear matter [35]. On the contrary, the pairing strengths $v^{\pi q}[n_n,n_p]<0$ of the Brussels–Montreal functionals were determined so as to reproduce the ${}^1{S}_0$ pairing gaps in infinite homogeneous neutron matter and in symmetric nuclear matter at $T=0$ and in the absence of currents (these reference gaps will be denoted by ${\widehat{\Delta }}_{NM}$ and ${\widehat{\Delta }}_{SM}$ respectively), as obtained from many-body calculations using realistic potentials (see [35,36,37] for details). Very accurate analytical expressions for the pairing strengths were obtained in [38]:

$$\begin{array}{ccc}\hfill v^{\pi q}\left[n_n,n_p\right]& =& -\frac{4{\pi }^2{\hslash }^2{\mathsf{\Sigma }}_q}{m{k}_{Fq}}\left[log\left(\frac{64m_q^{\oplus }{\mathsf{\Sigma }}_q}{m}\frac{{\epsilon }_{Fq}{\epsilon }_{\mathsf{\Lambda }}}{{\widehat{\Delta }}_q^2}\right)\right.\hfill \\ & +& {\left.2\sqrt{1+\frac{{\epsilon }_{\mathsf{\Lambda }}}{{\epsilon }_{Fq}}\frac{m}{m_q^{\oplus }{\mathsf{\Sigma }}_q}}-2log\left(1+\sqrt{1+\frac{{\epsilon }_{\mathsf{\Lambda }}}{{\epsilon }_{Fq}}\frac{m}{m_q^{\oplus }{\mathsf{\Sigma }}_q}}\right)-4\right]}^{-1}.\hfill \end{array}$$

(73)

The parameters ${\mathsf{\Sigma }}_q$ are used here to distinguish Brussels–Montreal functionals BSk17-29 [18,36,39,40,41] which neglect self-energy corrections (${\mathsf{\Sigma }}_q=1$) from the most recent series BSk30-32 [32] which include them (${\mathsf{\Sigma }}_q=m/m_q^{\oplus }$). Since reference pairing gaps $\widehat{\Delta }(n_n,n_p)$ for arbitrary composition are needed, the following interpolation ansatz was adopted in [36] for BSk17 and subsequent functionals:

$${\widehat{\Delta }}_q(n_n,n_p)={\widehat{\Delta }}_{SM}\left(n\right)(1-|\eta \left|\right)\pm {\widehat{\Delta }}_{NM}\left(n_q\right)\eta \frac{n_q}{n},$$

(74)

where $\eta =(n_n-n_p)/n$ and the upper (lower) sign is to be taken for $q=n\left(p\right)$. Because this parametrization is empirical, we have found that ${\widehat{\Delta }}_q(n_n,n_p)$ may become negative depending on the composition and density n. In such cases, we merely set ${\widehat{\Delta }}_q(n_n,n_p)=0$. As for the nucleon mass, it is defined as $m=2{(1/m_n+1/m_p)}^{-1}$.

For numerical calculations, we will adopt the Brussels–Montreal functional BSk24 [18]. The reference pairing gaps were taken from the extended Brueckner–Hartree–Fock calculations of [42]. The associated parameters are indicated in

Table 1. Parameters of the functional BSk24 [18]. The unit of length is femtometer and the unit of energy is megaelectronvolt.

${\mathit{t}}_0$	−3970.29
$t_1$	395.766
$t_2$	${10}^{-5}$
$t_3$	22648.6
$t_4$	−100.000
$t_5$	−150.000
$x_0$	0.894371
$x_1$	0.0563535
$x_2$	−0.138961 ×${10}^9$
$x_3$	1.05119
$x_4$	2.00000
$x_5$	−11.0000
$\alpha$	1/12
$\beta$	1/2
$\gamma$	1/12
${\epsilon }_{\mathsf{\Lambda }}$	16.0
${\mathsf{\Sigma }}_q$	1

and

Table 2. Parameters of the reference gaps from [37]. The unit of length is femtometer and the unit of energy is megaelectronvolt. With kind permission of The European Physical Journal (EPJ).

	${\mathsf{\Delta }}_0$	${\mathit{k}}_1$	${\mathit{k}}_2$	${\mathit{k}}_3$	${\mathit{k}}_{\mathbf{max}}$
${\widehat{\Delta }}_{SM}$	133.779	0.943146	1.52786	2.11577	1.51
${\widehat{\Delta }}_{NM}$	14.9003	1.18847	1.51854	0.639489	1.52

. The reference gaps can be conveniently represented as

$${\widehat{\Delta }}_{SM}\left(n\right)=H(k_{\max }-k_F){\Delta }_0\frac{k_F^3}{k_F^2+k_1^2}\frac{{(k_F-k_2)}^2}{{(k_F-k_2)}^2+k_3^2},$$

(75)

$${\widehat{\Delta }}_{NM}\left(n_q\right)=H(k_{\max }-k_{Fq}){\Delta }_0\frac{k_{Fq}^2}{k_{Fq}^2+k_1^2}\frac{{(k_{Fq}-k_2)}^2}{{(k_{Fq}-k_2)}^2+k_3^2},$$

(76)

where $k_F={(3{\pi }^{2}n/2)}^{1/3}$, H is the Heaviside unit-step function, and $k_1$, $k_2$, $k_3$ and $k_{\max }$ are fitted parameters. The functional BSk24 has been recently employed for determining the composition and the equation of state of dense matter throughout all regions of a neutron star [19,20] including the pasta mantle [21] and allowing for strong magnetic fields [22]. More importantly, as shown in [23,43,44], this functional turns out to be in very good agreement with existing astrophysical observations including those from the binary neutron-star merger GW170817 [45] as well as from PSR J0740+6620 and PSR J0030+0451 by the Neutron star Interior Composition Explorer (NICER) [46,47,48,49]. Results for the entrainment matrix at finite temperatures but in the absence of superflows have been recently published in [11] within Landau’s theory using values for the Landau parameters calculated for the Brussels–Montreal functionals including BSk24 and setting arbitrarily the critical temperatures. We will present here consistent numerical results for the pairing gaps, chemical potentials and entrainment matrix for arbitrary temperatures and superfluid velocities in different regions of neutron-star cores.

3.2. Numerical Implementation

The TDHFB equations are solved as follows. We first compute the pairing gaps ${\Delta }_q^{\left(0\right)}$ at zero temperature and in the absence of currents by solving Equations (43) and (44) for $T=0$ and ${\mathsf{𝕍}}_{\mathit{q}}=0$ via a root-finding method with a precision of ${10}^{-8}$, searching around the approximate solutions ${\mu }_q\approx {\epsilon }_{Fq}$ and the following expression given by Equation (14) in [38]:

$$\begin{array}{c}\hfill {\displaystyle {\Delta }_q^{\left(0\right)}\approx \frac{8\sqrt{{\epsilon }_{Fq}{\epsilon }_{\mathsf{\Lambda }}}}{1+\sqrt{1+{\epsilon }_{\mathsf{\Lambda }}/{\epsilon }_{Fq}}}exp\left[\frac{1}{v^{\pi q}\breve{\mathcal{D}}\left(0\right)}+\sqrt{1+\frac{{\epsilon }_{\mathsf{\Lambda }}}{{\epsilon }_{Fq}}}-2\right].}\end{array}$$

(77)

In a second stage, we use this solution to determine iteratively the pairing gaps ${\Delta }_q$ and the reduced chemical potentials ${\mu }_q$ at finite temperature $T>0$ and for given effective superfluid velocities ${\mathbb{V}}_{\mathit{n}}$ and ${\mathsf{𝕍}}_{\mathit{p}}$. An initial guess for ${\Delta }_q$ is obtained by solving Equation (59) using for ${\breve{\Delta }}_q^{\left(0\right)}$ the gap obtained previously. With this first estimate of the gap, Equation (44) is solved using the Newton-Raphson method and ${\mu }_q\approx {\epsilon }_{Fq}$ as the initial guess. Substituting these first estimates for ${\Delta }_q$ and ${\mu }_q$ in the right-hand side of Equation (43) leads to a new estimate for the pairing gap ${\Delta }_q$, which is injected in Equation (44) to refine the chemical potential ${\mu }_q$. The process is repeated until the difference in the pairing gaps between two successive iterations lies below ${10}^{-4}{\Delta }_q^{\left(0\right)}$. Having found ${\Delta }_q$ and ${\mu }_q$, the functions ${\mathcal{Y}}_q$ are calculated from Equation (38). The entrainment matrix can be easily inferred from Equations (19)–(25) together with Equations (69) and (70).

3.3. ${}^1$S${}_0$ Pairing Gaps

The ${}^1$S${}_0$ neutron and proton pairing gaps ${\Delta }_q^{\left(0\right)}$ for $npe\mu$ matter in beta-equilibrium at $T=0$ and ${\mathsf{𝕍}}_{\mathit{q}}=0$ are displayed in Figure 1 and Figure 2 at densities relevant for the outer core of neutron stars above the crust-core transition at density $n_{\mathrm{cc}}=0.08076$ fm${}^{-3}\approx 0.5n_0$, where $n_0=0.1578$ fm${}^{-3}$ is the nuclear saturation density with the corresponding mass density ${\rho }_0=m{n}_0=2.654\times {10}^{14}$ g cm${}^{-3}$. We have made use of the composition calculated in [19].

The approximate formula (77) is found to be in excellent agreement with the exact results, the deviations being contained within the thickness of the solid lines. With neutron-star matter containing only a few percents of protons, the reference pairing gaps for neutrons (74) are approximately given by that in pure neutron matter ${\widehat{\Delta }}_n(n_n,n_p)\approx {\widehat{\Delta }}_{NM}\left(n_n\right)$, as obtained from the many-body calculations of [42] using realistic potentials. On the contrary, the reference pairing gaps for protons are mainly determined by the interpolation ${\widehat{\Delta }}_p(n_n,n_p)\approx {\widehat{\Delta }}_{NM}\left(n_p\right)n_p/n$. This explains why the proton gaps ${\Delta }_p$ are significantly smaller than the neutron ones ${\Delta }_n$ unlike those usually employed in neutron-star studies, as e.g., in [4]. This result could reveal a deficiency of the interpolation (74). On the other hand, the proton pairing gaps remain highly uncertain (see, e.g., [25,50,51]). Recent many-body calculations [52] taking into account medium-polarization effects through self-energy and vertex corrections lead to very small proton pairing gaps in neutron-star matter of comparable magnitudes to those plotted in Figure 2. This study also shows that the three-body interactions, especially those between two protons and one neutron, reduce considerably the domain of temperatures and densities over which protons are superfluid (see also [53,54]).

The variations of the neutron and proton pairing gaps with temperature and effective superfluid velocity are found to be essentially independent of density when considering the dimensionless ratios ${\Delta }_q/{\Delta }_q^{\left(0\right)}$, $T/T_{cq}^{\left(0\right)}$ and ${\mathsf{𝕍}}_{\mathit{q}}/{\mathbb{V}}_{cq}^{\left(0\right)}$, with

$$\begin{array}{c}\hfill T_{cq}^{\left(0\right)}=\frac{{\mathrm{e}}^{\gamma }}{k_{\mathrm{B}}\pi }{\Delta }_q^{\left(0\right)},\end{array}$$

(78)

$$\begin{array}{c}\hfill {\mathbb{V}}_{cq}^{\left(0\right)}=\frac{\mathrm{e}}{2}\frac{{\Delta }_q^{\left(0\right)}}{\hslash k_{Fq}}.\end{array}$$

(79)

As shown in Figure 3 and Figure 4, the gaps for both neutrons and protons decrease monotonically with increasing temperature and effective superfluid velocity due to the excitation of quasiparticles. For vanishing effective superfluid velocities ${\mathsf{𝕍}}_{\mathit{q}}=0$ (i.e., in the absence of mass flow ${\mathit{\rho }}_{\mathit{q}}=\mathit{0}$), the temperature dependence of the pairing gaps is well fitted by the following expression [55]:

$$\begin{array}{c}\hfill \frac{{\Delta }_q(T\le T_{cq}^{\left(0\right)},{\mathsf{𝕍}}_{\mathit{q}}=0)}{{\Delta }_q^{\left(0\right)}}=\frac{{\mathrm{e}}^{\gamma }}{\pi }\sqrt{1-\frac{T}{T_{cq}^{\left(0\right)}}}\left(1.456\frac{T}{T_{cq}^{\left(0\right)}}-0.157\sqrt{\frac{T}{T_{cq}^{\left(0\right)}}}+1.764\right).\end{array}$$

(80)

This same formula was applied in [13] to evaluate the entrainment matrix. At zero temperature, the pairing gap remains equal to ${\Delta }_q^{\left(0\right)}$ until the effective superfluid velocity ${\mathsf{𝕍}}_{\mathit{q}}$ reaches Landau’s critical velocity ${\mathbb{V}}_{Lq}$, which for BCS condensates is given by [56]

$$\begin{array}{c}\hfill {\mathbb{V}}_{Lq}=\frac{{\Delta }_q^{\left(0\right)}}{\hslash k_{Fq}}.\end{array}$$

(81)

Beyond this point, the pairing gap decreases with increasing effective superfluid velocity and vanishes for ${\mathsf{𝕍}}_{\mathit{q}}={\mathbb{V}}_{cq}^{\left(0\right)}$. We find that this behavior is well reproduced by the following interpolating formula:

$$\begin{array}{c}\hfill \frac{{\Delta }_q(T=0,{\mathbb{V}}_{Lq}\le {\mathsf{𝕍}}_{\mathit{q}}\le {\mathbb{V}}_{cq}^{\left(0\right)})}{{\Delta }_q^{\left(0\right)}}=0.5081\sqrt{1-\frac{{\mathsf{𝕍}}_{\mathit{q}}}{{\mathbb{V}}_{cq}^{\left(0\right)}}}\left(3.312\frac{{\mathsf{𝕍}}_{\mathit{q}}}{{\mathbb{V}}_{cq}^{\left(0\right)}}-3.811\sqrt{\frac{{\mathbb{V}}_{cq}^{\left(0\right)}}{{\mathsf{𝕍}}_{\mathit{q}}}}+5.842\right).\end{array}$$

(82)

The maximum relative error does not exceed 0.13%.

The critical temperature and critical effective superfluid velocity delimiting the superfluid and normal phases, plotted in Figure 5 is well fitted by the following expression:

$$\begin{array}{c}\hfill {\displaystyle \frac{T_{cq}}{T_{cq}^{\left(0\right)}}({\mathsf{𝕍}}_{\mathit{q}}\le {\mathbb{V}}_{cq}^{\left(0\right)})\simeq {\left[1-{\left(\frac{{\mathsf{𝕍}}_{\mathit{q}}}{{\mathbb{V}}_{cq}^{\left(0\right)}}\right)}^2\right]}^{2/5}.}\end{array}$$

(83)

This interpolation is valid for both neutrons and protons. The errors are contained within the thickness of the lines in Figure 5.

The universality observed in the superfluid properties of both neutrons and protons (after a suitable choice of normalizations) is the consequence of the weak-coupling regime, as discussed in Section 2.5. Indeed, as shown in Equation (65), the normalized pairing gaps ${\breve{\Delta }}_q/{\breve{\Delta }}_q^{\left(0\right)}$ are independent of the pairing strength $v^{\pi q}$ (hence also of the associated pairing cutoff ${\epsilon }_{\mathsf{\Lambda }}$) and depend only on the rescaled temperature $T/T_{cq}^{\left(0\right)}$ and effective superfluid velocity ${\mathsf{𝕍}}_{\mathit{q}}/{\mathbb{V}}_{cq}^{\left(0\right)}$. Estimating the exact pairing gaps ${\Delta }_q^{\left(0\right)}$ from Equation (77) and substituting in Equation (65) lead to a very good approximation for the exact pairing gaps ${\Delta }_q/{\Delta }_q^{\left(0\right)}$ at finite temperatures and arbitrary effective superfluid velocities. The largest absolute deviations are found at the crust-core interface: they are of order of ${10}^{-3}$ for neutrons and lie within the numerical errors for protons.

3.4. Reduced Chemical Potentials

The TDHFB theory allows the determination of the chemical potentials consistently with the pairing gaps. Let us recall that in Landau’s theory adopted in previous studies [13,14,15], the reduced chemical potential ${\mu }_q$ was approximated by the corresponding Fermi energy ${\epsilon }_{Fq}$; effects induced by pairing, temperature, and currents were therefore ignored. To assess the precision of this approximation, we have computed ${\mu }_q$ numerically by solving simultaneously Equations (43) and (44) varying the temperature and the neutron effective superfluid velocity. The largest relative errors between ${\mu }_q$ and the Fermi energy ${\epsilon }_{Fq}$ we have found (at the crust-core interface) are 0.14% for neutrons and 0.052% for protons. Such errors have been obtained for low temperatures and small effective superfluid velocities for which pairing effects are the most important. Focusing on these conditions, we have plotted in Figure 6 the ratio ${\mu }_q/{\epsilon }_{Fq}$ as a function of density. As expected, the higher the density, the more precise are Landau’s approximations. To a large extent, the small deviations between ${\mu }_q$ and ${\epsilon }_{Fq}$ stem from the rather small pairing gaps predicted by the functional BSk24. Larger deviations cannot be excluded if another functional is adopted. In any case, let us recall that both Equations (43) and (44) should be solved simultaneously to obtain fully consistent pairing gaps and chemical potentials.

3.5. Functions ${\mathcal{Y}}_q$

Having computed the pairing gaps ${\Delta }_q$ as well as the reduced chemical potentials ${\mu }_q$ at finite temperatures and for arbitrary effective superfluid velocities by solving Equations (43) and (44), we can now evaluate the functions ${\mathcal{Y}}_q$ from Equation (38). When expressed in terms of the dimensionless ratios $T/T_{cq}^{\left(0\right)}$ and ${\mathsf{𝕍}}_{\mathit{q}}/{\mathbb{V}}_{cq}^{\left(0\right)}$, results are found to be essentially independent of density and are summarized in Figure 7 and Figure 8.

The functions ${\mathcal{Y}}_q$ are well approximated by the functions ${\mathsf{\Phi }}_q$ defined by Equation (62) where the pairing gaps ${\breve{\Delta }}_q$ are computed from Equation (59) and provided ${\breve{\Delta }}_q^{\left(0\right)}$ are evaluated from Equation (77). Since the deviations decrease with increasing density, we have focused on the crust-core interface. The absolute errors are found to be at most of order ${10}^{-3}$ for ${\mathcal{Y}}_n$ and ${10}^{-4}$ for ${\mathcal{Y}}_p$. It follows from Equation (67) that the function ${\mathcal{Y}}_q$ is universal. In the absence of superflow ${\mathsf{𝕍}}_{\mathit{q}}=0$, the temperature dependence of the functions ${\mathcal{Y}}_q$ can be well fitted by the following expression [57] (errors not exceeding 2.6%):

$$\begin{array}{cc}\hfill {\mathcal{Y}}_q(T\le T_{cq}^{\left(0\right)},{\mathsf{𝕍}}_{\mathit{q}}=0)& \simeq {\left[0.9443+\sqrt{0.{0557}^2+{\left(0.1886\frac{\pi }{{\mathrm{e}}^{\gamma }}\frac{{\Delta }_q\left(T\right)}{{\Delta }_q^{\left(0\right)}}\frac{T_{cq}^{\left(0\right)}}{T}\right)}^2}\right]}^{1/2}\hfill \\ \hfill & \times exp\left[1.753-\sqrt{1.{753}^2+{\left(\frac{\pi }{{\mathrm{e}}^{\gamma }}\frac{{\Delta }_q\left(T\right)}{{\Delta }_q^{\left(0\right)}}\frac{T_{cq}^{\left(0\right)}}{T}\right)}^2}\right]\hfill \end{array}$$

(84)

with ${\Delta }_q\left(T\right)$ computed using the interpolation (80).

3.6. Effective versus True Superfluid Velocities

The results we have presented so far have been conveniently expressed in terms of the effective superfluid velocities ${\mathsf{𝕍}}_{\mathit{q}}$, which are related to the original superfluid velocities ${\mathsf{𝕍}}_{\mathit{q}}$ by Equations (29) and (32). These relations are highly nontrivial, recalling that the coefficients ${\mathcal{I}}_{q{q}^{\prime }}$, defined by (21)–(25), depend on ${\mathsf{𝕍}}_{\mathit{q}}$ through the functions ${\mathcal{Y}}_q$.

So far, we have treated the effective superfluid velocities as free parameters. In reality however, ${\mathsf{𝕍}}_{\mathit{n}}$ and ${\mathsf{𝕍}}_{\mathit{p}}$ are determined by the dynamics of the star, as pointed out in the previous analysis of entrainment effects in [15]. In particular, in the study of low-frequency oscillations, it is a very good approximation to assume that the electric current in the normal frame vanishes, as shown in the classical work of [58]. Considering that leptons are co-moving with quasiparticle excitations, the previous condition reads ${\mathit{v}}_{\mathit{p}}=\mathit{0}$ (in the normal frame). It immediately follows from Equation (53) that ${\mathsf{𝕍}}_{\mathit{p}}=\mathit{0}$. In the following, we will restrict to this case as in [15] since it is of most physical interest. Under such condition, the vectors ${\mathit{V}}_{\mathit{n}}$ and ${\mathit{V}}_{\mathit{p}}$ are aligned, and are given by

$$\begin{array}{ccc}\hfill {\mathit{V}}_{\mathit{n}}& =& \left[1-{\displaystyle \frac{2}{{\hslash }^2}}\left({\displaystyle \frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}}+{\displaystyle \frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}}\right)\left(m_p^{\oplus }{n}_p+m_n^{\oplus }{n}_n{\mathcal{Y}}_n\right)+{\displaystyle \frac{16}{{\hslash }^4}}{\displaystyle \frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}}{\displaystyle \frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}}{m}_p^{\oplus }{n}_p{m}_n^{\oplus }{n}_n{\mathcal{Y}}_n\right]\hfill \\ & & \times \mathsf{\Xi }{\displaystyle \frac{m_n^{\oplus }}{m}}{\mathsf{𝕍}}_{\mathit{n}},\hfill \end{array}$$

(85)

$$\begin{array}{c}\hfill {\mathit{V}}_{\mathit{p}}=\left[\frac{2}{{\hslash }^2}\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\right)\left(1-{\mathcal{Y}}_n\right)m_n^{\oplus }{n}_n\right]\mathsf{\Xi }{\displaystyle \frac{m_p^{\oplus }}{m}}{\mathsf{𝕍}}_{\mathit{n}},\end{array}$$

(86)

$$\begin{array}{c}\hfill \mathsf{\Xi }={\left[1-{\displaystyle \frac{2}{{\hslash }^2}}\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}+\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\right)\left(m_p^{\oplus }{n}_p+m_n^{\oplus }{n}_n\right)+\frac{16}{{\hslash }^4}\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}{m}_p^{\oplus }{n}_p{m}_n^{\oplus }{n}_n\right]}^{-1}.\end{array}$$

(87)

These superfluid velocities depend on the baryon density n, the temperature T and the neutron effective superfluid velocity ${\mathsf{𝕍}}_{\mathit{n}}$. Please note that under Landau’s approximations, the norm of (85) reduces to Equation (79) of [15] (these authors adopted the notation ${\tilde{\mathsf{\Phi }}}_q$ for ${\mathcal{Y}}_q$, $m_n^{\ast }$ for the neutron effective mass $m_n^{\oplus }$, ${\tilde{V}}_n$ for the neutron effective superfluid velocity ${\mathsf{𝕍}}_{\mathit{n}}$ and $v_n$ for the neutron true superfluid velocity $V_n$; the Landau parameters $F_1^{q{q}^{\prime }}$ are given by Equation (100) of [17]).

Results for the norms, considering $npe\mu$ matter in beta-equilibrium, are displayed in Figure 9 and Figure 10 for two different densities. These superfluid velocities are only defined in the superfluid phase, for effective superfluid velocities and temperatures lower than their associated critical values given by (83). Indeed, in the normal phase, the abnormal densities ${\tilde{n}}_q$ vanish identically and the associated superfluid velocities are therefore ill defined. However, this has no physical implications since the superfluid velocities are irrelevant in this case.

Although the neutron superfluid velocity is roughly equal to the effective superfluid velocity, $V_n\approx {\mathbb{V}}_{\mathit{n}}$, the proton superfluid velocity exhibits a more complicated behavior as a function of ${\mathsf{𝕍}}_{\mathit{n}}$. From Equation (86), we have $V_p\propto \left(1-{\mathcal{Y}}_n\right){\mathsf{𝕍}}_{\mathit{n}}$. For sufficiently low neutron effective superfluid velocities, ${\mathcal{Y}}_n\approx 0$ therefore $V_p$ increases linearly with ${\mathsf{𝕍}}_{\mathit{n}}$. However, ${\mathcal{Y}}_n$ increases with ${\mathsf{𝕍}}_{\mathit{n}}$ leading to a decrease of $V_p$ (for ${\mathsf{𝕍}}_{\mathit{n}}\simeq {\mathbb{V}}_{Ln}$), which vanishes when ${\mathsf{𝕍}}_{\mathit{n}}={\mathbb{V}}_{cn}^{\left(0\right)}$ corresponding to ${\mathcal{Y}}_n=1$.

3.7. Entrainment Matrix

Having computed the pairing gaps, chemical potentials, and functions ${\mathcal{Y}}_q$ we can now determine the entrainment matrix from Equation (19). To better see the influence of temperature and superflows, results will be compared to those obtained at zero temperature and in the limit of small currents (conditions for which pairing can be ignored) using the following expression that we have previously calculated within the time-dependent Hartree–Fock (TDHF) theory [16]:

$$\begin{array}{c}\hfill {\rho }_{np}^{\mathrm{TDHF}}={\rho }_{pn}^{\mathrm{TDHF}}=-\frac{2}{{\hslash }^2}{\rho }_n{\rho }_p\left(\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_0}-\frac{\delta E_{\mathrm{nuc}}^j}{\delta X_1}\right),\end{array}$$

(88)

$$\begin{array}{c}\hfill {\rho }_{nn}^{\mathrm{TDHF}}={\rho }_n-{\rho }_{np}^{\mathrm{TDHF}},\end{array}$$

(89)

$$\begin{array}{c}\hfill {\rho }_{pp}^{\mathrm{TDHF}}={\rho }_p-{\rho }_{pn}^{\mathrm{TDHF}}.\end{array}$$

(90)

These matrix elements are shown in Figure 11 for $npe\mu$ matter in beta-equilibrium at densities relevant for the outer core of neutron stars. Results within the TDHFB theory for finite temperatures and different neutron effective superfluid velocities (recalling that we set ${\mathsf{𝕍}}_{\mathit{p}}=0$ as discussed in Section 3.6) are plotted in Figure 12, Figure 13 and Figure 14 at the crust-core interface, and in Figure 15, Figure 16 and Figure 17 at the saturation density.

Quite remarkably, the entrainment matrix at $T=0$ remains independent of the neutron effective superfluid velocity provided the latter does not exceed Landau’s critical velocity. In other words, the expressions obtained in [16] in the limit of vanishing small effective superfluid velocities are actually exact for any effective superfluid velocity lying below Landau’s critical value. This can be traced back to the vanishing of the function ${\mathcal{Y}}_q$ for ${\mathsf{𝕍}}_{\mathit{q}}\le {\mathbb{V}}_{Lq}$, as can be seen in Figure 7. This also entails that the entrainment matrix does not depend on the pairing gaps in this regime, thus justifying a posteriori our application of the TDHF theory [16] instead of TDHFB [17] since the gaps can thus artificially be set to zero. However, the TDHFB theory is still required for the determination of the actual value for ${\mathbb{V}}_{Lq}$.

At finite but sufficiently low temperatures, the entrainment matrix remains weakly dependent on the neutron effective superfluid velocity provided ${\mathsf{𝕍}}_{\mathit{q}}\le {\mathbb{V}}_{Lq}$. For higher neutron effective superfluid velocities, the entrainment matrix elements ${\rho }_{nn}$ and ${\rho }_{np}={\rho }_{pn}$ are reduced, even at $T=0$. The element ${\rho }_{pp}$ is essentially independent of ${\mathsf{𝕍}}_{\mathit{n}}$. The influence of T and ${\mathsf{𝕍}}_{\mathit{n}}$ becomes increasingly important as these parameters approach their critical value. In particular, ${\rho }_{nn}$ and ${\rho }_{np}={\rho }_{pn}$ both vanish when ${\mathsf{𝕍}}_{\mathit{n}}={\mathbb{V}}_{cn}^{\left(0\right)}$: neutron superfluidity is destroyed and the neutron mass is thus entirely transported by the normal fluid, as can be seen from Equation (1). On the other hand, protons remain superconducting but are no longer entrained by neutrons: the two species are dynamically uncoupled. The proton mass current (2) thus reduces to the familiar expression for a single superfluid.

3.8. Chemical Potentials

Knowing the relation between the effective superfluid velocities ${\mathbb{V}}_{\mathit{q}}$ and the superfluid velocities ${\mathit{V}}_{\mathit{q}}$, given by Equations (85) and (86) respectively, and using Equation (71) for the potentials $U_q$ together with Equation (50) for the momentum densities ${\mathit{j}}_{\mathit{q}}$ and Equation (51) for the kinetic-energy densities ${\tau }_q$, we can compute the true chemical potentials ${\lambda }_q$ from Equations (48) and (49). Results for ${\mathsf{𝕍}}_{\mathit{p}}=0$ (as discussed in Section 3.6) are plotted in Figure 18 for $n=n_{\mathrm{cc}}$ and in Figure 19 for $n=n_0$ respectively. Although the chemical potentials ${\lambda }_q$ are found to be very weakly dependent on the temperature and on the neutron effective superfluid velocity (in the regime for which superfluidity exists), they deviate substantially from their corresponding Fermi energies ${\epsilon }_{Fq}$ due to the contribution from the potential $U_q$.

4. Conclusions

We have studied the neutron–proton superfluid mixture present in the outer core of a neutron star in the framework of the nuclear-energy-density functional theory. In particular, we have calculated consistently the ${}^1$S${}_0$ pairing gaps ${\Delta }_q$ of each nucleon species q, their chemical potentials ${\lambda }_q$, and the entrainment matrix elements ${\rho }_{q{q}^{\prime }}$ relating the mass currents ${\mathsf{\rho }}_{\mathit{q}}$ to the so-called superfluid “velocities” ${\mathit{V}}_{\mathit{q}}$ (actually representing superfluid momenta per unit mass) in the normal-fluid rest frame.

To this end, we have solved numerically the self-consistent TDHFB equations with the Brussels–Montreal functional BSk24 [18] for $npe\mu$ matter in beta-equilibrium over the whole range of temperatures and velocities for which nuclear superfluidity can exist using the composition published in [19]. We have considered the full TDHFB equations without any approximation. In particular, the vector potentials ${\mathit{I}}_{\mathit{q}}$ and the contributions from the momentum densities $\hslash {\mathit{j}}_{\mathit{q}}$ to the potentials $U_q$ and to the kinetic densities ${\tau }_q$ have been fully taken into account. We have shown that $\hslash {\mathit{j}}_{\mathit{q}}$ represents the total momentum density of a given nucleon species q, accounting not only for the superfluid momentum density ${\rho }_q{\mathit{V}}_{\mathit{q}}$ but also for the momentum density carried by the quasiparticles, as shown in Equation (54). Because the true velocity ${\mathit{v}}_{\mathit{q}}={\mathit{\rho }}_{\mathit{q}}/{\rho }_q$ of the nucleon species q in the normal frame is related to the effective superfluid velocity ${\mathsf{𝕍}}_{\mathit{q}}$ introduced in Equation (29) through Equation (53), ${\mathsf{𝕍}}_{\mathit{q}}$ appears as the natural variable to characterize the superflow of the nucleon species q.

The ${}^1$S${}_0$ proton pairing gaps ${\Delta }_p^{\left(0\right)}$ at zero temperature and in the absence of flows are found to be significantly smaller than the neutron gaps ${\Delta }_n^{\left(0\right)}$, unlike those generally considered in neutron-star simulations. Although proton gaps are mainly determined by the empirical interpolation (74) between the reference gaps in symmetric nuclear matter and pure neutron matter, they turn out to be consistent with recent diagrammatic calculations taking into account medium-polarization effects and considering both two- and three-body interactions [52]. We have shown that the gaps ${\Delta }_q^{\left(0\right)}$ are accurately reproduced by the approximate formula (77).

The normalized ${}^1$S${}_0$ pairing gaps ${\Delta }_q/{\Delta }_q^{\left(0\right)}$ and the fraction ${\mathcal{Y}}_q$ of quasiparticles are found to be universal functions of $T/T_{cq}^{\left(0\right)}$ and ${\mathsf{𝕍}}_{\mathit{q}}/{\mathbb{V}}_{cq}^{\left(0\right)}$, with the critical temperature and critical velocity given by Equations (78) and (79) respectively, in the sense that they depend neither on the composition nor on the density, and are the same for both neutrons and protons. This result can be understood from the fact that ${}^1$S${}_0$ nucleon superfluidity in the core of neutron stars is in the (weak-coupling) BCS regime. We have found that the temperature dependence of the normalized pairing gaps in the absence of flows is well fitted by Equation (80) proposed in [55]. We have obtained new accurate interpolating formulas for describing the velocity-dependence of the normalized pairing gaps at zero temperature (82), as well as for the critical temperatures (83). For arbitrary temperatures and velocities, the pairing gaps can be determined with a very good accuracy by solving numerically Equation (59) instead of the full TDHFB equations.

We have found that the approximations reducing the TDHFB equations to Landau’s theory provide accurate results for the entrainment matrix ${\rho }_{q{q}^{\prime }}$ provided the critical temperatures $T_{cq}^{\left(0\right)}$ in the absence of superflow are given. Moreover, the reduced chemical potentials ${\mu }_q$ defined by Equation (35) are well approximated by the corresponding Fermi energies ${\epsilon }_{Fq}$ given by Equation (41). However, this conclusion may change depending on the functional, especially if the adopted one predicts stronger pairing. Moreover, numerical solutions of the full TDHFB equations are still required for calculating the chemical potentials ${\lambda }_q$.

Together with the results published in [19,20,21,22] for the composition and the equation of state, our calculations provide consistent microscopic inputs for modeling the global structure and dynamics of superfluid neutron stars. Although we have considered the specific functional BSk24 because it has been shown to be in excellent agreement with existing nuclear data and astrophysical observations [23,44], we have also derived all the necessary equations to evaluate superfluid properties for any other functional of the Brussels–Montreal type (this includes all the functionals based on standard Skyrme effective nucleon-nucleon interactions). Extension of the TDHFB theory to account for ${}^3$PF${}_2$ neutron superfluidity in the inner core of massive neutron stars is left for future studies.