Mitigating astrometric bias in barycentric correction with PEXO

Extremely precise radial velocity is essential for the detection of sub-m/s radial velocity of stars induced by Earth-like planets. Although modeling of the barycentric correction of radial velocity could achieve 1 mm/s precision, the input astrometry could be biased due to nonlinear motions of stars caused by companions. To account for astrometry-induced bias in barycentric correction, we correct for astrometric bias by minimizing the scatter of reduced RV data with PEXO. In particular, we apply this method to the barycentric correction for 266 stars from HARPS data archive. We find that the RV scatter for 8 targets are significantly reduced due to correction of astrometric bias. Among these targets, 2 targets exhibit bias caused by known massive companions, while for the remaining 6 targets, the bias could be attributed to unknown companions or Gaia systematics. Furthermore, 14 targets have an astrometry-induced annual RV variation higher than 0.05 m/s, and 10 of them are closer than 10 pc. We show the results of Barnard's star as an example, and find that an annual RV bias of 10 cm/s is mitigated by replacing BarCor by PEXO as the barycentric correction code. Our work demonstrates the necessity of astrometric bias correction and the utilization of barycentric correction code within a relativistic framework in high-precision RV for the detection of Earth-like planets.


INTRODUCTION
Inspired by the Large UV/Optical/Infrared Surveyor (LUVOIR) (Team 2019) and the Habitable Exoplanet Observatory (HabEx) (Gaudi et al. 2020), the Habitable Worlds Observatory (HWO) is proposed to search for life on nearby Earth-like exoplanets (Mamajek & Stapelfeldt 2023).One essential task to archive this goal is to detect several nearby Earth-like planets.
Extremely precise radial velocity (EPRV) is important for the detection of small radial velocity (RV) signals caused by small planets.To detect the RV variation induced by an Earth-like planet (Earth-size planets in habitable zone), the noise caused by host star, instrument and data reduction pipeline should be reduced to a level of ∼0.1 m s −1 .While RV facilities such as Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) mounted on the European Southern Observatory's Very Large Telescope (VLT) (Pepe et al. 2021) can achieve an instrumental prediction of about 0.1 m s −1 , the data reduction pipeline and stellar activity introduce the main challenges to further improvement of long-term RV precision for the detection of Earth-like planets.
Many sources of stellar noise should be mitigated in order to reliably detecting Earth-like planets.For example, star spots and plages would induce RV signals of > 1 m s −1 (Saar & Donahue 1997); Star convection have the effect of ∼ 500 m s −1 for solar-type stars (Dravins et al. 1986); Effects of micro-telluric lines on precise radial velocities would lead to a noise of > 0.1 m s −1 (Cunha et al. 2014).Advanced noise modeling of time-correlated noise (e.g.Feng et al. 2016a;Rajpaul 2017), line-by-line analysis by Dumusque 2018 andLisogorskyi et al. 2019 as well as high-cadence RV measurements (Hall et al. 2018) are proposed to disentangle planetary signals from stellar noises.
Data reduction puts another challenge to EPRV.One of the main components in RV data reduction is barycentric correction, which subtracts the motion of the Earth from the measured RVs of a star to derive the RV variation of the star viewed from the heliocenter.This so-called "barycentric correction" has been well modeled to 1 mm s −1 precision in recently years (Wright & Eastman 2014;Feng et al. 2019a).However, these models require precise input astrometry, precise mid-exposure time, and solar system ephemerides.The first requirement is not satisfied when a star move nonlinearly due to companion perturbations, while the other two requirements could be easily satisfied by using color-dependent high-precision photon counters (Pepe et al. 2021;Jurgenson et al. 2016) and well updated JPL ephemerides (Folkner et al. 2014).
The nonlinear motion of a star (or reflex motion) caused by unknown companions would vary the line-of-sight direction and the projection of Earth's velocity along this direction.Because both the nonlinear stellar motion and the Earth's annual motion would change the line of sight, the biased astrometry inferred under the assumption of linear stellar motion (e.g.Gaia Collaboration et al. 2016Collaboration et al. , 2018Collaboration et al. , 2021Collaboration et al. , 2022a) ) would not only introduce an RV variation with a period of the orbital motion of the unknown companion but also induce an annual RV variation caused by the coupling of the reflex motion and the Earth's motion.Considering that this coupled RV variation is quite similar to the RV signal induced by an Earth-like planet in terms of periodicity (300-400 days) and amplitude (∼0.1 m s −1 ), it is crucial to mitigate the astrometric bias in barycentric correction.
To mitigate the bias in input astrometry for barycentric correction, we model the barycentric RV with offsets in position and proper-motion.We use the Precise EXOplanetology (PEXO) package to calculate the gradient of RV relative to these parameters, and infer these parameters through likelihood maximization.We apply this approach to the High Accuracy Radial velocity Planet Searcher (HARPS) RV data (named "RVbank") reduced by Trifonov et al. 2020 using the SERVAL pipeline (Zechmeister et al. 2018).To date, this is the first time that astrometric bias is systematically considered in barycentric correction.
This paper is structured as follows.In section 2, we introduce the data we use in this work; in section 3, we introduce the method for correction and testing; in section 4, we show the main results of the correction; in section 5, we conclude and discuss the results.

