# Insights into the physics of neutron star interiors from pulsar glitches

M. Antonelli1, A. Montoli3, P. Pizzochero3,4 1 CNRS/IN2P3, Laboratoire de
Physique Corpusculaire, 14050 Caen, France
2 Dipartimento di Fisica, Università degli Studi di Milano, via Celora 16,
20133 Milano, Italy
3 Istituto Nazionale di Fisica Nucleare, via Celoria 16, 20133 Milano, Italy

###### Abstract

The presence of superfluid phases in the interior of a neutron star affects
its dynamics, as neutrons can flow relative to the non-superfluid (normal)
components of the star with little or no viscosity. A probe of superfluidity
comes from pulsar glitches, sudden jumps in the observed rotational period of
radio pulsars. Most models of glitches build on the idea that a superfluid
component of the star is decoupled from the spin-down of the normal component,
and its sudden re-coupling leads to a glitch. This transition in the strength
of the hydrodynamic coupling is explained in terms of quantum vortices (long-
lived vortices that are naturally present in the neutron superfluid at the
microscopic scale). After introducing some basic ideas, we derive (as a
pedagogical exercise) the formal scheme shared by many glitch studies. Then,
we apply these notions to present some recent advances and discuss how
observations can help us to indirectly probe the internal physics of neutron
stars.

##  I Pulsar glitches: general concepts

The present work111 Chapter of the volume _Astrophysics in the XXI Century
with Compact Stars_, Eds. C.A.Z. Vasconcellos and F. Weber, World Scientific
(2022) Vasconcellos and Weber (2022), submitted in August 2021. Two recent
reviews on pulsar glitches, complementary to this one, are Zhou  _et al._
(2022); Antonopoulou  _et al._ (2022).  is an introduction to the current
interpretation and modelling of the glitch phenomenon, a sudden change in the
stable rotation period of pulsars (complementary reviews on the same subject
are (D’Alessandro, 1996; Haskell and Melatos, 2015), see also (Manchester,
2018) for a historical introduction). The field is far from being settled, as
most of the results are qualitative and based on (Newtonian) phenomenological
models. Nonetheless, the glitch phenomenon is probably the most clear
observational evidence for the existence of an extended (kilometre-sized)
superfluid region in the inner layers of a neutron star (NS). Nucleon
superfluidity in the core of NSs has found additional support from the real
time monitoring of the rapid cooling of the young NS in Cassiopeia A (Page
_et al._ , 2011; Shternin  _et al._ , 2011), but this observational result and
its interpretation still have to be firmly assessed (Posselt and Pavlov, 2018;
Wijngaarden  _et al._ , 2019).

Pulsar glitches are sudden increases in the rotation frequency of a spinning
down neutron star (Lyne  _et al._ , 2000; Espinoza  _et al._ , 2011; Fuentes
_et al._ , 2017; Manchester, 2018). Soon after the first glitch observations
in the Vela pulsar (Radhakrishnan and Manchester, 1969; Reichley and Downs,
1969), Baym  _et al._ (1969) associated the long post-glitch relaxation
timescales of order days to months to the presence of a superfluid component
in the internal layers.

The existence of superfluids in NSs was conjectured even before the discovery
of pulsars. In fact, neutron superfluidity is theoretically expected in NS
interiors as most of the star will be cold enough for neutrons to undergo
Cooper pairing (Sauls, 1989; Sedrakian and Clark, 2019). Depending on the age
of the NS, the interior temperatures are in the range $10^{7}-10^{9}\,$K,
which are low on the scale defined by the typical Fermi energies of about
$10-100\,$MeV (equivalent to $10^{11}-10^{12}\,$K), so that the ratio of
temperature to Fermi energy in a NS is comparable to that of laboratory
superfluid 3He in the millikelvin temperature range.

For what concerns the glitch phenomenon, the presence of a neutron superfluid
significantly modifies the internal dynamics of a NS, allowing for the
superfluid to flow relative to the “normal” (i.e., non superfluid) components.
This relative flow determines the angular momentum reservoir needed to explain
the observed jumps of the rotation period.

###  I.1 An analogy with type-II superconductors: pining and creep

To date, the most promising explanation of the glitch mechanism is based on
the seminal idea of Anderson and Itoh (1975) and is built on the analogy
between the superfluid current of dripped neutrons in the inner crust of a NS
and the electron current in a type-II superconductor (Anderson and Kim, 1964).
In fact, a type-II superconductor immersed in a magnetic field can sustain
electron currents in a non-dissipative state as long as the vortices carrying
a quantum of magnetic flux are anchored to the crystalline structure of the
specimen – i.e., the quantized flux-tubes are “pinned”. On the contrary, the
superconductor enters a resistive state when the flux-tubes unpin and start to
move (Blatter  _et al._ , 1994).

Similarly, in a NS a current of neutrons can flow without resistance as long
as the quantized vortex lines, that are naturally present in a rotating
superfluid, find energetically favorable to pin to some impurity in the
medium. In the non-homogeneous environment of the inner crust the pinning
impurities may be the ions that constitute the crustal lattice (Alpar, 1977;
Donati and Pizzochero, 2004; Bulgac  _et al._ , 2013; Seveso  _et al._ ,
2016). If the proton fluid in the outer core forms a type-II superconductor
(Baym  _et al._ , 1969), pinning of neutron vortices to the quantized flux-
tubes may be a viable possibility as well, see (Alpar, 2017; Sourie and
Chamel, 2021) and and references therein. However, the pinning in the outer
core is, possibly, even less studied than the already difficult issue of
pinning in the inner crust. In fact, the nature of the proton superconductor
may be more complex than that of a usual type-II (Haber and Schmitt, 2017;
Wood  _et al._ , 2020), and also the features of the vortices threading the
spin-triplet neutron superfluid expected in most of the outer core can differ
from the simpler spin-singlet Cooper pairing case realized in the crust (Brand
and Pleiner, 1981; Leinson, 2020).

Only when the quantized vortices unpin, the neutron superfluid and the rest of
the star can interact at the hydrodynamic level via a dissipative force (the
_vortex-mediated mutual friction_), which tends to drive the system towards a
zero current equilibrium state.

To complete the above qualitative picture, it is necessary to provide a
reasonable argument for how a relative current of neutrons should develop in
the first place. The seminal Anderson’s idea is based on the fact that a
rotating superfluid can spin down only by expelling part of its quantized
vorticity. However, if the vortex lines are somehow pinned to the nuclear
clusters in the crust, their natural outward motion is hindered and the
superfluid cannot exactly follow the spin-down of the crust (Pines and Alpar,
1985). In this way, the superfluid component naturally tends to rotate faster
than the spinning-down normal component, storing an excess of angular
momentum. This reservoir of momentum can then be suddenly released during a
glitch when several vortex lines ($\sim 10^{12}$ out of $\sim 10^{18}$ in a
typical pulsar) are expelled from the superfluid domain222 More precisely, the
unpinned vortices will tend to move outward – away from the NS’s rotation axis
– but their displacement is not necessarily uniform across the superfluid
region (see e.g. Alpar  _et al._ , 1984a; Antonelli and Pizzochero, 2017;
Howitt  _et al._ , 2020). . Clearly, the same set of ideas applies also in
reverse to explain anti-glitches in spinning-up pulsars, namely sudden spin-up
events during which the creep motion of vortices is directed inward (Ducci
_et al._ , 2015; Ray  _et al._ , 2019).

###  I.2 Collective vortex dynamics

The glitch phenomenon is probably the result of a competition between pinning
and the hydrodynamic lift – the so-called _Magnus force_ – felt by pinned
vortices because of the background flow of neutrons. Interestingly, vortex
avalanches occur routinely in two dimensional simulations of many vortices
(Howitt  _et al._ , 2020), when chains of unpinning events are triggered
collectively by vortex-vortex repulsion (Warszawski  _et al._ , 2012), in
accordance with smaller scale studies based on the Gross-Pitaevskii equation
(Warszawski and Melatos, 2011).

This process is sketched in Fig. 1: in principle, we can imagine to associate
a certain energy to each vortex configuration in the inner crust, by
accounting for the vortex-vortex and the vortex-lattice interactions. This
would generate an intricate energy landscape in the space of all possible
configurations of vortices (reduced to only two dimensions in the figure). If
no velocity lag between the components is present (as in the left panel of
Fig. 1), then each configuration relaxes to the local minimum by rearranging
the mutual positions and shape of the vortices. However, the external torque
slowly builds a velocity lag, which tends to tilt the energy landscape.
Whether this produces a continuous rearrangement of the vortex configuration,
or a sudden slip to a new minimum (as in the central panel) is still not clear
and it really depends on the properties of the landscape. The idea that
glitches may the byproduct of a self-organized criticality of the vortex
configuration (Morley and Garcia-Pelayo, 1993; Melatos and Peralta, 2007) is
based on the latter possibility. Whether or not this process would be a
manifestation of self-organization in the strict sense (Watkins  _et al._ ,
2016), Anderson and Itoh (1975) described this driven relaxation process as a
_noisy vorticity creep_ process, which may be thought of as a sequence of
“quakes” in the vortex configuration.

Despite being just a cartoon, this abstract point of view on the glitch
phenomenon allows one to recognize that the sequence of dynamical phases
(explored in Sec. IV) experienced by a pulsar during its slowly driven
evolution depends on the properties of an intricate energy landscape, so that
the observed probability distributions of the avalanche sizes and waiting
times are the manifestation of a stress-release model with multiple thresholds
(Carlin and Melatos, 2021). A fast driven system (for which the tilt of the
potential tends to be steep) is less likely to develop such intermittent fast-
slow dynamics; on the contrary a smoother and continuous creep of vortices has
to be expected (Alpar  _et al._ , 1984a; Baym  _et al._ , 1992; Antonelli and
Haskell, 2020).

Figure 1:  Cartoon of the relaxation dynamics between different metastable
vortex configurations. A particular configuration of vortex lines, here
represented by the blue (or green) point, is immersed in the space of all
possible configurations. The total energy associated to a configuration
defines a disordered landscape over the configuration space. The external
torque drives the system and continuously tilts such landscape. When a
critical tilt is reached, a metastable configuration can relax to a new
minimum (during the relaxation process the tilt angle decreases as the lag
between the components is reduced). Since the energy landscape is extremely
irregular, different configurations (like the green and the blue points) will
start to relax at different tilt angles, following different paths. The final
configuration depends on how quickly the tilt angle decreases during the
process and on the initial configuration itself. In principle, the same idea
applies in reverse for anti-glitches in spinning-up pulsars.

The nature of the trigger for vortex unpinning is still unclear and proposals
range from vortex avalanches (Anderson and Itoh, 1975; Pines, 1999; Melatos
and Peralta, 2007) to hydrodynamical instabilities (Andersson  _et al._ ,
2003; Glampedakis and Andersson, 2009; Khomenko  _et al._ , 2019) or crust
quakes (Ruderman, 1969; Alpar  _et al._ , 1994; Giliberti  _et al._ , 2020).
Recently, during the observation of the 2016 Vela glitch, changes in the
pulsar magnetosphere have been observed333 An even more dramatic change was
observed in the high-magnetic field pulsar J1119-6127, where the pulse profile
switched from single to double following large glitch (Weltevrede  _et al._ ,
2011).  (Palfreyman  _et al._ , 2018). This may be a consequence of one of the
trigger processes listed above, in particular a quake (Bransgrove  _et al._ ,
2020), or it may even be a clue of the involvement of the magnetic field
itself during the glitch phenomenon (Ruderman  _et al._ , 1998). This
interpretation seems to be also supported by the recent analysis of Ashton
_et al._ (2020), that reveals an increased flickering activity of the
magnetosphere around the glitch epoch. This result, if confirmed, is
interesting, as the flickering activity may be linked to precursors and
aftershocks, much like in terrestrial earthquakes. If this is the case, then
in Fig. 1 the relaxation to a new metastable vortex configuration can be
triggered before the local critical tilt of the energy landscape is reached,
because the quake provides an initial kick to initiate the process.

Whatever the trigger mechanism, one of the key ingredients in the above
picture is the pinning force that the crustal lattice can exert on a vortex.
Similarly to what happens in type-II superconductors where pinning enhances
the critical current value, this quantity determines the maximum amount of
angular momentum that could be exchanged during a glitch (Antonelli  _et al._
, 2018). The understanding of how much angular momentum can be stored in
different regions of the star would, in fact, allow a comparison with
observations of glitching pulsars, potentially constraining the internal
structure of a NS (Ho  _et al._ , 2015; Pizzochero  _et al._ , 2017) when used
in tandem with other methods based on the study of the average glitch activity
(Andersson  _et al._ , 2012; Chamel, 2013; Carreau  _et al._ , 2019; Montoli
_et al._ , 2021), dynamical fits of the spin-up phase (Graber  _et al._ ,
2018; Ashton  _et al._ , 2019; Montoli  _et al._ , 2020a) or the glitching
history of pulsars (Montoli  _et al._ , 2020b; Gügercinoğlu and Alpar, 2020).

In the following we review some notions on glitch modelling and then move to
discuss some recent advances in the field. In the last years a number of
studies have been devoted to demonstrating that the analysis of glitches
across the pulsar population has the potential to put constraints on the
dense-matter properties. These analyses are only first steps, as more high-
cadence pulsar timing observations are necessary before it will be possible to
falsify some of the theoretical predictions and test the microscopic inputs.

##  II Basic formalism for glitch modelling

As we have seen, glitch modelling is based on several fruitful analogies
between spinning neutron stars and superfluid, or superconducting, systems
that are studied in the laboratory (see also Graber  _et al._ (2017)). In
fact, the tools and ideas that are currently used to describe the rotational
dynamics of a NS are the byproducts of an attempt to adapt and extend the
Tisza-Landau two-fluid model for superfluid Helium to account for the several
fluid species, and their mutual interactions, present in the interior of a NS,
as it has recently been reviewed in (Haskell and Sedrakian, 2017; Chamel,
2017a; Andersson, 2021).

###  II.1 Phenomenological model with two rigid components

The pulsar rotation frequency $\nu(t)=1/P(t)$ can be measured at a given
reference time $t$ by means of pulsar timing techniques, namely the continuous
recording of the time of arrival of each pulse at the telescope. The precise
timing of rotation-powered pulsars reveals a steady and extremely slow
increase of the pulse period, indicating that these objects lose angular
momentum and kinetic energy due to some emission mechanism. However, the
rotational frequency of a pulsar occasionally undergoes sudden jumps of
amplitude $\Delta\nu>0$, which are followed by a period of slow recovery that
may last for days or months, the so-called glitches. The signature of a glitch
in the timing residuals is rather clear for relative jumps bigger than
$\Delta\nu/\nu\sim 10^{-7}$: when a glitch occurs, the model $P(t)$ for the
slow period evolution before the glitch can no longer predict the time of
arrivals, which have to be described by a different timing solution by
including a discontinuity in $\nu$ and $\dot{\nu}$ around the inferred epoch
of the glitch.

Since the star is an isolated object, the observed sudden spin-up can be
explained only in terms of a sudden change of the NS’s moment of inertia, or
by invoking the action of an internal torque (or both). In particular, the
subsequent slow post-glitch relaxation of $\nu(t)$ observed in several pulsars
is an indication that an internal torque is at work. In other words, there is
a loose component inside the star whose rotation is not directly observable,
but that is exchanging angular momentum with the observable component.

_A minimal linear model -_ Baym  _et al._ (1969) proposed what is probably the
simplest phenomenological description of a glitch: the observable component
(assumed to be rigidly coupled to the magnetosphere, where the pulsar signal
originates) and the loose component in the interior are both modelled as rigid
and non-precessing bodies rotating around a common axis and coupled by a
linear viscous-like internal torque (see, e.g., Gügercinoǧlu and Alpar (2017)
and Celora  _et al._ (2020) for recent nonlinear extensions),

$\displaystyle
I_{o}\dot{\Omega}_{o}+I_{s}\dot{\Omega}_{s}=-I|\dot{\Omega}_{\infty}|$ (1)
$\displaystyle
I_{s}\dot{\Omega}_{s}=-\frac{I_{s}\,I_{o}}{I\,\tau_{exp}}(\Omega_{s}-\Omega_{o})=-\frac{I_{s}}{\tau}(\Omega_{s}-\Omega_{o})$
(2)

