Characterising the rotational irregularities of the Vela pulsar from 21~yr of phase-coherent timing
Pulsars show two classes of rotational irregularities that can be used to understand neutron-star interiors and magnetospheres: glitches and timing noise. Here we present an analysis of the Vela pulsar spanning nearly $21$~yr of observation and including $8$~glitches. We identify the relative pulse number of all of the observations between glitches, with the only pulse-number ambiguities existing over glitch events. We use the phase coherence of the timing solution to simultaneously model the timing noise and glitches in a Bayesian framework, allowing us to select preferred models for both. We find the glitches can be described using only permanent and transient changes in spin frequency, i.e., no step changes in frequency derivative. For all of the glitches, we only need two exponentially decaying changes in spin frequency to model the transient components. In contrast to previous studies, we find that the dominant transient components decay on a common $\approx 1300$~d time scale, and that a larger fraction ( $\gtrsim 25\%$) of glitch amplitudes are associated with these transient components. We also detect shorter-duration transient components of $\approx 25$~d, as previously observed, but are limited in sensitivity to events with shorter durations by the cadence of our observations. The timing noise is well described by a steep power-law process that is independent of the glitches and subdominant to the glitch recovery. The braking index is constrained to be $< 8$ with 95\% confidence. This methodology can be used to robustly measure the properties of glitches and timing noise in other pulsars.
vela_timing_6.tex

%RMS 21 Oct 2015

\pdfminorversion=5

\documentclass[useAMS,usenatbib]{mn2e}
%\documentclass[iop,numberedappendix]{emulateapj}
%% use this for emulateapj
%\documentclass[letterpaper,11pt]{article}

%\documentclass[preprint]{aastex}




\usepackage{graphicx}
\usepackage{bm}
\usepackage{psfrag}
\usepackage{mathrsfs}
\usepackage{amssymb,amsmath}
%\usepackage{astronat}
\citestyle{aa}

\newcommand{\be}{\begin{eqnarray} }
\newcommand{\ee}{ \end{eqnarray} }

\newcommand{\nodata}{...}
\newcommand{\pasa}{PASA}
\newcommand{\apj}{ApJ}%                                         % Journal abbreviations
\newcommand{\apjs}{ApJS}
\newcommand{\apjl}{ApJL}
\newcommand{\aap}{A{\&}A}
\newcommand{\aaps}{A{\&}AS}
\newcommand{\mnras}{MNRAS}
\newcommand{\aj}{AJ}
\newcommand{\araa}{ARAA}
\newcommand{\pasp}{PASP}
\newcommand{\aplett}{ApL}
\newcommand{\nat}{Nature}
\newcommand{\iaucirc}{IAU circ.}
%\newcommand{\be}{\begin{align}}
%\newcommand{\ee}{\end{align}}

%\tightenlines
%\singlespace
%


%\renewcommand{\normalsize}{\fontsize{12}{12}\selectfont}
%\renewcommand{\labelfonts}{\fontsize{11}{11}\selectfont}
%\newcommand{\captionfonts}{\fontsize{10}{10}\selectfont}

%\newcommand{\captionfonts}{\fontsize{11}{11}\selectfont}

%\renewcommand{\figurename}{\fontsize{11}{11}\selectfont}


