# Global numerical simulations of the rise of vortex-mediated pulsar glitches in full general relativity

Aurélie Sourie, Nicolas Chamel, Jérôme Novak, Micaela Oertel

► To cite this version:

Aurélie Sourie, Nicolas Chamel, Jérôme Novak, Micaela Oertel. Global numerical simulations of the rise of vortex-mediated pulsar glitches in full general relativity. Monthly Notices of the Royal Astronomical Society, 2017, 464 (4), pp.4641-4657. 10.1093/mnras/stw2613. hal-01554210

HAL Id: hal-01554210

https://hal.science/hal-01554210v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Global numerical simulations of the rise of vortex-mediated pulsar glitches in full general relativity

A. Sourie,<fnref n="1" />★ N. Chamel,<fnref n="2" /> J. Novak<fnref n="1" /> and M. Oertel<fnref n="1" />

Accepted 2016 October 7. Received 2016 October 7; in original form 2016 July 27

## ABSTRACT

In this paper, we study in detail the role of general relativity on the global dynamics of giant pulsar glitches as exemplified by Vela. For this purpose, we carry out numerical simulations of the spin up triggered by the sudden unpinning of superfluid vortices. In particular, we compute the exchange of angular momentum between the core neutron superfluid and the rest of the star within a two-fluid model including both (non-dissipative) entrainment effects and (dissipative) mutual friction forces. Our simulations are based on a quasi-stationary approach using realistic equations of state (EoSs). We show that the evolution of the angular velocities of both fluids can be accurately described by an exponential law. The associated characteristic rise time $\tau_r$, which can be precisely computed from stationary configurations only, has a form similar to that obtained in the Newtonian limit. However, general relativity changes the structure of the star and leads to additional couplings between the fluids due to frame-dragging effects. As a consequence, general relativity can have a large impact on the actual value of $\tau_r$: the errors incurred by using Newtonian gravity are thus found to be as large as ~40 per cent for the models considered. Values of the rise time are calculated for Vela and compared with current observational limits. Finally, we study the amount of gravitational waves emitted during a glitch. Simple expressions are obtained for the corresponding characteristic amplitudes and frequencies. The detectability of glitches through gravitational wave observatories is briefly discussed.

**Key words:** gravitational waves – methods: numerical – stars: neutron – pulsars: general – pulsars: individual: PSR B0833–45.

## 1 INTRODUCTION

Pulsars are very compact stars rotating rapidly with exceptionally stable periods spanning from ~1.4 milliseconds to a few seconds. Nevertheless, some pulsars exhibit sudden increases in their observed angular velocity $\Omega$, with relative amplitude $\Delta\Omega/\Omega$ ranging between ~$10^{-11}$ and ~$10^{-5}$ (Wong, Backer & Lyne 2001; Espinoza et al. 2011). These spin-up events, known as *glitches*, are usually followed by a slow relaxation on time-scales up to months or years and are sometimes accompanied by abrupt changes in the pulsar spin-down rate, $\Delta\Omega/\Omega \sim 10^{-4} - 10^{-2}$ (we use a dot to denote time derivative). Presently, 472 glitches have been detected in 165 pulsars<fnref n="1" /> (Espinoza et al. 2011), with angular velocities ranging from 0.09 to 327 Hz (see e.g. the ATNF Pulsar Database;<fnref n="2" /> Manchester

et al. 2005). At least two distinct glitching behaviours have been identified (Espinoza et al. 2011; Yu et al. 2013): (i) quasi-periodic giant glitches with a very narrow spread in size around $\Delta\Omega/\Omega \sim 10^{-6}$, and (ii) smaller glitches of various sizes at random intervals of time. The most emblematic pulsar of the first kind is Vela (PSR B0833–45) with a rotation frequency $f = \Omega/(2\pi) \simeq 11.19$ Hz (Dodson, Lewis & McCulloch 2007) corresponding to a period $P = 1/f \simeq 89$ ms. Since its discovery in 1969, 19 glitches have been detected so far every ~2–3 years (Espinoza et al. 2011). The second type of glitching pulsars is exemplified by the Crab (PSR B0531+21) with a rotation frequency $f \simeq 29.95$ Hz ($P \simeq 33$ ms).

Since the first detections of glitches (Radhakrishnan & Manchester 1969; Reichley & Downs 1969), different mechanisms have been proposed to explain these events (see e.g. the review by Haskell & Melatos 2015). A glitch is nowadays commonly thought as the manifestation of an internal process, except possibly for highly magnetized neutron stars for which some evidence of magnetospheric activity has been found (e.g. Archibald et al. 2013; Keith, Shannon & Johnston 2013; Antonopoulou et al. 2015). The interior

<footnote n="1">¹Laboratoire Univers et Théories, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="2">²Institut d'Astronomie et d'Astrophysique, CP-226, Université Libre de Bruxelles, B-1050 Brussels, Belgium</footnote>

<footnote>* E-mail: aurelien.sourie@obspm.fr</footnote>

<footnote n="1">¹ http://www.jb.man.ac.uk/pulsar/glitches.html</footnote>

<footnote n="2">² http://www.atnf.csiro.au/research/pulsar/psrcat.</footnote>
of neutron stars can thus be probed using observations of pulsar glitches.

Glitches were first suggested to arise from crustquakes (Ruderman 1969; Baym & Pines 1971). Following this idea, the presence of a solid crust (which crystallized when the star was young and rapidly rotating) prevents readjustments of the stellar shape, as the star spins down due to electromagnetic emission. Crustal stresses thus build up, until the crust cracks and the star suddenly adopts a more spherical shape. The resulting reduction of the moment of inertia leads the pulsar to spin up, assuming conservation of angular momentum. This scenario can account for small glitches, such as those exhibited by the Crab pulsar. However, as pointed out by Ruderman (1969), this mechanism fails to predict the occurrence frequency of giant glitches, as observed in the Vela pulsar.

Giant glitches are generally thought to be the manifestation of superfluid matter inside neutron stars, as suggested by the very long time-scales observed during post-glitch relaxations (Baym, Pethick & Pines 1969a). From theoretical calculations, the interior of a neutron star is expected to contain an isotropic neutron superfluid in the inner crust, an anisotropic neutron superfluid in the outer core, and possibly other superfluid species in the inner core (see e.g. Page et al. 2013). In a seminal work, Anderson & Itoh (1975) proposed that glitches themselves could be triggered by the sudden unpinning of neutron superfluid vortices. The idea is the following. It is well known from laboratory experiments (see e.g. Yarmchuk, Gordon & Packard 1979; Abo-Shaeer et al. 2001; Zwierlein et al. 2005) that a superfluid can only rotate by forming an array of quantized vortices, each carrying a quantum $ħ$ of angular momentum, where $ħ$ is the Planck–Dirac constant. The neutron superfluid present in the core and the inner crust of neutron stars is thus expected to be threaded by a huge number of vortex lines, with a mean surface density given by

$$n_v (\text{cm}^{-2}) = \frac{4m_n \Omega_n}{h} \approx \frac{10^7}{P(\text{ms})}, \quad (1)$$

where $h = 2\pi\hbar$ is the Planck constant, $m_n$ is the neutron rest mass and the coarse-grained averaged angular velocity $\Omega_n$ of the neutron superfluid is approximated by that of the star (Ginzburg & Kirzhnits 1965). The neutron superfluid is supposed to be weakly coupled to the rest of the star by so-called mutual friction forces arising from the dissipative forces acting on individual vortices (Alpar, Langer & Sauls 1984c). Due to the spin down of the star induced by the electromagnetic torque, vortices tend to move away from the rotation axis. The key assumption of vortex-mediated glitch theories is that vortices can pin to nuclear clusters in the inner crust (Anderson & Itoh 1975) and/or to quantized magnetic flux tubes in the core if protons form a type II superconductor (Baym et al. 1969a; Sauls 1989; Ruderman, Zhu & Chen 1998). In such case, the neutron superfluid is decoupled from the rest of the star, and can rotate more rapidly, as schematically illustrated in Fig. 1. The lag $\delta\Omega = \Omega_n - \Omega > 0$ induces a Magnus force on the vortices. The larger the lag, the stronger the force. For some critical value $\delta\Omega_0$ of the lag, vortices will suddenly unpin, the superfluid will spin down and, by conservation of angular momentum, the rest of the star will spin up leading to the observed glitch (see Fig. 1). During the subsequent relaxation, these vortices are thought to progressively repin and a lag can grow anew. This vortex-mediated scenario is supported by laboratory experiments (Tsakadze & Tsakadze 1980) and the ability of the vortex creep model to reproduce the post-glitch relaxations in different pulsars (Alpar et al. 1984a,b, 1993, 1996; Güercinoğlu & Alpar 2014).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic time evolution of the observed angular velocity $\Omega$ of the pulsar and of the angular velocity $\Omega_n$ of the interior neutron superfluid over two glitch events.</figcaption></figure>

It is noteworthy to mention that the two mechanisms described above are not necessarily independent. Indeed, starquakes can be induced by the presence of superfluids in neutron-star interiors, whether superfluid vortices are pinned (Ruderman 1991) or not (Carter, Langlois & Sedrakian 2000; Chamel & Carter 2006). In turn, sudden motions of vortices can be triggered by quakes (Ruderman 1991; Chau & Cheng 1993; Alpar et al. 1996; Eichler & Shaisultanov 2010).

Although mesoscopic studies of large collections (~$10^2-10^4$) of vortices provide useful insight (e.g. Warszawski & Melatos 2011, 2013), simulating pulsar glitches requires to follow the dynamics of all the superfluid vortices contained in the star. Given their huge number, ~$10^{17}$ for Vela, the overall transfer of angular momentum between the neutron superfluid and the rest of the star can be studied using a smooth-averaged hydrodynamic approach, still involving microscopic parameters determined by the local dynamics of individual vortices (see e.g. Bulgac, Forbes & Sharma 2013, and references therein). Whereas the general relativistic framework for describing starquakes was developed a long time ago (Carter & Quintana 1975), the general relativistic formulation of the vortex-mediated glitch model is more recent (Langlois, Sedrakian & Carter 1998). As a matter of fact, most global numerical simulations of pulsar glitches have been performed within the Newtonian framework (e.g. Larson & Link 2002; Peralta et al. 2006; Sidery, Passamonti & Andersson 2010; Haskell, Pizzochero & Sidery 2012). Recently, Seveso, Pizzochero & Haskell (2012) and Antonelli & Pizzochero (2017) have developed a non-relativistic hydrodynamic model for describing the different stages of the glitch phenomenon based on the static structure of the neutron star computed in general relativity. However, general relativity could also play an important role for the global dynamics of glitches. Furthermore, general relativity is essential to determine the amount of gravitational waves associated with glitch events. Observations of gravitational waves are of particular interest since they could potentially provide additional information on the glitch phenomenon (see e.g. Stopnitzky & Profumo 2014; Haskell & Melatos 2015 and references therein).

In this paper, we present global numerical simulations of vortex-mediated pulsar glitches. We focus on the spin-up stage regardless of the glitch triggering mechanism. On the other hand, we study the glitch dynamics in full general relativity. We also derive the associated gravitational wave characteristic amplitudes and frequencies using the standard quadrupole formula. The paper is organized as follows. We start by presenting, in Section 2, the different assumptions on which our model is based. In Section 3, we introduce the evolution equations governing the transfer of angular momentum in the interior of a pulsar during a glitch. Results of stationary
rotating configurations are discussed in Section 4. In Section 5, we detail the numerical procedure underlying our simulations. Results for the glitch rise time are presented and discussed. We study the emission of gravitational waves in Section 6. Finally, we conclude in Section 7.

## 2 MODEL ASSUMPTIONS

### 2.1 Quasi-stationary approach

The glitch phenomenon can be decomposed into distinct stages (i.e. the pre-glitch evolution, the spin up, and the post-glitch relaxation), which can be modelled separately in view of the different associated time-scales suggesting different physical mechanisms. Focusing on the sudden spin up of the pulsar after the catastrophic unpinning of vortices, stellar dynamics are essentially governed by the mutual friction force between the superfluid and the rest of the star. This force acts on a characteristic time-scale corresponding to the glitch rise time $\tau_r$, which has not been fully observationally resolved yet. The most stringent observational constraint on $\tau_r$ comes from the 2000 and 2004 Vela glitch timing data: $\tau_r < 30–40$ s (Dodson, McCulloch & Lewis 2002; Dodson et al. 2007).

It is interesting to compare $\tau_r$ with the typical time $\tau_h$ for the star to go back to hydrodynamic equilibrium once being driven out of it by a change in its rotation rate. Sometimes referred to as hydrodynamic time-scale, this time is roughly given by (Shapiro & Teukolsky 1983)

$$ \tau_h \sim (G\bar{\rho})^{-\frac{1}{2}} \approx 0.1 \left(\frac{R}{12 \text{ km}}\right)^{\frac{3}{2}} \left(\frac{M}{1.4 \text{ M}_\odot}\right)^{-\frac{1}{2}} \text{ms,} \quad (2) $$

where G denotes the gravitational constant, M is the neutron-star mass, R is the stellar radius, and $\bar{\rho} \sim 3M/(4\pi R^3)$ is the average density of the star. The hydrodynamic time $\tau_h$ represents the time for a sound wave with speed $c_s$ to propagate throughout a star of radius R, i.e. $\tau_h \sim R/c_s$ (Epstein 1988).

In the following, we shall assume that $\tau_r \gg \tau_h$, as suggested by previous studies (Haskell et al. 2012), so that the dynamical evolution of the pulsar can be reasonably well described by a sequence of quasi-stationary equilibrium configurations.

### 2.2 Two-component model

