Nonlinear approach to the entrainment matrix of superfluid nucleon mixture at zero temperature

TL;DR

This paper derives nonlinear equations for the entrainment matrix in superfluid neutron-proton mixtures, revealing significant velocity-dependent effects that impact neutron star dynamics.

Contribution

It introduces a nonlinear approach to calculating the entrainment matrix considering velocity-dependent energy gaps, a novel aspect in superfluid hydrodynamics.

Findings

Entrainment matrix components depend nonlinearly on superfluid velocities.

Velocity dependence of energy gaps significantly affects neutron star models.

Numerical solutions show strong nonlinear behavior in typical cases.

Abstract

The superfluid drag effect, in hydrodynamics of pulsating neutron stars, is conventionally described with the aid of the entrainment matrix relating the mass currents with the velocities of superfluid flows in the system. Equations for the entrainment matrix of a superfluid mixture of neutrons and protons are derived with allowance for the strong dependence of the energy gaps on the velocities of superfluid flows. The calculations are carried out in the frame of the Fermi-liquid theory. The equations obtained are highly nonlinear. Numerical solutions to the equations for some typical cases demonstrate that the components of the entrainment matrix possess a highly nonlinear dependence on the velocities of the two superflows simultaneously. This effect, previously ignored, can greatly influence the dynamics of neutron stars.

Figures

Equations

$${\mathbf{j}}_{\mathsf{n}}$$

$${\mathbf{j}}_{\mathsf{p}}$$

$${\mathbf{j}}_{\alpha }=\sum _{\mathbf{k},\sigma }{m}_{\alpha }\frac{\partial {\tilde{\epsilon }}_{\mathbf{k}}^{(\alpha )}}{\partial \mathbf{k}}{\tilde{n}}_{\mathbf{k}}^{(\alpha )},$$

$$v_{F\alpha }=\frac{d{ϵ}_k^{(\alpha )}}{dk}=\frac{k}{m_{\alpha }^{\ast }},$$

$$p_{F\alpha }={\left(3{\pi }^2{n}_{\alpha }\right)}^{1/3}.$$

$${ϵ}^{(\alpha )}\left(k\right)={\mu }_{\alpha }+{\xi }_{\mathbf{k}}^{(\alpha )},$$

$${\xi }_k^{(\alpha )}=v_{F\alpha }\left(k-p_{F\alpha }\right).$$

$${\epsilon }_{\mathbf{k}}^{(\alpha )}={\xi }_k^{(\alpha )}+\sum _{{\mathbf{k}}^{\prime }{\sigma }^{\prime }\beta }{f}^{\alpha \beta }\left({\mathbf{k},\mathbf{k}}^{\prime }\right)\left(n_{{\mathbf{k}}^{\prime }}^{(\beta )}-{\Theta }_{{\mathbf{k}}^{\prime }}^{(\beta )}\right),$$

$$n_{\mathbf{k}}^{(\alpha )}\equiv ⟨\left|a_{\mathbf{k}\sigma }^{(\alpha )†}{a}_{\mathbf{k}\sigma }^{(\alpha )}\right|⟩$$

$${\tilde{\epsilon }}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}={\xi }_k^{(\alpha )}+\frac{{\mathbf{kq}}_{\alpha }}{m_{\alpha }^{\ast }}+\sum _{{\mathbf{k}}^{\prime }{\sigma }^{\prime }\beta }{f}^{\alpha \beta }\left({\mathbf{k},\mathbf{k}}^{\prime }\right)\left({\tilde{n}}_{{\mathbf{k}}^{\prime }{+\mathbf{q}}_{\beta }}^{(\beta )}-{\Theta }_{{\mathbf{k}}^{\prime }{+\mathbf{q}}_{\beta }}^{(\beta )}\right),$$

$$\delta {\epsilon }_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}=\frac{{\mathbf{kq}}_{\alpha }}{m_{\alpha }^{\ast }}+\sum _{{\mathbf{k}}^{\prime }\beta {\sigma }^{\prime }}{f}^{\alpha \beta }\left({\mathbf{k},\mathbf{k}}^{\prime }\right)\delta n_{{\mathbf{k}}^{\prime }{+\mathbf{q}}_{\beta }}^{(\beta )},$$

$$\delta n_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}\equiv \left({\tilde{n}}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}-{\Theta }_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}\right)-\left(n_{\mathbf{k}}^{(\alpha )}-{\Theta }_{\mathbf{k}}^{(\alpha )}\right).$$

$${\tilde{n}}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}\equiv ⟨\left|a_{{\mathbf{k}+\mathbf{q}}_{\alpha }\sigma }^{(\alpha )†}{a}_{{\mathbf{k}+\mathbf{q}}_{\alpha }\sigma }^{(\alpha )}\right|⟩,$$

$${\tilde{n}}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}={\tilde{u}}_{\mathbf{k}}^{(\alpha )2}{\mathcal{F}}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}+{\tilde{v}}_{\mathbf{k}}^{(\alpha )2}\left(1-{\mathcal{F}}_{-{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}\right),$$

$${\mathcal{F}}_{\pm {\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}=\frac{1}{1+e^{\frac{1}{T}{\mathfrak{E}}_{\pm {\mathbf{k}+\mathbf{q}}_{\alpha }}^{\left(\alpha \right)}}}$$

$${\tilde{u}}_{\mathbf{k}}^{(\alpha )2}=\frac{1}{2}\left(1+\frac{{\xi }_{\mathbf{k}}^{(\alpha )}}{{\tilde{E}}_k^{\left(\alpha \right)}}\right),$$

$${\tilde{v}}_{\mathbf{k}}^{(\alpha )2}=\frac{1}{2}\left(1-\frac{{\xi }_{\mathbf{k}}^{(\alpha )}}{{\tilde{E}}_k^{\left(\alpha \right)}}\right).$$

$${\tilde{E}}_k^{\left(\alpha \right)}\equiv \sqrt{{\xi }_{\mathbf{k}}^{(\alpha )2}+{\tilde{\Delta }}^{(\alpha )2}},$$

$$n_{\mathbf{k}}^{(\alpha )}=\frac{1}{2}\left(1-\frac{{\xi }_{\mathbf{k}}^{(\alpha )}}{E_k^{\left(\alpha \right)}}\right),$$

$$E_k^{\left(\alpha \right)}\equiv \sqrt{{\xi }_{\mathbf{k}}^{(\alpha )2}+{\Delta }_0^{(\alpha )2}}$$

$$\delta {\epsilon }_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}=\sum _{\beta }{\gamma }_{\alpha \beta }\left({\mathbf{v}}_{\mathsf{p}s},{\mathbf{v}}_{\mathsf{n}s}\right){\mathbf{kv}}_{\beta s},$$

$${\tilde{\epsilon }}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}={\xi }_k^{(\alpha )}+\sum _{\beta }{\gamma }_{\alpha \beta }\left({\mathbf{v}}_{\mathsf{p}s},{\mathbf{v}}_{\mathsf{n}s}\right){\mathbf{kv}}_{\beta s},$$

$$\delta n_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}$$

