text:
# The Crustal Rigidity of a Neutron Star, and Implications for PSR 1828-11 and other Precession Candidates

## 1 Introduction

Stairs et al. (2000) recently reported strong evidence for free precession in the radio pulsar PSR 1828-11. This pulsar has period $`P_s=0.405`$s, $`\dot{P}_s=6.00\times 10^{14}`$, a characteristic age of $`\tau _c=`$0.11 Myr and an inferred B-field of $`5.0\times 10^{12}`$G. Timing residuals for PSR 1828-11 show strong periodic modulation, with periods 511 and 256 days. A 1009-day periodicity is also claimed, but with lower confidence.

The 511- and 256-day timing residuals are matched by periodic modulations of the pulse shape, which strongly supports the interpretation that the modulations are due to NS precession, with a wobble amplitude of $`3^{}`$. In the model of Link & Epstein (2001), free precession sets the 511-d timescale, while coupling of the precession to the external electromagnetic torque produces variations in the NS’s spindown rate that give the main contribution to the timing residuals (due to the changing angle between the neutron star’s magnetic dipole axis and spin direction) and produce a harmonic at 256 d.

While periodic or quasi-periodic timing residuals have been observed in several other isolated pulsars (including B1642-03, the Crab, and Vela) and tentatively interpreted as evidence for torque-free precession (see Jones & Anderrson 2001 for a review), PSR 1828-11 certainly represents the most convincing case for free precession in a neutron star (NS). This pulsar has prompted us to re-examine theoretical estimates of a NS’s (body-frame) free-precession frequency, $`\nu _p`$. If the crust is essentially relaxed, the body-frame precession period is $`P_p=(\nu _s)^1I_C/\mathrm{\Delta }I_d`$, where $`\nu _s`$ is the spin frequency, $`I_C`$ is the moment of inertia of the crust (plus any portion of the NS strongly coupled to the crust on a timescale $`<\nu _s^1`$), and $`\mathrm{\Delta }I_d`$ is the residual oblateness the NS would have if it were spun down to zero frequency without the crust breaking or otherwise relaxing (Pines & Shaham 1972, Munk & MacDonald 1960). For a crust that is relaxed at its current spin rate, we can write

$$\mathrm{\Delta }I_d=b\mathrm{\Delta }I_\mathrm{\Omega },$$
(1)
where $`\mathrm{\Delta }I_\mathrm{\Omega }`$ is the piece of NS’s inertia tensor due to its spin (i.e., to centrifugal force), approximately given by (Jones 2000):

$$\mathrm{\Delta }I_\mathrm{\Omega }0.3I_0(\nu _s/\mathrm{kHz})^2,$$
(2)
where $`I_0`$ is the total moment of inertia for the star when non-rotating. The coefficient $`b`$ in Eq. (1) is sometimes called the NS’s “rigidity parameter” (Jones & Andersson 2000).

This paper is concerned with an accurate calculation of $`b`$. For the case of a star with constant density and constant shear modulus, the result is $`b=\frac{57}{10}\mu V/|E_g|`$, where $`\mu `$ is the shear modulus, $`V`$ is the star’s volume, and $`E_g=\frac{3}{5}GM^2/R`$ is star’s gravitational binding energy (Love 1944). As an estimate for a neutron star with a liquid core, Baym & Pines (1971) took:

$$b\frac{57}{10}|E_g|^1\mu 𝑑V10^5,$$
(3)
for a $`1.4M_{}`$ NS. The Baym-Pines estimate for $`b`$, Eq. (3), has been routinely adopted in the NS-precession literature.

In this paper, we mainly do two things. First, we calculate $`b`$ for realistic NS structure – a solid crust afloat on a liquid core. We solve for the strain field that develops as the NS spins down, and find that $`b`$ is smaller than found by Baym and Pines by a factor of $`40`$. A partial explanation of why Eq. (3) yields so large an overestimate is given in the Appendices.

Second, we use our improved $`b`$ value to show that PSR 1828-11 is precessing $`250`$ times faster than one would expect, assuming its crust is nearly relaxed. Under the hypothesis that crust rigidity is the dominant source of the deformation bulge $`\mathrm{\Delta }I_d`$ (and that the observed modulations are indeed due to precession), we conclude that the crust of PSR 1828-11 is not currently relaxed, and that its reference spin $`\nu _{s,ref}`$ (roughly, the spin at which the crust is most relaxed) is approximately $`40`$Hz; i.e., $`\sqrt{250}`$ times higher than the current spin. (In this sense, the precession of PSR B1828-11 is different from the Chandler wobble of the Earth, since the Earth is almost relaxed).

As we will explain in §5, the current strain tensor in the crust $`\sigma _{ab}`$ can be written as the sum of “spindown” (sd) and “reference” (ref) pieces: $`\sigma _{ab}=\sigma _{ab}^{sd}+\sigma _{ab}^{ref}`$. The $`\sigma _{ab}^{ref}`$ piece is basically the strain that would exist in the crust, even if it were spinning at $`\nu _{s,ref}`$; $`\sigma _{ab}^{ref}`$ depends on the detailed history of local crust quakes/relaxations, and is thus unknown. The spindown strain $`\sigma _{ab}^{sd}`$ is the strain induced as the NS spins down from $`\nu _{ref}`$ to its current value $`\nu _s`$. For PSR 1828-11, we show that its weighted-average is $`\overline{\sigma }_{\mathrm{ave}}^{sd}5\times 10^5`$.

Although $`\sigma _{ab}^{ref}`$ is unknown, it seems likely to us that the full strain is larger, on average, than the spindown strain. (This is because the full strain $`\overline{\sigma }\sqrt{0.5\sigma _{ab}\sigma ^{ab}}`$ contains contributions from all harmonics, while the spindown strain is entirely in the $`Y_{20}`$ harmonic.) *I*f this is true, then the NS’s crustal breaking strain $`\overline{\sigma }_{max}`$ must be larger than the average spindown value: $`\overline{\sigma }_{max}>\overline{\sigma }_{\mathrm{ave}}^{sd}5\times 10^5`$. Although $`\overline{\sigma }_{max}`$ is poorly constrained, both empirically and theoretically, a value this large seems quite reasonable: the usual assumption has been that $`\overline{\sigma }_{max}`$ for NS crusts is somewhere in the range $`10^510^1`$. (Unlike the crust’s shear modulus, which can be estimated from simple energetics, the breaking strain $`\overline{\sigma }_{max}`$ depends on the type, density, and propagation of dislocations in the lattice, which are highly uncertain.)

If significant pinning of superfluid vortex lines to the nuclear lattice or to flux tubes in the core occurs, the precession frequency is not set by the material properties of the crust (Shaham 1977). In this situation, the pinned superfluid behaves effectively as a large $`\mathrm{\Delta }I_d`$, equal to the moment of inertia $`I_p`$ of the pinned superfluid. For example, if the entire crust superfluid pins to the solid, $`I_p`$ is $`0.01I_0I_C`$. The (body-frame) precession frequency should then be of order the spin frequency (or at most $`100`$ times slower, if the crust-core coupling is much larger than usually presumed, so $`I_CI_0`$). Instead, in PSR 1828-11 and the other precession candidates reviewed in §6, the precession frequency is $`68`$ orders of magnitude lower than the spin frequency. Thus claims for long-period NS precession appeared at first to conflict with one explanation for large pulsar glitches, wherein superfluid vortices do pin to crustal nuclei, and large glitches represent large-scale unpinning events (see, e.g., Anderson & Itoh 1975). This conflict was recently resolved by Link & Cutler (2001), who showed that, in PSR 1828-11 and other precession candidates, the precessional motion itself exerts Magnus forces on the vortices that are $`>100`$ times larger than the Magnus forces responsible for giant glitches. (The fact that sufficiently large, precession-induced Magnus forces would unpin superfluid vortices had been noted previously by Shaham 1986). That is, in all likelihood the precessional motion itself immediately unpins the superfluid vortices, and keeps them unpinned. Accordingly, in this paper, we set $`I_p=0`$. (Of course, pinning can persist in precessing stars if the wobble angle is sufficiently small.)

This paper is concerned with whether crustal stresses alone can give a NS precession period of order a year (for spin rates of a few Hz), but we mention that magnetic field stresses can also give a precession period of this order if the core is a Type II supeconductor, as shown recently by Wasserman (2002), in an extension of work by Mestel and collaborators (Mestel & Takhar 1972; Mestel et al. 1981; Nittman & Wood 1981; see all Spitzer 1958).

The plan of this paper is as follows. In §2 we review the dynamics of NS precession, partly to establish notation. In §3 we present our formalism for calculating the rigidity parameter $`b`$ for a NS. In §4 we present the results of this calculation: the value of $`b`$ and the corresponding crustal stresses. Knowing $`b`$ allows us to estimate the precession frequency of a NS with relaxed crust. In §5 we extend this result to a NS with substantially strained crust. We apply these results to PSR 1828-11 in §6, where we show that the observed precession frequency is $`250`$ times higher than would be predicted for a relaxed crust. We calculate the spindown strain levels that must exist today in the crust of PSR 1828-11, assuming that crustal shear stresses are indeed responsible for the observed precession frequency.

In §7 we briefly discuss the implications of our work for some other NS’s that show evidence for precession (the Crab, Vela, PSR B1642-03, PSR 2217+47, PSR B0959-54, and the possible remnant in SN 1987A). Our conclusions are summarized in §8.

The Appendices represent our attempts to check and understand our results for $`b`$ – in particular to understand why the Baym-Pines estimate is so much larger than the value we obtain. In Appendix A we represent $`\mathrm{\Delta }I_d`$ as an integral over crustal shear stresses. The result for $`b`$ we obtain in this way is consistent with our result in §4, and gives more insight into how how different stress components and different layers in the crust contribute to the final answer. In Appendix B we treat a case where $`b`$ can be obtained analytically– that of a uniform-density, incompressible star with a thin crust of constant shear modulus afloat on a liquid core. For this case, we show the actual $`b`$ is $`5`$ times smaller than the Baym-Pines estimate, Eq. (3). In the realistic NS case, $`b`$ is decreased by additional factor $`8`$, because the contributions to $`\mathrm{\Delta }I_d`$ from the different stress components and from different layers in the crust tend to cancel each other much more nearly than in the case of a thin, uniform crust.

## 2 Model and Results

Here we derive the basic equations describing the precession of a NS. In this section, and in §§3-6, we will assume that the only sources of non-sphericity in the NS are centrifugal force and crustal shear stresses.

We idealize the precessing NS as having two components: C (for crust) and F (for fluid), having angular velocities, $`\mathrm{\Omega }_C^a`$ and $`\mathrm{\Omega }_F^a`$, respectively. $`\mathrm{\Omega }_C^a`$ and $`\mathrm{\Omega }_F^a`$ have the same magnitude, $`\mathrm{\Omega }`$, but different directions. The $`C`$ piece includes the crustal nuclei, but may also include some component of the fluid that is coupled to the crust sufficiently strongly that it effectively co-rotates with the crust over timescales shorter than the rotation period $`2\pi \mathrm{\Omega }^1`$. The F piece is that portion of the fluid that is essentially decoupled from the crustal nuclei on this timescale; conservation of vorticity then implies that the fluid angular velocity $`\mathrm{\Omega }_F^a`$ remains fixed, and so is aligned with the NS’s total angular momentum $`J^a`$:

$$\mathrm{\Omega }_F^a=\mathrm{\Omega }\widehat{J}^a,$$
(4)
where $`\widehat{J}^a`$ is the unit vector along $`J^a`$. The angular velocity of the crustal nuclei is $`\mathrm{\Omega }_C^a=\mathrm{\Omega }\widehat{\mathrm{\Omega }}_C^a`$, where $`\widehat{\mathrm{\Omega }}_C^a`$ is a unit vector.

