© The Authors 2024

1 Introduction

Binary asteroids have attracted the attention of the scientific community due to their interesting properties and the significant impact they have on our understanding of the Solar System. Unlike single asteroids, binary systems offer unique insights into many fundamental processes, including the formation and evolution of planetary bodies, collision dynamics, and gravitational interactions.

Numerical simulations suggest the existence of nearly equal-sized binaries as a byproduct of catastrophic collisions (Durda et al. 2004) but they remain essentially unknown. On the other hand, multiple craters on Mars demonstrate that asteroid binaries with properties not represented by the known sample should exist, or have existed (Vavilov et al. 2022). More generally, the formation of asteroid companions by fragment ejection and in-orbit re-accumulation (Walsh & Richardson 2006; Walsh et al. 2008; Ćuk 2007; Pravec et al. 2010; Jacobson et al. 2013; Madeira et al. 2023) seems to be common among near-Earth asteroids (and probably small main belt objects) but the collisional evolution is also expected to play a role (Doressoundiram et al. 1997; Michel et al. 2001; Durda et al. 2007; Jutzi 2019), with an undefined boundary between the two creation mechanisms. However, some open questions on the origin of binaries and the evolution of asteroids cannot be answered without a greater number of known objects and a better knowledge of the sample.

The presence of asteroid companions has been revealed by different techniques, such as high-resolution imaging from ground-based and space-based telescopes, photometry, radar ranging, and stellar occultations. When taken together, these different approaches are complementary and cover a wide range of separations and size ratios for the system components. However, the known population of binary asteroids is still strongly biased.

While imaging favours large separations of the primary and the secondary, and bright primaries in the case of adaptive optics, photometric studies, based on detecting mutual events, strongly favour compact systems (Merline et al. 2002; Pravec et al. 2006; Richardson & Walsh 2006). Radar techniques are strongly limited in range and are efficient for near-Earth asteroids (Ostro et al. 2002).

Some satellites have also been discovered by space probes and studied during close encounters. While they constitute a precious sample of detailed information, in statistical terms their contribution to the knowledge of the global population of binaries is limited. Astrometric detection of binaries alone cannot solve all these problems, but can very usefully fill an unexplored space of exploration (Pravec & Scheirich 2012). The astrometric signature of the presence of asteroid satellites is expected to show up in the residuals of the asteroid motion (on the sky or, equivalently, concerning the heliocentric orbit) due to the difference between the position of the photocentre (the one usually provided by the observations) and the centre of mass of the system. This difference, which in principle can also be due to the irregular shape of a single object, is enhanced by the presence of a satellite.

Anomalous residuals in the astrometry, linked to the presence of a satellite, have been revealed in the case of (90482) Orcus (Ortiz et al. 2011). For trans-Neptunian objects, the effect is amplified by the large separation of the components, particularly in cases of high size ratio (secondary/primary).

The astrometric signature is generally much harder to detect for main belt asteroids since it requires much more precise astrometry to reveal a perturbed motion due to a companion that is not spatially resolved. A rather solid detection, the first using Gaia astrometry, was obtained for (4337) Arecibo (Tanga et al. 2023). Also in this case, however, the presence of a satellite was already known. A systematic search of astrometry for new binary systems is a much more complex task. Only the exploration of a large data set of homogeneous measurements of extremely high quality can potentially be successful. The data produced by the Gaia mission of ESA are the main sources available today for this search. A first attempt on the very limited data set (in the time range and the number of objects) provided by Gaia DR2 was unsuccessful (Segev et al. 2023).

In the following we present a blind search of astrometric binaries in the Gaia DR3 asteroid astrometry. Our exploration represents a first step only where we focus on the detection of periodic signals (wobble, in the following) in the residuals obtained from orbit adjustment to the data. For the moment this approach limits our investigation to sets of astrometric positions collected in short time spans. As a consequence, short-period systems are favoured. The result of our search is also filtered on the basis of rather strict criteria. These choices produce results that are certainly not comprehensive of what can be obtained from Gaia data. Nonetheless, it provides a first consolidated set of candidate detections that can be fruitfully explored by other techniques for confirmation.

The article is organised as follows. In Sect. 2 we recall the main properties of the data set that we explore and the goal of the search. Section 3 presents the method of period search. The selection, filtering, and validation of the candidates found, by statistical and by physical criteria, are presented in Sect. 4. Section 5 presents the content of the final list of binary candidates. We put our findings in a broader context in Sect. 6, and resume the work, with future perspectives, in Sect. 7.

2 Astrometry and binary signatures: Data set and model

The asteroid astrometry in Gaia DR3 for a large number of asteroids provides an outstanding opportunity for a search of astrometric signatures produced by possible satellites. The properties of Gaia astrometry are illustrated in detail by Tanga et al. (2023).

Here we recall a few properties specifically relevant to our search.

The first property is the time distribution of observations. The two Gaia telescopes scan the sky continuously with a six-hour rotation period. The two fields of view (FOV) sweep the same area of the sky 106 minutes apart. The spin axis precesses in such a way that a large field overlap is granted every subsequent rotation. The combination of this rotation-precession with the revolution around the Sun results in a sparse series of detections (also called transits in the FOV) for any given source, fixed or moving. However, when the motion on the scanning path in the sky combines favourably with the position and displacement of a target source, a short consecutive series of detections over a few satellite rotations is possible.

Over the coverage of Gaia DR3 (34 months), it is rather common to obtain sequences of a few consecutive observations spanning a few rotations (12h, four detections) for each object. A few occurrences of longer sequences with several rotations up to 10–12 detections (five or six revolutions) can also happen, although less frequently.

As is discussed in the following section, the length of time windows of consecutive data and the presence of constant periods in the data sampling control, respectively, the longest detectable period in our analysis, and the possible emergence of spurious frequencies. At the same time, irregular sampling allows the exploration of a large range of frequencies (Eyer & Bartholdi 1999; VanderPlas 2018).

In addition, the nearly mono-dimensional astrometry provided for each transit brings accurate information only in the direction of the scan (along-scan, AL). The only measurable residual from orbital fitting is the projection of the total residual onto the AL direction, whose direction is known from the attitude data of the satellite. The direction of the total post-orbital fit residual in the sky, however, is not directly measurable from Gaia astrometry alone at a given epoch (or in an interval of time in which the geometry of the system concerning the scan direction does not evolve). Constraints on the total wobble in the astrometry, and the possible primary-secondary separation, are the result of applying some physical constraints to the measured wobble projection and period.

For this first search for astrometric binaries, we use the simplest first-order model for the wobbling motion, assuming systems composed of uniformly illuminated spheres, of equal bulk densities and albedos. The other important assumption is that the components are not resolved so that the wobble is generated by the difference in position between the barycentre of the system (to which the orbit is adjusted) and the photocentre. Given the size of the Gaia pixel (60 mas) in the AL direction, this should be the case for a wide range of binaries in a reasonable range of sizes.

Figure 1 illustrates the principles discussed, where at the top of the plot it is shown a representation of an ideal case where Gaia is observing a binary asteroid perpendicularly to the plane of rotation of the objects around the barycentre, the along-scan direction (AL) of the Gaia satellite is on the x-axis and the across-scan direction is on the y-axis. In principle, at any given epoch, the projection of the instantaneous wobble amplitude α for the satellite position on its orbit, projected along AL, is provided by the residual of Gaia astrometry to the orbital fit. The blue and green dashed lines represent the projection of the positions of the photocentre and the barycentre of the binary, respectively. The difference between these projections represents the AL residuals, as shown at the bottom of Fig. 1.

With the assumptions of uniform illumination and unresolved components, it is possible to derive the wobble amplitude as the difference between the position of the barycentre and the position of the photocentre. This is a function of the mass ratio q only (Hestroffer et al. 2010), scaled for the separation a:

$$\alpha =a\left|\frac{1}{1+q^{-1}}-\frac{1}{1+q^{-\frac{2}{3}}}\right|,q\in ]0,1]$$(1)

This function (Fig. 2) has a maximum at q_max ≈ 0.154 (corresponding to a size ratio $k=\sqrt[3]{q_{\max }}=0.536$) and tends to zero at the extremes of the range of q ∊]0, 1], that is, for a vanishingly small satellite or equal-sized components. During a single transit in the Astrometric Field of Gaia (about 40 seconds), up to nine consecutive positions are obtained. We fit all Gaia DR3 astrometry for each asteroid by a version of the ORBFIT soft-ware¹, optimised to exploit the full accuracy of Gaia, following the same procedure as in Gaia Collaboration (2018). From the orbital fit, for each astrometric point, a residual is obtained in the AL direction.

Given that the observation geometry over that time span is nearly constant (and a relative displacement of a satellite negligible) we average the residuals from all transits from one observation into a single transit. Their standard deviation also provides uncertainties on the measured residual itself. The average AL residuals per transit, and their uncertainty, are the data on which we base our analysis. To prevent future confusion, throughout the work we refer to the residuals from the orbital fitting of the average transit observations in the along-scan direction simply as AL residuals.

Fig. 1 Representation of a simplified case where the Gaia satellite would be observing a binary system perpendicularly to the plane of rotation of the bodies, along with the projections of the positions of the barycentre and the photocentre on the AL direction and the representation of the orbital fit residuals’ projection in the along-scan direction. AL direction is represented horizontally for simplicity, but in reality it can take any orientation, slowly rotating by some degrees on a scale of several hours.

3 Wobbling detection

For this first data exploration, our search targets intervals of time in which several consecutive astrometric measurements are available, sufficiently short to neglect any change of orientation of the asteroid with respect to Gaia. On the observations available over that interval, a period search is performed.

Fig. 2 Plot of the amplitude of the photocentre wobbling as a function of the binary mass ratio, as described in Eq. (1). In this example we use a separation a = 10 mas, which means that we would observe a maximum photocentre wobbling of 0.9 mas for a mass ratio q ≈ 0.1538.

3.1 Sample selection

To easily identify the longest sequences of consecutive astrometric measurements, we perform a convolution of the observation epochs with a rectangular function and find the peaks of large data density, as shown in Fig. 3. In the left plot, we show the density of AL residuals over time for asteroid (3457) Arnenordheim, and the peak close to 450 days represents the 2.75-day window of time with 17 observations shown in the right plot. The rectangular function has a unity value over an interval of ten days and zero elsewhere.