Due to the magnetic field, the electrically charged particles inside neutron stars are strongly coupled and essentially corotate with the crust and the magnetosphere at the observed angular velocity $\Omega$ (Glampedakis, Andersson & Samuelsson 2011). We do not account for other effects of the magnetic field on the global dynamics of the star, which could be important for the most strongly magnetized neutron stars (Bocquet et al. 1995; Chatterjee et al. 2015), but can be safely ignored for the ordinary pulsars considered here.

The simplest model of pulsars thus consists of at least two distinct dynamical components (Baym et al. 1969b): (i) a plasma of charged particles (electrons, nuclei in the crust, and protons in the core), and (ii) a neutron superfluid extending in the whole core. Because of (non-dissipative) mutual neutron-proton entrainment effects according to which the momenta are misaligned with the corresponding velocities (Andreev & Bashkin 1976), neutron superfluid vortices in the core of a neutron star carry a fractional magnetic quantum flux (Sedrakian & Shakhabasian 1980). Electrons scattering off the magnetic field of the vortex lines was shown to induce a strong coupling between the core superfluid and the crust

(Alpar et al. 1984c). For this reason, only the neutron superfluid permeating the inner crust of the star has been generally thought to be responsible for giant glitches. This scenario was also supported by the analysis of the glitch data, which suggested that the superfluid represents only a few per cent of the angular momentum reservoir of the star (Alpar et al. 1993; Datta & Alpar 1993; Link, Epstein & Lattimer 1999). On the other hand, this model has been recently challenged (Chamel & Carter 2006; Andersson et al. 2012; Chamel 2013; Delsate et al. 2016) by the realization that despite the absence of viscous drag the crust can still resist the flow of the neutron superfluid due to Bragg scattering (Carter, Chamel & Haensel 2005; Chamel 2004, 2005, 2012). It has been argued that crustal entrainment could be much weaker assuming that the superfluid coherence length is much smaller than the size of clusters (Martin & Urban 2016). However, as recognized by these authors, this condition is generally not satisfied. Even if crustal entrainment is ignored, the analyses of the 2007 glitch detected in PSR J1119–6127, as well as of the 2010 glitch in PSR B2334+61, indicate that the crust is not enough (Yuan et al. 2010; Alpar 2011; Akbal et al. 2015). These recent studies suggest that the core superfluid plays a more important role than previously thought. In particular, the core superfluid could be decoupled from the rest of the star due to the pinning of vortices to quantized magnetic flux tubes assuming protons form a type II superconductor (Gügercinoğlu & Alpar 2014).

In this work, we thus focus on the dynamics of the superfluid neutron-star core within a two-fluid model: a neutron superfluid coupled to the ‘normal’ fluid made of protons and electrons (simply referred to as ‘protons’ in the following). Quantities related to the two fluids will be labelled by indices ‘n’ and ‘p’, respectively. Note that, since in our model we do not consider the stellar crust, we do not account for any crust-core coupling mechanisms, such as Ekman pumping. These couplings could still have a strong impact on the glitch dynamics, especially during the post-glitch relaxation (van Eysden & Melatos 2010; Haskell & Melatos 2015).

### 2.3 Space-time symmetries

Our glitch simulations are based on the general relativistic equilibrium configurations of rotating superfluid neutron stars computed by Sourie, Oertel & Novak (2016). In this section, we recall the main assumptions on the space-time symmetries and the metric.

The star is supposed to be axisymmetric and stationarily rotating. Neglecting the very small non-circular motion of the neutron superfluid due to the radial displacement of the vortices during the glitch (Langlois et al. 1998), the two fluids are further assumed to rotate around a common axis with possibly different rotation rates. The space-time is thus stationary, axisymmetric, circular and asymptotically flat. While the angular velocity $\Omega_p = \Omega$ of the charged components can be reasonably assumed to be uniform, the angular velocity $\Omega_n$ of the neutron superfluid may vary throughout the star. We circumvent this difficulty by considering that both fluids are rigidly rotating as in the model of Sidery et al. (2010) in Newtonian gravity. In this case, $\Omega_n$ and $\Omega_p$ are to be understood as the angular velocities averaged over the whole star. More details on general relativistic equilibrium configurations of rotating superfluid neutron stars can be found in Prix, Novak & Comer (2005) and Sourie et al. (2016).

### 2.4 Chemical composition and equation of state

The dominant electroweak processes governing the composition of a neutron star are the direct (DU) and modified (MU) Urca
beta processes (Yakovlev et al. 2001). Within the assumption of rigid-body rotation with a non-vanishing lag $\delta\Omega$, beta equilibrium can only possibly be achieved on the rotational axis of the star (Andersson & Comer 2001). Assuming corotation ($\delta\Omega = 0$) and ignoring superfluidity, the relaxation times towards beta equilibrium are roughly given by (Yakovlev et al. 2001)

$$ \tau_{\beta}^{(\text{DU})} \approx 20 \left( \frac{T}{10^9 \text{ K}} \right)^{-4} \text{s} \qquad (3) $$

$$ \tau_{\beta}^{(\text{MU})} \approx \left( \frac{T}{10^9 \text{ K}} \right)^{-6} \text{months,} \qquad (4) $$

where T represents the interior temperature of the star. For glitching pulsars, whose characteristic 'ages' $\tau_{sd} = \Omega/(2|\dot{\Omega}|) > 10^3$ yr (Espinoza et al. 2011) correspond to temperatures below $\sim 10^9$ K (Gnedin, Yakovlev & Potekhin 2001), equations (3) and (4) lead to time-scales of the order of a few tens of seconds and a month for the DU and MU processes, respectively. Differential rotation could in principle change the chemical equilibrium (Langlois et al. 1998); however the resulting effects are presumably negligible in view of the very small lags $\delta\Omega \ll \Omega$. On the contrary, superfluidity can strongly reduce the rates of beta processes (Villain & Haensel 2005), making these time-scales even longer. Therefore, beta equilibrium may not be achieved during the spin up.

The DU and MU rates being poorly known, we assume that $\tau_{\beta}$ still remains much shorter than the interglitch time so that the star is in beta equilibrium on the rotational axis at the beginning of a glitch. During the glitch rise, in order to estimate the error induced by our lack of knowledge, we restrict to the two non-dissipative limiting cases:

(i) $\tau_{\beta} \gg \tau_{r}$, i.e. no reaction takes place during the glitch. The baryon masses of the two fluids $M_{n}^{B}$ and $M_{p}^{B}$ thus remain separately constant, which leads stellar matter on the rotational axis to be (slightly) out of beta equilibrium.

(ii) $\tau_{\beta} \ll \tau_{r}$, i.e. the stellar matter on the rotational axis goes back instantaneously to beta equilibrium so that the chemical potentials satisfy $\mu^{n} = \mu^{p}$, where $\mu^{n}$ ($\mu^{p}$) represents the neutron (proton) chemical potential. In our simulations, it is sufficient to impose this condition at the centre of the star, i.e. $\mu_{c}^{n} = \mu_{c}^{p}$, as discussed by Prix et al. (2005). Only the total baryon mass $M^{B} = M_{n}^{B} + M_{p}^{B}$ is conserved in this case.

We adopt the same equations of state (EoSs) as in Sourie et al. (2016). These EoSs referred to as DDH and DDH$\delta$ were calculated from density-dependent relativistic mean-field models, including $\sigma$, $\omega$, $\rho$ mesons for the former and in addition the $\delta$ meson for the latter (Avancini et al. 2009; Typel & Wolter 1999). They were adapted to a system of two fluids at zero temperature coupled by entrainment for arbitrary compositions.

# 3 EVOLUTION EQUATIONS

## 3.1 Angular momentum transfer

Let $J_{n}$ and $J_{p}$ be the neutron superfluid and proton fluid angular momenta, respectively (see Langlois et al. 1998; Sourie et al. 2016 for definitions and expressions). Neglecting any external torque, the dynamics of these two fluids during the spin up is simply governed by the following equations

$$ \begin{cases} \dot{J}_{n} = +\Gamma_{\text{mf}}, \\ \dot{J}_{p} = -\Gamma_{\text{mf}}, \end{cases} \qquad (5) $$

where $\Gamma_{\text{mf}}$ stands for the mutual friction torque and overdot for time derivative.

Langlois et al. (1998) derived a covariant expression for the (relativistic) mutual friction torque $\Gamma_{\text{mf}}$, considering straight vortices parallel to the rotation axis and arranged on a regular array. The dynamical evolution of the superfluid might potentially lead to the formation of a vortex tangle (Peralta et al. 2006; Andersson, Sidery & Comer 2007). However, the onset of superfluid turbulence remains highly speculative, and therefore we shall not consider this possibility here. In Langlois et al. (1998), the motion of unpinned vortex lines was assumed to be simply determined from the interplay between a Magnus force due to the neutron fluid and a dissipative drag force caused by the proton fluid.

Our numerical approach is based on the 3+1 formalism, in which the space-time is foliated by a family $(\Sigma_{r})_{r \in \mathbb{R}}$ of space-like hypersurfaces (see e.g. Gourgoulhon 2012). Neglecting the small contribution of the non-circular motion of the vortices, the 3+1 expression of the mutual friction torque derived by Langlois et al. (1998) reads

$$ \Gamma_{\text{mf}} = \int_{\Sigma_{r}} \mathcal{B} \Gamma_{n} n_{n} \omega_{n} \chi_{\perp}^{2} d^{3}\Sigma \times (\Omega_{p} - \Omega_{n}), \qquad (6) $$

in the absence of any dissipation related to chemical reactions [see equations (72) and (89) of Langlois et al. (1998)]. In this expression, $d^3\Sigma$ denotes the volume element on the hypersurface $\Sigma_{r}$ and $n_{n}$ is the density of the neutron fluid in its rest frame. The term $\Gamma_{n}$ stands for the Lorentz factor of the neutrons with respect to the so-called Eulerian observer or zero-angular-momentum observer (ZAMO), whose four-velocity corresponds to the unit future-oriented (time-like) vector normal to $\Sigma_{r}$ (e.g. Gourgoulhon 2012). The macroscopic neutron vorticity $\omega_{n}$ reads (we use Greek letters for space-time indices)

$$ \bar{\omega}_{n} = \sqrt{\frac{\bar{\omega}_{\mu\nu} \bar{\omega}^{\mu\nu}}{2}}, \qquad (7) $$

where the vorticity two-form $\bar{\omega}_{\mu\nu}$ is defined by

$$ \bar{\omega}_{\mu\nu} = \nabla_{\mu} p_{\nu}^{n} - \nabla_{\nu} p_{\mu}^{n}, \qquad (8) $$

$p_{\mu}^{n}$ denoting the conjugate superfluid momentum. We note here that, on length scales smaller than the intervortex separation $d_{v}$, typically of the order of $d_{v} \sim n_{v}^{-1/2} \simeq 10^{-3}$ cm (see equation 1), $\bar{\omega}_{\mu\nu}$ strictly vanishes because $p_{\mu}^{n}$ should be locally proportional to the gradient of a quantum scalar phase. Nevertheless, on the large scales we are interested in here, the neutron vorticity two-form is non-vanishing, along with its corresponding scalar amplitude $\bar{\omega}_{n}$. Moreover, the geometric quantity $\chi_{\perp}^{2}$ is defined by

$$ \chi_{\perp}^{2} = \chi_{\perp}^{\alpha} \chi_{\perp}^{\alpha}, \quad \chi_{\perp}^{\alpha} = \perp_{\beta}^{\alpha} \chi^{\beta}, \qquad (9) $$

where $\chi^{\alpha}$ is the Killing vector associated with axisymmetry and $\perp_{\alpha\beta}$ is the projection tensor orthogonal to the two-dimensional string-type world sheets representing the vortex cores [see Langlois et al. (1998) for details].

Finally, the mutual friction parameter $\mathcal{B}$ characterizes the efficiency of the angular momentum transfer through mutual friction. This parameter is given by (see e.g. Langlois et al. 1998; Carter 2001)

$$ \mathcal{B} = \frac{\mathcal{R}}{1 + \mathcal{R}^2} \qquad (10) $$

as a function of the positive dimensionless drag-to-lift ratio $\mathcal{R}$. Since the dissipative processes contributing to mutual friction are not the same in different stellar regions (Alpar et al. 1984c; Jones 1990, 1992; Epstein & Baym 1992; Sedrakian & Sedrakian 1995; Haskell,
Glampedakis & Andersson 2014), $\mathcal{R}$ is expected to vary throughout the star. However, the values of this coefficient remain very uncertain. Microscopic estimates differ by many orders of magnitude. Given the current lack of knowledge on the microscopic origin of the mutual friction force, and since we are interested in global models of neutron stars, we introduce the averaged coefficient

$$ \bar{\mathcal{B}} = \frac{\int_{\Sigma_i} \mathcal{B} \Gamma_n n_n \bar{\omega}_n \chi_{\perp}^2 d^3\Sigma}{\int_{\Sigma_i} \Gamma_n n_n \bar{\omega}_n \chi_{\perp}^2 d^3\Sigma}, \quad (11) $$

which we consider as a free input parameter in our numerical simulations. Although $\mathcal{R}$ is likely to vary in time during the glitch event (due to changes of the vortex velocity or the repinning of some vortices, for instance), $\bar{\mathcal{B}}$ is assumed to be time-independent for simplicity. The mutual friction torque (6) thus becomes

$$ \Gamma_{\mathrm{mf}} = -\bar{\mathcal{B}} \int_{\Sigma_i} \Gamma_n n_n \bar{\omega}_n \chi_{\perp}^2 d^3\Sigma \times \delta\Omega, \quad (12) $$

where $\delta\Omega = \Omega_n - \Omega_p$ is the lag between the fluids. As shown in Appendix A1, the Newtonian limit of equation (12) is in perfect agreement with the expression given by Sidery et al. (2010) (see their equation 58).

To describe any transfer of angular momentum, it is convenient to introduce the partial moments of inertia