$$+\frac{1}{2}\left({\mathcal{F}}_{{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}-{\mathcal{F}}_{-{\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}\right)+{\mathbf{q}}_{\alpha }\hat{\mathbf{k}}\delta \left(k-p_{F\alpha }\right),$$

$${\mathfrak{E}}_{\pm {\mathbf{k}+\mathbf{q}}_{\alpha }}^{\left(\alpha \right)}={\tilde{E}}_k^{\left(\alpha \right)}\pm \sum _{\beta }{\gamma }_{\alpha \beta }{\mathbf{kv}}_{\beta s}.$$

$${\tilde{E}}_k^{\left(\alpha \right)}\pm \sum _{\beta }{\gamma }_{\alpha \beta }{\mathbf{kv}}_{\beta s}\le 0.$$

$${\mathcal{F}}_{\pm {\mathbf{k}+\mathbf{q}}_{\alpha }}^{(\alpha )}=\Theta \left(\mp \sum _{\beta }{\gamma }_{\alpha \beta }{\mathbf{kv}}_{\beta s}-{\tilde{E}}_k^{\left(\alpha \right)}\right),$$

$$f^{\alpha \beta }\left({\mathbf{k},\mathbf{k}}^{\prime }\right)=\sum _l{f}_l^{\alpha \beta }{P}_l\left(\cos \theta \right),\cos \theta \equiv {\hat{\mathbf{k}}\hat{\mathbf{k}}}^{\prime }.$$

$$F_1^{\alpha \beta }\equiv f_1^{\alpha \beta }\sqrt{N_{0\alpha }{N}_{0\beta }},$$

$${\mathbf{V}}_{\alpha }\left({\mathbf{v}}_{\mathsf{p}s},{\mathbf{v}}_{\mathsf{n}s}\right)\equiv \sum _{\beta }{\gamma }_{\alpha \beta }\left({\mathbf{v}}_{\mathsf{p}s},{\mathbf{v}}_{\mathsf{n}s}\right){\mathbf{v}}_{\beta s}.$$

$$\sum _{\beta }{\gamma }_{\alpha \beta }{\mathbf{kv}}_{\beta s}$$

Full text

††pubyear: 2017††pagerange: Nonlinear approach to the entrainment matrix of superfluid nucleon mixture at zero temperature–LABEL:lastpage

1 Introduction

Pulsations of neutron stars consisting mostly of superfluid nucleons are of a great interest today (see Andersson et al., 1999; Andersson, 1998; Friedman & Morsink, 1998; Andersson et al., 2003; Andersson & Kokkotas, 2001; Gusakov et al., 2014; Arras et al., 2003; Gusakov et al., 2005; Sidery et al., 2010; Prix & Rieutord, 2002; Glampedakis & Jones, 2014; Sourie, Oertel & Novak, 2016; Alford & Pangeni, 2017). As is well known the theory of this phenomenon involves the drag effects leading to the mutual entrainment of superfluid flows of neutrons and protons. This effect is conventionally described with the aid of the so called "entrainment matrix" that relates the mass current density with a relative velocity of superfluid and normal components in the nucleon liquid. For the simplest case of the npe-matter the symmetric entrainment matrix $\rho_{\alpha\beta}$ $\left(\alpha,\beta=\mathsf{n,p}\right)$ can be defined by the Andreev-Bashkin relation (Andreev & Bashkin, 1975):

[TABLE]

where $\rho_{\alpha}=m_{\alpha}n_{\alpha}$. Hereafter $m_{\alpha}$ and $n_{\alpha}$ are the nucleon mass and number density of nucleon species $\alpha$, respectively; $\mathbf{j}_{\alpha}$ and $\mathbf{v}_{\alpha s}$ denote the mass current densities and the velocities of neutron and proton superfluid components. Finally, $\mathbf{v}_{\alpha\mathsf{q}}$ represent velocities of normal (non-superfluid) nucleons.111The normal component vanishes at zero temperature in a superfluid at rest. However, in superfluid flows, the normal component appears due to the superfluid motion. For definiteness we consider the mixture of superfluids in a coordinate frame, where $\mathbf{v}_{\mathsf{q}}=0$.

All the existing calculations of the entrainment matrix(see, e.g. Borumand et al., 1996; Chamel & Haensel, 2006; Gusakov et al., 2009a) do not take into consideration the change of the energy gap $\Delta^{\left(\alpha\right)}\left(\mathbf{v}_{\alpha s}\right)$ due to the motion of superfluid condensates and thus are justified only for small oscillation amplitudes, which are restricted by the condition $v_{\alpha s}\ll v_{cr}^{\left(\alpha\right)}$, where $v_{cr}^{\left(\alpha\right)}\sim\Delta_{0}^{\left(\alpha\right)}/p_{F\alpha}$ is the critical speed of the superfluid flow at which the superfluidity is destroyed ($\Delta_{0}^{\left(\alpha\right)}$ is the energy gap for Bogoliubov excitations in the superfluid at rest). A comparison of the critical velocity, which is necessary for the energy gap collapse with the known estimates of relative velocities between the superfluid and normal components inside oscillating neutron stars (Kantor & Gusakov, 2014; Gusakov & Kantor, 2013; Gusakov et al., 2005) clearly demonstrates that the speed of superfluid motion can exceed the critical values, which means that the considered gap suppression effect can be very important for the dynamics of superfluid pulsations. Our goal is to develop an approach that allows one to evaluate the entrainment matrix at arbitrary speed of the superflows in the npe-matter.

The suggested theory assumes that protons are paired in the spin-singlet 1S0 state while the neutron pairing occurs into the spin-triplet 3P2 state (Tamagaki, 1970; Takatsuka, 1972). We assume also that, in the npe matter, the coulomb interaction of protons is completely screened by relativistic degenerate electrons.

The paper is organized as follows. In Sec. 2 the problem is considered in the frame of the Fermi liquid theory, where the Fermi-liquid interactions are expressed through the Landau parameters. We derive the nonlinear equations for the entrainment matrix in a mixture of superfluid Fermi liquids, which are valid for any velocities of the superfluid flows. The equations obtained must be solved simultaneously with the equations for velocity-dependent energy gaps in the superfluid flows. In Sec. 3 we solve the corresponding gap equations for protons and neutrons in the case of zero temperature. In Sec. 4 we calculate the functions describing the balance between the superfluid and normal components in the superfluid flows. In Sec. 5 we incorporate the resulting solution into equations for the entrainment matrix. Section 6 contains a derivation of the expression for the superfluid mass currents in the system under consideration. In Sec. 7 we discuss possible applications of the nonlinear equations. We demonstrate the results of their numerical solution for some interesting cases and compare them with known results obtained in the linear approximation. The summary and conclusion are collected in Sec. 8. The appendix contains the derivation of the equation for the 3P2 energy gap in the moving superfluid flow.

We use the system of units in which $c=\hbar=1$, and the Boltzmann constant $k_{\mathrm{B}}=1$.

2 Nonlinear equations for the entrainment matrix

We begin with a general expression for the mass current density $\mathbf{j}_{\alpha}$ in a superfluid mixture of two species of baryons. For concreteness, we shall consider protons and neutrons, $\alpha=\mathsf{n,p}$, although the same theory can be applied in the case of superfluid hyperons and nucleons in the neutron star core. Following the Fermi-liquid theory, the mass current density can be evaluated from the same expression, which is normally used in the case of non-superfluid matter (Leggett, 1965, 1975)

[TABLE]

where $\mathbf{k}$ and $\sigma=\uparrow,\downarrow$ are the momentum of the quasiparticle and its spin projection onto the quantization axis, respectively. In this expression, the quasiparticle energy $\tilde{\varepsilon}_{\mathbf{k}}^{(\alpha)}$ depends on the distribution $\tilde{n}_{\mathbf{k}}^{(\alpha)}$ of all interacting particles in the system in the presence of superfluid motion. Here and below, the tilde above a letter indicates the values that depend on the velocity of the superfluid flow.

For a system at rest, the energy of protons and neutrons in a normal Fermi liquid $\epsilon_{k}^{(\alpha)}$ is renormalized by interactions only. According to the Landau theory the renormalized energy of a quasiparticle may be written with the aid of the effective mass $m_{\alpha}^{\ast}$, which is no longer equal to the bare mass $m_{\alpha}$ but is defined by

[TABLE]

where the derivative is evaluated at the appropriate Fermi surface given by the Fermi wave vector

[TABLE]

Since the nucleon matter is highly degenerate it is convenient to write the energy of a quasiparticle $\alpha$ as

[TABLE]

where $\mu_{\alpha}=\epsilon^{(\alpha)}\left(p_{F\alpha}\right)$ is the chemical potential (the Fermi energy) of a quasiparticle species $\alpha$, and $\xi_{k}^{(\alpha)}\equiv\epsilon^{(\alpha)}\left(k\right)-\mu_{\alpha}$. Near the Fermi surface one can approximate

[TABLE]

In the mixture of superfluid liquids at rest, the quasiparticle energy is renormalized as222This correction $\sim(\Delta/\mu)^{2}$ is small and normally ignored.

[TABLE]

where $\Theta_{\mathbf{k}}^{\left(\alpha\right)}\equiv\Theta\left(p_{F\alpha}-\left|\mathbf{k}\right|\right)$ is the step function, and

[TABLE]

is the distribution function of quasiparticles in the mixture of superfluid liquids at rest. For simplicity, we consider unitary states of the superfluid nucleons of both species. Since the uniform motion of a unitary superfluid liquid induces no spin polarization we shall use the spin-averaged Fermi-liquid interactions, which are parametrized by the functions $f^{\alpha\beta}\left(\mathbf{k,k}^{\prime}\right)$.

The supercurrent in the system arises when the Cooper condensate is formed by pairing of quasiparticles with momentum $(\mathbf{k+q}_{\alpha})$ and $(-\mathbf{k+q}_{\alpha})$. In this case, each pair acquires a nonzero momentum $2\mathbf{q}_{\alpha}$, and the condensate moves, as a whole, with a velocity $\mathbf{v}_{\mathsf{\alpha}s}=\mathbf{q}_{\alpha}/m_{\alpha}$ relative to non-superfluid components, if any.

The energy of quasiparticles participating in the superfluid motion is renormalized according to

[TABLE]

Hereafter we use the fact that inequalities $q_{\alpha}\ll p_{F\alpha}$, or equivalently, $v_{\alpha s}\ll v_{F\alpha}$ are well fulfilled, in superfluid Fermi liquids. Therefore we restrict ourselves to a linear in $q_{\alpha}/p_{F\alpha}$ terms to write $f^{\alpha\beta}\left(\mathbf{k\mathbf{+q}_{\alpha},k}^{\prime}\mathbf{+q}_{\beta}\right)\simeq f^{\alpha\beta}\left(\mathbf{k,k}^{\prime}\right)$.

Combining Eqs. (8) and (10) one finds that the renormalization of the quasiparticle energy owing to the motion of superfluids is of the form

[TABLE]

where the change of the distribution functions because of the supercurrents is given by

[TABLE]

All existing computations of the entrainment matrix (see, e.g. Borumand et al., 1996; Chamel & Haensel, 2006; Gusakov et al., 2009a) do not take into account the change in the energy gap $\Delta^{\left(\alpha\right)}\left(\mathbf{q}_{\alpha}\right)$ due to the motion of superfluid condensates and, consequently, are justified only for small amplitudes of oscillations that are restricted by the condition $v_{\alpha s}\ll v_{cr}^{\left(\alpha\right)}$, where $v_{cr}^{\left(\alpha\right)}\sim\Delta_{0}^{\left(\alpha\right)}/p_{F\alpha}$ is the critical velocity of the superfluid flow at which the superfluidity is destroyed ($\Delta_{0}^{\left(\alpha\right)}$ is the energy gap for Bogoliubov excitations in the superfluid at rest).

To develop the theory, which is valid in the case of near-critical supercurrents one needs to incorporate, into the distribution function, the elementary excitations arising in the system due to the motion of the superfluid condensate. To this end, it is necessary to express the quasiparticle operators $a_{\mathbf{k+q}_{\alpha}\uparrow}^{(\alpha)}$ and $a_{\mathbf{k+q}_{\alpha}\downarrow}^{(\alpha)}$, in the distribution function definition

[TABLE]

in terms of Bogoliubov excitation operators $b_{\pm\mathbf{k+q}_{\alpha}\sigma}^{(\alpha)}$ (see, e.g. Lifshitz & Pitaevskii, 1980). The corresponding conversion has been carried out by Gusakov & Haensel (2005). We shall only outline their result, referring the reader to their original work for more details. As it follows from Eqs. (26)-(30) of the work (Gusakov & Haensel, 2005) the distribution function (13) can be recast to the form

[TABLE]

where the functions

[TABLE]

represent the occupation numbers of Bogoliubov excitations of the energy $\mathfrak{E}_{\pm\mathbf{k+q}_{\alpha}}^{\left(\alpha\right)}$ in the system at the temperature $T$. The parameters $u_{\mathbf{k}}^{(\alpha)}$ and $v_{\mathbf{k}}^{(\alpha)}$ are defined as

[TABLE]

and

[TABLE]

Hereafter

[TABLE]

and $\tilde{\Delta}^{(\alpha)}$ denotes the energy gap, which depends on the superflow velocity.333Below it will be seen that the energy gap of each species $\alpha$ depends on the velocities of both superfluid flows, $\mathbf{v}_{\mathsf{p}s}$ and $\mathbf{v}_{\mathsf{n}s}$.

If the system is at rest, the Bogoliubov excitations are absent at zero temperature and the functions $\mathcal{F}_{\pm\mathbf{k}}^{(\alpha)}$ vanish, thus

[TABLE]

where

[TABLE]

is the energy of Bogoliubov excitations in the superfluid at rest.

In a superfluid flow, the functions $\mathcal{F}_{\pm\mathbf{k}}^{(\alpha)}$ does not vanish even at zero temperature. To demonstrate this let us expand the quasiparticle energy variation due to the superfluid motion up to the first order in the small parameter $q_{\alpha}/k_{F\alpha}\ll 1$.

[TABLE]

where $\gamma_{\alpha\beta}\left(\mathbf{v}_{\mathsf{p}s},\mathbf{v}_{\mathsf{n}s}\right)$ is some unknown matrix, which is taken at the Fermi surface of particle species (at $k=p_{F\beta}$) and depends on the flow velocities according to the relations $\mathbf{q}_{\alpha}\equiv m_{\alpha}\mathbf{v}_{\alpha s}$. Notice, the dependence of the matrix $\gamma_{\alpha\beta}$ on the flow velocities reflects the balance between the superfluid and normal components in the flow. This balance changes along with the flow velocity change and should be necessarily taken into account (see Sec. 4).

In the case of spin-singlet isotropic pairing of protons the form (21) is quite general because there are only two vectors $\mathbf{k}$ and $\mathbf{v}_{\mathsf{p}s}$ to form a scalar $\delta\varepsilon_{\mathbf{k+q}_{\alpha}}^{(\alpha)}$. The same form seems justified also in the case of anisotropic triplet pairing of neutrons if one assumes that the gap function is axial symmetric and that the preferred direction for the principal-axis of the gap matrix is specified by the direction of superflow motion, i.e. $Oz\parallel\mathbf{v}_{\mathsf{n}s}$. At least at zero temperature this choice of the principal axis seems reasonable because the direction of the flow is the only preferred direction which exists in the uniform system in the absence of external fields (see appendix).

Neglecting terms of the order of $(\Delta^{\left(\alpha\right)}/\mu_{\alpha})^{2}$ one can write $\varepsilon_{\mathbf{k}}^{(\alpha)}\simeq\xi_{k}^{(\alpha)}$. Than the quasiparticle energy (10) can be expanded as

[TABLE]

Equations for the unknown components of the matrix $\gamma_{\alpha\beta}$ can be obtained by combining of Eqs. (11) and (21). To this end one needs to calculate the change of the distribution functions owing to the supercurrents as is given in Eq. (12). We proceed to this calculation with the aid of Eqs. (14) - (18).

Inserting Eqs. (14), (16), (17), and (19) into Eq. (12) we get

[TABLE]

where $\mathbf{\hat{k}\,=k/}k$ denotes unit vector in the $\mathbf{k\,}$ direction.

In Eq. (23) the last term implements the fact that the relations $q_{\alpha}\ll p_{F\alpha}$, or equivalently, $v_{\alpha s}\ll v_{F\alpha}$ are well fulfilled, in superfluid Fermi liquids, and one can therefore keep only the terms linear in $\mathbf{q}_{\alpha}$, thus obtaining $\Theta_{\mathbf{k+q}}^{\left(\alpha\right)}-\Theta_{k}^{\left(\alpha\right)}\simeq-\mathbf{q}_{\alpha}\mathbf{\hat{k}\,}\delta\left(k-p_{F\alpha}\right)$.

General form for the energy of Bogoliubov excitations in the mixture of neutrons and protons with superfluid currents was found in Gusakov & Haensel (2005) by a minimization of the thermodynamic potential of the system. From Eqs. (29), (33), and (34) of that work one finds

[TABLE]

It is easily to see that this energy becomes negative provided

[TABLE]

The system becomes unstable when the flow velocity exceeds some critical value. From Eq. (25) it follows that if $\tilde{E}_{k}^{\left(\alpha\right)}-\sum_{\beta}\gamma_{\alpha\beta}\mathbf{kv}_{\beta s}\leq 0$ for $\mathbf{kv}_{\beta s}>0$ then $\tilde{E}_{k}^{\left(\alpha\right)}+\sum_{\beta}\gamma_{\alpha\beta}\mathbf{kv}_{\beta s}\leq 0$ for $\mathbf{kv}_{\beta s}<0$. When the energy of the Bogoliubov excitations vanishes, they are spontaneously created in pairs and accumulate in a superfluid system at zero temperature. As a result the distribution function for Bogoliubov excitations in the flow takes the form

[TABLE]

where $\Theta\left(x\right)$ is the Heaviside step-function.

As the velocities of superfluid flows are small in a scale of the Fermi velocities the arguments of the functions $f^{\alpha\beta}\left(\mathbf{k,k}^{\prime}\right)$ can be approximately put equal to their values at the corresponding Fermi surfaces. This allows one to write the interaction function in the form of expansion in Legendre polynomials parametrized by the Landau parameters $f_{l}^{\alpha\beta}$:

[TABLE]

For a nucleon matter the Landau parameters were calculated for various mean-field models in a series of papers (see Matsui, 1981; Henning & Manakos, 1987; Caillon et al., 2001, 2002, 2003), where it is shown that only first two spin-averaged Landau parameters are non-zero,444The Landau parameters with $l\geq 2$ are actually non-vanishing in nucleon matter, although they do not appear in simplified mean-field models. i.e., $f_{l}^{\alpha\beta}=0$ at $l\geq 2$. The same was found for nucleon-hyperon matter (Gusakov et al., 2009a). In view of this observation, one can employ only the parameters $f_{0}^{\alpha\beta}$ and $f_{1}^{\alpha\beta}$.

After summation over $\mathbf{k}^{\prime}$ in Eq. (11) the terms proportional to $f_{0}^{\alpha\beta}$ mutually cancel with the identical correction to the chemical potential.555In the system at rest, the quasiparticle energy is defined as $\xi_{k}^{(\alpha)}\equiv\epsilon^{(\alpha)}\left(k\right)-\mu_{\alpha}$, therefore, simultaneously with the energy correction $\delta\epsilon^{(\alpha)}$ one has to take into account the correction to the chemical potential $\delta\mu_{\alpha}$, caused by the same interactions. Since the $f_{0}^{\alpha\beta}$ contribution reduces to a constant mean field, it will give the identical additive to the energy and to the chemical potential, which mutually cancel in the expression for $\xi_{k}^{(\alpha)}$. Thus, only $f_{1}^{\alpha\beta}$ are to be taken into account. This contribution into the interaction (27) is odd with respect to the replacement $\mathbf{\hat{k}}^{\prime}\rightarrow-\mathbf{\hat{k}}^{\prime}$. Therefore, in Eq. (23) only odd terms, which are shown in the second line, will survive after substitution into Eq. (11) and subsequent integration over the solid angle $\mathbf{\hat{k}}^{\prime}$. Making use of this fact one can derive equations for the unknown matrix $\gamma_{\alpha\beta}\left(\mathbf{v}_{\mathsf{p}s},\mathbf{v}_{\mathsf{n}s}\right)$.

To this end we define the standard dimensionless Landau parameters

[TABLE]

where $N_{0\alpha}=m_{\alpha}^{\ast}p_{F\alpha}/\pi^{2}$, and introduce two auxiliary functions $\left(\alpha=\mathsf{p},\mathsf{n}\right)$:

[TABLE]

Substituting formulas (21) and (23) into Eq. (11) we get

[TABLE]

where

[TABLE]

In a standard way the summation over $\mathbf{k}^{\prime}$ can be converted into the integral

[TABLE]

The integral with respect to the angles on the right-hand side of Eq. (30) can be carried out using the addition theorem, which is valid for the Legendre polynomials:

[TABLE]

We get

[TABLE]

where the function $\Phi_{\alpha}\left(\tilde{\Delta}^{(\alpha)},V_{\alpha}\right)$ is given by the integral

[TABLE]

and the distribution function for Bogoliubov excitations is defined in Eq. (31).

To proceed we write $\mathbf{\hat{k}V}_{\alpha}$ in Eq. (34) in the form $\mathbf{\hat{k}V}_{\alpha}=\sum_{\beta}\gamma_{\alpha\beta}\mathbf{\hat{k}v}_{\beta s}$ and equate terms with the same $\mathbf{\hat{k}v}_{\alpha s}$ in the left- and right-hand sides thus obtaining the set of equations for $\gamma_{\alpha\beta}$:

[TABLE]

[TABLE]

Here $\alpha\neq\beta$ so that if $\alpha=\mathsf{n}$ then $\beta=\mathsf{p}$ and vice-versa. This formula must be supplemented by an expression for the effective mass (Sjöberg, 1973; Borumand et al., 1996):

[TABLE]

It is convenient to recast this relation in terms of the dimensionless Landau parameters:

[TABLE]

The Eqs. (36), (37) should be solved simultaneously with the gap equations in order to consistently take into account the gap dependence on the flow velocities. Below we consider these equations.

It is necessary to stress that Eqs. (36), (37) represent a set of nonlinear equations, where the integrals (35) depend on the unknown functions $\gamma_{\alpha\beta}$ and the flow velocities $\mathbf{v}_{\alpha s}$ according to Eqs. (31) and (29).

3 Gap equations

3.1 1S0 energy gap in a superfluid flow of protons

In the case of spin-singlet pairing the equation for dependence of the energy gap on the superflow velocity in a one-component superfluid is well known and repeatedly discussed in the literature (see, e.g. Bardeen, 1962; Alexandrov, 2003; Gusakov & Kantor, 2013). The 1S0 gap equation for a mixture of two interacting superfluids can be obtained in the same manner as the equation for a one-component superfluid. To take into account the interaction effects, it is only necessary to substitute the distribution functions for the Bogoliubov excitations with momentum $\pm\mathbf{k+q}_{\mathsf{p}}$, as indicated in the equation (31). Omitting the intermediate steps, we can write the gap equation in the form

[TABLE]

As in above, $\Delta_{0}^{\left(\mathsf{p}\right)}$ denotes the energy gap in the superfluid system at rest, that is at $\mathbf{q}=0$, while $\tilde{\Delta}^{(\mathsf{p})\,}$, as before, represents the energy gap in the presence of superfluid flows.

The integral over momenta in Eq. (40) may be carried out analytically. Then we obtain this equation in the form

[TABLE]

Here and below we use the notations

[TABLE]

The physical meaning of the introduced dimensionless parameters is quite understandable: $\eta_{\alpha}$ is the ratio of the gap amplitude of the species $\alpha=\mathsf{p,n}$ in the mixture of moving superfluids to the energy gap in the same superfluid at rest.

The parameter $\Lambda_{\alpha}$ is the effective flow velocity,

[TABLE]

in units of the characteristic velocity $\Delta_{0}^{\left(\alpha\right)}/p_{F\alpha}$. It is necessary to remind that the effective velocity $V_{\mathsf{p}}$ in Eq. (42) is defined by Eq. (43). Thus, in a mixture of two interacting superfluid flows, the energy gap of each species depends on the two velocities simultaneously.

Solution to Eq. (41) is depicted in Fig. 1, where the curve $\eta_{\mathsf{p}}$ is shown against the effective velocity $V_{\mathsf{p}}$ in units of $\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$.

One can see that at zero temperature there are two critical velocities. The first one,

[TABLE]

represents the flow velocity at which the Bogoliubov excitations appear in the superfluid system, and the energy gap begins to decrease.666At finite temperatures, $0<T<T_{c}$, the first critical velocity is smoothed because of the temperature breaking of the Cooper pairs. The number of broken Cooper pairs increases rapidly when the flow velocity builds up. Finally, all pairs are destroyed when the effective velocity reaches the second critical value

[TABLE]

This is the well known critical velocity at which the energy gap $\tilde{\Delta}^{\left(\mathsf{p}\right)}$ vanishes when the Fermi-liquid interactions are ignored. (see, e.g. Alexandrov, 2003; Gusakov & Kantor, 2013).

3.2 3P2 energy gap in a superfluid flow

We restrict our consideration to a non rotating neutron star and consider the unitary 3P2 states of the gap matrix.777The unitarity condition implies that the superfluid state under consideration retains time reversal symmetry and does not have, for example, spin polarization. In this case the equation for the amplitude of the energy gap $\tilde{\Delta}^{\left(\mathsf{n}\right)}$ in the superfluid neutron flow at zero temperature can be written as (see appendix)

[TABLE]

where $\Delta_{0}^{\left(\mathsf{n}\right)}$ is the energy gap in the superfluid at rest, and the distribution functions for the Bogoliubov excitations with momentum $\pm\mathbf{k+q}_{\mathsf{n}}$, are as indicated in Eq. (31):

[TABLE]

with

[TABLE]

Here $\mathbf{\bar{b}}\left(\mathbf{\hat{k}}\right)$ is some real vector in spin space which is normalized by the condition (82). It describes the angular dependence of the anisotropic energy gap which is given by $D_{\mathbf{\hat{k}}}=\tilde{\Delta}^{\left(\mathsf{n}\right)}\,\left|\mathbf{\bar{b}}\left(\mathbf{\hat{k}}\right)\right|$. A specific form of this vector depends on the phase state of the Cooper condensate. We restrict our consideration to the 3P2 pairing with a maximal projection $\left|M\right|=2$ of the total angular momentum onto the quantization axis which is directed along the superfluid flow motion (Arguments in favour of such a model, see appendix). In this case one has

[TABLE]

where $\theta$ is the angle between the quasiparticle momentum and direction of the superflow motion. It should be noted that, by virtue of Eq. (82), the amplitudes $\Delta_{0}$ and $\tilde{\Delta}$ are chosen as to represent the energy gap averaged over the solid angle. Defined in this way, the average energy gap furnishes an overall measure of the pairing correction to the ground-state energy in the preferred state.888It is necessary to notice that a definition of the gap amplitude is ambiguous in the literature. For example, in the case of $\left|M\right|=2$, our gap amplitude is $\sqrt{2/3}$ times larger than the gap amplitude in Ref. Yakovlev & Kaminker & Levenfish (1995). Ratio $\Delta_{0}^{\left(\left|M\right|=2\right)}/T_{c}=1.\,\allowbreak 6573$ differ in the same proportion from those reported in Yakovlev & Kaminker & Gnedin & Haensel (2001).

The integral over momenta in Eq. (46) may be carried out analytically. Then we obtain this equation in the form

[TABLE]

Solution to Eq. (50) is depicted in Fig. 2, where the curve $\eta_{\mathsf{n}}$ is shown against the effective velocity $V_{\mathsf{n}}$ in units of $\Delta_{0}^{\left(\mathsf{n}\right)}/p_{F\mathsf{n}}$.

Important to notice that the effective velocity $\mathbf{V}_{\alpha}$ is the auxiliary function, which is defined in Eq. (29). It coincides with a real velocity of the superfluid flow only in the absence of Fermi-liquid interactions, when $\gamma_{\alpha\alpha}=1$ and $\gamma_{\alpha\beta}=0$ for $\alpha\neq\beta$. The calculation of the entrainment matrix dependence on the realistic velocities of superfluid flows with allowance for Fermi-liquid interactions represents a complicated nonlinear problem.

4 Functions $\Phi_{\alpha}$

For the proton component of the superfluid mixture the integral in Eq. (35), can be easily calculated to get

[TABLE]

The curve $\eta_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)$ on the Fig. 1 clearly demonstrates that the condition $\Lambda_{\mathsf{p}}>\eta_{\mathsf{p}}$ is equivalent to the condition $\Lambda_{\mathsf{p}}>1$. Therefore one may replace $\Theta\left(\Lambda_{\mathsf{p}}-\eta_{\mathsf{p}}\right)\rightarrow\Theta\left(\Lambda_{\mathsf{p}}-1\right)$ in Eq. (51).

The simultaneous solving of Eqs. (36), (37) with the gap equation ensures that the value of the energy gap is taken along the curve $\eta_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)$ depicted in Fig. 1. This condition will be automatically fulfilled if to replace $\eta_{\mathsf{p}}\rightarrow\eta_{\mathsf{p}}\left(\Lambda_{\alpha}\right)$ in Eq. (51). To this end, we introduce the new function $\tilde{\Phi}_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)\equiv\Phi_{\mathsf{p}}\left[\Lambda_{\mathsf{p}},\eta_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)\right]$. Taking into account that $\Phi_{\mathsf{p}}=~{}1$ for $\Lambda_{\mathsf{p}}>e/2$ one can write