The considered threshold of ten days is a compromise between the length of the time series and the resulting number of candidates, as longer time series allow for better frequency analysis but reduce the number of targets. It should be noted that this selection method allows for gaps in the data sequence, where the observations are absent due to some technical reasons (missed observation, rejection for low quality). Over the whole DR3 data set of 156 801 asteroids, our selection extracts 30 030 objects with sequences satisfying the above-mentioned criteria.

Next, we identify the observations in each one of the windows and perform an outlier check, removing the points with nominal values or standard deviations larger than 3σ of the dispersion of AL residuals in the window of search.

Fig. 3 Plots of the density of observations (number in ten consecutive days) for asteroid (3457) Arnenordheim and selected sample of data. Left: resulting convolution of the observation time span in days vs the number of observations in each ten-day window. Right: sample of data at the peak of the 10-day window of search combining 17 observations. The x-axis is the observation date and the y-axis is the average AL residual from the orbital fit. Each dot represents one averaged set of transits observed in the FOV 1 in blue and FOV 2 in green.

3.2 Period search

The Generalised Lomb–Scargle Periodogram (GLSP) (Zechmeister & Kürster 2009; VanderPlas & Ivezic 2015; VanderPlas 2018) is a statistical technique for analysing irregularly sampled time series data and identifying periodic signals in the data. The GLSP method is a generalisation of the classical Lomb–Scargle Periodogram (LSP) (Lomb 1976; Scargle 1982).

As described in VanderPlas (2018), the GLSP corresponds to a weighted least-squares problem where a sinusoidal signal of unknown amplitude and phase plus an unknown constant are jointly estimated at each test frequency. Formally, using the notation from VanderPlas (2018), the GLSP P(f) can be defined from the following equations:

$$y_{\text{model }}(t;f)\equiv y_0(f)+A_{\text{f}}\sin \left(2\pi ft+{\varphi }_f\right),$$(2)

$${\chi }^2(f)\equiv {\displaystyle \sum _{n=1}^N{\left(\frac{y_n-y_{\text{model }}\left(t_n;f\right)}{{\sigma }_n}\right)}^2},$$(3)

$$P(f)\equiv \frac{{\widehat{\chi }}_0^2-{\widehat{\chi }}^2(f)}{{\widehat{\chi }}_0^2}.$$(4)

Here, y_model Eq. (2) is the assumed true model: a sinusoid (unknown frequency f, amplitude A_f, and phase ϕ_f) plus a constant offset term y₀ that may depend on the considered frequency. If the N data points of a time series y, y(t_n) ≡ y_n for n = 1, ⋯ , N, correspond to the sum of the true model of Eq. (2) plus uncorrelated Gaussian noise of variance ${\sigma }_n^2$ for each of the n samples, then the weighted sum in Eq. (3) is a random variable with N degrees of freedom.

The GLSP computes the residual sum of squares (RSS) of the weighted residuals (noted ${\widehat{\chi }}^2(f)$) obtained when adjusting a sinusoid with unknown amplitude and phase for each test frequency plus a constant and compares it to the RSS obtained when fitting only a constant $\left({\widehat{\chi }}_0^2\right)$. Precisely, ${\widehat{\chi }}^2(f)\left(\text{ resp., }{\widehat{\chi }}_0^2\right)$ is obtained by replacing y_model by the estimated sinusoid plus constant (resp., the estimated constant) in Eq. (3). The GLSP P(f) in Eq. (4) is thus the relative reduction in the sum of square brought by the adjusted model at considered frequency f(with some abuse of notation on f). The higher P(f), the smaller ${\widehat{\chi }}^2(f)$, which means that the model and the data are in greater agreement.

The GLSP is probably the most common technique to detect periodic signals in time series with irregular sampling. We opted for this approach because it has proven efficiency and is computationally simple enough to allow for large-scale Monte Carlo simulations (this is an important aspect of our study, see next section). The GLSP in Eq. (4) is implemented with the ASTROPY package (Astropy Collaboration 2013, 2018, 2022), with error bars σ_n estimated as explained in the last paragraph of Sect. 2.

Given the sampling frequency and the time range of ten days for our data sequences, we limit the search to periods between 3 h and 3 times the length of the observations sample and apply our period search to all the asteroids in the 30030 objects selected in the previous step.

The period search provides a periodogram for all the studied bodies, and they are very different depending on the time sampling of observations. The highest peak in the periodogram represents the period that best fits the sample. The left plot of Fig. 4 shows, as an example, the GLSP obtained for the asteroid (106016) 2000 SS293 from the sample of AL residuals shown in the right plot of the same figure. The red line is the GLSP from the fitting on the data of AL residuals in the current window of search. Using the frequency with the highest power from the GLSP on the data we obtain the sinusoidal fit with amplitude A and, along with the same time sampling and uncertainties, we run the GLSP to obtain the blue line that represents the profile of the periodogram only due to the sinusoidal signal. Finally, the spectral window (green line) is obtained by computing the Lomb–Scargle periodogram of a time series with constant unit amplitude².

It is interesting to compare the shapes of the GLSP of the data (red curve in the left panel of Fig. 4), of the main estimated sinusoid (blue curve) and spectral window (green curve). The spectral window indicates the specific frequency pattern related to the specific time sampling in a time series. Similarities in the structures of the spectral window and the GLSP are often a clear sign that some peaks are fake (sampling-induced) rather than genuine signal periodicities. Similarly, the pure sinusoid analysis (blue curve) encodes the same information but focuses on the analysis of a particular frequency. As an example, we see that the peak at almost 7 cycles per day in the data GLSP (red curve) is probably, from the blue curve, a spurious sidelobe created by the specific time sampling. This is also visible from the green curve, which shows that an alias peak appears at four rotations/days after the main peak.

From the considerations explained above, it is clear that a detection method based on the highest peak of the periodogram must evaluate the “statistical significance” of a given maximum peak’s height. One way to do this is to measure the probability of obtaining a value of the largest peak as high as the one observed in the data when there is only noise in the time series. This is investigated in the next section.

Fig. 4 Examples of periodograms obtained for asteroid (106016) 2000 SS293 (left plot), with residuals shown on the right plot. The y-axis shows the GLSP P(f) for each frequency shown on the x-axis. The red line represents the periodogram obtained for the data analysed. The blue line represents the GLSP obtained for a purely sinusoidal signal whose parameters correspond to the largest peak of P(f) (about 3 cycles per day here), with the same time sampling as the data and the same error bars for the weights. Finally, the green line represents the spectral window of the time sampling. The right plot shows the AL residual time series that reproduces the periodogram on the left plot.

4 Binary candidates selection method

4.1 Statistical selection

For each window of search of each object, we computed the GLSP as described in Sect. 3.2 and obtained the value (say, υ) of the highest peak. A common proxy for the “statistical significance” of υ is the p-value (see e.g. Efron 2010), defined as follows. Let V be the random variable defined as the maximum value of the periodogram, conditional on the time sampling grid and error bars of the considered data. The p-value of υ is the probability of obtaining a maximum peak value at least as large as υ under the signal-absent hypothesis, with the same time sampling and noise error bars (denoted by ℋ₀). Formally this probability is

$$p(\upsilon )\equiv \mathrm{Pr}\left(V>\upsilon \mid {\mathscr{H}}_0\right).$$(5)

In practice, the p-value p(υ) of a particular time series can be evaluated by Monte Carlo simulations by estimating empirically this probability as follows. We generate 10000 simulated time series sampled with the same time sampling as the time series under investigation. Each time series is composed of a Gaussian noise consistent with the error bars plus a random constant (an offset drawn randomly from a zero-mean Laplacian distribution, which is consistent with our data, see Appendix A). The GLSPs of the 10 000 time series are computed and we count how often their highest peak has an equal or larger power than the one from the GLSP under investigation (υ).

A p-value is computed in this way for each of the selected 30 030 time series. To be selected as a possible candidate, the time series must satisfy three conditions. The first condition is that the p-value estimated from the MC under ℋ₀ must be smaller than 0.05. The second condition is a check of the “robustness” of the periodicity which we call “quality factor” below, and which also provides a method to estimate uncertainties for the amplitude and frequency associated with the main peak of the GLSP. For this second condition, we perform MC simulations under the hypothesis that there is a sinusoidal signal in the data with the frequency and amplitude found from the GLSP of the candidate time series (plus noise with the considered error bars, plus an unknown offset, and with same time-sampling; we note this hypothesis ${\widehat{\mathscr{H}}}_1$). We then perform another series of 5000 MC simulations of synthetic time series under ${\widehat{\mathscr{H}}}_1$.

A robust period detection should repeat nearly identically several times with different realisations of the synthetic data. To obtain the quality factor of the signal, denoted by Q below, we estimate the fraction of time series for which the GLSP highest peak is found ±0.15 cycles per day around the injected frequency. Our second condition is that Q should be larger than 0.5, meaning that the initially injected period should be retrieved more than half of the time in nominal noise and sampling conditions.

We underline that a large maximum peak caused by noise in the GLSP (i.e. a false alarm) will lead to a high value of the quality factor, because the wobbling amplitude estimated from the data will be large, so any strong but ‘fake’ detection will appear as robust concerning the Q analysis. In the opposite direction, a maximum peak with low amplitude will lead to a low estimated wobbling amplitude and a low value of the quality factor. Hence, while a large value of Q is not a guarantee that the detection is robust, the cross-analysis brought by the quality factor essentially provides a supplementary flag that some detections may not be robust (this will be further investigated in Sect. 4.2).

Finally, we turn to the third condition. An important aspect is that we are looking for astrometric binaries, which means that Gaia satellite would not have been able to observe the asteroid and the satellite separately. A wobble amplitude of 20 mas translates to a minimum separation of about 220 mas, which is above the Gaia resolution limit of about 200 mas. Hence, the third condition is that the wobble amplitude estimated $\widehat{A}$ must be smaller than 20 mas.

In summary, if the signal detected in the current window has a p-value < 5% as estimated from MC simulations under ℋ₀, a quality factor Q > 0.5 as estimated from MC simulations under ${\widehat{\mathscr{H}}}_1$ and estimated wobble amplitude smaller than 20 mas, then the object is selected as a preliminary candidate and the values estimated in the current window are stored. The threshold of p-value < 5% is not very conservative since the expected number of false positives is about 1500 (5% of 30 030). In our case, the number of candidates that pass the triple test above is 3038. This number (larger than 1500) can be explained by two reasons. First, by a population of true binaries that lead to smaller p-values. Second, by possible modelling errors (for instance, some un-modelled, low-order polynomial trends would lead to GLSP with larger peaks at low frequencies)³. However, as we discuss in Sect. 4.2, the latter seems unlikely.

