{'fraction_non_alphanumeric': 0.03129016981259725, 'fraction_numerical': 0.011602733238617144, 'mean_word_length': 4.50512104283054, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Axially symmetric equations for differential pulsar rotation with superfluid entrainment

In this article we present an analytical two-component model for pulsar rotational dynamics. Under the assumption of axial symmetry, implemented by a paraxial array of straight vortices that thread the entire neutron superfluid, we are able to project exactly the 3D hydrodynamical problem to a 1D cylindrical one. In the presence of density dependent entrainment the superfluid rotation is non-columnar: we circumvent this by using an auxiliary dynamical variable directly related to the areal density of vortices. The main result is a system of differential equations that take consistently into account the stratified spherical structure of the star, the dynamical effects of non-uniform entrainment, the differential rotation of the superfluid component and its coupling to the normal crust. These equations represent a mathematical framework in which to test quantitatively the macroscopic consequences of the presence of a stable vortex array, a working hypothesis widely used in glitch models. Even without solving the equations explicitly, we are able to draw some general quantitative conclusions; in particular, we show that the reservoir of angular momentum (corresponding to recent values of the pinning forces) is enough to reproduce the largest glitch observed in the Vela pulsar, provided its mass is not too large.


INTRODUCTION
The detection of glitch events in many isolated radio pulsars indicates that a large amount of angular momentum must be first stored and then exchanged between an internal component and the rigid observable crust. Thus a neutron star has to be comprised of (at least) two components, spinning at slightly different velocities, that interact via some type of mutual friction in order to permit such an accumulation and exchange of angular momentum. Observations of very long post glitch relaxations indicate that one of these components should be superfluid, whereas the spin-up is so fast that it is still unresolved (less than one minute).
With this in mind, these recurring period instabilities are usually explained in terms of superfluid vortex dynamics (Anderson et al. 1975), a scenario where pinning of vortex lines to the crustal lattice freezes part of the angular momentum and glitches are the consequence of a vortex avalanche, a catastrophic event in which many vortices are suddenly expelled from the superfluid bulk and transfer their angular momentum to the crust. The nature of the event that triggers an avalanche is still debated: it could be a starquake E-mail: marco.antonelli@unimi.it (Ruderman 1976;Ruderman et al. 1998), a hydrodynamic instability (Glampedakis & Andersson 2009) or a manifestation of a self-organized critical system (Melatos et. al. 2008). Despite the variety of existing studies, however, glitch theory lacks consistent macroscopic models whose quantitative outputs can be compared with observational data: for example the need to take into account the spherical and non-uniform structure of the star in order to evaluate the angular momentum reservoir was pointed out only recently in the framework of the static "snowplow model" (Pizzochero 2011). Since the pioneering two-component model of Baym et al. (1969), detailed dynamical models have evolved in order to incorporate in the dynamical equations different subtle effects of multi-fluid hydrodynamics, such as superfluid entrainment (Carter 1989;Comer & Langlois 1994;Chamel 2008). Although the local hydrodynamical equations are well known, the great majority of studies of glitches have been carried out in the quite unphysical approximation of a cylindrical, uniform-density and rigidly rotating star: as pointed out by Andersson et al. (2012), there has been little progress on making the pulsar glitch models quantitative.
In this work we address the problem of how to construct a simple but fully consistent model for pulsar glitches. By consistent we mean that, starting from well stated hypothe-sis, we provide a coherent prescription of how to implement all the available microphysical information in our macroscopic set of equations. This permits to account for the fact that radially the star has a non-uniform structure which in turn implies that the rotation of the superfluid must be differential, even in the simplest situation of parallel and straight vortices.
Our study is in the framework of the two-component model: a charged, rigid p-component is spun down electromagnetically, while a non-observable n-component comprised of superfluid neutrons reacts to the spin-down via some type of mutual friction. The hydrodynamical formalization of this multi-fluid problem is nowadays quite clear and developed but leads to a computationally difficult 3D problem [see for example  and Haskell & Melatos (2015) for recent reviews and Peralta et al. (2006) for hydrodynamical simulations in Couette flow geometry]. We can reduce it to a simpler 1D model by making two widely used assumptions: axial symmetry and straight vortices. On the one hand, the assumption of axial symmetry around the rotational axis of the pulsar is a widespread working hypothesis and allows to reduce the dimensionality of the problem. In this way, however, we loose the possibility to study the nature of the trigger event, since probably it is related to some non-trivial effect of the fully three-dimensional problem. On the other hand, the assumption of a paraxial array of straight vortices which resist bending is justified by the collective rigidity of vortex bundles in coherent motion: it is consistent with and required by the assumption of axisymmetry in a rotating superfluid, as discussed in a seminal paper by Ruderman & Sutherland (1974). This allows explicit integration along the vortex lines, projecting the spherical star onto its equatorial plane. In this way, however, we rule out a priori turbulent motion and its consequences on the mutual friction between the two components ): this is a limit of the model to be kept in mind.
In this paper we develop a dynamical model for superfluid differential rotation with entrainment: this is a new scenario since in the literature we can only find rigid models with uniform entrainment [as for example Sidery et al. (2010)] or differential models without entrainment [as done in Haskell et al. (2012)]. The recent large values found for entrainment in the crust indicate that such a phenomenon cannot be ignored in the physics of pulsar glitches (Chamel 2012(Chamel , 2013Andersson et al. 2012). Moreover, entrainment turns out to be strongly density-dependent, as shown in section 4. Thus an important part of our work is devoted to dealing with the complication that in the presence of density-dependent entrainment the hydrodynamics cannot be columnar, even in the extremely simple case of straight and parallel vortices: because of non-uniform entrainment, the angular velocity of the n-component will depend on the coordinate along the axis of rotation.
In the next section we will describe our working hypothesis and how to handle the entrainment effect when there are macroscopically long vortices that pass trough spherical layers of different density and composition. In the absence of a layer of normal matter between the core and the inner crust, the normal vortex cores should pass continuously through the phase transition between the S-wave superfluid in the crust and the P-wave superfluid in the core (Zhou et al. 2004); therefore, we take vortices which thread the entire neutron superfluid. We stress, however, that this is only a possibility and that it is straightforward to generalize our model to account for two different superfluid components (as long as we keep the assumption of straight vortices). The main result of the paper are the equations of motion of the system presented at the beginning of section 3; they are written in a form that does not depend explicitly on the presence of entrainment or by the fact that vortices pass trough the crust-core interface. These dynamical equations generalize and extend the previous work of Haskell et al. (2012) to the case of density-dependent superfluid entrainment and (with a suitable choice of the mutual friction term) provide a dynamical realization of the static "snowplow model" introduced by Pizzochero (2011). Thus the present work sets a rigorous mathematical framework in which to test quantitatively the macroscopic consequences of the presence of a stable array of parallel vortices, a working assumption widely used in almost all existing glitch models. In the last two sections, we give a numerical example of our results, using consistently some of the microphysical information that can be found in the literature and employing two different equations of state to describe dense matter. This also permits us to give estimates of the superfluid angular momentum reservoir and of the spin-up timescales as a function of the stellar mass. A concluding section summarizes the physical relevance of the paper for future numerical studies of pulsar glitches.

DESCRIPTION OF THE HYDRODYNAMICAL MODEL
In this section we derive the dynamical equations needed to describe the differential rotation of a pulsar in the framework of the two component model with superfluid entrainment. As already discussed in the introduction, in the presence of density-dependent entrainment the superfluid velocity field cannot be columnar, even if boundary effects and turbulence are totally neglected. To overcome this, we will not consider directly the velocity of the superfluid as a dynamical variable, but rather an auxiliary variable directly associated to the vortex line density.

Model assumptions
In the following we describe the set of assumptions, some of which can be generalized at a later time, that are used to construct our model: (i) The star is spherical, spinning around the z-axis (namely there is no precession), with inner-to-outer crust interface at R d (the radius at which the neutron drip starts). In order to simplify the fully 3D hydrodynamical problem, we will use the hypothesis of axial symmetry around the rotation axis of the star. To implement this symmetry, the slow-down of the pulsar is given by an external braking torque with constant direction Text = −Textêz that acts on the star.
(ii) We consider only one superfluid component (the ncomponent) inside the spherical volume of radius R d . This superfluid component is threaded by straight vortices parallel toêz (as in Fig 1), each carrying a quantum of circulation κ = h/(2mn) (with mn the nucleon mass). We discuss further this point at the end of this section.
(iii) In the absence of a layer of normal matter between the core and the inner crust (Zhou et al. 2004) there is no clear distinction between the crustal superfluid and the superfluid in the core and we assume that normal matter vortex cores pass continuously through the phase transition between the S-superfluid to the P-superfluid. This is the simplest possibility but it can be eventually modified by considering superfluidity only inside a spherical shell of radii Rc and R d . Although the scenario of continuous vortex lines was already suggested by Ruderman (1976) as an alternative to the scenario of distinct vorticity in the crust and the core, in the subsequent literature only the latter assumption has been implemented until the "snowplow model" of Pizzochero (2011).
(iv) The rigid crust and the charged component (electron and proton fluids) are locked into the strong magnetic field of the neutron star on very short timescales (the p-component). This is justified as long as protons (or ions in the crust) and electrons are locked on timescales much shorter than the dynamical timescales we are interested in (Easson 1979). Using cylindrical coordinates (x, θ, z), the velocity field of this rigid p-component is vp = xΩpê θ , which is related to the observed pulsar period P = 2πΩ −1 p . (v) As usual in the two-fluid approach, the presence of vortices induces a mutual friction between the n-component and the p-component, due to the presence of internal dissipative channels that can be modeled as drag forces acting on vortex lines (Andersson et al. 2006). For each dissipative process ("d.p.") occurring at a given density ρ, we need a corresponding drag parameter η d.p. (ρ) that gives the intensity of the associated drag force on vortices.
(vi) Vortex lines can also pin to the different inhomogeneities present in the star, thus we need a quantity that (for different densities and pinning mechanisms) tells us to which extent the vortices can be blocked. We assume that pinning can be modeled by introducing a "pinning force per unit length (of vortex line)" fp(ρ), which expresses the strength of the vortex-inhomogeneity interaction, suitably mediated over a mesoscopic region at density ρ.
Given this set of assumptions, we are far from describing the real hydrodynamical problem of a spinning neutron star interior; nonetheless this is a first attempt to construct a differential model for the rotational dynamics that encodes the non-uniform structure of the star and the entrainment effect in a consistent way, being at the same time computationally easy to implement. In particular, the hypothesis of an array of straight vortex lines and its validity to describe neutron star rotation deserve some further discussion. On the one hand, the tension of a single line is much smaller than the hydrodynamical force acting on a vortex, which suggests bending and twisting of lines so that vortices would actually form a tangle: turbulence is thus expected to develop in neutron stars interiors and indeed, following a criterion for the onset of turbulence originally devised for liquid Helium, Andersson et al. (2007) find that the core is susceptible to become turbulent. On the other hand, as also explicitly remarked by the authors, this study does not account for the stabilizing effect of rotation in suppressing turbulence: the increased effective tension of a bundle of vortices as compared to a single vortex line resists bending at the macroscopic scale, so that the array of vortex lines can remain parallel to the rotation axis and thus implement a superfluid flow which is consistent with a generalisation of the Taylor-Proudman theorem to axisymmetric spinnig-down frictionless fluids (Ruderman & Sutherland 1974). Since superfluid turbulence is a well known phenomenon in laboratory experiments, it is likely to develop in some regime also in neutron stars, but its effect on macroscopic hydrodynamics and its relevance for glitch models are still debated: further effort along these lines is required to take the modelling of glitches to the next level.
As it stands, studies based on the widely used approximation of straight vortices, that is also likely to be valid in some dynamical regimes, are still meaningful and physically relevant.

Some preliminary definitions
The microphysical input that describe the stellar structure are functions of the density: to study the rotational dynamics of the two components we need the superfluid neutron fraction xn(ρ) (that is zero below the drip density), the normal fraction xp(ρ) = 1 − xn(ρ), the entrainment parameter n(ρ), the pinning force per unit length fp(ρ) and the drag parameter η(ρ) arising from the dissipative processes acting on vortex lines where, for example, a dissipative process can be the "electron scattering off vortex cores" (Alpar et al. 1984b;Bildsten & Epstein 1989). Given an equation of state (EOS) for the stellar matter and solving the Tolman-Oppenheimer-Volkov (TOV) equations, all these profiles are easily parametrized as functions of the spherical radius r.
It is useful to introduce the following integrations along a straight line placed at distance x from the rotation axis: where ρn,p(r) = xn,p(r)ρ(r) and z(x) = √ R 2 − x 2 . In the above integrations the spherical radius r is used as a shorthand for √ x 2 + z 2 . We will also use a "cylindrical mean", defined as with normalized weight The normalization factor reduces to the moment of inertia of the n-component if there is no entrainment ( n = 0). Since in all the above definitions the integrand functions are proportional to xn, the stellar radius R can be replaced by R d .

Conservation of vorticity
To derive the dynamical equations we could start from the well known equations for multifluid hydrodynamics, as in Haskell et al. (2012), but here we prefer to present a simple constructive derivation. The first part of the following argument is not totally new (Alpar et al. 1984a), but the effect of entrainment and of the radial stellar structure have not been considered previously. Following the variational formalism originally developed by Carter (1989), the momentum per particle is a linear superposition of the local kinematic velocities of the two components: where we adopted the notation for entrainment discussed by Prix et al. (2002). Now the Feynman relation takes the form of Bohr-Sommerfeld quantization rule: where N is the number of vortex lines enclosed by the contour of integration γ and the factor of two accounts for Cooper pairs. This expression can be made local by introducing an areal density of vortices nv at any point (x, θ) on the equatorial plane. This quantity is well defined because of assumption (ii) and depends only on the cylindrical radius x thanks to assumption (i). Thus, by taking the circulation around a ring of radius x, the Eq. (8) reads Translation of the integration contour along the z-axis always encircle the same number of vortex lines and thus the above relation implies that the azimuthal component of the momentum (denoted by pn) is a function of the cylindrical radius only, namely pn(x) ·ê θ = pn(x). On the one hand, it is convenient to introduce an auxiliary, columnar angular velocity defined by From Eq. (8) this quantity obeys the usual Feynman relation: On the other hand, the angular velocity of the n-component is non-columnar due to the density dependence of entrainment and is given by Ωp .
Stokes' theorem applied to Eq. (9) implies that the density of vortex lines nv is proportional to the component of the vorticity directed along the normal to the equatorial plane Exchange of angular momentum between the superfluid and the p-component is possible only if we allow for radial motion of vortices. Let vL = (v x L , v θ L , 0) the local mean velocity of vortex lines in the "laboratory" rest frame. The conservation of the number of intersections of the vortex lines with the equatorial plane implies the following continuity equation Taking the temporal derivative of Eq. (8) and using the above expressions (in particular the continuity equation), it's easy to obtain This is the first dynamical equation of the model: once the radial velocity v x L (x, t) is given, it expresses the conservation of the vorticity of the superfluid component, in terms of a first order partial differential equation for the auxiliary variable Ωv(x, t).

Total angular momentum balance
The isolated pulsar slows down under the action of the external braking torque due to radiation emission, namely With the substitution vn(x, z) = x Ωv(x) − n(r)Ωp 1 − n(r)ê θ and using the well known property ρn(r) n(r) = ρp(r) p(r), after a simple calculation Eq. (14) takes the form Here we used Eq.
(3) and Ip is defined as and in the absence of entrainment it reduces to the moment of inertia of the normal component; from equation (15) it is apparent that Ip is the moment of inertia associated to the p-component in the presence of entrainment. This quantity is well defined since the integrand is always positive, thanks to the physical constraints on the entrainment parameters valid in any layer of the star (Chamel & Haensel 2006). Obviously entrainment cannot change the total moment of inertia of the star and, since the relation holds, we conclude that Iv defined in Eq. (7) is indeed the correct moment of inertia associated to the auxiliary variable Ωv.

Relation between vortex velocity and lag
Equations (13) and (15) provide a close set of partial differential equations for Ωv and Ωp only after some reasonable functional dependence of the local velocity vL from the dynamical variables Ωv(x) and Ωp is given. We find it useful to introduce a "friction functional" B defined as The explicit form of B depends on how vortex dynamics is modeled; here we refer to a simple but widely used picture, where a viscous force acts to deviate the vortices from corotation with the bulk superfluid and this induces a hydrodynamical lift force (the Magnus force) on the vortices. As suggested in Epstein & Baym (1992) a vortex line experiences a viscous drag force per unit length fD and a Magnus force fM expressed by As discussed in Andersson et al. (2006), the presence of entrainment does not modify the form of the Magnus force, providing that the vector κ is given by In the present case of straight vortices this vector is obviously directed alongêz. Finally, the Kelvin's circulation theorem for ideal and barotropic fluids implies that free vortices in a background flow are transported with the fluid, thus behaving like massless particles. Since there is no inertia, the equation of motion for straight and rigid vortices is given by where the integration is performed over the line, as in Eqs (2)-(4). In components the above equation reads Replacing the variable Ωn with the aid of Eq. (11) and using the definitions given in Eqs (2) where the x-dependence is understood for notational simplicity. Here we have defined the angular velocity lag ω as In the absence of entrainment our lag ω reduces to Ωn − Ωp.
Of course if Ωv(x) → Ωp, vortices corotate locally with the two components and there is no outward motion (18) and (19) suggest to introduce the dimensionless parameters Comparing Eq. (17) and (21) we see that in the present approach the friction functional is only dependent on x, namely B[Ωv, Ωp, x] = B x (x). In the frame of reference of the crust, the velocity of vortex lines at radius x makes an angle θ d with the directionê θ (the dissipation angle) such that locally we have tan(θ d ) = d(x)/a(x). This generalizes the early result of Epstein & Baym (1992) in the presence of entrainment and for a non-uniform star.

Effective pinning
Vortex pinning to the crustal lattice is the fundamental mechanism to store the angular momentum needed to produce a glitch. One possibility to model pinning is to incorporate its effect in the general form of the functional B, in order to reproduce the local mean behavior of vortex lines. In Eq. (18) the mean outward velocity v x L is proportional to B x (x) and pinning can be addressed by saying that only a fraction Y of (unpinned) vortex lines is free to move under the action of the drag force. Formally this means that we are approximating the friction functional as Alternatively, by considering a large population of vortices around the radius x, we can also interpret Y as the "probability for unpinning" (Jahan-Miri 2006). To give an estimate of Y we can use the general scenario of the snowplow mechanism described in Pizzochero (2011): when the local lag ω(x) is greater than a critical lag for unpinning ωcr(x), the vortices break free and are able to move outward. For clarity we stress that an important feature of the original snowplow model is that the particular form of ωcr and the assumption of rapid expulsion of unpinned vortex lines suggest the formation of a thin sheet of accumulated vorticity that slowly creeps outward. However, since our model is dynamical, we don't need to incorporate this "vortex-sheet" assumption: here we just need to generalize the prescription for calculating the snowplow critical profile ωcr when entrainment is present.
The starting point is the mesoscopic vortex-lattice interaction, namely the average pinning force (per unit length) acting on a rigid vortex that passes trough many lattice sites. Donati & Pizzochero (2004, 2006 studied the microscopic vortex-nucleus interaction, whereas Seveso et al. (2016) provide the mesoscopic pinning force per unit length in the crust fP (ρ). The physical meaning of this quantity is that a segment of vortex line can locally unpin if the modulus of the forces (leaving out the pinning force) acting on it equals fP . In the following we present a simple extension of this to the rigid vortex motion described here, based on a "rigid" model of the pinning force: a vortex line unpins when the modulus of the total force acting on it is equal or greater than the integrated value of fP along the line.
Rigid pinning -When a pinning force per unit length fP (x) acts on a segment of vortex line at x , the general equations of motion of rigid vortices become The vectorial quantity fP (x) is not known, since no reliable estimate of the dependence of the microscopic pinning force on the vortex-nucleus separation is available. Note that, despite our notation, in general fP is not the modulus of fP : the mesoscopic quantity fP (ρ) studied by Seveso et al. (2016) represents the threshold for unpinning, namely it is the maximum value that the pinning force can sustain before letting the vortex segment free to move (analogous to a static friction force, which adjust itself to the external force up to a maximum value where motion sets in). Despite this difficulty, we can use the equation of motion to determine the critical lag for unpinning; the viscous drag force is such that |fD|v L =vp = 0, namely it does not act on pinned vortices comoving with the crust. Thus we can write the condition for unpinning in terms of the Magnus force alone, whatever the actual form of fP is: due to the macroscopic rigidity of the vortex line, the Magnus force acting on the line (at a cylindrical radius x) must be greater than in order to unpin the vortex. It follows that is the magnitude of the total Magnus force acting on a vortex when the critical velocity lag for unpinning is reached. Since in this case by using the definition of lag given in Eq. (20) and the definition of b(x) given in Eq. (3) it is straightforward to find that the critical lag profile is given by With no entrainment ( n = 0) this result reduces to the definition of the critical lag profile given in Pizzochero (2011). The simplest way to encode pinning in our equations is thus to choose When the absolute value of the local lag is smaller than ωcr(x) the fraction of unpinned vortices is zero and there is no radial motion of vortex lines; otherwise the vortices contained in a thin cylindrical shell of radius x can move inward or outward depending on the sign of ω(x). The introduction of such a non analytical term into the equations poses a mathematical complication. However we have in mind that our equations can be used to define a numerical simulation of pulsar glitches: discretizing the cylindrical radius in many cells, Eq. (25) just tells us whether vortices in that cell can be moved or not. Vortex-creep -By using Eq. (25) we have the advantage that there is no need to introduce new parameters into the model, even though it is an unfortunate choice from the analytical point of view. We expect that temperature, quantum effects and finite vortex tension all take part in smoothing the step-like form of Y , somewhat similarly to the early vortex-creep model of Alpar et al. (1984a) where the radial velocity is given by .
Here α is a dimensionless temperature that tunes the local rigidity of the pinning, v0 a typical velocity of the microscopic motion of vortex lines andωcr(x) is an opportune critical lag. In the framework of the original work of Alpar et al. (1984a) it is not simple to implement consistently the geometry of the problem (even with our strong assumption of straight and rigid vortices) and it is not clear which prescription we have to use to estimate the cylindrical profiles α(x) andωcr(x). However the expression for v x L of Alpar et al. (1984a) physically motivates the role of temperature and tells us that we can study the effect of a less rigid pinning by using a smooth Y functional; for example one of the simplest possibilities is where in the limit of "zero temperature and rigid vortices" α → 0 we recover Eq. (25). The disadvantage with respect to the original approach of Alpar et al. (1984a) is that now α has no direct physical interpretation. Viscous drag -Another way to encode pinning that is different from Eq. (23) is to consider a large value of the drag parameter η in the crust: suppose that η → ∞ in the interval Rc < r < R d , then d(x) → ∞ for every value of x. By means of Eqs (18) and (19), we have v x L = 0 and v θ L = xΩp, namely in the whole star all the vortices are pinned and co-move with the p-component. The idea is thus to use a viscous drag parameter η[ρ, ω] that has a local dependence on the lag, so that η assumes a very high value when |ω(x)| is well below ωcr. It is not surprising that a strong viscous drag is equivalent to pinning, but the practical problem is that for every value of ω(x) we have to calculate B x (x) using the opportune values η[ω(x), ρ(r)]. As an order of magnitude estimate, η 10 21 s cm g −1 in the crust is sufficient to give a friction parameter ranging from B x (0) ∼ 10 −9 to B x (R d ) ∼ 10 −12 ; this ensures that the vortices are almost blocked on timescales of many years (cf. equation (35)). On the other hand, a rapid expulsion of vortex lines occurs when η is such that d(x) ≈ a(x) and the dissipation angle is nearly π/4.

MACROSCOPIC EQUATIONS FOR THE TWO FLUID MODEL AND RELATION WITH A RIGID MODEL
In the previous section we derived the macroscopic equations of our cylindrical model, Eqs (13) and (15), where the physical (observable) quantity is Ωp. With the aid of the definition given in Eq. (5) we can rewrite this system for Ωv(x, t) and Ωp(t) in a more compact form as where q = Iv/Ip andΩ∞ = Text/I, while the definition of the cylindrical mean, weighted by dIv(x)/Iv = g(x)dx, is given in Eq. (5). The first remark is that the presence of entrainment does not modify the form of the above equations: all the effects of entrainment are contained in the renormalized quantities q and B, provided one uses as effective rotational variables Ωv(x, t) and Ωp(t), which are associated to the effective moments of inertia Iv and Ip respectively; in this sense we can interpret the v-component, directly related to the vortex density, as the actual reservoir of angular momentum when entrainment is present. Also, we can trigger a glitch by perturbing the B functional (or by suddenly increasing the Y fraction in some interval of the variable x): this mimics the sudden unpinning of many vortices. It is worth to note that if in Eq. (27) we constrain ∂tΩp = ∂tΩv for all the x values, it follows that the star responds to the braking torque as a whole, with a spindown rate of −Ω∞. As usual, we refer to this situation as "steady-state"; the value ofΩ∞ can be settled by observing the mean spin-down of the pulsar over a long period of time and of course is independent of entrainment.
As described in the previous section, the non uniform structure of the star and the presence of entrainment influence the friction functional B, the ratio q and in particular the cylindrical mean. On the other hand, the actual value of the stellar radius plays no role since the above equations do not change under the rescaling x → xR d and g → g/R d . In the following we will thus measure the radius in units of R d , namely x lies in the range 0 x 1.
As a direct consequence of the definition of g, we can calculate the angular momentum of the differential vcomponent as

Lv(t) = Iv Ωv(x, t) .
It follows that, given a generic lag profile ω(x), the reservoir of angular momentum associated to the lag is In the last section, we will describe how this fact can be used to set an upper limit to the mass of a large glitcher like the Vela.

Reduction to a rigid model
Most studies of pulsar glitches are based on body-averaged models with two rigid components: an "observable component" strongly coupled to the magnetic field and an "internal component". Some examples of this kind of approach can be found in the early work of Baym et al. (1969) or in the more recent studies of Andersson et al. (2012), Sidery et al. (2010) and Jones (2002). The form of the mutual friction is assumed to be linear in the velocity lag and the dynamical equations are where the label i (o) refer to the internal (observable) component and τ is the phenomenological timescale that describes the coupling between the components. In the original work of Baym et al. (1969), the aim was to study the response induced by a sudden change in the moments of inertia, thus the two components are the crust (plus the charged particles) and the totality of the superfluid neutrons; this gives Io ∼ 10 −2 I. On the other hand, under the assumption that glitches arise in the inner crust superfluid (taken as distinct from the core superfluid), the observable component is comprised by the normal crust plus the totality of the core. In this case the core turns out to be strongly coupled to the normal component thanks to the effect of electron scattering off magnetized vortex cores (Alpar et al. 1984b), another manifestation of entrainment; thus estimates of the moment of inertia of the superfluid in the crust give Ii ∼ 10 −2 I (Link et al. 1999), which is a quite small reservoir of angular momentum (Pizzochero 2011).
As argued by Andersson et al. (2012) and Chamel (2013), the fraction Ii/I is significantly lowered by entrainment, leading to serious difficulties in explaining large glitches. For this reason we decided to work under the assumption of vortex cores that pass trough the core-crust interface (as in the snowplow model), since this gives a larger angular momentum reservoir. It is comforting that this choice is so far favored microscopically.
It's interesting to study under which conditions the rigid model defined by Eqs (29) and (30) can be a reasonable approximation of the Eqs. (26) and (27). For this purpose we define the deviation δf of a generic function f as and an "angular acceleration" quantity A(x, t) as A = ( B + δB)(2Ωp + 2 ω + 2δω + x∂xδω)( ω + δω) .
We can rewrite the Eqs (26) and (27) as a system of three dynamical equations for two rigid variables (the mean lag ω and Ωp, that depend only on time) and a differential one δω(x, t): It's possible to approximate A in such a way that the Eqs (32) and (33) decouple from Eq. (31). A very general B functional makes the problem analytically difficult. On the other hand we can restrict our attention to a finite time interval during which we assume that the friction functional B does not vary rapidly over time and can be approximated with an effective cylindrical profile that is function of x only, namely B ≈ B x eff (x), where B x eff eventually accounts also for the fraction Y . With this assumption, we can decouple equation (31) from equations (32) and (33) by dropping the term x∂xδω and taking in A only the terms that are up to first order in the deviations: In this way the dynamics of Ωp and ω does not depend on δω. Now, if we identify Ωo with Ωp and Ωi with Ωp + ω , the Eqs (32) and (33) are equivalent to Eqs (29) and (30), where q = Ii/Io. We will see later that q ≈ 10 − 20 for realistic star parameters; thence we have a potentially large reservoir of angular momentum. It also follows that the relaxation time can be estimated as The relaxation time given here should be intended as the coupling timescale of a body-averaged model characterized by B x eff . In this way, thanks to our cylindrical mean, we can give a microphysical interpretation of the phenomenological parameter τ . We have also shown that the entrainment effect does not change the form of the rigid equations but only affects the values of Io, Ii (leaving I unchanged) and τ , without any need to introduce new terms.

NUMERICAL INTEGRATION OF THE INPUT PROFILES
To use the dynamical equations (26) and (27) we first need to set the ratio q and the normalized function g(x). This will fix the structure of the star. These parameters depend on the mass of the star and on the EOS of dense matter, as well as on entrainment. On the other hand, all the information about the dynamics of vortices is contained in B; if the prescription of Eq. (23) is used, we have to calculate a static drag-parameter profile B x (x) and the critical profile ωcr(x). These three functions of x and the number q are the macroscopic "input profiles" of the model and can be obtained from micro and mesoscopic physics using the integrations described in the previous sections; to construct them we need the density profile ρ(r) and the density dependent quantities xn(ρ), n(ρ), η(ρ) and fp(ρ).
Since in this preliminary work we are dealing with Newtonian moments of inertia, in order to keep the exposition simple we do not implement the general relativistic conversion between the baryon number density n and the energy density ρ, namely we take n = ρ/mn and nn = xnρ/mn. We will study this aspect, by incorporating the prescriptions of Haensel & Potekhin (2001) for a thermodynamically consistent EOS as well as the corrections for slow rotation in general relativity, in a dedicated future work.

Microphysical input
For a given central density and EOS we solve the TOV equations. This gives the total mass of the neutron star and the  Table 1 and Table  4 of Douchin & Haensel (2001). Here we show also the superfluid fraction consistent with the GM1 EOS in the core (dashed line). For densities below the drip point we have xn = 0 (no superfluid neutrons are present). In the inner-crust xn is the fraction of superfluid nucleons in the neutron gas outside nuclei.
density profile ρ(r) so that we can express all the input profiles as functions of the spherical radius. EOS -We use the SLy EOS proposed in Douchin & Haensel (2001) as a consistent description of the whole star. It provides the relation P (ρ) and the fraction of free superfluid neutrons xn(ρ) (shown in Fig. 2). This equation of state predicts a first order phase transition where nuclear clusters dissolve into homogeneous nuclear matter at baryon density n ∼ 0.76 fm −3 , corresponding to ρc ≈ 1.3 × 10 14 g cm −3 . To find the value of R d we assume that neutron drip appears at density ρ d = 4.3 × 10 11 g cm 3 .
To study the dependence from the EOS we also use the stiffer GM1 equation of state of Glendenning & Moszkowski (1991) in the core, while keeping the SLy to describe the crust. The corresponding xn fraction is shown in Fig 2. In all the figures we measure the density in units of ρ0 = 2.8 × 10 14 g cm −3 , the nuclear saturation density.
Entrainment -Superfluid entrainment in neutron stars was studied by Chamel and collaborators in a series of articles. Here we use the effective neutron mass m * n (ρ) given in Chamel (2012) and in Chamel & Haensel (2006) for the crust and the core respectively, as shown in Fig. 3. The entrainment parameter is n = 1 − m * n /mn where mn is the bare nucleon mass. It is evident from the figure that m * n depends strongly on density. Even though the neutron superfluid can flow without friction in the crust, it is entrained by the nuclei because of Bragg scattering of neutrons by the crustal lattice (Carter et al. 2005). The entrainment parameters that we use in this work are easily obtained as n = 1 − m * n /mn and p = xn n/(1 − xn).
Pinning force -Another important ingredient of the model is fp, the pinning force per unit length resulting from the vortex-lattice interaction in the inner crust. Recently Seveso et al. (2016) have proposed a numerical simulations to evaluate fp(ρ) at different densities in the inner crust, which accounts for finite vortex tension and random orientation with respect to the lattice. In particular we use one of the pinning profile corresponding to in-medium suppressed pairing gap (the case β = 3 and L = 5000 of Seveso et al. (2016), see table 3 therein), as shown in Fig. 4. Since theoretical estimates suggest that the protons in the core could form a type II superconductor, in principle the vortex motion could also be impeded by the interaction with the flux tubes. However, following the phenomenological study of Haskell et al. (2013) this pinning interaction should be small (we are currently studying this issue along the lines of Seveso et al. (2016)). For this reason we choose to constrain fp to be zero outside the crust.
Viscous drag -The drag parameter η is, probably, the most difficult and less studied microphysical input that we need: a realistic treatment of the drag exerted on vortices implies the study of many dissipative phenomena, although we can take into account only the dominant ones in each region of the star. Andersson et al. (2006) give an estimate of the drag on vortex lines due to the electron scattering off magnetized cores as (36) where xp = 1 − xn, p = xn n/xp (Chamel & Haensel 2006) and all the quantities are functions of density. The standard picture is that vortices in the core are magnetized, due to the entrainment effect between the neutron superfluid and the proton superconductor. For this reason the above estimate is valid only for ρ > ρc, thus we take ηes = 0 in the crust.
On the other hand, as pointed out by Jones (1990a) and Epstein & Baym (1992), for ρ < ρc the vortices excite phonons in the crustal lattice. This effect dominates the dissipation when the relative velocity between the vortex line and the lattice is small (|vL − vp| < 10 2 cm s −1 ). However, at present only rough estimates of the drag parameter in the crust can be found in the literature. For this reason, we decide to take as a title of example a fiducial constant value for the "phononic" drag parameter in the inner crust, namely η ph ∼ 10 2 g s −1 cm −1 . This choice is made on the basis of the inferred parameters presented in Haskell et al. (2012), see Table 1 therein.
When the velocity of vortex lines relative to the lattice nuclei is much higher, other dissipative processes come into the play. Following Jones (1992) in this case the dominant dissipative process involves the creation of quasiparticles (Kelvin waves) in the vortex cores. This phenomenon is still little studied and only rough estimates are available at the moment. However the phenomenological approach of Haskell et al. (2012) tells us that the drag associated to excitation of Kelvin waves is roughly seven orders of magnitude greater than the phononic drag, thus we consider η kw ∼ 10 9 g s −1 cm −1 . Figure 5 shows the two full drag profiles as a function of the density given by Eq. (1); in both cases ("phononic" and "kelvonic") the mechanism that dominates dissipation in the core is given by Eq. (36).
The fact that for different vortex velocities we have to interpolate between the phononic drag and a kelvonic drag is a direct manifestation of the approximation given in Eq. (23): the power dissipated per unit length by a vortex is not in general of the form η(vL − vp) 2 and the assumption of a velocity independent (or lag independent) parameter η  . The dimensionless profile B x ph (x) for three neutron stars with mass 1.1, 1.4 and 1.9 M , using η(ρ) = η ph (ρ)+ηes(ρ). Note that we started from a discontinuous η at ρ = ρc but the integration over the vortex length gives a continuous profile. Here we used the SLy equation of state.
forces the dissipated power to be quadratic in the vortex velocity. For this reason the prescription of a lag dependent viscous drag given in Sec. 2.6 seems a viable choice for a dynamical simulation. The difficulty is that to implement this "viscous drag prescription" in a physically reasonable way we need a better understanding of drag and pinning.

Numerical results
We use the physical input described in the previous section to construct the profiles B x (x), wcr(x) and g(x) by performing numerically the integrations over the z-axis for a given cylindrical radius x. To study how the input of our model are sensitive to a change of the EOS in the core, we decided to use SLy and the stiffer GM1. For both cases (SLy for the whole star or GM1 in the core plus SLy in the crust), some structural parameters of the star are summarized in Table   10 Figure 8. The dimensionless profile B x kv (x) for three neutron stars with mass 1.1, 1.4 and 1.9 M , using η(ρ) = η kv (ρ)+ηes(ρ). For vortex lines near the rotational axis the drag is almost entirely given by the electron scattering and there is nearly no difference with respect to the "phonon" case. Here we used the SLy equation of state.    Figure 10. The product g(x)ωcr(x) is plotted for three neutron stars with mass 1.1, 1.4 and 1.9 M using the SLy EOS. The area enclosed by each curves represents the average critical lag, ωcr , corresponding to that mass: they are listed in Table 2. In the insert we zoom the outermost cylindrical shell and the arrows indicate the surfaces x = Rc for the three masses.
1 for three different masses. As expected, the stiff equation of state GM1 gives rise to smaller values of q.
In the following we present the results of the cylindrical profiles for the SLy case only, since the profiles obtained using GM1 in the core are qualitatively similar. Since we want to compare stars with different masses, and we already remarked that the actual value of R d is irrelevant in our model, the profiles are shown rescaling the cylindrical radius to unity. In Fig 6 we plot the weight g(x), that is the fundamental quantity that permits to take into account the spherical geometry and the non-uniform structure of the star by projecting it onto a cylindrical model. Since our cylindrical averages are weighted by the moment of inertia, it is clear that the behavior of g near the origin is ruled by the factor x 3 , while the drop after the maximum depends on the decreasing density and on the spherical shape of the star; then the profile of g is corrected by the details of the EOS used and by the entrainment effect. It is worth to note that the mass of the star does not have a deep impact on the overall shape of the normalized profile g, thanks to the fact that in our model we are free to rescale the cylindrical radius. However by increasing the mass of the star we find that g tends to weight more the outermost cylindrical layers for x 0.9 R d . Moreover the cylindrical shells that are completely immersed in the crust (namely for x > Rc) contribute very little to the cylindrical averages and to the total moment of inertia of the v-component, as can be seen in the insert of Fig 6.  Figures 7 and 8 show the B x profile calculated with electron scattering in the core and either phononic or kelvonic drag in the crust. We label the total dimensionless drag B x in these two cases as B x ph and B x kv . The dominant contribution from electron scattering in the core makes these two profiles quite similar for x 0.6 R d . Notice that switching from phononic drag to kelvonic drag in the crust is relevant also for vortex lines that are not completely immersed in the crust (0.8 R d x < Rc).
The critical lag ωcr given by Eq. (24) is shown in Fig.  9 for the SLy EOS. When GM1 is used in the core the re-sult is qualitatively similar: the critical lag profiles, even in presence of entrainment still have the typical peak of the "snowplow-model". This peak gets closer to the cylindrical shell x = R d as we increase the mass of the star. When GM1 is used in the core, the main difference with respect to the SLy case is the higher value of the central plateau between 0.2 x/R d 0.8.
Finally, in Fig. 10 we show the product g(x)ωcr(x) for the SLy EOS; by definition, the area enclosed by each curve is the mean critical lag ωcr , listed in Table 2. It is worth noting that the main contribution to this comes from regions with x > 0.6 R d , which do not contain the possibly exotic inner core: the cylindrical angular momentum reservoir is spread across the inner crust and the outer core of the neutron star.

AN APPLICATION: ESTIMATES OF THE ANGULAR MOMENTUM RESERVOIR AND OF THE SPIN-UP TIMESCALES
Here we briefly show how some of the results developed in the previous section can be useful to make quantitative estimates even without solving explicitly the equations of motion (26) and (27). Firstly we give a rough estimate of the spin-up timescale: we assume that the outward velocity after triggering a large glitch event is given by Eq. (23) and that all vortex lines unpin so that Y ∼ 1 everywhere. Thus the rigid approximation of Sec. 3.1 implies that the spin-up timescale can be defined as in complete analogy with Eq. (35): the body averaged relaxation time is the weighted harmonic mean of the local relaxation times. This is reliable as long as repinning, or a significant change in the value of the friction functional, occurs for times greater than τ kv since the beginning of the avalanche. Similarly we can calculate τ ph , that is the same quantity but using B x ph . We list these values (in units of the pulsar period P ) in Table 2. These values are only global quantities, intimately related to the rigid approximation described in Sec 3.1. Moreover the listed spin-up timescales are referred to the extreme situation where all the vortex lines are unpinned (Y = 1). For the Vela pulsar, using a period P ∼ 0.09 s, all the listed spin-up timescales are less than a minute and are particularly fast for a soft EOS like SLy. They are definitely smaller than the present observational "blind window" of about 40 s (Dodson et al. 2002), thus pointing to fast dynamical transients that need a much better timing resolution to be revealed: the first few seconds of a glitch event may provide a direct insight into the dynamics of superfluid vortices inside a pulsar.
To give an upper limit to the angular momentum reservoir, we must estimate the theoretical maximal glitch for a given EoS and for different masses. In turn, this can be compared to observations in order to constrain the pulsar mass. Conservation of total angular momentum before and after the glitch reads IΩp + Iv ωpre = I(Ωp + δΩp) + Iv ωpost in terms of the initial and final average lags. We can define a maximal glitch amplitude by taking ωpre = ωcr and ωpost = 0, namely a glitch where initially the whole star is just subcritical (in the snowplow scenario, this maximizes the angular momentum reservoir) while after the glitch the average lag is reduced to zero and the reservoir is "empty". We thus obtain Note that we don't require that after the glitch the lag is everywhere null, but only that the lag is null in the mean at a given instant after the triggering of the avalanche has occurred. Although this maximal glitch does not correspond to a real situation and is probably never realized dynamically, it still sets an upper limit to the observed glitch amplitude. Indeed, a sketch (out of scale) of a typical large glitch event is shown in Fig 11; it is clear that δΩ max p is not really the maximum value of the frequency jump, but rather it is the difference between the value of Ωp when the mean lag is instantaneously zero and the initial Ωp(ttr). An overshoot of the p-component is expected since q 1. Now, the very short spin-up timescales that we find by putting Y = 1 suggest that for a large glitch (where all vortices have been unpinned thus maximizing the avalanche) the overshoot is realized on timescales comparable to τ kv , namely a few seconds for Vela. Thus if we observe a glitch amplitude δΩ obs about a minute after the triggering event (say at time t obs , corresponding to the observational blind window) we can be confident that we are observing the glitch amplitude after the condition of "null lag" has occurred (see Fig 11), so that δΩ max p is actually un upper limit for the observed glitch. Conversely, any dynamical simulation will give a glitch amplitude (measured at times well after the short-time transients) which is smaller than δΩ max p ; therefore only stellar structures that permit to obtain a δΩ max p greater than the largest observed glitch are physically acceptable for a given pulsar. Some values of δΩ max p are listed in Table 2 for different EOS and stellar masses. For example, looking at the listed values of δΩ max p , we can conclude that a mass of 1.4 M is enough to reproduce the maximum observed glitch of the Vela (≈ 2.1 × 10 −4 rad/s) when the SLy EoS is used. With more values than those listed we find the more stringent limit M V ela < 1.45 M . The stiffer GM1 allows for larger masses but still rules out the possibility that the Vela is a very massive neutron star (M V ela < 1.8 M ). These upper bounds to the mass of the Vela can only be lowered by the observation of future, even larger glitch events. Note that the identification of the Vela as an intermediatemass neutron star has also been proposed on the basis of its peculiar surface cooling (Yakovlev & Pethick 2004).
Finally we note that the maximum glitch amplitude δΩ max p does not depend on entrainment, as can be seen directly from its definition in Eq. (38). Our conclusion is thus robust against the uncertainties of the entrainment parameters.
Of course we can repeat this line of reasoning for all pulsars that display large glitches, but a detailed discussion and improvement of the estimate of the angular momentum reservoir will be given in a dedicated work. Here our aim is just to propose that the largest glitch observed in each pulsar can put a constraint on the pulsar mass. This approach to constrain some star properties through glitch observations is very different when compared to the estimate of the mo- Figure 11. Sketch (out of scale) of the expected behavior of a typical glitch, with the unresolved region where an overshoot can occur and the maximal glitch δΩ max p defined in the text. The time at which the glitch is triggered is labeled with ttr, whereas t obs is the earliest time that can be resolved by current observations. Table 2. Timescales and maximal glitches as a function of stellar mass are calculated using the SLy EOS (upper panel) and the GM1 EOS with its consistent xn(ρ) in the core (lower panel). The phonon and kelvon relaxation time are obtained using Eq. (37) and are given in units of the pulsar period P . ment of inertia fraction using the activity parameter (Link et al. 1999;Andersson et al. 2012;Chamel 2013); indeed, the activity parameter is a mean quantity over many decades of pulsar evolution, while here we take the largest observed event, which may be a better indicator of the pulsar's glitching strength.

CONCLUSIONS
In this paper we have developed a dynamical model for superfluid pulsars, that can take into account consistently the layered structure of the star, the differential rotation of the superfluid and the presence of strong density-dependent entrainment. In our treatment we tried to point out that all the complex behavior of the rearrangement of the vortex configuration is hidden into the functional B, introduced in Eq. (17). In our estimate of the Vela maximal mass and of the spin-up timescales we used the simple form of B given by Eq. (23). However it is always possible to refine the model with better proposals for B and to study the consequences on the dynamics of the observable variable Ωp. The main result of the paper is the set of dynamical equations, Eqs. (27) and (26), that constitute the mathematical framework for future simulations and modeling of pulsar glitches. They describe the exchange of angular momentum between two effective components, characterized by moments of inertia Iv and Ip, determined by entrainment, and whose angular velocities are Ωv(x, t) and the observable Ωp(t). Even without solving explicitly the equations, we were able to draw some general quantitative conclusions; in particular, we have sketched a new method to estimate an upper bound to the mass of a large glitcher using observations. Moreover our conclusion on the mass of the Vela pulsar is robust against the uncertainties of the entrainment estimates, as can be seen from Eq. (38).
A central feature of the quantitative example described in the previous section concerns the nature of superfluid vortices extending across the crust-core interface. As explained before, we consider that the change in the symmetry of the order parameter (from the S-superfluid to the P-superfluid) may affect the vortex core structure but doesn't decouple the core superfluid and the crustal superfluid. In most of the existing literature the reservoir of angular momentum is instead restricted to the crustal superfluid. However Chamel (2013) and Andersson et al. (2012) pointed out independently that with the inclusion of entrainment the superfluid into the crust has not enough mobility to fit the activity parameter of many pulsars (in particular the large glitchers, like the Vela). In our model with continuous vortices, also the vortex lines that are not totally immersed in the crust (x < Rc) play an important role and the reservoir of angular momentum is enough to reproduce the largest glitch observed in the Vela, provided its mass is not too large. A more detailed study of the parameter space for the different input quantities as well as simulations of glitches induced by a massive unpinning of vortex lines will be the issue of future work. Of course, this kind of simulation makes more sense if it can eventually be compared to better resolution timing data: the first few seconds of a glitch with their massive exchanges of angular momentum may be far more revealing and constraining for the star structure than the later post-glitch evolution.
Finally we remark that despite the early conclusion of Ruderman & Sutherland (1974), the analogy with laboratory experiments on superfluid Helium tell us that turbulence is likely to develop in some internal regions of neutron stars. However there is far from a consensus as to which regimes will lead to turbulence and how this will develop, and all the work on this topic is at the moment very exploratory. Thus the scope of our paper is to clarify how superfluid entrainment affects the macroscopic dynamics in the simplest situation of straight vortices and to test quantitatively the consequences of such a widespread assumption. In fact, all vortex-creep models (since the seminal work of Alpar et al. (1984a)) are derived under the simplification of straight vortices and do not include the effect of entrainment or stratification. Similarly, models based on the two-fluid formalism with mutual forces have either neglected entrainment  or have been studied under the assumption of uniform entrainment and rigid rotation of the neutron superfluid (Prix et al. 2002;Sidery et al. 2010). At this point it is quite simple to introduce general relativistic corrections for slow rotation, and to test different EoS and alternative input profiles for pinning, entrainment and drag.