Issue 
A&A
Volume 654, October 2021



Article Number  A47  
Number of page(s)  17  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202141499  
Published online  08 October 2021 
Revisiting neutron starquakes caused by spindown
^{1}
Instituto de Física, Pontificia Universidad Católica de Valparaíso, Av. Universidad 330, Curauma, Valparaíso, Chile
email: jsrencoret@gmail.com
^{2}
Núcleo de Astronomía, Universidad Diego Portales, Ejército 441, Santiago, Chile
^{3}
Departamento de Física, Facultad de Ciencias Básicas, Universidad Metropolitana de Ciencias de la Educación, Av. José Pedro Alessandri 774, Ñuñoa, Santiago, Chile
Received:
8
June
2021
Accepted:
18
July
2021
Context. Pulsars show a steady decrease in their rotational frequency, occasionally interrupted by sudden spinups called glitches, whose physical origin is still a mystery. One suggested explanation for at least the small glitches are starquakes, that is, failures of the solid neutron star crust, in which the progressive reduction in the centrifugal force deforms the star, stressing the solid until it breaks. This produces a spinup, dissipating energy inside the star.
Aims. We explore this suggestion by analyzing a mostly analytical model in order to understand the possible consequences of starquakes, particularly whether they can explain at least the small glitches.
Methods. We analyze the deformations and strains produced by the decreasing centrifugal force, modeling the neutron star with a fluid core and a solid crust, each with uniform density and with the core possibly denser than the crust, as a simple approximation to the strong density gradient present in real neutron stars.
Results. The deformation of a star with very different densities in the core and crust is qualitatively different from the previously studied case of equal densities. The former more closely resembles the behavior of a fluid star, in which the corecrust interface is a surface of constant gravitational plus centrifugal potential.
Conclusions. Regardless of the uncertain breaking strain, the glitch activity in this model is several orders of magnitude smaller than observed, even if only small glitches are considered. For a large breaking strain, suggested by simulations, glitches due to starquakes could be roughly of the correct size but much less frequent than observed glitches. The energy released in each such glitch is much larger than in the standard model of angular momentum transfer from a faster rotating superfluid in the inner crust. On the other hand, we cannot rule out that the heating produced by small starquakes could trigger glitches by allowing neutron superfluid vortices to move. We also confirm that stresses in the neutron star crust can in principle support an ellipticity much larger than some observational upper limits from pulsar timing and continuous gravitational wave searches.
Key words: stars: neutron / stars: rotation / dense matter / gravitational waves
© ESO 2021
1. Introduction
Pulsars are among the best clocks known in the universe, with very precisely periodic pulses due to the rotation of a neutron star with a strong, inclined dipolar magnetic field. These pulses slow down progressively as the star loses rotational energy. In addition, they exhibit timing irregularities, the most prominent of which are “glitches”, sudden spinup events, in which about 1% of the secular spindown is recovered (Lyne et al. 2000; Espinoza et al. 2011; Fuentes et al. 2017).
Since the detection of the first glitches (Hewish et al. 1968; Boynton et al. 1969), the scientific community has tried to understand their physical mechanism. Initially, Ruderman (1969) and Baym et al. (1969) explained them as “starquakes”, sudden fractures of the elastic neutron star crust. These authors argued that because of its fast rotation, the pulsar has an ellipticity ϵ_{0}, which, as the star spins down, would change toward a more spherical shape by Δϵ_{0} if the star were entirely fluid. In a neutron star, only the core is fluid, while the crust is an elastic solid; therefore, the crust resists the deformation and accumulates stress, reducing the change in ellipticity by Δϵ = bΔϵ_{0}, where b is the (very small) “rigidity parameter” of the neutron star (e.g., Cutler et al. 2003). When a critical strain is reached, the crust yields and relaxes, so the star becomes more spherical by an amount Δϵ, decreasing its moment of inertia and thus causing a sudden increase in the rotation rate.
In this scenario, the time between starquakes is proportional to the value of at which the crust breaks (set by its critical strain angle), and the size of the associated glitches is proportional to the corresponding value of (Baym & Pines 1971). Thus, the “glitch activity,” , the longterm average spinup rate due to glitches (Lyne et al. 2000), is independent of the rather uncertain critical strain of the neutron star crust but proportional to the rigidity parameter b.
Baym & Pines (1971) formulated the starquake problem by writing the total energy of the rotating, stressed star as
where E_{0} is the energy of the relaxed, nonrotating star, L and I are the angular momentum and moment of inertia, the third term is the correction to the gravitational potential energy due to the rotationinduced ellipticity ϵ, and the last term accounts for the elastic stress in the crust, where ϵ_{0} is an ellipticity at which the crust is relaxed. They evaluated the coefficients A and B as integrals over the displacement field for a completely solid, incompressible star with uniform density, uniform shear modulus, and Newtonian gravity (later reformulated by Carter & Quintana 1975 in general relativity). The presence of the fluid core was taken into account in a nonrigorous, heuristic fashion by restricting the integral for B to the volume of the crust, obtaining the rigidity parameter b = B/(A + B)≈10^{−5} for a star of crustal shear modulus μ = 10^{30} erg cm^{−3} and mass M = 1.36 M_{⊙} (and increasing toward smaller masses, for which the gravity is weaker and the crust is thicker).
At that point, with roughly one year of timing observations, only two glitches had been detected, namely a very small one in the Crab and a substantially larger one in the Vela pulsar, the latter of which spins down more slowly than the former. Assuming that these glitches are “typical” for each of these pulsars and occur roughly once per year, Baym & Pines (1971) concluded that starquakes could at best account for the small Crab glitch but by no means for the much larger glitch of Vela, unless it was an extremely unusual event. As a more viable alternative, Anderson & Itoh (1975) explained glitches as sudden transfers of angular momentum from a weakly coupled neutron superfluid in the inner crust to the rest of the star. Later, the detection of many more glitches in these and many other pulsars allowed us to infer that their size distribution is (at least) bimodal (Lyne et al. 2000; Espinoza et al. 2011; Fuentes et al. 2017), suggesting that different mechanisms might be responsible for small and large glitches. For instance, angular momentum transfer from a superfluid component could be invoked for the large glitches (Pizzochero 2011; Hooker et al. 2015; Howitt et al. 2016), perhaps still leaving room for starquakes to potentially explain the small ones (as mentioned, e.g., by Jones et al. 2017).
On the other hand, Link et al. (1998) and Franco et al. (2000) assumed that glitches are caused by starquakes that fracture the crust asymmetrically (biased by the structure of the stellar magnetic field), causing a damped precessional motion of the star, which would result in an increase in the angle between the rotational and magnetic axes and thus explain the observation that the spindown rate is increased by each glitch. The first of these papers used the model of Baym & Pines (1971), whereas the second improved on it, calculating the displacement field for a star with a fluid core and an elastic crust (still incompressible, with uniform density everywhere in the star, uniform shear modulus in the crust, and Newtonian gravity). From this model, one obtains a five times smaller rigidity parameter b (for the same assumed shear modulus).
The discovery of correlated, quasiperiodic variations in the pulse shape and spindown rate of PSR B1828−11 and its interpretation as free precession by Stairs et al. (2000) motivated Cutler et al. (2003) to consider a more realistic (though still Newtonian) model (following Ushomirsky et al. 2000) with compressible matter whose density and shearmodulus profiles are obtained from equations of state based on nuclear physics calculations, yielding an even smaller value, b ≈ 2 × 10^{−7}. The dependence on the assumed neutron star mass and equation of state was explored by Zdunik et al. (2008) (using an approximate, heuristic formula for B based on the results of Cutler et al. 2003), generally finding similar values (with a spread of a factor of a few).
The possibility of detecting continuous gravitational waves from neutron stars with an elastically deformed crust motivated further studies, such as that of JohnsonMcDaniel & Owen (2013), who calculated the maximum quadrupole moment sustainable by a stressed crust (as well as other, more exotic alternatives) in general relativity, finding that it is ∼6 times smaller than the Newtonian result, a correction that should also apply to the closely related rigidity parameter b.
Recently, Giliberti et al. (2020) numerically calculated the centrifugal deformation of a (Newtonian) polytropic star with different adiabatic indices and confirmed that although the value of the adiabatic index has some effect on the deformations, the strain accumulated between successive glitches of the Vela pulsar is in all cases much too small to break the crust; therefore, if starquakes are invoked to trigger glitches, the crust must always remain close to the breaking strain. The same authors (Giliberti et al. 2019) analytically calculated the centrifugal deformation for a star with uniform but possibly different densities in the fluid core and the elastic crust. Among other things, they analyzed the effect of the Cowling approximation (neglecting perturbations of the gravitational potential; Cowling 1941) for the case of equal densities in the crust and the core, finding that the displacements with this approximation are ≈2/5 of those for the exact solution. They do not evaluate the rigidity parameter b or the angular velocity changes resulting from eventual starquakes, which we do in the present paper.
Despite these theoretical advances and the much increased sample of glitches (now numbering ∼400; e.g., Fuentes et al. 2017), we are not aware of any papers reanalyzing the question first raised (and to some extent answered) by Baym & Pines (1971), namely whether starquakes are large and frequent enough to explain all the observed glitches or even a subset of these, such as the class of small glitches.
In the present paper, we consider the same model as Giliberti et al. (2019), namely a Newtonian star with uniform but possibly different densities in the fluid core and the solid crust^{1}, as a simple approximation to the strong (but continuous) density stratification expected to be present in real neutron stars. In order not to overcomplicate this already quite difficult problem, we apply the Cowling approximation, neglecting the gravitational potential perturbations in our calculations but making an approximate correction for it at the end. This approach allows us to analyze the displacements caused by the centrifugal force in the fluid core and the elastic crust of the star, and to show that there is a competition between the buoyancy force associated with the corecrust density difference (which tends to make the crustcore interface an equipotential surface for the combined gravitational plus centrifugal potential) and the shear stress (which tends to minimize its deformation from the original, unstressed shape).
We go on to show that although starquakes can qualitatively resemble glitches, the glitch activity associated with them is far too small to explain even the subclass of small glitches. The predicted glitch activity, and thus this conclusion, are independent of the assumed critical breaking strain σ_{c}, since the size of the glitches is proportional to σ_{c}, whereas their frequency is ∝1/σ_{c}.
On the other hand, for the large yield strains found in molecular dynamics simulations of neutron star crust matter (Horowitz & Kadau 2009), we find that spindown can produce starquakes at most a few times in the life of a neutron star. We also analyze the temperature increase that could be caused by starquakes, which might allow neutron star vortices to move, thus triggering glitches. Finally, we use the previous results to infer the possible ellipticity that might be sustained by crustal stresses in our model, confirming that it is much larger than observational bounds established for millisecond pulsars.
This article is organized as follows. Section 2 describes the displacement field caused in the star by a changing centrifugal force, we discuss the differences between fluid and elastic matter and evaluate the case of a fluid core and an elastic crust as a function of their density difference and the shear modulus of the crust. Section 3 shows that the predicted glitch activity due to starquakes is much lower than that observed even for small glitches, and we use estimates for the breaking strain to predict the likely size and frequency of glitches caused by starquakes. Other potentially observable effects of starquakes are discussed in Sect. 3.5. Section 4 evaluates the ellipticity that can be sustained by crustal stresses in our model and compare it to observational constraints based on gravitationalwave emission from millisecond pulsars. Section 5 contains the conclusions of this work.
2. Rotational deformations of a neutron star
As explained above, in the starquake scenario the decreasing centrifugal force deforms the neutron star toward a more spherical shape, stressing the elastic crust, which eventually yields, allowing an abrupt, additional deformation that reduces the moment of inertia and thus spins up the star. In an extreme quake (which we consider for simplicity), all elastic stresses are released in the quake. Thus, if, as illustrated in Fig. 1, the star starts from a relaxed state “a” (at birth or after a previous, large quake), we can regard the state “b” just before and the state “c” just after the quake as perturbations of “a” under the same applied force (the change in the centrifugal force), in which the process a → b is subject to the elastic stresses in the crust but a → c is not, that is, in the latter process the star behaves as if completely fluid. Therefore, we need to analyze the displacement fields characterizing the deformation of both an elastic solid and a fluid under a changing force field. For conceptual and notational clarity, we start by introducing the basic concepts and equations of elastostatics, followed by the calculation of the displacement fields for the elastic and the fluid matter.
Fig. 1. Schematic diagram of a starquake driven by spindown, with time t on the horizontal axis and rotation rate Ω on the vertical axis. In state “a”, the star is relaxed. From “a” to “b” it spins down, with its crust deforming and building up stress. At “b”, the crust reaches a critical stress and fractures, the moment of inertia decreases, and the star suddenly spins up to a new relaxed state “c”. 
2.1. Equilibrium elastic deformation under an arbitrary force
The elastostatic equilibrium condition in Newtonian gravity is
(e.g., Link et al. 1999; Cutler et al. 2003; Thorne & Blandford 2017), where T is the stress tensor^{2}, ρ is the mass density, g = −∇ϕ_{g} is the acceleration of gravity (ϕ_{g} is the Newtonian gravitational potential), and f_{c} is some other force density field acting on the stellar matter (which we later take to be the centrifugal force due to the stellar rotation but which could, in principle, represent other forces, such as a Lorentz force of the stellar magnetic field). The stress tensor can be decomposed as
where p is the pressure, g is the metric tensor of Euclidean, threedimensional space (not to be confused with the gravitational acceleration g defined above), and the symmetric, tracefree tensor S represents the shear stress, which vanishes in a fluid and in a relaxed solid. With this decomposition, Eq. (2) can be rewritten as
where
is the elastic force density produced by the shear stress.
If we consider a stressed state “b” (with f_{sh} ≠ 0) close to an unstressed state “a” (with f_{sh} = 0), the small Eulerian (local) changes of the variables (denoted by δ[.] ≡ [.]_{b} − [.]_{a}) are related by
We also need to consider Lagrangian (material) perturbations, which give the change in a physical quantity in a given matter element as it is displaced, e.g., , where ξ represents the displacement vector field from the position x of a matter element in state “a” to its position x + ξ in state “b”. In the linear limit, Eulerian and Lagrangian perturbations of any fluid variable are related by
In the elastic limit, the Lagrangian stress perturbation is proportional to the strain tensor (generalized Hooke’s law)^{3},
where the rankfour tensor Y^{ijkl} is known as the Young tensor (a property of the material) and
is the strain tensor, which can be decomposed into its irreducible tensorial parts as
where the trace
represents the expansion (or compression) of the matter element, the symmetric, traceless part
represents the shear deformations, and the antisymmetric part
represents rotations (see Thorne & Blandford 2017 for a pedagogical discussion).
In an isotropic solid, Hooke’s law (Eq. (8)) can be written separately for the two pieces of the stress tensor (Eq. (3)) as
and
where K is the bulk modulus, which quantifies the resistance of the matter to expansion or compression, and μ is the shear modulus, which accounts for the resistance of the solid to shear deformations (rotations do not produce stress; Thorne & Blandford 2017). Thus, the shear force in Eqs. (4) and (6) can be written as
with Σ given by Eq. (12).
In order to obtain the displacement field ξ(x) caused by the change δf_{c}(x) in the applied force field, we must solve Eq. (6) with the appropriate boundary conditions, to be discussed in Sect. 2.4.
2.2. Elastic deformations in a rotating star
For mathematical simplicity and to allow physical insight, we consider a simple model in which the neutron star matter is incompressible (K → ∞), with Θ = 0 = Δρ (but Δp = −KΘ ≠ 0). We also assume uniform density in each domain (core and crust), so the Eulerian density perturbation δρ vanishes as well (except at the boundaries). We also consider a uniform shear modulus μ throughout the crust (and of course μ = 0 in the fluid core). Up to this point, the model is identical to that of Franco et al. (2000). However, as recently done by Giliberti et al. (2019), we allow for different densities in the two domains, ρ_{s} in the elastic, solid crust, and in the fluid core, which it makes an important difference.
When the star spins down, it is deformed (toward a progressively more spherical shape) by the change in the centrifugal force density, δf_{c} = −ρ∇δϕ_{c}. Since we take ρ to be uniform and δρ = 0, Eq. (6) implies that the shear force can be derived from a potential,
where we introduced a total (effective) potential, ϕ_{T} ≡ ϕ_{g} + ϕ_{c}. Thus,
where the constant h needs to be obtained from the boundary conditions (which otherwise cannot be satisfied)^{4}.
Choosing a spherical coordinate system (r, θ, φ) aligned with the axis of rotation, the centrifugal potential perturbation can be written as
where Δ(Ω^{2}) is the change in the squared angular frequency, and P_{0}(cos θ) = 1 and are the Legendre polynomials of order ℓ = 0 (monopole) and ℓ = 2 (quadrupole), respectively.
The incompressibility condition, Θ ≡ ∇ ⋅ ξ = 0, forces the monopole component of the radial displacement field to vanish^{5}, . Thus, motivated by the form of the centrifugal force, we assume a purely quadrupolar displacement, with
Again applying the incompressibility condition, and using the Legendre equation, the polar component takes the form
Equation (17) implies that
where the nonzero components of the strain tensor are
For the displacement field given by Eqs. (21) and (22), Eq. (23) becomes an ordinary differential equation of one independent variable,
whose solution we conveniently write as
with the dimensionless radial displacement function
(see also Franco et al. 2000; Giliberti et al. 2019). Here, we have defined the dimensionless radial coordinate, , as well as the Keplerian (or “breakup”) frequency,
for a star of radius R_{⋆} and mass M_{⋆}. The coefficients A_{1}, A_{3}, A_{−2} and A_{−4} are dimensionless constants to be obtained from the boundary conditions. For later reference, we also write down the explicit expression for the shear potential,
2.3. Fluid displacements in a rotating star
We show below that, although the force balance equation for the fluid is the same as that for the elastic solid but with μ = 0 (and thus f_{sh} = 0), in some cases the displacement field for the fluid cannot be obtained by taking the limit μ → 0 in the results for the solid. One formal indication of this is that, for μ = 0, Eq. (23) is trivially satisfied and therefore does not imply Eq. (28) and its solution, Eq. (30).
On a conceptual level, we note that, in the elastic solid, neighboring matter elements cannot be displaced independently without causing large stresses, and therefore the vector field ξ(x) is well defined, with no ambiguity (except for a solidbody rotation around the symmetry axis). However, in a uniform, incompressible fluid, all matter elements are equivalent and interchangeable, so there is no welldefined correspondence of fluid elements in different equilibrium states, even if these states are very similar. Thus, any “small” displacement field satisfying the incompressibility condition ∇ ⋅ ξ = 0 and the appropriate boundary conditions is equally valid. This also implies that, although Eq. (18) with Ψ_{sh} = 0 can be used to determine the Eulerian pressure perturbation^{6}δp(r, θ), the ambiguity in the displacement field does not allow us to uniquely define a Lagrangian pressure perturbation, Δp = δp + ξ ⋅ ∇p.
Of course, the assumption of uniform, incompressible matter is an extreme idealization that does not apply to a real neutron star, in which the strong gravity causes a density gradient, and variations in the centrifugal force (or other forms of forcing) change the density. Furthermore, the composition of neutron star matter is not uniform. In equilibrium, the relative abundances of particle species (neutrons, protons, electrons, and possibly other particles in the core; different nuclear species, electrons, and neutrons in the crust) are functions of the local pressure. However, if the pressure in a certain matter element is changed, its composition (here denoted by a generic label Y that does not require a precise definition) does not respond instantaneously but on a finite timescale set by the required strong or weak interaction processes. Over short times (such as typical neutron star oscillation periods), the composition can be taken to be (approximately) conserved in each matter element, and the density must be treated as a function of (at least) two variables, namely pressure and composition, ρ(p, Y) (Finn 1987; Reisenegger & Goldreich 1992), corresponding to a nonbarotropic (or “baroclinic”) fluid. Over long timescales, such as the spindown time of a pulsar, the composition in each fluid element of the neutron star core adjusts, staying close to its equilibrium value at the local pressure (Reisenegger 1995), so the matter can be regarded as barotropic, with the density as a unique function of pressure, ρ(p). This may not be a good approximation for the neutron star crust, where successive equilibria are not a continuum, but they are separated by potential barriers (Gusakov et al. 2015).
In a nonbarotropic star with a changing centrifugal force, the displacement field ξ(x) is much better defined than in the uniform, incompressible idealization. The force balance equation,
implies that ∇p and ∇ϕ_{T} are everywhere parallel to each other. Therefore, in a given equilibrium state, we can regard the pressure as a function of the total potential, p = p(ϕ_{T}). Similarly, the curl of Eq. (33) implies that ρ = ρ(ϕ_{T}), and, from the equation of state, ρ = ρ(p, Y), we conclude Y = Y(ϕ_{T}). Thus, in any given equilibrium state, there is a family of nested equipotential surfaces (spheroids) on each of which all four quantities (ϕ_{T}, p, ρ, and Y) are constant, and, if Y is conserved in each fluid element, it allows to identify the equipotentials containing the same fluid elements to be identified in different equilibrium states. Special cases of these equipotentials are the stellar surface (as long as the outermost layer is fluid) or any phase transition or density discontinuity between two fluid regions within the star. This also means that, for a small perturbation connecting two equilibrium states, the Lagrangian changes Δp and Δρ must take the same value everywhere on a given equipotential (and of course ΔY = 0). Since the Eulerian changes are well defined, we can use any of the fluid variables to constrain the displacement field (e.g., ξ ⋅ ∇p = Δp − δp); therefore, the component of the displacement field perpendicular to the equipotentials is uniquely fixed, and the parallel component can be obtained from solving ∇ ⋅ ξ = −Δρ/ρ with the appropriate symmetry and boundary conditions (as it was done to obtain Eq. (22) from Eq. (21)).
The equipotentials still exist (with p, ρ, and ϕ_{T} constant on each of them) in a barotropic fluid (likely a good approximation for the core of a neutron star that spins down slowly) and in the idealized uniform and incompressible case. However, in these cases, there is no guarantee that the set of fluid elements lying on one equipotential in some equilibrium state still lie on one equipotential in another equilibrium, as they can be interchanged (allowing them to adjust to the local pressure) without changing the structure of the equilibrium state. Still, we can choose to identify fluid elements on equipotentials in two different equilibria, using a conserved integral quantity, such as the total baryon number enclosed by each equipotential, to match equipotentials in different equilibrium states. Although this choice is arbitrary (analogous to fixing the gauge for the electromagnetic potentials), it has the advantage of including the stellar surface and internal discontinuities as special cases of such equipotentials (if the adjacent matter is fluid).
Thus, in our incompressible, uniform density model, we proceed by writing Eq. (18) with Ψ_{sh} = 0 as
In what follows, we assume a slow rotation, Ω ≪ Ω_{K}, so the background equilibrium model can be taken as spherically symmetric, with p = p(r) satisfying the hydrostatic equilibrium condition dp/dr = −ρg(r), where g(r) = ∇ϕ_{g} is the local gravity, and with Δp = Δp(r). Since ξ^{r}(r, θ) is purely quadrupolar, the monopole component of Eq. (34) gives
whereas the quadrupole part yields
In the last equality, we also used the Cowling approximation, namely neglecting the gravitational potential perturbation, δϕ_{g}, so δϕ_{T} ≈ δϕ_{c}. Equation (36) determines the shape of all equipotentials in the fluid regions of the star, particularly the surface and the crustcore interface, as long as both the crust and the core can be modeled as fluids. For our toy model, with a uniform crust and a uniform core, Fig. 2 shows the radial dependence of this displacement field. The change in the slope at the interface is caused by the assumed density discontinuity, quantified by , through its effect on g(r). We note that, in the crustal region, , Eq. (36) implies
Fig. 2. Radial dependence of the normalized quadrupolar radial displacement field, , of the equipotential surfaces in a completely fluid star with uniform densities, ρ_{s} in the crust () and in the core (), with a dimensionless interface radius and different values of . 
which does not have the polynomial form of the elastic displacement (Eq. (30)), so these two functions can at best agree at a few discrete values of .
2.4. Boundary conditions
At the interfaces between domains with different properties (in the present case, the core and crust of the neutron star), we must impose continuity of the perpendicular displacement field, n ⋅ ξ, and of the stress across the boundary, n · T, where n is the normal vector of the interface. At external boundaries (the stellar surface, separating the crust from the exterior, approximated as a vacuum), the first of these conditions becomes meaningless, and the second implies n · T = 0 at the boundary.
If both the core and the crust of the star are uniform, incompressible fluids (of possibly different densities), the displacement of all equipotential surfaces (including the interface and the surface) is given by Eq. (36). Therefore the continuity condition on n ⋅ ξ is trivially satisfied, and the boundary conditions on n · T are only required for the calculation of the pressure perturbations. (In Eq. (35), the constant h could be different in the two domains.)
If, on the other hand, the crust supports elastic stresses, the displacement field does not follow equipotential surfaces and thus does not satisfy Eq. (36). This applies, in particular, to the displacement of the crustcore interface, where the fluid displacement field cannot be given by this equation and thus the continuity condition on n ⋅ ξ becomes inapplicable.
On the other hand, we do need to apply the boundary conditions on n · T. The displacement of the stellar surface clearly follows the matter displacement field evaluated at the stellar radius, so the condition of zero external stress implies a condition on the Lagrangian perturbation of the stress tensor
At the crustcore interface, the condition is less obvious. As the pressure at the interface increases, strong interaction processes quickly dissolve nuclei or “pasta” structures (Ravenhall et al. 1983) in the adjacent crust into free nucleons and electrons (and the opposite process if the pressure decreases), converting some solid to fluid (or vice versa). We follow the previous literature (e.g., Baym & Pines 1971; Franco et al. 2000; Giliberti et al. 2019) in ignoring this effect, which affects only a small fraction of the matter, and therefore we also assume that the interface moves with the local displacement field of the solid, yielding
at the interface radius R_{i} between the fluid core (labeled by f) and the solid crust (labeled by s)^{7}
In our model, these boundary conditions become
Equations (40) can be rewritten as
providing two independent linear homogeneous equations relating the four constants A_{k} in Eq. (30),
Equation (41) is most conveniently written as
This equation has a monopolar component, yielding , and a quadrupole component,
or, equivalently,
We note that in Eqs. (46) and (47) the lefthand side is proportional to ρ_{f} − ρ_{s} and the righthand side is proportional to μ. Thus, in the limit μ → 0 with ρ_{f} ≠ ρ_{s}, Eq. (46) becomes , that is, the interface becomes an equipotential surface (see Eq. (36)), because the solid behaves essentially like a fluid. On the other hand, for ρ_{f} − ρ_{s} → 0 with μ ≠ 0, the boundary condition is very different, namely , a homogeneous linear condition involving up to third derivatives of . Thus, as further discussed in Sect. 2.6, the two limits and μ → 0 do not commute, as they yield qualitatively different interface conditions.
Similarly, at the surface of the star, Eq. (42) yields
implying for the monopolar component and
or
for the quadrupolar component. In the limit μ → 0, the righthand side of Eqs. (49) and (50) vanishes, so the stellar surface also becomes an equipotential.
In all cases, the four linear Eqs. (44), (47), and (50) yield a unique solution for the four constants A_{k}, thus fully determining the displacement field in the neutron star crust.
2.5. Displacement field for a uniform star
Before presenting the full analytic solution for the displacement field in a twodomain star (crust and core with different properties), it is useful to consider the simple cases of a completely uniform (singledomain) star, both as a fluid and as a uniform elastic solid. For the twodomain model, with our convention of assuming that fluid elements move with the equipotentials (which they would, in the case of a nonbarotropic star), the dimensionless displacement field can be obtained from Eq. (36) as
where is the normalized radius. For the onedomain model (the case considered by Baym & Pines 1971), regularity of Eq. (30) at r = 0 implies A_{−2} = A_{−4} = 0, which, together with Eq. (43) (imposed only at the surface) and Eq. (50), gives
with the dimensionless parameter^{8}
where the numerical subscripts stand for the exponents in the reference magnitudes for neutron stars, expressed in cgs units, e.g., μ_{30} ≡ μ/(10^{30} erg cm^{−3}). Figure 3 compares these two cases, considering several values of β, which are chosen unrealistically large for the purpose of illustration.
Fig. 3. Dimensionless radial displacement in the uniform solid star for different values of the parameter β defined in Eq. (53), compared to the case of an uniform fluid, with the displacement field in the latter case defined in terms of equipotentials, as explained in Sect. 2.3. 
The radial displacement at the surface () of the elastic star is smaller by only a fractional amount 𝒪(β) than that of its fluid counterpart. On the other hand, the internal displacement is very different, as it increases linearly with r in the fluid case, whereas in the solid case it has a convex shape with a maximum at r/R_{⋆} = (8/9)^{1/2} ≈ 0.94, where the displacement is slightly larger (by a factor of ) than on the surface. As discussed in Sect. 2.3, the equilibrium displacement field within a uniform or barotropic fluid is a matter of convention. Our choice is based on the case of a nonbarotropic fluid, in which equipotential surfaces of one equilibrium must deform into equipotential surfaces of another.
If instead the matter is a nonbarotropic, elastic solid, there is a conflict between the tendency to follow equipotentials (as in Eq. (51)) and the tendency to minimize the shear stresses in the crust (as in Eq. (52)).
2.6. Displacement field for two uniform domains
The uniformdensity model of Franco et al. (2000) accounts for the shear stresses, but of course it does not include the buoyancy force associated with the density gradient present in a realistic neutron star. This issue is addressed in the present paper (as in Giliberti et al. 2019) by allowing for a density difference between the core and the crust, which makes the crust float on top of the core, so their interface would be an equipotential if shear forces could be ignored. This is a strong simplification, which reduces a gradient in density and composition (from the center to the surface of the star) to a single step. However, it has the advantage of being analytically calculable and introducing a single additional parameter, , thus making it easier to analyze than the more realistic numerical calculations of Cutler et al. (2003) or Giliberti et al. (2020). Using the average densities in the core and crust for the range of neutron star models shown in Fig. 1 of Fattoyev et al. (2018), one obtains a density contrast in the range .
The boundary condition at the corecrust interface (Eqs. (46) and (47)) implies that our model reduces to the uniformdensity model of Franco et al. (2000) when (except that we use the Cowling approximation, δϕ_{g} = 0, whereas Franco et al. 2000 calculate δϕ_{g}[r, θ] selfconsistently). On the other hand, considering also the surface boundary condition (Eqs. (49) and (50)), one can see that our model must approach the pure fluid case (with the displacement field of the crustcore interface and stellar surface following equipotentials) when .
The full solutions for the coefficients of the displacement field in the elastic crust with the boundary conditions presented in Sect. 2.4 are given in Appendix A, and specific cases are plotted in Figs. 4 and 5.
Fig. 4. Normalized radial displacements in the neutron star crust with interface radius , both for purely fluid matter (dashdotted line) and for an elastic crust with β = 10^{−4} (red dots). The upper panel corresponds to a very small density difference, , whereas the lower panel shows a star with a much denser core, . 
Fig. 5. Normalized radial displacement at the stellar surface (r = R_{⋆}; upper panel) and at the crustcore interface (r = R_{i} = 0.9R_{⋆}; lower panel) as a function of the normalized density difference between the core and the crust. The blue line corresponds to a completely fluid star, whereas the other lines correspond to stars with an elastic crust, for different values of β. 
Evaluating the full solution at the corecrust interface yields
and at the stellar surface,
where the coefficients n_{k}, d_{k}, d_{Δk} (for various integers k), defined in Appendix A, depend only on the dimensionless crust thickness , whereas . All these coefficients have been defined in such a way that they approach 1 as .
In these expressions, it is clear that the limits , , and β → 0 do not commute. Here we discuss them one by one.
In the simplest case , in which the crust disappears, so β and become irrelevant and , since both quantities correspond to the surface of a purely fluid star of uniform density.
In the uniformdensity case (), the exact solutions are
and
In this case, it can be seen that the displacement field depends on the shear modulus only through the constant factor of in the denominator, so its value affects only the amplitude but not the shape of the displacement field. Taking the limit β → 0, this solution, like the completely solid star of Baym & Pines (1971), approaches the fluid limit at the surface but not in the interior (particularly at the interface), where the fluid displacement field was chosen to follow equipotentials. This limit is illustrated by the upper panel of Fig. 4, where .
In the arguably more realistic case of , we can expand the normalized radial displacement at the interface as
and at the surface as
where we see that the term independent of β corresponds to the case of two fluids with different densities (Eq. (36)), and the crust elasticity comes in as a small correction ∝β. However, we can see in the lower panel of Fig. 4 that within the crust (), the formal solution in the limit β → 0 does not correspond exactly to our fluid displacement field (chosen to follow equipotentials).
3. Modeling a glitch as a starquake
3.1. Rigidity parameter
In the starquake model, illustrated in Fig. 1, the change in the rotation rate in a glitch, ΔΩ_{bc}, is due to conservation of angular momentum as the star’s moment of inertia I changes by ΔI_{bc} when its crust breaks and relaxes,
The change in the moment of inertia due to the starquake can be calculated as
where the change from the initial, relaxed state “a” to the stressed state “b” is resisted by the crustal shear stresses, whereas the change from “a” to the relaxed state “c” can be modeled in terms of a purely fluid star (with μ = 0). Both changes ΔI_{ab} and ΔI_{ac} can be written in terms of the change in the squared rotation rate between starquakes, Δ(Ω^{2}), and the induced normalized displacement fields as
where, for the stressfree (“fluid”) displacement (ΔI_{ac}), the dimensionless function must be evaluated at β = 0. Using Eqs. (61) and (62) and assuming slow rotation and small β, the “rigidity parameter”, defined by Cutler et al. (2003) as the ratio between the change in moment of inertia at a starquake and between starquakes, can be approximated as
where we have neglected a term ΔΩ_{bc}/ΔΩ_{ac} ∼ b(Ω/Ω_{K})^{2}, which is small in the limit of slow rotation, Ω ≪ Ω_{K}. We can evaluate
which, in the case without shear stresses in the crust (the purely fluid case; Eq. (37)) becomes
Replacing in Eq. (63), one obtains a general expression for the rigidity parameter in a star with a uniform crust and a uniform core with possibly different densities,
(A more explicit expression for b is given in Eq. (B.1).)
It can be useful to rewrite the dimensionless elasticity parameter β in two different ways:
where
is the volume of the crust,
is the gravitational binding energy of our model star, and is the gravitational binding energy of a uniform sphere of the same mass and radius.
Figure 6 shows that for all values of and for a density contrast in the range , the rigidity parameter b does not differ from the dimensionless combination μV_{c}/E_{g} by more than an order of magnitude.
Fig. 6. Rigidity parameter b/(μV_{c}/E_{g}) as a function of the difference in density between fluid core and solid crust in the horizontal axis and as a function of the core radius in the vertical axis. The fiducial parameters are fixed so that . 
The horizontal axis of Fig. 6, , corresponds to the case of a uniform elastic sphere studied by Baym & Pines (1971). In this case, we obtain and
where V_{T} is the total volume of the star, and the approximation corresponds to the realistic limit . The vertical axis, corresponding to , approaches the case studied by Franco et al. (2000), namely a star with a uniform fluid core and a uniform elastic crust, both with the same density, for which Eq. (66) with yields
where the latter approximation also assumes , making d_{1} = d_{2} = γ_{2} = 1.
A somewhat more general expression can be obtained by taking simultaneously , , and in Eq. (B.1). This allows us to set w_{0} = 1 and , keep only the lowestorder terms in and in the numerator and the denominator, and keep the lowestorder term in in each of their coefficients, which yields
Thus, as long as , we recover the result of Eq. (72), whereas for we have
The change of behavior near can be clearly seen in Fig. 6. To the right of this point, the density contrast is still small, but the elasticity parameter is even smaller, corresponding to taking first the limit and then , that is, in the opposite order than in Eq. (72). If the limit is taken first, the crustcore interface can be viewed as separating two fluids of different densities, so the interface is roughly an equipotential. If, instead, the limit is taken first, the crustcore interface separates a fluid and an elastic solid of the same density, so the magnitude of its deformation depends crucially on the shear modulus of the solid.
Since, in a real neutron star, is very small (∼10^{−5}, see Eq. (69)), the crust is relatively thin (), but the contrast between the average crust and core densities is in the range , we consider the limit of and but finite (roughly in the range ), for which we obtain
As the density contrast increases, this function increases from 9/5 = 1.8 at toward an asymptotic limit of 72/5 = 14.4 for . We note that these values exceed the one obtained for and (Eq. (72)) by a factor between 33/8 and 33. For , we obtain 12.7 < b/(μV_{c}/E_{g}) < 13.8, not far from the asymptotic limit. Taking b/(μV_{c}/E_{g}) ≈ 13, noting that in this regime , with given by Eq. (69), and applying a correction factor of 5/2 to account for the effect of gravitational potential perturbations (see Sect. 3.2), we obtain our best estimate for the rigidity parameter of a neutron star, based on our model with a uniform crust and uniform core of different densities,
3.2. Effect of the Cowling approximation
Both our results for an uniform, elastic star (Eq. (71)) and for a star with an elastic crust and a fluid core of equal, uniform densities (Eq. (72)) are lower by a factor of 2/5 than their counterparts found in the literature, namely the classical result of Lord Kelvin (reproduced by Baym & Pines 1971 and Cutler et al. 2003) for the first case, and the result of Cutler et al. (2003) for the model of Franco et al. (2000) for the second.
This factor can be understood from the calculations by Giliberti et al. (2019) of the displacement fields for a star with an elastic crust and a solid core. From their Eqs. (10) and (57), the quadrupolar radial displacement field can be written as
when gravitational potential perturbations are taken into account, and
when they are not, that is, in the Cowling approximation. In both equations, P and Q are functions whose precise form is unimportant for the present argument, and we have neglected higherorder terms in and (q and χ^{2}, in the notation of Giliberti et al. 2019), both of which are assumed to be small. From our Eqs. (62) and (63), we see that the rigidity parameter depends on the fractional difference between the displacements for a purely fluid star () and a star with an elastic crust (),
For the case with gravitational potential perturbations, the righthand side of this equation is , whereas in the Cowling approximation it is ; therefore, the ratio of the rigidity parameters is, to lowest order in and ,
confirming the factor found in our comparison.
We note that this ratio 2/5 does not come from the overall prefactor in Eq. (78) but from the ratio of the coefficients (24 and 60) of the small correction terms of order . It is by no means clear whether in the more general case with different densities this ratio is still 2/5 or whether it now depends on the density contrast . We can only guess that the ratio might not be very different, so our results for the rigidity parameter b in Sect. 3.1 should all be multiplied by a factor of ≈5/2 in order to account for gravitational potential perturbations (except the final Eq. (76), where it is already taken into account).
3.3. Recovery fraction and glitch activity
Combining Eqs. (60), (62), and (63), we obtain the ratio of the spinup in a glitch to the spindown between glitches,
valid for all the cases, where , given by Eq. (65), decreases from 1 for to for , with in the range , which we adopt as our standard value in what follows. The ratio given in Eq. (81) is particularly interesting, because both the numerator and the denominator on the lefthand side are measurable in observed glitches. If all glitches were equal, we could identify
where ν = Ω/(2π) is the rotation frequency, is its timederivative (the spindown rate of the pulsar),
is the glitch spinup rate, customarily called “glitch activity” (Lyne et al. 2000), which accounts for the cumulative spinup effect of glitches during the time of observation T of a given pulsar, and Δν_{i} is the size of its ith glitch observed during this time. Therefore, from Eq. (81), the glitch activity in the starquake model can be written as
where is the Keplerian frequency. We note that all glitching pulsars have spin frequencies substantially lower than 1 kHz (e.g., ν ≈ 11 Hz for Vela), whereas many (perhaps most) have measured glitch activities (Lyne et al. 2000; Espinoza et al. 2011; Fuentes et al. 2017). Thus, starquakes (crust breaking events) are clearly far from being able to account for all glitches. Figure 7 shows that the observed glitch spinup rates in most bins are 5–6 orders of magnitude higher than the prediction based on starquakes, clearly indicating that this model cannot account even for the subsample of small glitches (Δν < 10 μHz) identified from the observed bimodality in glitch sizes (Espinoza et al. 2011; Fuentes et al. 2017). We emphasize that this result is independent of the rather uncertain breaking strain of the neutron star crust material, which determines the size and frequency of starquakes but without affecting .
Fig. 7. Glitch activity as a function of . Black dots represent the glitch activity (sum of all glitch sizes divided by the sum of the observation times) for all pulsars within bins of , and blue dots represent the same quantity but restricted to small glitches (Δν < 10 μHz). (Data kindly provided by J. R. Fuentes, from the glitch sample discussed in Fuentes et al. 2017.) The green and red lines represent the glitch activity in the starquake model, as given by Eq. (84), for the values b = 10^{−5} and b = 10^{−6} of the rigidity parameter, corresponding to the cases and , respectively. 
3.4. Crust breaking, glitch size, and glitch frequency
The two most widely used criteria for the failure of solid materials are the von Mises criterion,
and the Tresca criterion,
where σ_{1} ≥ σ_{2} ≥ σ_{3} are the eigenvalues of the strain tensor σ, and σ_{M, c}, σ_{T, c} are dimensionless “critical strains”, a property of the material that is difficult to predict from theory (Christensen 2013). These two criteria can be compared by considering as a function of the intermediate eigenvalue σ_{2}. This function is represented by a parabola with minimum at σ_{2} = (σ_{1} + σ_{3})/2, at which σ_{M} = σ_{T}. On the other hand, at the extremes of the interval σ_{3} ≤ σ_{2} ≤ σ_{1}, one obtains , and therefore
making the two criteria essentially equivalent if σ_{M, c} ≈ σ_{T, c}, so we chose to work with the mathematically more tractable von Mises criterion.
In incompressible matter, the strain tensor is traceless; therefore σ^{ij} = Σ^{ij}, and we can also write
In our model, this expression can be evaluated using the equations in Sects. 2.2 and 2.6, and rewritten as
where is a dimensionless function of position in the crust, shown in Fig. 8.
By analogy with solids on Earth (at essentially zero pressure) and attempting to match the observed glitch phenomenology, it was originally assumed that the critical strains of neutron star crusts were very small, ≲10^{−4} (Smoluchowski 1970). However, molecular dynamics simulations (Horowitz & Kadau 2009; Hoffman & Heyl 2012) of solids subject to very high pressures, as in neutron star crusts, support much larger values. Horowitz & Kadau (2009) use a Cartesian box, applying an incompressible deformation^{9}ξ^{x} = αy/2, ξ^{y} = αx/2, ξ^{z} = 𝒪(α^{2}), and finding that the material yields when α ≈ 0.1, corresponding to σ_{M, c} = σ_{T, c} = α/2 ≈ 0.05 according to our definitions above^{10}. Baiko & Chugunov (2018) find that the breaking strain can vary by a factor of ∼3 with respect to the result of Horowitz & Kadau (2009), depending on the direction of the deformation. However, given the many other uncertainties involved in this problem, we use the constant value of Horowitz & Kadau (2009) as our reference.
Since the offdiagonal strain components Σ^{rθ} are forced to vanish at both the stellar surface () and the crustcore interface (), they are generally much smaller (by a factor of ) than the diagonal components, which can be written as
Thus, ignoring the offdiagonal terms, the squared von Mises strain becomes
which is a decreasing function of sinθ for small values of sinθ, and an increasing function for sinθ ∼ 1. Therefore, its maximum value, where the crust first breaks, must lie either at the equator (θ = π/2), which happens for , or at the pole (θ = 0), which happens outside this range.
For very small (the condition, for , is ), the radial displacement function has a maximum inside the star and decreases slightly from the interface to the surface (upper panel of Fig. 4), so , implying . Thus, as the star reduces its equatorial bulge under spindown, the crust becomes thicker at the equator and thinner at the poles. As shown in the first panel of Fig. 9, the strain is mostly due to the nonradial components and is maximum at the equator, where the crust is expected to break first.
Fig. 9. Strain components, von Mises strain σ_{M} (Eq. (85)), and Tresca strain σ_{T} (Eq. (86)) at the crustcore interface for and . From left to right: the dimensionless density contrast increases as labeled. 
For , the radial displacement is nearly independent of r, so Σ^{rr} ≈ 0, and the crust neither expands nor shrinks radially at any latitude. The total strain is due to the polar and azimuthal components, , whose sin^{2}θ dependence implies that the crust will first break at the equator (see Fig. 9, second panel).
The condition , for which the largest strain occurs at the poles, is satisfied for , which includes the arguably most realistic case, , in which the interface and the surface are roughly equipotentials, as for a completely fluid star. In this situation, the radial displacement increases almost linearly from the interface to the surface (lower panel of Fig. 4), so , implying that the crust compresses radially at the equator and expands at the poles. Figure 10 shows that the maximum strain, hence the first breaking point, occurs at the poles.
Fig. 10. Normalized von Mises strain on the crustcore interface at the poles (green line) and at the equator (red line) as a function of the dimensionless density contrast , with and dimensionless crust thickness . 
Figures 8 and 10 show that the maximum value of for ρ_{f} ≈ 10ρ_{s} is , obtained at the poles of the crustcore interface, where the criterion indicates that fractures will begin. Since β is small (≲10^{−4}), its precise value has no relevant influence on the shear strain. From Eq. (89), the change in the squared rotation rate required for an initially relaxed star to produce a starquake is
where σ_{c} is a generic critical strain, for whatever criterion is applied. Thus, for the high critical strains obtained by Horowitz & Kadau (2009), starquakes could only occur in stars whose initial rotation rate is a substantial fraction of Ω_{K} (see also Fattoyev et al. 2018). Assuming full relaxation in each starquake, at most a few such events could occur along the lifetime of the star, much less frequently than observed glitches. Using Eq. (81), the fractional size of the associated glitches would be
within the range of the observed glitches (Lyne et al. 2000; Espinoza et al. 2011). Of course, the frequency of the glitches could be increased by assuming either a smaller value of σ_{c} or only partial relaxation in each starquake, but this would make the size of the associated glitches proportionally smaller, so the glitch activity obtained in Eq. (84) is kept constant.
3.5. Energy released by a starquake
As shown by Baym & Pines (1971), the changes in the kinetic and gravitational potential energy in a starquake nearly cancel each other, so their sum is much smaller than the stored strain energy (roughly by a factor of b, as given in Eq. (66)). Thus, nearly all of the energy released as the crust breaks comes from the strain energy,
which can be written as
with shown in Fig. 11. Thus, we can approximate the maximum strain energy released when the crust breaks by rotation as
Fig. 11. Elastic energy released when the crust breaks as a function of the density contrast . It is normalized by the shear modulus μ, crust volume V_{c}, and critical von Mises strain σ_{M, c}. 
If all of this energy is deposited as heat into an initially cold star, its final temperature T_{f} is given by , where C(T) is the star’s heat capacity. Considering only nonsuperfluid neutrons, C(T)≈1.2 × 10^{39}(T/10^{9} K) erg K^{−1} (with some dependence on the stellar parameters and the assumed equation of state; Ofengeim et al. 2017), and the resulting interior temperature is T_{f} ≈ 1.7 × 10^{8}(E_{strain}/10^{46} erg)^{1/2} K. If, at the opposite extreme, all the neutrons and protons are paired and far below their respective superfluid or superconducting transition temperatures, the heat capacity is dominated by the electrons. In this case, the heat capacity is reduced by a factor of ∼8, and the final interior temperature is increased to T_{f} ≈ 3 × 10^{8}(E_{strain}/10^{46} erg)^{1/2} K. These temperatures are in the range where photon luminosity takes over from neutrino luminosity as the main cooling mechanism of neutron stars, which is reached at an age ∼10^{5} yr in standard cooling models (e.g., Potekhin et al. 2015). In younger, hotter neutron stars, the heating due to a starquake is not noticeable, but in older, cooler ones (such as millisecond pulsars) it could have a substantial effect lasting ∼10^{5} yr for each starquake.
It is interesting to compare the energy released in such a starquake with the energy released in the classical model of a glitch due to the sudden coupling transfer of angular momentum from the internal superfluid, of moment of inertia I_{s} and initial angular velocity Ω_{s}, to the observable crust and the rest of the star, which rotate together, with moment of inertia I_{c} and initial angular velocity Ω_{c} < Ω_{s}. The total rotational energy of both components before the glitch is , whereas after the glitch both components rotate together, with angular velocity Ω_{f} = (I_{c}Ω_{c} + I_{s}Ω_{s})/(I_{c} + I_{s}) and rotational energy . Thus, for a glitch of a given observed size ΔΩ_{gl} = Ω_{f} − Ω_{c}, the energy released is
which can be compared to the energy of a starquake of the same size, that is, with ΔΩ_{bc} = ΔΩ_{gl}. For this purpose, we use Eqs. (81) and (89) to relate the critical strain and the change in angular velocity in a starquake,
With this relation, the ratio of Eqs. (99) and (97) becomes
where we used , I_{c}/I_{s} ∼ 10^{2}, μV_{c}/E_{g}∼10^{−5}, and E_{rot}/E_{g}∼(Ω/Ω_{K})^{2}. We see that for a given change in the rotation rate, the energy released in a starquake (crust yielding) is much larger than in the classical glitch model (superfluidcrust coupling).
On the other hand, we note that a release of even a very small fraction of the strain energy (∼10^{42} erg) could heat the neutron star crust enough to induce a glitch by increasing the mobility of the superfluid neutron vortices (Link & Epstein 1996; Larson & Link 2002). Thus, although starquakes by themselves are far from able to explain all the observed glitches (or even the small ones), it might be that starquakes trigger glitches by allowing neutron vortices to move. Since small starquakes releasing only ∼10^{−4} times the total strain energy are sufficient to cause this effect, they could occur ∼10^{4} times more frequently than estimated above, coming close to explaining the total observed glitch activity.
4. Ellipticity and gravitational wave emission
As another application of the calculations presented in the previous sections, we note that rotating neutron stars with a deformation that is not axially symmetric around the axis of rotation emit continuous gravitational waves. Equations (60) and (95) yield an explicit expression for ΔI_{bc}/I, a measure of the largest deformation that can be sustained by crustal stresses (and thus released in a starquake) in our model. Of course, this deformation is axially symmetric around the rotation axis (z). However, the same deformation could in principle exist for a star that is rotating around an orthogonal axis, say, the xaxis. This star would emit gravitational waves of an amplitude proportional to the standard ellipticity,
(Jaranowski et al. 1998), where I^{jk} are the Cartesian components of the inertia tensor. Defining ΔI^{jk} ≡ I^{jk} − I_{0}, where I_{0} is the moment of inertia of the undeformed (spherical) star, it can be shown that , with ΔI^{zz} ≡ ΔI_{bc}, so
In concordance with previous work (e.g., Gittins et al. 2021 and references therein), we find this theoretical upper limit to be much higher than observational constraints for millisecond pulsars, ϵ ≲ 10^{−8} from the nondetection of gravitational waves (Abbott et al. 2020) and in some cases as low as 10^{−9} from the small measured slowdown of their rotation (e.g., Woan et al. 2018). Thus, we confirm that the strength of their crust could allow millisecond pulsars to be much less axially symmetric than observed.
5. Conclusions
We have studied the possibility of a sudden breaking of the solid crust of a neutron star (“starquake”) as it is stressed by a decreasing rotation rate. In order to do this study in a mostly analytical, parametrized fashion, we considered a Newtonian model in which both the fluid core and the solid crust have uniform densities (ρ_{f} and ρ_{s}, respectively), whose difference is parametrized by , mimicking the density gradient in a real neutron star. Additional important parameters are and , where ΔR is the crust thickness, R_{⋆} is the stellar radius, M_{⋆} is the stellar radius, G is the gravitational constant, and μ is the crustal shear modulus (also assumed uniform). The conclusions obtained from this model are the following:

