11institutetext: Université de Strasbourg, CNRS, Observatoire astronomique de
Strasbourg, UMR 7550, 11 rue de l’Université, 67000 Strasbourg, France
22institutetext: FNRS, Institut d’Astronomie et d’Astrophysique, Université
Libre de Bruxelles, boulevard du Triomphe, 1050 Bruxelles, Belgium
33institutetext: GEPI, Observatoire de Paris, Université PSL, CNRS, 5 Place
Jules Janssen, 92190 Meudon, France 44institutetext: Université Côte d’Azur,
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Bd de
l’Observatoire, CS 34229, 06304 Nice Cedex 4, France 55institutetext: IMCCE,
Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Lille,
77 av. Denfert-Rochereau, 75014, Paris, France 66institutetext: Telespazio UK
S.L. for ESA/ESAC, Camino bajo del Castillo, s/n, Urbanizacion Villafranca del
Castillo, Villanueva de la Cañada, 28692, Madrid, Spain

#  Gaia Data Release 3. Astrometric binary star processing

Jean-Louis Halbwachs 11 Dimitri Pourbaix deceased on 14 November 202122
Frédéric Arenou 33 Laurent Galluccio 44 Patrick Guillout 11 Nathalie Bauchet
55 Olivier Marchal 11 Gilles Sadowski 22 David Teyssier 66 jean-
<EMAIL_ADDRESS>
(Received/Accepted)

###### Abstract

Context. The Gaia Early Data Release 3 contained the positions, parallaxes and
proper motions of 1.5 billion sources, among which some did not fit well the
”single star” model. Binarity is one of the causes of this.

Aims. Four million of these stars were selected and various models were tested
to detect binary stars and to derive their parameters.

Methods. A preliminary treatment was used to discard the partially resolved
double stars and to correct the transits for perspective acceleration. It was
then investigated whether the measurements fit well with an acceleration model
with or without jerk. The orbital model was tried when the fit of any
acceleration model was beyond our acceptance criteria. A Variability-Induced
Mover (VIM) model was also tried when the star was photometrically variable. A
final selection has been made in order to keep only solutions that probably
correspond to the real nature of the stars.

Results. At the end, $338\,215$ acceleration solutions, about $165\,500$
orbital solutions and 869 VIM solutions were retained. In addition, formulae
for calculating the uncertainties of the Campbell orbital elements from
orbital solutions expressed in Thiele-Innes elements are given in an appendix.

###### Key Words.:

binaries:general – Catalogs – Astrometry – Methods: data analysis –
Astronomical instrumentation, methods and techniques

††offprints: J.-L. Halbwachs,

##  1 Introduction

###  1.1 Context

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.

###  1.2 Selection of stars to be processed

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.

##  2 Treatment overview

Figure 1: Overall organisation of the astrometric treatment of binary stars,
as it was eventually applied to the DR3. The cascade on the left side of the
figure is the so-called “main processing” hereafter.

###  2.1 Overall organisation

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.

###  2.2 Criteria for the acceptance of solutions

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.

####  2.2.1 The goodness of fit

The goodness of fit indicates if the model provides predictions which are
compatible with the actual measurements. As in the preparation of the
Hipparcos catalogue (ESA, 1997, vol. 1, Sect. 2.1), we use the $F_{2}$
estimator, which is given by the formula:

$F_{2}=\sqrt{\frac{9\nu}{2}}\left[\left(\frac{\chi^{2}}{\nu}\right)^{1/3}+\frac{2}{9\nu}-1\right]$
(1)

where $\nu$ is the number of degrees of freedom, and $\chi^{2}$ is the sum of
the squared normalised residuals. For a solution derived through an adequate
model, $F_{2}$ is expected to obey the normal distribution ${\cal N}(0,1)$.
This property is attributed to linear models (Stuart & Ord, 1994), but we
verified by simulations that it also applies to the orbital model presented
hereafter, when the semi-major axis of the orbit is significantly larger than
its uncertainty. A large $F_{2}$ means that the model is inadequate, or that
the uncertainties used to derive the $\chi^{2}$ are underestimated. In
practice, we know that this last possibility widely applies to the astrometric
transits used in DR3. Therefore, only very large values of $F_{2}$ indicate
that the model cannot be accepted. Nevertheless, assuming that the model is
adequate, the uncertainties may be revised in order to obtain a zero $F_{2}$.
The simplest correction method, since it keeps the relative weights of the
transits and therefore the solution of the model, consists in multiplying the
uncertainties of the observations by the coefficient:

$c=\sqrt{\frac{\chi^{2}}{\nu\left(1-\frac{2}{9\nu}\right)^{3}}}$ (2)

It is worth noticing that the correction by the coefficient $c$ also applies
to the uncertainties of the model parameters. In addition to the correction of
uncertainties, $F_{2}$ was also taken into account for the selection of
solutions. In the main processing, solutions of $F_{2}>1000$, if any, were
rejected, but a stricter condition was applied in the post-processing:
whatever the model, purely astrometric solutions with uncorrected $F_{2}>25$
were considered as questionable, and were finally not retained.

Although similar in principle, this approach is slightly different from that
followed in the treatment of single star solutions for EDR3 (Lindegren et al.,
2021): Eq. 2 shows that the coefficient $c$ is equal to
$\rm{UWE}\times[1-(2/9\nu)]^{-3/2}$, which is approximately equal to
$\rm{UWE}$. On the other hand, we do not apply a renormalisation, but this
does not matter since we only consider stars brighter than $G=19$.

####  2.2.2 The significance

The significance is a dimensionless quantity that was already introduced in
the Hipparcos catalogue (ESA, 1997, vol. 1, Sect. 2.3.3) in order to decide
whether the use of a model including additional parameters is really
justified. When the additional parameters are the coordinates of a two-
dimensional vector, as is the case for the acceleration models in Sect. 4, the
significance, $s$, is defined as the module of this vector divided by its
uncertainty. Therefore, $s$ is calculated with the equation:

$s=\frac{1}{\sigma_{1}\sigma_{2}}\sqrt{\frac{p_{1}^{2}\sigma_{2}^{2}+p_{2}^{2}\sigma_{1}^{2}-2p_{1}p_{2}\rho_{12}\sigma_{1}\sigma_{2}}{1-\rho^{2}_{12}}}$
(3)

where $p_{1}$ and $p_{2}$ are the coordinates of the additional vector
characterising the model (for instance, the acceleration), $\sigma_{1}$ and
$\sigma_{2}$ are the uncertainties of $p_{1}$ and $p_{2}$, respectively,
corrected by the coefficient $c$ derived with Eq. 2, and $\rho_{12}$ the
correlation coefficient between $p_{1}$ and $p_{2}$. The values of
$\sigma_{1}$, $\sigma_{2}$ and $\rho_{12}$ are all taken from the variance-
covariance matrix.

When the solutions were calculated, a preliminary selection was made through
basic filtering: the solutions of $s>12$ and $F_{2}<25$ were considered good
enough to be retained, without trying other models. Otherwise, the solutions
with $s>5$ and $F_{2}<1000$ were considered as “alternative” solutions (see
Sect. 2.1), and the smaller $F_{2}$ alternative was provisionally accepted at
the end of the cascade. The limit of 5 was adopted empirically on the basis of
initial tests of the data. Thresholds of 5 and 12, which amount to as many
$\sigma$s, seem exceptionally high and are much higher than the thresholds of
3 and 4 that were adopted on the basis of simulations for alternative
solutions and for direct acceptance respectively. This is because, in reality
the uncertainties of the transits of a star are certainly not all
overestimated by the same coefficient; correcting the uncertainties by a
single coefficient is only a stopgap measure, as it is not possible to
separate transits with a correct original uncertainty from those with a very
underestimated uncertainty.

The set of solutions thus selected was then subjected to a more severe
filtering in order to keep only the most credible solutions, as explained
below for each model. The final selection criteria are presented in Table 1.
For the acceleration models and for the VIMF model, they are more severe than
the direct acceptance criteria, and they lead to the rejection of alternative
solutions. However, all alternative solutions contribute to the statistical
properties of the main processing solutions, which are presented hereafter,
and which are the basis for the choice of the final criteria.

Table 1: Properties of the different models and final conditions for the
selection of solutions.

