Gaia Data Release 3. Astrometric binary star processing

Since the start of the Gaia satellite observations in summer 2014, several data releases have been made available to the scientific community. The astrometric parameters of more than one billion stars have been determined three times with ever increasing precision, but these data still described the stars as single. The third Data Release (DR3) and the Early third Data Release (EDR3) are the first to have a significant number of observations covering a sufficiently long time span (about 33 months) to make it worthwhile to use Gaia observations to search for binary stars.

The extension of the exploitation of Gaia observations to binary stars reveals two more difficulties than the application of the single star model. As with its predecessor Hipparcos, the double star candidates are covered by different numerical models, without being certain that any of them really corresponds to a given star. Therefore, we have to solve the problem of choosing the model that best applies to a given star, as well as the final acceptance of the solution. Before that, we have to filter the input data to retain only those stars that may be able to get a suitable solution.

The astrometric processing of the Gaia Early Data Release 3 (EDR3, Gaia Collaboration, 2021) left some 36.5 million stars that were considered to be possibly double. These stars were all brighter than G=19, had a renormalized unit weight error ($RUWE$; see Lindegren et al., 2021) greater than 1.4, and were observed over at least 12 visibility periods. These objects included many pairs of stars that were resolved whenever the orientation of the scan direction allowed, and seen as a single photocentre otherwise. These so-called “partially resolved double stars” are likely to appear as unresolved binary stars of extremely rapid orbital motion, and should be discarded as far as possible. It is planned to apply a specific treatment to them in the next DRs, but for the moment we can only eliminate them when possible. For this purpose, two quantities published in Gaia EDR3 and coming from the astrometric binary detection processing (Lindegren et al., 2021) were used to establish two additional criteria: The percentage of CCD observations where image parameter determination (IPD) detected more than one peak, called $ipd\_frac\_multi\_peak$, must be less than or equal to 2, and $ipd\_gof\_harmonic\_amplitude$, which is the amplitude of the natural logarithm of the goodness-of-fit obtained in the IPD versus position angle of scan, must be less than 0.1. As a result, the number of targets to be treated was reduced to 10.9 million. In order to further reduce the number of partially resolved double stars, a final selection was based on photometric properties. As the pixels of the astrometric field are much smaller than those of the BP and RP windows, a component of a close binary at the resolution limit appears less bright in $G$ magnitude (from the astrometric field) than one would expect from its magnitude in the BP and RP fields. The corrected BP and RP flux excess factor $C^{*}$ and its uncertainty $\sigma_{C^{*}}$ as defined in Riello et al. (2021, Eq. 6 and 18) were taken into account. The drastic condition $|C^{*}|<1.645\;\sigma_{C^{*}}$, corresponding to a 90% confidence interval, was applied and a selection of $4\,115\,743$ stars was finally obtained.

The organisation of the astrometric binary star processing, as it was finally realised, is summarised in Fig. 1. The input data consist of the $4\,115\,743$ selected stars, for which the following information is used: the along-scan (AL) coordinates of each transit of the astrometric field of view, the partial derivatives of the single star astrometric model, the $G$ magnitudes of the transits, and the parameters of a preliminary single star solution, i.e. the coordinate offsets, the parallax and the proper motion. The production of the astrometric data is described in Lindegren et al. (2021), while the magnitude reduction is explained in Riello et al. (2021). The data are complemented by radial velocities, when available (Seabroke et al., 2021).

The calculation of the astrometric solutions is done in the so-called main processing, as explained hereafter. The data discussed above are ingested by the preprocessing, which prepares them for the calculation of different types of solutions. The solutions are calculated one after the other, starting with a new calculation of the single star solution, then the acceleration solutions, the orbital solution, and, finally, the variability-induced mover (VIM) solutions if the star is photometrically variable. As soon as a solution is considered acceptable, the processing cascade is interrupted and the star goes to post-processing. The post-processing essentially consists of re-filtering the solutions, applying more stringent criteria than those used after the solutions were calculated. When a solution is discarded by post-processing, it is replaced by the EDR3 single star solution.

