<figure><img src="image_40.png" /><figcaption>**Fig. 40.** Same as Fig. 37, but for the APOGEE sample. We note that a significant number of systems are off scale.</figcaption></figure>

<table><caption>**Table 14.** Same as Table 12, but for the APOGEE gold sample of PW20.</caption><thead><tr><th></th><th>e</th><th>γ</th><th>K</th><th>P</th></tr></thead><tbody><tr><td>ΔT ≥ P</td><td>87.6 (121)</td><td>73.6 (121)</td><td>85.1 (121)</td><td>69.4 (121)</td></tr><tr><td>ΔT < P</td><td>100.0 (5)</td><td>100.0 (5)</td><td>100.0 (5)</td><td>80.0 (5)</td></tr><tr><td>All</td><td>88.1 (126)</td><td>74.6 (126)</td><td>85.7 (126)</td><td>69.8 (126)</td></tr></tbody></table>

been discussed by Petrosky et al. (2021). Most of them are (near-)contact EBs. However, as will be shown below, a non-negligible amount are classical pulsators, especially RR Lyrae stars. The infrared photometric measurements are made through the W1 filter centred at about 3.4 µm. A subset of 49 465 systems can be unambiguously associated to a DR3 source<fnref n="12" />. This sample allows us to test the reliability of our orbital periods for very close binaries, although other orbital parameters are not available.

The breakdown of DR3 spectroscopic solution types is as follows for the SB1s in Petrosky et al. (2021):

- SB1: 704 systems,

- SB1C: 3 systems,

- TrendSB1: 12 systems,

- SB2: 69 systems,

- SB2C: 23 systems,

- Without any DR3 spectroscopic deterministic solutions: 48 654 systems.

We find a close correspondence in most cases with either the main WISE period or its second harmonic (Fig. 41). As discussed by Petrosky et al. (2021), the orbital period they found for the binaries is usually half the true value because of the nearly symmetrical nature of the WISE light curves. The WISE period is recovered to within 1% for 25.3% and 55.4% of the cases for the first and second harmonic, respectively. The overall recovery rate at this level is thus about 81%. It remains quite high for lower tolerances (e.g. ~78% for 0.1%). The deviations lie in the vast majority of cases within five times the mutual uncertainties (Fig. 42).

Most of the WISE sources were assigned a classification in DR3 based on their light variations (Mowlavi et al. 2023; Ripepi et al. 2023; Clementini et al. 2023). As can be seen in Fig. 43, a (nearly) circular solution is found for the majority of the EBs, as may be expected for such close binaries. In contrast, the systems with high eccentricities are mostly classical pulsators (Cepheids or RR Lyrae) mistaken as binaries. As described in

<figure><img src="image_41.png" /><figcaption>**Fig. 41.** Ratio between spectroscopic DR3 and photometric WISE periods. DR3 values that can be identified with the first and second WISE harmonics are indicated as vertical, dashed lines.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>**Fig. 42.** Deviations in $\sigma$ units between the DR3 and WISE periods when the latter (or its harmonic) are recovered to within 1%.</figcaption></figure>

Sect. 10.1.1, those pulsating variables were rejected from the final DR3 catalogue during the post-filtering process.

### 9.4.5. Comparison against ASAS-SN

The physical parameters of 35 464 detached eclipsing binaries (DEBs) observed in the framework of the All-Sky Automated Survey for Supernovae (ASAS-SN; Kochanek et al. 2017)<fnref n="13" /> have been determined by Rowan et al. (2022). A total of 35 307 systems have a Gaia DR3 ID. The data are based on the analysis of g- and V-band light curves. The orbital periods nicely bridge the gap between the WISE sample of close binaries and other more general catalogues, such as APOGEE or the SB9 (Fig. 36). Contrary to WISE, other orbital elements (e.g. eccentricity) are provided besides the period, although without any uncertainties. For P, we assume a typical relative uncertainty of $10^{-6}$ (e.g. Holanda et al. 2018).

The breakdown of DR3 spectroscopic solution types is as follows for the binaries in Rowan et al. (2022):

- SB1: 626 systems,

- SB1C: 3 systems,

- TrendSB1: 6 systems,

- SB2: 92 systems,

- SB2C: 35 systems,

- Without any DR3 spectroscopic deterministic solutions: 34 545 systems.

Most objects have $G \ge 13$ mag and are fainter than our magnitude cut-off.

The majority of the 626 objects are classified in DR3 as EBs based on their photometric behaviour (Mowlavi et al. 2023). Conversely, none are identified as classical pulsators

<footnote n="12">¹² http://cdn.gea.esac.esa.int/Gaia/gedr3/cross_match/</footnote>

<footnote n="13">¹³ https://www.astronomy.ohio-state.edu/asassn/</footnote>
<figure><img src="image_43.png" /></figure>

**Fig. 43.** Distribution of eccentricities for the WISE sources for which the period (filled symbols) or its harmonic (open symbols) is recovered to within 1%. Lower panel: eccentricities for the systems classified as EBs (orange; Mowlavi et al. 2023), Cepheids or RR Lyrae (purple; Ripepi et al. 2023; Clementini et al. 2023). Those without a classification are shown in black. The systems are displayed in no particular order along the y-axis. Upper panel: histogram of the distributions for the EBs, pulsating variables, and all SB1 solutions (the last one is divided by 250).

<figure><img src="image_44.png" /></figure>

**Fig. 44.** Same as Fig. 37, but for the ASAS-SN sample. We note that a significant number of systems are off scale.

(Clementini et al. 2023; Ripepi et al. 2023). For this sample with SB1-type solutions, the overall recovery rate of the photometric period to within $3\sigma$ amounts to ~83% (Fig. 44). This estimate includes one system where the ASAS-SN period is actually the harmonic (as discussed for WISE above) or for nine objects twice the true value because one of the two eclipses is too shallow to be detected (see figure 6 of Rowan et al. 2022 for an example). Visual inspection of the ASAS-SN light curves for the nine DEBs supports the latter explanation.

<figure><img src="image_45.png" /></figure>

**Fig. 45.** e-log *P* diagrams for the 626 systems in common between DR3 and ASAS-SN. Each panel corresponds to a different threshold on the solution significance. The dotted and dashed lines show the upper envelopes defined by Mazeh (2008) and Halbwachs et al. (2005), respectively. A cut-off period of one day was adopted for the latter. The position of the curve is only illustrative as this choice is somewhat arbitrary.

The relevance of filtering out short-period systems with a low solution significance was discussed in Sect. 9.2. We show in Fig. 45 the DR3 and ASAS-SN data in an e-log *P* diagram for various cuts in solution significance (see also Fig. 30). The upper envelopes defined by Mazeh (2008) and Halbwachs et al. (2005), modified to $P_{\text{cut-off}}$ at 1 d) are overplotted. The former is based on an empirical fit of the SB9 data, but the latter rests on theoretical arguments. While *P* is recovered in ~67% of the cases for $K/\sigma_K$ values up to 15, the proportion sharply increases for higher significances (success rate of ~90% on average). The differences are generally small in absolute terms, which explains the close match in the period distributions shown in Fig. 45. However, an excess of seemingly eccentric systems with respect to Rowan et al. (2022) is readily visible. The eccentricity distributions become more and more similar when the significance cut-off increases until they become statistically indistinguishable above about 50. Rather than a poor determination of the period for close binaries, the issue thus seems to rather lie in solutions that are too eccentric when they are ill-constrained (see also Jayasinghe et al. 2023).
### 9.4.6. Summary

As a preamble to summarising the main conclusions of our analysis, it is worth emphasising that some of the validation data sets discussed above are of limited usefulness in the sense that they do not offer stringent tests of our results. In particular, the APOGEE orbits from PW20 are based on fewer epochs (about 14 on average concerning the 126 objects) than most NSS sources (see Fig. 14). As a consequence, their orbital solutions are not necessarily more robust and precise. We indeed find a poorer match for the periods with respect to the other validation catalogues, even at high significance. It might be no coincidence that, except for $γ$, the highest recovery rates for systems monitored over a sufficiently long baseline are found for the high-quality Griffin solutions (Table 12).

As discussed above, the success rate for the main orbital parameters (*P*, *K*, *e*, and $γ$) is usually well above 80% for the whole set of validation catalogues considered here (Griffin, SB9, APOGEE, WISE, and ASAS-SN). As expected, the performance increases for relatively short-period systems for which at least one orbital period is covered by the observations (about 2.5 years for the RVS spectroscopic processing data set considered here). In this case, the recovery rate commonly exceeds 85–90%. The WISE sample of (near-)contact binaries allows us to investigate the performance of our period search at high frequencies. We recover the WISE photometric periods at a level of about 80–85%. It increases to about 90% for the DEBs with wider orbits discussed by Rowan et al. (2022). We can also note the excellent success rate for the pulsation periods of the RR Lyrae stars that are below one day (Sect. 9.3). The good agreement on the period with both EB data sets advocates for a deeper study using these data. Their combination will certainly result in an extensive analysis that could reveal various kinds of objects such as EW binaries and/or ellipsoidal variables. Such a task is beyond the scope of the present paper. Our analysis globally suggests that the orbital periods below about 10 days are generally robust provided the solution significance is not too low. We find a roughly constant success rate at the level of about 90–95% once the significance exceeds 15–20. However, although they can often be regarded as compatible with zero within the uncertainties, it appears that the eccentricities may be overestimated even if the significance is larger.

