text:
# Precession of Isolated Neutron Stars II: Magnetic Fields and Type II Superconductivity

## 1 Introduction

The convincing observation of free precession of PSR 1828-11 (Stairs, Lyne & Shemar 2000) poses challenges for theories of neutron stars. Shaham (1977; 1986) argued that vortex line pinning in the neutron star crust should prevent long term precession. Sedrakian, Wasserman & Cordes (1999) showed that precession is still prevented if vortex lines are not pinned perfectly but vortex drag is strong. They also showed that even if vortex drag is weak, precession is damped away. Link & Cutler (2002) estimated the strength of vortex line pinning forces, and argued that PSR 1828-11 may precess at large enough amplitude to unpin superfluid vortices in the crust. If vortex drag is small this would remove one impediment to free precession, although the free precession would still damp away eventually if it is not excited continuously.

Here, we consider an additional feature of radiopulsars like PSR 1828-11, namely, that they are strongly magnetized, with magnetic axes that are at an angle to their rotation axes. Based on earlier work on magnetic stars by Mestel and collaborators (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981; see also Spitzer 1958), we argue that precession may be required – there is no equilibrium corresponding to solid body rotation without precession for a rotating star with an oblique magnetic field. For a fluid star, though, we shall see that although the fluid must precess, the magnetic axis rotates uniformly. Although Mestel et al. (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981) showed that hydrostatic balance also requires fluid motions in addition to the precession which affect the stellar magnetic field, these are too slow to be important observationally.

Once we also take account of the solid crust of a neutron star in addition to its fluid interior, we show that if the stellar distortions due to the magnetic field are larger than the distortion of the crust, then steady state rotation is very unlikely (but not necessarily impossible). If the magnetic stresses inside PSR 1828-11 are simply due to the classical Maxwell stress tensor, evaluated with the inferred dipole magnetic field strength, then they are too weak to require precession at a period $`1`$ year. However, if the interior of the neutron star consists of a Type II superconductor, the effective stress tensor is larger than for a classical magnetic field, according to Jones (1975) and Easson & Pethick (1977); as was emphasized by these authors, and Cutler (2002), the magnetic distortion is correspondingly larger. For PSR 1828-11, we estimate that the distortion that would result in a core that contains a Type II superconductor can lead to precession at a period of order 1 year. A sufficiently strong toroidal magnetic field ($`B_t10^{14}`$ G) could also lead to precession at this period without Type II superconductivity (see e.g. Eq. \[2.4\] in Cutler  with $`B_c0`$).

Even if magnetic stresses are not strong enough to prevent steady rotation, magnetic distortion in an oblique rotator alters the physics of free precession qualitatively and quantitatively compared to what one would expect for precession due to axisymmetric crustal distortions. Even if the crust is axisymmetric, misalignment between the symmetry axis of the crust and the magnetic field make the effective stellar moment of inertia inherently triaxial. From a phenomenological viewpoint, one manifestation of this loss of axisymmetry is that although the angular velocity of the star is a periodic function of time in the frame rotating with the crust and magnetic field, it is not a sinusoidal function of time. This feature is consistent with observations of PSR 1828-11, which reveal behavior at several different harmonically related frequencies (Stairs, Lyne & Shemar 2000). The relative strengths of the harmonics depend on the degree of nonaxisymmetry, and, presuming an axisymmetric crust, on the distortions induced by the magnetic field. For the ordinary Maxwell stresses evaluated just with the inferred dipole field strength, the nonaxisymmetry would be small, but distortions resulting from a Type II superconductor, or from a normal core with a large toroidal field, could produce sufficient nonaxisymmetry to lead to comparable amplitudes for at least the first few harmonics of the fundamental precession period, as is observed.

Periodic, but not sinusoidal, precession would also arise if the neutron star crust were simply nonaxisymmetric even if magnetic stresses were negligible. Thus, the detection of harmonic behavior in PSR 1828-11 cannot, by itself, be taken to be evidence for amplified magnetic stresses in its core. However, standard results for triaxial precession, which are reproduced as a byproduct of the calculations we present, show that the oscillations at harmonics of the precession frequency are smaller in amplitude than the oscillation at the fundamental frequency by powers of the precession amplitude. We shall see that this is not the case in models where magnetic stresses predominate: the amplitudes of oscillations at the precessiom frequency and twice the precession frequency may be comparable even at small amplitude.

Moreover, it is well known that the minimum energy state for rotation of a nonaxisymmetric body is one in which the angular velocity and angular momentum are aligned with the principal axis of the body with largest eigenvalue of the moment of inertia tensor. There is no precession at all in this state. However, we shall see that even when magnetic stresses are not strong enough to require free precession, the minimum energy state for an oblique rotator does not correspond to alignment of the angular velocity with any principal axis of the effective moment of inertia tensor. Thus, the minimum energy state is one in which the star precesses.

We shall see that the timing residuals associated with the precession can account for observations of long term oscillations in PSR 1828-11, but that the explanation only works if the core of the neutron star has sufficiently strong magnetic stresses that magnetic distortions are 100-1000 times larger than crustal deformations. Thus, we argue that observations of free precession in PSR 1828-11 offer evidence that the pulsar is in the regime where magnetic stresses dominate, because either its core is a Type II superconductor, or, if it is a normal conductor, possesses a very large toroidal field.

We review the arguments given by Mestel and collaborators (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981) in §2. In §3, we consider the conditions under which precession of an oblique rotator is inevitable. In §4, we consider the modified precession problem for an axisymmetric crust and misaligned magnetic field. There we show that the minimum energy state at a given angular momentum is one in which the star precesses, provided that the magnetic field is not either along or perpendicular to the symmetry axis of the crust. There, we also review classical results for free precession of triaxial bodies; we shall find that there are three distinct cases of interest for a star where magnetic distortions are important. We present limiting results for the timing residuals expected in this model in §5, and obtain approximate results for the limit in which crustal distortions dominate in §5.1, and for the opposite limit in which magnetic distortions dominate in §5.2.

Here, our main purpose is to present arguments that an oblique rotator must precess. In §5.3, though, we present brief application of the model to observations of PSR 1828-11. There we argue that only models in which magnetic distortions dominate can account for all of the observed features of the long-term periodic timing residuals from this pulsar. We also suggest that the data favor prolate rather than oblate magnetic distortions.

Several appendices present cumbersome mathematical details needed to derive (and verify!) analytic results presented in the text.

## 2 Review of Structure of Fluid Equilibria with Oblique Magnetic Fields

Let us begin by reviewing the theory of rotating stars with oblique magnetic fields developed by Mestel and collaborators (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981; see also Spitzer 1958). Consider a (fluid) star with angular momentun $`𝑳=L\widehat{𝑳}`$. In the absence of a magnetic field, the angular velocity of the star would be $`𝛀_0=𝑳/I\mathrm{\Omega }_0\widehat{𝑳}`$. Let us work in a reference system rotating at a rate $`𝝎=\omega \widehat{𝑳}`$, where $`\omega `$ is unknown, but must be $`\mathrm{\Omega }_0`$ to lowest order in small quantities. Assume that the magnetic field is axisymmetric about an axis $`\widehat{𝒃}`$ that is fixed in this frame; we can verify that this works out correctly to lowest order later.

Solve for the density field of the star under the assumption that it is stationary; Mestel et al. (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981) discuss the sizes of time dependent correction terms, which are smaller than any we retain here. The result is a superposition of three components assuming slow rotation and weak magnetic fields:

$$\rho (𝒓)=\rho _0(r)+\rho _\omega (r,\widehat{𝒓}\mathbf{}\widehat{𝑳})+\rho _B(r,\widehat{𝒓}\mathbf{}\widehat{𝒃}).$$
(1)
Here, $`\rho _0(r)`$ is the spherical profile of the undistorted star, and $`\rho _\omega (r,\widehat{𝒓}\mathbf{}\widehat{𝑳})`$ and $`\rho _B(r,\widehat{𝒓}\mathbf{}\widehat{𝒃})`$ are the distortions due to rotation and magnetic fields, axisymmetric about $`\widehat{𝑳}`$ and $`\widehat{𝒃}`$, respectively. In absolute magnitude, the rotational distortion is of order $`\omega ^2R^3/GM`$ for a star of radius $`R`$ and mass $`M`$, and the magnetic distortion is of order $`BHR^4/GM^2`$, where $`H=B`$ for ordinary magnetic fields, but $`H`$ corresponds to the first critical field strength in a Type II superconductor ($`H_{c1}10^{15}`$G), as discussed by Jones (1975), Easson & Pethick (1977) and Cutler (2002).

The moment of inertia tensor corresponding to Eq. (1) is of the form

$$I_{ij}=d^3r\rho (𝒓)\left(r^2\delta _{ij}r_ir_j\right)=I_0\delta _{ij}+I_\omega \left(\delta _{ij}3\widehat{𝑳}_i\widehat{𝑳}_j\right)+I_B\left(\delta _{ij}3\widehat{𝒃}_i\widehat{𝒃}_j\right),$$
(2)
where ($`P_2(\mu )=\frac{1}{2}\left(3\mu ^21\right)`$)

$`I_0`$ $`=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle d^3r\rho (𝒓)r^2}`$
$`I_\omega `$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle d^3r\rho _\omega (r,\widehat{𝒓}\mathbf{}\widehat{𝑳})P_2(\widehat{𝒓}\mathbf{}\widehat{𝑳})}`$
$`I_B`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle d^3r\rho _B(r,\widehat{𝒓}\mathbf{}\widehat{𝒃})P_2(\widehat{𝒓}\mathbf{}\widehat{𝒃})}.`$ (3)
Let $`𝛀=𝑳/I+\delta 𝛀𝛀_0+\delta 𝛀`$; then the angular momentum of the star is

$$𝑳I𝛀_0=(I_0+I_\omega +I_B)𝛀3I_\omega \widehat{𝑳}\widehat{𝑳}\mathbf{}𝛀3I_B\widehat{𝒃}\widehat{𝒃}\mathbf{}𝛀(I_02I_\omega +I_B)𝛀_0+I\delta 𝛀3I_B\widehat{𝒃}\widehat{𝒃}\mathbf{}𝛀_0,$$
(4)
which implies that

$$\delta 𝛀\left(1\frac{I_0}{I}+\frac{2I_\omega }{I}\frac{I_B}{I}\right)𝛀_0+\frac{3I_B\widehat{𝒃}\widehat{𝒃}\mathbf{}𝛀_0}{I}.$$
(5)
We can decompose this into components along and perpedicular to $`\widehat{𝑳}`$ i.e. $`\delta 𝛀=\delta \mathrm{\Omega }_{}\widehat{𝑳}+\delta 𝛀_{}`$, where

$`\delta \mathrm{\Omega }_{}`$ $``$ $`\mathrm{\Omega }_0\left[1{\displaystyle \frac{I_0}{I}}+{\displaystyle \frac{2I_\omega }{I}}+{\displaystyle \frac{2I_B}{I}}P_2(\widehat{𝒃}\mathbf{}\widehat{𝑳})\right]`$
$`\delta 𝛀_{}`$ $``$ $`{\displaystyle \frac{3I_B\mathrm{\Omega }_0(\widehat{𝒃}\mathbf{}\widehat{𝑳})}{I}}\left(\widehat{𝒃}\widehat{𝑳}\widehat{𝒃}\mathbf{}\widehat{𝑳}\right)={\displaystyle \frac{3I_B\mathrm{\Omega }_0(\widehat{𝒃}\mathbf{}\widehat{𝑳})[\widehat{𝑳}\mathbf{\times }(\widehat{𝒃}\mathbf{\times }\widehat{𝑳})]}{I}};`$ (6)
$`\delta 𝛀_{}=0`$ only if $`\widehat{𝒃}`$ is either parallel to or perpendicular to $`\widehat{𝑳}`$. Since Euler’s equations imply time independent $`𝛀`$ only if $`𝑳`$ and $`𝛀`$ are parallel, in general a rotating star with an oblique magnetic field must precess.

To find the angular velocity of the reference system, require that the magnetic axis rotates with angular velocity $`𝛀\omega \widehat{𝑳}+3\widehat{𝒃}I_B\widehat{𝒃}\mathbf{}𝛀_0/I`$; thus

$$\omega =\mathrm{\Omega }_0\left(2\frac{I_0}{I}+\frac{2I_\omega }{I}\frac{I_B}{I}\right)$$
(7)
In this reference system, matter precesses about $`\widehat{𝒃}`$ with angular velocity

$$𝝎_p=\frac{3I_B\mathrm{\Omega }_0\widehat{𝒃}\widehat{𝒃}\mathbf{}\widehat{𝑳}}{I}.$$
(8)
However, even though the fluid precesses, the magnetic field does not, so we do not expect any observable effects to arise from the precession.

## 3 Conditions Under Which Steady Rotation is Impossible and Precession is Inevitable

We can look at the Mestel et al. (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981) results in a slightly different way, which is identical to the same order of approximation, but differs in approach slightly. Work in a reference system corotating with the matter. 111Mestel et al. (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981) show that maintaining hydrostatic balance also requires additional motions of the fluid. We shall discuss this feature in more detail below. The perfect conductivity condition demands that the magnetic axis is fixed in this reference system. Instead of Eq. (1), suppose that

$$\rho (𝒓)=\rho _0(r)+\rho _\mathrm{\Omega }(r,\widehat{𝒓}\mathbf{}\widehat{𝛀})+\rho _B(r,\widehat{𝒓}\mathbf{}\widehat{𝒃});$$
(9)
in this case, the moment of inertia tensor (formerly Eq. ) becomes

$$I_{ij}=d^3r\rho (𝒓)\left(r^2\delta _{ij}r_ir_j\right)=I_0\delta _{ij}+I_\mathrm{\Omega }\left(\delta _{ij}3\widehat{𝛀}_i\widehat{𝛀}_j\right)+I_B\left(\delta _{ij}3\widehat{𝒃}_i\widehat{𝒃}_j\right),$$
(10)
where the various factors are defined just as in Eq. (3). The angular momentum is therefore

$$𝑳=\left(I_02I_\mathrm{\Omega }+I_B\right)𝛀3I_B\widehat{𝒃}\widehat{𝒃}\mathbf{}𝛀.$$
(11)
The condition for a time independent $`𝛀`$ is $`𝑳𝛀`$. This condition is satisfied only if $`\widehat{𝒃}`$ is either parallel to or perpendicular to $`𝛀\widehat{𝑳}`$. The star precesses at an angular frequency $`𝝎_op`$ as before. However, the magnetic axis rotates at a uniform angular velocity, because, from Eq. (11)

$$𝛀\mathbf{\times }\widehat{𝒃}=\frac{𝑳\mathbf{\times }\widehat{𝒃}}{I_02I_\mathrm{\Omega }+I_B}𝛀_B\mathbf{\times }\widehat{𝒃},$$
(12)
where

$$𝛀_B=\frac{𝑳}{I_0+2I_\mathrm{\Omega }+I_B}𝛀_0\left(2\frac{I_0}{I}+\frac{2I_\mathrm{\Omega }}{I}\frac{I_B}{I}\right)𝝎,$$
(13)
which is independent of time.

The second approach makes contact with the Euler problem clearer. Angular momentum conservation is consistent with time independent rotation as long as the angular velocity vector is aligned with one of the principal axes of the moment of inertia tensor. For anisotropic density distribution, the moment of inertia tensor is of the form

$$I_{ij}=I_0\delta _{ij}+\mathrm{\Delta }I_{ij}$$

where $`\mathrm{Tr}(\mathrm{\Delta }I_{ij})=0`$. Even though $`\delta I_{ij}I_0`$ for small distortions, the principal axes of $`I_{ij}`$ are the principal axes of $`\mathrm{\Delta }I_{ij}`$. For the magnetic fluid, Eqs. (9) and (10) imply that

$$\mathrm{\Delta }I_{ij}=I_\mathrm{\Omega }\left(\delta _{ij}3\widehat{𝛀}_i\widehat{𝛀}_j\right)+I_B\left(\delta _{ij}3\widehat{𝒃}_i\widehat{𝒃}_j\right).$$

For there to be a principal axis of $`\mathrm{\Delta }_{ij}`$ along $`𝛀`$,

$$\mathrm{\Delta }_{ij}𝛀_j=\mathrm{\Lambda }𝛀_i.$$