[TABLE]

Combining Eqs. (41) and (51) one can get the gap equation (41) in terms of the function $\Psi_{\mathsf{p}}$

[TABLE]

The function $\tilde{\Phi}_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)$ obtained as a result of numerical solution of this equation is shown in Fig. 1.

For the neutron component of the superfluid mixture the integral in Eq. (35) gives

[TABLE]

The equation for $\tilde{\Phi}_{\mathsf{n}}\left(\Lambda_{\mathsf{n}}\right)\equiv\Phi_{\mathsf{n}}\left[\Lambda_{\mathsf{n}},\eta_{\mathsf{n}}\left(\Lambda_{\mathsf{n}}\right)\right]$ can be found combining Eqs. (50) and (54). It is convenient to convert the unknown function to the form

[TABLE]

which takes into account that $\Phi_{\mathsf{n}}=~{}1$ for $\Lambda_{\mathsf{n}}>\sqrt{\left(3/2\right)e}$ (see appendix). Then Eq. (50) can be written in terms of the function $\Psi_{\mathsf{n}}$

[TABLE]

The function $\tilde{\Phi}_{\mathsf{n}}\left(\Lambda_{\mathsf{n}}\right)$ obtained as a result of numerical solution of this equation is shown in Fig. 2.

5 Entrainment matrix

