Ensemble Pulsar Time Scale

CHINESE ASTRONOMY AND ASTROPHYSICS Chinese Astronomy and Astrophysics 42 (2017) 430–441

Ensemble Pulsar Time Scale†  YIN Dong-shan1,2,3 1 2

GAO Yu-ping1,2

ZHAO Shu-hong1,2

National Time Service Centre, Chinese Academy of Sciences, Xi’an 710600

Key Laboratory of Time and Frequency Primary Standards, National Time Service Center, Chinese Academy of Sciences, Xi’an 710600 3

University of Chinese Academy of Sciences, Beijing 100049

Abstract Millisecond pulsars can generate another type of time scale that is totally independent of the atomic time scale, because the physical mechanisms of the pulsar time scale and the atomic time scale are quite different from each other. Usually the pulsar timing observations are not evenly sampled, and the internals between two data points range from several hours to more than half a month. Further more, these data sets are sparse. All this makes it difficult to generate an ensemble pulsar time scale. Hence, a new algorithm to calculate the ensemble pulsar time scale is proposed. Firstly, a cubic spline interpolation is used to densify the data set, and make the intervals between data points uniform. Then, the Vondrak filter is employed to smooth the data set, and get rid of the high-frequency noises, and finally the weighted average method is adopted to generate the ensemble pulsar time scale. The newly released NANOGRAV (North American Nanohertz Observatory for Gravitational Waves) 9-year data set is used to generate the ensemble pulsar time scale. This data set includes the 9-year observational data of 37 millisecond pulsars observed by the 100-meter Green Bank telescope and the 305-meter Arecibo telescope. It is found that the algorithm used in this paper can reduce effectively the influence caused by the noises in pulsar timing residuals, and improve the long-term stability of the ensemble pulsar time scale. Results indicate that the long-term (> 1 yr) stability of the ensemble pulsar time scale is better than 3.4 × 10−15 . Key words time—pulsars—methods: observational—data analysis †

Supported by National Natural Science Foundation (11303032, 11473029), and National Key Labora-

tory of Geographic Information Engineering (SKLGIE2013-M-1-1) Received 2015–10–26; revised version 2015–12–21  

A translation of Acta Astron. Sin. Vol. 57, No. 3, pp. 326–335, 2016 [email protected]

0275-1062/17/$-see front matter © 2017 Elsevier B.V. All rights reserved. doi:10.1016/j.chinastron.2017.08.010

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

1.

431

INTRODUCTION

At present, International Atomic Time (TAI) is a time scale calculated in quasi-real time, its accuracy is continuously raising along with the emergence of a large number of high-precision atomic clocks and primary/standard clocks, and the improvement of the algorithm of atomic time. TAI, however, must be monthly calculated because of the requirement of its real-time realization, and the calculation is restricted by many constraint conditions so as to affect its accuracy to a certain degree. So, TAI is not the optimal realization of Terrestrial Time (TT). In order to provide a better atomic time scale, the International Bureau of Weights and Measures (BIPM) has given another time scale of posterior calculation, namely, TT (BIPM). A new version is annually calculated for TT (BIPM), by using all the available data of the primary standard clocks to obtain the posterior realization of TT. The blackbody radiation correction and other corrections are added into the adopted primary standard clocks in the period of calculation, and after smoothing and interpolation of the frequency the TT (BIPM) is obtained, so that the TT (BIPM) is more accurate and more stable than TAI, and therefore can be used as the reference time of the ensemble pulsar time scale. Millisecond pulsars can generate another type of time scale that is totally independent of the atomic time scale, because the physical mechanisms of the pulsar time scale and the atomic time scale are quite different from each other. It is also one of the measures to realize an ideal time scale, and equally characterized by its stability, accuracy and measurability[1−2] . In 2012, the Time Committee of the International Astronomical Union organized a working group of pulsar time scale to coordinate the international cooperation in this research field. It implies that the international academic world begins to attach importance to the study on the ensemble pulsar time scale. In general, the pulsar timing observation is not evenly sampled, the sampling rate is far lower than that of the measurement of atomic time comparison, and the model and timing noise of a pulsar clock are different from those of atomic clocks, all this makes it very difficult to calculate an ensemble pulsar time scale, and therefore the algorithm of ensemble pulsar time scale should be stressed in the study. In this paper, the ensemble pulsar time scale is calculated using a new algorithm based on the 9-year observational data recently released by NANOGRAV, the effect of noise in the pulsar time residual is analyzed, and the stability of the ensemble pulsar time scale is estimated. 2.