In practice, the software for managing the cascade of calculation and acceptance of solutions is much more sophisticated than the scheme shown in Fig. 1. In fact, the solution that entered the post-processing was sometimes a so-called “alternative” solution (not to be confused with the ”OrbitalAlternative” solutions in Holl et al., 2022): a solution that met poor acceptance criteria, and which was retained because it was the solution that best matched the observations at the end of the cascade. More information on this aspect of the process is given at the end of Section 2.2.2, but we can already highlight two things: First, that we did not systematically try all the models to retain the one that was the most suitable; this was motivated mainly by the concern to save computational time, but also by the fear of accepting complex solutions for objects that could fall under a simpler model (for example, giving an acceleration solution for an orbital binary is a less serious error than giving a false orbital solution for a binary that should only have an acceleration solution). Despite its apparent mathematical simplicity, the VIMF model has been placed at the end because it combines data of two kinds (astrometric and photometric), which undermines the reliability of the solutions. Secondly, post-processing consists only of rejecting solutions, without trying to replace them with more complex or alternative solutions. Ideally, the filters applied in post-processing should be included in the main processing, but this was not possible in the time available.

We define hereafter two quality criteria, the goodness of fit and the significance, that can distinguish acceptable solutions from those that are irrelevant or at least potentially wrong. In order to decide whether a solution is really plausible, we will also consider astrophysical criteria and other constraints, which will allow us to identify some solutions as not very plausible. These bad solutions will allow us to delimit a quality domain where they are absent or rare, and whose solutions will be kept.

Each transit of a star through the Gaia astrometric instrument comes down to the measurement of two coordinates giving its position relative to a reference point assigned to the star (see Lindegren et al., 2021). Of these, only the coordinate along the scan axis, called hereafter the AL abscissa, $w$, is taken into account, being considerably more accurate than the other. The AL abscissae were measured along the 9 CCD transits that constitute each field-of-view (FoV) transit. As the duration of the FoV transit is only a few seconds while we are looking for binary stars with periods of at least several days, the 9 CCD transits amount to 9 measurements of the same quantity. Outliers were therefore detected by intercomparison between the CCD transits. In practice, each AL abscissa was compared to the median of the 9 values, and the transit CCD was rejected when the difference exceeded 5 times the uncertainty of the abscissa.

This rejection process was not the only one that was applied, and subsequently transits were rejected after the calculation of some astrometric solutions. The largest residual transit was rejected when it exceeded 5 times the uncertainty, and the solution was recalculated; however, this operation was subject to three restrictions: (1) it was only done for linear models, i.e. the single star model and the two acceleration models, (2) the $\chi^{2}$ of the solution should exceed the number of transits used multiplied by 1.41 and (3) the proportion of rejected transits was limited to 5 %.

The distance from the Sun and the apparent position of a nearby star are varying during the time-span of the Gaia mission, inducing slight changes in astrometric parameters which are usually assumed to be constant. These so-called “perspective effects” were taken into account in the calculation of the astrometric parameters of a few single stars in the Hipparcos Catalogue (ESA, 1997). Dravins, Lindegren & Madsen (1999) and Lindegren & Dravins (2003) have shown that these effects are so closely related to radial velocity that they could be used to measure it. For Gaia, Halbwachs (2009) has shown that perspective effects do not significantly affect the astrometric orbits of binary stars; however, they do result in a “perspective acceleration” of stars that could be mistaken for an acceleration due to orbital motion.

The perspective acceleration is based on the radial proper motion, $\mu_{r}$, which is defined as the radial velocity, $v_{r}$, divided by distance. In practice, $\mu_{r}$, expressed in yr^-1, is related to $v_{r}$ and parallax, $\varpi$ by the following relationship:

There are then two possibilities to take into account the perspective acceleration: either multiply the partial derivatives by $\mu_{r}(t-T)$, as they are in Eq. 5, or correct the abscissae by subtracting the perspective acceleration contribution, $\Delta w$. Because of its simplicity, we followed the latter method: $\Delta w$ was calculated from Eq. 5, using preliminary values of $\varpi$, $\mu_{\alpha*}$ and $\mu_{\delta}$, and was subtracted to the AL abscissa. This correction was made for transits of all stars of known radial velocity when they were closer than 200 pc, perspective effects being negligible beyond that. The astrometric binary solutions obtained for stars whose measurements have been corrected for perspective acceleration are identified by a flag in the solution tables.

The above operations may have changed the transits sufficiently that the star should no longer be considered a potential binary. For this reason, the single star solutions were recalculated, and they were accepted when $F_{2}$ was zero or negative. Of the $4\,115\,743$ stars that entered the processing chain, only 28 left at this stage. The rest were passed through the acceleration models.

