Entrainment effects in neutron-proton mixtures within the nuclear energy-density functional theory. II. Finite temperatures and arbitrary currents.

The TDHFB method is discussed, e.g. in the classical textbook of Ref. Blaizot and Ribka (1986). The energy $E$ of a nucleon-matter element of volume $V$ is expressed as a function of the one-body density matrix $n_{q}^{ij}$ and pairing tensor $\kappa_{q}^{ij}$ defined by the following thermal averages of the creation and destruction operators, $c_{q}^{i\dagger}$ and $c_{q}^{i}$ (using the symbol $\dagger$ for Hermitian conjugation), for nucleons of charge type $q$ in a quantum state $i$ (using the symbol $*$ for complex conjugation):

$$n_{q}^{ij}=\langle c_{q}^{j\dagger}c_{q}^{i}\rangle=n^{ji*}_{q}\,,$$

(4)

$$\kappa_{q}^{ij}=\langle c_{q}^{j}c_{q}^{i}\rangle=-\kappa_{q}^{ji}\,.$$

(5)

Introducing the generalized Bogoliubov transformation¹¹1In Ref. Blaizot and Ribka (1986), the matrices were denoted by $X_{j}^{i}=\mathcal{U}_{ij}$ and $Y_{j}^{i}=\mathcal{V}_{ij}$.

$\displaystyle\begin{pmatrix}b_{q}^{i}\\ b_{q}^{i\dagger}\end{pmatrix}=\sum_{j}\begin{pmatrix}\mathcal{U}^{(q)*}_{ij}&\mathcal{V}^{(q)*}_{ij}\\ \mathcal{V}^{(q)}_{ij}&\mathcal{U}^{(q)}_{ij}\end{pmatrix}\begin{pmatrix}c_{q}^{j}\\ c_{q}^{j\dagger}\end{pmatrix}\,,$

(6)

such that $\langle b_{q}^{j\dagger}b_{q}^{i}\rangle=\delta^{ij}f^{(q)}_{i}$ ($\delta^{ij}$ is the Kronecker’s symbol) and $\langle b_{q}^{j}b_{q}^{i}\rangle=\langle b_{q}^{j\dagger}b_{q}^{i\dagger}\rangle=0$, where $b_{q}^{i\dagger}$ and $b_{q}^{i}$ are creation and destruction operators of a quasiparticle in a quantum state $i$, the TDHFB equations, which formally take the same form at any temperature, can be written as Blaizot and Ribka (1986)

$\displaystyle i\hbar\frac{\partial\mathcal{U}^{(q)}_{ki}}{\partial t}=\sum_{j}\bigl{[}(h_{q}^{ij}-\lambda_{q}\delta^{ij})\mathcal{U}^{(q)}_{kj}+\Delta_{q}^{ij}\mathcal{V}^{(q)}_{kj}\bigr{]}\,,$

(7)

$\displaystyle i\hbar\frac{\partial\mathcal{V}^{(q)}_{ki}}{\partial t}=\sum_{j}\bigl{[}-\Delta_{q}^{ij*}\mathcal{U}^{(q)}_{kj}-(h_{q}^{ij*}-\lambda_{q}\delta^{ij})\mathcal{V}^{(q)}_{kj}\bigr{]}\,,$

(8)

where $\lambda_{q}$ denotes the chemical potentials. The matrices $h_{q}^{ij}$ and $\Delta_{q}^{ij}$ of the single-particle Hamiltonian and the pair potential, respectively, are defined as

$\displaystyle h_{q}^{ij}=\frac{\partial E}{\partial n_{q}^{ji}}=h_{q}^{ji*}\,,$

(9)

$\displaystyle\Delta_{q}^{ij}=\frac{\partial E}{\partial\kappa_{q}^{ij*}}=-\Delta_{q}^{ji}\,.$

(10)

The fermionic algebra of the particle operators ($c^{i}_{q}$ and $c_{q}^{i\dagger}$) and the quasiparticle operators ($b^{i}_{q}$ and $b_{q}^{i\dagger}$) yields the following identities

