© The Authors 2023

1 Introduction

The study of binaries is a constantly expanding field of research and a combination of multiple observational methods is frequently used for their characterisation because this allows the masses of each component to be determined directly. Absolute astrometry has long been used to identify and study invisible stellar, and more recently planetary, companions to stars. However, until now it has only been applied to relatively small samples. In this respect, as in many others, the ongoing Gaia mission (Gaia Collaboration 2016b) represents a game changer, and is expected to astrometrically detect thousands of planets and stellar companions. Nevertheless, its epoch data will only become available with the fourth Gaia data release¹. Until then, the data of Gaia’s predecessor, HIPPARCOS (ESA 1997), provide a very valuable test bed, because both missions similarly obtain astrometry along a single direction. Moreover, combining HIPPARCOS and Gaia data extends the period of measurement to over three decades, which is essential for longer-period binaries.

Several methods have already been used to mine the multiplicity information brought by the combination of HIPPARCOS and Gaia. The earliest was the comparison of the ‘short-term’ proper motions returned by both missions with the ‘long-term’ proper motion derived from the difference between the HIPPARCOS and Gaia positions (see Kervella et al. 2019, for a review of this method). This produced both astrometric acceleration catalogues (Brandt 2018; Kervella et al. 2019, 2022; Brandt 2021) and approximate statistical determinations of companion masses (Kervella et al. 2019). More recently, Brandt et al. (2021a,b) developed tools that can be used to adjust Gaia information (Brandt et al. 2021b), radial velocities and relative astrometric data (Brandt et al. 2021b), and HIPPARCOS Intermediate Astrometric Data (IAD, ESA 1997, Vol. 1, Sect. 2.8). Much earlier, Söderhjelm (1999) combined the raw HIPPARCOS data from the original HIPPARCOS reduction (ESA 1997), referred to as Transit Data (TD, ESA 1997, Vol. 1, Sect. 2.9), with ground-based observations in order to adjust orbital elements, but no attempt has yet been made to combine those with Gaia information.

BINARYS (orBIt determiNAtion with Absolute and Relative astrometRY and Spectroscopy) is our new tool and can be used to simultaneously adjust the residual abscissae from HIPPARCOS data (IAD or TD); the astrometric parameters available from Gaia; and complementary observations from relative astrometry and radial velocity. BINARYS uses a gradient descent method implementing automatic differentiation thanks to the R package TMB (Kristensen et al. 2016), and rigorously uses the information from HIPPARCOS and Gaia with minimal assumptions or simplifications.

In the following, we first present the data classes used by BINARYS in Sect. 2, and then the tool itself in Sect. 3 along with its limits and validation. In Sect. 4, we illustrate the capabilities of BINARYS with three systems: G1 494 where we combine relative astrometry with the HIPPARCOS IAD and Gaia data; GJ 2060, which we analyse with relative astrometry and HIPPARCOS TD; and HIP 88745, where we combine HIPPARCOS TD with Gaia-resolved observations.

2 The data

2.1 HIPPARCOS

The HIPPARCOS mission (ESA 1997) was in operation from 1989 to 1993. The satellite scanned the sky continuously along great circles, and projected the image of pre-selected stars through an alternatively transparent and opaque grid with a grid step of s = 1207.4 mas, which modulated their light. The one-dimensional (1D) position of the object along the scanning great circle is therefore encoded into the observed phase of the corresponding quasi-periodic signal. Each star was observed during multiple satellite transits and the observations are published as residual abscissae (noted ∆y), which are the difference between the observed position of the star and the predicted position along the scanning circle for the published best model of the star. Two data reductions are available: the original reduction (ESA 1997) and the new reduction (van Leeuwen 2007b). The tool can handle both reductions, but in the following sections we only use data from the new reduction.

The residual abscissae are published in the Intermediate Astrometric Data (IAD)². We note that Brandt et al. (2021a) found an issue on the IAD which is that the astrometric solutions obtained from them are not exactly the same as the published ones when the number of observations (NOB) is lower than the number of residual records (NRES). We confirm this issue but it does not seem to be due to data corruption, as we do not find repeated sequences for the along-scan errors described in Brandt et al. (2021a). We checked that no source studied here is affected. The IAD can be used when the observed object is a point-source for HIPPARCOS. When the object is instead a resolved binary or multiple system, the observed phase no longer measures its photocentre but instead something specific to the HIPPARCOS scanning grid method, dubbed the Hippacentre (Martin et al. 1997). Using only the IAD, it is possible to use the Hippacentre to constrain the mass and the intensity ratio of the components, as shown by Martin et al. (1997). However, the TD contain the full signal modulation parameters and therefore provide more constraints on the HIPPARCOS observations of a resolved system (Quist & Lindegren 1999). Söderhjelm (1999) pioneered the use of the TD to derive the masses of visual binaries. While the TD are provided for only one-third of the sources in the original HIPPARCOS solution, in the reduction of van Leeuwen (2007b) the TDs are available for all HIPPARCOS stars. Those TD are extracted by the Java Tool³ from the HIPPARCOS calibrated raw data. An ASCII version of the TD is in preparation as well as an update of the tool to retrieve those data. The TD files provided in the DVD suffer from a factor 10 issue on the β₅ and associated error for some stars and these should therefore not be used. The new reduction did not reprocess the photometric signal, which should therefore be retrieved in the original reduction Epoch Photometric Annex and Extension (accessible via ESASky legacy TAP query⁴). Some transits do not have photometric information available in the original reduction and we ignore those transits in our code.

For the original reduction, BINARYS can handle either the IAD, which lists the residual abscissae relative to the published five-parameter astrometric model, or the TD for which we apply the method described in Quist & Lindegren (1999). For the new reduction, the residual abscissae are given relative to the model used for each star, which need not be the five-parameter model. BINARYS only uses the IAD when the solution for the star is a five-parameter solution, and reverts to using the TDs when it was analysed with a different model⁵.