We ascribe the poor correspondence between the DR3 and Griffin/SB9 $γ$'s to a slight zero-point offset at the ~1 km s⁻¹ level. There is indeed no evidence for a significant bias with respect to APOGEE (Fig. F.10). For the other parameters, no obvious systematic discrepancies are observed.

In order to obtain an order-of-magnitude estimate of the period whenever possible, a choice is made to systematically favour an SB1 solution over a TrendSB1 one when the period is estimated to be in the range 0.5 to 1.5 $ΔT$. As a result, an SB1 solution could have been assigned despite the fact that the length of the observations is insufficient. However, although the proportion strongly depends on the validation data set, it concerns less than 15% of the cases. Conversely, overall about 79% of the TrendSB1 solutions are associated to binaries with periods significantly longer than the observations.

Our ability to recover the orbital parameters quite significantly hinges on the retrieval of the true period. When only counting the binaries for which this condition is satisfied in the reference catalogues discussed in Sects. 9.4.1 to 9.4.3, we find an increase of ~5–10% in the recovery rate of *e*, $γ$, and *K*. For *e* and *K*, it typically reaches ~95%. It is evident that other limitations, either related to the characteristics of the observations (e.g.

number of RVS transits) or the very nature of the binary system (e.g. velocity semi-amplitude) can also potentially play a significant role in the quality of the solutions. Some clues about the sensitivity (or lack thereof) to the combination of some parameters can be gauged from plots such as Fig. 46. Besides suggesting a lack of dependency between the recovery rate of the period and some quantities (e.g. *ruwe*), it illustrates, as discussed above, the difficulty in successfully finding the period at low significance values. However, there may be more complex dependencies lurking in the data that cannot be easily and reliably quantified given the quite limited number of cross-matches with the validation catalogues (a few hundreds systems at most).

In this respect, the fact that only a relatively small fraction of the SB1s in the validation catalogues have a spectroscopic deterministic DR3 solution is the consequence of the various filtering steps, both at CU4 and CU6 levels, described previously. In addition, we note that for the SB9 and Griffin's sample, only a subset of the reference catalogue was used. The case of the SB9, and the steps that led to the relatively small overlap with respect to DR3 (see Sect. 9.4.2) are discussed in detail in Appendix G. Similar arguments apply to the other external catalogues and explain the limited number of benchmark binaries with a spectroscopic deterministic solution in DR3, although the proportion of systems lost at every filtering step may markedly vary. As an illustration, only the very brightest ASAS-SN eclipsing binaries get a DR3 solution. While the $T_{eff}$ range is not an issue, we find that ~90% have a (combined) *G* magnitude fainter than 12 mag. Similar figures apply to the WISE sample. Regardless of these intrinsic limitations, it is naturally expected that the completeness level will substantially improve in DR4 thanks to the larger amount of RVs collected and improvements in data analysis.

We finally note that a few systems securely reported as double-lined in the literature have received an SB1-like solution (either eccentric or trend). The causes for this are multiple. For instance, the spectral signature of the secondary could be too weak to be detected by *Gaia* or the lines in the composite spectrum cannot be fully deblended at the moderate resolution of the RVS, as discussed in Sect. 8.2.3. The various cases which can be encountered are discussed in Damerdji et al. (in prep.).

## 10. Additional considerations

### 10.1. Post-filtering

The results of the sole NSS-SB1 processing (see Sect. 8) are providing a plethora of orbital solutions covering a large domain of the values for the orbital parameters. The pipeline nature of the processing necessarily implies that part of the solutions could be wrong or spurious (whatever could go wrong will go wrong at least once). Therefore, we scrutinised a subset of the results in order to detect surviving problems. Since the DPAC policy with respect to the final production of the DR3 catalogue is intended to favour purity of the list rather than completeness, it is necessary to blacklist solutions that are most certainly spurious. We call this procedure post-filtering and it intervenes between the internally filtered output of the pipeline and the publication of the catalogue. This post-filtering can be based on information coming from outside our internal validation. A clear flaw of this procedure is that this analysis will never identify all the pitfalls. Conversely, some objects could be blacklisted whereas the solution was correct. In the following, we address the identified problematic cases in the context of the classes associated to the sole spectroscopic processing.
<figure><img src="image_46.png" /><figcaption>**Fig. 46.** Corner plot for the ASAS-SN sample, with the strong outliers in *P* deviating by more than 5σ highlighted in red. They constitute about 11% of the whole sample.</figcaption></figure>

### 10.1.1. Intrinsic variables

Since the NSS-SB1 processing is applied to time series of variable RVs, it is expected to also treat some intrinsic variables that are exhibiting well-defined RV variations. This could be the case for example for Cepheids (actually the bright ones). The only firewall preventing the processing from delivering spurious orbital solutions for intrinsic variables is the fact that the processing only considers an harmonic content of the RV curve that is Keplerian (or nearly so). Unfortunately, Cepheids RV curves, for example, are renowned for their harmonic content being quite similar to the Keplerian one. This has been pointed out by Imbert et al. (1985); Imbert (1987); Imbert et al. (1989). The problem generated by intrinsic variables has been pointed out in Sect. 9.4.4 and in Sect. 9.4.5 where some of them were detected in the binary catalogue, mainly at short periods and high eccentricities.

As part of DR3, Ripepi et al. (2023) published a complete analysis of the *Gaia* photometric curves for a little more than 15 000 Cepheids (of all classes), 799 having time series of RVs. This list was cross-correlated with our SB1/SB1C list. They found 338 objects in common. Among them, the Cepheid nature is confirmed from the literature for 276 objects. After an eye inspection of photometric and RV curves, five of them were

found doubtful Cepheids and have been rejected. Finally, 333 spurious SB1 solutions due to the Cepheid nature of the objects have been blacklisted from the SB subcatalogue.

Along these lines, we also considered the case of RR Lyrae stars (of all classes too). Clementini et al. (2023) gathered a list of a little more than 270 000 such objects, among which some 1096 have time series of RVs. The cross-match with our SB1 list yields 152 objects in common. Among these, 5 turned out to have a discrepancy between the period from our SB1 pipeline and the one derived by Clementini et al. (2023). Their RR Lyrae classification were thus considered as insecure. Therefore, we rejected (blacklisted), from the SB1 list, 147 objects considered as bona fide RR Lyrae. In summary, we rejected from our SB1 list some 480 sources as being intrinsic variables. Although a few pulsating stars (Cepheids, RR Lyrae) are rejected from our SB subcatalogue, they are nevertheless studied as pulsators elsewhere (see e.g. Eyer et al. 2023).

Another class of intrinsic variable objects exhibiting large photometric amplitudes are the long-period variables (LPVs). No similar secure list exists in DR3 for the LPVs. A tentative preliminary list exists that is based on the analysis of the photometry: the second *Gaia* catalogue of long-period variable candidates (Lebzelter et al. 2023). However, the related results
<figure><img src="image_47.png" /><figcaption>**Fig. 47.** Variability of *Gaia* DR3 5829160851462523008. Upper panel: *G* lightcurve of this object plotted as a function of the phase. Lower panel: RV curve of the same object. The red dotted line represents the systemic velocity. The zero of phase is arbitrary and has been chosen to improve the readability of the graph. Error bars are not given since they have sizes of the order of the symbols.</figcaption></figure>

were not entirely accessible at the time of the execution of our pipeline. In addition, from the later provisional analysis, it became clear that this preliminary list most probably also contains a lot of other objects and in particular ellipsoidal variables. We illustrate here our concerns by further studying the case of *Gaia* DR3 5829160851462523008. This object was not known to vary before *Gaia*. It appears in the catalogue of LPV candidates with a derived photometric period of 90.5 d. This object also appears in the SB subcatalogue, produced in the present work, as a binary with a period of 178.886 d. Figure 47 exhibits both the RV curve and the photometric light curve of this object with the phase computed according to the RV period. From this figure, it is clear that the RV curve is unimodal, as it should be, whereas the *G* light curve presents two cycles in good agreement with the period attributed by Lebzelter et al. (2023). If the two maxima are very similar, the two minima are significantly different. In addition, the minima happen at phases where the RV curve crosses the systemic velocity. All this strongly suggests that this star is an ellipsoidal variable. Therefore, at the time of DR3, it appeared hazardous to filter our solutions on the basis of the existing LPV catalogue. Objects like *Gaia* DR3 5829160851462523008 would have been unduly rejected and we cautiously decided to not reject LPVs on this basis. Certainly, a sound analysis necessitates the concomitant study of RV and photometry, which was not foreseen for DR3. This should be done in the subsequent FPR (see Gaia Collaboration 2023c) where some 20% of the LPV candidates are expected to be ellipsoidal variables. We encourage the user of the SB subcatalogue to filter out true intrinsic variables thanks to these subsequent results.