$\displaystyle\sum_{k}\left(\mathcal{U}^{(q)}_{ik}\mathcal{U}^{(q)*}_{jk}+\mathcal{V}^{(q)}_{ik}\mathcal{V}^{(q)*}_{jk}\right)=\delta_{ij}\,,\qquad\sum_{k}\left(\mathcal{U}^{(q)*}_{ki}\mathcal{U}^{(q)}_{kj}+\mathcal{V}^{(q)}_{ki}\mathcal{V}^{(q)*}_{kj}\right)=\delta_{ij},$

(11)

$\displaystyle\sum_{k}\left(\mathcal{U}^{(q)}_{ik}\mathcal{V}^{(q)}_{jk}+\mathcal{V}^{(q)}_{ik}\mathcal{U}^{(q)}_{jk}\right)=0\,,\qquad\sum_{k}\left(\mathcal{U}^{(q)*}_{ki}\mathcal{V}^{(q)}_{kj}+\mathcal{V}^{(q)}_{ki}\mathcal{U}^{(q)*}_{kj}\right)=0\,.$

(12)

The one-body density matrix and the pairing tensor can be expressed in terms of the quasiparticle components as

$\displaystyle n_{q}^{ij}=\sum_{k}\left[f^{(q)}_{k}\mathcal{U}^{(q)}_{ki}\mathcal{U}^{(q)*}_{kj}+(1-f^{(q)}_{k})\mathcal{V}^{(q)*}_{ki}\mathcal{V}^{(q)}_{kj}\right]\,,$

(13)

$\displaystyle\kappa_{q}^{ij}=\sum_{k}\left[f^{(q)}_{k}\mathcal{U}^{(q)}_{ki}\mathcal{V}^{(q)*}_{kj}+(1-f^{(q)}_{k})\mathcal{V}^{(q)*}_{ki}\mathcal{U}^{(q)}_{kj}\right]\,.$

(14)

The energy $E$ is generally further assumed to depend on $n_{q}^{ij}$ and $\kappa_{q}^{ij}$ only through the following local densities and currents²²2The energy may be a functional of other densities and currents. We consider here only those relevant for calculating the entrainment couplings in homogeneous nuclear matter using the most popular functionals.:

Since the energy $E$ is real, it can only depend on the pair density through its square modulus $|\widetilde{n}_{q}(\boldsymbol{r},t)|^{2}$. The pairing potential (28) can thus be written as

$$\Delta_{q}(\boldsymbol{r},t)=2\frac{\delta E}{\delta|\widetilde{n}_{q}(\boldsymbol{r},t)|^{2}}\widetilde{n}_{q}(\boldsymbol{r},t)=4\frac{\delta E}{\delta|\widetilde{n}_{q}(\boldsymbol{r},t)|^{2}}\Psi_{q}(\boldsymbol{r},t)\,.$$

(29)

The last equality shows that $\Delta_{q}(\boldsymbol{r},t)$ has the same phase as the order parameter $\Psi_{q}(\boldsymbol{r},t)$. Using Eq. (22), the matrix elements (25) of the pairing field (29) will thus generally take the following form:

$$\Delta_{q}^{ij}=\frac{1}{2}\sum_{k,l}v^{(q)}_{ijkl}\kappa^{kl}_{q}\,,$$

(30)

with

$$v^{(q)}_{ijkl}\equiv 4\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}\int{\rm d}^{3}\boldsymbol{r}\,\frac{\delta E}{\delta|\widetilde{n}_{q}(\boldsymbol{r},t)|^{2}}\varphi^{(q)}_{i}(\boldsymbol{r},\sigma)^{*}\varphi^{(q)}_{j}(\boldsymbol{r},-\sigma)^{*}\varphi^{(q)}_{k}(\boldsymbol{r},\sigma^{\prime})\varphi^{(q)}_{l}(\boldsymbol{r},-\sigma^{\prime})\,.$$

(31)

Let us note that $v^{(q)}_{ijkl}=v^{(q)*}_{klij}=-v^{(q)}_{jikl}=-v^{(q)}_{ijlk}$. In the traditional formulation of the TDHFB equations (see, e.g., Ref. Blaizot and Ribka (1986)), $v^{(q)}_{ijkl}$ represents the matrix elements of the effective (in-medium) two-body pairing interaction.

The TDHFB Eqs. (7) and (8) can be conveniently rewritten as Blaizot and Ribka (1986)