To obtain estimates of the period and amplitude uncertainties, we use the results from the MC simulations under ${\widehat{\mathscr{H}}}_1$. In the case of mean period determination, we look for frequencies that are no more than 0.25 cycles per day away from the best frequency found with the first GLSP and scan along the distribution to find the mean value, the 68% confidence interval, which corresponds to the mean ±1σ for a Gaussian distribution, and the 95% confidence interval. The value estimated from the first GLSP run is selected as the wobble period and the 95% CI is the uncertainty on the period. The procedure is the same for the amplitude, but in this case, we limit the interval for amplitudes up to 2 mas away from the first amplitude estimate. Figure 5 shows an example of uncertainty determination for the period on the top plot and the amplitude on the bottom plot.

If the object under the period search has more than one window of data, as discussed in Sect. 3.1, we repeat the same procedure described in all the other windows as well. Figure 6 shows a flowchart that summarises the statistical selection processes discussed above.

Fig. 5 Examples of uncertainty determination from MC simulations under ${\widehat{\mathscr{H}}}_1$. The top and bottom plots show, respectively, the distribution of periods and amplitudes obtained from the best fit in each MC simulation (in blue). The red lines show the value estimated from the first GLSP run only on Gaia data, and the solid black line represents the median obtained from the distribution of amplitude and period values from the MC under ${\widehat{\mathscr{H}}}_1$. The dotted vertical lines show the 68% distribution interval while the dashed lines represent the 95% confidence interval of the distribution. The final results are given as the first estimate from the Gaia data (red line), and the uncertainty is the 95% interval. For instance, in this example, the final amplitude value is 1.85 ± 0.08.

4.2 Analysis of the detection performance of the statistical selection process

To evaluate the detection performance of our method, in this section, we conduct a general study based on the analysis of two cases. In the first case, there is indeed a signal in each time series analysed and we look at the effectiveness of the method at detecting and characterising it. In the second case, the time series have no signal and we examine how incorrect the method can be in falsely selecting residuals with only noise. A supplementary investigation into the detection of eccentric binary systems is discussed in Appendix B.

4.2.1 Signal detection

In this analysis, we study the signal detection performance by grouping the simulated time series into classes of signal-to-noise ratio and looking at the results in terms of detections (at the right period or not), missed detections, associated estimated p-values and Q-factors. We use two signal-to-noise ratio definitions:

$$\text{“BaseS/N ”}:S/N_{\text{base }}\equiv A/\sigma \text{\hspace{1em} and }$$(6)

$$\text{“Post-fitS/N ”}:\widehat{S/N}\equiv \widehat{A}/\overline{\sigma },$$(7)

where A is the true amplitude of the time series and σ is the nominal average uncertainty of the data set, $\widehat{A}$ is the amplitude estimate from the data time series and $\overline{\sigma }$ is the estimated average uncertainty of the data set. To obtain the performance analysis we perform the period search and selection method shown in Fig. 6 for six different S/N_base values, with Monte Carlo simulations based on 1000 time series of synthetic AL residual.

First, we choose a set of different values of S/N_base between 0.5 and 3.0. We create each set of synthetic data from a sinusoid with a known amplitude A of 1 mas, a period of 23 hours (a frequency close to but not exactly 1 cycle per day to avoid possible bias) and a random phase between 0 and 2π. The nominal average noise σ for each set of synthetic samples is estimated using Eq. (6).

For the simulations of the time series, we draw 12 observation epochs, randomly distributed in 6 days. The small number of points, the randomness in the data and the time sampling are aimed not only at representing the data from Gaia, which is usually modestly sampled over a few days but also at providing conservative results that cover most of the non-ideal cases.

For each time series, the actual standard deviations of the errors σ_n are obtained by a random perturbation of σ (using a zero-mean, normal distribution of std 0.8; this setting is aimed at creating some scatter in the actual S/N within each S/N class), and we add one random Laplacian offset as described in Appendix A. Then, for each one of the 6000 synthetic time series simulated in this way, we run the statistical selection (Sect. 4.1).

Figure 7 shows the results for the (p, Q) pairs of 6000 time series simulated. We found that, as expected, for low values of S/N_base the parameters p and Q tend to be more scattered and present less correlation between themselves. But for higher S/N_base values, the couples (p, Q) tend to cluster around the “best conditions” represented by very low p-values and high Q.

Intuitively, the larger the S/N_base of a signal the easier it would be for the method to detect it. It is also expected that very noisy signals hamper the period fitting and provide incorrect detections more often, and our results follow that behaviour. For time series with S/N_base close to 0.5, the method achieves 11.3% of detections, of which 69% have incorrect periods. For the S/N_base class of 1.0, the detections represent 56% and about half present incorrect periods. For the S/N_base classes larger than 2.0, the method detects more than 93% of the signals with a low rate of incorrect detections.

It is also interesting to check the distribution of the quality factor (Q) and the p-value depending on the S/N_base, as presented in Fig. 8. It shows the distribution of best periods fitted on the top plots, the distribution of Q in the plots in the middle and the p-value distributions on the bottom plots for each different class of S/N_base studied. As discussed before, we notice that for small S/N_base there is a significant amount of detections with periods far from the real value, represented by the dotted black line. In addition, the corresponding distribution of Q shows that the quality of the detection tends to be mostly below 0.8, much smaller than for larger values of S/N_base where Q is usually closer to 1.0. Similar behaviour can also be seen in the p-value distributions where the higher the S/N_base, the larger the concentration of results with small p-values.

The p-value is usually the only parameter used in the selection of detections. In fact, for S/N_base > 1, we find that almost all missed detections have Q > 0.5, so the p-value is the parameter ruling the selection and Q has almost no influence in the selection process. However, for S/N_base ≤ 1 there are many cases with p-values smaller than 5% where Q is used to filter the weak detections, as seen in Fig. 7.

Of course, when applying the detection method to real Gaia data, we can only estimate the signal-to-noise ratio based on the sinusoidal fit from GLSP. This is the reason why, in addition to the analysis of the detection performance based on the base S/N_base values, it is also interesting to study the results based on the post-fit S/N $(\widehat{\text{S}/\text{N}})$, which is indeed different from the S/N_base. As shown in Fig. 9 we see that the number of detections increases for the largest values of S/N_base, as in Table 1.

By combining all the results from the 6000 time series, we can simulate the situation where we do not know the original parameters of the time series and have only the post-fit signal-to-noise ratio. Then, we separate them in intervals of mean $\widehat{S/N}\pm 0.1$ and check the statistics, as it was made previously for the S/N_base classes. The results are presented in Table 2. From the comparison between the analysis based on S/N_base classes presented previously, and the one based on the post-fit $\widehat{\text{S}/\text{N}}$, these results indicate that there are no major differences in the statistics of detections between both approaches, especially for $\widehat{S/N}\ge 1.5$.

It is interesting to note from these results that in the scenario where there are signals in all of the time series analysed, but we do not know the S/N_base, $\widehat{\text{S}/\text{N}}\le 1$ lead to p-values and Q factors that often fail to pass the selection test. Meanwhile, for $\widehat{S/N}\ge 2$ we find more than 85% of the detections around the correct period with a rate of incorrect/correct lower than 10% and we miss less than 3% of the detections.

We conclude that a $\widehat{\text{S}/\text{N}}$ of about one appears to be a limit for the detection regime in this setting. For increasingly stronger signals, with larger $\widehat{\text{S}/\text{N}}$, we have an increasingly smaller rate of incorrect detections. The method succeeds accordingly at classifying the sample as detection, at recovering the correct period in the data, and at providing a good approximation of the real S/N.

Fig. 6 Flowchart showing the period search and statistical selection processes discussed in Sects. 3 and 4. It represents part of the full detection process, which is coupled with physical characterisation tests that are discussed in the following sections.

Fig. 7 Distribution of pairs (p, Q) for each of the 6000 sets of synthetic data simulated. The partially transparent dots represent the cases where the detections do not meet the requirements and were disqualified, while the remaining dots represent the detections selected. Of the 6000 simulations, 1303 have p-value>0.1. Hence, their Q factors were not calculated and the corresponding points are not shown in this plot.

4.2.2 Noise-only scenario

The purpose of this analysis is to understand how the method performs when there is no signal in the data and to assess its robustness in handling data with no wobble at all. For that purpose, we generated 6000 sets of simulated AL residuals by replicating the error bars and time sampling of original AL residuals from Gaia asteroids. These synthetic noisy time series were generated by drawing random samples with Gaussian noise perturbation corresponding to the target error bars and we added in each time series a random offset consistent with the Gaia time series, cf. Appendix A.

As discussed in Sect. 4.1, with a p-value threshold of 0.05 it is by definition expected that 5% of the time series are selected (on average) if none of the data present a real signal. Our results showed that in the considered noise-only scenario, the statistical selection process selected as detections about 4.8% of the time series analysed.

The left plot in Fig. 10 shows the distribution of (p, Q) pairs for the 609 time series with p-value < 0.1. The remaining 5391 time series present p-value > 0.1 and therefore their Q was not calculated. We see that the resulting dispersion in the (p, Q) plane is similar to that observed in Fig. 7 for classes of S/N_base ≤ 1. This is in contrast with classes of larger S/N values for which the points tend to cluster around the top left corner.

Turning to the right panel of Fig. 10, we see that noise-only data lead to post-fit S/N values distributed around $\widehat{S/N}\approx 1$, with few cases having $\widehat{S/N}\ge 2$. As expected, this distribution is compatible with the lowest S/N classes of Fig. 9.

In a nutshell, we conclude from this two-fold performance study that, when a wobble is indeed present, our statistical selection process behaves as expected and increasingly succeeds at retrieving this wobble as its amplitude increases. In the situation when there is no wobble, the incorrect detection rate is controlled and does not exceed the prescribed value (set to 5% here).

Finally, the S/N value of one appears to be a turning point in the detection regime, post-fit S/N values larger than ≈2 being a serious hint of the presence of a population of real wobbles⁴.

4.3 Period search results overview