DATA
With a resolution of  = 115000, HARPS (Pepe et al. 2002) installed at the ESO 3.6 m telescope is one of the most productive high-resolution spectrographs in the detection of exoplanets using the radial velocity method.More than 6000 stars have been observed by HARPS according to the ESO archive 1 .Because the fiber of HARPS was changed in June 2015 (Curto et al. 2015), we use the RV data obtained before June 1st, 2015 (JD=2457174.5)to avoid the potential RV offset caused by the change of fiber.
The raw spectra obtained by HARPS are reduced using the official data reduction software to deliver the reduced RV data (Pepe et al. 2002).The HARPS data collected before 2020 are re-reduced by Trifonov et al. (2020) using the SpEctrum Radial Velocity AnaLyser (SERVAL) pipeline (Zechmeister et al. 2018), which integrates the barycentric correction code primarily based on BarCor (Hrudkova 2009) and routines specifically designed for various spectrographs.After removing the long-term drift caused by instrumental effects from the RV data, Trifonov et al. (2020) obtain the RVbank HARPS data for ∼ 3000 stars.In SERVAL, BarCor is revised to consider proper motion of stars.Meanwhile, the secular acceleration (SA) is calculated separately.To make the results consistent with other pipelines, e.g.PEXO or Barycorr, which consider SA in barycentric correction, we compare the combined results of barycentric Earth radial velocity (BERV) and SA of SERVAL with PEXO or Barycorr.
In Trifonov et al. 2020, Simbad (Wenger et al. 2000) astrometric data, which mainly comes from Gaia DR2, is used for the barycentric correction of the RVbank dataset.In this research, the input astrometric data used for most targets are from Gaia DR3.(Gaia Collaboration et al. 2022a).For some bright targets without Gaia data, we use revised Hipparcos (van Leeuwen 2007) data instead.The sources of reference astrometry are listed in our data table in the appendix section.
Simbad (Wenger et al. 2000) is used to crossmatch the RVbank targets with the Gaia and Hipparcos IDs.TOPCAT (Taylor 2005) is used to join the tables.We choose 266 targets according to the following criteria: (i) We remove the targets with RV variation > 30 m s −1 to avoid the contamination of binaries.Because we want to investigate the influence of wide companions on barycentric correction rather than the direct influence of binary companions on RV of the primary star; (ii) The number of RV data points exceeds 50.In order to correct the sub-m s −1 RV variation induced by biased astrometry, we need a large amount of RV data with ∼ 1 m s −1 precision; (iii) The RV data is collected in more than 30 nights.This is to guarantee that the RV data is relatively evenly sampled such that the potential annual bias is regularly sampled; (iv) The observational time span for RV is longer than 100 days.This is because the potential astrometric bias caused by proper motion bias increases linearly over time, resulting in a significant bias for RV if the observational baseline is sufficiently long.
The Hertzsprung-Russell diagram for these 266 targets is shown in Fig. 1.Most of the targets observed by HARPS are solar-type stars.Additionally, the majority of these targets are located within a distance of 100 parsecs.The selected targets show a gap in the vicinity of 4500 K, which is consistent with the distribution of nearby stars.