$\displaystyle i\hbar\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*}-\Delta_{q}^{ik}\kappa_{q}^{kj*}\right)\,,$

(32)

$\displaystyle i\hbar\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*}-\lambda_{q}\delta^{kj})-\Delta_{q}^{ik}n_{q}^{kj*}-n_{q}^{ik}\Delta_{q}^{kj}\right]+\Delta_{q}^{ij}\,.$

(33)

As shown in Appendix A, Eq. (32) can be translated in coordinate space by the same continuity equations as in the absence of pairing (a similar conclusion was previously drawn in Ref. Gusakov (2010) within Landau’s theory)

$\displaystyle\frac{\partial\rho_{q}(\boldsymbol{r},t)}{\partial t}+\boldsymbol{\nabla}\cdot\boldsymbol{\rho_{q}}(\boldsymbol{r},t)=0\,,$

(34)

where the nucleon mass current is given by Chamel and Allard (2019)

$\displaystyle\boldsymbol{\rho_{q}}(\boldsymbol{r},t)=\frac{m}{m_{q}^{\oplus}(\boldsymbol{r},t)}\hbar\boldsymbol{j_{q}}(\boldsymbol{r},t)+\rho_{q}(\boldsymbol{r},t)\frac{\boldsymbol{I_{q}}(\boldsymbol{r},t)}{\hbar}\,.$

(35)

The effective mass $m_{q}^{\oplus}(\boldsymbol{r},t)$ and the vector field $\boldsymbol{I_{q}}(\boldsymbol{r},t)$ are defined by Eq. (27). As shown in our previous work Chamel and Allard (2019), the nucleon mass current can be expressed solely in terms of the momentum densities as

	$\displaystyle\boldsymbol{\rho_{q}}(\boldsymbol{r},t)$	$\displaystyle=$	$\displaystyle\hbar\boldsymbol{j_{q}}(\boldsymbol{r},t)\Biggl{\{}1+\frac{2}{\hbar^{2}}\Biggl{[}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}(\boldsymbol{r},t)}-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}(\boldsymbol{r},t)}\Biggr{]}\rho(\boldsymbol{r},t)\Biggr{\}}$		(36)
			$\displaystyle-\hbar\boldsymbol{j}(\boldsymbol{r},t)\frac{2}{\hbar^{2}}\Biggl{[}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}(\boldsymbol{r},t)}-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}(\boldsymbol{r},t)}\Biggr{]}\rho_{q}(\boldsymbol{r},t)\,,$

where $E^{j}_{\rm nuc}$ represents the nuclear-energy terms contributing to the mass currents, and we have introduced the following fields

$$X_{0}(\boldsymbol{r},t)=n_{0}(\boldsymbol{r},t)\tau_{0}(\boldsymbol{r},t)-j_{0}(\boldsymbol{r},t)^{2}\,,$$

(37)

$$X_{1}(\boldsymbol{r},t)=n_{1}(\boldsymbol{r},t)\tau_{1}(\boldsymbol{r},t)-j_{1}(\boldsymbol{r},t)^{2}\,.$$

(38)

Let us recall that 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 so-called “superfluid velocity” of the nucleon species $q$ at position $\boldsymbol{r}$ and time $t$ is defined by (see, e.g., Eq. (2.4.21) of Ref. Leggett (2006))

$\displaystyle\boldsymbol{V_{q}}(\boldsymbol{r},t)=\frac{\hbar}{2m}\boldsymbol{\nabla}\phi_{q}(\boldsymbol{r},t)\,,$

(39)

where $\phi_{q}(\boldsymbol{r},t)$ is the phase of the associated condensate defined through the order parameter (23)

$\displaystyle\Psi_{q}(\boldsymbol{r},t)=|\Psi_{q}(\boldsymbol{r},t)|\exp(i\phi_{q}(\boldsymbol{r},t))\,.$

(40)

In this section, we consider an homogeneous neutron-proton mixture with stationary nucleon currents in the normal rest frame ($\boldsymbol{v_{N}}=\boldsymbol{0}$). The latter assumption ensures that the entropy densities $s_{q}$ are independent of time. Recalling that Blaizot and Ribka (1986)

