A minimal model for finite temperature superfluid dynamics

A realistic model of a mature neutron star must account for various extreme aspects of physics. In addition to the strong gravity, we need to consider the state of matter at supranuclear densities. We also need to ask how the composition of matter changes as the density increases. Many different options, involving varying levels of exotica, have been suggested. For example, at the present time we do not know if the deep neutron star core contains deconfined quarks, perhaps in the form of a colour superconductor [1]. Meanwhile, the composition of the outer core is better understood [2]. Just above the crust–core transition density (at around 60% of the nuclear saturation density), the matter is predominantly made up of neutrons with a small fraction of protons and electrons (to make the conglomerate charge neutral). As the density increases muons appear, eventually followed by (this is where uncertainties creep in) hyperons and/or quark deconfinement. As a result, a neutron star is a bit like a layer cake where the composition of the different layers affects the transport coefficients which in turn determine the effect on cooling and the viscous damping of oscillations. Key to modeling this physics in a realistic fashion is the fact that several regions in the star are likely to exhibit superfluidity/superconductivity.

There is growing observational support for the presence of both superfluid neutrons and superconducting neutrons in neutron stars. The standard model for the spin-up glitches seen in radio pulsars (and some magnetars) [3] involves the transfer of angular momentum from an unseen superfluid component (which rotates by forming vortices) and the star's crust [4, 5]. However, the strongest current constraint on the theory comes from the apparent real-time cooling of the young remnant from the Cassiopeia A supernova. This neutron star is undergoing an apparent fast cooling epoch that is best explained through a recent onset of core superfluidity (following a significantly earlier formation of core superconductivity) [6, 7].

If we focus our attention on neutron star oscillations, then the added degree of freedom associated with superfluidity leads to the presence of extra families of oscillation modes [8–11] (analogous to the second sound in laboratory superfluids). This aspect of the problem has already been considered at some level of detail. However, most previous models of superfluid neutron star dynamics ignore finite temperature effects (beyond the fairly trivial observation that the entropy is likely to be coupled to the electrons, which form a normal fluid anyway). This may seem natural given that a typical neutron star is cold (in the sense that the temperature is orders of magnitude below the neutron Fermi temperature). While this statement is certainly true, we also need to acknowledge that the various superfluid pairing gaps are density dependent [12]. The upshot of this is that there should always be regions in a mature neutron star that are near the superfluid transition, see figure 1. In such a region finite temperature effects will dominate. These transition regions may be (highly) localized, but it would nevertheless be wise to consider them in detail. They may, in fact, have an unexpectedly large influence on the dynamics of the system. Consider for example the analogy with the viscous boundary layer at the crust–core interface and the suggested effect on the gravitational-wave driven r-mode instability [13–16]. In this case it is the finite extent of the boundary layer that leads to it having a large effect on the dynamics. Thus, without considering the problem in more detail we cannot know if the superfluid transition region is dynamically important or not.

Zoom In Zoom Out Reset image size

In order to address this problem, we need a model that accounts for finite temperature effects in neutron star superfluids. Such models have been discussed, in particular by Gusakov and collaborators [17–19], but the problem needs further attention. Our aim is to develop a phenomenological model that (i) enables us to consider finite temperature superfluid dynamics, and at the same time is (ii) flexible enough to allow for future extensions involving (for example) the star's magnetic field and/or the elastic crust. For obvious reasons the discussion will focus on the matter components, especially the extra degrees of freedom that enter in a superfluid system. Still, one should keep in mind that the complete model involves the 'live' spacetime of Einstein's theory. We will not discuss the associated (gravitational-wave) aspects here, but it is important to remember that the matter dynamics that we consider is coupled to the spacetime geometry. In other words, the spacetime metric g_ab is not taken to be fixed in the following. In any realistic application, one also has to solve (at least a subset of) the Einstein field equations.

As we intend to develop a covariant model, we distinguish between co- and contravariant quantities using the spacetime metric to relate the two. Spacetime indices are given as italic letters from the beginning of the alphabet, a, b, c, .... In cases where the focus is on spatial components the corresponding indices are i, j, k, ... . We assume that the spacetime signature is +2 and the Einstein summation convention applies. The spacetime indices should be distinguished from the constituent indices that are used to identify the different fluid components in the problem. These are given as roman, generically x, y, z,... and specifically (in the present context) n, p, e and s. Repetition of these component labels in a mathematical expression does not imply summation.

In contrast to the deep core, the conditions in the outer core of a neutron star are relatively well understood. In this region a (triplet pairing) neutron superfluid is thought to coexist with a (singlet pairing) type II proton superconductor [20] and a sea of relativistic electrons (we will ignore muons in the following, even though they are quite easy to account for). In a hot system we also need to account for the entropy. In principle, all four components can flow relative to one another. A general model would therefore consider four distinct fluid components; neutrons (labeled n in the following), protons (p), electrons (e) and the entropy (s). In order to focus on the thermal excitations that arise at finite temperatures, we will ignore two of the relative degrees of freedom throughout most of the discussion. The first represents the charge current due to a relative flow between protons and electrons (as required to maintain the global magnetic field [21]), and the second would correspond to a distinct heat flux [22, 23]. In principle, both of these aspects can be included in the model (at the cost of making the involved relations less manageable [21]) but the added complexity detracts from the key points of the analysis. Hence, we first of all assume that the protons and electrons are locked electromagnetically. Unless otherwise stated, we also add the assumption that the entropy is carried along with the only 'normal' component in the system, the electrons. This leaves us with two fluid degrees of freedom. This should be the minimal model for the problem of interest, namely how the presence of superfluid phase transitions impact on the fluid dynamics.

We take as our starting point the convective variational principle for multi-fluid systems devised by Carter [24] (see [25] for a review). In this model, the fluid dynamics is obtained from a matter Lagrangian Λ, which is taken to be a function of the various scalars that can be constructed from the individual fluxes $n_{\mathrm{x}}^a$, where the constituent index x = n, p, e or s. Basically, this means that Λ is a function of

$$\begin{equation} n_{\mathrm{x}}^2 = - n^{\mathrm{x}}_a n_{\mathrm{x}}^a, \qquad \mbox{ and} \qquad n_{{\mathrm{x}}{\mathrm{y}}}^2 = - n^{\mathrm{x}}_a n_{\mathrm{y}}^a. \end{equation} \tag{ 1 }$$

As will become apparent later, the second of these definitions involves effects due to the relative flow between the different components. While this flow is generally expected to be small in magnitude, its contribution is nevertheless significant since it encodes the entrainment effect [26–28]. This effect relies on quadratic (and higher order) terms in the Lagrangian. Hence, it is important that we do not make immediate use of any low-relative velocity assumption that may apply to the real systems. Such approximations should not be incorporated until the relevant momenta have been defined.

Note that the formalism is general and can accommodate any reasonable form for the matter Lagrangian. We are not advocating a particular model. Rather, the aim is to build a flexible formulation that can be adapted to deal with whatever information nuclear physics provides.

Given the general form for Λ, the variational procedure determines the momenta that are canonically conjugate to the fluxes [25]. We have

$$\begin{equation} \mu _a^{\mathrm{x}}= \frac{\partial \Lambda }{\partial n^a_{\mathrm{x}}} = \mathcal {B}^{\mathrm{x}}n_a^{\mathrm{x}}+ \sum _{{\mathrm{y}}\ne {\mathrm{x}}} \mathcal {A}^{{\mathrm{y}}{\mathrm{x}}} n_a^{\mathrm{y}}, \end{equation} \tag{ 2 }$$

where

$$\begin{equation} \mathcal {B}^{\mathrm{x}}= -2\frac{\partial \Lambda }{\partial n_{\mathrm{x}}^2}, \end{equation} \tag{ 3 }$$

$$\begin{equation} \mathcal {A}^{{\mathrm{y}}{\mathrm{x}}} = -\frac{\partial \Lambda }{\partial n_{{\mathrm{y}}{\mathrm{x}}}^2}. \end{equation} \tag{ 4 }$$

The latter encodes the entrainment between the different components. We will account for two distinct entrainment mechanisms. First of all, we know that the entrainment between neutrons and protons arises due to the strong interaction. Hence, we will have $\mathcal {A}^{{\rm p}{\rm n}}\ne 0$. We also know that the entropy entrainment plays a central role in ensuring causality of the heat flux [22, 23, 29]. The entropy entrainment is crucial in the variational two-fluid model for Helium, as it translates to the 'normal fluid density' in the standard description of that problem [30] (see also [31] where an analogous model is constructed from the field-theory point-of-view). Hence, we assume that $\mathcal {A}^{{\rm s}{\mathrm{x}}}\ne 0$ (with x = n, p, after locking the protons to the electrons and the entropy), as well. We will not account for entrainment between electrons and baryons.

It is useful to distinguish the entropy flux $s^a = n_{\rm s}^a$ and the associated thermal momentum $\theta _a = \mu ^{\rm s}_a$ [22–24] from the other components. This makes the description slightly less 'symmetric' but it makes it easier to make contact with thermodynamics. We then have

