11institutetext: School of Astronomy and Space Science, Key Laboratory of
Modern Astronomy and Astrophysics (Ministry of Education), Nanjing University,
Nanjing 210023, P. R. China
11email<EMAIL_ADDRESS>22institutetext: School of Earth Sciences and
Engineering, Nanjing University, Nanjing 210023, P. R. China 33institutetext:
Institute of Astrophysics, Foundation for Research and Technology-Hellas,
Voutes, 71110 Heraklion, Greece
33email<EMAIL_ADDRESS>44institutetext: Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

# Systematics of planetary ephemeris reference frames inferred from pulsar
timing astrometry

N. Liu 1122 Z. Zhu 11 J. Antoniadis 3344 J.-C. Liu 11 H. Zhang 11

(Received; accepted)

###### Abstract

Aims. This study aims to investigate the systematics in planetary ephemeris
reference frames through pulsar timing observations.

Methods. We used the published data sets from several pulsar timing arrays and
performed timing analyses for each pulsar using different planetary
ephemerides retrieved from the Jet Propulsion Laboratory’s Development
Ephemeris (DE), Ephemeris of Planets and the Moon (EPM), and INPOP
(Intégration Numérique Planétaire de l’Observatoire de Paris). Then, we
compared the timing solutions and modeled the differences in position and
proper motion by vector spherical harmonics of the first degree. The timing
solutions were also compared with those determined by very long baseline
interferometry (VLBI) astrometry.

Results.  The orientation offsets between the latest editions of the DE, EPM,
and INPOP series do not exceed 0.4 milliarcseconds (mas), while the relative
spins between these ephemerides are less than 5 microarcseconds per year
($\mathrm{\mu as\,yr^{-1}}$). We do not detect significant glides in either
position or proper motion between these ephemerides. The orientation of the
pulsar timing frames deviates from that of the VLBI frame from zero by
approximately $\mathrm{0.4\,mas}$ when considering the formal uncertainty and
possible systematics.

Conclusions. The orientation of current planetary ephemeris frames is as
accurate as at least 0.4 mas, and the nonrotating is better than
$\mathrm{5\,\mu as\,yr^{-1}}$.

###### Key Words.:

Astrometry – Ephemerides – Reference systems – (Stars:) pulsars: general

##  1 Introduction

Modern numerical planetary ephemerides have been used for various purposes in
deep space missions (e.g., Yang et al. 2022), fundamental physics (e.g.,
Pitjeva et al. 2021; Bernus et al. 2019; Bernus et al. 2020, 2020, 2022), dark
matter (e.g., Pitjeva & Pitjev 2013), and planetary science (e.g., Liu &
Capitaine 2017; Fienga et al. 2020a, 2020b; Yuan et al. 2021). There are three
main ephemeris series, namely, Jet Propulsion Laboratory’s (JPL) planetary and
lunar ephemerides Development Ephemeris (DE) (Park et al. 2021), Ephemeris of
Planets and the Moon (EPM) from the Institute of Applied Astronomy of the
Russian Academy of Sciences (Pitjeva et al. 2022), and INPOP (Intégration
Numérique Planétaire de l’Observatoire de Paris) ephemeris from the IMCCE
(Institut de mécanique céleste et de calcul des éphémérides) at Paris
Observatory (Fienga et al. 2022). These planetary ephemerides realize the
International Celestial Reference System (ICRS; Feissel & Mignard 1998) in a
dynamical sense and are aligned onto the International Celestial Reference
Frame (ICRF; Charlot et al. 2020) with an accuracy of sub-milliarcsecond (mas)
by very long baseline interferometry (VLBI) and very long baseline array
(VLBA) observations of the planetary spacecraft. These observations mainly
include tracking measurements of MESSENGER (MErcury Surface, Space
ENvironment, GEochemistry, and Ranging) spacecraft orbiting Mercury (Verma et
al. 2014); Venus Express (VEX) and Magellan spacecraft orbiting Venus (Folkner
et al. 1993, 1994a, 1994b); Mars Global Surveyor, Mars Odyssey, and Mars
Reconnaissance Orbiter spacecraft orbiting Mars (Park et al. 2015); Galileo
and Juno spacecraft orbiting Jupiter (Jones et al. 2021); and Cassini
spacecraft orbiting Saturn (Jones et al. 2011, 2015, 2020). The orbits of
inner planets are tied to ICRF by VLBI measurements of Mars-orbiting
spacecraft (Park et al. 2021).

Understanding the systematic errors of the planetary ephemerides is essential
for their users. The largest uncertainty in planetary positions comes from the
lacks in the dynamical modelings (Standish & Fienga 2002). There are at least
two methods to determine the uncertainty of the planetary ephemerides (Fienga
et al. 2009). The first is to compare the extrapolated positions of planets
with the observations outside the fit interval of the ephemerides (e.g.,
Fienga et al. 2013). The second is to assess the internal agreements between
ephemerides from the same series that are constructed by identical or similar
dynamical models (e.g., Standish 1990a) or external agreements between the
state-of-the-art ephemerides from different series (e.g., Pitjeva 2005). These
agreements are usually evaluated by comparing the heliocentric or geocentric
coordinates of planets such as Earth or Mars (e.g., Tkachuk & Choliy 2013;
Deng et al. 2014).

The pulsar timing technique precisely records the pulse times-of-arrival
(TOAs) of pulsars, from which information about the astrometry, rotation, and
orbital motion can be inferred. One of the main objectives of pulsar timing
astrometry is to search for nanohertz (nHz) gravitational wave signals in the
timing observations of an ensemble of pulsars forming a pulsar timing array
(PTA; Burke-Spolaor et al. 2019). The planetary ephemerides are used to
compute the theoretical TOAs. It has been noted for a long time that the
pulsar timing positions obtained by using different planetary ephemerides in
the pulsar timing analysis vary, for instance, as found in Kaspi et al.
(1994). Several authors reported that using updated versions of planetary
ephemerides yielded an improved solution, leading to a more accurate
measurement of astrometric parameters such as the parallax (Splaver et al.
2005; Hotan et al. 2006; Verbiest et al. 2008). As the precision of the timing
astrometry is continuously improved to the level of a few tens of
microarcseconds ($\mathrm{\mu as}$), the systematic error in the planetary
ephemerides is recognized as one of the main error sources that may degrade
the sensitivity of the detection of gravitational waves (e.g., Chen et al.
2021).

Pulsar timing astrometry can also contribute to the investigation of
systematics in planetary ephemerides with the advantage that pulsar timing
observations are fully independent of the creation process of the planetary
ephemerides. The signal in the residual of TOAs due to the systematics in the
planetary ephemerides, such as the uncertainties in the orbits of outer
planets, can be modeled in the timing analysis (Caballero et al. 2018;
Vallisneri et al. 2020). The comparison of pulsar positions determined from
VLBI and pulsar timing observations provides a direct tie between the
planetary frames and the extragalactic frames (Bartel et al. 1996; Wang et al.
2017; Liu et al. 2023). It is also possible to compare different ephemerides
via pulsar timing astrometry to provide an independent check on the external
agreements between planetary ephemerides, as done in Fienga et al. (2011).

In this paper, we extend the related work of Fienga et al. (2011) by using
recent data releases from several PTA data sets with increased sample sizes
and improved astrometric precision and accuracy (Sect. 2.2), updated VLBI
solutions of pulsars (Sect. 2.3), as well as modeling using the vector
spherical harmonics (VSH; Mignard & Klioner 2012) of the ephemeris systematics
(Sect. 2.4). Three motivations drive us to perform this study: (i) to
investigate the realistic accuracy of current planetary ephemeris frames, (ii)
to provide a transformation relation between planetary ephemerides, especially
those produced by different groups, and (iii) to evaluate the implication of
using different planetary ephemerides in the timing solutions on the derived
astrometric parameters for pulsars. For these purposes, we compare the timing
solutions based on the same sets of timing observations using different
ephemerides (Sect. 3.1) and compare the timing solutions with the latest VLBI
astrometric solutions (Sect. 3.3). Throughout the paper, the planetary
ephemeris reference frame is referred to as the reference frame based on
Earth’s orbit provided by a given planetary ephemeris. The formal
uncertainties (including the error bars) are always provided at a confidence
level of 68 percent. The name of pulsars is given by their J2000 coordinates.

##  2 Materials and methods

###  2.1 Numerical planetary ephemerides

Figure 1:  Release timeline of the planetary ephemerides of the DE, EPM, and
INPOP series. The release time is chosen to be the relevant date found in the
online documents or publications accompanied by each ephemeris and may not be
rigorously correct.  Table 1:  VLBI measurements of Mars spacecraft and
asteroid modeling used for each planetary ephemeris considered in this work.

Ephemeris | VLBI observations of Mars | Dynamical Modeling | Ref
---|---|---|---
| Period | Number | Accuracy | Main belt | Kuiper belt |
|  |  | ($\mathrm{mas}$) |  |  |
DE200 | … | … | … | 3 | … | 1
DE405 | 1989 | 2 | 10–100 | 300 | … | 2
DE410 | 1989–2003 | 33 | 0.5 | 300 | … | 3
DE421 | 2001–2007 | 94 | ¡0.5 | 343 | … | 4–5
DE430 | 2001–2013 | 151 | ¡0.2 | 343 | … | 6
DE440 | 2001–2020 | 316 | 0.18–0.25 | 343 | 30 + 1 ring | 7
EPM2011 | 1989–2010 | 144 | 0.8 | 301 + 1 ring | 21 + 1 ring | 8
EPM2015 | 1989–2013 | 204 | 0.18–6.2 | 301 + 1 ring | 30 + 1 ring | 9
EPM2017 | 1989–2014 | 204 | 0.18–6.2 | 301 + 3 rings | 30 + 3 rings | 10
EPM2021 | 1989–2014 | 204 | 0.18–6.2 | 277 + 3 rings | 30 + 3 rings | 11
INPOP06c | 1989–2003 | 44 | 0.5 | 300 + 1 ring | … | 12
INPOP08a | 1989–2007 | 96 | 0.4 | 303 + 1 ring | … | 13
INPOP10a | 1989–2007 | 96 | 0.4 | 161 + 1 ring | … | 14
INPOP10b | 1989–2007 | 96 | 0.4 | 161 + 1 ring | … | 15
INPOP10e | 1989–2007 | 96 | 0.4 | 161 + 1 ring | … | 16
INPOP13c | 1989–2007 | 96 | 0.4 | 139 | … | 17
INPOP17a | 1989–2013 | 194 | 0.3 | 168 | … | 18
INPOP19a | 1989–2013 | 194 | 0.3 | 343 | 3 rings | 19
INPOP21a | 1989–2013 | 194 | 0.3 | 343 | 500 | 20
111 The accuracy of the VLBI observations is given by a priori accuracy or the
scatter (root-mean-squared or standard deviation) of postfit residuals in the
literature.

(1) Standish (1990b); (2) Standish (1998); (3) Standish (2003); (4) Folkner et
al. (2007); (5) Folkner et al. (2009); (6) Folkner et al. (2014); (7) Park et
al. (2021); (8) Pitjeva (2013); (9) Pitjeva (2017); (10) Pitjeva & Pitjev
(2018); (11) Pitjeva et al. (2022); (12) Fienga et al. (2008); (13) Fienga et
al. (2009); (14) Fienga et al. (2011); (15) Fienga et al. (2012); (16) Fienga
et al. (2013); (17) Fienga et al. (2014); (18) Viswanathan et al. (2017); (19)
Fienga et al. (2019); (20) Fienga et al. (2021).