Barycentric correction
To detect the perturbation of a star from its companion, the Earth's motion are typically subtracted from the observed RVs.It is often assumed that the target star moves linearly with respect to the Solar System barycenter.This barycentric correction is conducted in a relativistic framework (Wright & Eastman 2014;Feng et al. 2019a).Current barycentric correction packages such as Barycorr (Wright & Eastman 2014) and PEXO (Feng et al. 2019a) can model barycentric correction to a precision better than 1 cm s −1 .They convert the Doppler shift of stellar spectrum measured on the Earth ( meas ) to the Doppler shift that would be measured by a virtual observer at the barycenter of the solar system ( true ) according to (1) 4000 4500 5000 5500 6000 6500 7000  Hertzsprung-Russell diagram for selected targets. eff refers to the effective temperature from Gaia DR3;   is the absolute magnitude derived from the Gaia DR3 magnitude and parallax.Color of the dots represents the type of the stars derived from  eff (Pecaut & Mamajek 2013).Colored dots are the selected targets and grey dots are the Gaia DR3 targets closer than 20 pc.
where   is the barycentric Doppler shift due to the Earth's motion around the solar system barycenter and various relativistic effects;  is the speed of light.The radial velocity measured at the solar system barycenter is  =  true . (2) We conduct barycentric correction using PEXO, a global modeling framework for nanosecond timing, microsecond astrometry, and mm s −1 radial velocities (Feng et al. 2019a).It is a global modeling framework which considers all the second-order effects including high-order Roemer delay and relativistic Doppler shift both in the solar system and in the target system.The philosophy of PEXO is to recommend simultaneous modeling of the barycentric motion of the observer and the reflex motion of the target to avoid bias caused by decoupling the Earth's motion and the reflex motion.Nevertheless, PEXO is suitable for high precision barycentric correction when the reflex motion is small or over long timescale.To optimize the astrometric offset, we use PEXO to numerically calculate the gradient of radial velocity with respect to the parameters of astrometric offsets.

Astrometric bias
Unbiased barycentric correction is needed to obtain the precise motion of a target.Fig. 2 illustrates the impact of biased astrometry on barycentric correction.Here we approximate the target motion by the motion of the target system's barycenter (TSB) of the target system because the positional difference between the target and the target system's barycenter does not introduce long-term RV variation or significant annual bias.Compared with the biased target astrometry given by Gaia, the barycentric motion has a much smaller astrometric offset with respect to the motion of the target and the offset does not accumulate with time.By assuming that the target is a single star, Gaia employs a five-parameter astrometric model to fit its raw data.The parameters of the model are right ascension (R.A.) denoted by , declination (Dec.)denoted by , proper motion in R.A. denoted by   , proper motion in Dec. denoted by   , and parallax denoted by .This approach could lead to an astrometric bias in cases where an unknown companion is present, as shown in Fig. 3, or instrumental bias (Lindegren et al. 2021).Only a small fraction of Gaia targets are analyzed without the single-star assumption in DR3 (443,205 in total 1.46 billion;Gaia Collaboration et al. 2022b).
Here is an order-of-magnitude estimation for this effect.The bias in the direction from the solar system barycenter to the target system barycenter (TSB) n having the bias Δ n would lead to an RV () bias (Δ) of where ì  ⊕ is Earth orbital velocity relative to the solar system barycenter as shown in Fig. 2.This term would introduce an annual bias.To quantify the effect of astrometric bias, we conduct a simulation to assess the magnitude of it, as shown in Fig. 4. In this simulation, every target star is assumed to have a mass of 1  ⊙ and to have an unresolved companion; The observational baseline of RV is assumed to be 10 years; The orbits of the systems are assumed to be circular; The target system's barycenter is assumed to have no proper-motion or parallax.The light emitted from the companions is considered in the simulation of photocenter motion.The G-band luminosity of the companions is calculated using the mass-luminosity relation given by Pecaut & Mamajek 2013.Based on Gaia sampling frequency, 36 points are sampled randomly in a 3 year observational baseline to simulate the along-scan position (abscissa) of Gaia source (Gaia Collaboration et al. 2022a).For a companion with a period  and mass , we uniformly sample its angular parameters and simulate its measured abscissae based on the photocenter of the system, and then fit a linear model to get the biased astrometry.The bias of relative radial velocity comes from the difference between the biased astrometry and the true astrometry of target star.Because the astrometric bias induced by companions is inversely proportional to the distance, the corresponding RV bias is also inversely proportional to the distance according to equation (3).
It is shown in Fig. 4 that this effect is most significant for target stars with long-period massive companions.This type of companions would lead to an annual RV bias of about 10 cm s −1 if the distance to the target system is 10 pc.This effect is suppressed for short-period companions, because Gaia is more likely to sample evenly within the phase of a short-period companion, leading to a fit with smaller bias.Meanwhile, this effect is slightly suppressed for stellar-mass companions.It is because the photocenter of the system would have a smaller nonlinear motion due to the light from the companion.However, for long-period stellar-mass companions, the light from the companion would make the system's photocenter far from the target star itself, leading to a significant bias.Because the astrometric bias is proportional to the parallax of the system, this effect is only significant for nearby stars.
To model the impact of biased astrometry on the reduced RV, we calculate the Jacobian of RV relative to astrometry offset parameters.Through numerical calculation using PEXO, we obtain the Jacobian of the relative radial velocity of a given target with respect to the astrometric correction, which would be useful for correcting the input astrometry.The elements of the Jacobian ( Ĵ) can be written as where   ( ) is the radial velocity at time   given astrometry  where  goes from 1 to the number of observation of this target . denotes an astrometric parameter, which could be R. A., Dec., proper-motion in R. A. or proper-motion in Dec..The impacts of potential bias in parallax () and systematic radial velocity on barycentric correction is not considerred in this work.This is because the influence of parallax would not increase linearly with respect to time.The difference of secular acceleration caused by parallax bias is negligible for Earth-like planet detection.We assume that the scale of astrometric correction is small so that the change of radial velocity of every observation can be linearly approximated by the first order Taylor expansion.Our model of ì  can be written as where ì  = (Δ, Δ, Δ  , Δ  ) T , is the amount of correction for this target;  is the RV offset; ì 1 is a  × 1 vector which all entries equals to 1.