We write the star’s inertia tensor as

$$I_{ab}=I_{F,ab}+I_{C,ab},$$
(5)
where $`I_{F,ab}`$ and $`I_{C,ab}`$ are the contributions from the fluid interior and the solid crust, respectively. Following (and somewhat extending) Pines and Shaham (1972), we approximate each contribution as the sum of a spherical piece, plus a centrifugal bulge that follows that component’s angular velocity, plus a deformation bulge (sustained by crustal rigidity) that follows some principal crustal axis $`\widehat{n}^a`$:

$`I_{F,ab}`$ $`=`$ $`I_{F,0}\delta _{ab}+\mathrm{\Delta }I_{F,\mathrm{\Omega }}(\widehat{J}^a\widehat{J}^b{\displaystyle \frac{1}{3}}\delta _{ab})`$
$`+`$ $`\mathrm{\Delta }I_{F,d}(\widehat{n}_d^a\widehat{n}_d^b{\displaystyle \frac{1}{3}}\delta _{ab})`$
$`I_{C,ab}`$ $`=`$ $`I_{C,0}\delta _{ab}+\mathrm{\Delta }I_{C,\mathrm{\Omega }}(\widehat{\mathrm{\Omega }}_C^a\widehat{\mathrm{\Omega }}_C^b{\displaystyle \frac{1}{3}}\delta _{ab})`$
$`+`$ $`\mathrm{\Delta }I_{F,d}(\widehat{n}_d^a\widehat{n}_d^b{\displaystyle \frac{1}{3}}\delta _{ab}).`$
where $`\delta _{ab}`$ is the unit tensor $`[1,1,1]`$. The sum $`\mathrm{\Delta }I_d\mathrm{\Delta }I_{C,d}+\mathrm{\Delta }I_{F,d}`$ we call the NS’s deformation bulge; it is the non-sphericity the body would retain if it were slowed down to zero angular velocity, without the crust breaking or otherwise relaxing. Although $`\mathrm{\Delta }I_d`$ is ultimately due to crustal shear stresses, the term $`\mathrm{\Delta }I_{F,d}`$ is non-zero (even in the absence of pinning), because the crust exerts force on the fluid, both gravitationally and via pressure.

For a fully relaxed, rotating crust, we expect $`\mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }`$, where $`\mathrm{\Delta }I_\mathrm{\Omega }\mathrm{\Delta }I_{F,\mathrm{\Omega }}+\mathrm{\Delta }I_{C,\mathrm{\Omega }}`$ is the centrifugal piece of the entire star’s inertia tensor, and the coefficient of proportionality $`b\mathrm{\Delta }I_d/\mathrm{\Delta }I_\mathrm{\Omega }`$ is the rigidity parameter. For a physical interpretation of $`b`$, let us take the Earth as an example. The Earth’s crust is essentially relaxed. If we could slow the Earth down to zero angular velocity without cracking its crust, it would not settle directly into a spherical shape, but rather would remain somewhat oblate. This is because the Earth’s relaxed, zero-strain shape is oblate, and after centrifugal forces are removed, the stresses that build up in the crust will tend to push it back toward that relaxed shape. In fact, for the Earth, $`b0.7`$ (Munk and MacDonald 1960, p.40), so stopping the Earth from spinning would reduce its oblateness by only $`30\%`$. For a NS, as we will see, $`b2\times 10^7`$, so halting a NS’s rotation would decrease its oblateness by a factor of five million. The quantity $`b`$ is so small for a NS because the gravitational energy density far exceeds the crust’s shear modulus.

The total angular momentum $`J^a`$ of fluid plus crust is

$$J_a=\mathrm{\Omega }\left(I_{F,ab}\widehat{J}^b+I_{C,ab}\widehat{\mathrm{\Omega }}_C^b\right)$$
(8)
Note from Eqs. (2)-(8) that $`\widehat{\mathrm{\Omega }}_C^a`$, $`\widehat{n}^a`$, and $`\widehat{J}^a`$ are all coplanar. Let $`J^a`$ lie along the z-axis, and let $`\widehat{\mathrm{\Omega }}_C^a`$ and $`\widehat{n}^a`$ lie, at some instant, in the $`xz`$ plane, as illustrated in Figure 1. The angle between the crust’s symmetry axis and angular velocity is $`\theta `$. The angle between the angular velocity and the angular momentum is $`\stackrel{~}{\theta }`$. These angles are constants of the motion, and satisfy:

$`\widehat{n}^a=\mathrm{cos}\theta \widehat{z}^a\mathrm{sin}\theta \widehat{x}^a`$ (9)
$`\widehat{\mathrm{\Omega }}_C^a=\mathrm{cos}\stackrel{~}{\theta }\widehat{z}^a+\mathrm{sin}\stackrel{~}{\theta }\widehat{x}^a,`$ (10)
From Eqs. (9)-(10) we have, up to third-order terms in $`\theta `$ and $`\stackrel{~}{\theta }`$,

$$\widehat{\mathrm{\Omega }}_C^a=\mathrm{\Omega }\left(1+\frac{\stackrel{~}{\theta }}{\theta }\frac{1}{2}\theta ^2\frac{1}{2}\stackrel{~}{\theta }\theta \right)\widehat{J}^a\frac{\stackrel{~}{\theta }}{\theta }\mathrm{\Omega }\widehat{n}^a.$$
(11)
From Eq. (11), we see that the precession is a superposition of two motions: $`\widehat{n}^a`$ precesses around $`\widehat{J}^a`$ with (inertial frame) precession frequency $`\dot{\varphi }=\mathrm{\Omega }(1+\stackrel{~}{\theta }/\theta )`$, up to terms second order in $`\theta `$ and $`\stackrel{~}{\theta }`$. Likewise, the coefficient of $`\widehat{n}^a`$ in Eq. (11) gives the body-frame precession frequency: $`\mathrm{\Omega }_p=\mathrm{\Omega }(\stackrel{~}{\theta }/\theta )`$. For $`\mathrm{\Omega }>0`$ and $`\mathrm{\Delta }I_d>0`$, an observer situated above the pole $`\widehat{n}^a`$ and fixed in the body frame sees the angular velocity $`\widehat{\mathrm{\Omega }}^a`$ circle around $`\widehat{n}^a`$ in the counter-clockwise direction. The ratio $`\stackrel{~}{\theta }/\theta `$ is obtained directly by taking the $`x`$-component of Eq. (8); the result is

$$\stackrel{~}{\theta }/\theta =\frac{\mathrm{\Delta }I_d}{I_{C,0}+\frac{2}{3}\mathrm{\Delta }I_{C,\mathrm{\Omega }}\frac{1}{3}\mathrm{\Delta }I_{C,d}},$$
(12)
where $`\mathrm{\Delta }I_d\mathrm{\Delta }I_{C,d}+\mathrm{\Delta }I_{F,d}`$. Since $`\mathrm{\Delta }I_{C,\mathrm{\Omega }}`$ and $`\mathrm{\Delta }I_{C,d}`$ are tiny compared to $`I_{C,0}`$, we now ignore those terms and drop the subscript “0” from $`I_C`$. The (body frame) precession period $`P_p`$ is thus

$$P_p=P_sI_C/\mathrm{\Delta }I_d.$$
(13)
Since the pulsar’s magnetic dipole moment and emission region are presumably fixed with respect to the crust, the angle between the dipole direction and $`\widehat{J}^a`$ varies on the timescale $`P_p`$; i.e., the body-frame precession period $`P_p`$ is the modulation period for the pulsed radio signal. Thus for PSR1828-11, the measured value of $`\mathrm{\Delta }I_d/I_C`$ is

$$\frac{\mathrm{\Delta }I_d}{I_C}=(\nu _sP_p)^1=9.2\times 10^9\left(\frac{511\text{ d}}{P_p}\right).$$
(14)
(We have included the factor $`511d/P_p`$ explicitly in case the precession period is actually $`1009`$ or $`256`$ days.)

## 3 Elastic Deformation of the Slowing Crust: Formalism

Our goal in §3-4 is to determine the rigidity parameter $`b\mathrm{\Delta }I_d/\mathrm{\Delta }I_\mathrm{\Omega }`$ for a NS with a relaxed crust. We first show that $`\mathrm{\Delta }I_d`$ can be re-expressed as the difference between (i) $`\mathrm{\Delta }I_\mathrm{\Omega }`$ for a purely fluid NS; and (ii) $`\mathrm{\Delta }I_\mathrm{\Omega }`$ for the same NS, but now having an elastic crust whose relaxed state is spherical. We then describe how we solve for $`\mathrm{\Delta }I_\mathrm{\Omega }`$ for both cases (i) and (ii), using perturbation theory.

Consider a 2-parameter family of NS’s, all with the same mass and composition, and all spinning about the $`z`$-axis. The two parameters are the NS’s actual angular velocity, $`\mathrm{\Omega }`$, and the angular velocity $`\mathrm{\Omega }_r`$ at which its crust is relaxed. For convenience, we change variables to $`\mathrm{\Omega }_r^2`$ and $`\mathrm{\Delta }\mathrm{\Omega }^2\mathrm{\Omega }^2\mathrm{\Omega }_r^2`$. The star’s inertia tensor $`I_{ab}`$ is a function of both these variables: $`I_{ab}=I_{ab}(\mathrm{\Omega }_r^2,\mathrm{\Delta }\mathrm{\Omega }^2)`$. By axial symmetry, we can write

$$I_{ab}=I_0(\mathrm{\Omega }_r^2,\mathrm{\Delta }\mathrm{\Omega }^2)\delta _{ab}+I_2(\mathrm{\Omega }_r^2,\mathrm{\Delta }\mathrm{\Omega }^2)\left(\widehat{z}^a\widehat{z}^b\frac{1}{3}\delta _{ab}\right).$$
(15)
We define the function $`\mathrm{\Delta }I_d(\mathrm{\Omega })`$ to be the “residual” $`I_2`$ of a NS that is non-rotating, but whose crust is relaxed at angular velocity $`\mathrm{\Omega }`$. Then, by definition,

$$\mathrm{\Delta }I_d(\mathrm{\Omega })I_2(\mathrm{\Omega }^2,\mathrm{\Omega }^2),$$
(16)
since $`\mathrm{\Omega }_r^2=\mathrm{\Omega }^2`$ means that the crust is relaxed at angular velocity $`\mathrm{\Omega }`$, and $`\mathrm{\Delta }\mathrm{\Omega }^2=\mathrm{\Omega }_r^2`$ means that the star is non-rotating.

Note that $`I_2(0,0)=0`$, since $`\mathrm{\Omega }_r^2=\mathrm{\Delta }\mathrm{\Omega }^2=0`$ means both that the star’s crust is relaxed and that the star is non-rotating, and hence the star is spherical. Therefore we can write, to first order in $`\mathrm{\Omega }_r^2`$ and $`\mathrm{\Delta }\mathrm{\Omega }^2`$,

$`\mathrm{\Delta }I_d(\mathrm{\Omega })`$ $``$ $`I_2(\mathrm{\Omega }^2,\mathrm{\Omega }^2)`$ (17)
$`=`$ $`\mathrm{\Omega }^2{\displaystyle \frac{I_2}{(\mathrm{\Omega }_r^2)}}|_{0,0}\mathrm{\Omega }^2{\displaystyle \frac{I_2}{(\mathrm{\Delta }\mathrm{\Omega }^2)}}|_{0,0}`$ (18)
$`=`$ $`I_2(\mathrm{\Omega }^2,0)I_2(0,\mathrm{\Omega }^2).`$ (19)
This suggests our strategy for computing $`\mathrm{\Delta }I_d(\mathrm{\Omega })`$: we will compute both $`I_2(\mathrm{\Omega }^2,0)`$ and $`I_2(0,\mathrm{\Omega }^2)`$ and take the difference. Now, $`I_2(\mathrm{\Omega }^2,0)`$ is the $`I_2`$ of a star spinning at $`\mathrm{\Omega }`$, whose crust is completely relaxed at that spin. But if the crust is relaxed (i.e., if all shear stresses vanish), the shape of the star is the same as for a completely fluid star with the same mass and angular velocity. On the other hand, $`I_2(0,\mathrm{\Omega }^2)`$ is the $`I_2`$ of a star whose angular velocity is $`\mathrm{\Omega }`$, but whose whose crust is relaxed when spherical (i.e., when it is nonspinning).

