format: text

text:
# Correlation between Gamma-Ray bursts and Gravitational Waves

## I Introduction

Thirty years after their discovery, Gamma Ray Bursts (GRBs) are still the most mysterious objects in the Universe, as their properties have not yet been associated to any well-known object, while for instance pulsars have been almost immediately identified as final objects of the evolution of stars and quasars have been well framed among galaxy nuclei in particular in the Seyfert class. Unfortunately “standard” astrophysical objects do not exist for GRBs. The “non standard” compact astrophysical objects which might reproduce the experimental characteristic of GRBs involve black-holes (BH), neutron stars (NS) and massive stars.

The wide range of characteristics prevents the systematic classification of GRBs. Their energy spectrum is continuous and non-thermal, and covers the range between $`1`$ keV and $`1`$ MeV; their emission lasts from $`10`$ ms to $`10^3`$ s. The Burst And Transient Source Experiment (BATSE) has detected GRBs with an event rate of about one per day between 1991 and 2000.

One of the most important results of observations is that many GRBs are at cosmological distances. This fact, recently confirmed also by the satellite BeppoSAX, which has detected some optical counterpart with redshifts close to $`z=1`$, implies the emission of a huge amount of electromagnetic radiation: in few seconds the GRBs release an energy of $`10^{51}10^{54}`$ erg. On the other hand, this implies that the GRBs are rare events in a galaxy: in fact, as the rate of GRB events is of the order $`1/\text{day}`$, their rate amounts to $`1`$ event$`/10^6`$ years per galaxy.

The theoretical framework which gives a coherent explanation of the experimental data set of GRBs is the so-called “*fireball*” model: an “inner engine” produces a flux of relativistic energy of the order $`10^{52}`$ erg and it is extremely compact and clearly not visible; this relativistic flux of particles is a kind of “fireball” able to produce electromagnetic radiation in a optically thin shell of matter. The inner engine works for 1-3 days, producing the afterglow, but the bulk of the explosion lasts typically $`10`$ s (with peaks around $`0.5`$ and $`30`$ s and a variation of six orders of magnitude, $`10^3÷10^3`$ s) when the major part of the energy is emitted.

The inner engine progenitors are likely to be coalescing binary systems, such as NS-BH or NS-NS systems, or single stars that collapse into BH in supernova-like events (the so-called “*collapsar*” models ). If this is the correct explanation of the GRBs, then we would have compact astrophysical systems for which gravity plays an important role, and gravitational waves would be emitted.

This scenario has motivated us to investigate the behavior of the gravitational wave detector AURIGA during time spans which include the arrival time of a GRB burst. The AURIGA detector  is an Al$`5056`$ resonant bar of about $`2.3`$ *tons* with a typical noise temperature of $`7`$ mK and a bandwidth of about $`1`$ Hz. The detector is sensitive to gravitational waves (g.w.) signals over 1 Hz bandwidth around each one of its two resonant frequencies i.e. 913 and 931 Hz. The sensitivity of a g.w. detector can be conveniently expressed by the quantity $`h_{min}`$, representing the minimum g.w. amplitude detectable at $`SNR=1`$, which for the AURIGA detector is $`h_{min}2÷5\times 10^{19}`$. The AURIGA sensitivity is enough to detect NS-NS, BH-NS and BH-BH mergers which take place within the Galaxy.

Previous work on GRB-GWB association appeared in the literature, specifically about a single GRB trigger , and about a set of GRB triggers ; here we present an experimental upper limit on such an association obtained with the AURIGA-BATSE data.

The plan of the paper is as follows: in Sec. II we discuss the two methods, coincidence and statistical searches, which we use for the association of the GRBs in the BATSE catalog with GWBs. The results obtained with the two methods are presented in Sec. III. In Sec. IV we analize our results and discuss the possibilities opened by our methods for future searches.

## II Search methods

The aim of our analysis is the search of an association between GRB and GWB arrival times within a time window $`W`$. The analysis has been carried out by means of two different procedures which have been already discussed in the literature: $`i`$) the correlation method which is based on the coincidence between the arrival time of candidate gravitational events selected over a given threshold and the GRB trigger and $`ii`$) a statistical method  which relies on a hypothesis test on a statistical variable representing the mean energy of the gravitational detector at the GRB trigger.

The coincidence window plays an important role in both the analysis. In fact the coincidence window should be wide enough to hold the delay distribution between GRB and GWB. We notice that its “optimal” value depends on the astrophysical sources and on the g.w. detector properties. If we assume that the GRBs are generated by internal shocks in the fireball, the delay between GRB and GWB is less than $`1`$ second  but there are still some uncertainties. Moreover, the delay is widen by the cosmological redshift. To be as much conservative as possible on the distribution of the delays we choose $`W=5`$ s. This value turns out to be also consistent with the requirements of the filtering procedures of the AURIGA detector. In fact, the search for g.w. bursts requires a filtering of the data by a Wiener-Kolmogorov filter matched to $`\delta `$-like signals  and its characteristic time is the inverse of the detector bandwidth i.e. $`1`$ s. This time, as we shall see below, establishes the timescale of the noise correlation of the filtered data.

In what follows, “$`\delta `$-like signal” means any g.w. signal which shows a nearly flat Fourier transform at the resonance frequencies of the detector ($`913`$ and $`931`$ Hz) over a $`1`$ Hz bandwidth. Therefore the metric perturbations $`h(t)`$ sensed by the AURIGA detector can be a large class of short signals (of millisecond duration) including the latest stable orbits of inspiralling NS-NS or NS-BH, the subsequent merging and the final ringdown  and the collapsar bursts  which could be expected signals associated with GRBs.

### A Coincidence search

The coincidence method has been successfully applied to the search of coincident excitations of different g.w. detectors ; much work has been devoted to exploit its potentialities and to develop robust estimates of the background of accidental, even in the presence of non stationary event rates . A g.w. candidate event is a local maximum of the filtered data corresponding to any excitation of the detector. From the AURIGA filtered output we extract a list of candidate events setting an adaptive threshold in their signal-to-noise ratio ($`SNR=5`$). The event lists used in this analysis belong to periods of satisfactory performance of AURIGA as described in Ref. . It is worth noticing that the AURIGA event search checks each event against the expected signal template by means of a $`\chi ^2`$ test . The event lists contain the information needed to describe a $`\delta `$-like signal namely, its time of arrival (in UTC units), the amplitude of the Fourier transform, and the detector noise level at that time. The BATSE data are taken from the $`4\mathrm{B}`$ catalog by Meegan et al. available on Internet and includes trigger time (in UTC units), right ascension and declination, error box and other information about every GRB triggered.

We label with $`t_A^{(i)}`$ the estimate arrival time of the $`i`$-th candidate g.w. event detected by AURIGA and with $`t_\gamma ^{(k)}`$ the trigger time of the $`k`$-th GRB detected by BATSE. A coincidence between the $`i`$-th AURIGA event and the $`k`$-th GRB trigger time is observed if $`|t_A^{(i)}t_\gamma ^{(k)}|W`$, where $`W`$ is the coincidence window.

The coincidences found have to be compared with the accidental coincidence background due to chance. Two standard methods to evaluate the probability of accidentals have been applied to the pair BATSE-AURIGA: $`i`$) performing thousands time-shifts of the arrival time of one detector with respect to the other and looking for accidental coincidences at each shift; $`ii`$) assuming independent Poisson distribution of event times and using the mean measured rates to estimate the accidental rate $`n_a`$

$$n_a=N_AN_\gamma \frac{\mathrm{\Delta }T}{T},$$
(1)
where $`N_A`$ is the number of AURIGA events in the period $`T`$, $`N_\gamma `$ the number of GRB events in the same period and $`\mathrm{\Delta }T=2W`$ is the total amplitude of the coincidence window (see Appendix A for a proof of this basic relation). To avoid the problem of multiple coincidence within the same window we prefer to estimate the background with $`W=1`$ s and to scale our results with the help of Eq. (1) to $`W=5`$ s. The consistency of the two estimates ensures that the arrival times of the AURIGA-BATSE pair are distributed as Poisson random points even if the AURIGA event rate has been found to be not stationary.

### B Statistical search

