Modeling and Calibration of Gaia, Hipparcos, and Tycho-2 Astrometric Data for the Detection of Dark Companions

Hidden within the Gaia satellite’s multiple data releases lies a valuable cache of dark companions. To facilitate the efficient and reliable detection of these companions via combined analyses involving the Gaia, Hipparcos, and Tycho-2 catalogs, we introduce an astrometric modeling framework. This method incorporates analytical least-square minimization and nonlinear parameter optimization techniques to a set of common calibration sources across the different space-based astrometric catalogs. This enables us to discern the error inflation, astrometric jitter, differential parallax zero-points, and frame rotation of various catalogs relative to Gaia Data Release 3 (DR3). Our findings yield the most precise Gaia DR2 calibration parameters to date, revealing notable dependencies on magnitude and color. Intriguingly, we identify submilliarcsecond frame rotation between Gaia DR1 and DR3, along with an estimated astrometric jitter of 2.16 mas for the revised Hipparcos catalog. In a thorough comparative analysis with previous studies, we offer recommendations on calibrating and utilizing different catalogs for companion detection. Furthermore, we provide a user-friendly pipeline (https://github.com/ruiyicheng/Download_HIP_Gaia_GOST) for catalog download and bias correction, enhancing accessibility and usability within the scientific community.

1. INTRODUCTION While our solar system boasts multiple rocky planets positioned closer to the Sun and several giant planets situated on wider orbits, the more than 4000 known planetary systems exhibit a variety of different configurations.The pursuit of Earth-like planets, particularly those resembling Earth in size and located within the habitable zone of Sun-like stars, is often regarded as the "holy grail" of exoplanet science.Additionally, giant planets play a crucial role in shaping planetary system architectures (e.g., Gomes et al. 2005;Tsiganis et al. 2005) and influencing the potential habitability of Earth-like worlds (e.g., Lunine 2001;Horner et al. 2010Horner et al. , 2020)).Detecting Earth-like planets poses challenges due to their faint signals (e.g., Hall et al. 2018 andGe et al. 2022), while our current capabilities enable relatively straightforward detection of cold giants (Wittenmyer et al. 2020;Feng et al. 2022;Laliotis et al. 2023).
The emergence of the Gaia satellite (Gaia Collaboration et al. 2016a,b, 2018, 2021, 2023a) marks a golden age for the detection of cold giants.Leveraging Gaia's astrometric data, combined with other astrometric and radial velocity data, provides a powerful tool for this purpose (Snellen & Brown 2018;Kervella et al. 2019;Feng et al. 2019;Brandt et al. 2019).Specialized tools such as orvara (Brandt et al. 2021) and BINARYS (Leclerc et al. 2023) have been developed to Feng et al. rotation and zero-point parallax in these catalogs (Brandt 2018;Kervella et al. 2019;Brandt et al. 2023).However, in certain cases, calibration is integrated into the analysis itself, addressing potential biases alongside companion signals (Feng et al. 2021).
The techniques developed for exoplanet detection can also be applied in the search for dark companions such as black holes (e.g., El-Badry et al. 2023), neutron stars (e.g., Shahaf et al. 2023), white dwarfs (e.g., Ganguly et al. 2023), and brown dwarfs (e.g., Feng et al. 2022).Detecting these dark companions opens a new window to understand the formation and evolution of massive stars and binaries (Heger et al. 2003) and to test whether there is a 2 − 5 M⊙ mass gap (Kreidberg et al. 2012;Lam et al. 2022;Ye & Fishbach 2022) between neutron stars and black holes.
In the absence of Gaia epoch data, the proper motion anomaly derived from Gaia and Hipparcos proper motions and positions is commonly used to constrain stellar reflex motion with decade-long orbital periods (e.g.Kervella et al. 2019;Brandt et al. 2021).However, this method faces significant information loss in transforming astrometry into proper motion anomaly, making it inadequate for constraining the mass of Jupiter analogs 1 and distinguishing between retrograde and prograde orbits (Li et al. 2021).To address this limitation, a method proposed by Feng et al. (2023) utilizes Gaia DR2 and DR3 five-parameter astrometry, along with Hipparcos intermediate astrometric data (IAD), combined with radial velocity data.This approach successfully constrains the orbit and mass of nearby Jupiter analogs, such as ϵ Eridani b, with a mass of 0.76 +0.14  −0.11 MJup and an orbital period of 7.36 +0.04  −0.05 yr.Therefore, detecting nearby cold Jupiters is feasible, considering the present precision of Gaia DR2 and DR3 of about 0.1 mas for bright stars (Gaia Collaboration et al. 2023a), and the 24-year observational baseline established by Hipparcos and Gaia.This feasibility is contingent upon effectively mitigating biases and systematics in the astrometric catalogs to a comparable precision through calibration, either a priori or a posteriori.
Despite the vast dataset of one billion stars observed by Gaia, only a small fraction (fewer than 10,000) have been studied with high-precision spectrographs.To detect exoplanets and other dark companions solely through astrometric data, we have developed an efficient algorithm based on analytical χ 2 minimization.This method considers various combinations of astrometric catalogs, exploring combinations such as Gaia Data Release 1 (GDR1; Gaia Collaboration et al. 2016b), Gaia Data Release 2 (GDR2; Gaia Collaboration et al. 2018), Gaia Data Release 3 (GDR3; Gaia Collaboration et al. 2023a), Tycho-2 (TYC; Høg et al. 2000), the original Hipparcos catalog (HIP1; Perryman et al. 1997), and the revised Hipparcos catalog (HIP2; van Leeuwen 2007).
We do not consider Gaia Early Data Release 3 (GEDR3; Gaia Collaboration et al. 2021) since its astrometric data is the same as GDR3.Our Bayes factor (BF) is derived from χ 2 under the assumption of Laplace's approximation (Schwarz et al. 1978;Kass & Raftery 1995), leading to the development of a modeling framework for combined analyses of astrometric data.Additionally, we employ orbital solutions from the GDR3 non-single-star (NSS; Halbwachs et al. 2023;Holl et al. 2023;Gaia Collaboration et al. 2023b) catalog and stars with both Gaia and Hipparcos data as test sets to estimate the frame rotations of various catalogs relative to GDR3 and the error inflation within the astrometric catalogs.
This paper is organized as follows: Section 2 outlines the formulae for the modeling framework applied to various combinations of astrometric catalogs.Section 3 introduces the calibration procedure, while section 4 presents the calibration outcomes.Finally, our findings and conclusions are summarized in Section 5.