Yet another source of contaminants are main-sequence pulsators. The main periods of δ Scuti stars are typically below 0.2 d (Gaia Collaboration 2023b) and therefore lie outside the frequency range explored by our pipeline (even considering aliasing). However, γ Doradus stars are a more serious concern because they pulsate with periods in the range 0.3–1.4 d (Gaia Collaboration 2023b) and exhibit RV variations of the order of a few km s⁻¹ (e.g. Mathias et al. 2004).

The problem of the rejection of intrinsic variables is certainly a difficult one, and much progress is still to be made on this topic. As mentioned above, the measurement of RVs by cross-correlating the spectra with template prevents us from identifying the origin of the RV variations.

### 10.1.2. Double-line transits

Since the RVs are originally measured by the STA processing on the basis of a per transit approach, there is initially no constraint from transit to transit. In particular, the same object could exhibit a single-line spectrum at a given time and a composite spectrum at another moment. As part of the spectroscopic processing, a value representing the percentage of transits leading to double-line detection is computed. If the star is constant or assumed non-variable, the median velocity is computed, a cut-off is then applied on this percentage in order to avoid the propagation of the effect of suspicious secondary RVs. If no transit exhibits a composite spectrum, the object presenting variations is processed as a single star. On the opposite, if a sufficient number of transits show a double-line spectrum, the object is processed as an SB2 (see Damerdji et al., in prep.). If the assumed constant object presents only a minority of double-line transits, the latter are ignored and the object is analysed on the basis of the sole single-line data. The presence of 10% or more double-line transits is considered as a source of difficulties for the processing of single constant objects. It has been proposed a posteriori that the same rule (as the one applied to constant stars) should apply to the NSS-SB1 chain. This led to the rejection of 1475 SB1 objects and of one SB1C object.

### 10.1.3. Scan-angle dependent signals

The RVS being a slitless spectrometer, the knowledge of the exact position of a star is necessary to have a correct calibration of the wavelength for its observed spectrum (Sartoretti et al. 2018). The astrometric solution assumes that all the objects are point-like, and the objects presenting incompatibility with this hypothesis have potentially wrong RVs. They are thus associated to large values of the *ruwe* parameter. This point is also mentioned in Sect. 2. To these deviating objects is associated an excess noise compared to the global astrometric solution (labelled *astrometric_excess_noise*). One of the main categories of the non-point-like objects are close pairs (physical or not) unresolved or partially resolved. A blend of objects that presents some elongated image will induce the following effect. The measured position of the star will have a position dependent on how the elongated axis of the pair enters the astrometric CCD field. Thus, the resulting RVs are dependent on the orientation of the main axis of the elongated objects with respect to the scan. It depends on the scan angle in a rather smooth manner. This potentially could generate, from transit to transit, RV variations that are a simple function of the scan angle. In particular, solutions for RV-constant objects could be associated to solutions mimicking a spectroscopic orbit at a period of 62.97 d; these orbits are consequently spurious. During the precession cycle, the scan angles of the various transits are performing a complete round. The phenomenon manifests itself as an excess of solutions with small *K* and thus small mass function (because the period is fixed). It is visible in Fig. 31 and is further illustrated in Fig. 48. Clearly, most of the solutions are spurious and are associated to large values of the astrometric excess noise. The suspicious *K* span the approximate interval between 0.8 and 16 km s⁻¹. Another illustration can be found in Gaia Collaboration (2023a),
<figure><img src="image_48.png" /><figcaption>Fig. 48. Semi-amplitude K as a function of the period for the SB1 solutions. The data points are coloured as a function of the astrometric_excess_noise. A population of solutions very close to the 62.97 d (within ± 3 % as underlined by the vertical red lines) precession period exhibits small K and large such noise. The polluted domain spans values of K between 0.8 and 16 km s⁻¹.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>Fig. 49. astrometric_excess_noise as a function of the semi-amplitude, K, for the selected subsample. The data points are coloured according to the ruwe. The blue line indicates the separation above which the K is lower than three times the excess noise.</figcaption></figure>

where it is also shown that the phenomenon could be related to a particular behaviour of ipd_gof_harmonic_amplitude. A detailed investigation of this behaviour is to be found in Holl et al. (2023a).

Thus, the epoch RVs are affected by unmodelled astrometric shifts along the scan direction. In order to alleviate the problem, the objects in the class SB1 were first selected on the basis of an astrometric ruwe larger than 1.4 and a solution leading to a period within 3% of the precession period (62.97 d). In addition, we selected objects obeying the condition:

$$ \begin{gathered} \text{astrometric\_excess\_noise}(\text{km s}^{-1}) \\ = 0.146 \times \text{astrometric\_excess\_noise}(\text{mas}) > 0.5 \times \epsilon, \end{gathered} \qquad (20) $$

where ϵ is the mean of the uncertainties on the individual epoch RVs over the full time series. Figure 49 illustrates the distribution of the astrometric_excess_noise as a function of the parameter K for the selected sample. We can clearly see that all the K values are lower than three times

the astrometric_excess_noise suggesting that these orbital solutions are spurious and are due to a scan-angle effect. The semi-amplitude is too small relatively to the error induced by the astrometric solution. It has been decided to reject and blacklist the 87 objects presenting these characteristics.

A similar approach concerning the satellite revolution period around the Sun (i.e. within 10% of one year) led to the rejection of another 77 sources.

### 10.1.4. Further remarks

All in all, 2107 objects were blacklisted for the above-mentioned three reasons. The mere sum of the individual rejected objects amounts to 2119, but 12 objects appear in two lists. The post-filtering described here above does not concern classes that have been combined. It could happen that a successful combination leads to a removal of some objects from the purely spectroscopic classes. If later the global fit is not successful, the object will disappear from the catalogue and will not be re-injected in any purely spectroscopic classes. Thus, some countings could be difficult to perform.