where $`\mathrm{\Lambda }`$ is the associated eigenvalue. For this form of $`\mathrm{\Delta }_{ij}`$, the condition becomes

$$3I_B\widehat{𝒃}_i\widehat{𝒃}𝛀=(\mathrm{\Lambda }+2I_\mathrm{\Omega }I_B)𝛀_i,$$

which is only true in general provided that either $`𝛀`$ is along $`\widehat{𝒃}`$ or $`𝛀`$ is perpendicular to $`\widehat{𝒃}`$, neither of which will be the case generally.

Next, consider what happens when we consider a neutron star model that consists of a rigid solid and a fluid, with an oblique magnetic field. In that case, Eq. (10) is generalized to

$$I_{ij}=I_0\delta _{ij}+C_{ij}+I_\omega (\delta _{ij}3\widehat{𝛀}_i\widehat{𝛀}_j)+I_B\left(\delta _{ij}3\widehat{𝒃}_i\widehat{𝒃}_j\right),$$
(14)
where $`C_{ij}`$ is the trace-free part of the moment of inertia of the crust. 222More precisely, $`C_{ij}`$ is the portion of the moment of inertia tensor of the crust that is not aligned with the magnetic distortions. Any crustal distortions symmetric about $`\widehat{𝒃}`$ would just renormalize $`I_B`$. The angular momentum of the star is therefore

$$𝑳_i=\left(I_02I_\mathrm{\Omega }\right)𝛀_i+\left[C_{ij}+I_B\left(\delta _{ij}3\widehat{𝒃}_i\widehat{𝒃}_j\right)\right]\mathrm{\Omega }_j.$$
(15)
A steady state is only possible if the star rotates along a principal axis of

$$\mathrm{\Delta }I_{ij}C_{ij}+I_B\left(\delta _{ij}3\widehat{𝒃}_i\widehat{𝒃}_j\right).$$

But for this to be true, we must have

$$C_{ij}𝛀_j3\widehat{𝒃}_iI_B(\widehat{𝒃}\mathbf{}𝛀)=(\mathrm{\Lambda }I_B)𝛀_i,$$
(16)
where $`\mathrm{\Lambda }`$ is the associated eigenvalue. The components of this equation perpendicular to $`𝛀`$ are

$$C_{ij}\widehat{𝛀}_j\widehat{𝛀}_i\widehat{𝛀}_kC_{kj}\widehat{𝛀}_j+3I_B(\widehat{𝒃}\mathbf{}\widehat{𝛀})\left[\widehat{𝛀}_i(\widehat{𝒃}\mathbf{}\widehat{𝛀})\widehat{𝒃}_i\right]=0.$$
(17)
Eq. (17) must have a solution (with $`\widehat{𝛀}=\widehat{𝑳}`$) in order for there to be no precession.

To understand the significance of Eq. (17), first review the argument against precession when magnetic fields are ignored. The moment of inertia tensor of the solid is not known observationally for any neutron star, so we can specify $`C_{ij}`$ freely. Since it is a trace-free symmetric tensor, $`C_{ij}`$ has five independent components. We only need to choose two of these (the orientations of one of the three orthogonal principal axes) to get a steady state. We can make these choices without considering the magnitude of the distortion of the solid, since only the directions of the eigenvectors matters for finding a steady state. Moreover, from a physical viewpoint, we expect that at a fixed $`𝑳`$ the energy associated with rotation is minimized if $`𝑳`$ (and hence $`𝛀`$) is along the principal axis with the largest eigenvalue; dissipative processes tend to drive the neutron star toward this state.

Next, let us consider what happens when we restore the magnetic field. Expand

$$C_{ij}=\underset{\mu }{}C_\mu \widehat{𝝀}_i^\mu \widehat{𝝀}_j^\mu ,$$
(18)
where $`C_\mu `$ are the eigenvalues and $`\widehat{𝝀}^\mu `$ the eigenvectors of $`C_{ij}`$, and rewrite Eq. (17) as

$$\underset{\mu }{}C_\mu \left[\widehat{𝝀}^\mu \mathbf{}\widehat{𝛀}\widehat{𝝀}^\mu \widehat{𝛀}(\widehat{𝝀}^\mu \mathbf{}\widehat{𝛀})^2\right]+3I_B(\widehat{𝒃}\mathbf{}\widehat{𝛀})\left[\widehat{𝛀}(\widehat{𝒃}\mathbf{}\widehat{𝛀})\widehat{𝒃}\right]=0.$$
(19)
We know that if $`I_B=0`$, a steady state is possible irrespective of the magnitudes $`C_\mu `$, and that if all $`C_\mu =0`$ there is no steady state. When $`I_B`$ is nonzero, we can adjust the $`\widehat{𝝀}^\mu `$ so that Eq. (19) is satisfied, just as for $`I_B=0`$, but the adjustment depends on the magnitudes $`C_\mu `$ relative to $`I_B`$. In general, we would expect there to be no solution if $`|C_\mu |I_B`$, in which case the neutron star must precess. In § 4 we shall also see that even when $`I_B`$ is not large compared to $`C`$, the rotational energy is minimized, to first order in distortions, when the angular velocity is not precisely along a principal axis of the effective moment of inertia tensor, and therefore not precisely along $`\widehat{𝑳}`$.

To appreciate this result better quantitatively, let us focus on the simplest special case in which the crust is axisymmetric. If $`\widehat{𝒌}`$ is the symmetry axis of the crust, then we may write

$$C_{ij}=C\left(3\widehat{𝒌}_i\widehat{𝒌}_j\delta _{ij}\right)$$
(20)
so that the largest eigenvalue of the moment of inertia of the crust is along $`\widehat{𝒌}`$. Then the condition for a steady state, Eq. (19), may be written as

$$3C(\widehat{𝒌}\mathbf{}\widehat{𝛀})\left[\widehat{𝒌}\widehat{𝛀}(\widehat{𝒌}\mathbf{}\widehat{𝛀})\right]+3I_B(\widehat{𝒃}\mathbf{}\widehat{𝛀})\left[\widehat{𝛀}(\widehat{𝒃}\mathbf{}\widehat{𝛀})\widehat{𝒃}\right]=0.$$
(21)
If Eq. (21) has a solution, it must have coplanar $`\widehat{𝒌},\widehat{𝒃}\widehat{𝛀}`$. Thus, let us adopt $`\widehat{𝒃}=\widehat{𝒆}_\mathrm{𝟑}`$ and

$`\widehat{𝛀}`$ $`=`$ $`\widehat{𝒆}_\mathrm{𝟑}\mathrm{cos}\chi +\widehat{𝒆}_\mathrm{𝟏}\mathrm{sin}\chi `$
$`\widehat{𝒌}`$ $`=`$ $`\widehat{𝒆}_\mathrm{𝟑}\mathrm{cos}\theta +\widehat{𝒆}_\mathrm{𝟏}\mathrm{sin}\theta :`$ (22)
substitute into Eq. (21) to find

$$\left(\widehat{𝒆}_\mathrm{𝟑}\mathrm{sin}\chi \widehat{𝒆}_\mathrm{𝟏}\mathrm{cos}\chi \right)\left[C\mathrm{cos}(\chi \theta )\mathrm{sin}(\chi \theta )I_B\mathrm{sin}\chi \mathrm{cos}\chi \right]=0.$$
(23)
Therefore, for a steady state, we must require

$$\mathrm{sin}\left[2(\chi \theta )\right]=\frac{I_B\mathrm{sin}2\chi }{C},$$
(24)
which is only possible if $`|I_B\mathrm{sin}2\chi |C`$. If $`I_BC`$, then steady state rotation is rather unlikely, with a probability that decreases with increasing $`I_B/C`$. When $`I_B\mathrm{sin}2\chi /C1`$, Eq. (24) has the solutions

$`\mathrm{cos}2\theta `$ $`=`$ $`{\displaystyle \frac{I_B\mathrm{sin}^22\chi }{C}}\pm \mathrm{cos}2\chi \sqrt{1{\displaystyle \frac{I_B^2\mathrm{sin}^22\chi }{C^2}}}`$
$`\mathrm{sin}2\theta `$ $`=`$ $`\mathrm{sin}2\chi \left(\pm \sqrt{1{\displaystyle \frac{I_B^2\mathrm{sin}^22\chi }{C^2}}}{\displaystyle \frac{I_B\mathrm{cos}2\chi }{C}}\right),`$ (25)
as can be verified by direct substitution; the choice of signs depends on specific parameter values.

In order of magnitude, $`3I_B/I\beta BHR^4/GM^2`$, where $`\beta 1`$ is a structure constant, so that

$$\omega _B\frac{3I_B\mathrm{\Omega }_0\mathrm{cos}\chi }{I}=1.9\times 10^9\mathrm{\Omega }_0\beta \mathrm{cos}\chi (BH)_{27}R_6^4M_{1.4}^2,$$
(26)
for a neutron star mass and radius $`M=1.4M_{1.4}M_{}`$ and $`R=10R_6\mathrm{km}`$, a magnetic axis inclination angle $`\chi =\mathrm{cos}^1(\widehat{𝒃}\mathbf{}\widehat{𝑳})`$, and $`BH=10^{27}(BH)_{27}\mathrm{G}`$. ($`BH10^{27}\mathrm{G}`$ for $`B10^{12}\mathrm{G}`$ and $`H10^{15}\mathrm{G}H_{c1}`$ – see Easson & Pethick 1977.) The period associated with $`\omega _B`$ is $`P_B=2\pi /\omega _B16.6P_0(\mathrm{s})(\beta \mathrm{cos}\chi )^1M_{1.4}^2R_6^4(BH)_{27}^1\mathrm{y}`$ for a neutron star rotation period $`P_0=1P_0(\mathrm{s})\mathrm{s}`$. If magnetic effects are strong enough to require precession, then we should expect a precession period of order $`P_B`$. We consider the precession period more completely in the next section.

For PSR 1828-11, the spin period is $`P_0=0.405\mathrm{s}`$ and the dipole field strength deduced from the observed spindown is $`B5\times 10^{12}\mathrm{G}`$ (Stairs, Lyne & Shemar 2000); these values would imply $`P_B1.35(\beta \mathrm{cos}\chi )^1[5/(BH)_{27}]\mathrm{y}`$, or about $`492(\beta \mathrm{cos}\chi )^1(BH_{27}/5)\mathrm{d}`$, similar to the observed period for $`\beta \mathrm{cos}\chi 1`$. Note that the precession period would be far longer for $`BHB^2`$, by a factor of about 200. Thus, if magnetic effects are the reason for the observed “precession” then ether the neutron star interior must be a Type II superconductor (Jones 1975, Easson & Pethick 1977, Cutler 2002), or there must be a substantial toroidal field in the core if it is not a Type II superconductor (e.g. Eq. \[2.4\] in Cutler 2002 with $`B_c0`$).

## 4 Modified Euler Problem

### 4.1 Basic Equations, Principal Axes and Eigenvalues

Next, let us consider the modified Euler problem. To allow an analytic treatment, continue to assume that the crust is axisymmetric. Although this is a simplification, the effective moment of inertia of the star will still be triaxial because of the misaligned distortion introduced by the magnetic field, except for special orientations.

The Euler equation in the rotating frame of reference is

$$I_{ij}^{eff}\frac{d^{}\mathrm{\Omega }_j}{dt}+3ϵ_{ijl}\mathrm{\Omega }_j\left(C\widehat{𝒌}_l\widehat{𝒌}_mI_B\widehat{𝒃}_l\widehat{𝒃}_m\right)\mathrm{\Omega }_m=0,$$
(27)
where the effective moment of inertia tensor is

$$I_{ij}^{eff}=(I_02I_\mathrm{\Omega }+I_BC)\delta _{ij}+3\left(C\widehat{𝒌}_i\widehat{𝒌}_jI_B\widehat{𝒃}_i\widehat{𝒃}_j\right).$$
(28)
and $`d^{}/dt`$ is the time derivative in the rotating frame. The two vectors, $`\widehat{𝒃}`$ and $`\widehat{𝒌}`$ define a plane, and one of the eigenvectors of $`I_{ij}^{eff}`$ is along the unit vector perpendicular to that plane, $`\widehat{𝒆}_\mathrm{𝟐}`$, and has an eigenvalue $`I_2=I_02I_\mathrm{\Omega }+I_BC`$. The other two eigenvectors lie in the $`\widehat{𝒃}\widehat{𝒌}`$ plane. As in the previos section, take $`\widehat{𝒆}_\mathrm{𝟑}=\widehat{𝒃}`$ and define $`\widehat{𝒌}`$ by Eq. (22). Then the other two eigenvectors, $`\widehat{𝒆}_\mathbf{\pm }`$ are 333The following results imply that
$$\widehat{𝒆}_\mathbf{+}\mathbf{\times }\widehat{𝒆}_{\mathbf{}}=\widehat{𝒆}_\mathrm{𝟐}$$

$$\widehat{𝒆}_\mathrm{𝟐}\mathbf{\times }\widehat{𝒆}_\mathbf{+}=\widehat{𝒆}_{\mathbf{}}$$

$$\widehat{𝒆}_\mathrm{𝟐}\mathbf{\times }\widehat{𝒆}_{\mathbf{}}=\widehat{𝒆}_\mathbf{+}$$
so that these axes define a right handed coordinate system.

$$\widehat{𝒆}_\mathbf{\pm }=\widehat{𝒆}_\mathrm{𝟑}\mathrm{cos}\sigma _\pm +\widehat{𝒆}_\mathrm{𝟏}\mathrm{sin}\sigma _\pm ,$$
(29)
where

$`\mathrm{sin}\sigma _\pm `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(1\pm {\displaystyle \frac{\eta }{\sqrt{1+\eta ^2}}}\right)^{1/2}`$
$`\mathrm{cos}\sigma _\pm `$ $`=`$ $`\pm {\displaystyle \frac{1}{\sqrt{2}}}\left(1{\displaystyle \frac{\eta }{\sqrt{1+\eta ^2}}}\right)^{1/2}`$
$`\eta `$ $`=`$ $`{\displaystyle \frac{I_BC\mathrm{cos}2\theta }{C\mathrm{sin}2\theta }}.`$ (30)
The eigenvalues associated with these eigenvectors are $`I_\pm =I_2+\delta I_\pm `$, where

$$\delta I_\pm =\frac{3}{2}\left[CI_B\pm \sqrt{C^22CI_B\mathrm{cos}2\theta +I_B^2}\right].$$
(31)
In Appendix A we show that for $`I_B>0`$, $`I_+>I_2>I_{}`$, and for $`I_B<0`$, $`I_+>I_{}>I_2`$.

Approximate results are derived for $`|I_B|C`$ in Appendix B and for $`|I_B|C`$ in Appendix C. When $`|I_B|C`$, the eigenvectors are nearly aligned with the principal axes of $`C_{ij}`$, and the eigenvalues are nearly the eigenvalues of $`C_{ij}`$, apart from small corrections $`|I_B/C|`$. For $`|I_B|C`$, the eigenvectors are nearly along and perpendicular to $`\widehat{𝒃}`$ and the eigenvalues are almost determined by the magnetic distortions alone, apart from corrections $`|C/I_B|`$. We shall argue below that the case $`I_B<0`$, $`|I_B/C|1`$ may be especially relevant to observations of PSR 1828-11.

### 4.2 Solution of the Euler Equations

#### 4.2.1 Basic Equations and Minimum Energy State

The modified Euler equations are simply what one finds in general for a triaxial system,

$`I_+{\displaystyle \frac{d^{}\mathrm{\Omega }_+}{dt}}\delta I_{}\mathrm{\Omega }_2\mathrm{\Omega }_{}`$ $`=`$ $`0`$
$`I_{}{\displaystyle \frac{d^{}\mathrm{\Omega }_{}}{dt}}+\delta I_+\mathrm{\Omega }_2\mathrm{\Omega }_+`$ $`=`$ $`0`$
$`I_2{\displaystyle \frac{d^{}\mathrm{\Omega }_2}{dt}}(\delta I_+\delta I_{})\mathrm{\Omega }_+\mathrm{\Omega }_{}`$ $`=`$ $`0.`$ (32)
As is well-known, Eqs. (32) conserve both the magnitude of $`𝑳`$ and the rotational energy,

$$E_{rot}=\frac{1}{2}\left(I_+\mathrm{\Omega }_+^2+I_{}\mathrm{\Omega }_{}^2+I_2\mathrm{\Omega }_2^2\right).$$
(33)
Note that the magnitude of the angular velocity is not conserved, so that in actuality $`I_\mathrm{\Omega }`$ is not independent of time. However we will assume that the variation is slow enough, and of small enough amplitude, that its effect is only higher order than any others we consider here.

The condition for steady rotation, Eq. (24), can be rederived from Eqs. (30) under the assumption that $`\widehat{𝛀}`$ is along either $`\widehat{𝒆}_\mathbf{+}`$ or $`\widehat{𝒆}_{\mathbf{}}`$, with $`\chi <\pi /2`$ or $`\chi >\pi /2`$, respectively. We are interested in what happens when precession is required, so let us assume that Eq. (24) cannot be satisfied for any choice of $`\theta `$ given $`\chi `$, $`I_B`$ and $`C`$. Nevertheless, we expect that the system seeks a minimum energy state, even if that state is not steady. To be specific, let $`\chi `$ be the angle between $`\widehat{𝒃}`$ and $`\widehat{𝛀}`$ when they are coplanar with $`\widehat{𝒌}`$ (and therefore with $`\widehat{𝒆}_\mathbf{\pm }`$). Let us determine $`\theta `$ from the requirement that $`E`$ is minimum for a given $`L^2`$. Thus, let us consider the quantity $`2E/L^2`$ at the epoch when $`\widehat{𝛀}`$, $`\widehat{𝒌}`$ and $`\widehat{𝒃}`$ are coplanar. If we define

$`\overline{I}`$ $``$ $`{\displaystyle \frac{1}{2}}(I_++I_{})=I_2+{\displaystyle \frac{3(CI_B)}{2}}`$
$`\delta I`$ $``$ $`{\displaystyle \frac{1}{2}}(I_+I_{})`$
$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{\delta I}{\overline{I}}}={\displaystyle \frac{3\sqrt{C^22CI_B\mathrm{cos}2\theta +I_B^2}}{2I_2+3(CI_B)}},`$ (34)
then we find