Model | Dimension | Significance | Selection conditions
---|---|---|---
significance | $F_{2}$ | $\varpi/\sigma_{\varpi}$ | Other
Acceleration : |  |  |  |  |  |
Constant | 7 | $s_{7}=\frac{g}{\sigma_{g}}$ | $>20$ | $<22$ | $>1.2s_{7}^{1.05}$ | -
Variable | 9 | $s_{9}=\frac{\dot{g}}{\sigma_{\dot{g}}}$ | $>20$ | $<25$ | $>2.1s_{9}^{1.05}$ | -
Orbital∗ : |  |  |  |  |  |
eccentric | 12 | $\frac{a_{0}}{\sigma_{a_{0}}}$ | $>\mathrm{max}\left(5,\frac{158}{\sqrt{P_{\mathrm{days}}}}\right)$ | $<25$ | $>\frac{20\;000}{P_{\mathrm{days}}}$ | $\sigma_{e}<0.079\ln P_{\mathrm{days}}-0.244$
circular or | 10 | $\frac{a_{0}}{\sigma_{a_{0}}}$ | $>\mathrm{max}\left(5,\frac{158}{\sqrt{P_{\mathrm{days}}}}\right)$ | $<25$ | $>\frac{20\;000}{P_{\mathrm{days}}}$ | -
pseudo-circular |
VIMF | 7 | $\frac{D}{\sigma_{D}}$ | $>20$ | $<25$ | $>30$ | -

*: These selection conditions do not apply to orbital solutions from main processing when they have been confirmed by a SB1 spectrocopic solution, as explained in Gosset et al. (2022).

##  3 The preprocessing

###  3.1 Rejection of outliers

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 %.

###  3.2 Correction of the perspective acceleration

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 milliarcseconds (mas) per year, is related
to $v_{r}$ and parallax, $\varpi$ by the following relationship:

$\mu_{r}=v_{r}\;\varpi\times\frac{24\;\pi\times\mathrm{year_{days}}}{180\times\mathrm{au_{m}}}$
(4)

where $v_{r}$ is in km.s-1, $\varpi$ is in mas, $\mathrm{year_{days}}$ is the
duration of the year in days, and $\mathrm{au_{m}}$ is the length of the
astronomical unit (au) in meters. Perspective acceleration changes the AL
abscissa of a star, $w$, by adding the following quantity:

$\Delta w=-\mu_{r}\;(t-T)\times\;\left(\frac{\partial
w}{\partial\varpi}\;\varpi+\frac{\partial
w}{\partial{\mu_{\alpha*}}}\;\mu_{\alpha*}+\frac{\partial
w}{\partial{\mu_{\delta}}}\;\mu_{\delta}\right)$ (5)

where $t$ is the epoch in years, and $T$ is an epoch close to the middle of
the DR3 mission (in practice, $T=2016.0$). The partial derivatives of the
abscissa with respect to the parallax and to the proper motion coordinates are
those of the 5-parameter single star model. If $\mu_{r}$ is unknown, then the
abscissae are related to the parameters by a system of non-linear equations.
Fortunately, the radial velocities of the bright stars are measured by Gaia,
and a calculation based on the single star model, ignoring the perspective
acceleration, already gives a good approximation of the parallax and of the
proper motion. Thus, $\mu_{r}$ may be determined from Eq. 4.

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.

###  3.3 New calculation of the single star solution

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.

##  4 The acceleration solutions

Acceleration models were already introduced in the reduction of the Hipparcos
catalogue in order to describe binary stars that could not be covered by the
orbital model because their period was too long. Two models were considered:
the constant acceleration model, where the trajectory is a parabola, and the
variable acceleration model, which includes the time derivatives of the
acceleration. The calculation of the solutions from these models was done with
rejection of outliers, as explained at the end of Sect. 3.1.

###  4.1 The constant acceleration model

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:

$\begin{array}[]{ll}\frac{\partial w}{\partial
g_{\alpha*}}&=\frac{1}{2}\frac{\partial
w}{\partial\alpha*}\;\left[(t-T)^{2}-\frac{\Delta T^{2}}{3}\right]\\\ \\\
\frac{\partial w}{\partial g_{\delta}}&=\frac{1}{2}\frac{\partial
w}{\partial\delta}\;\left[(t-T)^{2}-\frac{\Delta T^{2}}{3}\right]\\\
\end{array}$ (6)

where $T$ is 2016.0 as in Sect. 3.2, and $\Delta T$ is half the time interval
covered by the observation of the star. For simplicity, the same value of
$\Delta T$ was assumed for all stars, and we have adopted half the time
between the start of the Gaia science mission and the last DR3 observation,
which, according to the information available when the software was finalised,
lead to $\Delta T=1035/2=517.5$ d, or 1.417 yr. This is an overestimation,
however: in fact, the last observation of the DR3 was made on 28 May 2017, 3
days later than we had assumed, but, on the other hand, the astrometric data
of the first month of the Gaia scientific mission were not taken into account
(Lindegren et al., 2021). Moreover, the actual duration covered by the
observations of a star is necessarily shorter than the duration covered by all
the observations. For all these reasons, $\Delta T$ should be rather of the
order of 470 days. However, the resulting shift in the mean position is very
small: 0.06 mas times the acceleration in mas.yr-2.

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.

###  4.2 The variable acceleration model

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:

$\begin{array}[]{ll}\frac{\partial
w}{\partial\dot{g}_{\alpha*}}&=\frac{1}{6}\frac{\partial
w}{\partial\alpha*}\;\left[(t-T)^{2}-\Delta T^{2}\right](t-T)\\\ \\\
\frac{\partial w}{\partial\dot{g}_{\delta}}&=\frac{1}{6}\frac{\partial
w}{\partial\delta}\;\left[(t-T)^{2}-\Delta T^{2}\right](t-T)\end{array}$ (7)

where $\Delta T$ is the same as in Sect. 4.1 above. We emphasise that the only
effect of introducing $\Delta T$ is to produce a position and proper motion
similar to that which would result from applying the single star model. Since
acceleration models cannot locate the barycentre of a binary, these parameters
are still affected by the orbital motion of the photocentre relative to the
barycentre.

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.

###  4.3 Selection of the acceleration solutions

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:

$\Gamma=\frac{\sqrt{g_{\alpha*}^{2}+g_{\delta}^{2}}}{\varpi}$ (8)

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:

$\Gamma_{\mathrm{au/yr^{2}}}\approx
4\pi^{2}\left(\frac{{\cal{M}_{\cal{M}_{\odot}}}}{P_{\mathrm{yr}}^{4}}\frac{q^{3}}{(1+q)^{2}}\right)^{1/3}$
(9)

where $P$ is the period, $\cal{M}$ the mass of the bright star and $q$ the
mass ratio, i.e. the ratio between the mass of the dark companion and
$\cal{M}$.

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.

Figure 2: Selection of the acceleration solutions. The left-hand panels
correspond to solutions with constant acceleration, and those on the right to
those with variable acceleration. Panels (a) and (b): Significance vs.
acceleration diagrams from a processing trial without the selection based on
parallax error. The solutions with $F_{2}$ below the final threshold are in
blue, while those above are in orange. The solutions of $5<s<12$ were accepted
in the absence of better solutions from another model, while the others were
immediately accepted when $s>12$ and $F_{2}$ below 22 (for constant
acceleration solutions) or 25 (for variable acceleration solutions). The
amount of solutions with a mean acceleration $\Gamma>3\;\mathrm{au/yr^{2}}$
falls off when the significance exceeds 20. Panels (c) and (d): the $F_{2}$
vs. $s$ diagrams of the solutions finally selected in the post-processing. The
black dots are the remaining solutions with $\Gamma>3\;\mathrm{au/yr^{2}}$;
they still correspond to the smallest significances, but this essentially due
to the filtering based on the parallax error.

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.

##  5 The orbital solutions

An orbital solution is calculated when the orbit has been observed either in
its entirety or over a sufficient portion to extrapolate the missing part.

###  5.1 Calculation of the orbital solutions

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

$\begin{array}[]{ll}w=&w_{B}+\frac{\partial w_{B}}{\partial\delta}(\cos
E-e)\times A+\frac{\partial w_{B}}{\partial\alpha*}(\cos E-e)\times B\\\
&+\frac{\partial w_{B}}{\partial\delta}\sqrt{1-e^{2}}\sin E\times
F+\frac{\partial w_{B}}{\partial\alpha*}\sqrt{1-e^{2}}\sin E\times
G\end{array}$ (10)

where $w_{B}$ is the abscissa of the barycentre, which is given by the single
star model, $A$, $B$, $F$ and $G$ are the Thiele-Innes (TI, hereafter)
elements, whose definition is given in Eq. 20, $e$ is the eccentricity of the
orbit, and $E$ is the eccentric anomaly. $E$ depends on the observation epoch,
but also on the orbital period, $P$, on $e$, and on the periastron epoch,
$T_{0}$, so the model is finally based on 12 unknowns: the 5 parameters of the
single star solution, and $A$, $B$, $F$, $G$, $P$, $e$, $T_{0}$; the partial
derivatives of $w$ with respect to the TI elements are given in Eq. 10, and
those with respect to $P$, $e$ and $T_{0}$ can be found in Goldin & Makarov
(2006).

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.