Our computational strategy is to start with two nonspinning NS’s, (i) and (ii), identical except that (i) is treated as a fluid throughout, while (ii) has an elastic crust whose relaxed shape is spherical. Now spin them both up to angular velocity $`\mathrm{\Omega }`$, calculate $`I_2`$ for each, and take the difference. The result is $`\mathrm{\Delta }I_d`$ for a NS whose crust is relaxed at spin $`\mathrm{\Omega }`$. (More precisely, we find $`\mathrm{\Delta }I_d(\mathrm{\Omega })`$, up to fractional corrections of order $`\mathrm{\Omega }^2R^3/GM`$.)

### 3.1 The perturbation equations

The advantage to the strategy outlined above is that both problems (i) and (ii) can be treated with a linear analysis about a spherical background for the stellar structure, with centrifugal force acting as a source term that drives the star away from sphericity. The effect of centrifugal force is described by the following potential in spherical coordinates (with $`\theta =0`$ corresponding to the rotation axis):

$$\delta \varphi ^c=\frac{1}{2}\mathrm{\Omega }^2r^2\left(1\mathrm{cos}^2\theta \right).$$
(20)
This potential has both an $`l=0`$ and $`l=2,m=0`$ piece, $`\delta \varphi ^c=\delta \varphi _0^c+\delta \varphi _2^c`$, where

$$\delta \varphi _0^c=\frac{1}{3}\mathrm{\Omega }^2r^2$$
(21)
and

$$\delta \varphi _2^c(r,\theta ,\varphi )=\widehat{\lambda }Y_{20}(\theta ,\varphi )\mathrm{\Omega }^2r^2,$$
(22)
with

$$\widehat{\lambda }\frac{1}{3}\sqrt{\frac{4\pi }{5}}.$$
(23)
The $`l=0`$ piece causes spherically-symmetric expansion of the star, while the $`l=2`$ piece generates the equatorial bulge. Only the $`l=2`$ piece contributes to $`I_2`$. However, displacements due to both $`l=0`$ and $`l=2`$ pieces generate strains in the crust and, hence, are important in determining the maximum $`I_2`$ spin-change sustainable by the crust without cracking or otherwise relaxing. Previous work on spindown-induced crustal strain (e.g., Baym & Pines 1971; Franco et al. 2000) considered incompressible NS models, where this issue does not arise, as the $`l=0`$ piece of the displacement vanishes for incompressible models. We now derive the equations governing the response of the star to any $`\delta \varphi ^c`$, and then specialize to the $`l=2`$ and $`l=0`$ pieces separately. Note that our treatment will be purely Newtonian.

#### 3.1.1 NS response to centrifugal potential $`\delta \varphi ^c`$

Inside the crust, the perturbed stress tensor (the sum of pressure and shear stresses) can be written as

$$\delta \tau _{ab}=\delta pg_{ab}+2\mu \sigma _{ab},$$
(24)
where $`\delta p`$ is the pressure perturbation, $`\mu (r)`$ is the crust’s shear modulus, and $`\sigma _{ab}`$ is the crust’s strain tensor, defined by

$$\sigma _{ab}=\frac{1}{2}\left(_a\xi _b+_b\xi _a\frac{2}{3}g_{ab}^c\xi _c\right).$$
(25)
Here $`\xi ^a`$ is the displacement of the crust away from its relaxed state. The perturbation $`\delta \tau _{ab}`$ is determined by the condition of hydro-elastic balance,

$$^a\delta \tau _{ab}=\delta \rho g\widehat{r}_b+\rho _b\left(\delta \varphi +\delta \varphi _2^c\right),$$
(26)
where $`\widehat{r}_b`$ is the radial unit vector. The perturbation $`\delta \varphi `$ of the gravitational potential obeys

$$^2\delta \varphi =4\pi G\delta \rho .$$
(27)
Finally, the density change arising from the displacement $`\xi ^a`$ obeys the continuity equation,

$$\delta \rho =^a(\rho \xi _a).$$
(28)
The following bundary conditions arise from Eq. (26): 1) $`\delta \tau _r`$ (where “$``$” refers to components orthogonal to the radial vector $`\widehat{r}^a`$) must go to zero at both the crust-core boundary and the stellar surface (using the fact that shear stresses vanish just outside the crust), and 2) $`\delta \tau _{rr}=0`$ for perturbations that break spherical symmetry (since $`\delta p`$ must vanish both in the liquid core–else fluid would flow–and above the star). These conditions can be succinctly expressed as

$$\delta \tau _{ab}\widehat{r}^b=0,$$
(29)
at both boundaries. (See Ushomirsky, Cutler & Bildsten 2000 for a more detailed discussion of these boundary conditions).

In the fluid core there are only two first-order equations to solve, for $`\delta \varphi `$ and $`d\delta \varphi /dr`$. To derive them, first note that equilibrium between pressure, gravitational, and centrifugal forces requires

$$_a\delta p=\delta \rho _a\varphi \rho _a\left[\delta \varphi +\delta \varphi ^c\right].$$
(30)
Projecting Eq. (30) along the horizontal and radial directions, we have

$`\delta p(r)=\rho \left[\delta \varphi +\delta \varphi ^c\right]`$ (31)
$`{\displaystyle \frac{d\delta p}{dr}}=\delta \rho g\rho {\displaystyle \frac{d}{dr}}\left[\delta \varphi +\delta \varphi ^c\right].`$ (32)
Taking $`d/dr`$ of Eq. (31) and subtracting Eq. (32), we have

$$\delta \rho =g^1\frac{d\rho }{dr}\left[\delta \varphi +\delta \varphi ^c\right]$$
(33)
so Poisson’s equation becomes

$$^2\delta \varphi =4\pi G(g^1\frac{d\rho }{dr})\left[\delta \varphi +\delta \varphi ^c\right].$$
(34)
#### 3.1.2 $`l=2`$ deformations

We now seperate the above equations into spherical harmonics, thereby obtaining a system of first-order ODEs. Our treatment follows McDermott et al. (1988); see also Ushomirsky, Cutler, & Bildsten (2000, hereafter UCB). Throughout this subsection, it is implicit that all scalar perturbations have $`Y_{20}`$ angular dependence; we will generally suppress explicit $`(l,m)`$ subscripts on perturbed quantities.

We begin by writing the displacement vector $`\xi ^a`$ as the sum of radial and tangential pieces:

$$\xi ^a\xi ^r(r)Y_{lm}\widehat{r}^a+\xi ^{}(r)\beta ^1r^aY_{lm},$$
(35)
where $`\widehat{r}^a`$ is the unit radial vector and $`\beta \sqrt{l(l+1)}`$. Then inside the crust, there are 6 variables:

$`z_1{\displaystyle \frac{\xi ^r}{r}},`$ $`z_2{\displaystyle \frac{\mathrm{\Delta }\tau _{rr}}{p}}={\displaystyle \frac{\delta \tau _{rr}}{p}}z_1{\displaystyle \frac{d\mathrm{ln}p}{d\mathrm{ln}r}}`$ (36)
$`z_3{\displaystyle \frac{\xi ^{}}{\beta r}},`$ $`z_4{\displaystyle \frac{\mathrm{\Delta }\tau _r}{\beta p}}={\displaystyle \frac{\delta \tau _r}{\beta p}}`$
$`z_5{\displaystyle \frac{\delta \varphi }{c^2}}{\displaystyle \frac{R^2}{r^2}},`$ $`z_6{\displaystyle \frac{dz_5}{d\mathrm{ln}r}}.`$
where $`\delta `$ and $`\mathrm{\Delta }`$ refer to Eulerian and Lagrangian variations, respectively. The definitions of $`z_1`$-$`z_4`$ and $`\delta \tau _{ab}`$ are the same as in UCB. The above definition of $`z_5`$, with $`r^2`$ in the denominator, factors out the dominant behavior of $`\delta \varphi /c^2`$ near $`r=0`$.

The perturbation equations are (see McDermott et al. 1988; UCB):

$`{\displaystyle \frac{dz_1}{d\mathrm{ln}r}}`$ $`=`$ $`\left(1+2{\displaystyle \frac{\alpha _2}{\alpha _3}}\right)z_1+{\displaystyle \frac{1}{\alpha _3}}z_2+6{\displaystyle \frac{\alpha _2}{\alpha _3}}z_3,`$ (37a)
$`{\displaystyle \frac{dz_2}{d\mathrm{ln}r}}`$ $`=`$ $`\left(\stackrel{~}{U}\stackrel{~}{V}4\stackrel{~}{V}+12\mathrm{\Gamma }{\displaystyle \frac{\alpha _1}{\alpha _3}}\right)z_1+\left(\stackrel{~}{V}4{\displaystyle \frac{\alpha _1}{\alpha _3}}\right)z_2`$
$`+`$ $`6\left(\stackrel{~}{V}6\mathrm{\Gamma }{\displaystyle \frac{\alpha _1}{\alpha _3}}\right)z_3+6z_4`$
$`+`$ $`{\displaystyle \frac{\stackrel{~}{V}c_3}{Z_{}}}(2z_5+z_6)+2{\displaystyle \frac{\stackrel{~}{V}c_3}{Z_{}}}\widehat{\lambda }{\displaystyle \frac{R^2\mathrm{\Omega }^2}{c^2}},`$
$`{\displaystyle \frac{dz_3}{d\mathrm{ln}r}}`$ $`=`$ $`z_1+{\displaystyle \frac{1}{\alpha _1}}z_4,`$ (37c)
$`{\displaystyle \frac{dz_4}{d\mathrm{ln}r}}`$ $`=`$ $`\left(\stackrel{~}{V}6\mathrm{\Gamma }{\displaystyle \frac{\alpha _1}{\alpha _3}}\right)z_1{\displaystyle \frac{\alpha _2}{\alpha _3}}z_2`$
$`+`$ $`{\displaystyle \frac{2}{\alpha _3}}\left\{11\alpha _1\alpha _2+10\alpha _1^2\right\}z_3`$
$`+`$ $`\left(\stackrel{~}{V}3\right)z_4+{\displaystyle \frac{\stackrel{~}{V}c_3}{Z_{}}}z_5+{\displaystyle \frac{\stackrel{~}{V}c_3}{Z_{}}}\widehat{\lambda }{\displaystyle \frac{R^2\mathrm{\Omega }^2}{c^2}},`$
$`{\displaystyle \frac{dz_5}{d\mathrm{ln}r}}`$ $`=`$ $`z_6,`$ (37e)
$`{\displaystyle \frac{dz_6}{d\mathrm{ln}r}}`$ $`=`$ $`5z_6+{\displaystyle \frac{\stackrel{~}{U}Z_{}}{c_3}}[({\displaystyle \frac{\stackrel{~}{V}}{\gamma }}4{\displaystyle \frac{\alpha _1}{\alpha _3}})z_1`$
$``$ $`{\displaystyle \frac{1}{\alpha _3}}z_2+12{\displaystyle \frac{\alpha _1}{\alpha _3}}z_3],`$
where