Correction of astrometric bias in a Bayesian framework
The likelihood of our model for one target star can be written as where ì  is a  × 1 column vector, which contains the radial velocity data of this target star; ì  = ( ì  T , ) T is a parameter set to be estimated in this model; Σ = diag( 2 1 ,  2 2 , ...,  2  ), where   is uncertainty of the th radial velocity observation provided by HARPS; |Σ| is the determinant of Σ.The prior of our model can be written as where is the error bar given by Gaia DR3.According to Bayes' theorem, the posterior can be written as  To estimate ì , we employ the Maximum a posteriori (MAP) approach.By calculating the derivative of the logarithm of the posterior and equating it to zero, we obtain the solution By marginalizing the posterior and using the property of multivariate normal distribution, we can obtain the covariance matrix of correction ì  by The square root of the diagnal term of Σ ì  is the uncertainty for each correction.By using equation ( 9) and ( 10), we can obtain the correction ì  and its covariance Σ ì  .Significance of the correction is measured by Bayes factor (BF) estimated using Bayesian information criteria (BIC).BIC is used for Bayesian model selection.It is defined as where L is the maximized log likelihood;  is number of parameters of the model.For the astrometric correction model which includes , ,   ,   , and , we have  = 5; For the model with no astrometric correction which only includes a constant offset , we have  = 1.Under the assumption of Laplace approximation (Kass & Raftery 1995), BF is derived from BIC as follows where BIC i and BIC j represents the BIC value of the compared model and the baseline model.In this research, we need to compare the following models: model 0 (GDR2BarCor), model 1 (GDR3PEXO) and model 2 (OptAstroPEXO).GDR2BarCor is the constant offset model based on HARPS RVbank data; GDR3PEXO is the constant offset model based on PEXO reduced RV data which use Gaia DR3 astrometry; OptAstroPEXO is the astrometric correction model which is described previously in this section, see Table 1.We need to notice that GDR2BarCor and GDR3PEXO do not optimize astrometric parameters but only use different barycentric correction code and input astrometry, and OptAstroPEXO is fitted to the RV data for optimizing astrometric parameters thus having 4 more parameters (change of R.A., Dec., proper motion offsets in R.A. and Dec.) than GDR2BarCor and GDR3PEXO.There are variety of intermediate models between GDR3PEXO and GDR2BarCor which consider different barycentric correction code and input astrometry.They are not included in the model comparison framework due to their minimal differences with the GDR3PEXO model, as will be discussed in section 4.1.
• No model is favored if none the above rules are satisfied.
These criteria are used because the condition lnBF>5 is a robust criterion for model comparison (Kass & Raftery 1995), and also see appendix C. Using these criteria, most targets don't have a favored model, see section 4.
To evaluate the reduction in Earth motion after the astrometric correction, Bayes factor periodogram (BFP) is used to detect periodic signals in RV data.The value of the BF of periodogram is defined by equation ( 11) and ( 12), where BIC j and BIC i represent the BIC of a pure noise model and a periodic model respectively.The criterion ln BF = 5 is often chosen to be the threshold for significant periodic signals (Feng et al. 2016b).The BFP is calculated using Agatha (Feng et al. 2017b), which is a framework of periodograms with different noise models.Circular orbital model (the eccentricity is assumed to be 0) with white noise baseline model is used in this work.The workflow of the paper is shown in Fig. 5.