###  5.2 Pseudo-circular and circular orbital solutions

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:

$\begin{array}[]{rl}\sigma_{G}=|G|&\left[\left(\frac{\sigma_{A}}{A}\right)^{2}+\left(\frac{\sigma_{B}}{B}\right)^{2}+\left(\frac{\sigma_{F}}{F}\right)^{2}\right.\\\
&\left.+2\left(\frac{\rho_{AB}\sigma_{A}\sigma_{B}}{AB}+\frac{\rho_{AF}\sigma_{A}\sigma_{F}}{AF}+\frac{\rho_{BF}\sigma_{B}\sigma_{F}}{BF}\right)\right]^{1/2}\end{array}$
(11)

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

$\begin{array}[]{rl}\rho_{AG}=&+1\\\ \rho_{BG}=&-1\\\
\rho_{FG}=&+1\end{array}$ (12)

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

###  5.3 Selection of the orbital solutions

Figure 3: Selection of the astrometric orbital solutions. The first row refers
to the solution obtained from a processing trial without the selection based
on parallax error. (a): Period vs. mass function diagram, showing a clump of
plausible solutions, bounded by $f_{\cal M}\approx 0.3\;{\cal M}_{\sun}$, and
various sequences, organised around particular periods, which cross the clump
and extend well beyond. (b): Period vs. “parallax significance”. The solutions
with $f_{\cal M}<0.3\;{\cal M}_{\sun}$ are in green, and the large $f_{\cal
M}$, which are well clustered around the same periods as in panel (a), are in
black. The line is the limit: $\varpi/\sigma_{\varpi}\;=\;20\;000/P_{d}$ which
was implemented in the main processing to discard the dubious solutions. The
second row shows the density map diagrams used to set up the post-processing
filters to remove the concentrations still existing over specific periods.
(c): Period–eccentricity uncertainty diagram; solutions above the red line
were rejected (d): Period–significance diagram; solutions below the red line
were rejected. The third row shows properties of the DR3 solutions eventually
selected, after applying the post-processing filters. (e): the
period–eccentricity diagram; the black curve is the maximum eccentricity
assuming that the periods shorter than 10 days are circularized. The black
dots are the solutions with $f_{\cal M}>0.3\;{\cal M}_{\sun}$. The
astrometric-only orbital solutions are in green, while the solutions confirmed
by a spectroscopic SB1 orbit from Gaia are in purple. (f): histogram of
periods; the proportion of solutions with mass function above 0.3 solar
masses, in black, increases with period, as the gap between these solutions
and the totality, in green, is reduced (in order to highlight possible
overdensities, the step size of the green histogram is very fine; to avoid
fluctuations in counts from small numbers, the black histogram was plotted by
merging the count intervals into clusters of 4, and giving the mean value for
each cluster).

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:

$f_{\cal M}=\frac{{a_{0}}^{3}\times 365.25^{2}}{P^{2}\varpi^{3}}$ (13)

where $P$ is the period in days, and where $a_{0}$ and $\varpi$ are in the
same units; $f_{\cal M}$ is thus obtained in solar masses. Taking into account
that $a_{0}$ refers to the orbit of the photocentre, the third Kepler’s law
gives the following expression:

$f_{\cal M}=\frac{\|{\cal F}_{1}{\cal M}_{2}-{\cal F}_{2}{\cal
M}_{1}\|^{3}}{({\cal F}_{1}+{\cal F}_{2})^{3}({\cal M}_{1}+{\cal M}_{2})^{2}}$
(14)

where ${\cal F}_{1}$ and ${\cal F}_{2}$ are the photometric fluxes of the
brightest and of the faintest component, respectively. ${\cal M}$ refers to
the masses. When ${\cal F}_{2}$ is negligible compared to ${\cal F}_{1}$, this
equation becomes:

$f_{\cal M}({\cal F}_{2}\ll{\cal F}_{1})=\frac{{\cal M}_{1}q^{3}}{(1+q)^{2}}$
(15)

where $q={\cal M}_{2}/{\cal M}_{1}$ and where the index “1” refers to the
brightest component.

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:

$\varpi/\sigma_{\varpi}\;=\;20\;000/P_{d}$ (16)

and it allowed to keep $179\,234$ astrometric orbits as a result of the main
processing.

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.

##  6 The variability-induced movers

The variability-induced movers (VIM hereafter) were first described by Wielen
(1996). They consist in unresolved binary systems containing one
photometrically variable component. The result is a photocentre which is
moving between the components in accordance with fluctuations in the total
brightness of the system. Simultaneously, it is following the orbital motion,
if this is perceptible over the duration of the mission. Several types of VIM
were expected from the Gaia astrometric and photometric measurements: when the
orbital motion is not perceptible, the system is a fixed VIM, or VIMF. When
the orbital motion looks like a motion of the components relative to each
other, in a straight line and at constant velocity, the system is a VIM with
linear motion, or VIML. This is followed by VIM with acceleration (constant or
variable) or VIMA, and then VIM with orbital motion, or VIMO. Only VIMFs are
described in detail below, as the others could not be included in the DR3.

###  6.1 The VIM with fixed components (VIMF) model

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:

$\begin{array}[]{ll}\frac{\partial w}{\partial D_{\alpha*}}&=\frac{\partial
w}{\partial\alpha*}\;\left(\frac{\bar{\cal F}}{\cal F}-1\right)\\\ \\\
\frac{\partial w}{\partial D_{\delta}}&=\frac{\partial
w}{\partial\delta}\;\left(\frac{\bar{\cal F}}{\cal F}-1\right)\end{array}$
(17)

where ${\cal F}$ is the photometric flux of the transit. Therefore, the
position of the system and its proper motion also apply to the photocentre
when the total flux is $\bar{\cal F}$.

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:

$\sigma_{\mathrm{mod}}=\sigma_{\cal F}\;\frac{\bar{\cal F}}{{\cal
F}^{2}}\;\left|\left(\frac{\partial
w}{\partial\alpha*}D_{\alpha*}+\frac{\partial
w}{\partial\delta}D_{\delta}\right)\right|$ (18)

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

$\sigma_{\mathrm{VIMF}}=\sqrt{\sigma_{w}^{2}+\sigma_{\mathrm{mod}}^{2}}$ (19)

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}$.

###  6.2 Selection of the VIMF solutions

Figure 4: The “parallax significance” vs. $D$ diagram of VIMF solutions
obtained from a processing trial; most false solutions of $D>100$ mas are
eliminated if the selection criterion $\varpi/\sigma_{\varpi}>30$ (to the
right of the vertical red line) is applied.

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.

Figure 5: Post-processing filtering of the VIMF solutions. Panel (a): The
$F_{2}$ vs. $D$ diagram of the 2508 solutions obtained from the main
processing. The proportion of large $D$, and therefore false solutions,
increases gradually beyond $F_{2}=25$, to the right of the vertical red line.
Panel (b): The significance vs. $D$ diagram of the 1660 solutions with
$F_{2}<25$; acceptable solutions cluster on a sequence of $D$ increasing with
significance; abnormally large $D$ solutions are rejected by restricting to
significances greater than 20 (to the right of the vertical red line), leaving
a final selection of 869 solutions.

###  6.3 The other VIM models

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.

##  7 Discussion and conclusion

This first publication of astrometric binaries of different types observed by
Gaia was made under difficult conditions, since the uncertainties of the
measurements were often very underestimated. However, we were able to
calculate hundreds of thousands of reliable solutions by keeping the weights
of the observations and correcting the uncertainties. Nevertheless, this was
done at the cost of a drastic selection, based on very high quality criteria,
and it would be difficult to use the selected solutions to deduce the
statistical properties of binary stars, as the sample is incomplete and
biased.