The method of the statistical search has been proposed for a pair of interferometric detectors by L. S. Finn and coworkers ; here we slightly modify their approach to the case of a single resonant detector. The data we use are the AURIGA filtered data obtained by means of the Wiener filter matched to a $`\delta `$-like signal, without setting any threshold on their amplitudes. If a signal enters the detector at time $`t_0`$, its output can be written as

$$y(t)=hf(tt_0)+\eta (t),$$
(2)
where $`f(tt_0)`$ is the normalized signal template of the Wiener filter, $`h`$ its amplitude (i.e. $`f(0)=1`$), $`t_0`$ its arrival time and $`\eta (t)`$ is a stochastic process with zero mean and correlation

$$\eta (t),\eta (t^{})=\sigma _h^2f(tt^{}).$$
(3)
Here $`\sigma _h^2`$ is the variance of the filter output in the absence of any signal and $`f(t)`$ is a superposition of two exponentially damped oscillating functions  which can be approximately expressed as

$$f(t)e^{|t|/t_w}\mathrm{cos}(\omega _0t)\mathrm{cos}(\omega _Bt)(\omega _0t_w1),$$
(4)
where $`t_w`$ is the Wiener filter decaying time (i.e. the inverse of the detector bandwidth), $`\omega _0`$ a center carrier frequency and $`\omega _B`$ is an amplitude modulation frequency. Typical values for the AURIGA detector are $`t_w1`$ s, $`\omega _0920`$ Hz and $`\omega _B20`$ Hz.

Let us define now the random variable $`X`$, which represents an averaged measure of the energy released by a g.w. signal impinging on the detector at time $`t_0`$

$$X(t_0)=\frac{1}{2W}_W^W𝑑t|y(tt_0)|^2,$$
(5)
where $`t_0`$ is the center of the time window. If an association between GRB and gravitational waves exists, the filtered output of the gravitational waves detector, in periods just prior the GRB (“on-source” population) will differ (statistically) from the output at other times (“off-source” population).

A statistically significant difference between on- and off-source populations clearly supports a GWB-GRB association. For each of the $`N_{\mathrm{𝑜𝑛}}`$ GRB trigger we compute $`X(t_\gamma ^{(k)})`$, $`k=1,\mathrm{},N_{\mathrm{𝑜𝑛}}`$ which forms the $`\chi _{\mathrm{𝑜𝑛}}`$ set of on-source events, and construct a complementary set $`\chi _{\mathrm{𝑜𝑓𝑓}}`$ with $`N_{\mathrm{𝑜𝑓𝑓}}`$ off-source events using windows before and after the trigger. The sets $`\chi _{\mathrm{𝑜𝑛}}`$ and $`\chi _{\mathrm{𝑜𝑓𝑓}}`$ are samples drawn from the populations whose distributions we denote $`p_{\mathrm{𝑜𝑛}}`$ and $`p_{\mathrm{𝑜𝑓𝑓}}`$. The off-source events are taken in periods not correlated with the trigger time, at a distance in time greater than $`10^3`$ seconds from GRB the trigger, both before and after it; this should be sufficient to have a fair sample of the off-source events as any GRB-GWB association is reasonably excluded.

For windows $`W`$ greater than the Wiener filter characteristic time, the central limit theorem implies that $`p_{\mathrm{𝑜𝑓𝑓}}`$ is a normal distribution.

Now suppose that the GWBs fall within the window $`W`$ opened around the GRB trigger time and that the SNR of the gravitational signal associated with the GRB (averaged over the source population) is smaller than one, then $`p_{\mathrm{𝑜𝑛}}`$ is also a normal distribution with mean

$`\mu _{\mathrm{𝑜𝑛}}`$ $`=`$ $`\mu _{\mathrm{𝑜𝑓𝑓}}+\mathrm{E}\left[{\displaystyle \frac{1}{2W}}{\displaystyle _W^W}𝑑t|hf(tt_\gamma )|^2\right]`$ (6)
$``$ $`\mu _{\mathrm{𝑜𝑓𝑓}}+\left({\displaystyle \frac{t_w}{8W}}\right)\mathrm{E}[h^2](Wt_w),`$ (7)
where $`\mathrm{E}[]`$ is the average over the astrophysical source population of GRBs.

The basic idea behind the statistical approach is a hypothesis testing where the null hypothesis $`_0`$ to test is the equivalence of the off-source and on-source distributions:

$$_0:p_{\mathrm{𝑜𝑓𝑓}}(X)=p_{\mathrm{𝑜𝑛}}(X).$$
(8)
The rejection of $`_0`$ clearly supports a GWB-GRB association. Since $`p_{\mathrm{𝑜𝑛}}`$ and $`p_{\mathrm{𝑜𝑓𝑓}}`$ are normal and could differ only in their mean values, we can test $`_0`$ by the Student’s $`t`$-test.

The $`t`$ statistic is defined from $`\chi _{\mathrm{𝑜𝑛}}`$ and $`\chi _{\mathrm{𝑜𝑓𝑓}}`$ by

$`t`$ $`=`$ $`{\displaystyle \frac{\widehat{\mu }_{\mathrm{𝑜𝑛}}\widehat{\mu }_{\mathrm{𝑜𝑓𝑓}}}{\mathrm{\Sigma }}}\sqrt{{\displaystyle \frac{N_{\mathrm{𝑜𝑛}}N_{\mathrm{𝑜𝑓𝑓}}}{N_{\mathrm{𝑜𝑛}}+N_{\mathrm{𝑜𝑓𝑓}}}}}`$ (9)
$`\mathrm{\Sigma }^2`$ $`=`$ $`{\displaystyle \frac{(N_{\mathrm{𝑜𝑛}}1)\widehat{\sigma }_{\mathrm{𝑜𝑛}}^2+(N_{\mathrm{𝑜𝑓𝑓}}1)\widehat{\sigma }_{\mathrm{𝑜𝑓𝑓}}^2}{N_{\mathrm{𝑜𝑛}}+N_{\mathrm{𝑜𝑓𝑓}}2}},`$ (10)
where $`\widehat{\mu }_{\mathrm{𝑜𝑛}}`$ and $`\widehat{\mu }_{\mathrm{𝑜𝑓𝑓}}`$ ($`\widehat{\sigma }_{\mathrm{𝑜𝑛}}^2`$ and $`\widehat{\sigma }_{\mathrm{𝑜𝑓𝑓}}^2`$) are the sample means (variances) of $`\chi _{\mathrm{𝑜𝑛}}`$ and $`\chi _{\mathrm{𝑜𝑓𝑓}}`$, respectively.

The expected value of $`t`$ averaged on the source population and on the filtered output of the detector is

$$\mu _t=\mathrm{E}[t]=\left(\frac{t_w}{8W}\right)\frac{\mathrm{E}[h^2]}{\sigma }\sqrt{\frac{N_{\mathrm{𝑜𝑛}}N_{\mathrm{𝑜𝑓𝑓}}}{N_{\mathrm{𝑜𝑛}}+N_{\mathrm{𝑜𝑓𝑓}}}},$$
(11)
where $`\sigma =\mathrm{E}[\mathrm{\Sigma }]`$.

Let us consider the upper limit on $`\mathrm{E}[h^2]`$ in the assumption $`_0`$ is true: in this case the most probable value of $`\mathrm{E}[h^2]`$ is zero. From Eq. 11 we get

$`\left({\displaystyle \frac{t_w}{8W}}\right){\displaystyle \frac{\mathrm{E}[h^2]}{\sigma }}\mu _{t,max}\sqrt{{\displaystyle \frac{N_{\mathrm{𝑜𝑛}}+N_{\mathrm{𝑜𝑓𝑓}}}{N_{\mathrm{𝑜𝑛}}N_{\mathrm{𝑜𝑓𝑓}}}}}=`$ (12)
$`=\{\begin{array}{cc}\mu _{t,max}\sqrt{2/N_\gamma }\hfill & (N_{\mathrm{𝑜𝑛}}=N_{\mathrm{𝑜𝑓𝑓}}=N_\gamma )\hfill \\ \mu _{t,max}/\sqrt{N_{\mathrm{𝑜𝑛}}}\hfill & (N_{\mathrm{𝑜𝑓𝑓}}N_{\mathrm{𝑜𝑛}})\hfill \end{array}`$ (15)
where $`\mu _{t,max}`$ is the upper limit. Assuming $`_0`$ true and $`N_{\mathrm{𝑜𝑓𝑓}}N_{\mathrm{𝑜𝑛}}`$, we have

