© The Authors 2026

1 Introduction

The mass of an astronomical object is one of its most fundamental properties. Masses can be measured in a variety of ways, but these all ultimately rely on a calibration via dynamical measurements. Observations of the mutual orbit of two gravitationally bound bodies provides the means to determine the masses of both components.

The Gaia mission provides the data from which the orbits of millions of systems can be determined (Gaia Collaboration 2016). As part of its third data release (Gaia Collaboration 2023b), the Gaia Data Processing and Analysis Consortium (DPAC) published solutions for about 840 000 spatially unresolved non-single stars (NSSs). Of these, around 443 000 have two-body Keplerian orbit solutions. This includes astrometric, spectroscopic, and eclipsing binaries, as well as combinations of these.

Orbital data alone only allow us to infer the masses of both components in a limited set of cases, namely when we have radial velocities of both components (double-lined spectroscopic binaries, SB2) together with either astrometric or eclipsing orbits (Southworth 2020). This is unfortunate, because SB2 data is only available for systems that are bright enough to have good single-epoch radial-velocity measurements and where the secondary is not too faint (mass ratio not too large). It would be helpful if we could obtain masses just from astrometric orbits, because these are available for many more systems, including fainter sources, sources in crowded fields, and highly inclined orbits. To achieve this, we need additional information on the component masses. This can be obtained from broad-band photometry.

For insight into how photometry helps infer mass, an unresolved binary comprising two identical stars is twice as bright as just one such star, and it will therefore lie 0.75 mag higher in the colour-absolute-magnitude diagram (CAMD). To identify this shift, we need to convert the observed apparent magnitude into the absolute magnitude. We can do this via the parallax, which Gaia determines self-consistently along with the NSS astromet-ric orbit. For an unequal-mass binary, the shift in the CAMD will be different. However, as there is a relationship between mass and absolute magnitude for each component of the binary, the magnitude shift in the CAMD contains useable information on the system's mass ratio. The absolute magnitude does not only depend on the mass, however. By using (PARSEC) stellar models, we can fit a relationship between stellar mass, metallicity, and age on the one hand and absolute magnitude on the other.

In this paper, we used the NSS astrometric orbit together with the parallax and three-band Gaia photometry (G, BP, RP) to estimate the masses of both components of an unresolved binary. For such a binary, Gaia can only observe the motion of the system's photocentre, so the measured semi-major axis (SMA) depends not only on the components' mass ratios, but also on their flux ratios (see Figure 1). We can predict this flux ratio for given stellar parameters from the stellar model, subject to the constraint that their summed flux agrees with the observed flux. Assuming a common metallicity and age for the two components, this gives us four stellar parameters to infer from the three-band photometry (each for the two components together), the parallax, the photometric SMA, and the orbital period.We performed this inference by defining the likelihood of the data given the stellar parameters, adopting a broad prior, and then sampling the resulting posterior probability density function (PDF). The age and metallicity are only weakly determined by the photometry – and not at all by the astrometry – but their presence is necessary to account for the broad scatter in the mass–absolute magnitude relationship.

Although the term NSS refers to its subject as stars, there is nothing in the Gaia measurement technique that precludes the detection of sub-stellar mass companions, namely brown dwarfs or exoplanets. Their lower mass perturbs the primary star less, however, making it harder to detect such companions. When the companion is also very faint, inference is made easier by assuming it contributes no light to the observed photometry. However, we cannot assume a priori that the companion is either low mass or faint. This is relevant because a faint low-mass object around a bright massive object has an identical astrometric orbit as a stellar binary comprising two stars of similar masses and brightnesses. Neglecting the light from the companion would force the former solution. We must instead permit a variable flux ratio in the inference.

Several other studies have estimated masses, or limits thereon, of unresolved binaries using the Gaia astrometry. As part of the validation of Gaia Data Release 3 (Gaia DR3), the Gaia Collaboration (2023a) explored estimating masses for one or both components of the system for several different types of NSS solutions. Similar in spirit to our work, Pérez-Couto et al. (2024) adopted a mass–luminosity relation to help estimate masses. This exploited the binary star results provided by the Multiple Star Classifier (MSC) as part of the Gaia processing (Fouesneau et al. 2023), which estimates stellar parameters of unresolved binaries using their combined BP/RP spectrum. Other studies have combined Gaia data with astrometry from Hipparcos and/or radial velocities from other surveys, often to extend the temporal baseline (e.g. An et al. 2025) The Gaia data have also been used to identify black-hole (El-Badry et al. 2023; Gaia Collaboration 2024) and brown-dwarf (Winterhalder et al. 2024; Wallace et al. 2025) companions to stars.

The paper is laid out as follows. In Section 2, we describe our inference method. We summarise the results in Section 3 and compare these to other estimates in Section 4. We explore the impact of adding infrared photometry and spectroscopic orbital information to the inference in Sections 6 and 7, respectively. We finish up with a discussion of some of the limitations and caveats in Section 8 and conclude in Section 9. Appendix A presents the application of our method to simulated data to give a lower bound on the performance. Appendix B investigates the statistical contamination of stellar companions when looking for exoplanet companions.

Fig. 1 Sketch of two objects (1 and 2) of different masses orbiting their barycentre, with a flux ratio that places the photocentre away from the barycentre. In an unresolved binary, we observe the motion of the photocentre at about the barycentre.

2 Methods

2.1 What astrometry tells us

Kepler’s third law relates the physical SMA a to the period T and the sum of the masses of the two components: $${\left(\frac{a}{\text{au}}\right)}^3={\left(\frac{T}{\text{yr}}\right)}^2\left(\frac{{ℳ}_1+{ℳ}_2}{{ℳ}_{\odot }}\right),$$(1)

where yr is approximately the Earth’s sidereal year. From the Gaia epoch astrometry, one infers not the physical SMA but the SMA of the photocentre, a_p, which is related to a via (van de Kamp 1956) $$a_p=a\left|\frac{1}{1+{ℳ}_1/{ℳ}_2}-\frac{1}{1+f_1/f_2}\right|,$$(2)

where f₁ and f₂ are the fluxes of the two components in the band in which the astrometry is measured (in our case, the Gaia G band). Note that this expression is invariant under a permutation of the subscripts.

The DPAC used the Gaia epoch astrometry to compute both T and a_p. The NSS catalogue actually lists the orbital parameters as Thiele–Innes parameters (and their covariances) instead of the more familiar Campbell parameters, but Halbwachs et al. (2023) showed how to compute a_p and its uncertainty from these (and supplied a Python tool to do so). They provided the Thiele–Innes parameters in angular units (mas), so we divided by the parallax to put this in linear units (au). This assumes the distance equals the inverse parallax, which, although generally not true (Bailer-Jones 2015), is a good approximation for the vast majority of the systems we analysed due to their high-parallax signal-to-noise ratios (S/Ns).

Combining Equations (1) and (2) gives $${\left(\frac{a_p}{\text{au}}\right)}^3{\left(\frac{T}{\text{yr}}\right)}^{-2}=\left(\frac{{ℳ}_1+{ℳ}_2}{{ℳ}_{\odot }}\right){\left|\frac{1}{1+{ℳ}_1/{ℳ}_2}-\frac{1}{1+f_1/f_2}\right|}^3.$$(3)

The Gaia astrometry gives the left side of this equation, which we call the ‘astrometric measurement’ and denote with symbol A. When f₂ ≪ f₁, the right side of this equation is independent of the fluxes, so we obtain an estimate of a function of ℳ₁ and ℳ₂. For this reason, the right side is sometimes called the ‘astrometric mass function’. If we additionally assumed ℳ₂ ≪ ℳ₁, then the equation becomes ${ℳ}_2^3\simeq A{ℳ}_1^2$. To obtain estimates of each mass individually, we need at least one more mass-sensitive measurement, such as absolute magnitude in a passband or the radial-velocity variation.

2.2 Introducing photometry to break the mass degeneracy

The simplification of neglecting the flux of the secondary is often made when looking for faint, low-mass companions to stars, as this leads to a much simpler approach, but this overlooks the degeneracy between mass and flux present in Equation (3). Two stars of near-equal masses and near-equal fluxes can produce an arbitrarily small astrometric measurement. Yet, two objects with the right combination of large mass and flux ratios, such as a star and a planet, can produce the same astrometric measurement. So, unless we can constrain the fluxes, a given astrometric measurement can correspond to any combination of masses.