$\displaystyle s_{q}=-\frac{k_{\rm B}}{V}\sum_{i}\left[f^{(q)}_{i}\ln f^{(q)}_{i}+(1-f^{(q)}_{i})\ln(1-f^{(q)}_{i})\right]\,,$

(41)

where $k_{\rm B}$ denotes Boltzmann’s constant, the quasiparticle distribution functions $f^{(q)}_{i}$ are therefore also independent of time.

The single-particle wave functions are given by plane waves:

$\displaystyle\varphi_{j}(\boldsymbol{r},\sigma)=\frac{1}{\sqrt{V}}\exp\left(i\boldsymbol{k}_{j}\cdot\boldsymbol{r}\right)\chi_{j}(\sigma)\,,$

(42)

where $\chi_{j}(\sigma)=\delta_{\sigma_{j}\sigma}$ denotes the Pauli spinor. As can be seen from Eqs. (21), and (24), the density matrix (4) and the single-particle Hamiltonian matrix (9) are both diagonal in this basis.

The superfluid velocities, which are necessarily spatially uniform and independent of time, are conveniently written as

$\displaystyle\boldsymbol{V_{q}}\equiv\frac{\hbar\boldsymbol{Q_{q}}}{m}\,.$

(43)

This corresponds to a superfluid phase $\phi_{q}(\boldsymbol{r})=2\,\boldsymbol{Q_{q}}\cdot\boldsymbol{r}$ (modulo some arbitrary constant term without any physical consequence). Inserting Eq. (42) in Eq. (23) using Eqs. (16) and (22), the order parameter (23) thus reduces to

$\displaystyle\Psi_{q}(\boldsymbol{r},t)=|\Psi_{q}(t)|\exp\left(2i\boldsymbol{Q_{q}}\cdot\boldsymbol{r}\right)=\frac{1}{4V}\sum_{i,j}\kappa_{q}^{ij}\exp\left[i\left(\boldsymbol{k}_{i}+\boldsymbol{k}_{j}\right)\cdot\boldsymbol{r}\right](\sigma_{j}-\sigma_{i})\,.$

(44)

One simple choice is to consider that $\kappa^{q}_{ij}$ is non-zero only if the states $i$ and $j$ have opposite spins and wave vectors $\boldsymbol{k}_{i}$ and $\boldsymbol{k}_{j}$ are such that $\boldsymbol{k}_{i}+\boldsymbol{k}_{j}=2\boldsymbol{Q_{q}}$. We can always arrange states such that $\boldsymbol{k}_{i}=\boldsymbol{k}+\boldsymbol{Q_{q}}$ and $\boldsymbol{k}_{j}=-\boldsymbol{k}+\boldsymbol{Q_{q}}$. For convenience, we introduce the following shorthand notation

$\displaystyle k\equiv(\boldsymbol{k}+\boldsymbol{Q_{q}},\sigma)\,,\qquad\qquad\bar{k}\equiv(-\boldsymbol{k}+\boldsymbol{Q_{q}},-\sigma)\,.$

(45)

These quantum numbers define the conjugate states that are paired. Indeed, it follows from the definition (10) that the only nonzero matrix elements of the pair potential are of the form $\Delta_{q}^{k\bar{k}}=-\Delta_{q}^{\bar{k}k}$. In the absence of current ($\boldsymbol{Q_{q}}=\boldsymbol{0}$), the conjugate state $\bar{k}$ is the time-reversed state of $k$. The presence of a non-vanishing current ($\boldsymbol{Q_{q}}\neq\boldsymbol{0}$) breaks the time-reversal symmetry. The nonzero elements of the $\mathcal{U}^{(q)}$ and $\mathcal{V}^{(q)}$ matrices satisfying Eqs. (11) and (12) are of the form $\mathcal{U}^{(q)}_{kk}=\mathcal{U}^{(q)}_{\bar{k}\bar{k}}$ and $\mathcal{V}^{(q)}_{k\bar{k}}=-\mathcal{V}^{(q)}_{\bar{k}k}$. Substituting $i\hbar\partial/\partial t$ by the quasiparticle energies $\mathfrak{E}^{(q)}_{k}$, the TDHFB Eqs. (7) and (8) finally reduce to