where ${\Omega}_{s}(t)$ and ${\Omega}_{o}(t)=2\pi\nu(t)$ are the angular
velocities of the loose component (identified with the superfluid in an
unknown internal region) and of the observable component, respectively. The
other parameters are the moments of inertia of the two components, $I_{s}$ and
$I_{o}$, while $|\dot{\Omega}_{\infty}|>0$ sets the intensity of the external
spin-down torque, usually linked to electromagnetic radiation (its exact
nature is unimportant here, but it is assumed to be constant on the timescale
of interest). The total moment of inertia is $I=I_{s}+I_{o}$. Since a real NS
is not uniform and the rotation of the loose component may be differential
(its rotation rate may depend on the position inside the star), the quantities
that appear in the above equations have to be understood as suitable spatial
averages over some region in the interior. Finally, there are two common ways
to parametrize the internal torque: in terms of the observed post-glitch
exponential timescale $\tau_{exp}$ (the _relaxation_ timescale), or in terms
of the _coupling_ timescale $\tau$, that is more directly linked to the
microscopic processes responsible for the interaction between the two
components.

It is interesting to notice that the first equation (1) – representing the
total angular momentum balance of the star – is a fundamental statement, while
the modelling enters into the second equation (2). The latter, in fact,
includes the interaction between the loose and the observable components which
is responsible for the angular momentum exchange between the two components.
For a given angular velocity lag $\Omega_{so}^{0}$ at $t=0$ so that
$\Omega_{o}(0)=\Omega_{o}^{0}$ and
$\Omega_{s}(0)=\Omega_{o}^{0}+\Omega_{so}^{0}$, the solution to the system is

$\displaystyle\Omega_{o}(t)$
$\displaystyle=\Omega^{\infty}_{o}(t)+\Delta\Omega_{o}(t)$ (3)
$\displaystyle\Omega_{s}(t)$
$\displaystyle=\Omega^{\infty}_{o}(t)+\Omega_{so}^{0}-(I_{o}/I_{s})\,\Delta\Omega_{o}(t)\,.$
(4)

where $\Omega^{\infty}_{o}(t)=\Omega_{o}^{0}-t|\dot{\Omega}_{\infty}|$
describes the steady spin-down, while the angular velocity residue
$\Delta\Omega_{o}(t)$ reads

$\Delta\Omega_{o}(t)=\frac{I_{s}}{I}(\Omega_{so}^{0}-\Omega_{so}^{\infty})\big{(}1-e^{-t/\tau_{exp}}\big{)}\,.$
(5)

Here $\Omega_{so}^{\infty}=\tau|\dot{\Omega}_{\infty}|$ is the constant value
of the lag in the steady state, that is the asymptotic situation in which both
components spin-down at the same rate $|\dot{\Omega}_{\infty}|$.

The general solution (3) may be used to fit the observed post-glitch
evolution: assuming that the steady-state spin-down rate
$|\dot{\Omega}_{\infty}|$ and the angular velocity $\Omega_{o}^{0}$ are known
(so that $\Omega^{\infty}_{o}(t)$ is given), we just have to fit the residue,
an operation that would require the fitting of 2 parameters (the timescale
$\tau_{exp}$ and the overall amplitude) that are related to the physical
properties of the NS444 If the model is used to fit the post-glitch
relaxation, then the initial conditions $\Omega_{so}^{0}$ and $\Omega_{o}^{0}$
refer to some time $t=0$ after the unresolved spin-up. From the practical
point of view, observations may not be frequent enough to pin the exact
instant of the spin-up event. .

###  II.2 Spin-up due to a starquake

In their original work, Baym and collaborators assumed that the unresolved
spin-up is due to a change in the moments of inertia after a starquake,
following the abrupt release of part of the stress accumulated in the crustal
solid during the whole NS’s life. The solid crust formed when the rotation
rate was higher, so that its subsequent relaxation towards a less oblate shape
(because of the spin-down) may be realized via a sequence of failures. Within
this picture, the loose component, whose location and extension within the
star is unknown, only serves the purpose to explain why the slow post-glitch
relaxation proceeds over the timescale $\tau_{exp}$.

We can use the linear system in (1) and (2) to study the relaxation after a
spin-up induced by a starquake (see also Pines  _et al._ , 1974; Shapiro and
Teukolsky, 1983). We already have the general solution, so we just have to
link the initial condition $\Omega_{so}^{0}$ with the angular velocity jump
due to a change in the moments of inertia. The starquake can be modelled as an
instantaneous (i.e. unresolved) global change at $t=0$, such that

$\begin{split}I_{s,o}\,&\rightarrow\,(1-\epsilon_{s,o})I_{s,o}\\\
\Omega^{0}_{s,o}\,&\rightarrow\,\Omega^{0}_{s,o}/(1-\epsilon_{s,o})\\\
|\dot{\Omega}_{\infty}|\,&\rightarrow\,(1+\epsilon_{\infty})|\dot{\Omega}_{\infty}|\,,\\\
\tau\,&\rightarrow\,(1+\epsilon_{\tau})\tau\end{split}$ (6)

where $\epsilon_{s,o}>0$ since the star should evolve towards a less oblate
shape, while the sign of $\epsilon_{\tau,\infty}$ is undetermined. Assuming
the star to be in the steady state before the starquake, the post-glitch
relaxation can be written in terms of the pre-glitch parameters and the
relative variations $\epsilon$,

$\Omega_{o}(t)=(1+\epsilon_{o})\Omega_{o}^{0}-(1+\epsilon_{\infty})|\dot{\Omega}_{\infty}|t+\frac{I_{s}}{I}\big{[}(\epsilon_{s}-\epsilon_{o})\Omega_{o}^{0}-(\epsilon_{\infty}+\epsilon_{\tau}-\epsilon_{s})\Omega_{so}^{\infty}\big{]}\big{(}1-e^{-t/\tau_{exp}}\big{)}+O(\epsilon^{2})$
(7)

where, again, $\Omega_{so}^{\infty}=|\dot{\Omega}_{\infty}|\tau$ is the pre-
starquake steady-state lag (all the quantities but $\epsilon$ in (7) refer to
some instant before the starquake, including $\Omega_{o}^{0}$).

The observed glitch jump $\Delta\nu$ allows us to fix the parameter
$\epsilon_{o}$ as $\epsilon_{o}=\Delta\nu/\nu$, but it is difficult to
estimate the other $\epsilon$ parameters. However, Baym and collaborators
noticed that it is possible to relate $\epsilon_{s}$ to the so-called healing
parameter $Q$, which may be obtained from observations: to the first order in
$\epsilon$, equation (7) may be expressed as

$\Omega_{o}(t)=\Omega_{o}^{0}+\Delta\Omega_{o}-(1+\epsilon_{\infty})|\dot{\Omega}_{\infty}|t-Q\Delta\Omega_{o}\big{(}1-e^{-t/\tau_{exp}}\big{)}$
(8)

where $\Delta\Omega_{o}\approx\epsilon_{o}\Omega_{o}^{0}$ is the observed
glitch amplitude. Hence,

$Q=\frac{I_{s}}{I}\left(1-\frac{\epsilon_{s}}{\epsilon_{o}}+\frac{\epsilon_{\infty}+\epsilon_{\tau}-\epsilon_{s}}{\epsilon_{o}}\,\frac{\Omega_{so}^{\infty}}{\Omega_{o}^{0}}\right)\approx\frac{I_{s}}{I}\left(1-\frac{\epsilon_{s}}{\epsilon_{o}}\right)\,>\,\frac{I_{s}}{I}$
(9)

where the approximated expression is obtained in the limit
$\Omega_{so}^{\infty}\ll\Omega_{o}^{0}$. Crawford and Demiański (2003)
collected the all the measured values of the healing parameter $Q$ for Crab
and Vela pulsars available in the literature at that time. They found that a
weighted mean of the $\sim 20$ measured values of $Q$ for Crab glitches the
values yields $Q\approx 0.7\pm 0.05$. On the other hand, a mean value for the
glitches in the Vela pulsar yields the smaller value $Q\approx 0.12\pm 0.07$.
To test the starquake model for the Crab and the Vela pulsars, they also
estimated the healing parameter as $Q\approx I_{s}/I$ – which is valid in the
limit of $I_{o}\ll I_{s}$ and/or $\epsilon_{s}I_{s}\ll\epsilon_{o}I_{o}$ – for
different equations of state dense nuclear matter covering a range of neutron
star masses, and interpreting $I_{s}$ as the moment of inertia of the neutrons
in the core. They conclude that the results are consistent with the starquake
model predictions for the Crab pulsar, but the much smaller values of $Q<0.2$
of the Vela pulsar are inconsistent with the constraint in (9), since the
implied Vela mass would be unreasonably small, about $\sim 0.5\,M_{\odot}$.

###  II.3 Some insights from seismology

Despite the starquake model is a plausible model for the glitches in the Crab
pulsar, although some difficulty arises in the case of the Vela. Apart from
the large amplitudes of the Vela glitches, that would require a consistent
change in the moments of inertia of the order of
$\delta\nu/\nu\approx\epsilon_{o}\sim 10^{-6}$, the problem is also due to the
frequency of its large glitches (roughly one in $2-3$ years): the spin-down is
too slow for the crust to break because of the change in the centrifugal force
between two glitches. Despite the situation is less severe for the Crab, the
important point is that starquakes may still be a viable trigger of large
Vela-like glitches (Akbal  _et al._ , 2017), but only if (see, e.g., the
discussion in (Giliberti  _et al._ , 2020)):

  1. -

the crust never relieves all the accumulated stresses (which is also the case
of earthquakes, see Fig. 2),

  2. -

and/or some internal loading force (possibly due to vortex pinning (Ruderman,
1976) or to the magnetic field (Ruderman  _et al._ , 1998)) that is orders of
magnitude stronger than the variation in the centrifugal force during the
inter-glitch time is operating in the crust.

Despite the fact that crustal motion after failure is also invoked as a
plausible mechanism to explain magnetar emission, it is possible that neutron
star matter cannot exhibit brittle fracture at any temperature or magnetic
field strength (Jones, 2003). Hence, the crust of pulsars can deform in a
gradual plastic manner in spite of its extreme strength, leading to rarer
starquakes than calculated for a rigid crust (Smoluchowski and Welch, 1970).
On the other hand, the presence of defects and of crystal domains (Ruderman,
1991) and the fact that stresses are applied continuously for a long time
(Chugunov and Horowitz, 2010) are all elements that point in the opposite
direction, namely toward the possibility of more frequent crust failures. The
situation is far from being settled and the very question of whether or not
crustquakes actually can happen (and with what frequency) in radio pulsars is
uncertain.

However, the link between quakes and glitches goes beyond the questions about
the nature of the trigger and the magnitude of the change in the moment of
inertia in (6). In fact, an earthquake is triggered when a mechanical
instability occurs and a fracture (the sudden slip of a fault) occurs in the
Earth’s crust. Likewise, a glitch is is triggered by the same mechanism – in
the case the trigger is a crustquake – or by a mechanical instability in the
vortex configuration.

Both glitches and earthquakes are the byproduct of a slowly driven system that
sporadically releases some sort of internal stress in a sequence of short
bursts. This is the kind of dynamics expected to emerge in many complex
systems, that have the tendency to exhibit rare and sudden transitions
occurring over intervals that are short compared with the timescales of their
prior evolution (Jensen, 1998). Such extreme events are expressions of the
underlying forces (like the competing Magnus and pinning forces), or stresses,
usually hidden by an almost perfect balance.

Despite there is an opinion that earthquakes – and so glitches, if the analogy
between the two phenomena holds – could be inherently unpredictable,
observations have shown that the process of seismogenesis is not completely
random. The fact that earthquakes seem to be clustered in time more than would
be expected for a random process (namely, the Omori’s law, which describes the
increased rate of small earthquakes after large earthquakes) is considered to
be expression of a degree of determinism and predictability in the properties
of the earthquake population (Sornette, 2002). In particular, Shimazaki and
Nakata (1980) identified the three idealized scenarios shown in Fig. 2:

  1. _Time-predictable -_ Assume that the stress $\sigma(t)$ along a fault builds up at a nearly constant rate $\dot{\sigma}$ and that a certain amount of stress – resulting in a slip of the fault – is released when a certain (constant in time) critical threshold $\sigma^{c}$ is topped. If we observe a drop in the stress $\Delta\sigma$, then the time to the next quake will be $\Delta t\sim\Delta\sigma/\dot{\sigma}$. Even if $\dot{\sigma}$ is not directly observable, it can be estimated by observation of the previous waiting times between quakes and the amplitudes.

  2. _Slip-predictable -_ The critical threshold $\sigma^{c}$ may not be constant in time: the fault itself is spatially non-homogeneous and it may also rearrange in time (because of small rearrangements during the continuous stress build up or modifications induced by past quakes). However, if the stress in a quake is released till it drops below a certain level $\sigma_{low}$ that is fairly constant in time, then it is possible to predict the size of the quake: $\Delta\sigma\sim\dot{\sigma}\,\Delta t$, where now $\Delta t$ is the observed time elapsed since the preceding quake.

  3. _Periodic -_ When both thresholds $\sigma^{c}$ and $\sigma^{low}$ remain almost constant in time, then the stress release works a clock, so that the distribution of both waiting times between consecutive quakes and sizes should be peaked around a well defined value.

What is important is that observations of a sequence of events may reveal a
tendency of the system for one of those scenarios (or for the most general
one, in which neither $\sigma^{c}$ or $\sigma^{low}$ is constant), as shown in
the right panel of Fig. 2.

Glitch forecast may not be as socially important as earthquake forecast, but,
given the fact that observing glitches requires a constant monitoring effort,
it is interesting to try to give at least rough estimates of the time till the
next glitch (Itoh, 1983; Middleditch  _et al._ , 2006; Akbal  _et al._ ,
2017). To date, a statistical procedure for predicting a glitch in a given
pulsar has not yet been developed, as we still have a lot to learn about the
statistics of glitch occurrence in different objects.

Figure 2:  The three idealized scenarios (left panel) of the temporal
evolution of tectonic stresses (Shimazaki and Nakata, 1980): _periodic_ (where
the initial $\sigma^{c}$ and final $\sigma^{low}$ stresses of faulting are
constant), _time-predictable_ (only $\sigma^{c}$ is constant, so the larger a
quake the longer the following quiet period), _slip-predictable_ (the final
stress $\sigma^{low}$ is constant, so the longer the quiet period the larger
the quake that follows). The same scheme can be applied to the stress
accumulated in the crust of a NS and/or in the vortex configuration because of
the competition between the Magnus and pinning forces. The right panel shows
the cumulative observed slip (i.e., the cumulative quake size, or the
cumulative glitch amplitude) as a function of time for the three scenarios.

While statistical analysis of the whole glitch population has been carried out
by aggregating data from different pulsars (Shemar and Lyne, 1996; Lyne  _et
al._ , 2000; Fuentes  _et al._ , 2017), studies of the glitch statistics for
individual pulsars are considerably more difficult because of small number
statistics. Nonetheless, as the number of detected glitches increases, it is
becoming possible to have a rough idea of the probability density functions
for the glitch sizes and the waiting times in single pulsars (Ashton  _et al._
, 2017; Howitt  _et al._ , 2018; Fuentes  _et al._ , 2019). Most pulsars are
roughly compatible with exponentially distributed waiting times and power-law
distributed event sizes, which is the expected behaviour arising from self-
organized critical systems555 The concept of self-organized criticality is
believed to provide a natural explanation of the statistics of earthquakes,
including the Gutenberg-Richter law for the distribution of earthquake
magnitudes (Bak and Tang, 1989; Sornette and Sornette, 1989; Ito and
Matsuzaki, 1990; Bak  _et al._ , 1994).  (Melatos  _et al._ , 2008). However,
the statistical properties of the Vela pulsar and J0537-6910 are unique among
glitching pulsars, as they typically have large glitches which occur fairly
periodically. Their observed activity is shown in Fig. 3.