Let us now return to the equations for the entrainment matrix. Using of the functions $\tilde{\Phi}_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)$ and $\tilde{\Phi}_{\mathsf{n}}\left(\Lambda_{\mathsf{n}}\right)$ allows one to eliminate the gap equations from the set of equations for the entrainment matrix thus solving only Eqs. (36), (37) but with the replacement $\Phi_{\alpha}\left(\Lambda_{\alpha},\eta_{\alpha}\right)\rightarrow\tilde{\Phi}_{\alpha}\left(\Lambda_{\alpha}\right)$. After the replacement the equations take the form

[TABLE]

[TABLE]

Here $\tilde{\Phi}_{\alpha}\equiv\tilde{\Phi}_{\alpha}\left(\Lambda_{\alpha}\right)$, and $\alpha\neq\beta$ so that if $\alpha=\mathsf{n}$ then $\beta=\mathsf{p}$ and vice-versa.

The equations (57), (58) are to be solved with taking into account that the effective velocity $V_{\alpha}\left(\gamma_{\alpha\beta},\mathbf{v}_{\mathsf{\alpha}s},\mathbf{v}_{\mathsf{\beta}s}\right)$ is a complicated function of the two flow velocities, as is given in Eq. (43). Thus, the set of equations (57), (58) should be supplemented also by the equations

