Comparison of Pulsar Positions from Timing and Very Long Baseline Astrometry

TL;DR

This study compares pulsar positions obtained from timing methods and VLBI to understand their reference frames, revealing significant differences and the influence of ephemerides and reference source uncertainties on positional accuracy.

Contribution

It provides a detailed comparison of pulsar positions from timing and VLBI, highlighting discrepancies and the impact of different ephemerides on positional measurements.

Findings

Timing and VLBI positions differ significantly.

A change in the obliquity of the ecliptic was detected for JPL DE405.

Position estimates vary with different timing data segments.

Abstract

Pulsar positions can be measured with high precision using both pulsar timing methods and very-long-baseline interferometry (VLBI). Pulsar timing positions are referenced to a solar-system ephemeris, whereas VLBI positions are referenced to distant quasars. Here we compare pulsar positions from published VLBI measurements with those obtained from pulsar timing data from the Nanshan and Parkes radio telescopes in order to relate the two reference frames. We find that the timing positions differ significantly from the VLBI positions (and also differ between different ephemerides). A statistically significant change in the obliquity of the ecliptic of $2.16\pm0.33$\,mas is found for the JPL ephemeris DE405, but no significant rotation is found in subsequent JPL ephemerides. The accuracy with which we can relate the two frames is limited by the current uncertainties in the VLBI reference source positions and in matching the pulsars to their reference source. Not only do the timing positions depend on the ephemeris used in computing them, but also different segments of the timing data lead to varying position estimates. These variations are mostly common to all ephemerides, but slight changes are seen at the 10$\mu$as level between ephemerides.

Tables

Table 1. Table 1: Published VLBI positions for millisecond pulsars corrected to ICRF2. The epoch of the observation is given in modified Julian days (MJDs). The values in parentheses are the 1 σ 𝜎 \sigma errors in the least significant digit.

PSR	Right Ascension	${\sigma }_{\mathrm{RA}}$	Declination	${\sigma }_{\mathrm{Dec}}$	Epoch	Ref
		(mas)		(mas)	(MJD)
J0437$-$4715	04:37:15.883205	0.18	$-$47:15:09.03213	0.16	54100.0	Deller et al. (2008)
J1022$+$1001	10:22:57.995715	1.5	+10:01:52.76475	1	56000.0	Deller et al. (2016)
J1713$+$0747	17:13:49.530594	1.5	+07:47:37.51915	2	52275.0	Chatterjee et al. (2009)
J1939$+$2134	19:39:38.561182	2.7	+21:34:59.13155	2.4	46853.0	Bartel et al. (1996)
J2145$-$0750	21:45:50.458788	1.5	$-$07:50:18.51295	2	56000.0	Deller et al. (2016)

Table 2. Table 2: Published VLBI and new timing positions for millisecond pulsars. Positions are differential with respect to the VLBI ICRF2 position in Table 1. Epochs are as in Table 1. All positions and uncertainties are in mas. The dimensionless normalized cross correlation between RA and Dec is given by C cor subscript 𝐶 cor C_{\rm cor} .

PSR	Type	RA	${\sigma }_{\mathrm{RA}}$	Dec	${\sigma }_{\mathrm{Dec}}$	$C_{\mathrm{cor}}$
J0437$-$4715	ICRF1	0.680	0.86	0.270	0.78	-
	DE405	$-$0.393	0.0092	$-$1.950	0.0065	0.017
	DE414	$-$1.547	0.0090	$-$0.250	0.0066	$-$0.032
	DE421	$-$1.401	0.0082	$-$0.310	0.0059	0.005
	DE430	$-$0.911	0.0082	$-$0.440	0.0059	0.005
	DE432	$-$0.942	0.0082	$-$0.450	0.0059	0.005
	DE435	$-$0.951	0.0082	$-$0.370	0.0059	0.005
J1022+1001	ICRF1	$-$0.224	1.500	0.250	1.000	-
	DE405	22.415	19.383	57.580	49.065	0.999957
	DE414	15.126	16.191	41.180	40.998	0.999973
	DE421	16.809	19.374	45.990	49.044	0.999957
	DE430	15.333	16.227	40.940	41.087	0.999973
	DE432	15.453	16.223	41.230	41.078	0.999973
	DE435	23.292	18.696	61.910	47.329	0.999965
J1713+0747	ICRF1	0.085	1.500	$-$0.153	2.000	-
	DE405	2.840	0.031	4.120	0.065	0.230
	DE414	1.616	0.031	2.450	0.065	0.230
	DE421	1.719	0.035	2.490	0.073	0.234
	DE430	2.063	0.035	2.620	0.074	0.230
	DE432	2.033	0.035	2.630	0.074	0.230
	DE435	1.988	0.035	2.540	0.074	0.230
J1939+2134	ICRF1	-0.276	2.700	-0.046	2.400	-
	DE405	$-$0.095	0.341	4.470	0.455	$-$0.031
	DE414	$-$1.689	0.341	3.300	0.455	$-$0.031
	DE421	$-$1.509	0.341	3.310	0.455	$-$0.031
	DE430	$-$1.236	0.341	3.270	0.455	$-$0.031
	DE432	$-$1.295	0.341	3.280	0.455	$-$0.031
	DE435	$-$1.307	0.341	3.190	0.455	$-$0.031
J2145$-$0750	ICRF1	0.179	1.500	-0.053	2.000	-
	DE405	1.796	0.578	0.860	1.555	$-$0.957
	DE414	0.720	0.578	0.530	1.555	$-$0.957
	DE421	0.651	0.660	0.640	1.773	$-$0.957
	DE430	1.004	0.660	0.540	1.773	$-$0.957
	DE432	1.001	0.660	0.560	1.772	$-$0.957
	DE435	0.917	0.660	0.480	1.773	$-$0.957

Table 3. Table 3: Published VLBI ICRF2 and new timing positions for young pulsars. Timing positions are with respect to VLBI positions and are given in mas.