We downloaded the successive editions of the planetary ephemerides from the DE
series222https://ssd.jpl.nasa.gov/planets/eph_export.html, the EPM
series333https://iaaras.ru/en/dept/ephemeris/epm/, and the INPOP
series444https://www.imcce.fr/recherche/equipes/asd/inpop/. Only the planetary
ephemerides that worked properly with the timing analysis software were
included in this work. We only considered the JPL DE ephemerides since DE403
because they were aligned onto the extragalactic reference frame, except for
DE200, which was constructed on its own dynamical reference frame of J2000
(Standish 1982). This ephemeris was considered because it was intensively used
in the previous timing analyses. For many pulsars (especially young,
nonrecycled pulsars), the most recent timing solutions were given under the
DE200 frame. Including DE200 in the comparison can thus provide a convenient
transformation relation from these timing astrometric solutions referred to
the DE200 frame to those referred to more recent planetary ephemerides. Only
these DE ephemerides for a general purpose were considered, that is, DE200,
DE405, DE421, DE430, and DE440, except for DE410, which was used for
comparison with previous results (Sect. 4.1.1). The EPM series contained
EPM2011, EPM2015, EPM2017, and EPM2021. For the INPOP series, we considered
the latest version of each generation, which are INPOP06c, INPOP08a, INPOP10e,
INPOP13c, INPOP17a, INPOP19a, and INPOP21a. In addition, INPOP10a and INPOP10b
were included to illustrate the influence of different factors on the
planetary ephemeris frames (Sect. 4.2). In total, we considered six DE
ephemerides, four EPM ephemerides, and nine INPOP ephemerides.

Figure 1 shows a rough timeline of the release dates of the ephemerides. One
can expect that the ephemerides created at similar times may have comparable
precision and accuracy, which has been validated in previous studies (e.g.,
Fienga et al. 2008). Table 1 lists the major differences in the VLBI
observations of Mars spacecraft and asteroid perturbation modeling of the Main
belt objects and Kuiper belt objects between planetary ephemerides. These
differences are supposed to be the most relevant to the differences between
planetary ephemerides and the alignment accuracy of their frame axes onto
ICRF, which are of interest in this work.

###  2.2 Pulsar timing data sets

Figure 2:  Sky distribution of the pulsars contained in different PTA data
sets in the equatorial coordinate system. The dotted curve indicates the
location of the ecliptic plane.  Table 2:  Astrometric precision of the timing
solutions based on different PTA data sets using DE440.

PTA | Subset | Pos. Epoch | Nb. PSR | $\sigma_{\alpha^{*}}$ | $\sigma_{\delta}$ | $\sigma_{\mu_{\alpha^{*}}}$ | $\sigma_{\mu_{\delta}}$
---|---|---|---|---|---|---|---
|  |  |  | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$)
PPTA DR2e | All | 2011.8 | 24 | 47 | 89 | 12 | 22
| $|\beta|\\!\leq\\!10^{\circ}$ | 2011.3 | 4 | 101 | 1469 | 29 | 427
| $|\beta|\\!>\\!10^{\circ}$ | 2012.0 | 20 | 40 | 74 | 9 | 16
EPTA DR1 | All | 2009.5 | 42 | 158 | 359 | 69 | 165
| $|\beta|\\!\leq\\!10^{\circ}$ | 2009.5 | 12 | 197 | 2255 | 120 | 1046
| $|\beta|\\!>\\!10^{\circ}$ | 2009.5 | 30 | 140 | 302 | 67 | 136
NANOGrav12.5-WB | All | 2013.6 | 47 | 82 | 156 | 44 | 89
| $|\beta|\\!\leq\\!10^{\circ}$ | 2012.5 | 6 | 212 | 977 | 84 | 245
| $|\beta|\\!>\\!10^{\circ}$ | 2013.8 | 41 | 65 | 140 | 39 | 80
NANOGrav12.5-NB | All | 2013.6 | 47 | 58 | 127 | 38 | 79
| $|\beta|\\!\leq\\!10^{\circ}$ | 2012.4 | 6 | 230 | 996 | 82 | 247
| $|\beta|\\!>\\!10^{\circ}$ | 2013.8 | 41 | 58 | 108 | 30 | 58
IPTA DR2A | All | 2010.2 | 64 | 161 | 320 | 60 | 154
| $|\beta|\\!\leq\\!10^{\circ}$ | 2009.8 | 15 | 254 | 2171 | 119 | 924
| $|\beta|\\!>\\!10^{\circ}$ | 2010.3 | 49 | 147 | 248 | 45 | 77
IPTA DR2B | All | 2010.2 | 64 | 145 | 314 | 59 | 136
| $|\beta|\\!\leq\\!10^{\circ}$ | 2009.8 | 15 | 281 | 1763 | 140 | 977
| $|\beta|\\!>\\!10^{\circ}$ | 2010.3 | 49 | 109 | 173 | 44 | 70
Combination | All | 2011.4 | 80 | 107 | 207 | 49 | 107
| $|\beta|\\!\leq\\!10^{\circ}$ | 2010.4 | 15 | 240 | 1465 | 102 | 581
| $|\beta|\\!>\\!10^{\circ}$ | 2011.6 | 65 | 86 | 146 | 39 | 70

555 $\sigma_{\alpha^{*}}\\!=\\!\sigma_{\alpha}\cos\delta$ and
$\sigma_{\mu_{\alpha^{*}}}\\!=\\!\sigma_{\mu_{\alpha}}\cos\delta$. The first
four columns tabulate the PTA name, subset, mean position epoch, and the
number of pulsars, followed by the typical precision (median) of the position
and proper motion measurements. The last three rows represent the combination
of all PTA data sets above.

We retrieved the pulsar timing data sets from several pulsar timing arrays,
including data release 1.0 of the European Pulsar Timing Array (EPTA DR1;
Desvignes et al. 2016), the extended second data release of the Parkes Pulsar
Timing Array (PPTA DR2e; Reardon et al. 2021), the narrowband and wideband
versions of the 12.5 yr data set of the North American Nanohertz Observatory
for Gravitational Waves (NANOGrav12.5-NB and NANOGrav12.5-WB; Alam et al.
2021a, 2021b), and two versions (A and B) of the second data release of the
International Pulsar Timing Array (IPTA DR2A and IPTA DR2B; Perera et al.
2019). All of these timing observations were made between 1986 and 2018. While
the data set in the EPTA DR1, PPTA DR2e, and NANOGrav12.5 are obtained from
independent observations of three regional PTAs in Europe, North America, and
Australia, the IPTA combines the data sets from these individual PTAs to
improve the sensitivity. The combination also leads to an increasing number of
pulsars and better sky coverage compared to individual PTAs. For IPTA DR2, the
input TOAs came from the previous data releases of these PTAs, which are EPTA
DR1, the 9 yr data set from NANOGrav (Matthews et al. 2016), and the first
data release of PPTA and its extended version (Reardon et al. 2016). The main
difference between the two versions of IPTA DR2 (i.e., IPTA DR2A and IPTA
DR2B) lies in the modeling of the dispersion measure and handling of the noise
properties of pulsars. In this work, EPTA DR1, PPTA DR2e, and NANOGrav12.5
were used in addition to IPTA DR2 for two reasons. The first one was that PPTA
DR2e and NANOGrav12.5 surpassed their previous data releases used to generate
IPTA DR2 in terms of astrometric accuracy and sample size. There were five
pulsars in PPTA DR2e and ten pulsars in NANAGrav12.5 that were not included in
IPTA DR2. The other reason was that different PTAs with independent
observations and analyses could be used to examine the consistency between
results of planetary ephemeris comparison based on pulsar timing observations.
The timing solutions from all the data sets were combined to derive the final
results.

Since the planetary ephemerides are mainly constrained by observations on the
ecliptic plane, it would be interesting to use only pulsars right on or close
to the ecliptic plane. However, the major axis of the error ellipse tends to
align with the ecliptic longitude for timing astrometry. The timing
astrometric error ellipse would become extremely elongated for pulsars right
on the ecliptic, making the measurement of right ascension and declination
highly correlated and leading to a less reliable timing astrometric solution.
One example in our sample is PSR J1022+1001, whose position is usually
provided in the ecliptic coordinate system in the timing solution. We failed
to reproduce timing solutions using several planetary ephemerides for this
pulsar. This pulsar might not be suitable for our analyses and thus was
removed from our sample. We believe that missing one pulsar should not affect
our results much because we studied the common features yielded in the timing
solutions of all pulsars. To compare the results based on pulsars near or far
from the ecliptic plane, we divided the sample into three subsets: all
pulsars, pulsars close to the ecliptic plane (i.e., ecliptic latitude
$|\beta|\\!\leq\\!10^{\circ}$), and pulsars far from the ecliptic plane (i.e.,
$|\beta|\\!>\\!10^{\circ}$). Figure 2 depicts the sky distribution of all
pulsars in the equatorial coordinate system, from which we found that the
sample was dominated by pulsars at high ecliptic latitudes.

The timing models for each pulsar given in each PTA data set were directly
used to rerun the timing analysis, except that different planetary ephemerides
were adopted. We processed the published TOAs of all PTA data sets by
TEMPO2666https://bitbucket.org/psrsoft/tempo2/src/master/ (version 2022.05.1;
Hobbs et al. 2006). Noting that the planetary masses used in TEMPO2 (defined
as constant variables) usually differed from those used in the planetary
ephemerides, we manually modified the values of the related variables in
TEMPO2 to be consistent with each planetary ephemeris to be studied. The
timing astrometric solutions of all pulsars in each PTA under the same
planetary ephemeris formed a pulsar catalog, which was considered as an
independent representation of the reference frame of the used ephemerides.

Since the DE ephemerides were widely used in the pulsar timing analyses, we
chose the timing solutions using DE440 as the reference solution for each PTA
data set. Table 2 presents a brief overview of these reference solutions. The
typical timing position precision (median) was approximately 50 ${\rm\mu
as}$–150 $\mathrm{\mu as}$ in right ascension and 100 $\mathrm{\mu as}$–350
$\mathrm{\mu as}$ in declination. For the proper motion measurements, the
typical precision was in the range of 10 $\mathrm{\mu as\,yr^{-1}}$ to 70
$\mathrm{\mu as\,yr^{-1}}$ in right ascension and twice as poor in
declination. Comparing two samples of $|\beta|\\!\leq\\!10^{\circ}$ and
$|\beta|\\!>\\!10^{\circ}$, we found that the astrometric precision in either
right ascension or declination for pulsars near the ecliptic plane was twice
or more worse than those for pulsars away from the ecliptic plane, which is a
typical feature of timing astrometry.

Assuming the timing precision merely consists of the measurement noise, the
sensitivity to the systematics of the ephemeris frames to be examined would be
improved by a factor of the root square of the number of pulsars in the
sample. Using the values presented in Table 2, we obtained 8 ${\rm\mu as}$ to
24 $\mathrm{\mu as}$ for right ascension, 18 ${\rm\mu as}$ to 55 $\mathrm{\mu
as}$ for declination, 2 $\mathrm{\mu as\,yr^{-1}}$ to 11 $\mathrm{\mu
as\,yr^{-1}}$ for proper motion in right ascension, and 4 $\mathrm{\mu
as\,yr^{-1}}$ to 26 $\mathrm{\mu as\,yr^{-1}}$ for proper motion in
declination. The detection sensitivity in the position system was thus well
below the orientation accuracy of the planetary ephemeris reference frames
(i.e., $\sim$0.25 mas), enabling us to perform meaningful investigations into
the possible systematics in the planetary ephemeris reference frames.

###  2.3 VLBI solutions of pulsars

Table 3: Overview of timing and VLBI solutions for seven pulsars used for
evaluating the rotation between ephemerides frames and extragalactic frame.