$$\mathrm{E}[h^2]h_{max}^2=\left(\frac{8W}{t_w}\right)\frac{\mu _{t,max}}{\sqrt{N_{\mathrm{𝑜𝑛}}}}\sigma ;$$
(16)
the value of $`\mu _{t,max}`$ can be deduced by the selected confidence level. In this way we are able to set an upper limit on the average amplitude of gravitational signals associated with GRBs

$`h_{\mathrm{RMS}}^2`$ $``$ $`[1.4\times 10^{18}]^2{\displaystyle \frac{W}{5s}}\left({\displaystyle \frac{t_w}{1s}}\right)^1{\displaystyle \frac{\mu _{t,max}}{1.96}}\times `$ (17)
$`\times `$ $`\left({\displaystyle \frac{N_{\mathrm{𝑜𝑛}}}{100}}\right)^{1/2}{\displaystyle \frac{\sigma }{[5\times 10^{19}]^2}}.`$ (18)
## III Results

The AURIGA data used in the two analysis are relative to the years 1997 and 1998, and the number of GRBs which fall into the AURIGA data taking periods is 120.

### A Coincidence search

Within the window of 5 seconds we have found 2 events in coincidence. This experimental result has to be compared with the number of coincidences due to chance.

The shifts method consists of $`10^4`$ time shifts of the coincidence window; for each one we compute the number of coincidences; if the candidate GWB arrival times and the GRB triggers can both be modeled as Poisson random points, the number of accidental coincidences is fitted to a Poisson curve; from the fit we get the expected number of coincidences due to chance. The results are summarized in Fig. 1. From the fit we obtain a value of mean expected accidental coincidences of $`n_a=2.57\pm 0.04`$.

Another approach to evaluate the number of accidental coincidences is given by Eq. 1, that holds in case of Poisson random points. The total number of AURIGA events and GRB triggers are respectively $`N_A=26816`$, $`N_\gamma =120`$, assuming $`\mathrm{\Delta }T=2W=10`$ s and $`T=1.32\times 10^7`$ s we get $`n_a=2.4\pm 0.2`$. The error on $`n_a`$ can be easily estimated assuming the Poisson statistics for the fluctuations on $`N_A`$ and $`N_\gamma `$ i.e. $`\sqrt{N_A}`$ and $`\sqrt{N_\gamma }`$ respectively. The two estimates of the accidental number of coincidences are in good agreement and demonstrate that the two event rates are uncorrelated. We conclude that the 2 coincidence found are due to chance.

### B Statistical search

We have first tested the method using a Monte Carlo simulation by adding $`N_{on}`$ signals at fixed SNR over a gaussian noise generated by means of a noise model with the same parameters of the AURIGA detector. The simulated detector output is then fed to the same data filtering procedures used for the AURIGA experimental data. We have generated 120 $`\delta `$-like signals with SNR 3, 2 and 1 and 1000 with SNR 1 storing their true arrival times $`t_0`$. The signals have been superimposed to stationary gaussian noise and we have then formed the $`\mathrm{𝑜𝑛}`$ and $`\mathrm{𝑜𝑓𝑓}`$ populations and calculated the $`t`$ value. The probability $`P(t)`$ that the $`t`$ value obtained is due to chance is reported in Table I. For signal with $`SNR=3`$ and $`SNR=2`$ this probability is very small and the value of $`t`$ is statistically significant, showing a GRB-GWB association. On the other hand, for signals with $`SNR=1`$, we have to increase the number of GRBs to 1000 to get a statistical significance.

The $`\chi _{\mathrm{𝑜𝑓𝑓}}`$ and $`\chi _{\mathrm{𝑜𝑛}}`$ sets relative to the AURIGA-BATSE data between 1997 and 1998 are given in the histograms of Fig. 2 and Fig. 3, where the non gaussian tails are due to the non stationarity of the AURIGA noise. The value of the Student’s $`t`$-test obtained is $`t=0.58`$, which corresponds to a probability of $`0.28`$ that it is due to chance.

As we conclude that there is no evidence of an association of GRB-GWB, we can put a constraint on gravitational radiation emitted from GRBs averaged over the source population, $`h_{\mathrm{RMS}}`$. The value of $`\mu _{t,max}`$ can be found setting a confidence level C.L. and solving the equation:

$$_{\mathrm{}}^{\mu _{t,max}}f_d(t)𝑑t=\text{C.L.},$$
(19)
where $`f_d(t)`$ is the distribution function for the Student’s $`t`$ with $`d`$ degrees of freedom. Notice that for high $`d`$, $`f_d(t)`$ tends to a normal curve with zero mean and unitary variance. If we choose C.L.=95% and set $`d=N_{\mathrm{𝑜𝑛}}+N_{\mathrm{𝑜𝑓𝑓}}2`$, we get $`\mu _{t,max}=1.65`$. Using Eq. 18, we set the following upper limit on the averaged g.w. amplitude associated with GRBs:

$$h_{\mathrm{RMS}}1.5\times 10^{18}.$$
(20)
## IV Discussion

The two search methods reported in this paper show no evidence of a correlation between $`\gamma `$-ray bursts and gravitational waves: First, the two coincidences found with the coincidence search have no statistical significance and can be reconduced to chance. Next, the statistical method doesn’t lead to a statistically significative value of the Student’s $`t`$. However we were able to put an upper limit on gravitational signals associated to the GRBs averaged over the source population, $`h_{\text{RMS}}1.5\times 10^{18}`$.

The existence of burst-like excitations (usually referred as non-gaussian noise) in the AURIGA data is well known, and here we deal with this noise selecting the periods of time when the detector was operating in a satisfactory way. The vetoing procedure of data dominated by non gaussian noise has been set by the AURIGA data analysis  and the resulting duty cycle is about 40 % of the total operating time. Care was also taken to cope with non stationarities of the AURIGA noise that give rise to the non gaussian tails in Fig. 2 and Fig. 3 excluding such tails from the fits we used to estimate the statistical variable $`t`$. However, it is important to notice that even with the above limitations the analysis has reached a level of sensitivity which is astrophysically of some interest. In order to get better upper limits we have also explored the possibility of using in our analysis the incoming direction of GRBs but we have found no significant improvements (see Appendix B).

To increase the confidence in our results we have applied the non-parametric Mann-Whitney $`u`$-test  to the sample sets of the statistical search, obtaining a second confirmation of the null hypothesis. Ranking the elements of the union set $`\chi _{\mathrm{𝑜𝑓𝑓}}\chi _{\mathrm{𝑜𝑛}}`$ in increasing order, we get a statistical parameter $`z=1.59`$, smaller than the critic value $`z=1.95`$ imposed by the fixed C.L. of 95 $`\%`$.

The upper limit we have set can be improved in the future simply by the increasing of common data taking periods of AURIGA ($`N_{on}`$ increases) and the new experiment HETE-II (High Energy Transient Explorer ) which is going to substitute by now the wasted BATSE satellite in the GRB search. Another possibility to reach more astrophysically interesting sensitivities is the upgrade of the AURIGA detector which is now in progress. The predicted sensitivity and bandwidth of the AURIGA detector equipped with the new read out system would be respectively $`h_{min}8\times 10^{20}`$ and $`t_w^110`$ Hz. This sensitivity, together with an enhancement of noise stationarity and duty cycle of the detector, would correspond to the lowering of the upper limit $`h_{\mathrm{RMS}}`$ of about 2 order of magnitude in one year of correlation analysis.

## A

To derive Eq. 1 let us consider $`n`$ points random distributed in a time interval $`[0,T]`$. The probability $`P\left(\{k_a,t_a\}\right)`$ that $`k_a`$ points lies in the time window $`t_a=t_2t_1`$ is given by the binomial distribution  with probability $`p=t_a/T`$ that a single point lie in $`t_a`$. If $`n1`$ and $`t_aT`$, using the Poisson theorem we get

$$P\left(\{k_a,t_a\}\right)=e^{nt_a/T}\frac{(nt_a/T)^{k_a}}{k_a!}.$$
(A1)
For $`m`$ not overlapping windows, it can be demonstrated that in the limit of $`n\mathrm{}`$, $`T\mathrm{}`$ and $`n/T`$ constant, the probability of {$`k_1`$ points in $`t_1`$}, …, {$`k_m`$ points in $`t_m`$} is