Pulsar	P0	VLBI		Epoch	Timing - VLBI (mas)
	(s)	RA${}_{\text{ICRF2}}$	Dec${}_{\text{ICRF2}}$	(MJD)	RA${}_{\text{DE430}}$	Dec${}_{\text{DE430}}$	Ref
J0034$-$0721	0.9429	00:34:08.8703(1)	$-$07:21:53.409(2)	52275.0	$-$446(73500)	$-$4609(154000)	1
J0108$-$1431	0.8075	01:08:08.347016(88)	$-$14:31:50.187139(1069)	54100.0	$-$195(360)	723(640)	2
J0139+5814	0.2724	01:39:19.7401(12)	+58:14:31.819(17)	52275.0	$-$223(181)	14(120)	1
J0332+5434	0.7145	03:32:59.3862(10)	+54:34:43.5051(150)	51544.0	$-$136(30)	103(27)	3
J0358+5413	0.1563	03:58:53.71650(330)	+54:13:13.7273(50)	51544.0	2.7(50)	12(11)	1
J0452$-$1759	0.5489	04:52:34.1057(1)	$-$17:59:23.371(2)	52275.0	49(153)	$-$31(221)	1
J0454+5543	0.3407	04:54:07.7506(1)	+55:43:41.437(2)	52275.0	$-$33(32)	23(32)	1
J0538+2817	0.1431	05:38:25.0572(1)	+28:17:09.161(2)	52275.0	3(225)	$-$1539(2100)	1
J0630$-$2834	1.2444	06:30:49.404393(43)	$-$28:34:42.778813(372)	54100.0	$-$33(50)	$-$33(64)	2
J0659+1414	0.3849	06:59:48.1472(7)	+14:14:21.160(10)	51544.0	32(315)	$-$840(1700)	4
J0814+7429	1.2922	08:14:59.5412(10)	+74:29:05.3671(150)	51544.0	-243(285)	-170(86)	2
J0820$-$1350	1.2381	08:20:26.3817(1)	$-$13:50:55.859(2)	52275.0	26(53)	$-$56(84)	1
J0953+0755	0.2530	09:53:09.3071(10)	+07:55:36.1475(150)	51544.0	$-$25(63)	$-$83(161)	3
J1136+1551	1.1879	11:36:03.1829(10)	+15:51:09.7257(150)	51544.0	50(21)	91(34)	3
J1239+2453	1.3824	12:39:40.3589(10)	+24:53:50.0194(150)	51544.0	12(28)	49(37)	3
J1509+5531	0.7396	15:09:25.6298(1)	+55:31:32.394(2)	52275.0	$-$66(69)	7(38)	1
J1543+0929	0.7484	15:43:38.8250(1)	+09:29:16.339(2)	52275.0	$-$45 (180)	209(320)	1
J1932+1059	0.2265	19:32:13.94970(30)	+10:59:32.4198(50)	51544.0	78(17)	$-$34(29)	5
J1935+1616	0.3587	19:35:47.8259(1)	+16:16:39.986(2)	52275.0	$-$6.6(127)	$-$17(19)	5
J2018+2839	0.5579	20:18:03.8332(10)	+28:39:54.1564(150)	51544.0	11(60)	137(68)	3
J2022+2854	0.3434	20:22:37.0712(10)	+28:54:23.0337(150)	51544.0	$-$14(26)	1(28)	3
J2022+5154	0.5291	20:22:49.8655(10)	+51:54:50.3811(150)	51544.0	$-$3(63)	$-$23(43)	3
J2048$-$1616	1.9615	20:48:35.640637(40)	$-$16:16:44.553501(147)	52275.0	460(405)	346(1400)	2
J2157+4017	1.5252	21:57:01.8495(1)	+40:17:45.986(2)	52275.0	$-$111(92)	$-$36(80)	1
J2313+4253	0.3494	23:13:08.6209(1)	+42:53:13.043(2)	52275.0	$-$7.5(825)	12.6(78)	1
(1) Chatterjee et al. (2009); (2) Deller et al. (2009); (3) Brisken et al (2002); (4) Brisken et al (2003); (5) Chatterjee et al (2004)

Table 4. Table 4: Rotation angles (in mas) between VLBI positions and timing positions using different JPL ephemerides for five millisecond pulsars.

Ephermeris	$A_x$	${\sigma }_{A_x}$	$A_y$	${\sigma }_{A_y}$	$A_z$	${\sigma }_{A_z}$	${\chi }_{norm}^2$
DE405	2.36	0.27	0.67	0.61	$-$1.19	0.70	0.93
DE414	0.76	0.26	1.26	0.58	$-$0.02	0.66	0.84
DE421	0.77	0.26	1.11	0.59	$-$0.02	0.67	0.86
DE430	0.87	0.26	0.99	0.58	$-$0.43	0.66	0.83
DE432	0.88	0.26	1.01	0.58	$-$0.42	0.66	0.83
DE435	0.76	0.28	0.92	0.63	$-$0.27	0.71	0.98

Table 5. Table 5: Rotation angles (in mas) after correction for a suspected bias in the right ascension of PSR J0437 − - 4715.

Ephermeris	$A_x$	${\sigma }_{A_x}$	$A_y$	${\sigma }_{A_y}$	$A_z$	${\sigma }_{A_z}$	${\chi }_{norm}^2$
DE405	2.16	0.33	0.13	0.78	-1.36	0.71	0.92
DE414	0.57	0.31	0.73	0.74	-0.19	0.67	0.83
DE421	0.58	0.32	0.59	0.76	-0.18	0.68	0.85
DE430	0.67	0.31	0.47	0.74	-0.59	0.67	0.82
DE432	0.69	0.31	0.49	0.74	-0.58	0.67	0.82
DE435	0.57	0.34	0.40	0.81	-0.44	0.73	0.97

Table 6. Table 6: Rotation angles (in mas) under DE430 and ICRF2 for five millisecond pulsars and all pulsars in the VLBI sample