2.2 Gaia EDR3

The Gaia mission (Gaia Collaboration 2016b) started observing in 2014 and is ongoing. The observation strategy of Gaia is similar to that of HIPPARCOS, except that the satellite records small images rather than a periodically modulated signal. The raw data are not available yet and the different data releases to date (DR1: Gaia Collaboration 2016a, DR2: Gaia Collaboration 2018 and EDR3: Gaia Collaboration 2021) only provide the two or five astrometric parameters that best match the observations, without taking into account a possible multiplicity. Multiplicity will be taken into account for the first time in the forthcoming DR3 release. The scanning law of the satellite, which contains the pointing direction and the scanning angle as a function of time, is also published⁶ and provides the conditions under which a given star was observed.

Before using the Gaia data, we also examine ancillary information such as the ruwe (Renormalised Unit Weight Error) and the multi peak flag (ipd_frac_multi_peak) present in EDR3 (Lindegren et al. 2021b). The ruwe evaluates the quality of the five-parameter solution, and a value above 1.4 indicates that the published solution may not describe the object well (Lindegren 2018). The multi-peak flag indicates the percentage of the windows used for the astrometric processing of the source which contain a double peak, and a high value is evidence of flux contamination. For us to use the Gaia data, the signal must originate from the source alone, because accounting for flux contamination would require a model of the line spread function fitting, which is not published at this point. When Gaia does not fully resolve the system and the secondary contributes non-negligible light, we cannot use the Gaia data, which correspond to separations of between 9 mas and 0.27″ depending on the magnitude difference (Lindegren 2022). For smaller separations, the photocentre can be used. If Gaia fully resolves the system and gives a separate solution for each component, then BINARYS can use those solutions, even when they are only two-parameter solutions, which ignore parallax and proper motion. Analysis of partially resolved systems with non-negligible light from companions will be undertaken when more detailed Gaia data are provided by the DR4 release.

To combine Gaia and HIPPARCOS data, we have to bring them into the same reference frame. We somewhat arbitrarily chose to convert the Gaia positions and proper motions to the HIPPARCOS proper motion reference frame (Lindegren et al. 2018, Brandt 2018, Kervella et al. 2019). As a consequence, the astrometric parameters that BINARYS adjusts to the data are in the HIPPARCOS reference frame at epoch: HIP_epoch = J1991.25.

The rotation to be used for Gaia EDR3 is (Fabricius et al. 2021), and the transformed Gaia astrometric parameters are given by:

(1)

with α* = α cos δ, , and the polar to Cartesian coordinates transformation matrix:

(2)

with Δ_GH being the difference between the Gaia (EDR3_epoch = J 2016) and HIP_epoch epochs. Similarly, we correct for the parallax zero-point difference: ϖ_new = ϖ − ϖ^shift, with mas (Fabricius et al. 2021) after applying the Gaia EDR3 parallax correction proposed by Lindegren et al. (2021a). The uncertainties on the five astrometric parameters are inflated according to the parallax error under-estimation factor derived by El-Badry et al. (2021).

2.3 Relative astrometry

Relative astrometry data can originate from either direct imaging or interferometric observations and consist of relative positions of the components at one or several epochs. The inputs for BINA-RYS are the date of the observation and the relative position of the two components in α and δ direction (ξ and η) with their associated uncertainties (σ_ξ and σ_η).

When the relative positions are published as a separation and position angle (ρ and θ) and unless the publication includes a full covariance matrix, we adopt as covariance matrix Σ_ξη:

(3)

with J being the Jacobian of the polar to Cartesian transformation and assuming that cor(ρ, θ) is null.

2.4 Radial velocity

The radial velocity (RV) inputs contain the date of the observation, the RV, its uncertainty, and optionally a code for the instrument that was used. The latter is needed when the radial velocity inputs were obtained with multiple instruments, and allows offsets to be adjusted to account for RV zero-point differences. A jitter can also be added, either to increase the instrument noise or to take into account an unmodelled stellar variability. The radial velocities can be adjusted for either the primary or the secondary stars, allowing both single-lined and double-lined binaries to be handled.

3 Method: Combination of absolute astrometry with relative astrometry and radial velocity

In the following, we describe how BINARYS estimates the orbital and astrometric parameters (OPs and APs) of a binary system. The adjusted OPs are expressed using the Campbell elements , where X can be either {a₂₁} or {M₁}; P is the period in years, defines the epoch of one periastron, counted from J2000.0 in units of the orbital period ; a₁ and a₂₁ are the semi-major axis of respectively the orbit of the primary and the relative orbit, in au; e is the eccentricity, and i is the inclination of the orbital plane to the tangent plane of the sky, oriented with the convention that 0 ≤ i ≤ 90° for a direct (defined by an increasing positional angle that is counted positive from north towards east direction) apparent motion and 90 ≤ i ≤ 180^° for a retrograde apparent motion. Also, ω₁ is the argument of periastron of the primary, counted from the ascending node and in the direction of the motion. The argument of periastron of the secondary is linked with that of the primary by ω₂ = ω₁ + π. Furthermore, Ω is the position angle of the ascending node, with the conventions used in HIPPARCOS and Gaia: it is the position angle – counted counterclockwise from the δ direction – of the intersection of the orbital and tangent planes. When radial velocities are available and resolve the ambiguity between the two nodes, Ω corresponds to the node where the primary star recedes from the observer; otherwise, we arbitrarily impose 0 ≤ Ω ≤ 180°. Finally, M₁ is the mass of the primary in units of solar masses. When combining absolute astrometry with relative astrometry, fitting a₂₁ is a natural choice, while combining absolute astrometry with radial velocities is easier using M₁, a parameter which has the advantage of having spectroscopic and/or photometric estimates. Both parameters lead in practice to the mass ratio information q = M₂/M₁ through the equations:

(4)