METHOD
In this section, we present the model framework designed for the integrated analysis of multiple astrometric catalogs.We introduce the circular reflex motion model and the simulation of Gaia epoch data using the Gaia Observation Forecast Tool (GOST) in subsection 2.1.The combined modeling of astrometric catalogs is categorized into three cases: GDR1+GDR2+GDR3, TYC+GDR2+GDR3, and HIP+GDR2+GDR3, which are detailed in subsections 2.2, 2.3, and 2.4, respectively.For each case, we elucidate the astrometric model and provide the analytical solution for circular reflex motion.Subsequently, in subsection 2.5, we augment the modeling framework with non-linear orbital parameters.The determination of the threshold for detecting a reflex motion signal induced by a companion is presented in subsection 2.6.

Astrometry model for circular orbits
In the case of unresolved binaries, Gaia observes the motion of the system's photocenter around the barycenter, focusing on the displacement of the combined center of light rather than tracking the reflex motion of individual stars within the binary system.Given a mass-luminosity function of F (m), the angular semi-major axis of the photocentric motion is where m1 and m2 are respectively the masses of the primary and the secondary, F1 and F2 are the observed fluxes of the primary and the secondary, a b is the binary semi-major axis, d is the distance derived from parallax, d = Au/ϖ where Au ≡ 1 au.For dark companions, F2/F1 ≈ 0, the photocentric motion is equivalent to the reflex motion, and thus the angular semi-major axis of the reflex motion is ar 1 A Jupiter analog is an exoplanet of Jupiter's mass (0.5 to 2 M Jup ) orbiting a FGK star in a configuration akin to Jupiter, with a semi-major axis in the range 2.5 to 10 au.
To derive analytical formulae for linear regression of astrometric model, we assume that the orbit of a star is circular and could be described by (3) where ∆ti = ti − t ref is the time difference between epoch ti and the reference epoch t ref (GDR3 referernce epoch, J2016, by default), ϕ is the orbital phase, and P is the orbital period.The coordinates in the orbital plane are transformed into the sky plane using where A ′ , B ′ , F ′ , G ′ are scaled Thiele Innes constants, and are functions of inclination I, argument of periastron ωT of the target star, longitude of ascending node Ω.We multiply A ′ , B ′ , F ′ , G ′ by ap to define the Thiele Innes constants A, B, F, G.
Here we assume that the reflex motion is equivalent to the photocentric motion as long as the mass function is not derived from the Thiele Innes constants, as in the case of GDR3 non-single-star catalog (NSS).
For circular orbit, ωT can be arbitrarily chosen such that ϕ = 0. Hence we choose ϕj = 0 when tj = t ref and we have where C t i = cos 2π P ∆ti and S t i = sin 2π P ∆ti .Relative to the reference epoch tj, the motion of the target system barycenter (TSB) viewed from the solar system barycenter is where (α b * j , δ b j , µ b α,j , µ b δ,j ) represents the TSB astrometry at the reference epoch tj.
The above propagation of the TSB does not account for second or higher order effects such as perspective acceleration.We consider these effects a priori by subtracting them from various catalogs.Because Gaia orbital solutions use the along-scan (AL) coordinates of the star (i.e., abscissa or IAD), we project the combined motion of the star onto the AL direction to simulate abscissa2 .With the Gaia scan angle θ and AL parallax factor f AL from GOST, the synthetic Gaia abscissae is where ∆ti = ti − t ref , S p i ≡ sin θi, C p i ≡ cos θi.We define the parameters of TSB astrometry at the Gaia DR3 ref- ).With these definitions, the synthetic abscissae could be written as ⃗ ξ = κ b ⃗ β b + κ r ⃗ β r .There are five astrometric parameters for the TSB of targets in GDR2 and GDR3.For a three-catalog combination, the catalog cross-matched with GDR2 and GDR3 falls into one of the following categories: • GDR1 five-parameter solutions (G1P5):.Stars without TYC data in GDR1 only have R.A. and Decl.given.In this case, only position is fitted to the GOST synthetic data, denoted as ⃗ β b GDR1 = (∆ αGDR1, ∆ δGDR1, 0, 0, 0) T .The model is described in subsection 2.2.
• TYC: Given that TYC proper motions are derived through combined analyses of TYC and previous photographic catalogs (Høg et al. 2000), we use only the TYC R.A. and Decl. to constrain the orbital parameters.Although the abscissae of TYC are not provided, we model the position of TYC, and the analytical solution is described in subsection 2.3.
• Hipparcos: For stars with Hipparcos data, GDR1 solutions are not independent of the Hipparcos data.We utilize the Hipparcos IAD in combination with GDR2 and GDR3 catalog data to constrain binary orbits.The detailed modeling of Hipparcos IAD is described in subsection 2.4.