Data	$A_x$	${\sigma }_{A_x}$	$A_y$	${\sigma }_{A_y}$	$A_z$	${\sigma }_{A_z}$	${\chi }_{norm}^2$
5 MSPs	0.67	0.31	0.47	0.74	-0.59	0.67	0.82
All pulsars	0.68	0.28	0.49	0.63	$-$0.63	0.71	0.86

Table 7. Table 7: Comparison of the rotation angles (in mas) between DE405 and ICRF2 with three different algorithms.

Using all five millsecond pulsars
Alg.	$A_x$	${\sigma }_{A_x}$	$A_y$	${\sigma }_{A_y}$	$A_z$	${\sigma }_{A_z}$	${\chi }_{norm}^2$
1	2.1558	0.2703	0.1296	0.6101	$-$1.3568	0.6898	0.9156
2	2.4728	0.3127	0.9878	0.7600	$-$2.3433	0.8651	0.7717
3	2.2045	0.1905	0.2881	0.4507	$-$1.5362	0.5056	0.6566
Excluding two millisecond pulsars near the ecliptic plane
1	2.6657	0.4304	1.5024	1.0990	$-$2.9308	1.2575	0.7794
2	2.6658	0.4305	1.5027	1.0993	$-$2.9312	1.2578	0.7797
3	2.5191	0.2716	1.1276	0.6964	$-$2.4937	0.7939	0.4832
Alg1=this work; Alg2=RA & Dec uncor.; Alg3=Madison et al. (2013)

Equations

$$x_{\mathrm{eq}}=\cos \alpha \cos \delta ;y_{\mathrm{eq}}=\sin \alpha \cos \delta ;z_{\mathrm{eq}}=\sin \delta .$$

$$\mathbf{D}_{\rm eq}=\left(\begin{array}[]{cc}-\sin{\alpha}\cos{\delta}&-\cos{\alpha}\sin{\delta}\\ \cos{\alpha}\cos{\delta}&-\sin{\alpha}\sin{\delta}\\ 0&\cos{\delta}\end{array}\right).$$

$${\hat{\mathbf{n}}}_{\mathrm{DE}}={\mathbf{\Omega }}_{\mathbf{IC}}^{\mathbf{DE}}{\hat{\mathbf{n}}}_{\mathrm{IC}}$$

$$\Lambda=\left(\begin{array}[]{ccc}0&A_{z}&-A_{y}\\ -A_{z}&0&A_{x}\\ A_{y}&-A_{x}&0\end{array}\right).$$

$$\mathbf{B}=\left(\begin{array}[]{ccc}0&-z&y\\ z&0&-x\\ -y&x&0\end{array}\right),$$

$$(d{\mathbf{\theta }}_{\mathrm{DE}}-d{\mathbf{\theta }}_{\mathrm{IC}})=({\mathbf{D}}_{\mathbf{eq}}^{\mathbf{T}}{\mathbf{D}}_{\mathbf{eq}}{)}^{-1}{\mathbf{D}}_{\mathbf{eq}}^{\mathbf{T}}\mathbf{BA}+{ϵ}_{\theta }.$$

$$\mathrm{RMS}({\theta }_p)=\sqrt{Dt(1000d)}c\zeta /[(1-\zeta )V_{p}1000d].$$

Full text

1 Introduction

Since the early 1980s (e.g., Bartel et al., 1985 and references therein) VLBI observations of pulsars have been providing pulsar positions with sub-milliarcsecond (mas) accuracy. The pulsar timing method is based on observations of pulse times-of-arrival (ToAs) from one or more pulsars. These ToAs are predicted using a model for the astrometric, rotational and orbital properties of each pulsar. If the position of the pulsar is not included with sufficient accuracy in the timing model then the timing residuals (the differences between the predicted and ToAs) will include a sine term with a period of a year. Pulsar timing packages, such as tempo2 (Hobbs, Edwards & Manchester 2006), can be used to fit for the amplitude and phase of the sinusoidal signal and, hence, improve the modelled position of the pulsar. As ToAs for pulsars with millisecond periods can often be measured with a precision of $\sim$100 ns the position of such pulsars can also be determined with sub-mas accuracy using these timing methods (see, e.g., Edwards, Hobbs & Manchester 2006).

Pulsar timing arrays use millisecond pulsars to search for the signatures of ultra-low-frequency gravitational waves (e.g., Shannon et al. 2015) and provide exceptionally precise source positions. These measurements (from the Parkes Pulsar Timing Array project; Manchester et al. 2013) form the basis of the work that we are describing here. Unfortunately there are not many VLBI measurements published for millisecond pulsars because such pulsars are relatively weak.

Younger pulsars with periods of order 1 s are more common and often brighter than the millisecond pulsars, but it is usually not possible to measure their pulse ToAs with such accuracy and the long-term stability of these pulsars is poorer (see, e.g., Hobbs et al. 2010 for a description of the timing noise in 366 young pulsars). However, there are many VLBI determinations of their positions. One of the goals of this project was to determine if inclusion of a much larger number of young pulsars would improve the results that could be obtained using many fewer millisecond pulsars. As we will show, results from the young pulsars were consistent with those from the millisecond pulsars, but they did not improve the frame-tie precision.

Fomalont et al. (1984) compared positions obtained using interferometry and the timing method. They noted that the differences in the pulsar positions obtained via these two methods could result from many causes including calibration errors in the interferometry and instabilities in the pulsar rotation or errors in the planetary ephemeris for the timing method. Bartel et al. (1996) computed rotation matrices to transform between various ephemerides, including the JPL planetary ephemeris DE200111The JPL ephemerides are described in ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/ExplSupplChap8.pdf.. Since the year 2003 (starting with the DE405 ephemeris) the reference frame for the ephemerides has been the International Celestial Reference Frame (ICRF) as VLBI measurements of planets and spacecraft were used in their construction. For completeness we note the following major changes in the ephemerides since DE405222based on ftp://ssd.jpl.nasa.gov/pub/eph/planets/README.txt:

• •

DE414: Created in 2005 this included an upgraded orbit of Pluto and ranging data from the MGS and Odyssey spacecraft.