Our model accommodates the general case where f₂ is not negligible. Clearly it would be useful to know the flux from each component, but because the system is unresolved we can only measure the sum of their fluxes. If absolute magnitude depended only on stellar mass in some known way, then the sum of the fluxes together with the parallax (to convert flux to absolute magnitude), along with the astrometric measurement, would allow us to solve for the two masses (two equations, two unknowns)¹. In reality, however, the absolute magnitude on the main sequence is also a function of age, metallicity, and other parameters, so there is not a one-to-one relationship between mass and absolute magnitude. This means we need the absolute magnitude in more than one band in order to estimate individual masses (and other parameters). To do this, we constructed a forward model of the absolute magnitude in each band as a function of mass, age, and metallicity. These, together with Equation (3), form a set of non-linear equations that can be solved numerically.

Real stars and noisy measurements do not conform precisely to a set of deterministic equations, however, so these will not provide a solution. We therefore adopted the likelihood approach in which we define the probability of the measurements given the masses (and other parameters) and then sought the most likely parameters. Rather than simply maximising the likelihood, we defined and sampled the posterior probability density function (PDF), which is the product of the likelihood with a prior over the parameters. By sampling the posterior, we obtained a point estimate of the parameters and a confidence interval on that estimate.

2.3 Photometric forward models from PARSEC

We used PARSEC 1.2S stellar models (Bressan et al. 2012) to define the relationship between absolute magnitude (M) in each of the three Gaia bands – G, B (=BP), and R (=RP) - and the stellar parameters mass, ℳ; age, τ; and metallicity, [M/H]. PARSEC is a stellar evolution code that solves the equations of stellar structure using information on nuclear-reaction rates, material opacities, and transport processes to predict the variation of stellar properties as a function of time. We used the public interface of PARSEC to generate a synthetic grid of stars over the mass range of 0.09-5.0 ℳ_⊙ with a variable mass step size, a log₁₀(τ/yr) range of 7.6–10.0 in steps of 0.15 dex, and a metallicity range of −2.0 to 0.25 dex in steps of 0.15 dex, only retaining evolutionary stages 1 (main sequence), 2 (sub-giant branch), and 3 (red giant branch). To fit the forward models, we only selected sources within the 2.6 < M_G < 15.4 and 0.0 < M_B – M_R < 4.5 box. This reduced the mass range to 0.09-1.61 ℳ_⊙ (the age and metallicity ranges were unchanged) and left us with a set of 7985 models to fit the forward models. These are shown in Figure 2. The majority (7333) are main-sequence stars.

The forward model for the absolute magnitude was defined for each band independently by the 3D function $$M_{i,X}=y_X\left({ℳ}_i,\tau ,[\text{M}/\text{H}]\right),$$(4)

where X is one of G, B, or R, and i = 1 or 2 to indicate the primary and secondary component, respectively. We experimented with various tree and neural network architectures for fitting these models, concluding on gradient-boosted trees as implemented by GradientBoostingRegressor in the scikit-learn Python package of Pedregosa et al. (2011). We used a squared error-loss function with 200 boosting stages. 85% of the PARSEC models were used for training and the rest for testing. The median absolute deviations of the fits are 0.013 mag, 0.014 mag, and 0.010 mag in G, B, and R, respectively. The bias of each fit is less than 0.001 mag in all cases.

To illustrate the models, we made 4000 stars by drawing log mass, log age, and metallicity from uniform distributions over the range of the fits and simulate photometry for each. The variation of M_G with mass colour-coded by metallicity is shown in Figure 3.

The forward models predict the absolute magnitudes for each component, from which we were able to compute the corresponding flux from $${\widehat{f}}_{i,X}={10}^{-0.4\ast \left(M_{i,X}+z_X-5\left[1+{\log }_{10}\varpi \right]\right)},$$(5)

where ϖ is the system parallax in arcseconds and z_X is the magnitude zeropoint of the Gaia photometric system². This assumes there is no interstellar extinction, which is a reasonable approximation for many stars within a couple of hundred of parsecs, but may be a source of error for more distant stars. We assess this assumption in Section 3.2. Equation (5) also assumes the parallax is the inverse of the distance, which, because the parallax S/N is large for most systems in this study, is not a limiting approximation.

The lower mass limit of the grid is 0.09 ℳ_⊙, but we would like to be able to infer the masses of lower mass objects. As such low-mass objects contribute little flux, it might be sufficient to simply adopt zero-flux contribution when the proposed mass falls below this grid limit during the inference. This is problematic, however, when both components in a system are below this limit (although Gaia sees such systems only when they are very nearby). In lieu of reliable grids for lower mass objects, we simply linearly extrapolated the absolute magnitude in each band as a function of log mass for masses below 0.09 ℳ_⊙, with no dependence on age or metallicity. This is not very accurate, but it is preferable to zero flux. For systems with primary masses above 0.09 ℳ_⊙, the exact extrapolation has little influence on the inferred masses. It turns out that none of the primaries in our Orbital300 sample (defined below) have a median mass below this limit.

Summing the fluxes for both components in a given band gives the flux we would expect to measure: $${\widehat{f}}_X={\widehat{f}}_{1,X}+{\widehat{f}}_{2,X}$$(6)

Adopting a Gaussian noise model for the Gaia flux measurements and assuming the forward model is true, the probability of measuring a flux, f_X, from an unresolved binary system is described by a Gaussian with mean of ${\widehat{f}}_X$ and standard deviation of Σ_fx. The latter is $${\text{Σ}}_{f_X}^2={\sigma }_{f_X}^2+0.848s_X^2\left({\widehat{f}}_{1,X}^2+{\widehat{f}}_{2,X}^2\right)$$(7)

which is the quadrature combination of the uncertainty in the measured flux, σ_fX, and the uncertainties in the forward-model fit for each component, s_X (in magnitudes). The factor (ln(10)/2.5)² = 0.848 comes from the flux-to-magnitude conversion. We set s_X equal to the median absolute deviation of the forward-model fits quoted above.

Fig. 2 CAMD of the PARSEC models used to fit the forward models of Equation (4).

Fig. 3 Variation of MG with mass predicted by the fitted model of Equation (4), colour-coded by median metallicity.

2.4 The posterior probability density function

The posterior PDF is a product of the prior over the parameters, P(θ), and the four likelihood terms: $$P\left(\theta \mid A,\varpi ,\left\{f_X\right\}\right)=P(\theta )P(A\mid \theta ){\displaystyle \prod }_{X}P\left(f_X\mid \theta ,\varpi \right),$$(8)

where $A=a_p^3/T^2$ is the astrometric measurement and θ = (ℳ₁, ℳ₂, τ, [M/H]) are the four sampling parameters. The last term on the right side is the product of the three Gaussian likelihoods in the Gaia fluxes described in Section 2.3 (X = G, B, R). The middle term is the probability of obtaining the astrometric measurement given the parameters. When the model and measurements are perfect, these quantities are linked by Equation (3). However, the right side of Equation (3) is model-predicted (from the parameters θ), whereas the left side is measured. Making the convenient (but not entirely accurate) assumption that deviations between these are Gaussian, the likelihood term P(A | θ) is thus a Gaussian distribution with mean equal to the right side of Equation (3) and standard deviation σ_A, where $${\left(\frac{{\sigma }_A}{A}\right)}^2={\left(\frac{3{\sigma }_{a_p}}{a_p}\right)}^2+{\left(\frac{2{\sigma }_T}{T}\right)}^2$$(9)

is obtained by combining the uncertainties in the SMA and the period in quadrature, which is reasonable for small fractional uncertainties. The posterior is therefore the product of four 1D Gaussian likelihoods and the prior. As mass and age must be positive, the sampling is done over the logarithm of the masses and age (metallicity is already a logarithmic expression of abundance).

We used a separable prior in all four parameters. The mass priors are uniform in log mass between 0.001 ℳ_⊙ and the maximum of the forward model grid (1.61 ℳ_⊙), which is minimally informative (Bailer-Jones 2017; Section 5.3.1), and zero outside. We later imposed absolute-magnitude and colour cuts on the data to exclude stars likely to be more massive than this upper limit. The age prior corresponds to a constant star formation rate from the maximum age in the grid, τ = τ_max, to the present date, τ = 0, so the prior is P(log τ) ∝ τ(τ_max – τ) for ages within the range of the grid and is zero outside it. The prior in metallicity increases linearly from the lowest metallicity in the grid, [M/H]_min, and is therefore P([M/H]) ∝ [M/H] – [M/H]_min for metallicities within the range of the grid and is zero outside it. This linear increase (triangle distribution) could reflect a steady enrichment over time. Neither this metallicity prior nor the age prior is particularly realistic, but they are arguably more realistic than uniform log priors. The Gaia fluxes are only weakly sensitive to age and metallicity, so a more complex and informative prior is not justified. We later demonstrate that adopting the least informative priors possible (uniform log) instead barely changes our results.