$$\begin{equation} \theta _a = \frac{\partial \Lambda }{\partial s^a} = \mathcal {B}^{\rm s}s_a + \sum _{{\mathrm{y}}\ne {\rm s}} \mathcal {A}^{{\mathrm{y}}{\rm s}} n_a^{\mathrm{y}}. \end{equation} \tag{ 5 }$$

Finally, a variation with respect to the spacetime metric tells us that the stress–energy tensor for the system takes the form

$$\begin{equation} T^a{}_{b} = \Psi \delta ^a{}_{b} + \sum _{\mathrm{x}}n_{\mathrm{x}}^a\mu _b^{\mathrm{x}}, \end{equation} \tag{ 6 }$$

where the generalized pressure is defined as

$$\begin{equation} \Psi = \Lambda - \sum _{\mathrm{x}}n_{\mathrm{x}}^a\mu _a^{\mathrm{x}}. \end{equation} \tag{ 7 }$$

The stress–energy tensor acts as source in Einstein's equations in the usual way. In the case where the fluxes are individually conserved, i.e., when the dynamical timescale is much faster than that associated with the relevant reactions, we also have

$$\begin{equation} \nabla _a n_{\mathrm{x}}^a = 0, \qquad \nabla _a s^a = 0. \end{equation} \tag{ 8 }$$

In this (the fully variational) case, each component satisfies an equation of motion of form

$$\begin{equation} 2n^a_{\mathrm{x}}\nabla _{[a}\mu ^{\mathrm{x}}_{b]} = 0, \end{equation} \tag{ 9 }$$

(where the square brackets indicate anti-symmetrization) and

$$\begin{equation} 2s^a\nabla _{[a}\theta _{b]} = 0. \end{equation} \tag{ 10 }$$

These Euler equations combine to ensure the conservation of energy and momentum, ∇^aT_ab = 0.

As in the case of relativistic heat flux [22, 23], it is useful to focus on a particular fluid frame. Basically, we need to decide which observer is measuring the temperature etcetera. In the present context, it seems natural to work in the entropy frame (associated with the flow of the 'normal' fluid). Hence, we introduce an observer with four velocity u^a such that s^a = su^a. In order to quantify the relative flow associated with the superfluid components, we define the relative velocities $v_{\mathrm{x}}^a$ such that

$$\begin{equation} n_{\mathrm{x}}^a = \gamma _{\mathrm{x}}n_{\mathrm{x}}\big(u^a + v_{\mathrm{x}}^a\big), \qquad u^av^{\mathrm{x}}_a = 0, \qquad \gamma _{\mathrm{x}}= \big(1-v_{\mathrm{x}}^2\big)^{-1/2}, \qquad v_{\mathrm{x}}^2 = v_{\mathrm{x}}^av^{\mathrm{x}}_a. \end{equation} \tag{ 11 }$$

In many cases of interest $v_{\mathrm{x}}^2$ will be small, and it makes sense to work with an expansion with these 'drift velocities' as small parameters. Moreover, rather than working with the individual velocities, it is often convenient to introduce the relative velocities $w_{{\mathrm{x}}{\mathrm{y}}}^a = v_{\mathrm{x}}^a - v_{\mathrm{y}}^a$. Then the entrainment parameter depends on

$$\begin{equation} n_{{\mathrm{x}}{\mathrm{y}}}^2 \approx n_{\mathrm{x}}n_{\mathrm{y}}\bigg( 1 + {w_{{\mathrm{x}}{\mathrm{y}}}^2 \over 2} \bigg), \end{equation} \tag{ 12 }$$

and we see that the equation of state, Λ, must contain information at order $w_{{\mathrm{x}}{\mathrm{y}}}^2$ in order for the momenta, $\mu ^{\mathrm{x}}_a$, to be correct at linear order. Problems that require higher order information must be considered on a case-by-case basis. For example, in the Hartle–Thorne slow-rotation approximation for a two-fluid system [32] (accurate to order $w_{{\mathrm{x}}{\mathrm{y}}}^2$) one can show that contributions originating from order $w_{{\mathrm{x}}{\mathrm{y}}}^4$ terms cancel exactly [33].

Let us, in the first instance, assume that the electrons are locked to the entropy, which should generally be the case. Then the number densities measured in the entropy frame (quantities in this frame will be denoted by * in the following) are

$$\begin{equation} n_{\mathrm{x}}^* = - u_a n_{\mathrm{x}}^a = \gamma _{\mathrm{x}}n_{\mathrm{x}}. \end{equation} \tag{ 13 }$$

We then have the momenta (associated with the three different fluxes)

$$\begin{equation} \mu _a^{\rm n}= \mathcal {B}^{\rm n}n_a^{\rm n}+ \mathcal {A}^{{\rm p}{\rm n}} n_a^{\rm p}+ \mathcal {A}^{{\rm s}{\rm n}} s_a = ( \mathcal {B}^{\rm n}n_{\rm n}^* + \mathcal {A}^{{\rm p}{\rm n}} n_{\rm p}^* + s \mathcal {A}^{{\rm s}{\rm n}} ) u_a + \mathcal {B}^{\rm n}n_{\rm n}^* v^{\rm n}_a + \mathcal {A}^{{\rm p}{\rm n}} n_{\rm p}^* v^{\rm p}_a, \end{equation} \tag{ 14 }$$

$$\begin{equation} \mu _a^{\rm p}= \mathcal {B}^{\rm p}n_a^{\rm p}+ \mathcal {A}^{{\rm n}{\rm p}} n_a^{\rm n}+ \mathcal {A}^{{\rm s}{\rm p}} s_a = ( \mathcal {B}^{\rm p}n_{\rm p}^* + \mathcal {A}^{{\rm n}{\rm p}} n_{\rm n}^* + s \mathcal {A}^{{\rm s}{\rm p}} ) u_a + \mathcal {B}^{\rm p}n_{\rm p}^* v^{\rm p}_a + \mathcal {A}^{{\rm n}{\rm p}} n_{\rm n}^* v^{\rm n}_a, \end{equation} \tag{ 15 }$$

and

$$\begin{equation} \theta _a = \mathcal {B}^{\rm s}s_a + \mathcal {A}^{{\rm n}{\rm s}} n_a^{\rm n}+ \mathcal {A}^{{\rm p}{\rm s}} n_a^{\rm p}= ( \mathcal {B}^{\rm s}s + \mathcal {A}^{{\rm n}{\rm s}} n_{\rm n}^* + \mathcal {A}^{{\rm p}{\rm s}} n_{\rm p}^* ) u_a + \mathcal {A}^{{\rm n}{\rm s}} n_{\rm n}^* v^{\rm n}_a + \mathcal {A}^{{\rm p}{\rm s}} n_{\rm p}^* v^{\rm p}_a. \end{equation} \tag{ 16 }$$

Note that we have not yet locked the electrons (which carry entropy and therefore move along with the entropy component) to the protons (we still have $v_{\rm p}^a\ne 0$). Given these expressions we see that the 'chemical potentials' (in the entropy frame) are

$$\begin{equation} \mu _{\rm n}= -u^a \mu ^{\rm n}_a = \mathcal {B}^{\rm n}n_{\rm n}^* + \mathcal {A}^{{\rm p}{\rm n}} n_{\rm p}^* + s \mathcal {A}^{{\rm s}{\rm n}}, \end{equation} \tag{ 17 }$$

$$\begin{equation} \mu _{\rm p}= -u^a \mu ^{\rm p}_a = \mathcal {B}^{\rm p}n_{\rm p}^* + \mathcal {A}^{{\rm n}{\rm p}} n_{\rm n}^* + s \mathcal {A}^{{\rm s}{\rm p}}, \end{equation} \tag{ 18 }$$

and

$$\begin{equation} \mu _{\rm e}= -u^a \mu ^{\rm e}_a = \mathcal {B}^{\rm e}n_e, \end{equation} \tag{ 19 }$$

while the temperature θ follows from

$$\begin{equation} \theta = - u^a \theta _a = \mathcal {B}^{\rm s}s + \mathcal {A}^{{\rm n}{\rm s}} n_{\rm n}^* + \mathcal {A}^{{\rm p}{\rm s}} n_{\rm p}^*. \end{equation} \tag{ 20 }$$

We can use these expressions to rewrite the stress–energy tensor and the pressure. This is not necessary but it clarifies the physical meaning of the variational coefficients. We need

$$\begin{equation} n_{\rm n}^a \mu ^{\rm n}_a = - n_{\rm n}^* \mu _{\rm n}+ n_{\rm n}^* \big( \mathcal {B}^{\rm n}n_{\rm n}^* v_{\rm n}^2 + \mathcal {A}^{{\rm p}{\rm n}} n_{\rm p}^* v_{\rm n}^a v^{\rm p}_a\big), \end{equation} \tag{ 21 }$$

$$\begin{equation} n_{\rm p}^a \mu ^{\rm p}_a = - n_{\rm p}^* \mu _{\rm p}+ n_{\rm p}^* \big( \mathcal {B}^{\rm p}n_{\rm p}^* v_{\rm p}^2 + \mathcal {A}^{{\rm n}{\rm p}} n_{\rm n}^* v_{\rm p}^a v^{\rm n}_a\big), \end{equation} \tag{ 22 }$$