Modeling Gaia catalog data
The five-parameter model for synthetic abscissae for the Gaia DRj is where ⃗ βj is the astrometry fitted for Gaia DRj.The astrometry of Gaia DRj relative to the reference astrometry (Gaia DR3 by default) is ⃗ βj.We define the synthetic abscissae corresponding to the three Gaia DRs as where κGDRn is similar to κGDRn but the reference time is now the reference epoch of GDRn, and For GDR1 targets without Tycho data, only the positions are measured.To be compatible with GDR2 and GDR3 solutions, the GDR1 proper motions and parallax are fixed at zero and the covariance matrix is , where ρ is the the correlation between α * and δ for GDR1, the parallax uncertainty σϖ is set to 1000 mas, and the proper motion uncertainties, σµ α and σµ δ , are set to 1000 mas yr −1 .The χ 2 of the above model of abscissae is where Σ ξ is the covariance of abscissae.Assuming that the synthetic abscissae have constant measurement errors σ ξ , χ 2 ξ could be simplified as Transposing both sides of the above equation leads to and thus ⃗ The coefficient matrices of ⃗ β b and ⃗ β r are respectively and Because we do not consider perspective acceleration in the above linear modeling of synthetic GOST data, we propagate the GDR3 astrometry to the reference epochs of GDR1 and GDR2, and define the astrometric data as ⃗ y = ( ⃗ βGDR1, ⃗ βGDR2, ⃗ βGDR3) T .
We calculate the residual vector by subtracting ⃗ The χ 2 for this model of catalog astrometry is where , and Σ GDR3 are respectively the covariances for Gaia DR1, DR2, and DR3.The minimization of χ 2 GDR leads to The above equation array could be simplified as where Then the optimal parameters are ⃗ β br = η −1 ⃗ b and the corresponding covariance is C br = η −1 .The baseline model is simply a covariance-weighted mean of the catalog astrometry, where ) −1 measures the parameter covariance.

Modeling TYC data
For targets with TYC positional data, the target position relative to the barycenter is The vectorized version is where λ TYC = 1 0 0 0 0 0 1 0 0 0 (33) Hence the residual is where ⃗ yTYC ≡ ⃗ βTYC = (∆αTYC, ∆δTYC) is the TYC position relative to the position propagated from Gaia DR3 to the TYC R.A. and Decl.reference epochs.With the above definitions, λ and γ in eq.20 become The optimal parameter vector for the corresponding baseline model is where covariance

Modeling Hipparcos intermediate data
Because the astrometric offset parameters shown in eq.14 are defined with respect to GDR3, the Hipparcos raw IAD ( ⃗ ξraw) are transformed into relative Hipparcos IAD following where ⃗ β b HG is the difference between Hippacos astrometry and the astrometry propagated from the GDR3 epoch to the Hipparcos epoch.After this correction, the Hipparcos abscissa model is the same as that in eq.14.The Gaia scan angle θ is complementary to the Hipparcos scan angle ψ, i.e. θ = π/2 − ψ.The residual is Here we use κ b HIAD (with GDR3 as the reference epoch) instead of κb HIAD because the Hipparcos abscissae are already modified to be compatible with the GDR3 solution in eq.37.
The χ 2 for δ ⃗ ξHIAD is where Σ −1 HIAD is the covariance matrix of Hipparcos IAD.This covariance is diag(σ 2 1 , σ 2 2 , ..., σ 2 N HIAD ), where NHIAD is the number of Hipparcos epochs.The minimization of χ 2 leads to where ⃗ β br = ( ⃗ β b , ⃗ β r ) T .To model Hipparcos IAD together with Gaia DRs, we insert eq.40 into the definition of λ in eq.20, γ in eq.21, Σ −1 β in eq.24 and ⃗ y in eq.23 as follows: With these new definitions, eq.27 is still valid to find solutions of ⃗ β b and ⃗ β r .The optimal parameter vector for the corresponding baseline model becomes where covariance

Modeling eccentric orbits
By sampling orbital periods and minimizing χ 2 for each period, we find astrometric signals with circular orbits.This linear least-square optimization provides initial parameters for more robust non-linear regressions to find solutions for eccentric orbits.The Keplerian motion of system photocenter in the orbital plane is where P is orbital period, e is eccentricity, Ei is the eccentric anomaly at time ti and can be determined by solving the Kepler's equation: where Mi is the mean anomaly at ti.The mean anomaly is where M0 is the mean anomaly at the reference epoch.Compared with circular orbits, a Keplerian orbit has two nonlinear parameters, e and M0, in addition to orbital period P .Hence there would be three nonlinear parameters for a Keplerian model.We define this set of nonlinear Keplerian parameters as ⃗ β n ≡ (P, e, M0) T .Hence the parameter vector for a full astrometry model is For a given ⃗ β n , C t i and S t i in eq.7 are respectively replaced by C tr i ≡ cos Ei − e and S tr i ≡ √ 1 − e 2 sin Ei.Following the procedures introduced in section 2.2, 2.4, and 2.5, the linear parameters can be analytically optimized for a given ⃗ β n .Thus we only need to optimize ⃗ β n nonlinearly using algorithms such as the Levenberg-Marquardt optimization algorithm (Levenberg 1944;Marquardt 1963).The Thiele Innes elements A, B, F , G, can be converted to the Campbell elements, ap, I,ωT , and Ω although ωT and Ω are respectively degenerated with ωT + π and Ω + π.

Detection threshold
The preceding subsections outline the analytical formulae for computing optimal parameters through χ 2 minimization.Assuming uniform priors for model parameters and employing Laplace's approximation (Schwarz et al. 1978;Kass & Raftery 1995), we transform the minimum χ 2 of two models into the logarithmic Bayes Factor (lnBF) using the following equations: Here, L1 and L2 represent the maximum likelihood for models 1 and 2, while χ 2 1 and χ 2 2 denote the respective minimum χ 2 values for models 1 and 2. The variable ν stands for the number of extra free parameters in model 2 compared to model 1, and N represents the number of data points.We consider lnBF> 5 as the threshold for identifying statistically significant signals.This is based on section 3.2 and table 2 of Kass & Raftery (1995), where the criteria for decisive evidence in favor of a hypothesis are stated as 2lnBF > 10 or lnBF > 5 or BF > 150.We also note that Feng et al. (2016) shows that lnBF > 5 is suitable for identifying signals across various noise properties, particularly in the context of exoplanet detection.However, we note that our study, not yet having developed the modeling framework into a periodogram for signal detection, is not highly sensitive to the choice of threshold.We include a detection threshold primarily for the sake of completeness.

