Gaia coordinates in J2000 as done by the CDS

Question

So I'm interested in knowing how the CDS calculates the J2000 equatorial coordinates from the Gaia DR3 equatorial coordinates (that are given in epoch 2016.0).

For example, for the star HD 308819 = Gaia DR3 5333580380647188736, the Gaia coordinates are (ra, dec) = 174.50345338328, -63.37448911466 (deg, deg) in J2016. The CDS calculates for Vizier the corresponding J2000 coordinates as (RAJ2000, DEJ2000) = 11:38:00.8430618017 -63:22:28.166022064.

I tried to do this in Python by using astropy like so:

import astropy.units as u
from astropy.coordinates import SkyCoord

# Gaia coordinates (ICRS | J2016):
ra = 174.50345338328
dec = -63.37448911466
gaia_coords = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame = 'icrs', equinox = 'J2016.0')

# Transform coordinates into J2000
j2000_coords = gaia_coords.transform_to('fk5')
raj2000 = j2000_coords.ra.to_string(unit = u.hourangle, sep = ':', precision = 10, alwayssign = False, pad = True)
dej2000 = j2000_coords.dec.to_string(unit = u.degree, sep = ':', precision = 9, alwayssign = True, pad = True)
print(raj2000+' '+dej2000)

But I get 11:38:46.0096374246 -63:27:47.407336922 instead of the J2000 coordinates calculated by the CDS. There is a difference of almost an arcminute.

Reading carefully, we see that Vizier states that proper motions [are] taken into account. So, I not only have to transform between J2016 and J2000 but also backpropagate the space motion of the star 16 years back. I can do that with Python like so:

import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord

# Gaia coordinates and proper motion (ICRS | J2016):
ra = 174.50345338328
dec = -63.37448911466
pmra = -5.987
pmdec = 0.326
gaia_coords = SkyCoord(ra = ra*u.deg, dec = dec*u.deg,
                       pm_ra_cosdec = pmra*u.mas/u.yr, pm_dec = pmdec*u.mas/u.yr,
                       frame = 'icrs', equinox = 'J2016.0')

# Transform coordinates into J2000 accounting for proper motion
delta_t = Time('J2000', format = 'jyear_str')-Time('J2016.0', format = 'jyear_str')
j2000_coords = gaia_coords.apply_space_motion(dt = delta_t)
raj2000 = j2000_coords.ra.to_string(unit = u.hourangle, sep = ':', precision = 10, alwayssign = False, pad = True)
dej2000 = j2000_coords.dec.to_string(unit = u.degree, sep = ':', precision = 9, alwayssign = True, pad = True)
print(raj2000+' '+dej2000)

But the result is 11:38:00.8430610377 -63:22:28.166028471, which is way closer the values given by the CDS. But still they differ in the microarcsecond. This might seem stupid but I actually need to understand those extra decimals.

Is it, perhaps, that Python is calculating with larger precision than the CDS (or viceversa)? Is it because the frame shouldn't be FK5? What is the issue here?

Answer 1

Parallax, proper motion and radial velocity should be taken into account. You can use the function EPOCH_PROP_POS for epoch propagation as follow:

from astroquery.gaia import Gaia

query = """
SELECT source_id, ra, dec,
COORD1(EPOCH_PROP_POS(ra, dec, parallax, pmra, pmdec, radial_velocity, ref_epoch, 2000)) AS ra2000,
COORD2(EPOCH_PROP_POS(ra, dec, parallax, pmra, pmdec, radial_velocity, ref_epoch, 2000)) AS dec2000
FROM gaiadr3.gaia_source
WHERE source_id = '5333580380647188736'
"""

job = Gaia.launch_job(query)

table = job.get_results()

print(table[['ra', 'dec', 'ra2000', 'dec2000']])

Output:

        ra                dec               ra2000           dec2000      
       deg                deg                                             
------------------ ------------------ ----------------- ------------------
174.50345338327568 -63.37448911466239 174.5035127575071 -63.37449056168441

If you want to do the conversion on your computer, use erfa:

import erfa
import astropy.units as u
from astropy.table import Table
from astropy.time import Time
from astroquery.gaia import Gaia


source_id = 5333580380647188736

job = Gaia.launch_job(
f"""
SELECT ra, dec, pmra, pmdec, parallax, radial_velocity
FROM gaiadr3.gaia_source_lite
WHERE source_id = '{source_id}'
"""
)

tbl = job.get_results()

# In case there is no radial_velocity
tbl['radial_velocity'].fill_value = 0
tbl = tbl.filled()


ra1 = tbl['ra'].quantity.to('rad').value[0]
dec1 = tbl['dec'].quantity.to('rad').value[0]
pmr1 = tbl['pmra'].quantity.to('rad/yr').value[0]
pmd1 = tbl['pmdec'].quantity.to('rad/yr').value[0]
px1 = tbl['parallax'].quantity.to('arcsec').value[0]
rv1 = tbl['radial_velocity'].quantity.to('km/s').value[0]

t1 = Time('J2016')
t2 = Time('J2000')

# Use starpm or pmsafe
ra2, dec2, pmr2, pmd2, px2, rv2 = erfa.pmsafe(ra1, dec1, pmr1, pmd1, px1, rv1, t1.jd1, t1.jd2, t2.jd1, t2.jd2)


print('RA_2016  :', (ra1*u.rad).to('deg'))
print('RA_2000  :', (ra2*u.rad).to('deg'))
print('DEC_2016  :', (dec1*u.rad).to('deg'))
print('DEC_2000  :', (dec2*u.rad).to('deg'))

Output:

RA_2016  : 174.50345338327568 deg
RA_2000  : 174.50347999226292 deg
DEC_2016  : -63.37448911466238 deg
DEC_2000  : -63.37449056169426 deg