PRINCIPLES OF THE ALGORITHM

2.1 Vondrak Filtering By means of the Vondrak filtering, the observational data, no matter whether they are evenly spaced or not, may be effectively smoothed without the need to know their variations. The Vondrak filter is virtually similar to a single-sideband filter, as a numerical filter used for the

432

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

data treatment, it can preserve useful signals as many as possible by selecting an adequate cut-off frequency. Let the observational data be x (ti ) (i = 1, 2, · · · , N ), the basic criterion of the Vondrak filtering is[3] : Q = F + λ2 S = min , (1) in which F =



S=

2

pi [x (ti ) − x (ti )] ,

N −3 



Δ3 x (ti )

2

,

(2) (3)

i=1

ε = 1/λ2 ,

(4)

in which xi (ti ) represents the smoothed values to be obtained; pi , the weights of the observational data; F , the degree of fitting of the Vondrak smoothing; S means the total smoothness of the smoothed curve to be obtained. In essence, the Vondrak filtering method is a compromise between the absolute smoothing and the absolute fitting of observational data. In such kind of smoothing method there exists a smoothing factor ε, which determines whether the signals and noises in the observational data can be properly separated from each other, hence the selection of the smoothing factor is critical in the data smoothing process. 2.2 Classical Weighting Algorithm Let there be n pulsar clocks in a clock set, and their indications be hi (t), i = 1, 2, · · · , n. According to the classical weighting algorithm, the time scale TA(t) calculated from them is generally: n n   TA (t) = ωi (t) = 1 , (5) ωi (t) hi (t) , i=1

i=1

in which ωi (t) is the weight of the clocki . The noise of the time scale TA(t) is the weighted sum of noises of all the clocks: δS (t) =

n 

ωi (t) δi (t) .

(6)

i=1

In order to make the noise δS (t) of the synthetical clock minimum, the weights are usually determined as follows: 1/σi2 ωi (t) =  , j = 1, 2, · · · , n , (7) n 1/σj2 j=1

in which σi2 indicates the Allan variances or standard deviations, which are calculated by using the data of the every clock relative to TA(t), and the weight of the every clock is calculated with Eq.(7).

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

433

2.3 Estimation of σz (τ ) Because the pulsar’s rotation phase, frequency, rotation slowdown rate and the related astrometric parameters have been precisely determined by the applied timing model, the cubic difference of timing residuals is the remained deviation of the lowest order in the pulsar timing residuals, hence σz (τ ) is the estimate suitable for the study of the stability of pulsar time [1] . The calculation process of σz (τ ) is as follows: (1) According the order of observational time, to make the sequences which are noted as ti , xi , σi , i = 1, 2, · · · , n, for the observational times, timing residuals and errors, respectively, and the whole time span of the selected data is T = tn − t1 . (2) To extract σz (τ ). These data are continuously divided into the subsequences with the intervals of τ according to the order of the observational time. Let t0 be an arbitrary reference time, a cubic polynomial is applied to each subsequence as follows: 2

3

X (t) = c0 + c1 (t − t0 ) + c2 (t − t0 ) + c3 (t − t0 ) ,

(8)

in which ci (i = 0, 1, 2, 3) are coefficients to be fitted. A least squares fitting is made on 2 Eq.(8) so that [(xi − X (ti ))/σi ] is minimum. σz (τ ) is defined as: τ 2  1/2 σz (τ ) = √ c23 , 2 5

(9)

in which   means the weighted average value for all the subsequences, and the weight is inversely proportional to the square of the error of c3 . 3.

ALGORITHM OF THE ENSEMBLE PULSAR TIME SCALE

Because the pulsar timing observation data are not evenly sampled, and the sampling rate is far lower than that of the measurement of atomic time comparison, all this makes it very difficult to calculate an ensemble pulsar time scale, and the algorithm of ensemble pulsar time scale must be stressed in the study. The calculation process of the new algorithm of the ensemble pulsar time scale is: (1) Using the cubic spline interpolation to densify the observational data, and to make their intervals become uniform. (2) By means of the Vondrak filtering to smooth the observational set, and to get rid of high-frequency noises. (3) Using the weighting algorithm to generate the ensemble pulsar time scale, and finally to analyze the stability of the ensemble pulsar time scale. The schematic diagram of the algorithm of the ensemble pulsar time scale is shown in Fig.1.