$$\begin{equation} n_{\rm e}^a \mu ^{\rm e}_a = - n_{\rm e}\mu _{\rm e}, \end{equation} \tag{ 23 }$$

and

$$\begin{equation} s^a \theta _ a = - s \theta . \end{equation} \tag{ 24 }$$

Given these results we have

$$\begin{equation} \Psi = \Lambda + n_{\rm n}^* \mu _{\rm n}+ n_{\rm p}^* \mu _{\rm p}+ n_e \mu _e + s \theta - \big[ \mathcal {B}^{\rm n}( n_{\rm n}^* v_{\rm n})^2 + 2 \mathcal {A}^{{\rm p}{\rm n}} n_{\rm n}^* n_{\rm p}^* v_{\rm n}^a v^{\rm p}_a + \mathcal {B}^{\rm p}( n_{\rm p}^* v_{\rm p}) ^2 \big]. \end{equation} \tag{ 25 }$$

This is the (extended) Gibbs relation for the system, accounting for the presence of three distinct flows.

Let us now explicitly lock the protons to the electrons, i.e. assume that $v_{\rm p}^a=0$. Then we are left with

$$\begin{equation} \Psi = \Lambda + n_{\rm n}^* \mu _{\rm n}+ n_{\rm p}\mu _{\rm p}+ n_{\rm e}\mu _{\rm e}+ s \theta - \mathcal {B}^{\rm n}(n_{\rm n}^* v_{\rm n})^2. \end{equation} \tag{ 26 }$$

Finally, in most situations of practical interest the deviation from co-motion will be small. When that is the case we can neglect terms of order $v_{\rm n}^2$, leaving us with

$$\begin{equation} \Psi = \Lambda + n_{\rm n}\mu _{\rm n}+ n_{\rm p}\mu _{\rm p}+ n_{\rm e}\mu _{\rm e}+ s \theta . \end{equation} \tag{ 27 }$$

In the following we will assume that these various approximations are valid for the neutron star problem. Should we be interested in a more general context, it is straightforward to reinstate the redshift factors, etcetera.

Given the above relations, we can work out the corresponding form for the stress–energy tensor. To do this, we need

$$\begin{eqnarray} &&\sum _{\mathrm{x}}n_{\mathrm{x}}^a \mu ^{\mathrm{x}}_b \approx \left( \Psi -\Lambda \right)u^a u_b + n_{\rm n}\mu _{\rm n}v_{\rm n}^au_b + n_{\rm n}( \mathcal {B}^{\rm n}n_{\rm n}+ \mathcal {A}^{{\rm n}{\rm p}} n_{\rm p}+ \mathcal {A}^{{\rm n}{\rm s}} s) u^a v^{\rm n}_b \nonumber \\ &&= ( \Psi -\Lambda )u^a u_b + n_{\rm n}\mu _{\rm n}\big(v_{\rm n}^au_b + v^{\rm n}_b u^a \big). \end{eqnarray} \tag{ 28 }$$

The final expression for the stress–energy tensor is

$$\begin{eqnarray} T^a_{\ b} &=& \Psi \delta ^a_{\ b} + ( \Psi -\Lambda )u^a u_b + n_{\rm n}\mu _{\rm n}\big(v_{\rm n}^au_b + v^{\rm n}_b u^a \big) \nonumber \\ &=& -\Lambda u^a u_b + \Psi \perp ^a_{\ b} + n_{\rm n}\mu _{\rm n}\big(v_{\rm n}^au_b + v^{\rm n}_b u^a \big), \end{eqnarray} \tag{ 29 }$$

where we have used the standard projection

$$\begin{equation} \perp ^{ab}= g^{ab} + u^a u^b. \end{equation} \tag{ 30 }$$

It is worth noting that the relative flux enters $T^a_{\ b}$ in the same way as the heat flux would [22, 23]. This is, of course, as expected.

Noting that the energy density measured by an observer moving along with the entropy–proton–electron conglomerate is

$$\begin{equation} \rho = u_a u_b T^{ab} = -\Lambda , \end{equation} \tag{ 31 }$$

we recognize the Gibbs relation (27) as the definition of the pressure Ψ = P, with θ = T the temperature. In other words, in the linear case our final expressions are

$$\begin{equation} T^a_{\ b} = \rho u^a u_b + P \perp ^a_{\ b} + n_{\rm n}\mu _{\rm n}\left(v_{\rm n}^au_b + v^{\rm n}_b u^a \right), \end{equation} \tag{ 32 }$$

and

$$\begin{equation} P+\rho = n_n \mu _{\rm n}+ n_{\rm p}\mu _{\rm p}+ n_{\rm e}\mu _{\rm e}+ sT. \end{equation} \tag{ 33 }$$

The simplified problem has two dynamical degrees of freedom (the remaining fluxes). To determine these we need two momentum equations, e.g. (9) and (10). However, if we want to compare to the standard single-fluid problem, then it makes sense to first of all consider

$$\begin{equation} \nabla _b T^{ab} = 0, \end{equation} \tag{ 34 }$$

representing the conservation of total energy and momentum. In addition, we will use the Euler equation for the neutrons

$$\begin{equation} 2n^a_{\rm n}\nabla _{[a}\mu ^{\rm n}_{b]} = 0, \end{equation} \tag{ 35 }$$

where, cf, equation (14),

$$\begin{equation} \mu ^{\rm n}_a =\mu _{\rm n}u_a + \mathcal {B}^{\rm n}n_{\rm n}v^{\rm n}_a. \end{equation} \tag{ 36 }$$

So far, the neutron superfluidity is only manifest from the fact that the corresponding component is allowed to flow relative to the other fluids in the system, i.e. there is no viscous coupling. This kind of model arises in situations where quantized vortices are present in the system. Basically, the macroscopic equations are obtained by averaging over the vortices. In order to make contact with the standard analysis of laboratory systems, it is useful to consider the fact that a 'pure' superfluid should be irrotational. In this case the flow is potential, such that the momentum follows from

$$\begin{equation} \mu ^{\rm n}_a = m_{\rm n}\nabla _a \phi _{\rm n}, \end{equation} \tag{ 37 }$$

where ϕ_n is a scalar potential (proportional to the order parameter of the condensate). In this case (35) is automatically satisfied. Later, we need to make contact with the 'orthodox' model for superfluidity, centered on the notion of a superfluid (spatial) 'velocity' $V_{\rm n}^i$. This quantity represents the scaled momentum of the superfluid component [25], and (in the entropy frame) we have

$$\begin{equation} V_{\rm n}^i \equiv {1\over m_{\rm n}} \nabla ^i \phi _{\rm n}. \end{equation} \tag{ 38 }$$

From (36) we then see that we must have

$$\begin{equation} \mu _{\rm n}= -u^a \mu ^{\rm n}_a = -m_{\rm n}u^a \nabla _a \phi _{\rm n}, \end{equation} \tag{ 39 }$$

and

$$\begin{equation} V_{\rm n}^a = {1 \over m_{\rm n}}\perp ^a_{\ b} \mu _{\rm n}^b = {\mathcal {B}^{\rm n}n_{\rm n}\over m_{\rm n}} v_{\rm n}^a. \end{equation} \tag{ 40 }$$

This shows that the 'velocity' $V_{\rm n}^a$ is parallel, but not identical, to the relative velocity $v_{\rm n}^a$. As we will see later, the difference is associated with the entrainment effect (arising when the neutron effective mass differs from the bare mass).

The model we have developed is complete, at least formally. However, in neutron star modeling there tends to be a void between theory formulations and applications. The main reason for this is that the true supranuclear equation of state remains elusive. To date, most discussions have focused either on bulk properties (as required to construct stationary self-gravitating models, leading to mass and radius [2]) or the transport coefficients associated with cooling [34] (reaction rates etcetera). Comparatively little attention has been given to the coefficients required to model superfluid neutron star dynamics, like the entrainment or additional dissipation coefficients [35, 36]. The notable exceptions are the pioneering work of Alpar et al [37] and Borumand et al [38] at the Newtonian level and Comer and Joynt [39] in the relativistic case. The most detailed recent work is that of Chamel and collaborators, see [40–42] for results for the neutron star crust and [43] for the core, and Gusakov et al for hyperon cores [44, 45]. It is also worth noting recent discussions of the additional bulk viscosities that arise in superfluid systems [46–48]. Another important point concerns the fact that many quantum equation of state calculations are carried out in a non-relativistic framework. This is natural since the nuclear physics scale is always small enough that one can work in a local inertial frame [49]. However, it means that we need to be able to implement non-relativistic equation of state parameters in our fully relativistic framework for superfluid dynamics if we want to make progress (given the available information).

To discuss this issue, it is useful to first focus on a zero temperature model (essentially following [50]). Then we are dealing only with two components (the electrons are locked to the protons and there is no entropy in the system) and many of the issues become conceptually clearer. Moreover, as a first stab at the finite temperature problem it makes sense to build on the cold model, simply assuming that the thermal effects of (say) the strong interaction are small. We develop this strategy in this section and then return to the finite temperature problem.