$\displaystyle\begin{pmatrix}\xi^{(q)}_{k}&\Delta^{(q)}_{k}\\ \Delta^{(q)*}_{k}&-\xi^{(q)*}_{\bar{k}}\end{pmatrix}\begin{pmatrix}\mathcal{U}^{(q)}_{kk}\\ \mathcal{V}^{(q)}_{k\bar{k}}\end{pmatrix}=\mathfrak{E}^{(q)}_{k}\begin{pmatrix}\mathcal{U}^{(q)}_{kk}\\ \mathcal{V}^{(q)}_{k\bar{k}}\end{pmatrix}\,,$

(46)

where $\xi^{(q)}_{k}\equiv\epsilon^{(q)}_{k}-\lambda_{q}$ ($\epsilon^{(q)}_{k}$ denoting the eigenvalues of the single-particle Hamiltonian matrix) are explicitly given by

$\displaystyle\xi^{(q)}_{k}=\frac{\hbar^{2}}{2m_{q}^{\oplus}}\left(\boldsymbol{k}+\boldsymbol{Q_{q}}\right)^{2}+U_{q}+\boldsymbol{I}_{q}\cdot\left(\boldsymbol{k}+\boldsymbol{Q_{q}}\right)-\lambda_{q}\,.$

(47)

The matrix elements $\Delta^{(q)}_{k}\equiv\Delta_{q}^{k\bar{k}}=-\Delta^{(q)}_{\bar{k}}$ are given by the following equation

$$\Delta^{(q)}_{k}=\frac{1}{2}\sum_{l}v^{(q)}_{k\bar{k}l\bar{l}}\,\kappa_{q}^{l\bar{l}}\,.$$

(48)

Using Eqs. (13), (15), (18), and (21), the nucleon number and momentum densities read

$\displaystyle n_{q}=\frac{1}{V}\sum_{k}n_{q}^{kk}=\frac{1}{V}\sum_{k}\left[|\mathcal{U}^{(q)}_{kk}|^{2}f^{(q)}_{k}+|\mathcal{V}^{(q)}_{\bar{k}k}|^{2}\left(1-f^{(q)}_{\bar{k}}\right)\right]\,,$

(49)

$\displaystyle\boldsymbol{j_{q}}=\frac{1}{V}\sum_{k}\boldsymbol{k}_{k}n_{q}^{kk}=\frac{1}{V}\sum_{k}\left(\boldsymbol{k}+\boldsymbol{Q_{q}}\right)\left[|\mathcal{U}^{(q)}_{kk}|^{2}f^{(q)}_{k}+|\mathcal{V}^{(q)}_{\bar{k}k}|^{2}\left(1-f^{(q)}_{\bar{k}}\right)\right]\,,$

(50)

respectively, where the quasiparticle distribution is given by Blaizot and Ribka (1986)

$$f^{(q)}_{k}=\left[1+\exp\left(\beta\mathfrak{E}^{(q)}_{k}\right)\right]^{-1}=\frac{1}{2}\left[1-\tanh\left(\frac{\beta}{2}\mathfrak{E}^{(q)}_{k}\right)\right]\,.$$

(51)

where $\beta\equiv\left(k_{\textrm{B}}T\right)^{-1}$.

The solutions of Eq. (46), readily obtained by diagonalizing the HFB matrix, are given by³³3Equation (46) actually admits two kinds of solutions corresponding to the eigenvalues $\mathfrak{E}^{(q)}_{k\pm}=(\xi^{(q)}_{k}-\xi^{(q)}_{\bar{k}})/2\pm\sqrt{\varepsilon^{(q)2}_{k}+|\Delta^{(q)}_{k}|^{2}}$. Those associated with $\mathfrak{E}^{(q)}_{k-}$ are such that the expressions (53) and (54) of $|\mathcal{U}^{(q)}_{kk}|^{2}$ and $|\mathcal{V}^{(q)}_{k\bar{k}}|^{2}$ are swapped. However, these solutions lead to the same values for the nucleon number density (49) and the momentum density (50) using the fact that $\mathfrak{E}^{(q)}_{k-}=-\mathfrak{E}^{(q)}_{\bar{k}+}$. Therefore, they give the same mass current (36) and, consequently, the same entrainment matrix.