In this section, we apply the detection method to the AL residuals data set. Figure 11 visually illustrates the results of the processes of selection and parameters determination for the case of an object that satisfies the requirements (at least one window with ten observations in less than 10 days, p-value < 0.05, Q > 0.5 and estimated wobble amplitude A smaller than 20 mas) and was selected as a possible candidate by the period search. In this case, for (3457) Arnenordheim, we have a window with 17 observations spread over 2.75 days (as shown in Fig. 3, right panel). The periodogram shows a large peak (Fig. 11 top left) that provides the best fit in the observations with a period of 45.95h (Fig. 11 top right). The MC analysis leads to a very low p-value (< 0.001%, Fig. 11 bottom left) and a quality factor Q of one (Fig. 11 bottom right).

The distribution of estimated periods and wobble amplitudes  for the 30 030 candidates (Fig. 12) shows that small amplitudes are preferentially discarded, as expected since the signal-to-noise ratio tends to be small, but the general distribution remains overall the same. We notice a kind of “tail” in the distribution for the large period and amplitude values, and they are mostly not selected as detections. This phenomenon might be caused by unmodelled low-frequency trends (see Sect. 4.1).

It is also noticeable in Fig. 12 a concentration of detections around 3, 4, 6, and 12 h with larger amplitudes, even though there is no over-representation of such periods. That happens due to the typical frequencies present in the Gaia scanning law (time interval between the two FOVs, rotation period of 6 h, precession of about 2 months, and combination of them). In some cases, the sampling of the data set has a dominant effect on the period search causing the GLSP to find the best period at, or very close to, values of 2, 4, 6, and 8 cycles per day. However, some of the detections found at such values might be real, so we leave at the next selection steps the task of cleaning the spurious detections that are not physically meaningful. In future improvements of the method, we intend to account for this effect. Finally, we note that the apparent vertical alignment of points for small values of A is due to the bin size in the determination of the mean value estimate, which is more noticeable due to the logarithmic scale of the plot.

We find that 83.8% of the windows of data analysed lead to non-detections due to a p-value larger than the 5% threshold and about 6.6% are disqualified because Q < 0.5 even though the p-value is within the limit. Eventually, only 9.6% of the sets of AL residuals comply with the requirements and are selected as possible binary candidates by the statistical selection, which represents 3038 asteroids, of which 40 have estimates from two (or more) different windows of observations.

Some results from different windows are very similar, but most of them show considerably divergent values. The difference between the solutions is typically up to ~35 h in the estimated period and up to ~2 mas in amplitude as can be seen in Fig. 13. At this point, we are not able to decide which window has the most physically meaningful results. Therefore, an asteroid that has been estimated with different periods and amplitudes in different windows at this stage will each have a pair of possible solutions (period, amplitude) considered in the next stage, where the least physically consistent results will be rejected.

Fig. 8 Distribution of periods from best fit (top row), quality factors (centre row), and p-values (bottom row) as percentages of the detections selected with p-value < 5% and Q > 0.5. The different coloured bars represent the different classes of S/Nbase values explored. The dotted line at 23 h represents the nominal (true) period of the injected signal.

Fig. 9 Distribution of time series selected as detections (coloured bars) and discarded (grey) as a function of the post-fit S/N, for the six different S/Nbase classes explored.

4.4 Selection based on physical properties

For the 3038 Gaia asteroids whose data passed the statistical detection test above, we move forward to check if the measured wobbling is compatible with the fundamental physical properties expected for an orbiting satellite. Under the hypothesis that the periodic signals correspond to the orbital period, we can exploit known properties measured by other surveys, such as the object size, to infer the range of size ratios and densities that could generate the observed signal. Density, in particular, should be within acceptable ranges, and compatible with the expected mineralogy for a corresponding spectral type.

We start with Kepler’s third law of planetary motion

$$\frac{T^2}{a^3}=\frac{4{\pi }^2}{M_{\text{total}}G},$$(8)

where M_total is the total mass of the studied system (M₁ + M₂), T is the orbital period of the secondary body around the primary, which we assume to be the same as the wobbling period, and G is the gravitational constant. We manipulate the equation to express it as a function of the mass ratio of the system $q=\frac{M_2}{M_1}=\frac{D_2^3}{D_1^3}$, the diameters of the primary and secondary component (D₁ and D₂, respectively), and the bulk density ρ.

Furthermore, we also consider that the currently observed diameter corresponds to a sphere of surface brightness equivalent to the sum of the binary components, $D^2=D_1^2+D_2^2$. Hence, the total mass of the system can be written as

$$M_{\text{total }}=\frac{\rho \pi }{6}\frac{D^3}{{\left(1+q^{2/3}\right)}^{3/2}}(1+q).$$(9)

Finally, substituting Eqs. (9) and (1) into Eq. (8) and with some manipulations, we obtain a relation describing the density ρ as a function of the mass ratio q:

$$\rho =\frac{24\pi {\alpha }^3}{G{T}^2{D}^3}\left|\frac{{\left(1+q^{\frac{2}{3}}\right)}^{{}^{\frac{9}{2}}}{(1+q)}^2}{q^2{\left(q^{\frac{1}{3}}-1\right)}^3}\right|,0<q<1.$$(10)

By this expression, the dependency of the density from the mass ratio q can be evaluated, by assuming an estimation of D, once T and α are measured from the astrometric wobbling. We retrieve D from literature using the SsODNet service⁵ (Berthier et al. 2023), along with the associated uncertainty. If SsODNet does not provide the diameter for a given object, we estimate it from the absolute magnitude H by taking the mean diameter for the geometric albedos p_υ spanning the interval [0.06, 0.2]. We assume a standard deviation of 15% for convenience, considering that without this assumption the errors in the diameter estimates are larger than 50%. This restrictive choice is partially compensated by the generous density interval that is allowed, as explained in the following.

While T is a direct result of our signal analysis, the measured amplitude α is only a minimum value, due to the unknown projection factor of the real wobble in the AL direction. In this respect then, having fixed all other parameters, Eq. (10) provides a minimum density as a function of q.

The Eq. (10), as shown in Fig. 14, is a concave function with a minimum corresponding to the maximum wobbling amplitude in Eq. (1). This minimum represents the smallest ρ of the system compatible with the observed wobbling. The uncertainties are obtained with the linear propagation of the errors in D, T and α.

For a first selection, we consider a large density range ρ ∊ [0.3, 7] g cm⁻³. We are interested in selecting the possible ranges of mass ratio q corresponding to densities in this interval. Therefore, when the minimum of the density function is <0.3 g cm⁻³, we obtain two separated intervals delimited by 4 different values of mass ratio ((q_min1, q_max1) and (q_min2, q_max2)). From them, and the associated values of ρ, the corresponding values of the total mass can be computed and used in Eq. (8) to derive two different intervals of possible separation of the binary ((s_max1, s_min1) and (s_min2, s_max2)). This is the case of (7616) Sadako in Fig. 14.

In cases where the minimum density estimate is never smaller than the inferior limit, as for (6261) Chione, we get a single interval of q and possible separations. Eventually, if the function has a minimum value larger than the adopted upper threshold, the wobble signature is considered not compatible with a binary object as this would require an excessive density. By adopting the density threshold mentioned above, 934 objects are selected and their possible ranges of q and s have been computed.

A further step consists in considering the smallest separation, for each candidate, compatible with an orbiting companion. We find that for some of the objects, the range of separations partially overlaps with the Roche limit⁶, computed by assuming the same density for the two components. Binary systems inside the Roche limit are unlikely to exist, and the detected wobble could instead be caused by a very irregular shape. Of course, more accurate physical data, in particular diameter measurements, could change the results for some objects and put them beyond the Roche limit. To be conservative, for this first search, we have decided to reject candidates that have separation intervals overlapping, even partially, with the Roche limit. With this selection, we remain with 362 wobble measurements for 358 objects.

As an attempt for a further verification, we considered the largest selected objects in the table having a known shape in the DAMIT⁷ database (Durech et al. 2010) to compute the distance between photocentre and centre of mass, namely for the asteroids (476) Hedwig (Marciniak et al. 2023), (487) Venetia (Marciniak et al. 2018), (516) Amherstia (Durech et al. 2009; Hanuš et al. 2011), (532) Herculina (Kaasalainen et al. 2002; Hanuš et al. 2017), (542) Susanna (Hanuš et al. 2021), (556) Phyllis (Marciniak et al. 2007), (699) Hela (Marciniak et al. 2012), (1200) Imperatrix (Ďurech et al. 2020), (1639) Bower (Ďurech et al. 2019), (2219) Mannucci (Ďurech et al. 2020), (70045) 1999 DA5 (Ďurech & Hanuš 2023). We simulated the photocentre shift due to the illuminated shape alone, at all epochs used for our period search, by using Lambert’s emmition and Lommel-Seeliger’s scattering laws (Lambert 1760; Pedrotti et al. 2017; Fairbairn 2005) and the corresponding illumination geometries. Then, its projection on the AL direction was used to look for correlations with the observed residuals. This test reproduces the one reported in Tanga et al. (2023) for (21) Lutetia.

We repeated the test with the multiple pole and shape solutions available for each object, also checking the possible role of a rotation phase shift concerning the nominal rotation origin. Our results show that the amplitude of the photocentre displacement could be correlated to the residuals for (487) Venetia, (1200) Imperatrix (for one of the two pole solutions available) and (2219) Mannucci (in one of the two observation windows that we consider), as shown in Fig. 15. It is also anti-correlated for (532) Herculina and (556) Phyllis. In the other cases, no correlation shows up.

This preliminary test is very sensitive to several details such as shape and pole coordinate errors, scattering law, presence of concavities, etc, and should be taken with caution. However, we notice that the order of magnitude of the wobble amplitude due to a (single) shape would be about the same as the signature produced by a satellite for large asteroids (separations of ~200 mas or more, spatially resolved by Gaia). However, we notice that the amplitude of the photocentre shift is not dissimilar from our residuals. In a synchronous system (primary rotation equal to secondary orbital period) they could even add up with the same or opposite sign, depending on the phase angle of illumination. Only in this case, when the sign is opposite, this could result in the observed anti-correlation which would then be evidence supporting the satellite’s presence.

All these considerations point in the direction of not excluding the largest objects remaining on our list. They certainly deserve more accurate characterisation by future observations and more extensive, complex modelling efforts when more data becomes available.