Let us consider the low-relative velocity version of the variational model. In principle, the Newtonian limit involves two key steps. First, we need to allow for the fact that the relativistic Lagrangian Λ includes the rest-mass energies of the nucleons. Secondly, we need to consider the equations in the limit where the speed of light goes to infinity. The detailed analysis would proceed in the same way as in, for example, [25]. However, for the present purposes a short-cut will suffice (and hence we will not make the speed of light factors explicit in the following). We can extract the information we need directly from the momenta and the chemical potentials. In the case of the neutrons, and to leading order in the drift velocities, the momentum takes the form

$$\begin{equation} \mu ^{\rm n}_a \approx ( n_{\rm n}\mathcal {B}^{\rm n}+ n_{\rm p}\mathcal {A}^{{\rm n}{\rm p}} ) u_a + n_{\rm n}\mathcal {B}^{\rm n}v^{\rm n}_a + n_{\rm p}\mathcal {A}^{{\rm n}{\rm p}} v^{\rm p}_a. \end{equation} \tag{ 41 }$$

Here and in the following we let the two components be the neutrons and the 'protons' (n and p, respectively, with the latter including the electrons). The chemical potential (measured by an observer moving with u^a) follows from

$$\begin{equation} \mu _{\rm n}= - {\partial \Lambda \over \partial n_{\rm n}} = - u^a \mu ^{\rm n}_a \approx n_{\rm n}\mathcal {B}^{\rm n}+ n_{\rm p}\mathcal {A}^{{\rm n}{\rm p}}. \end{equation} \tag{ 42 }$$

In order to compare to the corresponding non-relativistic quantity, we use the decomposition

$$\begin{equation} \Lambda = - \sum _{\mathrm{x}}m_{\mathrm{x}}n_{\mathrm{x}}- E, \end{equation} \tag{ 43 }$$

where E is the total energy in the Newtonian setting. Defining the Newtonian chemical potential in the standard way, we have

$$\begin{equation} \mu _{\rm n}^{\mathrm{N}}= {\partial E \over \partial n_{\rm n}} = \mu _{\rm n}- m_{\rm n}. \end{equation} \tag{ 44 }$$

Meanwhile, the Newtonian momentum can be obtained from a Lagrangian of the form

$$\begin{equation} \mathcal {L} = {1\over 2} \sum _{{\mathrm{x}}{\mathrm{y}}} g_{ij} K^{{\mathrm{x}}{\mathrm{y}}}n_{\mathrm{x}}^i n_{\mathrm{y}}^j+{U}(n_{\mathrm{x}}) \qquad {\mathrm{x}}={\rm n},{\rm p}, \end{equation} \tag{ 45 }$$

expressed in terms of the mobility matrix, with elements K^xy, used by Chamel [50] and the internal energy U. The fluxes $n_{\mathrm{x}}^i=n_{\mathrm{x}}v_{\mathrm{x}}^i$ are (obviously) spatial (relative to u^a) and U depends only on the number densities. This leads to [35]

$$\begin{equation} p^{\rm n}_i = \left( {\partial \mathcal {L} \over \partial n_{\rm n}^i }\right)_{n_{\mathrm{x}}} = n_{\rm n}K^{{\rm n}{\rm n}} v^{\rm n}_i + n_{\rm p}K^{{\rm n}{\rm p}} v^{\rm p}_i, \end{equation} \tag{ 46 }$$

and

$$\begin{equation} \mu ^{\mathrm{N}}_{\rm n}= {\partial U \over \partial n_{\rm n}}. \end{equation} \tag{ 47 }$$

In the low-velocity limit, we have $\mu ^{\rm n}_i \approx p^{\rm n}_i$ so (from the spatial components of (41)) we identify

$$\begin{equation} \mathcal {B}^{\rm n}\approx K^{{\rm n}{\rm n}}, \end{equation} \tag{ 48 }$$

and

$$\begin{equation} \mathcal {A}^{{\rm n}{\rm p}} \approx K^{{\rm n}{\rm p}}. \end{equation} \tag{ 49 }$$

Finally, we note that Galilean invariance requires [50]

$$\begin{equation} n_{\rm n}K^{{\rm n}{\rm n}} + n_{\rm p}K^{{\rm n}{\rm p}} = m_{\rm n}. \end{equation} \tag{ 50 }$$

leading to the observation that (42) implies (to leading order)

$$\begin{equation} n_{\rm n}\mathcal {B}^{\rm n}= \mu _{\rm n}- n_{\rm p}\mathcal {A}^{{\rm n}{\rm p}} \approx m_{\rm n}-n_{\rm p}K^{{\rm n}{\rm p}}. \end{equation} \tag{ 51 }$$

The entrainment has an intuitive interpretation in terms of an effective (dynamical) mass. This quantity follows by considering the neutron momentum in the proton frame (where $v_{\rm p}^i=0$), and vice versa [28]. That is, we have

$$\begin{equation} m_{\rm n}^*= m_{\rm n}-n_{\rm p}K^{{\rm n}{\rm p}} \quad \longrightarrow \qquad n_{\rm p}K^{{\rm n}{\rm p}} = m_{\rm n}- m_{\rm n}^*, \end{equation} \tag{ 52 }$$

which means that

$$\begin{equation} n_{\rm n}K^{{\rm n}{\rm n}} = m_{\rm n}^*. \end{equation} \tag{ 53 }$$

A similar analysis for the protons leads to

$$\begin{equation} n_{\rm p}K^{{\rm p}{\rm p}} = m_{\rm p}^*, \end{equation} \tag{ 54 }$$

and

$$\begin{equation} n_{\rm n}K^{{\rm n}{\rm p}} = m_{\rm p}- m_{\rm p}^*. \end{equation} \tag{ 55 }$$

From this we deduce that that the two effective masses are related [28]:

$$\begin{equation} m_{\rm n}- m_{\rm n}^* = {n_{\rm p}\over n_{\rm n}} ( m_{\rm p}- m_{\rm p}^* ). \end{equation} \tag{ 56 }$$

Before we proceed, it is worth noting that we can (obviously) define an analogous effective mass, $\bar{m}_{\rm n}^*$, in the relativistic model [50]. In this case, the argument leads to (using (42))

$$\begin{equation} \bar{m}_{\rm n}^* = n_{\rm n}\mathcal {B}^{\rm n}= \mu _{\rm n}+ (m_{\rm n}^* - m_{\rm n}). \end{equation} \tag{ 57 }$$

This relation shows that, in the absence of interactions (read, entrainment), the effective relativistic mass is equal to the chemical potential. This is as expected.

Summarizing the discussion, we have shown how a relativistic model can be 'inferred' from the Newtonian limit in the sense that we can obtain $\mathcal {B}^{\mathrm{x}}$ and $\mathcal {A}^{{\mathrm{x}}{\mathrm{y}}}$ from the rest masses m_x, the effective masses $m_{\mathrm{x}}^*$ and the chemical potentials $\mu _{\mathrm{x}}^\mathrm{N}$ (which ultimately derives from the energy functional U). The procedure is not 'unique', but we have to resort to it in lieu of relativistic equation of state models. In particular, the analysis provides the translation between the variational model and the non-relativistic microphysics, e.g. the information encoded in the entrainment. The relevant information is expressed in terms of a single parameter, which can be taken to be the effective mass $m_{\rm p}^*$.

We now want to build on this strategy, extending it to account for finite temperature effects. To do this, we need to make contact with the 'orthodox' framework for superfluidity, because the problem has so far only been considered in that context [17, 51]. In other words, we want to compare our results to ones based on splitting each condensate into a 'normal' part and a 'superfluid' part (each expressed in terms of a mass density). This language provides a convenient description of phenomena observed in laboratory measurements. The variational model contains exactly the same information [30], although the system is decomposed in a different fashion. The preference of one model over the other is essentially a matter of personal choice. After all, one cannot truly decouple the involved degrees of freedom (e.g. the entropy and the nucleons, or the superfluid fraction and the condensate). They are inextricably linked.

As already mentioned, the orthodox model focuses on the superfluid velocities (the rescaled momenta). In the Newtonian case, we have $V^{\mathrm{x}}_i = p^{\mathrm{x}}_i/m_{\mathrm{x}}$ which leads to a total momentum flux

$$\begin{equation} j_{\mathrm{x}}^i = \rho _{\mathrm{x}}v_{\mathrm{x}}^i = \rho _{{\mathrm{x}}{\mathrm{x}}} V_{\mathrm{x}}^i + \rho _{{\mathrm{x}}{\mathrm{y}}} V_{\mathrm{y}}^i\ , \qquad {\mathrm{y}}\ne {\mathrm{x}}, \end{equation} \tag{ 58 }$$

where the entrainment is now expressed in the elements of the mass density matrix ρ_xy [26, 27]. Comparing to the variational model, we see that (after a bit of algebra and setting m_n = m_p = m for simplicity)