Under the effect of the gravitational force, the two components of a binary are accelerated towards each other. If they cannot be resolved and the photocentre does not coincide with the barycentre, the apparent motion of their photocentre is accelerated; when the orbital period is much longer than the time interval covered by the observations, the acceleration is nearly constant in intensity and direction. These physical considerations are at the origin of the constant acceleration model. Each coordinate of the photocentre varies with time according to a parabola. However, fitting a parabola to the set of observations means setting the mid-mission position of the star as given by the first two terms of the solution, so that it is tangent to the parabola, and therefore far from the mean position of the star. To provide a position closer to the mean position of the star, we have followed the method explained in the Hipparcos catalogue (ESA, 1997, vol. 1, Sect. 2.3.3), and the partial derivatives of the abscissa with respect to the acceleration, $(g_{\alpha*},g_{\delta})$, are given by the following equations:

The significance of the constant acceleration model is defined from the acceleration vector. It is derived from Eq. 3, as explained in Sect. 2.2.2.

When the gravitational force is significantly changing over the duration of the mission, but when the orbital period is still much longer than the observation time span, the time derivative of the acceleration, $(\dot{g}_{\alpha*},\dot{g}_{\delta})$, is added to the parameters of the model. However, the introduction of these parameters could modify the proper motion of the star. In order to derive a proper motion corresponding to the average displacement of the photocentre during the mission, the partial derivatives of the abscissa with respect to $\dot{g}_{\alpha*}$ and $\dot{g}_{\delta}$ are given by the equations:

The significance of the variable acceleration model is defined from the $(\dot{g}_{\alpha*},\dot{g}_{\delta})$ vector. As above, it is derived from Eq. 3, as explained in Sect. 2.2.2.

It is worth noticing that the acceleration models are both linear models, and that the solutions are found by solving a system of linear equations by singular value decomposition.

In the Hipparcos reduction, the constant acceleration solution was considered as acceptable when the acceleration was significant, but the variable acceleration solution was always tried, and, when it was significant, it was preferred to the constant acceleration solution. We have kept this strategy, and, for that reason, the variable acceleration model was tried first. The variable acceleration solution was retained when it satisfied the conditions: $s>12$ and $F_{2}<25$. Otherwise, the constant acceleration solution was derived, and it was accepted when it satisfied the same conditions. Because of this acceptance, no other models were tried. Therefore acceleration solutions were retained for many binary stars with periods shorter than the duration of the Gaia mission. This choice was originally made to save computing time; this reason was no longer relevant at the time of processing the DR3, but the strategy remained unchanged because it was considered that the risk of producing an orbit with a wrong period was more serious than that of giving only a polynomial fit for a trajectory that could be described by an ellipse. This is a matter of judgement which may be considered questionable, and will probably be approached differently in the next data release.

To delimit a domain of acceptance of the solutions, we consider, for both acceleration models, the true acceleration projected on the sky, called $\Gamma$ in the following. $\Gamma$ is derived from the equation:

When the apparent acceleration is expressed in mas.yr^-2 and when the parallax is in mas, $\Gamma$ is obtained in au.yr^-2.

In order to estimate a maximum value for $\Gamma$, we consider a bright star with a dark companion, so that the acceleration of the photocentre is that of the bright star. Applying the law of gravitation, and assuming that the distance between the components is approximately equal to the semi-major axis of the orbit of the bright star around the companion, we deduce from Kepler’s third law that the acceleration is then:

An upper limit of $\Gamma$ is obtained when the bright star is a giant solar-mass star with a companion still on the main sequence or on the subgiant branch, with $q$ slightly less than 1, and when the period is twice the duration of the Gaia mission. Such a system is not frequent, but not exceptional in a sample bounded by a limit in magnitude. With a period of 5.7 yr, the acceleration may be as large as about 2.4 $\mathrm{au/yr^{2}}$ for a circular orbit; over half a period, the average acceleration can then be as high as 2.0 $\mathrm{au/yr^{2}}$. This result is however a rough estimate: We have neglected the eccentricity, which will produce a larger acceleration near the periastron, and we have assumed a period much longer than the duration of the mission, whereas we know that many solutions are for shorter periods. Considering the distribution of values obtained from the solutions (see Fig. 2, (a) and (b)), we see a concentration below $\Gamma=3\;\mathrm{au/yr^{2}}$, which we adopt as the limit for an acceptable acceleration solution. Stars whose acceleration is really beyond this limit are certainly very rare, but they are also objects of great scientific interest, such as binaries with a massive dark component, and for this reason we do not reject them directly.

The following conditions were introduced to restrict the proportion of $\Gamma$ larger than this limit:

1. 1.

   As mentioned in Sect. 4.2 above, the main processing was performed with the threshold $s>12$ for direct acceptance and $s>5$ for the alternative solutions. In addition, conditions based on the significance of the parallax were introduced for both levels: $\varpi/\sigma_{\varpi}\;>\;1.2\;s^{1.05}$ for the constant acceleration model, and $\varpi/\sigma_{\varpi}\;>\;2.1\;s^{1.05}$ for the variable acceleration model. These conditions have been set because they very effectively reduce the amount of high accelerations; it is worth noticing that, in fact, they are not criteria based on the quality of the solutions: reading them from right to left, they both mean that, for a given distance, the most significant solutions are rejected. This results in rejecting the largest accelerations because many of them are highly significant, and their occurence rate is thus artificially reduced. However, it was necessary to apply these criteria to avoid too much pollution by partially resolved double stars. Another beneficial consequence of this filter is that the short period binaries whose acceleration solution was rejected in this way were able to have an orbital solution. Including alternative solutions, the main processing yielded $808\,992$ constant acceleration solutions and $569\,022$ variable acceleration solutions. These solutions must be filtered in the post-processing to keep only those that seem most likely to be real.

2. 2.

   The significance, $s$, must be over 20 for both acceleration models.

3. 3.

   The goodness of fit, $F_{2}$, must be below 22 for the constant acceleration model, and below 25 for the variable acceleration model. The limit is more severe for the constant acceleration solutions so that the sequences of the HR diagrams deduced from the parallaxes of these solutions are always thinner than those deduced from the single star model.

Conditions (2) and (3) reduced the number of solutions to $246\,798$ for constant accelerations, and $91\,227$ for variable accelerations.

The two lower panels of Fig. 2 show the distributions of the selected solutions in the $F_{2}$ vs. $s$ diagrams. Thanks to the parallax filtering, item (1) above, the solutions with $\Gamma>3$ au.yr^-2 are very rare: 0.03 % of the constant acceleration solutions, and 0.3 % for the variable accelerations. If filters (2) and (3) are applied to a test sample unaffected by filter (1), the proportions of large $\Gamma$ are 3.1 and 5.6 % respectively. Unless one is prepared to sacrifice the vast majority of solutions, these high rates resist the quality criteria we have tried; this issue is still being discussed in Sect. 7, in the light of independent investigations. Several origins have been considered for the high acceleration solutions: short-period binaries, partially resolved binaries or outlier astrometric measurements. We can only conclude that the stars with acceleration solutions are most likely genuine binaries, but the physical interpretation of the accelerations is hazardous.

The abscissa of an unresolved astrometric binary can be written as in the following equation:

The derivation of the orbital solutions is based on the Levenberg-Marquardt algorithm, but this can only lead to the solution if it starts from a good starting point. In practice, a calculation starting from $e=0$ often leads to the right result if the period tested is close to the actual period, but scanning a grid covering the space of $P$, $e$ and $T_{0}$ is still necessary to be sure to find the solution. The period range tried is from 10 days to the duration covered by the observations of the star divided by 0.6. The solutions are calculated and published in terms of TI elements, but the calculation of the elements usually used to describe an orbit, which are the so-called Campbell elements (semi-major axis, inclination of the orbit, position angle of the ascending node and periastron longitude), is presented in Appendix A, where the calculations of their uncertainties are also given. The extension of this calculation to the Thiele-Innes elements introduced by both astrometric and spectroscopic orbits is presented in Appendix B.

The significance of the orbital model is defined as the semi-major axis of the orbit, $a_{0}$, divided by its uncertainty. This definition of the significance was adopted since a significant orbit is expected to have a large semi-major axis, when, on the contrary, an orbit with a semi-major axis close to or smaller than its uncertainty may be obtained for a single star. The semi-major axis is not explicitly a parameter of the orbital model that we have used, but it may be derived from the TI elements, thanks to Eq. 21. Its uncertainty is derived from Eq. 26 and Eq. 27.

Rare orbits (8 in total) have been directly calculated as circular and are given with $e=0$. These calculations required special arrangements, which have been also adopted for orbits with very small eccentricities, as explained below.