Fig. 10 Resulting distribution of parameters from the 6000 simulations in the noise-only scenario. Left: distribution of pairs (p, Q) for the time series simulated. 5391 of the cases have p-value > 0.1 (for these cases the Q factors were not calculated and the corresponding points are thus not shown). The grey dots represent the cases where the data was disqualified as a detection (cases for which p > 0.05, Q < 0.5, or both) and the red dots represent the detections. The two grey points in the top left square are detections discarded due to  > 20 mas. Right: post-fit S/N distribution for the pairs in the left plot.

Fig. 11 Visualisation of the output of the main steps of the selection process in the period search for one candidate example ((3457) Arnenordheim). The top left plot shows the periodogram from the fitting of the 17 observations in the window. The peak represents a period of 45.95h and the best fit is shown in the phased data in the top right plot. The distribution of data points not phased is shown in the right plot in Fig. 3. The bottom left plot shows the empirical distribution of GLSP power V (maximum value of the GLSP under ℋ0) and the vertical dashed line shows the position of the value of υ, the power of the peak from the periodogram with the real data, with p(υ) < 0.001%. Finally, the bottom right plot shows the distribution of best frequencies from the 5000 MC simulations under ${\widehat{\mathscr{H}}}_1$). The distribution tends to accumulate around the originally estimated frequency, showing that this detection is robust (here Q is close to 1.)

Fig. 12 Distribution of estimated wobble amplitudes and periods obtained from the period search for all objects. In grey are the time series disqualified due to a p-value > 0.05; in blue are the cases disqualified because Q < 0.5, but with p-value within the limit of 0.05; and in red are the cases selected as binary candidates from the period search for complying with p-value< 0.05 and Q > 0.5.

Fig. 13 Distribution of differences between the solutions for the 40 objects with multiple windows of consecutive observations.

Fig. 14 Examples of the bulk density ρ as a function of the mass ratio q for two asteroids, (7616) Sadako (top left) and (6261) Chione (top right). The red line indicates the nominal value. Shown in grey is the error range on the density. The blue and orange horizontal lines represent physically reasonable limits. qmin1, qmax1, qmin2 and qmax2 mark the corresponding extremes of the mass ratio values. The bottom plot shows the case for (930) Westphalia, eliminated in the density filtering because the minimum estimated density is higher than the chosen range of bulk density.

Fig. 15 Comparison between the average DR3 astrometric residuals and estimated average photocentre offset for candidates (487) Venetia, (1200) Imperatrix and (2219) Mannucci. The top plots show the variation of the data over time in the along-scan direction, while the bottom plots show the correlation between the two sets of data.

4.5 Constraints from taxonomy

The average density value for each taxonomic class can help to further constrain the physical properties of objects in our selection. We retrieved the taxonomic class for all objects in the list of preliminary candidates using SsODNet. Of the 358 candidates selected from the previous steps, 264 have their classes determined, corresponding to one of the complexes for which an average density is provided in Carry (2012).

By using the average density for the taxonomic class, we can repeat the analysis illustrated in Sect., recompute the permitted intervals of possible physical parameters, and refine the selection of the sample. We reached a total of 156 candidates. In addition, since large real wobbles tend to produce high $\widehat{\text{S}/\text{N}}$ (cf. Sect. 4.2) we flag as best candidates those with $\widehat{S/N}>1$.

The objects of our sample on which taxonomy-derived density constraints are applied are presented in Fig. 16, which shows the estimated radius of the primary as a function of the allowed separations (in units of the radius of the primary) for the 156 objects. The purple and pink points represent the first and second intervals of possible mass ratio values. Those with dark borders are the best candidates with $\widehat{S/N}>1$, and the points with grey borders are the candidates with $\widehat{S/N}<1$ that are not in the list of the best cases but still passed through all the filtering steps.

5 The final sample of candidate binary asteroids

The total outcome of 358 objects from this first search of astrometric binary candidates appears in the full version of Table 3 in CDS. It includes 6 known binaries, 67 of which (highlighted) present results with the density constrained by taxonomy and $\widehat{S/N}>1$.

Table 4 summarises the number of candidates after each step of our selection process. We select the 30 030 asteroids in DR3 with consecutive observations in a short time. From the statistical selection, we obtained 3038 preliminary candidates that presented p-value< 0.05 and Q > 0.5. The validation on the physical ground and the selection based on known taxonomic type further restricted this list to 156 objects presenting characteristics compatible with binary systems. Most of them have an estimated $\widehat{S/N}<1$, as the regime of the scarceness of the data and potentially low wobble amplitudes make detectability difficult. Therefore, we can consider that those with $\widehat{S/N}>1$ among the 156, are the best characterised. The distribution of their parameters concerning the statistical thresholds of period determination is shown in Fig. 17.

In the final list of 358 candidates, there are also 31 objects whose estimated period is not very reliable. In fact, their values reach is towards the highest ones in the allowed range of the GLSP (reaching three times the length of the observation window).

An example is (1509) Esclangona, a known binary with an estimated period of ~410 ± 120 h (assuming a bulk density of ρ = 2.0 ± 1.0 𝑔 cm⁻³, Marchis et al. 2012) well beyond the periodogram range (left panel of Fig. 18). A poorly defined best period is found, in a plateau of the periodogram where a whole range of values is equally possible. Going towards the extreme of the highest periods tested, the fitting function to the residuals (covering 84 hours right panel) tends to locally approach a straight line, resulting in a whole range of similar fitting quality.

It should be noted that the clear trend could indicate a genuine long period. Hence, the data may contain part of a periodic signal that we cannot consistently estimate, as this would require a longer data sequence. In this case, and the other ones resulting in a similar periodogram, we are only able to give a minimum boundary and estimate that the period should be larger than the length of the observation window. Candidate binaries in this category are identified by an orange background in the final list of candidates in Table 3.

Additionally, we notice that the method seems to favour the selection of objects presenting periods longer than 10 h, which is perhaps the sign that high frequencies not captured by the constant plus sinusoid model of the GLSP are present in the time series. By filtering the physically meaningful candidates we eliminate most of the cases with p-values > 2% and Q < 0.9 and the best candidates are mostly those with very low p-value and Q close to 1.

The comparison between Figs. 19 and 10, showing the heavier right tail of the $\widehat{\text{S/N}}$ distribution obtained from Gaia data, further strengthen the evidence that our candidates are indeed binary asteroids. Two of the best candidates in the list, (2219) Mannucci and (542) Susanna present matching wobble periods from two different time windows, separated in time by weeks or months. At such large intervals, the orientation of the wobbling motion relative to the scan direction can be different (as confirmed by the different wobble amplitude), thus making these measurements independent. The fact that the two estimated wobble periods of each object match is another sign that our method of period search is rather robust.

The bodies selected as best binary candidates, including the known binary asteroids, are in large majority main-belt asteroids (Inner MB: 11, Middle MB: 26, Outer MB: 22), but 3 Phocaea, 3 Mars-Crosser, 2 Hungaria, 1 NEA Amor, and 1 Trojan are also present. All of our best candidates present separations <200 mas, below the limit at which the onboard Video Processing Unit of Gaia treats sources separately.

Fig. 16 Plot of the radius of the primary as a function of the separation intervals normalised with the Roche limit for fluid satellites, for the sample of 156 asteroids with known taxonomy. The squares represent the mean value of the first mass ratio interval estimated and the circles represent the second interval. The coloured symbols with dark borders are the objects from the taxonomic filtering with $\widehat{S/N}>1$, and the grey symbols are those with $\widehat{S/N}<1$.

5.1 Known binaries

Our results contain 6 known binaries: (1509) Escanglona, (5817) Robertfrazer, (2871) Schober, (4337) Arecibo, (317) Roxane and (18301) Konyukhov. Such a small number of known binaries is not unexpected, since our filtering and statistical validation reduce a large sample to a very small population. For each binary, the probability of showing an astrometric signature in Gaia DR3 is not very high, as it requires the appropriate orientation of the wobble motion concerning the scan motion, during the periods of consecutive observations.

Table 5 contains the parameters described in the literature and obtained from observations of the six mentioned known binary systems. These six known binaries offer in principle the possibility of validating our approach, but in practice, other limitations appear. As mentioned previously, (1509) Esclangona cannot be exploited for the comparison, as it has a small satellite with an orbital period well beyond the reach of our data sample. The most probable periodicity given by our periodogram, 123.71 ±11.09 hr, is roughly between 2 and 4 times smaller and should be considered as a lower limit (Sect. 5).

In the case of (5817) Robertfrazer the data from Gaia DR3 presents two windows of observation that we use, but the result obtained from one of them does not pass the physical validation process and is therefore discarded. So, for the remaining window, we also find a matching size ratio (k = D₂/D₁), however, our period is more than 5 times larger than reported (~28.862 hr), and again close to the upper limit of the frequency range of the GLSP.

Interestingly, however, for both these incorrect matches the apparent wobble signal is strong, with a $\widehat{\text{S/N}}$ of 6.09 for Esclangona, and 2.64 for Robertfrazer, and the wobble amplitudes are significantly high (5.09 ± 0.54 mas for Esclangona and 1.61 ± 1.19 mas for Robertfrazer), especially for Esclangona. A possible alternative interpretation could be the presence of another satellite in the system, and future Gaia astrometry will probably help to clarify which scenario is more probable. In both cases, we estimate that the interest in selecting these objects by our procedure remains.

For (2871) Schober, the size ratio from the literature matches our results. We find solutions from two different time windows of observation, with wobble periods of 20±4 hr and 50 ± 11 hr: the first can be considered an alias emerging from the data at about double the frequency due to the irregular sampling from Gaia, while the second results in intervals of separation that overlap partially with the Roche limit and, therefore, it was rejected. Both windows are compatible with the value reported in the literature, obtained by photometry (Benishek et al. 2023a).

In the case of (317) Roxane, we find that our results also match well with the size ratio from the literature. We also estimated the orbital period of the secondary based on the separation, as was made for Esclangona, and we found that our period estimate is about half of the expected value. The observation window from Gaia DR3 that was selected to perform the search is shorter than the expected period of the secondary, which once again limits the extent of our period search and invisibilises the method to find the correct period.