Although the post-filtering has been applied with great caution, it also happens that some problems persist. One of these cases is the following. Some objects exhibit a good SB1 or SB1C solution with a well established significant orbital semi-amplitude for Gaia. However, these solutions could sometimes be doubtful because the same object observed in the framework of the APOGEE DR17 survey (Abdurro'uf et al. 2022) or LAMOST survey (Cui et al. 2012; Bashi et al. 2022) appears to be constant in RV as detailed in Sect. 8.2.4. This anomaly has not been fully explained but it could be related in some cases to the presence of a nearby neighbour (see below). These objects have not been rejected. We encourage the catalogue user to be very cautious when using the present SB subcatalogue despite the fact that persisting problems should only represent a minority. All the remaining and identified problems should be corrected in the future DR4.

## 10.2. An explanation for gaia DR3 3376949338201650112

In Sect. 8.2.4, we report the existence of SB1 solutions presenting a marked periodic variation, whereas ground-based data conclude to the constancy of the related object. Here below, we propose an explanation for this discrepancy by investigating in details the case of Gaia DR3 3376949338201650112. This example has the advantage of being enlightening and will allow us to correct the pipeline for future releases.

The astrometric solution for this object is characterised by a ruwe of 25.35, pinpointing to a problem with the fit. Referring to the IPD parameters (image parameter determination; see e.g. Holl et al. 2023a), we remark that ipd_gof_harmonic_amplitude is equal to 0.164 which is certainly different from zero. This is a clue that the object is badly treated as a single object. The value of ipd Frac_multi_peak turns out to be 84 suggesting that the object is expected to be double or multiple on the plane of the sky: most probably a resolved or, more precisely, a partially resolved pair (according to the definition of Holl et al. 2023a). As explained by these authors, such pairs (or higher order) have a predisposition to suffer from scan-angle effects that induce a bias on the measured RVs that is function of the relative orientation of the scanning direction on the sky at the time of the observations. This is confirmed for Gaia DR3 3376949338201650112 by Fig. 50 where it becomes clear that the measured RVs are function of the
<figure><img src="image_50.png" /><figcaption>**Fig. 50.** RVs of *Gaia* DR3 3376949338201650112 as a function of the scan angle at the time of the observations. The blue curve represents the function given by Eq. (2) in Holl et al. (2023a) with the derived *K* taken from the catalogue.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**Fig. 51.** Amplitude spectral window as a function of frequency as computed from the acquired time series for *Gaia* DR3 3376949338201650112. The insert is a zoom on the important second-highest peak, the first being at the origin as usual.</figcaption></figure>

scan angle. The parameter ipd_gof_harmonic_phase adopts a value of 103.4° that suggests that the line joining the two stars of the pair has a position angle on the sky of 13.4°, after the application of a 90° correction as explained in Sect. 3.1.2 of Holl et al. (2023a). The semi-amplitude *K* of the fake orbital solution is about 16.6 km s⁻¹ which translates according to Eq. (2) of Holl et al. (2023a) into a separation of the pair of 114 mas. The existence of such a pair is confirmed by a preliminary run of AGIS (AGIS 4.1) linked to the future DR4 that actually detects two objects *Gaia* DR4 3376949342501945344 and *Gaia* DR4 3376949342501945472 that are separated by some 200 mas and are exhibiting an orientation of 9°, thereby supporting our interpretation.

In principle, the bias due to the scan-angle effect is intended to produce variability with a period either of 62.97 d or of 1 y. This property has been used in Sect. 10.1.3 to apply post-filtering. However, in the present case, the period is 0.657582 d which does not appear to be compatible and could shed some doubt on our interpretation. The spectral window associated to the DR3 RV time series corresponding to this object is given in Fig. 51.

Besides the main peak at $v = 0$ d⁻¹ as usual, it exhibits an important peak at $v = 1.51784$ d⁻¹ that is exceeding a value of 0.8; this is an indication that isolated aliasing could be important for this particular time series. An aliasing with $\delta v = 1.51784$ d⁻¹ explains the ambiguity between a period of 0.657582 d and that of 1 year. The combination of the aliasing phenomenon and of the scan-angle effect explains the existence of such fake SB1 objects in our SB subcatalogue: such objects were not post-filtered, since the post-filtering was based on the 62.97 d or 1 y timescales. We conclude that it would be more appropriate to include a statistical test inspired by graphs like the one shown in Fig. 50 rather than on a criterion based on the period. However, a true periodicity with these timescales will imply a similar behaviour of the RV dependency on the scan angle even if no scan-angle effect is present.

The explanation given here does not apply to the other object cited in Sect. 8.2.4 (Gaia DR3 702159502868472832).

## 10.3. The case of astrometric binaries

In Sects. 10.1.3 and 10.2, we discussed the effect of astrometric close pairs on the determination of RVs. A close pair fixed on the sky could often generate a scan-angle effect and lead to corrupted RV values but also could generate fake orbital motions. An additional degree of complexity is reached if we consider astrometric binaries. Due to their binarity, these close pairs are not fixed on the sky and the main axis of the pair is rotating due to the orbital motion. Therefore, there is interaction between the above-mentioned scan-angle effect and the astrometric orbital motion itself. The situation still increases in complexity if the astrometric orbital motion and the spectroscopic RV changes correspond to two aspects of the very same orbital motion and thus have the same period. In this latter case, the correction presents a timescale equal to the phenomenon we are studying. The correcting scheme is explained in Sect. 3.3.1 of Holl et al. (2023a) and is also illustrated there with the (somewhat extreme) example of *Gaia* DR3 6631710607341412096. For details, we refer to this paper.

The effect addressed here was not anticipated and was not fully understood at the time of the final delivery of the NSS code for DR3 processing. As a consequence, some of our solutions are damaged by this phenomenon. In addition, those pertaining to the AstroSpectroSB1 class should be the most affected. For DR4, a preliminary run of the NSS astrometric channel has been included in order to compute a correction before the CU6 processing performing the RV measurements.

# 11. The SB subcatalogue

## 11.1. Description of the output

The present first catalogue of pure spectroscopic solutions is part of the *Gaia* Data Release 3. The orbital solutions (in class SB1 or SB1C as indicated in the entry nssSolution_id) are part of an auxiliary table named gaiadr3.nss_two_body_orbit. We output in this table the classical Campbell orbital parameters introduced in Sect. 3.1 for the SB1 class ($P, \gamma, K, e, \omega, T_0$) along with their estimated uncertainties ($1\sigma$). The same is done for circular orbit parameters ($P, \gamma, K, T_0$) introduced in Sect. 3.2 and corresponding to the SB1C class. Beyond the orbital parameters themselves, this table also contains some indicators of the quality of the solution (see Sect. 6.4). The flags associated to the solution (see Sect. 6.5) are also provided. The boolean mask bit_index defines the parameters that are actually fitted, and
corresponds for the SB1 class (respectively SB1C) to the value 127 (respectively 31). In addition, is also provided the upper triangle of the correlation matrix of the parameters. It is expressed as `corr_vec` that contains the cross-coefficients folded into a vector and corresponding to the ordering of the adopted list of parameters. The first term of the vector C(0) gives the element (0,1) of the matrix, the element C(1) gives the element (0,2) of the matrix and element C(2) the (1,2) one. Elements C(3), C(4), C(5) correspond respectively to the matrix elements (0,3), (1,3), (2,3), and so on. This is thus a column-major storage scheme. The details of the output are described in Table H.1 for the class SB1 and in Table H.2 for the SB1C one. The table contains 181327 objects (with solutions) in the SB1 class and 202 in the SB1C one. The marked difference in the counting is due to our decision to consider eccentric solutions down to as small eccentricities as possible before going circular. This was described in Sect. 6.3.

We also released the trend solutions as discussed in Sect. 3.3 and in Sect. 6.3.5. The objects belonging to the TrendSB1 class are part of the auxiliary table `gaiadr3.nss_non_linear_spectro`. The subtype of trend solutions (linear or second-order trends) are indicated in `nssSolution_id` (respectively `FirstDegreeTrendSB1` and `SecondDegreeTrendSB1`). They are characterised by a `bit_index` of 7 and 15, respectively. A few statistics about the quality of the fit are also given. The table contains 24 083 objects with solutions of type `FirstDegreeTrendSB1` and 32 725 objects with `SecondDegreeTrendSB1` ones, thus totalling 56 808 objects with solutions in the class `TrendSB1`. The details of the content of the latter table are described in Table H.3. We recall that the StochasticSB1 solutions are not delivered in DR3.

We should also point out that the auxiliary file `gaiadr3.nss_two_body_orbit` also contains the orbital parameters of the classes `AstroSpectroSB1` and `EclipsingSpectro` that result from the action of the combiner (see Sects. 3.5 and 7.2).

## 11.2. Weaknesses, caveats, and recommendations

The present section gathers the caveats dispersed all over the paper and delivers a summary of the characteristics of the catalogue that could be beneficial to the general user. The RVs studied here have the quality and shortcomings of the measurements of individual (epoch) RVs for constant stars. The properties and precision of these RVs for SB1 studies should be of the same order. Therefore, we refer the reader to the work of Katz et al. (2023) concerning the characteristics and condition of utilisation of these data. Seabroke et al. (in preparation) is also a work to be consulted.

The information extracted from variations of the RVs could be delicate to use. Firstly, it must be clear that the way the measurements are performed (cross-correlation with a template) makes the implicit assumption that the line is supposed to vary without changing its shape. This is not necessarily true. Two textbook cases at least come to mind. Some stars, notably the intrinsically variable ones, could exhibit dominant variations of their line shapes. In such a case, the variations do not represent a Doppler shift of the object as a whole and thus are not representative of an SB1 (see Sect. 10.1.1). In principle, the fact we are fitting models where the harmonic content is typical of a Keplerian motion is certainly a protection against variable stars appearing in the SB1 class. However, some variable stars, like Cepheids are well-known to display RV changes that

closely mimic the Keplerian motion (e.g. Imbert et al. 1989). As explained in Sect. 10.1.1, efforts have been made to post-filter the Cepheid and the RR Lyrae stars, but this post-filtering is certainly not complete for the time being. In addition, no such step was made for other intrinsic variables in DR3 (notably Long-Period Variables).

Another tricky case is that of true SB2 stars, where the two lines never fully deblend. They are appearing most probably in the SB1 channel. The RV variations are then due to the corresponding change of the full profile made by the blended lines. This has little effect on the period determination but the impact on the eccentricity and mainly on the semi-amplitude might be huge (Sect. 8.2.3). The related solutions are characterised by a too small value of *K* and of the mass function. They could be considered as suspicious but we cannot have any certitude that all of these types of solutions are necessarily wrong.

The algorithm used is intended to fit orbital solutions to the RV curves that it analyses. It is doing so but is not protected against some anomalies. The RV data set suffers from all the problems not attenuated by the CU6 processing and is not particularly protected against outliers, miss-classified double-line objects, intrinsic variables, emission-line objects and other problematic situations (e.g. badly normalised continuum, etc...). In addition, it is not (for the moment) able to properly deal with higher degree systems. It is possible that the code detects a binary motion, but this one could be quantitatively corrupted by the presence of additional components (see Sect. 8.2.2). In the future (DR4), we hope to include a cleaning (whitening) of the dominant signal followed by a new search for a second period or a trend. However, the presence of two signals renders the statistical analysis much more complicated.

Instrumental effects are also present. The scan-angle effect notably induces periodic RV variations for constant stars that do not fulfil the assumptions used in astrometry (objects are assumed to be isolated and thus single). Close astrometric pairs could thus mimic SB1 RV variations (see the work by Holl et al. 2023a), and are polluting the present sample (even if some post-filtering took place; see also Sect. 10.1.3). If, in addition, the close pair is rotating on the plane of the sky as is the case for astrometric binaries, the measured RV could be corrupted (see Sect. 10.3). This problem will be corrected for DR4.

In addition, it should be noted that all decisions are taken in a statistical way; this means that the output is statistically correct but that some incorrect individual fits can survive near the borders where the decisions are taken (e.g. between an orbital fit and a trend solution for example).

We recommend to the user of the catalogue to keep in mind the weaknesses and caveats described here above. We also strongly recommend, when possible, the use of the various quality factors accompanying each solution. If the catalogue is expected to be well-behaved from the statistical point of view, the situation could be more intricate for individual objects or particular set of objects. It remains that biases in the catalogue could be present from the point of view of stellar population studies. Such an analysis is beyond the scope of the present paper. For the time being, the genuine selection function of the catalogue is poorly known and any user should be aware of this and behave accordingly. The user is strongly encouraged to adapt the selection of solutions to their particular objective instead of using it blindly. They should combine the various delivered quality factors according to their aim. We advise a cautious approach. An example of interesting complex combinations is detailed in Bashi et al. (2022). To some extent, other information can be considered as quality factors (such as flags but, mainly, the number of
points in the time series $N_{good}$), which could also be included in the selection process.

The global weak point of the algorithms could be the determination of the period even if here above it is shown that the success rate is at the very least 80% and is more often located in the range 85–95%. There is evidently at a very low number of transits, a difficulty to correctly determine the period that forced us to neglect data with $N_{good}$ less than 10. However, at slightly larger $N_{good}$, the situation is not so simple. From the simulations described in Sect. 9.1, we conclude that a value of $N_{good}$ smaller than 17 is somewhat problematic (although it depends on the signal-to-noise ratio) for the period recovery. We also show, from the validation (see Sect. 9.2), that an overfitting phenomenon characterises part of the solutions associated to an $N_{good}$ in the range 10–15. We should further add that, at small $N_{good}$, the periodogram has a larger propensity to produce isolated aliasing peaks (generating ambiguity due to aliasing). While interpreting the simulations, it appears that it is unclear if the determining parameter is the $N_{good}$ or the significance. According to these simulations, significance lower than 10–15 can be associated with difficulties to have the proper recovery rate. This seems to be confirmed by the validation in Sect. 9.4, in particular in the corner plot of Fig. 46. It is also clear, from Fig. 30, that the circularisation phenomenon is not outstandingly visible for significance much lower than 30. Therefore, we consider that depending on the studies performed, the significance should be cut between 15 and 25, or more conservatively at 30. If 10 is still affordable in some cases, the solutions characterised by a significance between 5 (our cut-off) and 10 are most probably hazardous to accept.

The spectroscopic solutions could be merged with photometric or astrometric ones, and in these cases the resulting output is certainly very secure for the period. If the user of the catalogue only considers the sole results of the spectroscopy fits, it is advisable to restrict the sample to some kind of gold sample. This can be done by performing additional selections in the catalogue: for example, a larger number of transits, a higher significance, periods in a dedicated range (between ten days and $\Delta T$, or allowing it to go down to smaller periods but at the expense of a selection on the basis of a larger significance), periods with small uncertainties, solutions with negligible phase gaps (although the gap values are not provided in the DR3 catalogue), and so on. In Appendix J, we present indicative but simple examples of extraction of a gold sample from the catalogue.

## 12. Conclusions

Within the framework of the Gaia Data Release 3, we present the first version of the Gaia catalogue of spectroscopic orbital solutions (SB subcatalogue). We present 181 327 objects that have been found to be characterised by class SB1 and 202 objects by class SB1C. We also present objects exhibiting trends of variations of their RVs over time.

The present catalogue constitutes the largest collection of objects with a spectroscopic orbital solution by far, based on a comparison with the notorious SB9 compilation<fnref n="14" />. The improvement in terms of number of objects is more than one order (nearly two) of magnitude. The present SB subcatalogue has also the major advantage of being homogeneous from the point of view of both data acquisition (properties) and data processing, surpassing the characteristics of existing compilation. The only exception is the typical number of measurements that constitutes

the individual time series, which remains limited (at least for the present DR3 catalogue). The situation will undoubtedly improve for future data releases.

After a description of the data entering the NSS-SB1 processing, we describe in detail the logic and structure of the related code. We reveal the general philosophy behind, along with the mathematical formalism used. We also describe, following the main processing details, the filtering of the results performed to discard bad solutions and improve the catalogue. We give a number of examples illustrating the quality of the output, as well as a few examples of detected problematic cases. A full validation of the output is performed both on the basis of internal consistency and on the basis of external existing catalogues from the literature. We further analyse the characteristics of the catalogue and also describe the way the SB subcatalogue has been included in the Gaia Data Release 3. We additionally scrutinises the weaknesses of the approach and provide the necessary caveats and recommendations for scientific applications of the catalogue.

## Data availability

The epoch RVs constituting the time series analysed here have not been published in the context of the Gaia Data Release 3. The orbital solutions and trends we derived are available since 13 June 2022, in the ESA Gaia DR3 archive accessible at https://gea.esac.esa.int/archive. The same results can also be found at CDS Strasbourg via the link https://cdsarc.cds.unistra.fr/viz-bin/Cat?I/357. Appendices D to J are only available on Zenodo (see https://zenodo.org/records/13990211).

Acknowledgements. The authors are particularly thankful to Alain Detal who maintained over the years the Liège computer park with great efficiency. EG is also thanking F. Mélen for her consequent help in the manipulation of various computer codes. The authors are acknowledging the numerous scientific discussions that took place with Christos Siopis. The participation in this endeavour has been made possible through an Impulsion Credit from the University of Liège to EG. 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. The Gaia mission and data processing have, in Liège, been continuously supported by the Belgian Politique Scientifique Fédérale (BELSPO) through PRODEX contracts Gaia DPAC: Binaires, étoiles extrêmes et objets du système solaire (convention PEA 4200090289; years 2007–2013), Gaia DPAC Early Mission Belgian consolidation (Gaia-be) (convention PEA 4000110152; years 2014–1016), Gaia DPAC Early Mission Belgian consolidation (Gaia-be2) (convention PEA 4000119840; years 2017–1019; years 2020–2022). TM acknowledges financial further support from BELSPO for contract PLATO mission development. We acknowledge support from an ARC grant (13/17-06; Massive Stars: drivers of the evolution of the Universe) for concerted research actions financed by the French community of Belgium (Wallonia-Brussels Federation). We are also largely indebted to the Fonds de la Recherche Scientifique (F.R.S.-FNRS, Belgium) for various supports. YD acknowledges support to the Algerian Centre de Recherche en Astronomie, Astrophysique et Géophysique of Bouzareah Observatory. GS acknowledges support from BELSPO/PRODEX under various contracts and in particular convention PEA 40000119826. AS acknowledges support from the Italian Space Agency (ASI) under contract 2018-24-HH.0 “The Italian participation to the Gaia Data Processing and Analysis Consortium (DPAC)” in collaboration with the Italian National Institute of Astrophysics.

## References

* Abdurro'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35

* Andrae, R., Fouesneau, M., Sordo, R., et al. 2023, A&A, 674, A27

* Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360

<footnote n="14">¹⁴ https://sb9.astro.ulb.ac.be/</footnote>
Armstrong, D. J., Kirk, J., Lam, K. W. F., et al. 2015, A&A, 579, A19

Bashi, D., Shahaf, S., Mazeh, T., et al. 2022, MNRAS, 517, 3888

Bevington, P. R., & Robinson, D. K. 2003, Data Reduction and Error Analysis for the Physical Sciences, 3rd edn. (MA: McGraw-Hill)

Binnendijk, L. 1960, Properties of Double Stars; A Survey of Parallaxes and Orbits (Philadelphia: University of Pennsylvania Press)

Bloomfield, P. 1976, Fourier Analysis of Time Series: An Introduction (New York: Princeton University, John Wiley and Sons, Inc.)

Boulkaboul, A., Damerdji, Y., Morel, T., et al. 2022, MNRAS, 517, 1849

Brandt, T. D. 2021, ApJS, 254, 42

Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513

Clementini, G., Ripepi, V., Garofalo, A., et al. 2023, A&A, 674, A18

Cropper, M., Katz, D., Sartoretti, P., et al. 2018, A&A, 616, A5

Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Res. Astron. Astrophys., 12, 1197

Cumming, A. 2004, MNRAS, 354, 1165

David, M., Blomme, R., Frémat, Y., et al. 2014, A&A, 562, A97

Deeming, T. J. 1975, Ap&SS, 36, 137

Deeming, T. J. 1976, Ap&SS, 42, 257

Delchambre, L. 2010, Recherche d'une évaluation robuste de l'excentricité dans des systèmes binaires d'étoiles, Master Thesis Liège University (available at https://orbi.uliege.be/handle/2268/322803)

Eichhorn, H. K., & Xu, Y.-L. 1990, ApJ, 358, 575

Endl, M., Kürster, M., Els, S., et al. 2002, A&A, 392, 671

ESA, ed. 1997, ESA Special Publication, 1200, The HIPPARCOS and TYCHO catalogues. Astrometric and photometric star catalogues derived from the ESA HIPPARCOS Space Astrometry Mission

ESA, 1997, VizieR Online Data Catalog: I/239

Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4

Eyer, L., Audard, M., Holl, B., et al. 2023, A&A, 674, A13

Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5

Ford, E. B., & Gregory, P. C. 2007, in Astronomical Society of the Pacific Conference Series, 371, Statistical Challenges in Modern Astronomy IV, eds. G. J. Babu, & E. D. Feigelson, 189

Fraser, O. J., Hawley, S. L., & Cook, K. H. 2008, AJ, 136, 1242

Frémat, Y., Royer, F., Marchal, O., et al. 2023, A&A, 674, A8

Fukushima, T. 1996, AJ, 112, 2858

Gaia Collaboration (Brown, A. G. A., et al.) 2016a, A&A, 595, A2

Gaia Collaboration (Prusti, T., et al.) 2016b, A&A, 595, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2018a, A&A, 616, A1

Gaia Collaboration (Katz, D., et al.) 2018b, A&A, 616, A11

Gaia Collaboration (Arenou, F., et al.) 2023a, A&A, 674, A34

Gaia Collaboration (De Ridder, J., et al.) 2023b, A&A, 674, A36

Gaia Collaboration (Trabucchi, M., et al.) 2023c, A&A, 680, A36

Gaia Collaboration (Vallenari, A., et al.) 2023d, A&A, 674, A1

Golub, G., & Pereyra, V. 2003, Inverse Probl., 19, R1

Gosset, E., Royer, P., Rauw, G., Manfroid, J., & Vreux, J. M. 2001, MNRAS, 327, 435

Griffin, R. F. 1969, MNRAS, 145, 163

Griffin, R. F. 2007, The Observatory, 127, 113

Halbwachs, J. L., Mayor, M., & Udry, S. 2005, A&A, 431, 1129

Halbwachs, J.-L., Pourbaix, D., Arenou, F., et al. 2023, A&A, 674, A9

Heck, A., Manfroid, J., & Mersch, G. 1985, A&AS, 59, 63

Holanda, N., Barros, P. P., & Silva, A. N. 2018, New A, 65, 21

Holl, B., Fabricius, C., Portell, J., et al. 2023a, A&A, 674, A25

Holl, B., Sozzetti, A., Sahlmann, J., et al. 2023b, A&A, 674, A10

Imbert, M. 1987, A&A, 175, 30

Imbert, M., Andersen, J., Ardeberg, A., et al. 1985, A&AS, 61, 259

Imbert, M., Andersen, J., Ardeberg, A., et al. 1989, A&AS, 81, 339

Jayasinghe, T., Rowan, D. M., Thompson, T. A., Kochanek, C. S., & Stanek, K. Z. 2023, MNRAS, 521, 5927

Katz, D., Sartoretti, P., Cropper, M., et al. 2019, A&A, 622, A205

Katz, D., Sartoretti, P., Guerrier, A., et al. 2023, A&A, 674, A5

Kawata, D., Baba, J., Ciucä, I., et al. 2018, MNRAS, 479, L108

Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72

Khoperskov, S., Gerhard, O., Di Matteo, P., et al. 2020, A&A, 634, L8

Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502

Koppelman, H. H., Helmi, A., Massari, D., Roelenga, S., & Bastian, U. 2019, A&A, 625, A5

Kunder, A., Kordopatis, G., Steinmetz, M., et al. 2017, AJ, 153, 75

Lebzelter, T., Mowlavi, N., Lecoeur-Taibi, I., et al. 2023, A&A, 674, A15

Lester, K. V. 2020, PhD thesis, Georgia State University, United States

Lindegren, L. 2019, Gaia Data Processing and Analysis Consortium (DPAC) technical note GAIA-C3-TN-LU-LL-125, http://www.cosmos.esa.int/web/gaia/public-dpac-documents

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2

López-Corredoira, M., & Sylos Labini, F. 2019, A&A, 621, A48

López-Corredoira, M., Garzón, F., Wang, H. F., et al. 2020, A&A, 634, A66

Lucy, L. B., & Sweeney, M. A. 1971, AJ, 76, 544

Marchetti, T., Rossi, E. M., & Brown, A. G. A. 2019, MNRAS, 490, 157

Mathias, P., Le Contel, J. M., Chapellier, E., et al. 2004, A&A, 417, 189

Mazeh, T. 2008, in EAS Publications Series, 29, EAS Publications Series, eds. M. J. Goupil, & J. P. Zahn, 1

Monari, G., Famaey, B., Siebert, A., Wegg, C., & Gerhard, O. 2019, A&A, 626, A41

Mowlavi, N., Holl, B., Lecoeur-Taibi, I., et al. 2023, A&A, 674, A16

Nelson, R. H., & Robb, R. M. 2017, Inform. Bull. Variable Stars, 6203, 1

Nitschai, M. S., Cappellari, M., & Neumayer, N. 2020, MNRAS, 494, 6001

Perryman, M. 2009, Astronomical Applications of Astrometry: Ten Years of Exploitation of the Hipparcos Satellite Data (Cambridge University Press)

Petrosky, E., Hwang, H.-C., Zakamska, N. L., Chandra, V., & Hill, M. J. 2021, MNRAS, 503, 3975

Pourbaix, D., Arenou, F., Gavras, P., et al. 2022, Gaia DR3 Documentation. Chapter 7: Non-single stars, Gaia DR3 Documentation, European Space Agency; Gaia Data Processing and Analysis Consortium. Online at https://gea.esac.esa.int/archive/documentation/GDR3/index.html, 7

Pourbaix, D., Tokovinin, A. A., Batten, A. H., et al. 2004, A&A, 424, 727

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (Cambridge University Press)

Prevot, L., Andersen, J., Ardeberg, A., et al. 1985, A&AS, 62, 23

Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, ApJ, 837, 20

Price-Whelan, A. M., Hogg, D. W., Rix, H.-W., et al. 2020, ApJ, 895, 2

Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, A72

Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2023, A&A, 674, A29

Riello, M., De Angeli, F., Evans, D. W., et al. 2018, A&A, 616, A3

Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3

Ripepi, V., Clementini, G., Molinaro, R., et al. 2023, A&A, 674, A17

Rowan, D. M., Jayasinghe, T., Stanek, K. Z., et al. 2022, MNRAS, 517, 2190

Sartoretti, P., Katz, D., Cropper, M., et al. 2018, A&A, 616, A6

Sartoretti, P., Blomme, R., David, M., & Seabroke, G. 2022, Gaia DR3 Documentation. Chapter 6: Spectroscopy, Gaia DR3 Documentation, European Space Agency; Gaia Data Processing and Analysis Consortium. Online at https://gea.esac.esa.int/archive/documentation/GDR3/index.html, 6

Sartoretti, P., Marchal, O., Babusiaux, C., et al. 2023, A&A, 674, A6

Schlesinger, F. 1910, Publ. Allegheny Observ. Univ. Pittsburgh, 1, 33

Seabroke, G. M., Fabricius, C., Teyssier, D., et al. 2021, A&A, 653, A160

Sharma, S. 2017, ARA&A, 55, 213

Soubiran, C., Jasniewicz, G., Chemin, L., et al. 2018, A&A, 616, A7

Sterne, T. E. 1941, PNAS, 27, 175

Tatum, J. B. 1968, MNRAS, 141, 43

Trick, W. H., Coronado, J., & Rix, H.-W. 2019, MNRAS, 484, 3291

Wilson, E. B., & Hilferty, M. M. 1931, PNAS, 17, 684

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577
## Appendix A: Historical developments

Once the period has been evaluated, the general model to fit (Eq. 1) remains non-linear in some of the parameters (mainly the eccentricity, *e*). This implies that the $\chi^2$ of the fit could present a few or even numerous minima in the parameter space, and to find the global one is a true challenge. In particular, the derivation of the eccentricity is a notorious pitfall. The search for the global minimum might be performed empirically by guessing its approximate position and then by descending the $\chi^2$ valley by classical numerical methods (such as the Levenberg-Marquardt algorithm). However, in a pipeline as the one developed for Gaia, identifying the location of the global minimum is uncertain and a specific type of search is necessary. This is particularly frustrating because a well-trained expert eye is able to estimate a rough value for the eccentricity with great efficiency. In this context, we explored the possibilities to adopt alternative approaches. We acted along two main axes. On the one hand, we tried to invent a method based on pattern recognition in the phase diagram of the RVs folded with the candidate periods. Such an approach would allow us to take a first guess for the eccentricity, which could then be used for further classical derivation. In order to apply pattern recognition techniques, it is mandatory to define an interpolation function in the phase diagram to obtain a robust description of the shape of the RV curve. Several tools have been explored such that the nearest-neighbour interpolation, the linear interpolation, the interpolation by cosines, the local Akima periodic interpolation, the periodic interpolation by cubic splines, and the periodic interpolation by Catmull-Rom splines. We also explored the possibilities of performing classical adjustment by least-square methods. In particular, we considered classical polynomials, Fourier series, Bézier curves (using De Casteljan algorithm) and finally, cubic B-splines in their open uniform version. Once the interpolation curve is chosen and the shape is defined, we tested various pattern recognition methods such as algorithms for regression based on the nearest neighbourhood, algorithms for regression based on trees and in particular on extremely randomised trees and, finally, algorithms for regression based on neural networks. Full details of this exploration and detailed references are given in the work of Delchambre (2010).

On the other hand, we developed an approach based on genetic algorithms. These algorithms have the advantage to circumvent any preliminary guess of the solution. This is at the expense of being obliged to explore (even cleverly) the space of the parameters using the generation of random numbers (Monte Carlo methods). This work is described in details in Delchambre (2010). The use of genetic algorithms is a trade-off between the probability to effectively detect the global $\chi^2$ minimum and the consumption of computer time. The codes using genetic algorithms were written, delivered and tested. The detection of the global minimum was possible but was implying a prohibitively large amount of computer time despite strong efforts made to accelerate it. The large amount of objects to treat in the Gaia processing suggests using much more efficient algorithms. The genetic algorithm code was not activated for the operational runs.

Other methods are now available to derive orbital solutions. They became widespread in the domain of exoplanets, where the orbits are particularly incompletely covered and complex. The majority of these methods are based on Monte Carlo simulations. Very often, they use the Markov chain Monte Carlo (MCMC) method (see e.g. Sharma 2017) often utilised in conjunction with Bayesian approaches. These algorithms necessitate a dense coverage of the full space of priors, at the expense of the computer time. Comparing the works of Price-Whelan et al. (2017) and Price-Whelan et al. (2020) pinpoints the strong impact that could have the choice of priors on the resulting posteriors. The Bayesian methods have admittedly the big advantage to potentially provide marginal posterior probabilities for a series of models and thus estimate how much the observational data could favour a particular model. According to Ford & Gregory (2007), the challenge of estimating marginal posterior probabilities is still a crucial problem and obvious estimators are suffering from poor convergence properties. In addition, for the preliminary DR3 catalogue, the directives were that we should output one and only one model per treated object, thus mitigating the interest of using Bayesian methods. In addition, these Monte Carlo methods are computer time-consuming and were not acceptable for the DR3 processing. Finally, we adopted a method based on fully analytical and classical numerical developments.

## Appendix B: Practical details on the processing

The Data Processing Centre at CNES (DPCC<fnref n="15" />), located at the Centre National d'Etudes Spatiales in Toulouse (France), ran the spectroscopic processing chain (CU6) and the object processing chain (CU4). Globally, it ran the whole set of NSS codes and, in particular, the spectroscopic orbital solution chain described in the present paper. Our NSS-SB1 processing chain delivers the results to the combiner. The Main Data Base (MDB) is fed with solutions from this combiner. If the combiner finds no other solutions to combine with, the spectroscopic results are directly transmitted to the MDB. If a combination is possible, only the output of the combination is saved.

In the framework of the Gaia endeavour, the DPCC is hosting a cluster of 250 computers (6000 cores, 40 TBytes of RAM and 7 PBytes of HDFS disk storage). The resources are managed by a Hadoop system. The codes are integrated into a framework SAGA (System of Accomodation of Gaia Algorithms) developed by Thales. The spectroscopic processing was executed over 2500 cores and occupied 17 TBytes of RAM. Some $3 \times 10^6$ hours of CPU were necessary, equivalent in real time to 120 days. Concerning the NSS chain, the full treatment (Ingestion, Processing, and Post-Processing) was executed on the Hadoop cluster over 1400 cores and necessitated 16 105 days of CPU corresponding to an effective duration of 17 days. All the codes were written in Java 8.0. The version of the code tagged and executed for the operational run was 20.10.

<footnote n="15">¹⁵ https://www.cosmos.esa.int/web/gaia/dpac</footnote>
# Appendix C: Separable variable technique applied to the Keplerian model fit

## Appendix C.1. A reminder of the Levenberg-Marquardt method

Eq. (10) links the N dependent variables RV($t_i$) (the RVs of the model) to the independent ones $t_i$ (the observing times). Hereafter, we represent these two variables by the column arrays RV and t respectively. f represents the Keplerian function, such that

$$ \mathbf{RV} = f(\mathbf{t}, \mathbf{p}), \qquad (C.1) $$

with p being the array of 6 parameters [$\gamma, \alpha, \beta, v_k, e, T_0$]$^T$ where the superscript T denotes the transpose operation. The Newton-Raphson algorithm is based on the Taylor first order development

$$ \mathbf{RV} = f(\mathbf{t}, \mathbf{p}^{[n]}) + \mathbf{J}^{[n]} \Delta \mathbf{p}^{[n]}, \qquad (C.2) $$

where $\Delta \mathbf{p}^{[n]} = \mathbf{p}^{[n+1]} - \mathbf{p}^{[n]}$, and {$n$} denotes the iteration number. J is the Jacobian matrix which elements are

$$ J_{i,j}^{[n]} = \frac{\partial f(t_i, \mathbf{p}_j^{[n]})}{\partial \mathbf{p}_j^{[n]}} \quad \text{for } 1 \le i \le N, 1 \le j \le 6. \qquad (C.3) $$

The objective function to be minimised is the quadratic sum of residuals as

$$ \chi^2 = [\mathbf{RV}_{\text{obs}} - f(\mathbf{t}, \mathbf{p})]^T \mathbf{W} [\mathbf{RV}_{\text{obs}} - f(\mathbf{t}, \mathbf{p})], \qquad (C.4) $$

with W being the inverse of the data variance-covariance matrix. The $\chi^2$ is given at the n-th iteration by

$$ \chi^{2,[n]} = [\mathbf{RV}_{\text{obs}} - f(\mathbf{t}, \mathbf{p}^{[n]}) - \mathbf{J}^{[n]} \Delta \mathbf{p}^{[n]}]^T \mathbf{W} [\mathbf{RV}_{\text{obs}} - f(\mathbf{t}, \mathbf{p}^{[n]}) - \mathbf{J}^{[n]} \Delta \mathbf{p}^{[n]}], \qquad (C.5) $$

where RV$_{obs}$ contains the observed RVs. The optimal parameter step $\Delta \mathbf{p}^{[n]}$ is obtained by solving

$$ \frac{\partial \chi^2}{\partial \Delta \mathbf{p}^{[n]}} = 0. \qquad (C.6) $$

In the case of Newton-Raphson, this leads to the following algebraic system:

$$ (\mathbf{J}^{[n]T} \mathbf{W} \mathbf{J}^{[n]}) \Delta \mathbf{p}^{[n]} = \mathbf{J}^{[n]T} \mathbf{W} [\mathbf{RV}_{\text{obs}} - f(\mathbf{t}, \mathbf{p}^{[n]})]. \qquad (C.7) $$

The matrix $\mathbf{H}^{[n]} = \mathbf{J}^{[n]T} \mathbf{W} \mathbf{J}^{[n]}$ is a positive semi-definite matrix called the hessian matrix (it is not a positive definite matrix because it might be singular).

The Levenberg-Marquardt algorithm consists in rescaling the diagonal elements of $\mathbf{H}^{[n]}$ by the damping factor $1 + \lambda$. The system (Eq. C.7) becomes

$$ [\mathbf{H} + \lambda \text{diag}(\mathbf{H})]^{[n]} \Delta \mathbf{p}^{[n]} = \mathbf{J}^{[n]T} \mathbf{W} [\mathbf{RV}_{\text{obs}} - f(\mathbf{t}, \mathbf{p}^{[n]})]; \qquad (C.8) $$

where $\lambda$ is the damping factor, chosen to be large when H approaches the singularity.

## Appendix C.2. Separable variable technique

The parameters $\mathbf{p}_1 = [\gamma, \alpha, \beta]^T$ are linear because they satisfy the following property:

$$ f(\mathbf{t}, \mathbf{p}) = \sum_{k=1}^{3} \frac{\partial f(\mathbf{t}, \mathbf{p})}{\partial \mathbf{p}_1(k)} \mathbf{p}_1(k) = \mathbf{J}_1 \mathbf{p}_1. \qquad (C.9) $$

Hereafter, subscripts l and nl refer to linear and non-linear respectively. Then, $\mathbf{J}_1$ is the Jacobian matrix whose columns are the derivatives of f with respect to $\mathbf{p}_1$.

While in both Newton-Raphson and the (commonly used) Levenberg-Marquardt (LM) algorithms, the linear and non-linear parameters are treated the same way, the separable variable technique has the advantage that $\mathbf{p}_1$ can be obtained for every non-linear parameter $\mathbf{p}_{nl}^{[n]} = [v_k, e, T_0]^{[n]T}$ by a linear least-square regression, such that 16

$$ \mathbf{H}_1^{[n]} \mathbf{p}_1^{[n]} = \mathbf{Z}_1^{[n]}, \qquad (C.10) $$

<footnote n="16">16 The inversion of the normal equations is known to be numerically unstable and should be avoided in favour of orthogonal decomposition approaches.</footnote>
with

$$ \mathbf{H}_1^{[n]} = \mathbf{J}_1^{[n]T} \mathbf{W} \mathbf{J}_1^{[n]}, \qquad (C.11) $$

$$ \mathbf{Z}_1^{[n]} = \mathbf{J}_1^{[n]T} \mathbf{W} \mathbf{R} \mathbf{V}_{\text{obs}}. \qquad (C.12) $$

Therefore, Eq. C.2 becomes a function of $\mathbf{p}_{nl}^{[n]}$ only, such that

$$ \mathbf{R} \mathbf{V} = f(\mathbf{t}, \mathbf{p}^{[n]}) + \tilde{\mathbf{J}}_{nl}^{[n]} \Delta \mathbf{p}_{nl}^{[n]}, \qquad (C.13) $$

where the step of the non-linear parameters is $\Delta \mathbf{p}_{nl}^{[n]} = \mathbf{p}_{nl}^{[n+1]} - \mathbf{p}_{nl}^{[n]}$.

The Jacobian matrix $\tilde{\mathbf{J}}_{nl}^{[n]}$ encompasses both the derivatives of $f$ with respect to $\mathbf{p}_{nl}^{[n]}$ and $\mathbf{p}_1^{[n]}$

$$ \tilde{\mathbf{J}}_{nl}^{[n]} = \frac{\partial f(\mathbf{t}, \mathbf{p}^{[n]})}{\partial \mathbf{p}_{nl}^{[n]}} + \frac{\partial f(\mathbf{t}, \mathbf{p}^{[n]})}{\partial \mathbf{p}_1^{[n]}} \frac{\partial \mathbf{p}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}} \qquad (C.14) $$