$`P(\{k_1,t_1\},\mathrm{},\{k_m,t_m\})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{m}{}}}e^{nt_i/T}{\displaystyle \frac{(nt_i/T)^{k_i}}{k_i!}}`$ (A2)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{m}{}}}P\left(\{k_i,t_i\}\right)`$ (A3)
showing that the events {$`k_i`$ points in $`t_i`$} and {$`k_j`$ points in $`t_j`$} are independent for every $`i`$ and $`j`$. Substituting $`m=N_\gamma `$, $`t_i=2W=\mathrm{\Delta }T`$ and $`n=N_A`$ we get

$`P(\{k_1,\mathrm{\Delta }T\},\mathrm{},\{k_{N_\gamma },\mathrm{\Delta }T\})=`$ (A4)
$`=`$ $`e^{N_\gamma N_A\mathrm{\Delta }T/T}\left({\displaystyle \frac{N_A\mathrm{\Delta }T}{T}}\right)^k{\displaystyle \underset{i=1}{\overset{N_\gamma }{}}}{\displaystyle \frac{1}{k_i!}}.`$ (A5)
We can now evaluate the probability $`P(k)`$ to have $`k`$ coincidences, that can be obtained taking into account, at fixed $`k`$, all the possible sets $`[k_i]=\{k_1,\mathrm{},k_{N_\gamma }\}`$ with the constrain $`_{i=1}^{N_\gamma }k_i=k`$; we get

$`P(k)`$ $`=`$ $`{\displaystyle \underset{[k_i]}{}}P(\{k_1,\mathrm{\Delta }T\},\mathrm{},\{k_{N_\gamma },\mathrm{\Delta }T\})`$ (A6)
$`=`$ $`e^{N_\gamma N_A\mathrm{\Delta }T/T}\left({\displaystyle \frac{N_A\mathrm{\Delta }T}{T}}\right)^k{\displaystyle \underset{[k_i]}{}}{\displaystyle \underset{i=1}{\overset{N_\gamma }{}}}{\displaystyle \frac{1}{k_i!}}`$ (A7)
$`=`$ $`e^{N_\gamma N_A\mathrm{\Delta }T/T}\left({\displaystyle \frac{N_A\mathrm{\Delta }T}{T}}\right)^k{\displaystyle \frac{N_\gamma ^k}{k!}}.`$ (A8)
Rearranging the factors in Eq. A8 we get a Poisson distribution as in Eq. 1. The last equality can be easily demonstrated using the multinomial expansion relation

$$\left(\underset{i=1}{\overset{N}{}}x_i\right)^k=k!\underset{[k_i]}{}\underset{i=1}{\overset{N}{}}\frac{x_i^{k_i}}{k_i!}$$
(A9)
and substituting $`x_i=1`$, $`i=1\mathrm{}N`$.

## B

As the incoming direction of GRBs is known, one may wonder if a selection of the GRBs based on a cutoff on the AURIGA antenna pattern, averaged on the unknown polarizations of the GWB, could increase the sensitivity of our analysis.

The antenna pattern of a resonant bar detector is given by

$$F(\theta )=1\left(\widehat{𝐤}\widehat{𝐳}\right)^2=\mathrm{sin}^2(\theta ),$$
(B1)
where $`\widehat{𝐤}`$ and $`\widehat{𝐳}`$ are unitary vectors parallel to the GRB direction and the antenna bar axis and $`\theta `$ is the angle between $`\widehat{𝐤}`$ and $`\widehat{𝐳}`$. Sources that fall outside the two cones, which are defined by the equation $`\widehat{𝐤}\widehat{𝐳}\mathrm{cos}(\xi )`$, have a figure pattern $`F(\theta )\mathrm{sin}^2(\xi )`$ and therefore the average energy associated with these g.w. sources is

$`\mathrm{E}^\xi [h^2]`$ $``$ $`{\displaystyle \frac{15}{16}}{\displaystyle _{\mathrm{\Omega }_\xi }}F^2(\theta ){\displaystyle \frac{d\mathrm{\Omega }}{4\pi }}`$ (B2)
$`=`$ $`{\displaystyle \frac{75}{64}}\left[\mathrm{sin}(\xi )+{\displaystyle \frac{1}{6}}\mathrm{sin}(3\xi )+{\displaystyle \frac{1}{50}}\mathrm{sin}(5\xi )\right],`$ (B3)
where the solid angle $`\mathrm{\Omega }_\xi `$ is defined by $`F(\theta )F(\xi )`$. Moreover, as GRBs are isotropically distributed over the sky, the cutoff on the figure pattern decreases the number of available GRB: $`N_{on}^\xi =N_{on}\mathrm{sin}(\xi )`$. Therefore the net effect of this cutoff on the expected value of $`t`$ in Eq. 11 is

$`\mu _t^\xi `$ $`=`$ $`\mu _t{\displaystyle \frac{75}{64}}[\mathrm{sin}(\xi )+{\displaystyle \frac{1}{6}}\mathrm{sin}(3\xi )+{\displaystyle \frac{1}{50}}\mathrm{sin}(5\xi )]\times `$ (B4)
$`\times `$ $`\mathrm{sin}(\xi )^{1/2}.`$ (B5)
The function $`\mu _t^\xi /\mu _t`$ is a continuously increasing function in the range $`\xi [0,\pi /2]`$, and $`\mu _t^\xi /\mu _t=1`$ at $`\xi =\pi /2`$ (i.e. the whole solid angle). We must conclude that a cutoff on the GRB direction does not enhance the sensitivity of the statistical search.


text:
# Precession Interpretation of the Isolated Pulsar PSR B1828-11

## 1 Introduction

The monitoring of pulsar spin behavior through precision timing measurements provides a means for studying the dynamics and internal structure of neutron stars. In this connection, detection of precession in pulsars is of particular interest. Relativistic (geodetic) precession has been observed in the double neutron star binaries PSRs B1913+16 and B1534+12 (Weisberg, Romani & Taylor 1989; Arzoumanian 1995), and provides confirmation of General Relativistic orbital mechanics. Classical precession in isolated pulsars could probe the coupling of the stellar crust to the interior, but its detection is complicated by the presence of timing noise. Evidence for long-period precession is seen in the Crab pulsar (Lyne, Pritchard & Smith 1988), the Vela pulsar (Deshpande & McCulloch 1996) and PSR B1642-03 (Cordes 1993; Shabanova, Lyne & Urama 2001). The 35-d periodicity seen in the accreting system Her X-1 (Tannanbaum et al. 1972) has been interpreted as precession by many authors (e.g., Brecher 1972; Trümper et al. 1986; Cadez, Galicic & Calvani, 1997; Shakura, Postnov & Prokhorov 1998).

Recently, Stairs, Lyne & Shemar (2000) discovered the most compelling evidence to date for precession of an isolated neutron star. The period residuals of PSR B1828-11 have an amplitude of $`1`$ ns and have strong Fourier power at periods of $`1000`$ d, 500 d and 250 d. These Fourier components alone might suggest a Doppler effect through gravitational interaction with planets. However, correlations of the pulse shape and duration with variations in the period derivative $`\dot{p}`$ strongly suggest that the viewing angle with respect to the emission beam is varying, as one would expect if the star were precessing (see Fig. 1). There is no obvious damping of the observed fluctuations, though the timing behavior is not entirely stable. See Table 1 for the spin parameters of PSR B1828-11.

The purpose of this paper is to investigate if the timing behavior and pulse profile changes of PSR B1828-11 can be quantitatively understood in terms of stellar precession. We assume that the observed timing residuals represent true changes in the pulsar’s rotation rate. An immediate question is which of the three observed Fourier components corresponds to the actual precession frequency. The residuals of the period $`p`$ and period derivative $`\dot{p}`$ show the strongest Fourier power at $`500`$ d (see Fig. 1 and Stairs, Lyne & Shemar 2000 Fig. 3); we tentatively assume that this Fourier component represents the precession period. The wobble angle $`\theta `$ between the star’s symmetry axis and its angular momentum can be estimated from the pulse profile changes. The difference in widths between the “wide” and “narrow” beam profiles is $`2^{}`$ (see Fig. 2), implying a wobble angle of similar magnitude for plausible beam shapes.

This precession interpretation could be complicated by magnetospheric processes. Precession of the star could create changes in the magnetosphere that affect the geometry and location of the emission region. Changes in the height of the emission region would introduce light-travel time contributions to the pulse arrival times. Different regions of the magnetosphere might be viewed as the star precesses. We do not consider these effects in this paper.

In the next section, we give an overview of precession. In §3 we study how the electromagnetic torque, responsible for the star’s spin down (at least in part), modifies precession. In §4, we relate this torque effect to the observed timing. In §5, we propose a model of the the observed pulse duration variations. In §6 we discuss some of our chief results. We conclude with discussion of observational tests of the precession interpretation.

## 2 Overview

Most treatments of pulsar precession to date have considered mainly free precession (e.g., Goldreich 1970; Pines & Shaham 1974; Shaham 1977; Suto & Iso 1985; Orford 1987; Bisnovatyi-Kogan, Mersov & Sheffer 1989; Glendenning 1990; Nelson, Finn & Wasserman 1990; Bisnovatyi-Kogan & Kahabka 1993; Cadez, Galicic & Calvani 1997; Shakura, Postnov & Prokhorov 1998; Sedrakian, Wasserman & Cordes 1999; Jones & Andersson 2000). The observed variations in $`p`$ and $`\dot{p}`$ in PSR B1828-11, however, are far too large to be explained by free precession for a wobble angle of $`\theta 2^{}`$. In §3 we calculate the precession of a neutron star subject to a vacuum magnetic dipole braking torque and find a solution that is consistent with observations.

Consider first the free precession of a biaxial neutron star; the geometry is as shown in Fig. 3. At any instant, the symmetry axis $`\widehat{e}_3`$, angular momentum $`𝐋`$ and angular velocity $`𝝎`$ span a plane. Let the oblateness be $`ϵ(I_3I_1)/I_1`$, where $`I_3`$ and $`I_1`$ are the largest and smallest moments of inertia about the principal axes. For small oblateness $`ϵ1`$, the angle between $`𝐋`$ and $`𝝎`$ is $`\widehat{\theta }ϵ\theta \theta `$ (see, e.g., Jones & Andersson 2000). In the inertial frame, $`\widehat{e}_3`$ and $`𝝎`$ rotate in a right-handed sense about $`𝐋`$ at frequency $`\omega (1+ϵ)`$. Combined with this motion is a slow retrograde motion of the body about $`\widehat{e}_3`$ at precession frequency $`\omega _pϵ\omega `$.

For a precessing pulsar, the pulse emission times are defined by the instants at which the beam center passes closest to the direction of the observer; the observed arrival times are simply offset from the emission times (apart from Doppler shifts from translational motion). If the star is freely precessing, the pulse arrival times vary sinusoidally as $`\mathrm{\Delta }t(p_0/2\pi )\theta \mathrm{cot}\chi \mathrm{sin}(\omega _pt+\beta )`$ where $`p_0`$ is the spin period, $`\chi >\theta `$ is the inclination angle between the center of the pulsar beam and the body’s symmetry axis and $`\beta `$ is a phase (see Nelson, Finn & Wasserman 1990). This variation is a purely geometrical effect, as the magnitude of the angular velocity is unchanged by free precession. The corresponding changes in period derivative due to this geometrical effect are $`\mathrm{\Delta }\dot{p}(2\pi )^1(p_0\omega _p)^2\theta \mathrm{cot}\chi \mathrm{sin}(\omega _pt+\beta )`$ (see eq. 25 below). For a precession period of 500 d and $`\chi >30^{}`$, the residuals in $`\dot{p}`$ are less than 0.02 s s-1, a factor of at least 10 smaller than observed.

The electromagnetic torque on a neutron star modifies the precession and increases the magnitude of timing residuals for a given wobble angle (Jones 1988; Cordes 1993; Jones & Andersson 2000). In the body frame of a precessing object, the angular velocity vector $`𝝎`$ moves about the body’s symmetry axis $`\widehat{e}_3`$ at frequency $`\omega _p`$ in a right-handed sense. For a magnetized neutron star with a magnetic moment misaligned with respect to the symmetry axis, the angle between the dipole moment $`𝐦`$ and $`𝝎`$ changes by about $`\theta `$ over a precession period (see eq. 26 below). The external dipole torque, which depends on this angle, thus varies over a precession period in addition to its slow decay over an evolutionary timescale. For small wobble angle $`\theta `$, these torque changes produce changes in the star’s spin rate $`\omega `$ of (see eq. 16)

$$\frac{\mathrm{\Delta }\dot{\omega }}{\omega _0}=\frac{\theta }{2\tau }\left[\mathrm{sin}2\chi \mathrm{cos}(\omega _pt+\beta )+\frac{\theta }{2}\mathrm{sin}^2\chi \mathrm{cos}(2\omega _pt+2\beta )\right],$$
(1)
For simplicity we assume that the pulsed radiation is emitted along the axis of the magnetic dipole, i.e., both the beam and the dipole are at an angle $`\chi `$ with respect to the body’s symmetry axis. The torque-induced timing variations are much larger than the geometrical effects associated free precession (Cordes 1993; Jones 1988; see eq. 25 below). In addition, and as pointed out by Jones & Andersson (2000), the dependence of the torque on the angle between the star’s angular velocity and its dipole moment creates a strong harmonic of the precession frequency in the star’s spin rate if the star is a nearly orthogonal rotator. In §5 we apply this picture to PSR B1828-11 and find that $`\chi =89^{}`$ and a wobble angle of $`\theta 3^{}`$ can account for the 500-d Fourier component with a harmonic at 250 d. Precession does not give the whole story, however, as it does not explain the Fourier power at 1000 d.

As the star precesses we see different sweeps through the beam of the pulsar, allowing us to partially map the beam morphology. As we will show, the observed variations in pulse duration imply that the beam has an “hour glass” shape, i.e., the pulse duration is longer for sweeps farthest from the beam center.

## 3 Precession Under the Dipole Torque

Precession of a neutron star will occur if the symmetry axis of the crust is not aligned with the angular momentum. The elastic crust has two contributions to its oblateness: a centrifugal deformation associated with the star’s spin and a Coulomb deformation sustained by the rigidity of the lattice (Alpar & Pines 1985; Jones & Andersson 2000). The principal axis of inertia of the centrifugal deformation follows the instantaneous spin axis of the crust, and so does not affect the precession. In §6 we discuss possible origins of Coulomb deformation.

The period and amplitude of the crust’s precession are affected by the coupling of the crust to the liquid interior. The precession creates time-dependent velocity differences between the crust and liquid that vary over the star’s spin period. If the coupling time $`\tau _f`$ between the crust and the liquid interior is much longer than the crust’s spin period $`p`$, the precession will damp over $`2\pi \tau _f/p`$ precession periods (Bondi & Gold 1955; Sedrakian, Wasserman & Cordes 1999). All studies of damping of differential rotation between the crust and various parts of the liquid interior give $`\tau _fp`$ (e.g., Alpar & Sauls 1988; Epstein & Baym 1992; Jones 1992; Abney, Epstein & Olinto 1996; Mendell 1998). To a good approximation, therefore, the liquid interior can be treated as decoupled from the solid (Bondi & Gold 1955; Sedrakian, Wasserman & Cordes 1999); in this regime the crust precesses almost as if the liquid interior were not there.111A precessing rigid shell containing a non-spherical normal fluid is inertially coupled to the fluid (Lamb 1952). The fluid, which tends to a configuration symmetric about its rotation axis, exerts a reaction force on the solid. The relevance of this effect to an elastic neutron star containing superfluid neutrons and protons is unclear. Inertial coupling would not change the precession dynamics that we calculate in this paper, though it would change the dependence of the precession period on the crust oblateness.

In our calculations we treat the neutron star as consisting of a biaxial, rigid crust afloat on a liquid core. The solid portion of the star comprises less than $`1\%`$ of the star’s total moment of inertia and we regard the liquid of the inner crust as an extension of the core liquid. We neglect angular momentum exchange between the solid and the liquid. For simplicity we ignore the centrifugal bulge (which does not affect the precession) and treat the crust as perfectly rigid. For illustration we describe the spin dynamics of the crust with the vacuum dipole torque of Davis & Goldstein (1970). While plasma effects could contribute to the spin down, a torque similar to the vacuum torque will always be present. In a frame corotating with the crust, the equations of motion are (Goldreich 1970; Melatos 2000)

$$𝐈\frac{d𝝎}{dt}+𝝎\times 𝐋=\frac{2\omega ^2}{3c^3}(𝝎\times 𝐦)\times 𝐦+\frac{1}{Rc^2}(𝝎𝐦)(𝝎\times 𝐦),$$
(2)
where $`𝐈`$ is the moment of inertia tensor of the crust, $`𝝎`$ is the angular velocity, $`𝐋`$ is the angular momentum, $`𝐦`$ is the magnetic dipole moment (centered on the star) and $`R`$ is the radius of the star. The first contribution to the torque is due to the radiation far field and is responsible for spinning down the star. The second term represents the near-field radiation torque. The coefficient of this term depends on how the star’s internal magnetization is treated (Melatos 2000); its precise value does not affect our results. We will find that variations in the spin rate associated with precession are determined predominantly by the far-field torque (see. eq. 16).

We define a Cartesian coordinate system (1,2,3) with the 3 axis along the symmetry axis $`\widehat{e}_3`$ of the crust. Without loss of generality, we take the magnetic dipole to lie in the 1-3 plane, inclined with respect to the 3 axis by an angle $`\chi <\pi /2`$. At $`t=0`$ the spin rate is $`\omega _02\pi /p_0`$. The equations of motion are

$`I_1\dot{\omega }_1+\omega _2\omega _3(I_3I_1)`$ $`=`$ $`T_1={\displaystyle \frac{1}{2\tau }}I_3(\omega _3\mathrm{sin}\chi \omega _1\mathrm{cos}\chi )\mathrm{cos}\chi `$ (4)
$`+{\displaystyle \frac{3x_0^1}{4\tau }}I_3(\omega _1\mathrm{sin}\chi +\omega _3\mathrm{cos}\chi )\omega _2\mathrm{sin}\chi `$
$`I_1\dot{\omega }_2\omega _1\omega _3(I_3I_1)`$ $`=`$ $`T_2={\displaystyle \frac{1}{2\tau }}I_3\omega _2`$ (6)
$`+{\displaystyle \frac{3x_0^1}{4\tau }}I_3(\omega _1\mathrm{sin}\chi +\omega _3\mathrm{cos}\chi )(\omega _3\mathrm{sin}\chi \omega _1\mathrm{cos}\chi )`$
$`I_3\dot{\omega }_3`$ $`=`$ $`T_3={\displaystyle \frac{1}{2\tau }}I_3\mathrm{sin}\chi (\omega _1\mathrm{cos}\chi \omega _3\mathrm{sin}\chi )`$ (8)
$`{\displaystyle \frac{3x_0^1}{4\tau }}I_3(\omega _1\mathrm{sin}\chi +\omega _3\mathrm{cos}\chi )\omega _2\mathrm{sin}\chi ,`$
where $`\tau 3c^3I_3/4m^2\omega _0^2`$ is the characteristic time over which the far-field torque spins down the star, $`m|𝐦|`$ and $`x_0R\omega _0/c<<1`$. For the oblate star, $`I_3>I_2=I_1`$. Terms on the right-hand side of eqs. - proportional to $`x_0^1`$ arise from the near-field torque; the remaining terms are due to the far-field torque. The near-field torque is stronger than the far-field torque by a factor $`x_0^1`$, though its direction is always orthogonal to $`𝝎`$.

Identification of the small parameters in the problem allows perturbative treatment of the torque. Because the torque acts over timescales that are much longer than the observed precession period, we expect the spin solution to resemble that for free precession. Hence, the ratio $`\omega _p/\omega _0`$, where $`\omega _p`$ is the frequency of free precession, should be approximately $`ϵ(I_3I_1)/I_1`$. This ratio is $`10^8`$ for PSR B1828-11 if the precession period is $`500`$ d. Changes in the pulse duration (Fig. 2) through the precession period suggest a wobble angle of $`\theta 2^{}`$. For these small values of $`ϵ`$ and $`\theta `$, $`\widehat{\theta }ϵ\theta \theta `$, so $`\theta `$ and $`\alpha `$ are nearly equal (Jones & Andersson 2000). In terms of the components of the angular velocity, $`\alpha `$ is approximately $`(\omega _1^2+\omega _2^2)^{1/2}/\omega _01`$. We therefore regard $`\omega _1`$ and $`\omega _2`$ as perturbations about the secular spin down. The relative magnitudes of the small dimensionless parameters are given by

$$1\alpha ϵϵ\alpha >x_0^1\frac{1}{\omega _0\tau }\frac{1}{\omega _0\tau }.$$
(9)
To order $`ϵ\omega _1`$ and $`ϵ\omega _2`$, eqs.  and  are

$`\dot{\omega _1}+ϵ\omega _2\omega _0=0`$ (10)
$`\dot{\omega _2}ϵ\omega _1\omega _0=0.`$ (11)
For $`\alpha 1`$, the solutions are

$`\omega _1`$ $`=`$ $`\omega _0\alpha \mathrm{cos}(\omega _pt+\beta )`$ (12)
$`\omega _2`$ $`=`$ $`\omega _0\alpha \mathrm{sin}(\omega _pt+\beta )`$ (13)
$`\omega _p`$ $`=`$ $`ϵ\omega _0={\displaystyle \frac{I_3I_1}{I_1}}\omega _0,`$ (14)
where $`\beta `$ is a phase. The behavior of the component of the angular velocity orthogonal to the symmetry axis is essentially unaltered by the torque; the precession frequency is that of free precession and $`\alpha \mathrm{sin}\alpha =(\omega _1^2+\omega _2^2)^{1/2}/\omega _0=`$constant. In the body frame, both $`𝐋`$ and $`𝝎`$ (which are always coplanar with $`\widehat{e}_3`$) precess slowly at frequency $`\omega _p\omega _0`$ about $`\widehat{e}_3`$ and in the same sense. For the limit $`ϵ1`$ that we are considering, $`\alpha `$ and $`\theta `$ are nearly equal.

For biaxial free precession, $`\omega |𝝎|`$ is constant. Under the dipole torque, however, $`\omega `$ undergoes changes that are important for the observed timing behavior. Combining eqs. - gives

$$\frac{d\omega ^2}{dt}=\frac{2}{I_1}\left(𝝎𝐓ϵ\frac{I_1}{I_3}\omega _3T_3\right).$$
(15)
The $`𝝎𝐓`$ term does not depend on the near-field torque. The near-field torque changes the spin rate only through the small final term, which we will neglect. The secular spindown in the absence of precession is given by $`\dot{\omega }\omega \mathrm{sin}^2\chi /2\tau `$. Combining eq.  with eqs. -, setting $`\alpha =\theta `$, retaining terms to order $`\theta ^2(\omega _0\tau )^1`$, and subtracting the secular spindown gives (dropping constant terms):

$$\frac{\mathrm{\Delta }\dot{\omega }}{\omega _0}=\frac{\theta }{2\tau }\left[\mathrm{sin}2\chi \mathrm{cos}(\omega _pt+\beta )+\frac{\theta }{2}\mathrm{sin}^2\chi \mathrm{cos}(2\omega _pt+2\beta )\right].$$
(16)
Hence $`\dot{\omega }`$ undergoes small variations of order $`\theta (\omega _0\tau )^1`$ and $`\theta ^2(\omega _0\tau )^1`$, driven by changes in the far-field dipole torque. These torque variations are due to changes in the angle between $`𝝎`$ and $`𝐦`$ as $`𝝎`$ precesses through the body. The harmonic at twice the precession frequency arises from quadratic dependences of $`\omega _1T_1`$ and $`\omega _2T_2`$ on the components of the angular velocity; the harmonic dominates the fundamental for $`\theta \mathrm{tan}\chi >4`$. For a small $`\theta `$, the dipole moment must be nearly orthogonal to the rotation axis for the harmonic to be significant. The small oscillations of $`\dot{\omega }`$ largely determine the observed pulse timing; they give period variations $`\mathrm{\Delta }\dot{p}(p_0^2/2\pi )\mathrm{\Delta }\dot{\omega }`$.

## 4 Timing

The observed pulse timing at a given time in the precession cycle depends on the star’s spin rate and the orientation of the angular velocity $`𝝎`$ with respect to the observer. To relate the solution of the previous section to the observed pulse timing, we go to an inertial coordinate system $`x,y,z`$ with the angular momentum along the $`z`$ axis and the observer in the $`xz`$ plane with $`x>0`$ and $`z>0`$. The azimuthal and polar angles of the magnetic dipole $`𝐦`$, $`\mathrm{\Phi }`$ and $`\mathrm{\Theta }`$, are given in this frame by

$`\mathrm{tan}\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{m_y}{m_x}}={\displaystyle \frac{\left(\mathrm{cos}\psi \mathrm{sin}\varphi +\mathrm{cos}\theta \mathrm{cos}\varphi \mathrm{sin}\psi \right)\mathrm{sin}\chi \mathrm{sin}\theta \mathrm{cos}\varphi \mathrm{cos}\chi }{\left(\mathrm{cos}\psi \mathrm{cos}\varphi \mathrm{cos}\theta \mathrm{sin}\varphi \mathrm{sin}\psi \right)\mathrm{sin}\chi +\mathrm{sin}\theta \mathrm{sin}\varphi \mathrm{cos}\chi }}`$ (17)
$`\mathrm{cos}\mathrm{\Theta }`$ $`=`$ $`{\displaystyle \frac{m_z}{m}}=\mathrm{sin}\theta \mathrm{sin}\psi \mathrm{sin}\chi +\mathrm{cos}\theta \mathrm{cos}\chi ,`$ (18)
where $`\psi `$, the wobble angle $`\theta `$ and $`\varphi `$ are Euler angles following the definitions of Landau & Lifshitz (1976) and Goldstein (1980).

The pulsar beam is not necessarily in the same direction as the dipole moment. Nevertheless, for simplicity we define a pulse as occurring when the azimuthal angle $`\mathrm{\Phi }`$ of the magnetic dipole equals the azimuth of the observer, i.e., when $`m_y=0`$ and $`\mathrm{\Phi }=0`$. Hence $`\dot{\mathrm{\Phi }}`$ is the observed pulse frequency. For small $`\theta <\chi `$ and $`\theta \mathrm{cot}\chi 1`$, eqs.  and  give to second order in $`\theta `$,

$$\dot{\mathrm{\Phi }}=\dot{\varphi }+\dot{\psi }+\theta \mathrm{cot}\chi \dot{\psi }\mathrm{sin}\psi \frac{\theta ^2}{2}(1+2\mathrm{cot}^2\chi )(2\mathrm{cos}^2\psi 1)\dot{\psi }.$$
(19)
The angles $`\psi `$ and $`\varphi `$ are given by

$`\psi `$ $`=`$ $`\mathrm{tan}^1{\displaystyle \frac{L_1}{L_2}}=\mathrm{tan}^1{\displaystyle \frac{\omega _1}{\omega _2}}={\displaystyle \frac{\pi }{2}}\omega _pt\beta `$ (20)
$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{\omega _3\dot{\psi }}{\mathrm{cos}\theta }}\omega +\omega _p+{\displaystyle \frac{1}{2}}\theta ^2\omega _p,`$ (21)
The observed beam sweep rate is (dropping constant terms)