Figure 3:  The glitch activity of J0537-6910 and the Vela pulsar (J0835-6910).
The upper and lower points of the cumulative glitch amplitude have been fitted
with the standard linear regression procedure to test whether those two quasi-
periodic pulsars display a tendency for the slip-predictable or time-
predictable scenarios. The residuals of the fits are shown in the right
panels. For J0537-6910 the residuals relative to the slip-predictable fit line
are more dispersed than the ones relative to the time-predictable case,
indicating some slight preference of the system for the latter.

In particular, J0537-6910 is one of the most actively glitching pulsars we
know of, probably because of its very high spin-down rate (see (Ho  _et al._ ,
2020) and (Akbal  _et al._ , 2021) for a recent summary of its properties).
Its glitching behaviour seems to be particularly favourable to glitch
prediction due to the presence of some degree of correlation between the
waiting time from one glitch to the next and the amplitude of the first glitch
(Middleditch  _et al._ , 2006; Ferdman  _et al._ , 2018; Antonopoulou  _et
al._ , 2018). The analysis of its cumulative glitch amplitude in Fig. 3 also
reveals an imprint of such a time-predictable behaviour.

Apart from the case of J0537-6910 and to a less extent J1801-2304, only a few
statistically compelling correlations are found between the size of a glitch
and the time to the preceding or the succeeding glitch (Melatos  _et al._ ,
2018). This may sound counterintuitive, if glitches are threshold-triggered
events, as both the starquake and the vortex avalanches scenarios seems to
indicate. However, this absence of obvious correlations can be reproduced in
terms of a state-dependent Poisson process for the internal stress
$\sigma(t)$, where the rate with which a glitch is triggered diverges as the
internal stress reaches a threshold $\sigma_{c}$ (Fulgenzi  _et al._ , 2017;
Carlin and Melatos, 2019). Alternatively, Fig. (2) suggests that the lack of
obvious correlations may be due to a process that does not conserve neither
$\sigma^{c}$ nor $\sigma^{low}$. Within the vortex avalanche scenario the non-
constancy of the upper and lower thresholds is naturally explained because the
vortex configurations is expected to move from a local minimum to another in a
complex landscape as sketched in Fig. 1: each starting minimum would be
associated to a different $\sigma^{c}$, and every final position in the
landscape to a certain $\sigma^{low}$. Unsurprisingly, the exact sequence of
such values would depend on the past history of the system (also in view of
the possible presence of hysteresis) and only very long observations could
unveil some statistical properties of such landscape (Antonelli and Haskell,
2020)666 Note added to the arxiv version: besides (Antonelli and Haskell,
2020), hysteresis and properties the pinning landscape are also explored in
the recent paper (Link and Levin, 2022), thanks to simulations of vortex
filaments similar to those developed in (Link, 2009). .

As the number of glitches detected in single pulsars grows enough for allowing
systematic and reliable studies of the statistical properties of the glitching
behavior, it will be more and more important to ask why and how these
properties change across the pulsar population in relation with the estimated
pulsar age, magnetic field, temperature and spin-down rate.

##  III Two-fluid equations and mutual friction

Since the starquake model alone seems to fail to explain the observational
pattern of pulsar glitches (see also (Chamel and Haensel, 2008) and (Haskell
and Melatos, 2015) on this point), much attention in the last decades has been
devoted to understanding and refining fluid theories for the several
thermodynamic and hydrodynamic degrees of freedom that are necessary to
describe neutron star matter at the macroscopic scale. In fact, in the
presence of superfluid phases, the set of thermodynamic variables needed to
define a homogeneous equilibrium state must be enlarged to account for the
possible presence of long lived metastable states with persistent currents. It
follows that the system is a multifluid, namely a conducting medium where
different (chemical or even abstract, like the entropy) species can flow with
respect to each other (Carter, 1989).

The basic multifluid formalism necessary to build pulsar glitch models is
based on the seminal works of Carter and collaborators (e.g., (Carter and
Khalatnikov, 1992; Langlois  _et al._ , 1998)) and has already been reviewed
in (Haskell and Sedrakian, 2017; Chamel, 2017a; Andersson, 2021). An extended
presentation of relativistic multifluids can be found in (Andersson and Comer,
2007), and a comparison of different hydrodynamic languages used to describe
the internal dynamics of superfluid NSs can be found in (Gavassino and
Antonelli, 2020). Moreover, aspects of superfluidity and superconductivity in
NSs and nuclear matter with a focus on pairing phenomena have been discussed
by (Sauls, 1989) and (Sedrakian and Clark, 2019).

Since the neutrons in the crust are expected to form Cooper pairs in the
${}^{1}S_{0}$ channel of the neutron-neutron interaction potential (i.e., the
superfluid order parameter in the crust is a scalar), the superfluid in the
crust is expected to behave analogously to the laboratory superfluid Helium-4,
whose hydrodynamics is well described by the HVBK hydrodynamics (Barenghi  _et
al._ , 2001; Sonin, 2016). On the other hand, according to the current
understanding of nuclear interactions obtained from nucleon scattering data,
the neutron pairing most favourable at densities above nuclear saturation is
probably in the ${}^{3}P_{2}$ state rather than the ${}^{1}S_{0}$ pairing
(Fujita and Tsuneto, 1972; Yasui  _et al._ , 2019). In the following we will
neglect the complex issues that arise when the order parameter for neutron
superfluidity is not a scalar (Brand and Pleiner, 1981), like the possible
presence of domain walls (Yasui and Nitta, 2020), exotic vortices of half-
integer circulation and spin textures (Leinson, 2020). Despite these
difficulties, some formulations for the neutron superfluid and proton
superconductor mixture in the core have been proposed as well (see, e.g.,
Mendell, 1998; Andersson  _et al._ , 2006; Glampedakis  _et al._ , 2011;
Gusakov and Dommes, 2016; Rau and Wasserman, 2020; Sourie and Chamel, 2021).

Almost all glitch models to date are based on the assumption that (locally) it
is possible to phenomenologically model neutron conduction in terms of two
coupled components (or “chemical species”): a first fluid that should
represent the neutron superfluid (usually labelled by an index $x=n$) and an
effective charge-neutral fluid ($x=p$). This latter is a mixture of charged
particles (or nuclear clusters if we are in the crust), and all the normal
(i.e., non-superfluid) neutrons. All these species can be combined in the
p-component as long as charge neutrality is satisfied over macroscopic regions
(local charge imbalances are expected to be equilibrated by the electron fluid
on very short timescales). Here, we adopt this simplified zero-temperature
two-fluid effective description, where the generation of entropy and heat by
friction is not taken into account777An exception is the “thermally driven”
glitch model, where the sudden vortex motion generates heat that affects the
vortex creep rate and the related angular momentum exchange between the
components (Link and Epstein, 1996; Larson and Link, 2002)..

_Fluid equations in the non-transfusive limit -_ We adopt the Newtonian
multifluid formalism of (Prix, 2004), in the limit where there are no chemical
reactions involving the two species (see, e.g., (Chamel and Carter, 2006;
Gusakov and Andersson, 2006; Andersson and Comer, 2011) for other formulations
of the hydrodynamic equations and (Gavassino and Antonelli, 2020) for a
dictionary to translate between different formalisms). This simplified
description is the most commonly used setting when it comes to glitch
modelling. The set of hydrodynamic equations consist of two ($x=n,p$)
conservation equations,

$\partial_{t}\rho_{x}+\nabla_{i}(\rho_{x}v^{x}_{i})=0$ (10)

and of two momentum equations,

$\rho_{x}\left(\partial_{t}+L_{v_{x}}\right)p^{x}_{i}+\rho_{x}\nabla_{i}\left(\tilde{\mu}_{x}-|v_{x}|^{2}/2\right)=F_{i}^{x}\,,$
(11)

where $L_{v_{x}}p_{i}^{x}$ is the Lie derivative of the momentum per particle
$p_{i}^{x}$ along the kinematic velocity field $v_{i}^{x}$,

$L_{v_{x}}p_{i}^{x}=v_{x}^{k}\partial_{k}p_{i}^{x}+p_{j}^{x}\partial_{i}v_{x}^{j}\,.$
(12)

In (11), $\tilde{\mu}_{x}$ is the specific chemical potential of the relative
component, and $\rho_{x}=m_{x}n_{x}$ is the mass density, where $m_{x}$ is the
baryon mass and $n_{x}$ is the effective baryon number density. The effect of
gravity or external fields can be included in the force density $F_{i}^{x}$.

_Velocity, momentum and entrainment -_ The field $v_{x}^{i}$ is sometimes
referred to as kinematic velocity to stress the fact that it is the quantity
that appears in the conservation equations (10). The momentum of the
superfluid species $p_{i}^{n}$ is, at the inter-vortex scale, proportional to
the gradient of the phase of the superfluid order parameter, so that it is an
irrotational field. Hence, $p_{i}^{n}$ is analogous to the “superfluid
velocity” of the Tisza-Landau or HVBK models for superfluid Helium (see Carter
and Khalatnikov, 1992; Andersson and Comer, 2011; Gavassino and Antonelli,
2020, and references therein). At the macroscopic scale $p_{i}^{n}$ is the
averaged momentum in the fluid element (possibly containing several vortex
lines), and its curl locally reflects the average arrangement of the quantum
vortices.

To gain some intuition on the distinction between velocity and momentum we
just have to recall their role. The velocity field can be defined starting
from a well-defined conserved current for the particles labelled by $x$: this
field is the product of a counting procedure, that is exactly the same that is
used in simple fluids (just applied on a subset of the total ensemble of
particles). The momentum, on the other hand, is a dynamical quantity. Assume a
homogeneous system where there are no gradients and for which we have already
defined the velocities $\mathbf{v}_{n}$ and $\mathbf{v}_{p}$ with the counting
procedure. We can now make a thought experiment and apply an external force
$\mathbf{F}^{ext}_{n}$ only to the neutrons: clearly, we expect
$\mathbf{v}_{n}$ to change, but, thanks to the interaction between the two
species, also $\mathbf{v}_{p}$ may change by some amount. If there were only
neutrons, one would have
$\mathbf{F}^{ext}_{n}=n_{n}\dot{\mathbf{p}}_{n}=\rho_{n}\dot{\mathbf{v}}_{n}$,
but now

$\mathbf{F}^{ext}_{n}=n_{n}\dot{\mathbf{p}}_{n}=\rho^{*}_{nn}\dot{\mathbf{v}}_{n}+\rho^{*}_{pn}\dot{\mathbf{v}}_{p}\,,$
(13)

where $\rho^{*}_{nn}$ and $\rho^{*}_{pp}$ are some effective mass densities
that depend on the details of the interaction and that were first introduced
in the description of laboratory superfluid mixtures (Khalatnikov, 1957;
Andreev and Bashkin, 1976).

The last piece of the puzzle that is missing is the following: why does it
happen that $\mathbf{p}_{n}$, and not $\mathbf{v}_{n}$, is the quantity that
is related to the gradient of the phase of the order parameter? At the formal
level, this is a self-consistency requirement for this kind of hydrodynamic
theory at the inter-vortex scale (Carter and Khalatnikov, 1992), see also the
discussion in Sec. 3 of (Gavassino and Antonelli, 2020). The intuition behind
this fact was first discussed by Anderson (1966) on the basis of the Josephson
effect in the context of superfluid Helium.

In the end, the relation (13) between $p_{i}^{x}$ and $v^{i}_{x}$ turns out to
be a convex combination of the two velocity fields, that can be written in
terms of a dimensionless parameter $\epsilon_{x}$ (Prix, 2004),

$p_{x}^{i}=m_{x}\left[(1-\epsilon_{x})v_{x}^{i}+\epsilon_{x}v_{y}^{i}\right]\quad\text{for}\quad
x\neq y\,.$ (14)

where $\rho_{n}\epsilon_{n}=\rho_{p}\epsilon_{p}$ to guarantee some
fundamental properties of the hydrodynamic stress tensor $T^{ik}$ to be
discussed below (Carter  _et al._ , 2006).

In NS cores, where both superfluid neutrons and the charged component move as
fluids, the entrainment is due to nuclear interactions between the species
that endows a neutron with a cloud of protons, and vice versa (Gusakov and
Haensel, 2005; Chamel and Haensel, 2006; Leinson, 2017; Chamel and Allard,
2019).

On the other hand, in the inner crust the entrainment phenomenon is thought to
be due to the Bragg scattering of the superfluid neutrons with the nuclear
clusters in the lattice (Chamel, 2012). Here, $\rho_{n}$ is not the total
density of neutrons, but rather the density of the dripped neutrons that
leaked outside the heavy nuclear clusters that constitute the solid lattice.
Hence, $\epsilon_{n}$ also defines the mobility of those dripped neutrons and,
ultimately, the effective density

$\rho^{*}_{n}=\frac{\rho_{n}}{1-\epsilon_{n}}$ (15)

that has to be ascribed to the part of neutrons that are free to move (Carter
_et al._ , 2006). The quantity $\rho^{*}_{n}$ also defines the effective
moment of inertia of the neutrons in glitch models (Chamel and Carter, 2006;
Antonelli and Pizzochero, 2017) and is chemically gauge independent (Carter
_et al._ , 2006; Gavassino  _et al._ , 2021), namely it does not depend on the
ambiguity regarding how many neutrons should be counted as free in the inner
crust. In the end, $\rho^{*}_{n}$ can be considered to be the effective
“superfluid density”, the quantity that is analogous to the superfluid density
in the Tisza-Landau model (Prix, 2004; Carter  _et al._ , 2006; Chamel and
Carter, 2006; Andersson and Comer, 2011). The precise calculation of
$\epsilon_{n}$ in the inner crust is particularly delicate (Martin and Urban,
2016; Watanabe and Pethick, 2017), as its value can also be very dependent on
the arrangement of the crustal lattice and the presence of disorder (Sauls
_et al._ , 2020).

_Mutual friction -_ In glitch models, the force $F_{i}^{x}$ in (11) accounts
for the so-called _mutual friction_ (Hall and Vinen, 1956), which allows for
momentum transfer between the two components. The adjective “mutual”
emphasises the fact that, in the absence of external forces (i.e.
$F_{i}^{n}+F_{i}^{p}=0$), the total momentum density

$j_{i}=n_{n}p_{i}^{n}+n_{p}p_{i}^{p}=\rho_{n}v_{i}^{n}+\rho_{p}v_{i}^{p}$ (16)

is conserved, in the sense that

$\partial_{t}j^{i}+\nabla_{k}T^{ik}=0\qquad
T^{ik}=\sum_{x}n_{x}v_{x}^{i}p_{x}^{k}+\Psi\delta^{ik}\,,$ (17)

where $\Psi$ is a thermodynamic potential that reduces to the usual pressure
when the two species comove. The fact $T^{ik}=T^{ki}$ is guaranteed by the
before-mentioned property $n_{n}\epsilon_{n}=n_{p}\epsilon_{p}$.

We now introduce the vortex-mediated mutual friction force. First, we have to
understand why $F_{i}^{x}$ is related to the presence of vortices. The idea is
simple: given a vortex configuration, then also the phase of the order
parameter is assigned. If the vortex configuration is frozen, the phase field
can not change, nor its gradient. It follows that in order to change
$\mathbf{p}_{n}$ we have to displace some vortices: the faster the
displacement, the larger will be $\dot{\mathbf{p}}_{n}$. Hence, it is
reasonable to expect that $\mathbf{F}^{n}$, which defines
$\dot{\mathbf{p}}_{n}$ in the homogeneous system, can be written as a function
of the local averaged vorticity and of the local average vortex velocity. In
particular, we can expect $\dot{\mathbf{p}}_{n}$ (and so $\mathbf{F}^{n}$) to
be proportional to the average “vortex current” in the fluid element. The
problem is to define the concept of current for string-like objects, but it is
possible to do so if the vortices are locally parallel (which does not
necessarily mean that they have to be straight or form an array at the
macroscopic scale). Under this hypothesis, the macroscopic vorticity

$\bm{\omega}=\nabla\times\mathbf{p}_{n}/m_{n}=\kappa
n_{v}\hat{\bm{\omega}}\qquad\kappa n_{v}=|\bm{\omega}|$ (18)

