Profile Stochasticity in PSR J1909-3744
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 25 September 2018 (MN L A TEX style file v2.2)
Profile Stochasticity in PSR J1909 − L. Lentati (cid:63) , R. M. Shannon , Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge, CB3 0HE, UK CSIRO Astronomy and Space Science, Australia Telescope National Facility, Box 76 Epping, NSW, 1710, Australia International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia
25 September 2018
ABSTRACT
We extend the recently introduced Bayesian framework ‘Generative Pulsar Timing Analysis’to incorporate both pulse jitter (high frequency variation in the arrival time of the pulse) andepoch to epoch stochasticity in the shape of the pulse profile. This framework allows for a fulltiming analysis to be performed on the folded profile data, rather than the site arrival times asis typical in most timing studies. We apply this extended framework both to simulations, andto an 11 yr, 10 cm data set for PSR J1909 − − γ = / × − , consistent with previouslypublished limits. Key words: methods: data analysis, pulsars: general, pulsars:individual
The emission of electromagnetic radiation from pulsars in the radioband is thought to be the result of highly energetic charged particlesbeing accelerated, and then escaping, along open magnetic fieldlines in the pulsar’s magnetosphere. When the magnetic axis is notperfectly aligned with the rotational axis, the result is a beam ofradiation that sweeps across the sky. If an observer lies in the pathof this beam, they will observe a series of regularly timed pulses,giving the pulsar its trademark lighthouse e ff ect.The time of arrival (TOA) of these pulses is predicted by thepulsar’s ‘timing model’. This is a deterministic description of thepulsar, describing its rotational period and spin down, astrometricproperties, and if the pulsar is part of a binary system, additionalparameters such as its orbital period, and geometry of the orbit withrespect to the Earth.When a pulsar is observed for the purpose of timing, the in-dividual pulses from a single observing epoch are folded togetherusing some fiducial timing model to form an average pulse pro-file for that observation. This averaged data can then be comparedto a model for the pulse profile, a process that results in the cre-ation of a set of observed TOAs, referred to as site arrival times (cid:63) E-mail: [email protected] (SATs). These SATs can then be compared to those predicted bythe pulsar’s timing model, where the di ff erence between the twoare known as the ‘timing residuals’. These residuals carry physicalinformation about the unmodelled e ff ects in the pulse propagation,including intrinsic variation in the pulse emission time (e.g. Shan-non et al. 2014; Hobbs et al. 2010), or extrinsic factors, such asperturbations due to gravitational waves (GWs) (e.g. Sazhin 1978;Detweiler 1979).The most stringent 95% upper limit on the amplitude of anisotropic stochastic GW background (GWB) formed from the in-coherent superposition of GWs emitted from a large number ofmerging supermassive black-hole binaries is A < × − for areference frequency of one yr − (Shannon et al. 2015, henceforthS15). This is equivalent to a red noise process with an rms ampli-tude of ∼
70 ns in a 10 yr data set. In Fig. 1, we show the e ff ectof a 100ns deviation due to a GWB on the arrival time of the aver-age profile at a single observational epoch. In the PSR J1909 − ∼ one tenth of a phase bin.Several approaches to building a profile template have previ-ously been used. For example, a template can be formed by aver-aging over the individual profiles at each observational epoch tocreate a single high S / N profile that can then be used to form the c (cid:13) a r X i v : . [ a s t r o - ph . I M ] S e p L. Lentati et al. −0.05 0 0.05 0.1 0.15Time (ms) A r b i t r a r y F l u x −4 −2 0 2x 10 −3 Figure 1.
Noiseless model for the deterministic profile of PSR J1909 − ff ect of a passing GW from a 10 − isotropic GWB (red line). Thetimespan along the x-axis covers 7.5% of a full rotation for this pulsar. Inthe PSR J1909 − ∼ one tenth of a phase bin. TOAs directly. This, however, has the disadvantage that any noisepresent in the template can bias the resulting TOA estimates (Hotanet al. 2005), especially if the template noise is correlated with noisein the profiles. Alternatively, rather than use the averaged data, asmoothed version of the template can be used (e.g. Demorest et al.2013), or analytic functions can be fit to the averaged profile toform a noise free template (e.g. Manchester et al. 2013).Once a template has been developed it is then used to form theTOAs for each observational epoch. This is most commonly donevia the ‘Fourier phase-gradient method’ (Taylor 1992) in which thephase o ff set between the two is computed using the Fourier trans-form of both the template, and the profile at each epoch, and a crosscorrelation between the two performed. Alternative time domainapproaches have also been used (e.g. Hotan et al. 2005), howeverregardless of the approach, they all share a common assumption;that the profile is stable within radiometer noise from epoch toepoch.While long term stability of pulse profiles has been shownin some pulsars (e.g. Shao et al. 2013), epoch to epoch variationin the profile shape has also been observed. For example, a studyof morphological variability in PSR J1022 + − ff ects of these distortions have been modelled in the TOA do-main. This is done both by including additional white noise param-eters referred to as ‘EFAC’ and ‘EQUAD’, which scale and add inquadrature to the formal TOA uncertainties, and by incorporating amodel for low frequency timing noise into the analysis (e.g. Coleset al. 2011; van Haasteren & Levin 2013; Lentati et al. 2014). Thishas the disadvantage that in a single pulsar these distortions couldbe covariant with the GW signal (see Figure 2).In the top panel of Fig. 2, we compare simulated residuals in-duced by the GW signal from an isotropic stochastic backgroundwith an amplitude of 1 × − (black line), with those that resultfrom the passage of an additional Gaussian component throughthe profile, with an amplitude of 0.5% that of the observed pro-file, which is not appropriately modelled by the single averageprofile used to form the TOAs (see bottom panels). All TOAswere simulated using the highest signal to noise profile in the PSRJ1909 − − ff ersome concluding remarks in Section 9. The methods used in this analysis are drawn from those presentedin L15. Here, our pulsar timing analysis is performed entirely withthe profile data, rather than the TOAs formed from those profiles.Qualitatively, in each likelihood calculation, we construct a modelfor the deterministic (or average) profile using a shapelet basis, andgenerate a model time of arrival at each observational epoch forthat profile using the pulsars timing model. Both these steps oc-cur simultaneously, such that both the parameters that describe theshapelet model, and the timing model parameters are free to varywithin our analysis. c (cid:13) , 000–000 rofile Stochasticity in PSR J1909 − MJD R e s i dua l ( n s ) P r o f il e A m p li t ude ( % )
424 512 600−101 Profile Bin 424 512 600Profile Bin P r o f il e R e s i dua l ( % )
424 512 600Profile Bin
Figure 2. (Top) Black line - Simulated residuals due to a GW signal from an isotropic stochastic background with an amplitude of 1 × − , consistent withthe most stringent 95% upper limits set by Shannon et al. (2015). TOAs were simulated using the highest signal to noise profile in the PSR J1909 − ff erent positions in the main profile, and residuals from the profile fit. While a full description of the general framework we will useis available in L15, in this work we will be extending the method-ology significantly to incorporate the possibility of epoch to epochvariation in the profile. We include models for pulse jitter – a shiftin the arrival time of the deterministic profile – as well as shapevariation that could be of instrumental, or astrophysical origin. Assuch we will include below an overview of the basic framework,before providing details on the modifications required to incorpo-rate profile stochasticity.
A thorough description of the Shapelet formalism can be found inRefregier (2003), with astronomical uses being described in e.g,Kelly & McKay (2004); Lentati et al. (2013); Refregier & Bacon(2003). Here we give only an outline to aid later discussion.Shapelets are described by a set of dimensionless basis func-tions, which in one dimension can be written as: φ n ( x ) ≡ (cid:104) n n ! √ π (cid:105) − / H n (cid:32) x √ (cid:33) exp (cid:32) − x (cid:33) , (1)where n is a non-negative integer, and H n is the Hermite polyno-mial of order n . Therefore the lowest order shapelet is given by astandard Gaussian ( H ( x ) = Λ which is a free parameter to be fitted for, in order to construct dimensional basisfunctions: B n ( x ; Λ ) ≡ Λ − / φ n ( Λ − x ) . (2)These basis functions are orthonormal, i.e: (cid:90) ∞−∞ d x B n ( x ; Λ ) B m ( x ; Λ ) = δ nm , (3)where δ nm is the Kronecker delta, so that we can represent a func-tion s ( x ) as the sum: s ( x , ζ, Λ ) = n max (cid:88) n = ζ n B n ( x ; Λ ) , (4)where ζ n are the shapelet amplitudes, and n max the number ofshapelet basis vectors included in the model.In our analysis of pulsar timing data, we form a single profileshape, which will then be scaled from epoch to epoch. We thereforeredefine Eq. 4, such that we have a single global amplitude A , and n max − ζ n which are the amplitudes for the shapeletcomponents that have n >
0. These therefore represent the relativecontribution to the overall profile shape, relative to the zeroth-orderterm, which we take to have an amplitude of 1. Written in this wayEq. 4 becomes: s ( x , A , ζ, Λ ) = A n max (cid:88) n = ζ n B n ( x ; Λ ) . (5) c (cid:13) , 000–000 L. Lentati et al.
Finally, the total integrated flux in the model profile, F tot , isgiven by F tot = (cid:90) ∞−∞ d x s ( x ) = A (cid:88) n = even ζ n (cid:104) − n √ π Λ (cid:105) (cid:32) nn / (cid:33) . (6)While there is no natural basis in which to describe the shapeof the pulse profile, shapelets present several advantages whencompared to, for example, a set of von Mises functions. A shapeletmodel requires a single width parameter, and a set of n max am-plitudes which multiply the linear basis vectors given that width.Phase wrapped Gaussians, or von Mises functions, however, have n widths and n − We begin by considering our data, d , a set of N d integrated pulseprofiles, where the profile at observational epoch i consists of a setof N i values representing the flux density of the profile as measuredat a set of times t i . We represent d as a sum of both a deterministicand a stochastic component: d = s + n , (7)where d represents the (cid:80) N d i = N i total data points for the N d pulse pro-files for a single pulsar, with s and n the deterministic and stochas-tic contributions to the total respectively, where any contributionsto the latter will be modelled as random Gaussian processes.We first consider the stochastic component of the signal, n , tobe described solely by an uncorrelated random Gaussian processwith RMS σ i , determined individually at each epoch i . The deter-ministic component, s , consists only of (i) our shapelet model forthe pulse profiles with the centroid position for each model pro-file determined using the set of arrival times τ ( (cid:15) ) predicted by thepulsar’s m timing model parameters, (cid:15) , and (ii) an arbitrary base-line o ff set for each profile epoch. With the inclusion of the timingmodel, we can rewrite Eq. 5 to describe the shapelet model for aparticular epoch i as: s i ( t , A , ζ, Λ , (cid:15) ) = A n max (cid:88) n = ζ n B n ( t − τ i ( (cid:15) ); Λ ) , (8)with τ i the barycentric arrival time for pulse i predicted by the setof timing model parameters (cid:15) . In L15 it was assumed that the correction to the SSB for eachprofile could be computed at the SAT obtained previously usingtraditional techniques for that profile. In our analysis of the PSRJ1909 − ff sets between di ff erent ob-serving systems (referred to as ‘Jumps’) means that we must firstcorrect the SATs using the values for these o ff sets, and the overallphase o ff set. We then compute the correction to the Solar SystemBarycenter at these new values. This is simply because the correc-tion due to the Roemer delay can vary by ∼
100 ns across the set oftimes t i that span the pulse period of the pulsar, which is significantin the context of the highest precision observations available today.We can then write the likelihood that the data is described onlyby the shapelet parameters, θ ≡ ( A , ζ, Λ ), the timing model param-eters (cid:15) , and the baseline o ff set, γ i , for each epoch i as:Pr( d | θ, (cid:15) ) ∝ N d (cid:89) i = √ det N i (9) × exp (cid:34) −
12 ( d − s i ( θ, (cid:15) ) − γ i ) T N − i ( d − s i ( θ, (cid:15) ) − γ i ) (cid:35) , where N i is the white noise covariance matrix for epoch i , withelements ( N i ) jk = σ i δ jk . In Eq. (9) the amplitude of the radiometer noise in the profile data istaken to be a known quantity at each epoch. In principle we wouldwant to include σ i as a free parameter for each epoch, howeverthis is not currently computationally tractable for large data sets. Inour analysis we can obtain a good initial estimate of the radiometernoise from the o ff -model region of the profile data, where the am-plitude of the model profile is less than 0 . α , which scales our initial estimate of the radiometernoise for all epochs associated with a given observing system.Our value for the radiometer noise then becomes:ˆ σ ij = ( α i σ ij ) . (10) We define profile jitter in our model as a shift in the arrival timerelative to those predicted by the pulsar’s timing model of the deter-ministic profile that is uncorrelated between di ff erent observationalepochs. This could be the result of systematic e ff ects, or intrinsichigh frequency stochasticity in the emission time of the pulse (seeFig. 3, e.g. Shannon et al. (2014)). In the latter case, even if we as-sume that the mean profile is the same at each epoch, the observedarrival time could still be subject to a shift due to this stochasticityin emission (see Fig. 3). Regardless of its origin, in the TOA do-main this is typically modelled with an EQUAD parameter, whichadds in quadrature to the formal TOA uncertainty. In the profile do-main we can directly incorporate a model that describes a stochasticshift in the mean profile. In our model we assume that the jitter am-plitudes at each epoch can be described by a Gaussian distribution,where the variance of that distribution is a free parameter in ouranalysis. c (cid:13)000
12 ( d − s i ( θ, (cid:15) ) − γ i ) T N − i ( d − s i ( θ, (cid:15) ) − γ i ) (cid:35) , where N i is the white noise covariance matrix for epoch i , withelements ( N i ) jk = σ i δ jk . In Eq. (9) the amplitude of the radiometer noise in the profile data istaken to be a known quantity at each epoch. In principle we wouldwant to include σ i as a free parameter for each epoch, howeverthis is not currently computationally tractable for large data sets. Inour analysis we can obtain a good initial estimate of the radiometernoise from the o ff -model region of the profile data, where the am-plitude of the model profile is less than 0 . α , which scales our initial estimate of the radiometernoise for all epochs associated with a given observing system.Our value for the radiometer noise then becomes:ˆ σ ij = ( α i σ ij ) . (10) We define profile jitter in our model as a shift in the arrival timerelative to those predicted by the pulsar’s timing model of the deter-ministic profile that is uncorrelated between di ff erent observationalepochs. This could be the result of systematic e ff ects, or intrinsichigh frequency stochasticity in the emission time of the pulse (seeFig. 3, e.g. Shannon et al. (2014)). In the latter case, even if we as-sume that the mean profile is the same at each epoch, the observedarrival time could still be subject to a shift due to this stochasticityin emission (see Fig. 3). Regardless of its origin, in the TOA do-main this is typically modelled with an EQUAD parameter, whichadds in quadrature to the formal TOA uncertainty. In the profile do-main we can directly incorporate a model that describes a stochasticshift in the mean profile. In our model we assume that the jitter am-plitudes at each epoch can be described by a Gaussian distribution,where the variance of that distribution is a free parameter in ouranalysis. c (cid:13)000 , 000–000 rofile Stochasticity in PSR J1909 − A r b i t r a r y F l u x Mean Profile Phase δ t Phase
Figure 3. (Bottom) A series of simulated single pulses with a Gaussianprofile, and a high frequency stochastic element to the pulse emission time.(Top) The mean profile that results from averaging the single pulses. Evenif we assume that the mean profile is the same at each epoch, the observedarrival time could still be subject to a shift due to this stochasticity in emis-sion. In the TOA domain this is typically modelled with an EQUAD param-eter, which adds in quadrature to the formal TOA uncertainty. In the profiledomain we can directly incorporate a model that describes a stochastic shiftin the mean profile.
If we write the amplitude of the shift due to profile jitter atany epoch i as δ t i , we can modify Eq. 8 trivially to include this newparameter as: s i ( t , A , ζ, Λ , (cid:15), δ t i ) = A n max (cid:88) n = ζ n B n ( t − τ i ( (cid:15) ) − δ t i ; Λ ) . (11)If this was the only modification to our analysis method, clearlyintroducing a shift parameter per epoch would result in completecovariance with the timing model parameters, it would be equiv-alent to introducing a jump in between each TOA. We will laterinclude a prior on these shift amplitudes that constrains them tocome from a Gaussian distribution, and will show in the context ofsimulations, that when the jitter model is appropriate, this method produces results that are completely consistent with the simulatedjitter amplitudes, and constraints that are superior to those obtain-able in the TOA domain, with no loss in timing precision.Written as in Eq. 11 we have to include N d free parametersdescribing the jitter in our analysis, which, as for the overall pro-file amplitudes would rapidly become intractable for large datasets. An advantage of the shapelet basis, is that for any arbitraryshapelet model we can analytically expand Eq. 11 for small ampli-tudes δ t i relative to the width parameter Λ . This allows us to lin-earise the shift, and so marginalise analytically over the individualamplitudes, a process we will describe in Section 4. A completederivation of this linear translation operator, along with operatorsfor other transformations, can be found in Refregier (2003). Theresult is that, for any profile model s i ( t , ζ, Λ , (cid:15) ), we can write downa ‘jitter profile’ as: j i ( t , δ t i , ζ, Λ , (cid:15) ) = δ t i √ Λ n max (cid:88) n = ζ n (cid:16) √ nB n ( t − τ i ( (cid:15) ); Λ ) (12) − √ n + B n + ( t − τ i ( (cid:15) ); Λ ) (cid:17) . (13)For example, in the case where n max =
0, such that we onlyinclude the Gaussian term in our model, our jitter profile will begiven by: j i ( t , δ t i , ζ, Λ , (cid:15) ) = δ t i √ Λ ζ ( − B ( t − τ i ( (cid:15) ); Λ )) (14) = − δ t i Λ ζ ( t − τ i ( (cid:15) )) exp (cid:32) −
12 ( t − τ i ( (cid:15) )) Λ (cid:33) (15)which is equivalent to expanding exp (cid:18) −
12 ( t − τ i ( (cid:15) ) − δ t i ) Λ (cid:19) for small δ t i .We illustrate this example in Fig. 4.Our shifted shapelet model will then be given by: s (cid:48) i ( t , A , δ t i , ζ, Λ , (cid:15) ) = s i ( t , A , ζ, Λ , (cid:15) ) + j i ( t , δ t i , ζ, Λ , (cid:15) ) . (16)As mentioned previously, we must include a constraint on theamplitudes of the shift parameters δ t . We do this by defining a prioron the amplitude parameters that is equivalent to fitting for the vari-ance of the distribution. This variance can be defined across an ar-bitrary set of observational epochs, for example, to model the vari-ance separately as a function of observing system, or for all epochssimultaneously.The covariance matrix of the jitter amplitudes, which we de-note J , is written: J ij = (cid:68) δ t i δ t j (cid:69) = J i δ ij , (17)where the set of coe ffi cients J represent the theoretical variance ofthe jitter model at epoch i .We normalise the amplitudes δ t i at each epoch such that theycan be described by a single variance J i for all observations inunits of seconds. We cannot simply use the maximum likelihoodamplitude of the profile model however, as if we believe our data tobe a ff ected by jitter, this amplitude will not be correct, which willdecrease our sensitivity to the amplitude of jitter in the data set.We therefore define the scaling factor C i such that: C i = F d , i F tot , i , (18)where F tot , i is the total flux in the unscaled model profile at epoch i , and so is independent of the overall model amplitude, and F d , i is an estimate of the total signal flux in the data at epoch i , which c (cid:13) , 000–000 L. Lentati et al.
Figure 4.
Example of the linear translation operator for the lowest ordershapelet coe ffi cient. The black solid curve represents the model profile atthe arrival time predicted by the pulsar’s timing model. The red line is thejitter profile defined in Eq. (12). Adding the jitter profile to the deterministicprofile results in a shift by an amount δ t (black dotted line). We stress thatin addition to the shift amplitudes representing jitter, we enforce a strongadditional constraint in the form of our prior in Eq. (24). All jitter ampli-tudes are constrained to come from a Gaussian distribution, whose varianceis free parameter in our analysis. We will show using simulations in Section7, that when the stochastic component of the data really is well described bya shift in the arrival time, this profile domain model gives completely con-sistent constraints to those obtained in the TOA domain using an EQUADparameter. we calculate by integrating over the bins where the model profile isgreater than 0.1% of its total amplitude.We can then write our joint probability density for the individ-ual jitter amplitudes, and the variance parameters as:Pr( d | θ, (cid:15) , δ t , J ) = Pr( d | θ, (cid:15) , δ t ) × Pr( δ t |J ) (19)where Pr( d | θ, (cid:15) , δ t ) has the same functional form as in Eq. 9, butwith the the shapelet profile model being replaced with our shiftedmodel in Eq. (16), and Pr( δ t |J ) is given by:Pr( δ t |J ) ∝ √ det J exp (cid:34) − δ t T J − δ t (cid:35) . (20) In addition to our jitter model, we also include epoch to epoch vari-ation in the shape of the deterministic profile. This could be the re-sult of systematic e ff ects, such as errors in polarization calibration,or due to intrinsic shape changes in the profile, where the radiome-ter noise for a particular observing system is below the level of theintrinsic stochasticity.We define this stochasticity in two regimes, high and low fre-quency fluctuations in the profile shape, where by frequency, werefer to the scale of the fluctuations in phase space. Any stochasticity in the profile that occurs on scales of less than thebin width in phase space will present itself in any one epoch simply as an increase in the uncorrelated noise in the on pulse region ofthe profile. In the most general case, we can define a ‘stochasticenvelope’, which represents the increase in the radiometer noisedue to this process such that: E i ( t , θ e , (cid:15) ) = (cid:114) N i T i C i s p ( t , θ e , (cid:15) ) , (21)where θ e are the parameters that describe the shape of the envelope, T i is the integration time for the observation at epoch i , and C i is thenormalisation constant defined in Eq. 18. We then add this stochas-tic profile to our instrument noise σ i such that the new variance ina bin j for observational epoch i is given by:ˆ σ ij = ( α i σ ij ) + E i ( t , θ e , (cid:15) ) j . (22)While this is the most general case, such that the shape of the en-velope is free to vary in the analysis, we will assume a simplermodel so that the stochastic envelope has the same shape as themean profile, and is proportional to the amplitude of the profilein any one observation. We note that giant pulses observed in bothcanonical pulsars, and MSPs (e.g. Wolszczan et al. (1984); Sallmenet al. (1999); Romani & Johnston (2001)) show structure on ∼ β which represents theincrease in the variance, and redefine E such that: E i ( t , θ, (cid:15) ) = β (cid:114) N i T i C i s p ( t , θ, (cid:15) ) . (23) To model low-frequency stochastic shape changes in the profile weuse the same shapelet basis as for the mean profile itself. Our goal isthen to robustly determine the power spectrum of the shape changesas a function of the scale (or order) of the component in the model.We therefore define the set of shape variation power spectrumcoe ffi cients S , such that for a given scale i the covariance matrix,which we denote S , will be given by: S ij = (cid:68) ( ζ i − ¯ ζ i )( ζ j − ¯ ζ j ) (cid:69) = P d , i S i δ ij , (24)where ¯ ζ i is the mean value of the shapelet coe ffi cient ζ i , which inour analysis we take to be the values that describe the determinis-tic profile. The normalising factor 1 / P d , i , with P d , i the power in theon pulse region at epoch i , is included to give us the power spec-trum coe ffi cients S i in units of the fraction of the total power in theprofile.In any one observational epoch i , we can therefore write thesum of the deterministic, and stochastic components of the profileas: s (cid:48) i ( t , ¯ θ, (cid:15), ζ ) = ¯ s i ( t , ¯ θ, (cid:15) ) + s i ( t , θ, (cid:15) ) . (25)We can then write our joint probability density for the individ-ual jitter amplitudes, and the variance parameters as:Pr( d | ¯ θ, ζ, (cid:15) , S ) = Pr( d | ¯ θ, (cid:15) , ζ ) × Pr( ζ |S ) (26) c (cid:13)000
Example of the linear translation operator for the lowest ordershapelet coe ffi cient. The black solid curve represents the model profile atthe arrival time predicted by the pulsar’s timing model. The red line is thejitter profile defined in Eq. (12). Adding the jitter profile to the deterministicprofile results in a shift by an amount δ t (black dotted line). We stress thatin addition to the shift amplitudes representing jitter, we enforce a strongadditional constraint in the form of our prior in Eq. (24). All jitter ampli-tudes are constrained to come from a Gaussian distribution, whose varianceis free parameter in our analysis. We will show using simulations in Section7, that when the stochastic component of the data really is well described bya shift in the arrival time, this profile domain model gives completely con-sistent constraints to those obtained in the TOA domain using an EQUADparameter. we calculate by integrating over the bins where the model profile isgreater than 0.1% of its total amplitude.We can then write our joint probability density for the individ-ual jitter amplitudes, and the variance parameters as:Pr( d | θ, (cid:15) , δ t , J ) = Pr( d | θ, (cid:15) , δ t ) × Pr( δ t |J ) (19)where Pr( d | θ, (cid:15) , δ t ) has the same functional form as in Eq. 9, butwith the the shapelet profile model being replaced with our shiftedmodel in Eq. (16), and Pr( δ t |J ) is given by:Pr( δ t |J ) ∝ √ det J exp (cid:34) − δ t T J − δ t (cid:35) . (20) In addition to our jitter model, we also include epoch to epoch vari-ation in the shape of the deterministic profile. This could be the re-sult of systematic e ff ects, such as errors in polarization calibration,or due to intrinsic shape changes in the profile, where the radiome-ter noise for a particular observing system is below the level of theintrinsic stochasticity.We define this stochasticity in two regimes, high and low fre-quency fluctuations in the profile shape, where by frequency, werefer to the scale of the fluctuations in phase space. Any stochasticity in the profile that occurs on scales of less than thebin width in phase space will present itself in any one epoch simply as an increase in the uncorrelated noise in the on pulse region ofthe profile. In the most general case, we can define a ‘stochasticenvelope’, which represents the increase in the radiometer noisedue to this process such that: E i ( t , θ e , (cid:15) ) = (cid:114) N i T i C i s p ( t , θ e , (cid:15) ) , (21)where θ e are the parameters that describe the shape of the envelope, T i is the integration time for the observation at epoch i , and C i is thenormalisation constant defined in Eq. 18. We then add this stochas-tic profile to our instrument noise σ i such that the new variance ina bin j for observational epoch i is given by:ˆ σ ij = ( α i σ ij ) + E i ( t , θ e , (cid:15) ) j . (22)While this is the most general case, such that the shape of the en-velope is free to vary in the analysis, we will assume a simplermodel so that the stochastic envelope has the same shape as themean profile, and is proportional to the amplitude of the profilein any one observation. We note that giant pulses observed in bothcanonical pulsars, and MSPs (e.g. Wolszczan et al. (1984); Sallmenet al. (1999); Romani & Johnston (2001)) show structure on ∼ β which represents theincrease in the variance, and redefine E such that: E i ( t , θ, (cid:15) ) = β (cid:114) N i T i C i s p ( t , θ, (cid:15) ) . (23) To model low-frequency stochastic shape changes in the profile weuse the same shapelet basis as for the mean profile itself. Our goal isthen to robustly determine the power spectrum of the shape changesas a function of the scale (or order) of the component in the model.We therefore define the set of shape variation power spectrumcoe ffi cients S , such that for a given scale i the covariance matrix,which we denote S , will be given by: S ij = (cid:68) ( ζ i − ¯ ζ i )( ζ j − ¯ ζ j ) (cid:69) = P d , i S i δ ij , (24)where ¯ ζ i is the mean value of the shapelet coe ffi cient ζ i , which inour analysis we take to be the values that describe the determinis-tic profile. The normalising factor 1 / P d , i , with P d , i the power in theon pulse region at epoch i , is included to give us the power spec-trum coe ffi cients S i in units of the fraction of the total power in theprofile.In any one observational epoch i , we can therefore write thesum of the deterministic, and stochastic components of the profileas: s (cid:48) i ( t , ¯ θ, (cid:15), ζ ) = ¯ s i ( t , ¯ θ, (cid:15) ) + s i ( t , θ, (cid:15) ) . (25)We can then write our joint probability density for the individ-ual jitter amplitudes, and the variance parameters as:Pr( d | ¯ θ, ζ, (cid:15) , S ) = Pr( d | ¯ θ, (cid:15) , ζ ) × Pr( ζ |S ) (26) c (cid:13)000 , 000–000 rofile Stochasticity in PSR J1909 − where Pr( d | ¯ θ, (cid:15), ζ ) has the same functional form as in Eq. 9, butwith the the shapelet profile model being replaced by the sum ofthe mean, and stochastic profile models as in Eq. (25), and Pr( ζ |S )is given by:Pr( ζ |S ) ∝ √ det S exp (cid:34) − ζ T S − ζ (cid:35) . (27) Ignoring any additional stochastic parameters, our model for thepulse profile so far contains a set of n max − ff ectssuch as scintillation in the interstellar medium can cause both thee ff ective band-integrated shape, and the amplitude of the profile tovary significantly (Rickett 1969). As such, the overall amplitudemust be free to vary from one epoch to another, and in addition, thebaseline o ff set must also be treated as a separate parameter for eachepoch.When including stochastic parameters we then have an addi-tional amplitude parameter per epoch that describes the shift in thearrival time due to pulse jitter, as defined in Section 3.2, as wellas a set of amplitude parameters that describe the low frequencystochastic variation in the pulse shape for that epoch as describedin Section 3.3.2.In principle this introduces a large number of additional freeparameters, however we can marginalise analytically over all ofthese amplitude parameters without introducing large matrix op-erations and thus significantly decrease the size of parameter spacewe must sample over numerically.If the number of amplitude parameters we wish to marginaliseover analytically for any epoch i is N m we define a N m × N i ma-trix, which we denote P i such that, for example the profile overallamplitude parameter and o ff set amplitude would be included as:( P i ) , j = P i ) , j = ( s i ) j , (28)with the remaining basis vectors describing pulse jitter and profilevariation being included in rows ( P i ) .. N m . We then define the diago-nal, N m × N m matrix Ψ i , which combines our priors for the jitter andlow frequency stochastic parameters, and is zero for the elementscorresponding to the overall amplitude and baseline o ff set. Finallywe define the length- N m vector b i , which contains the amplitudeparameters that multiply our basis vectors in P i , which allows us towrite our likelihood as:Pr( d | ¯ θ, (cid:15), δ t , ζ )Pr( ζ |S )Pr( δ t |J ) ∝ N d (cid:89) i = √ det N i det Ψ i × (29)exp (cid:34) −
12 ( d i − P i b i ) T N − i ( d i − P i b i ) (cid:35) × exp (cid:34) − b i T Ψ i − b i (cid:35) . In order to perform the marginalisation over the amplitudes b i , we define P T i N − i P i + Ψ i − as Σ i and P T i N − i d i as ¯d i and then integrate over the elements in b i analytically to get the likelihoodthat the data is described by the remaining parameters alone.Our marginalised probability distribution for the remainingparameters alone is then given by:Pr( ¯ θ, (cid:15), S , J| d ) ∝ N d (cid:89) i = det ( Σ i ) − √ det ( N i ) × exp (cid:34) − (cid:16) d T i N − i d i − ¯d T i Σ − i ¯d i (cid:17)(cid:35) . (30) Bayesian Inference provides a consistent approach to the estimationof a set of parameters Θ in a model or hypothesis H given the data, D . Bayes’ theorem states that:Pr( Θ | D , H ) = Pr( D | Θ , H )Pr( Θ | H )Pr( D | H ) , (31)where Pr( Θ | D , H ) ≡ Pr( Θ ) is the posterior probability distributionof the parameters, Pr( D | Θ , H ) ≡ L ( Θ ) is the likelihood, Pr( Θ |H ) ≡ π ( Θ ) is the prior probability distribution, and Pr( D | H ) ≡ Z is the Bayesian Evidence.For model selection, the evidence is key, and is simply thefactor required to normalise the posterior over Θ : Z = (cid:90) L ( Θ ) π ( Θ )d n Θ , (32)where n is the dimensionality of the parameter space.As the average of the likelihood over the prior, the evidenceis larger for a model if more of its parameter space is likely andsmaller for a model where large areas of its parameter space havelow likelihood values, even if the likelihood function is very highlypeaked. Thus, the evidence automatically implements Occam’s ra-zor: a simpler theory with a compact parameter space will have alarger evidence than a more complicated one, unless the latter issignificantly better at explaining the data.The question of model selection between two models H and H can then be decided by comparing their respective posteriorprobabilities, given the observed data set D , via the posterior oddsratio R : R = P ( H | D ) P ( H | D ) = P ( D | H ) P ( H ) P ( D | H ) P ( H ) = Z Z P ( H ) P ( H ) , (33)where P ( H ) / P ( H ) is the a priori probability ratio for the twomodels, which can often be set to unity but occasionally requiresfurther consideration.The posterior odds ratio then allows us to obtain the probabil-ity of one model compared with the other, simply as: P = R + R . (34)Typically in our analysis we deal with the log odds ratio,which is then simply the di ff erence in the log Evidence for twocompeting models. Table 1 lists a common interpretation of theBayes factor given by Kass & Raftery (1995). In all our analysiswe will consider a di ff erence in the log Evidence between modelsof 3 to be the minimum required to justify the addition of extramodel components. c (cid:13) , 000–000 L. Lentati et al.
Phase A r b i t r a r y F l u x Phase A r b i t r a r y F l u x PDFB1PDFB2 PDFB4WBCORR
Figure 5. (Left) Post-fit residuals for PSR J1909 − χ of 1.7. Colours represent di ff erent observing systems: (black)WBCORR, (red) PDFB1, (green) PDFB2, (dark blue) PDFB4, (light blue) PDFB4a. (Right) Example profile data for PSR J1909 − ff erentobservational epochs, each with a di ff erent observing system. The di ff erent overall fluxes, o ff sets and noise levels in each are apparent. Table 1.
Interpretation of the evidence R log( R ) P Strength of evidence1 → → . → .
75 Not worth more than a bare mention3 →
20 1 → . → .
95 Positive20 →
150 3 → . → .
99 Strong > > > .
99 Very strong
The nested sampling approach (Skilling 2004) is a Monte-Carlomethod targeted at the e ffi cient calculation of the evidence, butalso produces posterior inferences as a by-product. Recently a newnested sampling algorithm, P olychord (Handley et al. 2015), wasintroduced that makes evidence calculation in hundreds of dimen-sions a tractable process. The P olychord algorithm makes use of‘slice sampling’ (Neal 2000). In one dimension, given a likelihood L , a point x is within the ‘slice’ if L ( x ) > L . Starting from a seedpoint x , sampling boundaries are set by expanding a random initialbound of width w until L ( x ± w ) < L . A new point x is then ob-tained within the slice by sampling uniformly within these bounds.If x is not in the slice, w is decreased and x is drawn again fromthese new bounds.In d dimensions, P olychord first ‘whitens’ the parameterspace, performing a linear skew transformation which turns degen-erate contours in the original parameter space into into ones withdimensions O(1) in all directions. An initial live point is then cho-sen randomly with coordinates in this whitened space given by x .A random initial direction is then chosen ˆn , and one dimensionalslice sampling is performed in that direction to generate a new point x . This process is repeated O ( d ) times to generate a new uniformlysampled point x N which is decorrelated from the initial point x .In the analysis presented in this work we will be dealing withmodels that have up to ∼
100 free parameters, as such P olychord isa vital tool for performing robust model comparison in the profiledomain.
We perform our analysis, and construct simulations that are basedo ff the 10cm data set for PSR J1909 − ff erent spectrometers.The improved computing power of the later spectrometers enabledthe implementation of polyphase digital filters, and an increasednumber of spectral channels and pulse phase bins. These data wereprocessed using procedures described in Manchester et al. (2013),which we will describe in brief below, in order to draw comparisonwith the methods presented in this work.Individual observations were first averaged fully in both timeand frequency to form a single pulse profile for each observation.Pulse TOAs were then measured by convolving model templates inthe Fourier Domain with the total intensity pulse profiles for eachobservational epoch. Model templates were constructed by sum-ming a series of von Mises functions, each of which is describedby 3-parameters, a phase, an amplitude, and a width. Initial guessesfor component locations, amplitudes and widths were added by eyeand then fitted using a non-linear algorithm. The residuals of thefit were inspected to determine if the fit was good or if additionalcomponents were needed, and the processes iterated upon until theresiduals appeared to have white-noise characteristics.Because of its relatively simple morphology (Dai et al. 2015)only three components were needed to model the pulse profile,however, as in this work, the weak interpulse was not included inthe model. Templates were formed for each set of systems that hadmarkedly di ff erent backend architectures, in order to deal with pro-file distortions that might arise from the time, and frequency sam-pling for the di ff erent systems. Pulse-profile features narrower thana pulse bin will be unresolved. Similarly, dispersive smearing canbroaden a pulse if frequency channels are wide (and, like the sys-tems used here, coherent dedispersion is not employed) relative tothe intra-channel dispersion, however, consistent timing solutions c (cid:13)000
We perform our analysis, and construct simulations that are basedo ff the 10cm data set for PSR J1909 − ff erent spectrometers.The improved computing power of the later spectrometers enabledthe implementation of polyphase digital filters, and an increasednumber of spectral channels and pulse phase bins. These data wereprocessed using procedures described in Manchester et al. (2013),which we will describe in brief below, in order to draw comparisonwith the methods presented in this work.Individual observations were first averaged fully in both timeand frequency to form a single pulse profile for each observation.Pulse TOAs were then measured by convolving model templates inthe Fourier Domain with the total intensity pulse profiles for eachobservational epoch. Model templates were constructed by sum-ming a series of von Mises functions, each of which is describedby 3-parameters, a phase, an amplitude, and a width. Initial guessesfor component locations, amplitudes and widths were added by eyeand then fitted using a non-linear algorithm. The residuals of thefit were inspected to determine if the fit was good or if additionalcomponents were needed, and the processes iterated upon until theresiduals appeared to have white-noise characteristics.Because of its relatively simple morphology (Dai et al. 2015)only three components were needed to model the pulse profile,however, as in this work, the weak interpulse was not included inthe model. Templates were formed for each set of systems that hadmarkedly di ff erent backend architectures, in order to deal with pro-file distortions that might arise from the time, and frequency sam-pling for the di ff erent systems. Pulse-profile features narrower thana pulse bin will be unresolved. Similarly, dispersive smearing canbroaden a pulse if frequency channels are wide (and, like the sys-tems used here, coherent dedispersion is not employed) relative tothe intra-channel dispersion, however, consistent timing solutions c (cid:13)000 , 000–000 rofile Stochasticity in PSR J1909 − were found if one template was used for all the backends, and at 10cm this pulsar is not a ff ected by dispersive smearing.Between backends discrete phase o ff sets (referred to as jumps)occur because of di ff erent propagation times through the digitalsystems. These o ff sets can also occur within a backend when it,or the observatory clock-distribution system experiences significanthardware changes. While it is possible in principle to measure rel-ative o ff sets between backends independent of the pulsar observa-tion, (by for example injecting a common signal into multiple back-ends and measuring the delay), these systems have yet to demon-strate the accuracy necessary to correct the highest-precision timingobservations. As a result, we include as free parameters in our anal-ysis o ff sets between all of the backends, and, for the Parkes DigitalFilterbank (PDFB) 4, at MJD 55319 when hardware changes re-defined the location within the signal processing at which a time-stamp was assigned, resulting in the definition of the PDFB4a sys-tem after this jump date. We construct four simulations using a fiducial timing modelfor PSR J1909 − ff set. We then include di ff erent levels of the radiometer noise,systematic pulse jitter, and additional non-stationary profile com-ponents to the simulations to test the e ffi cacy of our profile domainjitter and stochasticity models described in Sections 3.2 and 3.3.For each simulation we process the profile data in the same way asdescribed in Section 6 to form TOAs against which to compare ourprofile domain analysis.In simulations 1 and 2 we include radiometer noise consistentwith the levels determined from the real data set. In simulation 2we then include an additional systematic jitter term to the WideBand Correlator (WBCORR, Manchester et al. 2013) profiles withan amplitude of 5 × − s, chosen to be consistent with the value ob-served in the real data set when performing a TOA domain analysis.In simulations 3 and 4 we include a factor 10 times less radiometernoise than the real data set, and in simulation 4 we include an addi-tional Gaussian component to the profile, that has a width 0 .
1% thepulse period, and 1% the amplitude of the profile for any given ob-servational epoch. The component moves deterministically throughthe main pulse profile over the course of 4000 days. We use this forour simulation rather than including stochasticity in the same basisas our model in order to check the e ffi cacy of our method when itis attempting to correct for epoch to epoch profile variations thatare not well described by zero mean Gaussian fluctuations in theamplitude of the profile components.In Fig. 6 we show the residuals for all simulations after sub-tracting the fiducial timing model. The additional jitter in the WB-CORR system is apparent in the early data for simulation 2, as isthe timing noise induced by the deterministic passage of an addi-tional Gaussian component through the pulse profile for simulation4. In order to compare the profile and TOA domain analysis for sim-ulations 1 and 2, we define the following models: • TOA Model 1: Timing model only. • TOA Model 2: Timing model and additional EFAC andEQUAD parameters for each observing system. • Profile Model 1: Deterministic profile and timing model only. • Profile Model 2: Deterministic profile, timing model and ad-ditional EFAC and jitter parameters for each observing system.In both profile domain models we include only the n = ffi cient in our analysis. We find that the timing modelparameter estimates and uncertainties for equivalent models in bothsimulations are completely consistent, as we would expect, as theassumptions about the stochastic properties in the profile domainwere propagated correctly into the TOA analysis.In simulation 2 we find that the profile domain jitter modelresults in superior constraints on the high frequency variation inpulse arrival times compared to the TOA domain EQUAD model.In Fig. 7 (Top) we show a simulated profile for a single observa-tional epoch from simulation 2, for the WBCORR system, to whichwe applied the additional systematic jitter signal. In the left panelwe show the simulated profile (blue line), the model profile, evalu-ated at the arrival time predicted by the maximum likelihood timingmodel from our analysis (green line), and the maximum likelihoodjitter profile (red line). In the right panel we have subtracted themodel profile from the data, and the residual is clearly consistentwith our jitter profile.In Fig. 7 (bottom left) we compare the one and two-dimensional posterior parameter estimates for the EFAC andEQUAD parameters for the WBCORR system obtained in the TOAdomain analysis (solid lines in the one-dimensional plots, solidcontours in the two-dimensional plot), with the EFAC and jitter pa-rameter estimates from our profile domain analysis (dashed linesin the one-dimensional plots, filled dashed contours in the two-dimensional plot). The most striking feature of this comparisonis the di ff erence in constraints placed on the EFAC parameters inboth domains. In the profile domain we still have the o ff pulse re-gion to constrain our estimates of the radiometer noise, informa-tion that is no longer available when we move to the TOA domain,leading to an order of magnitude decrease in precision when theEFAC parameter is estimated as part of the analysis. As a by prod-uct the constraints on the jitter amplitude in the profile domain areimproved over the TOA domain EQUAD counterpart. This is be-cause in the TOA domain, the EFAC and EQUAD parameters arecorrelated (indicated by the ellipsoidal contours). In the profile do-main this is not the case, as seen in the bottom right panel, the jitteramplitude, and EFAC parameters are completely decorrelated, as ashift in the arrival time of the profile does not look like an increasein the radiometer noise, however in the TOA domain, a multiplica-tive prefactor (EFAC), and an additional term added in quadrature(EQUAD), can be completely covariant depending on how similarthe error bars are across di ff erent observations. For simulations 3 and 4 we consider the following models: • TOA Model: Timing model, EQUAD parameters for each ob-serving system, and additional red timing noise. • Profile Model 1: Deterministic profile, timing model and ad-ditional red timing noise. • Profile Model 2: Deterministic profile, timing model, high andlow frequency stochastic profile parameters, and additional red tim-ing noise. c (cid:13) , 000–000 L. Lentati et al.
Figure 6.
Post-fit residuals for PSR J1909 − × − s, chosen to beconsistent with the value observed in the real data set when performing a TOA domain analysis. Simulations 3 and 4 include a factor 10 times less radiometernoise than the real data set, and in simulation 4 we include an additional Gaussian component to the profile, that has a width 0 .
1% the pulse period, 1% theamplitude of the profile for any given observational epoch, and which moves deterministically through the main pulse profile over the course of 4000 days.
Additionally, for both simulations and for all 3 models we ob-tain upper limits on the amplitude of a steep red noise process witha spectral index of γ = /
3, consistent with that expected for agravitational wave background formed from a population of inspi-ralling super-massive black hole binaries, (Phinney 2001).As for simulations 1 and 2, we find that for both simulations3 and 4 that the timing model parameters are consistent betweenall 3 models, indicating that, as for the jitter parameters includedin simulations 1 and 2, including the stochastic profile parametersin our analysis does not adversely a ff ect the precision that we canachieve when no such stochasticity is present in the data. In addi-tion, in Fig. 8 (left plot) we find the limits obtained on a steep rednoise process within the data set are also consistent between thethree models, with 95% upper limits of A < × − in all cases.For simulation 4, however, we find that the profile domainanalysis results in a factor 2 improvement in upper limit when in-cluding profile stochasticity, compared to either the standard TOAdomain analysis or the equivalent analysis using Model 1 in the pro-file domain. We show this explicitly in Fig. 8 (right plot), in whichwe plot the upper limits for the same 3 models applied to simulation 4. While TOA Model 1 (black line), and Profile Model 1 (red line)are consistent, with 95% upper limits of A < × − , we find afactor 2 improvement in the upper limit when incorporating profilestochasticity in our analysis (blue line), obtaining A < × − .The only di ff erence between simulations 3 and 4 was the inclusionof an additional non-stationary Gaussian component to the profilemodel, that moved through the profile over a period of 4000 days.Despite having an amplitude of only 1% that of the profile inany observational epoch, the e ff ect on our TOA domain analysiswas to increase the upper limit by an order of magnitude. Whileour stochastic profile domain model does not fully recover the up-per limit obtained from simulation 3, we purposefully simulatednon-stationary behaviour in the pulse profile that was not well de-scribed by a random Gaussian process. Despite this, we still obtaina significantly improved limit compared to the TOA domain analy-sis. The key result both from this comparison, and from the resultsof the jitter analysis in simulations 1 and 2, is that by performingthe analysis using the profile data, rather than the SATs as is typ-ical in pulsar timing analysis, we can decorrelate profile variation c (cid:13)000
3, consistent with that expected for agravitational wave background formed from a population of inspi-ralling super-massive black hole binaries, (Phinney 2001).As for simulations 1 and 2, we find that for both simulations3 and 4 that the timing model parameters are consistent betweenall 3 models, indicating that, as for the jitter parameters includedin simulations 1 and 2, including the stochastic profile parametersin our analysis does not adversely a ff ect the precision that we canachieve when no such stochasticity is present in the data. In addi-tion, in Fig. 8 (left plot) we find the limits obtained on a steep rednoise process within the data set are also consistent between thethree models, with 95% upper limits of A < × − in all cases.For simulation 4, however, we find that the profile domainanalysis results in a factor 2 improvement in upper limit when in-cluding profile stochasticity, compared to either the standard TOAdomain analysis or the equivalent analysis using Model 1 in the pro-file domain. We show this explicitly in Fig. 8 (right plot), in whichwe plot the upper limits for the same 3 models applied to simulation 4. While TOA Model 1 (black line), and Profile Model 1 (red line)are consistent, with 95% upper limits of A < × − , we find afactor 2 improvement in the upper limit when incorporating profilestochasticity in our analysis (blue line), obtaining A < × − .The only di ff erence between simulations 3 and 4 was the inclusionof an additional non-stationary Gaussian component to the profilemodel, that moved through the profile over a period of 4000 days.Despite having an amplitude of only 1% that of the profile inany observational epoch, the e ff ect on our TOA domain analysiswas to increase the upper limit by an order of magnitude. Whileour stochastic profile domain model does not fully recover the up-per limit obtained from simulation 3, we purposefully simulatednon-stationary behaviour in the pulse profile that was not well de-scribed by a random Gaussian process. Despite this, we still obtaina significantly improved limit compared to the TOA domain analy-sis. The key result both from this comparison, and from the resultsof the jitter analysis in simulations 1 and 2, is that by performingthe analysis using the profile data, rather than the SATs as is typ-ical in pulsar timing analysis, we can decorrelate profile variation c (cid:13)000 , 000–000 rofile Stochasticity in PSR J1909 − −0.7 −0.3 0.1 −6.4 −6.3 −6.2 −6.1Log EQUADLog EFAC
Log E Q U A D −0.7 −0.3 0.1−6.4−6.3−6.2−6.1 −0.005 0 0.005 −6.4 −6.3 −6.2 −6.1Log Jitter AmpLog EFAC
Log J i tt e r A m p −0.005 0 0.005−6.4−6.3−6.2−6.1 Figure 7. (Top left) Simulated profile (blue line) for a single observational epoch from simulation 2, for the WBCORR system, to which we applied theadditional jitter signal. Over plotted are the maximum likelihood model profile, evaluated at the arrival time predicted by the maximum likelihood timingmodel from our analysis (green line), and the maximum likelihood jitter profile (red line). (Top right) As in the previous panel, however we have subtractedthe model profile from the data, leaving a residual that is clearly consistent with our jitter model. (Bottom left) One and two-dimensional posterior parameterestimates for the EFAC and EQUAD parameters for the WBCORR system obtained in the TOA domain analysis (solid lines in the one-dimensional plots, solidcontours in the two-dimensional plot), with the EFAC and jitter parameter estimates from our profile domain analysis (dashed lines in the one-dimensionalplots, filled dashed contours in the two-dimensional plot). The uncertainty in the EFAC parameter in the profile domain is an order of magnitude smallerthan in the TOA domain, as we can use the statistical information present in the o ff pulse region to help constrain our estimates of the radiometer noise. Theconstraints on the jitter amplitude in the profile domain are also improved compared to the TOA domain EQUAD counterpart. This is because in the TOAdomain, the EFAC and EQUAD parameters are correlated (indicated by the ellipsoidal contours). In the profile domain this is not the case, as seen in thebottom right panel, the jitter amplitude, and EFAC parameters are completely decorrelated in our profile domain analysis. from the GW signal and significantly improve both our sensitivityto GWs, and our understanding of the stochastic processes presentin the data. As models for profile variation improve, so too willthe upper limits on GWs obtainable from a profile domain analy-sis. Given the covariance between the signal induced in the timingresiduals from GWs, and from long-term temporal profile variation(see Fig. 2), however, only by performing a simultaneous fit can ro-bust estimates of both signal components be obtained, without oneabsorbing the other.We find the amplitude of both the high and low frequencystochastic models are consistent with the properties of the addi-tional profile component added to simulation 4. In Fig. 9 we showthe one dimensional marginalised posterior parameter estimates forboth the high frequency (top panel) and low frequency (bottompanel) stochastic profile parameters from our analysis of simula-tion 4. We included separate scaling factors, β , describing the highfrequency profile variation for each observing system. Values are quoted as the fraction of the total profile amplitude added to theuncorrelated variance for each bin. For all systems we find the frac-tional increase in the uncorrelated noise is < ff erent systems varies by a factor ∼ ff erent period of time (See Fig. 5)this could be the result of the non-stationary nature of the additionalGaussian component added to the profile data.For the low frequency profile variations we show either themean parameter estimates and 1 σ confidence intervals (points witherror bars), or 95% upper limits (arrows). Upper limits are plottedwhere the posterior for that scale is consistent with zero at greaterthan 5% probability. We find that only stochastic variations in the n = (2 , ,
4) profile coe ffi cients are significant. Given the widthof the profile, Λ = .
007 in phase, which corresponds to ∼ n = (cid:54) c (cid:13) , 000–000 L. Lentati et al. P r obab ili t y GWB Amplitude
0 5e-16 1e-15 1.5e-15 2e-15 P r obab ili t y GWB Amplitude
Figure 8.
One-dimensional marginalised posterior parameter estimates for the amplitude of a steep red noise process with a spectral index of γ = / × − . In simulation 4 we find that the TOA domain analysis, and a profile domain analysis that includes a red noise signal but no profile stochasticityresults in consistent 95% upper limits of 1 × − . Despite the fact that the additional Gaussian component was only 1% of the profile amplitude, and 0.1%of the pulse width, it has a significant impact on the upper limit, increasing it by an order of magnitude. Our profile domain analysis that includes profilestochasticity (blue dashed line) significantly improves upon the TOA domain analysis, with a 95% upper limit of 5 × − . While this is still much greaterthan for the data set that included no profile variation it still demonstrates clearly that by modelling these e ff ects in the profile domain they can be betterisolated from the GWB signal of interest. high frequency stochastic parameters, and thus less significant inour analysis. − We now apply our profile domain framework to the PSRJ1909 − • TOA Model 1: Timing model only. • TOA Model 2: Timing model, and additional EFAC andEQUAD parameters for each observing system. • Profile Model 1: Deterministic profile, and timing model only. • Profile Model 2: Deterministic profile, timing model, and ad-ditional EFAC and jitter parameters for each observing system. • Profile Model 3: Deterministic profile, timing model, and ad-ditional EFAC, jitter and high frequency stochastic profile parame-ters for each observing system. • Profile Model 4: Deterministic profile, timing model, EFACparameters for each observing system and high and low frequencystochastic profile parameters.For each model we find the optimal number of shapelet coe ffi -cients to include by increasing the number of included coe ffi cientsin steps of 5 until the change in the log evidence is less than 3. Table2 lists the optimal number of coe ffi cients included in each model,along with the log evidence values for that model.We find that Profile model 1, which includes no additionalstochastic parameters requires the largest number of components(40), while subsequent models that include either profile jitter orprofile stochasticity require 25% fewer components to describe thedata. Initially this may seem like a large number compared to the3 von Mises functions used in the standard analysis, however in Table 2.
Evidence values for di ff erent models fit to the PSR J1909 − c log Z TOA Model 1 - 0TOA Model 2 - 35.7Profile Model 1 40 0Profile Model 2 30 154.2Profile Model 3 30 1302.7Profile Model 4 30 1280.5 this case a separate profile model is fit to the WBCORR, PDFB1,PDFB2, and PDFB4 systems. Given a set of n von Mises functionsrequire 3 n -1 parameters ( n amplitudes, n widths and n − c (cid:13)000
Evidence values for di ff erent models fit to the PSR J1909 − c log Z TOA Model 1 - 0TOA Model 2 - 35.7Profile Model 1 40 0Profile Model 2 30 154.2Profile Model 3 30 1302.7Profile Model 4 30 1280.5 this case a separate profile model is fit to the WBCORR, PDFB1,PDFB2, and PDFB4 systems. Given a set of n von Mises functionsrequire 3 n -1 parameters ( n amplitudes, n widths and n − c (cid:13)000 , 000–000 rofile Stochasticity in PSR J1909 − -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 P r obab ili t y Log [ High Frequency Stochastic Amplitude ] WBCORRPDFB1PDFB2PDFB4PDFB4a -4.5-4-3.5-3-2.5-2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Log F r a c t i ona l P r o f il e S t o c ha s t i c i t y Log Profile Component
Figure 9. (Top) One-dimensional posterior parameter estimates for the am-plitude of the high frequency stochastic element to the profile for Simulation4. Separate amplitudes were fit to each observing system. The units of theabscissa are fraction of the total profile amplitude added to the uncorrelatedvariance for each bin. The amplitudes are consistent with <
1% of the pro-file amplitude for all systems however they show a factor ∼ ff erent period of time (SeeFig. 5) this could be the result of the non-stationary nature of the additionalGaussian component added to the profile data. (Bottom) Mean values and1 σ confidence intervals (red points with error bars), and 95% upper limits(black arrows) for low frequency stochastic variations in the profile. Upperlimits are plotted where the posterior for that scale is consistent with zeroat greater than 5% probability. related instrumental noise for the WBCORR system, however, isnot consistent with a value of 1 in profile model 3, implying thatthe residuals are not consistent with our initial estimate of the in-strumental noise for this backend. This can be seen more clearlyin the bottom panel of Fig. 10, which shows a comparison of theone-dimensional marginalised posterior parameter estimates for theWBCORR EFAC parameter from profile models 2 (black line) and3 (red line). Assuming that the mean profile su ff ers only a shift fromepoch to epoch as in profile model 3 (and in the standard timinganalysis paradigm) results in an EFAC that is systematically greaterthan 1, implying that the model used is not su ffi cient to describe the −0.7 −0.3 0.1 −7.5 −7 −6.5WBCORR Log EQUADWBCORR Log EFAC W B C O RR Log E Q U A D −0.7 −0.3 0.1−7.4−7.2−7−6.8−6.6−6.4 -0.01 0 0.01 0.02 0.03 0.04 P r obab ili t y WBCORR Log EFAC
Figure 10. (Top) One and two-dimensional posterior parameter estimatesfor the EFAC and EQUAD parameters for the WBCORR system obtainedin the TOA domain analysis (solid lines in the one-dimensional plots, solidcontours in the two-dimensional plot), with the EFAC and jitter parame-ter estimates from our profile domain analysis (dashed lines in the one-dimensional plots, filled dashed contours in the two-dimensional plot). Asin simulation 2, we find that the inferred amplitude of the profile domainjitter model was consistent with the TOA domain EQUAD model. We sim-ilarly see that as before the EFAC parameters are better constrained inthe profile domain by approximately an order of magnitude, with no ob-served correlation with the jitter parameters, leading to better overall con-straints being placed on their amplitudes. (Bottom) Comparison of the one-dimensional marginalised posterior parameter estimates for the WBCORREFAC parameter from profile models 2 (black line) and 3 (red line). As-suming that the mean profile su ff ers only a shift from epoch to epoch as inprofile model 3 (and in the standard timing analysis paradigm) results in anEFAC that is systematically greater than 1, implying that the model usedis not su ffi cient to describe the on pulse region. This is confirmed by theBayesian evidence, which is significantly greater ( ∆ log Z > on pulse region. This is confirmed by the Bayesian evidence, whichis significantly greater ( ∆ log Z > c (cid:13) , 000–000 L. Lentati et al.
Bin Number A r b i t r a r y F l u x Bin Number A r b i t r a r y F l u x Figure 11.
Example of two WBCORR profiles that su ff er from non-Gaussian ‘sawtooth’ noise. Approximately 10% of WBCORR points showthis e ff ect by eye, however the amplitude varies from profile to profile. Inthe TOA domain this manifests itself as an EQUAD term – high frequencyvariation in the arrival times of the pulses. In the profile domain we find theBayesian evidence is significantly in support ( ∆ log Z > ff ectsand separate them from ‘true’ shifts in the arrival time due to pulse jitter. ff er from non-Gaussian ‘sawtooth’ noise. Approxi-mately 10% of WBCORR points show this e ff ect by eye, howeverthe amplitude varies from profile to profile. While clearly the op-timal solution would be to simply fix such issues at the level offorming the folded data, where that is not possible by working inthe profile domain one can better model such systematic e ff ects andseparate them from ‘true’ shifts in the arrival time due to pulse jit-ter. Finally, we find no evidence for low frequency stochasticity inthe PSR J1909 − As a final comparison, we take TOA model 2, and profile model3 and include a red noise model for each. In both cases we ob-tain upper limits on the amplitude of a steep red noise process witha spectral index of γ = /
3. Fig. 13 shows the one-dimensionalmarginalised posterior parameter estimates for the amplitude of thered noise process for the TOA domain (black line), and profile do-main (red line) analysis. Both are consistent with one another, with95% upper limit of 1 × − . Given that we found no evidence forlow frequency profile stochasticity in the data, and only the olderWBCORR system, which contributes relatively little weight to the -2.2-2-1.8-1.6-1.4-1.2-1-0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Log F r a c t i ona l P r o f il e S t o c ha s t i c i t y Log Profile Component
Figure 12.
95% upper limits on the fractional profile stochasticity in thePSR J1909 − ∼ . P r obab ili t y Red Noise Amplitude
Figure 13.
One-dimensional marginalised posterior parameter estimates forthe amplitude of a red noise process in the PSR J1909 − / × − , consistent with values reported in S15. In this data set we foundno evidence of low frequency profile stochasticity, and only the older WB-CORR system was significantly a ff ected by high frequency systematics, soit is not surprising that upper limits using the optimal models for both theTOA and profile domain analyses are consistent with one another. full data set, was significantly a ff ected by high frequency system-atics, it is not surprising that upper limits using the optimal modelsfor both the TOA and profile domain analyses are consistent withone another. c (cid:13)000
One-dimensional marginalised posterior parameter estimates forthe amplitude of a red noise process in the PSR J1909 − / × − , consistent with values reported in S15. In this data set we foundno evidence of low frequency profile stochasticity, and only the older WB-CORR system was significantly a ff ected by high frequency systematics, soit is not surprising that upper limits using the optimal models for both theTOA and profile domain analyses are consistent with one another. full data set, was significantly a ff ected by high frequency system-atics, it is not surprising that upper limits using the optimal modelsfor both the TOA and profile domain analyses are consistent withone another. c (cid:13)000 , 000–000 rofile Stochasticity in PSR J1909 − In this paper we have extended the Bayesian framework ‘Genera-tive Pulsar Timing Analysis’ to incorporate both pulse jitter (highfrequency variation in the arrival time of the pulse) and epoch toepoch stochasticity in the shape of the pulse profile. This frame-work allows for a full timing analysis to be performed using thefolded profile data, rather than the site arrival times as is typical inmost timing studies.We first used simulations to show that while in a standard TOAdomain analysis of pulsar timing data the e ff ect of profile stochas-ticity and the signal induced by a gravitational wave backgroundare highly covariant, these two signals can be decorrelated whenperforming the analysis in the profile domain, resulting in a sig-nificantly improved upper limit on a stochastic gravitational wavebackground. We showed that this was also true for high frequencyvariation in the arrival time of the pulse. In the TOA domain thisis typically modelled using an ‘EQUAD’ parameter, that adds inquadrature to the TOA uncertainties. This is highly covariant with ascaling of the error bars (using ’EFACs’), which is typically used tomodel uncertainty in the value of the radiometer noise. In the pro-file domain a model for pulse jitter is completely decorrelated fromthe model for the thermal radiometer noise, resulting in improvedconstraints on the magnitude of the e ff ect in the pulsar timing data.Our simulation included temporal variation in the pulse profilein the form of an additional Gaussian component with an amplitude1% that of the main profile, and a width 0.1% of the pulse period.Despite the relatively small amplitude of the additional signal, thisresulted in an order of magnitude increase in the upper limit ob-tained in the TOA domain analysis. The techniques described inthis paper allow for the simultaneous evaluation of both contribu-tions to the data, resulting in robust upper limits, and eventual de-tection of a GWB.We also applied our methodology to a PSR J1909 − ff erent profile model components. Usingour profile domain framework we obtain upper limits on a red noiseprocess with a spectral index of γ = / × − , consistentwith previously published limits.With current upper limits for a stochastic GWB now begin-ning to rule out existing models for black hole binary populations,profile variability at the <
1% level can become an important limit-ing factor in the sensitivity achievable in pulsar timing experiments.By moving away from the standard timing paradigm, where a tim-ing model is fit to a set of site arrival times, and instead workingdirectly with the profile data, these e ff ects can be modelled appro-priately, and the covariance with a possible background, or othersignals of interest, can be significantly decreased. As upcoming ra-dio telescopes such as the Square Kilometre Array come on-line,improvements in sensitivity will increase the significance of thesee ff ects on the pulse arrival times as measured through traditionalmethods, making a profile domain analysis the natural choice in or-der to improve the prospects for the detection of gravitational waveswith pulsars.Finally, the methodology we have described can be easily ex-tended to account for shape variations associated with geodetic pre- cession in relativistic systems, such as the Hulse-Taylor Binary andthe double pulsar system, in which pulsars show dramatic secularshape variability. While these shape variations are not stochastic,by accounting for shape variability simultaneously to timing, a ro-bust characterization of the relativistic orbit can be obtained, andthe precision of tests of general relativity will be increased.
10 ACKNOWLEDGEMENTS
The Parkes radio telescope is part of the Australia Telescope Na-tional Facility which is funded by the Commonwealth of Australiafor operation as a National Facility managed by CommonwealthScience and Industrial Research Organization (CSIRO). We thankall of the observers, engineers, and Parkes observatory sta ff mem-bers who have assisted with the observations reported in this paper.LL was supported by a Junior Research Fellowship at Trin-ity Hall College, Cambridge University. RMS acknowledges travelsupport from the CSIRO through a John Philip award for excellencein early career research. REFERENCES
Coles W., Hobbs G., Champion D. J., Manchester R. N., VerbiestJ. P. W., 2011, Mon. Not. R. Astron. Soc., 418, 561Dai S., et al., 2015, Mon. Not. R. Astron. Soc., 449, 3223Demorest P. B., et al., 2013, The Astrophys. J., 762, 94Detweiler S., 1979, The Astrophys. J., 234, 1100Handley W. J., Hobson M. P., Lasenby A. N., 2015, ArXiv e-printsHankins T. H., Cordes J. M., 1981, The Astrophys. J., 249, 241Hobbs G., Lyne A. G., Kramer M., 2010, Mon. Not. R. Astron.Soc., 402, 1027Hotan A. W., Bailes M., Ord S. M., 2004, Mon. Not. R. Astron.Soc., 355, 941Hotan A. W., Bailes M., Ord S. M., 2005, Mon. Not. R. Astron.Soc., 362, 1267Hotan A. W., Bailes M., Ord S. M., 2006, Mon. Not. R. Astron.Soc., 369, 1502Jenet F. A., Anderson S. B., 1998, PASP, 110, 1467Kass R. E., Raftery A. E., 1995, Journal of the American Statisti-cal Association, 90, 791Kelly B. C., McKay T. A., 2004, The Astrophys. J., 127, 625Lentati L., Alexander P., Hobson M. P., 2015, Mon. Not. R. As-tron. Soc., 447, 2159Lentati L., Alexander P., Hobson M. P., Feroz F., van HaasterenR., Lee K. J., Shannon R. M., 2014, Mon. Not. R. Astron. Soc.,437, 3004Lentati L., Carilli C., Alexander P., Maiolino R., Wang R., Cox P.,Downes D., McMahon R., Menten K. M., Neri R., Riechers D.,Wagg J., Walter F., Wolfe A., 2013, Mon. Not. R. Astron. Soc.,430, 2454Lentati L., Hobson M. P., Alexander P., 2014, Mon. Not. R. As-tron. Soc., 444, 3863Liu K., Desvignes G., Cognard I., Stappers B. W., VerbiestJ. P. W., Lee K. J., Champion D. J., Kramer M., Freire P. C. C.,Karuppusamy R., 2014, Mon. Not. R. Astron. Soc., 443, 3752Liu K., Karuppusamy R., Lee K. J., Stappers B. W., Kramer M.,Smits R., Purver M. B., Janssen G. H., Perrodin D., 2015, Mon.Not. R. Astron. Soc., 449, 1158Lyne A., Hobbs G., Kramer M., Stairs I., Stappers B., 2010, Sci-ence, 329, 408 c (cid:13) , 000–000 L. Lentati et al.
Manchester R. N., et al., 2013, PASA, 30, 17Neal R. M., 2000, ArXiv Physics e-printsPennucci T. T., Demorest P. B., Ransom S. M., 2014, The Astro-phys. J., 790, 93Phinney E. S., 2001, ArXiv Astrophysics e-printsRefregier A., 2003, Mon. Not. R. Astron. Soc., 338, 35Refregier A., Bacon D., 2003, Mon. Not. R. Astron. Soc., 338, 48Rickett B J., 1969, 221, 158Romani R. W., Johnston S., 2001, The Astrophys. J., 557, L93Sallmen S., Backer D. C., Hankins T. H., Mo ff ett D., Lundgren S.,1999, The Astrophys. J., 517, 460Sazhin M. V., 1978, Soviet Astronomy, 22, 36Shannon R. M., et al., 2014, Mon. Not. R. Astron. Soc., 443, 1463Shannon R. M., Ravi V., Lentati L. T., Lasky P. D., Hobbs G., KerrM., Manchester R. N., Coles W. A., Levin Y., Bailes M., Bhat N.D. R., et. al. 2015, Science, 349, 1522Shao L., Caballero R. N., Kramer M., Wex N., Champion D. J.,Jessner A., 2013, Classical and Quantum Gravity, 30, 165019Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds,American Institute of Physics Conference Series Vol. 735 ofAmerican Institute of Physics Conference Series, Nested Sam-pling. pp 395–405Taylor J. H., 1992, Royal Society of London Philosophical Trans-actions Series A, 341, 117van Haasteren R., et al., 2011, Mon. Not. R. Astron. Soc., 414,3117van Haasteren R., Levin Y., 2013, Mon. Not. R. Astron. Soc., 428,1147Wolszczan A., Cordes J., Stinebring D., 1984, in Reynolds S. P.,Stinebring D. R., eds, Birth and Evolution of Neutron Stars: Is-sues Raised by Millisecond Pulsars a Single Pulse Study of theMillisecond Pulsar 1937 +
214 (Poster). p. 63 c (cid:13)000