$$\frac{2\overline{I}E}{L^2}=\frac{1+\mathrm{\Delta }(2\widehat{𝛀}_+^21)}{1+\mathrm{\Delta }^2+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)}1\mathrm{\Delta }(2\widehat{𝛀}_+^21),$$
(35)
where the approximation holds for $`\mathrm{\Delta }1`$, and

$`\widehat{𝛀}_+\widehat{𝛀}\mathbf{}\widehat{𝒆}_\mathbf{+}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\{\mathrm{cos}\chi [1{\displaystyle \frac{(I_BC\mathrm{cos}2\theta )}{\sqrt{C^22CI_B\mathrm{cos}2\theta +I_B^2}}}]^{1/2}`$
$`+\mathrm{sin}\chi [1+{\displaystyle \frac{(I_BC\mathrm{cos}2\theta )}{\sqrt{C^22CI_B\mathrm{cos}2\theta +I_B^2}}}]^{1/2}\}`$
$`\widehat{𝛀}_{}\widehat{𝛀}\mathbf{}\widehat{𝒆}_{\mathbf{}}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\{\mathrm{sin}\chi [1{\displaystyle \frac{(I_BC\mathrm{cos}2\theta )}{\sqrt{I_B^22CI_B\mathrm{cos}2\theta +C^2}}}]^{1/2}`$ (36)
$`\mathrm{cos}\chi [1+{\displaystyle \frac{(I_BC\mathrm{cos}2\theta )}{\sqrt{I_B^22CI_B\mathrm{cos}2\theta +C^2}}}]^{1/2}\}.`$
Assuming that $`\theta \pi /2`$ we find

$$\mathrm{\Delta }(2\widehat{𝛀}_+^21)=\frac{3\left\{C\mathrm{cos}[2(\chi \theta )]I_B\mathrm{cos}2\chi \right\}}{2I_2+3(CI_B)}.$$
(37)
To first order in $`\mathrm{\Delta }`$, the energy is minimized when $`\mathrm{\Delta }(2\widehat{𝛀}^21)`$ is maximized, which happens when $`\theta =\chi `$.

Note that in the minimum energy configuration, $`\widehat{𝛀}_{}`$ need not be very small, although it vanishes in the minimum energy state for $`I_B0`$. For small values of $`I_B/C`$,

$$\widehat{𝛀}_{}\mathrm{sin}(\chi \theta )\frac{I_B}{2C}\mathrm{sin}2\theta \mathrm{cos}(\chi \theta )\frac{I_B\mathrm{sin}2\theta }{2C},$$
(38)
where the last result assumes $`\theta \chi `$. Thus, for small $`I_B`$, we find small but nonzero $`\widehat{𝛀}_{}`$ in the minimum energy state. For $`|I_B/C\mathrm{cos}2\theta |1`$, $`\widehat{𝛀}_{}I_B\mathrm{sin}2\chi /2C`$. On the other hand, for $`|I_B/C\mathrm{cos}2\theta |1`$, $`\widehat{𝛀}_{}\mathrm{sin}\chi `$ or $`\mathrm{cos}\chi `$, depending on the sign of $`I_B`$.

When $`I_B0`$, $`\theta =\chi `$ minimizes $`2E\overline{I}/L^2`$ exactly, and we find, as usual, $`\widehat{𝛀}_+=1`$, corresponding to an angular velocity aligned with the principal axis with the largest moment of inertia. When $`I_B0`$, $`\theta =\chi `$ does not correspond to exact alignment of the angular velocity and the principal axis with the largest moment of inertia. This is because the eigenvalues of $`I_{ij}^{eff}`$ depend on $`\theta `$, so $`\mathrm{\Delta }`$ depends on $`\theta `$. If $`\mathrm{\Delta }`$ were independent of $`\theta `$, then the energy would be minimized for $`\widehat{𝛀}_+=1`$, the largest possible value of $`\widehat{𝛀}_+`$. The value of

$$\widehat{𝑳}\mathbf{}\widehat{𝛀}=\frac{1+\mathrm{\Delta }(2\widehat{𝛀}_+^21)}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}12\mathrm{\Delta }^2\widehat{𝛀}_+^2\widehat{𝛀}_{}^2$$
(39)
is only exactly one when either $`\widehat{𝛀}_+=1`$ or $`\widehat{𝛀}_{}=1`$; the angle between $`\widehat{𝑳}`$ and $`\widehat{𝛀}`$ is $`\pm 2\mathrm{\Delta }\widehat{𝛀}_+\widehat{𝛀}_{}`$ more generally.

The condition that $`\theta =\chi `$ means that the angle between the angular velocity vector and the magnetic field is the same as the angle between the symmetry axis of the crust and the field. Thus, it corresponds to angular velocity that is parallel to the symmetry axis of the crust at the epoch when all of the vectors lie in a plane. Since the angle between the angular velocity vector and the magnetic axis is fixed in this case, the portion of the rotational energy due to the rotating magnetic distortion is fixed. Minimizing the energy then amounts to minimizing the portion of the energy associated with rotation of the crustal distortions. This is achieved if the angular velocity is along the symmetry axis of the crust. The resulting minimum is only a local minimum, at a given value of $`\chi `$. The global minimum is achieved for $`\chi =0`$ or $`\pi /2`$. We assume that the star can evolve on a slow disspative timescale toward the local minimum, and on a longer timescale, probably the spindown timescale (e.g. Goldreich 1970) toward the global minimum.

The fact that $`\mathrm{\Delta }=\mathrm{\Delta }(\mathrm{cos}2\theta )`$ when $`I_B0`$ implies that even when alignment of $`𝛀`$ and $`𝑳`$ is possible, the minimum energy state is not steady rotation, but rather precession. When $`I_BC`$, the minimum energy state corresponds to an angle of approximately $`|\widehat{𝛀}_{}||I_B\mathrm{sin}2\chi |/2C`$ between $`\widehat{𝛀}`$ and $`\widehat{𝒆}_\mathbf{+}`$, and an angle $`3|I_B\mathrm{sin}2\chi |/2\overline{I}`$ between $`\widehat{𝛀}`$ and $`\widehat{𝑳}`$.

#### 4.2.2 Nonlinear Solution: $`I_B>0`$

The complete, nonlinear solution to Eqs. (32) is given in Landau & Lifshitz, §37, in terms of elliptic functions. We will have to consider separately the two cases $`I_B>0`$ and $`I_B<0`$. Here, we consider $`I_B>0`$, which implies $`I_+>I_2>I_{}`$ according to Eq. (A2) and the ensuing discussion. Adapting the solution in Landau & Lifshitz, §37, to this situation (and our notation) we have (see their Eqs. \[37.8\]-\[37.12\]) 444There is actually a sign ambiguity in the solution given by Landau & Lifshitz, which we resolve by always choosing $`\mathrm{\Omega }_{}\widehat{𝛀}_{}`$, rather than $`\mathrm{\Omega }_{}\sqrt{1\widehat{𝛀}_+^2}`$. The Euler equations Eq. (27) require that if we choose $`\mathrm{\Omega }_{}\widehat{𝛀}_{}`$, then we should also choose the sign of the coefficient of $`\mathrm{\Omega }_2`$ to be the same as the sign of the coefficient of $`\mathrm{\Omega }_{}`$.