• •

DE421: Created in 2008 this ephemeris included new fits to planetary and lunar laser ranging data

• •

DE430: From 2013, this was referenced to the ICRF 2.

• •

DE432: This was created in 2014 specifically for the New Horizons mission for Pluto. It did not include effects relating to nutation.

• •

DE435: Created in 2016 for the Cassini project and primarily provides an updated orbit for Saturn333see ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/de435.iom.pdf for details..

Pulsar timing observations analyzed with these ephemerides therefore should agree with VLBI observations. The most recent analysis was carried out by Madison et al. (2013) who did not detect a significant rotation between the frames. Here we extend their work and we find statistically significant differences in the obliquity of the ecliptic for the more recent JPL solar system ephemerides DE405 through 435 compared to ICRF2, but show that these differences are dominated by the VLBI and timing position of a single pulsar, PSR J0437$-$4715. When we update the position of this pulsar with more recent information, we only find a statistically significant difference for DE405.

As described by Madison et al. (2013), this type of work is of special importance for pulsar timing. The use of VLBI measurements in pulsar timing arrays would improve the sensitivity of the data sets both for gravitational wave searches and for searches for unknown solar-system bodies (such as the postulated planet X; see e.g., Batygin & Brown 2016).

In order to ensure that the work presented here is reproducible we have made our data files and MATLAB code publically available. Details are provided in Appendix B.

2 Observations

The timing data sets were obtained from the Nanshan 25m-diameter telescope (see Wang et al. 2001 for an introduction to the Nanshan observing system) and the Parkes 64m-diameter telescope (see Reardon et al. 2016 for a description of the Parkes data sets that we used here). All the VLBI pulsar positions are from published work (Bartel et al. 1996; Brisken et al. 2002, 2003; Chatterjee et al. 2001, 2004, 2009; Deller et at. 2008, 2009, 2012, 2013, 2016). Throughout this paper, uncertainties are denoted as $\sigma$ and represent 68% confidence intervals.

2.1 Timing data analysis

Regular timing observations of 74 pulsars using the Nanshan telescope of the Xinjiang Astronomical Observatory commenced in 2000 January in the 18 cm observing band. A cryogenic receiver system has been used since July 2002. The total bandwidth of 320 MHz is divided into 128 channels for each circular polarization and digitized at 1 ms intervals. The observation time for each pulsar is usually between 4 and 16 minutes. We selected from these pulsars the 25 for which VLBI positions have been published.

The Parkes pulsar timing array project (PPTA) has, since 2005, been regularly timing around 20 millisecond pulsars in the 10, 20 and 50 cm observing bands (Manchester et al. 2013; Reardon et al. 2016). We selected the five for which VLBI positions have been published. We have used the PPTA observations after they have been combined with earlier timing observations in the 20 cm band for some of the pulsars (Verbiest et al. 2009). Our timing analysis is based on the tempo2 software package. As described by Hobbs, Edwards & Manchester (2006) this package corrects all known phenomena at the $\sim$1 ns level. This requires accounting for variations in the Earth’s rotation. These are accounted for within tempo2 using the C04 series of the Earth Orientation Parameters (EOP) that are published by the International Earth Rotation Service (IERS). Details are available in McCarthy & Petit (2004). tempo2 also requires the telescope positions on the Earth (which are determined with high accuracy in geocentric coordinates) and many other effects. Details are provided in the Hobbs, Edwards & Manchester (2006) and Edwards, Hobbs & Manchester (2006) series of papers about the tempo2 software.

The earlier observations cannot be corrected for fluctuations in the interstellar medium (ISM) because they are at a single band but they can be analyzed together with the new ISM-corrected observations using the “split-Cholesky” method discussed by Reardon et al. (2016). We followed their method and used their noise models to re-run the timing analysis using different planetary ephemerides. We also extracted the covariance matrix of each pulsar’s right ascension (RA) and declination (Dec), which was not published by Reardon et al. (2016) and changed the timing position epoch to match that of the VLBI observation. It should be noted that the proper motion of many pulsars greatly exceeds the uncertainty in position over the data span, so the measurement of the proper motion is important. When the VLBI measurement is at the edge of the timing span over even outside it, the uncertainty in the proper motion can be the primary limitation.

The Nanshan pulsars were reanalyzed and new noise models obtained. The white noise was estimated using the efacEquad plugin in the tempo2 package to rescale the ToA uncertainties so that they better represent the observed scatter in the residuals (Wang et al. 2015). The red noise was estimated using the spectralModel plugin (Coles et al. 2011). The epoch of the timing position for each pulsar was set to be the same as that of the VLBI measurement. As before, we obtained the astrometric parameters using a range of different planetary ephemerides and we obtained the covariance matrix of RA and Dec for each analysis.

2.2 VLBI data

The most recent available observations of PSR J1939+2134 (Bartel et al. 1996) also provide a weighted average of their position with previous observations. We use this weighted average as it reduces the uncertainties by about 50% and improves the overall analysis by about 10%. In most cases the authors provide the formal RA and Dec uncertainties including the known positional uncertainties in the calibration sources used. However in the case of PSR J0437$-$4715 the uncertainties were given independently. In this case we summed the quoted differential uncertainty and the uncertainty in the calibrator source in quadrature. However it was completely dominated by the uncertainty in the calibrator position. In fact all the VLBI position uncertainties are dominated by the ICRF reference sources. In most cases the VLBI positions were referenced to ICRF1 (Ma et al. 1998). The observations of PSR J1939+2134 were referenced to the International Earth Rotation Service 1993 reference frame. Since this (IERS, 1993) was the immediate precursor to ICRF1 we have treated them as the same. The same calibration sources are now available in ICRF2 (Fey et al. 2009) so we have corrected the source positions to ICRF2. We did not correct the published uncertainties to those in ICRF2 except for PSR J0437$-$4715 because we cannot separate the calibration error from other errors. It is clear that the VLBI position uncertainties are likely to be underestimated, but by exactly how much is difficult to quantify and will vary between pulsars. When converting from ICRF1 to ICRF2 we typically halve the error in the position of the reference source, but the remaining errors remain unchanged. For our work we assume that the position uncertainties may be underestimated by as much as 50% for pulsars other than PSR J0437$-$4715, but we do not attempt to further quantify such errors. We note that the root-mean-square (rms) difference between ICRF1 and ICRF2 is $\sim$0.2  mas for our millisecond pulsars. The rotation between ICRF1 and ICRF2 has been estimated by Fey (2009) who found it consistent with zero with an uncertainty of $\sim$0.02 mas.