$`\stackrel{~}{U}{\displaystyle \frac{d\mathrm{ln}g}{d\mathrm{ln}r}}+2,\stackrel{~}{V}{\displaystyle \frac{\rho gr}{p}}={\displaystyle \frac{d\mathrm{ln}p}{d\mathrm{ln}r}},\alpha _1\mu /p,`$ (38)
$`\mathrm{\Gamma }={\displaystyle \frac{\mathrm{ln}p}{\mathrm{ln}\rho }}|_{\mu _\mathrm{e}},\alpha _2\mathrm{\Gamma }{\displaystyle \frac{2}{3}}{\displaystyle \frac{\mu }{p}},\alpha _3\mathrm{\Gamma }+{\displaystyle \frac{4}{3}}{\displaystyle \frac{\mu }{p}},`$
$`\gamma ={\displaystyle \frac{d\mathrm{ln}p}{d\mathrm{ln}\rho }}|_{},c_3=\left({\displaystyle \frac{r}{R}}\right)^3\left({\displaystyle \frac{M}{M_r}}\right),Z_{}={\displaystyle \frac{GM}{Rc^2}}.`$
In the fluid, we rewrite Eq. (34) in terms of our dimensionless variables $`z_5`$ and $`z_6`$:

$`{\displaystyle \frac{dz_5}{d\mathrm{ln}r}}`$ $`=`$ $`z_6,`$ (39)
$`{\displaystyle \frac{dz_6}{d\mathrm{ln}r}}`$ $`=`$ $`5z_5{\displaystyle \frac{\stackrel{~}{U}\stackrel{~}{V}}{\gamma }}\left(z_5+\widehat{\lambda }{\displaystyle \frac{R^2\mathrm{\Omega }^2}{c^2}}\right).`$ (40)
The boundary conditions are as follows. At the center of the star, $`d\delta \varphi /dr`$ must vanish, giving

$$z_6=0,$$
(41)
while $`z_5`$ takes on an (unknown) finite value. At both the crust-core boundary ($`r=r_c`$) and at the top of the crust ($`r=R`$), we have, from Eq. (29),

$`z_2`$ $`=`$ $`\stackrel{~}{V}\left\{{\displaystyle \frac{\rho _l}{\rho _s}}z_1{\displaystyle \frac{c_3}{Z_{}}}\left(z_5+\widehat{\lambda }{\displaystyle \frac{R^2\mathrm{\Omega }^2}{c^2}}\right)\right\},`$ (42)
$`z_4`$ $`=`$ $`0.`$
In addition, at the top of the crust ($`r=R`$),

$$5z_5+z_6=0.$$
(43)
Finally, at the crust-core boundary $`\delta \varphi `$ and $`d\delta \varphi /dr`$ must be continuous, requiring that $`z_5`$ and $`z_6`$ are as well. We thus have 8 equations (two in the core and six in the crust) and 8 boundary conditions (one at the center, three at the surface, and four at the crust-core boundary). For simplicity, we set the outer boundary of the NS at the top of the crust. For convenience, we approximate $`\mathrm{\Gamma }`$ as $`\gamma `$ in Eqs. (37) (i.e., we do not distinguish $`d\mathrm{ln}p/d\mathrm{ln}\rho `$ of the background model from the derivative at constant composition). Solutions of these equations are described in § 4.

The $`I_2`$ of the NS’s inertia tensor is related to the NS’s quadrupole moment $`Q_{20}`$ by

$$I_2=\sqrt{\frac{4\pi }{5}}Q_{20},$$
(44)
where $`Q_{20}`$ is defined by

$$Q_{20}\delta \rho (r)r^4𝑑r.$$
(45)
$`Q_{20}`$ is related to $`z_5`$ at the NS surface by

$$Q_{20}=\frac{5c^2}{4\pi G}z_5|_{r=R}.$$
(46)
As a check on the accuracy of our solutions, we also calculate $`Q_{20}`$ directly from the following integrals. In the crust, $`\delta \rho =(\rho \stackrel{}{\xi })`$, so

$`Q_{20}^{\mathrm{crust}}`$ $`=`$ $`2{\displaystyle _{r_c}^R}\rho (z_1+3z_3)r^4𝑑r\left[r^5\rho z_1\right]_{r=r_c}^{r=R}`$
$`=`$ $`{\displaystyle \frac{MR^2}{4\pi }}\{2{\displaystyle _{r_c}^R}{\displaystyle \frac{\stackrel{~}{U}}{c_3}}(z_1+3z_3)\left({\displaystyle \frac{r}{R}}\right)^5d\mathrm{ln}r`$
$``$ $`\left[\left({\displaystyle \frac{r}{R}}\right)^5{\displaystyle \frac{\stackrel{~}{U}}{c_3}}z_1\right]_{r=r_c}^{r=R}\}.`$ (48)
In the liquid part of the star, we use Eq. (33) to obtain

$$Q_{20}^{\mathrm{fluid}}=\frac{MR^2}{4\pi Z_{}}_0^{r_c}\frac{\stackrel{~}{U}\stackrel{~}{V}}{\gamma }\left(z_5+\widehat{\lambda }\frac{R^2\mathrm{\Omega }^2}{c^2}\right)\left(\frac{r}{R}\right)^5d\mathrm{ln}r.$$
(49)
In a star with both crust and core, we take $`Q_{20}=Q_{20}^{\mathrm{fluid}}+Q_{20}^{\mathrm{crust}}`$. In practice, we refine our numerical solution until the sum of (48) and (49) agrees with (46) to within required tolerance–about one part in $`10^{10}`$.

To find $`Q_{20}`$ for a purely fluid star, we solve the same fluid equations, but instead of matching at the bottom of the crust, we impose the boundary condition Eq. (43) at $`r=R`$. Neither solution (i) nor (ii) is really accurate to one part in $`10^{10}`$, because we approximate the star’s outer boundary as being at the top of the crust (ignoring the thin ocean on top), and because we approximate $`\mathrm{\Gamma }`$ by $`\gamma `$ in Eqs. (37). However since we use the same approximations in the solutions to both problems (i) and (ii), the difference is due solely to the presence of an elastic crust in model (ii). Therefore we claim we do know difference between (i) and (ii) to three significant figures despite the fact that the two quantities agree to seven significant figures. This claim is checked in Appendix A, where we obtain an approximate analytical expression for $`b`$.

#### 3.1.3 $`l=0`$ deformations

The $`l=0`$ component of the centrifugal potential, $`\delta \varphi _0^c`$, generates a spherically symmetric radial displacement field $`\xi ^a=\xi ^r(r)\widehat{r}^a`$, as well as pressure, density, and gravitational potential perturbations. In this subsection, all scalar perturbed quantities are purely functions of $`r`$.

In order to compute $`\xi ^r`$ for the $`l=0`$ case, we neglect the small shear modulus of the crust (effectively taking the entire star to be fluid). Of course, it is the small shear modulus that determines $`\mathrm{\Delta }I_d`$, and so in the previous subsection our equations included it. However, our present purpose is only to determine (roughly) the $`l=0`$ part of the crustal strain resulting from NS spindown, and for this we can neglect $`\mu `$, since it has only a small influence on the radial displacement. This approximation is justified by the smallness of the crust’s shear modulus compared to the pressure. With this approximation, the perturbations in pressure, density, and gravitational potential obey

$`{\displaystyle \frac{d\delta p}{dr}}`$ $`=`$ $`\delta \rho g\rho {\displaystyle \frac{d}{dr}}\left(\delta \varphi +\delta \varphi _0^c\right)`$ (50)
$`{\displaystyle \frac{\delta p}{p}}`$ $`=`$ $`\mathrm{\Gamma }{\displaystyle \frac{\delta \rho }{\rho }}+\left(\gamma \mathrm{\Gamma }\right)\xi ^r{\displaystyle \frac{d\mathrm{ln}\rho }{dr}}`$ (51)
$`^2\delta \varphi `$ $`=`$ $`4\pi G\delta \rho .`$ (52)
Eq. (51) comes from setting $`\frac{\mathrm{\Delta }p}{p}=\mathrm{\Gamma }\frac{\mathrm{\Delta }\rho }{\rho }`$.

We define dimensionless variables

$`y_1`$ $``$ $`{\displaystyle \frac{\xi ^r}{r}},y_2{\displaystyle \frac{\delta p}{p}}`$ (53)
$`y_3`$ $``$ $`{\displaystyle \frac{\delta \varphi }{c^2}},y_4{\displaystyle \frac{dy_3}{d\mathrm{ln}r}}`$
and obtain

$`{\displaystyle \frac{dy_1}{d\mathrm{ln}r}}`$ $`=`$ $`\left({\displaystyle \frac{\stackrel{~}{V}}{\gamma }}3\right)y_1{\displaystyle \frac{1}{\gamma }}y_2`$ (54a)
$`{\displaystyle \frac{dy_2}{d\mathrm{ln}r}}`$ $`=`$ $`\stackrel{~}{V}\left(1{\displaystyle \frac{1}{\gamma }}\right)y_2{\displaystyle \frac{\stackrel{~}{V}c_3}{Z_{}}}\left\{y_4\left({\displaystyle \frac{R}{r}}\right)^2{\displaystyle \frac{2}{3}}{\displaystyle \frac{R^2\mathrm{\Omega }^2}{c^2}}\right\}`$ (54b)
$`{\displaystyle \frac{dy_3}{d\mathrm{ln}r}}`$ $``$ $`y_4`$ (54c)
$`{\displaystyle \frac{dy_4}{d\mathrm{ln}r}}`$ $``$ $`{\displaystyle \frac{\stackrel{~}{U}Z_{}}{c_3\gamma }}\left({\displaystyle \frac{r}{R}}\right)^2y_2y_4.`$ (54d)
Boundary conditions at $`r=0`$ are obtained from the requirement that the solution be regular there. This yields

$$\left(4\frac{\mathrm{\Gamma }}{\gamma }\right)y_1+\frac{1}{\mathrm{\Gamma }}y_2=0\mathrm{and}y_4=0.$$
(55)
At $`r=R`$, $`\mathrm{\Delta }p`$ must vanish and the gravitational potential must match onto its solution in empty space ($`\delta \varphi 1/r`$), so there

$$y_2\stackrel{~}{V}y_1=0\mathrm{and}y_3+y_4=0,$$
(56)
while $`y_3`$ and $`y_4`$ are continuous at the crust-core interface.

## 4 Elastic Deformation of the Slowing Crust: Results

In this section, we describe the solutions to the elastic perturbation equations derived in § 3, and present our results for the value of $`b=\mathrm{\Delta }I_d/I_\mathrm{\Omega }`$ and compare it to the Baym-Pines estimate, Eq. (3).

In order to solve the perturbation equations, we need to specify the background NS model. We build NS models using several representative equations of state (EOS), summarized in Table 1. These models vary somewhat in their treatment of the pressure-density relation in the core and the crust. In particular, model FPS uses the equation of state described in Lorenz, Ravenhall & Pethick (1993); model WS uses EOS UV14+TNI of Wiringa et al. (1988) Wiringa\_88  in the core, matched to model FPS in the crust. Models AU and UU use equations of state AV14+UVII and UV14+UVII, respectively, matched to EOS of Negele & Vautherin (1973) in the crust. In practice, we use tabulations of pressure versus density provided with the code RNS<sup>1</sup><sup>1</sup>1N. Stergioulas, http://www.gravity.phys.uwm.edu/rns/. Finally, model n1poly uses pressure-density relation of an $`n=1`$ polytrope throughout the star. In all cases, we assume zero temperature, take $`M=1.4M_{}`$, and construct the NS model by solving the purely Newtonian equations of hydrostatic balance, $`dp/dr=\rho g`$, and continuity, $`dM_r/dr=4\pi r^2\rho `$. The resulting NS parameters (radius, $`I_0`$, $`I_C`$, and $`\mathrm{\Delta }I_\mathrm{\Omega }/\mathrm{\Omega }^2`$) are shown in Table 1.