[TABLE]

which reflect the fact that the gamma matrix depends actually on the specific flow velocities $\mathbf{v}_{\mathsf{\alpha}s},\mathbf{v}_{\mathsf{\beta}s}$. Moreover, the gamma matrix depends on the relative direction of the two flows.

To find a solution to Eqs. (57) and (58) we convert them to the form

[TABLE]

At first glance, these equations are similar to the equations obtained by Gusakov and Haensel in Gusakov & Haensel (2005), but the meaning of the functions $\tilde{\Phi}_{\alpha}\left(\Lambda_{\alpha}\right)$ is different. Actually, Eqs. (60)-(62) recover the result obtained in Gusakov & Haensel (2005) only when the flow velocities are small and $\tilde{\Phi}_{\alpha}\rightarrow 0$, see below.

6 Superfluid mass currents

We now turn to a calculation of the mass current caused by the superfluid flows. Substituting Eqs. (22) and (14) into Eq. (3) and performing simple integrations we find

[TABLE]

where the functions $\tilde{\Phi}_{\alpha}\left(\Lambda_{\alpha}\right)$ are defined in Eqs. (52), and (59)

[TABLE]

and $\gamma_{\alpha\beta}\,\left(\mathbf{v}_{\mathsf{\alpha}s},\mathbf{v}_{\mathsf{\beta}s}\right)$ denotes the solution to equations (57)-(59) for fixed $p_{F\alpha}$ and $p_{F\beta}$. As follows from Eqs. (63) and (64) the entrainment matrix can be written in the form

[TABLE]

where $\rho_{\alpha}=m_{\alpha}n_{\alpha}$. It represents a complicated nonlinear function of the velocities of superfluid flows $\mathbf{v}_{\mathsf{p}s}$ and $\mathbf{v}_{\mathsf{n}s}$.

It is important to notice that the right-hand side of Eqs. (60)-(62), as well as the entrainment matrix $\rho_{\alpha\beta}$ depend on $\Lambda_{\alpha}$ only. This greatly simplifies the problem. In practice one needs to solve numerically the set of equations

[TABLE]

to find $\Lambda_{\mathsf{p}}\left(\mathbf{v}_{\mathsf{p}s},\mathbf{v}_{\mathsf{n}s}\right)$ and $\Lambda_{\mathsf{n}}\left(\mathbf{v}_{\mathsf{p}s},\mathbf{v}_{\mathsf{n}s}\right)$ as functions of realistic flow velocities $\mathbf{v}_{\mathsf{p}s}$ and $\mathbf{v}_{\mathsf{n}s}$. Than, the entrainment matrix can be calculated with the aid of Eqs. (60)-(62) and (65).

7 Limiting cases and applications

7.1 Limiting case of small velocities

If the flow velocities are small as compared to the critical velocity, given in Eq. (44) we obtain $\Lambda_{\alpha}\left(\mathbf{v}_{\mathsf{\alpha}s},\mathbf{v}_{\mathsf{\beta}s}\right)\ll 1$ and thus $\tilde{\Phi}_{\alpha}=0$. Substituting this into Eqs. (57), (58) we get

[TABLE]

in agreement with the results obtained in Borumand et al. (1996); Gusakov & Haensel (2005).

In the opposite limiting case, when $V_{\alpha}\left(\mathbf{v}_{\mathsf{\alpha}s},\mathbf{v}_{\mathsf{\beta}s}\right)$ tends to the critical velocity, we obtain $\tilde{\Phi}_{\alpha}\rightarrow 1$ which leads to ${\rho}_{\alpha\beta}=0$ because the superfluidity disappears.

7.2 Mixture of superfluid protons and non-superfluid neutrons

We now turn to solving of Eqs. (57)-(59) in some particular cases typical for the nucleon matter in superfluid cores of neutron stars. Such a calculation requires a knowledge of the Landau parameters for a nucleon matter in beta equilibrium. This is a separate complicated problem, which is beyond the scope of our consideration. Since our goal is to illustrate the importance of the nonlinear effects, we shall employ the parameters derived in Gusakov et al. (2009b). Namely we consider the npe matter with a density of nucleons $n_{b}=n_{\mathsf{n}}+n_{\mathsf{p}}=3n_{0}$, where $n_{0}=0.16$ fm*-3* is the normal nuclear density, suggesting that the asymmetry parameter $\delta=\left(n_{\mathsf{n}}-n_{\mathsf{p}}\right)/n_{b}=0.837$. For this case it was obtained $F_{1}^{\mathsf{nn}}=-1,F_{1}^{\mathsf{pp}}=-0.55,F_{1}^{\mathsf{np}}=F_{1}^{\mathsf{pn}}=-0.3$. The superfluid transition temperatures for neutrons and protons are not reliably known. We adopt typical values used in models of neutron stars $T_{c\mathsf{p}}=5\times 10^{9}\,\mathsf{K}$ and $T_{c\mathsf{n}}=6\times 10^{8}\,\mathsf{K}$. For the energy gaps in the proton and neutron superfluids at rest we take $\Delta_{0}^{\left(\mathsf{p}\right)}=1.764T_{c\mathsf{p}}$, $\Delta_{0}^{\left(\mathsf{n}\right)}=1.657T_{c\mathsf{n}}$, respectively (see footnote8).

