Hydrodynamical instabilities in superfluid neutron stars with background flows between the components

TL;DR

This paper investigates hydrodynamical instabilities in superfluid neutron stars caused by background flows, analyzing oscillation modes and their potential link to observable phenomena like pulsar glitches.

Contribution

It introduces a detailed plane-wave analysis of oscillation modes considering background flows, entrainment, and different mutual friction models in superfluid neutron stars.

Findings

Standard mutual friction leads to unstable inertial and sound waves.

Entrainment causes quantitative differences between crust and core instabilities.

Isotropic mutual friction shows no unstable modes, suggesting turbulence may be linked to glitching.

Abstract

The interiors of mature neutron stars are expected to be superfluid. Superfluidity of matter on the microscopic scale can have a number of large scale, potentially observable consequences, as the superfluid component of the star can now flow relative to the `normal' component that is tracked by electromagnetic emission. The most spectacular of such phenomena are pulsar glitches, sudden spin-ups observed in many pulsars. A background flow of a normal-fluid component with respect to the superfluid is also known to lead to a number of instabilities in laboratory superfluids, possibly leading to turbulence and modifying the nature of the mutual friction coupling between the two fluids. In this paper we consider modes of oscillation in the crust and core of a superfluid neutron star, by conducting a plane-wave analysis. We explicitly account for a background flow between the two components (as would be expected in the presence of pinning) and the entrainment, and we consider both standard (Hall-Vinen) and isotropic (Gorter-Mellink) forms of the mutual friction. We find that for standard mutual friction there are families of unstable inertial and sound waves both in the case of a counter-flow along the superfluid vortex axis and for counterflow perpendicular to the vortex axis and find that entrainment leads to a quantitative difference between instabilities in the crust and core of the star. For isotropic mutual friction we find no unstable modes, and speculate that instabilities in a straight vortex array may be linked to glitching behaviour, which then ceases until the turbulence has decayed.

Figures

Tables

Table 1. Table 1: The two prototype cases, for the crust and core of the star, that have been tested in the numerical analysis of the determinant. For each of these five cases we consider nine different relative orientations of the local lag 𝐯 n subscript 𝐯 𝑛 \mathbf{v}_{n} and 𝐤 𝐤 \mathbf{k} and the three mutual friction scenarios listed in table 2 . Each of these cases has been investigated by considering two relative velocity speeds between neutrons and the normal component: | 𝐯 n | = 1 subscript 𝐯 𝑛 1 |\mathbf{v}_{n}|=1\, cm/s and | 𝐯 n | = 10 4 subscript 𝐯 𝑛 superscript 10 4 |\mathbf{v}_{n}|=10^{4}\, cm/s

Case	${\rho }_n$	$c_n$	$c_p$	${\epsilon }_n$	${\epsilon }_p$	$\Omega$	$\left|𝝎\right|$	$\widehat{𝝎}$	$A_{n,p}$	$C_{n,p}$
Crust	$4{\rho }_p$	${10}^9$cm/s	$0.51c_n$	-10	-40	$100$rad/s	$200$rad/s	$\widehat{𝐳}$	0	0
Core	$4{\rho }_p$	${10}^9$cm/s	$0.51c_n$	0.15	0.6	$100$rad/s	$200$rad/s	$\widehat{𝐳}$	0	0

Table 2. Table 2: The mutual friction scenarios considered in the numerical analysis: for each case listed in table 1 we consider these four mutual friction regimes.

Description	$B^{\prime }$	$B$	$ℛ$
Pinned vortices	1	0	$\infty$
Free vortices	0	0	0
Strong drag	0.5	0.5	1
Weak drag	${10}^{-4}$	${10}^{-2}$	$\approx {10}^{-2}$

Equations

$$\begin{array}{rl}& \left(\frac{\partial }{\partial t}+v_{\mathrm{x}}^j{\nabla }_j\right){\tilde{p}}_i^{\mathrm{x}}+{\epsilon }_{\mathrm{x}}{w}_j^{\mathrm{yx}}{\nabla }_i{v}_{\mathrm{x}}^j+\\ & +{\nabla }_i\left({\Phi }_R+{\tilde{\mu }}_{\mathrm{x}}\right)+2{ϵ}_{ijk}{\Omega }^j{v}_{\mathrm{x}}^k=f_i^{\mathrm{x}}.\end{array}$$

$${\Phi }_R=\Phi -\frac{1}{2}{\Omega }^2{r}^2\sin {\theta }^2,$$

$${\nabla }^2\Phi =4\pi G\sum _{\mathrm{x}}{\rho }_{\mathrm{x}}.$$

$${\tilde{p}}_i^{\mathrm{x}}=v_i^{\mathrm{x}}+{\epsilon }_{\mathrm{x}}{w}_i^{\mathrm{yx}},$$

$$\frac{\partial {\rho }_{\mathrm{x}}}{\partial t}+{\nabla }_j({\rho }_{\mathrm{x}}{v}_{\mathrm{x}}^j)=0.$$

$$f_i^{\mathrm{x}}=\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}{n}_{\mathrm{v}}{\mathcal{B}}^{{}^{\prime }}{ϵ}_{ijk}{\kappa }^j{w}_{\mathrm{xy}}^k+\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}{n}_{\mathrm{v}}\mathcal{B}{ϵ}_{ijk}{\hat{\kappa }}^j{ϵ}^{klm}{\kappa }_l{w}_m^{\mathrm{xy}},$$

$${\kappa }^i{n}_{\mathrm{v}}={\omega }^i={ϵ}^{ijk}{\nabla }_j{\tilde{p}}_k^{\mathrm{n}}+2{\Omega }^i,$$

$${\kappa }_i{n}_{\mathrm{v}}=2{\Omega }_i^{\mathrm{n}}+2{\epsilon }_{\mathrm{n}}({\Omega }_i^{\mathrm{p}}-{\Omega }_i^{\mathrm{n}}).$$

$$\mathcal{B}=\frac{\mathcal{R}}{1+{\mathcal{R}}^2},\text{mbox}and{\mathcal{B}}^{{}^{\prime }}=\frac{{\mathcal{R}}^2}{1+{\mathcal{R}}^2},$$

$$f_{\mathrm{x}}^i=-\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}{A}_{GM}{w}_{\mathrm{xy}}^2{w}_{\mathrm{xy}}^i,$$

$$Q(\mathbf{r},t)\approx Q_B(\mathbf{r})+\delta Q(\mathbf{r},t),$$

$$\frac{\partial \delta {\rho }_{\mathrm{x}}}{\partial t}+{\nabla }^j(v_j^{\mathrm{x}}\delta {\rho }_{\mathrm{x}}+{\rho }_{\mathrm{x}}\delta v_j^{\mathrm{x}})=0,$$

$$\begin{array}{rl}& ({\partial }_t+v_{\mathrm{x}}^j{\nabla }_j)(\delta v_i^{\mathrm{x}}(1-{\epsilon }_{\mathrm{x}})+w_i^{\mathrm{yx}}\delta {\epsilon }_{\mathrm{x}}+{\epsilon }_{\mathrm{x}}\delta v_i^{\mathrm{y}})+\delta v_{\mathrm{x}}^j{\nabla }_j{\tilde{p}}_i^{\mathrm{x}}+\\ & \delta ({\epsilon }_{\mathrm{x}}{w}_j^{\mathrm{yx}}{\nabla }_i{v}_{\mathrm{x}}^j)+{\nabla }_i(\delta {\tilde{\mu }}_{\mathrm{x}}+\delta {\Phi }_R)+2{ϵ}_{ijk}{\Omega }^j\delta v_{\mathrm{x}}^k=\delta f_i^{\mathrm{x}},\end{array}$$

$$\begin{array}{rl}& ({\partial }_t+v_{\mathrm{x}}^j{\nabla }_j)(\delta v_i^{\mathrm{x}}(1-{\epsilon }_{\mathrm{x}})+{\epsilon }_{\mathrm{x}}\delta v_i^{\mathrm{y}}+w_i^{\mathrm{yx}}\delta {\epsilon }_{\mathrm{x}})+\delta v_{\mathrm{x}}^j{\nabla }_j{v}_i^{\mathrm{x}}+\\ & {\epsilon }_{\mathrm{x}}\delta v_{\mathrm{x}}^j{\nabla }_j{w}_i^{\mathrm{yx}}+{\epsilon }_{\mathrm{x}}(\delta w_j^{\mathrm{yx}}{\nabla }_i{v}_{\mathrm{x}}^j+w_j^{\mathrm{yx}}{\nabla }_i\delta v_{\mathrm{x}}^j)+\delta {\epsilon }_{\mathrm{x}}{w}_j^{\mathrm{yx}}{\nabla }_i{v}_{\mathrm{x}}^j+\\ & {\nabla }_i\delta {\tilde{\mu }}_{\mathrm{x}}+2{ϵ}_{ijk}{\Omega }^j\delta v_x^k=\delta f_i^x,\end{array}$$

$$\delta w_{\mathrm{xy}}^i=\delta v_{\mathrm{x}}^i-\delta v_{\mathrm{y}}^i.$$

$$dE={\tilde{\mu }}_{\mathrm{x}}d{\rho }_{\mathrm{x}}+{\tilde{\mu }}_{\mathrm{y}}d{\rho }_{\mathrm{y}}+\alpha d{w}_{\mathrm{yx}}^2,$$

$$\displaystyle\tilde{\mu}_{\mathrm{x}}\,=\,\frac{\partial{E}}{\partial\rho_{\mathrm{x}}}\Big{|}_{\rho_{\mathrm{y}}w_{\mathrm{y}\mathrm{x}}^{2}}\,\quad\text{and}\quad\alpha\,=\,\frac{\partial{E}}{\partial w_{\mathrm{y}\mathrm{x}}^{2}}{\Big{|}}_{\rho_{\mathrm{x}}\rho_{\mathrm{y}}}\,.$$

$${\epsilon }_{\mathrm{x}}=\frac{2\alpha }{{\rho }_{\mathrm{x}}}.$$

$$\delta {\epsilon }_{\mathrm{x}}=\frac{2\delta \alpha }{{\rho }_x}-\frac{{\epsilon }_{\mathrm{x}}\delta {\rho }_x}{{\rho }_x}.$$

$$\delta {\tilde{\mu }}_{\mathrm{x}}=\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial {\rho }_{\mathrm{x}}}\delta {\rho }_{\mathrm{x}}+\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial {\rho }_{\mathrm{y}}}\delta {\rho }_{\mathrm{y}}+\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial w_{\mathrm{yx}}^2}\delta (w_{\mathrm{yx}}^2).$$

$$\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial w_{\mathrm{yx}}^2}=\frac{{\partial }^{2}E}{\partial w_{\mathrm{yx}}^2\partial {\rho }_{\mathrm{x}}}=\frac{\partial }{\partial {\rho }_{\mathrm{x}}}\left(\frac{\partial E}{\partial w_{\mathrm{yx}}^2}\right)=\frac{\partial \alpha }{\partial {\rho }_{\mathrm{x}}},$$

$$\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial w_{\mathrm{yx}}^2}=\frac{1}{2}{\epsilon }_{\mathrm{x}}+\frac{1}{2}{\rho }_{\mathrm{x}}\frac{\partial {\epsilon }_{\mathrm{x}}}{\partial {\rho }_{\mathrm{x}}}.$$

$$\begin{array}{rl}{c}_{\mathrm{x}}^2& ={\rho }_{\mathrm{x}}\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial {\rho }_{\mathrm{x}}}\\ C_{\mathrm{x}}& ={\rho }_{\mathrm{y}}\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial {\rho }_{\mathrm{y}}}={\rho }_{\mathrm{y}}\frac{\partial {\tilde{\mu }}_{\mathrm{y}}}{\partial {\rho }_{\mathrm{x}}}\\ A_{\mathrm{x}}& =\frac{4}{{\rho }_{\mathrm{x}}}\frac{\partial \alpha }{\partial w_{\mathrm{np}}^2}\\ {\alpha }_{\mathrm{x}}& =2\frac{\partial \alpha }{\partial {\rho }_{\mathrm{x}}}=2\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial w_{np}^2}={\epsilon }_{\mathrm{x}}+{\rho }_{\mathrm{x}}\frac{\partial {\epsilon }_{\mathrm{x}}}{\partial {\rho }_{\mathrm{x}}}.\end{array}$$

$${\epsilon }_{\mathrm{x}}+{\rho }_{\mathrm{x}}\frac{\partial {\epsilon }_{\mathrm{x}}}{\partial {\rho }_{\mathrm{x}}}\approx {\epsilon }_{\mathrm{x}}+n_B\frac{\partial {\epsilon }_{\mathrm{x}}}{\partial n_B}$$

$$\begin{array}{rl}\delta {\tilde{\mu }}_{\mathrm{x}}& =\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial {\rho }_{\mathrm{x}}}\delta {\rho }_{\mathrm{x}}+\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial {\rho }_{\mathrm{y}}}\delta {\rho }_{\mathrm{y}}+\frac{\partial {\tilde{\mu }}_{\mathrm{x}}}{\partial w_{\mathrm{yx}}^2}\delta w_{\mathrm{yx}}^2\\ & =c_{\mathrm{x}}^2\frac{\delta {\rho }_{\mathrm{x}}}{{\rho }_{\mathrm{x}}}+C_{\mathrm{x}}\frac{\delta {\rho }_{\mathrm{y}}}{{\rho }_{\mathrm{y}}}+{\alpha }_{\mathrm{x}}{w}_i^{\mathrm{yx}}\delta w_{\mathrm{yx}}^i.\end{array}$$

$$\begin{array}{rl}\delta \alpha & =\frac{\partial \alpha }{\partial {\rho }_{\mathrm{x}}}\delta {\rho }_{\mathrm{x}}+\frac{\partial \alpha }{\partial {\rho }_{\mathrm{y}}}\delta {\rho }_{\mathrm{y}}+\frac{\partial \alpha }{\partial w_{\mathrm{yx}}^2}\delta w_{\mathrm{yx}}^2\\ & =\frac{{\alpha }_{\mathrm{x}}{\rho }_{\mathrm{x}}}{2}\frac{\delta {\rho }_{\mathrm{x}}}{{\rho }_{\mathrm{x}}}+\frac{{\alpha }_{\mathrm{y}}{\rho }_{\mathrm{y}}}{2}\frac{\delta {\rho }_{\mathrm{y}}}{{\rho }_{\mathrm{y}}}+\frac{A_{\mathrm{x}}{\rho }_{\mathrm{x}}}{2}{w}_i^{\mathrm{yx}}\delta w_{\mathrm{yx}}^i.\end{array}$$

$$\delta {\epsilon }_{\mathrm{x}}=A_{\mathrm{x}}{w}_{\mathrm{yx}}^i\delta w_i^{\mathrm{yx}}+({\alpha }_{\mathrm{x}}-{\epsilon }_{\mathrm{x}})\frac{\delta {\rho }_{\mathrm{x}}}{{\rho }_{\mathrm{x}}}+{\alpha }_{\mathrm{y}}\frac{\delta {\rho }_{\mathrm{y}}}{{\rho }_{\mathrm{x}}}.$$

$$\begin{array}{rl}& \delta f_i^{\mathrm{x}}=\delta \left(\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}\right)[n_{\mathrm{v}}{\kappa }^j{\mathcal{B}}^{{}^{\prime }}{ϵ}_{ijk}{w}_{\mathrm{xy}}^k+n_{\mathrm{v}}{\kappa }_l\mathcal{B}{ϵ}_{ijk}{\hat{\kappa }}^j{ϵ}^{klm}{w}_m^{\mathrm{xy}}]+\\ & \frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}{n}_{\mathrm{v}}{\kappa }^j{\mathcal{B}}^{{}^{\prime }}{\epsilon }_{ijk}\delta w_{\mathrm{xy}}^k+\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}{n}_{\mathrm{v}}{\kappa }_l\mathcal{B}{\epsilon }_{ijk}{\hat{\kappa }}^j{\epsilon }^{klm}\delta w_m^{\mathrm{xy}}+\\ & \frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}\delta (n_{\mathrm{v}}{\kappa }^j){\mathcal{B}}^{{}^{\prime }}{\epsilon }_{ijk}{w}_{\mathrm{xy}}^k+\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}\delta (n_{\mathrm{v}}{\kappa }_l)\mathcal{B}{\epsilon }_{ijk}{\hat{k}}^j{ϵ}^{klm}{w}_m^{\mathrm{xy}}+\\ & \frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}{n}_{\mathrm{v}}{\kappa }_l\mathcal{B}{\epsilon }_{ijk}\delta ({\hat{\kappa }}^j){\epsilon }^{klm}{w}_m^{\mathrm{xy}}.\end{array}$$

$$\delta f_i^{\mathrm{x}}=(\delta f_i{)}_{sound}+(\delta f_i{)}_w+(\delta f_i{)}_{n_{\mathrm{v}}\kappa }+(\delta f_i{)}_{\hat{\kappa }}.$$

$$\begin{gathered} (\delta f_i{)}_{sound}=\left(\frac{\delta {\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}}-\frac{{\rho }_{\mathrm{n}}}{{\rho }_{\mathrm{x}}^2}\delta {\rho }_{\mathrm{x}}\right)(n_{\mathrm{v}}{\kappa }^j{\mathcal{B}}^{{}^{\prime }}{ϵ}_{ijk}{w}_{\mathrm{xy}}^k+ \\ +n_{\mathrm{v}}{\kappa }_l\mathcal{B}{ϵ}_{ijk}{\hat{\kappa }}^j{\epsilon }^{klm}{w}_m^{\mathrm{xy}}). \end{gathered}$$

Full text

††pagerange: Hydrodynamical instabilities in superfluid neutron stars with background flows between the components.–A.2††pubyear: 2019

1 Introduction

Neutron stars (NSs) are unique objects composed of matter at extremes of temperature, density and pressure that cannot be reproduced terrestrially. Describing their complex inner structure requires theoretical methods consistent with observations, which to date have been based on the analysis of the electromagnetic radiation emitted by the star, which carries mostly information on the processes in the magnetosphere or outermost, low density, layers. Despite the recent detection of gravitational waves emitted during a NS merger (Abbott et al., 2017), which is likely to deliver interesting constraints on the equation of state of matter at high densities in the near future (Abbott et al., 2018), direct data from the deep layers remains essentially inaccessible.

Nevertheless the state of matter in the stellar interior has a strong impact on the dynamics and observed astrophysical signals. In particular at high densities in the inner crust of a NS, neutrons can pair and form a superfluid, with protons also expected to be superconducting in the core (see Haskell & Sedrakian (2018) for a recent review). Superfluidity and superconductivity profoundly alter the dynamics of the components that can now flow relative to each other. The most striking phenomena connected to superfluidity are glitches, sudden spin-ups of the star, that are instantaneous to the accuracy of the data, and are generally assumed to be due to the sudden re-coupling of the interior superfluid neutrons and the normal crust, which is tracked by the electromagnetic emission (Anderson & Itoh, 1975; Haskell & Melatos, 2015). Glitches allow for an indirect probe of the NS interior, and can be used to obtain constraints on physical parameters of the star (Alpar et al., 1984a; Ho et al., 2015; Pizzochero et al., 2017; Akbal et al., 2017; Haskell et al., 2018). The superfluid component, however, also has a strong impact on the spectrum of oscillation of the NS, as the increase of degrees of freedom leads to an increase in the number of modes. This in turn will have observational consequences, as modes of oscillation can lead to gravitational wave emission (Andersson & Kokkotas, 1998) and may be at the heart of the observed quasi-periodic oscillations in the tail of magnetar giant flares (Watts & Strohmayer, 2007; Passamonti & Lander, 2014; Gabler et al., 2016).

The hydrodynamics of coupled superfluid and normal fluid systems plays a key role in the analysis of not only NSs, but also laboratory superfluids. In particular it is well known that thermal counterflow of the normal fluid with respect to the superfluid can generate turbulence, and that the vortex array will be disrupted by the so called Donnelly-Glaberson (DG) instability (Glaberson et al., 1974). Many experimental and theoretical studies have recently investigated counterflow turbulence, and the link between the motion of the normal fluid and that of the superfluid (Guo et al., 2010; Kivotides, 2011; Gao et al., 2017; Kivotides, 2018). The instability is a general feature of superfluids with vorticity, and is thus likely to operate in NS interiors (Tsubota et al., 2003; Sidery et al., 2008), as phenomena such as precession and Ekman flow can lead to large scale motions along the rotation axis (Peralta et al., 2005). The DG instability is not stabilized by the magnetic field, and could play an important role in a NS.(van Eysden & Link, 2018).

Very few experiments (see e.g. Swanson et al. 1983; Finne et al. 2003), however, have dealt with instabilities due to counterflow coupled to rotation, which is a setup reminisicent of the neutron star interior, in which the normal component can rotate at a different rate with respect to the neutron superfluid. A notable exception is represented by the experiment of (Tsakadze & Tsakadze, 1980), who considered rotating spherical containers of superfluid He II, that abruptly changed angular velocity, and studied the relaxation of the fluid (see van Eysden & Melatos (2011) for a recent analysis of the experiment).

To make progress in the neutron star problem we adopt a hydrodynamical approach similar to that used in the study of superfluid He II, the so called HVBK two fluid model (Hall & Vinen, 1956; Bekarevich & Khalatnikov, 1961), and based on the multi-fluid neutron star hydrodynamics proposed by Andersson & Comer (2006). This allows us to study modes of oscillation of the coupled fluids and investigate not only the mode spectrum, but also which modes may be driven unstable by the counterflow of the superfluid and normal components. Such instabilities may signal the transition to a turbulent state, and may both trigger a glitch and affect the response of the fluid after the event. An example of such trigger is the the Kelvin-Helmholtz instability on the isotropic-anisotropic interface between the ${}^{1}S_{0}$ and ${}^{3}P_{2}$ neutron superfluids close to the crust-core transition (Mastrano & Melatos, 2005). Another group of instabilities, namely two stream instabilities, may also be related to glitches since the entrainment effect could provide a sufficiently strong coupling for the instability to be astrophysically plausible (Prix et al., 2002; Andersson et al., 2003; Glampedakis & Andersson, 2009; Andersson et al., 2013). Instabilities in pinned superfluids have been studied by Link (2012b, a), who considered a regime in which vortices are pinned either to nuclei in the crust, or to superconducting fluxtubes in the core of the star, with only a small number of them ‘creeping’ out. On a larger scale, the global flow in a NS may also be susceptible to various instabilities in spherical Couette flow (Peralta & Melatos, 2009; Peralta et al., 2006a, 2008), and such instabilities may also be related to timing noise in radio pulsars (Melatos & Link, 2014), and precession of the rotational axis of NSs (Glampedakis et al., 2008, 2009).

In the present work a local three dimensional plane wave analysis of equations of motion for a multiconstituent fluid (Andersson & Comer, 2006) is performed. We consider two fluids, the superfluid neutrons and a normal charge neutral fluid that consist of electrons and protons (Mendell, 1991). Dissipative coupling between the fluids is given by mutual friction, while a non-dissipative coupling is due the entrainment effect, that accounts for the reduced mobility of neutrons especially in the crust (Prix, 2004; Chamel, 2017). We start from the analysis of Sidery et al. (2008), who analysed the oscillations of such a system, assuming the fluids to be locked in the background configuration, and in our calculation also allow for a velocity difference between the constituents in the background, and explicitly account for entrainment. Physically, such a background velocity difference will be built up if vortices are pinned in the crust or core of the star, and is crucial for our analysis.

2 Two fluid hydrodynamics

Our starting point will be the multifluid formalism of Andersson & Comer (2006). We consider the hydrodynamical equations of motion for two dynamical degrees of freedom, the superfluid neutrons and a charge neutral fluid consisting of protons and electrons locked together by electromagnetic interactions on time scales shorter than those of interest for our problem. The momentum equations take the form:

[TABLE]

Indices $\mathrm{x}$ and $\mathrm{y}$ label the constituents, and the inequality $\mathrm{x}\neq\mathrm{y}$ is always understood to be true. The proton-electron fluid will be denoted as $\mathrm{p}$ (and often referred to as the ‘proton’ fluid, as electrons ensure charge neutrality, but only carry a small fraction of the inertia of the fluid) and superfluid neutrons, labeled as $\mathrm{n}$. Indices $i,j,k$ label the spacial coordinates. Summation over repeated indices is implied (excluding the summation over constituent $\mathrm{x}$ and $\mathrm{y}$ indices). The velocity of constituent $\mathrm{x}$ is $v_{i}^{\mathrm{x}}$ while $w_{i}^{\mathrm{y}\mathrm{x}}=v_{i}^{\mathrm{y}}-v_{i}^{\mathrm{x}}$ is the difference between velocities of components. The angular velocity of the star in the background configuration is given by the vector $\Omega^{i}$, and we have included the centrifugal term in the potential $\Phi_{R}$, that in spherical coordinates can be written as

[TABLE]

where $\Phi$ is the gravitational potential which is given by the Poisson equation:

[TABLE]

$\tilde{\mu}_{\mathrm{x}}=\mu_{\mathrm{x}}/m_{\mathrm{x}}$ is the chemical potential per unit mass, and $\rho_{\mathrm{x}}$ is the density of the $\mathrm{x}$ constituent. In the following we make the approximation $m_{\mathrm{p}}=m_{\mathrm{n}}=m$. Finally the momentum per unit mass $\tilde{p}_{i}^{\mathrm{x}}$, which due to entrainment is not aligned with the velocity $v^{i}_{\mathrm{x}}$, is:

[TABLE]

where $\varepsilon_{\mathrm{x}}$ is the entrainment coefficient. Assuming that individual species are conserved the continuity equations are:

[TABLE]

The force $f_{i}^{\mathrm{x}}$ on the right hand side of the equations (1) is the vortex mediated mutual friction (Hall & Vinen, 1956), and it represents an average of individual interactions between vortices and the normal fluid on the sub-hydrodynamical scale. As such its form depends strongly on the properties of the vortex configuration within the fluid. For straight vortices and laminar flow, vortex mediated mutual friction force takes the form (Andersson et al., 2006):

[TABLE]

where hats represent unit vectors, and the vector $\kappa^{i}=\kappa\hat{\Omega}^{i}$ points along the vortex array, which is co-linear with the rotation axis, and $\kappa=h/2m_{\mathrm{n}}$ is the quantum of circulation. The vortex density per unit area, $n_{\mathrm{v}}$, can be linked to the average large scale vorticity $\omega^{i}$ of the superfluid as

[TABLE]

which for two fluids rotating rigidly around a common axis with angular velocities $\Omega_{\mathrm{p}}$ and $\Omega_{\mathrm{n}}$, reduces to

[TABLE]

For vanishing entrainment, or if the two fluids co-rotate, the above expression reduces to the standard Feynman-Onsager relation for rotating superfluids, $\kappa^{i}n_{\mathrm{v}}=2\Omega^{i}_{\mathrm{n}}$.

The mutual friction parameters $\mathcal{B}$ and $\mathcal{B}^{{}^{\prime}}$ in (6) can be expressed in terms of dimensionless drag parameter $\mathcal{R}$ as:

[TABLE]

where $\mathcal{R}$ encodes the microphysics of the dissipation processes that give rise to mutual friction in the stellar interior.

Theoretical calculations provide values for the drag parameter $\mathcal{R}\approx 10^{-4}$ for electron scattering on the vortex cores in the NS core (Alpar et al., 1984b) and $\mathcal{R}\approx 10^{-10}$ for phonon scattering in the crust (Jones, 1990). Higher value $\mathcal{R}\approx 1$ are expected due to Kelvons in the crust (Jones, 1992; Epstein & Baym, 1992), if vortices are moving rapidly past pinning sites. Using the two fluid model of Khomenko & Haskell (2018), that allows for vortex accumulation and differential rotation, Haskell et al. (2018) recently derived values of $\mathcal{R}$ from observations of glitches in Vela and Crab pulsar. They lies in the ranges $\mathcal{B}\approx 10^{-4}-10^{-3}$ for the Vela core and $\mathcal{B}\approx 10^{-5}-10^{-4}$ for the Crab crust. Similar results were obtained by (Graber et al., 2018), who also calculated the density dependence of the Kelvon mutual friction parameter.

In the presence of turbulence, however, the vortex array is disrupted, and a vortex tangle is likely to develop. In this case the form of the mutual friction in (6) will no longer be appropriate. The form of the mutual friction to be used in the case of a polarized turbulent tangle in a neutron star is highly uncertain (Andersson et al., 2007), however in analogy to the cause of homogeneous isotropic turbulence in He-II, we will consider the form proposed by Gorter & Mellink (1949)

[TABLE]

where the parameter $A_{GM}$ has the dimensions of the inverse of a circulation and $w_{\mathrm{x}\mathrm{y}}^{2}=w_{\mathrm{x}\mathrm{y}}^{i}w^{\mathrm{x}\mathrm{y}}_{i}$. As discussed by Andersson et al. (2007), this form is essentially phenomenological and its relevance for neutron stars is still unclear. However, it is interesting to use this form of mutual friction as a completely different, physically motivated, alternative to (6), to investigate the effects of fully developed turbulence on the mode structure of the star.

3 Perturbed equations of motion

We begin our analysis by linearising the equations of motion in a frame rotating with angular velocity $\Omega=\Omega_{\mathrm{p}}$. Formally we expand a generic quantity $Q$ in a neighbourhood of a certain position $\mathbf{r}$, as

[TABLE]

where $Q_{B}$ is the background (time independent) value of the field, and the perturbation is assumed to be such that $|\delta Q|\ll|Q_{B}|$. With this in mind, the perturbed continuity equation follows from (5) and is:

[TABLE]

while from (1) we obtain:

[TABLE]

where $\delta f_{i}^{\mathrm{x}}$ is a perturbed mutual friction contribution. In order to keep our equations tractable we apply the Cowling approximation, taking $\delta\Phi_{R}=0$. Using equation (4), we can thus expand (13) to write

[TABLE]

where

[TABLE]

3.1 Thermodynamical relations

As it is clear from the form of the perturbation equations in the previous section, to compute the effect of a velocity perturbation we need to calculate perturbations of the chemical potential and of the entrainment. For a charge neutral two-fluid system it is known that a complete thermodynamic description in the zero temperature limit can be achieved by considering the two densities $\rho_{\mathrm{x}}$, the squared velocity lag between the two species $w_{\mathrm{y}\mathrm{x}}^{2}$ and the related intensive variables, namely the two chemical potentials per unit mass $\tilde{\mu}_{\mathrm{x}}$ and a dimensionless parameter $\alpha$ which regulates the strength of the conservative entrainment coupling between the two species (Prix, 2004). In particular, the first law of thermodynamics can be expressed as (Andersson & Comer, 2001a; Andersson et al., 2006)

[TABLE]

where ${E}$ is the energy density of the fluid mixture, while the intensive variables are given by

[TABLE]

The parameter $\alpha$ has the dimensions of a mass density and is related to the usual dimensionless entrainment parameters $\varepsilon_{\mathrm{x}}$ via (Prix, 2004; Carter et al., 2006)

[TABLE]

Therefore, the variation $\delta\varepsilon_{\mathrm{x}}$ is related to the variations of the more basic thermodynamic quantities $\alpha$ and $\tilde{\mu}_{\mathrm{x}}$ as

[TABLE]

While for the chemical potential we need to calculate

[TABLE]

The derivatives of chemical potential with respect to the squared background velocity may be computed using the first law of thermodynamics in (16), which leads to:

[TABLE]

which, from (18) gives:

[TABLE]

We are now ready to calculate perturbations of the entrainment and chemical potential. To do this we introduce the sound velocities $c_{\mathrm{x}}$, the chemical couplings $C_{\mathrm{x}}$ (Andersson & Comer, 2001b) and two additional parameters related to the dependence on the lag of $\alpha$ and the chemical potentials,

[TABLE]

Note that our definition of $c_{\mathrm{x}}^{2}$ coincides with the one used by Sidery et al. (2008). In our more general analysis, we let the entrainment vary according to the thermodynamic relations, which forces us to introduce also $A_{\mathrm{x}}$ and $\alpha_{\mathrm{x}}$, at least at this formal level. There are no estimates in the literature for the parameters $A_{\mathrm{x}}$ in a neutron star, so in the following we will neglect it and set $A_{\mathrm{x}}=0$ in practical examples. The parameter $C_{\mathrm{x}}$ is also highly uncertain and generally computed from phenomenological models (Prix et al., 2002), while for $\alpha_{\mathrm{x}}$ one can obtain an estimate in the inner crust, shown in Fig 1, by assuming that the ratio between the densities of the two species is frozen at its chemical equilibrium value. In this case one has that

[TABLE]

where $n_{B}$ represents the total baryon number density. This rough argument gives a plausibility interval for the unknown values of $\alpha$.

Thanks to the fact that $\delta w_{\mathrm{x}\mathrm{y}}^{2}=2(v^{\mathrm{x}}_{i}-v^{\mathrm{y}}_{i})(\delta v_{\mathrm{x}}^{i}-\delta v_{\mathrm{y}}^{i})$, the variation of the chemical potentials can be conveniently written as

[TABLE]

Following the same procedure, we can write for the entrainment:

[TABLE]

which combined with equations (25) and (26) with (19), leads to

[TABLE]

3.2 Perturbing the mutual friction

The next step is to perturb the mutual friction. Let us start by considering the anisotropic form valid for a straight vortex array, given in (6). In this case it is

[TABLE]

where we have assumed that $\delta\mathcal{B}=\delta\mathcal{B}^{{}^{\prime}}=0$. In principle one could vary also these quantities, but given the large uncertainties on these coefficients, and to simplify our analysis, we choose to ignore their perturbations. In order to discuss the terms in a more tractable form, and keep track of their physical origin, we separate the expression above in the following way:

[TABLE]

The first term, $(\delta f_{i})_{sound}$ connects sound waves and the mutual friction and takes the form

[TABLE]

The term, containing the perturbations of the background velocity, which would exist also if the fluids co-rotate in the background, is:

[TABLE]

The final two terms, $(\delta f_{i})_{n_{\mathrm{v}}\kappa}$ and $(\delta f_{i})_{\hat{\kappa}}$ depend on variations of the vorticity, both in magnitude and direction, and to calculate them we need to perturb the vorticity $\omega^{i}=\kappa^{i}n_{\mathrm{v}}$, so that

[TABLE]

which using

[TABLE]

gives

[TABLE]

We also need to perturb the vorticity unit vector $\hat{\kappa}^{i}=\hat{\omega}^{i}$, which can be done readily in terms of a projection operator $\perp$ that projects orthogonally to the vortex lines:

[TABLE]

We can thus write

[TABLE]

and

[TABLE]

which, together with (33) completes our calculation of the perturbed mutual friction. Before moving on, let us remark that we could have obtained our expression directly in terms of the perturbed vorticity $\delta\omega^{i}$, by noting that the term $\epsilon_{ijk}\hat{\omega}^{j}\epsilon^{klm}\omega_{l}{w}_{m}^{\mathrm{y}\mathrm{x}}$ in the expression in (6) for the mutual friction, can be written in terms of the projection operator in (35) as

[TABLE]

and that

[TABLE]

3.2.1 Isotropic mutual friction

We consider also perturbations of the isotropic mutual friction force in (10), proposed for the case in which we have fully developed turbulence and an isotropic vortex tangle. In this case the perturbations of the vorticity play no role and from equation (10) we have:

[TABLE]

4 Plane wave analysis

We are now ready to study the oscillation spectrum of our neutron star model, and to do so we will make the plane wave approximation. To do this we assume that our oscillations can be described as plane waves in a neighbourhood $\mathbf{r}+\mathbf{x}$ of the point $\mathbf{r}$ at which we calculated our background quantities $Q_{B}(\mathbf{r})$, so that

[TABLE]

with

[TABLE]

We assume that $|\mathbf{x}|\ll|\mathbf{r}|$, so that all gradients of background quantities are ignored in our perturbation equations. We use a local cartesian coordinate system, as sketched in figure 2, which is taken to be co-rotating with the normal component, such that $\Omega=\Omega_{\mathrm{p}}$ and $v_{\mathrm{p}}=0$ in the following. This naturally implies that

[TABLE]

Each perturbed quantity $\delta Q$ can be expressed in terms of the two fields $\delta v^{i}_{\mathrm{n}}$ and $\delta v^{i}_{\mathrm{p}}$ by means of the continuity equations and of the thermodynamic relations and it will be convenient to express the amplitude $\bar{Q}$ as:

[TABLE]

Following this procedure, the perturbed version of the Euler-like equations (1) can be used to define a linear system of six equations for the six independent amplitudes $\bar{v}_{\mathrm{x}}^{i}$ which has the form

[TABLE]

The only way to have non trivial solutions of the above equation is to impose that $\det{(M)}=0$. The determinant of (45) is in general a high degree polynomial in the components of the wave vector $\mathbf{k}$ and of the frequency $\omega$ and requires a numerical analysis. Particular cases that can be studied analytically will, however, be presented in the following sections.

To calculate the components of the matrix in (45) first of all we insert the plane wave ansatz in the continuity equation (12), to obtain

[TABLE]

We also need the plane wave perturbations of the chemical potential and entrainment, which from (25) and (27), are:

[TABLE]

and

[TABLE]

With the above expressions at hand we can now compute the plane wave perturbations of the Euler equations in (14). The full result is complex, but a procedure to obtain all the variations in the presence of a background lag and with no approximations is given in Appendix A. It is, however, instructive to study the case in which we neglect perturbations of the entrainment and set $\delta\varepsilon_{\mathrm{x}}=0$ and also neglect terms involving $\partial\varepsilon_{\mathrm{x}}/\partial\rho_{\mathrm{x}}$, so that one has $\alpha_{\mathrm{x}}\approx\varepsilon_{\mathrm{x}}$ in (47). In this case the Euler equations can be written as

[TABLE]

where defining $\bar{p}_{\mathrm{n}}^{k}=\bar{v}_{\mathrm{n}}^{k}+\varepsilon_{\mathrm{n}}\bar{w}_{\mathrm{p}\mathrm{n}}^{k}$, we have

[TABLE]

5 Waves in the absence of mutual friction

In this section we consider some simple examples where progress can be made analytically, before moving on to the numerical results for the full set of equations. In particular we will first consider sound waves and inertial waves in the limit of vanishing mutual friction, and show how mutual friction can damp or drive some of the modes instable. The approximations will be outlined for each case as we proceed.

5.1 Sound waves

Let us begin by considering pure sound waves. If we set $\Omega=\varepsilon_{\mathrm{x}}=w_{\mathrm{n}\mathrm{p}}=v_{\mathrm{n}}=0$ in (49) and (50) the dispersion relation follows from

[TABLE]

we recover simply two branches of sound waves, one for each fluid

[TABLE]

Including a background flow, still with $\Omega=\varepsilon_{\mathrm{x}}=0$ does not couple the fluids, and simply introduces a correction:

[TABLE]

5.1.1 Chemical coupling and entrainment

We now move on to considering the chemical coupling between the two fluids in the regime in which $\varepsilon_{\mathrm{x}}=w_{\mathrm{n}\mathrm{p}}=v_{\mathrm{n}}=0$ and $\Omega=0$. In this case the dispersion relation follows from

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

as well as the dimensionless parameter

[TABLE]

The roots of (57) are organized in four branches that have the form

[TABLE]

where the constants $c_{1}$ and $c_{2}$ are two velocities for entrainment-coupled sound waves,

[TABLE]

with

[TABLE]

It is immediate to see that, in the limit of no entrainment and no chemical coupling, $c_{1}$ and $c_{2}$ are exactly equal to $c_{n}$ and $c_{p}$ of pure one-component sound waves. To grasp the physical effect of chemical coupling we furthermore require that $c_{\mathrm{x}}^{2}\gg C_{\mathrm{x}},C_{\mathrm{y}}$. In this limit it is useful to define $\Delta_{0}$ as the limit of $\Delta$ when $C_{\mathrm{x}}=0$,

[TABLE]

so that the the limit of weak chemical coupling can be equally seen as $C_{x}/\Delta_{0}\ll 1$. In this case, the velocities in (69) can be expanded up to the second order in the chemical couplings as

[TABLE]

In the limit of $\varepsilon_{\mathrm{n}}=\varepsilon_{\mathrm{p}}$, and using the relation $C_{\mathrm{p}}=\frac{\rho_{\mathrm{p}}}{\rho_{\mathrm{n}}}C_{\mathrm{n}}$ we recover the result (Sidery et al., 2008):

[TABLE]

However, retaining the entrainment terms, the lowest order correction to the dispersion relation comes the term that is linear in $C_{\varepsilon}$: when entrainment is considered, the effect of the chemical coupling on sound waves is enhanced. Note that this result is independent of the assumption $\delta\varepsilon_{\mathrm{x}}=0$, and remains the same also in the more general case in which perturbations of the entrainment are considered.

5.1.2 Including rotation

We now consider the effect of rotation on the sound waves, in the limit of no chemical coupling, i.e. we take $C_{\mathrm{p}}=C_{\mathrm{n}}=w_{\mathrm{n}\mathrm{p}}=v_{\mathrm{n}}=0$ in the background. In this case the dispersion relation is still too complex to solve for a general case. However if we consider sound waves with a dispersion relation of the form $\omega\sim c_{\mathrm{x}}k$, we expect that $\omega\gg\Omega$ even for the smallest possible $k$, which is of the order of the inverse of the stellar radius $R^{-1}$ (i.e. the sound velocity in nuclear matter is such that $c_{\mathrm{x}}\gg\Omega R$). Therefore, it is a very good approximation to expand the two branches of $\omega^{2}(k)$ to the lowest order in $\Omega$:

[TABLE]

where

[TABLE]

where $c_{1}$ corresponds to the upper sign choice. The quantities $f_{1}$ and $f_{2}$ are

[TABLE]

where $f_{1}$ corresponds to the upper sign choice and $\theta$ is the angle between $\boldsymbol{\Omega}$ and $\mathbf{k}$. Finally, combining the above expressions, we find that the correction to the four branches of the dispersion relation for sound waves produced by the Coriolis force is

[TABLE]

which, neglecting entrainment ($\varepsilon_{\mathrm{n}}=\varepsilon_{\mathrm{p}}=0$) reduces to (Sidery et al., 2008):

[TABLE]

5.2 Inertial waves

We now continue to focus on the effect of rotation, but search for solutions that are inertial waves, for which the dispersion relation is not linear in $k$. If we set $c_{\mathrm{n}}=c_{\mathrm{p}}=C_{\mathrm{n}}=C_{\mathrm{p}}$ we find, as expected, families of inertial modes such that

[TABLE]

In this case we still have two families of modes, but these do not correspond to oscillations in the single fluids as in the sound wave case, but rather represent an inertial mode where the two fluids flow together, which represents the standard inertial mode of the system, and a mode in which the two fluid oscillate relative to each other, the frequency of which is affected by entrainment (Haskell et al., 2009).

6 Including mutual friction

So far it has been assumed that the two fluids can oscillate without any dissipation mechanism. The mutual friction will couple the two fluids and generally tend to damp any relative motion. To study this problem we consider a selection of analytically tractable cases, and defer the reader to section 7 for the full analysis. We begin with a discussion of inertial waves, as it is known that inertial modes may be dynamically unstable in the presence of a background flow (Sidery et al., 2008; Prix, 2004; van Eysden & Link, 2018).

Let us start our analysis of inertial waves by considering the weak drag regime, such that $\mathcal{R}\ll 1$, which results in $\mathcal{B}^{{}^{\prime}}\ll\mathcal{B}$. We thus take $\mathcal{B}^{{}^{\prime}}=0$, and also set $c_{\mathrm{n}}=c_{\mathrm{p}}=C_{\mathrm{n}}=C_{\mathrm{p}}=0$ as we are interested in inertial waves. In the case in which there is no background flow, i.e. $w^{i}_{\mathrm{n}\mathrm{p}}=v^{i}_{\mathrm{n}}=0$, the dispersion relation follows from

[TABLE]

where $\tilde{\rho}=\rho_{\mathrm{n}}/\rho_{\mathrm{p}}$. The solutions are still an undamped inertial mode, which we identify with the co-moving mode, as in this case there is no relative motion and no mutual friction (note that in a full spherical analysis rotation couples the two motions at higher order, and results in mutual friction damping even for the standard co-moving mode, see Haskell et al. 2009), and the damped counter-moving mode, modified by entrainment.

[TABLE]

[TABLE]

Note that the denominator is always positive, $1-\varepsilon_{\mathrm{n}}-\varepsilon_{\mathrm{p}}>0$, as required by stability on a microscopic scale (Carter et al., 2006). The mode is therefore always stable, but we see here that our study of large scale hydrodynamical instabilities, such as the one that would arise if $1-\varepsilon_{\mathrm{n}}-\varepsilon_{\mathrm{p}}<0$, successfully capture a more general physical instability of the system.

6.1 Background flow

We now consider the physically realistic case in which the fluids are not co-rotating in the background, for example due to vortex pinning in the crust. The analysis is now more involved, so we will consider two cases separately: counterflow along the axis of the vortex array (the DG instability), and counterflow perpendicular to the axis.

6.1.1 The Donnelly-Glaberson instability

To study the instabilities that arise due to counterflow along the vortex axis we specialize our setup, and consider only modes propagating along the vortex axis, which is taken to be alligned with the rotation axis along $z$. We thus have $k_{x}=k_{y}=0$. To begin our analysis we consider the simplified case with no entrainment, i.e. $\varepsilon_{\mathrm{x}}=\varepsilon_{\mathrm{y}}=0$ and also start from the assumption $\tilde{\rho}=\rho_{n}/\rho_{p}=0$, which corresponds to assuming that the mutual friction only acts on the neutrons. This not only simplifies the calculations, but it corresponds to assuming that the protons and electrons are ’clamped’ to the normal fluid (e.g. their motion is entirely dictated by the magnetic field). In this case the dispersion relation is:

[TABLE]

Which, apart from the trivial solutions corresponding to the choice of coordinate system, has solutions:

[TABLE]

The first solution corresponds to standard inertial waves (the factor $\cos\theta=1$ due to our choice of direction of $\mathbf{k}$), the second corresponds to the DG instability and is generally unstable as long as $k_{z}v^{z}>0$, i.e. $|v_{\mathrm{n}}|>2\Omega/|k_{z}|$. If we consider a standard pulsar such that $2\Omega\approx 100$ rad/sec and hydrodynamical scales such that $|k_{z}|<1$ cm*-1*, then we have an instability if $|v_{\mathrm{n}}|>10^{2}$ cm/s, which corresponds to a lag of $\Delta\Omega\approx 10^{-4}$, in the outer layers of the star, which is easily sustainable by pinning forces in the crust (Seveso et al., 2016). Such an instability is thus always likely to be present, as vortex bending and large scale flows in the normal fluid, coupled to the superfluid, will always induce counterflow along the vortex array. The fact that the instability exists even if $\tilde{\rho}=0$ and mutual friction is not acting on the protons, suggests the instability will exist even in the presence of strong magnetic fields, as suggested by the analysis of van Eysden & Link (2018).

If we relax our approximations, and take $\tilde{\rho}\neq 0$ and include entrainment in our analysis the full solution is intractable. We thus consider the weak drag case, for which $\mathcal{R}\ll 1$, and we can take $\mathcal{B}^{{}^{\prime}}\approx 0$. Furthermore we consider two limiting cases, the small entrainment limit, which is relevant for the core of the star, and the large entrainment limit, relevant for the crust. For small entrainment $\varepsilon_{\mathrm{n}}<\varepsilon_{\mathrm{p}}\ll 1$, we have

[TABLE]

and

[TABLE]

We see that in the presence of a background flow entrainment modifies also the co-moving mode, which is now potentially unstable. For weak drag, if we expand to first order in $\mathcal{B}$, we see that the mode is always unstable for

[TABLE]

which means that in the presence of a lag, due to core pinning (e.g. of neutron vortices to proton fluxtubes), the vortex array is generally unstable on all hydordynamical lengthscales for standard parameters ($\mathcal{B}\approx 10^{-4}$, $\Omega\approx 100$, $\varepsilon_{\mathrm{p}}=0.6$, $|v_{\mathrm{n}}|=10^{4}$ cm/s).

For the crust of the neutron star we can make the large entrainment approximation, we to simplify our calculation we take to be the limit $\varepsilon_{\mathrm{x}}\longrightarrow\infty$. In this case the dispersion relation is:

[TABLE]

with:

[TABLE]

from which we can see that one of the modes is always generically unstable in the limit $|k_{z}v^{z}_{\mathrm{n}}|/\Omega\ll 1$, which we expect to be satisfied in most cases in a NS interior.

6.1.2 Two stream instability

Let us now consider the case of a two stream instability, in which there is a counterflow perpendicular to the vortex axis, which we take to be aligned with the rotation axis such that $\hat{\kappa}^{i}=\hat{\Omega}^{i}=\hat{z}^{i}$. We thus have $v_{x}^{\mathrm{n}}\neq 0$, with $v_{y}^{\mathrm{n}}=0,v_{z}^{\mathrm{n}}=0$, and to keep the calculation tractable we consider the case in which $\varepsilon_{\mathrm{x}}=0$ and take the protons to be clamped to the normal fluid($\tilde{\rho}=0$). The dispersion relation is:

[TABLE]

Apart from the trivial solution we now have solutions of the form

[TABLE]

This mode can be unstable due to the induced flow in the $y$ direction caused by the background flow in the $x$ direction. An expansion for $\mathcal{B}\ll 1$ reveals that the criterion for instability is:

[TABLE]

This gives the condition for instability to develop:

[TABLE]

Note, this condition doesn’t contain dependence on $k_{y}$. To obtain this we need expand to the second order in small parameter $|v_{\mathrm{n}}|/\Omega$, which leads to:

[TABLE]

As we can see mutual friction doesn’t affect the critical velocity for the instability to set in, but affects its growth rate, with higher mutual friction causing a faster rise of the instability with increasing velocity.

6.2 The strong drag regime

We now move on to consider the case of a strong drag, that corresponds to $\mathcal{R}\gg 1$, in which case the parameter $\mathcal{B}^{{}^{\prime}}$ can no longer be neglected. This case may be applicable in regions with strong pinning (Ruderman et al., 1998; Link, 2006, 2003), and may have important consequences for the development of the r-mode instability (Haskell et al., 2009; Glampedakis & Andersson, 2009). Let us begin our analysis by the equivalent of the DG instability, in which the relative flow is along the vortex axis, taken to coincide with rotational axis, or $v_{x}^{\mathrm{n}}=v_{y}^{\mathrm{n}}=0$ and $v_{z}^{\mathrm{n}}\neq 0$, and we also take $\varepsilon_{\mathrm{x}}=0$. We consider the limit $\mathcal{R}\ll 1$ for which $\mathcal{B}\approx 0$ and $\mathcal{B}^{{}^{\prime}}\approx 1$. In this case the spectrum is given by:

[TABLE]

with

[TABLE]

For standard neutron star parameters these modes are not unstable, and in fact if we expanding far $B^{{}^{\prime}}\ll 1$ we obtain:

[TABLE]

Which are stable, modified inertial modes.

6.2.1 Sound waves and mutual friction

To study the effect of mutual friction on sound waves we can consider an expansion in $\Omega$ around the $\Omega=0$ solution $\omega=\pm c_{\mathrm{n}}k$. Let us consider the modes of the superfluid, in the clamped proton approximation. Defining $\theta$ to be the angle between the wave vector and the rotation axis, we find, for $\varepsilon_{\mathrm{x}}=C_{\mathrm{x}}=w_{\mathrm{x}\mathrm{y}}=v_{\mathrm{n}}=0$, we have:

[TABLE]

The effect of the mutual friction disappears when the wave vector is aligned with the vorticity, i.e. $\theta=0$: in this case the velocity perturbation is parallel to the vortex array and there cannot be any mutual friction effect (all the cross products in 28 are zero). Things are very different if a background velocity lag is present (even if this lag is parallel to the vortex array) because vorticity perturbations enter the game: this will be investigated numerically in the next sections.

In the weak drag limit, we may set $B^{\prime}\approx\mathcal{R}^{2}$ and $\mathcal{B}\approx\mathcal{R}$ and the relevant terms become

[TABLE]

6.3 Isotropic (Gorter-Mellink) mutual friction

We now consider the Gorter-Mellink form for the mutual friction in (10) to investigate if and how the presence of an isotropic vortex tangle modifies the previous results. In the limit $kc_{\mathrm{n}}\gg A_{GM}|\mathbf{v}_{n}|^{2}$ it is possible to write the exact implicit form of the dispersion relation, which is:

[TABLE]

for inertial and sound waves respectively. Here $f^{\pm}_{I}\sim 1$ and $f^{\pm}_{S}\sim 1$ are four involved functions of the angles between the vectors $\mathbf{v}_{n}$, $\mathbf{k}$ and $\boldsymbol{\Omega}$, as well as the two positive functions $\alpha^{\pm}_{I}>0$ and $\alpha^{\pm}_{S}>0$. Irrespectively of the mutual orientation and magnitude of the three vectors (provided that $kc_{n}$ is always much bigger than $\Omega$ and $A_{GM}|\mathbf{v}_{n}|^{2}$), we are always in a damped regime in which the damping timescale is of the order of $A_{GM}^{-1}|\mathbf{v}_{n}|^{-2}$ for both inertial and sound waves.

It would thus appear that a rectilinear vortex array is rapidly destabilized as unstable inertial modes develop due to counterflow. If an (approximately) isotropic vortex tangle develops this remains stable, at least until the turbulent tangle decays. Once a rectilinear array is restored the system is again unstable, leading to a recurrent mechanisms that may be linked to the trigger of pulsar glitches. A more detailed analysis in the case of a polarized turbulent tangle should be the focus of future work.

Note that as the modes in (98) are stable, we will not consider isotropic mutual friction in the numerical analysis in the following section.

7 Numerical results

All the machinery needed to derive the general dispersion relations for the various modes of oscillation of the two-fluid system is presented in Appendix A. The final result is the matrix $M$ in (125), which depends on some basic local quantities that define how the system responds to a perturbation in the velocity fields: the full matrix (or its determinant) can be seen as a function of $\omega$ and $\mathbf{k}$ and depends on several parameters, namely

[TABLE]

Leaving aside the obvious dependence on the pulsation $\omega$ and on $\mathbf{k}$, the other parameters (that describe completely the hydrodynamic state of the background configuration when there is no turbulence and magnetic field) have been formally divided into three sets, which helps us to discuss their role.

The first set comprises the variables $\mathbf{v}_{n}$, $\boldsymbol{\omega}$ and $\boldsymbol{\Omega}$: they define the state of motion of the background configuration in which the normal component rotates rigidly. We allow for a stationary non-zero local velocity lag $\mathbf{v}_{n}=\mathbf{w}_{\mathrm{n}\mathrm{p}}$, which is expected to be non constant on the stellar radius length scale. Therefore, in a purely local analysis, the background vorticity field $\boldsymbol{\omega}$ is not expressible in terms of the local (constant) value of $\mathbf{w}_{np}$: only a true solution of the unperturbed equations of motion would lead to locally consistent values for $\mathbf{w}_{\mathrm{n}\mathrm{p}}$, $\boldsymbol{\omega}$ and $\boldsymbol{\Omega}$, but would prevent us to test the more general case in which the local direction and magnitude of these vectors are chosen at will.

The second set of parameters describes the local state of matter in the stationary configuration. It can be further divided into two subsets. The four parameters $\varepsilon_{\mathrm{x}}$ and $\rho_{\mathrm{x}}$, are defined by considering the first law of thermodynamics (16) and their unperturbed value can oscillate: indeed, we introduced the linear combinations (46) and (48) to express their amplitudes in terms of the perturbed velocities (the same is valid for the chemical potentials $\tilde{\mu}_{\mathrm{x}}$, with the only caveat that the unperturbed value of the chemical potentials do not enter into the explicit expression of $M$). The second subset comprises $\alpha_{\mathrm{x}}$, $A_{\mathrm{x}}$, $C_{\mathrm{x}}$ and $c_{\mathrm{x}}$, that are related to second order derivatives of the internal energy $E$: in a first-order analysis these quantities are fixed to their unperturbed value and allow, together with the unperturbed values $\rho_{\mathrm{x}}$ and $\varepsilon_{\mathrm{x}}$, to express the fundamental thermodynamic perturbations $\delta\rho_{\mathrm{x}}$, $\delta\varepsilon_{\mathrm{x}}$ and $\delta\tilde{\mu}_{\mathrm{x}}$ as described in section 3.1.

Finally, the two parameters $\mathcal{B}$ and $\mathcal{B}^{\prime}$ define the “state” of the vortex array (note again, that we only consider the standard, anisotropic form of the mutual friction, and not the Gorter Mellink form, as we have found modes to be stable for this form of the mutual friction). We do not consider variations of these parameters and their value has to be fixed according to the mutual friction scenario that we are interested to test. In the subsequent numerical analysis we consider different possibilities, listed in table 2.

From the analytic point of view, writing down the full matrix $M$ is of little interest but its complete and explicit form in components can be easily obtained with the aid of any software for symbolic computation, as outlined in appendix A. Although extremely complex, also the full and exact determinant of $M$ can be obtained as well. At this point we substitute a particular choice for the background parameters into the complete expression of the determinant; the roots of $\det M(\omega,\mathbf{k})=0$, which is a high degree polynomial equation in the pulsation $\omega$, can be computed numerically for different values of $\mathbf{k}$. By varying $\mathbf{k}$ and the parameters, some branches of the dispersion relation $\omega(\mathbf{k})$ may result in a positive imaginary part. The associated instability timescale is then defined as

[TABLE]

Since we worked within a purely hydrodynamic framework, it is important to stress that not all the values of $|\mathbf{k}|$ are physically meaningful. As also discussed in van Eysden & Link (2018), the wave vector should be much smaller than the inverse of the stellar radius, otherwise the effect of stratification is expected to modify the dispersion relation (namely, the background quantities have spatial dependence on such length scales and cannot be considered uniform, as in the present local analysis). Since the stellar radius is about $R\sim 10\,$km, we should consider $|\mathbf{k}|\gg 10^{-6}\,$cm*-1*: therefore, we expect the present analysis to be valid for $|\mathbf{k}|\sim 10^{-3}\,$cm*-1* since significant density changes are expected to occur on the length scale of about ten meters, especially at the core-crust interface.

On the other hand, the definition of the superfluid momentum and vorticity needs a suitable average on a macroscopic sample of matter containing many vortex lines, which average separation is expected to be of the order of $\sim 10^{-3}\,$cm. Therefore, we expect the present analysis to break down when the wave vector reaches the critical value $|\mathbf{k}|\sim 10^{3}\,$cm*-1*. For this reason, the region outside the physically interesting range $10^{-3}<|\mathbf{k}|\,$cm$\,<10^{3}$ is shaded in the plots of the instability timescales.

Considering the relative orientation of four different vectors ($\mathbf{v}_{n}$, $\boldsymbol{\Omega}$, $\boldsymbol{\omega}$ and $\mathbf{k}$) leads to a huge parameter space; therefore, in the numerical analysis we stick to the case in which $\boldsymbol{\Omega}=\Omega\hat{\mathbf{z}}$ and $\boldsymbol{\omega}=2\boldsymbol{\Omega}$, while $\mathbf{v}_{n}$ and $\mathbf{k}$ are chosen to be aligned with one of the three directions $\hat{\mathbf{x}}$, $\hat{\mathbf{y}}$ or $\hat{\mathbf{z}}$, for a total of nine combinations111 The presence of a non zero background lag may locally modify the vorticity direction. However, in our local analysis the background lag is treated as locally uniform, while the local direction of the vorticity depends on the large scale, global arrangement of the background lag in the stationary configuration: we are therefore forced, due to the fact that the analysis is local, to assume three independent directions for the wave vector, the velocity lag and the vorticity. However, in the following analysis, we keep the vorticity fixed along the z-axis (i.e. aligned with $\mathbf{\Omega}$). This assumption, as well as the fact that we will use $|\mathbf{\omega}|=2\Omega$ in the background configuration, is not severe when the lag in angular velocity is small compared to the rotation rate of the star. Note that only the background vorticity is fixed, we account for variations in $|\boldsymbol{\omega}|$ in the mutual friction. . However, for simplicity, only the interesting cases that allow for unstable modes are shown in figures and discussed.

First of all we consider the DG instability, for strong and weak drag and in the ‘pinned’ scenario, which we mimic by taking $\mathcal{B}^{{}^{\prime}}=1$, but $\mathcal{B}=0$. In Figures 3 and 5 we plot the instability timescale for strong drag (weak drag is qualitatively similar) and for pinned vortices, for varying values of the background lag. In general we confirm the results of the analysis in section 6.1.1, there exists a family of mixed inertial and sound waves, that are unstable on dynamical timescales, both in the core and crust. Mutual friction has little effect on damping the instability (in fact it is at the heart of driving it), and rather it is the lag that plays a role in determining the onset of the instability. For small values of the lag of $v_{\mathrm{n}}^{z}=1$ cm/s, corresponding to $\Delta\Omega\approx 10^{-5}$ in the outer core or inner crust, assuming a radius of $R\approx 10$ km, only the shortest length-scales are unstable. For higher lags ($v_{\mathrm{n}}^{z}=10^{4}$ cm/s, $\Delta\Omega\approx 10^{-2}$) all the dynamical range is unstable, and we can thus assume that once pinning allows for a significant enough lag to develop, the array will go unstable, possibly playing a role in triggering pulsar glitches or contributing to timing noise. In Figures 4 and 6 we plot the dispersion relation for both the strong drag and pinned case. We see that the nature of the unstable modes changes: in the pinned case inertial modes are unstable, while in the strong drag case (and we have verified that the same is true for weak drag) the unstable mode is a modified sound wave.

In Figure 7 we consider the two-stream instability, in which the background counterflow is taken perpendicular to the vortex (and rotation) axis. Here again we confirm our analytical results from section 6.1.2. There are mixed sound-inertial waves (as can be verified from the dispersion relation in Figure 8) that are unstable on dynamical timescales and the unstable range depends strongly on the magnitude of the background counterflow velocity. For high values, as may be expected from strong pinning, $v_{\mathrm{n}}^{x}=10^{4}$ cm/s, the whole dynamical range is unstable, with the instability developing on similar timescales at all scales. For $v_{\mathrm{n}}^{x}=1$ cm/s, on the other hand, only the smallest length-scales are unstable, and the instability develops on longer timescales, that may be affected by other viscous mechanisms, such as shear viscosity.

8 Conclusions

We have studied the modes of oscillation of a superfluid neutron star, accounting for background counterflows and including entrainment. We have considered the effect for both standard HVBK mutual friction that arises if the vortex array is straight, and the Gorter-Mellink isotropic form, relevant for fully developed turbulence.

We find that for standard mutual friction there is always a fast instability in the case of counterflow along the vortex axis (which may arise, for example if the star is precessing, or if vortices bend on a large scale), which is the neutron star analogue of the Donnelly-Glaberson instability that is observed in laboratory studies of superfluid helium. We also confirm the existence of a rapid two-stream instability when the background counterflow is perpendicular to the vortex axis (as is expected if a difference in angular velocity between the superfluid and the normal fluid arises due to pinning).

We find that entrainment plays an important role in the instability, and its effect is twofold. On the one-side it extends the dynamical range over which inertial modes are unstable in the crust, and but on the other it significantly shortens the growth time in the core of the star.

When we consider isotropic Gorter-Mellink mutual friction we find that the previous instabilities are stabilized. This suggests that once a large enough lag sets in, the instabilities we have presented will disrupt the straight vortex array and a turbulent tangle will develop. At this point the system is stable until the turbulence decays and the process will start again. It is thus very likely that transitions to turbulence play a role in triggering pulsar glitches and in timing noise (Peralta et al., 2006b; Melatos & Link, 2014).

An important effect that we have not considered is that of the magnetic field. The problem is complicated, as the interior field configuration of a neutron star is generally unknown, but van Eysden & Link (2018) found that for simple field geometries the two-stream instability can be stabilized, while the Donnelly-Glaberson instability is always present. We have mimicked the effect of the magnetic field by studying a setup in which the proton fluid is ‘clamped’ and not affected by the mutual friction (to imitate the situation in which it is held in place by magnetic stresses) and find that the two-stream instability is still present in the weak drag regime when entrainment is present. Future work should focus on the full problem, accounting also for superconductivity in the core and the interaction between neutron vortices and superconducting flux tubes (Glampedakis et al., 2011).

Acknowledgments

We acknowledge support from the Polish National Science Centre grant SONATA BIS 2015/18/E/ST9/00577. Partial support comes from PHAROS, COST Action CA16214. We thank the Institute for Nuclear Theory at the University of Washington for its kind hospitality and stimulating research environment. This research was supported in part by the INT’s U.S. Department of Energy grant No. DE-FG02- 00ER41132.

Appendix A Derivation of the dispersion matrix

In this appendix we provide a more self-contained description of how to perform the local variations by using directly the plane wave ansatz and find the dispersion matrix of the two-fluid system. In the following, the indexes $\mathrm{x}$ and $\mathrm{y}$ are always meant to be different ($\mathrm{x}\neq\mathrm{y}$), while $\mathrm{x}$ and $\mathrm{y}^{\prime}$ can assume the same value.

For clarity, let us list all the quantities that appear into the Euler equations (1) as the entries of a vector $Q$, namely

[TABLE]

Note that the gravitational potential $\Phi$ is not included into the set of variables that we perturb, as well as the mutual friction parameters $\mathcal{B}$, $\mathcal{B}^{\prime}$ and $A_{GM}$. The angular velocity $\Omega$ is also kept constant and defines the rotating frame in which $\bf k$ and the pulsation $\omega$ of the plane wave ansatz (42) are measured. Since the quantities in $Q$ are not all independent, it is convenient to calculate their variation following the order in the list. Once the quantities in $Q$ have been expressed in terms of the velocity fields only, the vector is expanded at the frist order as

[TABLE]

Hence, the generic amplitudes $\bar{Q}$ introduced in equation (44) are obtained by means of the coefficients

[TABLE]

For $\rho_{x}$, $\tilde{\mu}_{\mathrm{x}}$ and $\varepsilon_{\mathrm{x}}$ the result has already been given in section 4. From the contintuity equation we find the diagonal elements $\rho^{\mathrm{x}\mathrm{x}}_{j}$ of the three $2\times 2$ matrices $\rho^{\mathrm{x}\mathrm{y}^{\prime}}_{j}=\delta_{\mathrm{x}\mathrm{y}^{\prime}}\rho^{\mathrm{x}\mathrm{x}}_{j}$, see equation (46). Now, the thermodynamic relations of section 3.1 allow us to calcualte the explicit form of

[TABLE]

and

[TABLE]

where the matrices $\mu^{\mathrm{x}\mathrm{y}^{\prime}}_{j}$ and $\varepsilon^{\mathrm{x}\mathrm{y}^{\prime}}_{j}$ have been given in equations (47) and (48) respectively. Using the result for $\rho^{\mathrm{x}\mathrm{y}^{\prime}}_{j}$, the final expressions read

[TABLE]

for the specific momentum and

[TABLE]

for the entrainment parameters.

At this point it is immediate to address also the variations of the momenta and of the vorticity. To give a concrete example of how the linear combination (44) looks like when applied to vectorial quantities, we start by formally writing the momentum amplitude as

[TABLE]

Performing the expansion (102) on the momenta defined in (4) immediately gives

[TABLE]

The last term can be expanded thanks to (105), and (102) tells us that

[TABLE]

When the entrainment coefficients (107) are inserted into the above formulas, we see that the momentum perturbations have a correction which depends on the coefficient $A_{\mathrm{x}}$. Not surprisingly, these terms are of the second order in the background velocity lag; it is not particularly convenient to neglect them here, but from the numerical point of view these terms are expected to have little effect on the calculated dispersion relation.

The last ingredient that we need is to write the vorticity perturbation in the same fashion of (44), namely

[TABLE]

Starting from equation (34), or applying directly (102) on the definition of vorticity together with (105), it is straightforward to check that

[TABLE]

Alternatively, the above result can be found directly from (109) by observing that $\bar{\omega}_{i}=i\,\epsilon_{iab}\,k^{a}\,\bar{p}_{n}^{b}$. Similarly, the variation of the vorticity versor can be easily obtained by means of the projector in equation (35).

A.1 Left hand side of the Euler equations

Let us define the right hand side of equations (1) as

[TABLE]

where the gravitational potential $\Phi$ and the centrifugal term have been ignored: we are interested in the variation $\delta\mathcal{E}_{\mathrm{x}}^{i}$ and these two terms only contribute to define the background configuration. We extend the notation (42) to $\mathcal{E}^{\mathrm{x}}_{i}$, namely

[TABLE]

where

[TABLE]

Using systematically the results obtained so far, a small amount of algebra gives

[TABLE]

and

[TABLE]

In a dynamical regime in which the mutual friction is zero (i.e. in the rather ideal situation in which the vortex lines can be considered free and no drag nor pinning interactions act upon them), this result is sufficient to build the matrix $M$ of equation (45). More explicitly, the six equations of the full linear system (45) can be also be written as

[TABLE]

where the indexes $i$ and $a$ run over the three spatial indexes (the summation over $a$ is understood). Therefore, it should be clear that each block of the full matrix $M$ can be obtained as

[TABLE]

in the free vortex limit.

A.2 Perturbing the mutual friction

In order to treat consistently the mutual friction term it is convenient to start from the case $\mathrm{x}=n$ in (6), namely

[TABLE]

Once the variation $\delta\mathbf{f}^{n}$ has been obtained, the reaction on the normal component will just be given by

[TABLE]

Note that this last equation is valid also for the Gorter-Mellink mutual friction (10). In complete analogy with what has already been discussed for equation (115), we may write

[TABLE]

The calculation leading to an explicit form of the coefficients ${f}^{\mathrm{x}\mathrm{y}^{\prime}}_{ia}$ is laborious and the final expression is complex enough to be useless in an analytic approach. However, it is possible to implement an exact procedure that can be solved by any software for symbolic calculus. First, the amplitude in (123) can be formally obtained as

[TABLE]

where the derivative is evaluated on the background configuration. Therefore, the matrices in (123) are computed as

[TABLE]

In this way, it is possible to find the complete matrix $M$ of equation (99): not surprisingly, its form is defined in terms of the four $3\times 3$ blocks

[TABLE]

The determinant of this $6\times 6$ matrix is a huge rational function of complex coefficients, which numerator defines a high degree polynomial in $\omega$ and in the components of $\mathbf{k}$. The roots of this polynomial define the dispersion relation $\omega(\mathbf{k})$ of the oscillation modes of the two-fluid system. In section 7 the dispersion relations relative to some physically interesting cases are studied numerically within the symplifing assumptions that $C_{\mathrm{x}}=\alpha_{x}=A_{\mathrm{x}}=0$. Notice that, according to (107), this does not imply that entrainment variations are null (as it has been assumed in particular in (49) and (50)).