$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_{}}{L}}`$ $`=`$ $`{\displaystyle \frac{\widehat{𝛀}_{}\mathrm{cn}(\tau )}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}`$
$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_2}{L}}`$ $`=`$ $`\widehat{𝛀}_{}\mathrm{sn}(\tau )\sqrt{{\displaystyle \frac{2\mathrm{\Delta }(1\mathrm{\Delta })}{[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2](1\mathrm{\Delta }_0)(\mathrm{\Delta }+\mathrm{\Delta }_0)}}}`$
$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_+}{L}}`$ $`=`$ $`{\displaystyle \frac{\widehat{𝛀}_+\mathrm{dn}(\tau )}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}\widehat{𝛀}_+\mathrm{dn}(\tau ),`$ (40)
where $`\mathrm{\Delta }_03(CI_B)/2\overline{I}`$, $`\mathrm{cn}(\tau )\sqrt{1\mathrm{sn}^2(\tau )}`$, $`\mathrm{dn}(\tau )\sqrt{1k^2\mathrm{sn}^2(\tau )}`$, $`\mathrm{sn}(\tau )`$ is defined by

$$\tau =_0^{\mathrm{sn}(\tau )}\frac{dx}{\sqrt{(1x^2)(1k^2x^2)}},$$
(41)
with

$`k^2`$ $`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }\mathrm{\Delta }_0)(1\mathrm{\Delta })(1\widehat{𝛀}_+^2)}{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1+\mathrm{\Delta })\widehat{𝛀}_+^2}}{\displaystyle \frac{(\mathrm{\Delta }\mathrm{\Delta }_0)(1\widehat{𝛀}_+^2)}{(\mathrm{\Delta }+\mathrm{\Delta }_0)\widehat{𝛀}_+^2}}`$
$`{\displaystyle \frac{d\tau }{dt}}`$ $`=`$ $`{\displaystyle \frac{L\left|\widehat{𝛀}_+\right|}{\overline{I}}}\sqrt{{\displaystyle \frac{2\mathrm{\Delta }(1+\mathrm{\Delta })(\mathrm{\Delta }+\mathrm{\Delta }_0)}{(1\mathrm{\Delta }_0)(1\mathrm{\Delta }^2)[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}`$
The motion is periodic, with a dimensionless period $`4K(k^2)`$ where

$$K(k^2)=_0^1\frac{dx}{\sqrt{(1x^2)(1k^2x^2)}}.$$
(43)
Eqs. (40) and (LABEL:ktaudefpos) also imply that

$$\frac{\overline{I}^2|𝛀|^2}{L^2}=\frac{1[2\mathrm{\Delta }(\mathrm{\Delta }\mathrm{\Delta }_0)(1\widehat{𝛀}_+^2)/(1+\mathrm{\Delta })(1\mathrm{\Delta }_0)]\mathrm{sn}^2(\tau )}{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2},$$
(44)
which is independent of time up to terms $`\mathrm{\Delta }^2`$, thus validating the approximation of time independent $`I_\mathrm{\Omega }`$ used above.

One of the distinguishing features of the timing model we will develop below is that we shall not demand that $`|\widehat{𝛀}_{}|`$ be small. In particular, we shall see that when $`|I_B|C`$, $`\widehat{𝛀}_{}`$ will not be small in general, but the observable effect of the precession on pulse arrival times could still be small. This is because in the limit where magnetic distortions are far larger than crustal distortions the star tends to precess about its magnetic axis. If there were no crust at all, as in the magnetic fluid stars considered by Mestel and collaborators (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981), the star would precess exactly around its magnetic axis, which would therefore rotate uniformly. The crust breaks this symmetry, and allows the precession to be observable. We shall see this emerge in some detail when we consider timing residuals in §5.2.

For the opposite case, where the crustal deformations dominate, the precession amplitude is set by $`\widehat{𝛀}_{}`$, according to Eqs. (40). In that case, we also see that Eq. (LABEL:ktaudefpos) shows that $`k^2\widehat{𝛀}_{}^2`$. Thus, if the precession amplitude is small, so is $`k^2`$ irrespective of how triaxial the crust might be. Since $`k^2`$ governs the importance of oscillations at harmonics of the fundamental precession frequency, we see that small amplitude will imply oscillations predominantly at the fundamental if crustal deformations dominate. (This will also be true for the solution given in §4.2.3 for $`I_B<0`$, where Eq.  will also imply that oscillations at harmonics of the precession frequency are suppressed for small precession amplitude.) This is not the case when magnetic deformations dominate.

The better-known results for free precession of an axisymmtric star are recovered for $`\mathrm{\Delta }\mathrm{\Delta }_0=0=k^2`$. However, note that in the triaxial case, $`k^2`$ need not be small. We have already noted that $`1\widehat{𝛀}_+^2`$ does not have to be small, even in the minimum energy state, and from Eqs. (34) and the definition of $`\mathrm{\Delta }_0`$ we see that

$$\frac{\mathrm{\Delta }\mathrm{\Delta }_0}{\mathrm{\Delta }+\mathrm{\Delta }_0}=\frac{\sqrt{C^22CI_B\mathrm{cos}2\theta +I_B^2}(CI_B)}{\sqrt{C^22CI_B\mathrm{cos}2\theta +I_B^2}+CI_B},$$
(45)
which is not necessarily small either. For $`0<I_B/C1`$, we note that $`k^2I_B^3\mathrm{sin}^22\chi \mathrm{sin}^2\chi /4C^3`$ to lowest order in $`I_B/C`$ in the minimum energy state, so $`k^21`$ in this case. However, for $`I_BC`$, we shall see that $`k^2I_B/C1`$, and in general, since $`\mathrm{\Delta }_0<0`$ for $`I_B>C`$, it is possible for $`k^2`$ to exceed one rather generally in that regime. The solution to the Euler problem is still given in terms of elliptic equations when $`k^2>1`$, but we have to make the replacements

$`\mathrm{sn}(\tau )`$ $``$ $`{\displaystyle \frac{\mathrm{sn}(k\tau )}{k}}`$
$`\mathrm{cn}(\tau )`$ $``$ $`\mathrm{dn}(k\tau )`$
$`\mathrm{dn}(\tau )`$ $``$ $`\mathrm{cn}(k\tau )`$ (46)
in Eqs. (40), and $`\mathrm{dn}^2(k\tau )=1k^2\mathrm{sn}^2(k\tau )`$. The explicit solution for $`k^2>1`$ is

$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_{}}{L}}`$ $`=`$ $`{\displaystyle \frac{\widehat{𝛀}_{}\mathrm{dn}(\stackrel{~}{\tau })}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}`$
$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_2}{L}}`$ $`=`$ $`\widehat{𝛀}_+\mathrm{sign}(\widehat{𝛀}_{})\mathrm{sn}(\stackrel{~}{\tau })\sqrt{{\displaystyle \frac{2\mathrm{\Delta }(1+\mathrm{\Delta })}{(1\mathrm{\Delta }_0)(\mathrm{\Delta }\mathrm{\Delta }_0)[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}`$
$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_+}{L}}`$ $`=`$ $`{\displaystyle \frac{\widehat{𝛀}_+\mathrm{cn}(\stackrel{~}{\tau })}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}},`$ (47)
where $`\mathrm{sign}(\widehat{𝛀}_{})`$ is the sign of $`\widehat{𝛀}_{}`$, which may be positive or negative, $`\mathrm{dn}^2(\stackrel{~}{\tau })1\stackrel{~}{k}^2\mathrm{sn}^2(\stackrel{~}{\tau })`$, and

$`{\displaystyle \frac{d\stackrel{~}{\tau }}{dt}}`$ $`=`$ $`{\displaystyle \frac{L}{\overline{I}}}\sqrt{{\displaystyle \frac{2\mathrm{\Delta }(\mathrm{\Delta }\mathrm{\Delta }_0)(1\widehat{𝛀}_+^2)}{(1\mathrm{\Delta }_0)(1+\mathrm{\Delta })[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}`$
$`\stackrel{~}{k}^2`$ $`=`$ $`{\displaystyle \frac{1}{k^2}}={\displaystyle \frac{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1+\mathrm{\Delta })\widehat{𝛀}_+^2}{(\mathrm{\Delta }\mathrm{\Delta }_0)(1\mathrm{\Delta })(1\widehat{𝛀}_+^2)}}.`$ (48)
This somewhat unfamiliar case is not treated in Landau & Lifshitz, but is found eaasily from the equations given there, and corresponds simply to the transformations in Eqs. (46). Physically, it turns out to be important for a substantially prolate figure, which is what happens when $`I_B`$ is large and positive.

An important difference between axisymmetric and triaxial precession is that even though both are periodic, only the axisymmetric case is precisely sinusoidal. The functions $`\mathrm{cn}(\tau )`$ and $`\mathrm{sn}(\tau )`$ contain all odd harmonics of $`\pi \tau /2K(k^2)`$, whereas the function $`\mathrm{dn}(\tau )`$ contains all even harmonics. When $`k^2`$ is small, the expansions are dominated by their leading terms, but, as we have seen, the general case does not demand small values of $`k^2`$. The key parameter in the expansions is

$$q(k^2)\mathrm{exp}\left[\frac{\pi K^{}(k^2)}{K(k^2)}\right],$$
(49)
where

$$K^{}(k^2)=K(1k^2)$$
(50)
(Abramowitz & Stegun, Eqs. \[16.1.1\], \[16.23.1-3\]). Since $`q(k^2)<1`$ by its definition, only the first few harmonics ought to be prominent in the solution, but as long as $`k^2`$ is not especially small compared to one, the amplitudes of the first few harmonics ought to be roughly comparable. This feature is consistent with the observed properties of PSR 1828-11 (Stairs, Lyne & Shemar 2000). If this interpretation of the observations is correct, then the fundamental precession period must be 1000 (or perhaps 2000) days. 555The longer period appears in the Fourier analysis of the timing residuals, but not in those of the period, period derivative, or shape variations, given in Fig. 2 of Stairs, Lyne & Shemar 2000. In any event, the data shown there only span about 2000 days, so this feature may not be so well-established.

#### 4.2.3 Nonlinear Solution: $`I_B<0`$

For $`I_B<0`$, the solution is analogous to what was given in §4.2.2, except that, according to Eq. (A4) and ensuing discussion, $`I_+>I_{}>I_2`$. Instead of Eqs. (40) we have

$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_2}{L}}`$ $`=`$ $`\widehat{𝛀}_{}\left\{{\displaystyle \frac{2\mathrm{\Delta }(1\mathrm{\Delta })}{[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2](1\mathrm{\Delta }_0)(\mathrm{\Delta }+\mathrm{\Delta }_0)}}\right\}^{1/2}\mathrm{cn}(\tau )`$
$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_{}}{L}}`$ $`=`$ $`\left[{\displaystyle \frac{\widehat{𝛀}_{}\mathrm{sn}(\tau )}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}\right]`$
$`{\displaystyle \frac{\overline{I}\mathrm{\Omega }_+}{L}}`$ $`=`$ $`\left\{{\displaystyle \frac{[\mathrm{\Delta }_0+\mathrm{\Delta }^2+\mathrm{\Delta }(1+\mathrm{\Delta }_0)(2\widehat{𝛀}_+^21)]}{[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2](1+\mathrm{\Delta })(\mathrm{\Delta }+\mathrm{\Delta }_0)}}\right\}^{1/2}\mathrm{dn}(\tau ),`$ (51)
and instead of Eq. (LABEL:ktaudefpos) we have

$`k^2`$ $`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }_0\mathrm{\Delta })(1\mathrm{\Delta })(1\widehat{𝛀}_+^2)}{\mathrm{\Delta }_0+\mathrm{\Delta }^2+\mathrm{\Delta }(1+\mathrm{\Delta }_0)(2\widehat{𝛀}_+^21)}}`$
$`{\displaystyle \frac{d\tau }{dt}}`$ $`=`$ $`{\displaystyle \frac{L}{\overline{I}}}\sqrt{{\displaystyle \frac{2\mathrm{\Delta }[\mathrm{\Delta }_0+\mathrm{\Delta }^2+\mathrm{\Delta }(1+\mathrm{\Delta }_0)(2\widehat{𝛀}_+^21)]}{(1\mathrm{\Delta }_0)(1\mathrm{\Delta }^2)[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}.`$ (52)
The definitions of $`\mathrm{sn}(\tau )`$, $`\mathrm{cn}(\tau )`$ and $`\mathrm{dn}(\tau )`$ are the same as before. The magnitude of the angular velocity is now

$$\frac{\overline{I}^2|𝛀|^2}{L^2}=\frac{1+[2\mathrm{\Delta }(\mathrm{\Delta }_0\mathrm{\Delta })(1\widehat{𝛀}_+^2)/(1+\mathrm{\Delta })(1\mathrm{\Delta }_0)]\mathrm{cn}^2(\tau )}{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}.$$
(53)
Eq. (53) implies that $`|𝛀|`$ is independent of time up to terms $`\mathrm{\Delta }^2`$, just as we found from Eq. (44) for $`I_B>0`$. In this case, $`k^2<1`$, so these solutions suffice.

### 4.3 Effect of the Spindown Torque

The solutions given in §4.2.2 and 4.2.3 are for free precession. Radiopulsars spin down as a result of electromagnetic radiation, on a timescale $`t_{sd}P/2\dot{P}`$ that is long compared with the precession periods of interest: for PSR 1828-11, $`t_{sd}1.1\times 10^5`$ years. A careful examination of the Euler equations with spindown included shows that the precession is described by the torque-free solutions up to small corrections (just as was found by Link & Epstein 2001). We also note here that the quantity $`2E\overline{I}/L^2`$ only changes on a still longer timescale, $`t_{sd}/\mathrm{\Delta }`$. However, sinusoidal variations of the spindown torque result from the precession, since the angle between $`\widehat{𝒃}`$ and $`𝛀`$ varies with time. The amplitude of these variations need not be small, and the associated timing residuals can dominate (Cordes 1993). In fact, we shall see that they are dominant, just as was found by Link & Epstein (2001) for precession of an axisymmetric neutron star.

For vacuum magnetic dipole radiation,

$$\frac{d𝑳}{dt}=K\widehat{𝒃}\mathbf{\times }(\widehat{𝒃}\mathbf{\times }𝛀),$$
(54)
where $`K`$ is a constant to sufficient accuracy, and the magnitude of the angular velocity changes at a rate

$$\frac{d\mathrm{\Omega }}{dt}=\frac{1}{\overline{I}}\frac{dL}{dt}\frac{\mathrm{\Omega }_0}{2t_{sd}}\left(1+\frac{\zeta }{\mathrm{sin}^2\chi }\right),$$
(55)
where $`\mathrm{\Omega }_0`$ is the angular frequency at some reference time $`t=0`$; in the absence of precession, $`\zeta =0`$. We use the conventional definition of the spindown time as $`t_{sd}=P/2\dot{P}=\mathrm{\Omega }/2\dot{\mathrm{\Omega }}`$. We can use the solutions to the Euler problem found above to find the time dependent $`\zeta `$ case by case:

$`\zeta `$ $`=`$ $`\widehat{𝛀}_{}^2\mathrm{sn}^2(\tau )\left[{\displaystyle \frac{(\mathrm{\Delta }\mathrm{\Delta }_0)\mathrm{cos}^2\sigma _+}{\mathrm{\Delta }+\mathrm{\Delta }_0}}+\mathrm{sin}^2\sigma _+\right]2\widehat{𝛀}_+\widehat{𝛀}_{}\mathrm{cos}\sigma _+\mathrm{sin}\sigma _+[1\mathrm{dn}(\tau )\mathrm{cn}(\tau )]`$
$`\zeta `$ $`=`$ $`\widehat{𝛀}_+^2\mathrm{sn}^2(\stackrel{~}{\tau })\left[{\displaystyle \frac{(\mathrm{\Delta }+\mathrm{\Delta }_0)\mathrm{sin}^2\sigma _+}{\mathrm{\Delta }\mathrm{\Delta }_0}}+\mathrm{cos}^2\sigma _+\right]2\widehat{𝛀}_+\widehat{𝛀}_{}\mathrm{sin}\sigma _+\mathrm{cos}\sigma _+[1\mathrm{dn}(\stackrel{~}{\tau })\mathrm{cn}(\stackrel{~}{\tau })]`$
$`\zeta `$ $`=`$ $`\widehat{𝛀}_{}^2\mathrm{cn}^2(\tau )\left[{\displaystyle \frac{(\mathrm{\Delta }\mathrm{\Delta }_0)\mathrm{cos}^2\sigma _+}{\mathrm{\Delta }+\mathrm{\Delta }_0}}+\mathrm{sin}^2\sigma _+\right]`$ (56)
$`+2\widehat{𝛀}_{}\mathrm{sin}\sigma _+\mathrm{cos}\sigma _+\left[\mathrm{dn}(\tau )\mathrm{sn}(\tau )\sqrt{{\displaystyle \frac{\mathrm{\Delta }_0+\mathrm{\Delta }(2\widehat{𝛀}_+^21)}{\mathrm{\Delta }+\mathrm{\Delta }_0}}}\widehat{𝛀}_+\right],`$
for, respectively, $`I_B>0,k^2<1`$, $`I_B>0,\stackrel{~}{k}^2=1/k^2<1`$, and $`I_B<0`$. Here, we do not restrict the solution to small values of $`\widehat{𝛀}_{}^2`$, and Eqs. (55) and (56) imply oscillations of the spindown torque at both even and odd harmonics of the precession frequency. (There are also evidently zero frequency corrections but these can always be combined with $`\mathrm{sin}^2\chi `$ and factored out.) Note that for simple triaxial precession with a small tilt of the angular velocity away from the principal axis with maximum moment of inertia, the amplitude of the oscillations at $`2\omega _p`$ would be smaller, by a factor $`|\widehat{𝛀}_{}|`$, than the amplitude of the oscillations at $`\omega _p`$. For the model developed here, $`\widehat{𝛀}_{}`$ can become substantial, particularly if $`I_B<0`$ and $`|I_B|C`$.

## 5 Pulse Arrival Times

To determine pulse arrival times, we need to represent $`\widehat{𝒃}(t)`$ in the inertial frame of reference of the observer. We can define the angular momentum vector $`𝑳`$, which is conserved apart from spindown in this reference frame, to lie along the $`z`$ axis, and we can further choose to place the observer in the $`xz`$ plane. Pulses arrive when $`\widehat{𝒃}_y(t)=0`$ (actually, only half of the solutions correspond to pulse arrival – the other half might be an interpulse or else unobserved). The general problem is addressed in Appendix D, where the three different types of solutions are treated separately in Appendices D.1 and D.2 for positive and negative $`I_B`$, respectively. Approximate solutions are also derived in those appendices (Eqs. \[D17\], \[D28\] and \[D41\]), but those results only apply when the star is nearly axisymmetric and $`|\widehat{𝛀}_{}|1`$.

Here, we shall investigate two different nearly axisymmetric cases, $`|I_B/C|1`$ and $`|I_B/C|1`$. For $`|I_B|C`$, it will turn out that $`|\widehat{𝛀}_{}|1`$, and the results of Appendices D.1 and D.2 will be directly applicable. However, for $`|I_B|C`$, we have already mentioned that $`\widehat{𝛀}_{}`$ need not be small (see, for example, discussion following Eq. ).

### 5.1 Pulse Arrival Times for $`|I_B|C`$

When $`|I_B|C`$, we can apply the results in Eqs. (D17) and (D41) directly. Using the approximations given in Appendix B, and accounting for pulsar spindown using the results of § 4.3, we find that the oscillating parts of the phase residuals are

$$\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}\widehat{𝛀}_{}\mathrm{cot}\chi \left(\mathrm{sin}\tau +\frac{\mathrm{\Omega }_0\mathrm{cos}\tau }{t_{sd}\omega _p^2}\right)$$
(57)
for $`I_B>0`$, and

$$\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}\widehat{𝛀}_{}\mathrm{cot}\chi \left(\mathrm{cos}\tau +\frac{\mathrm{\Omega }_0\mathrm{sin}\tau }{t_{sd}\omega _p^2}\right)$$
(58)
for $`I_B<0`$, where

$$\widehat{𝛀}_{}\widehat{𝛀}_{}^{(0)}\frac{\widehat{𝛀}_+^{(0)}I_B\mathrm{sin}2\theta }{2C}$$
(59)
and $`\widehat{𝛀}_{}^{(0)}\mathrm{sin}(\chi \theta )`$, which is zero in the minimum energy configuration. (See Appendix B.) Thus, to lowest order in the (presumed) small quantity $`\widehat{𝛀}_{}\mathrm{cot}\chi `$, phase residuals oscillate only at $`\omega _p`$. Oscillations at higher harmonics, such as $`2\omega _p`$, have amplitudes that are smaller by additional factors of $`\widehat{𝛀}_{}`$. It is possible for these to be comparable in magnitude to the terms $`\widehat{𝛀}_{}`$ retained in Eqs. (57) and (58), but only if $`\mathrm{cot}\chi `$ is very small i.e. if $`|\widehat{𝛀}_{}\mathrm{cot}\chi |\widehat{𝛀}_{}^2`$, or $`|\mathrm{cot}\chi ||\widehat{𝛀}_{}|1`$. This was also found by Link & Epstein (2001), who required $`\chi `$ very close to $`\pi /2`$ in order for their axisymmetric precession model to account for the observed precession of PSR 1828-11. (See also Rezania 2002.) Here, we also note that $`k^2e^2\widehat{𝛀}_{}^2`$, where $`e^2`$ represents the deviation of the star from axisymmetry. Thus, nonaxisymmetric effects alone cannot introduce substantial harmonic structure in the phase residuals for small $`\widehat{𝛀}_{}^2`$ either.

The importance of spindown in the timing residuals is measure by the nondimensional parameter

$$\mathrm{\Gamma }_{sd}\frac{\mathrm{\Omega }_0}{\omega _p^2t_{sd}}=\frac{P_p^2}{2\pi P_0t_{sd}}376P_{p,1000}^2P_0^1t_{sd,5}^1$$
(60)
where $`P_p=1000P_{p,1000}`$ days is the precession period, $`P_0`$ is the pulsar period in seconds, and $`t_{sd}=10^5t_{sd,5}`$ years. For PSR 1828-11, we have $`\mathrm{\Gamma }_{sd}844`$, so the pulsar spindown dominates the oscillatory terms.

### 5.2 Pulse Arrival Times for $`|I_B|C`$

When $`|I_B|C`$, the moment of inertia tensor is once again approximately axisymmetric, but neither $`\widehat{𝛀}_{}`$ nor $`\widehat{𝛀}_+`$ has to be small. Thus, the expansions in Appendices D.1 and D.2 are not applicable, and we shall have to solve the timing equation in a different way. For doing so, Eqs. (C2) and (C4) prove to be useful. To keep the notation compact, we define the nondimensional parameter $`\stackrel{~}{C}C/|I_B|`$.

Let us consider the two possible sign choices separately. For $`I_B>0`$, we use Eq. (D24) to evaluate $`\widehat{𝒃}_y`$. The problem is simplified since, from Eq. (C2), $`\widehat{𝒃}_+\stackrel{~}{C}\mathrm{sin}2\theta /21`$. Consequently, to first order in $`\stackrel{~}{C}`$, we find