First we consider the case when only protons in the mixture of nucleons are superfluid and condensed into a spin-singlet state. Since the neutron superfluidity sets on at the temperature which is well below the critical temperature for Cooper pairing of protons, i.e., $T_{c\mathsf{n}}\ll T_{c\mathsf{p}}$, the limit of zero temperature can serve as a good approximation for calculation of the superfluid mass density of protons at $T_{c\mathsf{n}}<T\ll T_{c\mathsf{p}}$.

For normal (non-superfluid) neutrons one should put $\tilde{\Phi}_{\mathsf{n}}=1$ in Eqs. (57) and (58), thus obtaining $\rho_{\mathsf{nn}}=0$ and$\ \rho_{\mathsf{np}}=\rho_{\mathsf{pn}}={0}\,$. The remaining component of the entrainment matrix is given by

[TABLE]

where $\tilde{\Phi}_{\mathsf{p}}\left(\Lambda_{\mathsf{p}}\right)$ is given in Eqs. (52) and (53), and $\gamma_{\mathsf{pp}}\left(\Lambda_{\mathsf{p}}\right)$ obeys the equation

[TABLE]

The function $\Lambda_{\mathsf{p}}\left(v_{\mathsf{p}s}\right)$ is to be found from the equation

[TABLE]

If we denote

[TABLE]

then the partial density of superfluid protons can be written as999Notice, Eqs. (69) and (72) are similar to the equations obtained by Gusakov and Haensel in Gusakov & Haensel (2005), but the meaning of the functions $\tilde{\Phi}_{\mathsf{p}}$ and $\bar{\Phi}_{\mathsf{p}}$ is different.

[TABLE]

The result of numerical calculation is represented in Fig. 3, where the partial density of superfluid protons $\rho_{\mathsf{pp}}/\rho_{\mathsf{p}}$ is shown versus the superfluid velocity of protons. As it follows from Fig. 3 the mass density of superfluid protons drops very rapidly when the flow speed overcomes the value $v_{\mathsf{p}s}\simeq 0.65\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$ and completely disappears when the speed reaches $v_{\mathsf{p}s}\simeq 0.73\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$. Notice that this value is smaller than the known critical superfluid velocity (45) at which the energy gap collapses in a superfluid Fermi gas at zero temperature.

7.3 The case when protons and neutrons are superfluid

It is believed that protons and hyperons in the inner core of neutron stars condense into the isotropic 1S0 state, while the neutron pairing occurs in the triplet 3P2 state with an anisotropic energy gap (Tamagaki, 1970; Takatsuka, 1972).

Even at rest, the 3P2 superfluid neutron liquid has a complicated structure, which depends on the direction of the total angular momentum of a Cooper pair with respect to a quantization axis $\mathbf{\hat{n}}$. In a uniform system without external fields and at absolute zero, the orientation of the quantization axis in the superfluid at rest is arbitrary. The states of different orientation are degenerate. For equilibrium at a non-zero temperature this leads to the formation of a loose domain structure (Brueckner et al., 1960), where each domain has a preferred axis and these domains are randomly oriented. This fact is normally used for the superfluids at rest in order to simplify the calculations (see,e.g. Baiko et al., 2001). In the work (Gusakov & Haensel, 2005) this trick is used in the calculation of the entrainment matrix of a superfluid neutron-proton mixture. It is assumed that on the average the entrainment matrix is isotropic and can be calculated in the same way as in the case of the isotropic pairing but with some effective energy gap. Such an approach can be helpful only in the limit of small velocities of the superfluid flows, because the velocity dependence of the energy gap in the neutron superfluid is very different in the cases of isotropic and anisotropic pairing. This is clearly seen from a comparison of Fig. 1 and Fig. 2. Therefore, we will consider a more realistic model in which the pairing of neutrons in the moving superfluid liquid occurs in an anisotropic 3P2 state with a maximum value of the total angular momentum $\left|M\right|=2$, and the quantization axis is directed along the motion of the flow (see appendix).

We again consider the npe matter in beta equilibrium with a density of nucleons $n_{b}=n_{\mathsf{n}}+n_{\mathsf{p}}=3n_{0}$ with the interaction parameters, as in the previous example. The result of this calculation is demonstrated in Fig. 4, where the entrainment matrix components are shown as functions of the two velocities of superfluid flows of protons and neutrons. The computation was carried out for the case when superflows of neutrons and protons move in opposite directions.

Notice that the velocities of the two superfluid flows are scaled differently in the plots. Namely, the superfluid velocity of the proton flow $v_{\mathsf{p}s}$ is indicated in units of $\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$, while the superfluid velocity of the neutron flow $v_{\mathsf{n}s}$ is given in units of $\Delta_{0}^{\left(\mathsf{n}\right)}/p_{F\mathsf{n}}$. In practice, these two scales can be very different. For the case under consideration one can easily find that $\Delta_{0}^{\left(\mathsf{n}\right)}/p_{F\mathsf{n}}\ll\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$.

As can be seen, the $\rho_{\mathsf{pp}}$ component is almost independent of the velocity of superfluid neutron flow but strongly depends on the superfluid velocity of the proton flow. The $\rho_{\mathsf{pp}}$ value drops rapidly after the superfluid velocity of the proton flow exceeds approximately $0.65\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$ and vanishes completely when the velocity reaches the value $0.73\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$. Notice that this value is substantially smaller than the known critical superfluid velocity (45) at which the energy gap vanishes in a superfluid Fermi gas at zero temperature. Such a simple form of the $\rho_{\mathsf{pp}}$ component is due to the fact that, according to Eq. (59), protons give the dominant contribution to the dimensionless effective velocity $\Lambda_{\mathsf{p}}$.

The off-diagonal mixed components of the entrainment matrix coincide to one another, $\rho_{\mathsf{pn}}=\rho_{\mathsf{np}}$, in accordance with the general theory (Andreev & Bashkin, 1975). Negative values of these components vanish rapidly when the superfluid velocity of the proton flow exceeds $0.65\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$ and/or the superfluid velocity of the neutron flow exceeds $0.8\Delta_{0}^{\left(\mathsf{n}\right)}/p_{F\mathsf{n}}$.

The effective velocity $V_{\mathsf{n}}\left(\mathbf{v}_{\mathsf{n}s},\mathbf{v}_{\mathsf{p}s}\right)$ is governed by the magnitudes of both the superfluid velocities. Therefore the component $\rho_{\mathsf{nn}}$ rapidly drops along with growing of both the velocities. This is true, however, as long as both nucleon constituents are in superfluid states. When the superfluid velocity of the proton flow reaches $v_{\mathsf{p}s}\simeq 0.7\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$ protons become normal (non-superfluid), and their contribution to the effective velocity $V_{\mathsf{n}}$ is eliminated. The dependence of the entrainment matrix on the superflow velocity becomes very similar to the case, when only one of the components is superfluid in the mixture of nucleons. Namely, $\rho_{\mathsf{pn}}=\rho_{\mathsf{np}}=0$, $\rho_{\mathsf{pp}}=0$, and the $\rho_{\mathsf{nn}}$ value becomes independent of the proton velocity $v_{\mathsf{p}s}$ but drops rapidly along with growing of $v_{\mathsf{n}s}$.

Let us assume now that the superflows of neutrons and protons move in the same direction. Results of the calculation for this case are demonstrated in Fig. 5, where the components of the entrainment matrix are depicted as functions of the dimensionless velocities of the superfluid flows of protons and neutrons. We see that the $\rho_{\mathsf{pp}}$ component of the entrainment matrix is only slightly distinct from that component in Fig. 4, while remaining components are modified substantially due to the change of the relative direction of superfluid flows.

The domain of velocities where the mixed components, $\rho_{\mathsf{pn}}=\rho_{\mathsf{np}}$, differ from zero is a little larger than in Fig. 4. The explanation for this is the negative values of the off-diagonal elements of the $\gamma$-matrix. In this case, for the same directions of velocities $\mathbf{v}_{\mathsf{n}s}$ and $\mathbf{v}_{\mathsf{p}s}$, the effective velocity $V_{\mathsf{n}}$ is less than in the case of their opposite directions.

The behavior of the $\rho_{\mathsf{nn}}$ component can be interpreted in the same way as in the previous case. The velocity dependence of the component reduces to a more simple form when the proton flow speed $v_{\mathsf{p}s}$ exceeds $0.73\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$, and the proton superfluidity disappears. In that domain the component $\rho_{\mathsf{nn}}$ becomes independent of the proton flow velocity. Furthermore, $\rho_{\mathsf{pn}}=\rho_{\mathsf{np}}=0$, $\rho_{\mathsf{pp}}=0$, and $\rho_{\mathsf{nn}}$ drops rapidly when the superfluid velocity of the neutron flow $v_{\mathsf{n}s}$ exceeds $1.6\Delta_{0}^{\left(\mathsf{n}\right)}/p_{F\mathsf{n}}$.

