© The Authors 2025

1 Introduction

Finding and characterizing exoplanets has become one of the most active areas in astronomy. So far, most exoplanets have been found by the transit and the radial velocity (RV¹) techniques, as seen in the few publicly available exoplanet catalogs. Notably, Gaia absolute astrometry is expected to identify (tens of) thousands of new exoplanets and brown dwarfs (BD) in the near future (Perryman et al. 2014; Sahlmann et al. 2015; Holl et al. 2022; Gaia Collaboration 2023a; Holl et al. 2023).

The current number of exoplanet candidates identified with third Gaia data release (GDR3) astrometry (72; Gaia Collaboration 2023a) is still much below expectations. Therefore, a major challenge is to exploit the Gaia data currently made public in the online catalogs in its most recent data release (DR3; Gaia Collaboration 2021) to detect unknown exoplanet candidates, as nicely illustrated with the discovery of AF Lep b (Mesa et al. 2023; Franson et al. 2023; De Rosa et al. 2023). Incidentally, Gaia's astrometry can also help validate (or reject) candidate exoplanets detected by other means (RV, transit, imaging) and further characterize them (Kiefer et al. 2019, 2021; Kiefer 2019; Kervella et al. 2019, 2022; Brandt et al. 2019; Dalal et al. 2021; Brandt 2021; Feng et al. 2021, 2022; Xiao et al. 2023; Philipot et al. 2023b,a) or aid in assessing the existence of a companion (possibly supplementary) of a given star or set of stars of interest.

With this objective, we set up a tool called GaiaPMEX for Gaia DR3 proper motion anomaly and astrometric noise excess based on the original works of Kiefer et al. (2019); Kervella et al. (2019); Kiefer (2019); Kiefer et al. (2021); Kervella et al. (2022) that allows for determination of the mass of possible candidate companions and their relative semi-major axis in relation to their primary star (abbreviated to sma hereafter) from consideration of, individually or in combination, the constraints from the proper motion anomaly (hereafter PMa; Kervella et al. 2019; Brandt 2021; Kervella et al. 2022), the astrometric excess noise (AEN; see Kiefer et al. 2019; Kiefer 2019; Kiefer et al. 2021), and the renormalized unit weight error (ruwe; see Lindegren et al. 2018, 2021). This tool models, within a Bayesian framework, the observed AEN, ruwe, and PMa through simulated outcomes of Gaia's observations of a source if it had a companion of a given mass and sma. It leads to a 2D confidence map of the companion mass and sma. Introducing this tool is the purpose of the present paper; a series of further papers will report the results of its application on other systems.

In Sect. 2, we recall the definitions of AEN, ruwe, and PMa. In Sect. 3, we describe our reverse-engineering method to determine the noise levels of Gaia's observations of individual sources. In Sect. 4, we explain the modeling of any star’s orbital motion due to a companion and the simulations of Gaia astrometric measurements of that star. In Sect. 5.1 we define the PMa, ruwe, and AEN astrometric signatures. In Sect. 6, we present the GaiaPMEX tool in detail. In Sect. 8, we show illustrative examples of the application of GaiaPMEX on a few chosen sources. Finally, in Sect. 9, we discuss the perspectives opened by the application of this tool regarding the detection of exoplanets and brown dwarfs using Gaia.

2 Astrometric excess noise, RUWE, and proper motion anomaly

2.1 The astrometric excess noise

The AEN of a source, as introduced in Gaia Collaboration (2016), is the excess of scatter in the residuals of along scan angle measurements compared to the astrometric displacement of the source modeled as a single-star, that includes position, linear proper motion and parallactic motion. At each epoch of transit of a source along one of the detectors, there is a specific scan direction, the along scan direction (AL), along which the source image is moving during the rotation of the spacecraft. The position of the source on the detector can be determined in 2D, since there is also an across scan (or AC) direction, but it is much less precisely measured along the AC than along the AL direction. Therefore, in all Gaia data releases, only the AL angles are used as astrometric measurements to determine the main astrometric data of a source (Lindegren et al. 2016, 2018, 2021).