reflects completely also the vortex arrangement in the local fluid element;
$\kappa$ is the quantized circulation around a single vortex, $n_{v}$ is the
density of vortices measured in a plane orthogonal to $\hat{\bm{\omega}}$, the
local direction of the vortices. In other words, the hypothesis that the
quantum vortices are locally straight and parallel is so strong that there is
practically no loss of information in performing the average procedure over
the fluid element to obtain the macroscopic momentum $\mathbf{p}_{n}$. If, on
the other hand, vortices are locally arranged in a more complex, possibly
tangled, way, then the single vector $\bm{\omega}$ cannot be used to describe
the vortex arrangement. How to obtain the mutual friction in such a case is
still an open problem, even though progress can be made (Andersson  _et al._ ,
2007; Mongiovì  _et al._ , 2017; Haskell  _et al._ , 2020) by using concepts
of quantum turbulence imported from the study of vortex tangles in superfluid
Helium (Barenghi  _et al._ , 2001).

Now, we can derive the phenomenological form of $\mathbf{F}^{n}$, that must be
a Galilean invariant. Since we have at our disposal $\hat{\bm{\omega}}$ and
the relative velocity $\mathbf{v}_{np}=\mathbf{v}_{n}-\mathbf{v}_{p}$, we can
always project $\mathbf{F}^{n}$ on a (non-normalized) right-handed orthogonal
basis:

$\perp_{\hat{\bm{\omega}}}\mathbf{v}_{np}=-\hat{\bm{\omega}}\times(\hat{\bm{\omega}}\times\mathbf{v}_{np})\,,\qquad\hat{\bm{\omega}}\times\mathbf{v}_{np}\,,\qquad(\hat{\bm{\omega}}\cdot\mathbf{v}_{np})\hat{\bm{\omega}}$
(19)

where the operator $\perp_{\hat{\bm{\omega}}}$ is the projector that kills the
components parallel to $\hat{\bm{\omega}}$. Therefore, the force per unit
volume $\mathbf{F}^{n}$ can be decomposed as

$\mathbf{F}^{n}/(\rho_{n}\,|\bm{\omega}|)=\mathcal{B}_{d}\,\hat{\bm{\omega}}\times(\hat{\bm{\omega}}\times\mathbf{v}_{np})+\mathcal{B}_{c}\,\hat{\bm{\omega}}\times\mathbf{v}_{np}+\mathcal{B}_{k}\,\hat{\bm{\omega}}(\hat{\bm{\omega}}\cdot\mathbf{v}_{np})\,,$
(20)

where the factors $\rho_{n}$ and $|\bm{\omega}|$ have been introduced to make
the $\mathcal{B}$ coefficients dimensionless. If we want $\mathbf{F}^{n}$ to
be a force that drives the system towards the zero-lag state
$\mathbf{v}_{np}=0$, then we have to require that the dissipative part is
directed in the opposite direction to
$\perp_{\hat{\bm{\omega}}}\mathbf{v}_{np}$, so that $\mathcal{B}_{d}>0$.

The $\mathcal{B}_{k}$ coefficient may have a role when the vortices develop
Kelvin waves (helical displacements of a vortex’s core), as discussed by
Barenghi  _et al._ (1983). For strictly straight vortices $\mathcal{B}_{k}=0$
is usually neglected in glitch and Helium studies but it should be included in
the case a certain amount of quantum turbulence is expected888 The
phenomenological expression (20) is general enough that it can be supposed to
be valid also if vortices are not parallel to each other in the fluid element.
In particular, (20) also contains the Gorter-Mellink result for the mutual
friction in the isotropic turbulent regime (see (Andersson  _et al._ , 2007))
as the limiting case $\mathcal{B}_{c}=0$,
$\mathcal{B}_{d}=\mathcal{B}_{k}\propto|\mathbf{v}_{np}|^{2}$. The difference
with the non-turbulent case is that $\bm{\omega}$ alone is not sufficient to
describe the vorticity in the fluid element, so more involved additional terms
may be needed in (20).. Relativistic extension of the geometric construction
in (19) can be found in the seminal work of Langlois  _et al._ (1998), see
also (Gusakov, 2016; Gavassino  _et al._ , 2021).

Finally, the phenomenological $\mathcal{B}$ coefficients in (20) can be
functions of the lag $\mathbf{v}_{np}$ and the vorticity, making the mutual
friction non-linear in the lag (Antonelli and Pizzochero, 2017; Celora  _et
al._ , 2020). The hard task is to find their expressions from the physics at
the microscopic and mesoscopic scale.

###  III.1 The relation between mutual friction and vortex dynamics

The $\mathcal{B}$ coefficients in (20) are phenomenological macroscopic
coefficients that have been introduced on the basis of a purely geometrical
argument. We can shift the problem of their determination by assuming a
certain phenomenological model for the dynamics of vortex lines. This will not
solve all the problems, as vortex’s dynamics will be governed by a new set of
phenomenological parameters (again, introduced geometrically), that should
then be derived from the microscopic properties of the system at the scale of
the vortex’s core (roughly defined by the coherence length scale of the
superfluid).

We can now make clear the relation between the average vortex velocity and the
mutual friction outlined in the previous section (see (Sedrakian  _et al._ ,
1999), (Antonelli and Haskell, 2020) and references therein for further
details). For locally parallel vortex lines, we can introduce the average
vortex velocity $\mathbf{v}_{L}$ by writing down the equation of transport for
the macroscopic vorticity,

$\partial_{t}\bm{\omega}+\nabla\times(\bm{\omega}\times\mathbf{v}_{L})\,=0\,.$
(21)

In the above equation, $\mathbf{v}_{L}$ can be taken to be orthogonal to
$\bm{\omega}$ with no loss of generality. In this way, the current of vortex
lines in the fluid element is simply $n_{v}\mathbf{v}_{L}$. A geometrical
definition of $\mathbf{v}_{L}$ in general relativity is discussed in
(Gavassino  _et al._ , 2020). Equation (21) must be consistent with the
vorticity equation derived by taking the curl of (11) for $x=n$,

$\partial_{t}\bm{\omega}+\nabla\times(\bm{\omega}\times\mathbf{v}_{n})=\nabla\times(\rho_{n}^{-1}\,\mathbf{F}_{n})\,.$
(22)

Equations (21) and (22) are consistent if

$\mathbf{F}_{n}=-\rho_{n}\bm{\omega}\times(\mathbf{v}_{L}-\mathbf{v}_{n})=-\rho_{n}|\bm{\omega}|\mathbf{f}_{M}$
(23)

that tells us that the friction between the components is zero when the
vortices are advected in such a way that $\mathbf{v}_{L}=\mathbf{v}_{n}$. In
the above equation we also defined the quantity $\mathbf{f}_{M}$, the _Magnus
force_ acting on a straight vortex segment999 The quantity
$\mathbf{f}_{M}=\bm{\omega}\times(\mathbf{v}_{L}-\mathbf{v}_{n})$ is just a
rescaling of what in the literature is typically called Magnus force per unit
length (of vortex line), namely
$\rho_{n}\bm{\omega}\times(\mathbf{v}_{L}-\mathbf{v}_{n})$. Note that our
$\mathbf{f}_{M}$ is (apart from the overall sign) the same quantity appearing
in (20). .

If the vortices stick to the normal component because of pinning, then
$\mathbf{v}_{L}=\mathbf{v}_{p}$ and
$\mathbf{F}_{n}=\rho_{n}\bm{\omega}\times\mathbf{v}_{np}$, namely
$\mathcal{B}_{c}=1$ and $\mathcal{B}_{d}=0$, so that there is no dissipation
in this limit. More generally, if we have a model that allow us to calculate
the average vortex velocity in a fluid element
$\mathbf{v}_{L}(\mathbf{v}_{np})$ for a given value of the lag, then we can
substitute this expression back in (23) to obtain the expression of the
$\mathcal{B}$ coefficients in (20).

_Magnus and drag forces -_ In both Helium (Barenghi  _et al._ , 1983) and NSs
(Mendell, 1991; Andersson  _et al._ , 2006), the simplest phenomenological
model for vortex dynamics is in terms of balance of the Magnus force and a
drag force (proportional to a dimensionless drag parameter $\mathcal{R}$)
acting on a straight segment of vortex line,

$\mathbf{f}_{M}-\mathcal{R}\perp_{\hat{\bm{\omega}}}(\mathbf{v}_{L}-\mathbf{v}_{p})=\hat{\bm{\omega}}\times\left[\mathbf{v}_{L}-\mathbf{v}_{n}+\mathcal{R}\hat{\bm{\omega}}\times(\mathbf{v}_{L}-\mathbf{v}_{p})\right]=0\,,$
(24)

that can be inverted in the plane orthogonal to $\hat{\bm{\omega}}$ to find
$\mathbf{v}_{L}$ as a function of $\mathbf{v}_{n}$ and $\mathbf{v}_{p}$.
Plugging the result for $\mathbf{v}_{L}$ back into (23) allows to find

$\mathcal{B}_{c}\,=\,\frac{\mathcal{R}^{2}}{1+\mathcal{R}^{2}}\qquad\qquad\mathcal{B}_{d}\,=\,\frac{\mathcal{R}}{1+\mathcal{R}^{2}}\qquad\qquad\mathcal{B}_{k}\,=\,0\,.$
(25)

This tells us that, as long as $\mathcal{R}$ can be treated as a constant, the
mutual friction is linear in the lag $\mathbf{v}_{np}$.

In the outer core, the presence of entrainment causes the magnetisation of the
neutron superfluid vortices, generated by the current of protons that are
entrained by the local circulation associated with the vortex (Alpar  _et al._
, 1984b; Andersson  _et al._ , 2006). The scattering of the relativistic free
electrons with these magnetised vortex cores is the main process that is taken
into account in the estimates of $\mathcal{R}$.

In the inner crust, the drag depends also on the relative velocity between the
superfluid vortices and the normal component. For low values of this relative
velocity the drag is caused by phonon excitation in the crustal lattice
(Jones, 1990), while for higher values of this parameter the drag force is
caused by excitation of Kelvin waves on the superfluid vortices (Jones, 1992;
Epstein and Baym, 1992). A recent review on the calculation of the drag
parameter due to Kelvin waves excitation has been presented in Graber  _et
al._ (2018), while dynamical aspects related to the fact that vortices can
bend are discussed in (Link, 2009; Wlazlowski  _et al._ , 2016).

It is also worth mentioning that in the inner crust the complication induced
by the _chemical gauge_ problem (i.e., the problem of defining which neutrons
can be considered to be free, in the sense that they are not bounded into the
lattice clusters, see (Carter  _et al._ , 2006)) demands a more general
treatment than the one outlined above: in addition to the Magnus and the drag
forces in (24) also other contributions should be included in the vortex
equation of motion in order to ensure that the final hydrodynamic theory is
well behaved under redefinitions of the effective density of free neutrons
(Gavassino  _et al._ , 2021).

_Magnus, drag and pinning forces -_ If vortices are immersed in a non-
homogeneous medium, then their interaction with the substrate must be taken
into account. The vortex creep model of Alpar  _et al._ (1984a) deals with
this problem by assuming that the dependence of the jump rate from one pinning
center to another one is given by the typical Arrhenius formula, where the
activation energy for the jump is corrected by the presence of the background
superfluid flow. This makes $\mathbf{F}_{n}$ nonliner in the velocity lag
(Gügercinoǧlu and Alpar, 2017).

In fact, when an additional pinning force is included in equation (24) for the
motion for the vortices (as done in, e.g., (Sedrakian, 1995; Link, 2009;
Haskell and Melatos, 2016)), the $\mathcal{B}$ coefficients are not constant
anymore. For small $\mathbf{v}_{np}$ the parameter $\mathcal{B}_{d}$ is highly
suppressed, while $\mathcal{B}_{c}\approx 1$. When the lag reaches a certain
value (dependent on the details of the pinning landscape), the vortices start
to move and $\mathcal{B}_{d}$ tends to grow. For very high values of the lag,
the pinning interaction will be just a small correction to the Magnus and drag
forces, so that we can expect that the linear regime in (25) is eventually
recovered for very high values of the lag (remember, however, that the
parameter $\mathcal{R}$ may not be constant in the first place).

Pinning can also induce hysteresis in the pinning-depinning transition as the
lag is first increased and then decreased in a cyclic transformation, making
the behaviour of the parameters $\mathcal{B}_{c}$ and $\mathcal{B}_{d}$
history-dependent (Antonelli and Haskell, 2020). This may provide a natural
explanation for the behaviour of most pulsars, which is neither time-
predictable nor slip-predictable (as sketched in Fig. 2).

###  III.2 Circular models

In general, the fluid equations (10) and (11) can be solved only numerically
(Peralta  _et al._ , 2005; Howitt  _et al._ , 2016). Thus, in order to obtain
a glitch model that is more easily tractable, it is necessary to simplify the
problem by restricting ourselves to a specific subset of all possible fluid
motions. In this section, we assume a circular (not necessarily uniform)
velocity field and derive the equations of motion for the angular velocity of
the two components: this constitutes the the backbone of many glitch models.

Let us work with the usual cylindrical coordinates $(x,\varphi,z)$, where $x$
is the distance from the rotation axis of the star (directed along $z$). The
spherical radius is $r=\sqrt{x^{2}+z^{2}}$. In the absence of precession and
neglecting meridional flows101010 In General Relativity this would be the
circular spacetime approximation (Sourie  _et al._ , 2016; Antonelli  _et al._
, 2018; Gavassino  _et al._ , 2020; Sourie and Chamel, 2021). , the two
species can only rotate as

$\begin{split}\mathbf{v}_{p}&=x\,\Omega_{p}(x,z,t)\,\mathbf{e}_{\varphi}\\\
\mathbf{v}_{n}&=x\,\left(\Omega_{p}(x,z,t)+\Omega_{np}(x,z,t)\right)\,\mathbf{e}_{\varphi}\,.\end{split}$
(26)

It is also useful to define an auxiliary angular velocity variable directly
from the momentum $\mathbf{p}_{n}$,

$\Omega_{v}(x,z,t)=p_{n}^{\varphi}/x=\Omega_{p}+(1-\epsilon_{n})\Omega_{np}\,.$
(27)

Therefore, the macroscopic vorticity (18) reads

$\bm{\omega}\,=\,-x\partial_{z}\Omega_{v}\,\mathbf{e}_{x}\,+\,(2\Omega_{v}+x\partial_{x}\Omega_{v})\,\mathbf{e}_{z}\,,$
(28)

that reduces to the standard results $\bm{\omega}=2\Omega_{v}\,\mathbf{e}_{z}$
and $n_{v}=2\Omega_{v}/\kappa$ if the superfluid is rigidly rotating at the
macroscopic scale.

We can use (28) and (20) to obtain (we are working with locally parallel
vortices, so $\mathcal{B}_{k}=0$)

$\mathbf{F}_{n}=\rho_{n}\,x\,\Omega_{np}\Big{(}\mathcal{B}_{c}\left[x\partial_{z}\Omega_{v}\mathbf{e}_{z}+\left(2\Omega_{v}+x\partial_{x}\Omega_{v}\right)\mathbf{e}_{x}\right]-\mathcal{B}_{d}\,|\bm{\omega}|\,\mathbf{e}_{\varphi}\Big{)}\,.$
(29)

We are interested in specifying the equations (11) directed along the
azimuthal $\varphi$ direction for both the normal and superfluid components.
We are not interested in the other two directions, namely the components of
(11) along $\mathbf{e}_{x}$ and $\mathbf{e}_{z}$, because they do not contain
any time derivative: their role is to define the instantaneous configuration
of the rotating NS, whose deviation from the spherical hydrostatic equilibrium
must be constantly adjusted as $\Omega_{n}$ and $\Omega_{p}$ vary in time. In
the slowly rotating limit (Andersson and Comer, 2001), it is not a bad
approximation to neglect the components of (11) along $\mathbf{e}_{x}$ and
$\mathbf{e}_{z}$ and assume some given stratification of the star, i.e.,
$\rho_{n,p}=\rho_{n,p}(t,x,z)\approx\rho_{n,p}(r)$. Hence, we will restrict
ourselves to the dynamics along $\mathbf{e}_{\varphi}$, that is particularly
simple since the Lie derivative and the gradient term (because of the axial
symmetry hypothesis) do not have any component directed along
$\mathbf{e}_{\varphi}$. A direct substitution in (11) gives:

$\displaystyle
x\rho_{n}(\partial_{t}\Omega_{n}-\epsilon_{n}\partial_{t}\Omega_{np})=x\rho_{n}\partial_{t}\,\Omega_{v}\,=\,{F}^{\varphi}_{n}$
(30) $\displaystyle
x\rho_{p}(\partial_{t}\Omega_{p}+\epsilon_{p}\partial_{t}\Omega_{np})\,=\,-{F}^{\varphi}_{n}+F^{\varphi}_{Bp}\,,$
(31)
$\displaystyle{F}^{\varphi}_{n}=-\rho_{n}x\mathcal{B}_{d}|\bm{\omega}|\Omega_{np}$
(32)
$\displaystyle|\bm{\omega}|=\sqrt{(2\Omega_{v}+x\partial_{x}\Omega_{v})^{2}+(x\partial_{z}\Omega_{v})^{2}}$
(33)

where we have also added an extra external force $F^{\varphi}_{Bp}$ to model
the braking effect of dipole radiation on the normal component (an external
torque due to gravitational wave emission would act on both components). Note
that, despite $\mathcal{B}_{c}\neq 0$, this parameter only appears in the two
equations that define the stellar structure. If we sum the two dynamical
equations in (30) and (31) we have ($\rho=\rho_{n}+\rho_{p}$)

$x(\rho\,\partial_{t}\Omega_{p}\,+\,x\rho_{n}\partial_{t}\Omega_{np})\,=\,{F}^{\varphi}_{Bp}\,.$
(34)

that, in the limit ${F}^{\varphi}_{Bp}=0$ is the $\varphi$ component of
$\eqref{prezzemolo}$. The above equation may be used in place of, say, (31).

Another widely used assumption in glitch modelling is that the normal
component rotates rigidly. This assumption is beyond the fluid equations (11),
as it is motivated on some piece of physics that is not implemented in the
simple dynamical equations we are using here. Hence, when doing this we have
to ensure the consistency of the fluid equations: if we simply put
$\Omega_{p}=\Omega_{p}(t)$ into the above equation, then inconsistencies may
arise. Assuming rigid motion for the normal component forces us to substitute
one of the local equations (say (34) for convenience) with an integral one,

$\int
d^{3}x\,x^{2}\,(\rho\,\partial_{t}\Omega_{p}\,+\,\rho_{n}\,\partial_{t}\Omega_{np})\,=\,\int
d^{3}x\,x\,\,{f}^{\varphi}_{Bp}\,.$ (35)

It is only at this point that we can really restrict the motion of the charged
component to be perfectly rigid: assuming $\Omega_{p}=\Omega_{p}(t)$ in the
above equations allows us to bring it outside the integral sign to obtain

$I\dot{\Omega}_{p}\,+\,\int
d^{3}x\,x^{2}\rho_{n}\partial_{t}\Omega_{np}\,=\,-I|\dot{\Omega}_{\infty}|$
(36)

where $I$ is the total moment of inertia of the star for the given
stratification $\rho_{n}(x,z,t)+\rho_{p}(x,z,t)\approx\rho(r)$ and
$|\dot{\Omega}_{\infty}|$ is a phenomenological parameter (possibly a function
of time and/or the angular velocities) that sets the strength of the total
braking torque,

$|\dot{\Omega}_{\infty}|\,=\,\frac{1}{I}\,\left|\,\int
d^{3}x\,x\,\,{f}^{\varphi}_{Bp}\,\right|\,.$ (37)

The other equation that we need can be taken to be (30),

$\partial_{t}\,\Omega_{v}(t,x,z)\,=\,-\mathcal{B}_{d}\,|\bm{\omega}|\,\Omega_{np}\,\,=\,-\mathcal{B}_{d}\,|\bm{\omega}|\,\frac{\Omega_{vp}}{1-\epsilon_{n}}\,,$
(38)

where

$\Omega_{vp}(x,z,t)=\Omega_{v}-\Omega_{p}=(1-\epsilon_{n})\,\Omega_{np}$ (39)

is the lag between the superfluid angular velocity111111 We recall that the
variable $\Omega_{v}$ is defined by the vortex configuration and is,
basically, the _superfluid velocity_ of Tisza and Landau (not to be confused
with the _velocity of the species that is superfluid_ , $\Omega_{n}$). Hence,
$\Omega_{v}$ is not an exotic variable, but, rather, it is the standard choice
for laboratory superfluids. It has also the nice property of being chemical-
gauge independent, see (Carter  _et al._ , 2006; Gavassino  _et al._ , 2021).
$\Omega_{v}$ and the velocity of the normal component.

If we want to use $\Omega_{v}$ in place of $\Omega_{n}$ as primary variable,
then (36) can be equivalently written as

$I\dot{\Omega}_{p}\,+\,\int
d^{3}x\,x^{2}\,\rho^{*}_{n}\,\partial_{t}\Omega_{vp}\,=\,-I|\dot{\Omega}_{\infty}|\,,$
(40)

where $\rho^{*}_{n}$ is the entrainment corrected superfluid density
introduced in (15). In this way the equations are formulated in terms of the
chemically gauge invariant quantities $\Omega_{p}$, $\Omega_{v}$,
$\rho^{*}_{n}$ and $\rho$.

_Axially symmetric dynamical equations -_ To summarize, for a given
stratification of the NS defined by $\rho_{n,p}\approx\rho_{n,p}(r)$ and
$\epsilon_{n}\approx\epsilon_{n}(r)$, the axially symmetric dynamical
equations are (36) and (38). Due to the form of the vorticity in (33), the two
equations can be conveniently written by selecting as primary variables
$\Omega_{p}(t)$ and $\Omega_{v}(t,x,z)$. By introducing an average operation
over the local moment of inertia defined as

$dI_{v}=d^{3}x\,x^{2}\,\rho^{*}_{n}\qquad\quad I_{v}=\int
dI_{v}\qquad\quad\langle\,f\,\rangle=I_{v}^{-1}\\!\\!\int dI_{v}\,f(t,x,z)\,$
(41)

the system defined by (40) and (38) is equivalent to

$\displaystyle
I\partial_{t}\,{\Omega}_{p}(t)\,+\,I_{v}\,\partial_{t}\langle\,\Omega_{vp}\,\rangle\,=\,-I|\dot{\Omega}_{\infty}|$
(42)
$\displaystyle\partial_{t}\Omega_{v}(t,x,z)\,=\,-\mathcal{B}_{d}(r,\Omega_{vp})\,|\bm{\omega}|\,\frac{\Omega_{vp}(t,x,z)}{1-\epsilon_{n}(r)}$
(43)

where $\Omega_{vp}(t,x,z)=\Omega_{v}(t,x,z)-\Omega_{p}(t)$,
$\langle\Omega_{vp}\rangle=\langle\Omega_{v}\rangle-\Omega_{p}$ and
$|\bm{\omega}|$ is given in (33). Since the stratification is given, the drag
parameter $\mathcal{B}_{d}$ can be assumed to be fixed and almost spherically
stratified, i.e. $\mathcal{B}_{d}\approx\mathcal{B}_{d}(r)$, or to contain
also an extra dependence on the lag $\Omega_{vp}$, see the discussion below
(20).

Notice that to derive (42) and (43) we have used only the momentum equations
(11) because the conservation equations (10) are automatically satisfied when
assuming axial symmetry and no meridional flows.

Finally, (42) and (43) have been derived in the limit of slack vortices, in
which the effect of bending a quantum vortex is neglected (Antonelli  _et al._
, 2018). If the vortex tension (Andersson  _et al._ , 2007; Link, 2014;
Haskell and Melatos, 2016) is taken into account (as it is sometimes done in
the HVBK equations, see (Barenghi  _et al._ , 2001)), then an additional term
containing derivatives along $z$ and $x$ should appear in the left hand side
of (43). Opposite to the slack limit, we have the extreme scenario in which
$\Omega_{v}$ is cylindrical (i.e., it has no $z$ dependence, so that
$\bm{\omega}$ is everywhere parallel to the rotation axis), which is formally
recovered if the vortex tension is infinite. Such an unrealistic scenario is,
however, still useful from the practical point of view as it allows to further
simplify the equations: several glitch models are actually _cylindrical_ ,
meaning that they are built on this simplifying assumption, see e.g. (Alpar
_et al._ , 1984a; Pizzochero, 2011; Haskell  _et al._ , 2012; Antonelli and
Pizzochero, 2017; Graber  _et al._ , 2018; Khomenko and Haskell, 2018; Haskell
_et al._ , 2018).

_Models with rigid components -_ The phenomenological model of Baym and
collaborators defined by equations (1) and (2) lacks direct connection with
the local properties of neutron star matter. However, it is possible to link
its parameters with the quantities appearing in (42) and (43). The rigid model
can be recovered by averaging (43) over the superfluid region by using (41),

$\langle\partial_{t}\Omega_{v}\rangle\,=\,\partial_{t}\langle\Omega_{v}\rangle\,=\,-\left\langle\mathcal{B}_{d}\,|\bm{\omega}|\,\frac{\Omega_{vp}}{1-\epsilon_{n}}\right\rangle$
(44)

Clearly, the approximate reduction to a rigid model by performing the average
of the equations will work well if stratification and spatial correlations
between the different quantities are weak, that is

$\left\langle\mathcal{B}_{d}\,|\bm{\omega}|\,\frac{\Omega_{vp}}{1-\epsilon_{n}}\right\rangle\approx\left\langle\,\frac{\mathcal{B}_{d}\,|\bm{\omega}|}{1-\epsilon_{n}}\right\rangle\langle\Omega_{vp}\rangle\approx
2\langle\Omega_{v}\rangle\left\langle\,\frac{\mathcal{B}_{d}\,}{1-\epsilon_{n}}\right\rangle\langle\Omega_{vp}\rangle\,.$
(45)

The factor $\langle\Omega_{p}\rangle\approx\Omega_{p}$ is almost constant, but
the above expression can still be highly non-linear if $\mathcal{B}_{d}$
depends on the local lag. The linear model of Baym and collaborators can only
be found if $\mathcal{B}_{d}$ is locally constant over the timescale of
interest. In this case, a direct comparison of (42) and (44) with (1) and (2)
allows us to make the following identifications:

$\begin{split}&I_{s}=I_{v}\qquad\qquad\,\,I_{o}=I-I_{v}\\\
&\Omega_{s}=\langle\Omega_{v}\rangle\qquad\quad\Omega_{o}=\Omega_{p}\\\
&\frac{1}{\tau}=\frac{I-I_{v}}{I\,\tau_{exp}}\,\approx\,2\Omega_{p}\left\langle\,\frac{\mathcal{B}_{d}\,}{1-\epsilon_{n}}\right\rangle\end{split}$
(46)

For example, this tells us how a fit of the timescale $\tau_{exp}$ from
observations can be interpreted in terms of the physical input at smaller
scales. General relativistic corrections should be also taken into account,
which affect both the expressions for the moments of inertia and for $\tau$
(Antonelli  _et al._ , 2018; Sourie  _et al._ , 2016; Gavassino  _et al._ ,
2020).

A similar reduction to a rigid model can be performed by favouring
$\Omega_{n}$ in place of $\Omega_{v}$. In this case, extra internal torques
due to the entrainment coupling appear in the final equations (Chamel and
Carter, 2006; Sidery  _et al._ , 2010).

Finally, the superfluid region (i.e., the domain of $\Omega_{v}$, that extends
in the inner crust and outer core) can also be split into an arbitrary number
of different, non-overlapping, sub-regions to obtain the multi-component
generalization of the model of Baym and collaborators (see (Montoli  _et al._
, 2020a) for the formal derivation of the three-components analogue of Baym’s
model).

##  IV Dynamical phases and the Vela’s 2016 glitch precursor

It is interesting to make some general statements about the dynamics defined
by (42) and (43). Those equations do not allow to track the stresses in the
solid components and to directly model starquakes and changes of the moments
of inertia $I$ and $I_{v}$ should be implemented by hand, as in Baym’s model
of Sec. II.1. Hence, we assume here the Anderson and Itoh (1975) paradigm, in
which the moments of inertia $I$ and $I_{v}$ can be taken to be constant.

The angular momentum balance equation (42) and equation (43) have two
different roles. The momentum balance (42) allows to identify some general
dynamical phases, but the ability of a real pulsar to switch from one phase to
another depends on the details of friction (i.e., on the details of the
$\mathcal{B}_{d}(\Omega_{vp})$ function) or, more generally, on the
statistical properties of the landscape in Fig. 1. In other words, the
possibility to explore certain phases critically depends on how the angular
momentum can be exchanged between the components and transported between
different regions during the pulsar’s life.

Despite we still do not have a complete and satisfactory understanding of
these issues, it is still instructive to consider each dynamical phase as a
separate (theoretical) possibility. The following analysis is valid for a
system where the superfluid component can rotate non-uniformly (this system
has an infinite number of degrees of freedom, so that each phase is highly
degenerate and can be realized in a number of ways), but the same reasoning
can be applied to simpler models with two or more rigid components, that have
a finite number of degrees of freedom and are defined by ordinary differential
equations.

We plot the different phases defined by the angular momentum balance (42) in
Fig. 4.

Figure 4:  The evolution of the total angular momentum (42) can be regarded as
a constraint (the blue line) in the
$(\partial_{t}\langle{\Omega}_{vp}\rangle,\partial_{t}{\Omega}_{p})$ plane,
whatever the number of internal “components” (formally infinite if the
superfluid can rotate non-uniformly) or the details of the mutual friction.
The numbers indicate the different dynamical phases described in the main
text. The slope of the curve depends on the ratio $I_{v}/I$.

Going from negative to positive values of $\partial_{t}\Omega_{p}$ we have the
following phases:

  1. 1.

_Inward vortex creep -_ From (42) we see that
$\partial_{t}\Omega_{p}<-I/(I-I_{v})|\dot{\Omega}_{\infty}|$ can be realized
only if the internal torque contributes to slowing down the observable
component (the internal torque has the same sign as the external one). Since
$\mathcal{B}_{d}$ is a positive function (or functional) of the local lag, as
demanded by the second law of thermodynamics (Gusakov, 2016; Gavassino  _et
al._ , 2021), we should have that this phase can be realized only if there is
a sufficiently extended region in which the lag is negative (in fact,
$\text{sign}\Omega_{vp}=\text{sign}\Omega_{np}$). This also implies that the
vortex velocity is, on average, directed inward.

  2. 2.

_Average pinning -_ In this dynamical phase, the superfluid does not spin-
down. Mathematically, we have $\partial_{t}\langle\Omega_{v}\rangle=0$. The
lag between the two components builds up as
$\partial_{t}\langle\Omega_{vp}\rangle=-\partial_{t}\Omega_{p}$. Equation (40)
implies that $\partial_{t}\Omega_{p}=-I/(I-I_{v})|\dot{\Omega}_{\infty}|$.
This phase can be achieved with perfect pinning (the superfluid is _completely
decoupled_ from the normal component) or by a balance of inward and outward
creep (unlikely unless the actual lag has very big gradients and is partially
reversed).

  3. 3.

_Moderate outward creep -_ The average lag satisfies
$0<\partial_{t}\langle\Omega_{vp}\rangle<-\partial_{t}\Omega_{p}$, meaning
that the superfluid is only partially decoupled. This state is likely to occur
during the long term post-glitch relaxation: the reduced post-glitch lag
slowly increases, i.e. $\partial_{t}\langle\Omega_{vp}\rangle>0$. Consequently
the superfluid gradually recouples with the observable component (the system
is expected to reach the steady state if no other glitches occur in the
meanwhile).

  4. 4.