$\displaystyle\mathfrak{E}^{(q)}_{k}=\frac{\xi^{(q)}_{k}-\xi^{(q)}_{\bar{k}}}{2}+\sqrt{\varepsilon^{(q)2}_{k}+|\Delta^{(q)}_{k}|^{2}}\,,$

(52)

$\displaystyle|\mathcal{U}^{(q)}_{kk}|^{2}=\frac{1}{2}\left(1+\frac{\varepsilon^{(q)}_{k}}{\sqrt{\varepsilon^{(q)2}_{k}+|\Delta^{(q)}_{k}|^{2}}}\right)\,,$

(53)

$\displaystyle|\mathcal{V}^{(q)}_{k\bar{k}}|^{2}=\frac{1}{2}\left(1-\frac{\varepsilon^{(q)}_{k}}{\sqrt{\varepsilon^{(q)2}_{k}+|\Delta^{(q)}_{k}|^{2}}}\right)\,,$

(54)

where

$$\varepsilon^{(q)}_{k}\equiv\frac{\xi^{(q)}_{k}+\xi^{(q)}_{\bar{k}}}{2}=\varepsilon^{(q)}_{\bar{k}}\,.$$

(55)

Using Eq. (47), we have

$$\varepsilon^{(q)}_{k}=\frac{\hbar^{2}}{2m_{q}^{\oplus}}\left(\boldsymbol{k}^{2}+\boldsymbol{Q_{q}}^{2}\right)+U_{q}+\boldsymbol{I}_{q}\cdot\boldsymbol{Q_{q}}-\lambda_{q}\,,$$

(56)

$$\frac{\xi^{(q)}_{k}-\xi^{(q)}_{\bar{k}}}{2}=\hbar\boldsymbol{k}\cdot\boldsymbol{\mathbb{V}_{q}}\,,$$

(57)

where we have introduced the effective superfluid velocity

$$\boldsymbol{\mathbb{V}_{q}}=\frac{\hbar}{m_{q}^{\oplus}}\boldsymbol{Q_{q}}+\frac{\boldsymbol{I_{q}}}{\hbar}=\frac{m}{m_{q}^{\oplus}}\boldsymbol{V_{q}}+\frac{\boldsymbol{I_{q}}}{\hbar}\,.$$

(58)

For standard Skyrme functionals, the effective mass $m_{q}^{\oplus}$ depends only on the nucleon densities, whereas the potential $U_{q}$ depends on the pairing gaps $\Delta^{(q)}_{k}$ (which in turn depend also on $\boldsymbol{Q_{n}}$ and $\boldsymbol{Q_{p}}$). For the extended Skyrme functionals proposed in Ref. Chamel et al. (2009), the potential $U_{q}$ depends also explicitly on $j_{n}^{2}$, $j_{p}^{2}$ and $\boldsymbol{j_{n}}\cdot\boldsymbol{j_{p}}$. In general, the dependence of the energies (56) on the superfluid velocities may be therefore quite complicated.

Using Eqs. (31) and (42), it can be shown that the effective pairing interaction $v^{(q)}_{k\bar{k}l\bar{l}}$ is independent of the wave vectors, and depends only on the spins. The only nonzero elements are given by (denoting the spins with arrows for clarity)

$$v^{(q)}_{\downarrow\uparrow\downarrow\uparrow}=\frac{4}{V}\frac{\delta E}{\delta|\widetilde{n}_{q}|^{2}}<0\,,$$

(59)

and any other element obtained by permutation of the spin indices. It thus follows from Eq. (48), that $\Delta^{(q)}_{k}$ is also independent of the wave vectors $\boldsymbol{k}$ (but $\Delta^{(q)}_{k}$ still depends on the given wave vectors $\boldsymbol{Q_{q}}$). Dropping the wave vector $\boldsymbol{k}$ as a subscript, introducing the pairing gap

$$\Delta_{q}\equiv\Delta^{(q)}_{\downarrow\uparrow}=\frac{1}{V}\int d^{3}\boldsymbol{r}\,|\Delta_{q}(\boldsymbol{r})|\geq 0\,,$$

(60)