Pulsar | Timing | VLBI | Ref.
---|---|---|---
Epoch | $\sigma_{\alpha^{*}}$ | $\sigma_{\delta}$ | $\sigma_{\mu_{\alpha^{*}}}$ | $\sigma_{\mu_{\delta}}$ | Eopch | $\sigma_{\alpha^{*}}$ | $\sigma_{\delta}$ | $\sigma_{\mu_{\alpha^{*}}}$ | $\sigma_{\mu_{\delta}}$
|  | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$) |  | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$) |
J0437$-$4715 | 2009.5 | 5 | 5 | 1 | 2 | 2007.0 | 1015 | 994 | 50 | 90 | 1
J1012+5307 | 2009.5 | 49 | 64 | 10 | 14 | 2016.9 | 900 | 1000 | 90 | 140 | 2
J1640+2224 | 2009.5 | 46 | 65 | 13 | 19 | 2016.3 | 971 | 1000 | 80 | 140 | 3
J1713+0747 | 2009.5 | 5 | 10 | 1 | 2 | 2002.0 | 1486 | 2000 | 170 | 160 | 4
J2010$-$1323 | 2009.5 | 148 | 614 | 56 | 237 | 2012.2 | 1459 | 4000 | 329 | 303 | 5
J2145$-$0750 | 2009.5 | 97 | 268 | 25 | 69 | 2012.2 | 1486 | 4000 | 52 | 90 | 5
J2317+1439 | 2009.5 | 88 | 173 | 25 | 49 | 2012.2 | 1451 | 1000 | 465 | 704 | 5

777The timing solutions were obtained through a timing analysis of the IPTA
DR2B data set using DE440. The last column gives the reference to the VLBI
position.

(1) Deller et al. (2009); (2) Ding et al. (2020); (3) Vigeland et al. (2018);
(4) Chatterjee et al. (2009); (5) Deller et al. (2019).

We found the VLBI astrometric solutions for seven pulsars in our sample (PSR
J1022+1001 excluded), which mostly came from the PSR$\pi$ and MSPSR$\pi$
campaigns (Deller et al. 2019). These sources were all included in the IPTA
DR2 sample. Among these seven pulsars, six sources were common to the EPTA DR1
and NANOGrav samples, and two pulsars were found in PPTA DR2e. There are two
pulsars with $|\beta|\\!\leq\\!10^{\circ}$, which are PSR J2010$-$1323 and PSR
J2145$-$0750.

We noted that the positional uncertainty in Deller et al. (2009) did not
include the systematics due to the calibrator position uncertainty, core
shift, and phase referencing errors from the primary calibrator to the in-beam
calibrator. Therefore, we inflated the positional uncertainty therein using
the same method as described in Liu et al. (2023). That is, we added in
quadrature to the positional uncertainty given in Deller et al. (2009) an
empirical value of 0.8 mas in each coordinate for the core-shift, 1.26 mas in
right ascension and 0.59 mas in declination accounting for the phase
referencing error, and the error in the absolute position of the calibrator.
The positional precision of the timing and VLBI astrometric solutions is
displayed in Table 3. The VLBI positional uncertainty now reaches
$\sim\mathrm{1\,mas}$, while the timing positional uncertainty is at least ten
times smaller. The correlations between right ascension and declination in the
VLBI solutions were not available and were assumed to be zero.

The Gaia data also provided the astrometric solutions for pulsars in the Gaia
reference frame (Gaia-CRF3; Gaia Collaboration et al. 2022a). We followed the
method given in Antoniadis (2021) to search for pulsars in Gaia Data Release 3
(Gaia Collaboration et al. 2022b). We only found two pulsars, preventing us
from performing a robust comparison between the planetary ephemeris frames and
the Gaia reference frame.

###  2.4 Modeling of systematic differences in position and proper motion

The pulsar timing position and proper motion differences between using
different ephemerides and using DE440 were computed as

$\displaystyle\Delta\alpha^{*}$ $\displaystyle=(\alpha_{\rm e}-\alpha_{\rm
r})\cos\delta_{\rm r},$ (1) $\displaystyle\Delta\delta$
$\displaystyle=\delta_{\rm e}-\delta_{\rm r},$ (2)
$\displaystyle\Delta\mu_{\alpha^{*}}$ $\displaystyle=\mu_{\alpha^{*},{\rm
e}}-\mu_{\alpha^{*},{\rm r}},$ (3) $\displaystyle\Delta\mu_{\delta}$
$\displaystyle=\mu_{\delta,{\rm e}}-\mu_{\delta,{\rm r}},$ (4)

where the subscripts (as well as superscripts shown below) “r” and “e”
represent the reference solution using DE440 and the timing solutions using
ephemerides other than DE440, respectively. The notation
$\mu_{\alpha^{*}}\\!=\\!\mu_{\alpha}\cos\delta$, together with other similar
notations (e.g., $\sigma_{\alpha^{*}}\\!=\\!\sigma_{\alpha}\cos\delta$), was
used throughout the paper. We used the four-dimensional vector in the form of

$\bm{y}=(\Delta\alpha^{*},\,\Delta\delta,\,\Delta\mu_{\alpha^{*}},\,\Delta\mu_{\delta})^{\rm
T}$ (5)

as the input observable. The full covariance matrix of the observable was
expressed as

$\bm{C}_{\bm{y}}=\left[\begin{array}[]{cccc}\sigma^{2}_{\Delta\alpha^{*}}&C_{\Delta\alpha^{*},\Delta\delta}&C_{\Delta\alpha^{*},\Delta\mu_{\alpha^{*}}}&C_{\Delta\alpha^{*},\Delta\mu_{\delta}}\\\
C_{\Delta\alpha^{*},\Delta\delta}&\sigma^{2}_{\Delta\delta}&C_{\Delta\delta,\Delta\mu_{\alpha^{*}}}&C_{\Delta\delta,\Delta\mu_{\delta}}\\\
C_{\Delta\alpha^{*},\Delta\mu_{\alpha^{*}}}&C_{\Delta\delta,\Delta\mu_{\alpha^{*}}}&\sigma^{2}_{\Delta\mu_{\alpha^{*}}}&C_{\Delta\mu_{\alpha^{*}},\Delta\mu_{\delta}}\\\
C_{\Delta\alpha^{*},\Delta\mu_{\delta}}&C_{\Delta\delta,\Delta\mu_{\delta}}&C_{\Delta\mu_{\alpha^{*}},\Delta\mu_{\delta}}&\sigma^{2}_{\Delta\mu_{\delta}}\end{array}\right],$
(6)

where

$\displaystyle\sigma_{\Delta\alpha^{*}}$
$\displaystyle=\sqrt{\sigma^{2}_{\alpha^{*},{\rm
e}}+\sigma^{2}_{\alpha^{*},{\rm r}}},$ (7)
$\displaystyle\sigma_{\Delta\delta}$
$\displaystyle=\sqrt{\sigma^{2}_{\delta,{\rm e}}+\sigma^{2}_{\delta,{\rm
r}}},$ (8) $\displaystyle\sigma_{\Delta\mu_{\alpha^{*}}}$
$\displaystyle=\sqrt{\sigma^{2}_{\mu_{\alpha^{*}},{\rm
e}}+\sigma^{2}_{\mu_{\alpha^{*}},{\rm r}}},$ (9)
$\displaystyle\sigma_{\Delta\mu_{\delta}}$
$\displaystyle=\sqrt{\sigma^{2}_{\mu_{\delta},{\rm
e}}+\sigma^{2}_{\mu_{\delta},{\rm r}}},$ (10) $\displaystyle
C_{\Delta\alpha^{*},\Delta\delta}$ $\displaystyle=\rho^{\rm
e}_{\alpha,\delta}\sigma_{\alpha^{*},{\rm e}}\sigma_{\delta,{\rm e}}+\rho^{\rm
r}_{\alpha,\delta}\sigma_{\alpha^{*},{\rm r}}\sigma_{\delta,{\rm r}},$ (11)
$\displaystyle C_{\Delta\alpha^{*},\Delta\mu_{\alpha^{*}}}$
$\displaystyle=\rho^{\rm e}_{\alpha,\mu_{\alpha^{*}}}\sigma_{\alpha^{*},{\rm
e}}\sigma_{\mu_{\alpha^{*}},{\rm e}}+\rho^{\rm
r}_{\alpha,\mu_{\alpha^{*}}}\sigma_{\alpha^{*},{\rm
r}}\sigma_{\mu_{\alpha^{*}},{\rm r}},$ (12) $\displaystyle
C_{\Delta\alpha^{*},\Delta\mu_{\delta}}$ $\displaystyle=\rho^{\rm
e}_{\alpha,\mu_{\delta}}\sigma_{\alpha^{*},{\rm e}}\sigma_{\mu_{\delta},{\rm
e}}+\rho^{\rm r}_{\alpha,\mu_{\delta}}\sigma_{\alpha^{*},{\rm
r}}\sigma_{\mu_{\delta},{\rm r}},$ (13) $\displaystyle
C_{\Delta\delta,\Delta\mu_{\alpha^{*}}}$ $\displaystyle=\rho^{\rm
e}_{\delta,\mu_{\alpha^{*}}}\sigma_{\delta,{\rm
e}}\sigma_{\mu_{\alpha^{*}},{\rm e}}+\rho^{\rm
r}_{\delta,\mu_{\alpha^{*}}}\sigma_{\delta,{\rm
r}}\sigma_{\mu_{\alpha^{*}},{\rm r}},$ (14) $\displaystyle
C_{\Delta\delta,\Delta\mu_{\delta}}$ $\displaystyle=\rho^{\rm
e}_{\delta,\mu_{\delta}}\sigma_{\delta,{\rm e}}\sigma_{\mu_{\delta},{\rm
e}}+\rho^{\rm r}_{\delta,\mu_{\delta}}\sigma_{\delta,{\rm
r}}\sigma_{\mu_{\delta},{\rm r}},$ (15) $\displaystyle
C_{\Delta\mu_{\alpha^{*}},\Delta\mu_{\delta}}$ $\displaystyle=\rho^{\rm
e}_{\mu_{\alpha^{*}},\mu_{\delta}}\sigma_{\mu_{\alpha^{*}},{\rm
e}}\sigma_{\mu_{\delta},{\rm e}}+\rho^{\rm
r}_{\mu_{\alpha^{*}},\mu_{\delta}}\sigma_{\mu_{\alpha^{*}},{\rm
r}}\sigma_{\mu_{\delta},{\rm r}},.$ (16)

We used $\sigma$ and $\rho$ to represent the formal uncertainty and
correlation coefficient, respectively. To quantify the significance of the
position and proper motion difference, we adopted a normalized quantity
$X_{\bm{y}}$ defined as

$X_{\bm{y}}^{2}=\left[\begin{array}[]{llll}\Delta\alpha^{*}&\Delta\delta&\Delta\mu_{\alpha^{*}}&\Delta\mu_{\delta}\end{array}\right]\bm{C}_{\bm{y}}^{-1}\left[\begin{array}[]{l}\Delta\alpha^{*}\\\
\Delta\delta\\\ \Delta\mu_{\alpha^{*}}\\\
\Delta\mu_{\delta}\end{array}\right].$ (17)

We considered the global differences between both the position and proper
motion systems of the planetary ephemerides, which would manifest themselves
in the systematic changes of the pulsar timing positions and proper motions
when different planetary ephemerides were used in the timing analyses. These
global differences were modeled by the VSH of the first degree, including a
rotation vector and a glide vector.

For the position system, the rotation vector $\bm{R}\\!=\\!(R_{\rm X},R_{\rm
Y},R_{\rm Z})^{\rm T}$ mainly represents the orientation offsets between the
planetary ephemeris reference frames at the coordinate epoch (i.e., J2000). On
the other hand, the glide vector $\bm{G}\\!=\\!(G_{\rm X},G_{\rm Y},G_{\rm
Z})^{\rm T}$ models a dipolar positional offset field that can be caused by
the large-scale deformation in the planetary ephemeris reference frames and
relative displacements and motions of the Solar System Barycenter (SSB)
defined implicitly by the planetary ephemerides. The pulsar positional offsets
at J2000 modeled by these two vectors can be expressed as