RESULTS
In this section, we first compare intermediate models for the 266 HAPRS targets to investigate the impact of pipeline difference and input astrometry on the reduced RV data.We then compare GDR2BarCor (model 0), GDR3PEXO (model 1), and OptAstroPEXO (model 2) to study whether the optimization of input astrometry would improve the RV data reduction precision.Finally, we discuss the results for eight targets that show significantly improved RV precision after optimizing input astrometry or changing reduction pipeline.

Model comparison
We calculate the lnBFs of intermediate models and describe the model comparison shown in Fig. 6 as follows: • GDR2BarCor and GDR2Barycorr.These two models share the same input astrometry but with different barycentric correction codes.As seen in the right panel of Fig. 6, the lnBF distribution (blue lines) is quite asymmetric, with 31 targets favoring GDR2Barycorr and 17 favoring GDR2BarCor.Such asymmetry is more than 3 away from Poisson noise and the preference for Barycorr over BarCor is likely caused by the fact that the BarCor does not account relativistic effects as Barycorr does.The RV difference between these two models exceed 0.1 m s −1 for 4 targets and exceed 0.01 m s −1 for 55 targets, indicating the importance of pipeline selection in detecting Earth twins.
• GDR2Barycorr and GDR2PEXO.The RV difference between these models are less than 1 cm s −1 , consistent with the 1 cm s −1 RV precision of GDR2Barycorr.Such difference does not lead to asymmetry in the lnBF distribution (orange line) shown in the right panel of Fig. 6.Hence both PEXO and Barycorr are reliable in terms of conduct barycentric correction at the precision level of 1 cm s −1 .
• GDR2PEXO and GDR3PEXO.Only three targets shows RV difference of more than 1 cm s −1 due to change of Gaia DRs.The distribution of lnBF is symmetric, indicating little difference in barycentric correction due to the choice of Gaia DRs.After the comparison of intermediate models, we move on to investigate whether the optimized astrometry would improve the RV data reduction precision.Specifically, we compare OptAstroPEXO with GDR2BarCor and GDR3PEXO in the Bayesian framework to assess the impact of optimized astrometry.The distribution of inferred parameters of OptAstroPEXO combined with the parallax of the targets are shown in Fig. 7.It is apparent that the nearer targets have larger astrometric corrections.Meanwhile, there is a positive correlation between the magnitude of correction for positions and proper motions.
We calculate the lnBFs for each target and model.The results of model comparison and the magnitude of correction are shown in Table 2.Among the 266 targets, 13 of them favor GDR2BarCor; 19 of them favor GDR3PEXO; 8 of them favor OptAstroPEXO, where the "favored model" is defined in section 3.3.Considering the significance of RV bias correction, we find 141 targets having bias with slope of linear trend with signal-to-noise ratio (SNR) greater than 3 when using GDR3PEXO, and 143 targets have this property when using OptAstroPEXO.Among them, 3 targets have slope of linear trend ( 10 ) higher than 0.01 m s −1 yr −1 using GDR3PEXO, and 9 targets have slope of linear trend ( 20 ) higher than 0.01 m s −1 yr −1 using OptAstroPEXO.Meanwhile, we find 259 targets having an annual bias with an SNR greater than 3 when using GDR3PEXO or OptAstroPEXO.Among them, 11 targets have annual variation semi-amplitude ( annual 10 ) larger than 0.05 m s −1 when using GDR3PEXO, which is comparable with the RV signal of Earth-like planets, and 14 targets have annual variation semi-amplitude ( annual 20 ) larger than 0.05 m s −1 when using OptAstroPEXO.Among this 14 targets, 10 of them are within the distance of 10 pc.
Fig. 8 shows the distance and the number of nights with RV observation for selected targets.The 8 targets that have significant astrometric correction (favor OptAstroPEXO) are closer than 20 pc, and most of their RV data are obtained in more than 100 nights.The targets without significant astrometric correction do not display such preference.By calculating the BFPs of different models, we find that the annual signals of 96 targets are attenuated through the utilization of GDR3PEXO or OptAstroPEXO.The criterion for attenuation is defined as a decrease in the BF of the periodogram for signals corresponding to periods of 365.25 days and 182.625 days.The decrease of BF of periodogram for these targets is shown in Fig. 9.It is shown that the decrease of the periodic signal is larger for nearby targets, which agrees with the intuition that astrometric bias is more significant for nearby targets.
According to equation ( 3), the annual bias in barycentric correction is proportional to the line-of-sight direction.Because the bias of proper-motion would increase the bias of line-of-sight direction linearly with respect to the observational baseline Δ, it is expected that the correction for proper-motion should be larger than the correction for position, i.e.   account for this astrometric bias.The Gaia DR3 RUWE for this target is larger than 1.4, indicating the existance of unknown companions (Gaia Collaboration et al. 2022a).The RV bias of it might thus come from the bias of Gaia astrometry.The semi-amplitude of the annual variation of this target is  annual 20 = 5.00 ± 0.02 cm s −1 , which is comparable to the signal of Earth-like planets.
• 61 Vir (GJ 506).61 Virginis (GJ 506) is a G-type star 8.5 pc away from the Earth.It is included in the target star list of the HWO (Mamajek & Stapelfeldt 2023) ).However, these planets can not account for the significant RV bias.The G-band magnitude of it is 4.5 mag.Such a bright target would lead to pixel saturation which can not be fully mitigated by reducing integration time controlled by blocking gates.This issue may lead to astrometric bias (Lindegren 2020).
The correction results for this target is shown in Fig. 10.The top left panel shows that the RV data is well sampled in the observational baseline of ∼ 10 yr; The bottom left panel shows that the pattern for RV correction is an annual signal superposed on a linear trend, which is a typical pattern for correction for all the targets; The top right panel shows that the signals of all its confirmed planets are significant; The bottom right panel shows that the annual signals are attenuated after the correction.Meanwhile, the significance of all confirmed exoplanets is enhanced after using the corrected astrometry.This is because the reduction in annual signals increase the relative significance of non-annual signals, including the signals of planets.
• GJ 616.GJ 616 is a G-type star 14 pc away from the Earth.It is included in the target star list of the HWO (Mamajek & Stapelfeldt 2023).Laliotis et al. 2023 find a super Earth candidate using radial velocity method.The mass for this exoplanet candidate is  ≥ 6.77 ± 0.86 ⊕ and its period is  = 19.8777± 0.0062 d.The G-band magnitude of it is 5.33 mag, which is probably suffering from the systematics of bright stars.Meanwhile, there exists possibility that the astrometric bias is attributed to an unidentified long-period massive companion.
• GJ 691.GJ 691 ( Arae) is a G-type star 15.6 pc away from the Earth.It is included in the target star list of the HWO (Mamajek & Stapelfeldt 2023).It has 4 confirmed exoplanets detected using radial velocity (Pepe et al. 2007).Its heaviest planet GJ 691 e is a planet with about 1.97 Jupiter mass and a period of 4019 days, which is too small to account for the astrometric bias and radial velocity mitigation.The G-band magnitude for this object is 4.94 mag, which could potentially introduce astrometric bias.The semi-amplitude of the annual variation of this target is  annual 20 = 8.37 ± 0.06 cm s −1 .• GJ 785.GJ 785 is a K-type star 8.8 pc away from the Earth.It is included in the target star list of the HWO (Mamajek & Stapelfeldt 2023).It has two confirmed Neptune-sized planets whose orbital periods are  b = 74.278± 0.035 d and  c = 549.1 ± 4.5 d, and mass is  b ≥ 14.28 +0.64  −0.63  ⊕ and  c ≥ 14.96 +1.21 −1.18  ⊕ respectively (Pepe et al. 2011).These companions can not explain the astrometric bias.The G-band magnitude of it is 5.5 mag, which indicates that the bias may due to the Gaia astrometry.
• GJ 845.GJ 845 is a K-type star 3.6 pc away from the Earth.It is included in the target star list of the HWO and it is considered to be a good (tier A) target (Mamajek & Stapelfeldt 2023).It have a binary brown dwarf companion and a planet of about 3 Jupiter mass (Feng et al. 2019b).According to Fig. 4, these companions would lead to an RV bias of ∼ 0.1 m s −1 , which is comparable with the annual signals of Earth-like planets.Meanwhile, the semi-amplitude of the annual variation of this target is  annual 20 = 15.89 ± 0.4 cm s −1 , which is consistent with the results of Fig. 4.This indicates that the bias of GJ 845 is caused by its companions.
•  Ceti (HD 10700). Ceti (HD 10700) is a G-type star 3.7 pc away from the Earth with 4 planets (Feng et al. 2017a).It is included in the target star list of the HWO (Mamajek & Stapelfeldt 2023).No known planets could induce such bias.This target have RUWE=2.6343,which is larger than 1.4, indicating strong systematics of the Gaia astrometry.Meanwhile, the semi-amplitude of the annual variation of this target is  annual 20 = 10.77± 0.15 cm s −1 , which is comparable with the annual signals of Earth-like planets.• Barnard's star (GJ 699).Barnard's star is an M-dwarf 1.8 pc away from the Earth.Ribas et al. 2018 claims the existence of a super-Earth planet orbiting near the snow line of Barnard's star which period is  = 233 d.The input astrometry of Gaia DR3 does not have any significant bias.The corrected RV and the BFP is shown in Fig. 11.The RV difference is mainly caused by the utilization of PEXO.It is shown that the periodic signal of the exoplanet candidate is more significant while the annual alias due to the Earth motion is reduced when using PEXO.
Table 3. Inferred quantities for selected targets. is the preferred model;  is the linear trend slope of the change of radial velocity corrected by the preferred model, i.e. the slope of the bottom left plot of Fig. 10;  annual represents the annual variation semi-amplitude of the preferred model; lnBF i0 is the lnBF comparing model i with the baseline model using RVbank.The results for the remaining targets are accessible to the public.A comprehensive description of the data product is shown in the appendix.