$$\widehat{𝒃}_y\mathrm{sin}\chi \mathrm{cos}\varphi \frac{\stackrel{~}{C}\mathrm{sin}2\theta }{2}\left(\mathrm{sin}\tau \mathrm{sin}\varphi +\mathrm{cos}\chi \mathrm{cos}\varphi \right),$$
(61)
so pulses arrive when

$$\varphi \left(2n+\frac{1}{2}\right)\pi \frac{\stackrel{~}{C}\mathrm{sin}2\theta \mathrm{sin}\stackrel{~}{\tau }}{2\mathrm{sin}\chi };$$
(62)
mapping from phase to time implies

$$_0^t\frac{dtL}{\overline{I}}\left(2n+\frac{1}{2}\right)\pi \frac{\stackrel{~}{C}\mathrm{sin}2\theta \mathrm{sin}\stackrel{~}{\tau }}{2\mathrm{sin}\chi }+\frac{\stackrel{~}{C}\mathrm{sin}^2\theta \mathrm{sin}2\stackrel{~}{\tau }}{4\mathrm{cos}\chi }.$$
(63)
Taking account of pulsar spindown we find that the oscillatory part of the timing residuals is

$$\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}\stackrel{~}{C}\left[\frac{\mathrm{sin}^2\theta \mathrm{sin}2\stackrel{~}{\tau }}{4\mathrm{cos}\chi }\frac{\mathrm{sin}2\theta \mathrm{sin}\stackrel{~}{\tau }}{2\mathrm{sin}\chi }+\mathrm{\Gamma }_{sd}\left(\frac{\mathrm{sin}^2\theta \mathrm{cos}2\stackrel{~}{\tau }}{16}+\frac{\mathrm{cos}\chi \mathrm{sin}2\theta \mathrm{cos}\stackrel{~}{\tau }}{2\mathrm{sin}\chi }\right)\right].$$
(64)
Note that the phase residuals vanish as $`\stackrel{~}{C}0`$ even though the star still precesses. Moreover, there are oscillations at both $`\omega _p`$ and $`2\omega _p`$ whose amplitudes may be comparable, in agreement with observations of PSR 1828-11.

For $`I_B<0`$, we find, to order $`\stackrel{~}{C}`$, we find that pulses arrive when

$$\frac{dtL}{\overline{I}}\left(2n+\frac{1}{2}\right)\pi \frac{\stackrel{~}{C}\mathrm{sin}2\theta \mathrm{cos}\tau }{2\mathrm{sin}\chi }\frac{\stackrel{~}{C}\mathrm{sin}^2\theta \mathrm{sin}2\tau }{4\mathrm{cos}\chi },$$
(65)
and taking account of spindown results in oscilating timing residuals

$$\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}\stackrel{~}{C}\left[\frac{\mathrm{sin}2\theta \mathrm{cos}\tau }{2\mathrm{sin}\chi }\frac{\mathrm{sin}^2\theta \mathrm{sin}2\tau }{4\mathrm{cos}\chi }+\mathrm{\Gamma }_{sd}\left(\frac{\mathrm{sin}^2\theta \mathrm{cos}2\tau }{16}\frac{\mathrm{sin}2\theta \mathrm{cos}\chi \mathrm{sin}\tau }{2\mathrm{sin}\chi }\right)\right].$$
(66)
Once again, this involves oscillations at both $`\omega _p`$ and $`2\omega _p`$ which can have comparable magnitudes. We see again that as $`\stackrel{~}{C}0`$, the oscillatory phase residuals disappear.

### 5.3 Application to PSR 1828-11

For PSR 1828-11, timing residuals appear to oscillate at both $`\omega _p`$ and $`2\omega _p`$ with similar amplitude. The results of §5.1 show that this situation is incompatible with small values of $`|\widehat{𝛀}_{}|`$ if $`|I_B|C`$. Moreover, we note that the results of §5.1 continue to hold as $`|I_B|0`$, so we see that equal amplitudes at $`\omega _p`$ and $`2\omega _p`$ cannot arise from a model without magnetic stresses, but with a triaxial crust, unless there is some fine-tuning of parameters (as in the axisymmetric model of Link & Epstein , which requires $`\chi `$ very close to $`\pi /2`$).

Thus, we focus on the strongly magnetic case, $`|I_B|C`$. In this case, Eqs. (64) and (66) show that it is possible for the phase residuals to oscillate with comparable amplitudes. The difference between the large and small $`|I_B|/C`$ limits is that for small $`|I_B|/C`$, the amplitude of the observed timing residuals is determined solely by $`|\widehat{𝛀}_{}|`$, but at large $`|I_B|/C`$ the amplitude is determined by $`C/|I_B|`$ primarily. Thus, by contrast to what we found for $`|I_B|C`$, small amplitude need not suppress the oscillations at $`2\omega _p`$.

For PSR 1828-11, we also know that $`\mathrm{\Gamma }_{sd}8441`$, so let us approximate the oscillatory timing residuals further as

$`\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}`$ $``$ $`{\displaystyle \frac{\stackrel{~}{C}\mathrm{\Gamma }_{sd}\mathrm{sin}^2\theta }{16}}\left(\mathrm{cos}2\stackrel{~}{\tau }+{\displaystyle \frac{16\mathrm{cos}\stackrel{~}{\tau }}{\mathrm{tan}\theta \mathrm{tan}\chi }}\right)(I_B>0)`$
$`\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}`$ $``$ $`{\displaystyle \frac{\stackrel{~}{C}\mathrm{\Gamma }_{sd}\mathrm{sin}^2\theta }{16}}\left(\mathrm{cos}2\tau {\displaystyle \frac{16\mathrm{sin}\tau }{\mathrm{tan}\chi \mathrm{tan}\theta }}\right)(I_B<0).`$ (67)
The parameter

$$u\frac{16}{\mathrm{tan}\chi \mathrm{tan}\theta }$$
(68)
governs the relative strengths of the two harmonics in the timing residuals. If the precession is in the minimum energy state, $`\theta =\chi `$, then $`u=1`$ for $`\chi 76^{}`$.

The timing residuals in Eq. (67) vary between different minimum and maximum values. For $`I_B>0`$, $`\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}`$ is minimum when $`\mathrm{cos}\stackrel{~}{\tau }=u/4`$, provided that $`|\mathrm{tan}\chi \mathrm{tan}\theta |<1`$; for $`\chi =\theta `$, this is so as long as $`\chi >63^{}`$. Presuming this to be so, the minimum value is

$$(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})_{min}=\frac{\stackrel{~}{C}\mathrm{\Gamma }_{sd}\mathrm{sin}^2\theta }{16}\left(1+\frac{u^2}{8}\right).$$
(69)
The maximum is at $`\mathrm{cos}\stackrel{~}{\tau }=1`$, where we find

$$(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})_{max}=\frac{\stackrel{~}{C}\mathrm{\Gamma }_{sd}\mathrm{sin}^2\theta }{16}\left(1+u\right).$$
(70)
The ratio of maximum to minimum timing residual is

$$\frac{(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})_{max}}{|(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})_{min}}|=\frac{1+u}{1+u^2/8}(I_B>0);$$
(71)
for $`u=1`$, this ratio is $`16/91.8`$, and the ratio is two for $`u=2`$. Observationally, the timing residuals for PSR 1828-11 appear skewed toward positive values, with a maximum about twice the magnitude of the minimum (see Fig. 1 in Stairs, Lyne & Shemar ). Similarly, for $`I_B<0`$, the maximum value of $`\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}`$ occurs when $`\mathrm{sin}\tau =u/4`$, and the minimum occurs when $`\mathrm{sin}\tau =1`$; in this case, the ratio of minimum to maximum values is

$$\frac{(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})_{max}}{|(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})_{min}|}=\frac{1+u^2/8}{1+u}(I_B<0);$$
(72)
for $`u=1`$, this ratio is $`9/160.56`$, and the ratio is two for either $`u16.9`$ or $`u0.89`$. To the extent that we expect $`u>0`$ (manifestly so for $`\chi =\theta `$), and $`u1`$ (but not $`1`$) the observed residual arrival times appear to favor $`I_B>0`$, i.e. prolate magnetic distortions.

Fig. 1 illustrates the residuals in the arrival times, $`\mathrm{\Delta }t`$ (dotted, left panel), period derivatives $`\mathrm{\Delta }\dot{P}`$ (solid, dotted) and period $`\mathrm{\Delta }P`$ (solid, right panel) for $`\theta =\chi `$ and $`u=5`$ for the prolate (upper) and oblate (lower) cases. The period an period derivatives are computed by differentiating $`\mathrm{\Omega }_0\mathrm{\Delta }t_{osc}`$ with respect to time:

$`\mathrm{\Delta }P(t)`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{\Omega }_0^2}}{\displaystyle \frac{d(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})}{dt}}={\displaystyle \frac{P_0^2}{P_p}}{\displaystyle \frac{d(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})}{d\tau }}`$
$`\mathrm{\Delta }\dot{P}(t)`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{\Omega }_0^2}}{\displaystyle \frac{d^2(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})}{dt^2}}=2\pi \left({\displaystyle \frac{P_0}{P_p}}\right)^2{\displaystyle \frac{d^2(\mathrm{\Omega }_0\mathrm{\Delta }t_{osc})}{d\tau ^2}}.`$ (73)
The agreement between these evaluations and the results plotted in Fig. 2 of Stairs, Lyne & Shemar (2000) is good, superficially, for the prolate model. Better agreement is seen for the variable period and period derivative than for the arrival times themselves; this may have been expected (e.g. Cordes 1993).

Using the results plotted in Fig. 1, we can estimate the magnitude of $`\stackrel{~}{C}`$ and the angles involved, even though these results do not constitute a true fit to the data, but just a plausible model. The curves for $`\mathrm{\Delta }\dot{P}`$ and $`\mathrm{\Delta }P`$ resemble the observational results better, so let us focus on those. Fig. 1 was prepared for $`\chi =\theta `$ and $`u=5`$, which corresponds to $`\chi =\theta =60.8^{}`$. From Fig. 2 in Stairs, Lyne & Shemar (2000), we see that the maximum values of $`\mathrm{\Delta }P`$ and $`\mathrm{\Delta }\dot{P}`$ are about 1 ns and $`0.2\times 10^{15}`$ for PSR 1828-11. For the prolate model, which resembles the observations better, Fig. 1 implies a peak value of $`460\stackrel{~}{C}`$ ns for $`\mathrm{\Delta }P`$, and $`2.6\stackrel{~}{C}\times 10^{14}`$ for $`\mathrm{\Delta }\dot{P}`$. We therefore estimate $`\stackrel{~}{C}0.002`$ from $`\mathrm{\Delta }P`$, and $`\stackrel{~}{C}0.008`$ from $`\mathrm{\Delta }\dot{P}`$. Since we have not attempted true curve fits (i.e. by varying the parameters $`u`$, $`\theta `$ and $`\stackrel{~}{C}`$) we regard this as acceptable agreement, provisionally.

## 6 Discussion

Here, we have extended previous studies of precession of neutron stars to incorporate the effects of oblique magnetic fields. We have shown that if the magnetic stresses are large enough, then steady rotation is unlikely, and the neutron star must precess. Moreover, even when the magnetic stresses are relatively weak, so that steady rotation is possible irrespective of the obliquity of the magnetic field, the minimum energy state is not the steady state. Thus, even in this case, the neutron star will precess. We argued, in §4, that the minimum energy, precessing state is a local energy minimum that applies at fixed angle between the magnetic and rotational axes. On a longer timescale, we would expect the star to seek its global energy minimum, which should correspond to either aligned or perpendicular magnetic and rotation axes, and no precession. We might expect short timescale dissipative effects to drive the system toward its local minimum, and that the global minimum is only achieved on somewhat longer timescale, perhaps as a result of electromagnetic spindown torques (e.g. Goldreich 1970).

The effective moment of inertia tensor of a neutron star with an inclined magnetic field is inherently triaxial. Consequently, the precession is periodic but not sinusoidal in time. In general, the solution for the rotational angular velocity of the star can be expanded in a Fourier series involving harmonics of the precession frequency. We have shown that at least the first few terms in such an expansion can have comparable magnitudes provided that the interior magnetic stresses are not very small.

The condition that magnetic stresses play an important role is that the magnetic-induced distortions are comparable to or larger than the distortions of the stellar crust. For precession periods of order years, the implied magnetic stresses exceed those expected from the classical Maxwell stress tensor, evaluated using the inferred dipole magnetic field strength, by a couple of orders of magnitude. However, if the interior of a neutron star contains a Type II superconductor, or else is a normal conductor but possesses large toroidal magnetic fields, the magnetic stresses are larger, and the implied distortions can be of the right order of magnitude (Jones 1975, Easson & Pethick 1977, Cutler 2002). Thus, the observation of neutron star precession can be taken as indirect evidence for enhanced magnetic stresses, due to either Type II superconductivity, or large toroidal fields.

We postpone detailed application of the ideas set forth here to PSR 1828-11 to another paper (Akgun, Epstein & Wasserman 2002). However, in §5.3 we argued that only a model with $`|I_B|C`$ can lead to time residuals that oscillate with comparable amplitude at both $`\omega _p`$ and $`2\omega _p`$. In this case, the amplitude of the observed time residuals is set by the dimensionless ratio $`\stackrel{~}{C}C/|I_B|`$, not by the tilt of the angular velocity vector away from any principal axis, which may be large.

By contrast, if the stellar distortions associated with magnetic stresses are smaller than those associated with the crust, then Eqs. (57) and (58) show that the timing residuals oscillate predominantly at $`\omega _p`$. Oscillations at $`2\omega _p`$ would be down by factors $`|\widehat{𝛀}_{}|`$, where $`\widehat{𝛀}_{}`$ is given by Eq. (59), and is $`|I_B|/C`$ in the minimum energy configuration. Moreover, we argued, in §5.1 and §5.3, that precession of a triaxial crust alone would probably not, at small precession amplitude, be capable of producing oscillations of comparable magnitude at both $`\omega _p`$ and $`2\omega _p`$, because the precession amplitude is proportional to $`|\widehat{𝛀}_{}|`$ and the oscillations at harmonics of $`\omega _p`$ are suppressed by factors $`|\widehat{𝛀}_{}|`$. Thus, a solution in which the crustal distortion is responsible for precession is unlikely to explain the data on PSR 1828-11. The model in which magnetic stresses dominates is not merely precession of a triaxial body because small amplitude phase residuals can arise even when $`|\widehat{𝛀}_{}|`$ is not small.

For a precessing star in which the distortion is due to magnetic stresses primarily, the precession period is

$$P_p\frac{P_0\overline{I}}{3I_B\mathrm{cos}\chi }\frac{492\mathrm{days}}{\beta \mathrm{cos}\chi (BH)_{27}/5},$$
(74)
where the value given is for PSR 1828-11. This is similar to the observed period, about 1000 days, provided that $`\beta \mathrm{cos}\chi `$ is not very small. Thus, if $`\chi 60^{}`$ the timing residuals could oscillate with about the right period, amplitude and relative importance of the fundamental precession period and its first harmonic. In this case, we note that the condition that there is no steady state, $`|\mathrm{sin}2\chi |\stackrel{~}{C}`$ can be satisfied: for $`\chi =60^{}`$, for example, $`\mathrm{sin}2\chi =\sqrt{3}/2`$ whereas we estimated that $`\stackrel{~}{C}0.0010.01`$ in §5.3. This possible explanation of the timing residuals for PSR 1828-11 only works if the magnetic stresses in this star are $`200`$ times larger than would be indicated by its dipole magnetic field. Thus, we conclude that either the interior is a Type II superconductor, or is a normal conductor with a toroidal field whose strength is $`10^{14}`$ G. Otherwise, the expected magnetic stresses are far smaller than is needed for this solution to apply.

We also noted, in §5.3, that the ratio of the minimum and maximum timing residuals, and the shape of the variation of $`\mathrm{\Delta }P`$ and $`\mathrm{\Delta }\dot{P}`$ seen in PSR 1828-11 appear to favor a model in which $`I_B>0`$, so that magnetic distortions are prolate. Prolate distortions would arise naturally from the stresses due to a toroidal field, with or without Type II superconductivity (Cutler 2002), but may also result if magnetic flux tubes have been transported outward in the core and accumulate at its outer boundary (Ruderman, Zhu & Chen 1998, Ruderman & Chen 1999), which could “pinch” the interior.