2.3 The effects of ephemerides and ICRF frames

Timing positions depend on the planetary ephemeris used. The strongest millisecond pulsar, PSR J0437$-$4715, shows this clearly. In Figure 1, we plot the apparent position according to VLBI referenced to both ICRF1 and ICRF2. We also plot the apparent positions of this pulsar from pulsar timing using planetary ephemerides DE405 through DE435 (Standish 1998; Folkner et al. 2008). DE405 was the first ephemeris to be referenced to ICRF1 and DE430 was the first referenced to ICRF2. The error bars on the timing positions are extremely small and have been increased in the figure by a factor of 10 to be visible. All timing positions determined with different ephemerides differ by much more than their statistical uncertainties, but positions determined using DE405 are particularly inconsistent. This pattern also applies to the other millisecond pulsars in our data set. In general the timing position uncertainty is dominated by ephemeris differences just as the ICRF uncertainty dominates the VLBI position uncertainty.

Since ICRF2 positions are available for all VLBI calibrators we use them as the reference throughout the rest of this paper. These VLBI positions for the five millisecond pulsars we used are listed in Table 1. The differential positions (with respect to ICRF2) are given in Table 2 for VLBI under ICRF1 and for timing with ephemerides DE405 through 435. The timing measurements of RA and Dec for PSR J1022+1001 are highly correlated and the uncertainties reflect this. It was necessary to carry six significant figures in the cross correlation to analyze this pulsar. We note that the ecliptic longitude is well-defined and this pulsar does contribute significantly to the frame tie.

The young pulsars are listed in Table 3. Even though we calculated positions with respect to all the planetary ephemerides, for brevity we only tabulate our results with respect to JPL DE430, the most recent major adjustment of the JPL ephemerides. Position determinations for these young pulsars are not sufficiently precise to show differences between ephemerides.

3 Estimation of the Frame-Tie Rotation

To find the rotation that relates the celestial and ecliptic frames we follow the procedure described by Madison et al. (2013), keeping the same notation (their section 5.1). Since the rotation angles and the measurement errors are small ($<10^{-8}$ rad) the transformation can be linearized with no sacrifice in accuracy. We assume that the observations, both timing and VLBI, are provided in equatorial coordinates, $\bm{\theta}_{\rm eq}=(\alpha,\delta)$ (where $\alpha$ and $\delta$ are the right ascension and declination respectively), and we convert these into direction cosines $\hat{\bm{n}}_{\rm eq}=(x_{\rm eq},y_{\rm eq},z_{\rm eq})$

[TABLE]

The differentials can be written as $d\hat{\bm{n}}_{\rm eq}=\mathbf{D_{\rm eq}}\ d{\bm{\theta}}_{\rm eq}$ where

[TABLE]

The position vectors in the ICRF (IC) and timing (DE) frames (both column vectors) can be related by a rotation matrix, $\mathbf{\Omega}$, as shown below.

[TABLE]

For the very small angles considered here the rotation matrix can be written as $\mathbf{\Omega^{DE}_{IC}=I+\Lambda}$ where $\mathbf{I}$ is the identity matrix and

[TABLE]

Here $A_{x}$, $A_{y}$ and $A_{z}$ are the (right-handed) Euler rotation angles about the x, y, and z axes respectively. Note that the x axis is towards the vernal equinox so $A_{x}$ represents an offset in the obliquity of the ecliptic.

In matrix form we can rewrite the offset between the timing and VLBI positions as $d\hat{\bm{n}}_{\rm DE}-d\hat{\bm{n}}_{\rm IC}=\mathbf{BA}+{\bm{\epsilon}}_{n}$. Here

[TABLE]

The direction cosines ($x$,$y$,$z$) in the equation above could be either the timing or VLBI coordinates as they agree to $\sim$9 significant figures. $\mathbf{A}$ is the column vector $(A_{x},A_{y},A_{z})$ and ${\bm{\epsilon}}_{n}$ is the vector of the uncertainties in the difference between the two position vectors. The covariance matrix of $d\hat{\bm{n}}$ is $\mathbf{D_{\rm eq}}\rm{Cov}({\bm{\theta}}_{\rm eq})\ \mathbf{D_{\rm eq}^{T}}$ and $\Sigma$, the covariance matrix of ${\bm{\epsilon}}_{n}$, is the sum of the covariance matrix of $d\hat{\bm{n}}_{\rm DE}$ and that of $d\hat{\bm{n}}_{\rm IC}$.

This is a linear least squares problem and is easily solved if $\Sigma$ is known. Madison et al. (2013) made the simplifying assumption that $\Sigma$ is diagonal and computed the variances on the diagonal by propagating the errors from ${\bm{\theta}}_{\rm eq}$ to $\hat{\bm{n}}_{\rm eq}$ using our Equation (1). This assumption creates two problems. $\Sigma$ is not diagonal and would not be diagonal even if ${\rm Cov}({\bm{\theta}}_{\rm eq})$ were diagonal. Furthermore $\Sigma$ is singular and cannot be inverted as required by the solution given by Equation (12) in Madison et al. (2013).

These problems can be minimized by extracting the full ${\rm Cov}({\bm{\theta}}_{\rm eq})$ from the tempo2 fit, computing $\Sigma$ from ${\rm Cov}({\bm{\theta}}_{\rm eq})$ and using the Moore-Penrose pseudo-inverse in place of $\Sigma^{-1}$ (see equation 18 in Zhu, et al., 2014). Unfortunately we do not have the full ${\rm Cov}({\bm{\theta}}_{\rm eq})$ for the VLBI measurements, but these are never as highly correlated as are the timing positions.