The displacement field and shear strain in the crust depend quantitatively and qualitatively on . If (which includes the case , previously studied by Franco et al. 2000), the radial displacement has a maximum at an intermediate radius inside the star, beyond which it decreases toward the surface, and the shape of the crustcore interface does not approach an equipotential in the limit μ → 0. In the opposite limit, , the radial displacement increases monotonically toward the surface, and the corecrust interface does approach an equipotential shape as μ → 0 (as expected in a fluid star).

Once the shear strain at some point in the crust reaches a critical value σ_{c}, the crust breaks, its moment of inertia is suddenly reduced, and an abrupt spinup occurs, qualitatively similar to the observed pulsar glitches. The size of the glitches caused by this mechanism is proportional to σ_{c}, but their frequency is proportional to 1/σ_{c}. Thus, the resulting glitch activity (sum of glitch sizes over a fixed time interval covering many glitches) is independent of the uncertain value of σ_{c} and turns out to be much smaller than the observed glitch activity in pulsars, even if only the subclass of small glitches is considered. Thus, even this subclass cannot be accounted for by the “starquake” model.

For the large values σ_{c} ∼ 0.1 found in molecular dynamics simulations of the neutron star crust and assuming that each starquake releases most of the strain energy accumulated since the previous event (or since the birth of the star), the associated glitches would be of similar size as their observed counterparts, but they would be extremely rare, happening at most a few times in the life of a neutron star born with a fast initial rotation. Smaller values of σ_{c} or partial relaxation of the crust could account for more frequent but proportionally smaller glitches.