$$\dot{\mathrm{\Phi }}=\omega \theta \omega _p\mathrm{cot}\chi \mathrm{cos}(\omega _pt+\beta )\theta ^2\omega _p(1+2\mathrm{cot}^2\chi )\mathrm{cos}^2(\omega _pt+\beta )$$
(23)
Hence, the sweep rate is the spin rate plus a modulation at the free precession frequency and its first harmonic. This modulation is a geometrical effect due to variations in the rate at which the beam sweeps past the observer as the symmetry axis precesses about the angular momentum; it is the only modulation that occurs if the precession is free. The sweep rate also changes as $`\omega `$ changes through secular spin-down and the variations $`\mathrm{\Delta }\omega `$ associated with torque variations as the star precesses. We henceforth focus on behavior of the period derivative since, as we will show, the torque changes evaluated in the previous section largely determine the observe spin behavior. Differentiating eq.  in time and combining with eqs. ,  and \[LABEL:phi\], gives the period derivative residuals with respect to the secular spin down:

$`\mathrm{\Delta }\dot{p}{\displaystyle \frac{p_0^2}{2\pi }}\ddot{\mathrm{\Phi }}`$ $``$ $`{\displaystyle \frac{p_0}{t_{\mathrm{age}}}}\theta \left[\mathrm{cot}\chi \mathrm{cos}(\omega _pt+\beta )+{\displaystyle \frac{\theta }{4}}\chi \mathrm{cos}(2\omega _pt+2\beta )\right]`$ (24)
$``$ $`{\displaystyle \frac{1}{2\pi }}(p_0\omega _p)^2\theta \left[\mathrm{cot}\chi \mathrm{sin}(\omega _pt+\beta )+\theta (1+2\mathrm{cot}^2\chi )\mathrm{sin}(2\omega _pt+2\beta )\right].`$ (25)
Here $`t_{\mathrm{age}}\tau /\mathrm{sin}^2\chi `$ is the characteristic spin-down age. This equation gives the complete expression for the star’s residuals in $`\dot{p}`$. The first set of square brackets gives the residuals due to torque variations from eq. . The second set of square brackets gives the residuals due to geometrical effects (see Nelson, Finn & Wasserman 1990; Jones & Andersson 2000). For PSR B1828-11, the first-order torque effects always dominate the geometrical effects by a factor of $`250`$ if the precession period is $`500`$ d. For the harmonic at frequency $`2\omega _p`$ to play a role in the timing (terms proportional to $`\theta ^2`$), $`\chi `$ must be nearly $`\pi /2`$. The near-field torque, because it does not significantly change $`\omega `$, is unimportant for the small-$`\theta `$ precession that appears to be taking place in PSR B1828-11.