We tested this procedure, but we also found a more straightforward method in which we altered the least squares equations $d\hat{\bm{n}}_{\rm DE}-d\hat{\bm{n}}_{\rm IC}=\mathbf{BA}+{\bm{\epsilon}}_{n}$ to avoid the singularity in the covariance matrix $\Sigma$. This was done by making the substitution $\mathbf{D_{eq}}\ (d{\bm{\theta}}_{\rm DE}-d{\bm{\theta}}_{\rm IC})=\mathbf{BA}+{\bm{\epsilon}}_{n}$, then pre-multiplication of both sides by $\mathbf{D_{eq}^{T}}$. New least squares equations are found by premultiplication of both sides by $\mathbf{(D_{eq}^{T}D_{eq})^{-1}}$. The new least squares problem applies directly to the observations rather than the transformed observations:

[TABLE]

Here $\mathbf{D_{eq}^{T}D_{eq}}$ is diagonal and is trivially inverted. The product $\mathbf{(D_{eq}^{T}D_{eq})^{-1}D_{eq}^{T}}$ is the pseudo inverse of $\mathbf{D_{eq}}$. The results from this method are numerically identical to the original method.

One of the pulsars in our set, PSR J1022+1001, is very close to the ecliptic, and so $\alpha$ and $\delta$ are highly correlated. This pulsar would be better analyzed in ecliptic coordinates because the ecliptic latitude and longitude are almost uncorrelated. However we found that we could solve the least squares problem in equatorial coordinates using double precision. We did this to avoid potential confusion in defining the ecliptic reference frame.

4 Results for the frame-tie rotation

Prior to applying our method to our data, we tested our procedure with the positions of the two millisecond pulsars (PSRs J0437$-$4715 and J1713+0747) that were processed by Madison et al. (2013). We found, as expected, negligible difference between our rotation matrix and theirs. In this case, with only two pulsars, there only four data points and there are three parameters to be obtained. Since there is only one degree of freedom, the problem is more like an inversion than a least squares fit and the weighting of the observations has very little effect on the resulting estimate.

Below we first present our results using only the five millisecond pulsars and then demonstrate the effect of including the young pulsars.

4.1 Millisecond Pulsars

In Table 4, we present the results for all five millisecond pulsars, using six different JPL ephemerides. There is no significant rotation about the y or z axes in any ephemeris, but there is a consistently significant rotation about the x axis, $A_{x}$, i.e., an offset in the obliquity of the ecliptic. This is $2.4\pm 0.3$ mas in DE405, but was significantly corrected in DE415 and has remained at $\sim 0.8\pm$ 0.3 mas in more recent ephemerides. The normalized $\chi^{2}$ is consistently slightly below unity, probably because the VLBI positions were slightly improved by re-referencing them to ICRF2, but, as noted earlier, the published errors were not reduced when making this correction. Overall the $\chi_{norm}^{2}$ indicates that $\rm{Cov}({\bm{\theta}}_{\rm eq})$ describes the uncertainty in the observations ${\bm{\theta}}_{\rm eq}$ reasonably well.

The positions before and after rotation are shown in Figure 2 for the most recent ephemeris, DE435. The timing position marked with an ‘+’ is defined to be the origin. The original VLBI position in ICRF2 is marked with an ‘x’, and after rotation it is marked with a square. One can see that the fitting process puts the transformed VLBI position of PSR J0437$-$4715 very close to its timing position. The VLBI position of the other pulsars are also transformed to be consistent, within an error bar, with their timing positions, but the error bars are much larger. PSR J0437$-$4715 is the brightest millisecond pulsar in the sky so both its timing and VLBI positions are exceptionally well-measured and it dominates the fit.

To further understand the dominance of PSR J0437$-$4715 we removed it from the array and repeated the fit. This increased both $A_{x}$ and $\sigma_{Ax}$, to 2.8 and 1.5 mas respectively. However when we remove the other pulsars from the fit, one at a time, we find that each contributes significantly to reducing the uncertainty, $\sigma_{Ax}$. Indeed the rotation $A_{x}$ did not exceed $3\sigma_{Ax}$ comfortably unless all five pulsars were included. Because the result depends so strongly on PSR J0437$-$4715 we discussed this VLBI measurement with the lead author of the corresponding publication (Deller, private communication, 2017). He noted that continuing unpublished VLBI observations suggest that there is a calibration error of $\sim+0.8$ mas in the published observation of PSR J0437$-$4715. Accordingly we recomputed the frame tie assuming this bias and increasing the uncertainty by 0.8 mas in quadrature because the bias is not yet well understood. The revised results, presented in Table 5, no longer show a significant rotation (except for the DE405 ephemeris). We will use this correction to the VLBI position of PSR J0437$-$4715 in the remainder of this work.

4.2 Inclusion of Young Pulsars

We were able to include all the young pulsars although their positional uncertainties are as much as $10^{5}$ times larger than those of the millisecond pulsars. We did not expect all to contribute, but including the very noisy pulsars provided a test of the numerical stability of the algorithm. We did find a small problem in estimating the uncertainty. As is the normal practice we scale the covariance matrix of the estimated parameters by the normalised $\chi^{2}$-value ($\chi_{norm}^{2}$). However adding young pulsars sometimes reduced the $\chi_{norm}^{2}$ below unity and caused the covariance matrix to be underestimated. We modified the scaling so that it is only applied when the $\chi_{norm}^{2}$ is greater than unity which solved this problem. The result of including all the young pulsars is given in Table 6. One can see that they are completely consistent with the millisecond pulsars.

We also searched for outliers among the young pulsars by appending them, one at a time, to the five millisecond pulsars. The resulting estimates of $A_{x}$ had a mean and standard deviation of $0.68\pm 0.01$ mas. As the standard deviation of $A_{x}$ is 0.3 mas the variation of $A_{x}$ when including different young pulsars is completely negligible. We conclude that young pulsars would not bias a frame tie, but they will not be useful unless the timing error can be greatly improved. Simply including more young pulsars will not help.