In this paper, we generally assume that $`I_C`$ equals the moment of inertia of the crustal nuclei (i.e., excluding the dripped neutrons in the inner crust), which is roughly $`1\%`$ of $`I_0`$, since estimates of frictional coupling give coupling times between the crust and core well in excess of the spin period (see, e.g., Alpar & Sauls 1988). If the stellar magnetic field penetrates the core, magnetic stresses will also play a role in the coupling (Abney, Epstein & Olinto 1996; Mendell 1998), but the magnetic coupling timescale (estimated as the Alfvén crossing time for waves in the core) is again much longer than the rotation period for most neutron stars. However it is possible that some process does strongly couple the core to crust (in which case $`I_C`$ could be a substantial fraction of $`I_0`$), so we shall indicate how our main results scale with $`I_C/I_0`$.

In all models, we place the crust-core boundary at the fiducial density of $`\rho _b=1.5\times 10^{14}`$ g cm<sup>-3</sup> (Lorenz, Ravenhall & Pethick 1993). The crust-ocean boundary is placed at $`\rho =10^9`$ g cm<sup>-3</sup> for all models except for the $`n=1`$ polytrope one, where it is placed at $`10^{11}`$ g cm<sup>-3</sup> in order to ensure numerical stability of our integration of the elastic equations<sup>2</sup><sup>2</sup>2In the $`n=1`$ polytrope model, the density falls of very rapidly with radius in the outer layers–much faster than in models with a more realistic EOS. This very rapid variation over many orders of magnitude causes numerical problems. However, this truncation does not affect the final results, as the outer crust contains a negligible fraction of the mass..

We use the shear modulus $`\mu `$ as computed by Strohmayer et al. (1991), via Monte-Carlo simulations of both bcc crystals and quenched solids. Their results can be conveniently rewritten in terms of the pressure $`p_e`$ of degenerate relativistic electrons,

$$\frac{\mu }{p_e}=\frac{6\times 10^3}{1+0.595(173/\mathrm{\Gamma }_{\mathrm{Coul}})^2}\left(\frac{Z}{8}\right)^{2/3},$$
(57)
where $`\mathrm{\Gamma }_{\mathrm{Coul}}=Z^2e^2/akT`$ is the ratio of the Coulomb energy of the lattice to the thermal energy and $`Z`$ is the ionic charge. Throughout most of the crust (i.e., except near the top) $`\mathrm{\Gamma }_{\mathrm{Coul}}>10^3`$, so we for simplicity we ignore the slight dependence of the shear modulus on temperature, effectively setting $`1/\mathrm{\Gamma }_{\mathrm{Coul}}`$ equal to $`0`$. In order to compute the composition ($`A`$ and $`Z`$ of nuclei, as well as the free neutron fraction) we use the fits of Kaminker et al. (1999) to the formalism of Oyamatsu (1993). The resulting run of $`\mu `$ with density is shown in Figure 2. For the $`n=1`$ polytrope model we use a fiducial value of $`\mu /p=10^2`$.

With the background model specified as above, we solve the perturbation equations of § 3 using relaxation with adaptive mesh allocation (see, e.g., Press et al. 1992 Press\_92  for a description). For the case of $`l=2`$ perturbations in a purely fluid star, as well as for $`l=0`$ perturbations, the solution is straightforward. However, for $`l=2`$ perturbations in a star with a crust, there are internal boundary conditions at the crust-core interface (Eq. 42), which makes application of the standard relaxation method more complicated. In this case, to deal with this internal boundary, we transform the independent variable from $`x=\mathrm{ln}r`$ to $`t`$, such that for the fluid part of the star $`x=(x_cx_0)t+x_0`$, while in the crust $`x=x_1(x_1x_c)t`$, where $`x_0`$, $`x_1`$, and $`x_c`$ are the values near the center, at the crust-ocean boundary, and at the crust-core boundary, respectively. With this change of the independent variable, both the fluid and crust equations can be solved simultaneously on the interval $`t[0,1]`$: the internal boundary conditions are eliminated.

In Figure 3 we show the radial and transverse displacements induced when a relaxed, nonrotating star is spun up to $`\mathrm{\Omega }=c/R`$. (This is the solution of our perturbation equations with $`(R\mathrm{\Omega }/c)^2`$ set equal to $`1`$ in our source term Eq. 20.) Of course, $`\mathrm{\Omega }=c/R`$ is an unreasonably large spin value, but our solutions scale linearly with the source term, i.e., quadratically in $`\mathrm{\Omega }`$, so in all applications we simply scale the results down to low $`\mathrm{\Omega }`$. Recall that, when solving the elastic perturbation equations, we imagine taking a spherically symmetric star and spinning it up. In this case, the $`l=0`$ part of the centrifugal force acts radially outward, and hence the $`l=0`$ component of the radial displacement (lower panel of Figure 3) is positive throughout the crust. On the other hand, the $`l=2`$ part of the centrigual force squeezes the crust at the poles and pushes it out at the equator. Hence, $`z_1`$ (top panel of Figure 3), which can be thought of as the radial displacement at the poles (times $`R^1\sqrt{4\pi /5}`$), is negative. The transverse displacement, $`z_3`$ (middle panel Figure 3), changes sign in order to satisfy mass continuity.

Adding the $`l=0`$ and $`l=2`$ perturbations, we can rewrite the strain tensor (Eq. 25) as

$$\sigma _{ab}=(\sigma _0+\sigma _{rr}Y_{20})(\widehat{r}_a\widehat{r}_b\frac{1}{2}e_{ab})+\sigma _rf_{ab}+\sigma _\mathrm{\Lambda }(\mathrm{\Lambda }_{ab}+\frac{1}{2}e_{ab}).$$
(58)
Here the tensors $`e_{ab}`$, $`f_{ab}`$, and $`\mathrm{\Lambda }_{ab}`$ are defined by

$`e_{ab}`$ $``$ $`g_{ab}\widehat{r}_a\widehat{r}_b,`$ (59a)
$`f_{ab}`$ $``$ $`\beta ^1r(\widehat{r}_a_bY_{lm}+\widehat{r}_b_aY_{lm})),`$ (59b)
$`\mathrm{\Lambda }_{ab}`$ $``$ $`{\displaystyle \frac{r^2}{l(l+1)}}_a_bY_{lm}+f_{ab},`$ (59c)
the $`l=2`$ strain components $`\sigma _{rr}`$, $`\sigma _r`$, and $`\sigma _\mathrm{\Lambda }`$ are given by

$`\sigma _{rr}(r)`$ $`=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{dz_1}{d\mathrm{ln}r}}+{\displaystyle \frac{1}{3}}\beta ^2z_3`$ (60a)
$`\sigma _r(r)`$ $`=`$ $`{\displaystyle \frac{\beta ^2z_4}{2\mu }}`$ (60b)
$`\sigma _\mathrm{\Lambda }(r)`$ $`=`$ $`\beta ^2z_3`$ (60c)
and the $`l=0`$ piece is

$$\sigma _0=\frac{2}{3}\frac{dy_1}{d\mathrm{ln}r}$$
(61)
In Figure 4 we show the radial functions $`\sigma _{rr}`$ (top panel), $`\sigma _r`$ (second panel), $`\sigma _\mathrm{\Lambda }`$ (third panel), and $`\sigma _0`$ (bottom panel), for equations of state AU, WS, and n1poly. Again, these results are normalized to $`(R\mathrm{\Omega }/c)^2=1`$.

The crust breaks (or deforms plastically) when

$$\overline{\sigma }\left(\frac{1}{2}\sigma _{ab}\sigma ^{ab}\right)^{1/2}>\overline{\sigma }_{\mathrm{max}},$$
(62)
where $`\overline{\sigma }_{\mathrm{max}}`$ is the yield strain. (This is the von Mises criterion; see §6 of UCB for a discussion.) From Eq. (58), we get

$$2\overline{\sigma }^2=\frac{3}{2}\left(\sigma _0+\sigma _{rr}Y_{20}\right)^2+\sigma _r^2\frac{5}{64\pi }\left(\mathrm{sin}2\theta \right)^2+\sigma _\mathrm{\Lambda }\frac{5}{32\pi }\mathrm{sin}^4\theta ,$$
(63)
where $`Y_{20}=(5/16\pi )^{1/2}(3\mathrm{cos}^2\theta 1)`$. $`\overline{\sigma }`$ is independent of $`\varphi `$ by symmetry. In Figure 5 we plot contours of $`\overline{\sigma }`$ on a meridional plane (i.e., plane that slices through the crust of the star at a constant longitude $`\varphi `$), and in Figure 6 we show $`\overline{\tau }=2\mu \overline{\sigma }`$, again on a meridional plane. Clearly, the stresses are largest at the base of the crust, where $`\mu `$ is largest, while the strains are highest at the top of the crust (where they cost the least, energetically).

In Table 1 we show the value of $`b`$ computed as described in § 3 using the solutions to the perturbation problem and compare them to the estimates of Baym & Pines (1971) given by Eq. (3). Eq. (3) overestimates $`b`$ by a factor of $`1550`$. To emphasize this difference, and also to explore its dependence on $`\rho _b`$ (the assumed density at the crust-core boundary), we show in Figure 7 the ratio of our calculated $`b`$ to the Baym-Pines estimate as a function of $`\rho _b`$. In the two Appendices, we explore why the estimate Eq. (3) is inaccurate. Briefly, when a uniform elastic sphere spins down, its strain energy is strongly concentrated toward the center of the star; thus this case (on which Eq. 3 is based) is unsuitable for estimating the rigidity of a realistic NS, where all the strain energy is near the surface, in the thin crust. Moreover, we show that in a realistic NS crust, contributions to $`\mathrm{\Delta }I_d`$ from different stress components and different depths cancel each other to a much higher degree than one might expect from consideration of a uniform, incompressible crust.

## 5 $`\mathrm{\Delta }I_d`$ for non-relaxed crusts

So far in this paper we have calculated the residual oblateness $`\mathrm{\Delta }I_d`$ that a spinning NS with a relaxed crust retains if it is gently torqued down to zero spin frequency. But a real NS that is spinning down (or spinning up, due to accretion) is probably not completely relaxed. Indeed, we shall see below that the precession frequency of PSR1828-11 is inconsistent with the assumption of a relaxed crust (assuming shear stresses are indeed responsible for this NS’s $`\mathrm{\Delta }I_d`$). The crust, because it has finite rigidity, becomes strained as its spins changes. (As shown in Fig. 4 and §6, if no cracking or relaxation occurs, then the strain near the top of the crust grows to rather large values; e.g., for PSR 1828-11, the strain near the top would be $`10^3`$ in this case.) In some portions of the crust, the strain could be relaxed in localized episodes of crust cracking (starquakes) or plastic flow, while in other regions the strain may simply build up. Precisely because the NS crust may have a complex history of strain build-up and release, its current strain pattern is impossible to predict exactly.