2.5 Sampling and summarising the posterior

We sampled the posterior using a Markov chain Monte Carlo (MCMC) method (emcee; Foreman-Mackey et al. 2013) with 32 samplers (chains) for 1000 iterations. The initial samples were drawn uniformly (in the sampling parameters) over the full range of the grid. For each of the 32 initial samples, the higher mass is assigned to the primary and the lower mass to the secondary. This was done just to provide a clear initial distinction between the primary and the secondary. During the sampling, there was no constraint for one mass to be higher than the other, and samples were never swapped. This is essential in order to retain detailed balance in the MCMC method and to allow the mass posteriors to overlap.

The last 500 iterations of the sampling, thinned by a factor of five, represent the 4D posterior. We summarised each marginal posterior using three quantiles: the median, the 15.87th percentile, and the 84.12th percentile. These are named – for the primary mass, for example – ℳ_1,med, ℳ_1,lo, and ℳ_1,up respectively. The median is taken as the point estimate. The difference between the median and the 15.87th percentile is the lower 1σ-like uncertainty ('lower error bar'), and the difference between the 84.13th percentile and the median is the upper 1σ-like uncertainty ('upper error bar'). If the median mass of the secondary was higher than the median mass of the primary, they were relabelled at this point, although this is rarely necessary.

Our Monte Carlo chains are relatively short, although the number of chains (walkers) compensates for this to some degree. As we report below, longer chains do not change the inferred masses significantly.

2.6 Adding spectroscopic (radial-velocity) information

In spectroscopic binaries, we observe variations in the radial velocity of the brighter component due to its orbit around the barycentre (type SB1). These measurements permit the determination of two more Thiele–Innes parameters (published in Gaia DR3; Gosset et al. 2025). Combining these with the astrometric data, we can determine the spectroscopic SMA a_s (Halbwachs et al. 2023, Appendix B). This quantity is related to the physical SMA via Equation (3) again, but this time with f₁/f₂ = ∞ (and so 1/(1 + f₁/f₂) = 0), because in an SB1 we only observe the motion of the primary. We therefore define the 'spectroscopic measurement' as the left side of the following equation: $${\left(\frac{a_s}{\text{au}}\right)}^3{\left(\frac{T}{\text{yr}}\right)}^{-2}=\frac{{ℳ}_2^3}{{\left({ℳ}_1+{ℳ}_2\right)}^2,}$$(10)

with masses in solar units again.

2.7 Comparison statistics

To quantify the overall difference between two mass estimates for a set of stars, we used two metrics. The first is the median absolute residual (MAR): $$\text{MAR}=\text{median}\left(\left|y_i-x_i\right|\right),$$(11)

where y_i and x_i are the two mass estimates for object i. This is a measure of the scatter between the estimates. The bias between the estimates is the median residual: $$\text{bias}=\text{median}\left(y_i-x_i\right).$$(12)

Note that the MAR includes the impact of a bias, in the sense that if x and y differ just by a large constant factor, both the bias and the MAR will be large. When comparing our results to those by other authors, y are our results and x are the comparison results.

3 Results using astrometry and Gaia photometry

We selected all systems with nss_solution_type= ‘Orbital’ from the nss_two_body_orbit table³ in Gaia DR3 that lie within our Gaia photometry limits, have matches in the infrared 2MASS (Skrutskie et al. 2006) and A11WISE (Wright et al. 2010; Mainzer et al. 2011) point source catalogues, and that lie within 300 pc of the Sun. These Gaia photometry limits are the same as those we used to build the forward models as described in Section 2.3, with the absolute G-band magnitude computed neglecting extinction and with distance equal to inverse parallax. We only make use of the infrared photometry later in Section 6. This selection provides 20334 systems, which we refer to as the 'Orbital300' sample. Its CAMD is shown in Figure 4. The sample appears to consist almost entirely of main-sequence stars (plus a few white dwarfs), although secondaries could be other stellar or sub-stellar types.

The distributions of the astrometric data are shown in Figure 5. The 1st and 99th percentiles of the period distribution are at 79 days and 2012 days (5.5 years), respectively. Shorter periods are harder to detect, because they have smaller astrometric signal. Periods much longer than the data span of Gaia DR3, which is 1038 days (2.8 years), are, similarly, harder to detect. The period distribution shows a sharp dip at around 365.25 days, which is the orbital period of Gaia around the Sun. The S/N of the astrometric measurement $a_p^3/T^2$ ranges from 1.8 to 55.8 (1st and 99th percentiles), with a median of just 9.4. The median parallax S/N is 210; 99% of sources have a parallax S/N above 49.

We ran our parameter inference and summarised the posteriors as described in the previous section. The distribution of the median of the inferred masses is shown in Figure 6. The median of the mass distribution is 0.44 ℳ_⊙ for the primaries and 0.07 ℳ_⊙ for the secondaries. There is a weak positive correlation between the primary and secondary masses. The lowest median primary mass is 0.093 ℳ_⊙, and the lowest median secondary mass is 0.0023 ℳ_⊙ (the limit set by the prior is 0.001 ℳ_⊙ in both cases). We discuss the lowest mass secondaries in Section 5.

Figure 7 shows the distribution of the mass ratio for three different primary mass ranges. As we are particularly interested in the lowest mass companions, we used the parameter r = ℳ₁/ℳ₂ for the mass ratio, rather than its inverse. For the log-mass ratio, this simply changes the sign. The peak of log₁₀ r is at 0.2–0.4, corresponding to mass ratios of 1.6–2.5 (this is not the maximum of r, because a maximum is not invariant under a log transformation). Thus, for lower mass primaries we do not see as many large mass ratios. This could be in part a consequence of the lower sensitivity to lower mass secondaries.

We refrained from fitting models to the mass ratio distribution and from inferring population statistics, because the NSS sample suffers from various selection effects, as emphasised by Halbwachs et al. (2023). It is notable, however, that for mass ratios below about 1.6, there is a steep drop in the number of binaries, meaning we find few near-equal-mass binaries. Gaia's high-precision astrometry makes it suitable for identifying high-mass-ratio systems. Nonetheless, the larger the mass ratio, the smaller the astrometric measurement (Equation (3)) in general, and so the more difficult it is to detect high-mass-ratio systems. This observational selection effect could explain some of the slope in Figure 7.

The variation of the mass ratio with the orbital period is shown in Figure 8. We see a slight increase in the median mass ratio at intermediate periods.

The CAMD of the sample colour-coded by the median component masses is shown in Figure 9. The median over a set of stars is shown at each point. In both panels, the magnitude and colour are those observed for the system. In the upper panel, the primary mass varies fairly smoothly and continuously, with the mass decreasing for fainter and redder systems. This is what we expect when the primary dominates the system photometry. The colour-coding by the secondary mass in the lower panel shows a roughly similar trend. However, for systems brighter than M_G ≃ 8.0 mag we see a locus of systems in the central part of the main sequence with higher average secondary masses, and thus lower mass ratios.

So far, we have only considered the point mass estimates, which are medians of the posterior distributions. However, the posteriors can have significant widths, something we characterise with the lower and upper 1σ-like uncertainties (see Section 2.5). Dividing each of these two widths by the median gives a lower and an upper fractional uncertainty, respectively. These are plotted for the primary and secondary mass estimates in Figure 10. Note that the lower fractional uncertainty cannot be larger than 1.0, whereas the upper one is unbounded. The uncertainties for the primaries are relatively small, with median lower and upper fractional uncertainties of 0.07 and 0.05, respectively. 90% have lower and upper fractional uncertainties less than 0.23 and 0.12, respectively. The fractional uncertainties for the secondaries, in contrast, are often much larger; only 50% have lower and upper fractional uncertainties less than 0.25. A sizeable fraction have lower fractional uncertainties that are nearly as large as 1.0. This means that the posterior permits an arbitrarily low mass (down to 0.001 ℳ_⊙, the lower limit of the inference). There is no obvious trend of the size of this fractional uncertainty with the size of the astrometric measurement.

Some example posteriors for the primaries and the secondaries are plotted as histograms in Figure 11. The first panel shows a typical case. It is quite common for a small fraction of the MCMC samples to be 'disconnected' from the main set of samples when plotted as a histogram in this way, but the lower 'peaks' are often either near to the main one or contain few samples. Significant bimodality is very rare in the full sample. We do, however, see a few cases where there is a long, low tail to low masses for the secondary, such as in panels 2 and 3 of Figure 11. Panel 4 shows an example where some MCMC chains moved between the primary and secondary during the sampling. However, the number of samples attributed to the other component is small, so the median is unaffected.