434

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

Raw data of all pulsars (not even sampled) Cubic spline interpolation Raw data of all pulsars (even sampled) Vondrak filter Processed data of all pulsars (even sampled) Weighted average Ensemble pulsar time scale Stability analysis Analyse and evaluate of result Fig. 1 The schematic diagram of the algorithm of the ensemble pulsar time scale

4.

EXPERIMENTAL RESULTS AND ANALYSIS

4.1 Introduction to the Experimental Data The timing observation data of 37 millisecond pulsars observed in 9 years and released recently by NANOGRAV are chosen, these data were obtained with the Green Bank Telescope (GBT) and the ARECIBO telescope. The data of PSR B1937+21 are rejected, because a precise model can not be established due to its timing noises, which are obviously higher than the white noise level and unable to be corrected. Fig.2 and Table 1 show the information in detail about the observational data of other 36 pulsars[4] . Table 1 lists their fundamental parameters, including Name, Period, Dispersion Measurement (DM), Post-fit residual, Span of data, Number of observational data (Nobs ), First epoch, and Last epoch. The post-fit timing residuals of all the pulsars are obtained from the treatment of TEMPO2, and the ephemeris in option is the DE421 ephemeris of planetary barycenters of the solar system[5−7] . 4.2 Calculated Results and Analysis In order to maintain the stability of pulsar time scale, a procedure of cubic spline interpolation is first applied to interpolating the raw observational data into the equally spaced data, which are then filtered with the Vondrak filtering method so as to remove the high-frequency noises as far as possible, finally, the ensemble pulsar time scale is calculated with the new algorithm, and an analysis of its stability is carried out.

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

435

4.2.1 Cubic spline interpolation It is apparent from Fig.2 and Table 1 that the data sets of NANOGRAV are generally rather sparse, and the datum points are unevenly distributed. To increase the datum points and to make the sampling intervals of datum points uniform, after trying various interpolation methods, including the linear interpolation, binomial interpolation and cubic spline interpolation, the result shows that the cubic spline interpolation can provide the optimal data sequence. This procedure is common-used for data treatment, the details will not be described here. What should be particulary mentioned is that the various observational data sets differ in the time span and the number of data points, the cubic spline interpolation is performed only on the observational data in their own time span, so as to keep the reliability and accuracy of the observational data. The long duration extrapolation will affect the accuracy of pulsar observational data, and degrade the stability of the generated ensemble time scale. 4.2.2 Analysis of denoise effect of the Vondrak Filtering High-frequency noises are included in the pulsar observational data because of the influences of the cosmic background radiation, artificial radio interference, interstellar medium, and the thermal noise of devices, etc. on the pulsar observation. Therefore, a denoise treatment on the pulsar observational data is carried out before the calculation of the ensemble pulsar time scale. In this paper, the high-frequency noises in the observational data are reduced by means of Vondrak filtering, taking the pulsar J1012+5307 as an example, the post-fit timing residual curves before and after the Vondrak filtering with the smoothing factor of 1 × 10−9 are shown in Fig.3. It is evident from Fig.3, Fig.4 and Fig.5 that the curve after filtering can mirror the details of the pulsar data, and the residuals of the raw observational data with respect to the filtered data fluctuate up and down around the zero value, and the residuals have a Gaussian distribution, which tallies with the nature of phase white noise or of high-frequency noise. 4.2.3 Ensemble pulsar time scale After making the cubic spline interpolation and Vondrak filtering procedure on the observational data of the 36 pulsars, a weighting algorithm is adopted for the post-treatment data, and the ensemble pulsar time scale is finally obtained. Because the observational time span differs with the particular pulsar, the number of millisecond pulsars differs with the observational time, and this is also its difference from the conventional weighting algorithm. The finally-obtained ensemble pulsar time scale is shown in Fig.6[8] .

436

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

B1855+09 B1953+29 J0023+0923 J0030+0451 J0340+4130 J0613+0200 J0645+5158 J0931-1902 J1012+5307 J1024-0719 J1455-3330 J1600-3053 J1614-2230 J1640+2224 J1643-1224 J1713+0747 J1738+0333 J1741+1351 J1744-1134 J1747-4036 J1832-0836 J1853+1303 J1903+0327 J1909-3744 J1910+1256 J1918-0642 J1923+2515 J1944+0907 J1949+3106 J2010-1323 J2017+0603 J2043+1711 J2145-0750 J2214+3000 J2302+4442 J2317+1439 Fig. 2 The post-fit residuals of 36 pulsars included in NANOGRAV, within an observation duration from 53216 (MJD) to 56599 (MJD)