$$\begin{equation} K^{{\rm n}{\rm n}} = {m^2 \rho _{{\rm p}{\rm p}} \over \mathrm{det}\ \rho }, \qquad K^{{\rm p}{\rm p}} = {m^2 \rho _{{\rm n}{\rm n}} \over \mathrm{det}\ \rho }, \qquad \mathrm{and} \qquad K^{{\rm n}{\rm p}} = -{m^2 \rho _{{\rm n}{\rm p}} \over \mathrm{det}\ \rho }, \end{equation} \tag{ 59 }$$

where

$$\begin{equation} \mathrm{det}\ \rho = \rho _{{\rm p}{\rm p}} \rho _{{\rm n}{\rm n}}- \rho _{{\rm n}{\rm p}}^2. \end{equation} \tag{ 60 }$$

For future reference, it is worth noting that we can relate the mass density matrix to the effective masses as well. We then have

$$\begin{equation} \rho _{{\rm n}{\rm n}} = \rho _{\rm n}m_{\rm p}^* \bigg[ m_{\rm p}^* - {n_{\rm p}\over n_{\rm n}} (m-m_{\rm p}^*)\bigg]^{-1} \approx \rho _{\rm n}, \end{equation} \tag{ 61 }$$

$$\begin{equation} \rho _{{\rm p}{\rm p}} = \rho _{\rm p}\bigg[ m - {n_{\rm p}\over n_{\rm n}} (m-m_{\rm p}^*) \bigg] \bigg[ m_{\rm p}^* - {n_{\rm p}\over n_{\rm n}} (m-m_{\rm p}^*)\bigg]^{-1} \approx \rho _{\rm p}\bigg( {m \over m_{\rm p}^*} \bigg), \end{equation} \tag{ 62 }$$

and

$$\begin{equation} \rho _{{\rm n}{\rm p}} = - \rho _{\rm p}(m- m_{\rm p}^*) \bigg[ m_{\rm p}^* - {n_{\rm p}\over n_{\rm n}} (m-m_{\rm p}^*)\bigg]^{-1} \approx -\rho _{\rm p}\bigg( {m-m_{\rm p}^* \over m_{\rm p}^*}\bigg), \end{equation} \tag{ 63 }$$

where the approximations are valid whenever n_p/n_n ≪ 1. That is, they should be accurate for neutron star cores where the proton fraction tends to be below 10%. We will make use of these approximations in our discussion below.

So far, most of the literature on superfluid neutron star dynamics has ignored finite temperature effects. The exception is a series of paper by Gusakov and collaborators [17–19]. Their approach builds on the orthodox view of superfluidity and makes use of the relativistic hydrodynamics prescription of Son [52]. We want to make contact with this model in order to infer the key finite temperature effects. In particular, we need to be able to use the results of Gusakov and Haensel [51], who express the finite temperature effects in terms of the mass density matrix. In order to make use of these results we first of all need to translate them to the non-relativistic limit of our framework, and then elevate the model to relativity. This will not lead to a 'complete' model, but it is the best that we can do given the lack of microphysics information relating to the problem we are considering.

The starting point for the discussion is a three-component model, which accounts for the superfluidity of both neutrons and protons, and where the thermal excitations are expressed in the form of quasi-particles (index qp below). These are, in turn, locked to the electrons. Then the total non-relativistic fluxes $j_{\mathrm{x}}^i$ can be expressed as;

$$\begin{equation} j^i_{\rm n}= \rho _{\rm nn} V^i_{\rm n} + \rho _{\rm np} V^i_{\rm p} + (\rho _{\rm n} - \rho _{\rm nn}- \rho _{\rm np}) v^i_{\rm qp}, \end{equation} \tag{ 64 }$$

$$\begin{equation} j^i_{\rm p} = \rho _{\rm pp} V^i_{\rm p} + \rho _{\rm pn} V^i_{\rm n}+(\rho _{\rm p} - \rho _{\rm pp}- \rho _{\rm pn}) v^i_{\rm qp}, \end{equation} \tag{ 65 }$$

(here, and in the following, we adapt the notation in such a way that the comparison with the variational description in section 2 becomes transparent). These expressions account for the presence of entrainment through the ρ_xy coefficients in the manner discussed in the previous section. It has also been assumed that the 'normal velocities' of protons and neutrons are identical, represented by $v^i_\mathrm{qp}$ (this is a true transport velocity, in contrast to $V_{\mathrm{x}}^i$). The quasi-particles carry the entropy in the model.

The flux expressions are quite intuitive. The density of neutrons flowing with $V_{\mathrm{x}}^i$ is ρ_nx which leaves ρ_n − ρ_nn − ρ_np to flow with $v_\mathrm{qp}^i$ (and similarly for the protons). The quasi-particles, representing the normal fluid, are only present at finite temperatures. In the T = 0 limit we have [37]

$$\begin{equation} \rho _{\rm n}= \rho _{{\rm n}{\rm n}} + \rho _{{\rm n}{\rm p}}, \end{equation} \tag{ 66 }$$

$$\begin{equation} \rho _{\rm p}= \rho _{{\rm p}{\rm p}} + \rho _{{\rm n}{\rm p}}, \end{equation} \tag{ 67 }$$

bringing us back to the expressions discussed in the previous section. At a finite temperature these relations no longer hold, the difference being the presence of thermal excitations. This effect was quantified by Gusakov and Haensel [51], working in the normal-fluid frame. Basically, this means that we let $v_\mathrm{qp}^i=0$ which leads to an entrainment matrix with the same functional form as in the T = 0 case. The similarity is, however, deceptive as the coefficients are now temperature dependent. In particular, the three coefficients ρ_xy are independent, so we have a three-parameter entrainment model.

In order to compare the variational model to this description, it is useful to first consider a general frame not tied to any of the components. Identifying the quasi-particles with the entropy carrying component (the normal fluid) we have $v_\mathrm{qp}^i=v_{\rm s}^i$. Assuming that all relative velocities are small enough that we can linearize, the momenta will be given by

$$\begin{equation} p_{\rm n}^i =n_{\rm n}\tilde{K}^{{\rm n}{\rm n}} v_{\rm n}^i + n_{\rm p}\tilde{K}^{{\rm n}{\rm p}} v_{\rm p}^i + s \mathcal {A}^{{\rm s}{\rm n}} v_{\rm s}^i, \end{equation} \tag{ 68 }$$

and

$$\begin{equation} p_{\rm p}^i = n_{\rm n}\tilde{K}^{{\rm n}{\rm p}} v_{\rm n}^i + n_{\rm p}\tilde{K}^{{\rm p}{\rm p}} v_{\rm p}^i + s \mathcal {A}^{{\rm s}{\rm p}} v_{\rm s}^i. \end{equation} \tag{ 69 }$$

Here we have introduced the entropy entrainment [22, 29], encoded in $\mathcal {A}^{{\rm s}{\mathrm{x}}}$. We have assumed that this quantity is trivially related to the corresponding relativistic parameter (this makes sense since $\mathcal {A}^{{\rm n}{\rm p}} = K^{{\rm n}{\rm p}}$ in the cold model). In general, the coefficients of the mobility matrix will be affected by temperature, so it will not be the case that $\tilde{K}^{{\mathrm{x}}{\mathrm{y}}} = K^{{\mathrm{x}}{\mathrm{y}}}$. Of course, this depends on the physics involved. One may, for example, assume that the standard entrainment between neutrons and protons will be essentially the same as in the cold model because it originates from the strong interaction, which is only weakly temperature dependent [51]. The coupling to the quasi-particles (represented by the entropy carrying component above) complicates the situation. We need an argument that links the three parameters $\tilde{K}^{{\rm n}{\rm p}}$, $\mathcal {A}^{{\rm s}{\rm n}}$ and $\mathcal {A}^{{\rm s}{\rm p}}$ to the ρ_xy elements of the mass density matrix.

If we want to compare (68) and (69) to the previous relations, (64) and (65), we need to invert them to get the respective fluxes. Thus, we arrive at

$$\begin{equation} j_{\rm n}^i = m n_{\rm n}v_{\rm n}^i = {m \tilde{K}^{{\rm p}{\rm p}} \over \mathrm{det} \ \tilde{K}} p_{\rm n}^i - {m \tilde{K}^{{\rm n}{\rm p}} \over \mathrm{det} \ \tilde{K}} p_{\rm p}^i + {m\over \mathrm{det} \ \tilde{K}} ( \mathcal {A}^{{\rm s}{\rm p}} \tilde{K}^{{\rm n}{\rm p}} - \tilde{K}^{{\rm p}{\rm p}} \mathcal {A}^{{\rm s}{\rm n}}) s v_{\rm s}^i, \end{equation} \tag{ 70 }$$

where

$$\begin{equation} \mathrm{det} \ \tilde{K}= \tilde{K}^{{\rm n}{\rm n}} \tilde{K}^{{\rm p}{\rm p}} - ( \tilde{K}^{{\rm n}{\rm p}})^2. \end{equation} \tag{ 71 }$$

We also have