For values of in the range ∼10−30 suggested by the average densities of the core and crust in stateoftheart neutron star models, the “rigidity parameter” b, namely the ratio between the change of the moment of inertia during a starquake and its change between starquakes, is about an order of magnitude larger than for .

If σ_{c} ∼ 0.1 and a starquake leads to complete relaxation of the crustal stresses, the energy released can heat an initially cold neutron star up to an internal temperature of a few times 10^{8} K, comparable to the value where thermal photon emission replaces neutrino emission as the main cooling mechanism, at an age ∼10^{5} yr. In older neutron stars, such as millisecond pulsars, this would be an appreciable effect.

Even much smaller starquakes could heat the crust enough to cause a large glitch by increasing the mobility of superfluid neutron vortices (e.g., Link & Epstein 1996; Larson & Link 2002), leaving open the possibility that glitches could be triggered by starquakes.

In our model, the maximum neutron star ellipticity that can be sustained by crustal stresses is ϵ_{max} ∼ 10^{−6}(σ_{c}/0.1), much larger than the constraints ∼10^{−8} − 10^{−9} obtained from timing or gravitationalwave constraints for millisecond pulsars. Thus, we confirm that these stars are much more axially symmetric than required by theory.
We stress that our model is quite simplified, assuming that both the core and the crust of the neutron star are uniform, only roughly correcting for gravitational potential perturbations, and not including the effects of general relativity. However, none of these should change our main conclusion, namely that the observed glitches are not starquakes, although even small quakes could release enough energy to trigger thermally induced glitches.
For the stress tensor, we follow the sign convention of Thorne & Blandford (2017), which is opposite to that of Franco et al. (2000) and Cutler et al. (2003).
Throughout this paper, we express tensor components in a coordinate basis (Schutz 2009) (in the following sections specifically spherical coordinates), distinguishing between “covariant” and “contravariant” components (with indices as subscripts or superscripts), and interpreting ∇_{k} as a covariant derivative with respect to the coordinate x^{k}. We also use the “Einstein convention”, implicitly assuming a sum over indices that appear as subscripts and superscripts in the same expression, here k and l.
It cannot in general be set to zero, as done in Eq. (31) of Cutler et al. (2003), where it should in fact be an arbitrary function of r. However, this does not appear to affect the equations actually solved by those authors, presented in their Sects. 3.1.2 and 3.1.3.
Clearly δp must in general be nonzero and nonspherical, contrary to statements in Sect. 3.1.1 of Cutler et al. (2003). However, these do not appear to affect the equations actually solved in that paper and thus its final results (Cutler & Link, priv. comm.).
Since ΔT^{rr} = δT^{rr} − ρgξ^{r}, these conditions are equivalent to imposing continuity of the Eulerian stress perturbation (as done in Ushomirsky et al. 2000 and Cutler et al. 2003) if the density is continuous across the respective boundary, but not if there is a density discontinuity, as in the general case considered here.
We note that β = (c_{t}/v_{K})^{2}, where is the (roughly constant) speed of transverse shear waves in the elastic crust and is the Keplerian rotation velocity at the stellar surface. For a star of uniform density ρ_{s}, this expression is equivalent to β = μ/(2p_{c}), where p_{c} is the pressure at the center of the star.
For ξ^{z} there is a misprint in Horowitz & Kadau (2009), but this does not affect our results. The correct expressions are given in Horowitz & Hughto (2008).
α is the “strain angle”, also used by Franco et al. (2000).
Acknowledgments
We thank Cristóbal Espinoza, J. Rafael Fuentes, and Luis Rodríguez for providing important data, as well as Enrique Cerda, Curt Cutler, Felipe Espinoza, Elia Giliberti, Charles Horowitz, and Bennett Link for useful discussions and communications. We are also grateful to the anonymous referee for comments that helped improve this paper. This work was funded by FONDECYT Regular Project 1201582.
References
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJ, 902, L21 [CrossRef] [Google Scholar]
 Anderson, P. W., & Itoh, N. 1975, Nature, 256, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Baiko, D. A., & Chugunov, A. I. 2018, MNRAS, 480, 5511 [NASA ADS] [CrossRef] [Google Scholar]
 Baym, G., & Pines, D. 1971, Ann. Phys., 66, 816 [NASA ADS] [CrossRef] [Google Scholar]
 Baym, G., Pethick, C., Pines, D., & Ruderman, M. 1969, Nature, 224, 872 [NASA ADS] [CrossRef] [Google Scholar]
 Boynton, P. E., Groth, E. J. I., Partridge, R. B., & Wilkinson, D. T. 1969, Int. Astron. Union Circ., 2179, 1 [Google Scholar]
 Carter, B., & Quintana, H. 1975, Ann. Phys., 95, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Christensen, R. 2013, The Theory of Materials Failure (OUP Oxford), https://books.google.es/books?id=7pdoAgAAQBAJ [CrossRef] [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Cutler, C., Ushomirsky, G., & Link, B. 2003, ApJ, 588, 975 [NASA ADS] [CrossRef] [Google Scholar]
 Espinoza, C. M., Lyne, A. G., Stappers, B. W., & Kramer, M. 2011, MNRAS, 414, 1679 [NASA ADS] [CrossRef] [Google Scholar]
 Fattoyev, F. J., Horowitz, C. J., & Lu, H. 2018, ArXiv eprints [arXiv:1804.04952] [Google Scholar]
 Finn, L. S. 1987, MNRAS, 227, 265 [NASA ADS] [Google Scholar]
 Franco, L. M., Link, B., & Epstein, R. I. 2000, ApJ, 543, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Fuentes, J. R., Espinoza, C. M., Reisenegger, A., et al. 2017, A&A, 608, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giliberti, E., Antonelli, M., Cambiotti, G., & Pizzochero, P. M. 2019, PASA, 36, e036 [NASA ADS] [CrossRef] [Google Scholar]
 Giliberti, E., Cambiotti, G., Antonelli, M., & Pizzochero, P. M. 2020, MNRAS, 491, 1064 [NASA ADS] [Google Scholar]
 Gittins, F., Andersson, N., & Jones, D. I. 2021, MNRAS, 500, 5570 [Google Scholar]
 Gusakov, M. E., Kantor, E. M., & Reisenegger, A. 2015, MNRAS, 453, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., & Collins, R. A. 1968, Nature, 217, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Hoffman, K., & Heyl, J. 2012, MNRAS, 426, 2404 [NASA ADS] [CrossRef] [Google Scholar]
 Hooker, J., Newton, W. G., & Li, B.A. 2015, MNRAS, 449, 3559 [NASA ADS] [CrossRef] [Google Scholar]
 Horowitz, C. J., & Hughto, J. 2008, ArXiv eprints [arXiv:0812.2650] [Google Scholar]
 Horowitz, C. J., & Kadau, K. 2009, Phys. Rev. Lett., 102, 191102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Howitt, G., Haskell, B., & Melatos, A. 2016, MNRAS, 460, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Jaranowski, P., Królak, A., & Schutz, B. F. 1998, Phys. Rev. D, 58, 063001 [Google Scholar]
 JohnsonMcDaniel, N. K., & Owen, B. J. 2013, Phys. Rev. D, 88, 044004 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, D. I., Ashton, G., & Prix, R. 2017, Phys. Rev. Lett., 118, 261101 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, M. B., & Link, B. 2002, MNRAS, 333, 613 [CrossRef] [Google Scholar]
 Link, B., & Epstein, R. I. 1996, ApJ, 457, 844 [NASA ADS] [CrossRef] [Google Scholar]
 Link, B., Franco, L. M., & Epstein, R. I. 1998, ApJ, 508, 838 [NASA ADS] [CrossRef] [Google Scholar]
 Link, B., Epstein, R. I., & Lattimer, J. M. 1999, Phys. Rev. Lett., 83, 3362 [NASA ADS] [CrossRef] [Google Scholar]
 Lyne, A. G., Shemar, S. L., & Smith, F. G. 2000, MNRAS, 315, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Ofengeim, D. D., Fortin, M., Haensel, P., Yakovlev, D. G., & Zdunik, J. L. 2017, Phys. Rev. D, 96, 043002 [NASA ADS] [CrossRef] [Google Scholar]
 Pizzochero, P. M. 2011, ApJ, 743, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Pons, J. A., & Page, D. 2015, Space Sci. Rev., 191, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Ravenhall, D. G., Pethick, C. J., & Wilson, J. R. 1983, Phys. Rev. Lett., 50, 2066 [NASA ADS] [CrossRef] [Google Scholar]
 Reisenegger, A. 1995, ApJ, 442, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Reisenegger, A., & Goldreich, P. 1992, ApJ, 395, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. 1969, Nature, 223, 597 [CrossRef] [Google Scholar]
 Schutz, B. 2009, A First Course in General Relativity [CrossRef] [Google Scholar]
 Smoluchowski, R. 1970, Phys. Rev. Lett., 24, 923 [NASA ADS] [CrossRef] [Google Scholar]
 Stairs, I. H., Lyne, A. G., & Shemar, S. L. 2000, Nature, 406, 484 [NASA ADS] [CrossRef] [Google Scholar]
 Thorne, K., & Blandford, R. 2017, Modern Classical Physics: Optics, Fluids, Plasmas, Elasticity, Relativity, and Statistical Physics (Princeton University Press), https://books.google.cl/books?id=U1S6BQAAQBAJ [Google Scholar]
 Ushomirsky, G., Cutler, C., & Bildsten, L. 2000, MNRAS, 319, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Woan, G., Pitkin, M. D., Haskell, B., Jones, D. I., & Lasky, P. D. 2018, ApJ, 863, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Zdunik, J. L., Bejger, M., & Haensel, P. 2008, A&A, 491, 489 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Coefficients of the radial displacement function