Because of the cascade structure that leads to the testing of different
models, we may have retained solutions by direct acceptance that would have
been superseded by another model, had it been tried. This choice has the
advantage of preserving the quality of the orbital solutions, by assigning
acceleration solutions to partially resolved double stars that would have been
processed despite the filters set up to reject them. However, it was done at
the expense of the completeness of the orbital solutions, and it is likely
that some acceleration solutions actually correspond to binaries with periods
short enough to have an orbital solution. Moreover, the quality criteria alone
could not bring the acceleration solutions to a reasonable rate of high
acceleration. Ulrich Bastian studied this problem long after the post-
processing was done, and he estimated that constant acceleration solutions are
only highly reliable when the parallax is greater than 5 mas. This criterion
leaves only $15\,000$ solutions, but none with acceleration greater than 3
au.yr-2. When only quality criteria are applied, ie. when criterion (1) of
Sect. 4.3 is not applied, the restriction to parallaxes greater than 5 mas
reduces the $\Gamma>3$ au.yr-2 rate from 3.1 % to less than 1 %; this confirms
the effectiveness of this extremely drastic selection. For future DRs, it is
hoped that the treatment of partially resolved double stars, as well as
testing all astrometric models and choosing the one most compatible with the
data, will allow a wider selection of acceleration solutions.

The selection of orbital solutions was particularly severe to reject orbits
that are concentrated in certain periods, which are likely to be produced by
the satellite’s scanning law. Despite these difficulties, about $165\,500$
astrometric orbits were finally selected in DR3. The periods of the false
orbits should be different in the following data releases, as the precession
spin has been reversed during the last 6 months of observations that will be
part of the DR4; however, if these noise-induced orbits do not disappear at
this time, they could be more difficult to detect, being distributed over a
larger number of periods. Furthermore, the lack of one-year periods, which
comes from a degeneracy between parallactic and orbital motion, is expected to
remain.

The acceleration and orbital models apply to unresolved binaries with either
components that do not vary photometrically, or one component always much
brighter than the other. This limits their interest in the study of young
stars, or more generally, pulsating stars. VIM models extend the scientific
scope of astrometric binaries, provided that only one of the components is
photometrically variable. Unfortunately, the number of VIM solutions appears
to be very modest, numbering only a few hundred, limited to the VIMs with
fixed components. This rarity is partly due to the fact that the model has
only been tried at the end of the cascade, and we do not know how many
acceptable VIM solutions were lost due to the acceptance of a solution of
another type.

When both components are variable, the VIM solution can be very degraded: When
the brightnesses of the two stars vary in the same direction, the displacement
of the photocentre can be very small while the brightness of the system varies
greatly; for the VIMF model, this is interpreted as a short D distance.
Conversely, if the brightness of one component increases while that of the
other decreases, the fluctuation in total brightness that may be small while
the photocentre displacement is large. As the luminosities of the components
do not generally vary in phase, one can expect VIMF solutions of large
$\chi^{2}$, and therefore large $F_{2}$, with $D$ values that can be minimal
or, on the contrary, very large. A large $\chi^{2}$ leads to a decrease in
significance, thanks to the coefficient $c$ of Eq. 2, and, as a result, the
pollution of VIMF solutions by systems whose two components are
photometrically variable will be limited thanks to the high acceptance
thresholds that have been adopted.

Another difference between fixed-component VIMs and other solutions is that
the selection could only be adjusted on the basis of the median distance
between the photocentre and the photometric variable component, excluding any
dynamic criteria. Therefore, these solutions should be considered with
caution. In turn, it is possible that some acceleration or orbital solutions
have been assigned to objects that should have received a VIM solution.
However, these cases are too rare to constitute detectable pollution: Arenou
et al. (2022) have shown that, statistically, acceleration solutions have
proper motions affected by orbital motion, which would not be the case for
VIMFs. Similarly, if VIMFs receive an orbital solution, the inclination of the
solution must be 90 degrees, and the distribution of inclinations does not
show an excess for this value.

Despite these limitations raised above, the DR3 is a major advance in the
field of binary stars. The hundreds of thousands of new astrometric orbits
will have considerable importance in the future, compared to the less than
3000 visual binary star orbits in the Sixth Catalog of Orbits of Visual Binary
Stars111https://www.usno.navy.mil/USNO/astrometry/optical-IR-prod/wds/orb6.
Such a change of scale could open a new era in the statistical study of binary
stars, especially if the next DRs could be less affected by selection effects.
Instead of looking at the properties of pairs with a primary component of a
given type, it may be possible to consider pairs with a given total mass. Such
an approach would shed new light on the star formation process, provided that
the many biases affecting binary selection are taken into account.

###### Acknowledgements.

This work is part of the reduction of the Gaia satellite observations
(https://www.cosmos.esa.int/gaia). The Gaia space mission is operated by the
European Space Agency, and the data are being processed by the Gaia Data
Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium)). The Gaia archive
website is https://archives.esac.esa.int/gaia. Funding for the DPAC is
provided by national institutions, in particular the institutions
participating in the Gaia MultiLateral Agreement (MLA). We acknowledge for the
financial support of the french “Centre National d’études spatiales” (CNES)
and of the BELgian federal Science Policy Office (BELSPO) through a PROgramme
de Développement d’Expériences scientifiques (PRODEX). Over the two decades
that it took to develop the calculation methods and write the software, we
benefited from the contributions of Sylvie Jancart, Rémy Onsay, Didier Pelat
and Patrick Gempin. We also received constant support from François Mignard
and Alain Jorissen. The acceptability of acceleration solutions was discussed
with Ulrich Bastian and Lennart Lindegren, using data produced by Sergei
Klioner and Hagen Steidelmüller. Finally, we are grateful to Michael Biermann
and to an anonymous referee for their careful reading of the manuscript and
the relevance of their comments.