$$\begin{equation} j_{\rm p}^i = m n_{\rm p}v_{\rm p}^i = {m \tilde{K}^{{\rm n}{\rm n}} \over \mathrm{det} \ \tilde{K}}p_{\rm p}^i - {m\tilde{K}^{{\rm n}{\rm p}} \over \mathrm{det} \ \tilde{K}} p_{\rm n}^i + {m\over \mathrm{det} \ \tilde{K}} ( \mathcal {A}^{{\rm s}{\rm n}} \tilde{K}^{{\rm n}{\rm p}} - \tilde{K}^{{\rm n}{\rm n}} \mathcal {A}^{{\rm s}{\rm p}}) s v_{\rm s}^i. \end{equation} \tag{ 72 }$$

If we introduce the superfluid 'velocities' then these relations take the form;

$$\begin{equation} j_{\rm n}^i = {m^2 \tilde{K}^{{\rm p}{\rm p}} \over \mathrm{det} \ \tilde{K}}V_{\rm n}^i - {m^2 \tilde{K}^{{\rm n}{\rm p}} \over \mathrm{det} \ \tilde{K}}V_{\rm p}^i + {m\over \mathrm{det} \ \tilde{K}} ( \mathcal {A}^{{\rm s}{\rm p}} \tilde{K}^{{\rm n}{\rm p}} - \tilde{K}^{{\rm p}{\rm p}} \mathcal {A}^{{\rm s}{\rm n}}) s v_{\rm s}^i, \end{equation} \tag{ 73 }$$

and

$$\begin{equation} j_{\rm p}^i = {m^2 \tilde{K}^{{\rm n}{\rm n}} \over \mathrm{det} \ \tilde{K}} V_{\rm p}^i - {m^2 \tilde{K}^{{\rm n}{\rm p}} \over \mathrm{det} \ \tilde{K}} V_{\rm n}^i + {m\over \mathrm{det} \ \tilde{K}} ( \mathcal {A}^{{\rm s}{\rm n}} \tilde{K}^{{\rm n}{\rm p}} - \tilde{K}^{{\rm n}{\rm n}} \mathcal {A}^{{\rm s}{\rm p}}) s v_{\rm s}^i. \end{equation} \tag{ 74 }$$

At this point, the correspondence with the fluxes in the orthodox formulation is obvious. In particular, we see that the quasi-particle contribution vanishes as s → 0, as expected. It is also straightforward to establish that (59) still holds, although now for the $\tilde{K}^{{\mathrm{x}}{\mathrm{y}}}$ coefficients. By comparing to (64) and (65), and solving for the two entropy entrainment coefficients we find

$$\begin{equation} \mathcal {A}^{{\rm s}{\rm n}} = {m \over s \det \rho }( \det \, \rho - \rho _{{\rm p}{\rm p}}\rho _{\rm n}+ \rho _{{\rm n}{\rm p}}\rho _{\rm p}), \end{equation} \tag{ 75 }$$

and

$$\begin{equation} \mathcal {A}^{{\rm s}{\rm p}} = {m \over s \det \rho }( \det \, \rho - \rho _{{\rm n}{\rm n}}\rho _{\rm p}+ \rho _{{\rm n}{\rm p}}\rho _{\rm n}). \end{equation} \tag{ 76 }$$

This essentially completes the translation between the two models. In particular, we see how the entropy entrainment coefficients are related to the mass density matrix elements from Gusakov and Haensel [51]. That is, we have all information required for a variational finite temperature model for superfluid neutron stars.

In order to use these results to create a workable relativistic model, it first of all makes sense to assume that $\mathcal {A}^{{\rm s}{\mathrm{x}}}$ is not affected by relativistic effects (which seems reasonable given the discussion in the previous section). Next, we extend the analysis from section 3 to include the entropy component. This leads to the chemical potential being given by

$$\begin{equation} \mu _{\rm n}= n_{\rm n}\mathcal {B}^{\rm n}+ n_{\rm p}\mathcal {A}^{{\rm n}{\rm p}} + s \mathcal {A}^{{\rm s}{\rm n}} = m_{\rm n}c^2 + \mu _{\rm n}^\mathrm{N} , \end{equation} \tag{ 77 }$$

and (to leading order in the c → ∞ limit)

$$\begin{equation} n_{\rm n}\mathcal {B}^{\rm n}= n_{\rm n}\tilde{K}^{{\rm n}{\rm n}}. \end{equation} \tag{ 78 }$$

As in the cold case this identification (together with the analogous result for the protons) allows us to infer the parameters required in the relativistic case. The model is not unique, but it has the correct Newtonian limit by construction.

We have shown how we can infer the entropy entrainment needed to complete the variational model from the temperature dependent mass density matrix. This means that we can draw on the results of Gusakov and Haensel [51]. They show that the (non-relativistic) matrix elements take the form;

$$\begin{equation} \rho _{{\rm n}{\rm n}}(n,T) = [ 1- f_{\rm n}(T) ] \bar{\rho }_{{\rm n}{\rm n}}, \end{equation} \tag{ 79 }$$

$$\begin{equation} \rho _{{\rm p}{\rm p}}(n,T) = [ 1- f_{\rm p}(T) ] \bar{\rho }_{{\rm p}{\rm p}}, \end{equation} \tag{ 80 }$$

and

$$\begin{equation} \rho _{{\rm n}{\rm p}}(n,T) =[ 1- f_{\rm n}(T) ][ 1- f_{\rm p}(T) ]\bar{\rho }_{{\rm n}{\rm p}}, \end{equation} \tag{ 81 }$$

where $\bar{\rho }_{{\mathrm{x}}{\mathrm{y}}}$ is only weakly temperature dependent⁴. The main temperature behavior is encoded in the functions f_x, which can be taken from [53]. These functions vanish in the zero-temperature limit and approach unity as one approaches the critical temperature for superfluidity, when T → T_cx. This obviously means that ρ_xn → 0 as T → T_cn, reflecting the fact that all neutrons are 'normal' in this limit (similarly for the protons).

Given this, we may consider the behavior of the coefficients in the variational model as one of the critical temperatures is approached. For this exercise, let us assume that T_cn < T_cp (this would be true near the peak of the neutron triplet pairing gap in the example illustrated in figure 1). Introducing the shorthand notation

$$\begin{equation} R = \bar{\rho }_{{\rm n}{\rm n}} \bar{\rho }_{{\rm p}{\rm p}} - (1-f_{\rm n})(1-f_{\rm p}) \bar{\rho }_{{\rm n}{\rm p}}^2, \end{equation} \tag{ 82 }$$

we have

$$\begin{equation} \tilde{K}^{{\rm n}{\rm n}} = {m^2 \bar{\rho }_{{\rm p}{\rm p}} \over (1-f_{\rm n}) R}, \end{equation} \tag{ 83 }$$

$$\begin{equation} \tilde{K}^{{\rm p}{\rm p}} = {m^2 \bar{\rho }_{{\rm n}{\rm n}} \over (1-f_{\rm p}) R}, \end{equation} \tag{ 84 }$$

and

$$\begin{equation} \tilde{K}^{{\rm n}{\rm p}} = - {m^2 \bar{\rho }_{{\rm n}{\rm p}} \over R}, \end{equation} \tag{ 85 }$$

showing that $\tilde{K}^{{\mathrm{x}}{\rm p}}$ is regular as T → T_cn but $\tilde{K}^{{\rm n}{\rm n}}$ (and hence $\mathcal {B}^{\rm n}$ in the relativistic model) diverges. This behavior is in sharp contrast with the model of Gusakov et al [17, 51], where all parameters remain finite.

However, the divergence of $\tilde{K}^{{\rm n}{\rm n}}$ has a natural explanation in terms of the entrainment. As the system approaches the critical temperature, the dynamics will become dominated by the thermal excitations (an increasing fraction of the neutrons are 'normal'). Expressed in terms of the entrainment, this has the effect that it is increasingly difficult to move the neutrons relative to the entropy (the normal fluid). Essentially, their effective thermal mass increases. To see this, let us assume that the entropy is locked to the proton/electron component (as in section 2). Then the two momenta take the form

$$\begin{equation} p_{\rm n}^i = n_{\rm n}\tilde{K}^{{\rm n}{\rm n}} v_{\rm n}^i + ( n_{\rm p}\tilde{K}^{{\rm n}{\rm p}} + s \mathcal {A}^{{\rm s}{\rm n}} ) v_{\rm p}^i, \end{equation} \tag{ 86 }$$

and

$$\begin{equation} p_{\rm p}^i = n_{\rm n}\tilde{K}^{{\rm n}{\rm p}} v_{\rm n}^i + ( n_{\rm p}\tilde{K}^{{\rm p}{\rm p}} + s \mathcal {A}^{{\rm s}{\rm p}} ) v_{\rm p}^i, \end{equation} \tag{ 87 }$$

and we see that we can define effective masses $\tilde{m}_{\mathrm{x}}^*$ that include finite temperature effects;

$$\begin{equation} \tilde{m}_{\rm n}^* = n_{\rm n}\tilde{K}^{{\rm n}{\rm n}}, \end{equation} \tag{ 88 }$$

and