Target
Model

DISCUSSION AND CONCLUSION
We introduce a PEXO-based method to mitigate astrometric bias in barycentric correction for the detection of Earth-like planets.By applying this method to 266 HARPS targets, we find that 8 targets, including GJ 105 A, GJ 17, GJ 506, GJ 616, GJ 691, GJ 785, GJ 845 and HD 10700, have significantly biased astrometry.All of them are within a distance of 20 pc.The RV bias in GJ 105 A and GJ 845, characterized by a strong linear trend and an annual bias of ∼ 0.1 m s −1 , is likely attributed to their known companions.For GJ 17, GJ 506, GJ 616, GJ 691, GJ 785 and HD 10700, they are bright targets without known massive companions.Their astrometric bias could be attributed either to the Gaia instrumental systematics or to unidentified long-period massive companions.We show the details of radial velocity correction and periodogram analysis for GJ 506.The significance of all its confirmed exoplanets is enhanced after the astrometric correction, because the astrometric correction reduce the RV systematis as well as the annual variation.
For the remaining 226 targets, neither input nor optimized astrometry is preferred probably due to lack of well-sampled high precision RVs over long time span.The periodogram analyses reveal that the annual RV variation of 96 targets is decreased when PEXO in combination with GDR3 or corrected astrometry.14 among the 266 targets have annual bias of at least 0.05 m s −1 using corrected astrometry, and 10 of them are within the distance of 10 pc.Our successful correction of these bias demonstrate the importance of unbiased barycentric correction in detecting Earth-like planets.
Based on a comparison of intermediate models that share the same input astrometry but have different barycentric correction code, the Barycorr and PEXO pipelines have at most 1 cm s −1 difference in barycentric correction and both are precise enough for detecting Earth-twin signals.Both of them are favored over the BarCor pipeline that is used in SERVAL by Trifonov et al. (2020) to process RVbank data.Because BarCor does not account for relativistic effects, it can introduce sub-m s −1 RV bias in barycentric correction.In particular, we show that the barycentric correction for Barnard's star can be significantly improved by replacing BarCor with PEXO.Therefore, we highly recommend using relativistic pipelines such as Barycorr and PEXO for barycentric correction when detecting Earth twins.
In future work, a simultaneous modeling of barycentric and reflex motion is needed to totally remove the assumptions of barycentric correction in the RV analyses, as implemented in tools such as TEMPO2 (Edwards et al. 2006) and PEXO (Feng et al. 2019a).