$\displaystyle\Delta{\alpha^{*}_{0}}=$ $\displaystyle-R_{\rm X}\cos\alpha_{\rm
r}\sin\delta_{\rm r}-R_{\rm Y}\sin\alpha_{\rm r}\sin\delta_{\rm r}+R_{\rm
Z}\cos\delta_{\rm r}$ $\displaystyle-G_{\rm X}\sin\alpha_{\rm r}+G_{\rm
Y}\cos\alpha_{\rm r},$ (18) $\displaystyle\Delta{\delta}_{0}=$
$\displaystyle+R_{\rm X}\sin\alpha_{\rm r}-R_{\rm Y}\cos\alpha_{\rm r}$
$\displaystyle-G_{\rm X}\cos\alpha_{\rm r}\sin\delta_{\rm r}-G_{\rm
Y}\sin\alpha_{\rm r}\sin\delta_{\rm r}+G_{\rm Z}\cos\delta_{\rm r}.$ (19)

For the proper motion system, the rotation vector $\bm{r}\\!=\\!(r_{\rm
X},r_{\rm Y},r_{\rm Z})^{\rm T}$ characterizes the relative spin of the
planetary ephemeris reference frames, which can be attributed to the mean
motion difference in the orbital determinations of the Earth. Meanwhile, the
glide vector $\bm{g}\\!=\\!(g_{\rm X},g_{\rm Y},g_{\rm Z})^{\rm T}$ can be
related to the relative accelerations of the SSB defined by different
planetary ephemerides (if they exist). Similar to Eqs.(18)–(19), we obtained
the pulsar proper motion offsets as

$\displaystyle\Delta\mu_{\alpha^{*},{\rm C}}=$ $\displaystyle-r_{\rm
X}\cos\alpha_{\rm r}\sin\delta_{\rm r}-r_{\rm Y}\sin\alpha_{\rm
r}\sin\delta_{\rm r}+r_{\rm Z}\cos\delta_{\rm r}$ $\displaystyle-g_{\rm
X}\sin\alpha_{\rm r}+g_{\rm Y}\cos\alpha_{\rm r},$ (20)
$\displaystyle\Delta\mu_{\delta,{\rm C}}=$ $\displaystyle+r_{\rm
X}\sin\alpha_{\rm r}-r_{\rm Y}\cos\alpha_{\rm r}$ $\displaystyle-g_{\rm
X}\cos\alpha_{\rm r}\sin\delta_{\rm r}-g_{\rm Y}\sin\alpha_{\rm
r}\sin\delta_{\rm r}+g_{\rm Z}\cos\delta_{\rm r}.$ (21)

In addition, the resulting proper motion offsets can also lead to a positional
displacement if the reference epoch $t$ of the pulsar position differs from
the coordinate epoch $t_{0}$ of the planetary ephemeris frame (i.e., J2000).
Therefore, the total positional offsets due to the global differences of the
aforementioned planetary ephemeris reference frame can be written as

$\displaystyle\Delta\alpha^{*}_{C}$
$\displaystyle=\Delta\alpha^{*}_{0}+\Delta\mu_{\alpha^{*},C}\left(t-t_{0}\right),$
(22) $\displaystyle\Delta\delta_{C}$
$\displaystyle=\Delta\delta_{0}+\Delta\mu_{\delta,C}\left(t-t_{0}\right).$
(23)

The VSH parameters (components in vectors $\bm{R}$, $\bm{G}$, $\bm{r}$, and
$\bm{g}$) were determined by a least-squares fit of the model given by Eqs.
(20)–(23) to the observables described by Eqs. (1)–(4). The observables were
weighted by the inverse of the covariance matrix $\bm{C}_{\bm{y}}$ given in
Eq. (6). To improve the robustness of the fitting against the potential
outliers, we computed the normalized difference $X_{\bm{y}}$ based on the
postfit residuals following Eq. (17). The sources were considered outliers
when the normalized differences exceeded the median value of the normalized
differences for all sources by a chosen factor $\kappa$ (called the clip
limit), as used in Gaia Collaboration et al. (2021). The fit was iterated
until no new outlier was detected. We set $\kappa\\!=\\!3$ in this work.
Different values of $\kappa$ were examined: the estimates of the VSH
parameters changed little. We also adopted a bootstrap fit to cross-check the
results of the least-squares fit. We generated 1000 trial samples, for which
we randomly selected the same number of pulsars from the input sample with
replacements. The least-squares fit was performed for each trial sample. The
median value and the median absolute deviation (MAD) multiplied by a factor of
1.4826 (i.e., normalized MAD, which is equivalent to 1-sigma in a Gaussian
distribution) were used as the final estimate and the associated uncertainty
for each unknown.

For the comparison between timing and VLBI astrometry, we first propagated the
VLBI positions from their own epoch to those of the timing positions using the
VLBI proper motions and then computed the offsets of the timing positions with
respect to the propagated VLBI positions. The uncertainties of these
positional offsets in either right ascension or declination were computed as
the square root of the quadratic sum of both uncertainties in VLBI and timing
positions as well as the propagation errors induced by the formal error in
VLBI proper motion. Since there were a limited number of pulsars in common, we
only modeled these offsets with the orientation offset, that is, the
contribution from the vector $\bm{R}$, as given in Eqs. (18)–(19).

We performed the fit separately for each PTA data set to check their
consistencies. The final estimates, together with their uncertainties, were
obtained by a least-squares fit to the combination of the timing solutions
from all PTAs.

##  3 Results

###  3.1 Offsets of the pulsar timing position and proper motion due to
different ephemerides

Figure 3:  Offsets of the pulsar timing positions in the EPM2021 frame with
respect to those in the DE440 frame as a function of right ascension (left)
and declination (right).

Figure 4:  Offsets of the pulsar timing positions in the INPOP21a frame with
respect to those in the DE440 frame as a function of right ascension (upper)
and declination (lower).

Figure 5:  Offsets of the pulsar timing proper motions in the EPM2021 frame
with respect to those in the DE440 frame as a function of right ascension
(upper) and declination (lower).

Figure 6:  Offsets of the pulsar timing proper motions in the INPOP21a frame
with respect to those in the DE440 frame as a function of right ascension
(upper) and declination (lower).

Figures 3–4 depict the offsets of the timing positions derived using EPM2021
and INPOP21a with respect to those using DE440 at the original timing position
reference epochs for each pulsar. The comparisons of the pulsar timing
positions between other ephemerides and DE440 are shown in Appendix A. The
differences in the right ascension are approximately $\mathrm{100\,\mu as}$
between EPM2021 and DE440 and $\mathrm{200\,\mu as}$ between INPOP21a and
DE440, while the respective offsets in the declination are $\mathrm{50\,\mu
as}$ and $\mathrm{300\,\mu as}$. The weighted root-mean-squares (WRMS) in
right ascension and declination are $\mathrm{118\,\mu as}$ and
$\mathrm{55\,\mu as}$ for EPM2021 versus DE440 and $\mathrm{232\,\mu as}$ and
$\mathrm{305\,\mu as}$ for INPOP21a versus DE440. These positional offsets are
not significant over 1-sigma compared to their uncertainties for approximately
half of the pulsars. When using older planetary ephemerides, the positional
offsets become more pronounced and significant; they are almost all confident
at 1-sigma for DE200 and DE405. In addition, we found obvious patterns in
these plots, especially in the plots of $\Delta\alpha^{*}$ versus $\delta$ and
$\Delta\delta$ versus $\alpha$, which implied clear dependencies of the
positional offsets on the right ascension and declination. The amplitude of
these signatures exceeded the detection sensitivity of the PTA data sets
mentioned in Sect. 2.2, suggesting that the global difference in the position
system between planetary ephemerides could be determined with the timing
solutions.

Similarly, the offsets of the pulsar timing proper motions using EPM2021 and
INPOP21a with respect to those using DE440 are presented in Figs. 5–6, while
the related plots for the remaining ephemerides can be found in Appendix B.
For most cases, we found a dominant proper motion bias in the right ascension,
while the proper motion in the declination was much less affected. Except for
a small fraction of pulsars, the differences in the proper motion in the right
ascension are approximately $\mathrm{5\,\mu as\,yr^{-1}}$ for INPOP21a and
EPM2021; they are only a few $\mathrm{\,\mu as\,yr^{-1}}$ in declination.
These proper motion differences, however, are at least one order-of-magnitude
smaller than their uncertainties and thus not statistically significant. We
observed similar results between other planetary ephemerides and DE440, except
those between DE200 and DE440 (Fig. 31).

We also compared the timing astrometric solutions using the same ephemerides
but based on different PTA data sets. These timing solutions based on the IPTA
DR2B data set were used as the reference solutions. The positional agreements
(WRMS) are between $\mathrm{30\,\mu as}$ and $\mathrm{60\,\mu as}$ in right
ascension and between $\mathrm{60\,\mu as}$ and $\mathrm{180\,\mu as}$ in
declination. For proper motion measurements, the agreements are $\mathrm{\leq
20\,\mu as\,yr^{-1}}$ for right ascension and between $\mathrm{20\,\mu
as\,yr^{-1}}$ and $\mathrm{100\,\mu as\,yr^{-1}}$ for declination.

###  3.2 Modeling the differences

Figure 7:  Rotation parameters around the $X$-, $Y$-, and $Z$-axes in the
position system of the planetary ephemeris frames with respect to DE440 in the
equatorial coordinate system (from top to bottom). $R_{\rm Z}$ between DE200
and DE440 is approximately $-12$ mas, which is beyond the range of the
vertical axis.  Figure 8:  Glide parameters around the $X$-, $Y$-, and
$Z$-axes in the position system of the planetary ephemeris frames with respect
to DE440 in the equatorial coordinate system (from top to bottom).  Figure 9:
Rotation parameters around the $X$-, $Y$-, and $Z$-axes in the proper motion
system of the planetary ephemeris frames with respect to DE440 in the
equatorial coordinate system (from top to bottom).  Figure 10:  Glide
parameters around the $X$-, $Y$-, and $Z$-axes in the proper motion system of
the planetary ephemeris frames with respect to DE440 in the equatorial
coordinate system (from top to bottom).

Figure 11:  WRMS of the postfit residuals as a function of the maximum degree
of the VSH model. Left: EPM2021 versus DE440; Right: INPOP21a versus DE440.
Table 4:  Rotation and glide of the position system of ephemerides with
respect to that of DE440 in the equatorial coordinate system.

Ephemeris | $R_{\rm X}$ | $R_{\rm Y}$ | $R_{\rm Z}$ | $G_{\rm X}$ | $G_{\rm Y}$ | $G_{\rm Z}$
---|---|---|---|---|---|---
| ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$)
DE200 |  $-2291$(22) |  $-12432$(19) |  $-12428$(13) |  $8$(14) |  $-1$(17) |  $-17$(25)
DE403 |  $-1838$(13) |  $597$(13) |  $-2309$(8) |  $10$(9) |  $-17$(12) |  $-6$(16)
DE405 |  $-1623$(13) |  $248$(12) |  $891$(9) |  $6$(9) |  $-8$(11) |  $-5$(16)
DE410 |  $-66$(9) |  $-665$(8) |  $-219$(5) |  $2$(5) |  $-1$(7) |  $-1$(11)
DE421 |  $-33$(7) |  $-227$(7) |  $-184$(4) |  $-2$(3) |  $-3$(5) |  $-0$(8)
DE430 |  $-106$(7) |  $-14$(7) |  $101$(3) |  $1$(3) |  $-1$(6) |  $2$(9)
EPM2011 |  $-113$(8) |  $132$(8) |  $327$(4) |  $-1$(4) |  $-6$(7) |  $2$(9)
EPM2015 |  $-204$(8) |  $97$(7) |  $258$(5) |  $-2$(5) |  $-2$(7) |  $-1$(10)
EPM2017 |  $-192$(8) |  $140$(8) |  $223$(4) |  $-1$(4) |  $-4$(7) |  $0$(10)
EPM2021 |  $-58$(6) |  $8$(6) |  $120$(3) |  $1$(3) |  $-1$(5) |  $2$(7)
INPOP06c |  $3$(11) |  $-275$(10) |  $-308$(7) |  $1$(8) |  $-3$(10) |  $-5$(13)
INPOP08a |  $244$(9) |  $32$(9) |  $-572$(5) |  $2$(4) |  $-6$(7) |  $0$(10)
INPOP10a |  $116$(10) |  $-66$(8) |  $-116$(5) |  $-3$(6) |  $6$(9) |  $-9$(12)
INPOP10b |  $38$(12) |  $388$(10) |  $29$(7) |  $2$(7) |  $-5$(10) |  $1$(16)
INPOP10e |  $-229$(8) |  $-74$(7) |  $-202$(5) |  $0$(5) |  $-5$(6) |  $-1$(10)
INPOP13c |  $-107$(7) |  $-86$(6) |  $-218$(4) |  $-3$(4) |  $1$(6) |  $0$(8)
INPOP17a |  $-176$(9) |  $208$(10) |  $-81$(4) |  $1$(4) |  $-4$(7) |  $1$(10)
INPOP19a |  $-83$(6) |  $246$(6) |  $4$(3) |  $-2$(3) |  $-2$(5) |  $-1$(7)
INPOP21a |  $-237$(4) |  $405$(5) |  $3$(3) |  $-3$(3) |  $-2$(5) |  $2$(6)
888 The uncertainties for the estimates given in parentheses were obtained
from the least-squares fit and could be underestimated several times, as
suggested by the bootstrap fit.