As it follows from Figs. 4 and 5, at zero temperature, the linear approximation is accurate in a wide range of the superfluid flow velocities only for the $\rho_{\mathsf{pp}}$ component of the entrainment matrix that does not depend on the neutron flow velocity and has a plateau up to the superfluid proton velocity of about $0.65\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$. Notice, this magnitude is model dependent of the values of the Landau parameters used. For the remaining components the domain near the origin of the velocity axes, where the linear approximation can be employed is very small. These components of the entrainment matrix, can be calculated in the linear approximation only for $\mathbf{v}_{\mathsf{p}s}\ll\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$ and $\mathbf{v}_{\mathsf{n}s}\ll\Delta_{0}^{\left(\mathsf{n}\right)}/p_{F\mathsf{n}}$.

8 Summary and conclusion

In the frame of the Fermi-liquid theory we have derived equations for the entrainment matrix in a superfluid mixture of nucleons at arbitrary velocities of superfluid flows. The obtained equations account for the velocity dependence of the energy gaps and, therefore, are highly nonlinear. The problem is solved for the case of zero temperature in the assumption of an isotropic 1S0 pairing of protons and anisotropic 3P2 pairing of neutrons. We have considered the model, where the pairing of neutrons in a moving superfluid liquid occurs in the state with a maximum value of the total angular momentum $\left|M\right|=2$, and the quantization axis is directed along the motion of the flow. Arguments in favour of this model are discussed in appendix.

It is found, that in a mixture of two interacting superfluid flows, the energy gap of each species $\alpha=\mathsf{p,n}$ depends on the two superfluid speeds simultaneously.101010This property was first predicted qualitatively in (Gusakov & Kantor, 2013). This dependence is realized by means of the effective velocities for each species $V_{\mathsf{\alpha}}\left(\mathbf{v}_{\mathsf{n}s},\mathbf{v}_{\mathsf{p}s}\right)$ which are complicated functions of the realistic flow velocities $\mathbf{v}_{\mathsf{n}s},\mathbf{v}_{\mathsf{p}s}$. It is shown that the motion of superfluid condensate leads to instability of the superfluid state manifested in the appearance of Bogoliubov excitations at zero temperature. For the case of 1S0 pairing we have shown that, at zero temperature, there exist two critical flow velocities. The first one, which is given in Eq. (44), represents the velocity of the flow at which the system becomes unstable. At this velocity the minimum energy of Bogoliubov excitations vanishes. As a result such excitations are spontaneously created in pairs and accumulate in the superfluid system. The number of broken Cooper pairs increases rapidly when the flow velocity builds up. All the pairs are destroyed when the effective velocity reaches the second critical value, which is given in Eq. (45). This is the well known critical velocity at which the energy gap collapses in the superfluid Fermi gas. In the considered case of triplet pairing of neutrons creation of the Bogoliubov excitations starts immediately with the motion beginning. In both cases the excitations result in appearance of the normal component and the energy gap decrease.

Although the equations for the entrainment matrix are highly nonlinear they admit a simple analytical solution in the limiting case of small velocities of superfluid flows. In this case we get the results obtained in Borumand et al. (1996); Gusakov & Haensel (2005).

For arbitrary velocities of the superfluid flows the analytical solution of the nonlinear equations is impossible. Therefore, to demonstrate the importance of the nonlinear approach, we have solved the equations numerically for some particular cases typical for the nucleon matter in superfluid cores of non-rotating neutron stars. Such a calculation requires a knowledge of the Landau parameters for a nucleon matter in beta equilibrium. This is a separate problem, which is beyond the scope of our consideration. Since our goal is to illustrate the importance of the nonlinear effects, we have adopted the parameters derived in Gusakov et al. (2009b). Namely, we have considered the npe matter with the density three times larger the normal nuclear density, suggesting that the asymmetry parameter $\delta=0.837$. Numerical calculation was carried out for two cases.

In the first case, it was assumed that only protons are superfluid while the neutron liquid remains in a normal (non-superfluid) state. In this situation the limit of zero temperature can serve as a good approximation for calculation of the superfluid mass density of protons in the case $T_{c\mathsf{n}}<T\ll T_{c\mathsf{p}}$. Our computation has shown that the mass density of superfluid protons drops very rapidly when the flow speed overcomes the value $v_{\mathsf{p}s}\simeq 0.65\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$ and completely disappears when the speed reaches $v_{\mathsf{p}s}\simeq 0.73\Delta_{0}^{\left(\mathsf{p}\right)}/p_{F\mathsf{p}}$.

Next we have considered the case when both species of nucleons are superfluid. In this case the limit of zero temperature can be a reasonable approximation for the entrainment matrix at temperatures much below the critical temperature at which the neutron superfluidity sets on $T\ll T_{c\mathsf{n}}$. We found that the components of the entrainment matrix possess a highly nonlinear dependence on speeds of the two superflows simultaneously. Vanishing of the superfluidity in the superfluid mixture takes place at the flow velocities of the order of (but not equal to) the critical values known for single-component superfluids.

Vanishing of the superfluidity in a moving superfluid mixture can be crucial for oscillating neutron stars. As it follows from Figs. 4 and 5 the domain near the origin of the velocity axes, where the linear approximation can be employed is very small. On the other hand, theoretical estimates of the neutron star pulsations (Kantor & Gusakov, 2014; Gusakov & Kantor, 2013; Gusakov et al., 2005) show that even at a sufficiently small amplitude of oscillations, the relative velocity of superfluid and normal components of the neutron matter can exceed several times the critical value at which the superfluidity dies out. This can affect the oscillations themselves, making them highly nonlinear. Detailed discussions of the role of an entrainment in neutron star dynamics can be found in Prix (2004); Carter et al. (2005); Chamel & Carter (2006).

The results obtained in the present paper can be generalized to the case of finite temperature. This can be done making use of the temperature dependent distribution functions (15) in Eqs. (35), (40), and (46). The corresponding specific calculations require more numerical work. Evidently, at finite temperatures the entrainment matrix should vanish at smaller velocities of the superfluid flows.

Appendix A 3P2 energy gap in a superfluid flow of neutrons

Let us consider the gap equation for a system of neutral fermions interacting through an attractive spin-orbit potential. Let us consider a uniform superfluid flow of neutrons of a velocity $\mathbf{v}_{s}=\mathbf{q}/m$ in the rest frame of the normal component. For simplicity we consider the gap equation only for the case when the order parameter $\tilde{\Delta}_{\alpha\beta}\left(\mathbf{p}\right)$ is a unitary matrix in spin space, $\left(\alpha,\beta=\uparrow,\downarrow\right)$, which satisfies the condition

[TABLE]

where $\tilde{D}^{2}(\mathbf{p})$ is real. Here and below, the tilde above a letter indicates values that depend on the velocity of the superfluid flow. The unitarity condition implies that the superfluid state under consideration retains time reversal symmetry and does not have, for example, spin polarization. For this important case the general equation for the order parameter in a superfluid flow valid for any type of the pairing was derived in Fujita & Tsuneto (1972) in the BCS-Gor’kov approximation (see Eq. (12) of that work). We consider the case of 3P2 pairing, caused by spin-orbit interactions. In other words we assume that the effective interaction is most attractive for $S=1,L=1,J=2$ and that there is no tensor coupling involving the states with different $S$ and $L$.

Since the spin-orbit interaction conserves the total momentum, its matrix element can be written in the form

[TABLE]

It is adopted (see e.g. Takatsuka (1972); Fujita & Tsuneto (1972); Khodel & Clark & Zverev (2001)) to write $V_{\alpha\gamma,\beta\delta}\left(\mathbf{p,p}^{\prime}\right)$ in the 3P2 channel in the form of an expansion

[TABLE]

where $\Phi_{L,\alpha\beta}^{JM}\left(\mathbf{\hat{p}}\right)$ with $\mathbf{\hat{p}=p}/p$ are the standard spin-angle matrices

[TABLE]

In our case it is more convenient to use vector notations. Let $\mathbf{b}_{M}\left(\mathbf{n}\right)$ be the vectors in the spin space that generate the standard spin-angle matrices in accordance with

[TABLE]

where $\bm{\sigma}=\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)$ are Pauli spin matrices. There are five 3P2 states with a projection of the total angular momentum to the quantization axis $M=0,\pm 1,\pm 2$. The corresponding vectors are of the following explicit form

[TABLE]

The angular dependence of $\mathbf{b}_{M}(\mathbf{\hat{p}})$ is represented by Cartesian components of the unit vector $\mathbf{\hat{p}=p}/p$ which involves the polar angles on the Fermi surface, $\hat{p}_{x}=\sin\theta\cos\varphi$, $\hat{p}_{y}=\sin\theta\sin\varphi$, $\hat{p}_{z}=\cos\theta$. The vectors are normalized by the condition