Preliminary calculations showed that the eccentricities below 0.006 were all at least three times smaller than their uncertainties. For such solutions, the uncertainties of the TI elements are also often very large, as is the uncertainty of $T_{0}$, which can be larger than the period. These anomalies arise from the indeterminacy of the position of the periastron, affecting the TI elements through their dependence on the longitude of the periastron, $\omega$ (see Eq. 20). They do not affect the validity of the solution, for which the Campbell elements calculated as explained in Appendix A have correct uncertainties, except for $\omega$, whose uncertainty far exceeds $2\pi$; however, they may confuse the user.

In order to limit these effects, orbits with eccentricity less than 0.0005 have been “pseudo-circularised” in the following way: The Campbell elements were derived, and the periastron longitude, $\omega$, was set to 0 in order to place the periastron on the ascending node. Therefore, the periastron epoch, $T_{0}$, has been advanced by subtracting the time interval corresponding to the orbit section between the ascending node and the initial periastron, which for a circular orbit is $P\times\omega/(2\pi)$. As the initial orbit is not circular, this formula is an approximation that is only possible for very small eccentricities; this is why it was only done for $e<0.0005$. The TI elements were then recalculated, keeping the semi-major axis, orbit inclination and line-of-nodes position angle of the initial solution.

After these transformations, the solution is based on only 10 unknowns: $e$ and $T_{0}$ are fixed, and only 3 of the 4 TI elements remain, since $\omega=0$ implies that $G=-AF/B$, or $F=-BG/A$. Therefore, $G$ was no longer considered as an unknown provided that $|B|>10^{-5}$ mas. Otherwise, it was $F$ that had to be removed, but this did not happen in practice. Thus, the solution variance-covariance matrix was recalculated for the 10 selected unknowns. This $10\times 10$ matrix was inserted into the $12\times 12$ matrix of orbital solutions, setting the unnecessary terms to 0.

In the end, this gives quite reasonable uncertainties for the TI elements and for $T_{0}$. However, the solution will give an uncertainty and correlation terms that are zero for $G$, which is wrong and needs to be corrected in order to deduce the uncertainties of Campbell’s elements as explained in Appendix A. The uncertainty of $G$, $\sigma_{G}$, is derived from the equation:

The correlation coefficients relating $G$ to the other TI elements are:

These modifications must also be made to circular orbits if the uncertainties of their Campbell elements are to be calculated.

The orbital solutions were first selected according to the thesholds set out in Sect. 2.2.2, the significance being calculated as explained in the previous section. In order to identify doubtful solutions, we have taken into account the mass function, which is derived from:

Panel (a) of Fig. 3 shows the period vs. mass function diagram. There are two kinds of false solutions:

• •

  Solutions of $f_{\cal M}>0.3{\cal M}_{\sun}$. This limit looks rather large for dwarf stars, but roughly acceptable for binaries with a giant component and a component still close to the main sequence, or for a triple system with a massive pair as a secondary component.

• •

  Solutions with particular periods, which probably depend on the 63-day precession period of the satellite. These solutions have usually very large mass function,but they descend sufficiently into the small $f_{\cal M}$ range to contaminate plausible solutions.

A large part of the last category of false solutions was discarded by using a filter based on the parallax significance, as indicated in Fig. 3b. This filter consists of selecting only the solutions that satisfy the following condition:

Although the filter of Eq. 16 has rejected most of the large mass function solutions, a few concentrations around particular periods still persist in the processing results (Fig. 3c and d). Most of these remaining false solutions were eliminated in the post-processing, where only solutions satisfying the following two filters were selected:

• •

  A filter based on the quality of the eccentricity: $\sigma_{e}<0.079\ln P-0.244$, where $P$ is in days.

• •

  A filter based on the significance: $s\;=\;a_{0}/\sigma_{a_{0}}>158/\sqrt{P}$, where $P$ is again in days.

After these operations, about $166\,500$ orbits were retained, and the concentrations along particular periods have all disappeared, as shown in Fig. 3e. Moreover, very few orbits have eccentricities beyond the 10-day circularisation limit, calculated according to Halbwachs, Mayor & Udry (2005). The proportion of $f_{\cal M}>0.3{\cal M}_{\sun}$ is now about 1 %, but it can be seen in Fig. 3f that it increases considerably with the period: for example, it is only 0.2 % for periods of less than 1 year, but reaches 3.3 % for periods longer than 1000 days.

Subsequently, all the main processing astrometric orbits were compared to the SB1 spectroscopic orbits derived from Gaia radial velocities (Gosset et al., 2022). The matching of the orbits of the two types resulted in nearly $30\,950$ combined orbits, of which about 2470 were recovered after having been rejected by the post-processing. At the same time, almost 3500 post-processing orbital solutions were rejected by the validation (Babusiaux et al., 2022), or discarded as redundant with solutions obtained with Markov Chain Monte Carlo and Genetic Algorithms (Holl et al., 2022). In the end, our contribution to the DR3 astrometric orbital solutions amounts to nearly 165500 orbits.