Calibration sources
The models introduced in section 2 assume unbiased catalog data, a presumption that might not hold true, especially in combined analyses of different catalogs.While it's possible to model bias concurrently with stellar reflex motion on a case-by-case basis, as demonstrated in Feng et al. (2019), such bias models often involve nonlinear parameters, impeding analytical optimization and diminishing the efficiency of periodogram computation.Consequently, we proactively address potential astrometric bias through calibration, leveraging both GDR3 NSS and single-star (SS) catalogs.While we opt for both NSS and SS as representations of distinct astrometric solutions, it is important to note that our choice of calibration sources lacks comprehensiveness in terms of capturing various magnitudes, colors, sky positions, and other pertinent parameters.In subsection 4.3, we will evaluate the extent of contamination from undetected binaries in the SS catalogs by examining the correlation between calibration parameters and the renormalized unit weight error (RUWE; Lindegren et al. 2018).
Recognizing that Gaia systematics may vary with solution types (Lindegren et al. 2021a), we cross-match the GDR3 singlestar catalog (G3SS)4 with GDR2 and the Hipparcos catalog to obtain a sample of 77,138 single-star targets (termed "G3SS") for the calibration of bright and nearby stars.This sample is further cross-matched with GDR2, G1P2, G1P5, TYC, HIP1, and HIP2, resulting in SS2, SS2P2, SS2P5, SS2T,SS2H2,and SS2H1 catalogs,encompassing 77,138,6,763,70,676,74,216,77,136,and 77,138 targets,respectively.Because the Hipparcos and TYC sources were observed by Gaia while most of the Gaia NSS sources were not observed by Hipparcos or TYC, the SS calibration catalogs have similar sample size while the NSS calibration catalogs have quite different sample size.

Calibration models
We employ the calibration sources introduced in the previous section to ascertain the relative frame rotation and offset between GDR3 and other catalogs.Given the primary objective of this series of papers-to identify companions using multiple astrometric catalogs-our focus lies solely on the differential frame rotation and parallax zero-point.Considering that our typical detection involves companions closer than 1 kpc, the absolute parallax zero-point of GDR3, approximately 0.05 mas (Lindegren et al. 2021a), does not significantly impact the detection process.
The relative difference between two reference frames is measured by a rotation vector ⃗ ω and a constant offset vector ⃗ ϵ: where tGDR3 is the reference epoch of GDR3.We denote ⃗ ϵ(tGDR3) by ⃗ ϵ0 ≡ (ϵx, ϵy, ϵz) T .Additionally, the zero point of parallax is found to be non-zero for Gaia catalogs (Lindegren et al. 2018;Arenou et al. 2018;Lindegren et al. 2021a).Without using quasars to calibrate parallax zero-point like Lindegren et al. (2021a), we derive the differential zero-point parallax (δϖ) of a catalog relative to GDR3.The astrometric bias caused by the differential frame roration and zero-point parallax is where˜denotes astrometric parameters measured in a given frame while astrometric parameters without this subscript denotes parameters measured in the GDR3 frame.Note that we define the differential rotation and parallax zero-point in a convention such that the bias is subtracted from the given catalog.We vectorize the above equations and get the differential astrometric bias (⃗ y d ) caused by differential frame rotation and parallax zero-point, where The differential frame rotation and parallax zero-point should be subtracted from the original data ⃗ y to get the corrected data in the GDR3 reference frame.
Because we only consider differential frame rotation and parallax zero-point, GDR3 astrometry is treated as "unbiased" and thus the last seven elements in ⃗ β d and the last five rows in κ d are all zeros.To model the systematics of GDR1 and GDR2 relative to GDR3 simultaneously, we can respectively add ⃗ β GDR13 and κ GDR13 to ⃗ β d and κ d .Hence eq.49 becomes When CAT1 is TYC, κ GDR13 and ⃗ β GDR13 are respectively replaced by κ TYC and ⃗ β TYC while the forms of eqs.52 and 53 do not change.If CAT1 is HIP1 or HIP2, the Hipparcos abscissae (HIAD) induced by frame rotation are and where κ HIP has the same form as κ GDR23 .Then eq.23 becomes When G3NST is used for calibration, β r is given by the catalog and is subtraced from the catalog astrometry, ⃗ y ′ = ⃗ y − γ ⃗ β r , and eq.57 becomes When G3NSB is used, the barycentric astrometry is already given, and thus the GDR3 part of the data vector would not change, i.e. ⃗ y ′ GDR3 = ⃗ y GDR3 .If G3SS is used, the data vector would not change, i.e. ⃗ y ′ = ⃗ y.
We infer the differential calibration parameters β d as well as the barycentric astrometry ⃗ β b for a number of calibration sources simultaneously.The data vector for a sample of n calibration sources is where Similar to eq. 27, the linear regression of the astrometric model leads to where where CAT1 could be G1P2, G1P5, TYC, HIP1, or HIP2.If no CAT1 is provided, Considering that the covariance given by a catalog is probably underestimated, we model the error inflation and jitter using a and b, respectively.The covariance of a given target becomes , where ρ is the correlation matrix.
Finally, by minimizing we get the optimal calibration parameters ⃗ β d .The corresponding formulae for the above χ 2 minimization is similar to that shown in eq.27.
The lnBF for models with (model 2) and without (model 1) parameters of frame rotation and parallax zero-point (i.e.⃗ β d ) is where n b ≡ n bd − n b , n d , and n bd are respectively the number of parameters in ⃗ β b sample , ⃗ β d sample , and ⃗ β bd sample , N sample = n × N single is the total number of data points for a sample of n calibration sources, and N single is the number of data points for a single source.If CAT1 in eq.3.2 is not given, the differential frame rotation and parallax zero-point of GDR2 relative to GDR3 will be modeled by 7 parameters (i.e., n d = 7).If CAT1 is G1P5, HIP1, or HIP2, there would be 7 calibration parameters in addition to the GDR2 calibration parameters.If GDR2 calibration parameters are fixed at their optimal values, only 7 parameters need to be optimized, n d = 7.If CAT1 is G1P2 or TYC, the proper motion of CAT1 is unknown.Hence only the offset between frames and the differential parallax zero-point at the reference epoch (i.e., ϵx, ϵy, ϵz, and δϖ) would be considered, and n d = 4.