4.3 Effect of Assuming that the Covariance Matrix is Diagonal

We use the positions determined with the JPL DE405 ephemeris to provide a basis for comparing different assumptions about the covariance matrix of the observations. We chose this ephemeris are there is a detectable frame-tie rotation in $A_{x}$ that provides a “test signal”. We compare the results obtained using three different assumptions: (1) our algorithm, (2) our algorithm neglecting the cross correlation between RA and Dec and (3) the Madison et al. (2013) algorithm which assumes that the matrix $\Sigma$ is diagonal. We then exclude the two pulsars near the ecliptic (PSRs J1022+1001 and J2145$-$0750) and repeat the comparison. The results are in Table 7.

When we exclude the near-ecliptic pulsars cases (1) and (2) are identical because the RA and Dec are almost uncorrelated. In this comparison case (3) is similar to cases (1) and (2) but the $\chi_{norm}^{2}$ and all the other errors are underestimated because there are only 3 degrees of freedom in the residuals and the algorithm assumes six degrees of freedom. When we include the two near-ecliptic pulsars the results are more complex. Case (2) gets the signal $A_{z}$ reasonably well but overestimates $A_{y}$ and $A_{z}$. In fact $A_{z}$ becomes marginally significant. Case (3) is worse, as one might expect. It gets the signal reasonably well but underestimates the error, with $A_{z}$ becoming a 3$\sigma$ detection.

In summary, the need for including the off-diagonal terms increases when they become larger. In this particular data set it became important for two of the five millisecond pulsars. This is easily done for the timing measurements. We urge VLBI observers to make an effort to specify the correlation between their right ascension and declination measurements.

4.4 Effect of the Number of Pulsars

As adding more young pulsars will not improve the frame tie, here we estimate the improvement that could be gained by adding more millisecond pulsars. If we had $N_{p}$ pulsars with equal, uncorrelated errors in RA and Dec, then we would have $2N_{p}$ independent measurements from which we need to estimate the three independent parameters (the three Euler angles). Accordingly we would expect the variance of each Euler angle to decrease as $1/(2N_{p}-3)$. We tested this by placing $N_{p}=$ 2, 5, 10, 20, 50 and 100 pulsars at uniformly distributed random positions on the sky. As we have assumed that the noise in RA and Dec is independent we explicitly ignore the possibility that some pulsars may be near the ecliptic. We assumed a 10-mas rms independent Gaussian error in the differential RA and Dec. We then calculated the rotation angles using our method. We simulated 20 realizations of the sky for each $N_{p}$, and 1000 realizations of the noise for each realization of the sky. For each sky realization we found the rms rotation angles over the 1000 simulations, $\sigma_{A_{x}}$ etc. From these we determined the rms rotation angle for each sky realization $\sigma_{A}=(\sigma_{A_{x}}^{2}+\sigma_{A_{y}}^{2}+\sigma_{A_{z}}^{2})^{0.5}$. The mean of $\sigma_{A}$ over the 20 sky realizations is plotted as a function of $N_{p}$ in Figure 3. The error bars on Figure 3 are the rms of $\sigma_{A}$ over the sky realizations. The line drawn through the simulated measurements is $\sigma_{A}=30{\rm mas}/(2N_{p}-3)^{0.5}$ and matches well with our simulations.

We can use this equation to estimate how many pulsars would be required to halve the uncertainty in our estimate of the obliquity. Two of our five millisecond pulsars are close to the ecliptic so they together effectively only contribute one independent measurement. Accordingly we currently have four “typical” millisecond pulsars. We would need roughly 11 such millisecond pulsars to halve our current rms uncertainty.

5 Time variations in the timing positions

Our ability to perform the frame tie is limited by the precision of the VLBI positions, but the timing positions are far more precise and clearly show the effects of different ephemerides. This suggests that a search for time variations in the apparent timing positions might provide a useful measure of the accuracy of the ephemerides because the timing positions (at a reference epoch) should not vary with time. We have done this for the two bright pulsars, PSRs J0437$-$4715 and J1713+0747. Obtaining accurate timing positions requires several years of observation so we break our data into blocks of four years overlapping by two years. These data have been corrected for fluctuations in ISM dispersion where possible, but the early data could not be corrected because there was only data available in the 20 cm observing band.

The mean positions for both pulsars vary significantly with the ephemeris used (DE405, 414, 421, and 430) but the time variation in PSR J1713+0747 is not statistically significant and that in PSR J0437$-$4715 is only marginally significant. The time variations are highly correlated between the ephemerides, so we have taken an average of the most recent (DE414, 421 and 430) as a reference. This reference position is plotted in the top panel of Figure 4 for PSR J0437$-$4715 and similarly in Figure 5 for PSR J1713+0747. We have estimated how much position variation could occur in the early observations of PSR J0437$-$4715 because of uncorrected fluctuations in the ISM and it is about a factor of 10 too weak to explain the observed variation (our calculation is presented in Appendix A). On the other hand the points between 2000 and 2004 which exceed 2$\sigma$ are not statistically independent because the blocks overlap by half. We conclude that the observed time variation is probably simply measurement error.

However the measurement error is the same for all ephemerides, so we can search for smaller differences between the ephemerides. We compare the positions determined by each ephemeris with the reference position and plot these in the lower 4 panels of Figures 4 and 5. These show much smaller variation, of order 10 $\mu$as. As the variations between ephemerides are highly correlated the errors are also highly correlated and the differences are much more significant that would be suggested by the error bars on the top panel. The variation in DE405 in both pulsars appears to be real and this would be reasonable because we can easily detect frame tie errors in this ephemeris. It is not clear that the variations in the more recent ephemerides are real although the fluctuations in DE421 appear to be smallest in both pulsars.

The variations between pulsar positions determined with different ephemerides, both in the mean and the temporal variations, demonstrate the sensitivity of pulsar timing and its potential for improving planetary ephemerides.