To test the impact of the prior on our results (Section 2.4), we reran the inference using uniform log priors on age and metallicity, which are minimally informative but also less physical. Averaged over the sample, this increased the primary masses by 0.007% and the secondary masses by 0.1% (Equation (12)), which is negligible. There is a slightly larger scatter (Equation (11)) of 1% for primaries and 3% for secondaries – meaning a single source selected at random would differ by these amounts, on average – but these differences are much less than the widths of the posteriors, and so they are accounted for by the error bars.

As MCMC sampling is a stochastic process, the posteriors and their summaries vary if the MCMC sample is extended for more iterations or rerun with a different intialisation. When we extended the MCMC sampling to five times as many iterations (also for the burn-in), we found that the mass estimates changed by only 1.5% for the primaries and 3.9% for the secondaries (scatter; Equation (11)). Rerunning with a different initialisation caused less scatter. Convergence noise therefore has only a small impact on our results.

We provide our mass estimates as an electronic table with the entries shown in Table 1. The log-flux ratio is the ratio computed from the forward-model-predicted G-band fluxes (Equation (4)) at the median inferred values of the masses, age, and metallicity for each system. This flux ratio is not equal to the value computed from Equation (3) using the median masses and measured a_p and T, because the median stellar parameters are found by posterior sampling, and not by solving that equation.

There is an obvious set of 14 white dwarfs (primaries) in the bottom left corner of the CAMD of Figure 4. Our photometric forward models are not appropriate for these stars, so we do not expect to infer appropriate masses from them. Indeed, we obtained median primary and secondary masses of around 0.1 ℳ_⊙ and 0.008 ℳ_⊙, respectively, for these systems. We left them in the published catalogue because they were not excluded from our analysis (they have negligible impact on the later comparison statistics).

Fig. 4 CAMD for the Orbital300 sample assuming zero extinction.

Fig. 5 Distribution of the data for the Orbital300 sample. Top row: the photometric SMA ap and period T. Middle row: the astrometric measurement $a_p^3/T^2$(with ap in au and T in yr, from the left side of Equation (3)) (left) and its S/N (right). Bottom row: parallax S/N and distance (limited by construction to less than 300 pc).

Fig. 6 Inferred masses (median of the posteriors) for the Orbital300 sample for the primary (ℳ1, top) and secondary (ℳ2, bottom).

Fig. 7 Distribution of mass ratio r = ℳ1/ℳ2 for the Orbital300 sample for three different primary mass ranges.

Fig. 8 Variation of the mass ratio ℳ1/ℳ2 with period for the Orbital300 sample. The red points are the median values of the mass ratio in equally spaced period bins.

Fig. 9 System CAMD for the Orbital300 sample colour-coded by the primary mass (top) and secondary mass (bottom).

Fig. 10 Distribution of the inferred fractional 1σ-like uncertainties in the masses for the Orbital300 sample for the primaries (left panel) and secondaries (right panel). The lower fractional uncertainty cannot be larger than 1.0, whereas the upper one is unbounded. There are a small number of secondaries with upper fractional uncertainties beyond the limit of the plot.

Fig. 11 Examples of marginalised mass posteriors for the primary (blue) and secondary (orange) for four sources (one per panel). Panel 1 is typical, panels 2 and 3 show examples of long tails to low masses for the secondaries, and panel 4 shows where some of the MCMC chains moved between the primary and secondary during the sampling. The vertical dotted lines show the median masses.

3.1 Relevance of flux ratio on the mass estimates

Our forward model (Equation (4)) predicts the G-band flux for each component of the unresolved binary from the mass, age, and metallicity. From this we computed the flux ratio at the median inferred values of the parameters. Figure 12 shows how this flux ratio varies with the inferred mass ratio for the Orbital300 sample. This is broadly what we would expect for stellar binaries – the flux ratio increasing with mass ratio, plus additional variance due to metallicity and age variations – and is a consequence of our forward model.

The impact of the flux from the secondary on the determination of the masses via the astrometry is determined by the size of the cubed quantity on the right side of Equation (3), or more conveniently by $$R=\left|1-\frac{1+{ℳ}_1/{ℳ}_2}{1+{\widehat{f}}_1/{\widehat{f}}_2}\right|,$$(13)

which we call the 'flux relevance ratio'. When the secondary contributes negligible flux, such that the flux ratio ${\widehat{f}}_1/{\widehat{f}}_2$ is very large, R is close to one. When, instead, this flux ratio is closer to unity for a given mass ratio, R is less than one, and the smaller R is, the more relevant the flux from the secondary is in the interpretation of the astrometric measurement. The distribution of R is shown in Figure 13. The majority of systems have R near to unity, but there is tail of systems with R as low as 0.4, and a handful are at even lower values. The figure also shows how R varies with the primary mass. Unsurprisingly, it is larger for lower primary masses, because this is where the effect of the secondary on both the astrometry and photometry is more noticeable.

We tested the importance of accounting for the flux of the secondary more directly by repeating the inference with this flux fixed to zero. In this case, f₂ = 0 both in Equation (3) and when computing the photometric likelihoods. We find that the estimated masses of the secondaries are then 9% lower on average (median), and 10% of the sources have masses that are 40% or more lower. This confirms that we need to accommodate for the unknown flux contribution of the secondary when trying to estimate its mass.

Fig. 12 Predicted G-band flux ratio versus inferred mass ratio for the Orbital300 sample.

Fig. 13 Flux relevance ratio (Equation (3)) plotted against the inferred mass ratio (left) for the OrbitaBOO sample shown on a logarithmic density scale. The right panel shows the distribution of the flux relevance ratio.

3.2 Correction for interstellar extinction

So far, we have ignored interstellar extinction. Such extinction would dim and redden the observed photometry of some stars to a certain degree. These stars would, in reality, be brighter and bluer than measured, and thus more massive than inferred. Individual extinction estimates for some of our target stars are available in Gaia DR3 from the GSP-Phot algorithm (Andrae et al. 2023; Fouesneau et al. 2023). This algorithm self-consistently estimated the extinction parameter A₀ together with the intrinsic stellar parameters via stellar models. From these and the photometric passbands, the algorithm computed the extinction in each of the G, BP, and RP bands. Two-thirds of our Orbital300 sample – 14 173 sources – have extinction estimates, 13 710 of which have extinction-corrected magnitudes that lie within the bounds of our forward model. Their median G-band extinction A_G is 0.21 mag (10th and 90th percentiles are 0.004 and 0.56 mag). The GSP-Phot authors noted that their extinction estimates may not be very accurate; in particular, some may be overestimated. However, these extinction values seem plausible for sources over the whole sky within 300 pc, and so we can use them at least to assess the impact of neglected extinction.

We used the GSP-Phot extinctions to adjust the forward-model predicted fluxes in Equation (5) for all three bands during the inference procedure. We then compared the inferred masses (posterior median) source-by-source with the extinction-free masses inferred earlier. When including the extinction, the primary masses are estimated to be 2.2% higher on average (median), and the secondary masses 0.4% higher. The ℳ₁/ℳ₂ ratio is 1.6% higher. Hence, there is very little systematic increase in the mass estimates when including plausible extinctions. The scatter between the estimates is somewhat larger than these biases, indicating that some individual stars' mass estimates vary by more, although some of this is also MCMC convergence noise. Overall, we can say that for stars within 300 pc, neglecting extinction has very little impact on the mass estimates with our method.

4 Comparison to other estimates

As part of the Gaia DR3 validation, Gaia Collaboration (2023a) examined various combinations of Gaia NSS data to estimate the mass of the primary and/or secondary. Which data were used to estimate which masses is indicated by the flag combination_method in their table⁴. In several cases, they had to adopt a value for the mass of the primary in order to estimate the mass of the secondary. Often, this primary mass was estimated using Gaia photometry and isochrones. This is similar to what we did in the present study, except that we estimated the masses of both components simultaneously using the astrometry and photometry. Here, we refer to the mass estimates from Gaia Collaboration (2023a) as the 'comparison set'. Some of their ℳ₁ masses are lower than their ℳ₂ masses, which we then swapped to be consistent with our results (our primary always has the higher median mass). For the Orbital300 sample analysed in the previous section, the comparison set provides estimates of the primary mass in 19 603 systems, but of the secondary mass in only 1574 systems, the latter being a subset of the former. They have further cases where the secondary has a lower and/or upper mass limit, but we did not use those here.