and using Eqs. (53) and (54), Eq. (48) reduces to the gap equation

$$\Delta_{q}=\frac{2}{V}\frac{\delta E}{\delta|\widetilde{n}_{q}|^{2}}\sum_{\boldsymbol{k}}\frac{\Delta_{q}}{\sqrt{\varepsilon_{k}^{(q)2}+\Delta_{q}^{2}}}(f_{k}+f_{\bar{k}}-1)\,.$$

(61)

The summation is only over the wave vectors $\boldsymbol{k}$, the summation over the spins has been already carried out. Note that the right-hand side of this equation explicitly depends on the wave vectors $\boldsymbol{Q_{q}}$ through Eqs. (51), (52), (56) and (57).

Using the solution of the TDHFB equations, the momentum density (50) can be alternatively written as

$\displaystyle\boldsymbol{j_{q}}=n_{q}\,\boldsymbol{Q_{q}}+\frac{1}{V}\sum_{\boldsymbol{k}}\boldsymbol{k}(f^{(q)}_{k}-f^{(q)}_{\bar{k}})\,.$

(62)

The first term in the right-hand side coincides with the expression (29) adopted in our previous work Chamel and Allard (2019), thus demonstrating the validity of this identification. The second term accounts for quasiparticle excitations (note that the summation over spin states has been already carried out). Remarking that the component of $\boldsymbol{k}$ orthogonal to $\boldsymbol{\mathbb{V}_{q}}$ does not contribute to the sum (since $f^{(q)}_{k}=f^{(q)}_{\bar{k}}$ in this case), the momentum density can be expressed as

$$\boldsymbol{j_{q}}=\frac{mn_{q}}{\hbar}\boldsymbol{V_{q}}-\frac{m_{q}^{\oplus}n_{q}}{\hbar}\mathcal{Y}_{q}\boldsymbol{\mathbb{V}_{q}}\,,$$

(63)

with the $\mathcal{Y}_{q}$ function defined as

$\displaystyle\mathcal{Y}_{q}(T,\boldsymbol{\mathbb{V}_{q}})\equiv-\frac{\hbar}{m_{q}^{\oplus}n_{q}\mathbb{V}_{q}^{2}}\frac{1}{V}\sum_{\boldsymbol{k}}\boldsymbol{k}\cdot\boldsymbol{\mathbb{V}_{q}}(f^{(q)}_{k}-f^{(q)}_{\bar{k}})\,.$

(64)

Substituting Eq. (51) yields

$\displaystyle\mathcal{Y}_{q}(T,\boldsymbol{\mathbb{V}_{q}})=\frac{\hbar}{m_{q}^{\oplus}n_{q}\mathbb{V}_{q}^{2}}\left[\frac{1}{V}\sum_{\boldsymbol{k}}\frac{\boldsymbol{k}\cdot\boldsymbol{\mathbb{V}_{q}}\sinh\left(\beta\hbar\boldsymbol{k}\cdot\boldsymbol{\mathbb{V}_{q}}\right)}{\cosh\left(\beta\mathfrak{E}^{(q)}_{\boldsymbol{k}}\right)+\cosh\left(\beta\hbar\boldsymbol{k}\cdot\boldsymbol{\mathbb{V}_{q}}\right)}\right]\,.$

(65)

In terms of the effective superfluid velocity (58), the mass current (35) reduces to

$\displaystyle\boldsymbol{\rho_{q}}=\rho_{q}\left(1-\mathcal{Y}_{q}\right)\boldsymbol{\mathbb{V}_{q}}\,,$

(66)

while the momentum density $\boldsymbol{j_{q}}$ becomes

$\displaystyle\boldsymbol{j_{q}}=\frac{\rho_{q}}{\hbar}\left(1-\mathcal{Y}_{q}\right)\boldsymbol{V_{q}}-\frac{m_{q}^{\oplus}n_{q}}{\hbar^{2}}\mathcal{Y}_{q}\boldsymbol{I_{q}}\,.$

(67)

Using this expression of $\boldsymbol{j_{q}}$ combined with the general expression of the vector field $\boldsymbol{I_{q}}$ Chamel and Allard (2019)