For (18301) Konyukhov we find that the observed size ratio of 0.26 (≈0.017 in mass ratio, Benishek et al. 2023b) is very close to the upper limit of the first interval of mass ratio of 0.013, translated into an upper value for size ratio of approximately 0.24. Regarding the period, our Monte Carlo procedure to estimate the confidence interval on this value reveals that the distribution of the estimated periods (cf. Sect. 4.1) is nearly bimodal, with two modes on 70.1 h and 140.6 h. As a consequence, the procedure provides a very large confidence interval (±70.02h), which reflects the possibility that the true period can actually be close to either of the two values.

For (4337) Arecibo we find that our estimated period is equivalent within the uncertainty to other sources (Tanga et al. 2023; Ďurech & Hanuš 2023). This binary has two different size ratio values from the literature. Stellar occultations Gault et al. (2022) provide k_occ = 0.533 ± 0.063, while the fit to both Gaia astrometry and occultations Tanga et al. (2023) results in k_Tanga = 0.35 k_occ, with indications of a probable ellipsoidal shape of the components. Using the best value of surface–equivalent diameter from SsODNet (19.746±0.324 km) and our estimates for wobble amplitude (0.61±0.07 mas) and wobble period (35.9284±6.3600 h), we obtain results close to the parameters of the system previously published (lower size ratio solution in Fig. 20). The remaining difference with respect to Tanga et al. (2023) arises from the different estimates of the wobble amplitude and period, and from the fact that we do not consider here additional constraints such as those coming from the observation of stellar occultations.

5.2 Comparison to the known population

Our candidate binaries can be compared to the existing population. As our observational biases are different from other techniques, the goal is not to reproduce the same distribution. In principle, we can discover asteroids with satellites whose dynamical properties are unusual, in particular in a domain of intermediate sizes (that we can conventionally define between 20 and 100 km) where traditional techniques have very limited capabilities of discovery. On the other hand, we can also reasonably expect that some of our candidates are similar, in physical and/or dynamical properties, to the binaries known so far.

A first comparison between the astrometric binary candidates and the known binary population is based on taxonomy. Our goal is to verify if an excess of S-type binaries is present in our sample, by an approach similar to Minker & Carry (2023). We start by grouping the taxonomic classes in large complexes, as shown in Table 6. Notably, P-type asteroids are in the C complex, not the X complex. The reference sample of asteroids was constructed based on a magnitude limit of absolute magnitude H < 16 and within ranges of a ∊ [2, 3.5] au, i ∊ [0°, 30º] and e ∊ [0, 0.5] roughly matching the distribution of the sample of binary candidates.

The reference sample is selected by choosing 30 asteroids around each binary candidate, from a partition of the 4D (a,e,i,D) parameter space (see Minker & Carry 2023 for details).

Figure 21 shows the comparison between the reference sample and our binary candidates with D < 20 km, which represent most of our candidates. The S-complex taxonomic types present a fraction above 1, which means that there is an over-representation of objects of this complex in our sample of small binary candidates. On the other hand, our small candidates show a deficit of C-complex taxonomic types in comparison with the background asteroid population.

The results from this analysis are in great agreement with Minker & Carry (2023), where it is shown that even though chondritic asteroids are the most common type in the Solar System, there is a preferential presence of S-type over C-type, among binary objects at small diameters. Thus, we can conclude that the taxonomic distribution of our candidates is comparable to the one of the previously known population.

Taking a closer look at the sub-sample with 30<D< 50 km, we find no dominance of any taxonomic type (Fig. 22), suggesting a disconnect between the small- and medium-sized binary asteroid populations. However, statistical comparisons to a known population are not possible, as the only known binary asteroid in this size range is the system (243) Ida – Dactyl. Interestingly, it was discovered by a spacecraft during a fly-by (Chapman et al. 1995). This could indicate that un-conventional techniques such as the astrometric method are necessary to find similar binaries. The limited presence of large (D>50 km) asteroids in our sample of binary candidates eliminates the possibility of continuing this taxonomic analysis at larger sizes.

In order to compare the dynamical properties of our candidates to the known population, we first present (Fig. 23) the distribution of the primary rotation period as a function of the estimated primary diameter. The rectangles in the plot represent a raw approximation of the known binary asteroid groups (Pravec & Harris 2007; Pravec et al. 2010; Margot et al. 2015). Group L is for fast-rotating and large objects. Group A represents small primary asteroids with fast rotation; B, systems with small primaries but size ratios D₂/D₁ ~ 1 that rotate slower as the primary diameter increases.

The first striking feature of our sample is the partial overlap with the pre-existing population discovered by photometry (Group A), which clearly appears below ~ 20 km. Some candidates also partially overlap group B, but their mass ratios are small in general (of the order of 0.1). At the other extreme of the size range, binary candidates from our search do not match group L (besides one). This is expected, as the wobbling signal from large objects should not be dominated by a small satellite, but by the photocentre displacements on the shape of the primary alone. Coherently with this interpretation, some potential candidates in this region are excluded due to the separation being too small relative to the Roche radius. Eventually, intermediate sizes (20–100 km), clearly underrepresented in the known sample, are present for our astrometric candidates, directly confirming the new potential of this technique.

The comparison of satellite orbits and sizes is more difficult and intrinsically ambiguous, as our analysis introduces two families of solutions, corresponding to two different possible ranges of mass ratio, as shown in Fig. 14. The estimated mass ratio intervals for the 67 selected best binary candidates, and their respective primary rotation periods, are shown in Fig. 24. The two mass ratio intervals are more or less separated depending on the minimum estimated density for the system. At this stage, we have no criteria to favour a solution over the other. However, we notice that objects of intermediate size populate preferentially extreme values of the mass ratio. This could imply that binaries of similar-sized components could exist in that population.

The analysis of the ratio of the wobble (orbital) period to primary rotation period, reveals the presence of synchronous objects and, to a minor extent, of objects in a 1:2 period ratio (Fig. 25). Among the known main belt binaries, almost all those with D_s/D_p > 0.5 are synchronous, as the tidal coupling is more efficient for large satellites. Applied to our sample, this evidence would imply that large mass ratios have a higher probability of showing orbit-spin coupling. Arecibo, at D_s/D_p > 0.53 and synchronous, clearly falls in this category. Periods on the 1:2 ratio could simply represent an alias and could be considered synchronous too. Interestingly, intermediate-size objects are in large majority in this category.

Fig. 17 Distribution of $\widehat{\text{S/N}}$ vs. p-value (top); estimated wobble period vs wobble amplitude (center); and p-value vs. quality factor (bottom) for the 3038 objects obtained from the statistical selection (in grey), the 156 objects that survived the taxonomic filtering (in yellow), and the 67 best candidates that passed all tests and have $\widehat{S/N}>1$ (in pink).

Fig. 18 (1509) Esclangona residuals distribution (right) and the periodogram from the fitting of the 21 data points in the window (left).

Fig. 19 Distribution of $\widehat{\text{S/N}}$ for our list of 156 binary candidates selected from the taxonomic filtering. The 67 best candidates are those with $\widehat{S/N}>1$

Fig. 20 Comparison between the parameters for the Arecibo system. The triangle with black edges and grey error bars represents the values from literature (Table 5), while the circles represent the first (transparent) and second intervals of size ratio estimated from our method using the same wobble amplitude as in Tanga et al. (2023).

6 Discussion

As explained above, this first attempt of binary search in Gaia DR3 astrometry is performed by adopting some conservative criteria and an approach that does not take into account all the observations available for each asteroid. Other choices are possible, such as those described in Segev et al. (2023), where a solution of the complete inverse problem leading to satellite orbits, and involving all the AL residuals of each object, is proposed. Their search is not successful and the main explanation given is the rejection mechanism of astrometric outliers implemented by the validation pipeline in DR2 (Gaia Collaboration 2018).

However, our results show that periodic fluctuations in the residuals are compatible with the expected amplitude and orbital periods of satellites. We thus think that the missed detection in the cited article is due to three main reasons: the strong increase in the time range and number of asteroids of Gaia DR3 concerning DR2 (by factors 1.5 and 11, respectively); the improved astrometric quality of DR3; and, of lesser importance, the possible presence of systematic errors in isolated astrometric measurements that are included in the analysis.

Concerning the first argument, by simple statistical scaling we can estimate that, if DR2 astrometry were good enough to reveal our best candidates, only a few of them (<10) would have been found. However, the probability of detection is not a linear function of the time range, as a reduction in the observed arc can have a strong adverse impact on the quality of the orbital fit and, in turn, on its residuals.

The fact that DR3 improves astrometric quality in a significant way over DR2, is well visible not only in the overall statistics but also in specific cases such as (4337) Arecibo, whose residuals show a clear trend in DR3 (Tanga et al. 2023) but less clearly in DR2 (Fig. 19 in Segev et al. 2023), over the same time span. Hence, this is the first work that provides a large list of binary candidates obtained from astrometry and crosses dedicated statistical and physical selection criteria. At this stage, without independent confirmation of the binary status and a clear discrimination between the two solutions that we obtain for each candidate, the impact of our possible discoveries on formation scenarios cannot be very specific. However, some general considerations can be attempted.

First, the intermediate size range (20–100 km) where very few binaries are currently known (only ~5%n of the total population) is probably the most interesting, and the one having a stronger potential impact. If they are synchronous and have large mass ratios, a well-known prototype at large sizes can be considered the asteroid (90) Antiope (size ratio ~0.95, diameter ~88 km, Merline et al. 2000). Even though there are several theories for the formation of the Antiope system (Pravec & Harris 2007; Weidenschilling et al. 2001), observations support the hypothesis that the Antiope system objects were formed simultaneously from the disruption of a parent body(Descamps et al. 2009; Marchis et al. 2011).

With a similar mechanism, as shown in numerical simulations (Doressoundiram et al. 1997; Durda et al. 2004), Escaping Ejecta Binaries (EEBs) result from the creation of collisional fragments at low relative velocity, that remain gravitationally bound. They preferentially have highly eccentric orbits, that could be circularised by tidal forces in their secular evolution. They should not have particularly fast primary rotation as they are not created by fission or fragment ejection. However, the current sample of known binaries, if simulations are correct, has a clear under-representation of this category. The only good candidate is (317) Roxane with its satellite Olympias (Drummond et al. 2021). The population of EEBs could also be compatible with binary Mars impactors revealed by Vavilov et al. (2022). Once validated, our sample could clearly contribute to revealing many details about the presence, origin and evolution of EEBs.