Crustal distortions are still needed in order for the precessional amplitude to be nonzero. In fact, to be more precise, the precessional amplitude depends on the component of the moment of inertia tensor of the crust that is not symmetric about the magnetic axis. This can be seen directly from Eqs. (64) and (66), which show that the oscillating time residuals vanish as $`\stackrel{~}{C}\mathrm{sin}\theta 0`$, where $`\theta `$ is the angle between the magnetic axis $`\widehat{𝒃}`$ and, in this axisymmetric distortion model, the symmetry axis of the relevant crustal deformation, $`\widehat{𝒌}`$. Thus, although all neutron stars precess as a consequence of their magnetic stresses for $`|I_B|C`$ in the picture advanced here, only those with sufficiently large nonaligned crustal deformations would have discernible oscillations of their timing residuals. In this sense, PSR 1828-11 may be special.

Although we have treated the basic physics of precession of an oblique rotator in some detail, we have not treated several effects that might play significant roles. We have not explicitly included either vortex line pinning or vortex drag. Link & Cutler (2002) have argued that vortex lines can unpin globally at large enough precession amplitude. For $`|I_B|C`$, as is required to explain the timing residuals in PSR 1828-11, the precession amplitude is $`\mathrm{sin}\chi `$ or $`\mathrm{cos}\chi `$ (depending on the sign of $`I_B`$), which is not small, so global unpinning is expected. Thus, we may expect that vortex lines are unpinned in neutron stars with magnetic fields that are strong enough to have $`|I_B|C`$ except for $`\chi `$ very close to either zero or $`\pi /2`$, depending on the sign of $`I_B`$. For these, vortex line drag, if weak enough, may simply serve to bring the neutron star superfluid into corotation with its crust, and drive the rotating star toward its minimum energy state. For large $`|I_B|/C`$, we have seen that precession is required, so weak vortex drag or other forms of dissipation need not prevent precession. For small values of $`|I_B|/C`$, the precession amplitude is given by Eq. (59) and can be very small; in the minimum energy state, Eq. (59) implies $`|\widehat{𝛀}_{}||I_B\mathrm{sin}2\theta |/C`$. In this case, it is possible that precession cannot overcome pinning forces, as discussed in Link & Cutler (2002), and so precession does not occur.

Theories of pulsar glitches involve the pinning, unpinning and repinning of crustal superfluid vortex lines (e.g. Anderson and Itoh 1975; Alpar et al. 1981, 1984a,b, 1993; Link, Epstein and Baym 1993). As we have seen, for large $`|I_B|/C`$, precession amplitudes are large, and vortices may be expected to unpin, but it is possible for vortex lines to remain pinned in neutron stars with $`C|I_B|`$. Thus, there could be a dichotomy between pulsars that glitch ($`C|I_B|`$) and those that precess ($`C|I_B|`$). If, in the course of a glitch, all crustal superfluid vortices were to unpin, then the star might precess briefly. Perhaps that explains the detection of damped, quasisinusoidal timing residuals in the Vela pulsar after and perhaps before its Christmas 1988 glitch (McCulloch et al. 1990).

We have kept the problem of precession of an oblique rotator as simple as possible by considering what happens when the magnetic field is axisymmetric about some axis, and the crustal distortions are also axisymmetric, but about a different axis. More realistically, both of these simplifying assumptions are likely to be violated. Most likely, the crust is not axisymmetric. When magnetic stresses dominate, we do not expect including intrinsic crustal asymmetry to alter the results found here qualitatively, since the effective moment of inertia is already triaxial here. Triaxiality of the crust, in the limit of rather small crustal distortions, would simply rotate the principal axes slightly. Furthermore, the magnetic field may have a more complicated structure than we have assumed. A substantial quadrupolar component would presumably render the contribution to the inertia tensor from magnetic stresses alone triaxial. We shall consider these complications elsewhere.

Although we have included the spindown torque in our evaluations of timing residuals, we did not include near zone electromagnetic torques (Good & Ng 1985, Melatos 1997, 1999, 2000). The principal effect of such torques would be to renormalize the moment of inertia tensor of the star. Near zone torques can play a role similar to the magnetic distortions considered here, but are smaller by a factor $`(H/B)(Rc^2/GM)`$, which is non-negligible even if $`H=B`$. However, we note here that the large magnetic distortions we propose would presumably apply to the spindown of magnetars and anomalous X-ray pulsars, in much the same fashion as proposed by Melatos (1999, 2000). We shall pursue this idea elsewhere.

We have also ignored motions of the fluid and crust of the star apart from rigid rotation. Mestel and collaborators (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981) have pointed out that the equilibrium in a fluid star with oblique magnetic field must involve fluid motions with velocities $`(\mathrm{\Omega }^2R^3/GM)\omega _pR`$. These distort the stellar magnetic field by a fractional amount $`\mathrm{\Omega }^2R^3/GM`$ over the precession period $`2\pi /\omega _p`$ (Mestel & Takhar 1972, Mestel et al. 1981, Nittman & Wood 1981). Similar motions might arise in the crust, but with magnitudes $`(\mathrm{\Omega }^2R^2/c_t^2)\omega _pR`$, where $`c_t`$ is the sound speed for transverse waves. The expected amplitude of the resulting magnetic wander is $`\mathrm{\Omega }^2R^2/c_t^2`$. Although slow, these displacements cause the magnetic field of the star to oscillate about its undisturbed, axisymmetric state, and might influence the long-term behavior of the observed spindown. We shall investigate whether there is any long-term observational signature of these motions elsewhere. In addition, it is likely that magnetic and rotational deformations of the core must also deform the crust, as they exert pressure on its inner boundary. We would expect the magnetic deformation of the core to promote crustal deformation symmetric about $`\widehat{𝒃}`$, which would not lead to observable precession, but the rotation-induced deformation need not be symmetric about $`\widehat{𝒃}`$, and should be substantial. An important question left unanswered here is whether the crustal distortion, $`C_{ij}`$, that leads to detectable precession is relatively steady, or is simply due to a seismic fluctuation. We shall explore the important issue of crustal deformations elsewhere.

Finally, we emphasize that any neutron star with strong enough core magnetic stresses ought to precess, but we may not be able to detect their precession because their magnetic axes can still rotate more or less uniformly. This is because, at small values of $`C/|I_B|`$, the neutron star precesses almost exactly about its magnetic axis, which therefore rotates almost uniformly as seen in the inertial frame. Although the precession may not be detectable readily from timing residuals for most pulsars, gravitational radiation amplitudes would be larger than would arise without enhanced internal magnetic stresses (e.g. Cutler & Thorne 2002, Cutler 2002). The distortions required for PSR 1828-11 are still smaller than would be needed for detection by LIGO, even if it were spinning faster (Brady et al. 1998). If there are young, highly magnetized neutron stars rotating rapidly, they would be the brightest emitters of gravitational radiation. Such objects have been hypothesized to be the sources of the highest energy cosmic rays (Blasi, Epstein & Olinto 2000, Arons 2002).

Partial support for this work was provided by a grant from IGPP at LANL. I thank T. Akgun, J. Cordes, R. Epstein and B. Link for comments.

## Appendix A Inequalities Among Eigenvalues

For $`I_B>0`$, we can rewrite Eq. (31) in the form

$$\delta I_\pm =\frac{3}{2}\left[CI_B\pm \sqrt{(CI_B)^2+2CI_B(1\mathrm{cos}2\theta )}\right],$$
(A1)
from which it follows that

$`\delta I_+`$ $``$ $`{\displaystyle \frac{3}{2}}\left(CI_B+|CI_B|\right)`$
$`\delta I_{}`$ $``$ $`{\displaystyle \frac{3}{2}}\left(CI_B|CI_B|\right).`$ (A2)
Eq. (A2) implies that $`\delta I_+3(CI_B)>0`$ and $`\delta I_{}<0`$ if $`I_B<C`$, and $`\delta I_+>0`$ and $`\delta I_{}<3(CI_B)<0`$ if $`I_B>C`$. Thus, if $`I_B>0`$, then $`\delta I_+>0`$ and $`\delta I_{}<0`$ irrespective of whether $`I_B>C`$ or $`I_B<C`$. Thus, for $`I_B>0`$, $`I_+>I_2>I_{}`$.

For $`I_B<0`$, we can rewrite Eq. (31) in the form

$$\delta I_\pm =\frac{3}{2}\left[C+|I_B|\pm \sqrt{(C|I_B|)^2+2C|I_B|(1+\mathrm{cos}2\theta )}\right],$$
(A3)
from which it follows that

$`\delta I_+`$ $``$ $`{\displaystyle \frac{3}{2}}\left(C+|I_B|+|C|I_B||\right)`$
$`\delta I_{}`$ $``$ $`{\displaystyle \frac{3}{2}}\left(C+|I_B||C+|I_B||\right)=0.`$ (A4)
Eq. (A4) implies that $`\delta I_+>3C`$ if $`|I_B|<C`$ and $`\delta I_+>3|I_B|`$ if $`|I_B|>C`$. Thus, we see that for $`I_B<0`$, $`I_+>I_{}>I_2`$.

## Appendix B Approximate Results for $`|I_B|C`$

For $`|I_B|C`$, we can expand these results to find

$`\mathrm{sin}\sigma _+=\widehat{𝒃}_{}`$ $``$ $`\mathrm{sin}\theta \left[1+{\displaystyle \frac{I_B\mathrm{cos}^2\theta }{C}}+{\displaystyle \frac{I_B^2(5\mathrm{cos}^4\theta 3\mathrm{cos}^2\theta )}{2C^2}}+\mathrm{}\right]`$
$`\mathrm{cos}\sigma _+=\widehat{𝒃}_+`$ $``$ $`\mathrm{cos}\theta \left[1{\displaystyle \frac{I_B\mathrm{sin}^2\theta }{C}}+{\displaystyle \frac{I_B^2(5\mathrm{sin}^4\theta 3\mathrm{sin}^2\theta )}{2C^2}}+\mathrm{}\right]`$
$`\delta I_+`$ $``$ $`3C\left[1{\displaystyle \frac{I_B(1+\mathrm{cos}2\theta )}{2C}}+{\displaystyle \frac{I_B^2\mathrm{sin}^22\theta }{4C^2}}+\mathrm{}\right]`$
$`\delta I_{}`$ $``$ $`{\displaystyle \frac{3I_B(1\mathrm{cos}2\theta )}{2}}{\displaystyle \frac{3I_B^2\mathrm{sin}^22\theta }{4C^2}}+\mathrm{}.`$ (B1)
Using Eq. (B1) we also find that