$$ = \mathbf{J}_{nl}^{[n]} + \mathbf{J}_1^{[n]} \mathbf{A}^{[n]}. \qquad (C.15) $$

$\mathbf{A}^{[n]} = \frac{\partial \mathbf{p}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}}$ is a 3×3 matrix containing the derivatives of the linear parameters $\mathbf{p}_1^{[n]}$ with respect to the non-linear ones $\mathbf{p}_{nl}^{[n]}$. It is obtained by computing the derivative of Eq. (C.10) with respect to $\mathbf{p}_{nl}^{[n]}$

$$ \mathbf{A}^{[n]} = [\mathbf{H}_1^{[n]}]^{-1} \left[ \frac{\partial \mathbf{Z}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}} - \frac{\partial \mathbf{H}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}} \mathbf{p}_1^{[n]} \right]. \qquad (C.16) $$

Using Eqs. (C.11) and (C.12), we deduce the derivatives as

$$ \frac{\partial \mathbf{H}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}} = \frac{\partial \mathbf{J}_1^{[n]T}}{\partial \mathbf{p}_{nl}^{[n]}} \mathbf{W} \mathbf{J}_1^{[n]} + \mathbf{J}_1^{[n]T} \mathbf{W} \frac{\partial \mathbf{J}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}} \qquad (C.17) $$