We also notice the compatibility of our highest mass ratio solutions with the findings by Scheeres (2007) and Pravec et al. (2010), that determine (theoretically and by observations) the conditions for splitting of pairs created by fission. Their main finding is a general tendency to create separate pairs when the mass ratio is <0.2, which is about the inferior threshold for our distribution of solutions for high values of mass ratio when small asteroids <20 km are considered (magenta circles in Fig. 24). Some of the systems close to (and above) that critical splitting threshold, could have been produced by fission. Considering that mass shedding and re-accumulation in orbit is the most likely formation process for small asteroids (D ≲20 km) (Pravec & Harris 2007; Walsh et al. 2008), we expect that future observations characterising our sample should find results compatible with the known population in this size range.

Fig. 21 Comparison of the fraction of small objects with D < 20 km in each taxonomic complex between the reference sample of background asteroids and the population of our astrometric binary candidates. A taxonomic type with no column indicates that asteroids of that type are present in the reference population but not in the group of binary candidates.

Fig. 22 Fraction of binary candidates with 30< D <50 km in each taxonomic complex, with respect to the reference sample of background asteroids.

Fig. 23 Plot of the primary rotation period vs. the estimated primary diameter for the best candidates. The black dots are known binaries in the shown intervals and the magenta circles are the best binary candidates obtained.

Fig. 24 Plot of the primary rotation period vs. the estimated mass ratio intervals for the 67 selected candidates. The lighter and darker coulored circles represent, respectively, the first and second intervals of possible mass ratio values for the best binary candidates.

Fig. 25 Plot of the ratio of the estimated orbital period (or wobble period) of the secondary over the photometric rotation period of the primary vs. the estimated size ratio of the binary system. The black dots represent all of the known binaries and the red dots highlight those with intermediate sizes (20 km < Dprimary < 100 km). The light and dark purple circles represent the average value of the best-estimated size ratio intervals for objects smaller than 20 km and the light and dark blue circles represent the larger objects. The blue line shows the condition where the estimated orbital period of the secondary is almost equal to the rotation period of the primary object, representing a 1:1 period ratio, while the orange line represents the 1:2 period ratio.

7 Conclusions

Our work presented a simple but robust method to detect periodic signals in the residuals of the orbital fit in Gaia astrometry, which is compatible with the presence of binary asteroids. We explored all astrometric data of asteroids available in the DR3 catalogue and, after a careful selection process combining statistical and physical analyses, we obtained a consolidated list of candidates that are likely to be binary asteroids.

Among the strengths of our results, we can highlight that the presence fraction of our candidates has physical properties similar to the known binary population (small asteroids with companions discovered by photometry), but also spans a parameter range that tends to fill the gap where other techniques have been failing to find binaries, as expected from preliminary studies (Mignard et al. 2007; Pravec & Scheirich 2012). The distribution of the taxonomy of our candidate binaries also provides interesting matches to the known population and some potentially new evidence. We also find a large fraction of possible synchronous binaries, which are notoriously difficult to discover by photometry, especially if the secondary rotation is also synchronised (doubly synchronous binaries).

Some limitations of our approach are also clear. As we deal with an astrometric accuracy of the order of 1 mas, we cannot exclude that some signals are not due to companions, but to the photocentre displacement occurring during the rotation of irregular, single asteroids. This could be the case, especially for the largest objects approaching the 100 km diameter, although the relatively small phase angles should mitigate the effect. However, our preliminary analysis of the detection of photocentre variation shows that even though for some cases it could be the source of the signal detected, the results are not significant enough to discard the possibility of the astrometric detection of a satellite.

We also note that our comparison with known satellites is limited to a very small sample, but seems to point out that in some cases we detect companions, although not with the parameters derived by previous surveys. The parameters estimated by making use of taxonomy must be considered with caution since they could be misclassified, which would significantly change the separations. Additionally, we face problems due to the limitations of the data. For instance, the length of the observation window can be shorter than the wobble period, the irregular sampling that causes the rise in aliases in the period determination, and the usually small number of data points along with low S/N that complicates the detections.

Future improvements may involve many aspects: the exploitation of the astrometry could adopt a more refined error model at the transit level; the inclusion of observations more distant in time, with an increase in the accessible period range; a more complex model for the astrometric signal, involving flattening of the components and non-uniform illumination; a more detailed analysis of the spectral properties and the choices of sizes and plausible densities. In addition, a full simulation of the whole process starting from a synthetic binary population and simulated Gaia astrometry is a complex task that should provide a better view of the biases present in our method, and its limitations. However, such extensive investigation is beyond the scope of this work.

Other Gaia data releases, such as the Focused Product Release (which appeared in October 2023) and, in particular, the DR4 (including all the asteroid astrometry from the nominal mission) can certainly offer opportunities to consolidate our approach and refine the results. We expect the astrometric quality to increase further, and the number of consecutive data sequences (also for the same asteroid) to be more frequent.

Despite these limitations, we think that we have revealed the capabilities of Gaia in discovering asteroid satellites by the astrometric method, and we have set the foundations for future, more refined searches. Stellar occultation campaigns focused on our candidates may pin down the orbits, sizes, and shapes of the companions. In the absence of our Gaia binary candidates, a blind search of asteroid companions by stellar occultation observations would be unlikely to succeed. In addition, without the stellar occultation technique, it is hard to better assess some of the candidate system properties.

For some of the candidates with similar characteristics to the known binary population discovered by photometry, light curves could be sufficient to determine the presence of a companion. Therefore, we warmly invite the community to collaborate on the campaigns that we are organising, by photometry and stellar occultations, to validate our findings and physically characterise the Gaia astrometric binaries.

Acknowledgements

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/Gaia. The Gaia archive website is https://archives.esac.esa.int/Gaia. This work was supported by the project Gaia Moons of the Agence Nationale de Recherche (France), grant ANR-22-CE49-0002. It was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) -– Finance Code 001, also by CAPES-PRINT Process 88887.570251/2020-00, by the French Programme National de Plan-etologie, and by the BQR program of Observatoire de la Côte d’Azur. P.B. acknowledges funding through the Spanish Government retraining plan “María Zambrano 2021–2023” at the University of Alicante (ZAMBRANO22-04). We made use of the software products: SsODNet VO service of IMCCE, Observatoire de Paris(Berthier et al. 2023); Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018, 2022); Matplotlib (Hunter 2007); Multiprocess package (McKerns et al. 2012).

Appendix A Amplitude offset

To evaluate the statistical significance of the highest GLSP peak, we explained in Sect. 4.1 that we use Monte Carlo simulations. Their purpose is to generate a large number of synthetic times series that are consistent with our model when there is no sinusoid, that is, an unknown constant sampled at the same time sampling grid as the data set under investigation plus noise consistent with the data error bars.

Therefore, to generate unknown constants (offsets) that are consistent with our Gaia residual time series, we estimated this offset from a large number of time series and plotted its distribution. The result is in Fig. A.1, limited between −2 and 2 mas for visualisation purposes.

The figure compares the empirical distribution with a Laplacian (magenta, parameters µ and b) and Gaussian (red, parameters µ and σ) best fit. Clearly, the Laplacian fit is better and we use this distribution to generate offsets in the Monte Carlo simulations.

Fig. A.1 Histogram of distribution of the estimated offset ${\widehat{y}}_0$ for each window of search explored. The grey bars represent the empirical distribution, while the coloured lines represent different distributions fitted to the data and their respective best-fit parameters.

Appendix B Detection of eccentric system

Highly eccentric binary asteroid systems are rare and not expected to be frequently found in our sample. However, it is interesting to evaluate how efficiently we detect eccentric binaries’ astrometric signal. For such analysis, we start by creating samples of data from synthetic eccentric binary asteroid systems in one dimension, which mimics the projection of the residuals from the orbital fit of the Gaia DR3 astrometric data in the AL direction.

The generic binary system’s parameters were chosen arbitrarily but consistent with what we expect to detect. It consists of a primary with a diameter of D₁ = 40 km and a secondary with D₂= 8km resulting in a mass ratio of q = 0.008. The bulk density of both objects was chosen to be ρ = 2.5 g cm⁻³, a separation of a = 120 km equivalent to 3 radii of the primary and eccentricities of e = 0, 0.2, 0.5, and 0.8. Assuming that it was observed at a distance of 2.8 au, the resulting signal has a nominal wobble amplitude of α =1.8 mas and a period of T=30.5h, as shown in Fig. B.1.

Fig. B.1 Nominal signal from the generic binary system with extreme eccentricity values.

We randomly select 20 data points from the nominal signal over 3.75 days and create 1000 different samples of data with a process similar to the one described in Sect. 4.2, as it is also done for the white noise samples.

For each one of the 5 situations analysed (4 different values of eccentricity + noise only), we then calculate the empirical probability of detection $\left({\widehat{P}}_{det}\right)$ and of false alarm $\left({\widehat{P}}_{FA}\right)$. When plotting ${\widehat{P}}_{det}$ as a function of ${\widehat{P}}_{FA}$ we obtain Fig. B.2, which provides a visual representation of how well our test distinguishes between the two classes (noise only, or noise plus wobble). The plot shows that the detection method is efficient since all curves are concentrated on the top-left corner.

For instance, if we choose to accept a false alarm probability of 5% we find that about 98% of the signals from the circular case are detected against approximately 85% from the highly eccentric systems. Turning to the estimated $\widehat{S/N}$ and estimated wobble amplitude for the case where ${\widehat{P}}_{FA}=5\%$, we obtained the results of Fig. B.3. Here, the considered wobble amplitudes correspond to a nominal S/N of approximately 1.7. We see that the estimated S/N are consistent with this value. We notice that as we increase the eccentricity, the estimated amplitude decreases and the uncertainty grows. Since the periodogram uses a sinusoidal model to fit the signal from the eccentric systems, we lose some precision on the amplitude estimate when the signal comes from an eccentric system, which also impacts the $\widehat{S/N}$ estimate as shown. Still, the statistical detection process succeeds at selecting at least 90% of the samples of data from the most extreme eccentricity case. Note finally that when there is no wobble in the data, the estimated S/N tends to be much smaller (below) than in the wobble-present case.

Fig. B.2 ROC curve for the different eccentricity values considered.