$$\begin{equation} \tilde{m}_{\rm p}^* = n_{\rm p}\tilde{K}^{{\rm p}{\rm p}} + s \mathcal {A}^{{\rm s}{\rm p}}, \end{equation} \tag{ 89 }$$

The first of these clearly diverges as the critical temperature T_cn is approached. The second, however, remains finite. This becomes clear once we note that $\tilde{K}^{{\rm p}{\rm p}}$ is regular and

$$\begin{equation} \mathcal {A}^{{\rm s}{\rm n}} = {m\over s} \left[ 1 - {\bar{\rho }_{{\rm p}{\rm p}}\rho _{\rm n}- (1-f_{\rm n})\bar{\rho }_{{\rm n}{\rm p}}\rho _{\rm p}\over (1-f_{\rm n})R }\right], \end{equation} \tag{ 90 }$$

and

$$\begin{equation} \mathcal {A}^{{\rm s}{\rm p}} = {m\over s} \left[ 1 - {\bar{\rho }_{{\rm n}{\rm n}}\rho _{\rm p}- (1-f_{\rm p})\bar{\rho }_{{\rm n}{\rm p}}\rho _{\rm n}\over (1-f_{\rm p})R }\right], \end{equation} \tag{ 91 }$$

where the former diverges as T → T_cn but the latter remains regular.

As a result of this behavior, the neutrons lock to the entropy (the 'normal' component). That is, we must have $v_{\rm n}^i\rightarrow v_{\rm s}^i$ in order to ensure that the dynamics remain regular. This is apparent from the fact that the neutron momentum (86) takes the form

$$\begin{equation} p_{\rm n}^i \rightarrow { \rho _{\rm n}\over (1-f_{\rm n}) \bar{\rho }_{{\rm n}{\rm n}}} \big(v_{\rm n}^i-v_{\rm p}^i\big) + \bigg[ m- {m\rho _{\rm p}\bar{\rho }_{{\rm n}{\rm p}} \over \bar{\rho }_{{\rm n}{\rm n}}\bar{\rho }_{{\rm p}{\rm p}}} \bigg] v_{\rm p}^i, \end{equation} \tag{ 92 }$$

as T → T_cn. Demanding that the momentum remains regular simply requires the neutrons to lock to the entropy+protons in the limit. This reflects the anticipated physics at the critical temperature, and a reduction in the number of distinct dynamical degrees of freedom.

In order to explore the main features of the finite temperature superfluid model that we have developed, let us consider a simplified dynamical problem. The main aim of this exercise is to shed light on the transition as the critical temperature is approached.

We take as our starting point the two-fluid model that results from locking the protons (electromagnetically) to the electrons and the entropy, cf, section 2. Meanwhile, the neutrons are superfluid and may move relative to the other components. We consider linear deviations from a uniform equilibrium where all fluids move together (the simplest situation). We assume that the fluids are in dynamical, chemical and thermal equilibrium in this state. This means that all background quantities (including the temperature) are constant (in time), and we also have μ_n = μ_p + μ_e. The equilibrium is associated with a frame such that u^a = [1, 0, 0, 0]. Perturbing this system, the normal fluid will move with u^a + δu^a, where u^aδu_a = 0, and the superfluid neutrons move with relative velocity $\delta v_{\rm n}^a$, where $u^a \delta v_a^{\rm n}= 0$. For simplicity, we will work in the local inertial frame associated with u^a. That is, we assume that spacetime is flat, which means that all derivatives are partial.

Let us first consider the conservation law for entropy, which leads to

$$\begin{equation} \partial _t \delta s + s \nabla _i \delta u^i = 0, \end{equation} \tag{ 93 }$$

where we have defined the co-moving time-derivative ∂_t = u^a∇_a and we have used the fact that the perturbed velocities are spatial (i = 1 − 3). Next we perturb the stress–energy tensor and work out its divergence. This leads to

$$\begin{eqnarray} &&\nabla _a \delta T^{ab} = (P+\rho ) u^b \nabla _a \delta u^a + (P+\rho ) \partial _t \delta u^b + u^b \partial _t (\delta P + \delta \rho ) \nonumber \\ &&+\, \nabla ^b \delta P + n_{\rm n}\mu _{\rm n}u^b \nabla _a \delta v_{\rm n}^a + n_{\rm n}\mu _{\rm n}\partial _t \delta v_{\rm n}^b = 0. \end{eqnarray} \tag{ 94 }$$

Contracting with u_b we arrive at

$$\begin{equation} \partial _t \delta \rho + (P+\rho ) \nabla _i \delta u^i + n_{\rm n}\mu _{\rm n}\nabla _i \delta v_{\rm n}^i = 0. \end{equation} \tag{ 95 }$$

Meanwhile, we have from the projection orthogonal to u^b;

$$\begin{equation} (P+\rho )\partial _t \delta u_i + \nabla _i \delta P + n_{\rm n}\mu _{\rm n}\partial _t \delta v^{\rm n}_i = 0 . \end{equation} \tag{ 96 }$$

In these two equations we recognize the standard perfect fluid terms, but also see the influence of the superfluid through the terms involving $\delta v_{\rm n}^i$. Next consider the neutron conservation law, which leads to

$$\begin{equation} \partial _t \delta n_{\rm n}+ n_{\rm n}\left( \nabla _i \delta u^i + \nabla _i \delta v_{\rm n}^i \right) = 0 , \end{equation} \tag{ 97 }$$

while the corresponding momentum equation can be written (the component along u^b vanishes identically)

$$\begin{equation} \mu _{\rm n}\partial _t \delta u_i + n_{\rm n}\mathcal {B}^{\rm n}\partial _t \delta v_i^{\rm n}+ \nabla _i \delta \mu _{\rm n}= 0. \end{equation} \tag{ 98 }$$

Based on these equations, it is easy to show that only the longitudinal (in the spatial sense) degrees of freedom are non-trivial. Hence, we take the divergence of (96) and (98) in order to obtain scalar equations. Then we can use (95) and (97) to remove the explicit dependence on the velocities, ending up with equations relating perturbed scalar variables. The relevant equations take the form;

$$\begin{equation} \partial ^2_t \delta \rho - \nabla ^2 \delta P = 0, \end{equation} \tag{ 99 }$$

and

$$\begin{equation} \left[ {\mu _n^2 - (P+\rho )\mathcal {B}^{\rm n}\over n_{\rm n}\mu _{\rm n}- (P+\rho ) }\right] \partial ^2_t \delta n_{\rm n}- \nabla ^2 \delta \mu _{\rm n}- \left[ {\mu _n - n_{\rm n}\mathcal {B}^{\rm n}\over n_{\rm n}\mu _{\rm n}- (P+\rho ) }\right] \partial ^2_t \delta \rho = 0. \end{equation} \tag{ 100 }$$

The first is, obviously, the usual wave equation for sound and the second equation is similar (even though the coefficients are a bit messier). If we consider the case where the neutrons are not interacting with the other particles (no entrainment), then we have $\mu _{\rm n}=n_{\rm n}\mathcal {B}^{\rm n}$ and the second equation collapses to

$$\begin{equation} \partial _t^2 \delta n_{\rm n}- {n_{\rm n}\over \mu _{\rm n}} \left( {{\rm d}\mu _{\rm n}\over {\rm d}n_{\rm n}}\right) \nabla ^2 \delta n_{\rm n}= 0, \end{equation} \tag{ 101 }$$

i.e., a decoupled wave equation with wave-speed [25]

$$\begin{equation} c_{\rm n}^2 = {n_{\rm n}\over \mu _{\rm n}} \left( {{\rm d}\mu _{\rm n}\over {\rm d}n_{\rm n}}\right). \end{equation} \tag{ 102 }$$

The general relations provide a neat demonstration of the existence of two identifiable wave degrees of freedom in the problem, leading to the existence of unique 'superfluid' oscillations modes in a neutron star core [8–11].

Before moving on, it is worth noting that the entropy is slaved to the other perturbations. Solving (95) and (97) for the divergence of δuⁱ and using the result in (93) we see that (after time integration)

$$\begin{equation} \delta s = {s \over P+\rho - n_{\rm n}\mu _{\rm n}} \left( \delta \rho - \mu _{\rm n}\delta n_{\rm n}\right). \end{equation} \tag{ 103 }$$

In order to progress further, we need to introduce a three-parameter equation of state. Opting to work with the density, neutron number density and entropy density as our primary variables, we have

$$\begin{equation} \delta P = \left( {\partial P\over \partial \rho } \right)_{n_{\rm n},s}\delta \rho + \left( {\partial P\over \partial n_{\rm n}} \right)_{\rho ,s}\delta n_{\rm n}+ \left( {\partial P\over \partial s } \right)_{\rho ,n_{\rm n}}\delta s, \end{equation} \tag{ 104 }$$

and

$$\begin{equation} \delta \mu _{\rm n}= \left( {\partial \mu _{\rm n}\over \partial \rho } \right)_{n_{\rm n},s}\delta \rho + \left( {\partial \mu _{\rm n}\over \partial n_{\rm n}} \right)_{\rho ,s}\delta n_{\rm n}+ \left( {\partial \mu _{\rm n}\over \partial s } \right)_{\rho ,n_{\rm n}}\delta s. \end{equation} \tag{ 105 }$$