The constants A_{k} in the radial displacement function (Eq. (30)) that satisfy the boundary conditions given in Eqs. (40), (41) and (42) are
where
and
depend only on the dimensionless core radius (or the complementary crust thickness ). The coefficients in Eqs. (A.5) to (A.19) are all defined so that their value is 1 in the limit of a very thin crust (), and their expansion up to linear order in is given by the approximate expressions at the end of each equation. On the other hand,
are functions of the thickness of the crust and of the dimensionless density difference .
Appendix B: Explicit expression for the rigidity parameter
The expression for the rigidity parameter given in Eq. (66) can be written more explicitly as
where is defined in Eq. (69), and
All these quantities satisfy η_{k} → 0 for (or ), and η_{k} → 1 for .
All Figures
Fig. 1. Schematic diagram of a starquake driven by spindown, with time t on the horizontal axis and rotation rate Ω on the vertical axis. In state “a”, the star is relaxed. From “a” to “b” it spins down, with its crust deforming and building up stress. At “b”, the crust reaches a critical stress and fractures, the moment of inertia decreases, and the star suddenly spins up to a new relaxed state “c”. 

In the text 
Fig. 2. Radial dependence of the normalized quadrupolar radial displacement field, , of the equipotential surfaces in a completely fluid star with uniform densities, ρ_{s} in the crust () and in the core (), with a dimensionless interface radius and different values of . 