Furthermore, when analysing the frequencies at which the maximum of the periodogram is found, we find that for the circular case about 97% of the detections are within the correct signal frequency± 10% against and 79% for the case with e=0.8. Even in extreme eccentricity cases, when a signal is correctly classified as a detection, the period is correctly estimated 63% of the time. We can conclude that the detection of eccentric binaries is indeed more difficult, but this study shows that the method is still efficient under these conditions. We finally note that highly eccentric binaries are not commonly found and they are also not expected to survive for long periods of time under these extreme configurations (Marchis et al. 2008; Walsh & Richardson 2008; Walsh & Jacobson 2015), so they must be exceptions and thus rarely found in our data.

Fig. B.3 Density histograms for estimated $\widehat{S/N}$ (top) and estimated wobble amplitude (bottom) for the cases with ${\widehat{P}}_{FA}=5\%$.

References

All Tables

All Figures

Fig. 1 Representation of a simplified case where the Gaia satellite would be observing a binary system perpendicularly to the plane of rotation of the bodies, along with the projections of the positions of the barycentre and the photocentre on the AL direction and the representation of the orbital fit residuals’ projection in the along-scan direction. AL direction is represented horizontally for simplicity, but in reality it can take any orientation, slowly rotating by some degrees on a scale of several hours.

In the text

	Fig. 2 Plot of the amplitude of the photocentre wobbling as a function of the binary mass ratio, as described in Eq. (1). In this example we use a separation a = 10 mas, which means that we would observe a maximum photocentre wobbling of 0.9 mas for a mass ratio q ≈ 0.1538.
In the text

Fig. 3 Plots of the density of observations (number in ten consecutive days) for asteroid (3457) Arnenordheim and selected sample of data. Left: resulting convolution of the observation time span in days vs the number of observations in each ten-day window. Right: sample of data at the peak of the 10-day window of search combining 17 observations. The x-axis is the observation date and the y-axis is the average AL residual from the orbital fit. Each dot represents one averaged set of transits observed in the FOV 1 in blue and FOV 2 in green.

In the text

Fig. 4 Examples of periodograms obtained for asteroid (106016) 2000 SS293 (left plot), with residuals shown on the right plot. The y-axis shows the GLSP P(f) for each frequency shown on the x-axis. The red line represents the periodogram obtained for the data analysed. The blue line represents the GLSP obtained for a purely sinusoidal signal whose parameters correspond to the largest peak of P(f) (about 3 cycles per day here), with the same time sampling as the data and the same error bars for the weights. Finally, the green line represents the spectral window of the time sampling. The right plot shows the AL residual time series that reproduces the periodogram on the left plot.

In the text

Fig. 5 Examples of uncertainty determination from MC simulations under ${\widehat{\mathscr{H}}}_1$. The top and bottom plots show, respectively, the distribution of periods and amplitudes obtained from the best fit in each MC simulation (in blue). The red lines show the value estimated from the first GLSP run only on Gaia data, and the solid black line represents the median obtained from the distribution of amplitude and period values from the MC under ${\widehat{\mathscr{H}}}_1$. The dotted vertical lines show the 68% distribution interval while the dashed lines represent the 95% confidence interval of the distribution. The final results are given as the first estimate from the Gaia data (red line), and the uncertainty is the 95% interval. For instance, in this example, the final amplitude value is 1.85 ± 0.08.

In the text

	Fig. 6 Flowchart showing the period search and statistical selection processes discussed in Sects. 3 and 4. It represents part of the full detection process, which is coupled with physical characterisation tests that are discussed in the following sections.
In the text

Fig. 7 Distribution of pairs (p, Q) for each of the 6000 sets of synthetic data simulated. The partially transparent dots represent the cases where the detections do not meet the requirements and were disqualified, while the remaining dots represent the detections selected. Of the 6000 simulations, 1303 have p-value>0.1. Hence, their Q factors were not calculated and the corresponding points are not shown in this plot.

In the text

	Fig. 8 Distribution of periods from best fit (top row), quality factors (centre row), and p-values (bottom row) as percentages of the detections selected with p-value < 5% and Q > 0.5. The different coloured bars represent the different classes of S/Nbase values explored. The dotted line at 23 h represents the nominal (true) period of the injected signal.
In the text

	Fig. 9 Distribution of time series selected as detections (coloured bars) and discarded (grey) as a function of the post-fit S/N, for the six different S/Nbase classes explored.
In the text

Fig. 10 Resulting distribution of parameters from the 6000 simulations in the noise-only scenario. Left: distribution of pairs (p, Q) for the time series simulated. 5391 of the cases have p-value > 0.1 (for these cases the Q factors were not calculated and the corresponding points are thus not shown). The grey dots represent the cases where the data was disqualified as a detection (cases for which p > 0.05, Q < 0.5, or both) and the red dots represent the detections. The two grey points in the top left square are detections discarded due to  > 20 mas. Right: post-fit S/N distribution for the pairs in the left plot.

In the text

Fig. 11 Visualisation of the output of the main steps of the selection process in the period search for one candidate example ((3457) Arnenordheim). The top left plot shows the periodogram from the fitting of the 17 observations in the window. The peak represents a period of 45.95h and the best fit is shown in the phased data in the top right plot. The distribution of data points not phased is shown in the right plot in Fig. 3. The bottom left plot shows the empirical distribution of GLSP power V (maximum value of the GLSP under ℋ0) and the vertical dashed line shows the position of the value of υ, the power of the peak from the periodogram with the real data, with p(υ) < 0.001%. Finally, the bottom right plot shows the distribution of best frequencies from the 5000 MC simulations under ${\widehat{\mathscr{H}}}_1$). The distribution tends to accumulate around the originally estimated frequency, showing that this detection is robust (here Q is close to 1.)

In the text

	Fig. 12 Distribution of estimated wobble amplitudes and periods obtained from the period search for all objects. In grey are the time series disqualified due to a p-value > 0.05; in blue are the cases disqualified because Q < 0.5, but with p-value within the limit of 0.05; and in red are the cases selected as binary candidates from the period search for complying with p-value< 0.05 and Q > 0.5.
In the text

	Fig. 13 Distribution of differences between the solutions for the 40 objects with multiple windows of consecutive observations.
In the text

Fig. 14 Examples of the bulk density ρ as a function of the mass ratio q for two asteroids, (7616) Sadako (top left) and (6261) Chione (top right). The red line indicates the nominal value. Shown in grey is the error range on the density. The blue and orange horizontal lines represent physically reasonable limits. qmin1, qmax1, qmin2 and qmax2 mark the corresponding extremes of the mass ratio values. The bottom plot shows the case for (930) Westphalia, eliminated in the density filtering because the minimum estimated density is higher than the chosen range of bulk density.

In the text

	Fig. 15 Comparison between the average DR3 astrometric residuals and estimated average photocentre offset for candidates (487) Venetia, (1200) Imperatrix and (2219) Mannucci. The top plots show the variation of the data over time in the along-scan direction, while the bottom plots show the correlation between the two sets of data.
In the text

Fig. 16 Plot of the radius of the primary as a function of the separation intervals normalised with the Roche limit for fluid satellites, for the sample of 156 asteroids with known taxonomy. The squares represent the mean value of the first mass ratio interval estimated and the circles represent the second interval. The coloured symbols with dark borders are the objects from the taxonomic filtering with $\widehat{S/N}>1$, and the grey symbols are those with $\widehat{S/N}<1$.

In the text

	Fig. 17 Distribution of $\widehat{\text{S/N}}$ vs. p-value (top); estimated wobble period vs wobble amplitude (center); and p-value vs. quality factor (bottom) for the 3038 objects obtained from the statistical selection (in grey), the 156 objects that survived the taxonomic filtering (in yellow), and the 67 best candidates that passed all tests and have $\widehat{S/N}>1$ (in pink).
In the text

	Fig. 18 (1509) Esclangona residuals distribution (right) and the periodogram from the fitting of the 21 data points in the window (left).
In the text

	Fig. 19 Distribution of $\widehat{\text{S/N}}$ for our list of 156 binary candidates selected from the taxonomic filtering. The 67 best candidates are those with $\widehat{S/N}>1$
In the text

	Fig. 20 Comparison between the parameters for the Arecibo system. The triangle with black edges and grey error bars represents the values from literature (Table 5), while the circles represent the first (transparent) and second intervals of size ratio estimated from our method using the same wobble amplitude as in Tanga et al. (2023).
In the text

	Fig. 21 Comparison of the fraction of small objects with D < 20 km in each taxonomic complex between the reference sample of background asteroids and the population of our astrometric binary candidates. A taxonomic type with no column indicates that asteroids of that type are present in the reference population but not in the group of binary candidates.
In the text

	Fig. 22 Fraction of binary candidates with 30< D <50 km in each taxonomic complex, with respect to the reference sample of background asteroids.
In the text

	Fig. 23 Plot of the primary rotation period vs. the estimated primary diameter for the best candidates. The black dots are known binaries in the shown intervals and the magenta circles are the best binary candidates obtained.
In the text

	Fig. 24 Plot of the primary rotation period vs. the estimated mass ratio intervals for the 67 selected candidates. The lighter and darker coulored circles represent, respectively, the first and second intervals of possible mass ratio values for the best binary candidates.
In the text

Fig. 25 Plot of the ratio of the estimated orbital period (or wobble period) of the secondary over the photometric rotation period of the primary vs. the estimated size ratio of the binary system. The black dots represent all of the known binaries and the red dots highlight those with intermediate sizes (20 km < Dprimary < 100 km). The light and dark purple circles represent the average value of the best-estimated size ratio intervals for objects smaller than 20 km and the light and dark blue circles represent the larger objects. The blue line shows the condition where the estimated orbital period of the secondary is almost equal to the rotation period of the primary object, representing a 1:1 period ratio, while the orange line represents the 1:2 period ratio.

In the text

	Fig. A.1 Histogram of distribution of the estimated offset ${\widehat{y}}_0$ for each window of search explored. The grey bars represent the empirical distribution, while the coloured lines represent different distributions fitted to the data and their respective best-fit parameters.
In the text

	Fig. B.1 Nominal signal from the generic binary system with extreme eccentricity values.
In the text

	Fig. B.2 ROC curve for the different eccentricity values considered.
In the text

	Fig. B.3 Density histograms for estimated $\widehat{S/N}$ (top) and estimated wobble amplitude (bottom) for the cases with ${\widehat{P}}_{FA}=5\%$.
In the text