[TABLE]

In the vector notations one can write the matrix element of the 3P2 interaction as

[TABLE]

The general form of the 3P2 order parameter is known to be of the form (see, e.g. Leinson, 2010)

[TABLE]

where $\tilde{\Delta}$ is a scalar amplitude of the energy gap which depends on the temperature and velocity of the superfluid flow, and $\mathbf{\bar{b}}\left(\mathbf{\hat{p}}\right)$ is some real vector in spin space which we normalize by the condition

[TABLE]

When substituting expressions (80) and (81) into Eq. (12) of the work Fujita & Tsuneto (1972), after performing spin traces one gets the equation for the 3P2 energy gap in the moving flow:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Here $\mu\simeq p_{F}^{2}/\left(2m\right)$ is the chemical potential of degenerate particles. In writing Eq. (84) we have neglected the terms quadratic in $q\ll p_{F}$. Since the volume of the system $\Omega$ drops out of the final results we hereafter set $\Omega=1$ for brevity.

We identify the energy gap as

[TABLE]

where $\mathbf{\bar{b}}\left(\mathbf{\hat{p}}\right)$ is chosen to be real in accordance to the unitarity condition (73), and note that the gap $\tilde{D}\left(\mathbf{p}\right)$ depends on the direction of the quasiparticle momentum and, in general, has nodes. The amplitude of the energy gap $\tilde{\Delta}$ must be real up to an arbitrary overall phase factor $e^{-i\mathbf{qR}}$. We, therefore, may adopt that the gap amplitude $\tilde{\Delta}\left(p\right)$ is a real function

Making use of the orthonormality condition (79) and the identity

[TABLE]

one can recast Eq. (83) to its final form

[TABLE]

In Eq. (89), the pairing interaction ${\large\upsilon}\left(k,p\right)$ depends on the absolute values of the relative particle momenta and the integration goes over infinite momentum space. However, theory of superfluidity needs the specific form of the pairing interaction only for a calculation of the order parameter (the energy gap) in the superfluid at rest. The latter completely defines all the properties of the superfluid. Let us demonstrate this for a dependence of the energy gap on the velocity of the superfluid flow. To this end we write Eq. (89) in the form

[TABLE]

where the functions $\mathcal{F}_{\pm\mathbf{k+q}}$ represent represent the occupation numbers of Bogoliubov excitations in the system, which are given by the Fermi distribution functions with a chemical potential equal to zero

[TABLE]

For a superfluid at rest (i.e. for $\mathbf{q}=0$) and temperature $T=0$ the gap equation (90) takes the form

[TABLE]

where

[TABLE]

Further we multiply Eq. (90) by $\Delta_{0}\,\left(k\right)$ and Eq. (92) by $\tilde{\Delta}\,\left(k\right)$. Combining the obtained equations we get

[TABLE]

The sums on both sides of this equation rapidly converge at a distance $\sim\Delta_{0}$ near the Fermi surface. Indeed, the functions $\mathcal{F}_{\pm\mathbf{p+q}}$ on the right-hand side decrease exponentially for $\xi_{p}\gg\Delta_{0}$. The function $(\tilde{E}_{\mathbf{p}}^{-1}-E_{\mathbf{p}}^{-1})$ also decreases rapidly with increase of a distance from the Fermi surface, since for $\xi_{p}\gg\Delta_{0}$ this function can be expanded as

[TABLE]

In the narrow vicinity of the Fermi surface the smooth functions may be replaced with their values at the Fermi surface, $\Delta_{0}\,\left(p\right)\rightarrow\Delta_{0}\,\left(p_{F}\right)\equiv\Delta_{0}$,$~{}\tilde{\Delta}\,\left(p\right)\rightarrow\tilde{\Delta}\,\left(p_{F}\right)\equiv\tilde{\Delta}$, ${\large\upsilon}\left(k,p\right)\rightarrow{\large\upsilon}\left(p_{F},p_{F}\right)$, and moved out from under the summation signs. This yields

[TABLE]

The summation on the left-hand side can be done with the aid of the formula

[TABLE]

This gives

[TABLE]

where $\Delta_{0}$ is uniquely related to the temperature of the superfluid transition $T_{c}$, which is assumed to be known, and $\tilde{\Delta}$ is the energy gap, which for a fixed density depends on the temperature and velocity of the superfluid flow. The distribution functions for the Bogoliubov excitations are given by

[TABLE]

with

[TABLE]

The solution to Eq. (98) allows one to find the gap $\tilde{\Delta}$ as a function of the superfluid velocity $\mathbf{v}_{s}$ at temperatures below the critical value $T_{c}$ for the superfluidity onset.

The vector $\mathbf{\bar{b}}$ can be generally written as $\bar{b}_{i}=A_{ij}\hat{p}_{j}$. In the case of a unitary 3P2 condensate the $3\times 3$ matrix $A_{ij}$ must be a real symmetric traceless tensor. It may be specified by giving the orientation of its principal axes and its two independent diagonal elements in its principal-axis coordinate system.

Since the ground state should be invariant under time reversal the states with magnetic quantum numbers $\pm M$ must be populated with equal likelihood. On the other hand, the gap tensor $A_{ij}$ should be diagonal. This second requirement excludes the possibility of populating states with $M=\pm 1$.  Therefore within the preferred coordinate system, there exist two simple solutions of Eq. (89).

For the first solution, which represents a condensation of the pairs into the state with $M=0$, we have $\mathbf{\bar{b}}(\mathbf{\hat{p}})=\mathbf{b}_{0}(\mathbf{\hat{p}})$ and $\mathbf{\bar{b}}^{2}(\mathbf{\hat{p}})=1/2\left(1+3\cos^{2}\theta\right)$. The second solution of Eq. (89) corresponds to $\left|M\right|=2$. In this case $\mathbf{\bar{b}}(\mathbf{\hat{p}})=1/\sqrt{2}\left[\mathbf{b}_{2}(\mathbf{\hat{p}})+\mathbf{b}_{-2}(\mathbf{\hat{p}})\right]$, and $\mathbf{\bar{b}}^{2}(\mathbf{\hat{p}})=3/2\left(1-\cos^{2}\theta\right)$. As is well known the two solutions, with $M=0$ and $\left|M\right|=2$, are almost degenerate in the neutron superfluid at rest.

We now turn to study these states in the moving condensate. In the superfluid at rest all directions for the principal-axis of the gap matrix are equivalent because the corresponding states are degenerate. However, the superfluid motion removes the degeneration. The terms proportional to $\mathbf{pv}_{s}$, in the quasiparticle energy (85), make the direction of the quantization axis along the motion more favourable. Therefore one may consider that at zero temperature the preferred direction for the principal-axis of the gap matrix is specified by the direction of superflow motion, i.e. $Oz\parallel\mathbf{q}$. This choice seems reasonable because the direction of the flow is the only preferred direction which exists in the uniform system in the absence of external fields.

Now we need to choose between the states $M=0$ and $\left|M\right|=2$. To this end, we can compare the maximum values of the superfluid flow velocity $\mathbf{v}_{s\max}=\mathbf{q}_{\max}/m$ at which the energy gap vanishes in each of the two cases. To find it we consider Eq. (98) for $\tilde{\Delta}/\Delta_{0}\ll 1$. In the limit $T\rightarrow 0$, this equation becomes of the form

[TABLE]

where $\xi_{0}=v_{F}/\Delta_{0}$ is the coherence length of correlated fermions at zero temperature, and $\Theta\left(z\right)$ is a unit-step function, $\Theta\left(z\right)=1$ if $z\geq 0$ and $\Theta\left(z\right)=0$ otherwise. Since in the rest frame the states $M=0$ and $\left|M\right|=2$ are almost degenerate one can put $\xi_{0}^{\left(M_{J}=0\right)}\simeq\xi_{0}^{\left(\left|M_{J}\right|=2\right)}=\xi_{0}$.

Doing the last integral in the limit $\tilde{\Delta}/\Delta_{0}\rightarrow 0~{}$ we get the relation

[TABLE]

For a 3P2 condensation with $M_{J}=0$ one has $\mathbf{\bar{b}}^{2}=\left(1/2\right)\left(1+3\cos^{2}\theta\right)$. In this case, from Eq. (102) we find limiting value of $q_{\max}^{\left(M_{J}=0\right)}$ at which the energy gap vanishes at zero temperature

[TABLE]

For a condensation with $\left|M_{J}\right|=2$ one has $\mathbf{\bar{b}}^{2}=\left(3/2\right)\sin^{2}\theta$. In this case Eq. (102) gives

[TABLE]

A comparison of Eqs. (103) and (104) allows one to prefer the mode, where the moving 3P2 condensate exists in the state with $\left|M_{J}\right|=2$. It is obviously true at least at $q>q_{\max}^{\left(M_{J}=0\right)}$.