In the text 
Fig. 3. Dimensionless radial displacement in the uniform solid star for different values of the parameter β defined in Eq. (53), compared to the case of an uniform fluid, with the displacement field in the latter case defined in terms of equipotentials, as explained in Sect. 2.3. 

In the text 
Fig. 4. Normalized radial displacements in the neutron star crust with interface radius , both for purely fluid matter (dashdotted line) and for an elastic crust with β = 10^{−4} (red dots). The upper panel corresponds to a very small density difference, , whereas the lower panel shows a star with a much denser core, . 

In the text 
Fig. 5. Normalized radial displacement at the stellar surface (r = R_{⋆}; upper panel) and at the crustcore interface (r = R_{i} = 0.9R_{⋆}; lower panel) as a function of the normalized density difference between the core and the crust. The blue line corresponds to a completely fluid star, whereas the other lines correspond to stars with an elastic crust, for different values of β. 

In the text 
Fig. 6. Rigidity parameter b/(μV_{c}/E_{g}) as a function of the difference in density between fluid core and solid crust in the horizontal axis and as a function of the core radius in the vertical axis. The fiducial parameters are fixed so that . 

In the text 
Fig. 7. Glitch activity as a function of . Black dots represent the glitch activity (sum of all glitch sizes divided by the sum of the observation times) for all pulsars within bins of , and blue dots represent the same quantity but restricted to small glitches (Δν < 10 μHz). (Data kindly provided by J. R. Fuentes, from the glitch sample discussed in Fuentes et al. 2017.) The green and red lines represent the glitch activity in the starquake model, as given by Eq. (84), for the values b = 10^{−5} and b = 10^{−6} of the rigidity parameter, corresponding to the cases and , respectively. 

In the text 
Fig. 8. Normalized von Mises strain, (see Eqs. (85) and (89)), for , , and . 

In the text 
Fig. 9. Strain components, von Mises strain σ_{M} (Eq. (85)), and Tresca strain σ_{T} (Eq. (86)) at the crustcore interface for and . From left to right: the dimensionless density contrast increases as labeled. 

In the text 
Fig. 10. Normalized von Mises strain on the crustcore interface at the poles (green line) and at the equator (red line) as a function of the dimensionless density contrast , with and dimensionless crust thickness . 

In the text 
Fig. 11. Elastic energy released when the crust breaks as a function of the density contrast . It is normalized by the shear modulus μ, crust volume V_{c}, and critical von Mises strain σ_{M, c}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.