We model the data for the 2000-day span beginning at MJD 49,000 as the observations were most closely spaced during this period and these data show the observed periodicities most clearly. Fitting eq.  to the data, neglecting the small geometrical terms, gives the parameters $`\theta `$, $`\omega _p`$, $`\beta `$ and $`\chi `$ featured in Table 2. Comparison of the fit to the data is shown in Fig. 1.

## 5 Pulse Duration Variations

The duration of PSR B1828-11’s pulses is seen to change with the same periodic structure as the spin residuals, suggesting different sweeps through the beam as the star precesses. To study this further, consider the behavior of the polar angle $`\mathrm{\Theta }`$ of the beam with respect to $`𝐋`$ (eq. 18). Expanding eq.  to first order in $`\theta `$, using eq. , gives

$$\mathrm{\Theta }=\chi \theta \mathrm{cos}(\omega _pt+\beta ).$$
(26)
As the star precesses, $`\mathrm{\Theta }`$ varies sinusoidally about $`\chi `$. The variations in the polar angle are in phase with the dominant contribution to the torque variations (the first term in eq. ).

To quantify this further, suppose the constant angle between the observer and the angular momentum is $`\gamma \xi +\chi `$ (see Fig. 4). We define the sweep angle $`\mathrm{\Delta }\mathrm{\Theta }`$ as the difference in polar angle between the observer and the dipole at the time of the pulse: $`\mathrm{\Delta }\mathrm{\Theta }=\xi +\theta \mathrm{cos}(\omega _pt+\beta )`$. As a simple model of the variations in pulse duration, we suppose that the pulse duration $`w`$ is a function of $`\mathrm{\Delta }\mathrm{\Theta }`$ only. We assume that $`w`$ has an extremum (it will turn out to be a minimum) for $`\mathrm{\Delta }\mathrm{\Theta }=0`$, corresponding to a viewing angle co-linear with the dipole axis. We take the expansion