$`\widehat{𝛀}_+`$ $``$ $`\widehat{𝛀}_+^{(0)}+{\displaystyle \frac{\widehat{𝛀}_{}^{(0)}I_B\mathrm{sin}2\theta }{2C}}+\mathrm{}`$
$`\widehat{𝛀}_{}`$ $``$ $`\widehat{𝛀}_{}^{(0)}{\displaystyle \frac{\widehat{𝛀}_+^{(0)}I_B\mathrm{sin}2\theta }{2C}}+\mathrm{}`$
$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{3C}{2\overline{I}}}\left(1{\displaystyle \frac{I_B\mathrm{cos}2\theta }{C}}+\mathrm{}\right)`$
$`\mathrm{\Delta }_0`$ $`=`$ $`{\displaystyle \frac{3C}{2\overline{I}}}\left(1{\displaystyle \frac{I_B}{C}}\right),`$ (B2)
where $`\widehat{𝛀}_+^{(0)}\mathrm{cos}(\chi \theta )`$ and $`\widehat{𝛀}_{}^{(0)}\mathrm{sin}(\chi \theta )`$.

## Appendix C Approximate Results for $`|I_B|C`$

For $`|I_B|C`$, we find

$`1\pm {\displaystyle \frac{\eta }{\sqrt{1+\eta ^2}}}`$ $``$ $`1\pm s_B\left[1{\displaystyle \frac{C^2\mathrm{sin}^22\theta }{2I_B^2}}+\mathrm{}\right]`$
$`\delta I_+`$ $``$ $`{\displaystyle \frac{3|I_B|}{2}}\left[1s_B+{\displaystyle \frac{C(1s_B\mathrm{cos}2\theta )}{|I_B|}}+{\displaystyle \frac{C^2\mathrm{sin}^22\theta }{2I_B^2}}+\mathrm{}\right]`$
$`\delta I_{}`$ $``$ $`{\displaystyle \frac{3|I_B|}{2}}\left[1s_B+{\displaystyle \frac{C(1+s_B\mathrm{cos}2\theta )}{|I_B|}}{\displaystyle \frac{C^2\mathrm{sin}^22\theta }{2I_B^2}}+\mathrm{}\right].`$ (C1)
where $`s_BI_B/|I_B|`$ is the sign of $`I_B`$. For $`I_B>0`$ these results imply

$`\mathrm{sin}\sigma _+=\widehat{𝒃}_{}`$ $``$ $`1{\displaystyle \frac{C^2\mathrm{sin}^22\theta }{8}}+\mathrm{}`$
$`\mathrm{cos}\sigma _+=\widehat{𝒃}_+`$ $``$ $`{\displaystyle \frac{C\mathrm{sin}2\theta }{2I_B}}+\mathrm{}`$
$`\widehat{𝛀}_+`$ $``$ $`\mathrm{sin}\chi +{\displaystyle \frac{C\mathrm{cos}\chi \mathrm{sin}2\theta }{2I_B}}+\mathrm{}`$
$`\widehat{𝛀}_{}`$ $``$ $`\mathrm{cos}\chi +{\displaystyle \frac{C\mathrm{sin}\chi \mathrm{sin}2\theta }{2I_B}}+\mathrm{}`$
$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{3I_B}{2\overline{I}}}\left(1{\displaystyle \frac{C\mathrm{cos}2\theta }{I_B}}+\mathrm{}\right)`$
$`\stackrel{~}{e}^2`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }+\mathrm{\Delta }_0}{\mathrm{\Delta }\mathrm{\Delta }_0}}{\displaystyle \frac{C\mathrm{sin}^2\theta }{I_B}}`$
$`\omega _p`$ $``$ $`{\displaystyle \frac{2\mathrm{\Delta }L\mathrm{cos}\chi }{\overline{I}}},`$ (C2)
and for $`I_B<0`$ the same results imply

$`\mathrm{sin}\sigma _+=\widehat{𝒃}_{}`$ $``$ $`{\displaystyle \frac{C\mathrm{sin}2\theta }{2|I_B|}}+\mathrm{}`$
$`\mathrm{cos}\sigma _+=\widehat{𝒃}_+`$ $``$ $`1{\displaystyle \frac{C^2\mathrm{sin}^22\theta }{I_B^2}}+\mathrm{}`$
$`\widehat{𝛀}_+`$ $``$ $`\mathrm{cos}\chi +{\displaystyle \frac{C\mathrm{sin}\chi \mathrm{sin}2\theta }{2|I_B|}}+\mathrm{}`$
$`\widehat{𝛀}_{}`$ $``$ $`\mathrm{sin}\chi {\displaystyle \frac{C\mathrm{cos}\chi \mathrm{sin}2\theta }{2|I_B|}}+\mathrm{}`$
$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{3|I_B|}{2\overline{I}}}\left(1+{\displaystyle \frac{C\mathrm{cos}2\theta }{|I_B|}}+\mathrm{}\right)`$
$`e^2`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }_0\mathrm{\Delta }}{2\mathrm{\Delta }}}{\displaystyle \frac{C\mathrm{sin}^2\theta }{|I_B|}}+\mathrm{}nonumber`$ (C3)
$`\omega _p`$ $``$ $`{\displaystyle \frac{2\mathrm{\Delta }L\sqrt{e^2+\widehat{𝛀}_+^2}}{I}}{\displaystyle \frac{2\mathrm{\Delta }L\mathrm{cos}\chi }{\overline{I}}}.`$ (C4)
## Appendix D Timing Solution

To find pulse arrival times, we need to determine the motion of $`\widehat{𝒃}`$ in the inertial reference frame of the observer. To do this, we need the Euler angle rotation from the rotating frame of reference to the inertial frame; for an arbitrary vector $`𝑽`$ this is (see e.g. Goldstein Eq. \[4-47\])

$`𝑽_x`$ $`=`$ $`𝑽_{\widehat{1}}(\mathrm{cos}\psi \mathrm{cos}\varphi \mathrm{cos}\alpha \mathrm{sin}\psi \mathrm{sin}\varphi )𝑽_{\widehat{2}}(\mathrm{sin}\psi \mathrm{cos}\varphi +\mathrm{cos}\alpha \mathrm{cos}\psi \mathrm{sin}\varphi )+𝑽_{\widehat{3}}\mathrm{sin}\alpha \mathrm{sin}\varphi `$
$`𝑽_y`$ $`=`$ $`𝑽_{\widehat{1}}(\mathrm{cos}\psi \mathrm{sin}\varphi +\mathrm{cos}\alpha \mathrm{sin}\psi \mathrm{cos}\varphi )𝑽_{\widehat{2}}(\mathrm{sin}\psi \mathrm{sin}\varphi \mathrm{cos}\alpha \mathrm{cos}\psi \mathrm{cos}\varphi )𝑽_{\widehat{3}}\mathrm{sin}\alpha \mathrm{cos}\varphi `$
$`𝑽_z`$ $`=`$ $`𝑽_{\widehat{1}}\mathrm{sin}\alpha \mathrm{sin}\psi +𝑽_{\widehat{2}}\mathrm{sin}\alpha \mathrm{cos}\psi +𝑽_{\widehat{3}}\mathrm{cos}\alpha ,`$ (D1)
where $`\alpha ,\varphi ,\psi `$ are the Euler angles defined in Fig. 47 of Landau & Lifshitz, §35, except that, to avoid confusion with our definition of $`\theta `$ as the angle between $`\widehat{𝒌}`$ and $`\widehat{𝒃}`$, we label their Euler angle $`\theta `$ as $`\alpha `$. We assume that $`𝑳`$ is along the $`\widehat{𝒆}_𝒛`$ direction in the inertial frame. We can then determine the two angles $`\alpha `$ and $`\psi `$ from the equations

$`L_{\widehat{3}}`$ $`=`$ $`L\mathrm{cos}\alpha `$
$`L_{\widehat{1}}`$ $`=`$ $`L\mathrm{sin}\alpha \mathrm{sin}\psi `$
$`L_{\widehat{2}}`$ $`=`$ $`L\mathrm{sin}\alpha \mathrm{cos}\psi ;`$ (D2)
the third Euler angle $`\varphi `$ is not determined by these relations, but can be found from

$$\frac{d\varphi }{dt}=L\left(\frac{I_{\widehat{1}}\mathrm{\Omega }_{\widehat{1}}^2+I_{\widehat{2}}\mathrm{\Omega }_{\widehat{2}}^2}{I_{\widehat{1}}^2\mathrm{\Omega }_{\widehat{1}}^2+I_{\widehat{2}}^2\mathrm{\Omega }_{\widehat{2}}^2}\right).$$
(D3)
The choice of axes $`\widehat{𝒆}_{\widehat{1},\widehat{2},\widehat{3}}`$ in the rotating frame of reference is somewhat arbitrary, and we shall make three different choices below, as the situation demands.

### D.1 Pulse Arrival Times for $`I_B>0`$

#### D.1.1 $`I_B>0`$ and $`k^2<1`$

For $`I_B>0`$ and $`k^2<1`$, we choose $`𝑽_{\widehat{1}}=𝑽_{}`$, $`𝑽_{\widehat{2}}=𝑽_2`$, and $`𝑽_{\widehat{3}}=𝑽_+`$; then Eqs. (D1) imply (Landau & Lifshitz Eq. \[37.15\])

$`\mathrm{cos}\alpha `$ $`=`$ $`{\displaystyle \frac{\widehat{𝛀}_+(1+\mathrm{\Delta })\mathrm{dn}(\tau )}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}\widehat{𝛀}_+\mathrm{dn}(\tau )`$
$`\mathrm{tan}\psi `$ $`=`$ $`{\displaystyle \frac{\mathrm{cn}(\tau )}{\mathrm{sn}(\tau )}}\sqrt{{\displaystyle \frac{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1\mathrm{\Delta })}{2\mathrm{\Delta }(1\mathrm{\Delta }_0)}}}{\displaystyle \frac{\mathrm{cn}(\tau )}{\mathrm{sn}(\tau )}}\sqrt{{\displaystyle \frac{(\mathrm{\Delta }+\mathrm{\Delta }_0)}{2\mathrm{\Delta }}}},`$ (D4)
where the approximations are to leading order in $`\mathrm{\Delta }`$. The remaining Euler angle $`\varphi `$ evolves according to (Landau & Lifshitz, Eq. \[37.16\])

$`{\displaystyle \frac{d\varphi }{dt}}`$ $`=`$ $`{\displaystyle \frac{L}{\overline{I}}}\left[{\displaystyle \frac{\mathrm{\Delta }+\mathrm{\Delta }_0+(\mathrm{\Delta }\mathrm{\Delta }_0)\mathrm{sn}^2(\tau )}{(1\mathrm{\Delta })(\mathrm{\Delta }+\mathrm{\Delta }_0)+(\mathrm{\Delta }\mathrm{\Delta }_0)(1+\mathrm{\Delta })\mathrm{sn}^2(\tau )}}\right]`$ (D5)
$``$ $`{\displaystyle \frac{L}{\overline{I}}}\left[1+{\displaystyle \frac{\mathrm{\Delta }[\mathrm{\Delta }+\mathrm{\Delta }_0(\mathrm{\Delta }\mathrm{\Delta }_0)\mathrm{sn}^2(\tau )]}{\mathrm{\Delta }+\mathrm{\Delta }_0+(\mathrm{\Delta }\mathrm{\Delta }_0)\mathrm{sn}^2(\tau )}}\right],`$
which is uniform up to corrections $`\mathrm{\Delta }`$. Note that in perfect axisymmetry, $`𝛀_z`$ and $`d\varphi /dt`$ are independent of time, another difference between the magnetic case, which is inherently triaxial, and precession with an axisymmetric crust. When the star is nearly axisymmetric,

$$\frac{d\varphi }{dt}\frac{L}{\overline{I}}\left[1+\mathrm{\Delta }(\mathrm{\Delta }\mathrm{\Delta }_0)\mathrm{sn}^2(\tau )\right];$$
(D6)
thus, $`d\varphi /dt`$ oscillates at twice the precession frequency in this limit. Associated with the time development of $`d\varphi /dt`$ would also be variability of the pulsar spindown, at even harmonics of the precession frequency. The amplitude of the main variation would be of order $`(\mathrm{\Delta }\mathrm{\Delta }_0)P_0/P_pϵ(P_0/P_p)^2`$, where $`ϵ(\mathrm{\Delta }\mathrm{\Delta }_0)/2\mathrm{\Delta }`$ is the fractional deviation from axisymmetry; successive harmonics would be smaller by powers of $`ϵ`$. For PSR 1828-11, we would have $`(\mathrm{\Delta }\mathrm{\Delta }_0)P_0/P_p3\times 10^{16}ϵ`$. For comparison, electromagnetic spindown produces oscillations with an amplitude $`|\widehat{𝛀}_{}|P_0/\tau _{sd}1.3\times 10^{13}|\widehat{𝛀}_{}|`$, where $`\tau _{sd}`$ is the spindown timescale, which is considerably larger (Link & Epstein 2001).

Applying Eqs. (D1) to $`\widehat{𝒃}=\widehat{𝒃}_+\widehat{𝒆}_\mathbf{+}+\widehat{𝒃}_{}\widehat{𝒆}_{\mathbf{}}`$, and using Eqs. (D4) we find

$`\widehat{𝒃}_y`$ $`=`$ $`{\displaystyle \frac{\widehat{𝒃}_{}}{\sqrt{1+\frac{\mathrm{sn}^2(\tau )(\mathrm{\Delta }\mathrm{\Delta }_0)(1+\mathrm{\Delta })}{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1\mathrm{\Delta })}}}}\left[\mathrm{sn}(\tau )\mathrm{sin}\varphi \sqrt{{\displaystyle \frac{2\mathrm{\Delta }(1\mathrm{\Delta }_0)}{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1\mathrm{\Delta })}}}+{\displaystyle \frac{\widehat{𝛀}_+(1+\mathrm{\Delta })\mathrm{dn}(\tau )\mathrm{cn}(\tau )\mathrm{cos}\varphi }{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}\right]`$
$`{\displaystyle \frac{\widehat{𝒃}_+\widehat{𝛀}_{}(1\mathrm{\Delta })\mathrm{cos}\varphi }{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}\sqrt{1+{\displaystyle \frac{\mathrm{sn}^2(\tau )(\mathrm{\Delta }\mathrm{\Delta }_0)(1+\mathrm{\Delta })}{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1\mathrm{\Delta })}}}`$
No approximations have been made in Eq. (LABEL:bhatinertialpos); in fact, there is also no explicit reference to magnetic distortions here, and so these results apply to triaxial stars in general. We note here that

$`\widehat{𝒃}_+`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[1+{\displaystyle \frac{(C\mathrm{cos}2\theta I_B)}{\sqrt{C^22I_BC\mathrm{cos}2\theta +I_B^2}}}\right]^{1/2}`$
$`\widehat{𝒃}_{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[1{\displaystyle \frac{(C\mathrm{cos}2\theta I_B)}{\sqrt{C^22I_BC\mathrm{cos}2\theta +I_B^2}}}\right]^{1/2},`$ (D8)
and in the minimum energy state we can take $`\theta \chi `$ to lowest order in distortions. When $`C`$ is large compared with $`I_B`$, these reduce to $`\widehat{𝒃}_+\mathrm{cos}\chi `$ and $`\widehat{𝒃}\mathrm{sin}\chi `$.

Pulse arrival times are found from the condition that $`\widehat{𝒃}_y=0`$ if we assume that the observer is in the $`xz`$ plane. Let us simplify the notation by introducing the parameters

$`e^2`$ $``$ $`{\displaystyle \frac{(\mathrm{\Delta }\mathrm{\Delta }_0)(1+\mathrm{\Delta })}{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1\mathrm{\Delta })}}`$
$`q_\pm `$ $`=`$ $`{\displaystyle \frac{1\pm \mathrm{\Delta }}{\sqrt{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}};`$ (D9)
$`e^2`$ measures the departure from axisymmetry, and $`q_\pm `$ are one to lowest order in the distortions. We can simplify further by defining $`\widehat{𝚺}_\pm =q_\pm \widehat{𝛀}_\pm `$; then it follows that $`\widehat{𝚺}_{+}^{}{}_{}{}^{2}+\widehat{𝚺}_{}^{}{}_{}{}^{2}=1`$. In terms of these, the pulse arrival times are the solution to

$$0=\frac{\widehat{𝒃}_{}[\sqrt{1+e^2}\mathrm{sn}(\tau )\mathrm{sin}\varphi +\widehat{𝚺}_+\mathrm{dn}(\tau )\mathrm{cn}(\tau )\mathrm{cos}\varphi ]}{\sqrt{1+e^2\mathrm{sn}^2(\tau )}}\widehat{𝒃}_+\widehat{𝚺}_{}\mathrm{cos}\varphi \sqrt{1+e^2\mathrm{sn}^2(\tau )};$$
(D10)
we can rewrite this as

$`\mathrm{cos}(\varphi \tau )`$ $`=`$ $`[\mathrm{sin}\tau \sqrt{1+e^2}\mathrm{sn}(\tau )]\mathrm{sin}\varphi `$ (D11)
$`+\left\{\mathrm{cos}\tau \widehat{𝚺}_+\mathrm{cn}(\tau )\mathrm{dn}(\tau )+{\displaystyle \frac{\widehat{𝒃}_+}{\widehat{𝒃}_{}}}\widehat{𝚺}_{}[1+e^2\mathrm{sn}^2(\tau )]\right\}\mathrm{cos}\varphi .`$
Define the functions

$`e^2F_S(\tau )`$ $`=`$ $`\mathrm{sin}\tau \sqrt{1+e^2}\mathrm{sn}(\tau )`$
$`e^2F_C(\tau )`$ $`=`$ $`\mathrm{cos}\tau \mathrm{cn}(\tau )\mathrm{dn}(\tau );`$ (D12)
then the pulse arrival times are the solutions of the equation

$`\mathrm{cos}(\varphi \tau )=e^2F_S(\tau )\mathrm{sin}\varphi +\{e^2F_C(\tau )\widehat{𝚺}_++(1\widehat{𝚺}_+)\mathrm{cos}\tau `$
$`+{\displaystyle \frac{\widehat{𝒃}_+}{\widehat{𝒃}_{}}}\widehat{𝚺}_{}[1+e^2\mathrm{sn}^2(\tau )]\}\mathrm{cos}\varphi .`$ (D13)
In general, we can solve Eq. (D13) numerically, but a good analytic approximation can be found if the RHS is small.

The terms on the RHS of Eq. (D13) will be small if the deviation from axisymmetry is small and the rotation is nearly aligned with one of the principal axes of the inertia tensor. If we assume that RHS is small, then we can develop an approximate solution by first substituting $`\varphi =\varphi \tau +\tau `$. If we let $`\eta \varphi \tau `$, then we find the equation

$$\mathrm{cos}\eta =\mathrm{sin}\eta \left[e^2F_S(\tau )\mathrm{cos}\tau F(\tau )\mathrm{sin}\tau \right]+\mathrm{cos}\eta \left[e^2F_S(\tau )\mathrm{sin}\tau +F(\tau )\mathrm{cos}\tau \right],$$
(D14)
where

$$F(\tau )e^2F_C(\tau )\widehat{𝚺}_++(1\widehat{𝚺}_+)\mathrm{cos}\tau +\frac{\widehat{𝒃}_+}{\widehat{𝒃}_{}}\widehat{𝚺}_{}[1+e^2\mathrm{sn}^2(\tau )].$$
(D15)
We see that $`F(\tau )`$ is small only if both $`e^21`$ and $`1\widehat{𝚺}_+1\widehat{𝛀}_+1`$, so the expansion implied here involves a pair of small parameters. To zeroth order the solution to Eq. (D14) is $`\eta =(2n+1/2)\pi `$; writing the solution in general as $`\eta =(2n+1/2)\pi +\delta `$ we find

$$\mathrm{tan}\delta =\frac{[e^2F_S(\tau )\mathrm{cos}\tau F(\tau )\mathrm{sin}\tau ]}{1[e^2F_S(\tau )\mathrm{sin}\tau +F(\tau )\mathrm{cos}\tau ]},$$
(D16)
or, for small values of $`\delta `$,

$$\delta [F(\tau )\mathrm{sin}\tau e^2F_S(\tau )\mathrm{cos}\tau )][1+e^2F_S(\tau )\mathrm{sin}\tau +F(\tau )\mathrm{cos}\tau +\mathrm{}].$$
(D17)
It should be apparent from Eq. (D17) that arrival times should exhibit oscillations of decreasing magnitude with frequencies $`r\omega _p`$, where $`r`$ is an integer. Rigorously, we should expand the functions $`F(\tau )`$ and $`F_S(\tau )`$ in powers of $`e^2`$ as well, and the hiden functions $`\mathrm{dn}(\tau )`$, $`\mathrm{cn}(\tau )`$ and $`\mathrm{sn}(\tau )`$ should be expanded in $`k^2`$ (which is small if $`1\widehat{𝛀}_+^2`$ is small: see Eq. \[LABEL:ktaudefpos\]).

For the triaxial case we have to consider, in addition to $`\delta `$, the mapping between $`t`$ and $`\varphi `$, which is complicated by the fact that $`d\varphi /dt`$ is time dependent. Eq. (D5) implies that, to lowest nontrivial order in $`e^2`$,

$$\frac{\overline{I}}{L}\frac{d\varphi }{dt}1+\mathrm{\Delta }(1e^2)+e^2\mathrm{\Delta }\mathrm{cos}2\tau ,$$
(D18)
which integrates to

$$\varphi \left[1+\mathrm{\Delta }(1e^2)\right]_0^t\frac{dtL}{\overline{I}}+\frac{e^2\mathrm{sin}2\tau }{4|\widehat{𝛀}_+|}.$$
(D19)
and therefore