$$ \frac{\partial \mathbf{Z}_1^{[n]}}{\partial \mathbf{p}_{nl}^{[n]}} = \frac{\partial \mathbf{J}_1^{[n]T}}{\partial \mathbf{p}_{nl}^{[n]}} \mathbf{W} \mathbf{R} \mathbf{V}_{\text{obs}}. \qquad (C.18) $$

Therefore, the system (Eq. C.8) becomes

$$ [\tilde{\mathbf{H}}_{nl} + \lambda \text{diag}(\tilde{\mathbf{H}}_{nl})]^{[n]} \Delta \mathbf{p}_{nl}^{[n]} = \tilde{\mathbf{J}}_{nl}^{[n]T} \mathbf{W} [\mathbf{R} \mathbf{V}_{\text{obs}} - f(\mathbf{t}, \mathbf{p}^{[n]})] \qquad (C.19) $$

with $\tilde{\mathbf{H}}_{nl} = \tilde{\mathbf{J}}_{nl}^{[n]T} \mathbf{W} \tilde{\mathbf{J}}_{nl}^{[n]}$ being the new hessian matrix of the system. We have then split the 6×6 system (Eq. C.8) into two smaller systems, each one of 3×3, in Eqs. (C.10) and (C.19), respectively.

## Appendix C.3. Analytical expressions of the involved matrices