Figure 14 compares the two mass estimates for the primaries and secondaries. The agreement for the primaries is reasonable. In contrast, we see some significant differences for many of the secondaries, with our estimates extending to lower masses than the comparison ones. We remind the reader that our estimates are the medians of a posterior which for the secondaries can have a long tail to low masses, as shown earlier in Figure 10.

Figure 10 compares the median mass estimates, yet the uncertainties (posterior widths) in both estimates are non-negligible. To accommodate the uncertainties in a comparison, we define the 'deviation' as the difference between the mass estimates divided by their combined error bar, using the lower error bar (subscript 'lo') of one estimate and the upper error bar (subscript 'up') of the other, according to which estimate is higher (to obtain the 'neighbouring' error bar): $$\text{deviation\hspace{0.17em}=}\{\begin{array}{cc}\frac{y-x}{\sqrt{{\sigma }_{y,lo}^2+{\sigma }_{x,up}^2}}& \text{if}y-x\ge 0\\ \frac{y-x}{\sqrt{{\sigma }_{y,up}^2+{\sigma }_{x,lo}^2}}& \text{otherwise}.\end{array}$$(14)

These deviations – with y as our estimates and x as the comparisons – are plotted against our estimates in Figure 15. While the deviations for the primary masses are mostly just 'one or two sigma', in some cases the deviations for the secondary masses are considerably larger, showing that some differences are statistically significant.

We quantified the difference in the mass estimates using the statistics in Section 2.7. For the primaries, the bias and scatter (MAR) in log₁₀ ℳ₁ are −0.024 dex and 0.028 dex, respectively. Differences in log base 10 quantities correspond to ratios 2.3 times larger, so we can also interpret these values as a bias of −5.5% and a MAR of 6.4% in ℳ₁. In absolute terms, the bias is −0.037 ℳ_⊙ and the MAR is 0.044 ℳ_⊙. The comparison set does include a few very high mass estimates (up to 123 ℳ_⊙), but as our statistics are based on the median, these have little impact.

For the secondaries, the bias and MAR in log₁₀ℳ₂ are −0.16 dex (−38%) and 0.17 dex (39%), respectively. In this case, the MAR is dominated by the bias: the MAR around the bias is 0.11 dex. This systematic difference is unlikely to be caused by us neglecting extinction, because that should equally effect the primary mass estimates. In absolute terms, the bias and MAR are −0.18 ℳ_⊙ and 0.18 ℳ_⊙, respectively, although these values make little sense given that ℳ₂ spans two orders of magnitude (see Figure 6).

Among the primary masses compared, 96% have the flag combination_method in the comparison table equal to Orbital+M1. This means that ℳ₁ was estimated from Gaia photometry. Among the secondary masses, 53% have combination_method set to Orbital+M1, which means that ℳ₂ was estimated from the photometrically inferred ℳ₁ together with the astrometric measurement (the astrometric orbit). The bias and MAR of our log₁₀ ℳ₂ inferences compared to this subset are higher than for the complete set, at −0.22 dex and 0.22 dex, respectively. The complementary subset – the other 47% – mostly contains solutions with combination_method equal to Orbital+SBl+M1. The bias and MAR of our log₁₀ ℳ₂ inferences for these are lower, at −0.10 dex and 0.13 dex, respectively. The comparison mass estimates in this complementary subset have used the spectroscopic orbit in addition to the astro-metric orbit (and photometry), so they should be more accurate than the other comparison mass estimates. Our inferences, in contrast, just used the photometry and astrometry in both subsets.

To investigate this further, we repeated the comparison on a different sample, namely those that have combination_method set to AstroSpectroSBl+M1 and that lie within 300 pc, a set of 9674 systems we call the AstroSpectro300 sample (defined fully in Section 7). These are all single-lined spectroscopic binaries, and so the Gaia Collaboration (2023a) used the combined spectroscopic and astrometric orbit along with an adopted value of ℳ₁ to estimate ℳ₂. We inferred the masses with our method using just the astrometry and photometry (but not the spectroscopy) and compared results. The bias and MAR for log₁₀ ℳ₁ are −0.001 dex and 0.019, respectively, and for log₁₀ ℳ₁₂ they are +0.007 dex and 0.062 dex, respectively. This is a much better agreement than we had above, as can be see in Figure 16, which compares the masses, and Figure 17, which shows the deviations. Part of the reason for the better agreement may be that the AstroSpectro300 sample is brighter – compare Figures 4 and 19 – implying better quality photometry and astrometry than the Orbital300 sample. Indeed, with the Orbital300 sample, there is a trend towards increasing discrepancy in the secondary masses towards fainter magnitudes. There is no such trend with the S/N of the astrometric measurement, however. In any case, the improved agreement with the AstroSpectro300 sample suggests that our method can achieve reasonable mass results for both components even without using the spectroscopic data.

Pérez-Couto et al. (2024) also estimated the masses of astro-metric binaries from Gaia DR3 data, taking advantage of the results from the MSC algorithm in Gaia DR3. MSC attempted to estimate the stellar parameters of both components of an unresolved binary assuming the target was a binary with a flux ratio less than five, using the composite BP/RP spectrum (Fouesneau et al. 2023). The validation undertaken by the DPAC showed that these stellar parameters are not very reliable, however, and should be used with caution (Fouesneau et al. 2023; Babusiaux et al. 2023). Adopting a one-to-one mass-luminosity relation, Pérez-Couto et al. (2024) derived the component masses from the astrometric measurement and the system's G-band photometry – corrected for extinction using the MSC estimate – by adjusting the G-band magnitude difference between the components in discrete steps and selecting the one in best agreement with the data.

A total of 10 321 of their sources overlap with our Orbital300 sample, although 21% of their primary masses lie in the 10–90 ℳ_⊙ range, which seems an implausibly large fraction. Excluding these, we nonetheless find a good agreement between our primary mass estimates, with a bias of +0.006 dex (1.4%) and a MAR of 0.018 dex (4.2%) in log₁₀ℳ₁. The secondary masses show very poor agreement, however. Specifically, our secondary masses are invariably lower. Most of their systems have mass ratios (ℳ₁/ℳ₂) close to 1.0; 90% have mass ratios below 1.5. Their method is simpler than ours, but it also uses fewer photometric data and appears not to allow for the broad intrinsic scatter in the stellar mass-luminosity relation (e.g. Figure 3; see the discussion in Section 2.2).

Fig. 14 Comparison of our mass estimates on the Orbital300 sample using astrometry and photometry with those from Gaia Collaboration (2023a), which adopted a mass for the primary and then inferred the mass of the secondary from a variety of data (in many cases there is no secondary mass estimate). The two panels show the primaries (left, 19 603 stars) and secondaries (right, 1574 stars). A few sources lie outside the range of the plots.

Fig. 15 Deviations, as defined by Equation (14), for the results shown in Figure 14.

Fig. 16 Comparison of our mass estimates on the AstroSpectro300 sample when not using the spectroscopic measurement with mass estimates from Gaia Collaboration (2023a), which adopted a primary mass and then used the astrometry and the spectroscopy to estimate the secondary mass. The two panels show the primaries (left) and secondaries (right), both with 9674 stars. The horizontal and vertical axes do not span a common range in each panel.

Fig. 17 Deviations, as defined in Equation (14), for the results shown in Figure

5 Planetary-mass companions

The Orbital300 sample contains 155 systems in which the secondary has a median mass below 0.012 ℳ_⊙ (about 13 ℳ_J). Many of these have posteriors extending to much higher masses, however, so they should not be taken as good planetary-mass candidates. We instead identified the 22 systems where ℳ_2,up < 0.012 ℳ_⊙, i.e. the 84th percentile of their mass posterior ('median plus one sigma'), is below 0.012 ℳ_⊙. Five of these systems lie on the white-dwarf sequence, so their primary and secondary masses are unreliable. This leaves 17 objects that we refer to as exoplanet candidates. These are listed in Table 2.

The Gaia Consortium published a list of 72 astrometric objects of interest (ASOIs) that may be binaries with an exo-planet secondary⁴. Of these, 52 are in the Orbital300 sample. Of those, ten are in our exoplanet candidate list.

Stefánsson et al. (2025) followed up on 28 of the ASOIs using radial velocities, five of which are in our exoplanet candidate list. Of these, four are designated as SB2 by Stefánsson et al. (2025). We would not expect to see the lines of an planetary-mass secondary in the spectrum, so this suggests these four objects are not planets (false positives in our list). A fifth object is designated 'single', meaning only one set of lines (the primary) was seen. The remaining five ASOIs are designated 'unknown', indicating there has been no known follow-up.

In total, Stefánsson et al. (2025) found 23 SBs. Of these, 22 are in our Orbital300 sample. As just noted, four are in our exo-planet candidate list. It is relevant, however, that the other 18 are not, because it means we did not erroneously identify them as exoplanets. Thus, while our method may not be very accurate in estimating the masses of exoplanet companions - and it is not designed to do so – it shows some promise in identifying stellar or brown-dwarf companions. It can therefore help to clean exoplanet target lists for follow-up.

The above caveats notwithstanding, our mass estimates for one new Gaia exoplanet and one new Gaia brown dwarf identified by Stefánsson et al. (2025) agree well with their estimates, which used high-precision (10–15 m s⁻¹) multi-epoch radial velocities and multi-band photometry in addition to the Gaia astrometry. For Gaia-4b (GDR3 1457486023639239296), they reported a mass of 11.8 ± 0.7 ℳ_J, which compares to our estimate of ${11.8}_{-8.2}^{+2.8}{ℳ}_{\text{J}}$. For Gaia-5b (GDR3 2074815898041643520), they reported a mass of 20.9 ± 0.7 ℳ_J, which compares to our estimate of ${21.4}_{-1.9}^{+1.3}{ℳ}_{\text{J}}$. Neither of these are in our exoplanet candidate list because ℳ_2,up > 0 12 ℳ_⊙ in both cases. The other three Gaia objects with exo-planet designations, Gaia-1b, −2b, and −3b (Panahi et al. 2022; Sozzetti et al. 2023), are not in our Orbital300 sample.

6 Inclusion of infrared photometry

As we used photometry to help estimate the masses of the two components, it is reasonable to expect that more photometric bands could improve the mass estimates, either in accuracy or in precision. In this section, we add the J and K_s bands from 2MASS and the W1 and W2 bands from AllWISE and include these as four extra terms in the photometric likelihood (Equation (8)). The astrometric likelihood is unchanged as we still use the Gaia G-band flux to compute the flux ratio in Equation (3). We applied this model to the Orbital300 set used in Section 3 and compare to the infrared-free results obtained there.

The overall results are very similar, as we see in Figure 18. There is, however, a slight tendency for the inclusion of the infrared photometry to increase the estimated masses of both the primaries and the secondaries, by 2.2% and 3.9%, respectively (median of the mass differences, i.e. the bias in Equation (12)). This could be a consequence of neglecting interstellar extinction. If a non-zero extinction is ignored, the system has a brighter intrinsic magnitude than we measured, which generally corresponds to a higher mass. So, by ignoring extinction we expect to underestimate the mass. As infrared photometry is less affected by extinction than the visual Gaia bands, it is plausible that its inclusion delivers higher mass estimates. It is partly for this reason that we limited our analysis to systems within 300 pc. We do expect some extinction in some directions within this horizon, and this interpretation is perhaps confirmed by a slight increase in the difference in the mass estimates with distance, albeit with a weak correlation (Pearson correlation coefficients of 0.01 and 0.10 for ℳ₁ and ℳ₂ differences, respectively). Yet, the above mass increases are small, and the analysis in Section 3.2 also indicates that neglected extinction has only a small effect on the mass estimates.

To address whether adding infrared photometry improves the precision of the mass estimates, we compared the widths of the posteriors using the lower and upper fractional uncertainties. We recall that Figure 10 shows their distributions without infrared photometry. Including infrared photometry reduces the fractional uncertainties on ℳ₁ by 7% and 34% for the lower and upper error bars, respectively (median reduction). For ℳ₂, they are smaller by −1 % for the lower error bars (i.e. larger, so slightly less precise) and by 14% for the upper error bar. In summary, we find that including infrared photometry gives more precise mass estimates, but it does not change the estimates by more than a few percent.

Fig. 18 Comparison of our mass estimates using infrared photometry (vertical axis) with that not using infrared photometry (horizontal axis) on the OrbitaBOO subset for the primary (left) and secondary (right).

7 Inclusion of spectroscopic orbits

Some of the stars in the Gaia DR3 NSS table have spectroscopic orbital solutions in addition to astrometric orbital solutions. To analyse the impact on the mass inference when also using the spectroscopic measurement (see Section 2.6), we created a new sample of single-lined spectroscopic binaries by selecting all those sources lying within 300pc that have nss_solution_type set to AstroSpectroSB1 and combination_method flag (Section 4) set to AstroSpectroSB1+M1. Retaining only sources within the photometric limits of our forward models (Section 2.3), this gives a set of 9674 systems, which we call the Astro-Spectro300 sample (used in Section 4). The CAMD of this sample is shown in Figure 19. Because these sources must have (good) Gaia radial-velocity spectra, this sample is brighter than the Orbital300 sample: the 95th percentile in the apparent G-magnitude distribution is 12.7 mag, compared to 15.8 mag for the Orbital300 sample.

We next added an extra term for the spectroscopic measurement to the likelihood (Equation (10)). Everything else remained unchanged. Infrared photometry was not used. The left panel of Figure 20 compares the sizes of the spectroscopic and astrometric measurements. The former are larger by 11% on average (median). This is what we expected, because the non-zero flux contribution from the secondary reduces the astro-metric measurement relative to the spectroscopic one (compare Equations (3) and 10),

Figure 21 compares the results of the mass estimates when including spectroscopy with those not including it. Overall, we see very little difference. The primary masses are very similar, differing by less than 0.1% (bias; Equation (12)). The scatter (MAR; Equation (11)) is ten times larger, but it is still small. The secondary masses are slightly smaller on average (bias) when using spectroscopy, but only by 0.7%, and the scatter is only 2.7%. The fractional uncertainties of these mass estimates are smaller when including spectroscopy, by up 40% on average for the upper fractional uncertainty for the secondaries, but only 2–14% for the other three fractional uncertainties.

Overall, we find that including the Gaia DR3 spectroscopic information for this sample does not significantly change the mass estimates for most systems, and it only modestly improves the precision of the estimates. This is an important result, because it implies that Gaia-precision radial-velocity follow-up of stars that lack spectroscopic orbits in Gaia will not help identify the lowest mass companions. However, this conclusion only applies to this sample, which is quite bright and nearby. Here, the S/N of the spectroscopic measurement is on average four times lower than that of the astrometric measurement, as we can see in Figure 20, with median values of 16.4 versus 4.1. More distant stars will generally have poorer astrometric orbits (as the astrometric measurement decreases in size with increasing distance). The Gaia epoch radial velocities are not very precise -which leads to the low S/N in the spectroscopic measurement - so higher precision radial velocities would help in follow-up studies.

Fig. 19 CAMD for the AstroSpectro300 sample assuming zero extinction.

Fig. 20 Histograms of the distributions of the astrometric and spectroscopic measurements (left) and their S/Ns (right) for the AstroSpectro300 sample.

Fig. 21 Comparison of our mass estimates using spectroscopy (vertical axis) with those not using spectroscopy (horizontal axis) on the Astro-Spectro300 sample for the primary (left) and secondary (right), shown on a logarithmic density scale.

8 Limitations of the method

To determine the intrinsic uncertainties in our inference method, we used it to estimate the masses of simulated binary systems. This is described in Appendix A. The systems and the S/N of the data are broadly similar to those of the Orbital300 sample. We found that the primary masses can be estimated with an accuracy (MAR; Equation (11)) of 3% and a bias of less than −0.2% (i.e. a near-negligible underestimation). For the secondary mass the MAR is 5% and the bias is −3%. This is much better than what we achieved on the Orbital300 sample.

There are a number of reasons why the performance on real data is probably worse than that yielded by these simulations. These are as follows:

1. A mismatch between the PARSEC models and real stars.

2. An assumption of zero interstellar extinction.

3. Non-Gaussian likelihoods and/or incorrectly estimated uncertainties.

The first of these issues is a known problem and could well be the dominant issue. In particular, our extrapolation of the fit below the lower mass limit of the PARSEC models should ideally be replaced by a brown-dwarf and exoplanet model. The second issue is a problem in principle, but our analyses using the GSP-Phot extinction values (masses systematically overestimated by less than 2%; see Section 3.2) and the results from including infrared photometry (see Section 6) suggest neglected extinction has little impact on our nearby sample (within 300 pc). We chose not to adjust the photometry using one of the various published 3D dust maps, because these tend not to fit very well in the solar neighbourhood, and they are not designed to give accurate results for individual stars. Concerning the third issue, we found that inflating the Gaia photometric uncertainties by a factor of five changes the inferred masses by much less than 1% on average and changes the S/N in these masses by less than 1% on average. On the other hand, non-Gaussianity in the sometimes noisy astrometric measurements may be an issue (Figure 5), not least because it is computed from the ratio of just two quantities (Equation (9)), so central limit theorem convergence to a Gaussian may not apply. This uncertainty estimate has also neglected the contribution from the scatter in the forward models, which could be a limiting factor when the astrometric measurements are more precise.

Some studies, such as Winterhalder et al. (2024) and Wallace et al. (2025), adopted a simplified version of the method in this paper. By assuming f₂ ≪ f₁, the observed photometry allows the mass of the primary to be determined independently of the secondary, and the simplified version of Equation (3) then yields the mass of the secondary. However, this assumes we already know that the secondary is faint. If the single-star main sequence were an arbitrarily thin locus, the absolute photometry would betray the presence of a bright secondary. However, due to age and metallicity variations the main sequence is not arbitrarily thin, and if the photometry is not dominated by the primary then the single star inference from the photometry will give an incorrect primary mass. Our method avoids this assumption to produce a self-consistent solution for the masses and fluxes of both components.

9 Summary and conclusions

We developed a method to infer the masses of both components of an unresolved binary using the photocentre semi-major axis, orbital period, parallax, and photometric flux measurements. This involved using a forward model to map the relationship between the fluxes and the two masses and two latent variables (system age and metallicity). We fitted this model using PAR-SEC stellar models. We sampled the resulting four-dimensional posterior using a Monte Carlo method.

We tested our method on Gaia astrometry and photometry (three bands) for binaries listed in the NSS table in Gaia DR3. Our main sample (OrbitaBOO) is a set of 20 334 systems within 300 pc that have a system photometry consistent with that of the stellar main sequence. We can estimate primary masses with median precisions (lcr-like uncertainties) of 6% and secondary masses with median precisions of 25%. For a significant fraction of the systems, however, the posterior for the secondaries extends to zero with non-negligible probability density, indicating a large lower uncertainty on the secondary mass.

The numbers above do not refer to deviations from the true masses. Dynamical-only 'true' mass estimates are rarely available, because they require additional information only accessible for a small number of objects (in particular, double-lined spectroscopic binaries). We could, however, compare our results to mass estimates obtained by Gaia Collaboration (2023a). For the OrbitaBOO sample, the primary masses are quite similar: the scatter (median absolute residual in the log masses) is 6.4%, which includes a bias of −5.5% (our estimates are lower). The secondary masses show larger deviations: a scatter of 39%, which includes a bias of −38%. However, on a sample of brighter stars (AstroSpectro300), the agreement is much better. For the secondaries, the scatter is 14%, which includes a bias of just +1.6%. For this, Gaia Collaboration (2023a) also used the spectroscopic orbits, whereas we did not.

Our method accommodates the initially unknown flux contribution of the secondary that affects the size of the observed semi-major axis. We found that neglecting this underestimates the mass of the secondary by 9% on average (median), but by 40% or more for 10% of the sources. Neglecting this could therefore be a source of false-positive low-mass companions.

We did not account for interstellar extinction, and this could in principle lead to mass underestimates. However, we found that correcting for the extinction using values from Gaia DR3 had very little impact on our mass estimates. This is supported by the fact that including infrared photometry (J, K_s, W1, W2) in the inference does not change the mass estimates by more than a few percent. This additional photometry does, however, produce more precise mass estimates (narrower posteriors).

Including Gaia's spectroscopic orbit in the inference hardly changes the mass estimates. This suggests that a Gaia-like spectroscopic follow-up of these systems would not improve the mass estimates. More precise radial velocities at the appropriate epochs would help, however.

Using our method on the larger and more complete set of astrometry that will be published in Gaia DR4, we plan to estimate the masses of both components of many more binaries using only the astrometric orbit and photometric data. This should allow us to identify sub-stellar companions without too much contamination from stellar companions that have the same astrometric signal.

Data availability

Full Table 1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/708/A249

Acknowledgements

This work uses 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 publication also uses data products from the Two Micron All Sky Survey (2MASS), a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation, and from the Wide-field Infrared Survey Explorer (WISE), a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. We also made use of CDS's Simbad (Wenger et al. 2000), NASA's ADS, and matplotlib (Hunter 2007) in this work.

References

Appendix A Performance on simulated astrometric binaries

Fig. A.1 Performance of mass estimation on the simulated data set using astrometry and Gaia photometry. Note that in each panel, the horizontal and vertical axis do not span a common range.

To determine the limiting performance of our inference method, we use it here to infer the masses of simulated astrometric binaries. To generate our sample, we first draw parameters of individual stars at random from uniform distributions in log₁₀ ℳ, log₁₀(τ/yr), and [M/H]. The limits of these distributions are the range of the forward model grids (defined at the beginning of section 2.3), except for the lower mass limit, which is set to a lower mass of 0.01 ℳ_⊙. We then use the fitted forward models described in section 2.3 to predict the photometry for each star, including the extrapolation for the lower mass stars.

We generate binaries by randomly pairing stars. This is done by drawing a period at random for each system from a distribution uniform in log period between 100 and 1500 days, then computing the SMA from equations 1 and 2 and the true (noise-free) astrometric measurement from equation 3. We assign a parallax to each system at random by drawing a distance from a uniform space density distribution (P(r) ∝ r²) up to 1 kpc with a lower limit set for each system to ensure that the photocentre SMA is below 100 mas and so would be unresolved by Gaia.

To generate noisy data for the inference we add Gaussian noise to all inputs. For simplicity we assume fixed S/Ns for each input with the following values: 9.9 for the astrometric measurement (this is dominated by the uncertainty in the SMA); 200 for the G band flux; 100 for each of the BP and RP band fluxes; 150 for the parallax. This S/N in the astrometric measurement is close to the median of the Orbital300 sample of 9.4. The high, constant S/Ns for the Gaia fluxes essentially assume that the uncertainties are dominated by systematic errors rather than photon noise, which is consistent with the NSS sources being bright. The uncertainties corresponding to these S/Ns are also used as the standard deviations σ_A (equation 9) and σ_fx (equation 7) in the likelihoods.

We compute and summarize the posteriors as done for the real data (section 2.5). The median of the posteriors for the primary and secondary masses are compared to the true masses in figure A.1. For the primaries, the bias and MAR log₁₀ ℳ₁ are −0.0008 dex and 0.013 dex respectively, that is, negligible bias and a median scatter of 3%. For the secondaries the bias and MAR are −0.0136 dex and 0.0236 dex respectively. The former means the inference underestimates the masses by 2% on average, it being generally lower for low mass systems and higher for higher mass systems.

If we include the four infrared photometric bands as additional measurements, the results are very similar. The biases for the primaries are secondaries are quite similar, as is the MAR on the secondaries. But the MAR on the primaries is three times smaller when using infrared photometry. The S/Ns on the mass estimates of the primaries are also about twice as large, for both the lower and upper uncertainties, when using the infrared photometry.

This performance on simulated data is considerably better than the consistency we see in section 4 between our real results and other literature estimates, in particular for the secondaries. Possible reasons for this are discused in section 8. The simulations obviously present an ideal situation, in which the measurements - one astrometric measurement and the three fluxes - differ from the true values by only random Gaussian noise. This nonetheless indicates that in principle our method of determining the four system parameters – two masses, age, and metallicity – from the four inputs can infer accurate masses, even for the secondary. These results on simulated data represent an approximate limit on the mass estimation accuracy that we can achieve on real data.

Appendix B Contamination of stellar companions in detecting exoplanet companions

A follow-up goal of this work when Gaia DR4 data become available is to estimate the masses of exoplanet companions in unresolved astrometric star-planet binaries. As the planet contributes negligible light to the integrated photometry, the estimation of the mass of both the stellar and planet mass from the astrometric measurement and the flux is easier than considered in this paper, because we can then assume f₁/f₂ = ∞ in equation 3 and need the forward model to map mass to luminosity for the stellar primary only.

However, we do not know a priori that the secondary is a planet. As discussed in section 2.2, a near equal-mass and equal-luminosity star-star binary can have an arbitrarily small photometric SMA and thus be confused as a star-planet binary. Although the method outlined in this paper is designed to overcome this, the mass uncertainties on the secondary can nonetheless be large (figure 10). But how often do 'contaminating' star-star binaries occur in practice?

Here we make some simple assumptions about the expected populations of star-star and star-planet binaries to investigate statistically how likely it is that the former have photometric SMAs similar to the latter. We look at stars with masses between 0.8 and 1.5 ℳ_⊙, which corresponds roughly to F and G dwarfs.

To model both types of binaries, we draw the mass of the primary from a truncated power law $P\left({ℳ}_1\right)\propto {ℳ}_1^{-1}$. For the star-star binaries we then draw the mass ratio from a flat distribution and use this to assign the mass of the secondary. The resulting distributions are shown in the upper panel of figure B.1. For the star-planet binaries we draw the planet mass from a gamma distribution truncated between 0.5 and 5.0 ℳ_J. The resulting distributions are shown in the lower panel of figure B.1 We simulate ten times as many star-star as star-planet binaries to reflect a higher occurence of stellar companions over massive planet companions around F and G stars (Pass et al. 2023).

The orbital period for each system is drawn from a uniform distribution between 30 days and 5 years. From the period and masses we compute the physical SMA from Kepler's third law (equation 1). To compute the SMA of the photocentre (a_p) for the star–star binaries, we need to take into account the flux contributions from both components. To do this we first assign a log age and a metallicity for the system by drawing these quantities from uniform distributions with limits equal to those of the forward model grid (section 2.3). We then use the forward models to predict the flux for both components and finally equation 2 to compute a_p. The resulting distributions are shown in the upper panel of figure B.2.

Fig. B.1 Distribution of the masses and mass ratio in the simulated star-star systems (top), and the simulated star-planet systems (bottom).

Fig. B.2 Distribution of the period and SMAs in the simulated star-star systems (top), and the simulated star-planet systems (bottom).The two left panels span a common range, as do the two centre ones, but the right ones do not. The blue curve in the upper right panel shows the distribution if we erroneously neglect the flux from the secondary in computing the photocentre of the SMA.

For the star-planet binaries we assume the planet contributes no flux, so we don't need to simulate its photometry. The SMA of the photocentre is given by equation 3 with f₁/f₂ = ∞. The period and SMA distributions for these binaries are shown in the lower panel of figure B.2.

It is already apparent from the above figures that the photometric SMAs of the star-planet binaries are usually much smaller than those of the star-star binaries. This is more obvious in figure B.3, which plots them together. The bottom row shows the ratio of the number of systems of each type with photometric SMA below some limit. We see that below 0.01 au, there are only 0.037 star-star binaries for every star-planet binary. So although is possible for a near-equal mass star-star binary to have a very small photometric SMA, it is rare that it is small enough to mimic a star-planet binary. This result is quite robust to changes both in the shapes of the mass distributions adopted and to the lower period limit. Even if we now assume stars-star binaries outnumber star-planet ones by a factor of 50, and so scale up the star-star curve by a factor of five, the above fraction becomse 0.18, a 'contamination' rate of 1 in 6.5.

Fig. B.3 Comparison of the photometric SMA distributions for the simulated star-star and star-planet systems. The right panel is a zoom of the left panel in both rows. The top row shows the actual number of systems per bin (taking into account that star-planet binaries are ten times rarer). The bottom row shows the ratio of the number of star-star to star-planet systems that have a photometric SMA less than the value shown on the horizontal axis.

A more detailed analysis of the expected stellar binary contamination in the Gaia DR4 exoplanet yield is given by Lammers & Winn (2026), and is also studied in Marcussen & Albrecht (2023).

All Tables

All Figures

	Fig. 1 Sketch of two objects (1 and 2) of different masses orbiting their barycentre, with a flux ratio that places the photocentre away from the barycentre. In an unresolved binary, we observe the motion of the photocentre at about the barycentre.
In the text

	Fig. 2 CAMD of the PARSEC models used to fit the forward models of Equation (4).
In the text

	Fig. 3 Variation of MG with mass predicted by the fitted model of Equation (4), colour-coded by median metallicity.
In the text

	Fig. 4 CAMD for the Orbital300 sample assuming zero extinction.
In the text

	Fig. 5 Distribution of the data for the Orbital300 sample. Top row: the photometric SMA ap and period T. Middle row: the astrometric measurement $a_p^3/T^2$(with ap in au and T in yr, from the left side of Equation (3)) (left) and its S/N (right). Bottom row: parallax S/N and distance (limited by construction to less than 300 pc).
In the text

	Fig. 6 Inferred masses (median of the posteriors) for the Orbital300 sample for the primary (ℳ1, top) and secondary (ℳ2, bottom).
In the text

	Fig. 7 Distribution of mass ratio r = ℳ1/ℳ2 for the Orbital300 sample for three different primary mass ranges.
In the text

	Fig. 8 Variation of the mass ratio ℳ1/ℳ2 with period for the Orbital300 sample. The red points are the median values of the mass ratio in equally spaced period bins.
In the text

	Fig. 9 System CAMD for the Orbital300 sample colour-coded by the primary mass (top) and secondary mass (bottom).
In the text

	Fig. 10 Distribution of the inferred fractional 1σ-like uncertainties in the masses for the Orbital300 sample for the primaries (left panel) and secondaries (right panel). The lower fractional uncertainty cannot be larger than 1.0, whereas the upper one is unbounded. There are a small number of secondaries with upper fractional uncertainties beyond the limit of the plot.
In the text

	Fig. 11 Examples of marginalised mass posteriors for the primary (blue) and secondary (orange) for four sources (one per panel). Panel 1 is typical, panels 2 and 3 show examples of long tails to low masses for the secondaries, and panel 4 shows where some of the MCMC chains moved between the primary and secondary during the sampling. The vertical dotted lines show the median masses.
In the text

	Fig. 12 Predicted G-band flux ratio versus inferred mass ratio for the Orbital300 sample.
In the text

	Fig. 13 Flux relevance ratio (Equation (3)) plotted against the inferred mass ratio (left) for the OrbitaBOO sample shown on a logarithmic density scale. The right panel shows the distribution of the flux relevance ratio.
In the text

Fig. 14 Comparison of our mass estimates on the Orbital300 sample using astrometry and photometry with those from Gaia Collaboration (2023a), which adopted a mass for the primary and then inferred the mass of the secondary from a variety of data (in many cases there is no secondary mass estimate). The two panels show the primaries (left, 19 603 stars) and secondaries (right, 1574 stars). A few sources lie outside the range of the plots.

In the text

	Fig. 15 Deviations, as defined by Equation (14), for the results shown in Figure 14.
In the text

Fig. 16 Comparison of our mass estimates on the AstroSpectro300 sample when not using the spectroscopic measurement with mass estimates from Gaia Collaboration (2023a), which adopted a primary mass and then used the astrometry and the spectroscopy to estimate the secondary mass. The two panels show the primaries (left) and secondaries (right), both with 9674 stars. The horizontal and vertical axes do not span a common range in each panel.

In the text

	Fig. 17 Deviations, as defined in Equation (14), for the results shown in Figure
In the text

	Fig. 18 Comparison of our mass estimates using infrared photometry (vertical axis) with that not using infrared photometry (horizontal axis) on the OrbitaBOO subset for the primary (left) and secondary (right).
In the text

	Fig. 19 CAMD for the AstroSpectro300 sample assuming zero extinction.
In the text

	Fig. 20 Histograms of the distributions of the astrometric and spectroscopic measurements (left) and their S/Ns (right) for the AstroSpectro300 sample.
In the text

	Fig. 21 Comparison of our mass estimates using spectroscopy (vertical axis) with those not using spectroscopy (horizontal axis) on the Astro-Spectro300 sample for the primary (left) and secondary (right), shown on a logarithmic density scale.
In the text

	Fig. A.1 Performance of mass estimation on the simulated data set using astrometry and Gaia photometry. Note that in each panel, the horizontal and vertical axis do not span a common range.
In the text

	Fig. B.1 Distribution of the masses and mass ratio in the simulated star-star systems (top), and the simulated star-planet systems (bottom).
In the text

	Fig. B.2 Distribution of the period and SMAs in the simulated star-star systems (top), and the simulated star-planet systems (bottom).The two left panels span a common range, as do the two centre ones, but the right ones do not. The blue curve in the upper right panel shows the distribution if we erroneously neglect the flux from the secondary in computing the photocentre of the SMA.
In the text

Fig. B.3 Comparison of the photometric SMA distributions for the simulated star-star and star-planet systems. The right panel is a zoom of the left panel in both rows. The top row shows the actual number of systems per bin (taking into account that star-planet binaries are ten times rarer). The bottom row shows the ratio of the number of star-star to star-planet systems that have a photometric SMA less than the value shown on the horizontal axis.

In the text