© The Authors 2024

1. Introduction

Gaia is a space mission of the European Space Agency (ESA) that was launched in 2013. The mission opens a window to explore the unprecedented high-precision astrometric data for a large population of asteroids. Its latest data release, the Focused Product Release (Gaia FPR) published in October 2023, contains 66 months of data for about 160,000 asteroids. By covering a main-belt asteroid’s typical orbital period, it has been shown that the Gaia data alone can provide very precise heliocentric orbits (Gaia Collaboration 2023). As shown in previous data releases, the Gaia data is precise along its scanning direction (the AL axis), in contrast to the perpendicular across-scan direction. In addition, the standard deviation at the CCD level is of the order of 0.2 mas (milliarcsecond) for an object of magnitude G ≈ 13 (Gaia Collaboration 2023). As a result, all data analyses and comparisons are done on this AL axis in the rest of this Letter.

Of more than 500 known binary and multiple systems in the Solar System, to date Gaia has observed approximately 300 of them in different regions. Thanks to the unprecedented precision of Gaia, it is possible to reveal the astrometric signature of binary asteroids. This is the case for the recently discovered binary (4337) Arecibo system (hereafter Arecibo), analysed with the Gaia DR3 data (Tanga et al. 2023), where the astrometric wobble was clearly detected in a time window of several days covering successive transits (see their Fig. 15 and Sect. 5.2.1).

In the present work we report on our continuing research on the Arecibo binary system by taking into account all FPR observations, corresponding to a total of 691 observations, or 90 transits. We performed a fit on both the heliocentric orbit and the relative orbit of the binary over approximately five years of data. By adding the positions issued from two stellar occultations (Gault et al. 2022), we were able to constrain the orbital and physical parameters of the binary. We thus derived the complete relative orbit, the flux ratio, and the mass of each component, and also give an estimate of their bulk densities and the corresponding ratio.

The second section presents the method and data model implemented. The third and forth sections are dedicated to the results and discussions.

2. Method

2.1. Heliocentric orbit

The heliocentric orbit is fitted to the FPR data by making use of the Numerical Integration of the Motion of Asteroids (NIMA) programme, presented in Desmars (2015). The orbit propagation and the computation of the partial derivatives includes the planetary ephemerides with four additional massive asteroids as perturbers, as well as relativistic corrections (on the observations and for the equation of motion). NIMA applies the weighted least-squares method to compute a differential correction to the initial state vector fitting the observations. The weighting scheme follows the Gaia error model for systematic and random error provided in the Gaia catalogue (Muinonen et al. 2022, Chapter 8). At this stage the post-fit residuals reflect, among other things, systematic effects due to the offset between the system’s computed centre of mass (CoM) and the observed centre of light (photocentre).

2.2. Binary photocentre–barycentre offset

In the case of a binary system, the total photocentre displacement is the result of the offset due to the presence of the satellite, added to the offsets from the components’ apparent discs and solar phase angle. Following Pravec & Scheirich (2012), this offset δη along AL can be written as

$$\begin{array}{c}\hfill \delta \eta =(\frac{I_2}{I}-\frac{M_2}{M})\mathrm{\Delta }\eta +\delta r_1+\delta r_2,\end{array}$$(1)

where Δη is the angular separation between the primary and the secondary projected onto the AL direction; I = I₁ + I₂ and M = M₁ + M₂ are the total intensity and total mass, respectively; and O_p = δr₁ + δr₂ is the projected photocentre offset due to the apparent discs of the components. This last term O_p can be neglected in a first approximation.

From the point of view of the Gaia telescope and its separation power, we had to consider different observational cases where a binary system can be either resolved or unresolved. An analytical model has been proposed in Lindegren (2022) to describe the photocentre versus CoM shift (δη) as function of the flux ratio (f = I₂/I₁), and mass ratio (q = M₂/M₁) of the binary components:

$$\begin{array}{c}\hfill \delta \eta =\{\begin{array}{cc}(\frac{f}{1+f}-\frac{q}{1+q})\mathrm{\Delta }\eta \hfill & \mathrm{if}|\mathrm{\Delta }\eta /u|\le 0.1\hfill \\ uB(f,\mathrm{\Delta }\eta /u)-\frac{q}{1+q}\mathrm{\Delta }\eta \hfill & \mathrm{if}0.1<|\mathrm{\Delta }\eta /u|\le 3-f\hfill \\ -\frac{q}{1+q}\mathrm{\Delta }\eta \hfill & \mathrm{if}3-f<|\mathrm{\Delta }\eta /u|\hfill \end{array}\end{array}$$(2)

Here u = 90 mas is the ‘resolution unit’ of Gaia (3/2 pixel size in AL).

In the first case, the binary system is not resolved. The recorded position corresponds to the binary’s centre of light, and the offset refers to the difference between the CoM and the photocentre of the whole system. Under the approximation of two spherical bodies, of the same albedo and the same density, this unresolved astrometric wobble amplitude can be expressed as function of the mass ratio alone: w = (1 + q^−2/3)⁻¹ − (1 + q⁻¹)⁻¹.

The second case is considered partially resolved. In this case the difference between the CoM and the photocentre depends on the bias function B(f, Δη/u), which is the root of the equation x + (x − p)fexp(p.x − p²/2) = 0 (Lindegren 2022). The bias is linked to the flux difference f and scaled projected separation p = Δη/u.

In the third case the system is resolved, and the offset corresponds to the position of the primary body with respect to the system’s centre of mass (assuming that only the primary is observed).

2.3. Relative orbit

We also developed a program to compute the Cartesian coordinates of the satellite with respect to the primary. This was done, under the assumption of two-body Keplerian motion, by taking into account the complete relative orbit (with seven orbital parameters: semi-major axis a, longitude of ascending node Ω, inclination i, argument of pericentre ω, eccentricity e, pericentre passed time t0, and period T), the flux ratio f, and the mass ratio q. This provides the separation Δη projected along the AL direction to be compared to the residuals from Sect. 2.1.

For each transit we computed a normal point by taking the mean AL residuals $\overline{x}$ for the N CCDs, and an associated standard deviation ${\sigma }_{\overline{x}}=\sqrt{{\displaystyle \sum _{i=1}^N}{\sigma }_i^2}/(N-1)$. We put a three-sigma threshold on these normal points to remove the lower quality data and the outliers. This resulted in a final set of 89 data points to perform the relative orbit fitting at transit level. This orbit fitting was done by making use of a Markov chain Monte Carlo (MCMC) method, minimizing the weighted mean square error $(\overline{x}-\delta \eta )/{\sigma }_{\overline{x}}$ on the offset in Eq. (2). The initial value was obtained from a simplex fit of the data; the prior followed a uniform distribution. A few hundred burn-in steps and a few thousand iterations were performed for each experiment.

As expected, we found that a strong correlation exists between the component masses and their flux ratio; we thus introduced a correcting factor s such that q = s.f^3/2, and eventually fit for f and s as free parameters. We also note that if both binary components are spherical bodies, with the same density and the same albedo, then s = 1. In this case, the flux and mass ratio are uniquely related to only the size ratio.

It appears that the Arecibo system was not resolved from the Gaia observations (the projected separation in the AL direction is always significantly smaller than the resolution unit). As a result, the semi-major axis is not constrained. To further fit this parameter, we included into the fit the two relative position data points measured from ground-based stellar occultation (Gault et al. 2022). The target function is now

$$\begin{array}{c}\hfill lnL=-\frac{1}{2}[{\displaystyle \sum {\left[\frac{\overline{x}-\delta \eta }{{\sigma }_{\overline{x}}}\right]}^2}+{\displaystyle \sum {\left[\frac{l-l_c}{{\sigma }_l}\right]}^2}+{\displaystyle \sum {\left[\frac{pa-p{a}_c}{{\sigma }_{\mathit{pa}}}\right]}^2}],\end{array}$$(3)

where l is the separation and pa the position angle, with their associated uncertainties.

3. Results

We obtain a Gaussian-like posterior probability distribution function in the orbital parameters; the fitted orbital parameters are given in Table 1. The left panel of Fig. 1 shows the computed photocentre shifts from the binary model as a function of the AL residuals from NIMA, which has a Pearson correlation coefficient of 0.67. On the right, the residuals are displayed after correcting for the binary astrometric wobble. The goodness of fit is thus globally improved: the RMS of the post-fit is decreased by approximately 25%. The total mass of the system, obtained from the third Kepler law, is well constrained by a and T, resulting in M = 4.3 ± 0.5 × 10¹⁵ kg. We note that the orbital period is very well defined, thanks to the long time-frame of the FPR data. The precision on the total mass hence remains mostly limited by the knowledge of the semi-major axis.

Fig. 1. Comparison between the average heliocentric transit residuals and the computed binary photocentre–barycentre offsets. The right panel shows the equivalent plot after correcting for the computed shift. The horizontal lines represent the displacement between the data points before and after correction. A typical error bar, not shown in the figure, is of the order of 0.6 mas.

In addition, two solutions of the spin pole and period of the primary are derived from the Gaia photometric data (Ďurech & Hanuš 2023). We include them in Table 1 with their uncertainties. The similarity between the orbital pole and one of the spin poles and the corresponding spin period, together with the relatively small eccentricity, suggests a tidally locked system on a quasi-circular equatorial prograde orbit.

The corner plot of all estimated parameters is given in Appendix A. This figure shows that the inclusion of occultation data, and the relative positions, is essential for deriving the semi-major axis and also narrowing down the other orbital parameters. Moreover, one can determine that there is no correlation between s, f, and the other parameters. We show in Fig. 2 the posterior distribution of s and f, with the global (unresolved) astrometric wobble amplitude w = f/(1 + f)−q/(1 + q). One can see that w is well constrained, while the flux ratio and the correcting factor are correlated (through mass ratio).

Fig. 2. Posterior distribution of the parameters s, f, and w. The astrometric wobble amplitude w is computed from the equation introduced in the text. The red line is plotted solving the relationship of f and s with w = 0.053 (median of its distribution). The blue line is taken at the mode of the distribution of f = 0.11, and we take it as one of the most probable solutions of this system. The corresponding value of s is 1.33.

4. Discussion

Several assumptions are made in this preliminary work: The photocentre offsets O_p in Eq. (1), due to the phase effect of each component, were neglected; we checked that they remain smaller than the astrometric wobble δη. Similarly, possible offsets in the relative position from the occultation data modelling were neglected. The orbit computations were performed sequentially, starting with the heliocentric orbit for the centre of mass. We did not consider any precession of the binary orbit over the five years of Gaia observations due to external perturbations. We also did not consider a possible dynamical flattening of the primary; assuming a theoretical c/a ratio as high as 1/2, the mass and density determination could be underestimated by approximately 15%. All these effects remain small, but will be included in a future work.

The relative separation of the Arecibo binary on the AL axis is small, and the system is never resolved by the Gaia telescope. All Arecibo data fell into the first two cases described in Eq. (2), which are almost identical given the small separation Δη (Lindegren 2022). As seen in Fig. 2, the astrometric wobble amplitude is converging to w = 0.053 ± 0.005, while f, s, and hence q are correlated. The correlation curve for w = 0.053 is plotted as a red line in the figure with flux ratio range from 0 to 1. If we take the mode of f as the most probable solution we have f = 0.11, s = 1.33, and hence q = 0.05. It should be noted that the correcting factor s differs significantly from unity, which suggests a non-spherical shape, a different density, or a different albedo, or all three combined.

The thermal infrared survey WISE (Mainzer et al. 2019) gives a surface equivalent diameter D_eff = 19.7 ± 0.2 km for Arecibo. If we consider a pseudo system composed of two spheres of diameter D₁ and D₂, and same albedo, we obtain from $f=(D_2^2/D_1^2)$ the absolute sizes D₁ = 18.7 km and D₂ = 6.2 km. Combining further with the mass ratio q, a bulk density of each component is derived, resulting in ρ₁ = 1200 kg/m³ and ρ₂ = 1600 kg/m³. These density values are in line with the fact that Arecibo, together with (90) Antiope and (379) Huenna, is a member of the Themis family, an icy family in the outer belt (Castillo-Rogez & Schmidt 2010; Hargrove et al. 2015; Marsset et al. 2016). Furthermore, irrespective of the components’ absolute sizes, but still assuming similar albedos, the density ratio is ρ₂/ρ₁ = s ≈ 1.33. Thus, the secondary is more dense, or less porous than the primary. In general, for two spheres, (ρ₂/ρ₁).(A₁/A₂)^3/2 = s, so that a more reflective primary would lower this density ratio. In any case, the shapes and sizes of the components appear to be poorly known, with inconsistent values (by a factor of 1.5) from either occultation data or other radiometric data from Akari. For this reason, no error bars are given here for the density estimates.

5. Conclusion

By analysing Gaia data from the FPR release, completed by two ground-based occultations, we derived the orbital and physical parameters of the binary (4337) Arecibo from an analysis of its relative orbit.

In previous studies, for example Margot et al. (2002), Ostro et al. (2006), Naidu et al. (2012), and the summary in Margot et al. (2015), radar observations have proved to be a powerful method for detecting and inferring the individual physical properties of near-Earth binary systems. In our work, we presented another method for obtaining such results for a small main-belt binary, by analysing high-precision optical astrometric data. The procedure can be applied to other known systems observed as unresolved and resolved binaries by Gaia and ground- or space-based telescopes. Moreover, it can be implemented on the binary candidates from Liberato et al. (2024), and can further provide the ephemeris of the mutual orbit and predicted position of the secondary. We will also include perturbed two-body problem for the dynamical modelling in future work.

The accurate estimation of the bulk density of each component is complicated by the limited knowledge of their sizes and shapes (and hence volume). More observations are needed to constrain their absolute or relative sizes. This can be done with additional occultations, and more chords distributed along the two bodies. Several occultations in the coming years can be of interest to this purpose. For instance, there are two stellar occultation predictions for the system in December 2024 and January 2025 (see Fig. B.1). It is also possible to predict the seasons of mutual phenomena (occultations, eclipses) in such a binary system (Emelyanov et al. 2023). The next season in this system will occur in November 2025 (see Appendix B), while at solar elongation of approximately 60°, with a predicted magnitude drop of ΔG ≈ 0.1 in the photometric lightcurve. Such mutual phenomena or stellar occultation observations will be very valuable for better constraining the components sizes, and hence the respective density. Additional data at longer elongations, or with resolved observations, would help in constraining the orbit, and better separating the mass and flux ratios.

Acknowledgments

The authors would like to thank F. Arenou (Paris observatory) for fruitful discussions. 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. We made use of the emcee (Foreman-Mackey et al. 2013) and corner (Foreman-Mackey 2016) python package for data fitting and visualization. In addition, we made use of ephemerides from IMCCE https://ssp.imcce.fr.

References

Appendix A: Full results

Fig. A.1. Posterior distribution of all nine fitted parameters. All angles are in radians, in ecliptic J2000 reference frame; the semi-major axis is in kilometres; and the time-related parameters are in days. In grey: MCMC sample fit from Gaia data only; in blue: MCMC samples, including the data from the two occultations. The text and dashed line are given at the 16th, 50th, and 84th percentiles for the distribution including occultation.

Fig. A.2. Residuals of the satellite position between the two observed occultations (blue dot) and the prediction (black cross). Left: 19 May 2021; Right: 9 June 2021. The primary’s position is at the centre (0,0). The shaded area is plotted following the error bar in Gault et al. (2022).

Appendix B: Event predictions

B.1. Stellar occultations

Fig. B.1. Two stellar occultation predictions of Arecibo system in the near future. Left: Event on 21 December 2024, with the star magnitude of 13.3 and a drop in magnitude of 2.5. Right: Event on 18 January 2025, with the star magnitude of 13.1 and a drop in magnitude of 2.5.

B.2. Mutual events

Fig. B.2. Prediction of the Earth’s elevation (sub-Earth point latitude βSEP) above the orbital plane. The zoomed-in plot shows the plane-of-sky view at the time of one predicted mutual event; the blue line is the projected relative orbit.

All Tables

All Figures

	Fig. 1. Comparison between the average heliocentric transit residuals and the computed binary photocentre–barycentre offsets. The right panel shows the equivalent plot after correcting for the computed shift. The horizontal lines represent the displacement between the data points before and after correction. A typical error bar, not shown in the figure, is of the order of 0.6 mas.
In the text

Fig. 2. Posterior distribution of the parameters s, f, and w. The astrometric wobble amplitude w is computed from the equation introduced in the text. The red line is plotted solving the relationship of f and s with w = 0.053 (median of its distribution). The blue line is taken at the mode of the distribution of f = 0.11, and we take it as one of the most probable solutions of this system. The corresponding value of s is 1.33.

In the text

Fig. A.1. Posterior distribution of all nine fitted parameters. All angles are in radians, in ecliptic J2000 reference frame; the semi-major axis is in kilometres; and the time-related parameters are in days. In grey: MCMC sample fit from Gaia data only; in blue: MCMC samples, including the data from the two occultations. The text and dashed line are given at the 16th, 50th, and 84th percentiles for the distribution including occultation.

In the text

	Fig. A.2. Residuals of the satellite position between the two observed occultations (blue dot) and the prediction (black cross). Left: 19 May 2021; Right: 9 June 2021. The primary’s position is at the centre (0,0). The shaded area is plotted following the error bar in Gault et al. (2022).
In the text

	Fig. B.1. Two stellar occultation predictions of Arecibo system in the near future. Left: Event on 21 December 2024, with the star magnitude of 13.3 and a drop in magnitude of 2.5. Right: Event on 18 January 2025, with the star magnitude of 13.1 and a drop in magnitude of 2.5.
In the text

	Fig. B.2. Prediction of the Earth’s elevation (sub-Earth point latitude βSEP) above the orbital plane. The zoomed-in plot shows the plane-of-sky view at the time of one predicted mutual event; the blue line is the projected relative orbit.
In the text