© 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 spin-up events, in which about 1% of the secular spin-down 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 ϵ₀, which, as the star spins down, would change toward a more spherical shape by Δϵ₀ 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Δϵ₀, 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 long-term average spin-up 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

(1)

where E₀ 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 rotation-induced ellipticity ϵ, and the last term accounts for the elastic stress in the crust, where ϵ₀ 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⁻⁵ for a star of crustal shear modulus μ = 10³⁰ erg cm⁻³ 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 spin-down 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, quasi-periodic variations in the pulse shape and spin-down 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 shear-modulus profiles are obtained from equations of state based on nuclear physics calculations, yielding an even smaller value, b ≈ 2 × 10⁻⁷. 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 Johnson-McDaniel & 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¹, 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 core-crust density difference (which tends to make the crust-core 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 spin-down 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 gravitational-wave 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 spin-down, 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

(2)

(e.g., Link et al. 1999; Cutler et al. 2003; Thorne & Blandford 2017), where T is the stress tensor², ρ 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

(3)

where p is the pressure, g is the metric tensor of Euclidean, three-dimensional space (not to be confused with the gravitational acceleration g defined above), and the symmetric, trace-free 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

(4)

where

(5)

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

(6)

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

(7)

In the elastic limit, the Lagrangian stress perturbation is proportional to the strain tensor (generalized Hooke’s law)³,

(8)

where the rank-four tensor Y^ijkl is known as the Young tensor (a property of the material) and

(9)

is the strain tensor, which can be decomposed into its irreducible tensorial parts as

(10)

where the trace

(11)

represents the expansion (or compression) of the matter element, the symmetric, traceless part

(12)

represents the shear deformations, and the antisymmetric part

(13)

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

(14)

and

(15)

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

(16)

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,

(17)

where we introduced a total (effective) potential, ϕ_T ≡ ϕ_g + ϕ_c. Thus,

(18)

where the constant h needs to be obtained from the boundary conditions (which otherwise cannot be satisfied)⁴.

Choosing a spherical coordinate system (r, θ, φ) aligned with the axis of rotation, the centrifugal potential perturbation can be written as

(19)

where Δ(Ω²) is the change in the squared angular frequency, and P₀(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⁵, . Thus, motivated by the form of the centrifugal force, we assume a purely quadrupolar displacement, with

(21)

Again applying the incompressibility condition, and using the Legendre equation, the polar component takes the form

(22)

Equation (17) implies that

(23)

where the nonzero components of the strain tensor are

(24)

(25)

(26)

(27)

For the displacement field given by Eqs. (21) and (22), Eq. (23) becomes an ordinary differential equation of one independent variable,

(28)

whose solution we conveniently write as

(29)

with the dimensionless radial displacement function

(30)

(see also Franco et al. 2000; Giliberti et al. 2019). Here, we have defined the dimensionless radial coordinate, , as well as the Keplerian (or “break-up”) frequency,

(31)

for a star of radius R_⋆ and mass M_⋆. The coefficients A₁, A₃, A₋₂ and A₋₄ are dimensionless constants to be obtained from the boundary conditions. For later reference, we also write down the explicit expression for the shear potential,

(32)

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 solid-body rotation around the symmetry axis). However, in a uniform, incompressible fluid, all matter elements are equivalent and interchangeable, so there is no well-defined 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⁶δ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 spin-down 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,

(33)

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

(34)

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

(35)

whereas the quadrupole part yields

(36)

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 crust-core 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

(37)

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 crust-core 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

(38)

At the crust-core 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

(39)

at the interface radius R_i between the fluid core (labeled by f) and the solid crust (labeled by s)⁷

In our model, these boundary conditions become

(40)

(41)

(42)

Equations (40) can be rewritten as

(43)

providing two independent linear homogeneous equations relating the four constants A_k in Eq. (30),

(44)

Equation (41) is most conveniently written as

(45)

This equation has a monopolar component, yielding , and a quadrupole component,

(46)

or, equivalently,

(47)

We note that in Eqs. (46) and (47) the left-hand side is proportional to ρ_f − ρ_s and the right-hand 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

(48)

implying for the monopolar component and

(49)

or

(50)

for the quadrupolar component. In the limit μ → 0, the right-hand 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 two-domain star (crust and core with different properties), it is useful to consider the simple cases of a completely uniform (single-domain) star, both as a fluid and as a uniform elastic solid. For the two-domain 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

(51)

where is the normalized radius. For the one-domain model (the case considered by Baym & Pines 1971), regularity of Eq. (30) at r = 0 implies A₋₂ = A₋₄ = 0, which, together with Eq. (43) (imposed only at the surface) and Eq. (50), gives

(52)

with the dimensionless parameter⁸

(53)

where the numerical subscripts stand for the exponents in the reference magnitudes for neutron stars, expressed in cgs units, e.g., μ₃₀ ≡ μ/(10³⁰ erg cm⁻³). 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 uniform-density 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 core-crust interface (Eqs. (46) and (47)) implies that our model reduces to the uniform-density model of Franco et al. (2000) when (except that we use the Cowling approximation, δϕ_g = 0, whereas Franco et al. 2000 calculate δϕ_g[r, θ] self-consistently). 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 crust-core 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 (dash-dotted 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 crust-core interface (r = Ri = 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 core-crust interface yields

(54)

and at the stellar surface,

(55)

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 uniform-density case (), the exact solutions are

(56)

and

(57)

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

(58)

and at the surface as

(59)

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,

(60)

The change in the moment of inertia due to the starquake can be calculated as

(61)

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, Δ(Ω²), and the induced normalized displacement fields as

(62)

where, for the stress-free (“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

(63)

where we have neglected a term ΔΩ_bc/ΔΩ_ac ∼ b(Ω/Ω_K)², which is small in the limit of slow rotation, Ω ≪ Ω_K. We can evaluate

(64)

which, in the case without shear stresses in the crust (the purely fluid case; Eq. (37)) becomes

(65)

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,

(66)

(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:

(67)

(68)

where

(69)

is the volume of the crust,

(70)

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/(μVc/|Eg|) 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

(71)

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

(72)

where the latter approximation also assumes , making d₁ = d₂ = γ₂ = 1.

A somewhat more general expression can be obtained by taking simultaneously , , and in Eq. (B.1). This allows us to set w₀ = 1 and , keep only the lowest-order terms in and in the numerator and the denominator, and keep the lowest-order term in in each of their coefficients, which yields

(73)

Thus, as long as , we recover the result of Eq. (72), whereas for we have

(74)

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 crust-core 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 crust-core 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⁻⁵, 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

(75)

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,

(76)

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

(77)

when gravitational potential perturbations are taken into account, and

(78)

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 higher-order terms in and (q and χ², 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 (),

(79)

For the case with gravitational potential perturbations, the right-hand 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 ,

(80)

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 spin-up in a glitch to the spin-down between glitches,

(81)

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 left-hand side are measurable in observed glitches. If all glitches were equal, we could identify

(82)

where ν = Ω/(2π) is the rotation frequency, is its time-derivative (the spin-down rate of the pulsar),

(83)

is the glitch spin-up rate, customarily called “glitch activity” (Lyne et al. 2000), which accounts for the cumulative spin-up effect of glitches during the time of observation T of a given pulsar, and Δν_i is the size of its i-th glitch observed during this time. Therefore, from Eq. (81), the glitch activity in the starquake model can be written as

(84)

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 spin-up 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,

(85)

and the Tresca criterion,

(86)

where σ₁ ≥ σ₂ ≥ σ₃ 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 σ₂. This function is represented by a parabola with minimum at σ₂ = (σ₁ + σ₃)/2, at which σ_M = σ_T. On the other hand, at the extremes of the interval σ₃ ≤ σ₂ ≤ σ₁, one obtains , and therefore

(87)

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

(88)

In our model, this expression can be evaluated using the equations in Sects. 2.2 and 2.6, and rewritten as

(89)

where is a dimensionless function of position in the crust, shown in Fig. 8.

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

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⁻⁴ (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⁹ξ^x = αy/2, ξ^y = αx/2, ξ^z = 𝒪(α²), and finding that the material yields when α ≈ 0.1, corresponding to σ_M, c = σ_T, c = α/2 ≈ 0.05 according to our definitions above¹⁰. 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 off-diagonal strain components Σ^rθ are forced to vanish at both the stellar surface () and the crust-core interface (), they are generally much smaller (by a factor of ) than the diagonal components, which can be written as

(90)

(91)

(92)

Thus, ignoring the off-diagonal terms, the squared von Mises strain becomes

(93)

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 spin-down, 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 crust-core 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²θ 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 crust-core 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 crust-core interface, where the criterion indicates that fractures will begin. Since β is small (≲10⁻⁴), 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

(94)

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

(95)

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,

(96)

which can be written as

(97)

with shown in Fig. 11. Thus, we can approximate the maximum strain energy released when the crust breaks by rotation as

(98)

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 Vc, 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³⁹(T/10⁹ K) erg K⁻¹ (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⁸(E_strain/10⁴⁶ 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⁸(E_strain/10⁴⁶ 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⁵ 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⁵ 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

(99)

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,

(100)

With this relation, the ratio of Eqs. (99) and (97) becomes

(101)

where we used , I_c/I_s ∼ 10², μV_c/|E_g|∼10⁻⁵, and E_rot/|E_g|∼(Ω/Ω_K)². 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 (superfluid-crust coupling).

On the other hand, we note that a release of even a very small fraction of the strain energy (∼10⁴² 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⁻⁴ times the total strain energy are sufficient to cause this effect, they could occur ∼10⁴ 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 x-axis. This star would emit gravitational waves of an amplitude proportional to the standard ellipticity,

(102)

(Jaranowski et al. 1998), where I^jk are the Cartesian components of the inertia tensor. Defining ΔI^jk ≡ I^jk − I₀, where I₀ is the moment of inertia of the undeformed (spherical) star, it can be shown that , with ΔI^zz ≡ ΔI_bc, so

(103)

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⁻⁸ from the nondetection of gravitational waves (Abbott et al. 2020) and in some cases as low as 10⁻⁹ 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 crust-core interface does not approach an equipotential in the limit μ → 0. In the opposite limit, , the radial displacement increases monotonically toward the surface, and the core-crust 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 spin-up 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 state-of-the-art 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⁸ K, comparable to the value where thermal photon emission replaces neutrino emission as the main cooling mechanism, at an age ∼10⁵ 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⁻⁶(σ_c/0.1), much larger than the constraints ∼10⁻⁸ − 10⁻⁹ obtained from timing or gravitational-wave 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.

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

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

(A.1)

(A.2)

(A.3)

(A.4)

where

(A.5)

(A.6)

(A.7)

(A.8)

(A.9)

(A.10)

(A.11)

(A.12)

(A.13)

(A.14)

(A.15)

(A.16)

(A.17)

and

(A.18)

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,

(A.19)

(A.20)

(A.21)

(A.22)

(A.23)

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

(B.1)

where is defined in Eq. (69), and

(B.2)

(B.3)

(B.4)

All these quantities satisfy η_k → 0 for (or ), and η_k → 1 for .

All Figures

Fig. 1. Schematic diagram of a starquake driven by spin-down, 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 (dash-dotted 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 crust-core interface (r = Ri = 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/(μVc/|Eg|) 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 crust-core 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 crust-core 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 Vc, and critical von Mises strain σM, c.
In the text