The additional parameters specific to the VIMF model are the coordinates of the photocentre, $(D_{\alpha*},D_{\delta})$, measured from the variable component when the total photometric flux of the system is equal to a reference flux, noted $\bar{\cal F}$. In practice, $\bar{\cal F}$ is the median flux of all transits of the system. The partial derivatives of the abscissa with respect to $D_{\alpha*}$ an $D_{\delta}$ are given by the following equations:

The parameters of the VIMF model are thus derived from a linear system that is solved by singular value decomposition, but the uncertainty of each transit deserves a dedicated calculation, since it depends on the astrometric uncertainty of the measured abscissa, $\sigma_{w}$, but also on the uncertainty of the photometric flux, $\sigma_{\cal F}$, on which an error from the model depends. For each transit, this error, $\sigma_{\mathrm{mod}}$, is calculated from the equation:

The total uncertainty of transit is then given by the equation:

The $\sigma_{\mathrm{VIMF}}$ uncertainties are used to derive any VIMF solution, as well as $F_{2}$ and the correction coefficient $c$ as explained in Sect. 2.2.1. As $\sigma_{\mathrm{VIMF}}$ depends on $D_{\alpha*}$ and $D_{\delta}$, and thus on the VIMF solution, the calculation is iterative and the VIMF model cannot be considered linear in the strict sense. The significance of the solution is defined from the modulus of the $(D_{\alpha*},D_{\delta})$ vector, as explained in Sect. 2.2.2.

As far as the VIML model is concerned hereafter, suffice it to say that it has two more parameters, which are the time derivatives of $D_{\alpha*}$ and $D_{\delta}$.

The main criterion we have for estimating whether a VIMF solution is plausible is the apparent distance $D$ that separates the median photocentre from the variable component. As the resolution of Gaia is at least 180 mas (Lindegren et al., 2021) and the components of a VIM should never be separated, $D$ must remain under a limit of no more than a hundred mas; we are therefore looking here for quality criteria such that this condition is met.

The $F_{2}$ vs. $D$ and significance vs $D$ diagrams obtained after a first processing trial do not immediately show any limits in $F_{2}$ or significance that would allow almost only $D<100$ mas. Fortunately, a selection becomes possible when we consider the “parallax significance”, i.e. the quantity $\varpi/\sigma_{\varpi}$. It is clear from Fig. 4 that the selection of solutions of parallax significance greater than 30 rejects most excessive values of $D$. Therefore, this criterion was applied in the main processing, in addition to the criteria $s>12$ and $F_{2}<25$. As before, the criteria for selecting alternative solutions were $s>5$ and $F_{2}$ less than that of the alternative solution that may have already been selected.

Two thousand five hundred and eight VIMF solutions were obtained from the main processing. The $F_{2}$ vs. $D$ diagram of these stars is shown in Fig. 5a. As for the previous models, solutions with $F_{2}$ higher than 25 seem doubtful, because the proportion of large $D$ values in them increases, and, even more so, small $D$ values become increasingly rare as $F_{2}$ increases. They are therefore rejected, which leaves 1660 solutions. The significance vs. $D$ diagram of these solutions shows a large scattering of $D$ values when significance is below 20, especially below 12. However, the threshold was set at 20 as a precaution. Beyond that, $D$ increases with $D/\sigma_{D}$ (Fig. 5b), which corresponds to a value of $\sigma_{D}$ that is approximately constant and close to one tenth of mas, which seems acceptable.

The 869 VIMF solutions thus retained in the end have distances $D$ below 40 mas, with three exceptions: three solutions of significance close to 100, when $D$ is 55, 56 and 101 mas. These values are not unrealistic, and it seems possible that all VIMF solutions are valid. Unfortunately, this may also be due to the low effectiveness of using $D$ to point out dubious solutions, as the application of other VIM models has shown below.

As explained at the beginning of this section, four other VIM models were tried, from VIML to VIMO. Contrarely to the VIMF model, the parameters of these models allow the estimation of a minimum value for the mass of the systems. The result was devastating: for all these models, the minimum masses are concentrated between 10 and $10\;000$ ${\cal M}_{\sun}$. We therefore decided to remove these models from the mass processing cascade.