$$w=w_0+w_2(\mathrm{\Delta }\mathrm{\Theta })^2.$$
(27)
Consistent with Stairs, Lyne & Shemar (2000), we define a shape parameter $`S`$ as

$$w\mathrm{max}(w)S+\mathrm{min}(w)(1S),$$
(28)
where max $`(w)`$ and min $`(w)`$ are the maximum and minimum values of the beam duration in the precession cycle. Combining eqs.  and  gives

$$S=\{\begin{array}{cc}1\left(\frac{\mathrm{\Delta }\mathrm{\Theta }}{|\xi |+\theta }\right)^2\hfill & w_2>0\hfill \\ \left(\frac{\mathrm{\Delta }\mathrm{\Theta }}{|\xi |+\theta }\right)^2\hfill & w_2<0.\hfill \end{array}$$
(29)
This shape parameter depends on $`\xi `$, $`\theta `$, $`\omega _p`$, $`\beta `$, and the sign of $`w_2`$ but not its magnitude.

With $`\theta `$, $`\omega _p`$ and $`\beta `$ given by the fit to the timing data, the shape parameter is determined entirely by $`\xi `$ and the sign of $`w_2`$. Interestingly, we find that $`w_2`$ must be positive, i.e., the beam duration is larger for beam sweeps farthest from the dipole axis. Such a beam pattern is not standard, but is required by the precession interpretation of PSR B1828-11. The best-fit value of $`\xi `$ is $`0.3^{}`$, as given in Table 2. Comparison of the fit to the data is shown in Fig. 1.

## 6 Discussion

We find an acceptable fit to the data with $`\chi =89^{}`$ between the magnetic dipole and the star’s symmetry axis, a precession period of 511 d and a wobble angle $`\theta =3.2^{}`$. The predicted shape changes are in qualitative agreement with the observations if the angle between the observer and the angular momentum is $`88.7^{}`$ and the radio beam pattern has the “hourglass” shape sketched in Fig. 5. The beam shape that we have inferred differs form the usual picture of a roughly circular pulsar beam, but might be explainable in terms of a patchy beam of the type found by Han & Manchester (2000) in their statistical study of pulsar beam morphology. If the emission is beamed both parallel and anti-parallel to the dipole axis, an inner pulse would be expected for the nearly-orthogonal inclination that we have inferred. None is seen, so the emission pattern must be more complex.

Our model does not account for 1000-d Fourier component seen in the data. However, significant aperiodic variations in the timing seen over timescales of years indicate that more than just precession is at work in PSR B1828-11. Slow changes in the emission region could be a contributing factor. Perturbations to the star’s figure brought about by relaxation of the star’s structure as it spins down could also play a role. Such starquakes would excite precession by perturbing the star’s angular velocity with respect to the symmetry axis (Link, Franco & Epstein 1998). The perturbations would introduce irregular contributions to the magnitude of the wobble angle and the precessional phase (Jones 1988).

Andersson & Jones (2000) have suggested that the 1000-d Fourier component represents the precession period, and that near orthogonality of the dipole moment with respect to the symmetry axis gives a harmonic at 500 d in the spin residuals. This scenario is unlikely for several reasons: 1) near perfect orthogonality ($`\chi \genfrac{}{}{0pt}{}{_>}{^{}}89.9^{}`$) is required to keep the first-order torque term in eq.  from giving too large a contribution to $`\dot{p}`$, 2) the strong Fourier component at 250 d is not accounted for, and, 3) the shape parameter, which would now be periodic at 1000 and 500 d, would bear little resemblance to the data.