To handle the photocentre motion, BINARYS also adjusts the fractional luminosity . The adjusted astrometric parameters are the usual .

3.1 Adjustment of relative astrometry data

For a given observation time t, the positions of the primary (1) and secondary (2) stars relative to the barycentre and along the α and δ direction (ξ and η) are computed as:

(5)

with the polar coordinates of the primary on its orbit and υ the true anomaly. Here, υ is related to the eccentric anomaly E through

(6)

and E is obtained by numerically solving Kepler’s equation,

(7)

over ten iterations (Heintz 1978). The relative positions between the two stars in au are obtained by calculating Δξ = ξ₂ − ξ₁ and Δη = η₂ − η₁, and these converted to angular separations by multiplying with the parallax (which is one of the astrometric parameters). We finally compute the residuals between the computed and observed relative positions.

3.2 Adjustment of radial velocity data

To predict the radial velocity, we first calculate from the orbital parameters the semi amplitude K (km s⁻¹) of the radial velocity signal as, for the primary and secondary:

(8)

with C = 29.78525 km and AU being the astronomical unit in kilometres. The predicted radial velocities of the primary and the secondary for a given epoch are then:

(9)

with RV₀ being the radial velocity of the barycentre to be adjusted. The predicted radial velocity is then compared to the observed radial velocity data.

3.3 Adjustment of HIPPARCOS data

To predict the residual abscissae Δv and compare them to the observed ones at each observing epoch, we must first project the separation of the two stars onto the HIPPARCOS scanning grid for the orientation of the grid at that epoch. The separations of the components along α and δ are the Δξ and Δη calculated in Sect. 3.1. To project them on the grid, we resort to the partial derivatives of the abscissa against the five astrometric parameters with , which are:

(10)

with ψ the position angle of the scanning direction, ϖ_factor the parallax factor (ESA 1997, Eq. (1.2.26)), and ΔT the observation epoch relative to HIP_epoch. These are available in the IAD and the TD of the new HIPPARCOS reduction, while the original reduction directly provides the five partial derivatives.

For each transit and for a given set of orbital parameters θ_OP, the projected separation ρ_p on the scanning direction is:

(11)

We then calculate ζ, the projected separation in units of the HIPPARCOS grid step:

(12)

Finally, we calculate Δv_B, the position shift along the scanning direction due to the binary motion:

(13)

with B the fractional mass, or in other words the mass of the secondary divided by the total mass of the system , and the phase ϕ = a tan 2(β sin ζ, 1 − β + β cos ζ) where β is the fractional luminosity . For a detailed explanation of these latter steps, see Martin et al. (1997), who used the function Angle(x, y), which is equivalent to atan2(y, x). If the flux of the secondary is negligible, β becomes zero and the Δv_B shifts are purely from the reflex orbit of the primary star. If the secondary star contributes light but the separation remains small compared to the grid step, the orbit is that of the photocentre.

We then calculate Δv by applying the partial derivatives to the difference between the five adjusted astrometric parameters (θ_AP) and the five published reference astrometric parameters :

(14)

These new computed residuals are then compared to the observed ones for each transit.

When neither the secondary flux nor the separation is negligible, we include three additional observational quantities in the adjustment, which are available in the HIPPARCOS transit data (van Leeuwen 2007a): β₄, β₅, and H_ACH_DC = Hp_ac − Hp_dc. Here, β₄ and β₅ describe the amplitude and the phase of the second harmonic of the grid-modulated signal and are closely related to an interferometric visibility measured at the corresponding angular frequency. Hp_dc and Hp_ac are magnitudes (in the Hp HIPPARCOS spectral bandpass) evaluated from the unmodulated DC and the modulated AC components of the HIPPARCOS transit signal, respectively. These three observable quantities can be computed from the adjustable parameters as:

(15)

where r is the ratio of the secondary and primary luminosities, which can be written as r = β/(1 − β), and nf = (1 + 2r cos ζ + r²)^3/2.

The uncertainties on Δv are modified according to the amplitude of the first harmonic of the modulated signal, which for a binary is decreased from its point-source value by a factor . When comparing the computed and observed Δv, we consequently have to increase the measurement errors σ_Δv by the factor .

When one of the components is variable, we adjust one value of the fractional luminosity β for each epoch rather than one common value. We note that β is then calculated as , with r = r₁ when the primary is variable, and r = r₂ when the secondary is variable:

(16)

where β₀ and are the mean fractional luminosity and total magnitude, respectively (van Leeuwen 2007b).

3.4 Adjustment of Gaia data

The Gaia observations also constrain the astrometric and orbital information, but despite the impressive size of the Gaia DR3 non-single star catalogue (Gaia Collaboration 2023), for the majority of the systems, we only have access to the set of five astrometric parameters, because those systems are not yet treated as binaries by Gaia DR3. For those systems, the Gaia orbital information is therefore encoded in the bias of the five published astrometric parameters away from their true barycentric values. We therefore proceed by computing what Gaia would have observed for a given set of orbital and barycentric astrometric parameters and comparing to the five ‘effective’ astrometric parameters published in the [E]DR3 catalogue.

To do that, we start by propagating the astrometric parameters of the barycentre to the mean epoch of Gaia DR3. For nearby stars, that transformation must take the radial velocity into account (ESA 1997: Sect. 1.5). Perspective acceleration during the HIPPARCOS and Gaia observations, which is only needed for the closest stars, is not yet taken into account. From the published Gaia scanning law, we then retrieve the epochs when one of the two Gaia fields of view passed over the target of interest, as well as the scanning angle ψ for each of those epochs. We assume in the following that all those epochs have contributed to the Gaia solution.

We proceed to compute, for each epoch, the orbital motion projected along the Gaia scanning direction. When Gaia does not resolve the system, that motion is that of the photocentre. If instead Gaia resolved the system and gives separate solutions for the two components, we compute individual offsets for the primary and the secondary stars. The positions relative to the barycentre for the photocentre (0), the primary (1), and the secondary (2) projected along the Gaia scanning direction are respectively:

(17)

The first equation in Eq. (17) for Gaia (also writable as is similar to Eq. (13) for HIPPARCOS except that the extra complication in the HIPPARCOS formulation comes from the fact that Eq. (13) corresponds to the Hippacentre instead of the photocentre because of the signal modulation (Martin et al. 1997). From those abscissa residuals along the Gaia scanning direction for each observing epoch and the astrometric parameters propagated to the mean Gaia epoch , we estimate the five astrometric parameters that Gaia would have observed (AP^G) and compare those with the published parameters:

(18)

where X is the variation in astrometric parameters that reflects the residuals due to binarity. X is obtained by solving the linear equation: R = D · X, where R is the matrix of residuals Δv_B and D is the matrix of partial derivatives (van Leeuwen & Evans 1998, Sect. 3.1). The partial derivative along-scan ϖ_factor is computed using the position of the system on the sky, the observation epoch, and the corresponding scanning angle, as well as the orbit of the Earth (ESA 1997, Eq. (1.2.26)). For targets where Gaia published only positions and no proper motion or parallax, we nonetheless compute all five astrometric parameters and then discard the parallax and the proper motion. In the Gaia processing, Galactic prior information is added to provide more realistic uncertainties, step that does not need to be reproduced here.

3.5 Template Model Builder: Source code, options, and limitations of the tool

Template Model Builder (TMB; Kristensen et al. 2016) is an open-source R package designed to quickly and robustly adjust non-linear models with a large number of parameters. The R code calls functions from a user-provided C++ file that compute the likelihoods, which we make available⁷. To briefly explore the uncertainties and degeneracies of the parameters, we post-process the TMB results with a short MCMC run using the companion R package tmbstan (Monnahan & Kristensen 2018). In the present paper, the orbits are well constrained and a single short MCMC chain of 3000 iterations – of which we discard the first 1500 as warm-up iterations – has been found to be enough to reach convergence.

To help TMB converge, we adopt starting values from the literature whenever available. For previously unstudied systems, we explore a large range of starting values and often initially fix some parameters to plausible values (e.g. starting with a circular orbit or fixing the primary mass). When adjusting to HIPPARCOS transit data, good starting values for the astrometric parameters that already take into account a preliminary orbit greatly help. An option to ignore for the first few iterations can also help to quickly obtain starting values for the astrometric parameters. When adjusting for a radial velocity jitter, its value is best determined through an MCMC run – which is then set fixed when running TMB – because MCMC is less disturbed by jitter than gradient descent algorithms. When the flux ratio of the two stars is available, whether from Gaia or from ancillary observations through similar filters, it can enter the adjustment as an observation with its uncertainty.

The TMB adjustment works for well-constrained orbits for which there exists sufficient available data. The tool takes into account the system’s perspective acceleration between the HIPPARCOS and Gaia epochs, but not along the Gaia mission as needed for very nearby and/or fast-moving stars. Also, for stars identified by HIPPARCOS as having a component solution, it is possible that the light that has contaminated the data does not come from a companion of the system but from another star, as can happen in clusters. Both issues are a matter for future developments.

3.6 The evaluation of the solution

We evaluate each adjustment through its goodness of fit F2 (Wilson & Hilferty 1931), which asymptotically follows a Gaussian distribution and which is defined as:

(19)

with k being the number of degrees of freedom (the number of observations minus the number of adjusted parameters) and the sum of the χ² contributions of the individual observational methods. To be qualified as good, the adjustment must have F2 below 3.

To test the improvement of our solution on the HIPPARCOS data, we compute, as in HIPPARCOS, an F2 using only the χ² contribution associated with the residual abscissa (with their uncertainties increased by the multiplicative factor described in Sect. 3.3), and the number of parameters adjusted for the published solution. For the global F2 of our adjustments, β₄, β₅, and H_ACH_DC also contribute to the HIPPARCOS χ².

4 Orbital study of benchmark systems

For illustration, we choose three binary systems, for which we use different data-type combinations and which also provide interesting astrophysical results. The Gl 494 system tests the combination of HIPPARCOS IAD with Gaia astrometric parameters and relative astrometry; GJ 2060 tests the combination of HIPPARCOS TD with relative astrometry (including a new GRAVITY observation, Table 1); and HIP 88745 tests the combination of HIPPARCOS TD with Gaia resolved observations. For these three stars, radial velocity data are also available but were not included in the adjustment: we only used them for independent verification of the results of the adjustments.

4.1 Gl 494: combination of absolute Intermediate Astrometric Data, Gaia, and direct imaging

The Gl 494 (HIP 63510, Ross 458) system around an M star contains a close binary, which was first detected astrometrically by Heintz (1994) and then resolved with adaptive-optics imaging (Beuzit et al. 2004; Mann et al. 2019; Bowler et al. 2020), as well as possibly a common proper-motion planetary mass T-dwarf (Gl 494c) at a projected distance of 1200 au (Goldman et al. 2010 and Scholz 2010). We study the inner pair using 16 relative positions compiled by Bowler et al. (2020) and obtained with adaptive-optics imagers PUEO on CFHT, NACO on the VLT, and NIRC2 on the Keck telescope, which cover the binary orbit well. There are also 65 radial velocity measurements from the HIRES spectrograph (Tal-Or et al. 2019), but these are highly impacted by the intrinsic variability of the magnetically active primary star: we only use them for validation purposes as well as to identify which of the two nodes of the orbit is the ascending one.

The K band contrast (Δm_K = 4.27 ± 0.02, Mann et al. 2019) and the much redder spectrum of the late-M secondary guarantee that the secondary star contributes negligible light in the HIPPARCOS and Gaia observing bands; we consequently neglect any light from the secondary in the following. The Gaia EDR3 information is indeed compatible with an unresolved source: while the ruwe value is a very high 4.19, the multi peak rate is low, and the Gaia signal is therefore compatible with the reflex motion of only one luminous star. We can therefore safely use the published astrometric parameters as representing the average motion of the primary during the first 34 months of the Gaia mission. The solution published in the original reduction of HIPPARCOS includes an acceleration (seven-parameter solution), and the F2 goodness of fit for the five-parameter solution in the new reduction of HIPPARCOS is 2.32, meaning that the reflex motion was already detected by HIPPARCOS itself. Kervella et al. (2019) also detected a proper-motion anomaly, therefore finding signal in the HIPPARCOS–Gaia difference.

We adjusted the orbital and astrometric parameters to the relative astrometric data, to the five Gaia astrometric parameters, and to the HIPPARCOS residual abscissae extracted from the IAD of the new reduction using the orbital parameters of Mann et al. (2019) as starting values for TMB. Figure 1 represents the relative astrometric observations together with the best orbit obtained with TMB, as well as 1500 MCMC orbit samples to illustrate the uncertainty. Figure 2 shows that this adjusted orbit also matches the (unused) radial velocity measurements well, which provides independent validation. The solution (Table 2) is also fully compatible with, but improves upon, the orbits published by Bowler et al. (2020) and Mann et al. (2019).

The goodness of fit of the TMB best solution is F2 = 3.84, and is dominated by two 3σ outliers amongst the relative astrometry observations; if we remove those two, the goodness of fit improves to F2 = 2.01. Our accounting for the reflex motion greatly improves the match to the HIPPARCOS residual abscissae, with a revised goodness-of-fit contribution of F2 = −0.17.

We note that, unsurprisingly for a system with both an orbital period that is approximately three times the length of the HIPPARCOS mission and a separation of several hundred milliarcseconds, and in agreement with the proper-motion anomalies previously detected (Makarov & Kaplan 2005; Frankowski et al. 2007; Kervella et al. 2019), the proper motion that we derive for the barycentre (Table 3) differs greatly from the published HIPPARCOS value ( mas yr⁻¹ and μ_δ = −13.6 ± 1.0 mas yr⁻¹) and is in full agreement with the long-term proper motion provided in Tycho-2 ( mas yr⁻¹, μ_δ = –25.1 ± 1.4 mas yr^–1, Høg et al. 2000). The revised proper motion is much less compatible than the HIPPARCOS value with the proper motion of the proposed third component C (Table 3 of Scholz 2010), with the χ² between the proper motions of AB and of c now corresponding to a p-value of 2.62 × 10⁻¹⁰ instead of 0.01. We conclude that Gl 494C does not co-move with Gl 494AB and is likely not gravitationally bound to it.

We also determine, for the first time, the masses of both components of Gl 494 purely from Newtonian physics and without having to adopt a mass of the primary from a mass–luminosity relation. Those mass values (also reported in Table 4) are M₁ = 0.584 ± 0.003 M_⊙ and M₂ = 87 ± 1 M_Jup, leading to a total mass for the system of M_tot = 0.667 ± 0.004 M_⊙, which is in agreement – with a smaller uncertainty – with the estimations given by Mann et al. (2019) (M_tot = 0.666 ± 0.035 M_⊙) and Bowler et al. (2020) (M_tot = 0.66 ± 0.02 M_⊙). The individual masses we derived are in agreement with the frequently used mass–luminosity relations of Delfosse et al. (2000) and Mann et al. (2019; Fig. 3), but the agreement is even better with the BT-Settl isochrones (Baraffe et al. 2015) in the age range of 150-800 Myr derived by Burgasser et al. (2010).

Fig. 1 Gl 494 best-fitting orbit from TMB (values in Table 2) as the green dotted line and sample orbits from the MCMC algorithm in black. The direct imaging observations (Bowler et al. 2020) with their associated error bars are in red.

Fig. 2 Radial velocity behaviour of Gl 494 predicted by the adjustment of the direct imaging, HIPPARCOS, and Gaia data. The best-fit orbit from TMB (values in Table 2) is the green dotted line and sample orbits from the MCMC are in black. The red dots represent the radial velocity observations with their associated error bars. The systemic velocity of 0.59 km s−1 was derived from a TMB adjustment of the radial velocities and the direct imaging data for the sole purpose of this visualisation.

4.2 GJ 2060: Transit data and direct imaging

The GJ 2060 (HIP 36349) M dwarf system is a member of the AB Doradus moving group, the study of which is essential for the age prediction of the group, which is not yet well constrained (100–150 Myr in the most recent studies: Barenfeld et al. 2013 and Bell et al. 2015). As such, GJ 2060 adds to the short list of young, tight binary systems amenable to dynamical measurements. It can be used to vet the evolutionary model predictions known to be impacted by several uncertainties at young ages and in the low-mass regime (see Mathieu et al. 2007, for a review).

Rodet et al. (2018) previously studied the system using direct imaging and radial velocity observations. We use the 17 known relative astrometric measurements gathered from multiple imaging instruments (VLT/NaCo, astralux, Gemini/NICI, VLT/SPHERE; Rodet et al. 2018), and one new and higher precision one obtained with the VLTI/GRAVITY instrument (Gravity Collaboration 2017). The point is reported in Table 1.

Rodet et al. (2018) also use ten radial velocity measurements from FEROS, which we choose not to use because a large jitter is present due to the stellar variability and the measurements are impacted by the flux of the secondary (~0.25 flux ratio in the FEROS bandpass). The spectra would ideally be reanalysed as double lined, but the velocities of the two components are not well separated at any of the FEROS epochs.

Both components of the system contribute significant flux in the HIPPARCOS (component solution), so we have to use the HIPPARCOS TD rather than the IAD. Gaia did not resolve the system in EDR3 – while both components are contributing to the signal (ipd_frac_multi_peak = 76) – so we cannot use the Gaia astrometric parameters. Whilst taking into account the flux of the secondary, the F2 = 2.07 in the new reduction of HIPPARCOS indicates a small astrometric signal.

For this system, we adjust the orbit to the relative astrometric observations and to the TD from the new HIPPARCOS reduction. The photometric variability (Messina et al. 2010) is of the same order of magnitude as the HIPPARCOS photometric data noise, and so we did not consider the variability of the primary star. We use the orbital solution of Rodet et al. (2018) as starting values for the TMB gradient descent. A full exploration of the parameter space was also tested, leading to the same solution. The solution is represented on the direct-imaging data in Fig. 4. Figure 5 shows that this solution is qualitatively consistent with the radial velocity data, which provides independent validation, and our orbital parameters (Table 2) are compatible with those of Rodet et al. (2018), with the same strong correlation between ω and Ω (Fig. B.2). The goodness of fit of the TMB best solution is F2 = 2.99, and the contribution to F2 of the TD is significantly better for our orbital solution (F2 = 0.76) than the published one (F2 = 2.07).

Rodet et al. (2018) reported the total system mass from the relative astrometry. Here, we directly determine the individual dynamical masses of both companions for the first time. The masses derived from our adjusted orbital parameters (Table 2) are M₁ = 0.61 ± 0.06 M_⊙ and (Table 4). The fractional mass deducted is and this is consistent with the estimation of Rodet et al. 2018 () from the SB2 assumption using the method proposed by Montet et al. (2015). Comparison of the masses and luminosities with stellar evolution models (Appendix A) points to an age greater than 100 Myr, which is consistent with the most recent age estimates for the AB Doradus moving group, which rely on kinematics and chemistry (Barenfeld et al. 2013), placement of group members on isochrones (Bell et al. 2015), or cosmochronology (Gagné et al. 2018).

Fig. 3 Gl 494 masses compared to the mass-luminosity relations of Delfosse et al. (2000) and Mann et al. (2019), and the BT-Settl isochrones (Baraffe et al. 2015).

Fig. 4 Orbits from the MCMC algorithm in black and of the best solution from TMB in green dotted line for GJ2060 (values in Table 2). The direct imaging observations with their associated error bars are in red and the new GRAVITY point with its error bars is in light blue.

Fig. 5 Radial velocity behaviour of GJ2060 predicted by the adjustment of direct imaging and HIPPARCOS data. The best solution from TMB is the green dotted line (values in Table 2) and the MCMC results are in black. The radial velocity observations are in red dots with the associated error bars and with RV0 = 28.8 km s−1 (Rodet et al. 2018).

Fig. 6 Orbits of HIP 88745 adjusted to the HIPPARCOS and Gaia data. The best-fit TMB solution is displayed as the green dotted line (values in Table 2, solution A), and a sampling of the MCMC solutions is shown in black. The direct imaging observations (in red) were not used in the adjustment and provide an independent validation of the orbit.

Fig. 7 Radial velocity behaviour of HIP 88745 predicted by the adjustment of the HIPPARCOS and Gaia data. The green dotted line represents the TMB best-fit solution (values in Table 2, solution A) and the black lines show the MCMC solutions. The radial velocity observations from SB9 are in red, the Gaia DR3 radial velocity is in orange, and the ELODIE archive is shown in blue. The systemic velocity of RV0 = 0.40 km s−1 was determined from an adjustment of the SB9 radial velocities together with TD and Gaia observations, which was performed for the sole purpose of this visualisation.

4.3 HIP 88745: HIPPARCOS transit data and resolved Gaia observation

HIP 88745 is known to be a binary system with a main sequence F star (Hutter et al. 2019) and a circumbinary polarised debris disk (Kennedy et al. 2012). Direct imaging studies of the system include Heintz (1972), Abt & Willmarth (2006), Kennedy et al. (2012), Malkov et al. (2012), and Jao et al. (2016). Söderhjelm (1999) studied the system using the TD of the original HIPPARCOS reduction, and it appears in his Table 4, which lists stars for which only the total mass of the system could be derived.

We chose to analyse this system because it has a component solution in HIPPARCOS and was analysed as resolved in Gaia EDR3. We therefore adjusted to the HIPPARCOS TD as well as to the Gaia EDR3 astrometric parameters for the primary (5AP_A) and secondary (2AP_B) components. The HIPPARCOS new reduction considered the secondary flux and an astrometric acceleration, but still contains a strong remaining signal (F2 = 9.77). The ruwe value of the 5AP_A Gaia EDR3 solution is 1.39.

Direct imaging and radial velocity data are also available for this star: 17 relative positions from the Fourth Catalog of Interferometric Measurements of Binary Stars (Hartkopf et al. 2004) and 19 radial velocities of the primary from SB9 Pourbaix et al. (2004) covering a small fraction of the orbital period. To complete those radial velocities, we also consider the Gaia DR3 radial velocity Katz et al. (2023) with the error inflation described in Babusiaux et al. (2023) and two observations available in the ELODIE archive⁸ (Baranne et al. 1996) – for which we consider a quite arbitrary 0.3 km s⁻¹ offset uncertainty – in our validation plots.

For this star, we performed two adjustments, one (A) to the absolute astrometry data from HIPPARCOS and Gaia, and one (AR) that additionally uses the available relative astrometry. The (A) adjustment tests what can be done with a pure absolute astrometry fit and provides parameters that are independent of the direct imaging and radial velocity, and the (AR) adjustment provides better constrained parameters. For both adjustments, we used the Jao et al. (2016) orbital parameters as starting values. Tables 2 and 3 present our astrometric and orbital parameters for both adjustments. Figures B.3 and B.4 show that, as expected, the correlations are reduced by the introduction of the direct imaging data in the fit. Figures 6 and 7 show that the orbital solution from the pure absolute astrometric adjustment matches the direct imaging and radial velocity data well, validating this solution. The fitted 3.81 ±0.02 mag contrast between the two components in the HIPPARCOS band is qualitatively consistent with the Gaia magnitude difference of Δm_G = 3.406 ± 0.005 mag given the bluer HIPPARCOS passband. Our adjusted parameters are compatible with Jao et al. (2016) within 3 σ.

The global F2 of the (A) and (AR) adjustments is strongly dominated by the HIPPARCOS TD, for which two 5σ outliers are removed. The global F2 values are F2 = 6.97 and 7.1 for the (A) and (AR) adjustments, respectively. The new F2 values for HIPPARCOS are F2 = 9.09 and 9.26 for the (A) and (AR) adjustments, respectively, which remain equivalent to the published F2 values that took into account an astrometric acceleration. We do not fully understand the reason for this high score, but a third component in the system is one possibility. One was previously listed in the Washington Double Star Catalog Mason et al. (2001) before being classified as a non-detection by Hutter et al. (2019). We note that both components are listed in Gaia EDR3 with a non-negligible multi-peak fraction (31 and 21% for the primary and the secondary). For the primary, those might be from transits where the secondary is not separately detected, as four times as many observations were used for the primary compared to the secondary.

We directly determine individual dynamical masses (Table 4) of both companions using the TD and Gaia astrometric parameters of the primary and secondary. The masses derived from the adjusted orbital parameters (Table 2) using only absolute astrometry are M₁ = 1.01 ± 0.04 M_⊙ and M₂ = 0.68 ± 0.04 M_⊙ and those using absolute astrometry and relative astrometry are M₁ = 0.96 ± 0.02 M_⊙ and . These values are at around 3σ from the estimations made by Jao et al. (2016) from the SB1 mass function and the relative and photocentre orbit, of namely M₁ = 0.89 ± 0.03 M_⊙ and M₂ = 0.51 ± 0.03 M_⊙. Figure 8 shows that our new masses are in closer agreement with the PARSEC isochrones (Bressan et al. 2012) for a metallicity of −0.6 dex (see references in Jao et al. 2016) than the Jao et al. (2016) masses. The primary star mass would lead to an age of the system of 5± 1.3 Gyr for a −0.6 dex metallicity.

Fig. 8 PARSEC isochrones with [M/H] = –0.6 dex colour-coded by age (in Gyr). The dark circles correspond to the masses derived in this work (AR solution) while the grey squares are from Jao et al. (2016).

5 Discussion and conclusion

We present our new BINARYS tool, which rigorously combines HIPPARCOS and Gaia observations of binary stars with relative astrometry and/or radial velocity observations. For systems where the secondary contributes significant light, BINARYS uses the raw HIPPARCOS transit data.

For illustration and validation, we present three systems studied with BINARYS. The adjustment of direct imaging, HIPPARCOS IAD, and Gaia EDR3 constrains the primary and secondary masses in the Gl 494 system to M₁ = 0.584 ± 0.003 M_⊙ and M₂ = 87 ± 1 M_Jup. That adjustment also indicates that Gl 494C is unlikely to co-move with Gl 494AB.

The adjustment of direct imaging and HIPPARCOS TD on the AB Doradus GJ 2060AB system determines the masses of its primary and secondary and , which in turn constrains the age of the system to older than 100 Myr, which is in good agreement with the most recent estimate of the moving group age.

Finally, the adjustment of HIPPARCOS TD and resolved Gaia observations of HIP 88745 gave masses for the primary and secondary of M₁ = 0.96 ± 0.02 M_⊙ and , with strong residuals in the HIPPARCOS TD. Those may reflect a potential remaining signal in the TD, which might become usable later with further information and which could be from a third component.

In the future, we plan to extend BINARYS to accommodate very nearby stars, which have significant perspective acceleration during the Gaia and HIPPARCOS missions, and stars in clusters, where light from a star outside the system can contaminate HIPPARCOS observations. BINARYS is also being extended for the study of triple systems (Lagrange et al. 2020) and to take into account the new non-single solutions (NSSs) that are provided by Gaia DR3 (Gaia Collaboration 2023).

Additionally, this tool prepares us for Gaia DR4, which will provide epoch observations. When DR4 is released, we will be able to combine HIPPARCOS and the Gaia equivalent of the TD. By that point, the tool will run mostly without HIPPARCOS constraints because of the huge sample size difference. Although non-single solutions will be provided by the Gaia-DPAC consortium, the combination of Gaia with external data will have to be done using the epoch data for an optimised solution, but also in order to derive solutions for systems with an insufficient Gaia signal to obtain a full NSS solution alone, and to handle specific cases such as multiple systems. The fine and accurate handling of the Gaia epoch data will be crucial for the study of the exo-planets expected to be discovered by Gaia, of which there may be up to ~70 000 for a ten-year mission (Perryman et al. 2014).

Acknowledgements

We thank the referee for his detailed comments that helped to improve the clarity of the manuscript. This work is supported by the French National Research Agency in the frame-work of the Investissements d’Avenir program (ANR-15-IDEX-02), in particular through the funding of the “Origin of Life” project of the Univ. Grenoble-Alpes. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Model comparison for GJ2060

Let us now use the dynamical masses obtained in Section 4.2 for the GJ2060 system to derive the age of the stars, and thus increase the constraints on the age of its young moving group ABDor. We retrieve the bolometric luminosities L of each star from Rodet et al. (2018). These were derived using a distance d = 15.69 ± 0.45 pc, which is compatible with the parallax that we obtain in this work (Table 3). As the binary is young, we use premain sequence (PMS) evolutionary models from the literature to relate mass, luminosity, and age.

Several evolutionary models for PMS stars rely on slightly different physics (e.g. atmospheric models, convection efficiency). We used models from Baraffe et al. (2015, hereafter BHAC15), D’Antona & Mazzitelli (1997, hereafter DM97), the PARSEC model (Bressan et al. 2012), the PISA model (Tognelli et al. 2011, 2012), the Darmouth model (Dotter et al. 2008; Feiden et al. 2015) and the one from Siess et al. (2000, hereafter Siess00). When the model requires stellar parameters (hydrogen, helium, or metal composition), we used the ones closest to the solar abundances (as given in Asplund et al. 2009). Such hypotheses are consistent with the solar-like metallicity derived for members of the ABDor moving group (McCarthy & Wilhelm 2014).

We plot the masses as a function of the system age for the given luminosity in Figure A.1. The shading shows the uncertainties associated with the error on the luminosity. Pre-main sequence low-mass stars are more luminous than their main sequence counterparts, meaning that a given luminosity can correspond to either a young low-mass star or an older more massive star. The plot diverges at the main sequence mass corresponding to the observed luminosity. Indeed, the luminosity evolves on much larger timescales when the star reaches the zero-age main sequence (at around 100 Myr old), meaning that all ages greater than 100 Myr are roughly compatible with the main sequence mass.

The discrepancy with the models is reduced compared to the study of Rodet et al. (2018), because of the slightly lower masses that we derived in this work. our values are now compatible with most of the models assuming the system is at least 100 Myr old. This age agrees with recent independent estimates of the ABDor moving group, arguing in favour of its similarity to the ~120-Myr Pleiades. However, our results are not compatible with the predictions from the Siess00 model, and are only marginally compatible with the predictions from DM97.

Fig. A.1 Primary (top) and secondary (bottom) mass of the GJ2060 system compared with mass–age relations coming from six evolutionary models. The red vertical lines correspond to the mass estimates derived in this work. The shades correspond to the uncertainty in the luminosities. The masses and luminosities suggest that the system is older than 100 Myr.

Appendix B Corner plots

Fig. B.1 Corner plot of all the MCMC iterations for Gl 494, showing the correlations between parameters and their density of probability diagonally. The delta values given for the five astrometric parameters are given with respect to the solution published in the IAD.

Fig. B.2 Same as Fig. B.1 but for GJ 2060.

Fig. B.3 Same as Fig. B.1 but for HIP 88745 adjustment (A) using only the absolute astrometry from HIPPARCOS and Gaia.

Fig. B.4 Same as Fig. B.1 but for HIP 88745 adjustment (AR) using both absolute and relative astrometric data.

References

All Tables

All Figures

	Fig. 1 Gl 494 best-fitting orbit from TMB (values in Table 2) as the green dotted line and sample orbits from the MCMC algorithm in black. The direct imaging observations (Bowler et al. 2020) with their associated error bars are in red.
In the text

Fig. 2 Radial velocity behaviour of Gl 494 predicted by the adjustment of the direct imaging, HIPPARCOS, and Gaia data. The best-fit orbit from TMB (values in Table 2) is the green dotted line and sample orbits from the MCMC are in black. The red dots represent the radial velocity observations with their associated error bars. The systemic velocity of 0.59 km s−1 was derived from a TMB adjustment of the radial velocities and the direct imaging data for the sole purpose of this visualisation.

In the text

	Fig. 3 Gl 494 masses compared to the mass-luminosity relations of Delfosse et al. (2000) and Mann et al. (2019), and the BT-Settl isochrones (Baraffe et al. 2015).
In the text

	Fig. 4 Orbits from the MCMC algorithm in black and of the best solution from TMB in green dotted line for GJ2060 (values in Table 2). The direct imaging observations with their associated error bars are in red and the new GRAVITY point with its error bars is in light blue.
In the text

	Fig. 5 Radial velocity behaviour of GJ2060 predicted by the adjustment of direct imaging and HIPPARCOS data. The best solution from TMB is the green dotted line (values in Table 2) and the MCMC results are in black. The radial velocity observations are in red dots with the associated error bars and with RV0 = 28.8 km s−1 (Rodet et al. 2018).
In the text

	Fig. 6 Orbits of HIP 88745 adjusted to the HIPPARCOS and Gaia data. The best-fit TMB solution is displayed as the green dotted line (values in Table 2, solution A), and a sampling of the MCMC solutions is shown in black. The direct imaging observations (in red) were not used in the adjustment and provide an independent validation of the orbit.
In the text

Fig. 7 Radial velocity behaviour of HIP 88745 predicted by the adjustment of the HIPPARCOS and Gaia data. The green dotted line represents the TMB best-fit solution (values in Table 2, solution A) and the black lines show the MCMC solutions. The radial velocity observations from SB9 are in red, the Gaia DR3 radial velocity is in orange, and the ELODIE archive is shown in blue. The systemic velocity of RV0 = 0.40 km s−1 was determined from an adjustment of the SB9 radial velocities together with TD and Gaia observations, which was performed for the sole purpose of this visualisation.

In the text

	Fig. 8 PARSEC isochrones with [M/H] = –0.6 dex colour-coded by age (in Gyr). The dark circles correspond to the masses derived in this work (AR solution) while the grey squares are from Jao et al. (2016).
In the text

	Fig. A.1 Primary (top) and secondary (bottom) mass of the GJ2060 system compared with mass–age relations coming from six evolutionary models. The red vertical lines correspond to the mass estimates derived in this work. The shades correspond to the uncertainty in the luminosities. The masses and luminosities suggest that the system is older than 100 Myr.
In the text

	Fig. B.1 Corner plot of all the MCMC iterations for Gl 494, showing the correlations between parameters and their density of probability diagonally. The delta values given for the five astrometric parameters are given with respect to the solution published in the IAD.
In the text

	Fig. B.2 Same as Fig. B.1 but for GJ 2060.
In the text

	Fig. B.3 Same as Fig. B.1 but for HIP 88745 adjustment (A) using only the absolute astrometry from HIPPARCOS and Gaia.
In the text

	Fig. B.4 Same as Fig. B.1 but for HIP 88745 adjustment (AR) using both absolute and relative astrometric data.
In the text