_Steady state -_ Here the coupling of the superfluid with the normal component
is stronger than in the previous phase, and such that the lag is on average
constant in time, namely $\partial_{t}\langle\Omega_{vp}\rangle=0$. Despite
both components spin-down at the same rate
$\partial_{t}\Omega_{p}=\partial_{t}\langle\Omega_{v}\rangle=-|\dot{\Omega}_{\infty}|$,
the local lag $\Omega_{vp}$ does not have to do so: for example, this
situation may be realized also if a sequence of packets of unpinned vortices
creep outward.

  5. 5.

_Speed up -_ This phase can be defined as
$-|\dot{\Omega}_{\infty}|<\partial_{t}\Omega_{p}<0$, meaning that the
observable component is spinning down but the average lag is decreasing in
such a way that
$-I/I_{v}|\dot{\Omega}_{\infty}|<\partial_{t}\langle\Omega_{vp}\rangle<0$: the
normal component is receiving some extra angular momentum from the superfluid
but at such a low rate that it is still spinning down. Such a phase has been
probably been observed in the Crab, and it has been speculated that it may be
the result of a small increase in the Crab pulsar’s internal temperature
(Vivekanand, 2017) that has the effect of increasing the outward vortex creep
rate (Alpar  _et al._ , 1984a).

  6. 6.

_Spin-up -_ In this phase the observable component spins up,
$\partial_{t}\Omega_{p}>0$ and the lag is decreasing with a high rate
$\partial_{t}\langle\Omega_{vp}\rangle<-I/I_{v}|\dot{\Omega}_{\infty}|$. In
real glitches the inequality is actually believed to be stronger,
$\partial_{t}\langle\Omega_{vp}\rangle\ll-I/I_{v}|\dot{\Omega}_{\infty}|$ but
a slow spin-up event has been observed in the Crab pulsar following a fast,
unresolved, glitch (Shaw  _et al._ , 2018). The particular value
$\partial_{t}\Omega_{p}=0$ can be realized only if the internal torque due to
the superfluid exactly balances the external one. This can happen
instantaneously when the observable component reaches the maximum amplitude in
a glitch. In the vicinity of this event the lag can be both positive or
negative (depending on the fact that an overshoot takes place or not, see also
Sec. VI).

This analysis is independent on the details of the mutual friction, but it
depends on the fact that the moments of inertia are assumed to be constant.

_The Vela precursor -_ Ashton  _et al._ (2019) observed a feature in the Vela
pulsar’s phase residuals around the Vela’s 2016 glitch recorded by Palfreyman
_et al._ (2018). This feature can be interpreted as a spin-down (an
antiglitch), occurring just before the main spin-up event and possibly linked
to the triggering of the glitch. Assuming that such observed behaviour can be
explained only in terms of angular momentum exchange between the superfluid
and the normal component, it is interesting to try to understand to which
dynamical phase this spin-down “precursor” may belong to.

According to the fit performed in (Ashton  _et al._ , 2019), this sudden
decrease of the Vela’s rotational frequency ($\nu\approx 11\,$Hz) has an
amplitude $\Delta\nu\sim-5\times 10^{-6}\,$Hz and it is likely to occur over a
time span of $\Delta t\sim 100$s or shorter, implying $\dot{\Omega}_{p}\approx
2\pi\Delta\nu/\Delta t\sim-10^{-7}\,$rad/s2. Since a spin-down of the normal
component would increase the lag, this precursor may also trigger the
subsequent glitch by causing a critical lag between crustal superfluid and the
crust (Ashton  _et al._ , 2019).

We should compare the above numbers with the secular spin-down parameter of
the Vela, $|\dot{\Omega}|_{\infty}\approx 10^{-10}$rad/s2: by looking at Fig.
4, we see that the Vela should be in either phase 1, 2 or 3 during the
precursor. Let us assume that phase 2 is realized, which means a complete
decoupling of the superfluid (or that in different parts of the superfluid the
creep is balanced between outward creep and inward creep). In this case, the
observable component could spin-down at a faster rate as the external torque
acts on a reduced moment of inertia, especially if $I_{v}\approx I$, which is
the extreme situation where the NS is almost entirely superfluid. Assuming
phase 2, we have $-\dot{\Omega}_{p}/|\dot{\Omega}|_{\infty}=I/(I-I_{v})\approx
10^{3}$, so that the the blue line in Fig. 4 would be practically vertical.
This scarcity of normal matter in a NS (i.e., $I-I_{v}\sim 10^{-3}I$) is not
predicted by any reasonable model of NS structure. Hence, we may assume that
the Vela was in phase 1 (phase 3 is excluded since already phase 2 is not a
viable option). However, this is also not a realistic possibility: vortices
can undergo inward creep (so that the mutual friction torque acts in the same
direction of the external torque) if the lag is momentarily reversed, which is
impossible as this event immediately precedes a standard spin-up glitch, where
the lag has to be positive.

In conclusion, such a huge “antiglitch” precursor cannot be explained in terms
of the standard cartoon of pulsar glitches without invoking the occurrence of
the _inward creep_ phase just before the main spin-up event. Since we only
used the general equation (42) for the dynamics of the total angular momentum,
which is insensitive to all the details of the internal dynamics, assuming
complex models for the mutual friction or adding additional (fluid or rigid)
components should not change this conclusion.

_Possible explanations of the precursor -_ Gügercinoğlu and Alpar (2020)
speculated that the precursor can be a consequence of the formation of high
vortex density regions surrounded by vortex depletion regions called _traps_
(Cheng  _et al._ , 1988). In this case, the displacement of vortices is a
byproduct of a starquake (that eventually is also the cause of the glitch but
does not significantly modify the moments of inertia). According to this
picture, the dislocations of the lattice induced by the quake provide stronger
pinning sites for the vortex lines. These new sites can attract vortices over
microscopic distances, hence the formation of traps. Since the velocity and
direction of the displaced vortices depend on the concentration of
dislocations, it is still not completely clear how the creation of traps can
induce the average inward creep needed to justify the precursor.

A tentative explanation in terms of stochastic fluctuations of the internal
torque has also been proposed by Ashton  _et al._ (2019) and Carlin and
Melatos (2020). In the case the continuous transfer of angular momentum from
the superfluid to the normal component is particularly noisy, then the
precursor may just be the effect of an abnormal fluctuation in the internal
torque. Such a fluctuation increases in a short time the average lag between
the components, possibly pushing it above the depinning threshold in some
regions, hence triggering the subsequent glitch. It follows that understanding
if such internal torque fluctuations – that, in turn, are fluctuations in the
state of motion of many vortices – can actually take place in the real system
is key to falsify this possibility. Moreover, the statistical occurrence of
such large fluctuations should be closely related to the observed statistical
occurrence of glitches.

Future observations of well-resolved glitches have, thus, the potential to
challenge the current understanding and, possibly, to inspire new physical
paradigms for glitch modelling. However, there is also the additional
complication that what we directly observe is not the normal component but
rather a signal generated in the magnetosphere of the pulsar, that has its own
dynamics. Since the antiglitch cannot be attained by invoking neither the
internal torque nor the external one (or both, as in phase 1), a more
conservative interpretation is that the precursor does not represent a a
change in $\Omega_{p}$. On the contrary, it is likely to reflect a change
pertaining the emission region of the Vela (Montoli  _et al._ , 2020a).
Moreover, a recent analysis revealed a flickering magnetospheric activity of
the Vela around the 2016 glitch epoch (Ashton  _et al._ , 2020). If confirmed,
this would provide a further clue that the magnetosphere is indeed a complex
dynamical region; see also (Bransgrove  _et al._ , 2020), who simulated the
magnetospheric perturbations induced by a crustquake.

##  V Constraints on neutron star structure

In this section we present some methods that can be used to extract
information on a glitching neutron star through its observed glitch behaviour.
As often occurs in astrophysics, extrapolation of information is an indirect
process, whose result depends on some state of the art microscopic input
parameters and on the ability to build models that reproduce the relevant
piece of physics. In the following, we describe two static constraints on the
pinning forces and on entrainment, coming, respectively, from the observation
of glitches of large amplitude and very active pulsars. Then, we discuss the
information that can be extracted from the dynamical modelling of the spin-up
phase.

###  V.1 Static constraint from the maximum glitch amplitude

Observations of glitches of large frequency jump $\Delta\nu$ can be used to
understand if the microscopically calculated pinning forces are enough to
explain it. Within the assumption of circular motion described in the previous
section, from (42) we have that (ignoring the external torque that acts on
much longer timescales)

$2\pi\Delta\nu=-\frac{I_{v}}{I}\langle\Delta\Omega_{vp}\rangle$ (47)

where $\Delta\Omega_{vp}$ is the difference between the lag after and before
the glitch, locally. The exact difference $\Delta\Omega_{vp}$ is impossible to
observe or estimate and it also depends on the past history of the pulsar.
However, we can give a theoretical upper limit $\Omega^{*}_{vp}$ to the local
lag $\Omega_{vp}$: defining as $f_{P}$ the maximum Magnus force per unit
length that pinning can sustain, from $|\mathbf{f}_{M}|=f_{P}$ with
$\mathbf{v}_{L}=\mathbf{v}_{p}$ we obtain

$\Omega_{vp}(t,x,z)<\Omega^{*}_{vp}(x,z)=\frac{f_{P}(r)(1-\epsilon_{n}(r))}{\kappa\,\rho_{n}(r)\,x}\,.$
(48)

At this point, an upper limit $\Delta\Omega_{\rm max}$ on the glitch amplitude
can be obtained by assuming that before the glitch the angular momentum
reservoir was full (its maximum capacitance being defined by the pinning
forces), while it is fully emptied after the glitch. It follows that

$2\pi\Delta\nu<\Delta\Omega_{\rm
max}=\frac{I_{v}}{I}\langle\Omega^{*}_{vp}\rangle\,,$ (49)

namely (we recall that $\nu$ is the frequency, not the angular velocity)

$2\pi\Delta\nu<\Delta\Omega_{\rm max}=\frac{\pi^{2}}{I\kappa}\int
dr\,r^{3}f_{P}(r)\,.$ (50)

Interestingly, the above result does not depend neither on the entrainment
parameter $\epsilon_{n}$ nor the density $\rho_{n}(r)$, but it only depends on
the radial profile $f_{P}(r)$ and on the total moment of inertia $I$, which in
turn depends on the total density profile $\rho(r)$ only. Clearly, the region
where $f_{P}(r)>0$, namely the region that contributes to the integral in
(50), should be contained within the region defined by $\rho_{n}(r)>0$.
Moreover, the constraint (50) has been derived in the limit of slack vortices,
but the exact same result is valid also for the extreme scenario of straight
vortices discussed below (43).

When general relativistic corrections in the slow rotation approximation
(Hartle, 1967; Andersson and Comer, 2001) are taken into account, the
constraint (50) reads (Antonelli  _et al._ , 2018)

$\Delta\Omega_{\rm
max}=\frac{\pi^{2}}{I\kappa}\int_{0}^{R_{d}}\mathrm{d}r\,r^{3}\,\mathrm{e}^{\Lambda(r)}\,\frac{\mathcal{E}(r)+P(r)}{m_{n}n(r)c^{2}}f_{P}(r)\,,$
(51)

where $\Lambda$ is the metric function which describes length dilatation,
$\mathcal{E}$ is the total mass-energy density and $P$ is the pressure and $n$
the total baryon density (the Equation of State, EOS, is $E=E(n)$ and
$P=P(n)$). Here, the radial profiles are obtained by integrating the TOV
equations, in the sense that the corrections to the stellar structure of the
second order in $\Omega_{p}$ can be safely neglected for slowly rotating two-
component neutron stars (Andersson and Comer, 2001).

Again, the result in (51) depends only on the pinning force profile, the mass
of the star and the EOS. Thus, once the microphysical parameters like the EOS
and the pinning forces are fixed, it provides a way to constrain the mass of a
pulsar via the observation of a large glitch.

Figure 5: The maximum glitch amplitude $\Delta\Omega_{\rm max}$ in equation
(51) calculated as a function of mass for different unified EOSs (Antonelli
_et al._ , 2018): SLy4 (Douchin and Haensel, 2001), BSk20 and BSk21 (Goriely
_et al._ , 2010). The largest glitch displayed by some pulsars is also shown:
the record-holder is J0205+6449 (Livingstone  _et al._ , 2009), that should
have a mass lower than $1.1-1.2\,M_{\odot}$, depending on the EOS employed.
The pinning forces employed are the ones derived in (Seveso  _et al._ , 2016)
for the inner crust.

This is represented in Fig. 5; not surprisingly, lighter pulsars can, in
principle, undergo larger glitches. If we observe a large glitch of amplitude
$\Delta\nu$ from a pulsar, all the input which predict a theoretical glitch
smaller than the observed one are ruled out, at least for that object. By
fixing the EOS and $f_{P}$, this translates into an upper limit on the mass of
a particular glitching object. To date, the pulsar that displayed the glitch
of largest amplitude is J0205+6449 (Livingstone  _et al._ , 2009), that,
therefore, should have a mass lower than $1.1-1.2\,M_{\odot}$ when the pinning
forces of Seveso  _et al._ (2016) are used in the inner crust. Since the
theoretical minimum mass of a neutron star estimated from supernovae
explosions is $\sim 1.17\,M_{\odot}$ (Suwa  _et al._ , 2018), this means that
the input used to calculate $\Delta\Omega_{\rm max}$ in Fig. 5 is only
marginally acceptable. Hence, the largest glitch of J0205+6449 tells us that
it may be necessary to extend the pinning region in the outer core, especially
if even larger glitches will be detected in the future or new and more refined
microscopic calculations will provide estimates of pinning forces in the crust
that are lower than the ones used here.

###  V.2 Static constraint from the activity parameter

Other type of information can be obtained by considering the so-called glitch
activity, that represents the mean angular acceleration due to glitches (Lyne
_et al._ , 2000; Fuentes  _et al._ , 2017). For a pulsar that has undergone
$N_{\rm gl}$ glitches of amplitude $\Delta\Omega_{i}$ at times $t_{i}$ with
size $\Delta\Omega_{i}$ in a time span $T_{obs}$, the absolute activity can be
defined as

$\mathcal{A}_{a}\,=\,\frac{1}{T_{obs}}\sum_{i=1}^{N_{gl}}\Delta\Omega_{i}\,.$
(52)

This definition works well under the assumption that the statistical
properties of the stochastic sequence $(t_{i},\Delta\Omega_{i})$ do not depend
on the window of observation, namely if the activities relative to sub-
intervals are practically constant. If $T_{obs}$ is not know and the only
information available is the sequence $(t_{i},\Delta\Omega_{i})$, a quick way
to estimate $\mathcal{A}_{a}$ would be

$\mathcal{A}_{a}\,\approx\,\frac{1}{t_{N_{\rm
gl}}-t_{1}}\left[\sum_{i=2}^{N_{gl}-1}\Delta\Omega_{i}+\frac{\Delta\Omega_{1}+\Delta\Omega_{N_{gl}}}{2}\right]\qquad\text{or}\qquad\mathcal{A}_{a}\,\approx\,\frac{N_{\rm
gl}-1}{N_{\rm gl}(t_{N_{\rm
gl}}-t_{1})}\sum_{i=1}^{N_{gl}}\Delta\Omega_{i}\,,$ (53)

in order to correct the bias induced by the fact that in the interval from
$t_{1}$ to $t_{N_{\rm gl}}$ starts and ends with the detection of a glitch.
The second expression in (53) is just a more democratic version of the first
one, where the temporal ordering of the observed sequence is completely lost
(the first and last glitch are not treated differently from the others), see
App. A of (Montoli  _et al._ , 2021).

Due to small-number statistics, this may be difficult to assess in practice,
apart for few obvious cases like the Vela and J0537-6910, where the cumulative
glitch series has a clear linear trend. For instance, too few glitches have
been detected in some pulsars, so that it may not be safe to conclude that the
inferred value of $\mathcal{A}_{a}$ corresponds to the value that would be
extrapolated if the observed sequence of glitches were longer121212 Moreover,
also having a high number of glitches $N_{gl}$ detected in a pulsar does not
necessarily guarantee that the standard linear fit procedure provides a good
estimate of $\mathcal{A}_{a}$. This can be intuitively understood by
considering a sequence in which the first (or last) glitch has an amplitude
much higher than the sum of all the others. A linear fit to the points of the
cumulative amplitude will return a value of $\mathcal{A}_{a}$ that is much
smaller than the one simply calculated via (52). Again, the problem is that in
this extreme example the activities calculated in subintervals are not
scattered around the value provided by (52). .