$\displaystyle\boldsymbol{I_{q}}=\frac{\delta E^{j}_{\rm nuc}}{\delta\boldsymbol{j_{q}}}=-2\left(\boldsymbol{j_{p}}+\boldsymbol{j_{n}}\right)\left(\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}\right)-4\boldsymbol{j_{q}}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}$

(68)

leads to a self-consistent system of equations for $\boldsymbol{j_{q}}$ and $\boldsymbol{V_{q}}$, whose solutions are

$\displaystyle\boldsymbol{I_{q}}=\sum_{q^{\prime}}\mathcal{I}_{qq^{\prime}}\boldsymbol{V_{q^{\prime}}}\,,$

(69)

$\displaystyle\mathcal{I}_{nn}$

$\displaystyle=\frac{2}{\hbar}\rho_{n}\left(1-\mathcal{Y}_{n}\right)\Theta\left[\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}\left(\frac{8}{\hbar^{2}}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}m_{p}^{\oplus}n_{p}\mathcal{Y}_{p}-1\right)-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}\right]\,,$

(70)

$\displaystyle\mathcal{I}_{pp}$

$\displaystyle=\frac{2}{\hbar}\rho_{p}\left(1-\mathcal{Y}_{p}\right)\Theta\left[\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}\left(\frac{8}{\hbar^{2}}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}m_{n}^{\oplus}n_{n}\mathcal{Y}_{n}-1\right)-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}\right]\,,$

(71)

$\displaystyle\mathcal{I}_{np}$

$\displaystyle=\frac{2}{\hbar}\rho_{p}\left(1-\mathcal{Y}_{p}\right)\Theta\left(\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}\right)\,,$

(72)

$\displaystyle\mathcal{I}_{pn}$

$\displaystyle=\frac{2}{\hbar}\rho_{n}\left(1-\mathcal{Y}_{n}\right)\Theta\left(\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}-\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}\right)\,,$

(73)

	$\displaystyle\Theta$	$\displaystyle\equiv\left[1-\frac{2}{\hbar^{2}}\left(\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}+\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}\right)\left(m_{n}^{\oplus}n_{n}\mathcal{Y}_{n}+m_{p}^{\oplus}n_{p}\mathcal{Y}_{p}\right)\right.$
		$\displaystyle\qquad\;\;\left.+\left(\frac{4}{\hbar^{2}}\right)^{2}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{0}}\frac{\delta E^{j}_{\rm nuc}}{\delta X_{1}}m_{n}^{\oplus}n_{n}m_{p}^{\oplus}n_{p}\mathcal{Y}_{n}\mathcal{Y}_{p}\right]^{-1}\,.$		(74)

Substituting Eq. (69) into (58), the entrainment matrix can be readily obtained from Eq. (66):

$\displaystyle\rho_{nn}(T,\boldsymbol{\mathbb{V}_{q}})=\rho_{n}\left(1-\mathcal{Y}_{n}\right)\left(\frac{m}{m_{n}^{\oplus}}+\frac{\mathcal{I}_{nn}}{\hbar}\right)\,,$

(75)

$\displaystyle\rho_{pp}(T,\boldsymbol{\mathbb{V}_{q}})=\rho_{p}\left(1-\mathcal{Y}_{p}\right)\left(\frac{m}{m_{p}^{\oplus}}+\frac{\mathcal{I}_{pp}}{\hbar}\right)\,,$

(76)

$\displaystyle\rho_{np}(T,\boldsymbol{\mathbb{V}_{q}})=\rho_{n}\left(1-\mathcal{Y}_{n}\right)\frac{\mathcal{I}_{np}}{\hbar}\,,\qquad\qquad\rho_{pn}(T,\boldsymbol{\mathbb{V}_{q}})=\rho_{p}\left(1-\mathcal{Y}_{p}\right)\frac{\mathcal{I}_{pn}}{\hbar}\,.$

(77)

Using Eqs. (72) and (73), one can notice that the entrainment matrix is manifestly symmetric (i.e. $\rho_{np}=\rho_{pn}$). Let us emphasize that Eqs. (75)-(77) give the exact expression for the entrainment matrix within the TDHFB theory in homogeneous nuclear matter. No approximation has been made at this point. In particular, the full dependence of the entrainment matrix elements on the currents is taken into account.