437

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

Table 1 The detailed parameters of 36 millisecond pulsars used in this paper Name

Nobs

Period

DM

Post-fit

Span

/ms

/(cm−3 ·pc)

residual/μs

/yr

5.362

13.30

1.347

8.9

B1953+29

6.133

104.50

6.437

7.2

1302

53945

56591

J0023+0923

3.050

14.32

1.339

2.3

4373

55758

56599

J0030+0451

4.865

4.33

1.438

8.8

2455

53394

56599

J0340+4130

3.299

49.58

5.145

1.7

3003

55972

56586

J0613−0200

3.062

38.78

1.28

8.6

7422

53448

56586

J0645+5158

8.853

18.25

0.57

2.4

2891

55704

56586

J0931−1902

4.638

41.49

3.671

0.6

712

56351

56586

J1012+5307

5.256

9.02

2.283

9.2

11597

53217

56586

J1024−0719

5.162

6.49

1.753

4.0

4830

55118

56586

J1455−3330

7.987

13.57

4.357

9.2

4983

53217

56585

J1600−3053

3.598

52.33

1.015

6.0

7804

54400

56585

J1614−2230

3.151

34.48

1.37

5.1

7323

54724

56585

J1640+2224

3.163

18.46

0.652

8.9

2503

53344

56598

B1855+09

4005

First epoch

Last epoch

(MJD)

(MJD)

53358

56598

J1643−1224

4.622

62.42

3.063

9.0

6921

53291

56585

J1713+0747

4.570

15.99

0.196

8.8

15257

53393

56598

J1738+0333

5.850

33.77

1.367

4.0

2623

55135

56591

J1741+1351

3.747

24.20

0.767

4.2

1588

55042

56593

J1744−1134

4.075

3.14

0.724

9.2

8665

53216

56585

J1747−4036

1.646

152.96

4.597

1.7

2771

55976

56585

J1832−0836

2.719

28.19

1.868

0.6

1131

56354

56585

J1853+1303

4.092

30.57

1.227

5.6

1369

54531

56577

J1903+0327

2.150

297.56

3.085

4.0

1802

55135

56591

J1909−3744

2.947

10.39

0.197

9.1

10259

53292

56598

J1910+1256

4.984

38.06

2.049

8.8

2631

53370

56598

J1918−0642

7.646

26.60

1.819

9.0

9664

53292

56585

J1923+2515

3.788

18.86

4.407

2.2

920

55791

56593

J1944+0907

5.185

24.34

3.190

5.7

1696

54505

56591

J1949+3106

13.138

164.13

4.6

1.2

1409

56139

56593

J2010−1323

5.223

22.16

2.582

4.1

7667

55094

56585

J2017+0603

2.896

23.92

0.718

1.7

1565

55989

56598

J2043+1711

2.380

20.71

0.603

2.3

1382

55758

56591

J2145−0750

16.052

9.03

1.579

9.1

7029

53267

56586

J2214+3000

3.119

22.56

1.656

2.1

2514

55844

56598

J2302+4442

5.192

13.73

5.791

1.7

3037

55972

56586

J2317+1439

3.445

21.90

0.857

8.9

2620

53355

56599

438

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

x 10 −6

4

Raw observational data (J1012+5307) After Vondrak filter

Clock Difference/s

3 2 1 0 −1 −2 −3

53000 53500 54000 54500 55000 55500 56000 56500 57000

MJD/d Fig. 3 The post-fit residuals of PSR J1012+5307 before (dotted line) and after (dot-dashed line) the Vondrak filtering

3

x 10−6

Residual/ns

2 1 0

−1 −2 −3

53000 53500 54000 54500 55000 55500 56000 56500 57000

MJD/d Fig. 4 The phase difference between the raw data and the filtered data

As for the ensemble pulsar time scale shown in Fig.6, its reference time is the Barycentric Dynamic Time (TDB), and its mean square root is 0.056 μs. In the first 4∼5 yr (before 2009), the fluctuation amplitude of the data is considerable, while after 2009, it is evidently reduced. This phenomenon is mainly caused by: (1) the renovation of observational equipment, mainly the treatment terminals GASP and ASP are replaced by GUPPI and PUPPI, respectively. (2) the addition of a large number of data of millisecond pulsar observations with a higher observational accuracy and smaller timing residual in the later period.

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