In any case, the absolute activity is typically inferred by performing an
ordinary linear regression on the cumulative glitch sequence, as in Fig. 3,
even though this may lead to an underestimation of the associated uncertainty
for cumulative data (Mandel, 1957; Beck, 1974; Montoli  _et al._ , 2021).

_Constraint on the moment of inertia of the pinning region -_ It is possible
to infer some constraints on the moment of inertia fraction for the angular
momentum reservoir by employing the activity parameter (Link  _et al._ ,
1999). Let us revise this argument, that is as robust as the one presented in
the previous subsection, in the sense that it does not depend on the poorly-
known details of the mutual friction in (43), but only on the general
statement (42). We formally divide $\Omega_{p}$ into the contributions due to
the smooth and slow evolution observed during the waiting-times ($R$) and the
fast one due to glitches ($G$), namely

$\dot{\Omega}_{p}=\dot{\Omega}_{p}^{R}+\dot{\Omega}_{p}^{G}\,.$ (54)

Now, we introduce an average over a long time interval of extension $T$,

$[\\![\,f\,]\\!]\,=\,T^{-1}\int_{T}dt\,f(t,...)\,$ (55)

so that, by definition,

$[\\![\,\dot{\Omega}_{p}^{G}\,]\\!]\,=\,\mathcal{A}_{a}$ (56)

Consider now the total angular momentum balance (42), and assume $I_{v}$
constant in time (possible readjustments of the stellar structure are
neglected). We have

$[\\![I_{v}\langle\partial_{t}\Omega_{vp}\rangle\,]\\!]+I\mathcal{A}_{a}+I[\\![\,\dot{\Omega}_{p}^{R}\,]\\!]=-I|\dot{\Omega}_{\infty}|$
(57)

The first term is always bounded due to the finiteness of the angular momentum
reservoir,

$[\\![I_{v}\langle\partial_{t}\Omega_{vp}\rangle\,]\\!]\,=\,\frac{I_{v}}{T}\bigg{(}\langle\Omega_{vp}\rangle_{t=T}-\langle\Omega_{vp}\rangle_{t=0}\bigg{)}<\frac{I_{v}}{T}\langle\Omega^{*}_{vp}\rangle\,=\,\frac{I}{T}\Delta\Omega_{\rm
max}$ (58)

where we have used the average in (55). In this way, if $T$ is long enough to
guarantee $\Delta\Omega_{\rm
max}/T\ll\mathcal{A}_{a}<|\dot{\Omega}_{\infty}|$, equation (57) can be
approximated as

$\mathcal{A}_{a}+[\\![\,\dot{\Omega}_{p}^{R}\,]\\!]=-|\dot{\Omega}_{\infty}|\,.$
(59)

We can use this expression to find a theoretical upper bound on the observed
activity. To do so we need to set a robust limiting value to
$\dot{\Omega}_{p}^{R}$. As discussed in Sec. IV, we can invoke the extreme
scenario in which the vortices are pinned between glitches
($\dot{\Omega}_{v}^{R}=0$). This idealized situation is probably never
realized in the real pulsar because of both thermal creep (Alpar  _et al._ ,
1984a) and disorder in the pinning landscape (Antonelli and Haskell, 2020), so
we can safely conclude that

$[\\![\,\dot{\Omega}_{p}^{R}\,]\\!]\,<\,-\,\frac{I\,|\dot{\Omega}_{\infty}|}{I-I_{v}}\,.$
(60)

Finally, plugging (59) in the inequality above gives the sought upper bound on
the observed activity,

$\frac{\mathcal{A}_{a}}{|\dot{\Omega}_{\infty}|}<\frac{I_{v}}{I-I_{v}}\,.$
(61)

Formally, the above constraint reduces to the one derived by Link  _et al._
(1999) for $\epsilon_{n}=0$. The entrainment corrections to the moments of
inertia were first derived by Chamel and Carter (2006): their result, albeit
derived by using a different – but equivalent – Newtonian formalism, is
consistent with (61). Relativistic corrections to $I_{v}$ can also be included
(Antonelli  _et al._ , 2018; Montoli  _et al._ , 2021).

The constraint (61) is usually derived by starting from a model with two rigid
components rotating with uniform angular velocities $\Omega_{p}$ and
$\Omega_{n}$. Within this framework and in the limit of small $I_{v}$, the
constraint turns out to be (Chamel and Carter, 2006; Andersson  _et al._ ,
2012; Chamel, 2013)

$\frac{\mathcal{A}_{a}}{|\dot{\Omega}_{\infty}|}<\frac{I_{n}}{(1-\langle\epsilon_{n}\rangle)\,I}\qquad\quad
I_{n}=\lim_{\epsilon_{n}\rightarrow 0}I_{v}\,,$ (62)

where $\langle\epsilon_{n}\rangle$ is the averaged entrainment in the
superfluid domain, see also (Sidery  _et al._ , 2010).

If we assume the pinned superfluid region to be limited in the crust of the
star, and we fix the microphysical parameters, namely the EOS and the
entrainment parameter, then the moments of inertia $I_{v}$ and $I$ in (61)
depend only on the stellar mass. Therefore, similarly to the case for maximum
glitch amplitude, also here we obtain a constraint on the mass of the pulsar
(see Fig. 6). Apart from the EOS, and in both cases the constraint depends on
a particular microphysical parameter: the maximum pinning force for the
constraint from the largest glitch and the entrainment for the activity
constraint. Thus, the two different methods would allow for two tests on
completely different microphysical parameters, if the real mass of the star
were known.

Figure 6:  The curve defining the activity constraint (61) for some EOSs.
Here, the expression of $I_{v}$ is the general relativistic one in the slow
rotation approximation, as defined in (Antonelli  _et al._ , 2018; Montoli
_et al._ , 2021). The superfluid reservoir is limited to the inner and the
entrainment is that of Chamel (2012). The ratio
$\mathcal{A}_{a}/|\dot{\Omega}_{\infty}|$ is plotted for the Vela pulsar,
along with the uncertainty as obtained via the statistical “bootstrap”
procedure described in (Montoli  _et al._ , 2021).

The entrainment parameter for the superflow in the crust of a neutron star has
been obtained in (Chamel, 2012), by calculating the effects of Bragg
scattering due to the presence of the crustal lattice. The main result is a
negative entrainment parameter $\epsilon_{n}\sim-10$ in the crust of the star,
which implies a severely reduced superfluid density $\rho_{n}^{*}$ and a
smaller $I_{v}$, making the requirement in (61) more difficult to be met. In
fact, studies on the activity parameter for the Vela yield an unrealistically
low constraint on the mass of the pulsar for many EOSs, if a crustal reservoir
and the entrainment calculated by Chamel (2012) are assumed (see Andersson
_et al._ , 2012; Chamel, 2013; Delsate  _et al._ , 2016). Several ways have
been proposed to overcome this impasse, including the use of stiff EOSs
(Piekarewicz  _et al._ , 2014), a Bayesian analysis of the EOS uncertainty
(Carreau  _et al._ , 2019) or an extension of the neutron superfluid
participating in the glitch beyond the crust-core transition (Ho  _et al._ ,
2015; Montoli  _et al._ , 2020b).

Furthermore, the exact values of the entrainment parameter in the inner crust
are still matter of debate. Apart from calculation based on the effects of
band structure on the effective neutron superfluid density in the crust
(Carter  _et al._ , 2005; Chamel, 2012), other alternative estimates have been
proposed, which include hydrodynamical approaches (Martin and Urban, 2016) or
Bragg scattering in a disordered crustal lattice (Sauls  _et al._ , 2020),
both of which yield milder entrainment effects in the crust. According to the
analysis of Watanabe and Pethick (2017), the effects of band structure – that
should lead to strong entrainment and to a smaller value of $I_{v}$ – are
suppressed when the pairing gap is of order or greater than the strength of
the crustal lattice potential – see the discussion in Chamel (2017b). In all
these cases, the reduction of the neutron superfluid density would be less
than estimated on the basis of band theory calculations and, consequently, it
may be premature to rule out models of glitches based on superfluidity in the
crust alone.

Finally, different ways of estimating the activity parameter yield different
values for the uncertainty of this parameter. For example, dropping the
homoscedasticity hypothesis for the cumulative glitch amplitude data that
define the glitch activity gives rise to larger uncertainties, see Fig. 6 and
(Montoli  _et al._ , 2021).

###  V.3 Stronger constraints

As soon as we tentatively extend the pinning region partially in the inner
core, the moment of inertia $I_{v}$ grows so much that the curves in Fig. 6
are shifted upwards in such a way that the mass of the Vela is consistent with
any value smaller than the maximum one allowed by the EOS. However,
observations of both the maximum glitch amplitude and the glitch activity may
be used in tandem to put stronger – but less robust – constraints on the
pulsar’s mass (Montoli  _et al._ , 2020b). In fact, while the constraints
coming from the maximum glitch amplitude and the activity are derived on the
basis of the angular momentum balance only, we now have to model the internal
dynamics, which amounts to select a particular mutual friction prescription.

Starting from corotation of the two components, forward integration in time of
the dynamical equations (42) and (43) allows to follow the evolution of the
lag $\Omega_{vp}(x,z,t)$ and of the associated quantity $\Delta\Omega_{m}(t)$
defined as

$\Delta\Omega_{m}(t)\,=\,\frac{I_{v}}{I}\langle\,\Omega_{vp}\,\rangle\,<\,\Delta\Omega_{\rm
max}\,,$ (63)

that is nothing but a theoretical upper bound on the amplitude of a glitch
triggered at time $t$ (assuming that at $t=0$ the star was in a state of
corotation). In other words, $\Delta\Omega_{m}(t)$ sets a theoretical upper
limit for the glitch amplitude at time $t$, given that a large glitch occurred
at $t=0$. Again, since this is just an upper bound, it is not really important
if glitches in real pulsars can empty the whole reservoir: this is just an
assumption to construct the theoretical upper bound $\Delta\Omega_{m}(t)$.

The information coming from the glitch activity is used to roughly estimate
the time $t$ needed to replenish the angular momentum reservoir to a level
that is sufficient for the largest glitch amplitude $\Delta\Omega_{\rm obs}$
observed in the pulsar under study, $t_{act}\,\sim\,\Delta\Omega_{\rm
obs}/\mathcal{A}_{a}$. Asking that
$\Delta\Omega_{m}(t_{act})>\Delta\Omega_{\rm obs}$ provides the sough
constraint on the pulsar’s mass, once the microphysical input and the mutual
friction prescriptions are fixed.

An alternative and interesting possibility is to adapt the activity constraint
in (61) by taking into account the thermal history of a pulsar. Ho  _et al._
(2015) used (62) to constrain the mass of a sample of active pulsars by
considering a superfluid reservoir that extends in the outer core, accounting
for the temperature dependence of superfluidity. They inferred the temperature
of pulsars from observations to constrain how much the superfluid reservoir
extends into the core. This can be done by choosing a model of critical
temperature of neutron (singlet-state) superfluidity as a function of baryon
number density (Lombardo and Schulze, 2001) and simulating the internal
temperature of a pulsar. If the age of the pulsar is known – or, at least, its
characteristic age $P/(2\dot{P})$ – this gives an estimate of where the
internal temperature of the cooling core is lower than the critical
temperature for superfluidity. Hence, this approach makes it possible to
obtain mass upper bounds by combining pulsar glitch data and the temperature
dependence of superfluidity131313 The cooling and rotational evolution in a NS
are coupled. For example, the late-time thermal evolution is affected by the
spin-down due to the deposition of heat due to the action of mutual friction:
besides the neutrino and photon emissions that are energy sinks, rotational
energy is converted into heat by friction due to differential rotation of the
superfluid (Page, 1997; Larson and Link, 1999). Conversely, the internal
temperature also regulates the thermal creep (Alpar  _et al._ , 1984a) and,
possibly, localized heat deposition may increase the creep rate to initiate a
glitch (Link and Epstein, 1996; Larson and Link, 2002). .

Differently from (63), this method does not assume anything about the mutual
friction, yet some extra information is needed. Instead of simulating the
internal rotational dynamics, the difficulty is to simulate the cooling, that
depends on several ingredients, most importantly the superfluid gaps (Page
_et al._ , 2006). Moreover, the characteristic age is not always an accurate
indication of the true age of a pulsar.

In the end, current mass estimates of glitching pulsar like the ones presented
in (Ho  _et al._ , 2015) and (Pizzochero  _et al._ , 2017) should be taken
more as an indication of how much is the influence of the assumed microscopic
parameters on the final results.

##  VI Dynamical constraints from spin-up observations

We already discussed in Sec. II.1 how observations of the slow relaxation
after a glitch is the main evidence for the presence of superfluid components
in NSs. Also, the short spin-up timescale associated to the 2000 and 2004
glitches in the Vela pulsar – of the order of less than a minute (Dodson  _et
al._ , 2002, 2007) – and the particular slow-rise detected in the Crab pulsar
(Shaw  _et al._ , 2018) can be interpreted as hints for the fact that the
mutual friction is not simply linear, but that, on the contrary, its strength
is modulated as the superfluid evolves across the various dynamical phases
(Haskell  _et al._ , 2018; Gügercinoǧlu and Alpar, 2019; Celora  _et al._ ,
2020). Despite this is the general picture at the basis of the vortex-mediated
friction, we anticipate that the mutual friction is typically taken to be
linear – as in the system (65) below – as long as the model is supposed to be
applied only to a single dynamical phase (see, e.g., Haskell  _et al._ , 2012;
Graber  _et al._ , 2018; Sourie and Chamel, 2020; Pizzochero  _et al._ ,
2020).

The recent detection of a glitch in the Vela pulsar in December 2016
(Palfreyman  _et al._ , 2018) has opened a window on the elusive spin-up
phase. The peculiarity of this detection lies in the fact that this is the
first pulse-to-pulse observation of a glitch in the act, allowing for
unprecedented timing resolution of the spin-up and the following phases (see
Sec. IV). In fact, analysis of the phase residuals revealed that the initial
fast spin-up occurs over a timescale $\tau_{exp}\lesssim 10\,$s, as well as
the presence of an _overshoot_ in the data in the first seconds of the glitch
(Ashton  _et al._ , 2019; Pizzochero  _et al._ , 2020), a transient dynamical
phase sometimes observed in glitch simulations that account for the non-
uniformity of the superfluid velocity lag and stellar stratification (Larson
and Link, 2002; Haskell  _et al._ , 2012; Antonelli and Pizzochero, 2017;
Graber  _et al._ , 2018).

It is possible to use the minimal Baym’s model in (5) to fit the data of
Palfreyman  _et al._ (2018) and extract a coupling timescale that, in turn,
may be used to infer the mutual friction strength $\mathcal{B}_{d}$ during the
spin-up. Taking the spin-up timescale to be $\tau\approx\tau_{exp}\sim 10\,$s
and assuming that the glitch is due to the superfluid in the crust – the
superfluid in the core is assumed to be corotating with the normal (Alpar  _et
al._ , 1984b), so that $I_{v}/I\sim 1\%$ – from (46) we conclude that141414
The estimate $\mathcal{B}_{d}\sim 10^{-5}$ given in equation (2) of Ashton
_et al._ (2019) is based on the assumption of a single superfluid component
that is limited to the inner crust ($I_{v}\sim 0.01I$). Since they do not
explicitly include entrainment, it should be reproduced by our
$\langle\epsilon_{n}\rangle\ll 1$ case. The two orders of magnitude
discrepancy is due to a missing $I_{v}/I$ factor.

$\begin{split}\langle\mathcal{B}_{d}\rangle\sim
10^{-2}\qquad&\text{if}\quad\langle\epsilon_{n}\rangle\sim 10\\\
\langle\mathcal{B}_{d}\rangle\sim
10^{-3}\qquad&\text{if}\quad\langle\epsilon_{n}\rangle\ll 1\end{split}$ (64)