## References

  * Arenou et al. (2022) Arenou, F., Babusiaux, C., Barstow, M. et al. 2022 A&A, accepted
  * Babusiaux et al. (2022) Babusiaux, C. et al. 2022 A&A, in preparation
  * Binnendijk (1960) Binnendijk, L. 1960, Properties of double stars; a survey of parallaxes and orbits, University of Pennsylvania Press, Philadelphia
  * Dravins, Lindegren & Madsen (1999) Dravins, D., Lindegren, L., & Madsen, S. 1999, A&A, 348, 1040
  * ESA (1997) ESA 1997, The Hipparcos and Tycho Catalogue, SP-1200
  * Gaia Collaboration (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., Prusti, T., de Bruijne, J. H. J., Babusiaux, C. & M. Biermann, M. 2021, A&A, 649, A1
  * Goldin & Makarov (2006) Goldin, A. & Makarov, V.V. 2006, ApJS, 166, 341
  * Gosset et al. (2022) Gosset, E. et al. 2022 A&A, in preparation
  * Halbwachs, Mayor & Udry (2005) Halbwachs, J.-L., Mayor, M. & Udry, S. 2005, A&A, 431, 1129
  * Halbwachs (2009) Halbwachs, J.-L. 2009, MNRAS394, 1075
  * Heintz (1978) Heintz,W.D. 1978, Double stars, Reidel, Dordrecht
  * Holl et al. (2022) Holl, B., Sozzetti, A., Sahlmann, J. et al. 2022 A&A, in preparation
  * Lindegren & Dravins (2003) Lindegren, L. & Dravins, D. 2003, A&A, 401, 1185
  * Lindegren et al. (2021) Lindegren, L., Klioner, S.A., Hernández, J., et al. 2021, A&A, 649, A2
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3
  * Seabroke et al. (2021) Seabroke, G., Fabricius, C., Teyssier, D., et al. 2021, A&A, submitted
  * Stuart & Ord (1994) Stuart A., Ord K., 1994, Kendall’s Advanced Theory of Statistics, vol. 1. Edward Arnold, London
  * Wielen (1996) Wielen, R. 1996, A&A, 314, 679

##  Appendix A Conversion of Thiele-Innes elements ($A,B,F,G$) into Campbell
elements ($a,i,\Omega,\omega$) and calculation of the uncertainties

The Thiele-Innes (TI) elements, $A$, $B$, $F$ and $G$ are given by the
following equations (Binnendijk, 1960; Heintz, 1978):

$\begin{array}[]{ll}A=&a\;(\cos\omega\cos\Omega-\sin\omega\sin\Omega\cos i)\\\
B=&a\;(\cos\omega\sin\Omega+\sin\omega\cos\Omega\cos i)\\\
F=&-a\;(\sin\omega\cos\Omega+\cos\omega\sin\Omega\cos i)\\\
G=&-a\;(\sin\omega\sin\Omega-\cos\omega\cos\Omega\cos i)\end{array}$ (20)

where $a$, $i$, $\Omega$ and $\omega$ are the Campbell elements. The
definitions of these elements and their calculations from $A$, $B$, $F$, $G$
are given hereafter:

  * •

$a$ (or $a_{0}$ for a photocentric orbit) is the semi-major axis of the orbit,
in the same unit as the TI elements.

  * •

$i$ is the inclination of the orbit, i.e. the angle between the orbital plane
and the plane of the sky. $i\in[0,\pi/2]$ when the orbital motion is in the
direct sense, and $i\in[\pi/2,\pi]$ when the photocentre is revolving around
the barycentre in a clockwise direction.

  * •

$\Omega$ is the position angle of the ascending node. For us, the ascending
node is the intersection between the orbital plane and the plane of the sky in
the direction where the right ascension is increasing. Therefore,
$\Omega\in[0,\pi]$.

  * •

$\omega$ is the periastron longitude, measured in the orbital plane from the
ascending node defined above. $\omega$ is measured in the sense ofthe orbital
motion.

Binnendijk (1960) gives the following method for converting TI elements into
Campbell elements:

The semi-major axis $a$, is derived from:

$\begin{array}[]{ll}u=&(A^{2}+B^{2}+F^{2}+G^{2})/2\\\ v=&AG-BF\\\
a=&\sqrt{u+\sqrt{(u+v)(u-v)}}\end{array}$ (21)

The angles $\Omega$ and $\omega$ are derived simultaneously. Preliminary
estimates, between $-\pi$ and $\pi$, are derived as follows:

$\begin{array}[]{ll}\omega+\Omega=&\arctan\frac{B-F}{A+G}\pmod{\pi}\\\ &\\\
\omega-\Omega=&\arctan\frac{B+F}{G-A}\pmod{\pi}\end{array}$ (22)

The exact values of $\omega+\Omega$ are obtained, modulo $2\pi$, taking into
account that $\sin(\omega+\Omega)$ and $(B-F)$ have the same sign. Similarly,
$\sin(\omega-\Omega)$ has the same sign as $(-B-F)$. $\omega$ and $\Omega$ are
then derived from the equations:

$\begin{array}[]{ll}\omega=&\frac{(\omega+\Omega)+(\omega-\Omega)}{2}\pmod{\pi}\\\
&\\\ \Omega=&\frac{(\omega+\Omega)-(\omega-\Omega)}{2}\pmod{\pi}\end{array}$
(23)

When $\Omega$ is negative, $\pi$ must still be added to both $\Omega$ and
$\omega$. After that, a correction of $2\pi$ may still be applied in order to
meet the conditions $\Omega\in[0,\pi]$ and $\omega\in[0,2\pi]$.

To derive the inclination $i$, we first need to calculate two intermediate
terms, called $d_{1}$ and $d_{2}$ hereafter:

$\begin{array}[]{ll}d_{1}=&\|(A+G)\cos(\omega-\Omega)\|\\\
d_{2}=&\|(F-B)\sin(\omega-\Omega)\|\end{array}$ (24)

In order to obtain the value of $i$, the calculation is different depending on
whether $d_{1}$ is greater than $d_{2}$, or the reverse:

$\begin{array}[]{ll}\mathrm{if}\;(d_{1}\geq
d_{2}),&i=2\arctan\sqrt{\|(A-G)\cos(\omega+\Omega)\|\;/\;d_{1}}\\\
\mathrm{else,}&i=2\arctan\sqrt{\|(B+F)\sin(\omega+\Omega)\|\;/\;d_{2}}\end{array}$
(25)

The uncertainties of $a$, $i$, $\Omega$ and $\omega$ are derived by
differentiating the equations above, and by taking into account error
propagation.

We call hereafter $\sigma_{p}$ the error on element $p$ and
$\mathrm{cov}$(X,Y) the term of the variance covariance matrix of the solution
linking the elements $X$ and $Y$, ie
$\mathrm{cov}(X,Y)=\rho_{XY}\sigma_{X}\sigma_{Y}$, where $\rho_{XY}$ is the
correlation coefficient between $X$ and $Y$. The uncertainties are then given
by the equations below.

To derive the uncertainty of $a$, the following intermediate terms must be
calculated from $u$ and $v$ derived in Eq. 21:

$\begin{array}[]{ll}t_{A}=&A\;+\;(A\;u-G\;v)\;/\sqrt{u^{2}-v^{2}}\\\
t_{B}=&B\;+\;(B\;u+F\;v)\;/\sqrt{u^{2}-v^{2}}\\\
t_{F}=&F\;+\;(F\;u+B\;v)\;/\sqrt{u^{2}-v^{2}}\\\
t_{G}=&G\;+\;(G\;u-A\;v)\;/\sqrt{u^{2}-v^{2}}\end{array}$ (26)

The uncertainty of the semi-major axis is then:

$\begin{array}[]{ll}\sigma_{a}=\frac{1}{2a}\;\times\;[&t_{A}^{2}\sigma_{A}^{2}+t_{B}^{2}\sigma_{B}^{2}+t_{F}^{2}\sigma_{F}^{2}+t_{G}^{2}\sigma_{G}^{2}\\\
&+\;2t_{A}t_{B}\mathrm{cov}(A,B)+2t_{A}t_{F}\mathrm{cov}(A,F)\\\
&+\;2t_{A}t_{G}\mathrm{cov}(A,G)+2t_{B}t_{F}\mathrm{cov}(B,F)\\\
&+\;2t_{B}t_{G}\mathrm{cov}(B,G)+2t_{F}t_{G}\mathrm{cov}(F,G)\;]^{1/2}\end{array}$
(27)

The uncertainties of $\omega$ and $\Omega$ are also derived using dedicated
intermediate terms. We first calculate the followings:

$\begin{array}[]{ll}k=&(A+G)^{2}+(B-F)^{2}\\\
l=&(G-A)^{2}+(B+F)^{2}\end{array}$ (28)

These terms are used to derive other intermediate terms dedicated to the
calculation of the uncertainty of $\omega$:

$\begin{array}[]{ll}w_{A}=&(F-B)/k+(B+F)/l\\\ w_{B}=&(A+G)/k+(G-A)/l\\\
w_{F}=&-(A+G)/k+(G-A)/l\\\ w_{G}=&(F-B)/k-(B+F)/l\end{array}$ (29)

The uncertainty of $\omega$, $\sigma_{\omega}$ is then derived with the
equation hereafter:

$\begin{array}[]{ll}\sigma_{\omega}=&[w_{A}^{2}\sigma_{A}^{2}+w_{B}^{2}\sigma_{B}^{2}+w_{F}^{2}\sigma_{F}^{2}+w_{G}^{2}\sigma_{G}^{2}\\\
&+\;2\;(\;w_{A}w_{B}\;\mathrm{cov}(A,B)+w_{A}w_{F}\;\mathrm{cov}(A,F)\\\
&+w_{A}w_{G}\;\mathrm{cov}(A,G)+w_{B}w_{F}\;\mathrm{cov}(B,F)\\\
&+w_{B}w_{G}\;\mathrm{cov}(B,G)+w_{F}w_{G}\;\mathrm{cov}(F,G)\;)\;]^{1/2}\;/\;2\end{array}$
(30)

The terms from Eq. 28 are also used to compute the following intermediate
terms, which are required to derive the uncertainty of $\Omega$:

$\begin{array}[]{ll}O_{A}=&(F-B)/k-(B+F)/l\\\ O_{B}=&(A+G)/k-(G-A)/l\\\
O_{F}=&-(A+G)/k-(G-A)/l\\\ O_{G}=&(F-B)/k+(B+F)/l\end{array}$ (31)

The uncertainty of $\Omega$ is then:

$\begin{array}[]{ll}\sigma_{\Omega}=&[O_{A}^{2}\sigma_{A}^{2}+O_{B}^{2}\sigma_{B}^{2}+O_{F}^{2}\sigma_{F}^{2}+O_{G}^{2}\sigma_{G}^{2}\\\
&+\;2\;(\;O_{A}O_{B}\;\mathrm{cov}(A,B)+O_{A}O_{F}\;\mathrm{cov}(A,F)\\\
&+O_{A}O_{G}\;\mathrm{cov}(A,G)+O_{B}O_{F}\;\mathrm{cov}(B,F)\\\
&+O_{B}O_{G}\;\mathrm{cov}(B,G)+O_{F}O_{G}\;\mathrm{cov}(F,G)\;)\;]^{1/2}\;/\;2\end{array}$
(32)

Like the derivation of $i$, the derivation of the uncertainty of the
inclination depends on a comparison between $d_{1}$ and $d_{2}$, derived in
Eq. 24 above.

$\sigma_{i}$ is derived as follows. We first calculate the terms:

$\begin{array}[]{ll}q=&\sin\Omega\cos\Omega\\\ r=&\sin\omega\cos\omega\\\
\end{array}$ (33)

if $d_{1}\geq d_{2}$, the continuation of the calculation is as explained
hereafter. An extra term, $p_{1}$, is added to $q$ and $r$:

$p_{1}=\cos(\omega-\Omega)\;\cos(\omega+\Omega)$ (34)

The intermediate terms dedicated to the calculation of $\sigma_{i}$ are then:

$\begin{array}[]{ll}h_{A}=&2Gp_{1}+(G^{2}-A^{2})\;(qw_{A}+rO_{A})\\\
h_{B}=&(G^{2}-A^{2})\;(qw_{B}+rO_{B})\\\
h_{F}=&(G^{2}-A^{2})\;(qw_{F}+rO_{F})\\\
h_{G}=&-2Ap_{1}+(G^{2}-A^{2})\;(qw_{G}+rO_{G})\end{array}$ (35)

and, still for $d_{1}\geq d_{2}$, $\sigma_{i}$ is derived from the equation:

$\begin{array}[]{ll}\sigma_{i}=&[h_{A}^{2}\sigma_{A}^{2}+h_{B}^{2}\sigma_{B}^{2}+h_{F}^{2}\sigma_{F}^{2}+h_{G}^{2}\sigma_{G}^{2}\\\
&+\;2\;(\;h_{A}h_{B}\;\mathrm{cov}(A,B)+h_{A}h_{F}\;\mathrm{cov}(A,F)\\\
&+h_{A}h_{G}\;\mathrm{cov}(A,G)+h_{B}h_{F}\;\mathrm{cov}(B,F)\\\
&+h_{B}h_{G}\;\mathrm{cov}(B,G)+h_{F}h_{G}\;\mathrm{cov}(F,G)\;)\;]^{1/2}\\\
&/\;[\;(\tan(i/2)+\tan^{3}(i/2))\;d_{1}^{2}\;]\end{array}$ (36)

When $d_{2}<d_{1}$, the additional term $p_{2}$ is derived:

$p_{2}=\sin(\omega-\Omega)\;\sin(\omega+\Omega)$ (37)

The intermediate terms dedicated the calculation of $\sigma_{i}$ are now:

$\begin{array}[]{ll}g_{A}=&(B^{2}-F^{2})\;(qw_{A}-rO_{A})\\\
g_{B}=&2Fp_{2}+(B^{2}-F^{2})\;(qw_{B}-rO_{B}\\\
g_{F}=&-2Bp_{2}+(B^{2}-F^{2})\;(qw_{F}-rO_{F})\\\
g_{G}=&(B^{2}-F^{2})\;(qw_{G}-rO_{G})\end{array}$ (38)

and, still for $d_{1}<d_{2}$, $\sigma_{i}$ is derived from the equation:

$\begin{array}[]{ll}\sigma_{i}=&[g_{A}^{2}\sigma_{A}^{2}+g_{B}^{2}\sigma_{B}^{2}+g_{F}^{2}\sigma_{F}^{2}+g_{G}^{2}\sigma_{G}^{2}\\\
&+\;2\;(\;g_{A}g_{B}\;\mathrm{cov}(A,B)+g_{A}g_{F}\;\mathrm{cov}(A,F)\\\
&+g_{A}g_{G}\;\mathrm{cov}(A,G)+g_{B}g_{F}\;\mathrm{cov}(B,F)\\\
&+g_{B}g_{G}\;\mathrm{cov}(B,G)+g_{F}g_{G}\;\mathrm{cov}(F,G)\;)\;]^{1/2}\\\
&/\;[\;(\tan(i/2)+\tan^{3}(i/2))\;d_{2}^{2}\;]\end{array}$ (39)

The four Campbell elements and their uncertainties are thus deduced from the
solutions calculated in TI elements. However, it is not possible to estimate
the correlation coefficients between the terms.

##  Appendix B Conversion of Thiele-Innes elements ($C_{1},H_{1}$) into
elements (${\rm a}_{1},\omega_{1}$) and calculation of the uncertainties

The $C_{1}$ and $H_{1}$ elements are functions of the ${\rm a}_{1}$ and
$\omega_{1}$ elements that come from the spectroscopic orbit:

$\begin{array}[]{ll}C_{1}=&{\rm a}_{1}\sin i\sin\omega_{1}\\\ H_{1}=&{\rm
a}_{1}\sin i\cos\omega_{1}\end{array}$ (40)

from which it appears that:

$\begin{array}[]{ll}{\rm a}_{1}=&\sqrt{C_{1}^{2}+H_{1}^{2}}/\sin i\\\
\omega_{1}=&\arctan C_{1}/H_{1}\end{array}$ (41)

knowing that sign$(\sin\omega_{1})={\rm sign}(C_{1})\textrm{ and
sign}(\cos\omega_{1})={\rm sign}(H_{1})$.

$\omega_{1}$ must be approximately equal to either $\omega$, or $\omega+\pi$.
This ambiguity is due to differences in definition: $\omega$ is the longitude
of the periastron of the photocentre orbit, which, in the absence of radial
velocity, is measured from the node whose position angle is between 0 and
$\pi$. For its part, $\omega_{1}$ is the longitude of the periastron of the
orbit of the spectroscopically most visible component, measured from the node
where this component moves away from the Sun. Therefore, the true orientation
of the orbital plane will depend on the position of the photocentre with
respect to to the barycentre and the spectroscopic component: if the
photocentre is assumed to be between the barycentre and the spectroscopic
component, then the true ascending node and periastron of the astrometric
orbit are marked by the same angles as those of the spectroscopic orbit. If
$\omega$ and $\omega_{1}$ are significantly different, it is necessary to add
$\pm\pi$ to $\omega$ to make them roughly coincide, and also to add $\pi$ to
$\Omega$ to take into account the change of reference node.

Alternatively, the photocentre may be opposite to the spectroscopic component
(this may occur, for example, when the spectral type of the brightest
component is very different from the templates used for spectroscopic
reduction, so that the spectroscopic orbit corresponds to the least bright
binary component). When this occurs, $\omega$ must be approximately equal to
$\omega_{1}\pm\pi$. If it is not, $\pm\pi$ must be added to $\omega$ and $\pi$
must be added to $\Omega$ for the orientation of the astrometric orbit to be
completely defined.

The calculation of the uncertainty of ${\rm a}_{1}$ depends on the values of
$d_{1}$ and $d_{2}$, given by Eq. 24. When $d_{1}\geq d_{2}$, we compute the
following terms:

$\left(\begin{array}[]{l}q_{A}\\\ q_{B}\\\ q_{F}\\\
q_{G}\end{array}\right)=\left(\begin{array}[]{l}h_{A}\\\ h_{B}\\\ h_{F}\\\
h_{G}\end{array}\right)\times\frac{{\rm a}_{1}}{\tan
i\tan\frac{i}{2}\left(1+\tan^{2}\frac{i}{2}\right)d_{1}^{2}}$ (42)

where $h_{A}$, $h_{B}$, $h_{F}$, and $h_{G}$ are taken from Eq. 35, and

$\left(\begin{array}[]{l}q_{C}\\\
q_{H}\end{array}\right)=\left(\begin{array}[]{l}C_{1}\\\
H_{1}\end{array}\right)\times\frac{{\rm a}_{1}}{C_{1}^{2}+H_{1}^{2}}$ (43)

The uncertainty $\sigma_{{\rm a}_{1}}$ is then calculated with equation:

$\begin{array}[]{ll}\sigma_{{\rm
a}_{1}}=&[q_{A}^{2}\sigma_{A}^{2}+q_{B}^{2}\sigma_{B}^{2}+q_{F}^{2}\sigma_{F}^{2}+q_{G}^{2}\sigma_{G}^{2}+q_{C}^{2}\sigma_{C_{1}}^{2}+q_{H}^{2}\sigma_{H_{1}}^{2}\\\
&+\;2\;(\;q_{A}q_{B}\;\mathrm{cov}(A,B)+q_{A}q_{F}\;\mathrm{cov}(A,F)\\\
&+q_{A}q_{G}\;\mathrm{cov}(A,G)+q_{A}q_{C}\;\mathrm{cov}(A,C_{1})\\\
&+q_{A}q_{H}\;\mathrm{cov}(A,H_{1})+q_{B}q_{F}\;\mathrm{cov}(B,F)\\\
&+q_{B}q_{G}\;\mathrm{cov}(B,G)+q_{B}q_{C}\;\mathrm{cov}(B,C_{1})\\\
&+q_{B}q_{H}\;\mathrm{cov}(B,H_{1})+q_{F}q_{G}\;\mathrm{cov}(F,G)\\\
&+q_{F}q_{C}\;\mathrm{cov}(F,C_{1})+q_{F}q_{H}\;\mathrm{cov}(F,H_{1})\\\
&+q_{G}q_{C}\;\mathrm{cov}(G,C_{1})+q_{G}q_{H}\;\mathrm{cov}(G,H_{1})\\\
&+q_{C}q_{H}\;\mathrm{cov}(C_{1},H_{1})\;)\;]^{1/2}\end{array}$ (44)

When $d_{2}>d_{1}$, these equations become:

$\left(\begin{array}[]{l}r_{A}\\\ r_{B}\\\ r_{F}\\\
r_{G}\end{array}\right)=\left(\begin{array}[]{l}g_{A}\\\ g_{B}\\\ g_{F}\\\
g_{G}\end{array}\right)\times\frac{{\rm a}_{1}}{\tan
i\tan\frac{i}{2}\left(1+\tan^{2}\frac{i}{2}\right)d_{2}^{2}}$ (45)

where $g_{A}$, $g_{B}$, $g_{F}$, and $g_{G}$ are taken from Eq. 38, and

$\left(\begin{array}[]{l}r_{C}\\\
r_{H}\end{array}\right)=\left(\begin{array}[]{l}C_{1}\\\
H_{1}\end{array}\right)\times\frac{{\rm a}_{1}}{C_{1}^{2}+H_{1}^{2}}$ (46)

The uncertainty $\sigma_{{\rm a}_{1}}$ is then calculated with the equation:

$\begin{array}[]{ll}\sigma_{{\rm
a}_{1}}=&[r_{A}^{2}\sigma_{A}^{2}+r_{B}^{2}\sigma_{B}^{2}+r_{F}^{2}\sigma_{F}^{2}+r_{G}^{2}\sigma_{G}^{2}+r_{C}^{2}\sigma_{C_{1}}^{2}+r_{H}^{2}\sigma_{H_{1}}^{2}\\\
&+\;2\;(\;r_{A}r_{B}\;\mathrm{cov}(A,B)+r_{A}r_{F}\;\mathrm{cov}(A,F)\\\
&+r_{A}r_{G}\;\mathrm{cov}(A,G)+r_{A}r_{C}\;\mathrm{cov}(A,C_{1})\\\
&+r_{A}r_{H}\;\mathrm{cov}(A,H_{1})+r_{B}r_{F}\;\mathrm{cov}(B,F)\\\
&+r_{B}r_{G}\;\mathrm{cov}(B,G)+r_{B}r_{C}\;\mathrm{cov}(B,C_{1})\\\
&+r_{B}r_{H}\;\mathrm{cov}(B,H_{1})+r_{F}r_{G}\;\mathrm{cov}(F,G)\\\
&+r_{F}r_{C}\;\mathrm{cov}(F,C_{1})+r_{F}r_{H}\;\mathrm{cov}(F,H_{1})\\\
&+r_{G}r_{C}\;\mathrm{cov}(G,C_{1})+r_{G}r_{H}\;\mathrm{cov}(G,H_{1})\\\
&+r_{C}r_{H}\;\mathrm{cov}(C_{1},H_{1})\;)\;]^{1/2}\end{array}$ (47)

The uncertainty of $\omega_{1}$, $\sigma_{\omega_{1}}$, is derived from:

$\sigma_{\omega_{1}}=\frac{\cos^{2}\omega_{1}}{H_{1}^{2}}[H_{1}^{2}\sigma_{C_{1}}^{2}+C_{1}^{2}\sigma_{H_{1}}^{2}-2C_{1}H_{1}\;\mathrm{cov}(C_{1},H_{1})]^{1/2}$
(48)


00footnotetext: Version 2015-07-03_13. See
https://shelah.logic.at/papers/1061/ for possible updates.

# On failure of 0-1 laws

Saharon Shelah  Einstein Institute of Mathematics
Edmond J. Safra Campus, Givat Ram
The Hebrew University of Jerusalem
Jerusalem, 91904, Israel
and
Department of Mathematics
Hill Center - Busch Campus
Rutgers, The State University of New Jersey
110 Frelinghuysen Road
Piscataway, NJ 08854-8019 USA<EMAIL_ADDRESS>http://shelah.logic.at

(Date: July 3, 2015)

###### Abstract.

Let $\alpha\in(0,1)_{{\mathbb{R}}}$ be irrational and
$G_{n}=G_{n,1/n^{\alpha}}$ be the random graph with edge probability
$1/n^{\alpha}$; we know that it satisfies the 0-1 law for first order logic.
We deal with the failure of the 0-1 law for stronger logics:
${\mathbb{L}}_{\infty,{\bf k}},{\bf k}$ a large enough natural number and the
inductive logic.

###### Key words and phrases:

finite model theory, zero-one laws, random graphs, inductive logic, infinitary
logic on finite structures

###### 2010 Mathematics Subject Classification:

Primary: 03C13, 05C80; Secondary: 03C10, 03C80

This work was partially supported by European Research Council grant 338821.
Publication 1061 on Shelah’s list. The author thanks Alice Leonhardt for the
beautiful typing.

_Dedicated to Yuri Gurevich on the occasion of his 75th birthday_

Let $G_{n,p}$ be the random graph with set of nodes $[n]=\\{1,\dotsc,n\\}$,
each edge of probability $p\in[0,1]_{{\mathbb{R}}}$, the edges being drawn
independently, (see $\boxplus_{1}$ below). On 0-1 laws (and random graphs) see
the book of Spencer [6] or Alon-Spencer [1], in particular on the behaviour of
the random graph $G_{n,1/n^{\alpha}}$ for $\alpha\in(0,1)_{{\mathbb{R}}}$
irrational. On finite model theory see Flum-Ebbinghaus [2], e.g. on the logic
${\mathbb{L}}_{\infty,\mathbf{k}}$ and on inductive logic, also called LFP
logic (i.e. least fix point logic). A characteristic example of what can be
expressed in this logic is “in the graph $G$ there is a path from the node $x$
to the node $y$”, this is closed to what we shall use. We know that $G_{n,p}$
(i.e. the case the probability $p$ is constant), satisfies the 0-1 law for
first order logic (proved independently by Fagin [3] and Glebskii-et-al [4]).
This holds also for many stronger logics like ${\mathbb{L}}_{\infty,{\bf k}}$
and the inductive logic. If $\alpha\in(0,1)_{{\mathbb{R}}}$ is irrational, the
0-1 law holds for $G_{n,(1/n^{\alpha})}$ and first order logic.

The question we address is whether this holds also for stronger logics as
above. Though our real aim is to address the problem for the case of graphs,
the proof seems more transparent when we have two random graph relations (with
appropriate probabilities; we make them directed graphs just for simplicity).
So we shall deal with two cases A and B. In Case A, the usual graph, we have
to show that there are (just first order) formulas $\varphi_{\ell}(x,y)$ for
$\ell=1,2$ with some special properties, (actually we have also
$\varphi_{0}(x,y)$). For Case B, those formulas are $R_{\ell}(x,y),\ell=1,2$,
the two directed graph relations. Note that (for Case B), the satisfaction of
the cases of the $R_{\ell}$ are decided directly by the drawing and so are
independent, whereas for Case A there are (small) dependencies for different
pairs, so the probability estimates are more complicated.

Recall

  1. $\boxplus_{1}$

a 0-1 context consists of:

    1. $(a)$

a vocabulary $\tau$, here just the one of graphs or double directed graphs,

    2. $(b)$

for each $n,K_{n}$ is a set of $\tau$-models with set of elements = nods
$[n]$, in our case graphs or double directed graphs,

    3. $(c)$

a distribution $\mu_{n}$ on $K_{n}$, i.e.
$\mu_{n}:K\rightarrow[0,1]_{{\mathbb{R}}}$ satisfying
$\Sigma\\{\mu_{n}(G):G\in K_{n}\\}=1$

    4. $(d)$

the random structure is called $G_{n}=G_{\mu_{n}}$ and we tend to speak on
$G_{\mu_{n}}$ or $G_{n}$ rather than on the context.

Note that in this work “for every random enough $G_{n}$ …” is a central
notion, where:

  1. $\boxplus_{2}$

for a given 0-1 context, let “for every random enough $G_{n}$ we have
$G_{n}\models\psi$, i.e. $G$ satisfies $\psi$” means that the sequence
$\langle{\rm Prob}(G_{n}\models\psi):n\in{\mathbb{N}}\rangle$ converge to 1;
of course, ${\rm Prob}(G_{n}\models\psi)=\Sigma\\{\mu_{n}(G):G\in K_{n}$ and
$G\models\psi\\}$.

But

  1. $\boxplus_{3}$

$G_{n,p}$ is the case $K_{n}={\rm graph}$ on $[n]$ and we draw the edges
independently,

    1. (a)

with probability $p$ when $p$ is constant, e.g. $\frac{1}{2}$, and

    2. (b)

with probability $p(n)$ or probability $p_{n}$ when $p$ is a function from
${\mathbb{N}}$ to $[0,1]_{{\mathbb{R}}}$.

In the constant $p$ case, the 0-1 law is strong: it is done by proving
elimination of quantifiers and it works also for stronger logics:
${\mathbb{L}}_{\infty,\mathbf{k}}$ and so also for inductive logic
${\mathbb{L}}_{{\rm ind}}$. Another worthwhile case is:

  1. $\boxplus_{4}$

$G_{n,1/n^{\alpha}}$ where $\alpha\in(0,1)_{{\mathbb{R}}}$; so
$p_{n}=1/n^{\alpha}$.

Again the edges are drawn independently but the probability depends on $n$.

The 0-1 law holds if $\alpha$ is irrational, but we have elimination of
quantifiers only up to (Boolean combination of) existential formulas. Do we
have 0-1 law also for those stronger logics? We shall show that not by proving
that for some so called scheme $\bar{\varphi}$ of interpretation, for any
random enough $G_{n}$, $\bar{\varphi}$ interpret an initial segment of number
theory, say up to $m(G_{n})$ where $m(G_{n})$ is not too small; e.g. at least
$\log_{2}(\log_{2}(n))$.

For the probabilistic argument we use estimates; they are as in the first
order case (see [1], so we do not repeat them).

For the full version see the author website or the mathematical arXive. The
statements for which we need more estimates will probably be further delayed;
those are the ones proving that:

  1. $\boxplus_{5}$

    1. $\bullet$

using $n^{\varepsilon}$ instead of $\log_{2}(\log_{2}(n))$ in the proof for
Case 1 so the value of “${\rm Prob}(G_{n,1/n^{\alpha}})=\psi$” may change more
quickly,

    2. $\bullet$

we can define “$n$ even” (i.e. ${\rm Lim}({\rm
Prob}(G_{n,1/n^{\alpha}}\models\psi$ iff $n$ is even) exists and is one; this
is done by defining a linear order on $G_{n,\bar{\alpha}}$.

    3. $\bullet$

we may formalize the quantification on paths, so getting a weak logic failing
the 0-1 law, but its naturality is not so clear.

A somewhat related problem asks whether for some logic the 0-1 law holds for
$G_{n,p}$ (for constant $p\in(0,1)_{\mathbb{R}}$, e.g. $p=\frac{1}{2}$) but
does not have the elimination of quantifier, see [5].

We now try to informally describe the proof, naturally concentrating on case
B.

Fix reals $\alpha_{1}<\alpha_{2}$ from $(0,\frac{1}{4})_{{\mathbb{R}}}$, so
$\bar{\alpha}=(\alpha_{1},\alpha_{2})$ letting $\alpha(\ell)=\alpha_{\ell}$;

  1. $\boxplus_{6}$

let the random digraph
$G_{n,\bar{\alpha}}=([n],R_{1},R_{2})=([n],R^{G_{n,\bar{\alpha}}}_{1},R^{G_{n,\bar{\alpha}}}_{2})$
with $R_{1},R_{2}$ irreflexive relations drawn as follows:

    1. $(a)$

for each $a\neq b$, we draw a truth value for $R_{2}(a,b)$ with probability
$\frac{1}{n^{1-\alpha_{2}}}$ for yes

    2. $(b)$

for each $a\neq b$, we draw a truth value for $R_{1}(a,b)$ with probability
$\frac{1}{n^{1+\alpha_{1}}}$ for yes

    3. $(c)$

those drawings are independent.

Now for random enough digraph $G=G_{n}=G_{n,\bar{\alpha}}=([n],R_{1},R_{2})$
and node $a\in G$ we try to define the set $S_{k}=S_{G,a,k}$ of nodes of $G$
not from $\cup\\{S_{m}:m<k\\}$ by induction on $k$ as follows:

For $k=0$ let $S_{k}=\\{a\\}$. Assume $S_{0},\dots,S_{k}$ has been chosen, and
we shall choose $S_{k+1}$.

  1. $\boxplus_{7}$

For $\iota=1,2$ we ask: is there an $R_{\iota}$-edge $(a,b)$ with $a\in S_{k}$
and $b$ not from $\cup\\{S_{m}:m\leq k\\}$?

If the answer is no for both $\iota=1,2$ we stop and let ${0pt}(a,G)=k$. If
the answer is yes for $\iota=1$, we let $S_{k+1}$ be the set of $b$ such that
for some $a$ the pair $(a,b)$ is as above for $\iota=1$., If the answer is no
for $\iota=1$ but yes for $\iota=2$ we define $S_{k+1}$ similarly using
$\iota=2$.

Let the height of $G$ be $\max\\{0pt(a,G):a\in G\\}$. Now we can prove that
for every random enough $G_{n}$, for $a\in G_{n}$ or easier- for most $a\in
G_{n}$, for not too large $k$ we have:

  1. $\boxplus_{8}$

$S_{G_{n},a,k}$ is on the one hand not empty and on the other hand with $\leq
n^{2\alpha_{2}}$ members.

This is proved by drawing the edges not all at once but in $k$ stages. In
stage $m\leq k$ we already can compute $S_{G_{n},a,0},\dots S_{G_{n},a,m}$ and
we have already drawn all the $R_{1}$-edges and $R_{2}$-edges having at least
one node in $S_{G_{n},a,0}\cup\dots\cup S_{G_{n},a,m-1}$; that is for every
such pair $(a,b)$ we draw the truth values of $R_{1}(a,b),R_{2}(a,b)$. So
arriving to $m$ we can draw the edges having a nod in $S_{m}$ and not dealt
with earlier, and hence can compute $S_{m+1}$.

The point is that in the question $\boxplus_{7}$ above, if the answer is yes
for $\iota=1$ then the number of nodes in ${S_{m+1}}$ will be small,
essentially smaller than in $S_{m}$. Further, if the answer for $\iota=1$ the
answer is no but for $\iota=2$ the answer is yes then necessarily $S_{m}$ is
smaller than say $n^{(\alpha_{1}+\alpha_{2})/2}$ but it is known that the
$R_{2}$-valency of of any nod of $G_{n}$ is near $n^{\alpha_{2}}$. So the
desired inequality holds.

By a similar argument, if we stop at $k$ then in $S_{0}\cup\dots\cup S_{k}$
there are many nodes-e.g. at least near $n^{\alpha_{2}}$ by a crud argument.
As each $S_{m}$ is not too large necessarily the height of $G_{n}$ is large.

The next step is to express in our logic the relation
$\\{(a_{1},b_{1},a_{2},b_{2}):$ for some $k_{1},k_{2}$ we have $b_{1}\in
S_{G_{n},a_{1},k_{1}},b_{2}\in S_{G_{n},a_{2},k_{2}},k_{1}\leq k_{2}\\}$.

By this we can interpret a linear order with $0pt(G_{n})$ members. Again using
the relevant logic this suffice to interpret number theory up to this height.
Working more we can define a linear order with $n$ elements, so can
essentially find a formula “saying” $n$ is even (or odd).

For random graphs we have to work harder: instead of having two relations we
have two formulas; one of the complications is that their satisfaction for the
relevant pairs are not fully independent.

## References

  * [1] Noga Alon and Joel H. Spencer.  The probabilistic method.  Wiley-Interscience Series in Discrete Mathematics and Optimization. John Wiley & Sons, Inc., Hoboken, NJ, third edition, 2008.  With an appendix on the life and work of Paul Erdős.
  * [2] Heinz-Dieter Ebbinghaus and Jörg Flum.  Finite model theory.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, enlarged edition, 2006.
  * [3] Ronald Fagin.  Probabilities in finite models.  Journal of Symbolic Logic, 45:129–141, 1976.
  * [4] Y.V. Glebskii, D.I. Kogan, M.I. Liagonkii, and V.A. Talanov.  Range and degree of reliability of formulas in restricted predicate calculus.  Kibernetica, 5:17–27, 1969.  translation of Cybernetics vol 5 pp 142-154.
  * [5] Saharon Shelah.  Random Graphs: stronger logic but with the 0-1 law. In preparation.
  * [6] Joel Spencer.  The strange logic of random graphs, volume 22 of Algorithms and Combinatorics.  Springer-Verlag, Berlin, 2001.