Applying clock comparison methods to pulsar timing observations
MMNRAS , 1–11 (2020) Preprint 9 November 2020 Compiled using MNRAS L A TEX style file v3.0
Applying clock comparison methods to pulsar timing observations
Siyuan Chen , , (cid:63) , Franc¸ois Vernotte , and Enrico Rubiola , , Station de Radioastronomie de Nanc¸ay, Observatoire de Paris, PSL University, CNRS, Universit´e d’Orl´eans, 18330 Nanc¸ay, France FEMTO-ST, Department of Time and Frequency, UBFC and CNRS, 25030 Besanc¸on, France Laboratoire de Physique et Chimie de l’Environnement et de l’Espace, LPC2E UMR7328, Universit´e d’Orl´eans, CNRS, 45071 Orl´eans, France Observatoire des Sciences de l’Univers THETA, UBFC and CNRS, 20510 Besanc¸on, France Istituto Nazionale di Ricerca Metrologica INRiM, Divsion of Quantum Metrology and Nanotechnology, 10135 Torino, Italy
Accepted . . . Received . . . ; in original form . . .
ABSTRACT
Frequency metrology outperforms any other branch of metrology in accuracy (parts in10 − ) and small fluctuations ( < − ). In turn, among celestial bodies, the rotation speedof millisecond pulsars (MSP) is by far the most stable ( < − ). Therefore, the precisemeasurement of the time of arrival (TOA) of pulsar signals is expected to disclose informationabout cosmological phenomena, and to enlarge our astrophysical knowledge. Related tothis topic, Pulsar Timing Array (PTA) projects have been developed and operated for thelast decades. The TOAs from a pulsar can be affected by local emission and environmentaleffects, in the direction of the propagation through the interstellar medium or universallyby gravitational waves from super massive black hole binaries. These effects (signals) canmanifest as a low-frequency fluctuation over time, phenomenologically similar to a red noise.While the remaining pulsar intrinsic and instrumental background (noise) are white. Thisarticle focuses on the frequency metrology of pulsars. From our standpoint, the pulsar is anaccurate clock, to be measured simultaneously with several telescopes in order to reject theuncorrelated white noise. We apply the modern statistical methods of time-and-frequencymetrology to simulated pulsar data, and we show the detection limit of the correlated red noisesignal between telescopes. Key words: pulsars: general – methods: data analysis
Millisecond pulsars (MSP) are considered extremely stable astro-nomical clocks because of their high energy and momentum in asmall size (Verbiest et al. 2009), albeit the observations can be af-fected by large observational white noise. Such noise can be dueto the low signal to noise ratio of the MSP observations, and de-pends on the radio-telescope (RT). On the other hand, red noisecould originate from the MSP, the propagation through the interstel-lar medium or from gravitational waves (GW) on the line of sight,hence it is common to all the RTs (Hellings & Downs 1983; Jenet& Romano 2015). The characterization of the spectral signature ofGWs is a major scientific challenge because it is expected to dis-close information about the astrophysical sources. As pulsars havebeen precisely timed for decades with up to weekly cadence, thefrequency range that can be probed by Pulsar Timing Array (PTA)(Foster & Backer 1990) projects is f ∈ [ − , − ] Hz, which isa window inaccessible with the LIGO/VIRGO interferometers andat the edge of the LISA band. The most interesting and likely to bedetected source is a cosmic population of super massive black hole (cid:63)
E-mail: [email protected] binaries (Taylor et al. (2016) and Perrodin & Sesana (2018) for arecent review), whose interaction with the MSP signals introducesa correlated red noise in the Time of arrival (TOA) series, with aphase power spectral density (PSD) proportional to 1 / f / (fora circular population, see eg. Phinney (2001); Chen et al. (2017)).The choice of the algorithm is therefore of paramount importanceto disentangle the pulsar specific small red noise (signal) from theeffects of GWs in the presence of large white noise (background) inthe shortest observation time. In this regard, PTAs (Desvignes et al.2016; Alam et al. 2020; Kerr et al. 2020; Perera et al. 2019) enablethe simultaneous measurement of the same pulsar with several RTs,especially the Large European Array for Pulsars (LEAP) project(Bassa et al. 2016).The extraction of a small signal from noise exploiting simulta-neous measurements is a well-known method for the measurementof oscillators and atomic clocks. Specific tools have been devel-oped over more than 50 years using the PSD and wavelet variances(Barnes et al. 1971; Rutman 1978; Gray & Allan 1974; Allan &Barnes 1981). Among them, the cross-spectrum method (Rubiola& Vernotte 2010) deserves separate mention in view of our applica-tion. Under certain assumptions, this method rejects the backgroundnoise (uncorrelated), and converges to the oscillator noise even if © 2020 The Authors a r X i v : . [ phy s i c s . d a t a - a n ] N ov Chen et al. it is significantly lower than the background. Wavelet variances arecommonly used in the time domain, the most known of which isthe 2-sample Allan variance (AVAR), after the pioneering work ofAllan and Barnes (Allan 1966). Such variances enable to distinguishdifferent types of noise defined by their exponent in the PSD, e.g. f for white phase, 1 / f for flicker phase modulation, 1 / f for whitefrequency modulation, 1 / f for flicker frequency modulation, 1 / f for random walk frequency modulation noise, etc. The parabolicvariance (PVAR) (Benkler et al. 2015; Vernotte et al. 2016) extendsthis concept, and exhibits (i) the highest rejection of white noise, and(ii) the efficient detection red noise underneath background with theshortest data record. Moreover, the wavelet covariance (Fest et al.1983; Lantz et al. 2019), i.e., Allan covariance (ACOV) or paraboliccovariance (PCOV), enables the rejection of the background usingtwo (or more) uncorrelated instruments.This article is intended to port the time-and-frequency metrol-ogy methods to pulsar astronomy. In this respect, the pulsar is theclock under test, observed simultaneously with two or more RTsplaying the role of phase meters. We compare the results of differentmethods applied to simulated time series representing the TOAsof millisecond pulsars. In section 1 we summarize the concept ofpulsar timing as well as the spectral and variance methods from thetime-and-frequency metrology. Section 2 describes the simulationprocess of the TOA time series, followed by the statistical methodsin section 3. Results and the conclusions are presented in sections 4and 5 respectively. Pulsars are spinning neutron stars that emit radio waves from theregion above their magnetic poles. If these are misaligned with therotation axis and happen to point towards Earth as they sweep space,we receive one radio pulse each rotation. The most stable pulsarshave rotation period of milliseconds and are stable on a very longtimescale. Through the decades that they have been timed very smallvariations have been detected, which translate to low frequency rednoise in the Fourier domain. A detailed review about pulsars can befound in Lorimer & Kramer (2012).In general the time of arrival (TOA) of the pulses can be de-scribed by a precise timing model (Hobbs et al. 2006): t obs = t model + t red + t white (1)where t obs is the observed and t model is the predicted TOA consid-ering all known pulsar properties and propagation effects. For thisstudy we also include the effects of dispersion measure variationover time into the perfect model. Ideally, we should be able to pre-cisely predict the TOAs, if all model parameters are perfectly knownand there are no unknown sources of noise left. However, this is notthe case in practise, the difference between the observed and modelobservation forms the residual series x = t red + t white . Therefore, weare interested in both white and red noise components.The pulse residual series can be transformed into Fourier fre-quency space. The noise components are then described by a PSD. The PTA noise description can be found in Lentati et al. (2016);Caballero et al. (2016); Arzoumanian et al. (2016). Following thenotation of the latter, the white noise is described by two parameters E k and Q k : S white = ( E k W ) + Q k (2)where W is the initial estimate of the white coming from the ra-diometer noise and template matching error, E k and Q k are the f [Hz]10 R M S [ s ] PSDARSCS posCS neg
Figure 1.
Average PSD (blue pluses), averaged residual spectrum (ARS)(green triangles) and CS (orange crosses/points) computed from one set of 5different residual series in the boundary case so-called EFAC and EQUAD parameters respectively. Note that S white is constant over the whole frequency range.The red noise is described by a power law with two parameters: S red ( f ) = A PTA f − γ PTA (3)where A PTA is the amplitude and γ PTA is the spectral index of thepower law. The total noise PSD S ( f ) becomes (Vernotte & Zala-mansky 1999) S ( f ) = S red ( f ) + S white (4) In time and frequency metrology, the cross spectrum is the mostoften used method to measure the phase noise of an oscillator. Themethod relies on the simultaneous measurement with two sepa-rate instruments, assuming that the background noise is statisticallyindependent. The same hypothesis holds for pulsar observationsfrom different radio telescopes. The telescopes are obviously inde-pendent, and sufficiently far away from each other for the tropo-spheric/ionospheric noise to be statistically independent. We summa-rize the main results from our tutorial on the cross spectrum method(Rubiola & Vernotte 2010).Given a residual series x from an instrument, we can write itsFourier transform as X = ℜ ( X ) + i ℑ ( X ) . The one-sided PSD canbe calculated as S x ( f ) = T a XX ∗ = T a (cid:16) | ℜ ( X ) | + | ℑ ( X ) | (cid:17) for f > , (5)where T a is the acquisition time, and the factor 2 is due to energyconservation after surpressing the negative frequencies.Measuring a clock (pulsar) signal c simultaneously with twoinstruments (telescopes), we get two residual series x and y , eachcontaining the white background noise of the measurement addedto c . The cross spectrum (CS) S yx = T a Y X ∗ (6)converges to S c = T a CC ∗ . Notice that, after compensating for thedifferential path, c goes only into ℜ ( S yx ) , while the ℑ ( S yx ) containsonly the background noise of the observation. Thus the fastest and MNRAS , 1–11 (2020) omparing pulsars like clocks [days]10 l o g VA R AVARPVARACOVPCOV
Figure 2.
Averaged data Allan variance (blue) and parabolic variance (or-ange), the corresponding covariances (light dotted lines) and the associateduncertainties computed from one set of 5 residual series in the boundary case unbiased estimator following Rubiola & Vernotte (2010) isˆ S c = T a ℜ ( S yx ) = T a (cid:16) ℜ ( X ) ℜ ( Y ) + ℑ ( X ) ℑ ( Y ) (cid:17) (7)and can give negative values. We average over all CS coming fromthe different RT pairs to produce one simultaneous observationseries for the analysis.Figure 1 shows an example comparison between the differentspectra. The PSD can be computed on the averaged residual series(ARS) or on each residual series individually and then averaged(PSD). As the ARS is very similar, but slightly worse than the CS,we opt to focus the analysis on the CS method, whose results shouldbe comparable to those from the ARS. The CS is better at rejectinguncorrelated white noise and produces a lower level than the averagePSD, which serves as a good reference.We also investigated the effects of the choice of the windowfunction on the Fourier transform by comparing the sinc functionused above against a Blackman window. Although there is spectralleakage of the red noise, the amount is only noticable in the low-est frequency and is within the uncertainty of the PSD/CS value.We found minimal effects on our analysis. In fact as the spectralvalue decreases, the recovered red noise parameters provide a worsematch to the injected values. We also found similar results with aprewhitening postdarkening procedure. Therefore we will only showresults using the equations above. Different variances have been proposed to describe the amount ofvariation in a given residual series x with a total time span of T fortime steps of τ . Two prominent ones are the Allan (Allan 1966) andmore recently the parabolic variance (Vernotte et al. 2016), whichwe summarize here.The AVAR is a simple computation comparing three observa-tions separated by m time steps from each other.AVAR ( τ ) = M τ M − ∑ i = (cid:104) ( − x i + x i + m − x i + m ) (cid:105) (8)where M = N − m + N in the residual series x . It has the advantage of quite well estimatingthe largest time step τ = T / ( τ ) = Mm τ M − ∑ i = (cid:104) m − ∑ k = ( m − − k )( x i + k − x i + m + k ) (cid:105) (9)where M = N − m + x and y to compute the ACOVACOV ( τ ) = M τ M − ∑ i = (cid:104) ( − x i + x i + m − x i + m )( − y i + y i + m − y i + m ) (cid:105) (10)and PCOVPCOV ( τ ) = Mm τ M − ∑ i = (cid:104) m − ∑ k = ( m − − k )( x i + k − x i + m + k )( m − − k )( y i + k − y i + m + k ) (cid:105) (11)A comparison of the two different (co)variances can be foundin figure 2. In general, variance plots have two prominent features:(i) at lower time steps the variance is dominated by white noise, thisleads to a power law decrease with fixed spectral index; (ii) at highertime steps red noise start to appear, this shows up as a flattening andthen rising of the variance depending on the amplitude and spectralindex of the red noise. The variance curve can be parametrized as VAR = A VAR τ γ VAR + B VAR τ β VAR (12)The relation between the PTA parameters EFAC, EQUAD, γ PTA and A PTA and the variance parameters is not trivial. As the whitenoise follows a fixed theoretical spectral index for a given variance β VAR = − / − S white fromthe PTA white noise in (2): S white = (cid:40) B VAR / B VAR / ( τ ) for PVAR (13)where τ is the sampling time.The two red noise spectral indices can be related via γ PTA = γ VAR + . (14)The red noise amplitudes also follow analytic relations for giveninteger spectral indices (see table 1 in Vernotte et al. 2016). Underthe condition that γ PTA ∈ [ . , . ] it is also possible to find ananalytic relation between A VAR and A PTA for non-integer γ PTA , seealso Vernotte et al. (2020): A VAR = A PTA π π yr + γ PTA A ( γ PTA ) (15)where A ( γ PTA ) is a function that relates to the coefficients of thevariance. Using the same approach as Walter (1994) we find thefolowing relationships for AVAR and PVAR: A ( γ PTA ) = ( γ PTA − − γ PTA − ) π γ PTA − Γ ( − γ PTA ) sin ( π − πγ PTA / ) for AVAR − γ PTA π γ PTA − ( − γ PTA ( − − γ PTA ) + + γ PTA − γ PTA ) Γ ( − − γ PTA ) sin ( π − πγ PTA / ) for PVAR(16) MNRAS000
Averaged data Allan variance (blue) and parabolic variance (or-ange), the corresponding covariances (light dotted lines) and the associateduncertainties computed from one set of 5 residual series in the boundary case unbiased estimator following Rubiola & Vernotte (2010) isˆ S c = T a ℜ ( S yx ) = T a (cid:16) ℜ ( X ) ℜ ( Y ) + ℑ ( X ) ℑ ( Y ) (cid:17) (7)and can give negative values. We average over all CS coming fromthe different RT pairs to produce one simultaneous observationseries for the analysis.Figure 1 shows an example comparison between the differentspectra. The PSD can be computed on the averaged residual series(ARS) or on each residual series individually and then averaged(PSD). As the ARS is very similar, but slightly worse than the CS,we opt to focus the analysis on the CS method, whose results shouldbe comparable to those from the ARS. The CS is better at rejectinguncorrelated white noise and produces a lower level than the averagePSD, which serves as a good reference.We also investigated the effects of the choice of the windowfunction on the Fourier transform by comparing the sinc functionused above against a Blackman window. Although there is spectralleakage of the red noise, the amount is only noticable in the low-est frequency and is within the uncertainty of the PSD/CS value.We found minimal effects on our analysis. In fact as the spectralvalue decreases, the recovered red noise parameters provide a worsematch to the injected values. We also found similar results with aprewhitening postdarkening procedure. Therefore we will only showresults using the equations above. Different variances have been proposed to describe the amount ofvariation in a given residual series x with a total time span of T fortime steps of τ . Two prominent ones are the Allan (Allan 1966) andmore recently the parabolic variance (Vernotte et al. 2016), whichwe summarize here.The AVAR is a simple computation comparing three observa-tions separated by m time steps from each other.AVAR ( τ ) = M τ M − ∑ i = (cid:104) ( − x i + x i + m − x i + m ) (cid:105) (8)where M = N − m + N in the residual series x . It has the advantage of quite well estimatingthe largest time step τ = T / ( τ ) = Mm τ M − ∑ i = (cid:104) m − ∑ k = ( m − − k )( x i + k − x i + m + k ) (cid:105) (9)where M = N − m + x and y to compute the ACOVACOV ( τ ) = M τ M − ∑ i = (cid:104) ( − x i + x i + m − x i + m )( − y i + y i + m − y i + m ) (cid:105) (10)and PCOVPCOV ( τ ) = Mm τ M − ∑ i = (cid:104) m − ∑ k = ( m − − k )( x i + k − x i + m + k )( m − − k )( y i + k − y i + m + k ) (cid:105) (11)A comparison of the two different (co)variances can be foundin figure 2. In general, variance plots have two prominent features:(i) at lower time steps the variance is dominated by white noise, thisleads to a power law decrease with fixed spectral index; (ii) at highertime steps red noise start to appear, this shows up as a flattening andthen rising of the variance depending on the amplitude and spectralindex of the red noise. The variance curve can be parametrized as VAR = A VAR τ γ VAR + B VAR τ β VAR (12)The relation between the PTA parameters EFAC, EQUAD, γ PTA and A PTA and the variance parameters is not trivial. As the whitenoise follows a fixed theoretical spectral index for a given variance β VAR = − / − S white fromthe PTA white noise in (2): S white = (cid:40) B VAR / B VAR / ( τ ) for PVAR (13)where τ is the sampling time.The two red noise spectral indices can be related via γ PTA = γ VAR + . (14)The red noise amplitudes also follow analytic relations for giveninteger spectral indices (see table 1 in Vernotte et al. 2016). Underthe condition that γ PTA ∈ [ . , . ] it is also possible to find ananalytic relation between A VAR and A PTA for non-integer γ PTA , seealso Vernotte et al. (2020): A VAR = A PTA π π yr + γ PTA A ( γ PTA ) (15)where A ( γ PTA ) is a function that relates to the coefficients of thevariance. Using the same approach as Walter (1994) we find thefolowing relationships for AVAR and PVAR: A ( γ PTA ) = ( γ PTA − − γ PTA − ) π γ PTA − Γ ( − γ PTA ) sin ( π − πγ PTA / ) for AVAR − γ PTA π γ PTA − ( − γ PTA ( − − γ PTA ) + + γ PTA − γ PTA ) Γ ( − − γ PTA ) sin ( π − πγ PTA / ) for PVAR(16) MNRAS000 , 1–11 (2020)
Chen et al. R e s i d u a l s [ s ] MJD
Figure 3.
Example set of 5 different simulated residual series of a 4000-daylong observing campaign with a monthly cadence in the boundary case
Clock comparison methods require several time series of the samepulsar observed simultaneously with different telescopes over a longtime span. The LEAP project uses the 5 major radio telescopes inEurope: Effelsberg in Germany, Lovell in UK, Westerbork in theNetherlands, Nanc¸ay in France and Sardinia in Italy as an interfer-ometer with monthly simultaneous observations since 2012 (Bassaet al. 2016; Smits et al. 2017).Thus we simulate 100 realizations of 5 TOA series over a 4000-day long observing campaign with a cadence of 30 days, this is ofthe order of the nominal LEAP observations up to today. These havebeen created using the libstempo package (Vallisneri 2015). Forsimplicity, we assume that all 5 RTs have comparable instrumentalwhite noise, and give all residuals the initial uncertainty of W = = . , EQUAD = S / N ∝ / W . This also translates into an EFAC injection of1 . / = .
21, while EQUAD remains at 0.In total we consider three cases:(i) white noise case: no red noise signal injected(ii) boundary case: (log A PTA = − , γ PTA = A PTA = − , γ PTA = The standard PTA analysis is done via the enterprise package(Ellis et al. 2017), details of the analysis can be found in Arzouma-nian et al. (2016, 2018). We use all 5 individual RT residual seriesfor the PTA analysis comparison. For the LSL analysis we only haveone residual series, but with a lower white noise level.We compute the PSD and CS through equations (5) and (7),see figure 1, and compare them directly to equation (4) to get con-straints on the noise parameters. The PSD can be averaged over the5 residual series, whilst the CS can be computed from the averageof the 10 different combinations of 2 out of the 5 residual series.The Allan and parabolic Variances are calculated through equa-tions (8) and (9), see figure 2, and matched to the parametrizationin equation (12) with the use of equations (13), (14) and (15). Un-like the PSD we use the averaged residual series to compute thevariances.We use Bayesian statistics and a nested sampling algorithm,implemented in the cpnest package (Del Pozzo & Veitch 2015), tolook for the posterior constraints on the white noise properties aswell as the spectral index and amplitude of the red noise parameters.A likelihood function can be derived from a weighted leastsquares fitting procedure: L = − σ i (cid:16) D − F ( φ ) (cid:17) , (17)where σ i are the weights, D are the computed PSD,CS, AVAR or PVAR data points and F ( φ ) are the cor-responding parametrizations with the model parameters φ =( EFAC , log EQUAD , γ PTA , log A PTA ) .For the spectral analyis, we note that the spread of the 100 real-izations give roughly equally broad distributions at each frequency,see figures 7 and 9 left panels. Thus, we model each frequency withas a logarithmic Gaussian centered around the computed value witha uncertainty of 0.3, ie. the weights are σ i = . ∗ M i , where M i isthe average value of the 100 realizations at a given frequency.For the (co)variance analysis, we compute the weights fromthe effective degrees of freedom for each time step. A more detaileddescription can be found in the appendix A.The choice of priors follows closely the standard PTA analysis(eg. Lentati et al. (2015)) wherever possible with uniform distribu-tions in: EFAC ∈ [ . , . ] , log EQUAD ∈ [ − , − ] , γ PTA ∈ [ , ] and log A PTA ∈ [ − , − ] . The only exception are the boundariesof γ PTA ∈ [ . , . ] for the variance analysis. This is due to the lim-ited range of validity of equation (16) (Vernotte et al. 2020). To ensure that our results are statistically robust, each of the threecases is based on 100 different realization sets and doing the analysisusing all the clock comparison methods and the standard PTA/LSLcomparisons. All following posterior distributions and recoveriesare based on the addition of all 100 individual posterior distributions,unless stated otherwise. They should thus be representative of theconstraints from the analyses of a single simulated set of 5 residualseries.
We first explore the performance of the clock methods in the case ofpure white noise in the residual series, and compare it to the tradi-tional PTA/LSL analysis. Figure 4 shows the PTA noise parameter
MNRAS , 1–11 (2020) omparing pulsars like clocks marginalized 1D and 2D posterior distributions in a corner plot. Asthere is no red noise present, A PTA and γ PTA return very uninfor-mative posteriors for all methods. γ PTA simply returns a nearly flatuniform prior. The quantity A PTA returns a posterior distribution thatis flat for low amplitudes, but has a cutoff at an higher amplitude.This can be interpreted as a upper limit case, as red noise with ahigher amplitude should have affected the residual series. Similarly,as there is no EQUAD injected, all analyses also show upper limits.The only parameter of interest in the case of white noise isEFAC. As expected the PSD and PTA methods both return distri-butions around the injected value, with the PTA distribution beingmore constraining than the PSD. The CS method shows a cleartrend to the lower range of white noise, see top panel of figure 4.This shows that there is no correlated white noise in the residualseries. In the middle panel of figure 4, both AVAR and PVAR showsimilar, but worse, contraints than the PTA analysis, where AVARis closer to the injected value and PVAR has a small bias towardslower values. The bottom panel shows that the covariances rejectmost of the white noise, pushing EFAC towards to lower end of theallowed values. This effect is intended as the covariances ’try’ tofind a correlated signal, but as there is only uncorrelated noise aupper limit on EFAC is placed.However, we also track the evolution of the performance ofthe different methods to recover white noise as the amount of rednoise gradually increases as we move through the 3 simulationcases. The main result in figure 5 is that as the red noise increases,all methods seem to progressively perform worse in constrainingthe white noise. For the spectral and variance analysis the posteriordistributions progressively widen. Particularly in the strong red noisecase the AVAR and PVAR analysis are completely unable to recoverthe injected white noise (bottom middle panel). For the covarianceanalysis the EFAC upper limit gradually increases, most notable inthe strong red noise case.A different visual aid is given in figure 6, which shows theconstraints on EFAC from the individual analyses of a subset of 10sets of realizations.
A table summarizing the constraints on the PTA parameters from thestandard PTA/LSL analysis and the clock comparison methods isgiven in table 1. The posterior distributions for all analysis methodswith the recovered spectra and variances can be found in figure 7.Examples from analyses of a single simulation set are shown infigure 8. Despite no injection of a EQUAD value, we still performthe all analyses with all 4 PTA noise parameters. However, therecovered EQUAD posterior distributions are all upper limits asin the white noise case above. Thus, we will omit the EQUADparameter from our discussion from here on, as the focus will be onthe red noise parameters.
The comparison of the marginalized 1D and 2D posterior distribu-tions of the PTA noise parameter for the boundary case from thespectral analysis can found in the top left corner plot in figure 7. Themiddle and bottom left panels in figure 7 show the injected red andwhite noise along with the recovered PSD and CS respectively.The red noise signal is not detectable for the PSD as only thelowest frequency bin in the left middle panel shows signs of deviat-ing away from white noise. The recovered spectral space is strongly l o g E Q U A D . . . . P T A . . . . . EFAC . . . . . l o g A P T A log EQUAD . . . . PTA . . . . . log A PTA
PSDPSDCSPSDCSPTAPSDCSPTALSL l o g E Q U A D . . . . P T A . . . . . EFAC . . . . . l o g A P T A log EQUAD . . . . PTA . . . . . log A PTA
AVARAVARPVARAVARPVARPTAAVARPVARPTALSLAVARPVARPTALSL l o g E Q U A D . . . . P T A . . . . . EFAC . . . . . l o g A P T A log EQUAD . . . . PTA . . . . . log A PTA
ACOVACOVPCOVACOVPCOVPTAACOVPCOVPTALSLACOVPCOVPTALSL
Figure 4.
Comparison of the posterior distribution corner plots for the whitenoise case: spectral (top), variance (middle) and covariance analyses (bottom)MNRAS000
Comparison of the posterior distribution corner plots for the whitenoise case: spectral (top), variance (middle) and covariance analyses (bottom)MNRAS000 , 1–11 (2020)
Chen et al.
EFAC
Injection: 1.05PSD: 1.05CS: 0.11PTA: 1.06LSL: 0.23
EFAC
Injection: 1.05AVAR: 1.00PVAR: 0.90PTA: 1.06LSL: 0.23
EFAC
Injection: 1.05ACOV: 0.18PCOV: 0.16PTA: 1.06LSL: 0.23
EFAC
Injection: 1.05PSD: 1.06CS: 0.11PTA: 1.06LSL: 0.22
EFAC
Injection: 1.05AVAR: 1.05PVAR: 0.91PTA: 1.06LSL: 0.22
EFAC
Injection: 1.05ACOV: 0.25PCOV: 0.21PTA: 1.06LSL: 0.22
EFAC
Injection: 1.05PSD: 1.06CS: 0.13PTA: 1.07LSL: 0.30
EFAC
Injection: 1.05AVAR: 0.90PVAR: 0.67PTA: 1.07LSL: 0.30
EFAC
Injection: 1.05ACOV: 0.53PCOV: 0.30PTA: 1.07LSL: 0.30
Figure 5.
Evolution of the recovery of the EFAC parameter from the white noise case (top row), boundary case (middle row) to the red noise case (bottomrow): spectral (left column), variance (middle column) and covariance methods (right coloumn). The numbers in the legends represent the median values of theposterior distributions.parameter EFAC γ PTA log A PTA injection 1 .
05 3 − . + . − . . + . − . − . + . − . LSL 0 . + . − . . + . − . − . + . − . PSD 1 . + . − . . + . − . − . + . − . CS 0 . + . − . . + . − . − . + . − . AVAR 1 . + . − . . + . − . − . + . − . PVAR 0 . + . − . . + . − . − . + . − . ACOV 0 . + . − . . + . − . − . + . − . PCOV 0 . + . − . . + . − . − . + . − . Table 1.
List of the average median posterior values and average 90% cen-tral region bounds for the 3 constrainable PTA noise parameters from theanalyses of 100 realizations in the boundary case with different analysismethods. dominated by the prior choice, ie. the median recovery is close towhite noise only. Thus, it is insufficient to put any constraints on the spectral index γ PTA and amplitude A PTA from the PSD method inthe top left panel.The decrease of the white noise level through the CS allowsthe lowest few frequency bins to contribute to the recovery of thered noise power law, see bottom left panel of figure 7. However, thelarge uncertainties hinder a clear detection and leaves the possibilityfor a long tail. The constraints on the two PTA red noise parametersare therefore slightly worse than the standard PTA analysis, whichshows a clearer hint of the injected values. In the LSL the signal isclearly detectable and the tightest constraints can be placed in thetop left panel.
The middle column of figure 7 shows the comparison red noiseparameter corner plot in the top panel and recovered Allan (middlepanel) and parabolic (bottom panel) variances with the injected rednoise signal. As equation (16) is only valid for 1 . ≤ γ PTA ≤ . γ PTA and amplitude A PTA can still be compared, however the cutoffs at 1.5 and 4.5 limitthe comparison and can thus influence the red noise recovery.The recovered variances in the middle and bottom panels of
MNRAS , 1–11 (2020) omparing pulsars like clocks
40 42 44 46 48 realization E F A C Injection: 1.05PSD: 1.05CS: 0.13 PTA: 1.06LSL: 0.23 40 42 44 46 48 realization E F A C Injection: 1.05AVAR: 1.02PVAR: 0.91 PTA: 1.06LSL: 0.23 40 42 44 46 48 realization E F A C Injection: 1.05ACOV: 0.22PCOV: 0.19 PTA: 1.06LSL: 0.2340 42 44 46 48 realization E F A C Injection: 1.05PSD: 1.06CS: 0.13 PTA: 1.06LSL: 0.23 40 42 44 46 48 realization E F A C Injection: 1.05AVAR: 1.06PVAR: 0.92 PTA: 1.06LSL: 0.23 40 42 44 46 48 realization E F A C Injection: 1.05ACOV: 0.29PCOV: 0.25 PTA: 1.06LSL: 0.2340 42 44 46 48 realization E F A C Injection: 1.05PSD: 1.06CS: 0.15 PTA: 1.07LSL: 0.30 40 42 44 46 48 realization E F A C Injection: 1.05AVAR: 0.91PVAR: 0.68 PTA: 1.07LSL: 0.30 40 42 44 46 48 realization E F A C Injection: 1.05ACOV: 0.57PCOV: 0.32 PTA: 1.07LSL: 0.30
Figure 6.
Examples of the recovery of the EFAC parameter from one set of realizations of the white noise case (top row), boundary case (middle row) and thered noise case (bottom row): spectral (left column), variance (middle column) and covariance methods (right coloumn). The numbers in the legends representthe average median values of the posterior distributions. the middle column of figure 7 illustrate the fact that the highestoctave of τ with its very large uncertainty is typically unreliableand contributes very little to the overall analysis. Since we arelooking at a boundary case, the injections are close or within theuncertainties of the variance computation. Therefore they do notnecessarily coincide with the median of the recovered variancevalues. Additionally, as with the spectral analysis above the priorchoice allows for a large distribution of small amplitude red noise,producing a median variance that is close to white noise only.For AVAR this results in it not detecting any red noise in theboundary case. However, as there a small turn-up in the highertime steps, PVAR is able to pick up some evidence of the injectedred noise. This can be seen in the top row middle panel cornerplot, especially in the A PTA parameter with PVAR and a smallimprovement in constraining the spectral index. The PTA analysis isalso able to pick up the red noise and is slightly better in performancethan PVAR. With the lowest white noise level, it is no surprise thatthe LSL analysis performs the best.
We can further focus on the red noise detection by sacrificing theconstraints on the white noise. This is done in a similar fashion aswith the CS by computing the covariances, such that uncorrelatedwhite noise will be rejected to allow for the correlated red noise to be more prominent. The right column of figure 7 shows the posteriordistributions, recovered ACOV and PCOV from top to bottom. Thesame limitation on the spectral index γ PTA for variances applies alsoto the covariances.Despite the lower white noise from the covariances the con-straints on the red noise parameters γ PTA and A PTA are actuallyworse than their variance counterparts, as it can be seen in the topright panel of figure 7. The white noise dominated lower time steps τ become upper limits with large uncertainties down to zero. Thisleaves only one or two time steps that contain information on the rednoise, which is shown in the right column middle and bottom panels.The overall recovery using all time steps thus favours low amplitudewhite and red noise. Using only time steps with small uncertaintiescould produce tighter, but less reliable, red noise constraints. The strong red noise case is of particular interest as it shows howwell the different methods can in principle divide different red noiserealizations in the strong signal regime. An overview on the con-straints on the PTA parameters for the red noise case is given bytable 2 for the standard PTA/LSL and clock metrology analysis. Theoverall posterior distributions for all methods are shown in figure 9.The results of some example analyses of a single simulation set canbe found in figure 10.
MNRAS000
MNRAS000 , 1–11 (2020)
Chen et al. . . . . PTA . . . . . l o g A P T A . . . . . log A PTA
PSDPSDCSPSDCSPTAPSDCSPTALSL . . . . PTA . . . . . l o g A P T A . . . . . log A PTA
AVARAVARPVARAVARPVARPTAAVARPVARPTALSL . . . . PTA . . . . . l o g A P T A . . . . . log A PTA
ACOVACOVPCOVACOVPCOVPTAACOVPCOVPTALSL f [Hz]10 R M S [ s ] PSD recoveryInjectionPSDMedian [days]10 l o g VA R AVAR recoveryInjectionAVARMedian [days]10 l o g VA R ACOV recoveryinjectionACOVMedian f [Hz]10 R M S [ s ] CS recoveryInjectionCS posMedian [days]10 l o g VA R PVAR recoveryInjectionPVARMedian [days]10 l o g VA R PCOV recoveryinjectionPCOVMedian
Figure 7.
The left, middle and right column show results from the spectral, variance and covariance analysis respectively for the boundary case. The top rowshows the posterior distribution corner plots. The middle row shows the recovered PSD (left), AVAR (middle) and ACOV (right). The bottom row shows therecovered CS (left), PVAR (middle) and PCOV (right). Each recovery figure contains the values from the 100 realizations in black dots, the injection in reddotted lines and median values and uncertainties in blue crosses/bands.
40 42 44 46 48 realization P T A Injection: 3.00PSD: 3.18CS: 3.24 PTA: 3.55LSL: 3.24 40 42 44 46 48 realization P T A Injection: 3.00AVAR: 2.87PVAR: 2.74 PTA: 3.55LSL: 3.24 40 42 44 46 48 realization P T A Injection: 3.00ACOV: 2.95PCOV: 2.97 PTA: 3.55LSL: 3.24
40 42 44 46 48 realization l o g A P T A Injection: -14.00PSD: -16.67CS: -15.93 PTA: -14.48LSL: -14.16 40 42 44 46 48 realization l o g A P T A Injection: -14.00AVAR: -16.37PVAR: -15.29 PTA: -14.48LSL: -14.16 40 42 44 46 48 realization l o g A P T A Injection: -14.00ACOV: -16.83PCOV: -16.46 PTA: -14.48LSL: -14.16
Figure 8.
Examples of the recovery of the γ PTA (top row) and A PTA (bottom row) parameters from one set of realizations of the boundary case: spectral (leftcolumn), variance (middle column) and covariance methods (right coloumn). The numbers in the legends represent the average median values of the posteriordistributions. MNRAS , 1–11 (2020) omparing pulsars like clocks . . . . PTA . . . . . l o g A P T A . . . . . log A PTA
PSDPSDCSPSDCSPTAPSDCSPTALSL . . . . PTA . . . . . l o g A P T A . . . . . log A PTA
AVARAVARPVARAVARPVARPTAAVARPVARPTALSL . . . . PTA . . . . . l o g A P T A . . . . . log A PTA
ACOVACOVPCOVACOVPCOVPTAACOVPCOVPTALSL f [Hz]10 R M S [ s ] CS recoveryInjectionCS posMedian [days]10 l o g VA R PVAR recoveryInjectionPVARMedian [days]10 l o g VA R PCOV recoveryinjectionPCOVMedian
Figure 9.
Top and bottom rows correspond those in figure 7 respectively in the same style. However, the results are for the red noise case.
40 42 44 46 48 realization P T A Injection: 3.00PSD: 2.97CS: 3.12 PTA: 3.04LSL: 3.09 40 42 44 46 48 realization P T A Injection: 3.00AVAR: 2.76PVAR: 2.73 PTA: 3.04LSL: 3.09 40 42 44 46 48 realization P T A Injection: 3.00ACOV: 3.03PCOV: 3.04 PTA: 3.04LSL: 3.09
40 42 44 46 48 realization l o g A P T A Injection: -13.00PSD: -14.01CS: -13.26 PTA: -13.02LSL: -13.03 40 42 44 46 48 realization l o g A P T A Injection: -13.00AVAR: -13.10PVAR: -13.01 PTA: -13.02LSL: -13.03 40 42 44 46 48 realization l o g A P T A Injection: -13.00ACOV: -13.29PCOV: -13.12 PTA: -13.02LSL: -13.03
Figure 10.
Examples of the recovery of the γ PTA (top row) and A PTA (bottom row) parameters from one set of realizations of the red noise case: spectral(left column), variance (middle column) and covariance methods (right coloumn). The numbers in the legends represent the median values of the posteriordistributions.
In general all methods can detect the injected red noise sig-nal very well. As the red noise is very strong, there is very littledifference in the constraints on the red noise parameters γ PTA andlog A PTA between the PTA and the LSL analysis. Both are centeredaround the injected values with tight constraints, see top panels of fig-ure 9. So does the CS recovery, although with a small tail. Whereasthe PSD gives broader distributions than the PTA comparison. Theuncertainties in the spectral methods allow for small values even atthe lowest frequencies. This leads to a tail of lower amplitudes being recovered, see table 2. Both variances and covariances perform verysimilarly to the PTA/LSL analysis in terms of the tightness of theconstraints on the red noise parameters, especially PVAR and PCOVare comparable to the LSL analysis. This is due to the fact that thehigher time steps are dominated by red noise and have therefore verysimilar VAR/COV values. Additionally, the absolute uncertainties σ i are the same for both variances and covariances.Consequently, the covariances of the lowest time steps τ arenot as well defined as the variances for the same time steps. These MNRAS000
In general all methods can detect the injected red noise sig-nal very well. As the red noise is very strong, there is very littledifference in the constraints on the red noise parameters γ PTA andlog A PTA between the PTA and the LSL analysis. Both are centeredaround the injected values with tight constraints, see top panels of fig-ure 9. So does the CS recovery, although with a small tail. Whereasthe PSD gives broader distributions than the PTA comparison. Theuncertainties in the spectral methods allow for small values even atthe lowest frequencies. This leads to a tail of lower amplitudes being recovered, see table 2. Both variances and covariances perform verysimilarly to the PTA/LSL analysis in terms of the tightness of theconstraints on the red noise parameters, especially PVAR and PCOVare comparable to the LSL analysis. This is due to the fact that thehigher time steps are dominated by red noise and have therefore verysimilar VAR/COV values. Additionally, the absolute uncertainties σ i are the same for both variances and covariances.Consequently, the covariances of the lowest time steps τ arenot as well defined as the variances for the same time steps. These MNRAS000 , 1–11 (2020) Chen et al. parameter EFAC γ PTA log A PTA injection 1 .
05 3 − . + . − . . + . − . − . + . − . LSL 0 . + . − . . + . − . − . + . − . PSD 1 . + . − . . + . − . − . + . − . CS 0 . + . − . . + . − . − . + . − . AVAR 0 . + . − . . + . − . − . + . − . PVAR 0 . + . − . . + . − . − . + . − . ACOV 0 . + . − . . + . − . − . + . − . PCOV 0 . + . − . . + . − . − . + . − . Table 2.
List of the average median posterior values and average 90% cen-tral region bounds for the 3 constrainable PTA noise parameters from theanalyses of 100 realizations in the red noise case with different analysismethods. large relative uncertainties and the lower COV values translate into asmall bias for the red noise to have steeper spectral indices γ PTA thanwith the variance methods. On the other hand the variance methodsare still a little bit impacted by the choice of priors, see table 2.This can be seen in the bottom middle panel of figure 9, where themedian recovered PVAR is lower than the median PVAR values ofthe 100 realizations. The same is also true for AVAR recovery.Table 2 also shows that both AVAR and ACOV have a longertail towards lower red noise amplitudes log A PTA when comparedto PVAR and PCOV. This can be the advantage of the parabolic(co)variance transitioning from a spectral index of β VAR = − γ VAR = β VAR = − γ VAR = We have shown that clock comparison methods, such as the crossspectrum and variances, can be applied to pulsar timing observa-tions. Initial tests show that these methods produce slightly worse,sometimes comparable constraints on the pulsar noise parameterswhen compared to the standard PTA and optimal LSL analysis.Which perform very well in recovering the injected white noise,regardless of the level of red noise. The red noise recovery itselfis broadly consistent with the injections, with different realizationsgiving similar constraints. In general, the metrology clock compar-ison methods perform well in the pure white noise case, worse inthe boundary case and comparable in the strong red noise case. Thespread of recoveries between different realizations is also larger thanthe comparison PTA/LSL analysis.The PSD method produces broader parameter posterior distri-butions against the PTA analysis. The CS method reduces the whiteto be negligible, however with large uncertainties. Therefore, thebenefits of having a large number of frequency bins contribute tothe red noise recovery is counteracted by the possibility to havelow A PTA and γ PTA being allowed by the uncertainties. Therefore, itperforms about as well as the PTA analysis, but worse than the LSLanalysis. Where the observations only have a fifth of the white noiselevel of the PTA observations, allowing for the tightest constraintson the red noise.Similar constraints can be put on the red noise parameters for both the variances and covariances. The (co)variance methodsperform particularly badly in recovering the white noise in the rednoise case. By design ACOV and PCOV should reject uncorrelatedwhite noise and only leave the correlated red noise signal. Thereforeno useful information can be gained for the white noise, whilesome extra power may be pushed into the red noise. This lack ofinformation may limit their performance in precisely constrainingthe red noise and manifest as a small bias in γ PTA or a long tail oflow red noise amplitudes A PTA .We continue the study and apply the methods described hereto the observations from the LEAP project, where we deal with thechallenges from real pulsar timing data. These include unevenlysampled data points, different levels of white noise between RTs andobservations and uncertainties in the computation of the spectraldensities and (co)variances.
DATA AVAILABILITY
The simulated residual series can be obtained from the authors uponrequest.
ACKNOWLEDGEMENTS
This work is funded by the ANR Programme d’Investissementd’Avenir (PIA) under the FIRST-TF network (ANR-10-LABX-48-01) project and the Oscillator IMP project (ANR-11-EQPX-0033-OSC-IMP), and by grants from the R´egion Bourgogne FrancheComt´e intended to support the PIA.The Nanc¸ay radio Observatory is operated by the Paris Obser-vatory, associated to the French Centre National de la RechercheScientifique (CNRS). We acknowledge financial support from”Programme National de Cosmologie and Galaxies” (PNCG)and ”Programme National Hautes Energies” (PNHE) funded byCNRS/INSU-IN2P3-INP, CEA and CNES, France.We acknowledge the support of our colleagues in the LargeEuropean Array for Pulsars and European Pulsar Timing Arraycollaboration.
REFERENCES
Alam M. F., et al., 2020, arXiv e-prints, p. arXiv:2005.06490Allan D. W., 1966, Proc. IEEE, 54, 221Allan D. W., Barnes J. A., 1981, in Proc. 35 FCS. Ft. Monmouth, NJ, pp470–474Arzoumanian Z., et al., 2016, Astrophysical Journal, 821, 13Arzoumanian Z., et al., 2018, Astrophysical Journal, 859, 47Barnes J. A., et al., 1971, IEEE Trans. Instrum. Meas., 20, 105Bassa C. G., et al., 2016, MNRAS, 456, 2196Benkler E., Lisdat C., Sterr U., 2015, Metrologia, 55, 565Caballero R. N., et al., 2016, MNRAS, 457, 4421Chen S., Sesana A., Del Pozzo W., 2017, MNRAS, 470, 1738Del Pozzo W., Veitch J., 2015, CPNest: Parallel nested sampling in python,https://github.com/johnveitch/cpnest, doi:10.5281/zenodo.835874Desvignes G., et al., 2016, MNRAS, 458, 3341Ellis J. A., Vallisneri M., Taylor S. R., Baker P. T., 2017, ENTERPRISE(Enhanced Numerical Toolbox Enabling a Robust PulsaR InferenceSuitE), https://github.com/nanograv/enterpriseFest D., Groslambert J., Gagnepain J.-J., 1983, IEEE Trans. Instrum. Meas.,32, 447Foster R. S., Backer D. C., 1990, Astrophysical Journal, 361, 300Gray J. E., Allan D. W., 1974, in Proc. Int’l Freq. Control Symp.. AtlanticCity, NJ, USA, pp 243–246, doi:10.1109/FREQ.1974.200027MNRAS , 1–11 (2020) omparing pulsars like clocks Hellings R. W., Downs G. S., 1983, Astrophysical Journal, 265, L39Hobbs G. B., Edwards R. T., Manchester R. N., 2006, MNRAS, 369, 655Jenet F. A., Romano J. D., 2015, American Journal of Physics, 83, 635Kerr M., et al., 2020, Publ. Astron. Soc. Australia, 37, e020Lantz E., Calosso C. E., Rubiola E., Giordano V., Fluhr C., Dubois B.,Vernotte F., 2019, IEEE Trans. Ultras. Ferroelec. Freq. Contr., 66, 1942Lentati L., et al., 2015, MNRAS, 453, 2576Lentati L., et al., 2016, MNRAS, 458, 2161Lesage P., Audoin C., 1973, IEEE Trans. Instrum. Meas., 22, 157Lorimer D. R., Kramer M., 2012, Handbook of Pulsar Astronomy. Cam-bridge University Press, 2012Perera B. B. P., et al., 2019, MNRAS, 490, 4666Perrodin D., Sesana A., 2018, Radio Pulsars: Testing Gravity and DetectingGravitational Waves. Springer, Cham, p. 95, doi:10.1007/978-3-319-97616-7˙3Phinney E. S., 2001, arXiv e-prints, pp astro–ph/0108028Rubiola E., Vernotte F., 2010, arXiv e-prints, p. arXiv:1003.0113Rutman J., 1978, Proc. IEEE, 66, 1048Smits R., et al., 2017, Astronomy and Computing, 19, 66Taylor S. R., Vallisneri M., Ellis J. A., Mingarelli C. M. F., Lazio T. J. W.,van Haasteren R., 2016, Astrophysical Journal, 819, L6Vallisneri M., 2015, libstempo - a Python wrapper for tempo2,https://github.com/vallis/libstempoVerbiest J. P. W., et al., 2009, MNRAS, 400, 951Vernotte F., Lantz E., 2012, IEEE Transactions on Ultrasonics, Ferroelectrics,and Frequency Control, 59, 523Vernotte F., Zalamansky G., 1999, IEEE Trans. Ultras. Ferroelec. Freq.Contr., 46, 1545Vernotte F., Lenczner M., Bourgeois P.-Y., Rubiola E., 2016, IEEE Trans.Ultras. Ferroelec. Freq. Contr., 63, 611Vernotte F., Rubiola E., Chen S., 2020, arXiv e-prints, p. arXiv:2005.13631Walter T., 1994, IEEE Trans. Instrum. Meas., 43, 69
APPENDIX A: WEIGHTS AND UNBIASING OF THEVARIANCES
The uncertainty σ i on each time step i , whether variances or covari-ances, can be calculated from the average values of the variances VAR : σ i = VAR / ν (A1)where ν are the effective degrees of freedom, whose expressionsdepend on the type of noise. The expressions used in the analysisfor AVAR can be found in Lesage & Audoin (1973): for white phasemodulation ( γ PTA = ν ≈ ( N + )( N − m ) ( N − m ) (A2)and flicker frequency modulation ( γ PTA = ν ≈ N m ( N + m ) (A3)The equation for PVAR is given in Vernotte et al. (2020) as ν ≈ Am / M − B ( m / M ) , (A4)where A and B are constants which depend on the type of the mod-ulation. They are given as A =
23 for γ PTA = A ≈
28 for γ PTA =
3, whilst B ≈
12 for all spectral indices. As pointed out inVernotte et al. (2020), the formula is unreliable for the highest timesteps, especially if there is red noise. Thus we adjust them usingthe 100 realizations by setting them to the ones from the AVARformulas. We also notice a mixture of white and red noise for ourrealizations. Consequently, we approximate the degrees of freedom τ
60 120 240 480 960 1920 γ PTA ν
66 65 63 8 3 1PVAR ν
103 52 27 8 3 1AVAR unbiasing factor 1 1 1 1 .
14 1 .
45 3 . .
14 1 .
45 3 . Table A1.
Effective degrees of freedom and unbiasing factors for the vari-ances in the red noise case. τ
60 120 240 480 960 1920 γ PTA ν
66 65 63 58 3 1PVAR ν
103 52 27 14 3 1AVAR unbiasing factor 1 1 1 1 1 .
45 3 . .
45 3 . Table A2.
Effective degrees of freedom and unbiasing factors for the vari-ances in the boundary case. τ
60 120 240 480 960 1920 γ PTA ν
66 65 63 58 46 6PVAR ν
103 52 27 14 11 1AVAR unbiasing factor 1 1 1 1 1 1 . . . Table A3.
Effective degrees of freedom and unbiasing factors for the vari-ances in the white noise case. by applying the corresponding equation only to the highest 3 octavesfor the red noise case and 2 octaves for the borderline case. Thedegrees of freedom used for our analysis can be found in tablesA1,A2 and A3Finally, Vernotte & Lantz (2012) have shown that the largesttime steps are statistically biased towards lower VAR/COV values.The factors that need to be applied for any time step τ with itsdegrees of freedom ν can be found in table II in Vernotte & Lantz(2012). We also show the factors used in this paper in tables A1,A2and A3. MNRAS000