Calibration procedure
We derive the calibration parameters for catalogs through the following sequential steps: 1) Data Collection: Gather GOST and catalog data for various calibration sources.
2) Matrix and Vector Calculations: Compute the matrices and vectors introduced in sections 2 and 3.2.
3) Random Subsampling: Randomly select 100 targets from the entire sample to create a subsample.Repeat this subsampling process at least 1000 times to generate an ensemble of subsamples; 4) Error Inflation and Jitter Sampling: Sample error inflation (1 ≤ a ≤ 2) and jitter (0 ≤ b ≤ 2 mas or mas/yr)5 parameters, creating a grid with a bin size of 0.02.

5)
Optimization and χ 2 Calculation: For each subsample and bin, infer optimal parameters and calculate the corresponding χ 2 values for both the combined model of barycentric astrometry and calibration (Model 2) and the model of barycentric astrometry alone (Model 1).Eliminate subsamples with χ 2 values exceeding the median χ 2 of the entire ensemble of subsamples for each bin.
6) Maximum lnBF Calculation: Calculate the maximum lnBF for each bin according to eq. 65, identify the optimal (a, b) at the globally maximum lnBF (lnBFmax), and determine the uncertainty of (a, b) based on the range that satisfies lnBFmax−lnBF < 0.5 (equivalent to a 1-σ confidence level, assuming Laplace's approximation).
7) Parameter Distribution Calculation: With a and b fixed at their optimal values, repeat step 5 at least 10 times to obtain the distribution of optimal ⃗ β d sample .Calculate the mode, 16%, and 84% quantiles for each parameter.
Figures 1 and 2 illustrate the lnBF distributions over (a, b) and the distribution of ⃗ β GDR23 based on analyses of the NST, NSB, and SS calibration sources, respectively.Notably, the preference for zero jitter and no error inflation holds true for both GDR2 and GDR3, and this conclusion remains robust across different choices of calibration sources, as demonstrated in Table 1.
The analysis of NST, NSB, and SS calibration sources reveals significant frame rotation and parallax offset for GDR2 relative to GDR3.The SS values align with those reported by Brandt (2018), Lindegren (2020)    Feng et al.
Having set (aGDR2, bGDR2), (aGDR3, bGDR3), and ⃗ β GDR23 to their optimal values, we proceed to optimize the calibration parameters for CAT1 using various combinations of CAT1, GDR2, and GDR3-based catalogs.The analysis, detailed in Table 1 and depicted in Fig. 7, reveals negligible error inflation and jitter for G1P2 and G1P5 across different calibration sources.The G1P2 frame exhibits an offset of -0.13 +0.05 −0.03 mas relative to GDR3 along the y-axis, as determined through NSTP2 analysis.This result is consistent with values obtained from NSTP2 and notably more precise than those derived from SSP2.This precision improvement is likely attributed to SSP2 having fewer calibration sources compared to NSTP2 and NSBP2.
The G1P5 frame exhibits a consistent offset of (ϵx, ϵy, ϵz) = (0.39 +0.04 −0.05 , −0.17 +0.05 −0.07 , 0.12 +0.05 −0.06 ) mas relative to GDR3, as determined through analysis of SSP5, as depicted in Fig. 2. Notably, the frame-rotation speed ω of G1P5 relative to GDR3 aligns with zero.The NSTP5 and NSBP5 calibration sources, however, do not ascertain these parameters with the same high precision, likely due to their limited sample size and the inappropriate subtraction of a photocentric motion from TGAS proper motions derived partially from TYC (Michalik et al. 2015).Given the potential "contamination" from TYC, we recommend utilizing only the reference position in G1P5 for detecting companions.
As indicated in Table 1, the error inflation and jitter for TYC are not in alignment with zero at the 1σ level based on analyses of NSTT and NSBT.However, upon inspecting the lnBF distributions for NSTT and NSBT depicted in Fig. 6, it becomes evident that these distributions exhibit significant non-Gaussian characteristics, with a seemingly random occurrence of high lnBF regions.In contrast, the lnBF distribution for SST appears smoother, with lnBF variations consistently below 3, indicative of negligible error inflation and jitter.The frame rotation of TYC relative to GDR3 is also found to be consistent with zero, as detailed in Table 1 and illustrated in Fig. 7.
HIP2 exhibits negligible error inflation, aligning with the less than 20% error inflation previously determined by Perryman et al. (1997) and Brandt et al. (2023).An analysis of NSTH2 reveals a jitter of 2.16 +0.12 −0.00 mas (or mas yr −1 ) for HIP2, a result consistent with the value of 2.25±0.04mas (or mas yr −1 ) reported by Brandt et al. (2023).Fig. 6 illustrates that NSBH2 and SSH2 also exhibit elevated lnBF around b = 2 mas (or mas yr −1 ), although not as prominently as NSTH2.As presented in Table 1 and Fig. 8, the differential calibration parameters of HIP2 deviate from zero at the 2σ level and align with the values provided by Lindegren et al. (2016) at the 1.5σ level.Additionally, our analysis indicates a negative differential zero-point parallax based on NSTH2 and NSBH2 calibrations, consistent with the value determined by Fabricius et al. (2021).
As outlined in Table 1, the error inflation of HIP1 is negligible, while the jitter of HIP1 is determined to be 0.08, 0.20, and 0.02 mas (or mas/yr) based on analyses conducted with NSTH1, NSBH1, and SSH1, respectively.Examining Fig. 6, it is evident that the (a, b) values for HIP1 are largely consistent with (1, 0).However, it's noteworthy that the differential calibration parameters for HIP1 are not as precisely determined as those for HIP2.Conversely, Figs. 2 and 8 reveal a pronounced correlation between ⃗ ω and ⃗ ϵ for G1P5, HIP1, and HIP2.This correlation is likely influenced by the 24-year positional difference between TYC, HIP2, HIP1, and GDR3, which imposes a more stringent constraint on ⃗ β d than the difference in proper motion does.This leads to a degeneracy between ⃗ ω and ⃗ ϵ (see eq. 48).In Fig. 3, we juxtapose the calibration parameters derived in this study with values reported in the literature.Given that literature values for frame rotation and zero-point parallax are referenced to specific realizations of the International Celestial Reference Frame (ICRF), we adjust the literature values by subtracting the GDR3 values determined by Lunz et al. (2023) for a meaningful comparison.Notably, due to the absence of literature calibration for the frame rotation of TYC and GDR1, our comparison in Fig. 3 focuses solely on calibration parameters for GDR2 and Hipparcos.
In the left panel, it is evident that the values of ⃗ β GDR23 obtained through analyses of NST, NSB, and SS calibration sources align well with those from previous studies within a 2σ range.Importantly, our parameters exhibit higher precision compared to previous studies, underscoring the robustness of our calibration methodology.Noteworthy consistency is observed between NST and NSB, as expected since both correct GDR3 astrometry to the binary barycenter albeit with different methods.While NST (or NSB) and SS yield similar values for ϵz, ωx, ωy, ωz, and δϖ, they diverge in ϵx and ϵy, indicating a potential dependence of frame rotation on calibration sources.
The right panel of Fig. 3 reveals significant scatter in ⃗ ω and δϖ for NSTH1, NSTH2, NSBH1, NSBH2, SSH1, and SSH2 in our study.In contrast, Lindegren et al. (2016) provide more precise values for these parameters by comparing Hipparcos with VLBI observations of 262 radio stars.Additionally, Brandt (2018) presents similar values for ⃗ ω by comparing the Hipparcos catalog with long-term proper motion derived from the positional difference between Hipparcos and Gaia, although the fitting does not include ⃗ ϵ.