$$ I_{XX} = \left( \frac{\partial J_X}{\partial \Omega_X} \right)_{\Omega_Y} \quad \text{and} \quad I_{XY} = \left( \frac{\partial J_Y}{\partial \Omega_Y} \right)_{\Omega_X}, \quad (13) $$

where the two different capital letters X and Y refer to protons or neutrons. Depending on the assumption on the chemical composition (see Section 2.4), these derivatives are taken either for fixed partial baryon masses $M_n^B$ and $M_p^B$ (case i) or for a constant total baryon mass $M^B$ with identical chemical potentials at the centre of the star (case ii). It is possible to show that $I_{np} = I_{pn}$ [see equation (3.10) from Carter (1975)]. Furthermore, the partial moments of inertia should obey the following conditions:

$$ I_{nn} > 0, \quad I_{pp} > 0 \quad \text{and} \quad I_{nn} I_{pp} > I_{np}^2; \quad (14) $$

see Appendix B for more details. We also define the neutron and proton moments of inertia, $\hat{l}_n$ and $\hat{l}_p$, as

$$ \hat{l}_n = I_{nn} + I_{np} \quad \text{and} \quad \hat{l}_p = I_{pp} + I_{np}, \quad (15) $$

and the total moment of inertia $\hat{l}$ by $\hat{l} = \hat{l}_n + \hat{l}_p$. We note here that these definitions are more general than the moments of inertia discussed in section II-D of Sourie et al. (2016) in the limiting case of corotating fluids. However, the two definitions coincide in the slow-rotation approximation.

Given the previous definitions, the time derivatives of the angular momenta can be expressed as

