Precession of the Isolated Neutron Star PSR B1828-11

Stairs, Lyne&Shemar have found that arrival time residuals from PSR B1828-11 vary periodically with a period of 500 days. This behavior can be accounted for by precession of the radiopulsar, an interpretation that is reinforced by the detection of variations in its pulse profile on the same timescale. Here, we model the period residuals from PSR B1828-11 in terms of precession of a triaxial rigid body. We include two contributions to the residuals: (i) the geometric effect, which arises because the times at which the pulsar emission beam points toward the observer varies with precession phase; (ii) the spindown contribution, which arises from any dependence of the spindown torque acting on the pulsar on the angle between its spin and magnetic axes. We use the data to probe numerous properties of the pulsar, most notably its shape, and the dependence of its spindown torque on the angle between its spin and magnetic axes, for which we assume a sum of a spin-aligned component (with a weight 1-a) and a dipolar component perpendicular to the magnetic beam axis (weight a), rather than the vacuum dipole torque (a=1). We find that a variety of shapes are consistent with the residuals, with a slight statistical preference for a prolate star. Moreover, a range of torque possibilities fit the data equally well, with no strong preference for the vacuum model. In the case of a prolate star we find evidence for an angle-dependent spindown torque. Our results show that the combination of geometrical and spin-down effects associated with precession can account for the principal features of PSR B1828-11's timing behavior, without fine tuning of the parameters.