439

40 35

Number

30 25 20 15 10 5 0 −3

−2

−1

0 Residual/s

1

2

3 x 10 −6

Fig. 5 The distribution of residuals between the raw data and the filtered data

3

x 10 −7

Clock Difference/s

2 1 0 −1 −2 −3 2004

2008

2012

2016

Year

Fig. 6 The ensemble pulsar time scale

The stability of the finally-generated ensemble pulsar time scale is analyzed by the method of σz (τ ) estimation. The obtained σz (τ ) of the ensemble pulsar time scale is given in Fig.7[9] . When the sampling interval is less than 100 day, the stability is worse than 1 × 10−14 , while the sampling interval is longer than 100 day, the stability is better than 1 × 10−14 . It is clear from the result that the long-term frequency stability of the pulsar time scale is very nice, while the short-term one is not so good. What accounts for it is that when the sampling interval is less than 80 day, the frequency flicker noise and random walk noise are the main factors which affect the frequency stability, and the two kinds of noises are of systematic ones, such as the irregularity of pulsar rotation, the noise of the observation system, and/or the interstellar scintillation noise. But when the sampling interval is longer than 80 day, the main component of noises is the frequency white noise, usually caused by the observational noise. In order to raise further the stability and accuracy we need a more

440

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

precise timing noise model, better observational equipment, and a longer observational data set, this will be the direction for the later development of pulsar time scale. It is shown from the data analysis that the timing residuals of several pulsars are obviously characterized by red noises, at present the mechanism of red noise formation is not clear, but some low-frequency observational noises remain to exist after all the pulsar parameters have been fitted[10−11] . It can be deduced from the previous studies that many factors may cause red noises, such as the ephemeris error of the earth’s barycenter adopted in observations, the parameter errors in the theoretical expressions applied for the transformation between TT and Barycentric Coordinate Time (TCB), the physical mechanism of pulsars themselves, and the errors of some adopted constants, etc. All these factors are responsible for the red noises in the timing residuals. In the follow-up work the fitting model will be further improved, and the systematic part of the red noises will be gradually detached. The red noises will be further analyzed so as to improve the fitting accuracy of all the parameters. −13

lg σz(τ)

−13.5

−14

−14.5

−15 5.5

6

6.5 lg(τ/s)

7

7.5

8

Fig. 7 The frequency stability of the ensemble pulsar time scale

5.

CONCLUSION

The stability and accuracy of the pulsar time scale generated from a single pulsar are reduced due to the influences caused by the unevenly sampling of the observational data, the difference of observational time span, and various noise sources. In this paper, a new algorithm to calculate the ensemble pulsar time scale is put forward in order to reduce all these influences, and the experimental verification is implemented by using the real observational data. The experimental results indicate that the new pulsar time scale can effectively reduce the influence of the noises in the pulsar timing residuals, and improve the long-term stability of the pulsar time scale.

YIN Dong-shan et al. / Chinese Astronomy and Astrophysics 41 (2017) 430–441

441

ACKNOWLEDGEMENTS Thanks to Dr. Scott Ranson and Dr. Paul Demorest of the National Radio Astronomical Observatory of the United States for their helps in the observational data. References 1

Zhong C. X., Algorithm of Ensemble Pulsar Time and the Application of Pulsar Time, Xi’an: National Time Service Center, Chinese Academy of Sciences, 2007

2

Hobbs G., Coles W., Machester R. N., et al., MNRAS, 2012, 427, 2780

3

Vondrak J., BAICz, 1977, 28, 84

4

Arzoumanian Z., Brazier A., Burke-spolaor S., et al., ApJ, 2015, 813, 65

5

Hobbs G., Edwards R. T., Manchester R. N., et al., MNRAS, 2006, 369, 655

6

Hobbs G., Jenet F., Lee K. J., et al., MNRAS, 2009, 394, 1945

7

Edwards R. T., Hobbs G., Manchester R. N., MNRAS, 2006, 372, 1549

8

Hobbs G., Archibald A., Arzoumanian Z., et al., CQGra, 2010, 27, 2006

9

Riley W. J., Handbook of Frequency Stability Analysis, Washington, U. S., Government printing office, 2008

10

Yang T. G., Gao Y. P., Tong M. L., et al., Acta Astronomica Sinica, 2015, 56, 370

11

Yang T. G., Gao Y. P., Tong M. L., et al., ChA&A, 2016, 40, 256