$$ \left\{ \begin{array}{l} \dot{J}_n = I_{nn} \dot{\Omega}_n + I_{np} \dot{\Omega}_p, \\ \dot{J}_p = I_{np} \dot{\Omega}_n + I_{pp} \dot{\Omega}_p. \end{array} \right. \quad (16) $$

Using expression (12) of the mutual friction torque, the angular momentum transfer (5) reads

$$ \left\{ \begin{array}{l} \dot{\Omega}_n = - \frac{\dot{l}_p}{I_{nn} l_{pp} - I_{np}^2} \times \bar{\mathcal{B}} \kappa \delta\Omega, \\ \dot{\Omega}_p = + \frac{\dot{l}_n}{I_{nn} l_{pp} - I_{np}^2} \times \bar{\mathcal{B}} \kappa \delta\Omega, \end{array} \right. \quad (17) $$

where we have introduced the quantity $\kappa$ defined by

$$ \kappa = \int_{\Sigma_i} \Gamma_n n_n \bar{\omega}_n \chi_{\perp}^2 d^3\Sigma. \quad (18) $$

It should be noted that, in view of the conditions (14), the denominator appearing in equation (17) never vanishes. The time evolution of the lag $\delta\Omega = \Omega_n - \Omega_p$ is thus governed by the simple equation

$$ \frac{\delta \dot{\Omega}}{\delta \Omega} = - \frac{\hat{l}}{I_{nn} I_{pp} - I_{np}^2} \times \bar{\mathcal{B}} \kappa. \quad (19) $$

## 3.2 Analytical estimate for the spin-up time-scale

We now focus on deriving an approximate analytical formula for the spin-up time-scale. Recalling that $\Delta\Omega/\Omega \ll 1$, where $\Delta\Omega$ represents the variation in the pulsar rotation rate during the spin up, it is a reasonable approximation to neglect the change in the different quantities appearing in the right-hand side of equation (19). Starting from an initial lag $\delta\Omega_0$ at the beginning of the spin up (see Section 5.1.2), the lag therefore approximately evolves as

$$ \delta\Omega(t) \approx \delta\Omega_0 \times \exp(-t/\tau_r), \quad (20) $$

where we have introduced the characteristic time-scale

$$ \tau_r = \frac{\hat{l}_n \hat{l}_p - \hat{l} I_{np}}{\hat{l} \bar{\mathcal{B}} \kappa}. \quad (21) $$

The time evolution of the two angular velocities is given by

$$ \Omega_n(t) \approx \Omega_n^0 + \frac{\dot{l}_p}{\hat{l}} \times \delta\Omega_0 \times (\exp(-t/\tau_r) - 1), \quad (22) $$

$$ \Omega_p(t) \approx \Omega_p^0 - \frac{\dot{l}_p}{\hat{l}} \times \delta\Omega_0 \times (\exp(-t/\tau_r) - 1), $$

where $\Omega_n^0$ and $\Omega_p^0$ are the fluid rotation rates when the glitch is triggered.

By analogy with the Newtonian limit (A10), we introduce the quantity $\zeta$ through the relation

$$ \kappa = 2\zeta \hat{l}_n \Omega_n. \quad (23) $$

Still, it should be remarked that general relativistic corrections are not only included in $\zeta$ but are also partially contained in $\hat{l}_n$. Using (23), the general relativistic rise time (21) now reads

$$ \tau_r = \frac{\hat{l}_p}{\hat{l}} \times \frac{1}{2\zeta \bar{\mathcal{B}} \Omega_n} \times \left( 1 - \frac{I_{np} \hat{l}}{\hat{l}_n \hat{l}_p} \right). \quad (24) $$

Considering slowly rotating stars, for which $\Omega_n$, $\Omega_p \ll \Omega_K$ where the Keplerian limit $\Omega_K$ is the maximum angular velocity above which mass-shedding occurs at the equator, the Newtonian limit of equation (24) is found to be in perfect agreement with results from Sidery et al. (2010) (see Appendix A2). Note that this approximation is quite reasonable, given the rotation frequencies of observed glitching pulsars and the estimated values of the Keplerian frequency $\Omega_K/(2\pi)$, which is of the order of ~900 to 1800 Hz for neutron stars with masses larger than 1.4 M$_\odot$ (e.g. Haensel et al. 2009; Fantina et al. 2013; Haensel et al. 2016).

The analytical expression (24) of the spin-up time-scale calls for several remarks. First, since the radial velocity of the vortex lines increases with the mutual friction parameter [see equation (36) of Carter (2001) in the Newtonian case and equation (85) of Langlois et al. (1998) in the relativistic framework], the larger $\bar{\mathcal{B}}$ is, the faster is the transfer of angular momentum, consistently with (24). As pointed out by Carter (2001), the parameter $\bar{\mathcal{B}}$ cannot take any arbitrary value: according to (10), we have $\bar{\mathcal{B}} \le 1/2$. This implies the following lower bound for the spin-up time-scale:

$$ \tau_r \ge \frac{\hat{l}_p}{\hat{l}} \times \frac{1}{\zeta \Omega_n} \times \left( 1 - \frac{I_{np} \hat{l}}{\hat{l}_n \hat{l}_p} \right), \quad (25) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Superfluid vorticity in a neutron star plotted with respect to the radial coordinate $r$, in the equatorial ($\theta = \pi/2$) and polar ($\theta = 0$) planes, for $f_n = f_p = 11.19$ Hz, assuming beta equilibrium. The neutron vorticity $\sigma_n$ is normalized to its expected value $2m_n\Omega_n$ in the Newtonian limit, if corotation is enforced (black dashed lines). Blue vertical lines represent the radii at which $n_n = 0$, for both polar and equatorial planes given that the star is approximately spherical at such a low angular velocity. Left: neutron vorticity profiles for a star with a gravitational mass $M_G = 1.4 M_{\odot}$, using DDH (blue lines) and DDH$\delta$ (orange lines) EoSs. Right: influence of the gravitational mass of the star on the vorticity profiles, using the DDH EoS. Results in the equatorial (polar) plane are plotted in solid (dashed) lines in both panels. See the online version for colours.</figcaption></figure>

which is only reached for $R = 1$. Furthermore, in the slow-rotation approximation, $\tau_r$ is also found to be inversely proportional to $\Omega_n$, which can be interpreted from the fact that the mutual friction torque is proportional to the surface density of vortex lines [see equation (65) of Andersson, Sidery & Comer (2006) in the Newtonian context] through the superfluid vorticity $\sigma_n$ (7), which in turn is roughly proportional to $\Omega_n$ (see equation 1). The roles of $\zeta$ and $I_{np}$ will be studied in details in the following section.

# 4 STATIONARY ROTATING CONFIGURATIONS

In this section, we present numerical results concerning the superfluid vorticity and the couplings between the fluids, which are both playing an important role in the angular momentum transfer during a glitch event. These quantities can be directly obtained from the equilibrium configurations computed by Sourie et al. (2016).

## 4.1 Superfluid vorticity

In Fig. 2, the superfluid vorticity $\sigma_n$ (7) is plotted as a function of the radial coordinate, for a star spinning at Vela's rotation frequency, i.e. $f_n = f_p = 11.19$ Hz, and assuming beta equilibrium. The vorticity is normalized to its Newtonian limit $\sigma_n^{\text{newt}} = 2m_n\Omega_n$. In the left panel, the results from both EoSs are compared, whereas in the right panel vorticity profiles are represented for different gravitational masses, using the DDH EoS.

The deviation from the Newtonian value can be simply interpreted in terms of the compactness parameter $\Xi$ of the star, defined as the dimensionless ratio of the gravitational mass $M_G$ of the star to its circumferential radius in the equatorial plane $R_{c,eq}$ (see Gourgoulhon 2010 for definitions), i.e.

$$ \Xi = \frac{GM_G}{R_{c,eq}c^2}. \qquad (26) $$

For a $1.4 M_{\odot}$ neutron star spinning at 11.19 Hz, the DDH$\delta$ EoS predicts a larger radius than the DDH one, which leads to $\Xi_{\text{DDH}} \gtrsim \Xi_{\text{DDH}\delta}$. The deviation from the non-relativistic case is therefore slightly stronger for DDH, as can be seen in the left panel

of Fig. 2. Similarly, since the compactness parameter increases with the mass of the star, the deviation from the expected value in the Newtonian limit is more important for more massive stars, and can reach a maximum of ~30 per cent (see the right panel of Fig. 2). It is interesting to note that, since the compactness parameter is smaller for higher rotation rates at fixed gravitational mass, the general relativistic correction is found to be less important for stars spinning more rapidly.

Moreover, the quantity $\zeta$ involved in the spin-up time-scale (24) can be determined from equations (18) and (23) using stationary configurations, by computing the superfluid vorticity profile and the appropriate moments of inertia. To reach high accuracy, the latter are calculated from equation (13) using a fourth-order finite difference method, either at given fluid baryon masses (case i) or for a fixed total baryon mass with chemical equilibrium at the centre (case ii), depending on the assumption on the composition (see Section 2.4).

For a given rotation frequency, neglecting the lag and assuming beta equilibrium, $\zeta$ is found to decrease with increasing mass and tends towards 1 for small compactness parameters (see Appendix A2). In particular, for a star spinning at 11.19 Hz, $\zeta$ thus changes from 0.869 (0.880) for $M_G = 1.4 M_{\odot}$ to 0.763 (0.792) for a $2 M_{\odot}$ neutron star with the DDH($\delta$) EoS. On the other hand, at fixed gravitational mass, $\zeta$ is approximately independent of the angular velocity in the slow-rotation approximation ($\Omega_n, \Omega_p \ll \Omega_K$). Since the moments of inertia are also nearly constant for low rotation rates, $\tau_r$ turns out to be inversely proportional to the angular velocity (see equation 24). However, for rotation frequencies substantially higher than that of Vela, the moments of inertia are found to increase more rapidly than $\kappa$, meaning that $\zeta$ decreases when the angular velocity gets higher. For instance, the DDH EoS leads to $\zeta = 0.834$ for a $1.4 M_{\odot}$ neutron star spinning at 327 Hz, which is the rotation frequency of the fastest glitching pulsar observed so far (see the ATNF Pulsar Database;<fnref n="3" /> Manchester et al. 2005).

<footnote n="3">³ http://www.atnf.csiro.au/research/pulsar/psrcat.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left: entrainment parameter $\tilde{\epsilon}_p$ and Lense-Thirring coefficients $\epsilon_{np}^{LT}$ and $\epsilon_{nn}^{LT}$ (see text and Appendix C for definitions) as functions of the gravitational mass $M_G$ of the star, for $f_n = f_p = 11.19$ Hz. Only results obtained from the DDH EoS are presented. Right: total coupling coefficient $\tilde{\epsilon}_p$ and entrainment parameter $\tilde{\epsilon}_p$ with respect to the gravitational mass, for a star spinning at 11.19 Hz, assuming corotation and beta equilibrium. Results obtained with the DDH($\delta$) EoS are shown in solid (dashed) lines.</figcaption></figure>

## 4.2 Entrainment and frame-dragging effects

The cross-moment of inertia $I_{np}$ introduced in equation (13) contains all the possible couplings between neutrons and protons. A first coupling is due to entrainment, which in the core of neutron stars comes from the strong interactions between nucleons. This effect has already been discussed in detail in section III of Sourie et al. (2016) for both DDH and DDH$\delta$ EoSs. In Newtonian gravity, entrainment is the main fluid coupling at low angular velocities (see Appendix A3). Nevertheless, as already mentioned in Sourie et al. (2016), a new coupling arises in the general relativistic context from the so-called Lense-Thirring or frame-dragging effect (see Carter 1975).

In what follows, we characterize the total coupling between the fluids through the following quantities

$$ \hat{\epsilon}_n = \frac{I_{np}}{\hat{l}_n} \quad \text{and} \quad \hat{\epsilon}_p = \frac{I_{np}}{\hat{l}_p}. \qquad (27) $$

From the stability conditions (14), these parameters are not arbitrary but must satisfy the following inequality:

$$ \left(\frac{1}{\hat{\epsilon}_p} - 1\right) \left(\frac{1}{\hat{\epsilon}_n} - 1\right) > 1. \qquad (28) $$

At first order in the lag $\delta\Omega = \Omega_n - \Omega_p$ and in the slow-rotation approximation, these coupling coefficients can be written as

$$ \hat{\epsilon}_X = \frac{\tilde{\epsilon}_X - \epsilon_{YX}^{LT}}{1 - \epsilon_{YX}^{LT} - \epsilon_{XX}^{LT}}, \qquad (29) $$

where X and Y ≠ X are n or p for neutrons or protons, respectively (see Appendix C for details). In this expression, the term $\tilde{\epsilon}_X$ characterizes entrainment effects averaged over the star (C6), whereas $\epsilon_{YX}^{LT}$ and $\epsilon_{XX}^{LT}$ represent, respectively, the frame-dragging effect on fluid X caused by the second fluid and fluid X itself (see equations C7 and C9). Quite remarkably, frame-dragging effects lead to similar fluid couplings as the entrainment arising from neutron-proton interactions (see equation C10). Since the different quantities involved in equation (29) are positive, the Lense-Thirring effect is found to act in an opposite way to entrainment in the core. The reason for the presence of minus signs in front of every quantity relative to frame-dragging effects is the following: a ZAMO will rotate in the same sense as the whole star, leading a zero-angular-velocity

observer to have an angular momentum with an opposite sign to the total angular momentum of the star, even if this observer is static (Carter 1975). As a consequence, in the absence of entrainment, the total coupling coefficient is still expected to be non-vanishing and negative. Although entrainment is likely to be small in the outermost regions of the core of neutron stars (Carter, Chamel & Haensel 2006a; Chamel & Haensel 2006), its overall effect on the whole star is not necessarily negligible and therefore $\hat{\epsilon}_X$ could be positive or negative.

To assess the relative importance of these two effects to the total coupling coefficient $\hat{\epsilon}_p$, the quantities $\tilde{\epsilon}_p$ and $\epsilon_{np}^{LT}$ are plotted in the left panel of Fig. 3 as functions of the gravitational mass $M_G$ of the star for the DDH EoS. Since typically $\delta\Omega \ll \Omega_n$, $\Omega_p$, we consider here that $\Omega_n = \Omega_p = 2\pi \times 11.19$ rad s$^{-1}$. The mean proton entrainment parameter $\tilde{\epsilon}_p$ is increasing with the mass, because higher densities are reached in the star (see Sourie et al. 2016, fig. 2). Since general relativistic effects are the strongest for the most massive stars, $\epsilon_{np}^{LT}$ increases with the stellar mass. It is interesting to note that $\tilde{\epsilon}_p$ and $\epsilon_{np}^{LT}$ are found to be roughly of the same order of magnitude, making the Lense-Thirring contribution to the total coupling be very important. A similar conclusion is reached with the DDH$\delta$ EoS. Since $\epsilon_{nn}^{LT}$ and $\epsilon_{np}^{LT}$ both characterize the contribution of the neutron fluid to frame-dragging effects, these two parameters are very close to each other, as can be seen in the left panel of Fig. 3. For the same reasons, we have $\epsilon_{pn}^{LT} \simeq \epsilon_{pp}^{LT}$. Moreover, we numerically find that $\epsilon_{pp}^{LT} \simeq \hat{l}_p/\hat{l}_n \times \epsilon_{nn}^{LT} \ll \epsilon_{nn}^{LT}$, which means that the relative contribution of the two fluids to the frame-dragging effect is mainly due to their relative proportion in mass, as expected. The typical values obtained for the Lense-Thirring parameters are found to be consistent with the rough estimates given by Carter (1975), within a factor of ~4 for both EoSs.

Using the DDH EoS, the equality of the cross-moments of inertia $I_{np}$ and $I_{pn}$ (see Section 3.1), is numerically verified with a precision better than ~0.1 per cent for $f < 327$ Hz. In the following, we mainly focus on the proton coupling parameter $\hat{\epsilon}_p$ (27). Indeed, the neutron parameter $\hat{\epsilon}_n$ can be simply deduced from equation (27), i.e. $\hat{\epsilon}_n = \hat{l}_p/\hat{l}_n \times \hat{\epsilon}_p \ll \hat{\epsilon}_p$. For the numerical results displayed in Fig. 3, we fixed the total baryon masses (case ii). However, for the low rotation frequency we considered, fixing the individual baryon masses (case i) would have led to essentially the same results.
As displayed in the right panel of Fig. 3 for a star spinning at 11.19 Hz, the total coupling coefficient $\hat{\epsilon}_p$ is found to decrease significantly when the gravitational mass increases. This means that, although entrainment effects become more important, the frame-dragging contribution is increasing even more rapidly (see the left panel of Fig. 3). The discrepancy between entrainment parameters and total coupling coefficients can be clearly seen in the right panel of Fig. 3. For slowly rotating stars (with $f \lesssim 100$ Hz), results obtained from equations (27) and (29) agree with a precision better than ~0.1 per cent for the DDH EoS and ~0.3 per cent for the DDH$\delta$ one. As DDH and DDH$\delta$ predict stars with similar compactness parameters, frame-dragging effects are nearly the same for the two EoSs. Nevertheless, since entrainment effects are much stronger with DDH (see fig. 2 of Sourie et al. 2016), the total coupling is higher for DDH than for DDH$\delta$.

Making use of the coupling parameters (27), the general relativistic spin-up time-scale (24) leads to a similar expression to that obtained within the Newtonian framework (see equation A12), namely

$$ \tau_r = \frac{\hat{I}_p}{\hat{I}} \times \frac{1 - \hat{\epsilon}_p - \hat{\epsilon}_n}{2\zeta \hat{B} \Omega_n} \qquad (30) $$

It should be stressed, however, that general relativistic effects are not only included in the coefficient $\zeta$ but can also change the other parameters (for more details, see Section 5.2.2). Furthermore, using equation (29) and considering that $\epsilon_{np}^{LT} \simeq \epsilon_{nn}^{LT}$ and $\epsilon_{pn}^{LT} \simeq \epsilon_{pp}^{LT}$, the spin-up time-scale is roughly given by

$$ \tau_r \approx \frac{\hat{I}_p}{\hat{I}} \times \frac{1 - \hat{\epsilon}_p - \hat{\epsilon}_n}{1 - \hat{\epsilon}_{pp}^{LT} - \hat{\epsilon}_{nn}^{LT}} \times \frac{1}{2\zeta \hat{B} \Omega_n} \qquad (31) $$

in the slow-rotation approximation. From this expression, we clearly see that the Lense-Thirring effect acts to slow down the angular momentum transfer, whereas entrainment effects in the core tend to make it more efficient. Still, note that entrainment in the inner crust, which comes from Bragg scattering of dripped neutrons by nuclei (Chamel 2005, 2012), is expected to lead to $\hat{\epsilon}_X < 0$. In such a case, both entrainment and frame-dragging effects shall contribute to slow down the glitch event.

So far, we only focused on the slow-rotation approximation, in which the two fluids composing a neutron star are only coupled through entrainment and frame-dragging effects. Although this assumption is totally justified for Vela, it may not remain valid for other glitch pulsars. In particular, for rotation frequencies higher than ~100 Hz, the deformation of each fluid due to rotation leads to additional couplings via gravity, similarly to the Newtonian case (see Appendix A3). Beyond ~100 Hz, the coupling coefficient $\hat{\epsilon}_p$ is found to decrease strongly when the rotation frequency increases (see Fig. 4), and also depends on the assumption made on chemical equilibrium.

# 5 NUMERICAL SIMULATIONS OF GLITCHES

## 5.1 Numerical procedure

### 5.1.1 Computational scheme

Starting from two angular velocities $\Omega_n^0$ and $\Omega_p^0$ at the beginning of the glitch event (see Section 5.1.2), the evolution of the fluid rotation rates is computed from a series of equilibrium configurations (see Section 2.1), either keeping fixed the total baryon mass $M^B$ with $\mu_c^n = \mu_c^p$ or for constant partial baryon masses $M_n^B$ and $M_p^B$ (see Section 2.4). This means that, for given angular velocities

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Coupling coefficient $\hat{\epsilon}_p$ computed from (27) as a function of the rotation frequency for a 1.4 M$_{\odot}$ neutron star, assuming corotation and beta equilibrium (case ii). Results are shown in red (blue) for DDH($\delta$) EoS. The maximum rotation frequency considered here is 327 Hz, which corresponds to the highest frequency of observed glitching pulsars.</figcaption></figure>

corresponding to the instant under consideration, the mutual friction torque on the right-hand side of equation (17) is calculated at equilibrium, using the code described in Sourie et al. (2016). The moments of inertia involved in equation (17) are computed either through a fourth-order finite difference method or from a spectral interpolation based on Chebyshev polynomials. In the former, the moments of inertia are taken as constants during the glitch and evaluated at the rotation frequencies at the end of the glitch. On the contrary, in the latter, the moments of inertia are calculated for the angular velocities corresponding to the instant under consideration.

The angular velocities are evolved in time employing a two-step explicit Adams-Bashforth method using a time-step $\delta t \ll \tau_r$. For the different results given in Section 5.2, we typically consider time-steps of the order of $\delta t \simeq \tau_r/10^4 - \tau_r/10^2$, where $\tau_r$ is estimated from equation (24). Note that, for the different simulations performed, the total baryon mass $M^B$ and the total angular momentum $J$ are conserved during the glitch with a precision better than $10^{-10}$.

### 5.1.2 Initial conditions

The lag $\delta\Omega_0 = \Omega_n^0 - \Omega_p^0$ at the beginning of the glitch should be determined from the pinning of vortex lines during the pre-glitch evolution (see e.g. Haskell et al. 2012). Still, it is possible to get the initial angular velocities $\Omega_n^0$ and $\Omega_p^0$ from basic considerations, regardless of the physical processes that build up the lag and trigger the glitch, as described in the following.

We denote by $\Omega_n^f$ and $\Omega_p^f$ the angular velocities at the end of the spin up, i.e. when the post-glitch relaxation starts. By assuming that the system relaxes completely during the glitch rise, as suggested by equation (20), we have

$$ \Omega_n^f = \Omega_p^f = \Omega^f. \qquad (32) $$

The initial proton rotation frequency $\Omega_p^0$ can be simply deduced from a given choice of the amplitude

$$ \frac{\Delta\Omega}{\Omega} = \frac{\Omega^f - \Omega_p^0}{\Omega_p^0} \qquad (33) $$

of the glitch to be modelled, which is typically of the order of ~$10^{-6}$ for Vela (Dodson et al. 2007). Furthermore, the total angular
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left: time evolution of the superfluid and proton angular velocities for realistic input parameters corresponding to Vela ($f = 11.19$Hz, $\Delta\Omega/\Omega = 10^{-6}$). The gravitational mass is $M_G = 1.4 M_\odot$. The mutual friction parameter is $\bar{B} = 10^{-4}$. The dynamical evolution of the two fluids is plotted in terms of $(\Omega_X(t) - \Omega_p^0)/\Omega_p^0$. Results obtained from the DDH($\delta$) EoS are shown in solid (dashed) lines. Right: lag $\delta\Omega = \Omega_n - \Omega_p$ as a function of time for the same input parameters. Straight lines correspond to fits with exponential decaying laws.</figcaption></figure>

momentum $J = J_n + J_p$ being conserved during the whole glitch event (see equation 5), the last unknown $\Omega_n^0$ can be determined from

$$J(\Omega_n^0, \Omega_p^0) = J(\Omega_n^f, \Omega_p^f). \quad (34)$$

Since the glitch amplitudes are extremely small, the final angular velocities $\Omega_X^f$ are very close to the initial values $\Omega_X^0$. Expanding $J$ to the first order in the angular velocities is thus sufficient to determine $\Omega_n^0$ with very high accuracy:

$$\Omega_n^0 \approx \Omega^f \left( 1 + \frac{\hat{I}_p}{\hat{I}_n} \frac{\Delta\Omega}{\Omega} \right), \quad (35)$$

where the moments of inertia, defined by (13), are computed at the end of the glitch, i.e. for $\Omega_p^f = \Omega_n^f$.

The initial lag $\delta\Omega_0$ is therefore given by

$$\delta\Omega_0 = \Omega_n^0 - \Omega_p^0 \approx \Omega^f \frac{\hat{I}}{\hat{I}_n} \frac{\Delta\Omega}{\Omega}, \quad (36)$$

which leads to $\sim 7.10^{-5}$ rad s$^{-1}$ for the Vela pulsar, taking $\Omega^f/(2\pi) = 11.19$ Hz. We deduce that the lag between the two fluids always verifies the condition

$$\frac{\delta\Omega}{\Omega} \le \frac{\delta\Omega_0}{\Omega} \approx \frac{\Delta\Omega}{\Omega} \ll 1, \quad (37)$$

meaning that the deviation from corotation remains very small during the glitch event.

To summarize, the numerical simulations require the following macroscopic ingredients:

(i) the rotation rate $\Omega^f$ of the star,

(ii) its gravitational mass $M_G$,

(iii) the glitch amplitude $\Delta\Omega/\Omega$,

which can be potentially directly obtained from observations. In addition, the following microscopic inputs need to be specified:

(i) the EoS used to describe the interior of the star (for the adopted composition as discussed in Section 2.4),

(ii) the mutual friction parameter $\bar{B}$.

Contrary to the total baryon mass $M^B$, the gravitational mass $M_G$, which corresponds to the observed mass of the pulsar, should vary

during the glitch spin up. Nevertheless, we note that the change in $M_G$ associated with the angular momentum transfer is found to be smaller than a few $10^{-10}$ for the different tests performed with $\Delta\Omega/\Omega = 10^{-6}$.

## 5.2 Numerical results

We now present various results obtained from the numerical simulations described in the previous section. In the following, we mainly consider slowly rotating pulsars like Vela for which the assumption on the composition is unimportant (for the actual numerical calculations, we consider case ii; see Section 2.4). The impact of the chemical equilibrium on the evolution of more rapidly rotating neutron stars is discussed at the end of Section 5.2.1.

### 5.2.1 Dynamical evolution

In Fig. 5, we show the temporal evolution of the two angular velocities and the lag $\delta\Omega = \Omega_n - \Omega_p$ for conditions corresponding to the Vela pulsar. The gravitational mass of this pulsar being not (well) known, we have chosen the canonical value of $M_G = 1.4 M_\odot$ in this example. Moreover, the mutual friction parameter is arbitrarily fixed to $\bar{B} = 10^{-4}$. Results are shown in Fig. 5 for both EoSs studied in Sourie et al. (2016), namely the DDH and DDH$\delta$ EoSs, with only small differences between both. As can be seen in the right panel, the evolution of the lag can be very well described by an exponential law of the form (20), as expected from Section 3.2. For the present example, we find $\tau_r = 4.23$ s with the DDH EoS and $\tau_r = 2.92$ s with the DDH$\delta$ one. It should be remarked that these characteristic times, obtained from the time evolution of the lag, correspond indeed to the spin-up time-scales that could be measured from precise timing observations of glitches (see equation 22).

In order to study the dependence of the rise time on the different input parameters and to compare with the results given in Section 3.2, we have performed a series of simulations, varying in particular the mutual friction parameter $\bar{B}$, the rotation rate $\Omega$ of the star as well as its gravitational mass $M_G$. In Fig. 6, the spin-up time-scale is plotted with respect to the pulsar angular velocity $\Omega^f$, for two different gravitational masses and both EoSs, assuming a glitch amplitude $\Delta\Omega/\Omega = 10^{-6}$. The rise time turns out to be inversely
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Glitch rise time $\tau_r$ as a function of the pulsar rotation rate $f = \Omega/(2\pi)$, for a glitch amplitude $\Delta\Omega/\Omega = 10^{-6}$. Results obtained from the DDH($\delta$) EoS are shown in solid (dashed) lines. Configurations with two different gravitational masses are displayed, using $\tilde{B} = 10^{-4}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Glitch rise time $\tau_r$ plotted with respect to the gravitational mass $M_G$, for a star spinning at 11.19 Hz. The glitch parameters $\tilde{B}$ and $\Delta\Omega/\Omega$ are respectively fixed to $10^{-4}$ and $10^{-6}$. Results obtained from the DDH($\delta$) EoS are displayed with blue dots (yellow crosses).</figcaption></figure>

proportional to the pulsar rotation rate with a high accuracy, consistently with equation (24). A very small deviation from this simple behaviour can be seen for $f \gtrsim 100$ Hz, due to the strong increase of the moments of inertia with the angular velocity in this range of values (see Section 4.1).

For the different cases considered in Fig. 6, the numerical results are found to agree with values inferred from equation (24) with a precision better than $\sim 5 \times 10^{-6}$. This limit of accuracy comes from the numerical errors associated with dynamical simulations, which are dominated by the discretization in time and the precision with which the moments of inertia are computed. The reason why equation (24) gives such a precise estimate for the spin-up time-scale comes from the extremely small glitch amplitudes that are observed. The spin-up time-scale can thus be very precisely estimated from equation (24) by merely computing stationary configurations and ignoring the change in the moments of inertia during the glitch.

All other parameters fixed, the rise time increases with the gravitational mass of the star, as highlighted in Fig. 7. Several reasons can be invoked to explain this fact. First, the proton fraction $x_p = n_p/(n_n + n_p)$ and therefore the ratio $\hat{l}_p/\hat{l}$ are strongly

increasing with the mass of the star (see Appendix D). As the neutron fraction decreases, the transfer of angular momentum becomes longer (see equation 30). Moreover, the coupling coefficient $\hat{e}_p$ and the quantity $\zeta$ are also found to decrease significantly as the gravitational mass increases (see Section 4). The transfer of angular momentum is thus slowed down (equation 30). It should also be noticed that, even if the coupling is much stronger for DDH (see the right panel of Fig. 3), the spin-up time-scale is systematically longer with the DDH EoS than with the DDH$\delta$ one because the proton fraction and thus the ratio $\hat{l}_p/\hat{l}$ predicted by this EoS are much higher.

Finally, a few tests have been also performed to study the influence of the assumption concerning the evolution of the chemical composition during the glitch rise (see Section 2.4). Whether considering constant individual baryon masses (case i) or fixed total baryon mass with beta equilibrium at the centre (case ii) leads to negligible differences for the glitch rise time at low rotation frequencies. For instance, the deviation is lower than $\sim 7 \times 10^{-5}$ for a star rotating at 11.19 Hz within both EoSs, assuming a small glitch amplitude and fixed moments of inertia. Nevertheless, the impact of the assumption on the chemical equilibrium increases sharply with the angular velocity: for both EoSs, the discrepancy between cases (i) and (ii) is of the order of $\sim 2 \times 10^{-3}$ for 65 Hz and $\sim 5 \times 10^{-2}$ for 327 Hz. Still, it is important to note that the influence of the assumption relative to chemical equilibrium is much smaller than the dependence of the rise time on any other input parameters of our model, such as masses or rotation rates. This is the reason why we consider only case (ii) in most of the results presented in the present paper.

### 5.2.2 Contribution of general relativity

To study the global contribution of general relativity to the spin-up time-scale, we compare the rise times obtained within both relativistic and Newtonian frameworks. For simplicity, we consider polytropic EoSs, as implemented by Prix et al. (2005).

In Fig. 8, the relative differences on $\tau_r$ are plotted with respect to the compactness parameter (26), obtained for a star rotating at 10 Hz by varying its (gravitational) mass. We consider two different EoSs (referred to as EoS I and II), with small and respectively strong entrainment effects, associated with the following Lagrangian densities $\Lambda = -\mathcal{E}$ with

$$ \mathcal{E} = \frac{1}{2} \kappa_n n_n^2 + \frac{1}{2} \kappa_p n_p^2 + \kappa_{np} n_n n_p + \kappa_{\Delta} n_n n_p \Delta^2 + (\rho c^2), \quad (38) $$

and

$$ \mathcal{E} = \frac{1}{2} \kappa_n n_n^{2.1} + \frac{1}{2} \kappa_p n_p^{2.3} + \kappa_{\Delta} n_n n_p \Delta^2 + (\rho c^2), \quad (39) $$

respectively, where $\rho = m_n n_n + m_p n_p$, $\Delta$ stands for the relative speed between the fluids and we have used the same notation as in Prix et al. (2005). Note that the rest mass energy density is only present in the general relativistic case (more details on the differences between Newtonian and relativistic computations can be found in section IV-A of Prix et al. 2005). The different parameters $\kappa_n$, $\kappa_p$ and $\kappa_{np}$ (see Table 1) are chosen in order to reproduce ‘realistic’ values for the mass, radius and proton fraction $x_p$ of the stars: for instance, EoS I leads to a constant proton fraction throughout the star with $x_p = 0.05$, whereas a varying proton fraction in the range $x_p \simeq 0.05-0.1$ is obtained with EoS II for a 1.4 M$_\odot$ (relativistic) neutron star spinning at 10 Hz. Both EoSs predict an external circumferential radius in the equatorial plane $R_{c,eq} \simeq 13$ km
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Relative differences between general relativistic and Newtonian rise times (blue solid lines) for a glitch amplitude $\Delta\Omega/\Omega = 10^{-6}$ as functions of the (relativistic) compactness parameter, for a star spinning at 10 Hz. Results are shown for the two polytropic EoSs described in the text: the left (right) panel corresponds to EoS I (EoS II). Other curves highlight the contribution of general relativity in the computation of the different terms involved in the spin-up time-scale (equation 30): results accounting only for the general relativistic corrections on the ratio $\hat{l}_p/\hat{l}$ are displayed in orange dash-dotted lines, whereas general relativistic corrections on both $\hat{l}_p/\hat{l}$ and $\zeta$ are considered in the red dashed lines. The green dotted lines include the additional corrections on the entrainment parameters $\tilde{\epsilon}_X$ but do not account for the frame-dragging contribution to the coupling coefficients $\hat{\epsilon}_X$ (29).</figcaption></figure>

<table><caption>**Table 1.** Polytropic parameters defining EoSs I (equation 38) and II (equation 39). For EoS I, $\kappa_n$, $\kappa_p$ and $\kappa_{np}$ are given in units of $\rho_0 c^2 n_0^{-2}$ where $n_0 = 0.1 \text{ fm}^{-3}$ and $\rho_0 = 1.66 \times 10^{17} \text{ kg.m}^{-3}$, whereas $\kappa_n$ and $\kappa_p$ are respectively expressed in units of $\rho_0 c^2 n_0^{2.1}$ and $\rho_0 c^2 n_0^{-2.3}$ for the second EoS. For both EoSs, $\kappa_\Delta$ is expressed in units of $\rho_0 n_0^{-2}$.</caption><thead><tr><th></th><th>κ<sub>n</sub></th><th>κ<sub>p</sub></th><th>κ<sub>np</sub></th><th>κ<sub>Δ</sub></th></tr></thead><tbody><tr><td>EoS I</td><td>0.05</td><td>0.5</td><td>0.025</td><td>0.02</td></tr><tr><td>EoS II</td><td>0.046</td><td>1.4</td><td>–</td><td>0.1</td></tr></tbody></table>

for a relativistic neutron star with $M_G = 1.4 M_\odot$. The entrainment contribution is included through the coefficient $\kappa_\Delta$. The forms of the EoSs are taken consistently with the fact that entrainment effects should vanish when one of the fluids disappears (see Sourie et al. 2016). For EoS I, we take $\kappa_\Delta = 0.02$ in order to satisfy all the required stability conditions (Chamel & Haensel 2006). This choice leads to $\tilde{\epsilon}_p \simeq 0.07$ for a (relativistic) neutron star spinning at 10 Hz, with $M_G = 1.4 M_\odot$. Although this value is quite small compared to realistic EoSs (Fig. 3), it still corresponds to the outer core of neutron stars, where entrainment effects are nearly vanishing. On the contrary, requiring stability for EoS II, a value of $\kappa_\Delta = 0.1$ leads to $\tilde{\epsilon}_p \simeq 0.29$ for the same mass and spin, which happens to be much more realistic.

General relativity is expected to play a role in determining the moments of inertia $\hat{l}_p$ and $\hat{l}$, the quantity $\zeta$ and the coupling coefficients $\hat{\epsilon}_X$ involved in the spin-up time-scale (see equation 30). To highlight the contribution of general relativity on these different terms, the following quantities are plotted in Fig. 8 for EoS I (left panel) and EoS II (right panel):

$$ (i) \quad \frac{\tau_{\hat{l}}^{\text{GR}} - \tau_{\hat{l}}^{\text{newt}}}{\tau_{\hat{l}}^{\text{GR}}} = 1 - \frac{\hat{l}_{\text{p}}^{\text{newt}}}{\hat{l}_{\text{p}}^{\text{newt}}} \frac{\hat{l}_{\text{GR}} - \hat{l}_{\text{p}}^{\text{newt}} - \hat{l}_{\text{n}}^{\text{newt}}}{1 - \hat{l}_{\text{GR}} - \hat{l}_{\text{p}}^{\text{GR}}} \zeta, \quad \text{denoting Newtonian and general relativistic quantities by 'newt' and 'GR' respectively (see equation 29 and Appendix C for the definitions of the total coupling coefficients $\hat{\epsilon}_X$ and the entrainment parameters $\tilde{\epsilon}_X$),} $$

$$ (ii) \quad 1 - \frac{\hat{l}_{\text{p}}^{\text{newt}}}{\hat{l}_{\text{p}}^{\text{newt}}} \frac{\hat{l}_{\text{GR}}}{\hat{l}_{\text{p}}^{\text{GR}}}, \quad \text{accounting only for the relativistic corrections on the ratio } \hat{l}_{\text{p}}/\hat{l}, $$

$$ (iii) \quad 1 - \frac{\hat{l}_{\text{p}}^{\text{newt}}}{\hat{l}_{\text{p}}^{\text{newt}}} \frac{\hat{l}_{\text{GR}}}{\hat{l}_{\text{p}}^{\text{GR}}} \zeta, \quad \text{considering both relativistic corrections on the moments of inertia and } \zeta, $$

$$ (iv) \quad 1 - \frac{\hat{l}_{\text{p}}^{\text{newt}}}{\hat{l}_{\text{p}}^{\text{newt}}} \frac{\hat{l}_{\text{GR}} - \hat{l}_{\text{p}}^{\text{newt}} - \hat{l}_{\text{n}}^{\text{newt}}}{1 - \hat{l}_{\text{GR}} - \hat{l}_{\text{p}}^{\text{GR}}} \zeta, \quad \text{taking into account relativistic corrections on the moments of inertia, } \zeta \text{ and the entrainment parameters } \tilde{\epsilon}_X. $$

As expected, general relativistic corrections tend to zero, when the compactness parameter decreases. Concerning EoS I, general relativistic corrections on the ratio $\hat{l}_p/\hat{l}$ are found to be extremely small (see the left panel of Fig. 8). This is due to the fact that $x_p$ is constant throughout the whole star within this EoS and $\hat{l}_p/\hat{l} = x_p$ in both Newtonian and general relativistic frameworks, for a slowly rotating star in beta equilibrium and in the limit of vanishing lag between the fluids. In Newtonian gravity, $\zeta$ is equal to 1. Since general relativity leads to $\zeta \lesssim 1$ (see Section 4.1), this quantity acts to lengthen the rise time. Moreover, the general relativistic entrainment parameters are found to be much higher than their Newtonian counterparts, because higher densities are reached when general relativity is considered. These general relativistic corrections on entrainment tend to lower the rise time. Finally, frame-dragging contribution to the coupling coefficients $\hat{\epsilon}_X$ (29) also leads to a longer spin-up time-scale (see equation 31). Similar remarks apply to EoS II but with two differences. First, general relativity slightly modifies the ratio of the moments of inertia (see the right panel of Fig. 8). Furthermore, entrainment effects are much more important within this EoS than in EoS I, leading to a much larger reduction of the rise time. To summarize, general relativistic corrections on the different terms involved in the spin-up time-scale are found to be roughly of the same order of magnitude but depend strongly on the EoS considered. In particular, the frame-dragging contribution to the fluid couplings is found to be important.

For values of the compactness parameter relevant for neutron stars, i.e. $\Xi \sim 0.15-0.20$, these two EoSs predict that an error of the order of ~20–40 per cent is made on the rise time by using Newtonian gravity instead of general relativity, as can be seen in Fig. 8. It is therefore necessary to account for general relativistic effects in order to get precise results on the spin-up time-scales. Furthermore, it should be mentioned that these errors also depend
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Lines of constant rise time $\tau_r$ in a $\tilde{B}-M_G$ diagram, assuming $f = 11.19$ Hz and $\Delta\Omega/\Omega = 10^{-6}$. Results obtained from the DDH($\delta$) EoS are shown in solid (dashed) lines. The $\tau_r = 30$ s line refers to the upper limit from Dodson et al. (2007).</figcaption></figure>

significantly on the rotation rate considered. For instance, the relative difference $(\tau_r^{GR} - \tau_r^{newt}) / \tau_r^{GR}$ obtained for a 1.4 M$_\odot$ neutron star varies from ~30 per cent at 10 Hz to ~10 per cent at 327 Hz, using EoS I.

## 5.3 Astrophysical considerations

In Fig. 9, lines of constant spin-up time $\tau_r$ are displayed in the $M_G-\tilde{B}$ plane, using DDH and DDH$\delta$ EoSs. These results are plotted for $f = 11.19$ Hz, which corresponds to the Vela pulsar. Considering the current upper limit $\tau_r < 30$ s (Dodson et al. 2007), the mutual friction parameter $\tilde{B}$ should be higher than ~$10^{-5}$ to explain Vela glitches. Using equation (10) with $\tilde{B} \approx \tilde{B}$, this limit implies that the averaged drag-to-lift ratio $R$ should verify $10^{-5} \lesssim R \lesssim 10^5$, which is not very constraining regarding the diversity of dissipative mechanisms that could give rise to mutual friction and the corresponding microscopic uncertainties. Note that similar conclusions were reached by Glampedakis, Andersson & Jones (2009) from the analysis of post-glitch relaxation data. Nevertheless, as the dependence of the spin-up time on $M_G$ is much less pronounced than on $\tilde{B}$, future more stringent observational limits on $\tau_r$ shall put interesting constraints on the process governing the angular momentum transfer during the spin up.

Finally, since $\tilde{B} \le 1/2$, the shortest possible rise time (see equation 25) associated with Vela glitches, $\tau_r \simeq 0.6-0.8$ ms for a 1.4 M$_\odot$ star (see Fig. 10), is found to be comparable with the hydrodynamical time-scale. Nevertheless, the actual value of $\tau_r$ is presumably much longer in view of current estimates of the mutual friction parameters. Therefore, the whole dynamical evolution of star during the spin up can be accurately computed by considering a sequence of stationary configurations only.

## 6 GRAVITATIONAL WAVES

In this section, we study the amount of gravitational waves emitted through the time evolution of the mass quadrupole of the star, as a consequence of the changes in the fluid angular velocities during the spin up. Because of the small deviation from spherical symmetry, the variations of the mass quadrupole are expected to be very small. Thus, we focus on the DDH EoS, with which the equilibrium

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Glitch rise time $\tau_r$ as a function of the drag-to-lift ratio $R$, for a star spinning at 11.19 Hz with $M_G = 1.4$ M$_\odot$. The glitch amplitude considered here is $\Delta\Omega/\Omega = 10^{-6}$. Results obtained from the DDH($\delta$) are shown in blue (red). The 30 s upper limit is also displayed.</figcaption></figure>

configurations obtained are sufficiently accurate to highlight these tiny variations (see Sourie et al. 2016).

### 6.1 Mass quadrupole

To compute the coordinate-independent mass quadrupole moment $Q$ of the star, we follow the prescription given by equation (11) of Pappas & Apostolatos (2012) (see also Friedman & Stergioulas 2013). The sign of the mass quadrupole moment is chosen such that $Q < 0$ for an oblate spheroid. Using the DDH EoS, $|Q|$ decreases as time evolves, which means that the shape of the star is getting less and less oblate.

For realistic glitch amplitudes, numerical results show that the time evolution of the mass quadrupole can be very well approximated by the formula

$$Q(t) = \Delta Q (1 - \exp(-t/\tau_r)) + Q_0, \quad (40)$$

where $\tau_r$ corresponds to the spin-up time-scale discussed in the previous sections. All other input parameters fixed, the variation $\Delta Q$ of the quadrupole moment during the glitch event is found to verify $\Delta Q \propto \Omega^2 \times \Delta\Omega/\Omega$. Whereas the dynamical evolution of $Q$ depends indeed on the mutual friction parameter through $\tau_r$, the variation $\Delta Q$ does not depend on $\tilde{B}$ because this latter does not play any role in determining the initial and final equilibrium configurations. On the other hand, the assumption on chemical equilibrium (Section 2.4) affects $\Delta Q$: assuming $\Delta\Omega/\Omega = 10^{-4}$ and $M_G = 1.4$ M$_\odot$, the relative differences on $\Delta Q$ between cases (i) and (ii) are of the order of ~1 per cent for $f = 65$ Hz and ~15 per cent for $f = 327$ Hz.

### 6.2 Gravitational wave amplitudes

At leading order in a multipolar expansion, the gravitational radiation field is given by the so-called quadrupole formula (see e.g. Bonazzola & Gourgoulhon 1996; Le Tiec & Novak 2016). Since the star remains axisymmetric during the spin up, the $\times$ polarization of the wave strain vanishes. On the other hand, the $+$ polarization reads

$$h_+(t) = -\frac{3}{2} \frac{G}{Dc^4} \sin^2 i \times \ddot{Q} (t - D/c), \quad (41)$$

where $D$ is the distance to the pulsar and $i$ denotes the angle between the rotation axis of the star and the direction from the star's centre
to the Earth. Using equation (40), the gravitational wave amplitude (equation 41) is given by

$$h_{+}(t) = h_{0} \sin^{2} i \times \exp\left(-\frac{(t - D/c)}{\tau_{r}}\right), \quad (42)$$

where the characteristic amplitude $h_0$ is defined as

$$h_0 = \frac{3}{2} \frac{G}{Dc^4} \frac{\Delta Q}{\tau_r^2}. \quad (43)$$

In the frequency domain, the corresponding characteristic strain $h_c(f_{GW})$ is given from the Fourier transform $\tilde{h}(f_{GW})$ of the signal $h_+(t)$ through the relation

$$h_c(f_{GW}) = 2f_{GW} \times |\tilde{h}(f_{GW})|, \quad (44)$$

see e.g. Moore, Cole & Berry (2015). Using (42), the characteristic strain reads

$$h_c(f_{GW}) = \frac{h_0 \sin^2 i}{\pi} \times \frac{f_{GW}/f_0}{\sqrt{1 + (f_{GW}/f_0)^2}}, \quad (45)$$

where we have introduced the characteristic frequency

$$f_0 = \frac{1}{2\pi\tau_r}. \quad (46)$$

Assuming $M_G = 1.4 M_\odot$ and beta equilibrium at the centre of the star, the gravitational wave characteristic amplitude (43) and frequency (46) obtained from the DDH EoS can be well-fitted by the following expressions:

$$h_0 \approx 1.0 \times 10^{-37} \left( \frac{D}{1 \text{ kpc}} \right)^{-1} \left( \frac{\bar{B}}{10^{-3}} \right)^2 \times \left( \frac{\Omega}{10^2 \text{ rads}^{-1}} \right)^4 \left( \frac{\Delta\Omega/\Omega}{10^{-6}} \right), \quad (47)$$

and

$$f_0 \approx 0.535 \left( \frac{\bar{B}}{10^{-3}} \right) \left( \frac{\Omega}{10^2 \text{ rads}^{-1}} \right) \text{Hz}, \quad (48)$$

provided that $\Delta\Omega/\Omega \ll 1$. For low rotation frequencies, typically $f \lesssim 65$ Hz, the two expressions (47) and (48) approximate $h_0$ and $f_0$ with a precision better than a few per cent and $\sim 0.1$ per cent, respectively. For a star spinning at 327 Hz, the precision of these estimates is reduced to $\sim 10$ per cent. Note that the prefactors involved in equations (47) and (48) are smaller for higher gravitational masses, mainly because the rise times are longer: for instance, we get $3.7 \times 10^{-38}$ and $0.357$ for $M_G = 2 M_\odot$.

It should be remarked here that the quadrupole formula (41) is only valid in the slow-motion approximation, meaning that the frequency $f_{GW}$ of the gravitational wave emitted should satisfy the condition $f_{GW} \times R \ll c$, where R is the characteristic size of the emitter. Taking $R \simeq 10^4$ m for the radius of the star, this leads to $f_{GW} \ll 3$ kHz, which in view of equation (48) is well verified for glitch pulsars. Furthermore, we can easily check that the energy lost by gravitational waves emission is completely negligible with respect to that associated with mutual friction, such that equation (5) is valid.

Applying equations (47) and (48) to the Vela pulsar, for which $\Delta\Omega/\Omega = 10^{-6}$, $f = 11.19$ Hz and $D \simeq 287$ pc (Dodson et al. 2003), the constraint on the mutual friction parameter $\bar{B}$ discussed in Section 5.3, i.e. $10^{-5} < \bar{B} < 0.5$, leads to

$$h_0 \sim 10^{-41} - 10^{-32} \quad \text{and} \quad f_0 \sim 4 \text{ mHz} - 200 \text{ Hz}. \quad (49)$$

Although the peak frequency $f_0$ is thus predicted to lie in the sensitivity bands of Advanced LIGO and Advanced Virgo (Acernese et al. 2015; LIGO Scientific Collaboration et al. 2015; Abbott et al. 2016), the corresponding gravitational-wave signal is too weak to be detectable with present detectors. However, other mechanisms associated with glitches such as Ekman pumping, which we have not taken into account in this work, might lead to a much stronger gravitational wave signal (van Eysden & Melatos 2008; Bennett, van Eysden & Melatos 2010).

## 7 CONCLUSION

In this paper, we have studied in detail the impact of general relativity on the global dynamics of giant pulsar glitches as observed in Vela. We have carried out numerical simulations of the spin up triggered by the sudden unpinning of quantized vortices. To this end, we have computed the exchange of angular momentum between the neutron superfluid in the core and the rest of the star within a two-fluid model including neutron-proton entrainment effects. Both fluids were assumed to be coupled by mutual friction arising from dissipative forces acting on individual vortices. Since the hydrodynamical time-scale is typically much smaller than the glitch rise time, we have described the time evolution of the two fluids by a sequence of quasi-stationary axisymmetric rigidly rotating configurations following Sourie et al. (2016). We have calculated the mutual friction torque considering straight vortices arranged on a regular array, following Langlois et al. (1998).

In order to get some physical insight, we first solved analytically the dynamical equations by expressing the change in the lag as $\delta\Omega/\delta\Omega \approx -1/\tau_r$, where the characteristic spin-up time-scale $\tau_r$ can be expressed in a form similar to that obtained in the Newtonian limit (see e.g. Carter 2001; Sidery et al. 2010). However, general relativity not only changes the structure of the star, but also impacts the fluid dynamics. In particular, frame-dragging effects induce additional fluid couplings of the same form as the entrainment arising solely from neutron-proton interactions. For all these reasons, general relativity can substantially change the glitch rise time.

To test the validity of this analytical approach and to assess the importance of general relativity, we have also numerically solved the equations governing the transfer of angular momenta. For this purpose, two different kinds of inputs are needed: macroscopic quantities (the rotation frequency of the star, the glitch amplitude and the neutron-star mass) and microscopic properties (the EoS and the stellar-averaged mutual friction coupling $\bar{B}$). We have explored in detail various stellar configurations, using two different relativistic mean-field EoSs and considering the observed properties of glitching pulsars. The results obtained by numerical simulations are found to be very well reproduced by the analytical approximation. In particular, the glitch rise time $\tau_r$ can thus be expressed in terms of the moments of inertia of the fluids, the stellar rotation rate and $\bar{B}$, which can be obtained from stationary configurations. Furthermore, we have studied the effects of general relativity on $\tau_r$ by using two different polytropic EoSs of the kind previously introduced by Prix et al. (2005). The effects of general relativity on both the structure of the star and the fluid couplings are found to be important and therefore realistic simulations of the global glitch dynamics should be carried out in full general relativity. Depending on the stellar compactness and on the rotation rate, the errors incurred by using Newtonian gravity instead of general relativity are found to be very sensitive to the adopted EoS, and amount to $\sim 20-40$ per cent. These errors, however, might not be the dominant source of uncertainties. In particular, neutron superfluid vortices may not be arranged on
a regular array parallel to the rotation axis, as assumed here. The dynamics of superfluid vortices and proton flux tubes remain highly uncertain, and warrant further studies.

Considering the current upper limit $\tau_r < 30$ s (Dodson et al. 2007), we have found that the mutual friction parameter $\bar{B}$ should be higher than $\sim 10^{-5}$ to explain Vela glitches. Since $\bar{B}$ represents the average over the whole star, the mutual friction coupling $\bar{B}$ might be locally much stronger ($\bar{B} \sim 10^{-4}-10^{-3}$) as discussed for instance by Sedrakian (2005) and Haskell et al. (2014). In any case, since the actual value of $\tau_r$ is found to be much longer than the hydrodynamical time-scale for current estimates of the mutual friction forces, the whole dynamical evolution of star during the spin up can be accurately computed by considering a sequence of stationary configurations only.

We have also determined the amount of gravitational radiation emitted by the star during the spin up. For this purpose, we have studied the time variation of the mass quadrupole moment of the star resulting from changes in the fluid angular velocities. Using the quadrupole formula, we have numerically computed the characteristic amplitudes and frequencies associated with glitch events. Although the peak frequencies are found to lie in the sensitivity bands of current interferometers like Advanced LIGO, the corresponding amplitudes are too small for the gravitational waves to be detected. Their observations would require to improve the sensitivity by orders of magnitude. In particular, the characteristic amplitude for Vela is estimated to be at most $\sim 10^{-32}$ for the (unrealistic) value $\bar{B} = 1/2$. If existing, the most promising sources would thus be pulsars rotating much more rapidly than Vela and undergoing high amplitude glitches.

Although glitches are unlikely to be detected through gravitational waves, the Low Frequency Array (LOFAR) radio telescope (Stappers et al. 2011) and the future Square Kilometre Array (SKA) (Watts et al. 2015) will be able to observe the spin up with unprecedented accuracy. It would thus lead to much more stringent constraints on the characteristic time $\tau_r$ and thereby on the underlying glitch mechanism. This calls for more realistic models of glitching pulsars including the crust magnetoelasticity and superfluidity (whose formalism has been already developed; see e.g. Carter, Chachoua & Chamel 2006b; Carter & Samuelsson 2006), and accounting for the local dynamics of quantized vortices.

## ACKNOWLEDGEMENTS

We would like to thank Ismaël Cognard for instructive discussions and Armen Sedrakian for interesting suggestions. This work has been partially funded by the ‘Gravitation et physique fondamentale’ action of the Observatoire de Paris (France), the Fonds de la Recherche Scientifique - FNRS (Belgium) under grant no CDR J.0187.16, the PHC Tournesol (no 35904ZJ) scientific cooperation program between France and Belgium, and the European COST action MP1304 ‘NewCompstar’.

## REFERENCES

* Abbott B. P. et al., 2016, Living Rev. Relativ., 19, 1

* Abo-Shaeer J. R., Raman C., Vogels J. M., Ketterle W., 2001, Science, 292, 476

* Acernese F. et al., 2015, Class. Quant. Gravity, 32, 024001

* Akbal O., Güercinoğlu E., Şaşmaz Muş S., Alpar M. A., 2015, MNRAS, 449, 933

* Alpar M. A., 2011, AIP Conf. Proc. Vol. 1379, ASTROPHYSICS OF NEUTRON STARS 2010: A Conference in Honor of M. Ali Alpar. Am. Inst. Phys., New York, p. 166

* Alpar M. A., Pines D., Anderson P. W., Shaham J., 1984a, ApJ, 276, 325

* Alpar M. A., Anderson P. W., Pines D., Shaham J., 1984b, ApJ, 278, 791

* Alpar M. A., Langer S. A., Sauls J. A., 1984c, ApJ, 282, 533

* Alpar M. A., Chau H. F., Cheng K. S., Pines D., 1993, ApJ, 409, 345

* Alpar M. A., Chau H. F., Cheng K. S., Pines D., 1996, ApJ, 459, 706

* Anderson P. W., Itoh N., 1975, Nature, 256, 25

* Andersson N., Comer G. L., 2001, Class. Quant. Gravity, 18, 969

* Andersson N., Sidery T., Comer G. L., 2006, MNRAS, 368, 162

* Andersson N., Sidery T., Comer G. L., 2007, MNRAS, 381, 747

* Andersson N., Glampedakis K., Ho W. C. G., Espinoza C. M., 2012, Phys. Rev. Lett., 109, 241103

* Andreev A. F., Bashkin E. P., 1976, Soviet J. Exper. Theor. Phys., 42, 164

* Antonelli M., Pizzochero P., 2017, MNRAS, 464, 721

* Antonopoulou D., Weltevrede P., Espinoza C. M., Watts A. L., Johnston S., Shannon R. M., Kerr M., 2015, MNRAS, 447, 3924

* Archibald R. F. et al., 2013, Nature, 497, 591

* Avancini S. S., Brito L., Marinelli J. R., Menezes D. P., de Moraes M. M. W., Providência C., Santos A. M., 2009, Phys. Rev. C, 79, 035804

* Baym G., Pines D., 1971, Ann. Phys., 66, 816

* Baym G., Pethick C., Pines D., 1969a, Nature, 224, 673

* Baym G., Pethick C., Pines D., Ruderman M., 1969b, Nature, 224, 872

* Bennett M. F., van Eysden C. A., Melatos A., 2010, MNRAS, 409, 1705

* Bocquet M., Bonazzola S., Gourgoulhon E., Novak J., 1995, A&A, 301, 757

* Bonazzola S., Gourgoulhon E., 1996, A&A, 312, 675

* Bulgac A., Forbes M. M., Sharma R., 2013, Phys. Rev. Lett., 110, 241102

* Carter B., 1975, Ann. Phys., 95, 53

* Carter B., 2001, in Blaschke D., Glendenning N. K., Sedrakian A., eds, Lecture Notes in Physics Vol. 578, Physics of Neutron Star Interiors. Springer-Verlag, Berlin, p. 54

* Carter B., Quintana H., 1975, Ann. Phys., 95, 74

* Carter B., Samuelsson L., 2006, Class. Quant. Gravity, 23, 5367

* Carter B., Langlois D., Sedrakian D. M., 2000, A&A, 361, 795

* Carter B., Chamel N., Haensel P., 2005, Nucl. Phys. A, 748, 675

* Carter B., Chamel N., Haensel P., 2006a, Int. J. Modern Phys. D, 15, 777

* Carter B., Chachoua E., Chamel N., 2006b, Gen. Relativ. Gravitation, 38, 83

* Chamel N., 2004, PhD thesis, Université Paris VI, France

* Chamel N., 2005, Nucl. Phys. A, 747, 109

* Chamel N., 2012, Phys. Rev. C, 85, 035801

* Chamel N., 2013, Phys. Rev. Lett., 110, 011101

* Chamel N., Carter B., 2006, MNRAS, 368, 796

* Chamel N., Haensel P., 2006, Phys. Rev. C, 73, 045802

* Chatterjee D., Elghozi T., Novak J., Oertel M., 2015, MNRAS, 447, 3785

* Chau H. F., Cheng K. S., 1993, Phys. Rev. B, 47, 2707

* Datta B., Alpar M. A., 1993, A&A, 275, 210

* Delsate T., Chamel N., Gürlebeck N., Fantina A. F., Pearson J. M., Ducoin C., 2016, Phys. Rev. D, 94, 023008

* Dodson R. G., McCulloch P. M., Lewis D. R., 2002, ApJ, 564, L85

* Dodson R., Legge D., Reynolds J. E., McCulloch P. M., 2003, ApJ, 596, 1137

* Dodson R., Lewis D., McCulloch P., 2007, Ap&SS, 308, 585

* Eichler D., Shaisultanov R., 2010, ApJ, 715, L142

* Epstein R. I., 1988, ApJ, 333, 880

* Epstein R. I., Baym G., 1992, ApJ, 387, 276

* Espinoza C. M., Lyne A. G., Stappers B. W., Kramer M., 2011, MNRAS, 414, 1679

* Fantina A. F., Chamel N., Pearson J. M., Goriely S., 2013, A&A, 559, A128

* Friedman J. L., Stergioulas N., 2013, Rotating Relativistic Stars. Cambridge Univ. Press, Cambridge, p. 429

* Ginzburg V. L., Kirzhnits D. A., 1965, Sov. Phys. JETP, 20, 1346

* Glampedakis K., Andersson N., Jones D. I., 2009, MNRAS, 394, 1908

* Glampedakis K., Andersson N., Samuelsson L., 2011, MNRAS, 410, 805

* Gnedin O. Y., Yakovlev D. G., Potekhin A. Y., 2001, MNRAS, 324, 725

* Gourgoulhon E., 2010, preprint (arXiv:1003.5015)
## APPENDIX A: ANGULAR MOMENTUM TRANSFERS IN TWO-FLUID NEWTONIAN MODEL

In Newtonian gravity, the fluid angular momenta read

$$ \left\{ \begin{aligned} J_n &= I_n (1 - \bar{\epsilon}_n) \Omega_n + I_n \bar{\epsilon}_n \Omega_p, \\ J_p &= I_p (1 - \bar{\epsilon}_p) \Omega_p + I_p \bar{\epsilon}_p \Omega_n, \end{aligned} \right. \qquad (A1) $$

see appendix A of Sourie et al. (2016). The moments of inertia $I_X$ involved in (A1) are given by the classical formula

$$ I_X = \int_{\Sigma_i} \rho_X r^2 \sin^2 \theta d^3\Sigma_i, \qquad (A2) $$

where $\rho_X$ is the mass density of fluid X and $d^3\Sigma_i$ stands for the volume element of flat space-time, while $r$ and $\theta$ refer to the radial and polar coordinates, respectively. The quantities $\bar{\epsilon}_X$, which characterize entrainment, are defined as

$$ \bar{\epsilon}_X = \frac{\int_{\Sigma_i} \epsilon_X \rho_X r^2 \sin^2 \theta d^3\Sigma_i}{\int_{\Sigma_i} \rho_X r^2 \sin^2 \theta d^3\Sigma_i}, \qquad (A3) $$

where the entrainment parameter $\epsilon_X$ is given by

$$ \epsilon_X = \frac{2\alpha}{\rho_X}, \qquad (A4) $$

and $\alpha$ is defined as in Prix et al. (2005). The mean entrainment parameters are related through the relation $I_n \bar{\epsilon}_n = I_p \bar{\epsilon}_p$. More details can be found in appendix A of Sourie et al. (2016).

### A1 Mutual friction torque

In the Newtonian limit, the geometric term $\chi^2_{\perp}$ (equation 9) appearing in equation (12) is simply given by

$$ \chi^2_{\perp} = r^2 \sin^2 \theta. \qquad (A5) $$

Moreover, in case of *constant* entrainment parameters throughout the star, the non-relativistic neutron vorticity $\omega_n$ (7) reads

$$ \omega_n = 2m_n [\Omega_n + \epsilon_n (\Omega_p - \Omega_n)], \qquad (A6) $$

where $m_n$ is the neutron mass. Considering that in the non-relativistic limit the volume element $d^3\Sigma$ is $d^3\Sigma_i = r^2 \sin\theta dr d\theta d\varphi$ and $\Gamma_n = 1$, the mutual friction torque (12) thus reduces to

$$ \Gamma_{\text{mf}}^{\text{newt}} = -\vec{B} I_n \omega_n \delta \Omega, \qquad (A7) $$

where the superfluid vorticity per unit mass $\omega_n$ reads

$$ \omega_n = \frac{\omega_n}{m_n} = 2 [\Omega_n + \epsilon_n (\Omega_p - \Omega_n)]. \qquad (A8) $$

The Newtonian limit (equation A7) corresponds to equation (58) of Sidery et al. (2010). The quantity $\kappa$ introduced in equation (18) is given by

$$ \kappa = I_n \omega_n. \qquad (A9) $$

To a very good approximation, the difference between $\Omega_n$ and $\Omega_p$ can be neglected, so that $\omega_n \approx 2\Omega_n$ and therefore

$$ \kappa \approx 2 I_n \Omega_n. \qquad (A10) $$

Gourgoulhon E., ed., 2012, 3+1 Formalism in General Relativity Lecture Notes in Physics Vol. 846. Springer-Verlag, Berlin

Güercinoğlu E., Alpar M. A., 2014, ApJ, 788, L11

Haensel P., Zdunik J. L., Bejger M., Lattimer J. M., 2009, A&A, 502, 605

Haensel P., Bejger M., Fortin M., Zdunik L., 2016, Eur. Phys. J. A, 52, 59

Haskell B., Melatos A., 2015, Int. J. Modern Phys. D, 24, 1530008

Haskell B., Pizzochero P. M., Sidery T., 2012, MNRAS, 420, 658

Haskell B., Glampedakis K., Andersson N., 2014, MNRAS, 441, 1662

Jones P. B., 1990, MNRAS, 243, 257

Jones P. B., 1992, MNRAS, 257, 501

Keith M. J., Shannon R. M., Johnston S., 2013, MNRAS, 432, 3080

Langlois D., Sedrakian D. M., Carter B., 1998, MNRAS, 297, 1189

Larson M. B., Link B., 2002, MNRAS, 333, 613

Le Tiec A., Novak J., 2016, in Auger G., Plagnol E., eds, An Overview of Gravitational Waves: Theory and Detection, World Scientific Press, Singapore

LIGO Scientific Collaboration et al., 2015, Class. Quant. Gravity, 32, 074001

Link B., Epstein R. I., Lattimer J. M., 1999, Phys. Rev. Lett., 83, 3362

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

Martin N., Urban M., 2016, preprint (arXiv:1606.01126)

Moore C. J., Cole R. H., Berry C. P. L., 2015, Class. Quant. Gravity, 32, 015014

Page D., Lattimer J. M., Prakash M., Steiner A. W., 2013, in Bennemann K. H., Ketterson J. B., eds, Novel Superfluids. Oxford Univ. Press, Oxford

Pappas G., Apostolatos T. A., 2012, Phys. Rev. Lett., 108, 231104

Peralta C., Melatos A., Giacobello M., Ooi A., 2006, ApJ, 651, 1079

Prix R., Novak J., Comer G. L., 2005, Phys. Rev. D, 71, 043005

Radhakrishnan V., Manchester R. N., 1969, Nature, 222, 228

Reichley P. E., Downs G. S., 1969, Nature, 222, 229

Ruderman M., 1969, Nature, 223, 597

Ruderman M., 1991, ApJ, 382, 587

Ruderman M., Zhu T., Chen K., 1998, ApJ, 492, 267

Sauls J., 1989, in Ögelman H., van den Heuvel E. P. J., eds, NATO Advanced Science Institutes (ASI) Series C Vol. 262, NATO Advanced Science Institutes, Kluwer Academic/Plenum Press, New York, p. 457

Sedrakian A., 2005, Phys. Rev. D, 71, 083003

Sedrakian A. D., Sedrakian D. M., 1995, ApJ, 447, 305

Sedrakian D. M., Shakhabasian K. M., 1980, Astrofizika, 16, 727

Seveso S., Pizzochero P. M., Haskell B., 2012, MNRAS, 427, 1089

Shapiro S. L., Teukolsky S. A., 1983, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects. John Wiley & Sons, New York

Sidery T., Passamonti A., Andersson N., 2010, MNRAS, 405, 1061

Sourie A., Oertel M., Novak J., 2016, Phys. Rev. D, 93, 083004

Stappers B. W. et al., 2011, A&A, 530, A80

Stopnitzky E., Profumo S., 2014, ApJ, 787, 114

Tsakadze J. S., Tsakadze S. J., 1980, J. Low Temperature Phys., 39, 649

Typel S., Wolter H. H., 1999, Nucl. Phys. A, 656, 331

van Eysden C. A., Melatos A., 2008, Class. Quant. Gravity, 25, 225020

van Eysden C. A., Melatos A., 2010, MNRAS, 409, 1253

Villain L., Haensel P., 2005, A&A, 444, 539

Warszawski L., Melatos A., 2011, MNRAS, 415, 1611

Warszawski L., Melatos A., 2013, MNRAS, 428, 1911

Watts A. et al., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14). Giardini Naxos, Italy, p. 43

Wong T., Backer D. C., Lyne A. G., 2001, ApJ, 548, 447

Yakovlev D. G., Kaminker A. D., Gnedin O. Y., Haensel P., 2001, Phys. Rep., 354, 1

Yarmchuk E. J., Gordon M. J. V., Packard R. E., 1979, Phys. Rev. Lett., 43, 214

Yu M. et al., 2013, MNRAS, 429, 688

Yuan J. P., Manchester R. N., Wang N., Zhou X., Liu Z. Y., Gao Z. F., 2010, ApJ, 719, L111

Zwierlein M. W., Abo-Shaeer J. R., Schirotzek A., Schunck C. H., Ketterle W., 2005, Nature, 435, 1047
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Coupling coefficient $\hat{\epsilon}_p$ and entrainment parameter $\tilde{\epsilon}_p$ as functions of the rotation frequency $f = \Omega/(2\pi)$ for a 1.4 M$_\odot$ neutron star, assuming corotation and beta equilibrium. The two EoSs considered are Newtonian versions of the EoSs used in Section 5.2.2.</figcaption></figure>

## A2 Glitch rise time

Using equations (A1), the moments of inertia defined by (13) lead in the Newtonian limit to

$$ 
\begin{align*}
I_{nn} &= I_n (1 - \tilde{\epsilon}_n) + \Omega_p \frac{\partial I_n}{\partial \Omega_n}, \\
I_{np} &= I_n \tilde{\epsilon}_n + \Omega_p \frac{\partial I_n}{\partial \Omega_p}, \\
I_{pn} &= I_p \tilde{\epsilon}_p + \Omega_p \frac{\partial I_p}{\partial \Omega_p}, \\
I_{pp} &= I_p (1 - \tilde{\epsilon}_p) + \Omega_p \frac{\partial I_p}{\partial \Omega_p},
\end{align*}
\tag{A11} $$

where the partial derivatives have been evaluated for vanishing lag. For slowly rotating stars, $I_n$ and $I_p$ are approximately independent of the rotation rates and the Newtonian moments of inertia should thus verify $\hat{I}_n \simeq I_n$, $\hat{I}_p \simeq I_p$, $\hat{I} \simeq I_n + I_p \equiv I$ and $I_{np} \simeq I_n \tilde{\epsilon}_n = I_p \tilde{\epsilon}_p$. Consequently, comparing (23) and (A10) yields $\zeta = 1$ in Newtonian gravity. As a result, the glitch rise time (equation 24) in the Newtonian limit is given by

$$ \tau_r = \frac{I_p}{I} \times \frac{(1 - \tilde{\epsilon}_p - \tilde{\epsilon}_n)}{2\tilde{B}\Omega_n}, \tag{A12} $$

which coincides with equation (69) from Sidery et al. (2010).

## A3 Fluid couplings

Using equations (A11), the Newtonian proton coupling coefficient $\hat{\epsilon}_p$ (27) is given by

$$ \hat{\epsilon}_p = \frac{\tilde{\epsilon}_p + \frac{\Omega_p}{I_p} \frac{\partial I_p}{\partial \Omega_n}}{1 + \frac{\Omega_p}{I_p} \left( \frac{\partial I_p}{\partial \Omega_p} + \frac{\partial I_p}{\partial \Omega_n} \right)}. \tag{A13} $$

We have computed Newtonian equilibrium configurations as discussed in Section 5.2.2. In Fig. A1, the coefficients $\hat{\epsilon}_p$ and $\tilde{\epsilon}_p$ are plotted with respect to the rotation frequency $f$, for a 1.4 M$_\odot$ neutron star, assuming corotation and beta equilibrium. As expected, at low angular velocities $\hat{\epsilon}_p \simeq \tilde{\epsilon}_p$ to a very good approximation since the moment of inertia $I_p$ is nearly constant.

For frequencies higher than ~100 Hz, the effects of rotation on the stellar structure become non-negligible, and are twofold. First, the proton entrainment parameter $\hat{\epsilon}_p$ is slightly decreased because the central density is lowered. More importantly, the moments of inertia change thus leading to large deviations between $\hat{\epsilon}_p$ and $\tilde{\epsilon}_p$.

## APPENDIX B: CONSTRAINTS ON THE PARTIAL MOMENTS OF INERTIA

From the application of the action principle to the asymptotically flat stationary states of an axisymmetric star composed of two fluids in circular motion with rigid angular velocities, the change in the energy of the star, between two nearby states, is given by

$$ \delta E = \sum_X \Omega_X \delta J_X, \tag{B1} $$

for a fixed total baryon mass (see equation 3.5 from Carter 1975). Using the definitions (equation 13) of the moments of inertia, the variation in the angular momenta simply reads

$$ \delta J_X = \sum_Y I_{XY} \delta \Omega_Y. \tag{B2} $$

Recalling that $I_{XY} = I_{YX}$, equation (B1) is thus given by

$$ \delta E = \sum_{X,Y} I_{XY} \Omega_X \delta \Omega_Y \tag{B3} $$

$$ = \frac{1}{2} \left( I_{nn} \delta (\Omega_n^2) + 2 I_{np} \delta (\Omega_n \Omega_p) + I_{pp} \delta (\Omega_p^2) \right). \tag{B4} $$

Considering very small rotation rates, the energy $E$ of a rotating state is therefore given by

$$ E = E_0 + \frac{1}{2} \left( I_{nn} \Omega_n^2 + 2 I_{np} \Omega_n \Omega_p + I_{pp} \Omega_p^2 \right), \tag{B5} $$

where $E_0$ stands for the energy of the static configuration. Rewriting this equation as follows

$$ E - E_0 = \frac{1}{2} I_{nn} \left( \Omega_n + \frac{I_{np}}{I_{nn}} \Omega_p \right)^2 + \frac{1}{2} \left( I_{pp} - \frac{I_{np}^2}{I_{nn}} \right) \Omega_p^2, \tag{B6} $$

the stability of the static state implies that the right-hand side should be strictly positive, leading to

$$ I_{nn} > 0 \text{ and } I_{nn} I_{pp} - I_{np}^2 > 0, \tag{B7} $$

which in turn gives $I_{pp} > 0$.

## APPENDIX C: RELATIVISTIC COUPLING PARAMETERS

Correcting a typo in equation (A1) of Sourie et al. (2016), the angular momentum of a fluid, say $X$, reads

$$ J_X = \int_{\Sigma_I} \left[ \Gamma_X^2 n_X \mu^X U_X + 2\alpha \frac{\Gamma_X^2}{\Gamma_\Delta^2} \left( \frac{\Gamma_Y}{\Gamma_\Delta \Gamma_X} U_Y - U_X \right) \right] \times \text{Br} \sin\theta \, d^3\Sigma. \tag{C1} $$

In this expression, $n_X$ and $\mu^X$ are, respectively, the particle density and the chemical potential of the fluid, as measured in its rest frame. The norms $U_X$ and $U_Y$ of the physical velocities of the fluids with respect to the ZAMO are given by

$$ U_X = \frac{B}{N} (\Omega_X - \omega) r \sin\theta, \tag{C2} $$
<figure><img src="image_13.png" /><figcaption>**Figure C1.** Left: proton fraction $x_p$ as a function of the total baryon density $n_B = n_n + n_p$, assuming corotation and beta equilibrium. Right: ratio of the proton moment of inertia $\hat{I}_p$ to the total one $\hat{I}$ with respect to the gravitational mass for a star spinning at 11.19 Hz, assuming beta equilibrium at the centre. In both figures, results obtained from the DDH($\delta$) EoS are plotted with red solid (green dashed) lines.</figcaption></figure>

where B, N and $\omega$ are different potentials involved in the space-time metric (see Sourie et al. 2016). The volume element is given by $d^3\Sigma = A^2Br^2 \sin\theta dr d\theta d\varphi$. $\Gamma_X$, $\Gamma_Y$ and $\Gamma_\Delta$ are Lorentz factors associated with $U_X$, $U_Y$ and $\Delta$, the relative speed between the fluids. The quantity $\alpha$, which characterizes the presence of entrainment, is linked to the more common entrainment parameter $\epsilon_X$ through

$$ \epsilon_X = \frac{2\alpha}{n_X \mu^X \Gamma_\Delta^2} \qquad (C3) $$

(see section III-B of Sourie et al. 2016).

In the slow-rotation approximation ($\Omega_n$, $\Omega_p \ll \Omega_K$) and to the first order in the lag $\delta\Omega = \Omega_n - \Omega_p$, equation (C1) becomes

$$ J_X \approx \int_{\Sigma_i} n_X \mu^X \frac{B^2}{N} r^2 \sin^2 \theta (\Omega_X - \omega) d^3\Sigma + \int_{\Sigma_i} n_X \mu^X \frac{B^2}{N} r^2 \sin^2 \theta \epsilon_X (\Omega_Y - \Omega_X) d^3\Sigma, \qquad (C4) $$

where the couplings by entrainment and Lense-Thirring effects are clearly visible, respectively through the terms $\epsilon_X$ and $\omega$. We now introduce the following quantities

$$ \tilde{I}_X = \int_{\Sigma_i} n_X \mu^X \frac{B^2}{N} r^2 \sin^2 \theta d^3\Sigma, \qquad (C5) $$

$$ \tilde{\epsilon}_X \tilde{I}_X = \int_{\Sigma_i} n_X \mu^X \frac{B^2}{N} r^2 \sin^2 \theta \epsilon_X d^3\Sigma, \qquad (C6) $$

and

$$ \tilde{\omega}_X \tilde{I}_X = \int_{\Sigma_i} n_X \mu^X \frac{B^2}{N} r^2 \sin^2 \theta \omega d^3\Sigma, \qquad (C7) $$

such that equation (C4) now reads

$$ J_X = \tilde{I}_X (\Omega_X - \tilde{\omega}_X) + \tilde{\epsilon}_X \tilde{I}_X (\Omega_Y - \Omega_X). \qquad (C8) $$

In the Newtonian limit, $\tilde{I}_X$ and $\tilde{\epsilon}_X$ are, respectively, given by equations (A2) and (A3) and $\tilde{\omega}_X$ simply vanishes, so that (C8) tends

towards (A1). We numerically find that $\tilde{\omega}_X$ can be approximated by a relation of the form

$$ \tilde{\omega}_X = \epsilon_{XX}^{LT} \Omega_X + \epsilon_{YX}^{LT} \Omega_Y, \qquad (C9) $$

with a precision better than 0.1 per cent for a star spinning at 65 Hz or less. In this equation, $\epsilon_{YX}^{LT}$ represents the frame-dragging contribution of fluid Y on fluid X and $\epsilon_{XX}^{LT}$ denotes the self-frame-dragging effect of fluid X. By making use of the different coupling parameters introduced so far, the angular momentum of fluid X is given by

$$ J_X = \tilde{I}_X (1 - \epsilon_{XX}^{LT} - \tilde{\epsilon}_X) \Omega_X + \tilde{I}_X (\tilde{\epsilon}_X - \epsilon_{YX}^{LT}) \Omega_Y. \qquad (C10) $$

To this level of approximation, the coupling parameter $\hat{\epsilon}_X$ (27) reads

$$ \hat{\epsilon}_X = \frac{\tilde{\epsilon}_X - \epsilon_{YX}^{LT}}{1 - \epsilon_{YX}^{LT} - \epsilon_{XX}^{LT}}. \qquad (C11) $$

## APPENDIX D: REALISTIC PROTON FRACTIONS

The proton fractions are plotted in left panel of Fig. C1 as functions of the total baryon density, for the DDH and DDH$\delta$ EoSs. For both EoSs, $x_p$ is strongly increasing with the density, until $n_B \gtrsim 1.1 \text{ fm}^{-3}$. The proton fraction profiles are very different between both EoSs because of the values considered for the symmetry energy and its slope at saturation density (see table I of Sourie et al. 2016). Consequently, as the mass of the stars increases, higher values of the proton fraction are reached and the quantity $\hat{I}_p/\hat{I}$ increases, as can be seen in the right panel of Fig. C1. This ratio is higher for the DDH EoS because, for densities found in neutron stars, the DDH proton fraction is more important than the DDH$\delta$ one and the central baryon densities are systematically higher in the case of DDH.

This paper has been typeset from a TeX/Latex file prepared by the author.