Table 5:  Rotation and glide of the proper motion system of ephemerides with
respect to that of DE440 in the equatorial coordinate system.

Ephemeris | $r_{\rm X}$ | $r_{\rm Y}$ | $r_{\rm Z}$ | $g_{\rm X}$ | $g_{\rm Y}$ | $g_{\rm Z}$
---|---|---|---|---|---|---
| ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$) | ($\mathrm{\mu as\,yr^{-1}}$)
DE200 |  $2.4$(1.6) |  $114.6$(1.5) |  $-259.0$(1.0) |  $-0.7$(1.0) |  $0.4$(1.3) |  $-0.3$(1.9)
DE403 |  $2.1$(1.0) |  $-4.8$(1.0) |  $20.1$(0.6) |  $-1.0$(0.7) |  $1.9$(0.9) |  $-0.7$(1.2)
DE405 |  $-0.7$(1.0) |  $-3.2$(0.9) |  $10.8$(0.7) |  $-0.8$(0.7) |  $1.1$(0.8) |  $-1.1$(1.2)
DE410 |  $-0.2$(0.7) |  $3.5$(0.6) |  $-4.7$(0.4) |  $-0.2$(0.4) |  $0.1$(0.5) |  $-0.1$(0.8)
DE421 |  $-0.2$(0.5) |  $0.3$(0.5) |  $-1.0$(0.3) |  $0.2$(0.3) |  $0.3$(0.4) |  $0.0$(0.6)
DE430 |  $-0.4$(0.6) |  $-3.2$(0.5) |  $6.1$(0.3) |  $-0.1$(0.3) |  $0.1$(0.4) |  $-0.3$(0.6)
EPM2011 |  $-0.1$(0.6) |  $-1.9$(0.6) |  $4.3$(0.3) |  $-0.0$(0.3) |  $0.7$(0.5) |  $-0.6$(0.7)
EPM2015 |  $0.0$(0.6) |  $-3.3$(0.6) |  $6.6$(0.4) |  $0.2$(0.4) |  $0.1$(0.6) |  $0.3$(0.8)
EPM2017 |  $-0.5$(0.6) |  $-5.1$(0.6) |  $11.2$(0.3) |  $0.1$(0.3) |  $0.4$(0.5) |  $-0.1$(0.7)
EPM2021 |  $0.0$(0.5) |  $-1.9$(0.5) |  $2.3$(0.2) |  $0.0$(0.3) |  $0.1$(0.4) |  $-0.2$(0.6)
INPOP06c |  $-0.5$(0.8) |  $-2.7$(0.8) |  $7.0$(0.5) |  $-0.3$(0.6) |  $0.7$(0.7) |  $-0.9$(1.0)
INPOP08a |  $-1.6$(0.7) |  $12.3$(0.6) |  $-30.9$(0.4) |  $-0.1$(0.3) |  $0.5$(0.5) |  $-0.2$(0.8)
INPOP10a |  $0.1$(0.7) |  $-10.0$(0.6) |  $21.5$(0.4) |  $0.3$(0.4) |  $-0.5$(0.6) |  $1.0$(0.9)
INPOP10b |  $0.3$(0.9) |  $-9.0$(0.7) |  $22.5$(0.5) |  $-0.1$(0.5) |  $0.4$(0.7) |  $0.1$(1.1)
INPOP10e |  $-0.1$(0.6) |  $-6.1$(0.5) |  $15.0$(0.4) |  $-0.0$(0.4) |  $0.4$(0.5) |  $0.2$(0.7)
INPOP13c |  $0.4$(0.5) |  $-5.9$(0.5) |  $13.3$(0.3) |  $0.2$(0.3) |  $-0.0$(0.4) |  $0.1$(0.6)
INPOP17a |  $-0.1$(0.7) |  $-3.0$(0.7) |  $6.8$(0.3) |  $-0.0$(0.3) |  $0.3$(0.6) |  $-0.2$(0.8)
INPOP19a |  $0.1$(0.5) |  $6.5$(0.5) |  $-16.0$(0.3) |  $0.2$(0.3) |  $0.1$(0.4) |  $0.1$(0.5)
INPOP21a |  $-0.0$(0.3) |  $1.5$(0.4) |  $-4.6$(0.2) |  $0.2$(0.2) |  $0.2$(0.4) |  $-0.2$(0.5)
999The uncertainties for the estimate given in parentheses were obtained from
the least-squares fit and could be underestimated several times, as suggested
by the bootstrap fit.

We modeled the pulsar timing position and proper motion offsets seen in the
previous section following methods introduced in Sect. 2.4. The least-squares
fits converged after two iterations or less. The WRMS of postfit residuals was
reduced to a few $\mathrm{\mu as}$ for the position and a few $\mathrm{\mu
as\,yr^{-1}}$ for the proper motion. This suggested that the variations in the
timing position and proper motion were well accounted for by our models.

However, we found a few large deviations from the bulk of the sample for a few
pulsars (e.g., PSR J1730–2304 and PSR J1949+3106), for which the precision of
the timing astrometric solutions was relatively low compared to the rest of
the pulsars (approximately 10 mas for position and $\mathrm{5\,mas\,yr^{-1}}$
for proper motions). This implied that the systematics in the planetary
ephemeris reference frames could be evinced in the pulsar timing data set only
when the timing astrometric precision is sufficiently high.

The estimates of the VSH parameters based on various PTA data sets are
depicted in Figs. 7–10, which show good agreement among different PTAs. We
observed a general trend of declining to zero for the successive DE
ephemerides, supporting a better global agreement between newer editions of
the DE ephemerides. Similar results can be found for the successive EPM and
INPOP ephemerides if one takes the results of the latest versions of each
series (i.e., EPM2021 and INPOP21a) as the reference. The positional glide
parameters (components in $\bm{G}$) were estimated to be $\mathrm{10\,\mu as}$
to $\mathrm{20\,\mu as}$ for most cases, obviously smaller than the
corresponding positional rotation parameters (components in $\bm{R}$). This
comment is also valid for the VSH parameters of the proper motion system. We
also noticed that the $X$-components of the rotation vector $\bm{r}$
($\sim\mathrm{1\,\mu as\,yr^{-1}}$) are much less pronounced than the $Y$\-
and $Z$-components.

Table 4 presents estimates of the rotation and glide parameters of the
position system among the planetary ephemerides based on the combination of
the timing solutions from all PTAs. Except for earlier DE planetary
ephemerides (DE200 and DE405), the orientation offsets were generally less
than $\mathrm{500\,\mu as}$, and almost no relative dipolar deformation was
effectively detected.

We tabulated the determinations of the VSH parameters among the planetary
ephemeris proper motion systems in Table 5. The spin around the $X$-axis is
much smaller than those around the $Y$\- and $Z$-axes. We converted the spin
vectors from the equatorial coordinate system to the ecliptic coordinate
system and found that the spins were mainly in the direction of the latitude,
suggesting that the spin was likely due to the differences in the mean orbital
motion of Earth. The glide parameters were less than $\mathrm{1\,\mu
as\,yr^{-1}}$ and were statistically insignificant for most cases.

The VSH fittings were also applied to the subsets of pulsars at high and low
ecliptic latitudes. The results of pulsars near the ecliptic plane showed some
discrepancies compared to the results based on subsets of all pulsars and only
pulsars far from the ecliptic plane, which usually occurred for parameters on
the Z-axis. The most pronounced differences occurred for the comparison
between DE200 and DE440, for which the sample of pulsars at low ecliptic
latitudes yielded a solution of $G_{\rm Z}\\!=\\!-634\,\pm\,375\,{\rm\mu as}$
and $g_{\rm Z}\\!=\\!36\,\pm\,32\,{\rm\mu as\,yr^{-1}}$. These discrepancies,
however, were not statistically significant with respect to the uncertainties
of the estimates for most cases.

The bootstrap fit gave consistent results with those from the least-squares
fit (i.e., reported in Tables 4 and 5), while the formal uncertainty was
usually several times larger than the formal uncertainty of the latter. To
cross-check the results, we processed the NANOGrav12.5 data set using
PINT101010https://pypi.org/project/pint-pulsar/ (version 0.9.1; Luo et al.
2021) developed by the NANOGrav team under the DE planetary ephemeris frames,
from which the results were consistent with those given by TEMPO2.

It was also possible to extend our VSH model to high degrees, for example, to
degree two to include the quadrupolar terms. However, we noticed that the
postfit WRMS in the timing positions and proper motions were already small.
Including higher orders of VSHs might not improve the fit much. The
quadrupolar signals, if they exist, should not exceed $5\,\mathrm{\mu as}$ in
the positions and $1\,\mathrm{\mu as\,yr^{-1}}$ in the proper motions.
Comparisons among timing solutions using DE440, EPM2021, and INPOP21a were
used to test this speculation. The VSH model was extended to different maximum
degrees $l_{\rm max}$ until $l_{\rm max}\\!=\\!5$, to which the differences
between timing solutions were fitted. We found that the estimates of the
rotation and glide parameters varied slightly when $l_{\rm max}\\!\leq\\!4$.
Figure 11 shows the evolution of the WRMS of postfit residuals as $l_{\rm
max}$ increases, which changes little compared to that of $l_{\rm
max}\\!=\\!1$. The results of this experiment suggest that the VSH model of
the first degree is sufficient to account for the global differences in
position and proper motions between planetary ephemerides.

###  3.3 Comparison between timing and VLBI solutions

Figure 12:  Offsets of the pulsar timing positions in the DE200 frame with
respect to those in the VLBI frame as a function of right ascension (left) and
declination (right).

Figure 13:  Offsets of the pulsar timing positions in the DE405 frame with
respect to those in the VLBI frame as a function of right ascension (left) and
declination (right).

Figure 14:  Offsets of the pulsar timing positions in the DE440 frame with
respect to those in the VLBI frame as a function of right ascension (left) and
declination (right).  Table 6:  Orientation offsets of the planetary ephemeris
frames with respect to the VLBI frame.