Figure 1 .
Figure1.Hertzsprung-Russell diagram for selected targets. eff refers to the effective temperature from Gaia DR3;   is the absolute magnitude derived from the Gaia DR3 magnitude and parallax.Color of the dots represents the type of the stars derived from  eff(Pecaut & Mamajek 2013).Colored dots are the selected targets and grey dots are the Gaia DR3 targets closer than 20 pc.

Figure 2 .
Figure 2. Effects of biased astrometry.E1, E2 and E3 represent the positions of the Earth at different epoch.n represents the direction from the Earth to the target system's barycenter (TSB).The biased target astrometry (T) provided by Gaia is expressed through five-parameter solutions.Δ n represents the bias in the direction of the target's barycenter due to Gaia astrometry.

Figure 3 .
Figure3.Astrometric bias caused by utilization of a five-parameter astrometric model to fit the data of a target with an unidentified companion.The orange curve is simulated by accounting for both proper motion and reflex motion due to a long-period companion without considering parallax.The dashed red line represents the best-fit position and proper motion solution.The solid blue line represents the barycentric motion of this target.The slope of a line is the proper-motion in this direction.The observation time comes from the Gaia GOST result of Barnard's star.

Figure 4 .
Figure 4. Magnitude of bias caused by an unresolved companion around a sun-like star. is the period of the massive companion;   is the mass of the companion;  is the scale-free maximum radial velocity bias normalized by distance.For a star at a distance of , its radial velocity bias caused by the companion is /.