Fortunately, however, many of our results can immediately be generalized to this case of a non-relaxed crust, by scaling them to the crust’s “reference spin” $`\nu _{s,ref}`$, defined as follows. We imagine adjusting the NS’s spin frequency $`\nu _s`$ while keepiing the crust’s preferred (zero-strain) shape fixed (i.e., without letting the crust crack or otherwise relax). While in general there will be no spin at which the crust is completely relaxed, if we plot the resulting quadrupole moment $`Q_{20}`$ as a function of $`\nu _s`$, there would be one spin value, $`\nu _{s,ref}`$, such that $`Q_{20}(\nu _{s,ref})`$ was identical to the $`Q_{20}`$ of a perfect-fluid NS with the same mass and spin, $`\nu _{s,ref}`$. (Here we are making the natural assumption that the crust’s preferred shape is roughly oblate, as opposed to prolate.) We define the crustal strain field at $`\nu _{s,ref}`$ to be $`\sigma _{ab}^{ref}`$. While $`\sigma _{ab}^{ref}`$ is in general non-zero, these strains make no net contribution to the NS’s quadrupole moment. In a rough way, one can think of $`\nu _{s,ref}`$ as the the spin at which the crust is most relaxed. As an example, consider a NS that is born with a relaxed crust at spin $`\nu _{s,init}`$ and never relaxes at all as it spins down to $`\nu _{s,final}`$. Then $`\nu _{s,ref}=\nu _{s,init}`$ and $`\sigma _{ab}^{ref}=0`$. If, on the other hand, the NS relaxes somewhat (but not completely) as it spins down, then $`\nu _{s,ref}`$ will lie somewhere between $`\nu _{s,init}`$ and $`\nu _{s,final}`$, and one would generally expect $`\sigma _{ab}^{ref}0`$.

The point of these definitions is the following. Imagine a NS with strained crust and spin $`\nu _s`$. To determine this NS’s precession frequency, we want to know what residual oblateness $`\mathrm{\Delta }I_d`$ it would have if it were spun down to zero frequency. We can imagine accomplishing the spindown in two steps. First, adjust the spin to $`\nu _{s,ref}`$. At this spin, the crust has some strain $`\sigma _{ab}^{ref}`$, but $`\sigma _{ab}^{ref}`$ has no net effect on the NS’s oblateness. Second, we spin the star down from $`\nu _{s,ref}`$ to zero angular velocity. We define $`\sigma _{ab}^{sd}`$ to be the extra strain that this spindown induces; i.e., $`\sigma _{ab}^{sd}\sigma _{ab}\sigma _{ab}^{ref}`$, where $`\sigma _{ab}`$ is the strain in the NS when it is nonrotating. Conceptually, we can consider both $`\sigma _{ab}^{ref}`$ and $`\sigma _{ab}^{sd}`$ as independent, small (i.e., linear) perturbations, and work to first order in both. To this order, $`\sigma _{ab}^{sd}`$ is just the strain induced when spinning down a relaxed NS from $`\nu _{s,ref}`$ to zero spin, which we solved for in §4. Also to this order, the oblateness $`\mathrm{\Delta }I_d`$ caused by $`\sigma _{ab}`$ is the sum of (i) the oblateness due to $`\sigma _{ab}^{sd}`$ and (ii) the oblateness due to $`\sigma _{ab}^{ref}`$, but (ii) vanishes, by definition. That is, $`\mathrm{\Delta }I_d`$ is just the oblateness induced by the spindown part of the strain field, Thus, generalizing Eq. (1), we can write

$$\mathrm{\Delta }I_d=b\mathrm{\Delta }I_\mathrm{\Omega }(\nu _{s,ref}),$$
(64)
where $`\mathrm{\Delta }I_\mathrm{\Omega }(\nu _{s,ref})0.3I_0\left(\nu _{s,ref}/\mathrm{kHz}\right)^2`$ is the oblateness of a fluid star with spin $`\nu _{s,ref}`$, and $`b`$ is the same coefficient (i.e., same numerical value) we calculated in § 3.

## 6 The precession frequency of PSR 1828-11

In §2 we showed that the measured precession period of PSR 1828-11 implies that for this NS,

$$\frac{\mathrm{\Delta }I_d}{I_C}|_{1828}=9.2\times 10^9\left(\frac{511\text{ d}}{P_p}\right).$$
(65)
For a NS with relaxed crust (i.e., relaxed except for the stresses induced by precessional motion itself), the predicted value is

$$\frac{\mathrm{\Delta }I_d}{I_C}=b\frac{\mathrm{\Delta }I_\mathrm{\Omega }}{I_0}\left(\frac{I_C}{I_0}\right)^1.$$
(66)
From the results of §4, we estimate (for our Newtonian NS models):

$`\mathrm{\Delta }I_\mathrm{\Omega }/I_0`$ $``$ $`0.3\left({\displaystyle \frac{\nu _s}{\mathrm{kHz}}}\right)^2\left({\displaystyle \frac{M}{1.4M_{}}}\right)^1\left({\displaystyle \frac{R}{12\text{ km}}}\right)^3,`$ (67)
$`I_C/I_0`$ $``$ $`0.01\left({\displaystyle \frac{M}{1.4M_{}}}\right)^2\left({\displaystyle \frac{R}{12\text{ km}}}\right)^4,`$ (68)
$`b`$ $``$ $`2\times 10^7\left({\displaystyle \frac{M}{1.4M_{}}}\right)^3\left({\displaystyle \frac{R}{12\text{ km}}}\right)^5.`$ (69)
The scaling with NS mass and radius are estimated by taking (i) $`\mathrm{\Delta }I_\mathrm{\Omega }/I_0MR^2/(M^2/R)=R^3/M`$ (i.e., to the ratio of the NS’s kinetic and potential energies), (ii) $`I_C/I_0M_C/M\rho _{bot}R^2\mathrm{\Delta }R/MR^2p_{bot}/(gM)R^4/M^2`$, and (iii) $`b\mu _{ave}V_C/(M^2/R)R^2\mathrm{\Delta }R/(M^2/R)R^2(p_{bot}/(g\rho _{bot}))/(M^2/R)R^5/M^3`$. Here $`M_C`$ is the mass of the crust, $`\mathrm{\Delta }R`$ is its thickness, $`V_C`$ its volume, $`\rho _{bot}`$ and $`p_{bot}`$ are the density and pressure at the bottom of the crust, $`\mu _{ave}`$ is a volume-weighted average of the crustal shear modulus, and $`gGM/R^2`$. Taking $`\nu _s=2.5`$Hz, we estimate for PSR 1828-11:

$$\frac{\mathrm{\Delta }I_d}{I_C}|_{\mathrm{relaxed}}3.8\times 10^{11}\left(\frac{I_C/I_0}{10^2}\right)\left(\frac{M}{1.4M_{}}\right)^2\left(\frac{R}{12\text{ km}}\right)^4.$$
(70)
Note that the measured value, Eq. (65), is $`250`$ times larger than our relaxed-crust estimate, for our fiducial NS parameters and precession period. And note that if some substantial fraction of the core is dynamically coupled to the crust (so $`I_C>>0.01I_0`$), then the discrepancy only gets worse. Even if we assume that PSR 1828’s precession period is really $`1000`$ days and take rather extreme values for the NS mass and radius, $`M=1.0M_{}`$ and $`R=20`$ km, we are still left with a factor $`8`$ discrepancy. We conclude that either our basic picture is wrong (i.e., either PSR 1828-11 is not actually precessing, or that some other mechanism besides crustal shear stress is responsible for its large $`\mathrm{\Delta }I_d`$), or that the crust is not relaxed. Here we adopt the latter explanation, and pursue its implications.

The 511-d precession period for PSR B1828-11 implies that the star’s reference spin is $`40`$ Hz. Subsequent spin down to its present spin of 2.5 Hz without significant structural relaxation (through quakes or plastic flow) would strain the crust and give the inferred $`\mathrm{\Delta }I_d/I_C=9.2\times 10^9`$. What are the current strain levels in the crust of 1828-11 in this scenario? Compared to $`40`$ Hz, its current spin of $`2.5`$ Hz is very slow. If one continued slowing PSR 1828-11 down to zero angular velocity, the spindown-induced crustal stresses would increase fractionally by only $`0.4\%`$. So (neglecting that $`0.4\%`$), the current strains are the sum of a) the reference strain $`\sigma _{ab}^{ref}`$ that this NS would still have if it were spun up to $`\nu _s=40`$Hz and b) the spindown strains $`\sigma _{ab}^{sd}`$ induced by spinning a relaxed NS down from $`\nu _s=40`$Hz to $`\nu _s=0`$. While $`\sigma _{ab}^{ref}`$ is practically unknowable, the spindown-induced strains $`\sigma _{ab}^{sd}`$ are those we found in §4.

We parametrize the strength of the spindown strain $`\sigma _{ab}^{sd}`$ by the scalar quantity $`\overline{\sigma }^{sd}`$, defined by

$$\overline{\sigma }^{sd}\left(\frac{1}{2}\sigma _{ab}^{sd}\sigma ^{sd,ab}\right)^{1/2},$$
(71)
and we define the $`\mu `$-weighted-average of the crustal spindown strain, $`\overline{\sigma }_{\mathrm{ave}}`$, by

$$\overline{\sigma }_{\mathrm{ave}}^{sd}\frac{\mu \overline{\sigma }𝑑V}{\mu 𝑑V}.$$
(72)
We have computed $`\overline{\sigma }_{\mathrm{ave}}^{sd}`$ for our five EOS; the results are in Table 1. They all give roughly the same result: $`\overline{\sigma }_{ave}^{sd}0.5`$ for $`R\mathrm{\Omega }/c=1`$. Scaling to the $`40`$Hz reference angular velocity of PSR1828-11 yields

$$\overline{\sigma }_{\mathrm{ave}}^{sd}=5\times 10^5\left(\frac{P_s}{0.4\mathrm{s}}\right)\left(\frac{P_p}{511\text{ d}}\right)^1\left(\frac{I_C/I_0}{0.01}\right)\left(\frac{b}{2\times 10^7}\right)^1.$$
(73)
Again, the total strain $`\sigma _{ab}`$ is the sum of $`\sigma _{ab}^{sd}`$ and $`\sigma _{ab}^{ref}`$, where $`\sigma _{ab}^{ref}`$ is determined by the detailed evolution of the crust, through structural adjustments, and so without knowing that history we can say nothing definitive about $`\sigma _{ab}^{ref}`$ (except that, by definition, these strains have no net effect on the NS’s quadrupole moment). However, since $`\sigma _{ab}^{ref}`$ will contain contributions from all harmonics (i.e., all $`Y_{lm}`$), while $`\sigma _{ab}^{sd}`$ is wholly $`l=2,m=0`$, it seems quite likely that the NS’s average total strain is larger than the average spindown strain. *I*f this is true, the average value of spindown portion of the strain places a lower limit on the crustal breaking strain: $`\overline{\sigma }_{max}>5\times 10^5(511\mathrm{d}/P_p)`$. In terrestrial solids, strains of order $`5\times 10^5`$ are rather easy to maintain, so it seems likely that the NS crust is strong enough to sustain the $`\mathrm{\Delta }I_d`$ implied by the precession frequency. If $`I_C/I_0`$ is of order unity, the implied average strain would be $`\overline{\sigma }_{ave}>5\times 10^3(511\mathrm{d}/P_p)`$, which is still plausibly below the crust’s breaking strain. We conclude that the neutron star crust is likely strong enough to deform PSR B1828-11 to the extent that it will precess with a period of $``$ 500 d.

The spindown strain field suggests that the evolutionary picture of crustal strain might be as follows. Suppose the star is born relaxed and rapidly-spinning. As the star spins down the strain becomes largest first in the upper crust (see Fig. 5). Cracking begins in the upper crust when the yield strain is reached there; however, because relatively little strain energy can be stored in this region of small shear modulus, these small quakes would not be expect to relax the strain that has developed deeper in the crust. When the yield strain is reached deeper, larger (i.e., more energetic) quakes occurs there. Eventually an equilibrium is reached in which the strain field is just sub-critical throughout the crust. If the Earth is any guide, then numerous small quakes should occur often, while large quakes occur rarely, but release most of the accumulated strain energy. These large events could excite the precession. We will explore this scenario further in a future publication.