Ephemeris | $R_{\rm X}$ | $R_{\rm Y}$ | $R_{\rm Z}$
---|---|---|---
| ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$) | ($\mathrm{\mu as}$)
DE200 | $-3141(221)$ | $-11806(216)$ | $-14803(189)$
DE405 | $-2478(219)$ | $-372(210)$ | $1200(187)$
DE410 | $-929(212)$ | $-1233(204)$ | $-152(183)$
DE421 | $-904(218)$ | $-820(211)$ | $4(187)$
DE430 | $-1011(220)$ | $-637(216)$ | $431(189)$
DE440 | $-896(221)$ | $-587(217)$ | $292(189)$
EPM2011 | $-1013(221)$ | $-472(217)$ | $668(190)$
EPM2015 | $-1106(221)$ | $-523(217)$ | $619(189)$
EPM2017 | $-1070(218)$ | $-505(212)$ | $529(188)$
EPM2021 | $-971(221)$ | $-600(216)$ | $412(189)$
INPOP06c | $-893(222)$ | $-882(218)$ | $61(190)$
INPOP08a | $-655(217)$ | $-429(211)$ | $-683(187)$
INPOP10a | $-769(217)$ | $-763(211)$ | $311(187)$
INPOP10b | $-900(221)$ | $-286(216)$ | $534(189)$
INPOP10e | $-1124(223)$ | $-726(217)$ | $242(190)$
INPOP13c | $-1002(221)$ | $-732(217)$ | $210(190)$
INPOP17a | $-1076(221)$ | $-413(217)$ | $281(190)$
INPOP19a | $-978(221)$ | $-276(217)$ | $130(189)$
INPOP21a | $-1170(221)$ | $-164(216)$ | $225(189)$
111111The uncertainties of these estimates given in parentheses were obtained
from the least-square fitting and could be underestimated several times.

The positional offsets between the timing and VLBI solutions are presented in
Figs. 12–14. Only the results of timing solutions using DE200, DE405, and
DE440 are shown here; the results of using other planetary ephemerides are
similar to those of using DE440 and thus not plotted. The positional offsets
between DE200 and VLBI reach $-$20 mas for some pulsars, while they decrease
to approximately 4 mas or less between DE440 and VLBI, with uncertainties of 1
mas to 4 mas in either right ascension or declination. The results based on
different PTA data sets are consistent within their astrometric precision.

The fittings of a frame rotation were performed to the positional offsets
using each PTA data set separately except PPTA DR2e (due to the small sample
size of pulsars in common) as well as their combination. The estimates of
rotation parameters based on the IPTA DR2 data sets always differed by 0.3 mas
to 0.6 mas from those based on the data sets from regional PTAs, although they
still agreed with each other with their uncertainties returned by the least-
squares fitting. These differences indicate possible systematics of the same
level in the input VLBI and/or timing astrometric solutions.

Table 6 reports the fitting results based on the combination of all data sets.
We found an orientation offset of $\mathrm{\sim-3\,mas}$ on the $X$-axis and
of $\mathrm{<\,-10\,mas}$ on the $Y$\- and $Z$-axes in the DE200 frame. We
also detected an orientation offset of $\mathrm{\sim-2.5\,mas}$ on the
$X$-axis for the DE405 frame. The rest of the planetary ephemerides yield a
common orientation offset of $\mathrm{\sim-1.0\,mas}$ on the $X$-axis.

Since there were only two pulsars near the ecliptic plane, we could not use
this subset to test the contribution of pulsars near the ecliptic plane to the
determination of the frame rotation parameters. Instead, we excluded these two
pulsars, leaving only pulsars with high ecliptic latitudes in the sample, and
reran the solution. We found excellent agreements in the estimates of $R_{\rm
X}$ and $R_{\rm Y}$, while there was an offset of $\mathrm{\sim-0.2\,mas}$ in
the estimate of $R_{\rm Z}$ with respect to the solution using all pulsars. We
speculated that this offset was caused by the removal of PSR J2145–0750, which
yielded an offset of $\mathrm{\sim+0.8\,mas}$ in right ascension between
timing and VLBI. An increase in the sample size might lead to a more robust
fit of the frame rotation.

##  4 Discussion

###  4.1 Comparison with previous results

####  4.1.1 Intercomparison among planetary ephemerides

Standish (2004a) compared DE200 and DE405 with DE409 (the newest edition of
the DE series at that time) to approximate the accuracy of DE200 and DE405.
Based on the heliocentric coordinates of the Earth-Moon barycenter, the author
gave the rotation relations between DE200 and DE409 and between DE405 and
DE409 at J2000. The rotation relation can be rewritten as

$\bm{R}_{\rm DE200}=\left[\begin{array}[]{c}-2200\\\ -11700\\\
-12300\end{array}\right]\,{\rm\mu as},\leavevmode\nobreak\ \bm{r}_{\rm
DE200}=\left[\begin{array}[]{c}+3\\\ +108\\\ -249\end{array}\right]\,{\rm\mu
as\,yr^{-1}},$ (24)

for DE200 and

$\bm{R}_{\rm DE405}=\left[\begin{array}[]{c}-1500\\\ +1000\\\
+1100\end{array}\right]\,{\rm\mu as},\leavevmode\nobreak\ \bm{r}_{\rm
DE405}=\left[\begin{array}[]{c}+1\\\ -7\\\ +15\end{array}\right]\,{\rm\mu
as\,yr^{-1}},$ (25)

for DE405 by converting to the same units used in this work. We computed the
relative orientation offsets and spins of DE200 and DE405 with respect to
DE410 by comparing the rotation parameters for DE200 and DE405 with those for
DE410 reported in Tables 4–5. Assuming that DE409 was very close to DE410, we
obtained consistent results with those of Standish (2004a) given in Eqs.
(24)–(25), which validated the effectiveness of our method. The relative spin
found in the DE200 frame mainly reflects the mean motion errors of DE200 for
the Earth-Moon barycenter (EMB), which was caused by the relatively short
time-span coverage of the input ranging data used for deriving DE200 and the
lack of asteroids in modeling their perturbations, as explained in Standish
(2004a). The relative spin between DE405 and DE410 suggested a difference in
the mean motion of Earth, as reported in Standish (2003), which was mainly
caused by the differences in the modeling of perturbation from the largest 20
asteroids and more accurate fitting of Mars’ orbit in DE410 compared to DE405.

Adopting a different computation of the asteroid perturbation compared to
DE405 (considering perturbations of all 300 asteroids on all the planets and
using an asteroid ring), the accuracy obtained with INPOP06 was improved to be
comparable to the accuracy of DE414 (Fienga et al. 2008). We observed in Table
4 that the rotation and glide estimates for INPOP06c were much smaller than
those of DE405 with respect to DE440. The accuracy of the Earth’s orbit of
DE421 was continuously improved by new VLBI tracking observations of
spacecraft (Standish 2004b; Folkner et al. 2007, 2009). In addition, the
asteroid perturbation was modeled as the effects of 67 asteroids with
individual masses and 276 asteroids whose masses were defined by different
density classes. These improvements were reflected in the smaller values of
the rotation and glide parameters between DE421 and DE440 compared to those
between DE410 and DE440. For INPOP08a, the ring model and asteroid selection
were revised. Fienga et al. (2009) showed that the heliocentric longitudes of
the EMB in the ecliptic coordinate computed using INPOP08 differed from those
computed using DE421 or INPOP06 by a general linear drift, resulting in an
ecliptic longitude offset of $\mathrm{-3\,mas}$ over 100 years. It
corresponded roughly to a difference of $\mathrm{-30\,\mu as\,yr^{-1}}$ in the
spin around the $Z$-axis in the ecliptic coordinate, which agreed with the
spin reported in Table 5.

The main differences between INPOP08a and INPOP10a were the input data sets
and the implementation of the fit. Fienga et al. (2011) used timing
astrometric solutions of 18 millisecond pulsars with a precision of better
than 10 mas to estimate the relative orientation offset between DE200, DE405,
DE414, DE421, INPOP08, and INPOP10a. The reported results agreed with those
obtained in this work within the formal uncertainties. Thanks to the
increasing numbers of millisecond pulsars, more timing observations
accumulated over more than 10 years and from more radio telescopes, and
improvements achieved in the timing observing system, our results were several
times more precise than those in Fienga et al. (2011), reaching a precision of
several tens of $\mathrm{\mu as}$. EPM2011 was the first edition of the EPM
series included in this work. Pitjeva (2013) reported that the maximum
differences in the heliocentric positions of Mars between EPM2011 and DE424
were approximately 0.7 mas in right ascension and 0.5 mas in declination. Our
fitting suggested that the frame rotation between EPM2011 and DE421 was
approximately 0.4 mas around the $Y$-axis and 0.5 mas around the $Z$-axis.

No ring model accounting for the perturbation of the Main belt objects was
considered in the INPOP series since INPOP13c. The comparison of the
geocentric coordinates of Mars between INPOP13c and DE430 provided by Fienga
et al. (2014) showed that the differences in right ascension and declination
were all below 0.5 mas. In this work, the agreement in the orientation of each
axis between INPOP13c and DE430 was found to be approximately 0.3 mas or
better, which was thus consistent with the results in Fienga et al. (2014).
The asteroid perturbation models were slightly modified in EPM2015 and EPM2017
compared to EPM2011. We found an orientation agreement of $\sim$0.1 mas among
these planetary ephemeris reference frames. The dynamical modeling of INPOP17a
was almost the same as that of INPOP13c. The differences in the heliocentric
positions of Mars between INPOP17a and DE430 given by the comparison in
Viswanathan et al. (2017) were approximately 0.2 mas in both right ascension
and declination, close to the differences in the rotation parameters between
these two planetary ephemeris reference frames. A new model for the Trans-
Neptunian objects was introduced in INPOP19a in addition to the increase in
the amount of main-belt asteroids in comparison to INPOP17a, leading to a
difference of $\sim$ 100 km in the barycentric position of the EMB that was
mainly in the direction perpendicular to the ecliptic (Fienga et al. 2019). It
seemed that, however, these improvements did not greatly affect the
orientation of the INPOP reference frame. Fienga et al. (2021) compared the
heliocentric position of Earth between DE440 and INPOP21a, which gave an
agreement of 0.2 mas in right ascension and 0.4 mas in declination. These
values agreed roughly with the orientation offset between DE440 and INPOP21a
frames determined from the timing astrometry.

Comparisons with previous results based on other data sets or methods have
validated the effectiveness of our method. Therefore, the transformation
relation between other planetary ephemerides that are not included in this
study might also be derived following the same procedure as described in Sect.
2.4.

####  4.1.2 Comparison between planetary ephemerides and extragalactic frames

DE200 was constructed on its own dynamical equinox of J2000.0 with an accuracy
of 1 arcsecond (Standish 1982). Folkner et al. (1994) estimated the
orientation offset between the DE200 frame and extragalactic frame to be
approximately $\mathrm{-2\,mas}$, $\mathrm{-12\,mas}$, and $\mathrm{-6\,mas}$
around the $X$-, $Y$-, and $Z$-axes in 1988 by a joint analysis of the lunar
laser ranging and VLBI measurements. These values roughly match our
determination of the frame rotation between the DE200 frame and the VLBI frame
at J2000.0, plus the propagation of the frame spin.

Since DE403, the tie between all the planetary ephemeris reference frames and
the ICRF is achieved through the adjustment of the planetary ephemerides to
the VLBI observations of the spacecraft relative to the extragalactic sources.
For the orbits of inner planets, this frame tie is dominated by VLBI
measurements of Mars-orbiting spacecraft. Therefore, the accuracy of the frame
tie between the planetary ephemeris reference frames and the ICRF is on the
same level as the accuracy of these VLBI measurements.

The DE405 frame was tied to ICRF1 (Ma et al. 1998) with an uncertainty of
$\sim\mathrm{1\,mas}$ (Standish 1998) based on VLBI observations of the
Magellan spacecraft in orbit around Venus. Later, in DE410, the frame tie was
improved by VLBI points of Mars-orbiting spacecraft (Standish 2003). The
improvement of frame-tie from DE405 to DE410 can be easily seen in Table 6,
although our estimations suggest larger orientation offsets between these two
ephemeris frames and the VLBI frame.

The DE421 frame was tied to the ICRF1 with an accuracy better than
$\mathrm{1\,mas}$ (Folkner et al. 2009). Fienga et al. (2013) reported that
the link between the INPOP10e reference frame and the ICRF was maintained with
an accuracy of $\mathrm{\sim 1\,mas}$ over 10 years. Pitjeva (2013) determined
the orientation offset of the EPM2011 frame with respect to the ICRF1, which
depended strongly on the time span of the input VLBI observations of various
spacecraft. Nevertheless, the frame tie between EPM2011 and ICRF1 is most
likely to also be better than $\mathrm{1\,mas}$. Our comparisons between the
timing and VLBI solutions agree with these results within the quoted
uncertainty.