6 Discussion and conclusions

Data contributing to the JPL ephemerides are dominated by ranging observations between solar-system bodies and have only recently begun to include VLBI measurements of the apparent position of solar-system bodies which can tie the celestial frame of VLBI observations to the ecliptic coordinates of the heliosphere (Folkner & Border 2015). The present accuracy of the frame tie is believed to be $\sim 0.25$ mas (Folkner, private communication 2017). So it is not surprising that we find an error in the obliquity of the ecliptic of 2.16$\pm 0.33$ mas in DE405, but that any errors in subsequent ephemerides are not statistically significant. However our work shows that pulsar timing observations can be used to detect differences between ephemerides of order $10\mu$as which suggests that they may be useful in refining future ephemerides.

We would like to use the VLBI position, proper motion and parallax in the pulsar timing models in order to remove the requirement to fit for these astrometric parameters. Currently such fitting removes all signals with a frequency near 1 yr*-1* and much of the power near 2 yr*-1*. Such signals potentially include ephemeris errors, variations in the relativistic correction from terrestrial to barycentric time and gravitational waves from individual supermassive binary black holes. However our work shows that present VLBI positions are not yet accurate enough to improve the frame tie or to improve the sensitivity of pulsar timing arrays.

In the near future improvements in the frame tie can be expected primarily from VLBI observations because there are already $\sim$50 precisely-timed millisecond pulsars in the International Pulsar Timing Array (IPTA, Verbiest et al. 2016). It will take many years to add more millisecond pulsars to this set. However we expect that significantly more VLBI observations of these pulsars will soon become available, e.g. through the PSR$\pi$ project (Deller et al. 2016). It is also possible that the Large European Array for Pulsars (LEAP) project could provide astrometric information on a some of these millisecond pulsars.

In the more distant future, more pulsars will be precisely timed with new instruments and, hopefully, some will be bright enough that precise VLBI positions can be determined. Pulsar timing programs are planned for both the Five Hundred Metre Spherical Telescope (FAST) currently being commissioned in China and the Square Kilometre Array (SKA).

7 Acknowledgements

The Parkes radio telescope is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by the Commonwealth Scientific and Industrial Research Organisation (CSIRO). This work is supported by West Light Foundation of CAS (No.XBBS201322) and National Natural Science Foundation of China (Nos. 11403086, U1431107, 1137300), the Strategic Priority Research Programme (B) of the Chinese Academy of Sciences (No. XDB230102000) and GH acknowledge support from the Australian Research Council (ARC) Future Fellowship. We acknowledge comments on early drafts of the manuscript from Dr A. Deller. This work used the astronomy & astrophysics package for Matlab (Ofek 2014).

Appendix A Expected position variations from turbulent phase gradients in the ISM

A transverse phase gradient in the ISM would cause a change in the apparent position of a pulsar. Such gradients will vary statistically due to the natural spectrum of turbulent variations in electron density. The resulting phase variations $\phi(\bm{r})$ can be described statistically by the phase structure function $D_{p}(\bm{s})=\langle(\phi(\bm{r}+\bm{s})-\phi(\bm{r}))^{2}\rangle$ where the angle brackets denote an ensemble average (You et al, 2007, Keith et al, 2013). Accordingly we can estimate the rms phase gradient over a spatial scale s, $\nabla\phi=\sqrt{D_{p}(\bm{s})}/|{\bm{s}}|$ and the resulting angular shift is given by $\bm{\theta}_{s}=\nabla\phi/k$ where $k=2\pi F/c$. The ISM models used by Keith et al. (2013) are given as a time delay structure function $D_{t}(\bm{s})=D_{p}(\bm{s})/(2\pi F)^{2}$ where $\bm{s}=\bm{V}T$ and $D_{t}(\bm{s}=\bm{V}\times 1000~{}{\rm days})$ is tabulated. Here $\bm{V}$ is the velocity of the line-of-sight with respect to the ISM.

The velocity $\bm{V}$ and the apparent position shift of the pulsar ${\bm{\theta}}_{p}$ depend on the location of the scattering region. We define the fractional distance from the pulsar to the scattering region as $\zeta$. Then we have $V=V_{p}(1-\zeta)$ and $\theta_{p}=\theta_{s}\zeta$. Here $V_{p}$ is the proper motion of the pulsar. Finally we have

[TABLE]

Putting in values for PSR J0437$-$4715 of $D_{t}(1000)=1.6\mu s^{2}$ (Keith et al, 2013), $V_{p}=105$ km/s (Reardon et al. 2016), and $\zeta=0.22$ (Bhat, Coles, private communication, 2015) we obtain $\theta_{p}=0.0024$ mas. Thus the probability that this process (with an rms of 0.0024 mas) produces the observed single oscillation of $\sim$0.1 mas peak to peak is low.

Appendix B Reproducing our work

Our software scripts and the input data files are available from the CSIRO Data Access Portal (DAP; https://data.csiro.au). The input pulsar timing data files for the millisecond pulsars are available from http://doi.org/10.4225/08/561EFD72D0409. The data were processed using tempo2 using the different solar system ephemerides described in this paper. The resulting tempo2 positions are recorded in a text file msp_timing. The VLBI positions are obtained from the literature and stored in msp_vlbi. The calibration source positions for the millisecond pulsars are stored in cal.txt. Similar files exist for the analysis of the young pulsars. We then use MATLAB to obtain the rotation angles. All these files can be obtained from the CSIRO DAP (http://doi.org/10.4225/08/58a0e1593c5be). Our MATLAB code requires the convertdms.m, coco.m, rotm_coo.m, cosined.m and obliquity.m files from the Astromatlab library and errorbare.m a plot routine for drawing error bars444The errorbare.m file can be obtained from https://au.mathworks.com/matlabcentral/fileexchange/23465-enhanced-errorbar-function/content/errorbare/errorbare.m and the Astromatlab library is described by Ofek (2014) and available from https://webhome.weizmann.ac.il/home/eofek/matlab/..