In the GDR3, as in previous releases, the process of fitting the astrometric data is iterative. At each iteration, individual errors, σ_AL, of AL angle measurements performed during a transit of a star on the detector are estimated or updated and then used to calculate a χ². Since DR2 (Lindegren et al. 2018), a spacecraft attitude excess noise σ_att is quadratically added to σ_AL in the calculation of the χ². Its amplitude is typically about 0.076 mas, while individual measurement errors are within 0.05−0.15 mas (Lindegren et al. 2021). Both form a “formal error” ${\sigma }_{\text{formal}}=\sqrt{{\sigma }_{\text{att}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}$. We give more details and estimation of their variations with respect to the magnitude, color, right ascension (RA) and declination (Dec) of targets in Sect. 3.1. Their time series will only be known upon the publication of the DR4. The monitoring of the residuals root-mean-square shows that the measurement and excess attitude errors are constant most of the time, with rare deviations (see Fig. A.2 in Lindegren et al. 2021). We thus assume in the following that the attitude excess noise of a time series for any given target remains relatively constant in time. With this assumption, the χ², as it appears in the archives (namely astrometric_chi2_al), written here ${\chi }_{\text{astro}}^2$ is $${\chi }_{\text{astro}}^2=\frac{\sum _{\ell =1,N}{R}_{\ell }^2}{{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2},$$(1)

where R_ℓ are the residuals of the N astrometric measurements (astrometric_n_good_obs_AL) after subtraction of the fitted model. If some additional calibration noise – that is, a non-subtracted residual instrumental jitter beyond the attitude excess noise – or real astrometric signal were to be present, it would not be accounted for in the formal errors used to calculate the χ² and the reduced χ² would be larger than 1. Deviations of the reduced χ² beyond 1 are accounted for in the AEN (astrometric_excess_noise). To calculate the final uncertainties of fitted parameters of a given target, the AEN is quadratically added to the formal error of any astrometric measurements such as to impose a reduced χ² of 1. Still assuming that the errors are uniform along the time series, the AEN is related to the residuals sum of square through $${\text{AEN}}^2+{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2=\frac{\sum _{i=1,N}{R}_i^2}{N-5},$$(2)

counting N − 5 degree of freedom (DOF), with five parameters fit to the astrometry. The exact definition of the AEN involves possibly non-uniform errors and it is fixed iteratively during the reduction. Its value might thus slightly deviate from this definition. The level of the additional calibration noise still present in the data, not accounted for in the formal error of Eq. (1) but contributing to the AEN in Eq. (2), strongly depends on the magnitude and the color of the observed targets (Lindegren et al. 2016, 2018, 2021). We invented a method to estimate it for any source from the whole Gaia catalog of bright sources with magnitude G < 16, as thoroughly explained in Sect. 3.4.

The identification of many zero-valued AEN for sources dimmer than G=13 led us to become aware of an issue with the estimation of the calibration noise in the GDR3’s reduction. When the χ² was smaller or equal to the 95 th-percentile of the χ² distribution with N_DOF degree of freedom, that is, when the reduced χ² was smaller than $1+\frac{N_{\text{DOF}}}{\sqrt{2N_{\text{DOF}}}}$, the AEN was almost always fixed to zero in the archives (Lindegren et al. 2012). For sources dimmer than G=13, the attitude excess noise, common to all sources observed at the same epoch on the detector, overestimates the calibration noise and thus the format error to compute the χ² (Lindegren, priv. comm.). This led to an AEN wrongly fixed to zero for many sources beyond G=13, thus erasing any information on supplementary signals. Below G=13 this problem did not arise, because the calibration noise was conversely underestimated by the attitude excess noise, leading always to strictly positive values of the AEN. Our present understanding is that the AEN can be used as a binarity indicator and even used to characterise orbital motion, as long as the calibration noise and the attitude excess noise are both well known, and that the zero-valued AEN are discarded. The renormalized unit weight error, discussed in the next section, being directly proportional to the reduced χ² will be less problematic in this regard because it is not cut off below some value.

2.2 The renormalized unit weight error

An alternative to overcome the above issue is to use the renormalized unit weight error, or ruwe, instead of the AEN. By definition (Lindegren et al. 2018), $$\mathtt{r}\mathtt{u}\mathtt{w}\mathtt{e}=\frac{1}{u_0}\times \sqrt{\frac{{\chi }_{\text{astro}}^2}{N-5}},$$(3)

where u₀ is a factor that depends on magnitude and color. It can be determined from the GDR3 database values of ${\chi }_{\text{astro}}^2$ or astrometric_chi2_al), ruwe and number points N or astrometric_n_good_obs_al. With the approximate Eqs. (1) and (2), the ruwe and the AEN are directly associated: $$\mathtt{r}\mathtt{u}\mathtt{w}\mathtt{e}\approx \frac{1}{u_0}\times \sqrt{\frac{{\mathrm{A}\mathrm{E}\mathrm{N}}^2+{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}{{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}}.$$(4)

The ruwe is a unit-less scalar, but by the use of this formula, it could be conveniently transformed to an AEN. With a unit of angle – expressed in milli-arcsecond (mas) in the catalog – the AEN is directly commensurate to any possible astrometric motion – in au if divided by the parallax. A large value (>1.4) for a source is often accepted as indicating binarity. In many cases, this is indeed true, but it is nevertheless a misinterpretation of the DR3’s documentation, rather cautiously indicating that well-behaved sources (single or not), that is, for which the five-parameter fit gives a reasonably good fit, should have ruwe<1.4. We noticed, indeed, that the deviation of the ruwe above 1 in the GDR3 catalog is sometimes unreliable as a binarity indicator. This is most frequent for sources whose Gaia data were fit using six parameters (astrometric_params_solved=95). The case of the star β Pictoris is an excellent counter-example, with a ruwe of 3.07, that, we show in Sect. 8.5, can be explained by noise only, for this very bright star.

2.3 The proper motion anomaly

The PMa, as initially introduced in Kervella et al. (2019), is the proper motion offset between the Hipparcos-Gaia average proper motion (with a baseline ∼24.5 years), and the GDR3 fitted linear proper motion (with a baseline of 36 months). It thus measures an acceleration of the primary star due to the presence of a long-period secondary companion. The most recent measurements of PMa can be found in Kervella et al. (2022) as well as in Brandt (2021) with a different treatment of the global reference frames matching between GDR3 and the Hipparcos International celestial reference system (ICRS for short). In brief, noting μ the 2D proper motion, with index HG for “Hipparcos-Gaia", $${\mu }_{\mathrm{P}\mathrm{M}\mathrm{a}}={\mu }_{\mathrm{G}\mathrm{D}\mathrm{R}3}-{\mu }_{\mathrm{H}\mathrm{G}}.$$(5)

The non-linear perspective acceleration is assumed to be corrected in μ_HG. In this sense, μ_HG is the average 3-D Hipparcos-Gaia linear proper motion projected on the tangent plane at GDR3 epoch. Moreover, the effect of perspective acceleration is taken into account in the GDR3 astrometric solution, and μ_GDR3 is thus already the proper motion of the star in the Cartesian tangent plane. With these definitions, we can thus consider that μ_PMa is the projected tangential PMa as measured from a reference frame co-moving with the system’s barycenter.

As illustrated in Fig. 1, for any mass and sma, there always exists a longitude of ascending node Ω that fits any PMa position angle, the mass and sma can only be constrained from ‖PMa‖. When referring to PMa in the rest of the text, we thus always refer to ‖PMa‖, that is: $$\mathrm{P}\mathrm{M}\mathrm{a}=\parallel {\mu }_{\mathrm{G}\mathrm{D}\mathrm{R}3}-{\mu }_{\mathrm{H}\mathrm{G}}\parallel .$$(6)

Even though the PMa measures a variability in the proper motion of a star, noise in the astrometric measurement may induce a non-zero PMa. We assess the significance of the PMa in Sect. 5.2.2. It turns out that this is different than comparing the value of the PMa to its error bar that is calculated from the published measurement errors from Gaia and Hipparcos.

3 Noises and errors in Gaia observations

One the issues with interpreting correctly the AEN, ruwe and PMa as indicators of binarity and even measurements of companion’s properties, is our ignorance, a priori, of the noise budget in those quantities. Indeed, measurement noise and instrumental calibration noise participate at a certain degree in the excess of residuals beyond “formal error” (see Sect. 3.1 for a definition), as well as in any excess of proper motion fitted to noisy astrometric data. To complicate the task further, the level of those noises and error in Gaia data for any given source is not published and thus unknown to the community. As a prerequisite to the functioning of GaiaPMEX, whose goal is to model the astrometric motion beyond noise in AEN, ruwe and PMa, we thus present, in the following sections, a method that we developed to determine the noises and error levels in Gaia data for any source with G<16.

3.1 The formal error

What we call the formal error, σ_formal, is the unknown error that appears in the denominator of ${\chi }_{\text{astro}}^2$ in Eq. (1), that is, $\sqrt{{\sigma }_{\text{att}}^2+{\sigma }_{\text{AL}}^2}$. Combining this equation with Eq. (2) led us to express a simple approximation of ${\chi }_{\text{astro}}^2$ with respect to AEN: $${\chi }_{\text{astro\hspace{0.17em}}}^2=(N-5)\times \frac{{\mathrm{A}\mathrm{E}\mathrm{N}}^2+{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}{{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}.$$(7)

The formal error could thus be guessed by inverting this formula for all the sources observed in the GDR3 that have an AEN not compatible with 0 mas, that is, with an astrometric_excess_noise_sig ≥2 (Lindegren et al. 2012): $${\sigma }_{\text{formal\hspace{0.17em}}}=\mathrm{AEN}{(\frac{{\chi }_{\text{astro\hspace{0.17em}}}^2}{N-5}-1)}^{-1/2}.$$(8)

This estimate of the typical errors used in the ${\chi }_{\text{astro}}^2$ for several million sources allowed us to study the impact of magnitude, color, RA and Dec on GDR3 astrometric errors, and, more specifically, as we show in Sects. 3.2 and 3.3 of the attitude excess noise and the AL measurement errors. We adopted the bins defined in Table 1. To adapt to more rapid variations of the errors between magnitudes of 10.5 and 13.5, we adopted a smaller bin size ∼0.1 between 10.5 and 12.5. Moreover, a strong discontinuity in the errors occurs at G=13. It is related to the change in window class (or WC) from G<13 (WC0) to G>13 (WC1). It goes with a different level of charge transfer inefficiency (CTI) that is increasing in WC0 up to G=13, but strongly decreasing in WC1 (Lindegren et al. 2021). Because of this, we had to adopt an even smaller bin size of 0.05 between 12.5 and 13.5 G-mag.

In each magnitude-color or RA–Dec bin, we calculate the median formal error of all sources in these bins, respectively σ_formal (mag, color) and σ_formal (RA, Dec). This gives the relationship between formal error and magnitude & color or RA & Dec. Figure 2 shows the variations of the median formal error for the Gaia sources with AEN>0 mas and brighter than G=16, with respect to those parameters. We consider separately the sources whose data were fit by a five parameters model (astrometric_params_solved = 31), hereafter called ‘5p’ dataset, and those whose data – astrometry plus photometry – were fit by a six parameters model (astrometric_params_solved = 95) that includes an astrometric estimate of the effective wavenumber, or pseudocolor, v_eff, hereafter called ‘6p’ dataset. The sources whose data were only fit by a two-parameter model were not considered. Interestingly, this shows that the most crowded regions of the Galaxy have a larger error on average, as well as the sources with a G-mag of about 7–9. This latter dependence on magnitude agrees well with Fig. A.1 of Lindegren et al. (2021).

3.2 Attitude excess noise

The Gaia spacecraft attitude is modeled during the data reduction. It converts a rigid solid-body motion within the ICRS reference frame into Gaia's own reference frame where the CCDs are fixed (Lindegren et al. 2021). It thus models part of the path followed by any source along the detectors during a transit. This model suffers from time-dependent attitude excess noise, due, for example, to micro-clanks, calibration errors, etc., and has a typical level of 76 μ as on average (Lindegren et al. 2021). The attitude excess noise varies with time but at a given epoch all stars observed share a common attitude excess noise (Lindegren, priv. comm.). Depending on the magnitude and the color, σ_att tends to over/under estimate the calibration noise (see also Sects. 2.1 and 3.4).

The time-dependency of formal errors or attitude excess noise is not available. We can only assume that for any source, those errors are relatively constant (see, e.g., the Fig. A.3 in Lindegren et al. 2021 for an example with time-dependent attitude excess noise variations). Nonetheless, for a given source with specific RA & Dec direction, the attitude excess noise is probed at more-or-less regularly spaced epochs because of the scanning law of the spacecraft. Sources in different directions might thus probe disjoint sets of attitude excess noise values, and the mean attitude excess noise might thus depend on the RA & Dec direction. Being fixed, by construction, for all stars observed at the same epoch, the attitude excess noise do not dependent on magnitude or color.

For any source in the GDR3 database that has AEN>0 mas, quadratically removing the AL measurement error from the formal error leads to the attitude excess noise. To do this computation, we need to know the AL measurement error for any source. By conversely quadratically subtracting the attitude excess noise from the formal error one in fact can estimate the AL measurement error. At any bin of magnitude & color, σ_formal (mag, color) is the median formal error among all sources in that bin, distributed on all directions of the sky. We thus expect that, at any magnitude & color, the median σ_att is close to 76 μ as. This led to a first estimation of the AL measurement error, with respect to the magnitude and the color of the source, by applying $${\sigma }_{\mathrm{A}\mathrm{L}}(\text{mag},\text{color})=\sqrt{{[{\sigma }_{\text{formal}}(\mathrm{m}\mathrm{a}\mathrm{g},\text{color})]}^2-{0.076}^2}.$$(9)

This estimation is refined in Sect. 3.3. The AL measurement error of a given source depends mainly on the optical properties associated with a CCD measurement of its point spread function (PSF) on the detector, thus related to the magnitude and the color of the source. Linearly interpolating through this magnitude-color relationship, we can estimate the AL measurement error for any source of given magnitude and color (within available convex hull), or σ_AL, mc. Our best guess of the attitude excess noise for any source can then be obtained by quadratically subtracting this σ_AL, mc from the formal error: $${\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}=\sqrt{{\sigma }_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}^2-{\sigma }_{\mathrm{A}\mathrm{L},\mathrm{m}\mathrm{c}}^2}.$$(10)

To allow for estimation of σ_att even if a source’s AEN is compatible with 0 mas, and to smooth out scatter among sources with a similar sky location, we calculated a median attitude excess noise in every RA-Dec bins described in Table 1. Those median attitude excess noises are given in Table H.1. Figure 3 shows the dependence of the median σ_att(RA, Dec) with the sky direction. It shows a strong dependence on this parameter, with more pronounced error, up to 0.13 mas, in crowded regions, such as the Magellanic clouds and the center of the Milky Way. For any given source, an estimation of the effective level of attitude excess noise is determined by linearly interpolating the RA-Dec relationship at the RA and Dec of the source. In the rest of the article, this interpolated value is called σ_att.

3.3 Along-scan angle measurement error

Once σ_att is estimated for any source given their RA-Dec direction, it is straightforward to determine the σ_AL for all the sources with AEN > 0 mas. We simply used $${\sigma }_{\mathrm{A}\mathrm{L}}=\sqrt{{\sigma }_{\text{formal}}^2-{\sigma }_{\text{att}}^2}.$$(11)

As for σ_att above, to allow for estimation of σ_AL, even if a source’s AEN is compatible with 0 mas and to smooth out scatter among sources with a similar magnitude and color, we calculated a median attitude excess noise in every magnitude–color bins are described in Table 1. The median AL measurement errors are given in Table H.2. It is available online with only an extract shown here at a Bp − Rp close to that of GJ 832, that is, Bp − Rp = 2.2. Figure 4 shows the dependence of the median σ_AL with the magnitude and color. For any given source, an estimation of the effective level of AL measurement error is determined by linearly interpolating the magnitude-color relationship at the G-mag and Bp − Rp color of the source. In the rest of the article, this interpolated value is called σ_AL.

3.4 Calibration noise

The level of the calibration noise truly present in the data depends mainly on the magnitude and the color of the observed sources (Lindegren et al. 2016, 2018, 2021). For any source, we calculated a normal model of the ${\chi }_{\text{astro}}^2\sim \mathcal{N}({\mu }_{{\chi }^2},{\sigma }_{{\chi }^2})$, as thoroughly detailed in Appendix D. It accounts for the correlations between the co-adjacent astrometric AL angle measurements performed at the same epoch. The mean μ_χ² of the distribution of the ${\chi }_{\text{astro}}^2$ is related to σ_att, σ_AL, and σ_calib by $${\mu }_{{\chi }^2}=\frac{N_{AL}}{{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}[(N_{\mathrm{F}\mathrm{o}\mathrm{V}}-5){\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}^2+N_{\mathrm{F}\mathrm{o}\mathrm{V}}{\sigma }_{\mathrm{A}\mathrm{L}}^2]$$(12)

where N_FoV is the number of field of view (FoV) transits on the detector, and N_AL is the average number of AL angles collected per transit, that is, ≈int(N/N_FoV). In the Gaia archives, the total number of AL angle measurements N is given by astrometric_n_good_obs_AL, while N_FoV is given by astrometric_matched_transit. This equation leads to an expression of the σ_calib for any source assumed single, that is, $${\sigma }_{\text{calib}}=\sqrt{\frac{{\chi }_{\text{astro}}^2\times ({\sigma }_{\text{att}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2)-N_{AL}{N}_{\mathrm{F}\mathrm{O}\mathrm{V}}{\sigma }_{\mathrm{A}\mathrm{L}}^2}{N_{AL}(N_{\mathrm{F}\mathrm{O}\mathrm{V}}-5)}}.$$(13)

In any of the magnitude and color bins (Table 1), the best estimation of σ_calib is thus that of single sources. The sources are separated into single and multiple stars, whose rate N (multiple)/N (sources) is unfortunately unknown. The distribution of σ_calib in a given bin is thus the combination of both populations. According to Gaia's DR2 documentation², rather than the median, one can more safely rely on the mode of the unit weight error (UWE) distribution³ to locate the median of single star’s distribution. Indeed, the mode is shown to be less affected by multiplicity than the median and is thus a better approximation of single star’s median. Conversely to what is adopted in the documentation, we have found that the 41st-percentile is not always a good approximation of the mode. We thus rather localized the mode in the σ_calib distribution by iteratively excluding sources with a σ_calib larger than twice 1.483 × MAD(σ_calib) above the median, where MAD is the median absolute deviation. We then defined the mode as the median of this reduced distribution. We found that 3 iterations were necessary and enough to localise the mode. Figure E.1 shows this mode localization in the cumulative density functions of the σ_calib distribution at some magnitude-color bins.

Figure 5 shows the distribution of σ_calib with respect to G-magnitude and Bp − Rp color for the 5p and 6p datasets. The σ_calib of the 6p dataset are systematically higher than the 5p dataset. This is an effect of the poorer-quality of the fit for those stars. For them, a large AEN or ruwe has to be interpreted with care. In both dataset, σ_calib and magnitude are strongly correlated, especially for bright stars with G-mag <6, and to a lesser extent for early types with Bp − Rp<0.5 and late-types with Bp − Rp > 2.5.

4 Modeling of Gaia and Hipparcos astrometry

For later use, we define in this section our process for modeling Gaia and Hipparcos astrometric data. Our aim is to model the key data, namely PMa, AEN, and ruwe, which allow us to characterize the presence of companions and determine their main parameters, such as mass and sma. For any set of fixed companion, star and orbital parameters, we modeled by simulation the system’s photocenter orbit as if it was observed by Gaia or Hipparcos. In doing so, we accounted for instrumental and measurement noises in Gaia and Hipparcos data, and then performed a five-parameter fit of those datasets. We obtained residuals, as well as proper motion and centroid simulated measurements at GDR3 and Hipparcos epochs, respectively 2016.0 (Gaia Collaboration 2021) and 1991.25 (van Leeuwen 2007). We explain our method and the technical details of the simulations and the fit procedures in the following sections.

4.1 Modeling photocenter orbits

The core of the orbit modeling is the same as the one used in the “Gaia astrometric noise simulation to derive orbital inclination” tool or GASTON for short (Kiefer et al. 2019, 2021). We always consider a 2-body system, with a primary A and a secondary B, possibly planetary, brown dwarf or stellar. We fix the reference frame of the orbit to be the system’s barycenter. With Keplerian parameters fixed for this system, we model the orbit of the photocenter of the system on the plane of the sky. The photocenter semi-major axis, a_phot is determined from the total system’s semi-major axis, that is, the relative semi-major axis of the companion to the primary star; here written sma, through (Kiefer et al. 2021): $$a_{\text{phot}}=\mathrm{sma}(\beta -B)\varpi$$(14)

with ϖ the parallax, β = q/(1 + q) the mass fraction, and B = L₂/(L₁ + L₂) the luminosity fraction. The relative luminosity of the secondary over the primary is determined from semi-empirical mass-luminosity relation on the main sequence, at a typical age of 5 Gyr (see Kiefer et al. 2021 for more details). By default, we consider that the secondary may contribute to the photocenter’s position. Depending on the case at hand, one may instead consider a dark companion, whose luminosity is thus not contributing to the photocenter’s displacement. An illustration of the result of assuming instead a dark companion is shown in Sect. 8 for the case of α CMaB, that is, Sirius B, and whose companion Sirius A is resolved by Gaia and thus not contributing to α CMa B’s photocenter’s displacement.

The modeled orbits are then sampled at specific epochs, along specific directions, according to the scan law of Gaia and Hipparcos during their observation campaign. Noise is finally added to the individual measurements in a way that is specific to each instrument. This is further explained in the next Sects. 4.2 and 4.3.

4.2 Gaia DR3 sampling, scan-law, and noise

We sampled the modeled orbits at the GDR3 FoV-passage epochs and along the AL direction. The Gaia spacecraft is composed of two FoV, the ‘preceding’ and the ‘following’. They are separated around the spin axis of the spacecraft by a basic angle of 106.5^∘ (Lindegren et al. 2012; see also Fig. 6). At a given epoch, the spin axis is moreover oriented in a certain direction conferring to the detector a certain orientation of its main axis, the AL axis, u_AL. The law of the position angle (or PA) of u_AL through time can be found in the GOST. Six to 9 astrometric measurements are performed at the same epoch during the transit of the source across the detector thanks to Gaia's spin, at a speed of 60^′′/min (Lindegren et al. 2012).

The date of passage of a star on the Gaia detector and the PA of u_AL can be predicted accurately using the Gaia observation forecast tool (GOST for short). However, this tool is only accessible online⁴. Instead, we built a code that performs the same predictions using the spacecraft scan law accessible from the commanded_scan_law database. As explained in Fig. 6, we calculate the angle between the direction of each detector (α_d, δ_d) with the direction of the star at the GDR3 epoch (α_s, δ_s). We used a gnomonic projection (see, e.g., Calabretta & Greisen 2002) to transform this angle into a vector on the plane of the detector. More specifically, there is a relationship between (η, ζ) the AL and AC coordinates on the detector, and the difference of coordinates between the pointing direction of the detector and the direction of the star. Moreover, given that the AL-direction is oriented at a PA=θ_AL, we found this relationship to be $$\begin{array}{rl}\eta =& \frac{\cos {\delta }_{\mathrm{s}}\sin ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}{\sin {\delta }_{\mathrm{s}}\sin {\delta }_{\mathrm{d}}+\cos {\delta }_{\mathrm{s}}\cos {\delta }_{\mathrm{d}}\cos ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}\cos {\theta }_{\mathrm{A}\mathrm{L}}\\ & -\frac{\sin {\delta }_{\mathrm{s}}\cos {\delta }_{\mathrm{d}}-\cos {\delta }_{\mathrm{s}}\sin {\delta }_{\mathrm{d}}\cos ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}{\sin {\delta }_{\mathrm{s}}\sin {\delta }_{\mathrm{d}}+\cos {\delta }_{\mathrm{s}}\cos {\delta }_{\mathrm{d}}\cos ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}\sin {\theta }_{\mathrm{A}\mathrm{L}}\end{array}$$(15) $$\begin{array}{rl}\zeta =& \frac{\cos {\delta }_{\mathrm{s}}\sin ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}{\sin {\delta }_{\mathrm{s}}\sin {\delta }_{\mathrm{d}}+\cos {\delta }_{\mathrm{s}}\cos {\delta }_{\mathrm{d}}\cos ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}\sin {\theta }_{\mathrm{A}\mathrm{L}}\\ & +\frac{\sin {\delta }_{\mathrm{s}}\cos {\delta }_{\mathrm{d}}-\cos {\delta }_{\mathrm{s}}\sin {\delta }_{\mathrm{d}}\cos ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}{\sin {\delta }_{\mathrm{s}}\sin {\delta }_{\mathrm{d}}+\cos {\delta }_{\mathrm{s}}\cos {\delta }_{\mathrm{d}}\cos ({\alpha }_{\mathrm{s}}-{\alpha }_{\mathrm{d}})}\cos {\theta }_{\mathrm{A}\mathrm{L}}.\end{array}$$(16)

Then we imposed that this vector should be contained within the used area of the detector. The zero origin of η and ζ is not located at the center of the detectors and is different in the two FoVs, as explained in (Lindegren et al. 2016). In terms of CCD (AL × AC) unit, compared to the center of the detectors, they are located at (−2.5,+0.5) for the preceding FoV and at (−2.5,−0.5) for the following FoV. The detectors have a common dimension of 0.66 × 0.74 degree² with a grid of 9×7 CCDs. Dead zones are the wavefront sensor WFS2 (Gaia Collaboration 2016) and the area exterior to the detector. We assumed that CCD regions at less than a quarter of a CCD-sized distance to a dead zone is also a dead zone. This led to a better match of the number of predicted transits with the actual number of transits for any given star. We rejected a detection if the star fell on the dead zones. Figure 7 shows a representation of the detector and the geometry of the assumed dead zones, with GJ 832’s predicted average positions and AL scan direction orientations.

The time sampling of the scan law is ∼11 s. Therefore, during any transit of a source on the detector, several epochs are found, whereas only one epoch is required per transit. The spacecraft rotates at 1^′′/s, and the largest of the diagonals of the detector have a dimension of 1^∘. We thus determined, for any transit, the average epoch and average PA from all the predicted transit epochs found within a 1 hour window. We thus obtained for a given star all its theoretical epochs of transits through any of the two detectors with their corresponding PA of the AL direction. As a final step, we removed the epochs that fell at known gaps published in the GDR3 catalog⁵. The position of GJ 832 on the detector and the PA of the scan directions during its transits as retrieved from the scan-law is shown in Fig. 7.

We verified that the retrieved FoV transits matched those predicted by the GOST. Moreover their number are always close to those given by the astrometric_matched_transit in the GDR3 catalog. We noted that our calculation, consistently with the GOST, sometimes overestimated the number of actual FoV transits retained to calculate the astrometric solution in the GDR3. This happens most often to bright stars, thus indicating an effect of saturation that led to removal of some of the transits in the solution. In those cases, we randomly selected the correct number of epochs effectively used by Gaia in the GDR3 among all retrieved epochs of FoV passages.

We then assumed that the Gaia's AL measurements along the AL direction are distributed according to a normal law with as standard deviation, the noise σ_AL determined by the G and Bp − Rp of the source in consideration (see Sect. 3.3). Adding to this error, we added an epoch-specific offset randomly drawn from a normal distribution with standard deviation σ_calib, determined with respect to the G and Bp − Rp of the source (see Sect. 3.4). One such simulation is shown in Fig. 8.

4.3 Hipparcos sampled epochs

The Hipparcos-2 intermediate astrometric data (IADs; van Leeuwen 2007) are necessary to model the PMa as determined by Kervella et al. (2022). We here only focus on the PMa between GDR3 and HG baselines. The location of the centroid of the Hipparcos-2 data in the source barycenter reference frame has to be determined for any modeled orbit to find the HG average proper motion between epochs 1991.25 and 2016.0. In the Hipparcos-2 database the published source centroid is located along the fitted solution. This is not adequate for us because Hipparcos-2 used, when possible, more elaborate models including acceleration or orbital motion. However, Kervella et al. (2022) only considers the result of a five-parameter fit of the Hipparcos-2 IADs to derive the Hipparcos-2 proper motion and the location of the Hipparcos-2 source centroid at epoch 1991.25.

Therefore, for any source also observed with Hipparcos, we downloaded the IADs residuals from the Hipparcos-2 Interactive Data Access Tool⁶. These data include the orbit number (IORB), the epoch, the cosine and sine of the Hipparcos scan angle (related to the Gaia scan angle convention by ψ=θ_AL-π/2; Brandt 2021), the residuals of the fitted model (RES) and the formal errors (SRES).

We removed the data with negative or zero SRES that are rejected observations. To model Hipparcos observations of the astrometric displacement of the photocenter due to an orbital motion, we only need the part of the residuals at each IORB that cannot be due to supplementary non-modeled displacement. We thus calculated corrected residuals (CRES) by removing the local average from common IORB residuals. Then when simulating an orbit, those CRES are added along the Hipparcos scan direction.

Besides, to have an estimation of the typical dispersion associated with a source centroid position, we also calculated the Hipparcos-2 positional error for the considered source from the RA and Dec positional error published in the Hipparcos-2 catalog: $${\sigma }_{\mathrm{p}\mathrm{o}\mathrm{s},\mathrm{H}\mathrm{I}\mathrm{P}}=\sqrt{{\mathrm{e}}_{-}{\mathrm{R}\mathrm{A}}^2{\cos }^2(\mathrm{D}\mathrm{E})+{\mathrm{e}}_{-}{\mathrm{D}\mathrm{E}}^2}.$$(17)

An illustration of the Hipparcos data modeled for an arbitrary orbit is shown in Fig. 8.

4.4 Gaia and Hipparcos five-parameter model fit

For each modeled orbit, we applied a five-parameter⁷ fit to Hipparcos-2 and GDR3 simulated data. It included the RA-Dec centroid of data points, the linear proper motions μ_α and μ_δ, and the parallax ϖ. Given that we placed ourselves in the barycenter reference frame of the considered system, we thus fit the excesses of (positive or negative) offset, proper motion, and parallax only due to the presence of an orbital motion. Because of the orbital motion, the parallax measured in GDR3 deviates from the true value. Assuming that the current orbit was the true one, we first estimated the parallax error Δ ϖ from the fitted parallax excess in a first simulation. We then performed a second simulation, correcting the parallax by ϖ → ϖ − Δϖ. The 2D-fitted linear model is $$M(t)=(\delta \mathrm{R}\mathrm{A},\delta \mathrm{D}\mathrm{e}\mathrm{c})+(\delta {\mu }_{\alpha },\delta {\mu }_{\delta })(t-t_0)+\delta \varpi \mathrm{\Pi }(t)$$(18)

where the δμ are proper motion in α=RA cos Dec and δ=Dec tangent plane directions, Π(t) is the parallax ellipse depending on the coordinate of the star, and t₀ the Hipparcos-2 or GDR3 epochs, respectively 1991.25 and 2016.0. We note that the effect of perspective acceleration (see, e.g., Michalik et al. 2014; Halbwachs et al. 2023) that mainly affect high proper motion targets close to Sun, is already corrected in the GDR3, so being a second-order effect we can ignore it here (Lindegren et al. 2021). To compare this linear model to the Hipparocs and Gaia measurements, we needed to project this model onto the along-scan directions, with position angle θ_AL, determined at the sampled epochs along the orbit: $$\begin{array}{rl}{M}_{\mathrm{A}\mathrm{L}}\left(t_i\right)=& (\delta \mathrm{R}\mathrm{A}+\delta {\mu }_{\alpha }(t_i-t_0)+{\delta \varpi \prod \left(t_i\right)|}_{\mathrm{R}\mathrm{A}})\sin {\theta }_{\mathrm{A}\mathrm{L}}\left(t_i\right)\\ & +(\delta \mathrm{D}\mathrm{e}\mathrm{c}+\delta {\mu }_{\delta }(t_i-t_0)+{\delta \varpi \prod \left(t_i\right)|}_{\mathrm{D}\mathrm{e}\mathrm{c}})\cos {\theta }_{\mathrm{A}\mathrm{L}}\left(t_i\right).\end{array}$$(19)

We separated components along Dec (north) and RA (east) directions. We subtracted this five-parameter model from the simulated data and calculate the residuals. For Gaia, they are further used in comparison to the tabulated AEN or ruwe published in the GDR3 catalog as explained in Sects. 5.1 and 6.

To calculate the PMa, the fit GDR3 proper motion, δμ_GDR3, was combined with the average proper motion between the fit positions of the photocenter at the Hipparcos-2 reference epoch and the GDR3 reference epoch, δμ_HG, through $$\mathrm{P}\mathrm{M}\mathrm{a}=\delta {\mu }_{\mathrm{G}\mathrm{D}\mathrm{R}3}-\delta {\mu }_{\mathrm{H}\mathrm{G}}$$(20) $$\text{with\hspace{0.17em}}\delta {\mu }_{\mathrm{H}\mathrm{G}}=\frac{(\alpha ,\delta {)}_{\mathrm{G}\mathrm{D}\mathrm{R}3}-(\alpha ,\delta {)}_{\mathrm{H}\mathrm{I}\mathrm{P}}}{24.75\mathrm{y}\mathrm{r}}.$$(21)

This modeled PMa is compared to the PMa published in Kervella et al. (2022), as explained in Sects. 5.1 and 6.

5 The non-singleness of stars observed with Gaia

5.1 Astrometric signatures

To assess the non-singleness of stars from AEN, ruwe, and PMa, we defined (and introduce here) the “astrometric signatures”, further written as α. They properly quantify the deviation of AEN, ruwe and PMa beyond the level that they must have had if the sources were single.

5.1.1 The residuals astrometric signature

We first introduced the residuals unbiased estimator of variance a posteriori (UEVA) related to the ${\chi }_{\text{astro}}^2$ by $$\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}=\frac{{\chi }_{\text{astro}}^2}{N-5}\times ({\sigma }_{\text{att}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2).$$(22)

This quantity’s square root, also known as regression standard error (RSE), is an unbiased estimator of the data typical error in the considered sample of measurements. The UEVA of any source, considering Eqs. (3), (7) and (22), can be estimated in two different ways, either using the AEN or the ruwe: $${\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}_{\mathrm{a}\mathrm{e}\mathrm{n}}={\mathrm{A}\mathrm{E}\mathrm{N}}^2+{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2$$(23) $${\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}_{\text{ruwe}}={(\mathtt{r}\mathtt{u}\mathtt{w}\mathtt{e}\times u_0)}^2({\sigma }_{\text{att}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2).$$(24)

We consider both in the rest of the study. The UEVA of a single source (UEVA_single) only accounts for calibration noise and AL astrometric noise. The distribution of UEVA_single is obtained from the normal model of the χ² distribution of the residuals, determined in Appendix D, Eqs. (D.8) and (D.9). It led to a normal distribution $\mathcal{N}$(μ_UEVA,single, σ_UEVA,single) of the UEVA_single with $${\mu }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A},\text{single}}=\frac{N_{AL}}{N_{AL}{N}_{\mathrm{F}\mathrm{o}\mathrm{V}}-5}[(N_{\mathrm{F}\mathrm{o}\mathrm{V}}-5){\sigma }_{\text{calib}}^2+N_{\mathrm{F}\mathrm{o}\mathrm{V}}{\sigma }_{\mathrm{A}\mathrm{L}}^2]$$(25) $$\begin{array}{rl}{\sigma }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A},\text{single}}^2=& \frac{2N_{AL}}{{(N_{AL}{N}_{\mathrm{F}\mathrm{o}\mathrm{V}}-5)}^2}[N_{AL}(N_{\mathrm{F}\mathrm{o}\mathrm{V}}-5){\sigma }_{\text{calib}}^4\\ & +N_{\mathrm{F}\mathrm{o}\mathrm{V}}{\sigma }_{\mathrm{A}\mathrm{L}}^4+2N_{\mathrm{F}\mathrm{o}\mathrm{V}}{\sigma }_{\mathrm{A}\mathrm{L}}^2{\sigma }_{\text{calib}}^2].\end{array}$$(26)

Ingredients needed to calculate μ_UEVA,single and σ_UEVA,single for any source are thus G-mag, Bp − Rp, RA and Dec to estimate the noises; and N_FoV, as well as N_AL=int(N/N_FoV). The variation of the single-star RSE with respect to magnitude is compared to the level of σ_AL in Fig. 9. It compares well and agrees with the same curves plotted in Fig. A.1 in Lindegren et al. (2021) determined directly from the unpublished time series and images. This shows that our estimation of UEVA_single in Eq. (25) provides a reliable estimation of the ground level of the UEVA for any source.

For a non-single star, when adding the orbital motion, the UEVA may positively deviate from UEVA_single. We define the residuals astrometric signature as the angular excess that has to be quadratically added to the ground level UEVA_single to recover the UEVA of the residuals of the given source. It thus measures, in units of milli-arcseconds, the strength of non-singleness of the source. It is written α_UEVA and is formally defined as $${\alpha }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}=\sqrt{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}-{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}_{\text{single}}}.$$(27)

Using μ_UEVA,single as the expectation value of UEVA_single and the AEN (Eq. (23)) or the ruwe (Eq. (24)) to estimate UEVA, we can determine α_UEVA for any source. For a single star, as further developed in Sect. 5.2, because of the diverse astrometric noises that depends on the star’s magnitude, color and sky coordinates, the UEVA follows a broadened distribution that extends around UEVA_single. It leads to α_UEVA=0 if the UEVA is smaller or equal to UEVA_single and positive values otherwise. For non-single sources, α_UEVA may become strongly positive if the astrometric motion dominates over the astrometric noises. In that sense, α_UEVA is indeed an astrometric signature.

Figure 10 compares the α_UEVA calculated from UEVA_aen or UEVA_ruwe to respectively AEN or ruwe for both 5p and 6p datasets. It shows that the α_UEVA is almost equal to AEN beyond AEN=2 mas. The AEN could thus be directly interpreted as an astrometric signature in this regime. While there is a clear linear correspondance between ruwe and α_UEVA beyond ruwe=1.4, with an approximative slope of ∼0.2, the thickness of the relation makes a direct astrometric interpretation more difficult. This is worse for the 6p dataset where the range of possible α_UEVA for a given ruwe is even more spread out, due to the larger levels of noise compared to the 5p dataset. Irrespective of the dataset, below AEN=2 mas and ruwe=1.4 the range of possible α_UEVA is significantly broader. Some values are as low as 10 μ as, thus dominated by noise and insignificant.

5.1.2 The PMa’s astrometric signature

Similar to α_UEVA, a PMa’s astrometric signature can determine the excess that has to be quadratically added to PMa_single to recover the PMa of the given source measured in Kervella et al. (2022). Indeed, for a given orbit, the PMa should be given by a constant vector (c), that is, the “noiseless” orbital contribution, plus a stochastic vector (ξ), that is, the pure noise contribution. The expectation value of the square-norm of the PMa is thus $$\begin{array}{rl}⟨\parallel \mathrm{P}\mathrm{M}\mathrm{a}{\parallel }^2⟩& =⟨\parallel \mathbf{c}+\xi {\parallel }^2⟩\\ & =\parallel \mathbf{c}{\parallel }^2+⟨\parallel \xi {\parallel }^2⟩\end{array}$$(28)

since 〈ξ〉=0 and c is constant. The first term is the pure orbital contribution to the PMa, null if there are no orbital motion, that is, the astrometric signature that we seek. The second rightmost term is the squared-norm of PMa for a single star, or ${\text{PMa}}_{\text{single}}^2$. Therefore, we introduce α_PMa the PMa’s astrometric signature, formally defined as $${\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}=\sqrt{{\mathrm{P}\mathrm{M}\mathrm{a}}^2-{\mathrm{P}\mathrm{M}\mathrm{a}}_{\text{single}}^2}.$$(29)

It is theoretically possible to determine the typical level of PMa_single. However, conversely to the UEVA, there are no theoretical formula for estimating this distribution. We thus needed to perform simulations to estimate its mean and standard deviation, depending on the RA, Dec, Bp-Rp and G of the sources. We determined the typical distribution followed by PMa_single in Sect. 5.2.2.

5.2 The significance of the astrometric signatures

5.2.1 Significance of α_UEVA

As introduced in Sect. 5.1.1, the UEVA_single could be modeled by a normal distribution $\mathcal{N}$(μ_UEVA,single, σ_UEVA,single) whose parameters are written in Eqs. (25) and (26). Rigorously speaking, as explained in Appendix D, the UEVA_single is a linear combination of χ² and normal distributions, with the main terms distributed according to the χ² law. Under the prescription of Wilson & Hilferty (1931) (see also Canal 2005), the cubic-root of the UEVA more closely resembles a normal distribution. We thus used UEVA^1/3 and assumed that it followed a normal distribution $\mathcal{N}$(μ_1/3, σ_1/3). We approximated its parameters by μ_1/3=μ^1/3 and σ_1/3=σμ^−2/3/3 by applying error propagation. Figure 11 shows the resulting distribution of 100 000 simulations of UEVA^1/3 obtained when assuming that a source – here, for example, HD 114762 – is a single star, and compares it to the normal model.

The significance of α_UEVA naturally corresponds to the p-value of the UEVA^1/3 as calculated from either the AEN or the ruwe from the GDR3, within the $\mathcal{N}$(μ_1/3, σ_1/3) single-star distribution. This p-value is converted to an N−σ significance, following the “normal law” relationship between the 1−2−3-σ levels and the 31.6−4.6−0.27% p-values. We defined that α_UEVA is significant at N_σ−σ if ${\chi }_{\mathrm{F}\mathrm{o}\mathrm{V}}^2$ > chi2.ppf(x_Nσ/100,N_FoV−5) where chi2.ppf is the function of the python’s scipy.stats-module, and x_Nσ corresponds to the usual percentage at N_σ−σ (N=1: x_Nσ=68.3%; N=2: x_Nσ=95.4%; N=3: x_Nσ=99.73%). An α_UEVA significant at 2−σ would thus imply that the AEN or ruwe would have a less than 4.6% chance to occur if the star was single.

We showed in Fig. 12 the distribution of α_UEVA among all datasets. Its full range goes from about 0.1 μ as to 10 mas, but the α_UEVA with a significance >2−σ rather span the range that is beyond 10 μ as.

The sample with an α_UEVA significance larger than 2−σ contains 19.1 × 10⁶ sources, that is, about 25% of the sample of 75.9 × 10⁶ Gaia sources with G<16 and ϖ>0. Close to 9% of the sources (6.7 × 10⁶) have a significance larger than⁸ 8−σ. Figure 13 shows d f_detec(N_σ) the number of detections per bin of (N_σ, N_σ + dN_σ)–σ significance. They are compared to the theoretical distribution for only single stars, that is, d f_single=N_single × exp $(-N_{\sigma }^2/2)$ dN_σ/ $\sqrt{2\pi }$. The number of single stars, N_single, is determined with respect to an assumed binary (and multiple) rate in the sample, Γ_b=N_binary/N_sample, using N_single=(1−Γ_b) × N_sample. Assuming a Γ_b of 0% (blue curve in Fig. 13), the single star distribution cannot explain the rate of detections beyond 1.5−σ. Moreover, since many systems must be non-single, considering 100% of single-stars in the sample obviously overestimates the number of detection below 1.5−σ. A more realistic value for Γ_b can be found by assuming that all the sources at 0−σ significance must be single. We then fixed N_single in such a way that df_single(N_σ=0) matches df_detec(N_σ=0). As illustrated with an orange line in Fig. 13, we found a single-star rate of 47% and thus a Γ_b of 53% in this sample of sources brighter than G=16. The pollution from false-positives (FP) clearly equates/dominates over true positives below 2−σ. Beyond 2−σ, with a single-star rate of conservatively 47−100%, we estimate that about 9−19% of the selected binary or planetary systems could be single-star false positives. Thus, more than 80% of the selected sample beyond 2−σ significance are bona-fide binary, multiple or planetary systems. Beyond 3−σ the sample reaches 16% of the 76 millions sources, leading to single-star false positives rate of 0.9−2%. Thus, virtually all >3− σ sources are true binary, multiple or planetary systems, but about 14−24% of the non-single star sample is lost compared to using the 2−σ threshold. In the perspective of identifying binary, multiple or planetary systems candidates, we thus recommend, and adopt, that an α_UEVA, whose significance is above the 2σ level, has to be considered as a strong candidate, with a 9−19% chance that it actually is a single-star FP.

We use HD 114762 as an illustration for the significance of astrometric signature. This system was known for a long time for hosting a candidate planet HD 114762 b with an m sin i of 11 M_J (Latham et al. 1989; Kane et al. 2011). Confirming suspicion (Cochran et al. 1991; Hale 1995; Halbwachs et al. 2000), it was further shown using the astrometry from Gaia DR1, then GDR3, that HD 114762 b was actually an M–dwarf (Kiefer et al. 2019; Winn 2022; Gaia Collaboration 2023a). This was first shown using the significant value of AEN of 1.09 mas (Kiefer et al. 2019); then a proper astrometric orbit solution was obtained using the private GDR3 timeseries (Gaia Collaboration 2023a). In the light of the new framework introduced above, we found indeed that its GDR3 residuals astrometric signature rejected the single star hypothesis.

If HD 114762 was a single star – or orbited by an undetectable companion – the Gaia measurements would lead to AEN_single=0.251 ± 0.035 mas, ruwe_single=1.11 ± 0.13, and UEVA_single=0.076 ± 0.018 mas². In comparison, using the values of AEN=0.708 mas and ruwe=3.16 found in the GDR3 for HD 114762, we determined that UEVA_aen=0.514 mas² and UEVA_ruwe=0.607 mas². The astrometric signature deduced in both cases by applying Eq. (27) are α_astro,aen=0.662 mas and α_astro,ruwe=0.729 mas. The p-value of ${\text{UEVA}}_{\text{ruwe}}^{1/3}$ and ${\text{UEVA}}_{\text{aen}}^{1/3}$ in the distribution of ${\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}_{\text{single}}^{1/3}$ corresponds to a significance >9−σ (see Fig. 11). A single star may explain the AEN and the ruwe in close to 0% of the simulations. Both the AEN and the ruwe thus indicate the presence of a companion around HD 114762.

Assuming that the 84-day companion is responsible for this astrometric signature and that the star’s semi-major axis a_*∼α_UEVA/ϖ with the parallax ϖ=26 mas, and crudely applying M_c=a_*(M_*/M_⊙)^2/3(P/year)^−2/3 with M_*=1.05 M_⊙ (Winn 2022), we find a possible mass of the companion of 86 M_J. The mass obtained is indeed much larger than 11 M_J but is less than its most recent estimation ∼ 0.3 M_⊙ (Winn 2022). This just shows that α_UEVA is not a measure of a_*, leading us to develop a more sophisticated framework for interpreting α_UEVA and infer main parameters of companions, as explained in Sect. 6.

5.2.2 Significance of α_PMa

The significance of α_PMa corresponds to the p-value of the PMa^2/3 within the distribution of ${\mathrm{P}\mathrm{M}\mathrm{a}}_{\text{single}}^{2/3}$. This p-value is converted to an N−σ significance (see Sect. 5.2 for more details). As mentioned in Sect. 5.1.2, the distribution of ${\mathrm{P}\mathrm{M}\mathrm{a}}_{\text{single}}^{2/3}$ is not known a priori because it strongly depends on the time sampling of the astrometric observations further fitted by a five-parameter model to measure the GDR3 proper motion, as well as the HG relative positions. Conversely to α_UEVA, determining the distributions of α_PMa thus requires performing many simulations of a single star observations given the main parameters, scan law, and noises of the given source.

We used the system of GJ 832 to illustrate the significance of PMa beyond the single star hypothesis. GJ 832 is an M0V star at a distance of 5pc and with a mass of 0.48 ± 0.05 M_⊙. Its planetary system was discovered by Bailey et al. (2009), reporting one Jupiter-like planet with a period of 3416 ± 131 days and minimum mass of 0.64 ± 0.06 M_J. A second Earth-like planet was proposed for detection (Wittenmyer et al. 2014) but finally identified as a stellar activity artifact (Gorrini et al. 2022). The sma and minimum mass of GJ 832 b were further updated to 3.6 ± 0.4 au and minimum mass of 0.74 ± 0.06 M_J (Gorrini et al. 2022). This system is illustrative for us because it harbors one of the smallest mass planets leading to significant astrometric acceleration detected by combining Hipparcos and Gaia, and astrometric excess noise in the GDR3.

Figure 14 shows the PMa^2/3 distribution obtained for GJ 832 in the hypothesis that it is a single star. We used the values of noises, including the Hipparcos position error, that are given in Table B.1. It leads to PMa_single=0.060 ± 0.032 mas yr⁻¹ and ${\mathrm{P}\mathrm{M}\mathrm{a}}_{\text{single}}^{2/3}$=0.153 ± 0.054 (mas yr⁻¹)^2/3, while Kervella et al. (2022) measures PMa=0.565 ± 0.027 mas yr⁻¹, and equivalently PMa^2/3=0.683 ± 0.022(mas yr⁻¹)^2/3. If comparing to a zero-point PMa offset of zero, that is, when neglecting noise, the PMa^2/3 would have an apparent signal-to-noise ratio (S/N) of 30. But the distribution of ${\mathrm{P}\mathrm{M}\mathrm{a}}_{\text{single}}^{2/3}$, because of noise, strongly departs from zero as shown in Fig. 14. It leads to rather consider a positive zero-point offset and a larger error, that imply a more modest S/N of 10 for the PMa of GJ 832. We thus stress that the noise brings a major contribution to the zero-point offsets and errors that are used to determinate the S/N. For GJ 832, the p-value of PMa^2/3 in the distribution of ${\mathrm{P}\mathrm{M}\mathrm{a}}_{\text{single}}^{2/3}$, implies a significance of the PMa >9−σ (see Fig. 14). Thus, it firmly indicates the presence of a companion in this system.

Since the orbital period of GJ 832 b is ∼6 yr, that is, smaller than the Gaia-Hipparcos baseline of 24.5 years, we would tend to crudely interpret the PMa here as the average orbital speed during GDR3 observations. Assuming thus that PMa/$\varpi \propto \sqrt{G{M}_c^2/\left(\text{sma}{M}_{\star }\right)}$, it leads to estimate the mass of GJ 832 b at 3.6 au to M_c ∼ 0.1 M_J. This mass is on the order of magnitude of the expected mass, though underestimated. The tool that we developed in Sect. 6 allowed us also to properly infer the main parameters of companions from the knowledge of PMa.

6 The GaiaPMEX method

We aim at calculating the confidence regions of possible mass and sma of a companion, for a Gaia source with a given α_UEVA – either determined from AEN or ruwe – and/or α_PMa. To do so, rather than using a Markov-chains Monte-Carlo approach, which is time consuming, we perform a mass-sma grid search within a Bayesian framework. At each node of the mass-sma grid, as summarized in Fig. 15, the values determined in the GDR3 and in Hipparcos-Gaia-(E)DR3 studies (Brandt 2021; Kervella et al. 2022) are compared to modeled PMa and UEVA (see Sect. 4). We defined a likelihood function in Sect. 6.3.1. The adopted Bayesian framework is explained in Sect. 6.3.2 and summarized in Fig. 16.

6.1 A uniform grid of log M_c and log sma

To probe different orders of magnitudes for the mass and sma, we define a uniform 2D-grid on log M_c and log sma with log-scaled Δ mass and Δ sma, as sketched in Fig. 15. In each bin, we draw log M_c and log sma within a uniform distribution bounded by the bin edges. We also draw Keplerian parameters of the possible companion orbits, namely e the eccentricity, ω the periastron longitude, Ω the longitude of ascending node, ϕ the phase, and I_c the orbit inclination, following the distributions summarized in Table 2. Parallax, ϖ, and stellar mass, M_*, are drawn from normal distribution defined by prior knowledge on those parameters. The parallax is taken from the GDR3 catalog. By default, if M_* is not given as input, it is first searched in the GDR3 Coordination unit 8 (CU8) catalog (Gaia Collaboration 2023b). We specifically looked for the final luminosity age and mass estimate (or FLAME) available for more than 140 millions sources with M_*>0.5 M_⊙, the mass-Flame that is determined from combining photometry, parallax and stellar models. Complying with studies of stellar masses reported in Gaia Collaboration (2023a) and Babusiaux et al. (2023), as well as the recommendations in the Gaia documentation⁹ and more specifically the “astrophysical parameters” section¹⁰, the uncertainty on the mass is assumed ∼ 10% if M_*>0.7 M_⊙ and 0.1 M_⊙ if the mass is <0.7 M_⊙. If the star is a giant we only use the CU8 stellar mass if it is within 1−2 M_⊙, and we assume a 30%-uncertainty. If missing in this catalog, the stellar mass is instead estimated from the two first letter and digit of the spectral type given in the SIMBAD database¹¹ assuming that the star is on the main sequence.

6.2 Modeling of PMa and UEVA

In each bin of the (mass, sma)-grid, we modeled by simulation many Gaia observations of photocenter orbits, due to the reflex motion of the source due to a companion at given mass – or log (mass) – and sma – or log (sma).¹² For each modeled observation, we added noise in Gaia data (see Sects. 2.1 and 3.1), and then performed a five-parameter fit. Technical details on this modeling by simulation are explained in Sect. 4.1. Each simulation leaded to a value of PMa and a value of UEVA.

For each bin of the mass-sma grid drawn in Sect. 6.1, the full set of simulations obtained at a given (log M_c, log sma)-bin leaded to distributions of possible UEVA and PMa that would be measured in the GDR3 if the companion had such a mass and sma. We found that at least 100 simulations per bin are necessary to lead to a fine quality map. Running 300 simulations per bin performed better, leading to cleaner noiseless maps, with a computation time that was still tractable, though 3 × longer¹³. The maps that are shown here were obtained with 100 simulations per bin.

6.3 A Bayesian scheme for the 2D posterior distribution

From the modeled distributions, we determined a likelihood of the actual values of UEVA – either determined from AEN or ruwe – and PMa given the (mass,sma)_n-model $\mathcal{L}$= p(PMa, UEVA∣mass_n, sma_n). It is derived in Sect. 6.3.1 below. This likelihood was further used for a Bayesian inversion to determine the posterior distributions on mass and sma of the hypothetical companion, as explained in Sect. 6.3.2.

6.3.1 A log-likelihood of PMa and UEVA

One of the important issues met when defining a log-likelihood for PMa and UEVA, was to determine their probability density function (PDF) with ln $\mathcal{L}$=ln PDF (data). They typically do not follow a normal-law, since both quantities are always positive. We show in Sect. 5.1 that, because the PDF of the squares PMa² and UEVA are similar to χ²-distributions, their transformations to PMa^2/3 and UEVA^1/3 follow close-to-normal laws (Wilson & Hilferty 1931; Canal 2005). This ensured that our data distributions had PDF that were more compact and symmetrical as shown in Fig. 17.

Figure 17 illustrates the typical differences in the distributions of UEVA and UEVA^1/3 as well as PMa, PMa² and PMa^2/3, for the case of a hypothetical 100−M_J companion at 1 au around GJ 832. Figure 18 also shows the difference in the KS-test of the normal law with the distributions of UEVA or PMa and respectively those of UEVA^1/3 or PMa^2/3 modeled at different values of sma. Generally for most sma from 0.5 to 100 au, the X^1/3−transformation is a closer match to the normal law than simply UEVA or PMa. We thus adopt using UEVA^1/3 and PMa^2/3 for calculating the likelihoods.

For any (mass,sma)-bin of the grid, a Gaussian kernel density estimation (with the gaussian_KDE library from scipy) was performed on the distributions of the modeled UEVA^1/3 and PMa^2/3. This gives a good approximation of the true PDF of those quantities, as long as the sampling is dense enough, which is why we preferred using distributions that are not too extended and long-tailed. They can be used to derive the log-likelihoods as: $$\begin{array}{rl}\ln {\mathcal{L}}_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}& =\ln \left[{\mathrm{P}\mathrm{D}\mathrm{F}}_{{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}^{1/3}}\left({\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}_{\mathtt{r}\mathtt{u}\mathtt{w}\mathtt{e}}^{1/3}\right)\right]\\ \ln {\mathcal{L}}_{\mathrm{P}\mathrm{M}\mathrm{a}}& =\ln \left[{\mathrm{P}\mathrm{D}\mathrm{F}}_{{\mathrm{P}\mathrm{M}\mathrm{a}}^{2/3}}\left({\mathrm{P}\mathrm{M}\mathrm{a}}_{\mathrm{H}\mathrm{G}}^{2/3}\right)\right].\end{array}$$(30)

Those log-likelihood give the probability of the data given a certain (mass,sma)-companion model and distributions of other Keplerian parameters. They are then used in GaiaPMEX in the inversion of the Bayesian formula to obtain the probability of the models given the data, and confidence regions on mass and sma, as explained in Sect. 6.3.2.

6.3.2 Bayesian inversion

Every bin n of the grid has thus a likelihood ${\mathcal{L}}_n$ and a corresponding log-likelihood $\ln {\mathcal{L}}_n$. The bin n_max at which $\log {\mathcal{L}}_n$ is maximized, reaching $\ln {\mathcal{L}}_{\text{max}}$, is the (mass,sma)-model for which the data are best matching the UEVA and PMa distribution. At any other bin, we measured a likelihood ratio (LR) through $\mathrm{\Delta }\ln {\mathcal{L}}_n=\ln {\mathcal{L}}_{\text{max}}-\ln {\mathcal{L}}_n$.

To derive a probability function for mass and sma, we needed to determine at each bin n, what is the p-value of $\mathrm{\Delta }\ln {\mathcal{L}}_n$. This is summarized in Fig. 16. In the ideal case of the likelihood-ratio test of some null hypothesis, the Wilks theorem states that with a large number of data, $2\mathrm{\Delta }\ln \mathcal{L}=2(max\ln \mathcal{L}-\ln {\mathcal{L}}_{\text{null}})$ should follow a χ² distribution with k degrees of freedom (Wilks 1938; Silvey 1970); k being the difference between the maximum number of degree of freedom (DOF_max=N_data−N_param), and the number of degree of freedom in the region constrained by the null-hypothesis where some of the parameters are fixed (DOF_null=N_data−N_{param,unfixed}). Here, fixing the only 2 parameters to vary, the mass and the sma, we have k=DOF_max-DOF_null=2. However, since the number of data is small (N_data=2) the Wilks theorem does not apply.

To convert the LR at bin n into a p-value, we must find the empirical distribution of the LR at that bin. This was done by assuming that the models (mass ± Δ mass, sma ± Δ sma)_n is the true one and draw many possible UEVA and PMa from these models, as if they were those measured by GDR3. For each drawn UEVA & PMa, we apply the same grid search as explained in the above paragraph, finding the likelihood optimum. At the considered bin n, this led to a distribution of $\mathrm{\Delta }\ln \mathcal{L}$. The corresponding percentile p_n of $\mathrm{\Delta }\ln {\mathcal{L}}_n$ within this distribution has a frequentist interpretation, as for Δ χ² inference on confidence limits (Press et al. 2002). It is the confidence level with which a confidence region r_n may contain the true (mass, sma) and whose boundary passes through model n¹⁴. Details on how the $\ln \mathcal{L}$ is calculated for a given bin are explained in Sect. 6.3.1.

We thus produced the GaiaPMEX constraints map that shows the shape of the posterior probability density functions on the Log-mass and Log-sma (or Log-period) with respect to the α_UEVA, HG PMa (Kervella et al. 2022), or both combined. The combination of both observables is most important in enabling the exclusion of either small period or long period companions, as discussed in Sect. 7.

7 The GaiaPMEX maps

7.1 Constraints from AEN, RUWE, and PMa on the mass and sma of a companion

We use the systems of GJ 832 and HD 114762 to illustrate the outcomes of the GaiaPMEX approach. Details on HD 114762 and GJ 832 are given respectively in Sects. 5.2 and 5.2.2 and their main parameters are summarized in Table B.1. Figs. 19 and 20 show the confidence regions, or maps, at confidence levels 68.3 (1−σ region), 95.4 (2−σ region), 99.73% (3−σ region), of the mass and sma of a candidate companion around GJ 832 and HD 114762 as calculated with GaiaPMEX, given either only ruwe or AEN – through calculating UEVA as explained in Sect. 5.1.1 – or only PMa, or combining the constraints from PMa and ruwe.

The maps give the possible sma and mass of a companion assuming it is responsible of the measured PMa, AEN or ruwe. The mass-sma degeneracy drawn by the confidence regions follows typical curves, essentially U-shaped for the PMa, and V-shaped for the AEN or the ruwe. They present mainly three features: a short-period (SP) branch, a long-period (LP) branch, and the near equal-mass binary solutions forming an horizontal branch around the mass of the primary (B-branch hereafter). The SP and LP semi-linear branches curve up into the B-branch because of the contribution of the secondary in the position of the photocenter. The GaiaPMEX combination of the constraints from PMa and ruwe, as shown on the green map in Figs. 19 and 20, combines the shape of the PMa and ruwe maps. Most importantly it lifts up a fraction of the degeneracy, leaving confidence regions of smaller extent.

In the case of GJ 832, as shown in Fig. 19, the α_UEVA (AEN) maps form a confidence region at 2−σ leading to a mass >0.1 M_J with 95.4% confidence and a mass upper-limits at any sma <10 au. However, the ruwe is only leading to the 2−σ upper-limit on mass at any sma. This is explained by the difference in the significance of the α_UEVA from either ruwe (1.7−σ) or AEN (2.5−σ). The PMa maps lead to a well-defined 3−σ confidence region with mass >0.5 M_J at 99.7% confidence and an upper-limit on mass within 1−40 au. Both 1−σ confidence regions of the α_UEVA and PMa maps encircles the minimum mass 0.74 ± 0.06 M_J at sma=3.6 ± 0.4 au for GJ 832 b (Gorrini et al. 2022), as well as the mass of ${0.99}_{-0.08}^{+0.09}$ M_J found by combining the RV and the Hipparcos-Gaia proper motions of GJ 832 in Philipot et al. (2023a). Nonetheless, the confidence regions still leave a lot of degeneracy in the sma and mass solutions. Combining the constraints from PMa and ruwe leads to reject most of the SP-branch. The remaining LP-branch spans ∼ 2−500 au with mass as low as 0.5 M_J. The 1−σ, or 68%, confidence region is restricted to within the LP-branch with sma=2-8 au and true mass within 0.5−2 M_J. Again, the 1−σ confidence region encircles the known possible mass and sma of GJ 832 b.

In the case of HD 114762, as shown in Fig. 20, the confidence regions drawn from AEN, ruwe and PMa agreed well on a companion mass >20 M_J for sma within 0.1−100 au. The ruwe and AEN even constrained sma to be lower than 20 au. Combining PMa and ruwe together led to still largely degenerated 2 and 3−σ confidence regions. The 1−σ region is now centered about 3 au and mass=50 M_J, and the 2−σ regions spans the stellar domain, embedding the known mass=225 ± 14 M_J at 0.38 au of HD 114762 b (Winn 2022). The wide-orbit companion at 130 au HD 114762 B cannot be compatible with neither Gaia nor Gaia-Hipparcos astrometry and thus did not contribute significantly to the photocenter’s motion as captured by Gaia.

The broken linear relationships drawn upon the U and V-shaped confidence regions in the PMa, AEN and ruwe maps are related to their physical interpretation. The semi-empirical laws behind those curves draw approximate constraints on the mass and sma for any system. We further elaborated on these insights in Sects. 7.2 and 7.3.

7.2 The mass-sma relationship constrained by the α_PMa

The typical broken U-shape of the mass-sma relationships that we found in the PMa map is shown in Fig. 19. It can be divided into four sections.

First, if P≫24.5 yr, then the HG baseline covers a part of the orbit. The contribution of the orbital motion in the HG PM vector and in the GDR3 PM vector are the mean orbital speed during respectively 24.5 and 3 yr. In this regime, the PMa is indeed an approximation of the orbital acceleration fixed by the relative sma of the companion, with PMa ∝ 1/sma². This approximation is increasingly accurate with increasing orbital period, as the orbital motion becomes more linear within 24.5 yr. If P<24.5 yr, for an orbital period sufficiently small, the contribution of the orbital motion in the 24.5-yr averaged PM becomes negligible.

In fact, the central point between those two main regimes is close to a period of 49 yr, that is, about P/2=24.5 yr. Indeed, for a star on a pole-on circular orbit with P/2=24.5 yr, the orbital contribution to the HG PM (PM_HG, orbit hereafter) that is, the average orbital speed along half the period, is ∼4 a_*/P; the orbital contribution to the GDR3 PM (PM_{GDR3, orbit} hereafter) is the instantaneous speed ∼ 2 π a_*/P; and they lead to a ratio PM_HG, orbit/PM_GDR3, orbit ∼0.6. For P/2<24.5 yr, the orbital motion contribution to the PM_HG, orbit becomes even smaller than half the orbital motion contribution to PM_GDR3, orbit. Conversely, for P/2>24.5 yr, the PM_HG, orbit and PM_GDR3, orbit strengths gradually reach equality with increasing P.

Below 49 yr, assuming that PM_HG, orbit becomes rapidly negligible, we find three different regimes that depend on GDR3 sampling and baselines at P<1 yr, 1<P<6 yr and 6<P<49 yr:

• P<1 yr: The Gaia data undersample the orbits, leading to strong aliasing of the orbital signal. The fit proper motion varies a lot from one value of sma to another, and almost nulling at resonances between the main sampling frequency and the orbital frequency (see also Kervella et al. 2019). There is nonetheless an average trend of the α_PMa that is linearly increasing with sma, and thus α_PMa ∝ sma;

• 1<P<6 yr: Single orbits are better phase-covered, and there is more than half an orbit cycle monitored during the 3 yr of the GDR3. In this regime, the smaller the period, the more cycles are monitored and the less sensitive the α_PMa is to the amplitude of the motion. Thus, α_PMa ∝ sma P ∝ sma^5/2;

• 6<P<49 yr: Less than half a phase of the orbital motion is covered by the Gaia observations and the fit proper motion is oriented along, and in intensity proportional to, the average orbital velocity of the star along this orbit segment. Thus, the α_PMa is approximately proportional to the instantaneous orbital speed, ${\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}\propto 1/\sqrt{\text{sma}}$.

Adding the correct dependence on the star mass M_*, the mass of the companion M_c, and the parallax ϖ, these relations lead to the four mass-sma log-linear relations observed in the PMa maps: $${\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}\propto \{\begin{array}{ll}{M}_c{M}_{\star }^{-1}\varpi \mathrm{s}\mathrm{m}\mathrm{a}& \text{in}(1):<1\mathrm{y}\mathrm{r}\Rightarrow M_c\propto {\mathrm{s}\mathrm{m}\mathrm{a}}^{-1}\\ M_c{M}_{\star }^{-3/2}\varpi {\mathrm{s}\mathrm{m}\mathrm{a}}^{5/2}& \text{in}(2):1-6\mathrm{y}\mathrm{r}\Rightarrow M_c\propto {\mathrm{s}\mathrm{m}\mathrm{a}}^{-5/2}\\ M_c{M}_{\star }^{-1/2}\varpi {\mathrm{s}\mathrm{m}\mathrm{a}}^{-1/2}& \text{in}(3):6-49\mathrm{y}\mathrm{r}\Rightarrow M_c\propto {\mathrm{s}\mathrm{m}\mathrm{a}}^{1/2}\\ M_c\varpi {\mathrm{s}\mathrm{m}\mathrm{a}}^{-2}& \text{in}(4):>49\mathrm{y}\mathrm{r}\Rightarrow M_c\propto {\mathrm{s}\mathrm{m}\mathrm{a}}^2\end{array}$$(31)

Those four regimes are all linear in log-log space but with different slopes. We note that the constraints stop following these linear models as the luminosity of the companion becomes comparable with the luminosity of the primary. When this occurs, as the luminosity of the companion increases, the photocenter semi-major axis shrinks more and more until reaching almost zero at a mass ratio¹⁵ q ∼ 1.

7.3 The mass-sma relationship constrained by the α_UEVA

We now explain the typical V-shape of the mass-sma relationship that we observed in the UEVA maps, that is, those of AEN and ruwe, shown in Fig. 19. Two regimes need to be distinguished. The regimes with P either shorter than 3 yr, when more than one orbital cycle is monitored, or longer than 3 yr, when only a part of a single orbit is covered by Gaia.

• If P<3 yr, several orbital cycles happen during the GDR3 observation campaign. The residuals amplitude, and thus α_UEVA, varies in proportion to the extent of the motion. The star’s semi-major axis a_* itself is proportional to sma, α_UEVA ∝ sma.

• If P>3 yr, only a part of a single orbit is covered by GDR3’s observations and is locally fit by a linear motion. The α_UEVA measures the non-linearity of the motion, that is, the acceleration. Therefore, in this regime, α_UEVA ∝ 1/sma².

Adding the correct dependence on the star mass M_*, the mass of the companion M_c and the parallax ϖ, these relations lead to the mass-sma log-linear relations observed in the UEVA maps: $${\alpha }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}\propto \{\begin{array}{ll}{M}_c{M}_{\star }^{-1}\varpi \mathrm{s}\mathrm{m}\mathrm{a}& \text{in}(1):<3\mathrm{y}\mathrm{r}\Rightarrow M_c\propto {\mathrm{s}\mathrm{m}\mathrm{a}}^{-1}\\ M_c\varpi {\mathrm{s}\mathrm{m}\mathrm{a}}^{-2}& \text{in}(2):>3\mathrm{y}\mathrm{r}\Rightarrow M_c\propto {\mathrm{s}\mathrm{m}\mathrm{a}}^2.\end{array}$$(32)

The log(mass) decreases linearly with the log(sma) down to about P=3 yr that corresponds to a range of sma within 1–3 au depending on the stellar mass. Then the log(mass) increases linearly with the log(sma) up to very long period, until it reaches stellar masses. We note that, as for the PMa, the constraints stop following these linear models as the luminosity of the companion becomes comparable with the luminosity of the primary, that is, beyond M_c∼0.5 M_*.

7.4 The mass-sma relationships and the minimum mass of companion

The relations expressed above between M_c and sma, at given ϖ and M_*, can be summarized as $$M_c=C_{\ell }\frac{\left({\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}\right)\text{or}\left({\alpha }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}\right)}{\varpi }{M}_{\star }^{\frac{2-n_{\ell }}{3}}{\text{sma}}^{n_{\ell }}$$(33)

with ℓ a number indexing the considered period regime, C_ℓ a multiplicative constant, and n_ℓ an exponent. To determine the C_ℓ’s, we considered an arbitrary source and modeled by simulation 1000 values of α_UEVA and α_PMa at a P of the corresponding regime – either 0.1, 2, 30 or 100 yr – and at fixed M_c=10 M_J, M_*=1 M_⊙ and ϖ=1000 mas. The inclination, eccentricity and ω were set randomly following the distributions of Table 2. This led, for given sma=P^2/3 $M_{\star }^{1/3},M_c,M_{\star }$, M_c, M_* and ϖ, to distributions of α_UEVA and α_PMa, and thus, inverting Eq. (33), to measurements of the constants, C_ℓ,simu ± σ_C,ℓ. Because this might lead to disconnected line segments at the transition between regimes, that is, at P=1, 3, 6, or 49 yr, we also imposed a continuity condition by noticing that from one regime ℓ to another contiguous ℓ + 1, there is a simple condition to fulfill: $$\log C_{\ell +1}-\log C_{\ell }=\frac{2}{3}(n_{\ell }-n_{\ell +1})\log P.$$(34)

To do so, we minimized a least-squares problem by varying for either the case of α_PMa or α_UEVA only C₁ and minimizing the following objective function f: $$f=\sum _{\ell }\frac{{(C_{\ell }-C_{\ell ,\text{simu}})}^2}{{\sigma }_{C,\ell }^2}$$(35)

where any C_ℓ≠1 is determined by using Eq. (34) and, for example, log C₃−log C₁=(log C₃−log C₂) + (log C₂−log C₁). We expected that the C_ℓ determined by this method would depend on the source considered due to the source-dependent noise levels, the shape of the unit parallax ellipse, and the Gaia and Hipparcos scan law. We thus used 20 arbitrary sources with diverse magnitudes, colors, and RA-DEC and determined their global average C_ℓ rounded at the significant digit, given standard deviations of about 10−20%: $${\alpha }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}⟹\{\begin{array}{ll}\ell =1:<3\mathrm{y}\mathrm{r},& n_1=-1,C_1=2300\\ \ell =2:>3\mathrm{y}\mathrm{r},& n_2=+2,C_2=260\end{array}$$(36) $${\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}⟹\{\begin{array}{ll}\ell =1:<1\mathrm{y}\mathrm{r},& n_1=-1,C_1=6800\\ \ell =2:1-6\mathrm{y}\mathrm{r},& n_2=-5/2,C_2=6800\\ \ell =3:6-49\mathrm{y}\mathrm{r},& n_3=+1/2,C_3=190\\ \ell =4:>49\mathrm{y}\mathrm{r},& n_4=+2,C_4=3.8\end{array}.$$(37)

Reporting these numbers in Eq. (33), the curve segments in these different regimes are added to the GaiaPMEX maps in Figs. 19 and 20. Given the actual mass uncertainty associated with the maps at any sma, the C_ℓ given here should only be considered as indicative, which is why we did not quote their uncertainties. These curves give an immediate approximate idea of what GaiaPMEX maps should look like without doing the full computation. They are not intended as a substitute for GaiaPMEX calculations. They do not produce, conversely to GaiaPMEX, exact constraints on the sma and mass of the companion given the observed AEN, ruwe and PMa. Nevertheless, these curves can allow for rapid estimation of the mass of the companion as a function of its possible sma given either AEN, ruwe, or PMa.

If the α_UEVA or the α_PMa are sufficiently significant, typically more than 2−σ, these relationships enable measuring of the minimum mass of the companion, located at the minimum of all curves. This is reached at an sma∼2–3 au: $${\alpha }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}⟹\{\begin{array}{l}{M}_{c,min}=1150M_{\star }^{2/3}\left({\alpha }_{\mathrm{U}\mathrm{E}\mathrm{V}\mathrm{A}}/\varpi \right)\\ {\mathrm{s}\mathrm{m}\mathrm{a}}_{min}=2.1M_{\star }^{1/3}(\mathrm{a}\mathrm{u})\end{array}$$(38) $$\begin{array}{}\text{(J)}& {\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}⟹\{\begin{array}{l}{M}_{c,min}=340M_{\star }^{2/3}\left({\alpha }_{\mathrm{P}\mathrm{M}\mathrm{a}}/\varpi \right)\\ {\mathrm{s}\mathrm{m}\mathrm{a}}_{min}=2.9M_{\star }^{1/3}(\mathrm{a}\mathrm{u}).\end{array}\end{array}$$(39)

These values can be used to identify planet candidates among catalogs of targets observed with Gaia.

8 Illustrative cases of GaiaPMEX results

8.1 HD 81040: Evidence for a short-period companion

HD 81040 is well-known as being the system of the first discovery from GDR3 data of an exoplanet companion using only an orbital fit of the astrometric data (Gaia Collaboration 2023a)¹⁶. HD 81040 is a solar-type G2/3V star of magnitudes V=7.72 and G=7.57, at a distance of 32.56 ± 1.31 pc from Earth. The planetary candidate HD 81040 b was first discovered by Sozzetti et al. (2006) with an orbital period of 1001 days, eccentricity of 0.53 and an msin i of 6.86 M_J. Using the PMa calculated from the Gaia EDR3 and Hipparcos data, there was further indications that this companion was indeed planetary with a mass of ${7.24}_{-0.37}^{+1.00}$ M_J (Li et al. 2021; Winn 2022). Ultimately, a joint fit of the astrometric GDR3 data and the RV data for this star fully confirmed a mass of ${8.04}_{-0.54}^{+0.66}$ M_J. At the orbital period, best constrained with RV, of 1001 days, the planet is predicted to have an sma=1.94±0.02 au (ref. Gaia ESA webpage).

Figure 21 shows the map constraining the mass and sma of a candidate companions around HD 81040 calculated with GaiaPMEX from combining ruwe and PMa constraints. Other individual maps using only AEN, ruwe and PMa are shown in the Appendix F, Fig. F.1. The parameters and results found for HD 81040 are summarized in Table B.1. The astrometric signatures of AEN and ruwe are both significant at >6.7−σ, with p-value∼0.00, while the α_PMa cannot reject the single star hypothesis at a 1.5−σ significance and a p-value of 0.13. Nevertheless, the constraints on mass-sma calculated by GaiaPMEX from the PMa stay compatible with the known planet in this system, but do not strongly reject an edge-on inclination.

Most interestingly, when combining the constraints from the ruwe and the PMa, the LP-branch disappeared from the GaiaPMEX map, only leaving an SP-branch. The case of HD 81040 is thus illustrative of the identification by GaiaPMEX of a short-period companion, typically detectable using RV. Conversely, in such a case, it may allow for putting constraints on the mass of a known RV companion, as done with GASTON in Kiefer et al. (2019); Kiefer (2019); Kiefer et al. (2021). Here, the GaiaPMEX maps agree well with the mass and sma of HD 81040 b found by fitting the astrometric time series and the RV.

8.2 AF Lep: Evidence for a long-period companion

AF Lep is an F8V star, initially characterized as an RS CVn. It has been shown that its variability is actually due to the presence of a companion, a super-Jupiter with mass ∼3 M_J at 8 au, detected through the combination of astrometric acceleration and direct imaging (Franson et al. 2023; De Rosa et al. 2023; Mesa et al. 2023). The main parameters of AF Lep are summarized in Table B.1. This system was selected for follow-up with direct imaging because it shows a significant astrometric acceleration, or PMa, compatible with a planet mass. This was a long-shot, because using PMa only, one cannot reject that this system is actually a short period binary that would remain unresolved by direct imaging, while lower mass or lower sma could remain plausible, also leading to a non-detection with direct imaging.

The GaiaPMEX maps obtained for AFLep are shown in Figs. 22 and F.2. In agreement with previous work on this source, the astrometric signature of the PMa is significant at 4−σ, that is, compatible with the single star hypothesis with a p-value of 6.334 × 10⁻⁵. Both astrometric signatures found for AEN and ruwe cannot reject the single star hypothesis for AFLep with a p-value of 0.96 and a significance <0.1−σ. It leads to strict upper-bounds on the companion mass at any sma, since otherwise the AEN and ruwe would have been more significant.

Interestingly, the mass-sma constraints from PMa are largely degenerate, but combining them with those from ruwe leads to rejection of most low-sma solutions, that is, the SP-branch.

This shows that combining the PMa and the ruwe would have led to much tighter constraints on the possible mass and sma of a companion around AFLep. The detected companion with ${4.3}_{-1.2}^{+2.9}$ M_J at sma=7.99+0.85 au using direct imaging (Mesa et al. 2023; Franson et al. 2023) falls indeed within the 1−σ bounds of the LP-branch.

8.3 HD 23596: An almost detection with astrometry

HD 23596 is a V=7.2-mag F8 star at 52 pc from the Sun that is known to host a companion planet. It was first discovered as an 8.2−M_J super-Jupiter with the ELODIE spectrograph (Perrier et al. 2003) and further characterized with a similar msini by Wittenmyer et al. (2009) (7.71±0.39 M_j) and Stassun et al. (2017) (9.03±0.74 M_J) using supplementary High Resolution Spectrograph (HRS) data and new M_* estimations. It was then re-established as a 14−M_J low-mass brown dwarf combining RVs and Hipparcos-Gaia PMa (Feng et al. 2022; Xiao et al. 2023). The orbital period of this companion is ${4.31}_{-0.055}^{+0.069}$ years, with an sma of 2.90 ± 0.08 au, a mass of ${14.6}_{-1.3}^{+1.5}$ M_J and an inclination being of either ${34.0}_{-2.9}^{+3.6\circ }$ (prograde) or ${146.0}_{-3.6}^{+2.9\circ }$ (retrograde) (Xiao et al. 2023). The inclination significantly non edge-on explains the lower 8.2-MJ m sin (i) initially found for this companion.

The GaiaPMEX analysis combining the constraints from both ruwe and PMa of HD 23596 led to the confidence regions on companion mass and sma shown in Fig. 23. The individual maps from either ruwe, AEN or PMa contraints are shown in the Appendix in Fig. F.3. The combination of PMa and ruwe leads to infer a companion in the brown dwarf domain, with a narrow constraint on mass within 10−30 M_J as well as on sma within 2–5 au at 68.3% confidence. This is in perfect agreement with the known companion of HD 23596. Surprisingly, this source was not identified as a non-single star and does not appear in the non-single star catalog (Gaia Collaboration 2022)¹⁷. This shows that Gaia in combination with Hipparcos, with the help of the GaiaPMEX tool, can detect and characterise a companion even without including RVs. More such planet candidates with strong constraints on the mass and the sma from only Gaia + Hipparcos astrometry will be presented in other forthcoming papers (paper II, Kiefer et al. 2025; paper III, Lagrange et al. 2024).

8.4 α CMa B: The effect of considering a dark companion

The case of α CMaB is an instructive opportunity to show the effect of considering either a dark or a non-dark companion in the GaiaPMEX analysis. Indeed, α CMa B is nothing else than Sirius B, the white dwarf companion of the massive 240 Myr old Sirius. It was first hypothesized by Bessel (1844), then officially discovered by A. G. Clark in 1862, and further characterized by Flammarion (1877). Even though its average separation with Sirius A is only 19.6 au, the close distance to the Sun of this system lead to an angular separation of 7.6^′′ (van den Bos 1960), and thus Sirius B was resolved in the GDR3. The respective mass of Sirius A & B are 2.14 M_⊙ and 1.05 M_⊙ (Gatewood & Gatewood 1978). Based on the orbital analysis of the Sirius AB system using astrometry, an additional companion Sirius C around either A or B was suspected to exist in this system, with a putative orbital period of 6 years (Volet 1932; Benest & Duvent 1995). Around Sirius B’s, this companion may have a mass <20 M_J and an orbit sma of 1–2.5 au (Bonnet-Bidaud & Pantin 2008). Using Keck/NIRC2 observations, Lucas et al. (2022) further excluded any companion of mass >10 M_J down to 0.2 au, >2.4 M_J down to 0.5 au and >0.7−1.2 M_J beyond 1 au around Sirius B.

In GDR3, as summarized in Table B.1, the α_UEVA, from both AEN and ruwe, of Sirius B is significant, >9−σ. Sirius B was not observed with Hipparcos, implying no PMa for this star. The map derived from ruwe is shown in Fig. 24. We considered, by default in the top figure, a luminous companion. Within Lucas et al. (2022) constraints, the ruwe allows at 2−σ for a companion with a mass as high as 10 M_J at less than 0.5 au. However, with a mass ∼2.1 M_⊙, the known companion of Sirius B, that is, Sirius A, at an average sma∼20 au, seems to explain the observed ruwe without the need to invoke a supplementary companion.

Nonetheless, for Sirius A, this approach is problematic for two reasons. First, the age of the system is overestimated since in the current version of GaiaPMEX, the mass-luminosity relation is determined only for an age of 5 Gyr. Second, and most importantly, even though it is much more luminous than B, Sirius A is well resolved by Gaia and cannot act in the position of the photocenter of Sirius B. The main effect of considering dark or luminous companion is to change the shape of the confidence region at companion mass on the order of the magnitude of the mass of the source, here 1.05 M_⊙. A luminous companion narrows down the photocenter’s orbit. For a given constant photocenter semi-major axis, at short periods it requires the companion’s sma to increase, while at long periods it requires the acceleration to increase and thus the sma to decrease. Considering instead a companion whose light does not contribute to the photocenter (dark or resolved), the confidence region rather more closely follows the curves derived in Sects. 7.3 and 7.2 that only describe the reflex astrometric motion of the main source under the gravitational pull of a companion.

The GaiaPMEX map derived from ruwe for Sirius B and considering a dark companion is shown in Fig. 24, bottom. We note that Sirius A is within the 3−σ region and upon the edge of the 2−σ region. Thus, we confirm that Sirius A may indeed explain the ruwe, although it would tend to generate, on average, a ruwe larger than the one published in the archives.

Finally, the GaiaPMEX maps of Sirius B cannot exclude any companion with a mass located below the confidence regions, since the perturbation from this other companion would be subdominant. It thus remains possible that another hidden companion exist around Sirius B. The IPD_frac_multi_peak and IPG_gof_harmonic_ampl of α CMaB are moderate, respectively 18% and 0.05, but according to our analysis performed in Appendix C, they imply that a source within 200–500 mas, that is, 0.5−1.3 au, with a flux ratio within 10⁻⁴−10⁻³ is acting on the shape of the PSF. Interestingly, the Lucas et al. (2022) analysis allows for a <1−2 M_J companion at 0.5−1.3 au below a limiting contrast of ∼10⁻³. This strongly suggest digging further within the 1-au surroundings of the white dwarf Sirius B in the quest of planets.

8.5 β Pictoris: An issue with the RUWE

We focus here on the system of β Pictoris to bring forward again that the ruwe and AEN cannot be interpreted directly as binarity indicators, but only through determining the astrometric signature within residuals and its significance. The main parameters of β Pictoris are summarized in Table B.1. β Pictoris is a south-hemisphere 20-Myr old A5V star (Mamajek & Bell 2014) located at 19.6pc from the Sun with a G-mag of 3.82. It is known for being the most furbished planetary system, after the Solar system, in the wealth of body types that planet formation can produce. It hosts a widely extended dusty and gaseous debris disk, exocomets, asteroids, and giant exoplanets with masses of ∼8 and ∼10 M_J. It is a system of great importance and a focus of attention from the community for understanding the running processes during the first hundred Myr of the Solar system and of planetary systems in general. In the GDR3 database, β Pictoris is announced with a ruwe of 3.07 and an AEN of 1.39 mas. It looks therefore at first sight that significant deviation to the five-parameter model has been detected in this system, if we follow the guidelines that ruwe>1.4 indicates a non well-behaved or non-single system (as recommended in the DR3 documentation). However, here it is not the case.

At a ${\chi }_{\text{astro}}^2=66,641.58$ for 231 good AL measurements, the UWE for β Pic is 17.2 and seems anomalously large, apparently corroborating the large ruwe. Considering the level of noises for this bright blue source of the 6p-dataset that we found in Sect. 3, that are σ_calib=1.548 mas, σ_AL=0.012 mas, and σ_AL=0.074 mas, the χ² expected for a single source according to Eqs. D.8 and D.9 is 80 000±24 000. And it corresponds to UWE=18.8 ± 2.8, in good agreement with the value found above, indicating that 17.2 is, in fact, not an anomalously large UWE. It is therefore surprising that the ruwe of β Pic is 3 times larger than 1.0, while the ruwe was introduced as a renormalized version of the UWE to recenter this goodness-of-fit indicator around unity for well-behaved sources. The renormalizing factor u₀ in Eq. (3) is thus ill-defined for sources in the 6p-dataset at the (G, Bp − Rp) corresponding to β Pic, that is, at G=3.82 and Bp − Rp=0.261. The AEN=1.39 mas is conversely well defined for β Pic. According to the relation of Eq. (7), the AEN implies ${\chi }_{\text{astro}}^2=77922$. This is in good agreement with the value expected for a single star determined above. The astrometric signature derived from AEN and ruwe is undefined because their corresponding UEVA are smaller than the single star’s expected UEVA_single in Eq. (27).

Figure 25 compares the UWE obtained for single stars from Eq. (D.8) to the u₀ found in the GDR3 auxiliary data¹⁸ over the whole sources database with G<16. The ratio of both quantities is ruwe_single, that should ideally be equal to 1. The distribution of this ratio peaks indeed at 1, but in the 6p-dataset, there are strong tails on both sides toward lower and higher values. In the 5p-dataset, some ruwe_single deviate from 1, although much tightly than in the 6p-dataset. In the 5p dataset, 462 over 71 042 992 sources (6.510⁻⁴%) have a ruwe >1.4, while in the 6p dataset there are 14 890 over 2 203 807 sources (0.68%) with a ruwe >1.4. In conclusion, regardless of the dataset, it is safer to interpret the values of AEN and ruwe only through calculating the astrometric signature of the UEVA, as defined in Sect. 5.1.

The GaiaPMEX maps of β Pictoris are shown in Fig. 26. Consistently with the above analysis, the AEN and ruwe both lead to sma and mass of a hypothetical companion to β Pic compatible with a single star, that is, mass=0 M_J or sma=0 au. The two planets (b: 10−11 M_J at 9.8±0.4 au; c: 7.8±0.4 M_J at 2.7±0.02 au; Lagrange et al. 2020) are compatible with the AEN and ruwe, as well as the PMa. There are no evidence in the GDR3 for any other yet unknown companion around β Pictoris.

9 The current sensitivity of Gaia to the presence of companions

We determined the detection rates of exoplanets of different mass and sma that one can expect from using either PMa, UEVA (determined from AEN or ruwe), or PMa ∪ UEVA¹⁹. They strongly depend on the stellar mass M_* and parallax ϖ. We call “detection” any value of α_UEVA or α_PMa more significant than 2−σ, that is, for which the UEVA^1/3 or the PMa^2/3 are above the 95.4th percentile of the respective single star’s distributions.

Considering pre-main to main-sequence (MS) stars, we explored a 2D-grid, with 30 bins per dimension uniformly spaced in log-scale, of M_* from 0.08 to 2.5 M_⊙ and ϖ from 1 to 1000 mas. At each bin with given (M_* ± ΔM_*, ϖ ± Δϖ), we modeled UEVA and PMa by simulation (see Sect. 4) of photocentric orbits as observed by Gaia and Hipparcos, due to companions with a mass M_c within bins delimited by 0.1, 0.2, 0.51, 2, 5, 10 and 20 M_J and with an sma within different orbital regimes: Mercury-Earth type (0.1−1 au), Earth–Mars type (1−3 au), Jupiter-Saturn type (3−10 au) and Uranus-Neptune type (10−30 au). The G-mag and the Bp − Rp color corresponding to a given M_* and ϖ are calculated from Pecaut et al. (2012) and Pecaut & Mamajek (2013)’s spectral type to flux conversion tables for pre-MS to MS stars²⁰. The absolute magnitude M_G is converted to apparent magnitude G using the distance modulus calculated from the parallax. For simplicity, and to draw the general picture, we assumed zero extinction. The results are thus susceptible to be only informative, especially beyond 100 pc. We fixed the other properties (such as RA, Dec, parallax unit ellipse, Hipparcos IADs) and the GDR3 epochs sampling, to those of GJ 832. At each bin, the noises σ_calib and σ_AL levels are fixed with respect to the median G and Bp − Rp. Since our initial sample covers 3 to 16 G-mags, our study of sensitivity in GDR3 is limited to this range.

For each bin, we simulated N=1000 orbits given the range of companion sma and M_c and host star’s M_* and ϖ, randomizing other orbital parameters according to the distributions defined in Table 1. We counted the percentage r – or detection rate – of UEVA^1/3 and PMa^2/3 exceeding the 95.4th-percentile in their respective single stars’ distributions. We determined when r>20, 50, 90 and 99%. We were peculiarly interested in the possibility to detect a planet by considering PMa, UEVA or both, that is, PMa ⋃ UEVA. For single stars, the frequency of FP beyond the 2−σ threshold is 4.6%, when PMa or UEVA are considered separately, and 9.2%, when considering PMa ⋃ UEVA. This is the worst case scenario (largest FP frequency) in which noise rather than orbital motion causes a significance larger than 2−σ. In the general case, some true companion detections might be serendipitous, that is, due to noise rather than orbital motion, and it is not possible to determine exactly the fraction of FP in this case. At best, a percentage r>4.6% (respectively 9.2%) indicates an increased sensitivity of Gaia to the detection of companions in the given range of mass and sma. At worst, the fraction of FP is 4.6/r (respectively 9.2/r). In particular, if r is close to 4.6% (respectively 9.2%), then the fraction of FP is close to 100%.

A map of the detection rates – or equivalently, Gaia's sensitivity to companion detection – with respect to star mass and parallax from using PMa, α_UEVA or both is shown in Fig. 27. Planets with mass <1 M_J around stars located farther than 100 pc (ϖ<10 mas) from the Sun lead to significant astrometric signal in less than 20% of the simulations. But Jupiter-mass planets (1−2 M_J) might be detected with a >20% chance around stars less massive than 1 M_⊙ and up to 100 pc distance. Similarly, planets of mass <0.1 M_J around solar-like stars (0.5−2 M_J) have a less than 20% chance of detection with Gaia whatever their sma and whatever the distance to the Sun. But, if their host is an M-type star closer than 10 pc and if their sma is within 1−10 au, planets of any mass >0.1 M_J have a more than 20% chance of being detected at 2−σ. Moreover, within 1−10 au and if their host star is an M-type star closer than 5pc, Jupiter-mass planets (1−2 M_J) were detected in >99% of simulations, and Neptune/Saturn-mass planets (0.1−0.2 M_J) in >50% of the simulations. Finally, super-Jupiter and brown dwarfs are easily detected at a rate >99% up to large distance, even beyond 100 pc and around A-type stars with M_*>2 M_⊙.

Focusing now on close-by low mass stars, we invoked again GJ 832, a 5-pc distant M-dwarf, and derived the maps of Gaia's sensitivity to companion detection around this star. We modeled the UEVA and PMa on a grid of sma and M_c, the same that was used for the GaiaPMEX constraints maps, and determined the detection rates r, as done above. The map for detection with UEVA ⋃ PMa is shown in Fig. 28 and the individuals maps for detection with either UEVA or PMa are shown in Fig. G.3. Around such a star, super-Jupiter and brown dwarfs with mass >10 M_J and with an sma within 0.2–20 au are detected in >99% of the simulations. Lighter planetary companions whose mass is within 2−10 M_J may lead to a significant astrometric signal in >99% of the cases, provided their sma is contained within a narrower range of 1−10 au. At masses in the Saturn-to-Jupiter regime (0.2−1 M_J) planetary companions have a lower detection rate, with a >50% chance of being detected for an sma within 2−10 au. Planets of mass <0.2 M_J are much less frequently detected (<50%) whatever their sma. The detection rates obtained when considering only PMa or only UEVA (Fig. G.3) show that the PMa is more sensitive to planets orbiting in the 2−20 au range, while the UEVA tends to perform the most efficiently at shorter separations from 0.2 to 2 au.

In summary, GDR3 currently performs best at detecting sub-stellar companions that are more massive than Jupiter on Earth-to-Saturn orbits (1−10 au) around any star closer than 100 pc from the Sun.

10 Conclusion

We have introduced GaiaPMEX, a tool that allows for characterization of the possible mass and sma of companions to stars observed with Gaia. It uses the proper motion anomaly (PMa) and the excess of the five-parameter model residuals, the unbiased estimator of variance a posteriori (UEVA), indicated by the values of renormalized unit weight error (ruwe) or astrometric excess noise (AEN). GaiaPMEX determines their significance with respect to the null hypothesis of a single star by calculating an astrometric signature for each indicator, α_PMa and α_UEVA, and then models them by the star’s reflex motion due to a companion with ranges of possible mass and sma.

Being otherwise highly degenerate in mass and sma, by combining the use of AEN or ruwe and PMa, the space of solutions is significantly reduced. This has never been done before and is highly promising for characterizing systems also observed by other means (radial velocities mostly). When α_UEVA or α_PMa are significant, the posterior maps follow mass-sma relationships that are given in Eqs. (36) and (37) and that vary with orbital periods. Thereby, the PMa could only be considered as a proxy for orbital acceleration when the orbital period of the companion is longer than 49 yr, while the UEVA varies positively with the photocenter semi-major axis only when the orbital period of the companion is shorter than 3 yr.

We determined companion detection rates and Gaia's sensitivity to exoplanet detection by comparing modeled UEVA and PMa for companions with a hypothesized mass and sma to expectations for single stars. We show that in the exoplanet domain, the best detection rates are obtained for M-dwarf sources that are closest to the Sun (<10 pc) and for exoplanet’s sma within Earth-to-Saturn orbits, that is, 1−10 au with masses down to 0.1 M_J.

In the present article, we introduced in greater detail the GaiaPMEX tool and the needed estimation of calibration, attitude and measurement noise levels. We plan on extensively exploiting GaiaPMEX for the search and characterization of exoplanets in further studies.

Data availability

Appendix E–H are available at the following url at https://zenodo.org/records/17192804?token=eyJhbGciOiJIUzUxMiJ9.eyJpZCI6ImFkYWV1NDg1LTQzMmQtNDczZC1hNzk3LWYzYzUwNmYxZjE1MSIsImRhdGEiOnt9LCJyYW5kb20iOiIyNjRlN2QyYjczOWRhODVlMjc1NmU1NjlhNjQ2MmFlMyJ9.n8FLZELuZ80ifu1hqDkJ6svxBxeBxmsDauL3LmlGmDVv50Jk3FxKYF4KuYMO_2ExMNdHpB4LkSpvN_urjrCLuw

Acknowledgements

We are very thankful to the anonymous referee for her/his thorough and courageous reading that led to significant improvements of this article. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data processing and analysis consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project EquipMeso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COBREX; grant agreement no. 885593). F.K. acknowledges funding from the initiative de recherches interdisciplinaires et stratégiques (IRIS) of Université PSL “Origines et Conditions d’Apparition de la Vie (OCAV)”, as well as from the Action Pluriannuelle Incitative Exoplanètes from the Observatoire de Paris – Université PSL. F.K. also acknowledges funding from the American University of Paris.

Appendix A Table of acronyms used in the text with their definitions and page references

Notation	Description	Page List
5p	five parameters	4, 6, 7, 10, 11, 22, 29, 30
6p	six parameters	4, 6, 7, 10, 11, 22
AC	across scan	2, 7
AEN	astrometric excess noise	1–4, 6, 9–14, 16–22, 25
AL	along scan	2–4, 6–10, 22, 29, 30
BD	brown dwarf	1
CU8	Coordination unit 8	13
DEC	declination	2–5, 8–11, 23
DOF	degree of freedom	2, 15
DPAC	Data processing and analysis consortium	25
DR	data release	1, 2, 6
FLAME	final luminosity age and mass estimate	13
FP	false-positives	12, 23, 24
FoV	field of view	6–11, 30
GASTON	Gaia astrometric noise simulation to derive orbital inclination	6, 20
GDR3	third Gaia data release	1–13, 15–19, 21–24, 28–30
GOST	Gaia observation forecast tool	7, 8
GaiaPMEX	Gaia DR3 proper motion anomaly and astrometric noise excess	1–3, 14–17, 19–25
HG	HIPPARCOS–Gaia	3, 8, 9, 12, 14–17
IADs	intermediate astrometric data	8
ICRS	International celestial reference system	3,4
IPD	Image parameter determination	29
LR	likelihood ratio	15
LSF	line spread function	29
MAD	median absolute deviation	6
MCMC	Markov-chains Monte-Carlo	13
MS	main-sequence	22, 23
PA	position angle	3,7,8
PDF	probability density function	14
PMa	proper motion anomaly	1–3, 6, 8, 9, 11–25
PSF	point spread function	4, 22, 29
RA	right ascension	2–5, 81, 23
RMS	root mean square	29
RSE	regression standard error	9, 10, 29
RSS	residuals sum of square	2
ruwe	renormalized unit weight error	1–3, 6, 9–14, 16–22, 24, 25, 29
RV	radial velocity	1, 16, 17, 19, 20
S/N	signal-to-noise ratio	13
sma	semi-major axis	1–3, 6, 12–25
UEVA	unbiased estimator of variance a posteriori	9–19, 21–25, 29
UWE	unit weight error	6, 22
WC	window class	4
WFS	wavefront sensor	7, 8

Appendix B Additional table

Appendix C The Gaia point and line spread function variations

The Image parameter determination (IPD) of Gaia makes a fit by a single point spread function (PSF; line spread function or LSF, if G>13) of the flux distribution within some defined window around any source transiting the detector. When polluting light is present and not masked out, it periodically affects the location of the measured photocenter compared to the theoretical photocenter at the light-barycenter of the main targeted source. The presence of unresolved sources (background or wide-orbit companions) nearby may thus cause undesirable variable shifts of the photocenter of the main source. The amplitude of the shifts is correlated with the angle between the AL scan direction (with position angle ψ) and the lines joining the polluting sources and the main source (with position angle θ).

The IPD_frac_multi_peak indicator published in the GDR3 archive gives the information on the fraction of exposures for which multiple modes have been detected. The window is then recalculated by masking out the parts containing the identified peaks (Holl et al. 2023). In theory, since some peaks may not be detected and masked out at all epochs, this can be a source of supplementary spurious variations in the 5p-fit residuals. These spurious variations are partly removed by the Gaia's reduction software. It fits a sinusoidal function of the angle of the AL scan direction ψ (a.k.a subpixel phase; Lindegren et al. 2021) to the astrometric time series of the targeted source (c + a sin ψ + b cos ψ; Eq. 9 in Lindegren et al. 2021). The shift function is not a perfect sinusoid due to the non-axisymetric shape of the PSF on the detector (Holl et al. 2023). Thus, a residual AL angle Δ η remains that may vary through time and mimick the signature of true orbital astrometric motion. The IPD produced another important indicator of PSF distortion, the IPD_gof_harmonic_amplitude. It measures the amplitude of the sinusoidal variation of the goodness-of-fit $(\text{GOF}=\ln {\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{e}\mathrm{d}}^2)$ of the PSF (LSF) fit. Fabricius et al. (2021) showed that sources with a ruwe indicative of possible binarity, and IPD_frac_multi_peak>2 or the IPD_gof_harmonic_amplitude>0.1, must be considered as resolved doubles.

We modeled the observation by Gaia, through time, of two Gaussian PSF of a wide-orbit binary, that is, with virtually no orbital motion during Gaia's 3-yr monitoring. We assumed separations ranging from 0 to 1000 mas and flux ratio of 0.001−1. An example is shown in Fig. C.1. We fixed the CCD noise to 1% of the flux, corresponding to AL measurement uncertainty of ∼0.03 mas. We arbitrarily fixed the flux at the tip of the PSF of the main source to 1. We fit the double star PSF by a single 2D-Gaussian varying only the scale and the centroid locus. A fit result is also shown in red in Fig. C.1. We obtained series of GOF and photocenter AL angle that we fit by a sinusoidal function following the AL scan rotation law, leading to GOF variation amplitude and residual Δ η beyond the five-parameter model. Figure C.2 shows the root mean square (RMS) of Δ η with respect to separation and flux ratio, as well as some isocontours of IPD_frac_multi_peak (at 2 and 98%) and IPD_gof_harmonic_amplitude (at 0.1 and 1%).

We found that significant residuals are associated with IPD_frac_multi_peak being different than 0 or 100% or with IPD_gof_harmonic_amplitude being larger than 0.1, while, insignificant variations on the order of the AL measurement uncertainty were mostly found when IPD_frac_multi_peak is close to 0 or 100% and IPD_gof_harmonic_amplitude<0.1. This agrees well with Fabricius et al. (2021) conclusions, but shows also that for fully resolved binary the IPD PSF fit remains unperturbed. As long as IPD_gof_harmonic_amplitude<0.1 and IPD_frac_multi_peak < 2 or >98%, the astrometric signature α_UEVA is only lightly affected by the binary astrometric bias and correctly measures supplementary astrometric motion. However, for the sources that do not comply to these conditions, the interpretation of α_UEVA, as indicative of orbital motion of the photocenter, is hazardous.

Appendix D Detailed calculation of the approximate probability density function followed by ${\chi }_{\text{astro}}^2$

We show here that for a single star, the ${\chi }_{\text{astro}}^2$ that can be found in the GDR3 archive does not actually follow a χ² distribution with N−5 degree of freedom. Instead we show that it more accurately follows a linear combination of normal distributions. The ${\chi }_{\text{astro}}^2$ is expressed with respect to the residuals of the five-parameter astrometric fit: $${\chi }_{\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{o}}^2=\sum _{\ell }^N\frac{R_{\ell }^2}{{\sigma }_{\ell }^2},$$(D.1)

where as assumed in the rest of this paper, the formal error σ_ℓ ≡ σ_f is approximately constant throughout all data points. The formal error is the quadratic sum ${\sigma }_f=\sqrt{{\sigma }_{\mathrm{A}\mathrm{L}}^2+{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2}$ (Lindegren et al. 2012) with σ_att the attitude excess noise (see Appendix 3.2) and σ_AL the along-scan measurement error (see Appendix 3.3). The formal error is explicitly calculated for each GDR3 source in Appendix 3.1. One residual R_ℓ is the sum of two contributions r_i and ${\xi }_j^{(i)}$, with r_i that varies randomly from one FoV transit (i) to another with normal distribution $\mathcal{N}$(0, σ_calib), and ${\xi }_j^{(i)}$ that varies for a given FoV transit (i) from one AL measurement (j) to another with normal distribution $\mathcal{N}(0,{\sigma }_{\mathrm{A}\mathrm{L}}^2)$. It follows that $$\sum _{\ell }^N{R}_{\ell }^2\approx N_{\text{AL}}\sum _i^{N_{\text{FoV}}}{r}_i^2+\sum _{\ell }^N{\xi }_{\ell }^2+2\sum _i^{N_{\text{FoV}}}\left(r_i\sum _j^{N_{\text{AL}}}{\xi }_j^{(i)}\right).$$(D.2)

We recall that N_FoV is the number of FoV transits on the detector, and N_AL is the average number of AL angle measurements per transit, that is, N_AL=int(N/N_FoV), given N is the astrometric_n_good_obs_AL. And thus ${\chi }_{\text{astro}}^2$ is the combination of three terms: $${\chi }_{\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{o}}^2\approx N_{AL}\frac{{\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}^2}{{\sigma }_f^2}X+\frac{{\sigma }_{\mathrm{A}\mathrm{L}}^2}{{\sigma }_f^2}Y+\frac{2}{{\sigma }_f^2}Z,$$(D.3)

with $$X\sim {\chi }^2(N_{\text{FoV}}-5)\sim \mathcal{N}(N_{\text{FoV}}-5,\sqrt{2(N_{\text{FoV}}-5)})$$(D.4) $$Y\sim {\chi }^2(N)\sim \mathcal{N}(N,\sqrt{2N})$$(D.5) $$Z\sim \mathcal{N}(0,\sqrt{N}{\sigma }_{\mathrm{A}\mathrm{L}}{\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}).$$(D.6)

The last term Z is obtained as the sum of the product of two normally distributed variables: $$\begin{array}{rl}Z& =\sum _i^{N_{\text{FoV}}}\sum _j^{N_{AL}}{r}_i{\xi }_j^{(i)}\approx \sum _i^{N_{\text{FoV}}}{r}_i{W}_i\\ & \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{W}_i=\mathcal{N}(0,\sqrt{N_{AL}}{\sigma }_{\mathrm{A}\mathrm{L}})\mathrm{a}\mathrm{n}\mathrm{d}{r}_i=\mathcal{N}(0,{\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}).\end{array}$$(D.7)

The standard deviation of the product of two normally distributed variables centered on zero is the product of their standard deviation, and the standard deviation of the sum of normally distributed variables is the root sum square of their standard deviation. This leads us to the formula of Z expressed in Eq. D.6, given moreover that the total number of astrometric points is N ≈ N_AL N_FoV. The two first terms X and Y dominate the spread, thus ${\chi }_{\text{astro}}^2$ should mainly be distributed according to a skewed χ² law. Thus, ${\chi }_{\text{astro}}^2$ follows a quasi-normal distribution $\mathcal{N}$(μ, σ) with $$\mu =\frac{N_{AL}}{{\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2}[(N_{\text{FoV}}-5){\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}^2+N_{\text{FoV}}{\sigma }_{\mathrm{A}\mathrm{L}}^2]$$(D.8) $$\begin{array}{rl}{\sigma }^2=& \frac{2N_{AL}}{{({\sigma }_{\mathrm{a}\mathrm{t}\mathrm{t}}^2+{\sigma }_{\mathrm{A}\mathrm{L}}^2)}^2}[N_{AL}(N_{\text{FoV}}-5){\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}^4\\ & +N_{\text{FoV}}{\sigma }_{\mathrm{A}\mathrm{L}}^4\\ & +2N_{\text{FoV}}{\sigma }_{\mathrm{A}\mathrm{L}}^2{\sigma }_{\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{b}}^2].\end{array}$$(D.9)

We compared the simple χ² distribution with N − 5 degrees of freedom and this non-trivial normal model to the distribution of 10,000 ${\chi }_{\text{astro}}^2$ values modeled for a single source from the 5p dataset in Fig. D.1. Models are obtained by simulation, as explained in Sect. 4 by generating noisy astrometric measurements of GJ 832 by Gaia, assuming zero orbital motion, and fitting-out a five-parameter model, leading to residuals R_ℓ and a χ² as defined in Eq. D.1. We indeed found that the non-trivial normal distribution is a more accurate model of the true probability law that is followed by ${\chi }_{\text{astro}}^2$, than a naïve χ²-law.

References

---

All Tables

All Figures