The DE430 frame was aligned with ICRF2 (Fey et al. 2015) with an accuracy of
$\mathrm{0.2\,mas}$ (Folkner et al. 2014; Folkner & Border 2015). Park et al.
(2015) also determined the orientation of the DE430 frame relative to the ICRF
by using VLBA measurements of Mars with an uncertainty of 0.23 mas. The link
accuracy of the INPOP13c frame to ICRF2 was improved to be better than
$\mathrm{0.5\,mas}$ compared to INPOP10e (Fienga et al. 2014). The EPM2015 and
EPM2017 frames were oriented to the ICRF2 with an accuracy better than 0.2 mas
at 3-sigma (Pitjeva 2017; Pitjeva & Pitjev 2018). A recent study presented by
Deram et al. (2022) used the observations of asteroids in the Gaia Data
Release 2 (Gaia Collaboration et al. 2018a) to determine the frame rotation
between the INPOP19a reference frame and Gaia-CRF2 (Gaia Collaboration et al.
2018b). The misalignment between these two reference frames was found to be
approximately $\mathrm{0.25\,mas}$; it decreased to a few $\mathrm{\mu as}$
when only very accurate observations of inner planets were considered.
Considering that the agreement between Gaia-CRF2 and ICRF3 is on the level of
$\mathrm{20\,\mu as}$–$\mathrm{30\,\mu as}$ (Charlot et al. 2020; Gaia
Collaboration et al. 2018b; Liu et al. 2020), the tie between the INPOP19a
frame and ICRF3 should be accurate at $\mathrm{0.25\,mas}$ or much better. The
average accuracy of the orientation of the inner planet orbits for DE440 tied
to ICRF3 is approximately 0.2 mas (Park et al. 2021). For these planetary
ephemeris frames, however, we detected a common frame rotation $R_{\rm
X}\\!\simeq\\!\mathrm{-1\,mas}$ that is confident at 4-sigma or more; this
frame rotation would decrease if only an individual PTA data set was used.
Therefore, we suspected that the common rotation was likely related to the
systematics in the input timing data sets.

Fienga et al. (2011) found a significant rotation of $\mathrm{\sim 10\,mas}$
in the planetary ephemerides with respect to the VLBI frame based on four
pulsars with both VLBI and timing astrometric measurements, which diminishes
in our results with improved timing and VLBI astrometric data. Therefore, the
rotation reported in Fienga et al. (2011) is likely caused by the unmodeled
systematics in the input timing or VLBI positions used therein. Wang et al.
(2017) studied the rotation between VLBI positions and timing positions using
different ephemerides (DE405, DE414, DE421, DE430, DE432, and DE435) for five
millisecond pulsars and reported a statistically significant rotation angle of
$\mathrm{\sim 2\,mas}$ around the $X$-axis for DE405. Similar results were
obtained in our comparison (see Table 6).

The orientation offsets of the DE440, EPM2021, and INPOP21a frames with
respect to the VLBI frame range from $\mathrm{0.2\,mas}$ to
$\mathrm{1.2\,mas}$ in an absolute sense. Taking the formal uncertainty of
$\mathrm{\sim 0.2\,mas}$ of the rotation parameters and the possible
systematics of $\mathrm{0.3\,mas}$–$\mathrm{0.6\,mas}$ (Sect. 3.3) into
consideration, the nonzero orientation offset of the ephemeride frame is
approximately $\mathrm{0.4\,mas}$. This deviation is similar to the accuracy
of the input VLBI measurements used for constructing the planetary
ephemerides.

We noted that the formal uncertainty of the frame rotation based on pulsar
astrometry was close to those based on the VLBA measurements, even though we
used only a small fraction of the sample from the MSPSR$\pi$ campaign. More
precise frame-tie between the planetary ephemeris frame and the VLBI frame via
pulsars can be expected with the increased sample size and reduced systematics
in the absolute VLBI positions in the future data release of this campaign.

###  4.2 Accuracy of state-of-the-art planetary ephemeris reference frames

The orientation offset and global spin of the ephemerides generally decrease
with respect to the latest ephemerides in each series (i.e., DE440, EPM2021,
and INPOP21a). These results can be expected and validate the improvements in
the error control for the successive ephemerides in the DE, EPM, and INPOP
series from the aspect of the pulsar timing astrometry.

According to Standish (2004a), there are three main factors determining the
accuracy of an ephemeris, including (i) the correctness and completeness of
the equations of motion, (ii) the method and algorithm of integration, and
(iii) the input observational data to create the ephemeris. The third factor
contributed mainly to the improvements of the ephemerides in the DE series
since DE405. Here, we use the INPOP series as an example to discuss the effect
of these factors on our results, for which INPOP21a is used as the reference
ephemeris. There is no modification in the dynamical modeling of INPOP10a
compared with INPOP08, but new observations, especially the normal points
deduced from flybys of the Messenger spacecraft around Mercury, were used to
create INPOP10a (Fienga et al. 2011). This may explain the improvements of
$\mathrm{0.4\,mas}$ in the orientation of the INPOP10a frame compared with the
INPOP08a frame. INPOP10a and INPOP10b were fitted on the same data sample,
while the modeling of the asteroid perturbations on planet orbits for INPOP10b
was improved and determined with an advanced algorithm (Fienga et al. 2012).
Our comparison suggests that the orientation offset of the INPOP10b frame was
reduced by $\mathrm{0.4\,mas}$ in the $Y$-axis compared to the INPOP10a frame,
which can be considered as a contribution from the first two factors.

The orientation agreement was found to be better than 0.15 mas between DE440
and EPM2021 frames and 0.40 mas between DE440 and INPOP21a frames (Sect. 3.1).
In addition, the discrepancy of orientation between the ephemerides frames and
VLBI frame was estimated to be $\mathrm{\sim 0.4\,mas}$ (Sect. 4.1.2). From
these results, we concluded that the orientation of the state-of-art planetary
ephemeris frames is accurate at 0.4 mas. The relative spin is approximately
$\mathrm{2\,\mu as\,yr^{-1}}$ between DE440 and EPM2021 and $\mathrm{5\,\mu
as\,yr^{-1}}$ between DE440 and INPOP21a, leading to the conclusion that the
nonrotating of the planetary frame is no greater than $\mathrm{5\,\mu
as\,yr^{-1}}$. In addition, there is no significant dipolar deformation in the
planetary frames at the level of $\mathrm{1\,\mu as}$ in the position and
$\mathrm{0.1\,\mu as\,yr^{-1}}$ in the proper motion.

##  5 Conclusions

In this work, we used the up-to-date pulsar timing observations from the data
releases of several pulsar timing arrays and astrometric VLBI measurements for
pulsars to investigate the systematics in the planetary ephemeris reference
frames. We reported the rotation and glide in the position and proper motion
systems of successive editions of the DE series, EPM series, and INPOP series
with respect to those of DE440, which is as precise as several $\mathrm{\mu
as}$ in position and $\mathrm{0.1\,\mu as\,yr^{-1}}$ in proper motion.
Depending on the sky location, the variations in the pulsar timing position
and proper motion due to the different choices are on the same order of
magnitude as the derived rotation and glide parameters. The orientation offset
between the ephemerides frames and VLBI frame is estimated to be approximately
$\mathrm{0.4\,mas}$ when considering the formal uncertainty. Our
determinations of rotation parameters are improved compared to previous
studies thanks to the precise pulsar timing observations of more than 60
millisecond pulsars, which were not available before. Based on the results
mentioned above, we conclude that the misalignment of the latest planetary
reference frames is no greater than 0.4 mas with a nonrotating of better than
$\mathrm{5\,\mu as\,yr^{-1}}$.

###### Acknowledgements.