In the following we will not, for the sake of clarity and brevity, state the variables that are held fixed explicitly. In these relations we can use (103) to replace δs, effectively reducing the problem to one with two parameters. Having done this, we can use the obtained expressions in (96) and (98), leading to two coupled wave equations for δρ and δn_n.

These equations are rather messy, with the first taking the form

$$\begin{equation} \partial _t^2 \delta \rho - a \nabla ^2 \delta \rho - b \nabla ^2 \delta n_{\rm n}= 0, \end{equation} \tag{ 106 }$$

where

$$\begin{equation} a = \left( {\partial P \over \partial \rho } \right) + {s \over P+\rho - n_{\rm n}\mu _{\rm n}} \left( {\partial P\over \partial s } \right), \end{equation} \tag{ 107 }$$

and

$$\begin{equation} b = \left( {\partial P\over \partial n_{\rm n}} \right) - {s \mu _{\rm n}\over P+\rho - n_{\rm n}\mu _{\rm n}} \left( {\partial P\over \partial s } \right). \end{equation} \tag{ 108 }$$

The second equation can be written

$$\begin{equation} c \partial ^2_t \delta n_{\rm n}- d \nabla ^2 \delta n_{\rm n}- e \nabla ^2 \delta \rho = 0, \end{equation} \tag{ 109 }$$

with

$$\begin{equation} c = {\mu _n^2 - (P+\rho )\mathcal {B}^{\rm n}\over n_{\rm n}\mu _{\rm n}- (P+\rho ) }, \end{equation} \tag{ 110 }$$

$$\begin{equation} d = \left( {\partial \mu _{\rm n}\over \partial n_{\rm n}} \right) - {s \mu _{\rm n}\over P+\rho - n_{\rm n}\mu _{\rm n}} \left( {\partial \mu _{\rm n}\over \partial s } \right) + \left[ {\mu _n - n_{\rm n}\mathcal {B}^{\rm n}\over n_{\rm n}\mu _{\rm n}- (P+\rho ) } \right] b, \end{equation} \tag{ 111 }$$

and

$$\begin{equation} e = \left[ {\mu _n - n_{\rm n}\mathcal {B}^{\rm n}\over n_{\rm n}\mu _{\rm n}- (P+\rho ) } \right] a + \left( {\partial \mu _{\rm n}\over \partial \rho } \right) + {s \mu _{\rm n}\over P+\rho - n_{\rm n}\mu _{\rm n}} \left( {\partial \mu _{\rm n}\over \partial s } \right). \end{equation} \tag{ 112 }$$

Given these expressions it is straightforward to work out the dispersion relation for plane waves in the system. Assuming that the perturbations behave as exp (iωt + ik_jx^j) and defining the phase-velocity σ = ω/k we find that

$$\begin{equation} ( \sigma ^2 - a) ( c \sigma ^2 - d) - be = 0. \end{equation} \tag{ 113 }$$

It is obviously easy to obtain the characteristic wave speeds in the system, once we introduce an explicit equation of state. We will, however, not do this. Instead, let us focus on two limiting cases. Consider first the low-temperature limit, i.e. when the system far below the critical temperature, deep into the superfluid region. Then we can take the s → 0 limit to obtain the dispersion relation;

$$\begin{equation} \left[ \sigma ^2 - \left({\partial P \over \partial \rho }\right) \right] \left[ \sigma ^2 - {n_{\rm n}\over \mu _{\rm n}} \left( {\partial \mu _{\rm n}\over \partial n_{\rm n}} \right)\right] - {n_{\rm n}\over \mu _{\rm n}} \left( {\partial P\over \partial n_{\rm n}} \right) \left( {\partial \mu _{\rm n}\over \partial \rho } \right) = 0. \end{equation} \tag{ 114 }$$

As a sanity check of this result, it is worth noting that in the zero-temperature limit we can rewrite the expression using n_n and n_p as the independent variables. Then defining

$$\begin{equation} c_{\mathrm{x}}^2 = {n_{\mathrm{x}}\over \mu _{\mathrm{x}}} \left( {\partial \mu _{\mathrm{x}}\over \partial n_{\mathrm{x}}} \right)_{{\rm n}_{\mathrm{y}}}, \end{equation} \tag{ 115 }$$

we arrive at the (more symmetric) dispersion relation

$$\begin{equation} \big( \sigma ^2 - c_{\rm n}^2\big) \big( \sigma ^2 - c_{\rm p}^2\big) - {n_{\rm n}\over \mu _{\rm n}} \bigg( {\partial \mu _{\rm n}\over \partial n_{\rm p}}\bigg) _{n_{\rm n}} {n_{\rm p}\over \mu _{\rm p}} \bigg( {\partial \mu _{\rm p}\over \partial n_{\rm n}}\bigg) _{n_{\rm p}} = 0. \end{equation} \tag{ 116 }$$

This result is identical to that derived in [25].

The other useful limit to consider corresponds to letting the temperature approach the critical temperature for the neutrons. We know from the discussion in the previous section that the parameters in the model are singular in this limit. In particular, we know that

$$\begin{equation} \mathcal {B}^{\rm n}\sim \tilde{K}^{{\rm n}{\rm n}} \rightarrow \infty , \end{equation} \tag{ 117 }$$

while the chemical potential μ_n remains finite. In order for the problem to remain regular, all terms proportional to $\mathcal {B}^{\rm n}$ in the equations must balance. Thus, we find that (109) becomes

$$\begin{equation} - {P+\rho \over n_{\rm n}} \partial ^2_t \delta n_n + b \nabla ^2 \delta n_{\rm n}+ a \nabla ^2 \delta \rho = 0. \end{equation} \tag{ 118 }$$

Adding this to (106) (and time integrating) we see that

$$\begin{equation} \delta n_{\rm n}= { n_{\rm n}\over P+ \rho } \delta \rho . \end{equation} \tag{ 119 }$$

That is, one of the dynamical degrees of freedom has been quenched. The remaining wave equation can be written

$$\begin{equation} \partial _t^2 \delta \rho - \left[ \left( {\partial P\over \partial \rho } \right) + {s \over P+\rho } \left( {\partial P\over \partial s } \right) + {n_{\rm n}\over P+\rho } \left( {\partial P\over \partial n_{\rm n}} \right) \right] \nabla ^2 \delta \rho = 0, \end{equation} \tag{ 120 }$$

which leads to the sound speed

$$\begin{equation} c_s^2 =\left( {\partial P\over \partial \rho } \right) + {s \over P+\rho } \left( {\partial P\over \partial s } \right) + {n_{\rm n}\over P+\rho } \left( {\partial P\over \partial n_{\rm n}} \right). \end{equation} \tag{ 121 }$$

To see that this is, indeed, the expected result we need to take a few more steps. First we note that the standard single fluid analysis is based on the total baryon number density n_b = n_n + n_p. Secondly, it is customary to work in terms of the entropy per baryon rather than the entropy number density. Thus we define $\bar{s} = s/n_\mathrm{b}$. We can then rewrite the thermodynamical derivatives in terms of ρ and $\bar{s}$. This leads to the final result

$$\begin{equation} c_s^2 = \left( {\partial P\over \partial \rho } \right)_{\bar{s}}, \end{equation} \tag{ 122 }$$

which is the standard expression for the (adiabatic) sound speed in the single-fluid case.

This model problem demonstrates that the finite temperature model we have developed connects the two-fluid region with the normal fluid region in a natural fashion.

Working in the context of the variational framework for relativistic multi-fluid systems [24, 25], we have developed a model that accounts for finite temperature effects in coupled superfluid systems. The model builds on recent progress on the analogous problem of heat flow [22, 23] and is based on treating the thermal excitations (e.g. phonons) in the systems as an additional entropy-carrying fluid. The approach breaks down in the zero-temperature limit, where the excitations propagate ballistically rather than collectively, but should be relevant for various systems of astrophysical interest. Our model is designed with the dynamics of the outer core of a mature neutron star (where superfluid neutrons are coupled to a conglomerate of protons and electrons) in mind, but the main ingredients are relevant for a range of analogous problems.

The developed formalism allows for detailed studies of the dynamics of neutron star superfluids near the critical density at which superfluidity sets in for a given (local) temperature. Setting the scene for more detailed work on this problem, we explored the key role of the entrainment between the entropy and the superfluid condensates, arguing that an increasing 'thermal effective mass' couples the components as the critical density/temperature is approached. We demonstrated this coupling in the specific case of the sound waves in the system. This exercise brings out the expected behavior; far below the critical temperature the systems exhibits second sound, but this phenomenon is quenched as one approaches the edge of the superfluid region. This result has implications for future work on the dynamics of neutron star superfluids. In particular it improves our understanding of the normal to superfluid transition, a key issue for studies of, for example, realistic neutron star seismology.

NA acknowledges support from STFC via grant number ST/J00135X/1. GLC acknowledges partial support from NSF via grant number PHYS-0855558.