In this picture, one expects that in the upper crust, $`\sigma _{ab}^{ref}`$ nearly cancels $`\sigma _{ab}^{sd}`$ (so that the sum is below the breaking strain), while near the bottom of the crust, where spindown strain in lower, the actual strain level $`\overline{\sigma }`$ is probably greater than $`\overline{\sigma }^{sd}`$. E.g., if the bottom of the crust has not significantly relaxed since the NS’s spin was $`60`$ Hz, then the total strain field near the bottom would presumably resemble $`\sigma _{ab}^{sd}`$ for spindown from $`\nu _{s,ref}=60`$Hz.

## 7 Other candidate precessing pulsars

So far we have focused our discussion on PSR 1828-11; we now apply our results to other pulsars showing modulations that may be due to precession. Here we discuss the Crab, Vela, PSRs B1642-03, B2217+47, B0959-54, and the possible remnant in SN 1987A. The strength of evidence for precession of these pulsars varies from more or less convincing to marginal; only one of these candidates (PSR B1642-03) shows strongly periodic and correlated changes in both pulse shape and phase similar to those observed in PSR 1828-11. In this section we merely summarize all published claims of precession, without attempting to assess their validity. Our discussion will be brief, since most these cases have recently been reviewed in depth by Jones & Anderson (2001). Note however that our conclusions are often different from Jones & Andersson (2001), because our more accurate value of $`b`$ is nearly two orders of magnitude smaller. Our results are summarized in Table 2.

Before proceeding, we note that two of the precession candidates that we discuss, the Crab and Vela pulsars, suffer glitches in spin rate. Glitches are thought to represent angular momentum transfer to the crust from the interior superfluid as the array of superfluid vortices, usually pinned to nuclei of the crust or magnetic flux tubes in the core, undergoes a sudden expansion (see, e.g., Anderson & Itoh 1975; Link & Epstein 1996; Ruderman, Zhu & Chen 1998). Even a small amount of vortex pinning is inconsistent with long-period precession (Shaham 1977). For the sake of this discussion, we ignore this objection. We refer the reader to Link & Cutler (2002) for a discussion of how precession could unpin vortices that are intially pinned to the crust.

### 7.1 PSR B1642-03

PSR B1642-03 has $`\nu _s=2.5`$ Hz (the same spin rate as PSR 1828-11), and shows a $`10^3`$ day periodicity in pulse shape, with modulation amplitude $`0.05`$ (Cordes 1993). This period is $`120`$ times shorter than the free-precession period, assuming a relaxed crust. If shear stresses are responsible for PSR B1642-03’s $`\mathrm{\Delta }I_d`$, then its reference angular velocity is $`27`$Hz.

### 7.2 PSR B0531+21 (the Crab pulsar)

The Crab pulsar has $`P_s=0.0331`$s; Lyne et al. (1988) observed a phase residual with period $`20`$ months, which Jones (1988) interpreted as evidence for precession. Thus $`P_s/P_p=6.4\times 10^{10}`$, while for our fiducial NS, we predict $`\mathrm{\Delta }I_d/I_C5.4\times 10^9`$ This factor $`8`$ discrepancy is not large considering the uncertainties in the NS EOS, crust thickness, and the mass and radius of the Crab; we conclude that a 20-month precession period is not in significant disagreement with our theoretical estimate. This conclusion is contrary to that of Jones & Anderson (2001), who concluded that the discrepancy was a factor $`700`$ instead of $`8`$.

### 7.3 PSR 2217+47

PSR 2217+47 has $`P_s=0.538`$ s, and was originally thought to have a single-component pulse profile. However, Suleymanova & Shitov (1994) reported the discovery of a weaker second component that varies with a period of $`68`$ years. Additionally, the braking index of 2217+47 has changed significantly between the original observation epoch (1974-1984) and the subsequent observations (Shabanova 1990). Suleymanova & Shitov (1994) interpret the change in the spindown rate as evidence for free precession with period at least comparable to the baseline of observations (i.e., $`P_p>20`$ years), and attribute the more rapid pulse shape variation to a patchy structure of the emission beam. Taking this interpretation at face value, the reference angular velocity of the crust is $`\nu _{s,ref}<12`$ Hz (compared to the current $`\nu _s=1.86`$ Hz) and $`\overline{\sigma }_{\mathrm{ave}}^{sd}<4\times 10^6`$.

### 7.4 PSR B0959-54

D’Alessandro & McCulloch (1997) show that the odd moments of the pulse profile of B0959-54 are negatively correlated with the variations in the timing residuals. The timing residuals appear quasi-periodic, and so they interpret this observation as evidence for free precession with $`P_p>2500`$ days (i.e., greater than or equal to the period of the timing residuals). This constraint on $`P_p`$ implies the crust’s reference angular velocity is $`\nu _{s,ref}<33`$ Hz (compared to spin frequency $`\nu _s=0.696`$ Hz), and $`\overline{\sigma }_{\mathrm{ave}}^{sd}<3\times 10^5`$.

### 7.5 PSR B0833-45 (the Vela pulsar)

The Vela pulsar has $`P_s=0.089`$s; Deshpande & McCulloch (1996) detected intensity variations with a $`165`$-d period, which they interpreted as possible evidence of precession. Thus $`P_s/P_p=6.2\times 10^9`$, while for our fiducial NS, we predict $`\mathrm{\Delta }I_d/I_C7\times 10^{10}`$. The factor $`9`$ discrepancy could be accounted for by Vela having a rather small mass and correspondingly large radius (e.g., $`M1.1M_{}`$ and $`R18`$ km), or by Vela having a crust whose reference angular velocity is $`\nu _{s,ref}32`$ Hz (i.e., $`3`$ times the current spin). The latter would entail an average crustal spindown strain of $`\overline{\sigma }_{ave}^{sd}3\times 10^5`$.

### 7.6 Remnant in SN1987A

Middleditch et al. (2001) have suggested that there is a precessing pulsar in the remnant of SN 1987A. They present evidence for the existence of a pulsar (now faded from view) with spin period $`P_s=2.14`$ms, and evidence for modulation (which they interpret as precession) on a timescale $`P_p10^3`$s, requiring $`\mathrm{\Delta }I_d/I_C2\times 10^6`$. Though failure to observe a pulsar in SN1987A since 1993 brings the existence of a pulsar into question, we apply our results assuming the pulsar exists. For our fiducial NS parameters, we would expect $`\mathrm{\Delta }I_d/I_C1.3\times 10^6`$. The factor $`1.5`$ difference is small compared to uncertainties in the NS EOS, mass, etc. We conclude (in contrast to Jones & Andersson 2001) that the $`10^3`$s precession timescale is quite reasonable, theoretically.

### 7.7 Discussion

Figure 8 summarizes the observational claims of precession presented in this section. In this figure, we plot the observed spin periods of neutron stars against their precession claimed periods (or lower limits on $`P_p`$ for PSRs B2217+47 and B0959-54). What limits can we place on the precession periods using the interpretation adopted in this paper, i.e., that the precession is due to a non-spherical shape of the crust? First, $`P_p=P_s(I_C/\mathrm{\Delta }I_d)=P_s(I_C/b\mathrm{\Delta }I_\mathrm{\Omega })\nu _{s,ref}^2`$ (for given $`\nu _s`$). Dashed lines in Figure 8 show the above relation for $`1/\nu _{s,ref}=1`$ ms and 20 ms, respectively. A pulsar lying below the 1 ms line would have a reference spin frequency faster than 1 kHz. Note that none of the pulsars with claims of precession require uncomfortably large $`\nu _{s,ref}`$, i.e., nowhere close to breakup spin of $`1.5`$ kHz for our fiducial mass and radius. (See Cook et al. 1994 for an exhaustive summary of breakup frequencies for various equations of state.) Second, if the crust of the neutron star is relaxed at its current spin frequency, then $`P_s=1/\nu _{s,ref}`$ and $`P_p\nu _{s,ref}^3`$. This relation is shown by the solid line in Figure 8. Since we presume that the NS had a relaxed crust in the past, and has since then spun down, theoretically feasible precession candidates must lie to the right of the solid line. All precession candidates except for the Crab pulsar satisfy this constraint. The Crab’s claimed precession frequency can also be explained in terms of its crustal rigidity, if one assumes that its mass and radius are somewhat different from our fiducial values (e.g., $`M=1.8M_{}`$ and $`R=8.0`$km).

## 8 Summary

A relaxed, self-gravitating object with a rigid crust has a portion of its spin-induced bulge that cannot follow changes in the direction of the solid’s rotation vector. Such an object can precess, at a frequency that is determined by the material properties of the crust and the strength of gravity, whose relative strengths are determined by the rigidity parameter $`b`$. In this paper showed that a NS’s rigidity parameter $`b`$ is a factor $`40`$ times smaller than the estimate of Baym & Pines (1971; see Eq. 3). Using this result, we showed that the precession frequency $`\nu _p`$ for PSR 1828-11 is $`250(511\mathrm{d}/P_p)`$ times faster than expected for a NS with a relaxed crust, indicating that its crust is significantly strained. The strain would naturally arise from the secular spin down of the star. Large quakes that partially relax the accumulated strain energy could also excite precession. Applying our results to PSR 1828-11 implies the average value of spindown strain in the crust is is $`\overline{\sigma }_{ave}^{sd}5\times 10^5`$ if the core and crust are effectively decoupled, and larger by a factor of $`100`$ if the crust and core are dynamically coupled over a spin period. It is not unreasonable to expect the NS crust to be able to sustain such modest strain values. Applying our model of a relaxed crust to other precession candidates, we found that only SN1987A and the Crab pulsar are consistent with the hypothesis of a relaxed crust. The other candidates can be reasonably explained as having strained crusts, though unfortunately this explanation has no predictive power; the one free parameter (the reference spin) is adjusted to fit the precession timescale.

We conclude here by mentioning one other application of our result for $`b`$. Cutler (2002) has shown that for rapidly rotating NS’s with relaxed crusts and a strong, interior toroidal magnetic field $`B_t`$, the prolate distortion of the star induced by $`B_t`$ dominates over the oblateness frozen into the crust for $`B_t>3.4\times 10^{12}\mathrm{G}(\nu _s/300\mathrm{H}\mathrm{z})^2`$. In this case, dissipation tends to drive the magnetic symmetry axis orthogonal to the spin direction, and the NS becomes a potent gravitational wave emitter. Had one used the Baym-Pines value for $`b`$, one would have arrived at a substantially larger requirement on the toroidal field (in order for it to dominate over crustal rigidity): $`B_t>1.4\times 10^{14}\mathrm{G}(\nu _s/300\mathrm{H}\mathrm{z})^2`$.

C.C.’s work was supported in part by NASA Grant NAG5-4093, B.L.’s work was supported in part by NSF Grant AST 00-98728, and G.U.’s work was supported by a Lee A. DuBridge postdoctoral fellowship at Caltech.

## Appendix A Approximate $`b`$ from integral of crustal stresses

Our solution for the NS’s rigidity parameter $`b`$ in § 3-4 involved calculating the exterior $`l=2`$ piece of the gravitational potential $`\delta \varphi `$ for two stars spun up to angular value $`\mathrm{\Omega }`$: i) a star that is completely fluid, and, ii) a star with relaxed, spherical crust that is otherwise identical to the first one. The relative difference in the exterior $`\delta \varphi `$ is $`b`$. Our result (Table 1) is that $`b2\times 10^7`$, so an accurate calculation of $`b`$ by the above subtraction method requires solving for $`\delta \varphi `$ to $`9`$ decimal places. Here, as an additional check on our work, we estimate $`b`$ in a way that does not require such high accuracy from our solutions.