Dependence of calibration parameters on magnitude and color
The relationship between calibration parameters in Gaia Data Releases (DRs) and stellar magnitude and color has been explored in previous studies (e.g., Lindegren et al. 2021a andCantat-Gaudin &Brandt 2021).Given the notably precise determination of calibration parameters for GDR2 relative to GDR3 ( ⃗ β GDR23 ) compared to other catalogs, our investigation focuses exclusively on exploring the dependence of ⃗ β GDR23 on the G magnitude and BP-RP color of Gaia sources.This dependence is illustrated for NST and SS calibration sources in Fig. 4.
The top two panels exhibit the correlation between calibration parameters and G magnitude.Notably, the faint NST calibration sources display more pronounced parameter variations than their bright SS counterparts.A distinct shift in GDR2-GDR3 frame offsets is observed around G = 13 mag, a phenomenon also noted by Brandt (2018), Lindegren et al. (2018), andCantat-Gaudin &Brandt (2021).This abrupt change is likely attributed to Hipparcos being unable to serve as a reference catalog for    calibrating GDR3 for targets with 11 < G < 13 mag (Cantat-Gaudin & Brandt 2021).Specifically, we identify a frame offset of ⃗ ϵ GDR23 = (0.14, 0.19, 0.03) ± (0.01, 0.02, 0.02) mas between the faint (G> 13 mag) and bright (G< 13 mag) frames of GDR2 relative to GDR3.However, the bright frame does not exhibit a sudden change in the rate of frame rotation (i.e., ⃗ ω) relative to the faint frame of GDR2, despite continuous variation across the G magnitude range.In comparison to NST calibration sources, the SS sources do not demonstrate a significant dependence of differential calibration on G magnitude.
The lower panels of Fig. 4 illustrate the color-dependent behavior of the differential calibration parameters for both SS and NST sources.Specifically, the color-dependent variations for SS and NST are respectively (0.02, 0.02, 0.02, 0.03, 0.02, 0.02, 0.03) mas (or mas yr −1 ) and (0.04, 0.05, 0.02, 0.02, 0.03, 0.02, 0.01) mas (or mas yr −1 ), exhibiting a significance level of approximately 3σ.Notably, we observe a consistent monotonic decrease in the zero-point parallax (δϖ) with BP-RP for the SS sources, aligning with the color-dependent trend of δϖ depicted in the top-left panel of Lindegren et al. (2021a).
Following a thorough comparison of various calibration results, we offer the following differential calibration recommendations relative to GDR3 for different catalogs: • For straightforward calibration of GDR2 data, employ SS-based values for sources with G <10.5 mag and NST-based values for sources with G >10.5 mag.
• For precise calibration of GDR2 data, utilize the Python scripts provided in the appendix to address bias on a case-by-case basis.
• Calibrate the 2-parameter GDR1 solutions using parameters determined with NSTP2.
• For GDR1 targets with 5-parameter solutions, use the SSP5 values to calibrate G < 10.5 mag stars and use the NSTP2 values to calibrate G > 10.5 mag stars.
• Use TYC positions without correction, in conjunction with GDR2 and GDR3, to constrain orbits.
• Calibrate HIP1 and HIP2 data using either non-zero ⃗ ϵ or non-zero ⃗ ω given by Lindegren et al. (2016) but not both, the ICRF2 could be transformed to GDR3 frame using the calibration parameters given by Charlot et al. (2020); Lunz et al. (2023).
• Quadratically add an astrometric jitter of 2.16 mas (or mas/yr) to HIP2 astrometry.These recommendations are automatically implemented in the Python script provided in the appendix.Investigating potential biases in the SS astrometric solutions arising from undetected binaries, we examine the impact of the RUWE parameter on calibration parameters.RUWE serves as an indicator of binarity or significant excess noise.Similar to the approach outlined in section 4.2, we categorize the SS calibration sources into subsamples based on different RUWE values.
The dependence of calibration parameters on RUWE is visually represented in Fig. 5. Remarkably, the parameters derived from the SS subsamples demonstrate minimal sensitivity to RUWE, indicating that our SS-based calibration remains largely unbiased even in the presence of undetected binaries.However, the uncertainty associated with calibration parameters increases, particularly for RUWE values exceeding 1.2.This outcome aligns with expectations, as sources with higher RUWEs or elevated astrometric jitters contribute substantial uncertainty to the calibration parameters.This reinforces the robustness of our approach as an effective calibration method.