The stellar oblateness inferred from the precession period of 511 d is $`ϵ=(I_3I_1)/I_1=9\times 10^9`$. For comparison, we estimate the maximum Coulomb oblateness. The maximum mass $`M_m`$ of a “mountain” that the rigid neutron star crust can support is given roughly by $`M_mgB\theta _c/R`$ where $`g`$ is the gravitational acceleration at the stellar surface, $`B`$ is the bulk modulus of the crust ($`10^{48}`$ erg) and $`\theta _c`$ is the critical strain angle at which the crust fractures. The maximum oblateness sustainable by crust rigidity is thus (see also Ushomirsky, Cutler & Bildsten 2000)

$$ϵ\frac{I_3I_1}{I_1}\frac{M_m}{M_c}10^5\frac{\theta _c}{10^2},$$
(30)
where $`M_c`$ is the mass of the crust. A mountain could form as the result of a starquake that perturbs the star’s principal axis of inertia away from its angular momentum axis (Link, Franco & Epstein 1998). However, even if the crust is relaxed (i.e., unstrained) before it is set into precession, rigidity prevents a portion of the star’s bulge from following the instantaneous rotation axis. The effective oblateness due to rigidity is in this case $`ϵ_r10^5ϵ_\mathrm{\Omega }`$, where $`ϵ_\mathrm{\Omega }`$ is the centrifugal oblateness (Alpar & Pines 1985; see also Munk & MacDonald 1960). For PSR B1828-11’s spin period $`ϵ_r2\times 10^{11}(p/0.4\text{ s})^2`$, giving a precession period of 700 yr. We conclude that PSR B1828-11’s crust is strained to the extent that the Coulomb oblateness is $`100`$ times larger than $`ϵ_r`$. A significantly strained crust can precess at a much higher frequency than the unstrained crust considered in some previous studies (e.g., Goldreich 1970 and Melatos 2000).

The precession of PSR B1828-11 seriously challenges the notion that superfluid vortices of the crust pin to nuclei (e.g., Anderson & Itoh 1975). Pinning would drastically alter the way in which the star precesses. Shaham (1977) showed that very effective pinning changes the precession period to $`(I_1/I_p)p`$, where $`I_p`$ is the portion of the star’s fluid moment of inertia that is pinned and $`I_1`$ is the moment of inertia of the part of the star that is precessing. Sedrakian, Wasserman & Cordes (1999) have shown that Shaham’s conclusion is essentially unaltered if the pinning is imperfect. Pulsar glitches, which might arise from variable coupling between the superfluid and the crust, indicate that $`I_p/I_{\mathrm{star}}>1.4`$% in stars that frequently glitch (Link, Epstein & Lattimer 1999); $`I_{\mathrm{star}}`$ is the total moment of inertia of the star. This degree of pinning in PSR B1828-11 would give a precession period of $`40`$ s, far shorter than the observed of 500 d precession period. It may be that vortex pinning is inhibited in stars that undergo precession with amplitudes as large as those in PSR B1828-11.

## 7 Observational Tests

Our proposed model for the fluctuations in the timing residuals of PSR B1828-11 implies that the pulsed emission should have several specific attributes that may be observationally tested. Fig. 5 illustrates these features that are characteristic of our model. During the 511 day precession cycle the observer sees both the broad upper and lower parts of the beam. One important test of our model would be evidence of these two different parts of the emission region.

We modeled the observed pulse durations with the parabolic form of eq. (27), giving the overall beam shape illustrated by the dashed lines in Fig. 5. The actual shape of the emission must be more complex, however, as we suggest with the irregular brightness contours also shown in the figure. Additionally, the linear polarization of the emission region is unlikely to be the same in the upper and lower portions of the emission region, as indicated by the double- headed arrows in the figure.

To indicate the possible signatures one may observe, we have marked certain viewing angles (A-C) in Fig. 5 and the corresponding observing times in Fig. 1. The dotted horizontal line in Fig. 5 is the shape equator (perpendicular to the largest principal moment). The horizontal dotted lines A and C give the range of the viewing angles seen from Earth. Line B represents the viewing angle at which the pulse profile is the narrowest. The pulse profile is wider at viewing angle A and widest at viewing angle C.

While the shape parameter is the same at phases above and below phase B, the precise pulse profiles and radio polarization should be noticeably different at these phases. The polarization angle will vary across the emission region so that the polarization direction will change continuously from phase A to C and then reverse going from C to A. Detection of these changes in pulse shape and polarization would indicate that the radio beam has the hour-glass shape we have inferred, and would support the interpretation that PSR B1828-11 is a precessing neutron star.

We thank A. G. Lyne, D. I. Jones, R. Bandiera, A. Melatos, J. Cordes, D. Chernoff and the other participants of the ITP workshop on Spin and Magnetism in Young Neutron Stars for helpful discussions. We are grateful to I. Stairs and A. G. Lyne for providing us with the timing data and beam templates for PSR B1828-11. The work was performed under the auspices of the U.S. Department of Energy, and was supported by NASA ATP grant NAG 53688, by IGPP at LANL and in part by the National Science Foundation under Grant No. PHY94-07194.