Figure 5 .
Figure 5. Workflow of the correction of the astrometric bias in RV data.

Figure 6 .
Figure 6.Comparison of intermediate models.Left: Distribution of the barycentric correction difference, where   =   represents the barycentric velocity obtained using the provided input astrometry and data reduction pipelines. Δ  corresponds to the scattering of the difference in   for a target star.The error bars indicate uncertainties caused by poisson fluctuation.Right: Distribution of lnBFs when utilizing different input astrometry and data reduction pipelines.The red dashed lines are the thresholds for significant difference.

Figure 7 .
Figure 7. Correlation between the astrometric correction and the parallax of the targets.Each point in the figure correspondence to a HARPS RVbank target.

Figure 8 .Figure 10 .
Figure 8. Distribution of selected targets on distance and number of nights with RV measurement. is the distance of the target stars obtained using parallax.The color of the points and histograms represents the favored model of the target.

Figure 11 .
Figure11.Change of radial velocity and the BFP when using GDR3PEXO.Configurations of this figure is similar to Fig.10, but for Barnard's star.In the top left panel, an offset of 16.3 m s −1 is applied for the corrected RV.In the top right panel, an offset of 22.8 is applied for the BFP after correction.

Table 1 .
Models used in this work

Table A1 .
Header for the new astrometry table .A. correction and R. A. proper-motion correction deg mas yr −1 cov_delta_ra_pmde Covariance between R. A. correction and Dec. proper-motion correction deg mas yr −1 cov_delta_de_pmra Covariance between Dec. correction and R. A. proper-motion correction deg mas yr −1 cov_delta_de_pmde Covariance between Dec. correction and Dec. proper-motion correction deg mas yr −1 cov_delta_pmra_pmde Covariance between R. A. proper-motion correction and Dec. proper-motion correction mas 2 yr −2