I N T RO D U C T I O N
The pulse arrival times of neutron stars can be found very accurately, which allows for the determination of the spin period and period derivative to very high precision. Normally, the arrival-time residuals, which are calculated by subtracting the period and the period derivative (and in some cases the period second derivative), are mostly white noise. However, residuals from a small number of rotating neutron stars are found to exhibit long-term cyclical, but non-oscillatory, variations with characteristic time-scales of the order of months to years (Cordes 1993). The variability may be temporary [e.g. the Vela pulsar during its Christmas glitch (McCulloch et al. 1990)] or persistent [e.g. the accreting neutron star Her X-1 (Tannanbaum et al. 1972), the Crab pulsar (Lyne, Pritchard & Smith 1988), and the pulsars PSR 1642−03 (Blaskiewicz 1992), PSR B0959−54 (D'Alessandro & McCulloch 1997) and PSR B1828−11 E-mail: akgun@astro.cornell.edu (TA); link@physics.montana.edu (BL); ira@astro.cornell.edu (IW) (Stairs, Lyne & Shemar 2000)]. The long time-scales that characterize the observed variations would arise naturally from precession, 1 when the principal axes of a body (defined through the moments of inertia, which we will take as I 1 I 2 I 3 ) revolve periodically around the angular momentum, as viewed in an inertial frame. An ellipticity = (I 3 − I 1 )/I 1 1 would be expected to produce variations in the timing residuals of an axisymmetric body with a period P p = P / ≈ 3.2P (s)(10 8 ) −1 yr, where P (s) is the rotation period in seconds. The arrival-time variations characteristic of precession would be strictly periodic, but not sinusoidal for a triaxial rotator.
There are two physical causes for arrival-time residuals ( t) in a precessing neutron star (Cordes 1993). (i) One is directly geometrical. As the rotating star precesses, the symmetry axis of its radiation L b z ϑ χ θ Figure 1. Definition of various angles: the wobble angle (θ ) is the angle between the angular momentum ( L) and the body z-axis (which is chosen as the principal axis corresponding to the largest moment of inertia, I 3 ); the beam swing angle (ϑ) is the angle between the angular momentum and the magnetic axis (b); and χ is the polar angle of the magnetic axis in the body frame. Note that the angles may not be coplanar. For an axisymmetric body the wobble angle remains constant; for a triaxial body it varies with time (see Appendix B).
beam crosses the plane defined by the angular momentum of the star and the direction to the observer at times that vary periodically over the precession cycle. The magnitude of the variability is set by the amplitude of the precession, which is roughly the wobble angle, θ (defined as the angle between the angular momentum and the principal axis corresponding to the largest moment of inertia, Fig. 1). Typically θ 1, and the amplitude of the arrival-time residuals is t geo ∼ θ P .
(ii) In addition, the dependence of the spin-down torque acting on the pulsar on the angle between its spin and magnetic axes produces a timing residual that can be comparable to and even exceed t geo . If we assume that the spin-down torque is proportional to Ω − a( Ω ·b)b, where Ω = Ω is the angular velocity,b is the magnetic axis and the dimensionless parameter a is a measure of the angular dependence (a ≡ 1 for a spinning magnetic dipole radiating into vacuum), then the spin-down rate varies over the precession cycle as well, producing a timing residual t sd ∼ aθ P 2 p t sd ∼ a P 2 p P t sd t geo , where t sd is the spin-down time-scale for the pulsar. The dimensionless parameter sd ≡ P 2 p P t sd ≈ 3.2P 2 p (yr) P (s) t sd /10 7 yr may be large. Associated with these arrival-time residuals are period residuals ( P/P ) geo ∼ θ P /P p ≈ 3.2 × 10 −8 θ P (s)/P p (yr) and ( P/P ) sd ∼ a sd ( P/P ) geo .
The best candidate to date for truly periodic long-term variations in arrival times is PSR B1828−11 (Stairs et al. 2000(Stairs et al. , 2003. Fourier analysis of these variations reveals harmonically related periodicities at approximately 1000, 500 and 250 d (Stairs et al. 2000; Fig. 2), The shape parameter is defined in terms of the weight of the narrow (A N ) and wide (A W ) standard pulse profiles that are present at every epoch, S = A N /(A N + A W ), so that S 1 for narrow pulses, and S 0 for wide ones (Stairs et al. 2000(Stairs et al. , 2003. The solid line uses a cubic spline to connect the points for P. with the latter two somewhat more pronounced than the first. The length of the time-scale of these variations implies that they are probably not of magnetospheric origin, since the natural time-scale in the magnetosphere is of the order of the spin period, which in this case is P = 0.405 s. Even the E × B drift of subpulses (e.g. Ruderman & Sutherland 1975) does not exceed ∼10 P . [However, Ruderman (2001) has suggested the possibility of drifts with periods of the order of a year.] As of the time of writing, there are no quantitative models for the data based solely on magnetospheric effects, but there are successful models based on precession (e.g. Jones & Andersson 2001;Link & Epstein 2001;Wasserman 2003). Link & Epstein (2001) previously modelled the timing residuals from this pulsar in terms of precession of an axisymmetric, oblate rotating rigid body slowing down according to the vacuum magnetic dipole radiation formula. They found that the observations could be accounted for in this model provided that the underlying pulsar is nearly an orthogonal rotator (magnetic obliquity χ ≈ 89 • to the body's symmetry axis) and has a nearly aligned angular momentum (wobble angle θ ≈ 3 • between angular velocity and symmetry axis). These are in accordance with the conclusions reached by Jones & Andersson (2001). Although the precession amplitude is small, it may suffice to unpin superfluid vortex lines (Link & Cutler 2002), thus avoiding a potential impediment to precession: pinning was shown to shorten the precession period to about 100 spin periods, and precession itself is dissipated over a time-scale of 100-10 000 precession periods (Shaham 1977(Shaham , 1986Sedrakian, Wasserman & Cordes 1999). Wasserman (2003) argued that the data could also be accounted for if the underlying neutron star has either a type II superconductor or a strong toroidal magnetic field in its core. In these models, the angle between the angular velocity and symmetry axis of the star could be larger than found by Link & Epstein (2001), and the star does not have to be a nearly orthogonal rotator; spin-down variations were found to dominate the timing residuals in this case as well. Link (2003) showed that the standard picture of the core of type II superconducting protons coexisting with superfluid neutrons is inconsistent with long-period precession; pinning of the neutron vortices to the proton flux tubes makes the precession frequency comparable to the rotation frequency of the star, a factor of 10 8 too fast. Possible implications include a normal core (both neutrons and protons), superfluid neutrons and normal protons, normal neutrons and superfluid protons (type I or type II), or superfluid neutrons and type I protons. Sedrakian (2005) studied the last possibility. He calculated the drag on neutron vortices moving in a type I superconducting core, and found that the drag is sufficiently small (that is, the vortices are sufficiently mobile with respect to the protons) that long-period precession is indeed possible in this scenario. Irrespective of the details, magnetic stresses in excess of the relatively weak ones that would arise from the pulsar's apparent dipole field strength, together with crustal stresses, would render the neutron star effectively triaxial in shape (Cutler 2002;Wasserman 2003;Cutler, Ushomirsky & Link 2003). Link & Epstein (2001) and Wasserman (2003) gave two alternative models that interpret the timing of PSR B1828−11 as precession. These models fit the data well, thus providing strong evidence that the observed timing variations do indeed represent precession. These two models, however, are special cases. The purpose of this paper is to do a thorough search of the parameter space to see what we can learn about the properties of the spin-down torque and the stellar figure. To this end, we analyse the period residuals from this pulsar in terms of a simple model in which the rotating neutron star is assumed to be a triaxial rigid body. Obviously, precise axisymmetry is a special case, and we do not expect it to hold generally, particularly if the crust of the star is not in a relaxed state (e.g. Cutler et al. 2003), or has substantial internal magnetic stresses that may not be axisymmetric to begin with. Thus, one of our goals is to see what the data from PSR B1828−11 reveal about the shape of the neutron star crust.
In this paper, we model a precessing neutron star as a single (rigid) body, rotating uniformly. Realistic neutron star modelling should take account of at least two different components -its solid crust and (super)fluid core. Bondi & Gold (1955) considered the precession of a body consisting of a solid crust coupled frictionally to a fluid core. Their work showed that the long-term precession of the composite system depends on the time-scale t cc on which the crust and core couple to one another. If t cc 1 then the crust and core are very tightly coupled to one another on time-scales smaller than a rotation period, and the moment-of-inertia tensor relevant to precession is that of the entire system, crust plus core. In this case, precession damps out slowly, on a time-scale ∼ ( t cc ) −1 precession periods. If t cc 1 the crust and core only couple on time-scales long compared to a rotation period, and the momentof-inertia tensor relevant to precession is that of the crust alone. In this case, precession also damps slowly, with a characteristic decay time-scale ∼ t cc precession periods. Estimates of the crust-core coupling time-scale vary, but the consensus is that the coupling is weak, with t cc ∼ 10 2 -10 4 1 (e.g. Alpar, Langer & Sauls 1984;Alpar & Sauls 1988;Sedrakian & Sedrakian 1995), so the precession dynamics are governed by the moment-of-inertia tensor of the crust alone. However, it also turns out that, as long as the crust and core couple on a time-scale short compared with a precession period, but still long compared to the spin period, the relevant moment of inertia for all spin-down effects, including those that vary periodically over a precession period, is the total stellar moment of inertia (Akgun, Link & Wasserman, in preparation); this is the appropriate regime as long as t cc 10 8 , which appears to be the case. Thus, our one-component model for precession is justified, apart from slow decay of the precession, which we neglect.
Another issue is that the neutron star crust is not perfectly rigid, but has a finite shear modulus. For a biaxial precessing star, the crust must be strained in order for the star to precess with a period of the order of a year (Cutler et al. 2003). In addition, the strain field will vary with time as the star precesses, making the star's moment-ofinertia tensor time-dependent in the body frame. For simplicity, we neglect these effects, and assume that the rotation of an imperfectly rigid, triaxial star is well described by the Euler equations for a rigid body, but with a moment-of-inertia tensor that is re-scaled to account for the finite shear modulus.
Since P = 0.405 s and t sd ≈ 10 5 yr for PSR B1828−11, sd ∼ 10 3 , and the spin-down contribution to the precession-induced timing residuals is particularly important. As a result, we may also hope to use these data to probe the value of a, that is, to probe the angular dependence of the spin-down torque. While it is common to assume that a ≡ 1 for rough analysis, on theoretical grounds we should not expect this to be true, for even an aligned rotator surrounded by a magnetosphere radiates energy, a process whose source is ultimately the rotational energy of the star, resulting in spin-down at a rate presumably not much different from its luminosity divided by its rotational frequency. One of our chief findings is evidence that the external torque that spins down a pulsar does indeed possess at least some angular dependence.
Our analysis uses the same segment of the data 2 that was the basis of Link & Epstein (2001); this facilitates direct comparisons between the results of the two studies. We focus on the period residuals because we can derive an analytic formula for them in terms of elliptic functions. [An analogous formula for an oblate axisymmetric rotator has been derived previously by Bisnovatyi-Kogan, Mersov & Sheffer (1990) and Bisnovatyi-Kogan & Kahabka (1993), and has been applied to the 35-d cycle of Her X-1.] Because the underlying triaxial model involves numerous parameters, using an analytic formula speeds up the computation considerably, which is a distinct advantage. By contrast, direct analysis of the timing residuals would require numerical integration of the model equations, a distinct disadvantage. Thus, for computational convenience, we analyse the period residuals rather than the arrival-time residuals. However, a straightforward analysis of the period residuals using their tabulated uncertainties yielded very large values of χ 2 (∼a few × 10 3 ) under the assumption that the residuals are due solely to precession. This indicated to us that there is extra noise in the period residuals, either because their estimated uncertainties are too small (which we regard as unlikely to account for all the noise) or because there is a physical source of period noise that smears out the smooth contribution from precession systematically. In order to account for this extra noise simply, we multiplied the tabulated uncertainties by a (single) factor F, and then marginalized over F to obtain posterior distributions of the (more interesting) parameters of the precession model. This method (Student's t-test) represents a computational realization of 'chi-by-eye' for data whose uncertainties may only be known incompletely. Details are given in Appendix D.
Section 2 contains basic features of our model; further details may be found in Appendices A, B and C. Section 3 contains results and implications of our analysis; statistical details, including the 'chi-byeye' method mentioned above, are found in Appendix D. Section 4 is a short digression on the pulse shape of PSR B1828−11, which is seen to vary systematically with precession phase (Stairs et al. 2000). Although we do not use this information in our statistical analysis, precession samples different regions in a pulsar's radio-emitting region, and offers the possibility of mapping out its shape [as has  Weisberg, Romani & Taylor (1989), and previously by Link & Epstein (2001) for PSR B1828−11].

M O D E L S
Here we review briefly the models that we use, highlighting some of the more important parameters. The derivations for the period residuals are lengthy and are left for the Appendices; we present the geometric model in Appendix A, and the spin-down model is derived in Appendix C. In what is next, we will follow the notation used in these Appendices. The parameters of our models are listed in Table 1.

Geometric model
What we refer to as the geometric model is the effect of triaxiality alone (i.e. torque-free precession). In this case, Euler's equation for the angular momentum can be solved analytically in terms of Jacobian elliptic functions (Landau & Lifshitz 1976). As we show in Appendix A, the period residuals are then found to be of the form P ge /P ≈ p f n , where P is the rotation period at a fiducial epoch, p is a dimensionless quantity of the order of P /P p ∼ , and f n is a complicated combination of the elliptic functions. Because of the inherent form of f n , the amplitude of the residuals is not trivial to predict in general, and they can exhibit very rich behaviour.
We denote the principal moments of inertia by I i ; the associated axes serve as the basis for the rotating (body) frame. We then define the following parameters: = (I 3 − I 1 )/I 1 , which measures the deviation from sphericity; e 2 = [I 3 (I 2 − I 1 )]/[I 1 (I 3 − I 2 )], which measures the degree of triaxiality; k 2 , which is the parameter of the Jacobian elliptic functions, and depends on the angular momentum and the moments of inertia; and λ, which determines the components of the angular momentum [and would be simply λ = L 1 /L 3 for an axisymmetric star, but is slightly different in the more general case (see equation A4)]. The last three are not independent: k = eλ. Note that k 2 does not depend on , but only on e 2 and λ.
Note that for e 2 = 0 the body is oblate and axisymmetric, and λ = 0 when there is no precession. In both cases k = 0, and the Jacobian elliptic functions reduce to the regular trigonometric functions. At k = 1 they become hyperbolic functions, and the angular momentum exponentially aligns with the principal axis corresponding to the intermediate moment of inertia, I 2 . Between these two extremes, the precession of the angular momentum takes place along the intersection of the sphere defined by the conservation of angular momentum (L 2 = L i L i ), and the ellipsoid defined by the conservation of energy (E = L 2 i /2I i ). The resulting shape of the trajectories is the Binet ellipsoid (Landau & Lifshitz 1976). In the limit e 2 → ∞ the body becomes prolate axisymmetric.

Spin-down model
The rotation of an isolated neutron star is not torque-free, but slows down, resulting in a gradual increase in the rotation period. It is thought that, because of the rotating magnetic field, angular momentum is lost to radiation. In the simple model of a rotating dipole in a vacuum, the pulses are emitted at the poles of the magnetic field and the torque has the form N where Ω is the instantaneous rotation axis andb is the pulse (and magnetic) axis.
This is a very crude model. The magnetic field may have nondipolar components of considerable amplitude, which we will not consider here. The pulsar is also not in a perfect vacuum, but is surrounded by a plasma-filled magnetosphere (Goldreich & Julian 1969). The vacuum torque vanishes when the angular velocity and the magnetic axis are aligned, while the presence of a magnetosphere would require a loss of angular momentum, no matter what the orientation is. Therefore, the vacuum dipole torque should give only an incomplete description at best. We adopt a general spin-down torque of the form where we introduce an additional parameter a. Loss of angular momentum mandates that a 1. It should also be positive, or we would have an angle-dependent torque that is opposite in sign to the dipole contribution. The vacuum case is retrieved by setting a = 1. The case of a = 0 corresponds to an external torque with no angular dependence, which would not produce periodic arrival-time residuals.
The torque-modified Euler's equation can be solved approximately for small , and a second contribution to the period residuals arises due to the torque, P sd /P ≈ sd ˜ . We will refer to this as the spin-down model; and the sum of both geometric and spin-down contributions will be referred to as the full model. Here, sd = I 3 N 0 / p L 2 ∼ P p /t sd ∼ 10 −5 and is determined by the spin-down properties of the neutron star (in particular, the period derivative,Ṗ, or the characteristic age, t sd ); ˜ is another complicated function of Jacobian elliptic functions and Legendre integrals. For an axisymmetric star, where ω p is now the precession frequency, and a i are some coefficients; the axisymmetric model thus has two harmonically related components. We take the two peaks in the spectrum of the period residuals of PSR B1828−11 with periods of ∼500 and ∼250 d to be the most significant. From the form of the axisymmetric model, we thus conclude that the precession period must be ∼500 d. It can also be shown that, for an axisymmetric star, in the region of interest, the geometric contribution is quite negligible compared to the spin-down contribution for a ∼ 1 (see Appendix C).
If the 1000-d period represented the precession period, we would expect to see variations of the pulsar beamwidth at the same period. While such changes were reported by Stairs et al. (2000), subsequent more detailed analysis has not shown a 1000-d period in the beamwidth data (Parry et al. 2005). We thus assume that the precession period is ∼500 d, and attribute the 1000-d component in the timing data to something unrelated to precession, such as timing noise. We note that our model cannot provide satisfactory fits to the data if the precession period is ∼1000 d.

Constraints and statistical analysis
The orientation of the angular momentum, L is fixed in the inertial frame. This is still true even in the presence of the spin-down torque, if is sufficiently small (see Appendix C). Then the requirement that the pulse beam, which we assume to be centred along the magnetic axis,b never precesses entirely out of our line of sight means that the angle between L andb should not vary by more than the angular width of the pulse itself. We will refer to the angle between L and b as the beam swing angle and denote it by ϑ. If the angular radius of the pulse is ρ, then the above constraint can be expressed as ϑ = ϑ max − ϑ min 2ρ; in general, we will require the beam swing variation to be less than some value ϑ max . We also will define the wobble angle, θ , as the angle between the angular momentum and the body z-axis (see Fig. 1 and Appendix B).
The duty cycle allows us to estimate the angular extent of the pulse, and for PSR B1828−11 this varies between 5 • and 7 • . For a circular pulse, this implies that the beam swing angle cannot be varying by more than a few degrees. Larger variations would require a more elongated pulse shape. Yet, at this time, there is not enough evidence to elaborate more on this. In particular, polarization data might be quite useful to determine the extent of the pulse. Stairs et al. (2000) also report periodic variations in the average pulse shape. We offer a possible explanation in a following section.
Another restriction may be that PSR B1828−11 does not have an interpulse. That can further restrict the relative orientation of the angular momentum and the magnetic axis. However, owing to uncertainties in the structure of the magnetic field, it is not clear that an interpulse will necessarily appear. Therefore, we do not impose this restriction. The observer's location is an additional parameter, and can be independently fixed.
We apply the two models -geometric (a = 0) and full (a = 0)under the given constraints to PSR B1828−11, using a Bayesian approach to obtain probability distribution functions (PDFs) for individual parameters. We assume specific priors in the full multidimensional parameter space (to be discussed next), but the effective priors exhibited in the projected (marginalized) 1D posterior PDFs shown in the figures below are integrals over the constraints. Because there appears to be systematic noise in the period residuals larger than their tabulated uncertainties, we scale the latter and then marginalize over the scaling factor as detailed in Appendix D. Once the likelihood is determined, the individual PDFs are obtained through integration over the remaining parameters and normalization.
Let { p k } denote the set of n parameters. Then the likelihood, L({ p k }), and the volume of integration, V({ p k }), are functions of this set. The latter also depends on the beam swing angle constraint, which itself is a function of a subset of the parameters. The priors, g i ( p i ) are functions only of the single parameter they refer to. Then the projected 1D posterior PDF for the ith parameter can be expressed as an integral of the likelihood over the remaining parameters, over the volume defined by the constraints, Similarly, the projected 1D prior can be expressed as It is these two quantities ( f i and h i ) that are plotted in Figs 4-9. Note that, if the volume of integration had not depended on the constraints, then we would simply have h i = g i .
Within the context of our precession model, we can use the PDF to compare how well different sets of model parameters fit the data. However, we cannot assess the extent to which the data demand explanation in terms of precession, as opposed to some other, completely different, physical model. Any model for the data will lead to a PDF with local maxima at certain values of the parameters of the model, and we can assess the relative significance of these peaks to quantify the extent to which the model parameters are determined by the data. Whether or not another model that is just as well motivated physically as our precession model can fit the data better is outside the scope of our analysis. Given a competitor model -of which we are unaware -Bayesian methods could be used for making model comparisons.

R E S U LT S A N D D I S C U S S I O N
The physical parameters that determine the form of the residuals are the two angles that specify the orientation of the magnetic axis in the body frame (the polar angle χ and the azimuthal angle ϕ; see Fig. A1), any two of e 2 , λ and k 2 , and a. There is also a τ 0 (measured in precession cycles) that determines the initial phase. Thus, the total number of parameters is six. The angle χ varies between 0 and π/2, and its prior is taken to be flat over cos χ ; and φ varies between 0 and 2π, and has a flat prior. In other words, we assume that orientations of the magnetic axis are equally likely over all solid angles. Priors for τ 0 and a are flat between 0 and 1. On the other hand, e 2 and λ can have any positive values, as long as the beam swing angle is constrained and k 2 < 1; therefore, we have to introduce cut-offs in their priors. The parameter λ is related to the wobble angle, and as a result of the beam swing angle constraint it cannot be too large; we take λ 0.2 with a flat prior.
The situation is slightly more complicated for e 2 . The crust of a neutron star (which in our model is the only component since we do not consider the liquid interior) can relax only through shearing motions as it spins down and so must be triaxial (Link, Franco & Epstein 1998;Franco, Link & Epstein 2000). Adding the magnetic stresses, which result from the multipolar field near the surface, would produce a very complicated figure. It is, therefore, quite unlikely that the star is oblate axisymmetric (e 2 = 0) or prolate axisymmetric (e 2 → ∞) to a very high precision. On theoretical grounds one might expect e 2 to be close to unity. Thus, we first take e 2 2 with a flat prior (see Figs 4-6). However, we find that within this range the PDF for e 2 is not confined. Because of that, we also consider a second case where we allow for larger values of e 2 (see Figs 7-9). The spin-down model we use is derived under the assumption that e 2 is not exceedingly large (see Appendix C). Therefore, we take e 2 4000, which seems to encompass the regions of interest, without violating our assumptions. Since the volume of integration is considerably larger at large values of e 2 than at small values, taking a flat PDF over e 2 in this case would greatly suppress the importance of small e 2 . Therefore, we need to incorporate a prior that is fair for both regimes: we take a flat prior over ln (1 + e 2 ), i.e. the prior for e 2 is 1/(1 + e 2 ). For both e 2 2 and e 2 4000, we calculate PDFs for three different values of the maximum beam swing angle constraint: ϑ max = 1 • , 3 • and 5 • .
In Fig. 3 we show the data that we use, together with some sample models. The parameters for these models are listed in Table 2. The best fit that we find is a purely geometrical model, which has a very large beam swing angle that is, in fact, outside our prior range (which was relaxed for determining an unconstrained, global 'favourite' model). The axisymmetric model given here is similar to that of Link & Epstein (2001) Table 2. Table 2. List of parameters for the models shown in Fig. 3.

Parameter
Geometric Axisymmetric Full constrained to be below 5 • ; in fact, as we show in Appendix C, any axisymmetric model is assured to yield quite similar results even when we introduce the additional torque parameter a.
Because of the large number of parameters, numerical integration for the PDFs is quite time consuming. The figures presented here typically have a resolution of about 100 points per parameter or less. This means that fine structure in the PDFs may have been missed. Nevertheless, the most notable structures in the PDFs are expected to remain.
In Figs 4-9 we show the projected 1D PDFs for each parameter, computed by integrating the multidimensional PDF over all other parameters. The prior is also a function of the entire set of parameters and is not separable for all except a. We define the 1D prior for a given parameter by integrating over the rest. A comparison with the full 1D PDF illustrates the importance of the period residuals in determining the PDF. Keep in mind that both the prior and the posterior PDFs also include the beam swing angle constraint. In these figures, the dotted lines are for the prior PDFs, the dashed lines are for the geometric model alone (equation A27), and the solid lines are for the geometric model and the spin-down model (equation C36) both combined.
We now discuss some of the main characteristics and implications of our analysis.

The torque parameter a
For e 2 2, Figs 4-6 show considerable probability over the whole range of acceptable values, with a peak at low a that becomes more prominent as ϑ increases. Another lower and wider peak appears for ϑ max = 5 • at larger values of a (Fig. 6). Nevertheless, neither of the peaks is highly significant, because they do not contain most of the probability. Therefore, we conclude that the data do not constrain     a strongly, and it can be quite different from the vacuum spin-down value a ≡ 1. As larger values of e 2 are permitted, Figs 7-9 show a peak at a → 1, but with a large tail extending over most of the parameter space. With increasing values of ϑ, lower a values become likelier, but most of the probability (>90 per cent) still lies at a 0.25. Thus, the value of a is not well determined, but there is evidence for an angle-dependent torque. The parameter a is truly a measure of the angle dependence of the spin-down torque; it does not depend on the geometric effect at all.

The magnetic inclination χ
The axisymmetric model discussed by Link & Epstein (2001) requires χ to be extremely close to 90 • . As discussed in Appendix C, this is true even when we allow a = 1. For a triaxial model, we find   Figure 7. PDFs for beam swing angle variation less than 1 • . In Figs 7-9 the prior for χ is flat over cos χ; the prior for e 2 is flat over ln(1 + e 2 ); all other priors are flat. The cut-off for e 2 is at 4000, and λ is confined to be below 0.2; χ is in degrees; a, e 2 and λ are dimensionless; e 2 is plotted on a log-log scale to reveal more detail. The percentages listed above the plots are the probabilities for the full model enclosed in the corresponding ranges. The prior probabilities are given for e 2 in parentheses for comparison with the posterior probabilities.
the range of acceptable χ values to be much larger. The geometric model has a peak at small χ , which moves on to higher values of χ and broadens with increasing beam swing angle. This trend is still seen at e 2 2 when a = 0 is turned on. The spin-down produces a narrow but strong peak in the vicinity of 90 • , which becomes more pronounced as ϑ is allowed to be bigger. This peak corresponds to the axisymmetric case, and implies that it requires larger ϑ values; in fact, the model discussed by Link & Epstein has ϑ 6.4 • . For e 2 4000 (Figs 7-9), the peak at large χ remains apparent, though now it is quite broad. The inclusion of points beyond e 2 ≈ 2 seems to favour a more important spin-down contribution, and the geometric effect is further suppressed. There is also a small cusp that appears in the PDF for ϑ < 5 • , at χ very near 90 • , corresponding     Figure 9. PDFs for beam swing angle variation less than 5 • , for the same set of priors as in Fig. 7. to the axisymmetric case. Yet, this peak is quite narrow, and the vast majority of the probability lies outside of it.

The triaxiality parameter e 2
For e 2 2, the PDF looks quite similar to the prior, implying that the data do not differentiate among values of e 2 (Figs 4-6). There is a narrow sharp peak at e 2 = 0 for the full model, corresponding to the oblate axisymmetric case, but the probability enclosed within the peak is very small. Note that, at the same time, the PDF for the geometric effect alone almost vanishes, i.e. for the axisymmetric case, spin-down is essential. When we allow e 2 > 2, another peak appears in the PDF (Figs 7-9). At these values of e 2 the star is once again almost axisymmetric, except that now I 1 < I 2 I 3 , i.e. the star is prolate. Qualitatively e 2 = 1 separates oblate and prolate shapes. Then, comparing the probabilities enclosed in the two regions, e 2 < 1 and e 2 > 1, we find that the prolate case contains, by far, most of the probability, even though we have adopted a prior that somewhat disfavours large e 2 values. Note that, as the beam swing angle constraint is relaxed, the probability becomes quite evenly distributed over a wide region in e 2 . This leads to the conclusion that there are many triaxial models with a wide range of e 2 that can fit the data. In other words, the data do not discriminate among values of e 2 , especially when ϑ is relatively large.

The wobble parameter λ
For both e 2 2 and e 2 4000, the PDF is contained in a region that seems to be mostly confined by the beam swing angle. The constraint results in a dramatic cut-off at the high end of λ. Beyond that point, we cannot find a value of e 2 for which ϑ will be smaller than the constraint. Within that region, there is considerable probability distributed over the whole range of λ; both models follow the same trend. We can conclude that the oblate case favours somewhat larger values of λ, which produces a peak that is partially visible for ϑ max = 5 • (Fig. 6). The prolate case, on the other hand, favours smaller values of λ, resulting in a second peak, which appears when we allow e 2 to be large (Fig. 9), but is absent when e 2 is confined to low values (Fig. 6). serves as an implicit way of evaluating the significance of the peaks in the PDF. Relatively large amplitude implies that there are no significant peaks in the PDF, meaning that no points or regions in the allowed parameter space are favoured strongly. As ϑ max is increased, two regions acquire prominence: a very narrow ridge at large χ, which extends over a wide region in a and includes the axisymmetric case, and a smaller peak at small χ and small a. Keep in mind that this second peak is further discriminated against by the prior, which is flat over cos χ, i.e. there is a factor of sin χ that also enters the PDF. Consequently, when the onedimensional PDFs are calculated, the second peak is considerably suppressed.

S H A P E O F T H E P U L S E
The pulse profile of PSR B1828−11 alternates between two different modes, one narrow and the other broad, and has Fourier power at both 250 and 500 d (Stairs et al. 2000(Stairs et al. , 2003. Stairs et al. (2003) describe how the pulse profile is determined. During a particular observing session, 16 pulse averages may be either broad or narrow, with a shape parameter S defined to be the fraction of the mean pulse shape for that session attributed to the narrow component. The shape parameter S varies systematically between ≈0 (all wide) and ≈1 (all narrow) over the precession cycle, with a strong Fourier component at the 'first harmonic' 1/250 d of the 'fundamental' precession frequency 1/500 d (Stairs et al. 2000(Stairs et al. , 2003. Link & Epstein (2001) Figure 11. Schematic of a pulse consisting of a bright core, and surrounded by smaller fainter conal blobs. Core emission is assured to be stationary, whereas conal emission could vary as the emitting blobs circulate about the beam axis.
must have an hourglass shape in order for S to exhibit substantial variability on the 250-d time-scale. [A similar elongated shape was inferred from studies of the geodetic precession of PSR B1913+16 by Weisberg & Taylor (2002).] However, they did not address the issue of mode switching during individual observing sessions at all. Here we present an alternative viewpoint centred around modelling profile mode switching within individual observing sessions, and argue that it may be able to produce some aspects of the required harmonic structure shown in Fig. 2.
The basic geometrical picture is shown in Fig. 11. We attribute the narrow component of the pulse to core emission centred around the beam axis. The broader profile is a superposition of core and conal emission. Thus, in the parlance of Rankin (1990Rankin ( , 1993, the pulsar alternates between presenting a core single (S t ) and triple (T) pulse profile. The relatively young spin-down age of PSR B1828−11 (about 0.11 Myr) and its large value of B 12 /P 2 are consistent with this categorization. However, the apparent pulsewidth in the narrow state appears to be anomalous: Rankin (1990) finds a full width at half-maximum (FWHM) pulsewidth W core = 2.45 • P −1/2 /sin χ, where χ is the angle between a pulsar's spin and magnetic axes. For PSR B1828−11, the FWHM of the narrow state is about 2.3 • , as opposed to W core ≈ 3.85 • /sin χ from Rankin's formula. We note that the bounding relationship W core 2.45 • /P −1/2 was derived from a set of 'interpulsars' thought to be nearly orthogonal rotators, so we would have expected Rankin's formula to work especially well if χ is near 90 • , as was suggested by Link & Epstein (2001); on the other hand, the discrepancy is also smallest for χ ≈ 90 • , which may be circumstantial evidence that PSR B1828−11 really is nearly an orthogonal rotator. [Moreover, the frequency dependence of the core width is relatively weak at high frequencies, so that the fact that the Rankin (1990) formula is for 1-GHz emission, whereas the Stairs et al. (2003) observations were at 1.6 or 1.7 GHz, is not responsible for the discrepancy.] We note that there are other exceptions to Rankin's (1990) bound, but not many [see e.g. fig. 23b in Graham-Smith (2003), adapted from Gould (1994)]; given uncertainties in χ there may be other pulsars with W core > 2.45 • / √ P but also W core < 2.45 • / √ P sin χ. If pulsar core emission intensity were Gaussian, we would expect an observed intensity of the form where β is the impact parameter of the observer's line of sight relative to the beam axis, φ is the pulse phase (centred on epoch of closest passage relative to the axis), and α = χ + β is the angle between the line of sight to the observer and the stellar spin axis. Here ρ 1 and ρ 2 define the extent and the shape of the beam, which would be elliptical when they are not equal. This formula assumes that β χ, and that emission is strongly beamed along magnetic field lines, but does not presume that the emission pattern is circularly symmetric with respect to the beam axis. The two directions 1 and 2 are relative to a coordinate system in whichê 3 =b coincides with the magnetic moment of the star, which is assumed to be the beam axis,ê 2 = L ×b/| L ×b| andê 1 = (b · Lb − L)/| L ×b|. For a Gaussian beam, equation (3) shows that the core component width is independent of the impact angle β, although the peak intensity is not (e.g. Rankin 1990). However, this is a unique property of a Gaussian profile. We can well imagine that the emission cuts off sharply (even discontinuously) for sufficiently large β, in which case the core width could be narrower than normal. Because the peak core intensity would also be lower in such cases, it would be harder to detect, which may account for the rarity of exceptions to Rankin's (1990) bound.
A sharper cut-off to the core emission beam would not only allow narrower core pulse profiles, but also introduce β dependence into the width. As a simple example, suppose that the beam profile is i.e. still a self-similar function but with a sharper cut-off than a Gaussian profile. The peak intensity is at φ = 0, where u = u min = β 2 /ρ 2 1 ; the FWHM is at phases ±φ 1/2 , where where the approximation is valid for small values of κ, irrespective of κu min . The cut-off becomes important once κu min ∼ 1 i.e. for β ρ 1 / √ κ. The core width decreases with increasing β, as does the peak intensity observed from the core. As we mentioned above, we ascribe the broader pulse profile state to a superposition of core and conal components. In keeping with the schematic Fig. 11, we assume that the conal emitting pattern is patchy and, as we discuss further below, probably only stationary in the mean. Consider an individual conal emitting region (hereafter 'blob') i; we assume that it is centred at (x 1 i , x 2 i ) = ρ i (cos σ i , sin σ i ). The emission pattern of blob i may be anisotropic in a complicated fashion, with possible preferred directions not only along theê 1,2 axes, but also along and perpendicular tox i = (cos σ i , sin σ i ). The observer sees an intensity that is a function of the two variables β − ρ i cos σ i and − where h i is the height of blob i above the core emitting region. The peak value of the intensity of radiation seen from any given blob is only a function of |β − ρ i cos σ i |, though, and we assume that blob i is detectable provided that where δ i may depend on σ i . Thus, the detectability of an individual blob varies through the precession cycle, and the probability of seeing any blobs at all also varies, thus affecting the observed beamwidth.
The problem of modelling the detectability of a given blob can be quite complex. To illustrate, suppose that each blob is at ρ i = ρ, and has the same anisotropic shape. Simplify even more by assuming that the emission profiles of the blobs have a characteristic length δ r along the (radial) direction from the beam axis to its centre, and a different length δ t in the direction tangential to it. Then we can detect the blob if (β − ρ cos σ i ) 2 δ 2 t sin 2 σ i + δ 2 r cos 2 σ i .
For convenience, we define δ 2 = 1 2 δ 2 t + δ 2 r and qδ 2 = 1 2 δ 2 t − δ 2 r . Note that q may be positive or negative, and |q| 1; q = 0 for a circular beam. Since we expect that, in general, δ ρ and, except for special values, δ β, we only expect blobs within a small range σ (β) to be visible. Fig. 12 shows the result of solving equation (7) for the range of observable σ /2π as a function of impact parameter β. (The solutions were not extended beyond β/ρ = 1 + δ √ 1 − q, where σ ≡ 0.) The results exhibit complicated behaviour even in this simple model. Given N conal blobs, the probability of seeing the broader pulse profile is large when 2 × N σ /2π 1, and is small when N σ /π 1. (The factor of 2 is because the observer's line of sight crosses the cone twice.) Thus, we may expect the shape parameter S to be small for impact parameters where σ /2π is large, and vice versa; during a precession cycle, both regimes may be sampled. Fig. 13 shows an example of how the probability of seeing only the narrow pulse would vary with precession phase in this model; this example captures the main features of the observed beamwidth variations shown in Fig. 2. For constructing the figure we adopted q = −0.8, δ = 0.1ρ, and assumed a sinusoidal variation of the impact parameter with precession phase, with β 0 = 0.95ρ and β 1 = 0.15ρ assumed for graphical purposes. The 'shape parameter' is taken to be S = (1 − σ /π) N , i.e. the probability that no blobs are detected; Fig. 13 assumes N = 6. Clearly, S varies periodically but not sinusoidally in this model, and also varies substantially in half a precession cycle. This distinctive 'doubly periodic' variation is only seen if the observer's line of sight crosses near β = ρ. This is consistent with our earlier discussion of core widths if the core emission is still visible but starting to cut off at such impact angles. Presumably, the peak intensity of the core emission must also far exceed that of any conal blob for this model to be viable; there are some indications that conal emission becomes more prominent as pulsars age (e.g. Rankin 1990). Although the range of variation of S in this example is smaller than in PSR B1828−11, extensions of the model, such as different assumptions about the conal emission [e.g. an hourglassshaped cone as in Link & Epstein (2001), or a more complicated version of blob anisotropy] may possibly yield a better account of the data.
If the conal blobs were stationary in the rotating frame, then the observer would see pulse profile variations as a function of precession phase, but would not see any variations at a given precession phase. However, it is likely that the conal blobs are not at fixed positions but rather circulate around the cone in a rotating carousel (Deshpande & Rankin 1999, 2001; see Rankin & Wright 2003 for a review). In this picture, which has empirical support (Deshpande & Rankin 1999, 2001, conal emission is from beams that circulate with a frequency d = f d relative to a reference frame rotating with the star; the circulation is probably the consequence of E × B drift (e.g. Ruderman & Sutherland 1975;Gil, Melikidze & Geppert 2003;Wright 2003), and f d 0.1 is a reasonable value. By contrast, the core emission is stationary, and from a much lower altitude than the conal emission (possibly from near the polar cap). In this picture, during a given observing session the core component is always visible, but the conal component fluctuates as emitting blobs periodically pop in and out of the observer's line of sight. The probability that the observer sees conal emission at all varies systematically during the precession cycle, and is fixed during any observing session lasting a day or so (i.e. far less than the precession period).
If this model is correct at least in broad outline, then the total beam swing during a precession cycle is 1-3 • , given expected core radii (Rankin 1993). Larger beam swings could be accommodated by a more complex model for the pulse shape (e.g. the hourglass shape of Link & Epstein 2001).

C O N C L U S I O N
We find a wide range of triaxial models that may explain the period residuals of PSR B1828−11 in terms of precession, even under the stringent constraints we have imposed. We find many fits that are as good as, or better than, the axisymmetric model considered before  (Link & Epstein 2001). In general, fits improve with larger beam swing angle variations ( ϑ), but if we assume that the pulse is confined to a region a few degrees in size, we have to rule them out. Prolate and oblate axisymmetric models seem to be favoured, especially for small ϑ, but that is not sufficient to rule out other triaxial models. Both the geometric and spin-down effects contribute to the fits. Oddly, if we relax our beam swing constraint completely, the data prefer a best fit that has a = 0 (no spin-down contribution), but ϑ for that model is unreasonably large (Fig. 3 and Table 1), so it is merely an unphysical curiosity.
In the oblate axisymmetric model, spin-down is the dominant effect, but it requires parameters (in particular, a large χ value) that could be expected to produce an interpulse, which is not seen in PSR B1828−11. If we were to impose the absence of an interpulse as a constraint, some of the models we have permitted in our analysis would be excluded, particularly those at large χ. Conceivably, the magnetic field and core beam structure of PSR B1828−11 are sufficiently complex that an interpulse would be absent even at χ → 90 • . We note that our model for shape variations suggests that we are only viewing the outskirts of the core emission in the component we detect, which may enhance the probability that emission from the opposite pole is undetectable. Thus, we do not impose the absence of an interpulse as a constraint on our analysis.
Our models do not require a = 1, so substantial deviations from the vacuum spin-down formula are allowed. In fact, rather small values of a are permitted for an oblate star (e 2 < 1). However, for a prolate star (e 2 > 1) we find that larger values are favoured (a > 0.25), thus providing evidence for an angle-dependent torque. To our knowledge, our analysis provides the first evidence that pulsars are spun down by a torque that depends on the angle between the magnetic moment and the instantaneous angular velocity.
The magnetic obliquity, χ , is no longer required to be extremely close to 90 • , and we find λ (which is related to the wobble angle) to be restricted mainly by the beam swing angle. Two peaks in e 2 are prominent, corresponding to the oblate axisymmetric (e 2 = 0) and prolate nearly axisymmetric (large e 2 ) cases. These are especially evident for small beam swing angle variations. For larger beam swing variations, the data do not discriminate among values of e 2 very much.
In summary, our precession model fits the data equally well for a broad range of parameters. We cannot constrain the shape of the star, but we do find evidence for angle-dependent spin-down torque. Overall, the ability of our physically motivated model to account for the principal features of the timing data of PSR B1828−11 without special choices of the parameters reinforces the idea that the pulsar is precessing. In particular, we have shown that the data can be fitted without resorting to a nearly orthogonal rotator with a vacuumlike dipole torque as Link & Epstein (2001) did in their preliminary work. Though we cannot strongly constrain the angular dependence of the spin-down torque with the present data, the potential remains for learning more about this important aspect of the neutron star magnetosphere from future observations. The parameters of our model are not very tightly constrained. Two possible reasons are that we have only three cycles of data, and that the data have a large degree of intrinsic scatter that our simple model cannot account for, creating wide PDFs. It will be interesting to see if the parameters can be more tightly constrained as more data become available over the next decade.
Our analysis does not employ the data on the shape parameter variations, because constructing a comprehensive mathematical description would require a reliable model for the pulsar beam. Although we do not possess such an accurate model, we offer an explanation using a compound pulse structure, with core and cone components.
Precession has interesting implications for pulsar observations, which so far have not been widely discussed. One immediate and very obvious effect would be disappearing pulsars, i.e. pulsars that, as a result of precession, would at some point leave the line of sight of the observer, but, excluding other effects, would eventually come back. The time-scale of such changes could be months to years.

AC K N OW L E D G M E N T S
We thank Ingrid Stairs for providing us with the data and for valuable discussion. This research is supported in part by NSF AST-0307273 (Cornell University), NSF AST-0098728 (Montana State University) and IGPP 1222R from LANL.

A P P E N D I X A : P E R I O D R E S I D UA L S F O R A T R I A X I A L R I G I D S TA R I N T H E A B S E N C E O F E X T E R NA L TO R Q U E S : T H E G E O M E T R I C C O N T R I B U T I O N
Euler's equation for a freely precessing rigid body is and can be solved analytically in terms of Jacobian elliptic functions (see Landau & Lifshitz 1976). Using the principal axes (I 1 I 2 I 3 ) as the basis for the body (rotating) frame, we can express the components of the angular momentum unit vector L as Here, the argument of the elliptic functions is and the parameter of the elliptic functions is where we make use of the following auxiliary definitions: = I 3 − I 1 I 1 and e 2 = I 3 (I 2 − I 1 ) The minus signs that we have explicitly included in our definitions are due to our choice of the initial orientation of axes. Note that ω p is not the precession frequency, since the elliptic functions do not have a period of 2π, or more precisely, ω p is not the time derivative of the angular displacement. Instead, the precession frequency is given through where 2π is the period of the elliptic functions and can be calculated using the Legendre elliptic integral of the first kind, τ = F(φ, k), with sin φ = sn τ , The values of the parameter k 2 are unrestricted, though different regimes require careful handling. Values k 2 < 1 correspond to precession around the body z-axis; k 2 = 1 corresponds to the unstable trajectories of the angular momentum, which decay exponentially towards the intermediate axis, y; and k 2 > 1 represent precession around the x-axis (see Binet ellipsoid). For now we will confine ourselves to the first case, and the other two will be left for a later section.
We are interested in an isolated neutron star, and we want to determine the arrival time of pulses (produced along the magnetic axis) for an inertial observer. Let the inertial z-axis be along the angular momentum vector, which remains constant; and let the inertial x-axis be defined by the orientation of the observer, whom we choose to locate in the first quadrant of the inertial xz plane. Let a unit vectorb denote the orientation of the magnetic axis, and b i be the rotating frame components. Then, whenever the inertial y-component (which we choose to denote by b y ) vanishes, while the inertial x-component (b x ) is positive, we get a pulse. The two frames are related through a rotation matrix constructed from the Euler angles θ , ψ and φ (see Goldstein 1980), whence the two conditions can be expressed as b y = b 1 (cos ψ sin φ + cos θ cos φ sin ψ) + b 2 (− sin ψ sin φ + cos θ cos φ cos ψ) − b 3 sin θ cos φ = 0 (A8) and b x = b 1 (cos ψ cos φ − cos θ sin φ sin ψ) − b 2 (sin ψ cos φ + cos θ sin φ cos ψ) + b 3 sin θ sin φ > 0.  Figure A1. Orientation of the magnetic axis in the body frame. We assume that the pulses are also emitted along the same axis.
Using the solution for the angular momentum (equation A2), the Euler angles are given through cos θ = 3 dn τ, sin θ = 1 1 + e 2 sn 2 τ , cos ψ = − √ 1 + e 2 sn τ √ 1 + e 2 sn 2 τ , The last equation can be written as Equation (A8) also implies that Pulses are seen when equations (A11) and (A12) are both satisfied. In the absence of precession ( 1 = 0 and k 2 = 0) the solution of equation (A12) for φ is simply given through where η = π/2 − ϕ (Fig. A1); and we can further restrict it to lie anywhere between 0 and 2π. Note that we have implicitly included the second requirement, equation (A9), by skipping every other possible solution for φ. This is a non-trivial assumption, and would break down if the magnetic axis b happens to lie between Ω and L. However, for a pulsar these two vectors are very nearly aligned since is extremely small. Therefore, we do not need to worry about such a case. We express the general solution for φ as where ζ is confined to lie within a period of tangent (i.e. π). Using tan η = b 1 /b 2 one can show that We now have two equations for φ (equations A11 and A14), which we can combine to get the arrival times of pulses: Here ζ 0 (for τ = 0) appears as a consequence of the fact that φ 0 = η + ζ 0 . We have Note that in the absence of precession (i.e. when 1 = 0) the arrival times reduce to the form which is the solution for pure rotation. The period (between two pulses) is given as If the precession period is much longer than the pulse period (as is the case for a neutron star), we can approximate the differences by derivatives: We will find it convenient to define the expression inside the parentheses as a new function: The derivative dζ /dτ is given through (from equation A15) where, from equation (A12), To evaluate dτ n /dn we will make use of the arrival-time equation (equation A16). Taking the derivative of both sides with respect to n, we get where we have defined a new dimensionless quantity, The pulse period will be given through whence the period residuals can be found to be where P = 2πI 3 /L is the rotation period of the star, and the last approximation results from our anticipation that will be sufficiently small. Indeed, from the above definitions we get, for small k 2 , The coefficient of the function f n in equation (A27) is For PSR B1828−11 the rotation period is 405.04 ms, and the precession period is about 511 d, so that B 0 3.8 ns.

A1 The axisymmetric body
For an axisymmetric body we have e 2 = k 2 = 0. We can also set b 2 = 0, which is equivalent to introducing some initial phase shift in the definition of τ . We thus get, after some rearrangement, Let us now assume that the angle θ between the symmetry axis and the angular momentum is small, i.e. 1 ∼ θ and 1 − 3 ∼ θ 2 /2, and working to second order compute the period residuals, Also let b 3 = cos χ. Then, (Here P = 2πI 3 /L and p = I 3 ω p /L = 3 .) Note that a harmonic arises from geometrical effects.

A2 Precession around the principal axis corresponding to the smallest moment of inertia (k 2 >1)
We will now look at this case in more detail. The solution given through equation (A2) is still valid. However, it will be mathematically and computationally convenient to carry out a transformation of the Jacobian elliptic functions with k 2 > 1 into functions with parameter 1/k 2 < 1. (That, in the limit 1 = 1 and 3 = 0, we get ω p = 0 whence τ = 0 and k 2 = ∞, provides further motivation.) Then, we can write the general solution as whereτ = τ k,ê 2 = 1/e 2 and the parameter of the elliptic functions is nowk 2 = 1/k 2 . We can also define a new frequency from equation (A3), Through an appropriate redefinition of axes, the solution can be expressed in a form identical to the k 2 < 1 case, except that nowÎ 1 >Î 2 >Î 3 . Define a new right-handed coordinate basis for the rotating frame, Letˆ 1 = 3 andˆ 3 = 1 . Then the components of the angular momentum can be expressed as The precession is now clockwise, as can also be verified from Euler's equation. Theˆ i have exactly the same form as before, in terms of the new moments of inertia, So doω p ,ê 2 andk 2 , as can be verified from the equations above. We have thus transformed the problem from a 'k 2 > 1 case for an I 3 > I 1 body' into a 'k 2 < 1 case for anÎ 1 >Î 3 body', which should not be surprising. The equations for the Euler angles (equations A10) remain of the same form, with the exception of cosψ. This is effectively a sign change, τ → −τ , in the argument, One must be careful with equation (A11) as well, where there is also a sign change due to the fact that nowÎ 1 >Î 3 , These two effects add up to modify the function f n defined through equation (A22),

A P P E N D I X B : T H E W O B B L E A N D B E A M S W I N G A N G L E S
We will define the wobble (θ) and beam swing (ϑ) angles according to (see Fig. 1) cos ϑ = L ·b = −b 1 1 cn τ − b 2 2 sn τ + b 3 3 dn τ.
By definition, the wobble angle is equivalent to Euler's angle θ and is constant for an axisymmetric star. The beam swing angle is related to the angle between the beam and the observer. It could exceed 90 • , but since the pulse will have a limited angular size, there is a restriction on how much it can vary throughout a precession period. Otherwise, the beam will leave the observer's line of sight. Therefore, the span of the beam swing angle serves as a constraint. The angle can be further restricted by imposing the conditions for an interpulse.
We define the widest span as ϑ = ϑ max − ϑ min . Then the constraint is that this be smaller than some value ϑ max , which is estimated based on information about the width and shape of the pulse. Note that the beam angle depends on four parameters: the two angles determining the orientation of the pulse, and any two of k 2 , e 2 and λ = 1 / 3 .

A P P E N D I X C : P E R I O D R E S I D UA L S F O R T H E S P I N -D OW N TO R Q U E
When external torques are present, Euler's equation becomes Taking the dot product with the angular momentum, we get the equation governing the evolution of its magnitude: If we now substitute L = L L in Euler's equation, we get, after some rearrangement, which governs the evolution of the orientation of the angular momentum. These two are the basic equations that need to be solved. Of course, only three (of the total of four) components are independent equations. In the classical rotating magnetic dipole model of pulsars, the angular momentum is lost to radiation. The electromagnetic torque for a spherical rigid star in vacuum is (Davis & Goldstein 1970) Note that the torque vanishes when Ω andb are aligned. However, the pulsar is not in a perfect vacuum, and is surrounded by a magnetosphere. Therefore, there should be loss of angular momentum even when these two vectors are aligned. We will therefore adopt a general spin-down torque of the form where a is a dimensionless parameter that measures the relative importance of the two components. The amplitude of the spin-down torque can be estimated from observed values of the spin-down time, and is small. We will be interested in a particular example (PSR B1828−11) where the spin-down time is Compare this with the observed precession period for the same pulsar, The ratio of the two gives The second term in equation (C3) has a magnitude of L 2 /I , and therefore the right-hand side (RHS) of that equation is quite negligible for the case of interest (as will be discussed below). The above form of the torque is true for spherical stars. This is nevertheless a good approximation, given how small N 0 and are. In fact, we will neglect all combinations of N 0 with . This is equivalent to taking Ω L within all torque terms, since the angle between these two vectors is of the order of . We therefore have Loss of energy (or angular momentum, given through the second equation) now clearly requires that a 1. Finally, we will also neglect any time dependence within N 0 itself. The third equation demands careful thought. In component form, where L i are the components of L, s = (I 2 − I 1 )/(I 3 − I 1 ), cos ϑ = L ·b, and we have already neglected second-order terms in . As long as the star is sufficiently non-spherical and the angular momentum is sufficiently misaligned with the body z-axis, we can neglect the RHS, as it causes changes in the orientation of the angular momentum smaller (by many orders of magnitude) than the second term. In other words, and 1 are small but not zero. (Keep in mind that there is no precession if either one is zero.) Also, s cannot be too close to unity, i.e. e 2 cannot be exceptionally large.
The same cannot be done on the RHS of the equation for the magnitude of the angular momentum, as it is the only term we have. Incidentally, setting N 0 = 0 would take us back to the torque-free precession case.
In order to write the equations in dimensionless form, let us divide all sides by a frequency ω p , defined in accordance with equation (A3), but where the magnitude of the angular momentum is no longer constant. Also define which, for a constant ω p , reduces to the familiar form of the torque-free case. Now, the differential equations become, after some rearrangement, Since L/ω p is time-independent, the second equation has exactly the same solution as before, except that τ is now different, and given through a differential equation on its own. In other words, the L i remain of the same form. Thus, we only need to solve equations (C14) and (C15). We define a new dimensionless constant, and write where ω p0 = p L 0 /I 3 . The differential equations now become It is worth noting that we make no assumptions in these substitutions.
If we finally define one more dimensionless constant, and let = sd˜ and δ = sdδ , then the two equations can be written as For the pulsar that we discuss here, we have sd ∼ This means that one may safely ignore the denominators of the two equations, thus further simplifying the results: Note that L · n = 1 − a( L ·b) 2 0, thus assuring that L is monotonically decreasing (as required by loss of angular momentum

C1 Arrival-time residuals
Since the L i remain of the same form, with the only difference being that τ is now determined through a differential equation (equation C14), the Euler angles remain the same (equations A10). However, owing to the time dependence of L, it is more convenient to express φ as a function of τ , and we need to replace equation (A11) by Thus, all we have to do is to replace Lt n /I 3 by τ n / p on the left-hand side of equation (A16). The period (given through equation A26) remains the same as well, as does the calculation of dτ n /dn. In fact, we run into trouble only with the period residuals, since the magnitude of the angular momentum is now changing. Define P res = P n − 2πI 3 L 0 and P n = P n − 2πI 3 L ; (C25) P n is formally the same as the torque-free case. However, observations give us only information about P res . In practice, one first determines the period (P ) at some epoch (t 0 ), and then finds the period derivative (Ṗ ) which is the secular term attributed to spin-down, and subtracts the two contributions, so that the residuals are then given through Consider the difference between the two definitions above, where we make use of equations (C19) and (C21). We thus get where P = 2π/ = 2πI 3 /L 0 . The first term is the geometric effect ( P ge /P ≈ p f n ) and the second term is the spin-down term ( P sd /P ≈ sd˜ ). There are still secular terms present in the spin-down term that need to be subtracted. This will be taken care of below. The relative amplitude of these two terms cannot be simply determined, and it is possible that either one is dominant, or that they are comparable.
Carrying out the integrals of the Jacobian elliptic functions, and substituting the values of k 2 and 2 , we get, after some rearrangement, where we have introduced the complementary parameter k 2 1 = 1 − k 2 . In equation (C31), E(am τ , k) is the Legendre elliptic integral of the second kind, and am τ is the Jacobi amplitude, am τ = sin −1 sn τ . Note that we have implicitly assumed that, at the zero of time, the angular momentum is in the xz plane of the body frame. This is a non-trivial assumption, and in general one does not have the freedom of randomly setting the initial orientation of the angular momentum. Therefore, in general we would have L = L(τ − τ 0 ), where τ 0 is the phase offset and is an additional parameter, and the overall result would be to replace˜ (τ ) above by˜ (τ − τ 0 ) −˜ (−τ 0 ), where˜ (−τ 0 ) is just a constant. For simplicity, we will continue to assume τ 0 = 0 in the rest of our derivations, but the general case should be kept in mind.
The oscillatory part of˜ , after secular terms have been removed, is given through The average is carried over a period, 2π = 4F(π/2, k), where, using equation (C31), we get ( L ·b) 2 = 1 2π 2π 0 ( L ·b) 2 dτ = 2 3 e 2 b 2 2 (1 + e 2 ) − b 2 1 k 2 1 + 2 3 e 2 b 2 1 − b 2 2 (1 + e 2 ) + b 2 3 e 2 E(π/2, k) F(π/2, k) , and we have made use of the relations am(2π) = 2π and E(2π, k) = 4E(π/2, k). We thus get which can now be used in equation (C28) to calculate the arrival-time residuals, We will find it convenient to express this equation in the following form: ˜ a = c 1 (1 − cn τ ) + c 2 sn τ + c 3 k 2 (1 − dn τ ) + c 4 k 2 E(π/2, k) F(π/2, k) τ − E(am τ, k) , where c 1 = 2 2 3 b 2 b 3 , c 2 = 2 1 3 b 1 b 3 , c 3 = −2 1 2 b 1 b 2 and c 4 = 2 1 b 2 1 − b 2 2 (1 + e 2 ) + b 2 3 e 2 . These coefficients are related to each other through It is also interesting to note that ˜ has a non-zero average over a precession period. The residuals may have non-zero average depending on when and how the period and its derivatives are calculated. This becomes particularly important when calculating the arrival-time residuals, which can be obtained by integrating the period residuals; and if the period residuals have a constant term, then the arrival-time residuals will have a linear term. Therefore, in calculating the arrival-time residuals, one will have to subtract any constant terms from the period residuals.

C2 Amplitude of the residuals
Consider the period derivative, which is given through the secular terms in˜ , where η = L · n = 1 − ac 0 and c 0 = cos 2 ϑ . The amplitude of the period residuals thus becomes, from equation (C36), and t sd = P /2Ṗ is the spin-down time. Recall that a measures the strength of the oscillating part of the spin-down torque, and must be 1. For PSR B1828−11 the period is 405.04 ms, the precession period is about 511 d, and the spin-down time is 0.11 Myr, so that we get A 0 409.95 ns.

C3 The axisymmetric body
For an axisymmetric star e 2 = k 2 = 0, but λ = k/e = 0. Owing to the symmetry we can set b 2 = 0 by shifting the zero of time through some phase τ 0 . (Note that the same cannot be done in a triaxial body, where we chose to fix the axes according to the principal moments of inertia.) In this case equation (C37) reduces to the form ˜ a = λ 1 + λ 2 sin 2χ sin(τ − τ 0 ) − λ 4 sin 2 χ sin 2(τ − τ 0 ) .
Note that the non-linearity of the dipole contribution of the torque naturally brings in a harmonic. The period residuals are then Here τ c is the characteristic time, It turns out that, for the axisymmetric case, the geometric term is quite negligible compared to the spin-down term, for the range of physical parameters of interest (Jones & Andersson 2001;Link & Epstein 2001). To convert our result for period residuals ( P/P 0 ) into residuals of the derivative of the angular velocity ( ˙ / 0 ) given by Link & Epstein (2001), we make use of which indeed gives the correct results, together with the initial phase difference of π between the two definitions, ˙ 0 = aλ 2τ c (1 + λ 2 ) − sin 2χ cos(τ − τ 0 ) + λ 2 sin 2 χ cos 2(τ − τ 0 ) .
There is one difference between the two derivations and that is the presence of the coefficient a, which measures the strength of the spin-down torque. With the addition of this new element, the number of unknowns increases to three (a, λ and χ ), while a fit to data will C 2005 The Authors. Journal compilation C 2005 RAS, MNRAS 365, 653-672 yield only two coefficients (a 1 and a 2 ; τ 0 does not contain any further information). This implies that there is a certain level of freedom in the choice of the physical coefficients. Let us denote the fitting function by f , f = a 1 sin(τ − τ 0 ) − a 2 sin 2(τ − τ 0 ). (C46) Then, the relations between the coefficients of this function and the physical parameters that we actually seek would be a 1 = aλ sin 2χ 1 + λ 2 and a 2 = aλ 2 sin 2 χ 4(1 + λ 2 ) .
Define λ = tan θ , and the ratio of the two coefficients gives tan χ tan θ = 8a 2 a 1 .
It is also possible to express χ and θ as functions of a. However, as it turns out, the range of the physical parameters is severely restricted by the beam swing angle constraint, which in the axisymmetric case is given through This forces one of the two angles to be small (which will have tan 1 < 0.09 even if we let ϑ max = 10 • ); while equation (C48) ensures that the other remains very close to 90 • . [The ratio of the coefficients is found to be a 2 /a 1 ∼ 0.4 for the data used by Link & Epstein (2001). This yields the condition tan 2 > 36, i.e. the second angle has to be larger than 88 • , in accordance with previous findings.]

A P P E N D I X D : S TAT I S T I C A L I N F E R E N C E
Denote the set of parameters by x. Then the PDF for the parameters can be calculated as, by Bayes' theorem, where D stands for data, M stands for the model and also takes into account any other information that is available on the problem apart from data (in this case the beam swing angle, which we impose as a restriction on the parameter space). P(x | M) is the prior probability for the parameters; P(D | x, M) is the likelihood; and P(D | M) is effectively a normalization constant. To find the PDF for a certain parameter, or a subset of parameters, we integrate equation (D1) over the remaining parameters. For that, we need to know the likelihood. For each data point y i at time t i , we have a theoretical prediction f i = f (t i | x, M). For well-known uncertainties with Gaussian distribution, we would then have If we assume that the error bars are not well determined and rescale them through some number F, the above equation becomes and we regard F as an additional parameter. We have to introduce a prior for F. Since we do not want it to depend much on the end-points, we take it to be flat over d lnF, i.e. proportional to dF/F, and integrate over all values of F. Define whence we get, for d data points, where we have dropped anything that does not depend on the remaining parameters x, including integrals that give constants, products of the original sigmas, and factors of √ 2π. Thus, our final result is where the first term is the prior probability, and the constant of proportionality can be computed from the condition that the final PDF is normalized to unity. This paper has been typeset from a T E X/L A T E X file prepared by the author.