UCB derived an identity relating a NS’s mass multipole moment $`Q_{22}`$ to an integral of stresses in the crust. The same calculation can be repeated for the $`Q_{20}`$ multipole moment. Defining $`Q_{20}`$ by

$$Q_{20}\delta \rho (r,\theta ,\varphi )Y_{20}(\theta ,\varphi )𝑑V$$
(A1)
we obtain

$`Q_{20}=`$ $`\left(1F\right)^1`$ $`{\displaystyle }{\displaystyle \frac{r^3}{g}}[{\displaystyle \frac{3}{2}}(4\stackrel{~}{U})t_{rr}+{\displaystyle \frac{1}{3}}(6\stackrel{~}{U})t_\mathrm{\Lambda }`$
$`+`$ $`\sqrt{{\displaystyle \frac{3}{2}}}(83\stackrel{~}{U}+{\displaystyle \frac{1}{3}}\stackrel{~}{U}^2{\displaystyle \frac{1}{3}}r{\displaystyle \frac{d\stackrel{~}{U}}{dr}})t_r]dr,`$
where $`g(r)=Gm(r)/r^2`$, $`\stackrel{~}{U}(r)4\pi \rho r^3/m(r)`$ (note $`\stackrel{~}{U}(r)<<1`$ in the crust), and the stress components of the stress tensor $`t_{ab}`$ are defined by

$`t_{ab}`$ $`=`$ $`t_{rr}(r)Y_{lm}(\widehat{r}_a\widehat{r}_b{\displaystyle \frac{1}{2}}e_{ab})`$
$`+`$ $`t_r(r)f_{ab}+t_\mathrm{\Lambda }(r)(\mathrm{\Lambda }_{ab}+{\displaystyle \frac{1}{2}}e_{ab}),`$
where the tensors $`e_{ab}`$, $`f_{ab}`$, and $`\mathrm{\Lambda }_{ab}`$ were defined in Eq. (59).

The factor $`(1F)^1`$ in Eq. (A2) requires some explanation. UCB derived Eq. (A2), without that factor, within the Cowling approximation, in which the self-gravity of the perturbation is neglected. So $`(1F)^1`$ is a correction factor that accounts the effect of the perturbation’s self-gravity. UCB derived an approximate expression for F (UCB, Eq. 72), leading to the estimate that $`F0.20.5`$, depending on the particular NS model.

We have calculated the right-hand side of Eq. (A) for the stresses obtained by solving our problem (ii), with the source term $`(R\mathrm{\Omega }/c)^2`$ set equal to $`1`$. The result is the residual $`Q_{20}`$ for a NS whose crust was relaxed at spin $`\mathrm{\Omega }=c/R`$, and which was then spun down to $`\nu _s=0`$. We find $`Q_{20}|_{\mathrm{residual}}=(1F)^14.24\times 10^{38}`$g-cm<sup>2</sup> for our fiducial $`1.4M_{}`$ NS, constructed with the AU EOS.

The rigidity parameter $`b`$ is the residual $`Q_{20}`$ of the non-rotating NS, divided by the $`Q_{20}`$ of the rotating model. The latter is given by

$$Q_{20}|_{\mathrm{rotating}}=\frac{5c^2R^3}{4\pi G}z_5|_{r=R},$$
(A4)
where $`z_5=\delta \varphi /c^2`$. For $`\mathrm{\Omega }`$ set equal to $`c/R`$, we find $`Q_{20}|_{\mathrm{rotating}}=2.8\times 10^{45}`$g-cm<sup>2</sup> for the same fiducial model. Thus integrating the stresses in the crust yields $`b=(1F)^1\times 1.5\times 10^7`$, or, for the $`F`$ in the range $`0.20.5`$, $`b1.93.0\times 10^7`$. This is in excellent agreement with the value $`b=2.47\times 10^7`$ in Table 1.

## Appendix B Analytic solution for $`b`$ for homogeneous sphere with thin crust

A classic problem solved by Lord Kelvin was the determination of $`b\mathrm{\Delta }I_d/\mathrm{\Delta }I_\mathrm{\Omega }`$ for the case of a constant-$`\rho `$, constant-$`\mu `$ sphere; the result can be written as $`b=\frac{57}{10}\mu V/|E_g|`$, where $`V`$ is the star’s volume and $`E_g`$ is its binding energy. (For the derivation, see the classic treatise by Love 1944.) In Fig. 9 we plot $`dE_{strain}/dr`$ for this case. We see that the strain energy is concentrated toward the center of the star, which already suggests that this model problem will not be a reliable guide for estimating the rigidity of a realistic NS (where all the shear stresses are near the surface.)

Here we solve the corresponding problem for a star with a thin crust. That is, we consider a constant-$`\rho `$ star that consists of two pieces: a fluid interior (where the shear modulus $`\mu `$ is zero), plus a thin crust where $`\mu `$ is constant. Our case is clearly much closer to a realistic NS, where the ratio (crust thickness)/(NS radius) is $`1/20`$. Most of the analysis we require has already been carried out by Franco, Link & Epstein (2001), who solved for the displacements and strain build-up that occur when such a star spins down. Franco, Link & Epstein (2001) work within the Cowling approximation–i.e., they neglect the gravitational perturbation $`\delta \varphi `$ induced by crustal distortions–so for convenience in this Appendix we neglect $`\delta \varphi `$ as well. As explained in Appendix A, the Cowling approximation underestimates $`b`$ by a factor $`5/3`$.

Franco, Link & Epstein (2001) did not restrict themselves to the case of a thin crust, but we do so, for convenience. Rather than repeat their derivation, we shall simply quote their formulae; e.g., we shall refer to their 39th numbered equation as (F39). Franco et al. define $`R^{}`$ and $`R`$ to be the radii at the bottom and top of the crust, respectively. We define $`\mathrm{\Delta }RRR^{}`$ and neglect terms that are of quadratic or higher order in $`\mathrm{\Delta }R`$. We just highlight the basic steps and refer the reader to Franco et al.(2001) for further details.

We will show that for the thin-crust case, and within the Cowling approximation,

$$b=\frac{12}{11}\frac{\mu V_c}{|E_g|}$$
(B1)
where $`V_c`$ is the volume of the crust and $`E_g`$ is still the binding energy of the entire star.

When the star’s spin-squared changes by $`\delta (\mathrm{\Omega }^2)`$, the crust undergoes a displacement $`\stackrel{}{u}(r,\theta )`$, whose radial piece can be shown to have the form $`u_r(r,\theta )=f(r)P_2(\theta )`$, where $`P_2`$ is the second Legendre polynomial, and where $`f(r=R)`$ has the form:

$$f(R)=\frac{5}{6}R\frac{\delta (\mathrm{\Omega }^2)}{v_k^2}\left(1b\right).$$
(B2)
Here $`v_k^2GM/R`$ and $`M`$ is the star’s mass. Eq. (F31) shows that the Eulerian change in the ($`l=2`$ piece of the) gravitational potential exterior to the star is $`f(R)`$, which makes it clear that the $`b`$ on the right-hand side of Eq. (B2) must also equal $`\mathrm{\Delta }I_d/\mathrm{\Delta }I_\mathrm{\Omega }`$. Expanding (F23) and (F24) to linear order in $`\mathrm{\Delta }R`$ and subsituting into (F19), we find

$`f(R)={\displaystyle \frac{1}{3}}AR^3{\displaystyle \frac{3}{16}}BR^2`$ (B3)
$`f^{}(R)={\displaystyle \frac{1}{4}}AR^2+{\displaystyle \frac{3}{4}}BR^3`$ (B4)
$`f^{}(R\mathrm{\Delta }R)=f^{}(R)\mathrm{\Delta }R[2AR{\displaystyle \frac{3}{4}}BR^4]`$ (B5)
where $`f^{}df/dr`$ and $`A`$ and $`B`$ are coefficients (to be solved for) appearing in the expansion of $`f(r)`$; see F19. Expanding F38 to linear order in $`\mathrm{\Delta }R`$ and comparing with Eq. (B5), we can solve for $`B`$ in terms of $`A`$:

$$BR^3=\frac{2}{3}AR^2+\frac{2}{5}AR\mathrm{\Delta }R.$$
(B6)
Plugging (B6) into Eqs. (B3)-(B4) yields

$`f(R)={\displaystyle \frac{11}{24}}AR^3{\displaystyle \frac{3}{40}}AR^2\mathrm{\Delta }R`$ (B7)
$`f^{}(R)={\displaystyle \frac{1}{6}}AR^2+{\displaystyle \frac{3}{10}}AR\mathrm{\Delta }R,`$ (B8)
and plugging Eqs. (B7)-(B8)into (F33) yields

$$f(R)=\frac{5}{6}R\frac{\delta (\mathrm{\Omega }^2)}{v_k^2}(1\frac{60}{11}\frac{\mathrm{\Delta }R}{R}\frac{c_t^2}{v_k^2}).$$
(B9)
where $`v_k^2\mu /\rho `$. The second term in parentheses on the righ-hand side of Eq. (B9) is $`b`$, and can be re-expressed as

$$b=\frac{12}{11}\mu V_c/|E_g|.$$
(B10)
Note that the coefficient $`12/11`$ in Eq. (B10) is smaller than the $`57/11`$ in the Baym-Pines estimate, Eq. (3), by $`209/405.23`$. So from our uniform, thin-crust model, we see that the Baym-Pines expression significantly overestimates the rigidity of a NS.

If we were to simply replace the $`57/11`$ in the Baym-Pines estimate by $`12/11`$, we would still overestimate the rigidity of a realistic NS by a factor $`8`$. To understand this factor, it is instructive to to evaluate $`Q_{20}|_{\mathrm{residual}}`$ for the uniform, thin-crust case using our integral expression Eq. (A2), and compare with the realistic case. For the uniform, thin-crust case, $`\stackrel{~}{U}=3`$, and using Eq. (76) in UCB we find that the correction factor $`(1F)^1`$ is exactly $`5/2`$. The results in Franco et al.(2001) and a few pages of algebra suffice to show that, to leading order in $`\mathrm{\Delta }R`$,

$$\sigma _{rr}=\frac{10}{33}\sqrt{\frac{4\pi }{5}}\frac{R^3\mathrm{\Omega }^2}{GM},\sigma _\mathrm{\Lambda }=\frac{10}{33}\sqrt{\frac{4\pi }{5}}\frac{R^3\mathrm{\Omega }^2}{GM},\sigma _r=0,$$
(B11)
for the case of star spun down from $`\mathrm{\Omega }`$ to zero angular velocity. Also, for a uniform star rotating at angular velocity $`\mathrm{\Omega }`$,

$$Q_{20}=\frac{5}{6}\sqrt{\frac{4\pi }{5}}\rho R^5\frac{R^3\mathrm{\Omega }^2}{GM}.$$
(B12)
Using these results and Eq. (A2), we find that the term $`\sigma _{rr}`$ in Eq. (A2) contributes $`(\frac{1}{2})\frac{12}{11}\mu V_c/|E_g|`$ to $`b`$, while the term $`\sigma _\mathrm{\Lambda }`$ contributes $`(\frac{3}{2})\frac{12}{11}\mu V_c/|E_g|`$.

Comparing these intermediate results to those for the realistic NS case (with compressible matter and a steep density gradient), we see that in the realistic case, there is much more cancellation within the integral for $`Q_{20}`$. As shown in Fig. 10, in the realistic case both $`\sigma _{rr}`$ and $`\sigma _\mathrm{\Lambda }`$ switch signs at different depths inside the crust (so that the contributions to $`b`$ from different layers tend to cancel each other, unlike in the uniform, thin-crust case). Also, the $`\sigma _{rr}`$ and $`\sigma _\mathrm{\Lambda }`$ contributions are clearly much closer in magnitude (though still of opposite sign) in the realistic case, and so cancel each other much more nearly.