%\makeatletter  % Allow the use of @ in command names
%\long\def\@makecaption#1#2{%
%  \vskip\abovecaptionskip
%  \sbox\@tempboxa{{\captionfonts #1: #2}}%
%  \ifdim \wd\@tempboxa >\hsize
%    {\captionfonts #1: #2\par}
%  \else
%    \hbox to\hsize{\hfil\box\@tempboxa\hfil}%
%  \fi
%  \vskip\belowcaptionskip}
%\makeatother 
\renewcommand{\footnotesize}{\normalsize}
%\renewcommand{\footnotesize}{\fontsize{10}{10}\selectfont}


\title[Phase-coherent timing of the Vela pulsar]{Characterising the rotational irregularities of the Vela pulsar from  21~yr of phase-coherent timing}

\author[R.~M.~Shannon et al. ]{ R.~M.~Shannon\thanks{E-mail: ryan.shannon@csiro.au}$^{1,2}$,  L.~T. Lentati$^3$, M.~Kerr$^1$,  S.~Johnston$^1$, G.~Hobbs$^1$, \newauthor \& R.~N.~Manchester$^1$\\
$^1$ CSIRO Astronomy and Space Science, Australia Telescope National Facility, Box 76 Epping, NSW, 1710, Australia. \\
$^2$ International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia. \\
$^3$ Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK.}



\begin{document}

\date{Accepted XX ; in original form \today}

\pagerange{1--10} \pubyear{2015}

\maketitle

\begin{abstract}

Pulsars show two classes of rotational irregularities that can be used to understand neutron-star interiors and magnetospheres: glitches and timing noise.
Here we present an  analysis of the Vela pulsar spanning nearly $21$~yr of observation and including $8$~glitches.
We identify  the relative pulse number of all of the observations between glitches, with the only  pulse-number ambiguities existing over glitch events.    
We use the  phase coherence of the timing solution to simultaneously  model the timing noise and glitches in a Bayesian framework, allowing us to select preferred models  for both. 
We find the glitches can be described using only permanent and transient changes in spin frequency, i.e., no step changes in frequency derivative.
For all of the glitches, we only need two exponentially decaying changes in spin frequency to model the transient components.
   In contrast to previous studies, we find that  the dominant transient components  decay on a common  $\approx 1300$~d time scale, and that a larger fraction ( $\gtrsim 25\%$) of  glitch amplitudes are associated with these  transient components.  
We also detect shorter-duration transient components of $\approx 25$~d, as previously observed, but are limited in sensitivity to events with shorter durations by the cadence of our observations.
 The timing noise is well described by a steep power-law process that is independent of the glitches and  subdominant to the glitch recovery. The braking index is constrained to be $< 8$  with 95\% confidence. 
This methodology can be used to robustly measure the properties of glitches and timing noise in other pulsars.

\end{abstract}



\begin{keywords}
 pulsars: general -- pulsars:~specific (PSR~B0833$-$45)  -- stars: neutron
\end{keywords}

\section{Introduction}

Pulsars are celebrated for the predictability of the arrival times of their pulses.
The power of the pulsar-timing method is realised when a phase-coherent model of the  pulse times of arrival  (TOAs) is  achieved;  that is, when there is a solution that unambiguously accounts for every rotation of the pulsar.
For example, in an analysis of $25$~yr of observations of the first millisecond pulsar, PSR~B1937$+$21, the last observed TOA in $2010$ is 518,595,058,289~periods after   first TOA in $1986$ \cite[][]{2013ApJ...766....5S}.
% and in $22$~yr of observation of the eclipsing binary pulsar PSR~B1259$-$63 the observations span 15,225,402,818 periods. 
Phase coherence enables the detection of subtle effects that only  slightly modify the pulse frequency to accumulate in the arriving phase of pulses.  These effects include variations in the orbit of binary pulsars associated with general-relativistic effects and potentially the  passage of gravitational waves
with frequencies in the nanohertz range.
However, nearly every pulsar shows evidence for intrinsic spin irregularities which also alter the TOAs.  Irregularities are phenomenologically bifurcated into two forms, timing noise and glitches, with glitches more common in younger pulsars. 

Timing noise manifests as a red-noise (time-correlated) process in  the TOAs, and is  typically described by a wide-sense stationary stochastic process \cite[][]{1975ApJS...29..453G},  modelled as either  random walks in the pulsar spin parameters, or  with a  more general  power spectrum. 
The origin of this noise is unclear.    While some of it may be caused by variable torques associated with changes in the pulsar magnetosphere state \cite[][]{2010Sci...329..408L},  a significant fraction is likely associated with rotational irregularities interior to the star \cite[][]{2014MNRAS.437...21M}.
The properties of timing noise vary markedly across the pulsar population, with its  strength depending on the pulsar spin frequency $\nu$ and  frequency derivative $\dot{\nu}$  \cite[][]{sc2010}, and likely other factors because of  the large dispersion levels of timing noise between pulsars.
 For example, the most stable millisecond pulsar (MSP), PSR~J1909$-$3744 shows no evidence for instabilities with phase variations limited  $\lesssim 100$~ns  ($\approx 0.3\%$ of pulse phase) over  $11$~yr \cite[][]{2015Sci...349..1522S}.   %Another MSP, PSR~B1937$+$21 shows phase variations of $10\%$ of pulse phase over $24$~yr \cite[][]{2013ApJ...766....5S}.    
 In contrast, the  timing noise in young pulsars and magnetars can contribute many cycles of pulse phase on week to month time scales \cite[][]{2011ApJ...730...66L},  both making  it difficult to find a phase coherent solutions  in poorly sampled data and  presenting challenges to TOA-modelling algorithms. 
 While finding  the origin of timing noise is important for understanding neutron stars,   it is also necessary to account for timing noise as  part of a general timing model for the pulsar.  This is necessary to eliminate, or at least mitigate, bias in  the estimation of other parameters in the model \cite[][]{2011MNRAS.418..561C,2013MNRAS.428.1147V,2014MNRAS.437.3004L,2015arXiv150706982K}.  
 
 In addition to exhibiting timing noise,  pulsars can  experience glitch events \cite[e.g.,][]{2011MNRAS.414.1679E,2013MNRAS.429..688Y}, in which they are observed to  suddenly change  spin state, with the most significant component being increases
  in spin frequency that can exceed $1:10^5$ \cite[][]{2011ApJ...736L..31M}. 
  Glitches have been modelled with  permanent changes in $\nu$ and $\dot{\nu}$, as well as transient components, in which changes in $\nu$  are modelled to decay (typically exponentially) on  time scales
  $\tau$.  
  Occasionally multiple glitch-decay components (with different time scales) are invoked, particularly when  high cadence (daily or higher) observations \cite[e.g.,][]{2002ApJ...564L..85D}.  
   Glitch components are often identified and characterised by searching for variation of spin frequency and frequency derivative in subsets of the data \cite[][]{2015MNRAS.446..857L}. 
   This method is suboptimal because it does not utilise the phase coherence of the pulse arrival times. 
  The study of glitches is further complicated by the presence of timing noise. 
   
    There are two prevailing theoretical models for glitches.   
 In the first, glitches are associated with the transfer of angular momentum between the superfluid interior and solid crust of the neutron star \cite[][]{1975Natur.256...25A,1984ApJ...276..325A}.
 A superfluid component is  present in the core of the neutron star and a portion pervades the inner crust \cite[][]{1969Natur.224..673B}.
The quantised angular velocity vortices can pin on the nuclear lattice of the crust. 
      As the crust and the normal fluid component of the NS spin down because of electromagnetic braking, differential angular momentum is built up between the pinned vortices and the other components. 
   Eventually an external trigger or the Magnus force \cite[e.g.,][]{2009ApJ...700.1524M}  causes the vortices to unpin and transfer angular momentum to the crust, spinning the star up.  
 In the second model, glitches are associated with star-quakes in the solid, crystalline crust of the neutron star.  The quakes are the result of changes in the equilibrium configuration as the oblate star relaxes toward a spherical state as it slows down and cools \cite[][]{1969Natur.223..597R}.     
 The latter model is presently disfavoured for most pulsars because it cannot account for  the amplitude distribution and event rate observed in glitches \cite[][]{2015IJMPD..2430008H}, though it might be suitable for a few pulsars such as the Crab pulsar. 
Empirical evidence suggests that glitches can be modelled using a deterministic signal in the pulsar timing model \cite[][]{2006MNRAS.372.1549E} with only a few parameters.  
  
% As glitches are discrete events, they can be included as a deterministic process with the pulsar timing model \cite[][]{2006MNRAS.372.1549E} b 
Theoretical models  predict that  glitch decay time scales should be constant for individual pulsars if the underlying physics driving glitch recovery is a linear process and that the transient components should have comparable magnitudes to the permanent components \cite[][]{2015IJMPD..2430008H}.   However, the transient components hitherto  measured have a small contribution relative to the permanent component and have variable decay time scales \cite[][]{2011MNRAS.414.1679E}.

Many young pulsars also show evidence for  measurable braking that is attributed to pulsar spin down.    This is a deterministic process that  primarily manifests as a  second derivative of spin frequency $\ddot{\nu}$, and is parametrized by a braking index $n$ (where $\dot{\nu} \propto \nu^n$) .  For electromagnetic braking associated with a dipolar magnetic field $n=3$;  however measured braking indices often depart from this markedly, with contributions likely arising from  angular momentum loss from particle winds and free precession.   


Here we study the  timing properties of the Vela pulsar \cite[][]{1968Natur.220..340L}, a relatively  young  (characteristic age of $\tau_c \approx 11$~kyr) pulsar that shows both large levels of timing noise and glitches at quasi-regular intervals.  
The Vela pulsar is the brightest known pulsar at decimetre wavelengths, with  period-averaged flux density  of $\sim 1$~Jy.   
The first glitch in any pulsar  was identified in the Vela pulsar \cite[][]{1969Natur.222..228R,1969Natur.222..229R} and  $16$ subsequent glitches have been detected \cite[][]{2013MNRAS.429..688Y}. 
Timing noise analysis of the pulsar \cite[][]{1980ApJ...239..640C,1988ApJ...330..847C} has typically been constrained to intervals between glitches.
Glitch analyses similarly have  been rarely conducted including the effects of timing noise. 
Similarly, the observation of a very low braking index of $n=1.4 \pm 0.2$  reported for the pulsar  \cite[][]{1996Natur.381..497L} did not account for timing noise. 


To analyse both the glitches and the timing noise simultaneously, we use a timing solution for the Vela pulsar that spans $21$~yr of observation and $8$ glitch events, enabling us to examine the stationarity of the timing noise and robustly estimate glitch parameters. In section \ref{sec:obs}, we present the observational data used.  In section \ref{sec:model}, we discuss the timing analysis and Bayesian methodology employed.  In section \ref{sec:results}, we compare phenomenological models for the timing noise and the glitches and select a preferred model.  In section \ref{sec:discuss}, we discuss the implications of this model. In section \ref{sec:conclusions}, we give our conclusions. 

\section{Observations}\label{sec:obs}

Our highly heterogeneous  data set comprises $1231$ TOAs obtained with the 64-metre Parkes radio telescope between $1992$~December~12 and $2014$~January~14. 
 While observations were conducted at frequencies between $0.4$ and $23$~GHz, most were made at a central frequency of $\approx 1.4$~GHz. 
Prior to $2003$ observations were made with a series of analogue-filterbank and digital-autocorrelation spectrometers; these observations are described in detail in \cite{2000MNRAS.317..843W} and \cite{2013MNRAS.429..688Y}.
   Most recently, the pulsar has been observed with digital polyphase filterbank spectrometers as part of a programme to monitor pulsars of interest to the {\it Fermi} gamma-ray observatory \cite[][]{2010PASA...27...64W}.   These observations have monthly cadence with a central observing frequency close to $1.4$~GHz and semi-annual cadence with a dual-band system capable of observing simultaneously at  central frequencies of $0.73$ and $3.0$~GHz.    

The primary data  in this analysis are TOAs, formed by correlating  observations that have been averaged in  frequency, time and (where recorded) polarisation with a template, using the commonly applied Fourier phase gradient method, described in \cite{1992RSPTA.341..117T}, and implemented in the pulsar analysis code {\sc psrchive} \cite[][]{2004PASA...21..302H}.   
Templates were produced individually for each backend/observing-band combination using an analytic model fitted to the average profile from that combination.  
Offsets between the backends were included in the timing model, as discussed in Section \ref{sec:model}.
The cross-correlation method assumes that the data can be described by the template and additive white noise.
  For our observations, this is not the case.  
Distortions of the pulse profile, especially prevalent in older observations, are introduced both by the high flux density  (in excess of system equivalent flux density) of the pulsar, and the  large dispersion sweep of the pulsar relative to the  pulse phase and frequency resolution of the observations.

%Digitisation can cause apparent negative flux density at the leading and trailing edges of pulse profile \cite[][]{1998PASP..110.1467J}.  
Saturation of the amplifiers, other non-linear effects in the receiver and downconversion chain, and low-bit digitisation can lead to artefacts in the pulse profile, such as  apparent negative flux density on the leading and trailing edges of the pulse \cite[][]{1998ApJ...498..365J}. 
Older observations were recorded  with analogue-filterbank spectrometers with single-bit digitisers and were especially susceptible to these artefacts. 
Additionally, the pulse profile can be artificially broadened if the dispersive delay across an individual channel bandwidth is larger than the pulse-phase resolution of the observation. 
Given the relatively narrow pulse ($2.1$~ms),  and  relatively high dispersion measure ($68$\,pc\,cm$^{-3}$), older observations conducted with wide channels at low frequency show this type of broadening.
Even in more recent observations where instrumental effects are minimised, stochasticity in the pulse shape   introduces additional timing error  \cite[referred to as pulse jitter, ][]{1985ApJS...59..343C} that limits the timing precision of the observations. 
The effects of all of these distortions are secondary to TOA variations induced by timing noise and glitch events.
 It is however necessary to account for these effects in the analysis,  in particular  when modelling transient glitch components in our sparsely sampled data set.
While we do not account for them while measuring TOAs \cite[][]{2015MNRAS.454.1058L}, we account for their effect in the pulsar timing model, as discussed in the next section. 




 
\section{Timing Analysis}\label{sec:model}


\begin{figure}
\begin{center} 
\begin{tabular}{c}
%\includegraphics[angle=-90,scale=0.5]{1937_mem.eps}
\includegraphics[scale=0.5]{residuals_all.eps}  \\
%\includegraphics[scale=0.5]{residuals_noise.eps}  \\
%\includegraphics[scale=0.5]{residuals_nb.eps} \\
%\includegraphics[scale=0.5]{residuals_brake.eps} \\
%\includegraphics[scale=0.5]{residuals_pl.eps} \\
\includegraphics[scale=0.5]{residuals_optimal.eps} \\
\includegraphics[scale=0.5]{glitch_model.eps}  \\
\includegraphics[scale=0.5]{residual_white.eps}\\
\includegraphics[scale=0.5]{dmvar.eps} \\
%\includegraphics[scale=0.5]{diff.eps} 
\end{tabular}
\caption{  \label{fig:residual_plot}  \footnotesize Residual arrival times for maximum-likelihood models  of the Vela pulsar, measured in units of time $\Delta t$ and cycles of phase $\Delta P$.   The stars show the epochs of the glitches.   {\em a}:  Only fitting for the spin frequency and frequency derivative. 
% Panel b:  Fitting for glitches in the presence of band-limited timing noise.     
{\em b}.    Fitting for the glitches but assuming power-law noise.   
{\em c}:  Modelled glitch signal from solution presented in panel b.    
{\em d}:  Whitened residuals for maximum-likelihood solution.
{\em e}:  Dispersion measure variations for the maximum-likelihood solution (solid line) .  The dashed line shows the $1\sigma$  uncertainties on the realisation.}  
\end{center} 
\end{figure}

The presence of strong timing noise and glitch events  make it difficult to produce phase-connected solutions  over long data spans for young, energetic pulsars like the Vela pulsar.  
As the data spans increase, the amplitude of the timing-noise signal increases rapidly \cite[with the timing noise having a power spectral density $P_r(f) \propto f^{-5 \pm 2}$  across the population, ][]{sc2010}, and relative to a spin period and period derivative at some fiducial epoch, the arrival times diverge.     Even if an initial phase-connected solution exists, it typically  has hitherto been difficult to fit the solution because small changes in model parameters  can change the residual arrival times by more than one cycle of pulse phase.

%A variety of methods have been used to produce phase connected timing solutions and to mitigate timing noise. 
%One commonly used method is to add an increasing number of polynomials or a Fourier basis of sinusoids to the timing model. 
%This can work for datasets with relatively weak timing noise. 
%Solutions with large numbers of sinusoids or polynomical  do not extrapolate well, limiting the effectiveness of the technique for long datasets. 

%Another method is to add phase wraps to the timing residuals.
%This can be done through visual inspection of the residual time series,  or automatically.
%Automatic phase-tracking capabilities are present in the  {\sc tempo} and  {\sc tempo2}  packages.   In this case the codes dynamically add integer phase offsets are automatically included in the timing model to locally minimize the derivative of the residuals.
%This requires that over short intervals (where short contains greater than a few TOAs), the variations to be smaller than a cycle of phase.  Over long data spans the timing noise can introduce many cycles of phase relative to a fiducial solution.

Within the pulsar-timing code {\sc tempo2}  \cite[][]{2006MNRAS.369..655H}, it is now possible to use the relative pulse numbers of the TOAs as a reference when  measuring the goodness of fit or the likelihood of a model. 
% The relative pulse numbers are independent of the model, so can be found for an unphysical model that correctly accounts for the phase variations.  
If the level of timing noise is large so that phase connecting the TOAs over the entire data set is difficult or impossible,   the relative pulse numbers for the entire data set can be determined after producing phase-connected time solutions for  subsets  of observations. 
 The overlapping intervals are chosen to be sufficiently long that the relative  pulse numbering can be checked for agreement 
% Intervals are chosen to cover epochs  where phase connection in unambiguous:  the noise in the TOA variations can be modelled with many fewer parameters than number of observations. 
 These solutions do not need to be physical, so timing noise can be whitened using sinusoids \cite[][]{2005MNRAS.360..974H} or derivatives of the pulse spin frequency.
The consistency of the solutions can be checked using TOAs common to different intervals, i.e.,  relative pulse numbers should be the same for overlapping observations.


Ambiguities in the pulse numbering potentially exist immediately following glitch events.
If the epoch of the glitch is poorly constrained, the glitch amplitude is large and the observation cadence is poor around the glitch epoch,  there may be uncertainty in the number of rotations of the pulsar between the glitch and the first post-glitch observation.  
%This will result in an offset in the pulse numbering if the glitch epoch is  mis-estimated by an amount 
Possible glitch epochs,  corresponding to changes in the pulse numbering by one unit, are separated by $\Delta \nu_g^{-1}$, where $\Delta \nu_g$ is the change in spin frequency associated with the glitch. 
%As it is highly unlikely that glitches introduce a phase offset in the neutron star rotation, any excess phase can be used to more precisely estimate the glitch epoch.
%  However, the phase offset is known only modulo the pulse rotation so uncertainties in  pulse numbering adding integer cycles of phase at the glitch epoch.
For most of the glitches in our data set,  the glitch epochs have been previously reported  to sufficient accuracy such that we can unambiguously identify the pulse rotations  through the glitch event.
For these glitches, the maximum-likelihood glitch phase parameter is found to be within  $< 0.1$ cycles of zero  offset. 
If the glitch epoch were perfectly determined, the glitch phase parameter would be zero.  
%If  the pulse numbering were incorrect (and the epoch correct), the glitch phase parameter would be an integer number of rotations of the pulsar.
  However, a few of the glitches epochs are poorly constrained. 
In this case we are unsure of the relative pulse numbering through the glitches.   
%In this case it is possible that the pulse numbering could be incorrect at the glitch epoch.
%For the glitch at MJD a55428, the closest phase offset was $0.16$ cycles, with possible glitch epochs of MJD $55427.9$ or $55428.5$. 
We tested the effects of having incorrect pulse numbers through the glitches, and found that no parameters in our analysis (beyond the glitch phase) were significantly affected, because the glitch epochs are constrained to a time much shorter than the fastest varying process in our model. 
This is to be expected, because we have included the glitch phase parameter in our timing analysis. 
Because we analytically marginalise over the glitch phase parameter, we are marginalising over the pulse-number uncertainty that exists at glitch epochs.  

Between glitches, we confirmed our pulse numbering is correct by whitening the entire dataset with a series of sinusoids \cite[][]{2005MNRAS.360..974H}.
Incorrect pulse numbering would result in TOAs that  have a random phase offset from the other arrival times, and  residuals  distributed across a full cycle of pulse phase in the whitened solution. 
%the sinusoidal time series or the residuals.   
We used a series of $40$ sinusoids with periods ranging from $23.2$ to $0.58$~yr.  This was sufficient to model most of the red timing noise with the residual TOAs constrained to  $\approx 0.05$ of pulse phase.  
We found no outlying points.  
In contrast, when we purposely introduced pulse numbering errors,  outlying TOAs were easily identified. 

%We cannot exclude ambiguities of integer numbers of rotations of the the pulsar 
%We conclude that our observations span  7,459,866,091 rotations of pulsar. 

Using this timing solution, we can directly calculate the absolute residual pulse phase of the TOAs  to the model, even if the difference between the model and data are $\gg 1$ cycle of phase.  
We used the code {\sc temponest} \cite[][]{2014MNRAS.437.3004L} to construct a complete timing model from the phase-connected solution.
The Bayesian framework implemented in {\sc temponest} enables us to simultaneously model stochastic parameters (e.g., timing noise) and deterministic parameters (e.g., glitch parameters) of interest, and marginalise over nuisance parameters of no interest to the analysis. 
For example we were able to search over a wide range of glitch decay times and glitch epochs, which is not possible with the fitting algorithm internal to {\sc tempo2}. 
 
 This  approach also enabled us  to select favoured models through the use of  Bayesian evidence (the integral of the likelihood over the parameter space weighted by the prior), which can be used to calculate the relative probabilities of  different models.
  {\sc temponest} was recently enhanced to use the {\sc PolyChord} algorithm \cite[][]{2015arXiv150201856H}  to sample posterior distributions and calculate the evidence, enabling more efficient searches and more robust calculation of evidences in high-dimensional parameter spaces than the alternative {\sc multinest} algorithm \cite[][]{2009MNRAS.398.1601F}. 

We also had to modify the code to incorporate higher precision ($128$-bit and $256$-bit) floating-point arithmetic\footnote{We implemented the arithmetic using the {\sc mlapack} library: {\tt http://mplapack.sourceforge.net/}.},   because of the  dynamic range required to model the arrival times.  The highest precision measurements have individual TOA errors (after accounting for systematics and pulse shape variations) of a few $\mu$s, while the plausible red-noise variations over the data span exceed $1000$~s.  This corresponds to a dynamic range  of $10^{18}$  for the  noise covariance.   
Matrices with this dynamic range need to be inverted as part of our analysis.  For our analysis, we found that $128$-bit precision was sufficient, and the $256$-bit precision was notably slower. 

Our model for the TOAs includes deterministic variations \cite[][]{2006MNRAS.372.1549E}  to account for the pulsar spin-down and astrometric terms.  There are  position, proper motion, and parallax measurements from long-baseline radio interferometry that exceed  the precision we can obtain through pulsar timing by factors of $\gg 10$ \cite[][]{2003ApJ...596.1137D}  and are consistent with our measurements. 
We therefore held the value fixed at the interferometrically determined position  in our analysis\footnote{While formally the VLBI astrometry should be incorporated into the prior in our analysis,  it is acceptable to fix the position at the VLBI position because it is of much greater precision, i.e., the Gaussian and delta-function priors are equivalent relative to the precision that can be measured from the data.}.
We analytically marginalised over deterministic parameters that are linear (or linearisable) in the timing model \cite[][]{2009MNRAS.395.1005V}, greatly reducing the time required to estimate non-linear and stochastic components of the timing model. 

We  also included terms to account for  stochastic time-independent (white-noise) and time-correlated (red-noise) contributions to TOAs. 
In addition to the white noise associated with the formal TOA uncertainty, terms are included to account for intrinsic shape variations and instrumental distortions, and other instrumental and astrophysical effects that are temporally uncorrelated between observations.  
We model the white noise by adjusting the uncertainty on individual TOAs to be
\be
\label{eqn:whitenoise}
\sigma^2 = F \sigma_r^2 +  \sigma_Q^2  %+ \sigma_J^2(T_{\rm obs}) ,
\ee
where $\sigma_r^2$ is the formal uncertainty derived from TOA fitting.
%, $\sigma_J$ is the contribution from single-pulse stochasticity. 
The factor $F$  (often referred to as EFAC\footnote{We have used the {\sc temponest} convention for defining EFAC and EQUAD which differs from the {\sc tempo2} definition.})  modifies $\sigma_R$ to account for instrumental distortions. 
 The term $\sigma_Q$ (often referred to as EQUAD)   accounts for additional observation-independent uncertainties.  
Both $F$ and $\sigma_Q$ are defined independently for each band-backend system.

%Assuming the shape variations are independent pulse-to-pulse, as the total number of averaged pulses   $N_p$  averaged increases , the jitter noise can be modelled to be $\sigma_J \propto N_p^{-1/2} \propto T_{\rm obs} ^{-1/2}$, where $T_{\rm obs}$ is the integration time \cite[][]{2014MNRAS.443.1463S}.   
%  Using high time resolution data we measured  $\sigma_J  = 25\, \mu s~(T_{\rm obs}/30~{\rm s})^{-1/2}$.    

Red-noise contributions to the TOAs include radio frequency-independent timing noise and dispersion-measure variations.  
Their contribution is often described by assuming the  amplitude of the fluctuations can be described using a power spectrum, which is suitable for wide-sense stationary processes.
The simplest model for timing noise we considered was a power-law power spectrum characterised with a spectral index $\beta$ and an amplitude $A$:
\be
\label{eqn:tnpl}
P_{r,{\rm PL}}(f) =  A  \left(\frac{f}{f_{\rm yr}}\right)^{\beta},
\ee
where  $f_{\rm yr}$ is a frequency of $1$ cycle per year.


In the second model, the power-law is modified  to include a spectral flattening
\be
\label{eqn:tnbl}
P_{r, {\rm BL}}(f) = \frac{A (f_c/f_{\rm yr})^{-\beta} }{\left[1 + (f/f_c)^{-\beta/2}\right]^2},
\ee
where $A$ again is the amplitude, and $f_c$ is a corner frequency.  With this definition, $P_r(f)$ is a power law when  $f \gg f_c$, and  when $f \ll f_c$, $P_r(f)$ is constant.  
This model is motivated by observations of non-power-law and apparently quasiperiodic timing  noise observed in many pulsars \cite[][]{2010MNRAS.402.1027H}.

We modelled the noise in the time domain using a harmonically related series of sinusoids constrained to have spectral density parametrized using the forms in Equations (\ref{eqn:tnpl}) or (\ref{eqn:tnbl}). 
It is necessary to mitigate spectral leakage of low-frequency power when modelling stochastic processes with Fourier series. 
However if the timing noise is band-limited or the power-law spectral index is relatively flat ($\beta >-6$), it is sufficient to start the series at $f = 1/T_{\rm span}$, because the inclusion of  $\nu$ and $\dot{\nu}$   act as a pre-whitening filter on the data set  \cite[][]{1984JApA....5..369B}.   
However, for young pulsars, timing noise can be very steep, e.g. PSR B1259$-$63 was measured to have $\beta \approx -9$ \cite[][]{2014MNRAS.437.3255S}.  
To mitigate spectral leakage in  this case, we included low frequencies $f < 1/T$ to model the lowest frequency timing noise at logarithmically spaced intervals, as described in  \cite{2015MNRAS.446.1170V}.   


We also searched for dispersion measure (DM) variations. 
Following a technique outlined in \cite{2015arXiv150403692L},  we include two components in the model.
The first is a stochastic component, in which the DM fluctuations (distinct from the  TOA fluctuations) are modelled as a time series constrained to have a power-law power spectrum, \be
\label{eqn:dmpower}
P_{\rm DM}(f) = D \left( f/f_{\rm yr} \right)^{\gamma},
\ee
where $D$ is the amplitude of the DM variations and $\gamma$ is the spectral index of the assumed power law.  
The second component is a quadratic polynomial in DM that accounts for secular trends observed in some pulsars \cite[][]{2013MNRAS.429.2161K} but also acts to mitigate spectral leakage of the first component. 
Together they account for stochastic variations  associated with the the turbulent ionised interstellar medium,  but the polynomial also models linear trends in DM that are also observed  by \cite{2013MNRAS.429.2161K}.
In the TOAs, the fluctuations induced by the DM variations  scale by the inverse-square of the observing frequency, enabling DM variations and timing noise to be distinguished. 


We also searched for pulsar braking,  a deterministic process that primarily induces low-frequency signals  in TOAs. 
The pulsar rotation rate  is expected to decelerate as the pulsar spins down because 
of the associated decrease in the magnetic torque.  
%For a rotating dipole, both the energy radiated and angular momentum lost is proportional to $\nu$.
 As a pulsar slows down, the rate of deceleration decreases, resulting in a positive second derivative.    
In general the braking can be approximated by a braking index $n$ which modulates the secular spin evolution of the pulsar:
\be
\dot{\nu} = K \nu^n.
\ee
 In  standard magnetic-dipole braking,  the pulsar  magnetic field strength and magnetic-dipole inclination angle are assumed to be constant and $n=3$.    Departures from this value are interpreted as being associated with magnetic field evolution, changes in the spin-magnetic axis orientation, or the effects of a pulsar wind. 

For young pulsars this braking  index induces a measurable second derivative of the spin frequency $\ddot{\nu}_b$,
\be
\label{eqn:brake_nudot}
\ddot{\nu}_{b} = n \frac{\dot{\nu}^2}{\nu}, 
\ee
 and, potentially a third frequency derivative $\dddot{\nu}$,   
\be
\label{eqn:brake_nuddot}
\dddot{\nu}_{b} = n (2n-1) \frac{\dot{\nu}^3}{\nu^2}.
\ee

   Measured values of $n$  disagree markedly with  that predicted from magnetic-dipole braking.  
   For the youngest pulsars  ($\tau_c < 10$~kyr) , the braking index generally has a small positive value ($1 < n < 5$).  
   For other young pulsars ($10 < \tau_c < 10^6$~kyr), the braking index is positive, but  typically large ($n \gtrsim 10$),  suggesting that $n$  is not  associated with braking but another  process,  such as recovery from previous glitches  \cite[][]{1999MNRAS.306L..50J,2001MNRAS.328..855W,2010MNRAS.402.1027H}.       
   For the oldest pulsars ($\tau_c > 10^6$~kyr), the braking index is measured to have  both  positive and negative values, again suggesting it is not associated with pulsar braking but another process.  In these pulsars,  $n$ is likely being masked by  red noise.  
   In  previous studies of braking indices, no attempt was made to model the timing noise simultaneously to estimating  $n$. 
    We searched for pulsar braking by including it as a parameter in our timing model and adding in its  contribution to  $\ddot{\nu}$ and $\dddot{\nu}$ as described in Equations (\ref{eqn:brake_nudot}) and  (\ref{eqn:brake_nuddot}). 











%In Figure \ref{fig:residual_plot}, we show the residuals for three different plo










%\begin{itemize}



\section{Results} \label{sec:results}  


In panel  {\em a} of Figure \ref{fig:residual_plot}, we show the residual arrival times, only fitting for the pulsar spin frequency and frequency derivative and therefore excluding terms that account for glitch events or the timing noise. 
In the absence of these terms, the residual arrival times $\Delta t$ show peak-to-peak variations of  $100$\,s, which is  $>1000$ cycles of pulse phase. 
At the epochs of the glitches (indicated by the stars at the bottom of the panel), the residuals show a discontinuous change in slope,  associated with sudden change in spin frequency. 
Between glitches, the residuals have positive curvature.
This process has been interpreted as a change in $\dot{\nu}$ and has been widely reported   \cite[][]{2011MNRAS.414.1679E,2013MNRAS.429..688Y}.  
Interestingly, the process apparently causes  the pulsar to return to comparable phase at each epoch, after accounting for the uncertainty in  $\nu$ and $\dot{\nu}$ through marginalisation.

%With the phase connected solution, we observe that  on only  connection in spin frequency lost in the timing noise and recovered in the glitch, but an apparent recovery of pulse phase as well.

\subsection{Models for glitches}

We  compare phenomenologically different scenarios for the glitches and timing noise, motivated by previously applied models and the observed residuals presented in panel {\em a} of Figure \ref{fig:residual_plot}:  

{\em  Permanent and  transient changes in $\nu$:}  As a minimal model, we assume that the glitches can be described by a permanent ($\Delta \nu_p$) and transient  ($\Delta \nu_t$)  change in the spin frequency. 
The permanent change introduces a change in the pulse phase at time $t>t_g$
\be
\Delta \phi(t) = \Delta \nu_p  (t-t_g),
\ee
where $t_g$ is the glitch epoch.
The transient component is assumed to exponentially  decay on a time scale $\tau$, so its contribution to pulse phase is
\be
%[1 ? e?(t?tg)/?d ]?d?d
\Delta \phi(t) =  \nu_d \tau \left[ 1-\exp\left( \frac{t-t_g}{\tau} \right) \right].
\ee

 We search over all possible values for the component amplitudes (both positive and negative) and decay times from  $1 < \tau < 10^4$~d.

 {\em Two decay time scales:}  We include an additional series of glitch decay times, to search for both long time-scale decays $\Delta \nu_\ell$ $\gg  80$~d and short time-scale glitch decay components $\Delta \nu_s$, which, respectively have decay time-scales $\tau_\ell$ and $\tau_s$. 

 {\em Common glitch-decay time scales:} In addition to assuming that the glitches can be described by permanent and transient changes in spin frequency, we assume that the transient components (long and short) all have the same time scales $\tau_s$ and $\tau_\ell$, respectively.  
 
 {\em   Permanent changes in $\dot{\nu}$:}  We assume that there are discrete changes in the spin frequency derivative $\Delta \dot{\nu}_p$ at the epoch of each glitch, as has been previously reported for the pulsar \cite[][]{2013MNRAS.429..688Y} and commonly modelled. 
The change in pulse phase is, for $t>t_g$,
\be
\Delta \phi (t)  = \frac{1}{2}  \Delta \dot{\nu}_p \left( t-t_g \right)^2
\ee
%

 {\em  Decay from glitch prior to our observations}:  We include exponential recovery associated with a glitch that occurred prior to the observation.   This would be recovery of a glitch occurring  $\approx 400$~d prior to our first observation \cite[][]{1991IAUC.5311....3F}, or to model long time-scale
 glitch recovery, which, as noted above, is hypothesised to be the source of timing noise in young pulsars.
 
{\em  Additional small glitches:}    In addition to the reported glitches, we search  the dataset for additional glitches.  This enables us to assess the apparent dichotomy between glitches and timing noise and search for micro-glitches that have been previously suggested \cite[][]{1988ApJ...330..847C}. 

In summary, the general  model for a glitch with permanent changes  $\Delta \nu_p$, and $\Delta \dot{\nu}_p$, and a single transient component $\Delta \nu_t$ that decays on a timescale $\tau$ is
\be
\Delta \phi(t) =  \Delta \Phi  +  \Delta \nu_p  (t-t_g), + \frac{1}{2}  \Delta \dot{\nu}_p \left( t-t_g \right)^2  \nonumber \\
  +\nu_d \tau \left[ 1-\exp\left( \frac{t-t_g}{\tau} \right) \right].
\ee
The phase $\Delta \Phi$ accounts for errors in  $t_g$  (if it is held fixed)  or errors in the pulse numbering at the glitch event. 

\subsection{Models for timing noise}

We consider a smaller number of models for the timing noise: 

{\em    No timing noise:}  We first consider a model containing no timing noise.  This model therefore assumes that the entirety of the time-correlated signal in the TOAs is associated with the glitches.

{\em    Power-law timing noise:}  We assume that the timing noise can be described by a wide-sense-stationary power-law process (Equation \ref{eqn:tnpl}) that is observed in many young pulsars. 

{\em   Band-limited timing noise:}  We assume that the timing noise is band limited and can be described by Equation (\ref{eqn:tnbl}), which would be the case if it was caused by state changing or a similar process. 

{\em    Non-stationary timing noise:}   We assume that in each inter-glitch period, the timing noise can be described by an independent power-law red-noise process. This model enables us to determine if the timing noise and glitch activity are correlated. 

{\em  Braking index:}  We additionally search for pulsar braking index $n$.
  
\subsection{Supported models}

For the models listed in the previous sections,
we calculated the Bayesian evidence  over the distribution of posterior parameters of interest while marginalising over the remaining nuisance parameters.

    In Table \ref{tab:evidence} we show a hierarchy of models that provided significantly improved evidence, culminating in the preferred model. 
In general, models that provide improvements  of $\Delta \log E >  3$ or have fewer parameters but comparable values of evidence are favoured. 

   
   
   \begin{table*}
\begin{centering}
\caption{Model comparison}
\label{tab:evidence}
\begin{tabular}{cllrr}
\hline
\hline
  \multicolumn{1}{c}{Model} & \multicolumn{1}{c}{TN Model}  & \multicolumn{1}{c}{Glitch Model} &  \multicolumn{1}{c}{log(E)} & \multicolumn{1}{c}{ $\Delta$log(E) }   \\
\hline
1 &  No red noise & $\Delta \nu_p + \Delta \nu_t$ 	&4559.4&	-5394.0\\
2 & Red  & $\Delta \nu_p +\Delta \nu_t$	&	9953.4	&0.0	\\
3 & Red  $+$ DM 	& $\Delta \nu_p + \Delta \nu_t$	&10057.5	&104.1	\\
4  & Red $+$ DM & $\Delta \nu_p+ \Delta \dot{\nu}$ & 10018.4 & 65.0 \\ 
5 & Red $+$ DM	&(3) $+$ prior glitch 	&10099.7	&146.3	\\
6 & Red $+$ DM	& (4) $+$ $\Delta \nu_s$  + $\Delta \nu_\ell$	&10143.9	&190.5	\\
7 & Red $+$ DM	&(5)  $+$ common glitch time scales &10144.7&	191.3 \\ 
\hline
\end{tabular}
\end{centering}
\begin{flushleft}
Models for timing of the Vela pulsar and their evidence.  For each model we show the timing noise model, and the glitch models, and the total evidence, and relative evidence.   .  Timing noise model:   RN:  Red noise;  DM:  dispersion-measure variations.   Glitch models:   $\Delta \nu_p$: permanent changes in spin frequency;  $\Delta \nu_t$: transient change in spin frequency (if only one transient component included), modelled by as an exponential variation with an amplitude and a time scale $\tau$ ; $\Delta \nu_\ell$:  long duration transient component (if two component are included), with time scale $\tau_\ell$;  $\Delta \nu_s$:  short-duration transient component (if two components included) with time scale $\tau_s$; $\Delta \dot{\nu}$:  permanent changes in spin-frequency derivative. 

\end{flushleft}
\end{table*}
     
{\em 1.  Power-law  stationary timing noise.}     In panel {\em b} of Figure \ref{fig:residual_plot} we show the residuals from the maximum-likelihood model  accounting for the glitches,  modelling the timing noise to be a power-law red noise process and marginalising over the unknown $\nu$ and $\dot{\nu}$. 
The dominant signal in the residuals is red noise, which induces peak-to-peak variations in the residuals of $12$~s. 
The residuals are well described by a power-law red noise process with a spectral index of $\beta \approx -6$. 
We measured a consistent  amplitude and spectral shape for the timing noise  using the spectral-modelling algorithm presented in \cite{2011MNRAS.418..561C}, after fixing the glitch parameters at their maximum-likelihood values.  We find no evidence that the level or amplitude of timing noise varies between glitches. 
 
 
{\em 2. Permanent frequency changes and two decaying component with  time-scales  (a short time scale of $\tau_s \approx 25$~d and a long time scale of $\tau_\ell \approx 1300$~d) that are common amongst the glitches.  }  The amplitude of the glitch components can be found in Table \ref{tab:glitch_params}.
For glitches $1$, $2$, $3$, and $5$, our data could not constrain the amplitudes of the short duration components $\Delta \nu_s$, likely because of poor sampling at these epochs.    
 In panel {\em c} of Figure \ref{fig:residual_plot}, we show the maximum-likelihood model of the glitch signal, after marginalising over the pulsar spin-down.     
The signal closely resembles the difference of panels {\em a} and {\em b}.   Once accounting for the timing noise (and marginalising over the uncertain $\nu$, $\dot{\nu}$,  DM, and instrumental jumps) the glitch signal appears to be relatively stationary in the residuals.     The maximum-likelihiood parameters for the glitches are displayed in Table \ref{tab:glitch_params}. 
We find no evidence for any variation in the glitch-decay time scales.
When modelled individually the long glitch decay time-scales show consistent posterior distributions. 
%While for most glitches we were able to measure the amplitudes of the short decay time scales, for a few the distributions were not constrained because the observing cadence around the glitches was poor.  
 %In Figure \ref{fig:post_plot}, we show the posterior distributions for the long-term glitch decay time scales ,fitted independently for each glitch.  
% While the distributions are broad, they show consistent time scales. 
 These long time-scale  components have not previously been observed and have been  modelled  previously as a combination of  inter-glitch timing noise and  permanent changes in slowdown ($\Delta \dot{\nu}$).
 Models containing only short time scale glitches and permanent changes in $\dot{\nu}$ (Model 4 in Table \ref{tab:evidence}) are significantly disfavoured over models containing only transient components (Model 5). 
We also find evidence for the decay from a previous glitch (Model 5 in Table \ref{tab:evidence}). 
 This term accounts for latent exponential decay from a glitch  $400$~d prior to our observations.    
In contrast to previous analyses, we find  a significant component of each glitch event is  associated with a transient component, with  between $ 30\%$ and $80\%$ of the total change in spin frequency  associated with a transient component. 


{\em 3.  No small glitches}. 
  In addition to modelling the $8$ large glitches in our data set, we searched for additional small glitches.    
We place limits on the amplitudes of other glitches of $\approx 10^{-7}$ averaged over the observing span.
Lower amplitude events are covariant with the timing noise.
% model. 

{\em 4.   No permanent changes in spin frequency derivative  $\Delta \dot{\nu}$.  }
This component instead is attributed to the long time-scale glitch decay.  The positive curvature in the residuals between epochs, usually modelled as a change in $\dot{\nu}$  is better modelled by long-term recovery. 
We verified that our methods could detect differences between permanent changes in $\dot{\nu}$  and long-term decays by simulating data sets that contained either significant $\Delta \dot{\nu}_p$ or long-term decays.      The simulated data set had either the maximum-likelihood values of  $\Delta \dot{\nu}$ (from a model that did not contain long-term decays) or those from the long time scale decays (as listed in Table \ref{tab:glitch_params}) obtained from our data set.  The simulated data sets also contained the maximum-likelihood realisation of red noise from our data set and simulated TOA uncertainties corrected for EQUAD and EFAC.  In both cases the correct model was selected with $\Delta \log E >   130$ (i.e., with probabilities of  $ > 1- e^{-130}$).  
%Evidences for GLF1 simulation 0.100153738566655957E+05-0.102941699E+005
%Evidences for decay simulations 10337.5-10222.53

In panel {\em d} of Figure \ref{fig:residual_plot}, we show the {\em whitened} maximum-likelihood residuals, after accounting for  glitches, timing noise, and dispersion-measure variations (discussed below).  The residuals show no systematic variations in arrival times, suggesting that the model is complete, and providing further confirmation that the pulse numbering is correct.    
 The plot also shows that the data quality has greatly improved in the most recent data, due almost entirely  to instrumental improvements. 
  

  



%\begin{figure}
%\begin{center} 
%\begin{tabular}{c}
%%\includegraphics[angle=-90,scale=0.5]{1937_mem.eps}
%\includegraphics[scale=0.5]{hist_tau.eps}  \\
%\end{tabular}
%\caption{  \label{fig:post_plot}  \footnotesize Posterior probability density $\rho$  for glitch decay times $\tau_g$.    We list the glitch number above the maximum-likelihood value for each curve.  The decay time for the glitch prior to our observations is identified as  {\em P} .   }  \end{center} 
%\end{figure}











\begin{table*}
\begin{centering}
\caption{Glitch parameters}
\label{tab:glitch_params}
\begin{tabular}{clrrrrrrc}
\hline
\hline
  & \multicolumn{1}{c}{MJD}  & \multicolumn{1}{c}{$\Delta \nu_p$} &  \multicolumn{1}{c}{$\Delta \nu_\ell$} & \multicolumn{1}{c}{$\Delta \nu_s$} &  \multicolumn{1}{c}{$\Delta \nu_g/ \nu$} &  \multicolumn{1}{c}{$Q$} & \multicolumn{1}{c}{$(\Delta \nu_g/\nu)_{\rm lit}$} &  \multicolumn{1}{c}{Ref.}\\
  & & ($\mu$Hz) &  \multicolumn{1}{c}{($\mu$Hz)} &\multicolumn{1}{c}{($\mu$Hz)}& \multicolumn{1}{c}{($10^{-6}$)}&& \multicolumn{1}{c}{($10^{-6}$)} &\\
\hline
P & 48457&	\nodata	&16(4)	&\nodata	&\nodata	&\nodata	&2.715(2)	&1\\
1 & 49559&	1.8(2)	&7.8(2)	&\nodata	&0.86(3)	&0.81(4)	&0.835(2)	&2\\
2& 49591 &	1.4(3)	&0.7(2)	&\nodata	&0.19(2)	&0.3(1)	&0.199(2) &	3\\
3 &50369	&     12.3(3)	&11.5(3)	&0.11(3)	&2.14(5)	&0.48(2)	&2.11(2)	&4\\
4& 51559 &	22.7(6)	&12.1(6)	&0.18(3)	&3.12(8)	&0.35(2)	&3.152(2)	&4\\
5 &53193&	10.6(3)	&12.3(3)	&0.1(1)	&2.06(4)	&0.54(2)	&2.100	&4\\
6 &53960&	19.5(4)	&9.2(3)	&0.3(1)	&2.59(5)	&0.83(3)	&2.62	&4\\
7 &55408&	9.1(3)	&12.1(3)	&\nodata	&1.89(4)	&0.33(2)	&1.94  &5\\
8 &56555 &	20.9(2)	&13.2(2)	&0.21(2)	&3.06(4)	&0.39(1)	&3.100	&6\\
\hline
\end{tabular}
\end{centering}
\begin{flushleft}
Maximum-likelihood glitch parameters.  For each glitch, we list the MJD, the permanent change in spin frequency $\Delta \nu_p$,  the long and short glitch-decay amplitudes (respectively,  $\Delta \nu_\ell$ and $\Delta  \nu_s$) , reported in absolute value.   For glitches $1$,$2$,$7$, our data could not significantly constrain the amplitudes of the short glitch recoveries.   We also show the total change in spin frequency at the glitch epoch $\Delta \nu_g$, which is the sum of the permanent and transient components,  relative to the pulsar spin frequency $\nu$.  The values in parentheses represent the nominal $1-\sigma$ uncertainties for our measurements and, where available, previous measurements.  $Q = (\Delta \nu_\ell + \Delta \nu_s)/(\Delta \nu_s + \Delta \nu_\ell + \Delta \nu_p)$ is the fraction of the glitch that is recovered. 
   Where available we also show previously reported total glitch-decay amplitudes $\Delta \nu_{g,{\rm lit}}$. 
   The references are (1) \cite{1991IAUC.5311....3F};  (2) \cite{1994IAUC.6038....2F}; (3)  \cite{1994IAUC.6064....2F} ;  (4) \cite{2013MNRAS.429..688Y} ;  (5) \cite{2010ATel.2768....1B}; and  (6) \cite{2013ATel.5406....1B}. 
\end{flushleft}
\end{table*}





%\section{Braking Index}





\subsection{Dispersion-measure variations}



As part of a larger study of dispersion-measure (DM) variations  young energetic pulsars, \cite{2013MNRAS.435.1610P} measured DM variations for the Vela pulsar, utilising some of the data presented here.  The DM variations  were modelled using a linearly interpolated time series described and as described in  \cite{2013MNRAS.429.2161K}.  
 We  searched for the amplitude and  power-law index associated with DM variations, as described above. 
 The results of our model for DM variations are presented in panel {\em e} of Figure \ref{fig:residual_plot} and  are consistent with those presented in  \cite{2013MNRAS.435.1610P}.
At $1.4$~GHz, the DM variations induce TOA fluctuations that are a factor of $10^4$ smaller than the timing noise.

\section{Discussion} \label{sec:discuss}

Compared to previous analyses of the glitches,  we find that a  much larger fraction of the changes in spin frequency is associated with transient decaying components. 
The dominant transient components are associated with a long time-scale recovery ($\tau \approx 1300$~d), with a time scale that is common to all glitches.   
Not surprisingly, the sum of our permanent and transient changes in spin frequency is comparable to previously published measurements of the permanent component as  displayed in Table \ref{tab:glitch_params}.
With only eight glitches in our data set, we have an insufficient sample size to determine if there is a correlation between the amplitude of the permanent and transient components to pulsar glitches. 
The two glitches observed on MJD 49559 and 49591 are significantly smaller in $\Delta \nu_p$ and $\Delta \nu_\ell$ than the other glitches but are also unusual because of their relative contemporaneity. 


Many first-principle models of glitches  predict that transient components of different glitches  should decay with the same time scale
\cite[e.g.,][]{2010MNRAS.409.1253V}, because the microphysics of the neutron star, which regulates glitch decay, does not change. 

%The decay time scales are often longer than the time between glitches.

The length of the decays suggests that the recovery  (and the glitches) could be associated with the non-linear regime in the vortex-creep model for glitches. 
It has been predicted that the transition from linear to non-linear creep regime for the Vela pulsars would occur on time-scales of $\approx 1000$~d, close to our glitch time-scale of $1300$~d \cite[][]{1993ApJ...409..345A}.
While the regime may be non-linear, the glitch recovery would remain linear because the perturbation is relatively weak \cite[][]{1989ApJ...346..823A}, with $\Delta \nu_g/|\dot{\nu}| \ll \tau_\ell$. 


Most of the inter-glitch TOA variations are associated with the long time-scale decay of the transient component (glitch recovery). 
This is consistent with observations of timing noise in other young pulsars  \cite[][]{2010MNRAS.402.1027H}, which show  $\ddot{\nu}>0 $, much like the Vela pulsar.
However, other pulsars show markedly different distributions of glitch waiting times.



We find that the timing noise can be described by a wide-sense stationary process. 
When we considered independent realisations of timing noise between glitches, we found that they all had consistent amplitudes and spectral indices.  
Furthermore, the evidence supported a single coherent process.
Once accounting for glitch recovery,  we exclude long time-scale  decay as being the origin of the timing noise.    The timing noise has a spectral index of $\beta = -6.0\pm0.5$~$(1\sigma)$, which is comparable to other young pulsars.  
    The stationarity of the timing noise suggests that  it is not related to the glitches. 
    Because the timing noise is subdominant and unaffected by the glitches, we have no evidence for any causal or correlated relationship between the two phenomena. 

We do not find any evidence for braking of the pulsar spin down and set an upper limit of $n<8$.  Our sensitivity to the braking index is limited by the large TOA variations induced by the glitch events, the presence of timing noise, and the covariance of braking index with these parameters. 
A previous measurement of braking index  relied only on measurements of $\nu$ and $\dot{\nu}$  at specific post-glitch epochs to measure the braking index     \cite[][]{1996Natur.381..497L}.
 Through $9$ glitches, they found measurements of $\dot{\nu}$  at a date  $\approx 150$~d after each glitch (and after apparent glitch recovery)  were consistent with a braking index of $1.4 \pm 0.2$.  We attempted to reproduce the results of \cite{1996Natur.381..497L} using our dataset, which span an independent set of glitches.
We whitened our dataset by including a series of sinusoids in our maximum-likelihood timing model (displayed in panel b of Figure \ref{fig:residual_plot}).   
We then used this model (which is the sum of the glitches, their recovery, and the timing noise) to to calculate $\dot{\nu}$ at epochs $150$~d after glitch events.
Like \cite{1996Natur.381..497L}, we find that there is a linear trend in ${\nu}$ at these epochs. 
However, we do not identify a linear trend in $\dot{\nu}$ at any epoch after the $8$~glitches in our observations. 
We derive an apparent braking index to be $0.0 \pm 0.9$ at epochs $120$~d after glitches. 
%These are an independent set of glitches to those presented in \cite[][]{1996Natur.381..497L}.   
We attribute the differences in the apparent braking index to the the noise realisations in the two independent data sets.
%These are an span an independent set of glitches to those presented in \cite[][]{1996Natur.381..497L}.



% Using observations $150$~d after these four glitches,  we measure  $ \ddot{\nu} 2.3 \pm 2 \times 10^{-23}$ Hz~s$^{-2}$.   There is large scatter in the measurements of post-glitch $\dot{\nu}$. 
 
Despite the high levels of rotational instability, the timing behaviour observed in Vela shares similarities to that observed in older pulsars which glitch less and have much lower levels of timing noise.  In a sample of $366$ pulsars monitored over $\gtrsim 30$~yr,   \cite{2010MNRAS.402.1027H} found a number of pulsars that showed cuspy timing events that occur quasiperiodically. These cuspy events represent increases in $\nu$,  followed by episodes of increased $\dot{\nu}$ that are interrupted (at the same residual phase) by another increase in $\nu$.  The best example of this is PSR~B1900$+$06 \cite[][]{2010MNRAS.402.1027H}.    If these cusps are interpreted as small glitches and the intervening episodes as glitch recovery, the rotational irregularities are analogous to what we observe in the Vela pulsar, but with a smaller magnitude. 

The robustness of the solution  presented here depends  on the underlying assumption that one of the model families considered is correct.
 The greatest uncertainty is  the assumption that the timing noise follows a power-law or broken power process.  
If the timing noise can be better constrained it may be possible to detect and characterise additional components or measure a significant braking index. 
These uncertainties affect both Bayesian and maximum-likelihood approaches. 
Searches for the most transient components are limited by the relatively poor (monthly) cadence of our observations.  

\section{Conclusions}\label{sec:conclusions}

We have presented a timing solution for the Vela pulsar that spans $\approx 21$~yr.    
The solution is nearly phase connected, with the only uncertainties associated with pulse numbering at glitch events. 
This solution, use of a full timing model, and Bayesian methodology has enabled the robust  parametrization  of the glitches and characterise the  spin noise. 

We have identified  dominant transient components to the glitches that decays on a common  $25$~d and $1300$~d time scales for all $8$ glitches in our dataset and  a subdominant steep red noise component.
 These methods can be applied to other young pulsars to identify long time-scale glitch decay components, characterise timing noise, and robustly measure braking indices.



\section*{Acknowledgments}

The Parkes radio telescope is part of the Australia Telescope National Facility, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO.    GH is the recipient of ARC Future Fellowship FT120100595.     LTL acknowledges support from  a junior research fellowship at Trinity Hall College, Cambridge University.  RMS acknowledges travel support for this work through a John Philip early-career research award from CSIRO. 

\bibliographystyle{mn2e}
%\bibliography{thesis}
\begin{thebibliography}{59}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi

\bibitem[{{Alpar} {et~al}\mbox{.}(1993){Alpar}, {Chau}, {Cheng}, \&
  {Pines}}]{1993ApJ...409..345A}
{Alpar} M.~A., {Chau} H.~F., {Cheng} K.~S., {Pines} D., 1993, \apj, 409, 345

\bibitem[{{Alpar}, {Cheng} \& {Pines}(1989){Alpar}, {Cheng}, \&
  {Pines}}]{1989ApJ...346..823A}
{Alpar} M.~A., {Cheng} K.~S., {Pines} D., 1989, \apj, 346, 823

\bibitem[{{Alpar} {et~al}\mbox{.}(1984){Alpar}, {Pines}, {Anderson}, \&
  {Shaham}}]{1984ApJ...276..325A}
{Alpar} M.~A., {Pines} D., {Anderson} P.~W., {Shaham} J., 1984, \apj, 276, 325

\bibitem[{{Anderson} \& {Itoh}(1975)}]{1975Natur.256...25A}
{Anderson} P.~W., {Itoh} N., 1975, \nat, 256, 25

\bibitem[{{Baym}, {Pethick} \& {Pines}(1969){Baym}, {Pethick}, \&
  {Pines}}]{1969Natur.224..673B}
{Baym} G., {Pethick} C., {Pines} D., 1969, \nat, 224, 673

\bibitem[{{Blandford}, {Romani} \& {Narayan}(1984){Blandford}, {Romani}, \&
  {Narayan}}]{1984JApA....5..369B}
{Blandford} R., {Romani} R.~W., {Narayan} R., 1984, Journal of Astrophysics and
  Astronomy, 5, 369

\bibitem[{{Buchner}(2013)}]{2013ATel.5406....1B}
{Buchner} S., 2013, The Astronomer's Telegram, 5406, 1

\bibitem[{{Buchner}(2010)}]{2010ATel.2768....1B}
{Buchner} S.~J., 2010, The Astronomer's Telegram, 2768

\bibitem[{{Coles} {et~al}\mbox{.}(2011){Coles}, {Hobbs}, {Champion},
  {Manchester}, \& {Verbiest}}]{2011MNRAS.418..561C}
{Coles} W., {Hobbs} G., {Champion} D.~J., {Manchester} R.~N., {Verbiest}
  J.~P.~W., 2011, \mnras, 418, 561

\bibitem[{{Cordes} \& {Downs}(1985)}]{1985ApJS...59..343C}
{Cordes} J.~M., {Downs} G.~S., 1985, \apjs, 59, 343

\bibitem[{{Cordes}, {Downs} \& {Krause-Polstorff}(1988){Cordes}, {Downs}, \&
  {Krause-Polstorff}}]{1988ApJ...330..847C}
{Cordes} J.~M., {Downs} G.~S., {Krause-Polstorff} J., 1988, \apj, 330, 847

\bibitem[{{Cordes} \& {Helfand}(1980)}]{1980ApJ...239..640C}
{Cordes} J.~M., {Helfand} D.~J., 1980, \apj, 239, 640

\bibitem[{{Dodson} {et~al}\mbox{.}(2003){Dodson}, {Legge}, {Reynolds}, \&
  {McCulloch}}]{2003ApJ...596.1137D}
{Dodson} R., {Legge} D., {Reynolds} J.~E., {McCulloch} P.~M., 2003, \apj, 596,
  1137

\bibitem[{{Dodson}, {McCulloch} \& {Lewis}(2002){Dodson}, {McCulloch}, \&
  {Lewis}}]{2002ApJ...564L..85D}
{Dodson} R.~G., {McCulloch} P.~M., {Lewis} D.~R., 2002, \apjl, 564, L85

\bibitem[{{Edwards}, {Hobbs} \& {Manchester}(2006){Edwards}, {Hobbs}, \&
  {Manchester}}]{2006MNRAS.372.1549E}
{Edwards} R.~T., {Hobbs} G.~B., {Manchester} R.~N., 2006, \mnras, 372, 1549

\bibitem[{{Espinoza} {et~al}\mbox{.}(2011){Espinoza}, {Lyne}, {Stappers}, \&
  {Kramer}}]{2011MNRAS.414.1679E}
{Espinoza} C.~M., {Lyne} A.~G., {Stappers} B.~W., {Kramer} M., 2011, \mnras,
  414, 1679

\bibitem[{{Feroz}, {Hobson} \& {Bridges}(2009){Feroz}, {Hobson}, \&
  {Bridges}}]{2009MNRAS.398.1601F}
{Feroz} F., {Hobson} M.~P., {Bridges} M., 2009, \mnras, 398, 1601

\bibitem[{{Flanagan}(1991)}]{1991IAUC.5311....3F}
{Flanagan} C., 1991, \iaucirc, 5311, 3

\bibitem[{{Flanagan}(1994)}]{1994IAUC.6064....2F}
{Flanagan} C., 1994, \iaucirc, 6064, 2

\bibitem[{{Flanagan} \& {McCulloch}(1994)}]{1994IAUC.6038....2F}
{Flanagan} C., {McCulloch} P.~M., 1994, \iaucirc, 6038, 2

\bibitem[{{Groth}(1975)}]{1975ApJS...29..453G}
{Groth} E.~J., 1975, \apjs, 29, 453

\bibitem[{{Handley}, {Hobson} \& {Lasenby}(2015){Handley}, {Hobson}, \&
  {Lasenby}}]{2015arXiv150201856H}
{Handley} W.~J., {Hobson} M.~P., {Lasenby} A.~N., 2015, arXiv:1502.01856

\bibitem[{{Haskell} \& {Melatos}(2015)}]{2015IJMPD..2430008H}
{Haskell} B., {Melatos} A., 2015, International Journal of Modern Physics D,
  24, 30008

\bibitem[{{Hobbs} {et~al}\mbox{.}(2005){Hobbs}, {Lorimer}, {Lyne}, \&
  {Kramer}}]{2005MNRAS.360..974H}
{Hobbs} G., {Lorimer} D.~R., {Lyne} A.~G., {Kramer} M., 2005, \mnras, 360, 974

\bibitem[{{Hobbs}, {Lyne} \& {Kramer}(2010){Hobbs}, {Lyne}, \&
  {Kramer}}]{2010MNRAS.402.1027H}
{Hobbs} G., {Lyne} A.~G., {Kramer} M., 2010, \mnras, 402, 1027

\bibitem[{{Hobbs}, {Edwards} \& {Manchester}(2006){Hobbs}, {Edwards}, \&
  {Manchester}}]{2006MNRAS.369..655H}
{Hobbs} G.~B., {Edwards} R.~T., {Manchester} R.~N., 2006, \mnras, 369, 655

\bibitem[{{Hotan}, {van Straten} \& {Manchester}(2004){Hotan}, {van Straten},
  \& {Manchester}}]{2004PASA...21..302H}
{Hotan} A.~W., {van Straten} W., {Manchester} R.~N., 2004, \pasa, 21, 302

\bibitem[{{Jenet} {et~al}\mbox{.}(1998){Jenet}, {Anderson}, {Kaspi}, {Prince},
  \& {Unwin}}]{1998ApJ...498..365J}
{Jenet} F.~A., {Anderson} S.~B., {Kaspi} V.~M., {Prince} T.~A., {Unwin} S.~C.,
  1998, \apj, 498, 365

\bibitem[{{Johnston} \& {Galloway}(1999)}]{1999MNRAS.306L..50J}
{Johnston} S., {Galloway} D., 1999, \mnras, 306, L50

\bibitem[{{Keith} {et~al}\mbox{.}(2013){Keith}, {Coles}, {Shannon}, {Hobbs},
  {Manchester}, {Bailes}, {Bhat}, {Burke-Spolaor}, {Champion}, {Chaudhary},
  {Hotan}, {Khoo}, {Kocz}, {Os{\l}owski}, {Ravi}, {Reynolds}, {Sarkissian},
  {van Straten}, \& {Yardley}}]{2013MNRAS.429.2161K}
{Keith} M.~J. {et~al.}, 2013, \mnras, 429, 2161

\bibitem[{{Kerr} {et~al}\mbox{.}(2015){Kerr}, {Johnston}, {Hobbs}, \&
  {Shannon}}]{2015arXiv150706982K}
{Kerr} M., {Johnston} S., {Hobbs} G., {Shannon} R.~M., 2015, arXiv:1507.06982

\bibitem[{{Large}, {Vaughan} \& {Mills}(1968){Large}, {Vaughan}, \&
  {Mills}}]{1968Natur.220..340L}
{Large} M.~I., {Vaughan} A.~E., {Mills} B.~Y., 1968, \nat, 220, 340

\bibitem[{{Lentati} {et~al}\mbox{.}(2014){Lentati}, {Alexander}, {Hobson},
  {Feroz}, {van Haasteren}, {Lee}, \& {Shannon}}]{2014MNRAS.437.3004L}
{Lentati} L., {Alexander} P., {Hobson} M.~P., {Feroz} F., {van Haasteren} R.,
  {Lee} K.~J., {Shannon} R.~M., 2014, \mnras, 437, 3004

\bibitem[{{Lentati} \& {Shannon}(2015)}]{2015MNRAS.454.1058L}
{Lentati} L., {Shannon} R.~M., 2015, \mnras, 454, 1058

\bibitem[{{Lentati} {et~al}\mbox{.}(2015){Lentati}, {Taylor}, {Mingarelli},
  {Sesana}, {Sanidas}, {Vecchio}, {Caballero}, {Lee}, {van Haasteren}, {Babak},
  {Bassa}, {Brem}, {Burgay}, {Champion}, {Cognard}, {Desvignes}, {Gair},
  {Guillemot}, {Hessels}, {Janssen}, {Karuppusamy}, {Kramer}, {Lassus},
  {Lazarus}, {Liu}, {Os{\l}owski}, {Perrodin}, {Petiteau}, {Possenti},
  {Purver}, {Rosado}, {Smits}, {Stappers}, {Theureau}, {Tiburzi}, \&
  {Verbiest}}]{2015arXiv150403692L}
{Lentati} L. {et~al.}, 2015, arXiv:1504.03692

\bibitem[{{Livingstone} {et~al}\mbox{.}(2011){Livingstone}, {Ng}, {Kaspi},
  {Gavriil}, \& {Gotthelf}}]{2011ApJ...730...66L}
{Livingstone} M.~A., {Ng} C.-Y., {Kaspi} V.~M., {Gavriil} F.~P., {Gotthelf}
  E.~V., 2011, \apj, 730, 66

\bibitem[{{Lyne} {et~al}\mbox{.}(2010){Lyne}, {Hobbs}, {Kramer}, {Stairs}, \&
  {Stappers}}]{2010Sci...329..408L}
{Lyne} A., {Hobbs} G., {Kramer} M., {Stairs} I., {Stappers} B., 2010, Science,
  329, 408

\bibitem[{{Lyne} {et~al}\mbox{.}(2015){Lyne}, {Jordan}, {Graham-Smith},
  {Espinoza}, {Stappers}, \& {Weltevrede}}]{2015MNRAS.446..857L}
{Lyne} A.~G., {Jordan} C.~A., {Graham-Smith} F., {Espinoza} C.~M., {Stappers}
  B.~W., {Weltevrede} P., 2015, \mnras, 446, 857

\bibitem[{{Lyne} {et~al}\mbox{.}(1996){Lyne}, {Pritchard}, {Graham-Smith}, \&
  {Camilo}}]{1996Natur.381..497L}
{Lyne} A.~G., {Pritchard} R.~S., {Graham-Smith} F., {Camilo} F., 1996, \nat,
  381, 497

\bibitem[{{Manchester} \& {Hobbs}(2011)}]{2011ApJ...736L..31M}
{Manchester} R.~N., {Hobbs} G., 2011, \apjl, 736, L31

\bibitem[{{Melatos} \& {Link}(2014)}]{2014MNRAS.437...21M}
{Melatos} A., {Link} B., 2014, \mnras, 437, 21

\bibitem[{{Melatos} \& {Warszawski}(2009)}]{2009ApJ...700.1524M}
{Melatos} A., {Warszawski} L., 2009, \apj, 700, 1524

\bibitem[{{Petroff} {et~al}\mbox{.}(2013){Petroff}, {Keith}, {Johnston}, {van
  Straten}, \& {Shannon}}]{2013MNRAS.435.1610P}
{Petroff} E., {Keith} M.~J., {Johnston} S., {van Straten} W., {Shannon} R.~M.,
  2013, \mnras, 435, 1610

\bibitem[{{Radhakrishnan} \& {Manchester}(1969)}]{1969Natur.222..228R}
{Radhakrishnan} V., {Manchester} R.~N., 1969, \nat, 222, 228

\bibitem[{{Reichley} \& {Downs}(1969)}]{1969Natur.222..229R}
{Reichley} P.~E., {Downs} G.~S., 1969, \nat, 222, 229

\bibitem[{{Ruderman}(1969)}]{1969Natur.223..597R}
{Ruderman} M., 1969, \nat, 223, 597

\bibitem[{{Shannon} \& {Cordes}(2010)}]{sc2010}
{Shannon} R.~M., {Cordes} J.~M., 2010, \apj, 725, 1607

\bibitem[{{Shannon} {et~al}\mbox{.}(2013){Shannon}, {Cordes}, {Metcalfe},
  {Lazio}, {Cognard}, {Desvignes}, {Janssen}, {Jessner}, {Kramer}, {Lazaridis},
  {Purver}, {Stappers}, \& {Theureau}}]{2013ApJ...766....5S}
{Shannon} R.~M. {et~al.}, 2013, \apj, 766, 5

\bibitem[{{Shannon}, {Johnston} \& {Manchester}(2014){Shannon}, {Johnston}, \&
  {Manchester}}]{2014MNRAS.437.3255S}
{Shannon} R.~M., {Johnston} S., {Manchester} R.~N., 2014, \mnras, 437, 3255

\bibitem[{{Shannon} {et~al}\mbox{.}(2015){Shannon}, {Ravi}, {Lentati}, {Lasky},
  {Hobbs}, {Kerr}, {Manchester}, {Coles}, {Levin}, {Bailes}, {Bhat},
  {Burke-Spolaor}, {Dai}, {Keith}, {Os{\l}owski}, {Reardon}, {van Straten},
  {Toomey}, {Wang}, {Wen}, {Wyithe}, \& {Zhu}}]{2015Sci...349..1522S}
{Shannon} R.~M. {et~al.}, 2015, Science, 349, 1522

\bibitem[{{Taylor}(1992)}]{1992RSPTA.341..117T}
{Taylor} J.~H., 1992, Royal Society of London Philosophical Transactions Series
  A, 341, 117

\bibitem[{{van Eysden} \& {Melatos}(2010)}]{2010MNRAS.409.1253V}
{van Eysden} C.~A., {Melatos} A., 2010, \mnras, 409, 1253

\bibitem[{{van Haasteren} \& {Levin}(2013)}]{2013MNRAS.428.1147V}
{van Haasteren} R., {Levin} Y., 2013, \mnras, 428, 1147

\bibitem[{{van Haasteren} {et~al}\mbox{.}(2009){van Haasteren}, {Levin},
  {McDonald}, \& {Lu}}]{2009MNRAS.395.1005V}
{van Haasteren} R., {Levin} Y., {McDonald} P., {Lu} T., 2009, \mnras, 395, 1005

\bibitem[{{van Haasteren} \& {Vallisneri}(2015)}]{2015MNRAS.446.1170V}
{van Haasteren} R., {Vallisneri} M., 2015, \mnras, 446, 1170

\bibitem[{{Wang} {et~al}\mbox{.}(2000){Wang}, {Manchester}, {Pace}, {Bailes},
  {Kaspi}, {Stappers}, \& {Lyne}}]{2000MNRAS.317..843W}
{Wang} N., {Manchester} R.~N., {Pace} R.~T., {Bailes} M., {Kaspi} V.~M.,
  {Stappers} B.~W., {Lyne} A.~G., 2000, \mnras, 317, 843

\bibitem[{{Wang} {et~al}\mbox{.}(2001){Wang}, {Manchester}, {Zhang}, {Wu},
  {Yusup}, {Lyne}, {Cheng}, \& {Chen}}]{2001MNRAS.328..855W}
{Wang} N., {Manchester} R.~N., {Zhang} J., {Wu} X.~J., {Yusup} A., {Lyne}
  A.~G., {Cheng} K.~S., {Chen} M.~Z., 2001, \mnras, 328, 855

\bibitem[{{Weltevrede} {et~al}\mbox{.}(2010){Weltevrede}, {Johnston},
  {Manchester}, {Bhat}, {Burgay}, {Champion}, {Hobbs}, {K{\i}z{\i}ltan},
  {Keith}, {Possenti}, {Reynolds}, \& {Watters}}]{2010PASA...27...64W}
{Weltevrede} P. {et~al.}, 2010, PASA, 27, 64

\bibitem[{{Yu} {et~al}\mbox{.}(2013){Yu}, {Manchester}, {Hobbs}, {Johnston},
  {Kaspi}, {Keith}, {Lyne}, {Qiao}, {Ravi}, {Sarkissian}, {Shannon}, \&
  {Xu}}]{2013MNRAS.429..688Y}
{Yu} M. {et~al.}, 2013, \mnras, 429, 688

\end{thebibliography}

\appendix


\end{document}


\begin{figure}
\begin{center} 
\begin{tabular}{c}
%\includegraphics[angle=-90,scale=0.5]{1937_mem.eps}
\includegraphics[scale=0.5]{compare_spectra.eps}  \\
\end{tabular}
\caption{  \label{fig:spectrum}  \footnotesize Power spectrum  for observations of PSR~B1259$-$63 (solid thick line).  The dotted line shows the power spectrum derived from a model that only includes Keplerian parameters.  The thin solid lines show models that include additional phase jumps across periastron passages.   These models show evidence for excess power, particularly at the fundamental and first harmonic of the orbital period, which are highlighted by vertical solid lines. 
The dashed line shows the $f^{-8}$ model for the noise, used both as a prewhitening filter to estimate pulsar parameters and the basis for simulations to test the robustness of parameter estimation.    }  \end{center} 
\end{figure}


\begin{table*}
\begin{centering}
\caption{Best-fitting parameters for PSR~B1259$-$63. The spin frequency and frequency derivative parameters presented here were fit while the other parameters were held fixed without use of the Cholesky prewhitening algorithm.  $\nu$ and $\dot{\nu}$ were held fixed when the other parameters were fit. }
\label{tab:params}
\begin{tabular}{ll}
\hline
\hline
\multicolumn{2}{c}{Measured Quantitites} \\ 
\hline
Right ascension, $\alpha$\dotfill &  13:02:47.6426(11) \\ 
Declination, $\delta$\dotfill & $-$63:50:08.665(8) \\ 
Pulse frequency, $\nu$ (s$^{-1}$)\dotfill & 20.93692420(15) \\ 
First derivative of pulse frequency, $\dot{\nu}$ (s$^{-2}$)\dotfill & $-$9.989(5)$\times 10^{-13}$ \\ 
Proper motion in right ascension, $\mu_{\alpha} \cos \delta$ (mas\,yr$^{-1}$)\dotfill & $-$6.6(1.8) \\ 
Proper motion in declination, $\mu_{\delta}$ (mas\,yr$^{-1}$)\dotfill & $-$4.4(1.4) \\ 
Orbital period, $P_b$ (d)\dotfill & 1236.724526(6) \\ 
First derivative of orbital period, $\dot{P_b}$\dotfill & 1.4(7)$\times 10^{-8}$ \\ 
Epoch of periastron, $T_0$ (MJD)\dotfill & 53071.2447290(7) \\ 
Projected semi-major axis of orbit, $x$ (lt-s)\dotfill & 1296.27448(14) \\ 
First derivative of $x$, $\dot{x}$  (lt-s s$^{-1}$)\dotfill & 2.100(8)$\times 10^{-11}$ \\ 
Second derivative of $x$, $\ddot{x}$ (lt-s s$^{-2}$) \dotfill & 2.0(7)$\times 10^{-20}$ \\ 
Longitude of periastron, $\omega_0$ (deg)\dotfill & 138.665013(11) \\ 
Orbital eccentricity, $e$\dotfill & 0.86987970(6) \\ 
Periastron advance, $\dot{\omega}$ (deg~yr$^{-1}$)\dotfill & 7.81(3)$\times 10^{-5}$ \\ 
Periastron acceleration $\ddot{\omega}$ (deg~yr$^{-2}$)\dotfill & 1.8(3)$\times 10^{-6}$ \\ 
Glitch epoch \dotfill & 50708.0(5) \\ 
Glitch frequency permanent frequency change \dotfill & 3.1(3)$\times 10^{-8}$ \\ 
Glitch transient frequency change \dotfill & 1.8(3)$\times 10^{-8}$ \\ 
Glitch transient decay time\dotfill & 82(18) \\ 
\hline
\multicolumn{2}{c}{Fixed Quantities} \\ 
\hline
Epoch of frequency determination (MJD)\dotfill & 50357 \\ 
Epoch of position determination (MJD)\dotfill & 55000 \\ 
Epoch of dispersion measure determination (MJD)\dotfill & 50357 \\ 
Dispersion measure, DM (pc~cm$^{-3}$)\dotfill & 146.8 \\ 
\hline
\multicolumn{2}{c}{Derived Quantities} \\
\hline
$\log_{10}$(Characteristic age, yr) \dotfill & 5.52 \\
$\log_{10}$(Surface magnetic field strength, G) \dotfill & 11.52 \\
\hline
\multicolumn{2}{c}{Assumptions} \\
\hline
Clock correction procedure\dotfill & TT(BIPM2011) \\
Solar system ephemeris model\dotfill & DE421 \\
\hline
\end{tabular}
\end{centering}
\end{table*}


\begin{table*}
\begin{centering}
\caption{Best-fitting parameters for PSR~B1259$-$63. The spin frequency and frequency derivative parameters presented here were fit while the other parameters were held fixed without use of the Cholesky prewhitening algorithm.  $\nu$ and $\dot{\nu}$ were held fixed when the other parameters were fit. }
\label{tab:params}
\begin{tabular}{ll}
\hline
\hline
\multicolumn{2}{c}{Measured Quantitites} \\ 
\hline
Right ascension, $\alpha$\dotfill &  13:02:47.6426(11) \\ 
Declination, $\delta$\dotfill & $-$63:50:08.665(8) \\ 
Pulse frequency, $\nu$ (s$^{-1}$)\dotfill & 20.93692420(15) \\ 
First derivative of pulse frequency, $\dot{\nu}$ (s$^{-2}$)\dotfill & $-$9.989(5)$\times 10^{-13}$ \\ 
Proper motion in right ascension, $\mu_{\alpha} \cos \delta$ (mas\,yr$^{-1}$)\dotfill & $-$6.6(1.8) \\ 
Proper motion in declination, $\mu_{\delta}$ (mas\,yr$^{-1}$)\dotfill & $-$4.4(1.4) \\ 
Orbital period, $P_b$ (d)\dotfill & 1236.724526(6) \\ 
First derivative of orbital period, $\dot{P_b}$\dotfill & 1.4(7)$\times 10^{-8}$ \\ 
Epoch of periastron, $T_0$ (MJD)\dotfill & 53071.2447290(7) \\ 
Projected semi-major axis of orbit, $x$ (lt-s)\dotfill & 1296.27448(14) \\ 
First derivative of $x$, $\dot{x}$  (lt-s s$^{-1}$)\dotfill & 2.100(8)$\times 10^{-11}$ \\ 
Second derivative of $x$, $\ddot{x}$ (lt-s s$^{-2}$) \dotfill & 2.0(7)$\times 10^{-20}$ \\ 
Longitude of periastron, $\omega_0$ (deg)\dotfill & 138.665013(11) \\ 
Orbital eccentricity, $e$\dotfill & 0.86987970(6) \\ 
Periastron advance, $\dot{\omega}$ (deg~yr$^{-1}$)\dotfill & 7.81(3)$\times 10^{-5}$ \\ 
Periastron acceleration $\ddot{\omega}$ (deg~yr$^{-2}$)\dotfill & 1.8(3)$\times 10^{-6}$ \\ 
Glitch epoch \dotfill & 50708.0(5) \\ 
Glitch frequency permanent frequency change \dotfill & 3.1(3)$\times 10^{-8}$ \\ 
Glitch transient frequency change \dotfill & 1.8(3)$\times 10^{-8}$ \\ 
Glitch transient decay time\dotfill & 82(18) \\ 
\hline
\multicolumn{2}{c}{Fixed Quantities} \\ 
\hline
Epoch of frequency determination (MJD)\dotfill & 50357 \\ 
Epoch of position determination (MJD)\dotfill & 55000 \\ 
Epoch of dispersion measure determination (MJD)\dotfill & 50357 \\ 
Dispersion measure, DM (pc~cm$^{-3}$)\dotfill & 146.8 \\ 
\hline
\multicolumn{2}{c}{Derived Quantities} \\
\hline
$\log_{10}$(Characteristic age, yr) \dotfill & 5.52 \\
$\log_{10}$(Surface magnetic field strength, G) \dotfill & 11.52 \\
\hline
\multicolumn{2}{c}{Assumptions} \\
\hline
Clock correction procedure\dotfill & TT(BIPM2011) \\
Solar system ephemeris model\dotfill & DE421 \\
\hline
\end{tabular}
\end{centering}
\end{table*}



\begin{table*}
\begin{centering}
\caption{Glitch parameters}
\label{tab:glitch_params}
\begin{tabular}{lrrrrrrc}
\hline
\hline
 \multicolumn{1}{c}{MJD}  & \multicolumn{1}{c}{$\Delta \nu$} &  \multicolumn{1}{c}{$\Delta \nu_t$} & \multicolumn{1}{c}{$\Delta \nu/ \nu$} &  \multicolumn{1}{c}{$\Delta \nu_t/ \nu$} &  \multicolumn{1}{c}{Q} & \multicolumn{1}{c}{$(\Delta \nu/\nu)_{\rm lit}$} &  \multicolumn{1}{c}{Ref.}\\
 & ($\mu$Hz)& ($\mu$Hz) &  \multicolumn{1}{c}{($10^{-6}$)} &\multicolumn{1}{c}{($10^{-6}$)}& & \multicolumn{1}{c}{($10^{-6}$)} &\\
\hline
49591 &	1.4(3)	&0.7(2)	&\nodata	&0.13(2)	&0.06(2)	&0.199(2) &	3\\
49559&	1.8(2)	&7.8(2)	&\nodata	&0.16(2)	&0.7(2)	&0.835(2)	&2\\
50369	&12.3(3)	&11.5(3)	&0.11(30	&1.10(2)	&1.03(20	&2.11(2)	&4\\
51559 &	22.7(6)	&12.1(6)	&0.18(3)	&2.03(5)	&1.08(5)	&3.152(2)	&4\\
56555 &	20.9(2)	&13.2(2)	&0.21(2)	&1.87(2)	&1.18(2)	&3.100()	&5\\
55428&	9.3(3)	&11.7(3)	&\nodata	&0.83(3)	&1.05(3)	&\nodata	&\nodata\\
53960&	19.5(4)	&9.2(3)	&0.3(1)	&1.74(3)	&0.83(3)	&2.62	&4\\
53193&	10.6(3)	&12.3(3)	&0.1(1)	&0.95(3)	&1.10(3)	&2.100()	&4\\
48957&	\nodata	&16(4)	&\nodata	&1.5(4)	&\nodata	&0.835	&2\\
\hline
\end{tabular}
\end{centering}
\begin{flushleft}
For each glitch we list the MJD, the permanent change in spin frequency $\Delta \nu$, transient changes in spin frequency  $\Delta \nu_g$, reported in absolute value and as a fraction of the pulsar spin frequency $\nu$.   We also report the fraction of the permanent change associated with the transient component  $Q = \Delta \nu_t/\Delta \nu$, and, where available, previously reported glitch-decay amplitudes, $\Delta \nu_{\rm lit}$. 
The values in parentheses represent the nominal $1-\sigma$ uncertainties for our measurements and, where available, previous measurements. 
References are (1) \cite{1991IAUC.5311....3F};  (2) \cite{1994IAUC.6038....2F}; (3)  \cite{1994IAUC.6064....2F} ;  (4) \cite{2013MNRAS.429..688Y} ; and  (5) \cite{2013ATel.5406....1B}. 
\end{flushleft}
\end{table*}



1604.02661.extracted_bib
[{{Alpar {et~al\mbox{.(1993){Alpar, {Chau, {Cheng, \& {Pines]{1993ApJ...409..345A {Alpar M.~A., {Chau H.~F., {Cheng K.~S., {Pines D., 1993, \apj, 409, 345
1604.02661.extracted_bib
[{{Alpar, {Cheng \& {Pines(1989){Alpar, {Cheng, \& {Pines]{1989ApJ...346..823A {Alpar M.~A., {Cheng K.~S., {Pines D., 1989, \apj, 346, 823
1604.02661.extracted_bib
[{{Alpar {et~al\mbox{.(1984){Alpar, {Pines, {Anderson, \& {Shaham]{1984ApJ...276..325A {Alpar M.~A., {Pines D., {Anderson P.~W., {Shaham J., 1984, \apj, 276, 325
1604.02661.extracted_bib
[{{Anderson \& {Itoh(1975)]{1975Natur.256...25A {Anderson P.~W., {Itoh N., 1975, \nat, 256, 25
1604.02661.extracted_bib
[{{Baym, {Pethick \& {Pines(1969){Baym, {Pethick, \& {Pines]{1969Natur.224..673B {Baym G., {Pethick C., {Pines D., 1969, \nat, 224, 673
1604.02661.extracted_bib
[{{Blandford, {Romani \& {Narayan(1984){Blandford, {Romani, \& {Narayan]{1984JApA....5..369B {Blandford R., {Romani R.~W., {Narayan R., 1984, Journal of Astrophysics and Astronomy, 5, 369
1604.02661.extracted_bib
[{{Buchner(2013)]{2013ATel.5406....1B {Buchner S., 2013, The Astronomer's Telegram, 5406, 1
1604.02661.extracted_bib
[{{Buchner(2010)]{2010ATel.2768....1B {Buchner S.~J., 2010, The Astronomer's Telegram, 2768
1604.02661.extracted_bib
[{{Coles {et~al\mbox{.(2011){Coles, {Hobbs, {Champion, {Manchester, \& {Verbiest]{2011MNRAS.418..561C {Coles W., {Hobbs G., {Champion D.~J., {Manchester R.~N., {Verbiest J.~P.~W., 2011, \mnras, 418, 561
1604.02661.extracted_bib
[{{Cordes \& {Downs(1985)]{1985ApJS...59..343C {Cordes J.~M., {Downs G.~S., 1985, \apjs, 59, 343
1604.02661.extracted_bib
[{{Cordes, {Downs \& {Krause-Polstorff(1988){Cordes, {Downs, \& {Krause-Polstorff]{1988ApJ...330..847C {Cordes J.~M., {Downs G.~S., {Krause-Polstorff J., 1988, \apj, 330, 847
1604.02661.extracted_bib
[{{Cordes \& {Helfand(1980)]{1980ApJ...239..640C {Cordes J.~M., {Helfand D.~J., 1980, \apj, 239, 640
1604.02661.extracted_bib
[{{Dodson {et~al\mbox{.(2003){Dodson, {Legge, {Reynolds, \& {McCulloch]{2003ApJ...596.1137D {Dodson R., {Legge D., {Reynolds J.~E., {McCulloch P.~M., 2003, \apj, 596, 1137
1604.02661.extracted_bib
[{{Dodson, {McCulloch \& {Lewis(2002){Dodson, {McCulloch, \& {Lewis]{2002ApJ...564L..85D {Dodson R.~G., {McCulloch P.~M., {Lewis D.~R., 2002, \apjl, 564, L85
1604.02661.extracted_bib
[{{Edwards, {Hobbs \& {Manchester(2006){Edwards, {Hobbs, \& {Manchester]{2006MNRAS.372.1549E {Edwards R.~T., {Hobbs G.~B., {Manchester R.~N., 2006, \mnras, 372, 1549
1604.02661.extracted_bib
[{{Espinoza {et~al\mbox{.(2011){Espinoza, {Lyne, {Stappers, \& {Kramer]{2011MNRAS.414.1679E {Espinoza C.~M., {Lyne A.~G., {Stappers B.~W., {Kramer M., 2011, \mnras, 414, 1679
1604.02661.extracted_bib
[{{Feroz, {Hobson \& {Bridges(2009){Feroz, {Hobson, \& {Bridges]{2009MNRAS.398.1601F {Feroz F., {Hobson M.~P., {Bridges M., 2009, \mnras, 398, 1601
1604.02661.extracted_bib
[{{Flanagan(1991)]{1991IAUC.5311....3F {Flanagan C., 1991, \iaucirc, 5311, 3
1604.02661.extracted_bib
[{{Flanagan(1994)]{1994IAUC.6064....2F {Flanagan C., 1994, \iaucirc, 6064, 2
1604.02661.extracted_bib
[{{Flanagan \& {McCulloch(1994)]{1994IAUC.6038....2F {Flanagan C., {McCulloch P.~M., 1994, \iaucirc, 6038, 2
1604.02661.extracted_bib
[{{Groth(1975)]{1975ApJS...29..453G {Groth E.~J., 1975, \apjs, 29, 453
1604.02661.extracted_bib
[{{Handley, {Hobson \& {Lasenby(2015){Handley, {Hobson, \& {Lasenby]{2015arXiv150201856H {Handley W.~J., {Hobson M.~P., {Lasenby A.~N., 2015, arXiv:1502.01856
1604.02661.extracted_bib
[{{Haskell \& {Melatos(2015)]{2015IJMPD..2430008H {Haskell B., {Melatos A., 2015, International Journal of Modern Physics D, 24, 30008
1604.02661.extracted_bib
[{{Hobbs {et~al\mbox{.(2005){Hobbs, {Lorimer, {Lyne, \& {Kramer]{2005MNRAS.360..974H {Hobbs G., {Lorimer D.~R., {Lyne A.~G., {Kramer M., 2005, \mnras, 360, 974
1604.02661.extracted_bib
[{{Hobbs, {Lyne \& {Kramer(2010){Hobbs, {Lyne, \& {Kramer]{2010MNRAS.402.1027H {Hobbs G., {Lyne A.~G., {Kramer M., 2010, \mnras, 402, 1027
1604.02661.extracted_bib
[{{Hobbs, {Edwards \& {Manchester(2006){Hobbs, {Edwards, \& {Manchester]{2006MNRAS.369..655H {Hobbs G.~B., {Edwards R.~T., {Manchester R.~N., 2006, \mnras, 369, 655
1604.02661.extracted_bib
[{{Hotan, {van Straten \& {Manchester(2004){Hotan, {van Straten, \& {Manchester]{2004PASA...21..302H {Hotan A.~W., {van Straten W., {Manchester R.~N., 2004, \pasa, 21, 302
1604.02661.extracted_bib
[{{Jenet {et~al\mbox{.(1998){Jenet, {Anderson, {Kaspi, {Prince, \& {Unwin]{1998ApJ...498..365J {Jenet F.~A., {Anderson S.~B., {Kaspi V.~M., {Prince T.~A., {Unwin S.~C., 1998, \apj, 498, 365
1604.02661.extracted_bib
[{{Johnston \& {Galloway(1999)]{1999MNRAS.306L..50J {Johnston S., {Galloway D., 1999, \mnras, 306, L50
1604.02661.extracted_bib
[{{Keith {et~al\mbox{.(2013){Keith, {Coles, {Shannon, {Hobbs, {Manchester, {Bailes, {Bhat, {Burke-Spolaor, {Champion, {Chaudhary, {Hotan, {Khoo, {Kocz, {Os{\lowski, {Ravi, {Reynolds, {Sarkissian, {van Straten, \& {Yardley]{2013MNRAS.429.2161K {Keith M.~J. {et~al., 2013, \mnras, 429, 2161
1604.02661.extracted_bib
[{{Kerr {et~al\mbox{.(2015){Kerr, {Johnston, {Hobbs, \& {Shannon]{2015arXiv150706982K {Kerr M., {Johnston S., {Hobbs G., {Shannon R.~M., 2015, arXiv:1507.06982
1604.02661.extracted_bib
[{{Large, {Vaughan \& {Mills(1968){Large, {Vaughan, \& {Mills]{1968Natur.220..340L {Large M.~I., {Vaughan A.~E., {Mills B.~Y., 1968, \nat, 220, 340
1604.02661.extracted_bib
[{{Lentati {et~al\mbox{.(2014){Lentati, {Alexander, {Hobson, {Feroz, {van Haasteren, {Lee, \& {Shannon]{2014MNRAS.437.3004L {Lentati L., {Alexander P., {Hobson M.~P., {Feroz F., {van Haasteren R., {Lee K.~J., {Shannon R.~M., 2014, \mnras, 437, 3004
1604.02661.extracted_bib
[{{Lentati \& {Shannon(2015)]{2015MNRAS.454.1058L {Lentati L., {Shannon R.~M., 2015, \mnras, 454, 1058
1604.02661.extracted_bib
[{{Lentati {et~al\mbox{.(2015){Lentati, {Taylor, {Mingarelli, {Sesana, {Sanidas, {Vecchio, {Caballero, {Lee, {van Haasteren, {Babak, {Bassa, {Brem, {Burgay, {Champion, {Cognard, {Desvignes, {Gair, {Guillemot, {Hessels, {Janssen, {Karuppusamy, {Kramer, {Lassus, {Lazarus, {Liu, {Os{\lowski, {Perrodin, {Petiteau, {Possenti, {Purver, {Rosado, {Smits, {Stappers, {Theureau, {Tiburzi, \& {Verbiest]{2015arXiv150403692L {Lentati L. {et~al., 2015, arXiv:1504.03692
1604.02661.extracted_bib
[{{Livingstone {et~al\mbox{.(2011){Livingstone, {Ng, {Kaspi, {Gavriil, \& {Gotthelf]{2011ApJ...730...66L {Livingstone M.~A., {Ng C.-Y., {Kaspi V.~M., {Gavriil F.~P., {Gotthelf E.~V., 2011, \apj, 730, 66
1604.02661.extracted_bib
[{{Lyne {et~al\mbox{.(2010){Lyne, {Hobbs, {Kramer, {Stairs, \& {Stappers]{2010Sci...329..408L {Lyne A., {Hobbs G., {Kramer M., {Stairs I., {Stappers B., 2010, Science, 329, 408
1604.02661.extracted_bib
[{{Lyne {et~al\mbox{.(2015){Lyne, {Jordan, {Graham-Smith, {Espinoza, {Stappers, \& {Weltevrede]{2015MNRAS.446..857L {Lyne A.~G., {Jordan C.~A., {Graham-Smith F., {Espinoza C.~M., {Stappers B.~W., {Weltevrede P., 2015, \mnras, 446, 857
1604.02661.extracted_bib
[{{Lyne {et~al\mbox{.(1996){Lyne, {Pritchard, {Graham-Smith, \& {Camilo]{1996Natur.381..497L {Lyne A.~G., {Pritchard R.~S., {Graham-Smith F., {Camilo F., 1996, \nat, 381, 497
1604.02661.extracted_bib
[{{Manchester \& {Hobbs(2011)]{2011ApJ...736L..31M {Manchester R.~N., {Hobbs G., 2011, \apjl, 736, L31
1604.02661.extracted_bib
[{{Melatos \& {Link(2014)]{2014MNRAS.437...21M {Melatos A., {Link B., 2014, \mnras, 437, 21
1604.02661.extracted_bib
[{{Melatos \& {Warszawski(2009)]{2009ApJ...700.1524M {Melatos A., {Warszawski L., 2009, \apj, 700, 1524
1604.02661.extracted_bib
[{{Petroff {et~al\mbox{.(2013){Petroff, {Keith, {Johnston, {van Straten, \& {Shannon]{2013MNRAS.435.1610P {Petroff E., {Keith M.~J., {Johnston S., {van Straten W., {Shannon R.~M., 2013, \mnras, 435, 1610
1604.02661.extracted_bib
[{{Radhakrishnan \& {Manchester(1969)]{1969Natur.222..228R {Radhakrishnan V., {Manchester R.~N., 1969, \nat, 222, 228
1604.02661.extracted_bib
[{{Reichley \& {Downs(1969)]{1969Natur.222..229R {Reichley P.~E., {Downs G.~S., 1969, \nat, 222, 229
1604.02661.extracted_bib
[{{Ruderman(1969)]{1969Natur.223..597R {Ruderman M., 1969, \nat, 223, 597
1604.02661.extracted_bib
[{{Shannon \& {Cordes(2010)]{sc2010 {Shannon R.~M., {Cordes J.~M., 2010, \apj, 725, 1607
1604.02661.extracted_bib
[{{Shannon {et~al\mbox{.(2013){Shannon, {Cordes, {Metcalfe, {Lazio, {Cognard, {Desvignes, {Janssen, {Jessner, {Kramer, {Lazaridis, {Purver, {Stappers, \& {Theureau]{2013ApJ...766....5S {Shannon R.~M. {et~al., 2013, \apj, 766, 5
1604.02661.extracted_bib
[{{Shannon, {Johnston \& {Manchester(2014){Shannon, {Johnston, \& {Manchester]{2014MNRAS.437.3255S {Shannon R.~M., {Johnston S., {Manchester R.~N., 2014, \mnras, 437, 3255
1604.02661.extracted_bib
[{{Shannon {et~al\mbox{.(2015){Shannon, {Ravi, {Lentati, {Lasky, {Hobbs, {Kerr, {Manchester, {Coles, {Levin, {Bailes, {Bhat, {Burke-Spolaor, {Dai, {Keith, {Os{\lowski, {Reardon, {van Straten, {Toomey, {Wang, {Wen, {Wyithe, \& {Zhu]{2015Sci...349..1522S {Shannon R.~M. {et~al., 2015, Science, 349, 1522
1604.02661.extracted_bib
[{{Taylor(1992)]{1992RSPTA.341..117T {Taylor J.~H., 1992, Royal Society of London Philosophical Transactions Series A, 341, 117
1604.02661.extracted_bib
[{{van Eysden \& {Melatos(2010)]{2010MNRAS.409.1253V {van Eysden C.~A., {Melatos A., 2010, \mnras, 409, 1253
1604.02661.extracted_bib
[{{van Haasteren \& {Levin(2013)]{2013MNRAS.428.1147V {van Haasteren R., {Levin Y., 2013, \mnras, 428, 1147
1604.02661.extracted_bib
[{{van Haasteren {et~al\mbox{.(2009){van Haasteren, {Levin, {McDonald, \& {Lu]{2009MNRAS.395.1005V {van Haasteren R., {Levin Y., {McDonald P., {Lu T., 2009, \mnras, 395, 1005
1604.02661.extracted_bib
[{{van Haasteren \& {Vallisneri(2015)]{2015MNRAS.446.1170V {van Haasteren R., {Vallisneri M., 2015, \mnras, 446, 1170
1604.02661.extracted_bib
[{{Wang {et~al\mbox{.(2000){Wang, {Manchester, {Pace, {Bailes, {Kaspi, {Stappers, \& {Lyne]{2000MNRAS.317..843W {Wang N., {Manchester R.~N., {Pace R.~T., {Bailes M., {Kaspi V.~M., {Stappers B.~W., {Lyne A.~G., 2000, \mnras, 317, 843
1604.02661.extracted_bib
[{{Wang {et~al\mbox{.(2001){Wang, {Manchester, {Zhang, {Wu, {Yusup, {Lyne, {Cheng, \& {Chen]{2001MNRAS.328..855W {Wang N., {Manchester R.~N., {Zhang J., {Wu X.~J., {Yusup A., {Lyne A.~G., {Cheng K.~S., {Chen M.~Z., 2001, \mnras, 328, 855
1604.02661.extracted_bib
[{{Weltevrede {et~al\mbox{.(2010){Weltevrede, {Johnston, {Manchester, {Bhat, {Burgay, {Champion, {Hobbs, {K{\iz{\iltan, {Keith, {Possenti, {Reynolds, \& {Watters]{2010PASA...27...64W {Weltevrede P. {et~al., 2010, PASA, 27, 64
1604.02661.extracted_bib
[{{Yu {et~al\mbox{.(2013){Yu, {Manchester, {Hobbs, {Johnston, {Kaspi, {Keith, {Lyne, {Qiao, {Ravi, {Sarkissian, {Shannon, \& {Xu]{2013MNRAS.429..688Y {Yu M. {et~al., 2013, \mnras, 429, 688