CONCLUSION
In order to identify dark companions and binaries across multiple Gaia DRs, we have developed an astrometric modeling framework, designed for the combined analysis of Gaia, Hipparcos, and Tycho-2 catalogs.To address biases inherent in these catalogs, we ascertain error inflation, astrometric jitter, differential frame rotation, and parallax zero-point relative to GDR3.This involves utilizing calibration sources selected from GDR3, both with and without orbital solutions.Through the simultaneous fitting of calibration parameters and barycentric astrometry to the calibration sources, our analysis reveals negligible error inflation across all catalogs.Notably, a significant jitter of 2.16 mas is identified for HIP2 IAD, a frame offset ranging from 0.12 to 0.26 mas, and a frame rotation of 0.07-0.15mas yr −1 for GDR2 relative to GDR3.Additionally, a substantial frame offset of approximately -0.12 mas along the y-axis is observed between G1P2 and GDR3.While our estimation of calibration parameters for HIP1 and HIP2 aligns with previous studies, it's worth noting that the precision of values determined in this work falls short compared to those derived from VLBI observations of radio stars, as demonstrated by Lindegren et al. (2016).
Given the higher precision with which we determine the differential calibration parameters for GDR2 in this study compared to prior research, we delve deeper into investigating the dependency of these parameters on G-magnitude and BP-RP color.Our findings reveal a magnitude-dependent frame offset of 0.08 ± 0.02 mas and frame rotation of 0.26 ± 0.02 mas yr −1 for faint stars (G> 11 mag) in GDR2.Additionally, we identify a noteworthy 0.24 mas offset between bright (G< 13 mag) and faint (G> 13 mag) frames, aligning closely with findings from earlier studies such as Brandt (2018), Lindegren et al. (2018), andCantat-Gaudin &Brandt (2021).Additionally, a significant dependence of zero-point parallax on color is found in the SS calibration sources, consistent with previous studies by Lindegren et al. (2021a).A pipeline to calculate the magnitude-andcolor dependence of calibration parameters is provided in the appendix.
After comparing the calibration conducted in this study with previous research, we propose the utilization of NST and SS values for the global calibration of GDR2.For the calibration of G1P2 and G1P5, we suggest employing NSTP2 values.Additionally, we recommend adopting the calibration values provided by Lindegren et al. (2016) and Lunz et al. (2023) for the calibration of HIP1 and HIP2.Between the two versions of Hipparcos data, we advocate for the use of HIP2, incorporating an astrometric jitter of 2.16 mas.To mitigate unknown systematics in constraining stellar orbits, our recommendation is to utilize TYC data without correction.Furthermore, for GDR1 sources, we propose using the position instead of five-parameter astrometry.
While we conduct global calibrations for the most widely used astrometric catalogs, our current focus is solely on the examination of the dependence of differential calibration parameters on magnitude and color specifically for GDR2.To enhance the a priori calibration of astrometric bias, future efforts should delve into a detailed investigation of how systematics vary with additional parameters such as position and distance across various catalogs.However, it is important to note that our ability to develop a comprehensive calibration function for different catalogs is constrained.Therefore, adopting an approach of analyzing astrometric data on a case-by-case basis is likely the optimal solution for bias mitigation, as demonstrated by studies such as Snellen &Brown 2018 andFeng et al. 2019.Nonetheless, our work ensures a robust calibration of the most prevalent all-sky astrometric catalogs, facilitating the efficient detection of dark companions and binaries in extensive datasets using our algorithm.
Feng et al.