in the inner crust crust, where we recall that $2\Omega_{p}\sim 100\,$rad/s
for the Vela. The two cases refer, roughly, to the strong (Chamel, 2012,
2017b) and weak (Martin and Urban, 2016; Sauls  _et al._ , 2020) entrainment
cases, respectively.

The rough estimate in (64) depends on the fact that we reduced a fluid model
to a much simpler simple model with only two rigid components. Clearly, such a
simplification cannot be valid if, for example, the superfluid in the core is
coupled to the normal component on a timescale that is comparable or longer to
the one of the superfluid in the crust. Moreover, it is not possible to
reproduce the overshoot feature since the response of the system in (5) can
only fit an exponential rise.

###  VI.1 Models with three rigid components

The minimal glitch model which allows for an overshoot requires at least three
rigid components. This can be formally obtained from the two-component fluid
equations (42) and (43) by splitting the whole superfluid domain in two non-
overlapping regions (Montoli  _et al._ , 2020a). Taking the average over these
two regions – see the discussion below (46) – gives us the three-component
analogue of the Baym’s model,

$\displaystyle\begin{split}&x_{p}\dot{\Omega}_{p}+x_{1}\dot{\Omega}_{1}+x_{2}\dot{\Omega}_{2}\,=\,-|\dot{\Omega}_{\infty}|\\\
&\dot{\Omega}_{1}=-b_{1}\left(\Omega_{1}-\Omega_{p}\right)\\\
&\dot{\Omega}_{2}=-b_{2}\left(\Omega_{2}-\Omega_{p}\right)\end{split}$ (65)

where $p$ indicates the normal (observable) component, $i=1,2$ labels the two
superfluid regions, $x_{i}=I_{i}/I$ is the moment of inertia fraction of the
region $i$ and $b_{i}$ is the inverse of the coupling timescale between the
$i$ region and the normal component. The normal fraction $x_{p}=1-x_{1}-x_{2}$
is not an independent parameter of the model. The system (65) can be solved
analytically. Following the same notation used in (5), the angular velocity
residue with respect to the steady-state spin-down law reads

$\Delta\Omega_{p}(t)\,=\,\Delta\Omega_{p}^{\infty}\left[1-\omega\,e^{-t\lambda_{+}}-(1-\omega)\,e^{-t\lambda_{-}}\right]\,,$
(66)

where $\omega$, $\Delta\Omega_{p}^{\infty}$ and $\lambda_{\pm}$ are functions
of the parameters $x_{1,2}$, $b_{1,2}$ and $|\dot{\Omega}_{\infty}|$, as well
as of the initial conditions for the two superfluids and the normal
components. In particular, $\Delta\Omega_{p}^{\infty}$ is the angular velocity
residue of the star at large $t$, after both exponential transients have fully
relaxed.

To extract information on the parameters $x_{1,2}$ and $b_{1,2}$, the formula
(66) must be supplemented with the exact dependence of its four parameters
$\omega$, $\Delta\Omega_{p}^{\infty}$ and $\lambda_{\pm}$ on the physical
parameters $x_{1,2}$, $b_{1,2}$ and the initial conditions. For example, the
two relaxation timescales $1/\lambda_{\pm}$ in (66) are functions of all the
four parameters $b_{1,2}$ and $x_{1,2}$,

$\lambda_{\pm}=\frac{1}{2}\bigg{[}\beta_{1}(1-x_{2})+\beta_{2}(1-x_{1})\pm\sqrt{[\beta_{1}(1-x_{2})+\beta_{2}(1-x_{1})]^{2}-4\beta_{1}\beta_{2}x_{p}}\bigg{]}\,,\quad\beta_{i}=b_{i}/x_{p}\,.$
(67)

The explicit formulae for $\omega$ and $\Delta\Omega_{p}^{\infty}$, as well as
the link of the averaged physical parameters $b_{1,2}$ and $x_{1,2}$ with the
local microscopic input (including entrainment), are given in (Montoli  _et
al._ , 2020a).

_Physical interpretation of the parameters -_ Exactly as the Baym’s model,
also the system (65) is agnostic about the location of the superfluid
components, so that the moment of inertia fractions $x_{1,2}$ of the two
superfluid components do not necessarily correspond to crust and core
superfluid. For example, the two superfluid components may be both located in
the core of the star (Sourie and Chamel, 2020). The hope is to fit (66) from
observations of a well-resolved spin-up, and then to infer the most likely
values for $x_{1,2}$ and $b_{1,2}$, so to constrain the superfluid regions and
the (averaged) friction parameters.

In the original Baym’s model defined in (1) and (2) the entrainment coupling
was not explicit, but hidden in the expressions of the phenomenological
parameters in terms of the local ones. We can reduce the full fluid problem to
the set of phenomenological ordinary differential equations in (65) by
following the same reasoning that led to (46). Now, instead of introducing a
single average over the whole superfluid region, we formally split it in two
non-overlapping domains (labelled by $i=1,2,$),

$dI_{v}=d^{3}x\,x^{2}\,\rho^{*}_{n}\qquad\quad
I_{i}=\int_{i}dI_{v}\qquad\quad\langle\,f\,\rangle=I_{i}^{-1}\\!\\!\int_{i}dI_{v}\,f\,$
(68)

In this way, the inverse coupling timescales and the moments of inertia
fractions can be expressed as (compare with equation 46)

$\begin{split}&x_{i}=I_{i}/I\qquad\qquad\,\,x_{p}=(I-I_{1}-I_{2})/I\\\
&\Omega_{1}=\langle\Omega_{v}\rangle_{1}\qquad\qquad\Omega_{2}=\langle\Omega_{v}\rangle_{2}\\\
&b_{i}\approx\,2\Omega_{p}\left\langle\,\frac{\mathcal{B}_{d}\,}{1-\epsilon_{n}}\right\rangle_{i}\approx\,2\Omega_{p}\frac{\langle\,\mathcal{B}_{d}\,\rangle_{i}}{\langle\,1-\epsilon_{n}\,\rangle_{i}}\end{split}$
(69)

Again, the presence of entrainment affects the interpretation of the
macroscopic parameters $x_{i}$ and $b_{i}$ in terms of the local input.

It is worth to remark that the idea behind the reduction to a model with rigid
components is that the correlations in (45) are negligible. Hence, it may be a
good strategy to identify the two regions in such a way that the resulting
parameters $b_{1,2}$ are as different as possible, so that the superfluid tend
to respond in a as uniform as possible way within a region. Because of this,
it is reasonable to assume that the two superfluid components represent
regions with different physical properties (like the superfluid in the crust
and the one in the inner core, but we will see that this is not the only
possibility).

Finally, the identities in (69) also tell us that the description of the
dynamical phases in Sec. IV is still valid, since

$I_{v}\,\langle\partial_{t}\Omega_{v}\rangle\,=\,I(x_{1}\dot{\Omega}_{1}+x_{2}\dot{\Omega}_{1})\qquad\qquad
I_{v}=I(x_{1}+x_{2})\,.$ (70)

In fact, the averaging procedure used to reduce a fluid model to a model with
two or more rigid components reproduces perfectly the the total angular
momentum balance (42).

###  VI.2 Overshoot condition

The maximum value $\Delta\Omega_{\rm over}$ of the angular velocity jump
during a glitch can be obtained by finding the maximum value that can be
attained in (66). There are two possible situations: the maximum of
$\Delta\Omega_{p}(t)$ is reached asymptotically – namely $\Delta\Omega_{\rm
over}=\Delta\Omega_{p}^{\infty}$ at large times – or there exists a finite
time for which $\Delta\Omega_{\rm over}>\Delta\Omega_{p}^{\infty}$. In the
latter case, the solution $\Delta\Omega_{p}(t)$ describes an overshoot and
this happens if $\omega>1$. However, the parameter $\omega$ is, in general, a
rather complex function of the initial lags and the physical parameters, so
that this overshoot condition is not of simple physical interpretation.

It is possible to simplify the problem by restring our attention to a
particular subset of solutions. The simplest overshoot condition can be found
when the unknown initial conditions rely on some physical assumptions on the
role of the two components (Pizzochero  _et al._ , 2020). One of the two
superfluid components, say 1, is assumed to be a “passive” one that does not
change its creep rate – meaning that $b_{1}$ is constant before and after the
triggering of the glitch. Hence, it is reasonable to set its initial lag to
the steady-state value, $\Omega_{1}-\Omega_{p}=|\dot{\Omega}_{\infty}|/b_{1}$
at $t=0$. On the other hand, the component 2 can be assumed to be the “active”
one, namely the component that undergoes a transition from a pinned to an
unpinned state. Since this is the component that is pinned before the glitch,
it is also the component that provides the extra angular momentum for the
observed spin-up.

In this setting, it is possible to show that the mathematical condition for an
overshoot is just $b_{1}<b_{2}$: the post-glitch timescale $1/b_{2}$ of the
active component (that in the pre-glitch state was only loosely coupled to the
rest of the star because of pinning) must be smaller than the coupling
timescale $1/b_{1}$ of the passive component (Pizzochero  _et al._ , 2020).
From the physical point of view, the overshoot occurs when the superfluid
region that stores the angular momentum for the glitch (the active component)
can transfer its excess of angular momentum to the normal component faster
than the typical timescale the rest of the superfluid reacts with. This idea
is further clarified by considering the limit $b_{1}\ll b_{2}$: the passive
component maintains its state of motion, being practically decoupled, while
the active superfluid imparts its momentum to the normal component only,
spinning it up. At this point, the lag between the passive superfluid and the
normal component is smaller than its positive steady-state value or even
possibly reversed. The post-overshoot relaxation, during which the lag
$\Omega_{1}-\Omega_{p}$ will tend to approach again the steady-state value, is
modulated by the slower response of the passive superfluid. This kind of
behaviour is also clear from the simulations of the Vela 2016 glitch performed
in (Graber  _et al._ , 2018), where the passive and active components are
assumed to be the superfluid in the core and the one in the inner crust,
respectively.

Sourie and Chamel (2020) applied the above picture to study the possible
impact of the pinning to the flux-tubes in the outer core, where the magnetic
field is expected to be predominantly toroidal and pinning to fluxoids may be
the most effective (Gügercinoglu and Alpar, 2014). In their setting the active
component is the superfluid in the outer core, while the passive one is the
superfluid in both the inner core and crust. In this way, they have been able
to translate the simplified overshoot condition $b_{1}<b_{2}$ into a statement
on the number of fluxoids that are, on average, pinned to a vortex. When the
number of pinned fluxoids exceeds a certain threshold, the parameter $b_{2}$
is suppressed (pinning tends to suppress mutual friction) and no overshoot is
possible.

###  VI.3 Constraint on the moment of inertia of the passive component

By using a three-component model of the kind in (65) and its particular
solution given in (Pizzochero  _et al._ , 2020), Sourie and Chamel (2020)
proposed a way to constrain the moment of inertia fraction $x_{1}$ of the
passive component from the observation of an overshoot in the angular velocity
residue. This result has been generalized in (Montoli  _et al._ , 2020b),
where the constraint is derived directly from the full solution to the problem
including entrainment corrections.

We have seen that an overshoot is likely to occur when $b_{1}<b_{2}$. By
expanding the value of $\Delta\Omega_{\rm over}$ in the limit $b_{1}\ll
b_{2}$, it is possible to show that

$x_{1}\,>\,1-\Delta\Omega_{p}^{\infty}/\Delta\Omega_{\rm
over}\qquad\quad\textrm{if}\quad b_{1}\ll b_{2}$ (71)

Hence, the a measurement of both the asymptotic amplitude
$\Delta\Omega_{p}^{\infty}$ and of the maximum amplitude at the peak of the
overshoot $\Delta\Omega_{\rm over}$ provides a lower bound to the superfluid
fraction $x_{1}$. For example, observing $\Delta\Omega_{\rm over}\approx
1.1\times\Delta\Omega_{p}^{\infty}$ would immediately imply $x_{1}>10\%$,
ruling out the possibility that the passive component is the inner crust
alone.

Both $\Delta\Omega_{p}^{\infty}$ and $\Delta\Omega_{\rm over}$ may be obtained
by fitting the data of a spin-up observation: we try to test the constraint
(71) with the data of the Vela 2016 glitch in Fig. 7.

Figure 7:  Corner plot for the posterior distribution of the quantities
$x_{1}$ and $1-\Delta\Omega_{p}^{\infty}/\Delta\Omega_{\rm over}$ obtained
from the Bayesian analysis of the Vela 2016 spin-up reported in (Montoli  _et
al._ , 2020a). The orange diagonal represents the constraint in equation (71)
of Sourie and Chamel (2020): most of the posterior points satisfy the
constraint, indicating the likely the presence of a glitch overshoot in the
data.

The median of the probability distribution of $x_{1}$ is larger than that of
$1-\Delta\Omega_{p}^{\infty}/\Delta\Omega_{\rm over}$. In the covariance plot
most of the points fall in the region below the diagonal, thus respecting the
constraint in (71). Since (71) is valid in the limit $b_{1}\ll b_{2}$ (but the
analysis in Fig. 7 is not), this also is a clue that it is likely that in the
2016 Vela glitch the two coupling parameters $b_{1,2}$ have values different
by at least an order of magnitude. In fact, posteriors for $b_{1,2}$ reported
by Montoli  _et al._ (2020a) indicate that the 16th-84th percentile intervals
are $b_{1}\sim 0.004-0.009\,$s-1 and $b_{2}\sim 0.08-24.64\,$s-1.

##  VII Conclusions

We discussed some of the fundamental ideas behind the current understanding
and modelling of pulsar glitches. Since the phenomenological minimal model of
Baym  _et al._ (1969) was developed, we now have now a better understanding of
how the hydrodynamic modelling of neutron star interiors works. While
phenomenological models are useful to fit observations, information on the
physics of the interior of a NS can be extracted only if we can reliably link
the parameters of the fits to a more fundamental physical description of the
system.

As a pedagogical exercise, we derived – starting from a local fluid
description that is more directly linked to the microphysics of the system –
the axially symmetric equations that constitute the backbone of several glitch
models. In this way, the phenomenological parameters of Baym’s model (or its
natural three-component generalization) can be expressed in terms of the local
properties of dense nuclear matter.

We have revised the current state of the art in constraining some aspects of
neutron star structure from timing observations of pulsar glitches. We
discussed two robust constraints coming, respectively, from the observation of
glitches of large amplitude (Antonelli  _et al._ , 2018) and from the average
glitch activity over an extended period of time (Link  _et al._ , 1999;
Chamel, 2013; Andersson  _et al._ , 2013; Delsate  _et al._ , 2016;
Piekarewicz  _et al._ , 2014; Carreau  _et al._ , 2019; Montoli  _et al._ ,
2021). Observation of large glitch amplitudes (the record holder is
J0205+6449) can be used to test microscopic estimates of pinning forces, while
the high activity of some pulsars (the Vela being the most active) provides a
way to test the entrainment coupling. The two tests are completely independent
(apart from the fact that both rely on the choice of an EOS) and seem to
indicate that part of the superfluid in the outer core should participate in
storing angular momentum: the upper limit on the inferred mass coming from the
constraints is, for some objects, close or smaller than the theoretical
minimum NS mass of $1.17\,M_{\odot}$ inferred from supernovae explosion
simulations (Suwa  _et al._ , 2018).

This is interesting, as pinning of vortices to flux tubes is actually a
possible mechanism to store angular momentum, as reviewed by Alpar (2017).
Those two static constraints may always be refined, but at the expense of
solving the hydrodynamic equations in the NS interiors, or invoking some
simplified model that accounts for the angular momentum exchange between the
components (Montoli  _et al._ , 2020b).

Then, we moved to consider a class of dynamical studies (as opposed to the two
static constraints from glitch amplitudes and activities). In particular, the
possibility of observe all the phases of a glitch in the act (Palfreyman  _et
al._ , 2018) justifies the efforts towards a better understanding and less
phenomenological modelling of the processes that mediate the angular momentum
exchange in a glitch. For instance, our ability to extract information on NS’s
structure and superfluid properties from glitch observations is directly