To perform the above computations, we need the expressions of the Jacobian matrices $\mathbf{J}_1^{[n]}$ and $\mathbf{J}_{nl}^{[n]}$. The former, relative to the linear system of Eq. (C.9) is given by

$$ \mathbf{J}_1^{[n]} = [1, \cos(\mathbf{v}) + e, \sin(\mathbf{v})]^{[n]}. \qquad (C.20) $$

Its derivatives with respect to the non-linear parameters are

$$ \frac{\partial \mathbf{J}_1^{[n]}}{\partial v_k} = [0, -\sin(\mathbf{v}), \cos(\mathbf{v})]^{[n]} \frac{\partial \mathbf{v}^{[n]}}{\partial v_k}, \qquad (C.21) $$

$$ \frac{\partial \mathbf{J}_1^{[n]}}{\partial e} = \left[ 0, 1 - \sin(\mathbf{v}) \frac{\partial \mathbf{v}}{\partial e}, \cos(\mathbf{v}) \frac{\partial \mathbf{v}}{\partial e} \right]^{[n]}, \qquad (C.22) $$

$$ \frac{\partial \mathbf{J}_1^{[n]}}{\partial T_0} = [0, -\sin(\mathbf{v}), \cos(\mathbf{v})]^{[n]} \frac{\partial \mathbf{v}^{[n]}}{\partial T_0}. \qquad (C.23) $$