$$\eta \left[1\mathrm{\Delta }(2|\widehat{𝛀}_+|1+e^2)\right]_0^t\frac{dt}{\overline{I}}+\frac{e^2\mathrm{sin}2\tau }{4|\widehat{𝛀}_+|}.$$
(D20)
The integral includes the effects of sinusoidal variations in pulsar spindown, which may be evaluated using the formulae in § 4.3. Thus, the final solution for arrival times is

$$\left[1\mathrm{\Delta }(2|\widehat{𝛀}_+|1+e^2)\right]_0^t\frac{dtL}{\overline{I}}\left(2n+\frac{1}{2}\right)\pi +\delta \frac{e^2\mathrm{sin}2\tau }{4|\widehat{𝛀}_+|}.$$
(D21)
We see that the variation of $`d\varphi /dt`$ with time may introduce variation in pulse arrival times at twice the precession frequency, a possibility that only arises when the moment of inertia tensor is nonaxisymmetric. Note that in the expansion in Eq. (D18) we only kept terms up to $`e^2`$. Higher order terms would introduce further variability (presumably at $`4\omega _p`$, $`6\omega _p`$ and so on).

#### D.1.2 $`I_B>0`$ and $`k^2>1`$

For $`k^2>1`$, we choose $`𝑽_{\widehat{1}}=𝑽_+`$, $`𝑽_{\widehat{2}}=𝑽_2`$ and $`𝑽_{\widehat{3}}=𝑽_{}`$. We then find that

$`\mathrm{cos}\alpha `$ $`=`$ $`(1\mathrm{\Delta })\mathrm{dn}(\stackrel{~}{\tau })\sqrt{{\displaystyle \frac{1\widehat{𝛀}_+^2}{1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2}}}`$
$`\mathrm{tan}\psi `$ $`=`$ $`\mathrm{sign}(\widehat{𝛀}_{}){\displaystyle \frac{\mathrm{cn}(\stackrel{~}{\tau })}{\mathrm{sn}(\stackrel{~}{\tau })}}\sqrt{{\displaystyle \frac{(1+\mathrm{\Delta })(\mathrm{\Delta }\mathrm{\Delta }_0)}{2\mathrm{\Delta }(1\mathrm{\Delta }_0)}}}`$
$`{\displaystyle \frac{d\varphi }{dt}}`$ $`=`$ $`{\displaystyle \frac{L}{\overline{I}}}\left[{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Delta }_0+(\mathrm{\Delta }+\mathrm{\Delta }_0)\mathrm{sn}^2(\stackrel{~}{\tau })}{(1+\mathrm{\Delta })(\mathrm{\Delta }\mathrm{\Delta }_0)+(1\mathrm{\Delta })(\mathrm{\Delta }+\mathrm{\Delta }_0)\mathrm{sn}^2(\stackrel{~}{\tau })}}\right].`$ (D22)
For this case, we define $`\widehat{𝚺}_\pm `$ as before, and

$$\stackrel{~}{e}^2\frac{(\mathrm{\Delta }+\mathrm{\Delta }_0)(1\mathrm{\Delta })}{(\mathrm{\Delta }\mathrm{\Delta }_0)(1+\mathrm{\Delta })};$$
(D23)
in terms of these variables, we find that

$$\widehat{𝒃}_y=\frac{\widehat{𝒃}_+[\mathrm{sign}(\widehat{𝛀}_{})\sqrt{1+\stackrel{~}{e}^2}\mathrm{sn}(\stackrel{~}{\tau })\mathrm{sin}\varphi +\widehat{𝚺}_{}\mathrm{dn}(\stackrel{~}{\tau })\mathrm{cos}\varphi ]}{\sqrt{1+\stackrel{~}{e}^2\mathrm{sn}^2(\stackrel{~}{\tau })}}\widehat{𝒃}_{}\widehat{𝚺}_+\mathrm{cos}\varphi \sqrt{1+e^2\mathrm{sn}^2(\stackrel{~}{\tau })},$$
(D24)
where $`\widehat{𝒃}_\pm `$ are given by Eq. (D8). Pulse arrival times are found by solving $`\widehat{𝒃}_y=0`$. The equations are exactly the same as for $`k^2<1`$ except for the replacements $`e^2\stackrel{~}{e}^2=1/e^2`$, $`\widehat{𝒃}_{}\mathrm{sign}(\widehat{𝛀}_{})\widehat{𝒃}_+`$, $`\widehat{𝒃}_+\widehat{𝒃}_{}`$. $`\widehat{𝚺}_+|\widehat{𝚺}_{}|`$ and $`\widehat{𝚺}_{}\widehat{𝚺}_+`$. Thus, if we define $`\eta =\varphi \stackrel{~}{\tau }`$, and also

$`\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })`$ $`=`$ $`\mathrm{sin}\stackrel{~}{\tau }\sqrt{1+\stackrel{~}{e}^2}\mathrm{sn}(\stackrel{~}{\tau })`$
$`\stackrel{~}{e}^2\stackrel{~}{F}_C(\stackrel{~}{\tau })`$ $`=`$ $`\mathrm{cos}\stackrel{~}{\tau }\mathrm{cn}(\stackrel{~}{\tau })\mathrm{dn}(\stackrel{~}{\tau })`$
$`\stackrel{~}{F}(\stackrel{~}{\tau })`$ $`=`$ $`\stackrel{~}{e}^2\stackrel{~}{F}_C(\stackrel{~}{\tau })|\widehat{𝚺}_{}|+(1|\widehat{𝚺}_{}|)\mathrm{cos}\tau +\mathrm{sign}(\widehat{𝛀}_{}){\displaystyle \frac{\widehat{𝒃}_{}}{\widehat{𝒃}_+}}\widehat{𝚺}_+[1+\stackrel{~}{e}^2\mathrm{sn}^2(\stackrel{~}{\tau })],`$ (D25)
then we find

$$\mathrm{cos}\eta =\mathrm{sin}\eta \left[\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })\mathrm{cos}\stackrel{~}{\tau }\stackrel{~}{F}(\stackrel{~}{\tau })\mathrm{sin}\stackrel{~}{\tau }\right]+\mathrm{cos}\eta \left[\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })\mathrm{sin}\stackrel{~}{\tau }+\stackrel{~}{F}(\stackrel{~}{\tau })\mathrm{cos}\stackrel{~}{\tau }\right].$$
(D26)
As before, the approximate solution is $`\eta =(2n+1/2)\pi +\delta `$, where

$$\mathrm{tan}\delta =\frac{[\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })\mathrm{cos}\stackrel{~}{\tau }\stackrel{~}{F}(\stackrel{~}{\tau })\mathrm{sin}\stackrel{~}{\tau }]}{1[\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })+\stackrel{~}{F}(\stackrel{~}{\tau })\mathrm{cos}\stackrel{~}{\tau }]},$$
(D27)
which becomes

$$\delta [\stackrel{~}{F}(\stackrel{~}{\tau })\mathrm{sin}\stackrel{~}{\tau }\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })\mathrm{cos}\stackrel{~}{\tau }][1+\stackrel{~}{e}^2\stackrel{~}{F}_S(\stackrel{~}{\tau })\mathrm{sin}\stackrel{~}{\tau }+\stackrel{~}{F}(\stackrel{~}{\tau })\mathrm{cos}\stackrel{~}{\tau }+\mathrm{}]$$
(D28)
for small $`\delta `$. This approximate solution is valid as long as $`|\widehat{𝛀}_{}|1`$ and $`\stackrel{~}{e}^21`$. The restriction to small $`|\widehat{𝛀}_{}|`$ may seem a bit strange, until we recall that for $`I_B>0`$ and $`\stackrel{~}{k}^2=1/k^2<1`$, the spinning star is somewhat prolate rather than oblate. Thus, simple solutions are expected for small amplitude around the axis of smallest moment of inertia in this case.

The mapping from pulse phase to pulse arrival times requires the connection between $`\eta =\varphi \stackrel{~}{\tau }`$ and $`t`$. From the last of Eqs. (D22) we find that

$$\frac{\overline{I}}{L}\frac{d\varphi }{dt}1\mathrm{\Delta }(1e^2)e^2\mathrm{\Delta }\mathrm{cos}2\stackrel{~}{\tau },$$
(D29)
which integrates to

$$\varphi \left[1\mathrm{\Delta }(1e^2)\right]\frac{dtL}{\overline{I}}\frac{e^2\mathrm{sin}2\tau }{4\sqrt{1\widehat{𝛀}_+^2}},$$
(D30)
when we use $`\omega _p2\mathrm{\Delta }L\sqrt{1\widehat{𝛀}_+^2}/\overline{I}`$. Thus we find that

$$\eta =\varphi \stackrel{~}{\tau }\left[1\mathrm{\Delta }\left(1+2\sqrt{1\widehat{𝛀}_+^2}e^2\right)\right]_0^t\frac{dtL}{\overline{I}}\frac{e^2\mathrm{sin}2\stackrel{~}{\tau }}{4\sqrt{1\widehat{𝛀}_+^2}},$$
(D31)
so that pulses arrive at times

$$\left[1\mathrm{\Delta }\left(1+2\sqrt{1\widehat{𝛀}_+^2}e^2\right)\right]_0^t\frac{dtL}{\overline{I}}\left(2n+\frac{1}{2}\right)\pi +\delta +\frac{e^2\mathrm{sin}2\stackrel{~}{\tau }}{4\sqrt{1\widehat{𝛀}_+^2}}.$$
(D32)
As before, we see that to lowest order in $`e^2`$, the mapping from phase to time introduces variation at $`2\omega _p`$. Retaining higher order terms would introduce further variability at $`4\omega _p`$, etc.

### D.2 Pulse Arrival Times for $`I_B<0`$

For $`I_B<0`$, $`𝑽_{\widehat{1}}=𝑽_2`$ and $`𝑽_{\widehat{2}}=𝑽_{}`$, and Eq. (D4) is replaced by

$`\mathrm{cos}\alpha `$ $`=`$ $`\mathrm{dn}(\tau )\sqrt{{\displaystyle \frac{[\mathrm{\Delta }_0+\mathrm{\Delta }^2+\mathrm{\Delta }(1+\mathrm{\Delta }_0)(2\widehat{𝛀}_+^21)](1+\mathrm{\Delta })}{(\mathrm{\Delta }+\mathrm{\Delta }_0)[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}`$
$`\mathrm{tan}\psi `$ $`=`$ $`{\displaystyle \frac{\mathrm{cn}(\tau )}{\mathrm{sn}(\tau )}}\sqrt{{\displaystyle \frac{2\mathrm{\Delta }(1\mathrm{\Delta }_0)}{(1\mathrm{\Delta })(\mathrm{\Delta }+\mathrm{\Delta }_0)}}},`$ (D33)
and Eq. (D5) is replaced by

$$\frac{d\varphi }{dt}=\frac{L}{\overline{I}}\left[\frac{2\mathrm{\Delta }+(\mathrm{\Delta }_0\mathrm{\Delta })\mathrm{sn}^2(\tau )}{2\mathrm{\Delta }(1\mathrm{\Delta }_0)+(\mathrm{\Delta }_0\mathrm{\Delta })(1+\mathrm{\Delta })\mathrm{sn}^2(\tau )}\right].$$
(D34)
In this case, we define

$`\widehat{𝚺}_+`$ $`=`$ $`\sqrt{{\displaystyle \frac{(1+\mathrm{\Delta })[\mathrm{\Delta }_0+\mathrm{\Delta }^2+\mathrm{\Delta }(1+\mathrm{\Delta }_0)(2\widehat{𝛀}_+^21)]}{(\mathrm{\Delta }+\mathrm{\Delta }_0)[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}`$
$`\widehat{𝚺}_{}`$ $`=`$ $`\widehat{𝛀}_{}\sqrt{{\displaystyle \frac{2\mathrm{\Delta }(1\mathrm{\Delta }_0)(1\mathrm{\Delta })}{(\mathrm{\Delta }+\mathrm{\Delta }_0)[1+2\mathrm{\Delta }(2\widehat{𝛀}_+^21)+\mathrm{\Delta }^2]}}}`$
$`e^2`$ $`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }_0\mathrm{\Delta })(1+\mathrm{\Delta })}{2\mathrm{\Delta }(1\mathrm{\Delta }_0)}};`$ (D35)
then

$$\widehat{𝒃}_y=\frac{\widehat{𝒃}_{}[\mathrm{cn}(\tau )\mathrm{sin}\varphi \widehat{𝚺}_+\sqrt{1+e^2}\mathrm{dn}(\tau )\mathrm{sn}(\tau )\mathrm{cos}\varphi ]}{\sqrt{1+e^2\mathrm{sn}^2(\tau )}}\widehat{𝒃}_+\widehat{𝚺}_{}\mathrm{cos}\varphi \sqrt{1+e^2\mathrm{sn}^2(\tau )},$$
(D36)
where $`\widehat{𝒃}_\pm `$ are given by Eq. (D8). Pulse arrival times are determined from the condition $`\widehat{𝒃}_y=0`$:

$`\mathrm{sin}(\varphi \tau )`$ $`=`$ $`[\mathrm{cos}\tau \mathrm{cn}(\tau )]\mathrm{sin}\varphi `$
$`+\left\{\widehat{𝚺}_+\sqrt{1+e^2}\mathrm{sn}(\tau )\mathrm{dn}(\tau )\mathrm{sin}\tau +{\displaystyle \frac{\widehat{𝒃}_+}{\widehat{𝒃}_{}}}\widehat{𝚺}_{}[1+e^2\mathrm{sn}^2(\tau )]\right\}\mathrm{cos}\varphi .`$
Define the functions

$`e^2G_C(\tau )`$ $`=`$ $`\mathrm{cos}\tau \mathrm{cn}(\tau )`$
$`e^2G_S(\tau )`$ $`=`$ $`\sqrt{1+e^2}\mathrm{sn}(\tau )\mathrm{dn}(\tau )\mathrm{sin}\tau `$
$`G(\tau )`$ $`=`$ $`\widehat{𝚺}_+e^2G_S(\tau )(1\widehat{𝚺}_+)\mathrm{sin}\tau +{\displaystyle \frac{\widehat{𝒃}_+}{\widehat{𝒃}_{}}}\widehat{𝚺}_{}[1+e^2\mathrm{sn}^2(\tau )];`$ (D38)
then the pulse arrival times are solutions of the equation

$$\mathrm{sin}\eta =[e^2G_C(\tau )\mathrm{cos}\tau G(\tau )\mathrm{sin}\tau ]\mathrm{sin}\eta +[e^2G_C(\tau )\mathrm{sin}\tau +G(\tau )\mathrm{cos}\tau ]\mathrm{cos}\eta ,$$
(D39)
where $`\eta \varphi \tau `$. The solution is $`\eta =2\pi n+\delta `$, where

$$\mathrm{tan}\delta =\frac{e^2G_C(\tau )\mathrm{sin}\tau +G(\tau )\mathrm{cos}\tau }{1[e^2G_C(\tau )\mathrm{cos}\tau G(\tau )\mathrm{sin}\tau },$$
(D40)
or, for small $`\delta `$,

$$\delta [e^2G_C(\tau )\mathrm{sin}\tau +G(\tau )\mathrm{cos}\tau ][1+e^2G_C(\tau )\mathrm{cos}\tau G(\tau )\mathrm{sin}\tau ].$$
(D41)
Once again, it is necesary that both $`e^21`$ and $`\sqrt{1\widehat{𝛀}_+}1`$ for the approximation to be valid.

For this case, the mapping from phase to arrival times that results from expanding Eq. (D34) leads to exactly the same result as for $`I_B>0`$ and $`k^2<1`$, namely

$$\eta \left[1\mathrm{\Delta }\left(2\widehat{𝛀}_+1+e^2\right)\right]_0^t\frac{dtL}{\overline{I}}+\frac{e^2\mathrm{sin}2\tau }{4|\widehat{𝛀}_+|}.$$
(D42)
In this case, then, pulses arrive when

$$\left[1\mathrm{\Delta }\left(2\widehat{𝛀}_+1+e^2\right)\right]_0^t\frac{dtL}{\overline{I}}2n\pi +\delta \frac{e^2\mathrm{sin}2\tau }{4|\widehat{𝛀}_+|}.$$
(D43)