We sincerely thank the referee, Dr. Agnès Fienga, for the constructive
comments and useful suggestions, which have greatly improve the work. We wish
to thank Prof. Xue-Mei Deng for her critical reading and comments on the first
draft. We also thank Dr. Jing Luo for his kind instructions on PINT. N. Liu
and Z. Zhu were supported by the National Natural Science Foundation of China
(NSFC) under grants Nos 11833004 and 12103026. N. Liu was also supported by
the China Postdoctoral Science Foundation (Grant Number: 2021M691530) and the
Yuxiu Postdoctoral Institute at Nanjing University (“Yuxiu Young Scholars
Program”). J. Antoniadis was supported by the Stavros Niarchos Foundation
(SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under
the 2nd Call of “Science and Society” Action Always strive for excellence –
“Theodoros Papazoglou” (Project Number: 01431). This research also made use of
Numpy (van der Walt et al. 2011), IPython (Perez & Granger 2007),
Astropy121212http://www.astropy.org (Astropy Collaboration et al. 2018), the
Python 2D plotting library Matplotlib (Hunter 2007), and NASA’s Astrophysics
Data System. All the results shown in this work are also presented in a series
of Python Jupyter Notebooks posted and publicly accessible
online131313https://git.nju.edu.cn/neo/pulsar-ephem-crf.

## References

  * Alam et al. (2021a) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021a, ApJS, 252, 4
  * Alam et al. (2021b) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2021b, ApJS, 252, 5
  * Antoniadis (2021) Antoniadis, J. 2021, MNRAS, 501, 1116
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bartel et al. (1996) Bartel, N., Chandler, J. F., Ratner, M. I., et al. 1996, AJ, 112, 1690
  * Bernus et al. (2019) Bernus, L., Minazzoli, O., Fienga, A., et al. 2019, Phys. Rev. Lett., 123, 161103
  * Bernus et al. (2020) Bernus, L., Minazzoli, O., Fienga, A., et al. 2020, Phys. Rev. D, 102, 021501
  * Bernus et al. (2022) Bernus, L., Minazzoli, O., Fienga, A., et al. 2022, Phys. Rev. D, 105, 044057
  * Burke-Spolaor et al. (2019) Burke-Spolaor, S., Taylor, S. R., Charisi, M., et al. 2019, A&A Rev., 27, 5
  * Caballero et al. (2018) Caballero, R. N., Guo, Y. J., Lee, K. J., et al. 2018, MNRAS, 481, 5501
  * Charlot et al. (2020) Charlot, P., Jacobs, C. S., Gordon, D., et al. 2020, A&A, 644, A159
  * Chatterjee et al. (2009) Chatterjee, S., Brisken, W. F., Vlemmings, W. H. T., et al. 2009, ApJ, 698, 250
  * Chen et al. (2021) Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970
  * Deller et al. (2019) Deller, A. T., Goss, W. M., Brisken, W. F., et al. 2019, ApJ, 875, 100
  * Deller et al. (2009) Deller, A. T., Tingay, S. J., Bailes, M., & Reynolds, J. E. 2009, ApJ, 701, 1243
  * Deng et al. (2014) Deng, X.-M., Fan, M., & Xie, Y. 2014, Chinese Astron. Astrophys., 38, 330
  * Deram et al. (2022) Deram, P., Fienga, A., Verma, A. K., Gastineau, M., & Laskar, J. 2022, Celestial Mechanics and Dynamical Astronomy, 134, 32
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341
  * Ding et al. (2020) Ding, H., Deller, A. T., Freire, P., et al. 2020, ApJ, 896, 85
  * Feissel & Mignard (1998) Feissel, M. & Mignard, F. 1998, A&A, 331, L33
  * Fey et al. (2015) Fey, A. L., Gordon, D., Jacobs, C. S., et al. 2015, AJ, 150, 58
  * Fienga et al. (2020a) Fienga, A., Avdellidou, C., & Hanuš, J. 2020a, MNRAS, 492, 589
  * Fienga et al. (2022) Fienga, A., Bigot, L., Mary, D., et al. 2022, in Multi-Scale (Time and Mass) Dynamics of Space Objects, ed. A. Celletti, C. Galeş, C. Beaugé, & A. Lemaître, Vol. 364, 31–51
  * Fienga et al. (2021) Fienga, A., Deram, P., Di Ruscio, A., et al. 2021, Notes Scientifiques et Techniques de l’Institut de Mecanique Celeste, 110
  * Fienga et al. (2019) Fienga, A., Deram, P., Viswanathan, V., et al. 2019, Notes Scientifiques et Techniques de l’Institut de Mecanique Celeste, 109
  * Fienga et al. (2020b) Fienga, A., Di Ruscio, A., Bernus, L., et al. 2020b, A&A, 640, A6
  * Fienga et al. (2011) Fienga, A., Laskar, J., Kuchynka, P., et al. 2011, Celestial Mechanics and Dynamical Astronomy, 111, 363
  * Fienga et al. (2009) Fienga, A., Laskar, J., Morley, T., et al. 2009, A&A, 507, 1675
  * Fienga et al. (2012) Fienga, A., Laskar, J., Verma, A., et al. 2012, Complementary version of INPOP planetary ephemerides, INPOP10b, online document available at https://www.imcce.fr/content/medias/recherche/equipes/asd/inpop/inpop10b.pdf
  * Fienga et al. (2008) Fienga, A., Manche, H., Laskar, J., & Gastineau, M. 2008, A&A, 477, 315
  * Fienga et al. (2013) Fienga, A., Manche, H., Laskar, J., Gastineau, M., & Verma, A. 2013, arXiv e-prints, arXiv:1301.1510
  * Fienga et al. (2014) Fienga, A., Manche, H., Laskar, J., Gastineau, M., & Verma, A. 2014, arXiv e-prints, arXiv:1405.0484
  * Folkner et al. (1994a) Folkner, W., Border, J., & Iijima, B. 1994a, JPL Interoffice memorandum, 335, 94
  * Folkner et al. (1993) Folkner, W., Kroger, P., & Iijima, B. 1993, JPL Interoffice memorandum, 335, 1
  * Folkner et al. (1994b) Folkner, W., Kroger, P., Iijima, B., & Border, J. 1994b, JPL Interoffice memorandum, 335, 94
  * Folkner & Border (2015) Folkner, W. M. & Border, J. S. 2015, Highlights of Astronomy, 16, 219
  * Folkner et al. (1994) Folkner, W. M., Charlot, P., Finger, M. H., et al. 1994, A&A, 287, 279
  * Folkner et al. (2007) Folkner, W. M., Standish, E. M., J., Williams, J. G., & Boggs, D. H. 2007, JPL Interoffice memorandum, 343R-07-005, 1
  * Folkner et al. (2009) Folkner, W. M., Williams, J. G., & Boggs, D. H. 2009, Interplanetary Network Progress Report, 42-178, 1
  * Folkner et al. (2014) Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., & Kuchynka, P. 2014, Interplanetary Network Progress Report, 42-196, 1
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1
  * Gaia Collaboration et al. (2022a) Gaia Collaboration, Klioner, S. A., Lindegren, L., et al. 2022a, A&A, 667, A148
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Klioner, S. A., Mignard, F., et al. 2021, A&A, 649, A9
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Mignard, F., Klioner, S. A., et al. 2018b, A&A, 616, A14
  * Gaia Collaboration et al. (2022b) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022b, arXiv e-prints, arXiv:2208.00211
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
  * Hotan et al. (2006) Hotan, A. W., Bailes, M., & Ord, S. M. 2006, MNRAS, 369, 1502
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Jones et al. (2021) Jones, D., Folkner, W., Park, R., & Dhawan, V. 2021, in American Astronomical Society Meeting Abstracts, Vol. 53, American Astronomical Society Meeting Abstracts, 142.05
  * Jones et al. (2015) Jones, D. L., Folkner, W. M., Jacobson, R. A., et al. 2015, AJ, 149, 28
  * Jones et al. (2020) Jones, D. L., Folkner, W. M., Jacobson, R. A., et al. 2020, AJ, 159, 72
  * Jones et al. (2011) Jones, D. L., Fomalont, E., Dhawan, V., et al. 2011, AJ, 141, 29
  * Kaspi et al. (1994) Kaspi, V. M., Taylor, J. H., & Ryba, M. F. 1994, ApJ, 428, 713
  * Liu & Capitaine (2017) Liu, J. C. & Capitaine, N. 2017, A&A, 597, A83
  * Liu et al. (2020) Liu, N., Lambert, S. B., Zhu, Z., & Liu, J. C. 2020, A&A, 634, A28
  * Liu et al. (2023) Liu, N., Zhu, Z., Antoniadis, J., et al. 2023, A&A, 670, A173
  * Luo et al. (2021) Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45
  * Ma et al. (1998) Ma, C., Arias, E. F., Eubanks, T. M., et al. 1998, AJ, 116, 516
  * Matthews et al. (2016) Matthews, A. M., Nice, D. J., Fonseca, E., et al. 2016, ApJ, 818, 92
  * Mignard & Klioner (2012) Mignard, F. & Klioner, S. 2012, A&A, 547, A59
  * Park et al. (2015) Park, R. S., Folkner, W. M., Jones, D. L., et al. 2015, AJ, 150, 121
  * Park et al. (2021) Park, R. S., Folkner, W. M., Williams, J. G., & Boggs, D. H. 2021, AJ, 161, 105
  * Perera et al. (2019) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS, 490, 4666
  * Perez & Granger (2007) Perez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21
  * Pitjeva et al. (2022) Pitjeva, E., Pavlov, D., Aksim, D., & Kan, M. 2022, in Multi-Scale (Time and Mass) Dynamics of Space Objects, ed. A. Celletti, C. Galeş, C. Beaugé, & A. Lemaître, Vol. 364, 220–225
  * Pitjeva (2005) Pitjeva, E. V. 2005, in IAU Colloq. 196: Transits of Venus: New Views of the Solar System and Galaxy, ed. D. W. Kurtz, 230–241
  * Pitjeva (2013) Pitjeva, E. V. 2013, Solar System Research, 47, 386
  * Pitjeva (2017) Pitjeva, E. V. 2017, Transactions of IAA RAS, 10
  * Pitjeva & Pitjev (2013) Pitjeva, E. V. & Pitjev, N. P. 2013, MNRAS, 432, 3431
  * Pitjeva & Pitjev (2018) Pitjeva, E. V. & Pitjev, N. P. 2018, Celestial Mechanics and Dynamical Astronomy, 130, 57
  * Pitjeva et al. (2021) Pitjeva, E. V., Pitjev, N. P., Pavlov, D. A., & Turygin, C. C. 2021, A&A, 647, A141
  * Reardon et al. (2016) Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751
  * Reardon et al. (2021) Reardon, D. J., Shannon, R. M., Cameron, A. D., et al. 2021, MNRAS, 507, 2137
  * Splaver et al. (2005) Splaver, E. M., Nice, D. J., Stairs, I. H., Lommen, A. N., & Backer, D. C. 2005, ApJ, 620, 405
  * Standish (2004a) Standish, E. M. 2004a, A&A, 417, 1165
  * Standish & Fienga (2002) Standish, E. M. & Fienga, A. 2002, A&A, 384, 322
  * Standish (1982) Standish, E. M., J. 1982, A&A, 114, 297
  * Standish (1990a) Standish, E. M., J. 1990a, A&A, 233, 272
  * Standish (1990b) Standish, E. M., J. 1990b, A&A, 233, 252
  * Standish (1998) Standish, E. M., J. 1998, JPL Interoffice memorandum, 312.F-98-048, 1
  * Standish (2003) Standish, E. M., J. 2003, JPL Interoffice memorandum, 312.N-03-009, 1
  * Standish (2004b) Standish, E. M., J. 2004b, JPL Interoffice memorandum, 343R-06-002, 1
  * Tkachuk & Choliy (2013) Tkachuk, V. V. & Choliy, V. Y. 2013, Advances in Astronomy and Space Physics, 3, 141
  * Vallisneri et al. (2020) Vallisneri, M., Taylor, S. R., Simon, J., et al. 2020, ApJ, 893, 112
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Verbiest et al. (2008) Verbiest, J. P. W., Bailes, M., van Straten, W., et al. 2008, ApJ, 679, 675
  * Verma et al. (2014) Verma, A. K., Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2014, A&A, 561, A115
  * Vigeland et al. (2018) Vigeland, S. J., Deller, A. T., Kaplan, D. L., et al. 2018, ApJ, 855, 122
  * Viswanathan et al. (2017) Viswanathan, V., Fienga, A., Gastineau, M., & Laskar, J. 2017, Notes Scientifiques et Techniques de l’Institut de Mecanique Celeste, 108
  * Wang et al. (2017) Wang, J. B., Coles, W. A., Hobbs, G., et al. 2017, MNRAS, 469, 425
  * Yang et al. (2022) Yang, P., Huang, Y., Li, P., et al. 2022, Advances in Space Research, 69, 1060
  * Yuan et al. (2021) Yuan, Y., Li, F., Fu, Y., & Chen, J. 2021, A&A, 654, A66

##  Appendix A Pulsar timing position differences using different ephemerides

Figures 15–30 presents the offsets of the pulsar position from the timing
analyses using successive versions of planetary ephemerides from the DE
series, EPM series, and INPOP series with respect to those from the timing
solution using DE440.

Figure 15:  Offsets of the pulsar timing positions in the DE200 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 16:  Offsets of the pulsar timing positions in the DE405 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 17:  Offsets of the pulsar timing positions in the DE410 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 18:  Offsets of the pulsar timing positions in the DE421 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 19:  Offsets of the pulsar timing positions in the DE430 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 20:  Offsets of the pulsar timing positions in the EPM2011 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 21:  Offsets of the pulsar timing positions in the EPM2015 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 22:  Offsets of the pulsar timing positions in the EPM2017 frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 23:  Offsets of the pulsar timing positions in the INPOP06c frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 24:  Offsets of the pulsar timing positions in the INPOP08a frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 25:  Offsets of the pulsar timing positions in the INPOP10a frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 26:  Offsets of the pulsar timing positions in the INPOP10b frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 27:  Offsets of the pulsar timing positions in the INPOP10e frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 28:  Offsets of the pulsar timing positions in the INPOP13c frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 29:  Offsets of the pulsar timing positions in the INPOP17a frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 30:  Offsets of the pulsar timing positions in the INPOP19a frame with
referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

##  Appendix B Pulsar timing proper motion differences using different
ephemerides

Figures 31–46 presents the offsets of the pulsar proper motion from the timing
solutions using successive versions of planetary ephemerides from the DE
series, EPM series, and INPOP series with respect to those from the timing
solution using DE440.

Figure 31:  Offsets of the pulsar timing proper motions in the DE200 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 32:  Offsets of the pulsar timing proper motions in the DE405 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 33:  Offsets of the pulsar timing proper motions in the DE410 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 34:  Offsets of the pulsar timing proper motions in the DE421 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 35:  Offsets of the pulsar timing proper motions in the DE430 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 36:  Offsets of the pulsar timing proper motions in the EPM2011 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 37:  Offsets of the pulsar timing proper motions in the EPM2015 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 38:  Offsets of the pulsar timing proper motions in the EPM2017 frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 39:  Offsets of the pulsar timing proper motions in the INPOP06c frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 40:  Offsets of the pulsar timing proper motions in the INPOP08a frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 41:  Offsets of the pulsar timing proper motions in the INPOP10a frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 42:  Offsets of the pulsar timing proper motions in the INPOP10b frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 43:  Offsets of the pulsar timing proper motions in the INPOP10e frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 44:  Offsets of the pulsar timing proper motions in the INPOP13c frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 45:  Offsets of the pulsar timing proper motions in the INPOP17a frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).

Figure 46:  Offsets of the pulsar timing proper motions in the INPOP19a frame
with referred to those in the DE440 frame as a function of the right ascension
(left) and declination (right).