Very similarly, the Jacobian matrix $\mathbf{J}_{nl}^{[n]}$ is given by

$$ \frac{\partial \mathbf{J}_{nl}^{[n]}}{\partial v_k} = [0, -\alpha \sin(\mathbf{v}), \beta \cos(\mathbf{v})]^{[n]} \frac{\partial \mathbf{v}^{[n]}}{\partial v_k}, \qquad (C.24) $$

$$ \frac{\partial \mathbf{J}_{nl}^{[n]}}{\partial e} = \left[ 0, \alpha \left( 1 - \sin(\mathbf{v}) \frac{\partial \mathbf{v}}{\partial e} \right), \beta \cos(\mathbf{v}) \frac{\partial \mathbf{v}}{\partial e} \right]^{[n]}, \qquad (C.25) $$
$$ \frac{\partial \mathbf{J}_{nl}^{[n]}}{\partial T_0} = [0, -\alpha \sin(\mathbf{v}), \beta \cos(\mathbf{v})]^{[n]} \frac{\partial \mathbf{v}^{[n]}}{\partial T_0}. \qquad (C.26) $$

The derivatives of the true anomalies v with respect to the three non-linear parameters are

$$ \frac{\partial \mathbf{v}}{\partial v_k} = 2 \pi (t - T_0) (1 + e \cos(\mathbf{v}))^2 / (1 - e^2)^{3/2}, \qquad (C.27) $$

$$ \frac{\partial \mathbf{v}}{\partial e} = \sin(\mathbf{v}) (2 + e \cos(\mathbf{v})) / (1 - e^2), \qquad (C.28) $$

$$ \frac{\partial \mathbf{v}}{\partial T_0} = -2 \pi v_k (1 + e \cos(\mathbf{v}))^2 / (1 - e^2)^{3/2}. \qquad (C.29) $$

## Appendix C.4. Computation steps

The parameter step $\Delta \mathbf{p}_{nl}^{[n]}$ in Eq. C.19 is computed as follows:

1. Given the preliminary non-linear parameters $\mathbf{p}_{nl}$, we compute the mean and eccentric anomalies given by Eq. 4. This transcendental equation is solved following the algorithm proposed by Fukushima (1996), and implemented in a common Gaia tool library.

2. The true anomaly trigonometric functions are computed through Eq. 2 and Eq. 3.

3. The Jacobian matrix $\mathbf{J}_1^{[n]}$ is filled given Eq. C.20. The linear parameters $\mathbf{p}_1^{[n]}$ are then derived by the linear least-square regression given by Eq. C.10.

4. The derivatives of the true anomalies are computed using Eq. C.27, Eq. C.28, and Eq. C.29.

5. We fill the Jacobian matrix $\mathbf{J}_{nl}^{[n]}$ following Eq. C.24, Eq. C.25, and Eq. C.26.

6. The matrix A is computed thanks to Eq. C.17, Eq. C.18, and Eq. C.16.

7. Finally, the total Jacobian matrix $\tilde{\mathbf{J}}_{nl}^{[n]}$ is obtained by Eq. C.15, and the linear system (Eq. C.19) is solved for $\Delta \mathbf{p}_{nl}^{[n]}$.