B. PIPELINE FOR DOWNLOADING DATA
The data used in this work is downloaded using several scripts, which are available in the GitHub repository: https://github.com/ruiyicheng/Download HIP Gaia GOST.The description of these scripts is as follows: • get hipIAD1997.pyThis code is dedicated to downloading IAD of the Hipparcos 1997 reductions.There are a total of 118,204 entries which have non-empty Hipparcos IAD information in the 1997 version.Using the get_hipIAD1997  function, the Hipparcos IAD of the corresponded source can be downloaded as a csv file.If the input HIP entry is non-empty, this csv file will contain columns of orbit number, source of abscissa (FAST or NDAC), partial derivatives of the abscissa with respect to five astrometric parameters (α * = α cos δ, δ, π, µα * , µ δ ), abscissa residual in mas, standard error of the abscissa in mas, correlation coefficient between abscissae, reference great circle mid-epoch in years, reference great circle mid-epoch in days, RA of the great circle pole in degrees, decl. of the great circle pole in degrees.Otherwise, The Hipparcos IAD of this HIP entry cannot be found will be printed.
• nssGDR123HIPTYC.sqlThis script is used to queue the Gaia data via ADQL (https://gea.esac.esa.int/archive/).It retrieves the five-parameter astrometric solution of Gaia DR1, DR2, and DR3 for non-single stars with two-body orbital solutions.Additionally, it obtains the cross-match results of these stars with the Hipparcos and Tycho catalogues.
• Obtain GOST.pyThis script is used to queue the Gaia GOST data using the interface in https://gaia.esac.esa.int/gost/GostServlet.It employs a similar pre-processing method as described in Brandt et al. (2021) to convert the server's returned results into a human-friendly pandas dataframe.The data is then filtered to match the GOST results obtained through the web page interface at https://gaia.esac.esa.int/gost/.The script requires the path of the csv file obtained using nssGDR123HIPTYC.sqlas an argument.The results are recorded in separate files named using the Gaia DR3 IDs.
• frame rotation correction.pyThis script is used to align the astrometry of different cataloges with Gaia DR3.The script requires the catalog astrometry, together with source and CAT1 as shown in Table 1.The results are recorded in a dictionary.An exception value would be given if there is no input proper motion or parallax.
• download single target.pyThis script offers a user-friendly interface for downloading and collaborating astrometric data.Users can queue the data using HIP id, TYC id, or Gaia DR1, DR2, DR3 source id.The Gaia astrometric data is acquired using an ADQL script similar to nssGDR123HIPTYC.sql.The HIP and TYC astrometry is gathered through cross-matching with the bulk data obtained from VizieR6 .The collaborated astrometry is obtained using frame rotation correction.py with the recommended startegy in section 4.2.Users have the option to download either the Hipparcos epoch data or Gaia GOST data for their target.The results are stored in the /results folder.

C. SYMBOLS AND ACRONYMS
The following table provides abbreviations used in the paper along with the meaning of acronyms and symbols used.
Table 2. Symbols used in this work.Note-Some of the variations of symbols are not shown and the meaning of symbol variations could be found in the corresponding text.

Symbols or Acronyms Meaning
where ∆tij = ti − tj, α b * j ≡ α b j / cos δ b j , δ b j , µ b α,j , µ b δ,j are respectively R.A, Decl., and proper motions in R.A. and Decl.at epoch tj.Here µα ≡ µα * = α cos δ.In this work, we only model the deviation from the motion determined by the reference astrometry.Hence the TSB relative motion is ∆α b The asymmetry in the uncertainty of parameters a and b is attributable to the restricted bin size of 0.02 employed during the sampling process.† † The first column displays the differential rotation and parallax zero-point between the RF0 and RF1 frames.For example, RF0-RF1 means that the astrometric offset induced by ⃗ β d should be added to the astrometry in the RF1 frame to derive the astrometry in the RF0 frame.In this work, the data (⃗ y) in the original reference frame of CAT1 is subtracted by κ d ⃗ β d to derive the astrometry in the GDR3 frame, which can be represented by CAT1-GDR3.a GDR3 Using optically bright radio stars observed by the very long baseline interferometry (VLBI), Lunz et al. (2023) measure the rotation of GDR3 frame (equivalent to GEDR3 frame) relative to the Third International Celestial Reference Frame (ICRF3; Charlot et al. 2020) at reference epoch J2016.0.b Lindegren et al. (2021b) conduct an ad hoc calibration using the difference between Hipparcos and GDR3 proper motions.c This calibration is conducted by Lindegren (2020) using the VLBI observations of 26 radio stars.Hence the parameters shown in this row measure the frame rotation of GDR2 bright sources relative to the frame determined by the 26 radio stars at epoch J2015.5.d The parameters are determined by Brandt (2018) by mixing HIP2 and HIP1 with a ratio of 60/40 (HIP12) and calibrating the proper motions of HIP12 and GDR2 using their positional differences (HG2).They claim negligible uncertainties on the parameters.e The value is obtained from table 1 of Fabricius et al. (2021).f The parameters measures the frame rotation of the original Hipparcos catalog relative to the ICRF2 (Fey et al. 2015) realized by radio observation of 262 sources Lindegren et al. (2016) at epoch J1991.25.However, ⃗ ϵ and ⃗ ω cannot be used together to correct for bias in the Hipparcos data because they are calculated under different assumptions.g Brandt et al. (2023) calibrate the HIP2 IAD using Hipparcos-GEDR3 (HG3) long-term proper motions and parallax differences.Note-Note that (a, b) for NST, NSB, and SS represent the error inflation and jitter for GDR2 while (a, b) for other calibration sources is for the corresponding CAT1.

Figure 1 .
Figure 1.Distribution of lnBF over the error inflation and jitter of GDR2 (first row), GDR3 (second row), G1P2 (third row) and G1P5 (fourth row) based on analyses of various calibration sources.The color encodes the value of lnBF of the calibration model relative to the model without calibration.The panels in the first and second rows are zoom-in versions of distributions over a bigger grid with 1 ≤ a ≤ 2 and 0 ≤ b ≤ 2, shown in the lower panels.

Figure 2 .
Figure 2. Distribution of calibration parameters for GDR2 and G1P5 based on 10,000 draws of 100 samples from various calibration sources.The parameter range is divided into 50 bins.The red dashed lines represent the modes of the distributions.The upper panels show the Pearson correlation value (Corr) and its significance indicated by the number of star-symbols.The p-values of less than 0.1, 0.05, 0.01, and 0.001 are respectively indicated by ".", "*", "**", and "***".No symbol is shown if the p-value is higher than 0.1.

Figure 3 .
Figure 3.Comparison of differential calibration parameters determined in this work and in literature.The left panel shows the values of β GDR23 while the right panel shows ⃗ β d of HIP1 and HIP2 relative to GDR3.Because Brandt (2018) determines the calibration parameters for a hybrid of HIP1 and HIP2, the same B18 values are presented for both HIP1 and HIP2.

Figure 4 .
Figure 4. Dependence of differential calibration parameters ( ⃗ β GDR23 ) on G magnitude (top) and BP-RP color (bottom) for the SS (left) and NST (right) calibration sources.Different colors represent different parameters.

Figure 5 .
Figure 5. Dependence of differential calibration parameters ( ⃗ β GDR23 ) on GDR3 RUWE for the SS calibration sources.Different colors represent different parameters.

Figure 7 .
Figure 7. Similar to Fig. 2, but for distribution of optimized calibration parameters for G1P2 (left) and TYC (right) based on 10,000 draws of 100 samples from various catalibration sources.

Figure 8 .
Figure 8. Similar to Fig. 2, but for distribution of optimized calibration parameters for HIP2 (left) and HIP1 (right) based on 10,000 draws of 100 samples from various catalibration sources.
axis of the photocenter relative to the mass center a b Binary semi-major axis ar Semi-major axis of the reflex motion of the target Au 1 au A, B, F , G Thiele Innes constants A ′ , B ′ , F ′ , G ′ Scaled Thiele Innes constants ⃗ b Vector of dependent variable in the normal equation η ⃗ β .offset, equivalent to ∆α * ∆ti ti − tj ∆tij Time difference between epoch ti and epoch tj ⃗ ϵ ≡ (ϵx, ϵy, ϵz) Offset between two reference frames at a reference epoch

Table 1 .
Lindegren et al. (2021a)While NST-based and NSB-based calibrations yield similar values for δ ⃗ β GDR23 , SS-based calibration results in slightly different values, suggesting a dependency of ⃗ β GDR23 on stellar parameters, as briefly mentioned inLindegren et al. (2021a).A detailed investigation of this dependence is presented in subsection 4.2.Calibration parameters determined for different calibration sources.

Table 1
continued on next page