Dynamical masses across the Hertzsprung–Russell diagram Open Access

ABSTRACT

We infer the dynamical masses of stars across the Hertzsprung–Russell (H–R) diagram using wide binaries from the Gaia survey. Gaia’s high-precision astrometry measures the wide binaries’ orbital motion, which contains the mass information. Using wide binaries as the training sample, we measure the mass of stars across the 2D H–R diagram using the combination of statistical inference and neural networks. Our results provide the dynamical mass measurements for main-sequence stars from 0.1 to 2 M_⊙, unresolved binaries, and unresolved triples on the main sequence, and the mean masses of giants and white dwarfs. Two regions in the H–R diagram show interesting behaviours in mass, where one of them is pre-main-sequence stars, and the other one may be related to close compact object companions like M dwarf-white dwarf binaries. These mass measurements depend solely on Newtonian dynamics with minimal assumptions on eccentricities, providing independent constraints on stellar evolutionary models, and the occurrence rate of compact objects.

1 INTRODUCTION

Mass is arguably the most fundamental stellar parameter of stars. Masses together with metallicities of stars determine their stellar structures, surface temperatures, luminosities, stellar and chemical evolution, lifetime, and their ultimate fates. However, individual stellar masses are often derived from stellar photometry using stellar models like isochrone fitting, and the foundation of stellar models relies on some benchmark stars that have model-independent mass measurements (Andersen 1991; Delfosse et al. 2000; Torres, Andersen & Gimenez 2010). Therefore, to validate and improve our understanding of stellar physics and models, it is critical to have model-independent mass measurements covering a wide range of the parameter space.

The orbital motions of binary stars directly reflect their mutual gravity and therefore their masses, providing the most reliable model-independent mass measurements (Popper 1980; Serenelli et al. 2021 and references therein). In particular, the individual masses can be derived for double-lined eclipsing binaries (e.g. Popper 1967; Burdge et al. 2022), where the radial velocities give the mass ratio and the eclipsing light curves provide inclination constraints. For binary stars where the components are resolved and therefore not eclipsing, if their orbital periods remain sufficiently short (≲10² yr) so that the observation can detect the orbital motion (i.e. the acceleration term), their component masses can also be derived to reasonable precision through the combination of astrometry and radial velocities (Bowler et al. 2018; Brandt, Dupuy & Bowler 2019).

Due to their rareness, the double-lined eclipsing binary method is challenging to cover a wide range of stellar parameter space. For example, since the eclipsing probability is proportional to the radii at a fixed semimajor axis, eclipsing binaries are rare among stars with smaller radii (e.g. white dwarfs and brown dwarfs), with only a few reported cases (e.g. Lodieu et al. 2015; Burdge et al. 2019; Martin et al. 2024). Hence, some other techniques have been developed to measure the masses. For instance, gravitational redshifts along with stellar radius measurements can derive the masses of white dwarfs (Falcon et al. 2010, 2012; Joyce et al. 2018; Pasquini et al. 2019; Chandra et al. 2020), where the white dwarfs are not necessarily in binary systems. A single white dwarf passing a background source can lead to an astrometric microlensing event, resulting in a precise individual mass measurement (Sahu et al. 2017). In recent years, asteroseismology has become another tool to measure the mass (see review by Chaplin & Miglio 2013), in particular for solar-type stars, giants, and pulsating white dwarfs (Hermes et al. 2017), although the model dependence of the mass measurements varies for different stellar types.

Different mass measurement methods have different systematics, and cross-validations among different methods are often difficult due to the limited availability of these measurements. Furthermore, while the double-lined eclipsing binary method provides reliable masses, it is challenging to decompose the photometry of unresolved binaries into individual components, adding hurdles in using them to constrain stellar models. Therefore, the goal of this paper is to develop a homogeneous method to measure the masses of stars for the entire H–R diagram, placing the mass measurements of different stellar populations on the same scale.

In this paper, we use resolved wide binaries to measure the dynamical masses across the H–R diagram. The high-precision astrometry from the Gaia survey (Gaia Collaboration 2016) enables the identification of millions of resolved wide binaries (El-Badry & Rix 2018; Hartman & Lépine 2020; Tian et al. 2020; El-Badry, Rix & Heintz 2021), critical to covering the entire H–R diagram. These resolved wide binaries have orbital periods longer than 10³ yr, and Gaia captures their current snapshot orbital motions without acceleration terms. Different from the double-lined eclipsing binary method that measures the mass of an individual binary, we measure the mass as a function of the H–R diagram using ensemble of wide binaries. Even with a single snapshot of the orbit, from the ensemble of wide binaries we can marginalize over the nuisance parameters, in particular the orbital orientation and eccentricity (Tokovinin 2020; Hwang et al. 2022a; Hwang, Ting & Zakamska 2022b), with proper statistical inference. Furthermore, since wide binaries are resolved, their photometry is directly associated with the component stars and can be used for stellar model comparison without photometry decomposition.

Our measurements are purely based on the Keplerian motions and do not rely on further astrophysical assumptions. Based on a similar motivation, Giovinazzi & Blake (2022) have measured dynamical masses as a function of Gaia’s absolute RP magnitudes using equal-mass (‘twin’) wide binaries (El-Badry et al. 2019). Here we use both equal-mass and non-equal-mass binaries to cover the 2D H–R diagram along the BP−RP colours and the absolute G-band magnitudes, enabling the investigation of the global mass structures in the H–R diagram from a homogeneous determination.

The paper is structured as follows. We detail our methodology and sample selection in Section 2. In Section 3, we present the mass measurements for the simulated wide binaries and the wide binaries observed by Gaia. We discuss interesting mass features in the H–R diagram in Section 4 and conclude in Section 5. We use G and M_G interchangeably for absolute G-band magnitudes. The term ‘wide binary’ and ‘wide pair’ are referred to systems where two component stars are resolved by Gaia (with typical physical separations ≳10² AU), regardless of whether any of the two components have unresolved companions or not. The term ‘single star’ refers to the system with only one star (i.e. no close companions) in Gaia’s spatial resolution unit. We use ‘component stars’ and ‘individual stars’ interchangeably for the unresolved system of a wide binary, and they can be either single stars or unresolved binaries. The component star with a brighter absolute G-band magnitude is referred to as the primary component, and the fainter one as the secondary. For wide binaries where both component stars are single stars, we call them ‘genuine wide binaries’. When specifically discussing wide binaries with an unresolved companion, we use ‘three-body system’ or ‘triple systems’. In this paper, we focus on triples consisting of an unresolved binary and a resolved wide tertiary, instead of the triples where Gaia resolves all three component stars.

2 METHODOLOGY AND SAMPLE SELECTION

2.1 Inferring mass from orbital velocities and separations

For a binary with a circular orbit, its orbital velocity v_3D (the relative velocity between two component stars) is

where a_3D is the semimajor axis of the circular binary, G is the gravitational constant, and m_tot = m₁ + m₂ is the total mass of the binary where the component masses are m₁ and m₂. The subscript ‘3D’ indicates that they are the magnitudes of full 3D vectors without any projection effects. Therefore, if v_3D and a_3D are known for a circular binary, the total mass of the binary can be derived by |$m_{\rm tot}=a_{\rm 3D} v^2_{\rm 3D} /G$|⁠.

In reality, full-3D magnitudes of v_3D and a_3D in equation (1) are observationally difficult to measure, and often only the components projected onto the plane of the sky can be well determined. Furthermore, binaries may have non-zero eccentricities, so the binary separation is a function of orbital phases, instead of a fixed value at the semimajor axis in the case of a circular orbit. Similar to equation (1), these projected quantities and phase-dependent separations are related to the system mass. Hence, we define u as:

where v is the 2D orbital velocity, s is the 2D orbital separation, and both v and s are the components projected in plane of the sky. In this paper, the orbital separations s are in units of AU, v in  km s⁻¹, and u in units of  km s⁻¹ AU^1/2. The definitions of variables and their units used in this paper are summarized in Table 1.

To infer the binary mass based on the observed distribution of u, we build the likelihood model for u, which depends on the binary mass, binary orientation, and eccentricity distribution. Since our Sun is not located at a special place for these wide binaries, the orientation of binaries should be randomly distributed with respect to the Sun. As a test, we find the direction of the projected separation vectors are randomly oriented relative to the Milky Way disc for wide binaries with angular separations >1 arcsec; below ∼1 arcsec, wide binaries’ resolvability depends on Gaia’s scanning direction and therefore the orientation of detected wide binaries depends on the scanning pattern. In this paper, we use wide binaries with angular separations >2 arcsec, and we adopt a random binary orientation in the inference.

The eccentricity distribution of wide binaries can be inferred from the ‘v–r angle’, which is the angle between the orbital velocity vector (⁠|$\vec{v}$|⁠) and the separation vector (⁠|$\vec{s}$|⁠). Using this method, Hwang, Ting & Zakamska (2022b) showed that the eccentricity distribution f_e(e) (i.e. a probability distribution of eccentricities) is close to uniform (f_e(e) = 1) for wide binaries at ∼10² AU, approximately thermal (f_e(e) = 2e) at 10³ AU, and then superthermal (f_e(e) ∝ e^1.3) at >10³ AU. Since our sample is dominated by wide binaries at ∼10³ AU, our mass inference assumes a thermal eccentricity distribution (f_e(e) = 2e), which is a theoretical distribution when the ensemble of binaries uniformly populates the phase space at a fixed energy (Jeans 1919; Ambartsumian 1937; Heggie 1975; Kroupa 2008; Hamilton 2022; Modak & Hamilton 2023). The dynamical effect of scattering and Galactic tides on wide binaries (>10³ AU) is nicely summarized in Hamilton & Modak (2023). It is possible that some populations among the wide binaries do not follow the thermal eccentricity distribution (e.g. twin wide binaries; Hwang et al. 2022a), but we reserve the modelling of different eccentricity distributions for future work.

Here we define

where the total mass m_tot is known from simulations or measurements. Throughout the paper, the values of |$\tilde{u}$| are in units of  km s⁻¹ AU^1/2  M_⊙^−1/2. Then we derive the probability distribution functions |$p(\tilde{u}|e)$| using simulations. For every eccentricity e from 0 to 1 with a step of 0.01, we simulate a large number of binaries where their orbital phases are randomly sampled in time (i.e. uniformly sampled in their mean anomaly). Their binary orientations are randomly sampled so that the directions of their angular momenta are uniformly distributed on a sphere. Our code is available on GitHub.¹|$p(\tilde{u}|e)$| is numerically normalized. Fig. 1 shows the resulting |$p(\tilde{u}|e)$|⁠, where the colour represents the value of e. For a circular orbit (purple line), its |$p(\tilde{u}|e=0)$| is non-zero up to its constant circular orbital velocity |$\tilde{u}=29.8$|⁠. Its sharp peak at |$\tilde{u}=18.4$| and the tail toward |$\tilde{u}=0$| is caused by the projection effect of s and v.

We numerically integrate |$p(\tilde{u}|e)$| to derive the case for the thermal eccentricity distribution, |$p(\tilde{u}|f_e=2e)$| (black line in Fig. 1). To accelerate the computation and to have a differentiable form for the later use of neural networks, we fit a functional form for |$p(\tilde{u}|f_e=2e)$| as

where A = 4.95 × 10⁻³, B = 2.24 × 10⁻³, C = 3.85, and |$\tilde{u}_0=36.09$|⁠. This function form is chosen simply to fit the distribution well with a minimal number of parameters. The functional fit is plotted in Fig. 1 as the dotted red line, showing that it is well fit to the simulated black line.

For binaries with different total masses, their p(u|e, m_tot) are

where the pre-factor |$1/\sqrt{m_{\rm tot}}$| ensures proper normalization. To simplify the notation, we use p(u|m_tot) ≔ p(u|m_tot, f_e = 2e) and |$p(\tilde{u}):= p(\tilde{u}|f_e=2e)$| for the case of the thermal eccentricity distribution.

One important property of u and |$\tilde{u}$| is that p(u|m_tot) and |$p(\tilde{u})$| are independent of the semimajor axis a. For an arbitrary binary, its |$u=v\sqrt{s}$| would follow p(u|m_tot) regardless of its a because the Keplerian motion is scale-free and the projection effect does not depend on a. Therefore, semimajor axes do not enter our likelihood later and there is no need to marginalize over the underlying semimajor axis distribution.

2.2 From total masses to individual masses

In the last section, we have derived the likelihood p(u|m_tot) (equation 5) so that we can infer the posterior of the total mass p(m_tot|u) for a wide binary using the Bayes’ theorem. However, how do we further constrain the masses of individual component stars? The key is that we have the photometry of individual component stars (in this paper, colours and absolute magnitudes), which provides the connection between the total mass (m_tot) and the component masses (m₁ and m₂).

For example, we can select a sample of equal-mass (twin) wide binaries where all individual components have identical colours and absolute magnitudes. Therefore, the binary mass ratios are 1, and then we can derive the individual masses from the total masses. This method is used in Giovinazzi & Blake (2022), where the authors measure the individual masses with respect to Gaia’s absolute RP magnitudes.

However, twin wide binaries only constitute <10 per cent of all wide binaries (El-Badry et al. 2019), and there are regions in the H–R diagram (e.g. white dwarfs) where twin systems are rare. Therefore, using non-equal-mass wide binaries is necessary to cover the entire H–R diagram. Indeed, non-equal-mass wide binaries can constrain the individual masses as well. For instance, say we have three types of individual stars (star |$\mathcal {A}$|⁠, |$\mathcal {B}$|⁠, and |$\mathcal {C}$|⁠), and there are three pairs of wide binaries, each consisting of star |$\mathcal {A}$| + |$\mathcal {B}$|⁠, |$\mathcal {B}$| + |$\mathcal {C}$|⁠, and |$\mathcal {A}$| + |$\mathcal {C}$|⁠. Once we obtain the total masses of each wide binary from their orbital velocities and separations, |$m_{\rm tot,1}=m_\mathcal {A}+m_\mathcal {B}$|⁠, |$m_{\rm tot,2}=m_\mathcal {B}+m_\mathcal {C}$|⁠, and |$m_{\rm tot,3}=m_\mathcal {A}+m_\mathcal {C}$|⁠, the three individual masses can be solved by |$m_\mathcal {A}=(m_{\rm tot,1}+m_{\rm tot,3}-m_{\rm tot,2})/2$|⁠, etc. This example demonstrates that non-equal-mass binaries are able to constrain the individual component masses. In reality, this procedure is more complicated because we do not deterministically measure the total mass of each wide binary. Instead, we have p(m_tot|u) for every wide binary, and we forward-model the component masses through proper statistical inferences that maximize the likelihood of the observations.

Now we are ready to formulate the statistical model. For each wide binary with an index i, we measure its |$u_i=v_i \times \sqrt{s_i}$| (equation 2), where v_i is the projected orbital velocity and s_i is the projected binary separation. Then for a sample of wide binaries, we have a set of {u_i} from observations. Our goal is to find a function |$\mathcal {M}(...)$| that maps the observables of an individual star to its individual mass. To derive the mass across the H–R diagram, we consider the BP−RP colour and the absolute G-band magnitude measured by Gaia as the observables in this paper. Therefore, |$m_{i,j} = \mathcal {M}(({\rm BP-RP})_{i,j}, {\rm G}_{i,j})$|⁠, where m_{i, j} is the individual mass of the j-th (j ∈ {1, 2}) component star of the i-th wide binary, and (BP−RP)_{i, j} and G_{i, j} are the colour and the absolute G-band magnitude of the component star. Then for the i-th wide binary, its total mass is m_{tot, i} = m_{i, 1} + m_{i, 2}. The component mass we are measuring is the mass sum within Gaia’s spatial resolution unit. For example, if the 1-st component star is an unresolved binary, then m_{i, 1} would be the total mass of the unresolved binary.

Since each wide binary is independent, we can write the likelihood as

where |$m_{\rm tot,i}=m_{i,1}+m_{i,2}=\mathcal {M}({\rm (BP-RP)_{i,1}}, {\rm G}_{i,1})+\mathcal {M}({\rm (BP-RP)_{i,2}}, {\rm G}_{i,2})$|⁠. p(u_i|m_{tot, i}) can be evaluated using equation (5).

Observationally there are measurement uncertainties associated with u. Therefore, we marginalize over the uncertainty:

where u is the observed value and u′ the true value, and the uncertainty distribution p(u_i|u′) is an assumed Gaussian distribution with a width of the associated uncertainty |$\sigma _{u_i}$|⁠. Here we assume a flat prior for p(u′) so we can approximate p(u_i|u′) using p(u′|u_i). We discuss |$\sigma _{u_i}$| from the data perspective in Section 2.4.

Our formulation of |$\mathcal {M}$| as a function of BP−RP and the absolute magnitude does not contain any astrophysical assumptions. For example, we do not assume any mass–temperature or mass–luminosity relation. One can include other observables of their interest, like other photometry bands and metallicity, in the parameters of |$\mathcal {M}$|⁠. What we assume in the formulation is that the mass is well-behaved as a function of BP−RP and the absolute magnitude, and the choice of these two observables is indeed astrophysically motivated. Other than this motivation, we emphasize that our likelihood model for mass measurements is purely based on Newtonian dynamics without inputs from astrophysical models.

Our approach is applicable if there is a mass distribution (instead of a single mass) at a certain location of the H–R diagram. Astrophysically, a particular location in the H–R diagram may have a mass distribution due to the crossing of multiple isochrones, e.g. different stellar ages, metallicities, unresolved stellar companions, and the presence of compact objects. In that case, the mass solution of the function |$\mathcal {M}$| would be close to the arithmetic mean [with some difference depending on the function |$f(\tilde{u})$|] of the underlying mass distribution so that the mass solution maximizes the likelihood in equation (6). In other words, when there are underlying mass distributions, then we are measuring the average mass for every point of the H–R diagram.

2.3 Neural network model for |$\mathcal {M}$|

There are several approaches to formulate the function |$\mathcal {M}($|BP−RP, G). For example, one can bin the H–R diagram and assign a mass to each bin, i.e. |$m[k]=\mathcal {M}({\rm BP-RP}, {\rm G})$| if (BP−RP, G) is in the k-th bin. We use this formulation along with the Markov Chain Monte Carlo (MCMC) method to derive the posterior distribution of the mass measurements in Appendix  A. However, this method requires pre-defined bins on the H–R diagram, and can only work with a small number of parameters due to its expensive computation.

Our goal is to map the mass across the H–R diagram without priori knowledge like pre-defined bins. We end up using a neural network to construct the target function |$\mathcal {M}({\rm BP-RP}, {\rm G})$|⁠. The advantage of using a neural network is that it does not require pre-defined bins on the H–R diagram. Instead, a neural network takes the two input parameters (BP−RP and G), and outputs the component mass after a series of algebraic transformations. The mass model from the neural network is trained by optimizing all the hyperparameters (weights and biases) in the network. Based on the Universal Approximation Theorem (e.g. Cybenko 1989; Funahashi 1989), a neural network can theoretically approximate any continuous functions, thus providing the most general formulation to map the observables to the mass.

We use the neural network implementation from PyTorch (Paszke et al. 2019; v2.0.1). We use a 5-layer fully connected neural network, with 128 neurons in each layer. The Gaussian error linear units (GELU) are used for the activation function. Dropout is used during training and prediction so that each neuron has a 20 per cent chance of being dropped out. Dropout avoids overfitting but also provides an estimate for the model prediction uncertainty based on the dropout variational inference (Kendall & Gal 2017; Leung & Bovy 2019). The neural network takes two inputs: the colour (BP−RP) and the absolute magnitude (G) of a star. To ensure that the prediction mass is always positive, the neural network’s output is the mass exponent x, and the mass is m = exp (x) M_⊙. During the training, there are 1024 wide binaries in each training batch, and the total number of the training sample is ∼0.1 million. Different from traditional neural networks which often use the L₂ (mean squared errors) loss function, we use the likelihood [equation (6) with the outlier mixture model introduced in the next section] as the loss function and then adjust the model parameters by computing the gradient of the loss function with respect to the parameters (‘back-propagation’). The Adam optimizer is used (Kingma & Ba 2015). The training typically goes through the entire sample 500–1000 times (‘epochs’). During the prediction, we use the trained model to generate 1000 predicted masses for a given star with dropout activated, and we use the median of the predicted masses as the reported mass and the standard deviation as the reported mass uncertainty. The training is done on the Google Colab platform with graphics processing units (GPU) to accelerate the computation. The training of 1000 epochs takes about 15 min using a V100 GPU.

2.4 Wide binaries from Gaia

We use the wide binary catalogue from El-Badry, Rix & Heintz (2021), where the resolved wide binaries were searched using the Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2016, 2018, 2021) out to 1 kpc from the Sun. These wide binaries were identified such that the component stars of a binary have parallaxes and proper motions consistent with being a gravitationally bound system. Wide binaries in clustered regions and resolved triples (where all three components are resolved by Gaia) were excluded from the catalogue. For each wide binary, El-Badry, Rix & Heintz (2021) estimated the probability of being a chance-alignment pair (R_chance) by doing the wide binary search at random coordinates where physical wide binaries are not expected. We adopt R_chance < 0.01 to avoid chance-alignment pairs, and the separations of the wide binaries of our interest are ∼10³ AU where the contamination from chance-alignment pairs is negligible.

For each Gaia wide binary, we measure its projected physical separation s_i (in AU) and projected orbital velocity v_i (in  km s⁻¹) to compute |$u_i=v_i\times \sqrt{s_i}$| [equation (2)]. The projected separation is s_i = 1000θ_i/ϖ_i, where the angular separation θ_i is in arcsec and the parallax ϖ_i is in mas. Following El-Badry, Rix & Heintz (2021), we use the parallax of the primary (the brighter component) of the wide binary to avoid possible parallax systematics in fainter sources. We require parallaxes over error >10. The orbital velocity is v_i = 4.74 × Δμ_i/ϖ_i, where the proper motion difference is |$\Delta \mu _i=|\vec{\mu }_{i,2}-\vec{\mu }_{i,1}|$| (in  mas yr⁻¹) and |$\vec{\mu } = (\mu _{\alpha }^{*}, \mu _{\delta })$|⁠. The uncertainty of proper motion difference is

where |$\sigma _{\mu ^{*}_{\alpha ,1}}$| and |$\sigma _{\mu _{\delta ,1}}$| are the uncertainty of |$\mu ^{*}_{\alpha ,1}$| and μ_{δ, 1}, |$\Delta \mu ^2_\alpha =(\mu ^{*}_{\alpha ,2}-\mu ^{*}_{\alpha ,1})^2$|⁠, and |$\Delta \mu ^2_\delta =(\mu _{\delta ,2}-\mu _{\delta ,1})^2$|⁠. The subscripts ‘1’ and ‘2’ indicate the two component stars of the wide binary. Our error analysis could be improved by considering the correlation between the two proper motion components. However, as most wide binaries in Gaia have a small correlation, we stick to this simpler error form.

The measurement of u depends on angular separations, parallaxes, and proper motion differences. Typically, angular separations can be determined better than 0.1 per cent precision, and parallaxes better than a few per cent precision. Therefore, the uncertainty of u is dominated by the proper motion difference, and we compute the uncertainty of u by σ_u = uσ_Δμ/Δμ. Based on Gaia’s current proper motion precision, we select wide binaries with s < 10^3.5 AU and parallaxes >2.5  mas (distances <400 pc). In particular, a circular wide binary with a total mass of 1 M_⊙ at a = 3000 AU has an orbital velocity of 0.55 km s⁻¹, corresponding to a proper motion difference of 0.3 mas yr⁻¹ at parallaxes of 2.5 mas. For typical proper motion uncertainties of 0.1 mas yr⁻¹ in Gaia EDR3, this sample selection ensures that most of the wide binaries can have signal-to-noise ratios (SNR) of proper motion difference |$\Delta \mu /\sigma _{\Delta _\mu }\gt 3$|⁠, and therefore u/σ_u > 3.

Observational data often exist with some outliers from various origins. For instance, if a star has a close stellar companion with orbital periods of several months, it may affect the astrometric measurements derived from the single-star model (Belokurov et al. 2020; Penoyre, Belokurov & Evans 2022a, b). Indeed, the majority of the resolved tertiaries of inner astrometric binaries (identified by Gaia Data Release 3; Halbwachs et al. 2023) are missed if the wide binary search uses the single-star solutions (Hwang 2023). In this paper, we only use the single-star astrometric solutions from Gaia EDR3 (which are unchanged in DR3). To ensure that the single-star astrometric measurements are reliable, we require that both component stars have renormalized unit weight errors ruwe<1.4 (Lindegren et al. 2018). Note that Gaia DR3 only processes the astrometric binary solutions for stars having ruwe>1.4 (Halbwachs et al. 2023), so our ruwe<1.4 selection effectively excludes the astrometric binaries identified in Gaia DR3. However, astrometric binaries may still be present at ruwe<1.4 (Penoyre, Belokurov & Evans 2022b); therefore our sample may still be subject to possible outliers due to unresolved companions.

Since |$p(\tilde{u})$| has a sharp cutoff at |$\tilde{u}\sim 40$| (Fig. 1), the likelihood model from equation (6) is vulnerable to high-u outliers. A high-u outlier would have an expensive penalty in the likelihood, thus overestimating the mass. To appropriately incorporate the high-u outliers, we adopt a mixture model (e.g. Hogg, Bovy & Lang 2010):

where p_i(u_i|m_{tot, i}) is equation (5) for well-behaved data, and F is the fraction of data that is better described by the outlier model |$p_{\rm outlier}(\tilde{u})$|⁠. We adopt a Gaussian (truncated at 0) for the outlier model:

where Z is the normalization factor (not exactly |$\sqrt{2\pi \sigma _{\rm outlier}^2}$| due to the truncation at 0). After trial-and-error with the Gaia data, we use |$\tilde{u}_{\rm outlier}=30$|⁠, σ_outlier = 20, and F = 0.20. We also use MCMC to constrain F from the data, ending up with a similar value. Our final log-likelihood [see equation (6)] becomes

The BP and RP fluxes in Gaia (E)DR3 are measured without deblending treatment, so a nearby source (wide binary companion or physically unrelated source) within ∼2 arcsec would affect one star’s BP and RP fluxes. Also, wide binaries with angular separations <2 arcsec have significantly higher ruwe likely due to centroiding errors (El-Badry, Rix & Heintz 2021). Therefore, we require wide binaries to have angular separations >2 arcsec, SNR of G/BP/RP fluxes >10, and each component star have bp_rp_color_excess_factor <1.8 to avoid unreliable BP−RP due to the nearby sources (Evans et al. 2018). We limit our sample to galactic latitudes |b| > 10 deg to avoid the dustier region in the thin disc, and we do not apply extinction correction for photometry. The photometric uncertainties play a minor role, and thus we do not include photometric uncertainties in our statistical model.

Fig. 2 (left) shows the distribution of primaries and secondaries from the resulting 99 680 pairs of Gaia wide binaries. The mean projected separation is 1.4 × 10³ AU. These wide binaries contain main-sequence (MS) stars, (sub-)giants, and white dwarfs, covering a wide range of the parameter space in the H–R diagram for mass measurements. The right panel of Fig. 2 shows the median u in each bin of the left panel. Since |$u\propto \sqrt{m_{\rm tot}}$| (equation 2), u is positively correlated with the mass in each bin, assuming the wide companion mass distribution is similar across the H–R diagram (which is a reasonable assumption but not strictly correct; El-Badry et al. 2019). Therefore, u in the right panel represents several features of mass, including the higher masses in blue main-sequence stars, the lower masses in red main-sequence stars, and the unresolved binaries parallel to the main sequence. This plot is a convenient way to present the data, but every u measurement is actually associated with two stars at different locations (unless it is a twin) in the H–R diagram. Therefore, it is necessary to use the mass inference model introduced in the earlier sections to recover the underlying masses across the H–R diagram.

3 MASS MEASUREMENTS ACROSS THE H–R DIAGRAM

3.1 Tests on mock wide binary samples

Here we demonstrate that our method can recover the underlying masses from the mock wide binary samples where their masses are known. Specifically, we simulate 100 000 mock wide pairs, similar to the number of Gaia wide binaries in our sample. Among the mock wide binaries, 50 per cent are genuine wide binaries without any unresolved companions, 40  per cent are triples consisting of one single star and one unresolved binary, and 10 per cent are quadruples consisting of two unresolved binaries. 20 per cent of the genuine wide binaries (i.e. 10 per cent of all pairs) are equal-mass binaries. The component masses are randomly drawn from the Salpeter initial mass function (Salpeter 1955) with a mass range between 0.2 and 1 M_⊙. We use brutus² (Speagle et al., in preparation) to generate their Gaia photometry (G, BP, RP) for single stars and unresolved close binaries based on the MIST isochrones (Choi et al. 2016; Dotter 2016), assuming solar metallicities and stellar ages of 10^9.5 yr. For this test, we only consider main-sequence stars. Then for each wide pair, we randomly draw a |$\tilde{u}_i$| from the distribution |$p(\tilde{u})$| assuming that they follow the thermal eccentricity distribution (equation 4), and scale it to |$u_i=\tilde{u}_i \sqrt{m_{{\rm tot},i}}$|⁠, where m_{tot, i} is the total mass (in  M_⊙) of the wide pair system. Here we assume no measurement uncertainties for u_i and no outlier measurements. We train the neural network for 1000 epochs with a learning rate of 0.001. The tests aim to infer the mass across the H–R diagram from the information of {u_i} and their photometry ((BP−RP)_i, G_i).

We start with a toy model where we manually set the masses of stars to discrete values such that their masses only depend on the BP−RP colours, regardless of whether they have unresolved companions or not. Specifically, the mass is 2.5–0.5 k M_⊙ for systems having 0.5 + 0.5k ≤BP−RP<1 + 0.5k where k ∈ [0, 1, 2, 3, 4]. The ground-truth masses of the mock sample are shown in Fig. 3, left panel of the first row. The double-track feature at BP−RP∼2.5 is because the minimum main-sequence mass prohibits some parameter space of unresolved non-equal-mass binaries. Then we apply our neural-network-based method in Section 2 (and set the outlier fraction F = 0) to derive the mass. The second panel shows that the inferred masses recover the underlying masses from the mock wide binaries. The third panel is the difference between inferred and true mass. Most stars have differences close to zero, except the boundaries of discrete masses have larger deviations up to ∼0.3 M_⊙. This is a general problem that a continuous function is difficult to fit discontinuous boundaries. In the fourth panel, stars close to the boundaries have larger mass uncertainties estimated using the dropout variation, suggesting that the uncertainty estimates capture the underlying uncertainty (see also the 5th panel about the fractional errors). The last panel shows the resulting |$\tilde{u}$| distribution using the mass measurements from the training (blue), in agreement with the likelihood model |$p(\tilde{u})$| from equation (4) (black line).

Fig. 4 is the second test using the astrophysical masses from the MIST isochrones. The mass here is the mass within the spatial resolution limit, which can be the mass of a single star or the mass of an unresolved binary. The first panel shows the ground-truth masses. Astrophysically, at a fixed metallicity, the masses of stars determine the surface temperatures and therefore BP−RP colours, and the dependence of mass on absolute magnitudes at a fixed BP−RP is due to the mass ratios of unresolved binaries. The second and the third panels show that the inferred masses agree well with the true mass. There are larger differences in the higher mass end because there are fewer high-mass stars due to the initial mass function. The fourth and fifth panels demonstrate that the uncertainty estimates also report reasonably higher uncertainties at the high-mass end. The last panel shows that the resulting |$\tilde{u}$| using the inferred masses (blue) agrees well with the likelihood model from equation (4) (black line).

The tests in Figs 3 and 4 suggest that our neural network method can measure the dynamical masses from wide binaries across the entire H–R diagram. The tests also show that the measurement may be less reliable at some discontinuous mass transition (if any), and our method would report a larger uncertainty there. The size (100 000) of our mock binary sample is similar to our Gaia sample, demonstrating that we can obtain rich information about mass from the actual Gaia sample.

3.2 Masses across the H–R diagram using Gaia wide binaries

Among the 99 680 pairs of Gaia wide binaries from Section 2.4, we randomly select 90 per cent of them as the training sample and 10 per cent as the validation sample. We then train our models using a learning rate of 1 × 10⁻⁴ and an outlier fraction of F = 0.20. During the training, the average loss of the validation sample stops decreasing at epochs around 400 and becomes flat afterwards. Therefore, we stop the training at 500 epochs to have a good performance in the validation sample and also avoid possible overfitting in longer epochs.

Fig. 5 shows the resulting inferred dynamical masses across the Hertzsprung–Russell Diagram, where each point is a component star in the wide binary sample. Here we use these component stars to present the model, but we emphasize that the model itself is a continuous function across the H–R diagram (function |$\mathcal {M}$| in Section 2.2), and there is no discrete binning involved. Only components with fractional errors <0.5 (99.98 per cent of the sample) are plotted. Fig. 6 is the inferred fractional error, showing that the masses can be measured better than 10 per cent in most regions of the H–R diagram. The Sun has M_G = 0.82 and BP−RP = 4.67 (Casagrande & VandenBerg 2018), where our measurement gives 0.96 ± 0.08 M_⊙ (the error represents the uncertainty of an individual object at Sun’s BP−RP and M_G). Our training result gives a mass of 0.67 ± 0.07 M_⊙ for white dwarfs at M_G = 12 and BP−RP = 0, including all spectral types of white dwarfs. This mass is well consistent with the gravitational redshift measurements, where the mean mass of hydrogen-atmosphere (DA) white dwarfs is |$0.647^{+0.013}_{-0.014}$| (Falcon et al. 2010) and that of helium-atmosphere (DB) white dwarfs is |$0.74^{+0.08}_{-0.09}$| M_⊙ (Falcon et al. 2012).

Due to the equation of state of degenerate electrons, more massive white dwarfs are smaller in size than the less massive ones, and thus more massive white dwarfs are dimmer at fixed colours. Therefore, white dwarfs formed from single-star evolution are expected to have masses ∼0.5 M_⊙ on the bright end and ∼1 M_⊙ at the faint end, which is observationally confirmed through gravitational redshift measurements (Chandra et al. 2020). In contrast, our result is consistent with a nearly constant mass around 0.70 M_⊙ for all white dwarfs. There is a weak mass trend that brighter white dwarfs are more massive by ∼0.07 M_⊙ than the faint white dwarfs, but this difference is only at 0.6-σ and thus is not statistically significant. The reason we do not recover the expected mass gradient in white dwarfs may be that the number of white dwarfs is insufficient in the sample, only ∼5 per cent of the component stars. Furthermore, the mass distribution of white dwarfs strongly peaks around 0.6 M_⊙(Tremblay et al. 2016), and the higher mass white dwarfs are rarer in wide binaries due to their significant mass-loss during the late stellar evolution (Hwang et al., in preparation). Hence, there are too few high-mass white dwarfs in the wide binary sample to detect their higher masses.

Fig. 7 shows the distribution of |$p(\tilde{u})$| from the training result (blue). Different from the theoretical |$p(\tilde{u})$| (orange line), the observed |$p(\tilde{u})$| has a main component and an extended high-|$\tilde{u}$| tail up to |$\tilde{u} \sim 80$|⁠. The main component agrees well with the |$p(\tilde{u})$| from the thermal eccentricity distribution (orange), supporting our observationally motivated assumption that the eccentricity distribution follows the thermal eccentricity distribution (Hwang, Ting & Zakamska 2022b). By construction of the model, the overall behaviour of the high-|$\tilde{u}$| tail is well described by our assumed outlier model |$p_{\rm outlier}(\tilde{u})$| with an outlier fraction of F = 0.20 (Section 2.4). The high-|$\tilde{u}$| tail up to |$\tilde{u} \sim 80$| cannot be explained by eccentricities because |$p(\tilde{u}|e)$| is strongly truncated above |$\tilde{u}=40$| for all e (Fig. 1) due to the small fraction of time spent at pericentre. The high-|$\tilde{u}$| tail remains similar if we only select wide binaries with smaller errors of u (e.g. u/σ_u > 20). If we select wide binaries where both components have apparent G-band magnitudes brighter than 16 mag, then the high-|$\tilde{u}$| tail is reduced to F ∼ 10 per cent, suggesting that a significant fraction of the high-|$\tilde{u}$| tail is associated with fainter sources.

The high-|$\tilde{u}$| tail can be due to either data stochastic outliers or physical causes. For example, it can be that the astrometric quality is worse in fainter sources and is not well characterized in Gaia’s reported uncertainties. Alternatively, our result assumes a one-to-one relation between the mass and the H–R diagram, and physical scenarios that do not follow this assumption may become outliers, for example the presence of compact object companions and the metallicity dependence of isochrone. Also, some outliers may be due to the astrometric noise caused by the orbital motion of the unresolved companion (Penoyre, Belokurov & Evans 2022b). Even though the nature of the high-|$\tilde{u}$| tail remains unclear, this paper aims to use a reasonable outlier model to take outliers into account, so that we can measure the mass through the main component reliably.

Several works have suggested that the orbital velocities of wide binaries at ∼10⁴ AU deviate from the expected Keplerian velocities, indicating the possible presence of non-Newtonian gravity in the low-acceleration regime (Banik & Zhao 2018; Pittordis & Sutherland 2019; Chae 2023; Hernandez 2023, but see alternative explanations in El-Badry 2019). In Fig. 7, we show that at ∼10³ AU wide binaries, there is already a tail of higher-than-Keplerian velocities. Therefore, they may be of the same origin (e.g. worse astrometric quality in fainter sources, astrometric noise from unresolved companions, etc.), and it is just that at ∼10⁴ AU binary separations, the Keplerian orbits become subdominant than the noise, making the observed velocities ‘non-Newtonian’. Future work is necessary to compare the difference in high-velocity tail behaviours at different wide binary separations.

Based on similar merit but quite different methods, Giovinazzi & Blake (2022) measure the mass with respect to absolute RP magnitudes using equal-mass wide binaries. In Appendix  A, we compare our results with theirs after photometry conversion. Our trained models and the interpolated grids are available on GitHub.³ Due to the stochastic nature of the training/validation sample selection and the neural network training, the mass values may not be exactly the same from different runs of training, but the differences are consistent within the inferred mass uncertainty. Therefore, the overall mass structures in the H–R diagram and our main results are unchanged in different training runs.

4 DISCUSSION

4.1 Wide binaries as representatives of field stars

We derive the dynamical mass across the H–R diagram from wide binaries, but is this relation applicable for general field stars? If wide binaries are a representative sampling of field stars, then the derived relation between mass and the H–R diagram applies to the field stars as well. Then the question becomes whether stars in wide binaries differ from the field stars that are not in wide binaries.

Metallicity plays an important role in determining an isochrone. We find that the wide (10³–10⁴ AU) binary fraction peaks around the solar metallicity (Hwang et al. 2021, 2022c). Therefore, the mean metallicity of wide binaries is similar to that of the field stars in the solar neighbourhood, and they are expected to follow a similar isochrone.

Field stars and stars in wide binaries have significantly different occurrence rates of close companions. Compared to stars without resolved wide companions, close binaries with orbital periods <10 d are more common among wide pairs (i.e. wide tertiaries) (Hwang et al. 2020; Fezenko, Hwang & Zakamska 2022; Hwang 2023). However, this difference does not affect our mass measurements as long as the close binaries with wide companions follow the same mass–HR relation as the close binaries without wide companions. There is a small subset of close binaries (e.g. contact binaries) that have interaction between two component stars and thus may change their photometry (Qian 2003) so they do not follow the same mass–HR relation; however, the occurrence rate of contact binaries is small (∼0.1 per cent; Paczynski et al. 2006; Kirk et al. 2016; Hwang & Zakamska 2020), leaving a negligible effect on our measurements.

Non-interacting close companions of compact objects (black holes, neutron stars, and faint white dwarfs) strongly alter the system’s total mass without changing the overall photometry. If the occurrence rate of compact object close companions is higher in wide binaries (thus forming a three-body system) than in single stars (i.e. those without wide companions), then our measured mass would be higher than that of a single star due to the mass contribution from compact objects. However, we expect the occurrence rate of compact objects to be low; hence, this potential difference in compact object occurrence rate likely has a small effect on our results. Furthermore, as discussed in Section 4.3 later, the mass measurement of single-like stars provides a unique constraint on the occurrence rate of compact objects.

Therefore, we consider wide binaries as reasonable representatives of field stars, and our derived mass measurements overall apply to general field stars. Our approach assumes an approximately one-to-one relation between mass and the H–R diagram. According to stellar physics, this assumption is only true for low-mass (≲ 1 M_⊙) main-sequence stars with the same metallicity and for white dwarfs. Given that the solar-neighbourhood is dominated by solar–metallicity stars, our mass measurements are mainly the mass relation at solar metallicities, although it is possible that some part of the H–R diagram is dominated by stars with different metallicities (e.g. the fainter part of the main sequence may have lower metallicities; Gallart et al. 2019; Sahlholdt, Casagrande & Feltzing 2019). In the regions of subgiants and giants in the H–R diagram, the isochrones with different ages are overlapped, and our measurements represent their masses in an average sense.

4.2 Masses of different stellar types

Gaia has opened a new era of understanding the global structures of the H–R diagram. Several features are directly identified based on the density of the H–R diagram, including the gap where MS stars transition from fully to partially convective (‘Jao’s gap’, Jao et al. 2018), and the anomalous accumulation of white dwarfs (‘Q track’) in their cooling track likely due to the crystallization at their cores (Cheng, Cummings & Ménard 2019; Tremblay et al. 2019). In this work, we provide a different view—the global mass structures of the H–R diagram (Fig. 5 and the annotated version in Fig. 12). For example, at (BP−RP, G)=(2.5, 10.1), approximately the Jao’s gap (Jao et al. 2018), our measured mass is 0.45 ± 0.03 M_⊙, close to (at 2.7σ) the theoretical transition mass of 0.37 M_⊙ for solar-metallicity stars to become fully convective (Spada et al. 2017). We note that the M-dwarf stellar models are strongly subject to the convection treatment, and thus the potential 2.7σ difference may provide an additional constraint on the stellar models.

Fig. 8 shows our measurements of single-MS stars compared with isochrone models. In the left panel, we fit a single-MS track (blue) which is a running median of absolute G-band magnitudes as a function of BP−RP (after excluding white dwarfs). Then we evaluate the mass measurements along this single-MS track in the right panel, with the width indicating ±1σ uncertainties. We then compare the results with the isochrones from PARSEC (version 1.2S, Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014; Chen et al. 2015) and from MIST (version 1.2 with rotation v_initial/v_critical = 0.4, Paxton et al. 2011, 2013; Choi et al. 2016; Dotter 2016). We consider solar-metallicity isochrones at stellar ages of 0.1 Gyr. Using isochrones with ages of 1 or 10 Gyr has negligible effects in the right panel, except that the high-mass (>1 M_⊙) stars are no longer on the main sequence. The right panel shows that our measurements are in good agreement with the isochrone model. At the low-mass end (BP−RP>2, or stellar masses <0.5 M_⊙), the PARSEC model is significantly more massive than MIST at a fixed BP−RP, and our measurements agree better with the PARSEC isochrone than MIST. Choi et al. (2016) discuss the difference between PARSEC and MIST in the low-mass dwarfs (their figs 15–17), suggesting that the difference may arise from the fact that PARSEC adopts an opacity-related boundary condition calibrated to observations (Chen et al. 2014). Therefore, the modelling of low-mass stars remains an open question, and our mass measurements provide further model-independent constraints on the stellar models.

For reference, we fit the mass-BP−RP relation of single-MS stars (blue line in the right panel of Fig. 8) using 15-order polynomials. The coefficients from the highest to lowest order are: (1.923, −1.918, 6.907, 0.839, −161.598, 648.392, −1278.028, 1542.543, −1240.595, 692.045, −272.217, 75.397, −14.405, 1.809, −0.134, 0.004). The input is BP−RP (mag) and the output is mass (M_⊙). The fit is valid for BP−RP between 0 and 4.

Fig. 9 presents our results at the unresolved binary (and triple) track. An equal-mass binary would have the same BP−RP colours as the single star, with a larger absolute magnitude by ΔG = 0.75 mag. Therefore, we obtain the orange line in the left panel of Fig. 9 by offsetting the single-MS line in Fig. 8 by ΔG = 0.75 mag. If unresolved equal-mass binaries dominate the sample at the ΔG = 0.75 mag track (orange line), then we would expect the measured mass to be a factor of two that of the single-MS track.

In the right panel of Fig. 9, we compare the mass measurements on the ΔG = 0.75 mag track with the single-MS track. The dark blue line is the mass of the single-star track and the light-blue line shows their double, the expected mass for the equal-mass binary. At BP−RP>2.5, the measured mass of the ΔG = 0.75 mag track (orange) is consistent with two times the single MSs (light blue), suggesting that the equal-mass binaries dominate the sample. At BP−RP<1, the masses at the ΔG = 0.75 mag track are consistent with those of the single-star track. This may be because the equal-mass binaries are less common in more massive stars and wide binaries (Moe & Di Stefano 2017; El-Badry et al. 2019), and also the G = 0.75 mag track overlap with the evolution paths of (sub-)giants, making our measurements more similar to single-star masses.

Interestingly, at BP−RP∼2 in the right panel of Fig. 9, the measured mass is larger than the expected equal-mass binaries by 0.6 M_⊙ at 3σ significance. An unresolved two-body system cannot explain the 0.6 M_⊙ excess and the flux excess of ΔG = 0.75 mag. Alternatively, the mass excess at this region of the H–R diagram may suggest that there is a significant fraction of unresolved triple systems.

Fig. 10 investigates more detail into the mass excess. The green line shows the measured mass versus absolute G-band magnitudes at a fixed BP−RP=1.6 (the green vertical line in Fig. 9, left panel). The horizontal axis of Fig. 10 is the inverted absolute G-band magnitude so that the right-hand side of the plot is brighter. The star symbols are the main-sequence masses from the PARSEC isochrones that intersect BP−RP = 1.6. We consider two metallicities, a solar metallicity [Fe/H] = 0 as the open star symbols and [Fe/H] = −0.35 as the filled star symbols. The single-star symbols indicate the mass and absolute magnitude for the single-star cases, the double-star symbols for unresolved equal-mass binaries, and the triple-star symbols for the unresolved equal-mass triple. At M_G = 6.6, our inferred mass reaches the highest value of 1.60 ± 0.15 M_⊙, exceeding the mass of equal-mass binaries for both metallicities. Since equal-mass binaries are the upper limit of mass for unresolved binaries at a fixed BP−RP (Fig. 4), the highest mass of 1.60 ± 0.15 M_⊙ can only be explained by unresolved triples. Having a compact object companion may explain the higher mass, but cannot explain the brighter absolute magnitude. Indeed, the highest mass and its M_G agree well with the unresolved equal-mass triples at [Fe/H] = −0.35, where the lower metallicity may be explained by the higher close binary fraction in lower metallicity stars (Raghavan et al. 2010; Yuan et al. 2015; Badenes et al. 2018; Moe, Kratter & Badenes 2019). Alternatively, the highest mass of 1.60 ± 0.15 M_⊙ may be non-equal-mass triples with solar metallicities, which can fill the parameter space between the solar-metallicity equal-mass binary (open double-star symbol) and equal-mass triple (open triple-star symbol) in Fig. 10.

Figs 9 and 10 suggest that there may be a population of unresolved triples near M_G = 6.6 and BP−RP=1.6. Recently, hundreds of compact triple systems (where the outer tertiaries have orbits smaller than a few AU) and compact quadruples have been identified (Borkovits et al. 2016, 2020a, b; Czavalinga et al. 2023; Kostov et al. 2023; Rappaport et al. 2023). The formation of compact triples is particularly challenging to theory because the tertiary orbit is smaller than the radius of an initial hydrostatic stellar core (Larson 1969). Therefore, our results may unveil the mysterious population of compact multiples and provide their total mass measurement.

Our mass measurements reveal two regions (A and B in Fig. 11) with interesting mass behaviours. Region A is located at the faintest part of the main sequence with BP−RP between 1.5 and 3. Intriguingly, stars in Region A have higher masses than the M dwarfs at similar BP−RP. Region A has a significantly higher u in Fig. 2 (right panel), in agreement with our higher inferred masses. PARSEC single-star isochrones with a low metallicity of [Fe/H] = −1.5 intersect Region A in the H–R diagram. However, a lower metallicity isochrone would have a lower mass than the solar-metallicity isochrone at the same BP−RP, inconsistent with the observed higher mass in Region A. The mass of unresolved metal-poor binaries is still insufficient to explain the observed higher mass. Furthermore, we expect all wide binaries in our sample to be co-natal, i.e. two component stars of a wide binary have the same stellar ages and metallicities (Hawkins et al. 2020; Nelson et al. 2021); therefore, the component stars of the same wide binary would be located on the same isochrone. However, several wide binaries associated with Region A (e.g. the blue solid line in Fig. 11; its primary’s Gaia DR3 source ID is 650350652407376768) are not aligned with a single metal-poor isochrone, suggesting that metallicity is not the cause for stars being in Region A. Lastly, metal-poor wide binaries are rare and are unlikely to dominate the sample in Region A (Hwang et al. 2021, 2022c). Therefore, the origin of Region A’s higher masses is not due to lower metallicities.

One possible explanation is that Region A is associated with close compact object companions. For instance, unresolved M dwarf-white dwarf binaries are located in a similar H–R diagram region (Rebassa-Mansergas et al. 2021). The unresolved white dwarf companion thus results in the higher inferred mass. This explains why some wide binaries in Region A are not aligned with any single-star isochrone, because the photometry is significantly affected by both white dwarf and M dwarf. We notice that many wide binaries associated with Region A have large u (e.g. the blue solid line in Fig. 11 has u = 117) and large |$\tilde{u}$|⁠, and Region A has a larger fraction of wide binaries associated with the high-|$\tilde{u}$| tail (Section 3.2 and Fig. 7). This may be because the unresolved M dwarf-white dwarf binaries in Region A have a large spread of the mass, and systems involving high-mass white dwarfs become the high-|$\tilde{u}$| tail. The fainter absolute magnitudes of Region A may also be the reason why the high-|$\tilde{u}$| tail is more related to the fainter sources, as discussed in Section 3.2. Future investigations are needed to confirm the nature of Region A and its connection with the high-|$\tilde{u}$| tail.

In Fig. 11, Region B has brighter absolute magnitudes and lower masses than the nearby unresolved binary/triple track. Unresolved 4-body systems may explain the bright absolute magnitudes of Region B but are inconsistent with its lower inferred masses. Region B is most likely associated with pre-MS stars (Zari et al. 2018; McBride et al. 2021), when the young stars are still contracting to the main sequence. The wide binary catalogue we use for mass inference has explicitly excluded clustered regions like star-forming regions (El-Badry, Rix & Heintz 2021), and our selection excludes systems at low Galactic latitudes (|b| < 10 deg). Interestingly, there are still sufficient remaining pre-MS wide binaries that dominate the sample in Region B, leaving a significant signal in our mass inference. In Fig. 11, the long solid red line shows an example of a double pre-MS wide binary (primary’s Gaia DR3 ID: 49366530195371392), where both of the components have an infrared excess of W1−W2 = 0.18 and 0.25 mag (Wright et al. 2010; Marrese et al. 2019; Marocco et al. 2021). This wide binary is considered a candidate member of the Taurus-Auriga star-forming complex (Kraus et al. 2017). The short solid red line is another double pre-MS wide binary (primary’s Gaia DR3 ID: 6016457297110757760), where both of the components have an infrared excess of W1−W2 = 0.39 and 0.19 mag. This wide binary is a candidate member of the Upper Centaurus-Lupus region (Kuruwita et al. 2018), and its TESS light curve presents the ‘dipper’ behaviours (Capistrant et al. 2022), which is a class of variability in young stellar objects that may be associated with circumstellar obscuration events, spots, or disc activities (Cody et al. 2014). These two examples have relatively low surrounding stellar densities, allowing them into our wide binary catalogue. These young wide binaries are particularly important for testing the formation theory of wide binaries (Kouwenhoven et al. 2010; Moeckel & Clarke 2011; Tokovinin 2017; Rozner & Perets 2023; Xu et al. 2023) and for understanding the multiplicity of young stars (Kraus & Hillenbrand 2009; Kraus et al. 2012; Kounkel et al. 2019).

Another application of our results is the mean mass of giants. At BP−RP = 1.2 and M_G = 1.0 (Fig. 12), our measured mean mass of giants is 1.31 ± 0.15 M_⊙, which coincides with the peak of the giant mass distribution measured from asteroseismology (Yu et al. 2018; Wu et al. 2023). The inferred mean mass of white dwarf is 0.67 ± 0.07 M_⊙, in good agreement with the literature measurements from gravitational redshifts (Falcon et al. 2010, 2012; Chandra et al. 2020). Unfortunately, we do not detect the mass gradient among white dwarfs. We discuss possible future improvements in the next section.

We summarize all the features of mass measurements in Fig. 12. With the wide binaries from the Gaia survey, we measure the dynamical masses homogeneously across the H–R diagram, covering main-sequence stars, giants, and white dwarfs.

4.3 Caveats and future prospects

One caveat of our work is that we do not recover the mass gradient in the white dwarfs detected by gravitational redshifts (Chandra et al. 2020). As discussed in Section 3.2, the non-detection may be due to the smaller sample sizes of high-mass white dwarfs among the wide binary sample. We have attempted to train a separate white dwarf model with no success. Future Gaia releases may be critical, improving the sample size of high-mass white dwarfs and the astrometric precision.

In recent years, there has been an exciting discovery of stellar mass black holes in binary orbits (Thompson et al. 2019; Chakrabarti et al. 2023; El-Badry et al. 2023a, b) and as a single system from microlensing (Sahu et al. 2022, but see Lam et al. 2022). One of them is a 9.62 ± 0.18 M_⊙ black hole surrounding a 0.93 ± 0.05 M_⊙ Sun-like star (El-Badry et al. 2023a), with BP−RP = 1.16, M_G = 5.35, and metallicity [Fe/H]=−0.2 ± 0.05. No comoving wide companion is found around the system. In principle, our mass measurements can constrain the occurrence rate of such ∼10 M_⊙ stellar black holes among wide binaries. For example, at BP−RP = 1 on the single-MS track in Fig. 8, our measured mass is 0.97 ± 0.08 M_⊙, while the PARSEC isochrone gives a mass of 0.89 M_⊙, assuming a solar metallicity. If we take the isochrone mass as the ground truth of a single star’s mass, we obtain an upper limit of |$\sim 1~{{\ \rm per\ cent}}$| of the occurrence rate of 10- M_⊙ black holes among these wide binaries. We are cautious that this rough estimate assumes that the black-hole companions are not associated with the high-|$\tilde{u}$| tail in the model, and future work is needed to investigate the nature of the high-|$\tilde{u}$| outliers. Nevertheless, the prospect of constraining compact objects’ occurrence rates is exciting, especially since this approach may be the only possible method to constrain such populations out to ∼10² AU separations from the host stars (as long as they remain a hierarchical, stable three-body system).

Here we discuss some other caveats and prospects for this work. We adopt a thermal eccentricity distribution in this work to perform the inference. However, it is known that the eccentricity distribution is a function of binary separations (Tokovinin 2020; Hwang, Ting & Zakamska 2022b), and also that some wide-binary populations like wide twins have distinct eccentricity distributions (Hwang et al. 2022a). Although Fig. 7 suggests that the final result is consistent with the thermal eccentricity distribution as the entire sample, it is more reliable to treat eccentricity distributions as free parameters in future work. Further investigation of the nature of the high-|$\tilde{u}$| tail and the underlying mass distribution are critical to robustly constraining compact object occurrence rates and to testing the Newtonian gravity at larger binary separations. In this work, we only consider two parameters (BP−RP and M_G), and it would be valuable to include other astrophysical parameters. In particular, metallicities are available from large spectroscopic surveys and Gaia itself (e.g. XP spectra). Mapping the metallicity-dependent mass measurements for the H–R diagram would provide unique insights into our understanding of stellar physics.

5 CONCLUSIONS

In this paper, we use wide binary stars from the Gaia survey to measure dynamical masses across the H–R diagram. Specifically, wide binaries’ orbital velocities and separations measured by Gaia provide critical information for their masses (Section 2.1 and 2.2), and we model the relation between mass and the H–R diagram using statistical inference and a neural network (Section 2.3). Our neural network-based method enables us to model the mass along two dimensions (colour and absolute magnitude), which is critical to cover the entire H–R diagram. Our mass measurements are purely based on Keplerian law, thus providing independent constraints for stellar evolution models. Using a mock wide binary sample, we show that our method can recover the underlying masses in the H–R diagram (Figs 3 and 4).

The resulting mass measurements provide a global view of mass structures across the H–R diagram (Figs 5 and 12). Our findings are summarized as follows:

• The inferred mass of the single main-sequence stars spans from 0.1 to 2 M_⊙ (Fig. 8). Our mass measurements overall agree with the PARSEC and MIST isochrone models at masses >0.6 M_⊙. At lower masses <0.6 M_⊙, our measurements are more consistent with the PARSEC model.

• Our results show an unresolved binary/triple track parallel to the main-sequence single stars (Fig. 9). An additional mass excess at BP−RP between 1.5 and 2 suggests a population of unresolved main-sequence triples (Fig. 10).

• Two regions in the H–R diagram show interesting mass behaviours (Fig. 11). Region A is located at the faint end of the main sequence at BP−RP between 2 and 3, with an unexpectedly higher mass of ∼0.8 M_⊙. The higher mass may be due to unresolved M dwarf-white dwarf binaries. Region B is located brighter than the nearby unresolved binary/triple track, with a smaller mass of ∼0.6 M_⊙ than the unresolved binaries/triples. Region B is mostly likely the pre-MS stars. We have identified two example wide binaries where both component stars are in Region B, have infrared excess, and are candidate members of nearby star-forming regions.

• The measured mean mass of giants and white dwarfs are 1.31 ± 0.15 M_⊙ and 0.67 ± 0.07 M_⊙, respectively. These measurements are consistent with the literature values. The non-detection of the mass gradient among white dwarfs may be due to the lack of high-mass white dwarfs in the wide binary sample.

• We discuss some future improvements, including modelling the eccentricity distribution during the mass inference and a better understanding of the high-|$\tilde{u}$| outliers. These improvements are critical to constraining the occurrence rate of compact object companions and to testing the non-Newtonian gravity theory using wide binaries.

ACKNOWLEDGEMENTS

The authors thank the referee, Hans-Walter Rix, for encouraging and constructive comments that significantly improved the manuscript. The authors appreciate discussions with Josh Winn, Jo Bovy, Chris Hamilton, and Nadia Zakamska. YST acknowledges financial support from the Australian Research Council through DECRA Fellowship DE220101520. HCH acknowledges the financial support of the Friends of the Institute for Advanced Study Member Fund at the Institute for Advanced Study. HCH thanks the constant support and encouragement for this work from Ting-Wei Young, Nae-Chyun Chen, and Ting-Wei Liao.

Facilities: Gaia.

Software:IPython (Pérez & Granger 2007), jupyter (Kluyver et al. 2016), Astropy (Astropy Collaboration 2013, 2018), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), PyTorch (Paszke et al. 2019).

DATA AVAILABILITY

The data underlying this article are available online. The data sets were derived from sources in the public domain: Gaia Data Archive https://gea.esac.esa.int/archive/

Footnotes

References

APPENDIX A: MARKOV-CHAIN MONTE-CARLO WITH A DISCRETE-MASS MODEL

Here we demonstrate how to use MCMC to measure the dynamical masses from wide binaries and discuss its limitations. We consider a model with five mass parameters m[k] and assign every star an integer mass index k ∈ [0, 4] based on their locations in the H–R diagram. The colour code in Fig. A1 represents the assigned mass indices, where m[4] is for all white dwarfs and m[0–3] are for main-sequence stars with different ranges of BP−RP colours. Then the total mass of wide binary i is m_{tot, i} = m[k_{i, 1}] + m[k_{i, 2}], where k_{i, 1} and k_{i, 2} are the mass indices of the component star 1 and 2. The same Gaia sample selection is used and described in Section 2.4. Among them, we randomly select 10 000 wide binary pairs to reduce the computation time.

We use the MCMC package emcee (Foreman-Mackey et al. 2013) to sample the posteriors of m[k]. The likelihood is equation (6) and equation (10), and we adopt flat priors for the mass parameters. The parameters in the outlier model is F = 0.20, u_outlier = 30, σ_outlier = 20. We marginalize the uncertainties of u using equation (7), where we approximate the integration by a discrete summation from u′ = 0 to u′ = 100 with a step of Δu′ = 1. Binaries with uncertainties σ_u < Δu′ have a narrow p(u_i|u′) that would be undersampled, so we manually set their σ_u to Δu′ = 1, which does not affect the main result because |$p(\tilde{u})$| is a smooth distribution on the scale of Δu′.

Fig. A2 shows the posterior distributions (made using corner; Foreman-Mackey 2016) of the mass parameters. This result consists of 32 walkers, each with 5000 steps. The computation time takes about 1.5 h on the Google Colab platform (no GPU is used). The 16-50-84 percentiles of masses are:

These measurements are consistent with the literature. For example, the mean white dwarf mass from MCMC is |$m[4] =0.631_{-0.059}^{+0.063}$| M_⊙, consistent with the mean mass of |$0.647^{+0.013}_{-0.014}$| M_⊙ DA white dwarfs using gravitational redshifts (Falcon et al. 2010). The left panel of Fig. A3 is the resulting distribution of |$p(\tilde{u})$|⁠, where |$\tilde{u}$| is computed using the 50- percentile masses, showing that the result is in good agreement with the model. Also, the high-|$\tilde{u}$| feature is similar to the one from the neural-network method (Fig. 7).

The right panel of Fig. A3 compares the MCMC results with the neural network results from Fig. 9. The horizontal values of the MCMC results (black crosses) are the median BP−RP in the sample, and their horizontal bars indicate the range of the BP−RP bins. The vertical bars are the 1-sigma mass uncertainties. Compared to the neural network’s single-star results (blue line), the MCMC masses are slightly overestimated because the MCMC sample includes both single stars and unresolved binaries. We also overplot the result from Giovinazzi & Blake (2022) as the red line, where the authors use twin wide binaries to derive the relation between mass and the absolute RP magnitudes. We apply their fit (their equation 8) to derive the masses for the field MS stars, and then the red line in Fig. 9 is the median mass binned by their BP−RP colours. The blue and red lines overall agree with some noticeable differences. We are cautious that this comparison is mainly for illustration, and the detailed difference is subject to the photometry conversion, the choice of the fitting function used in Giovinazzi & Blake (2022), and the different wide binary samples.

The advantage of MCMC is that it provides robust uncertainties for the masses. Using only 10 000 (8 per cent) out of 126 000 wide binaries, we can reach mass uncertainties down to 0.01–0.07 M_⊙. These uncertainties give us a sense of mass uncertainties in the neural-network method, which considers a larger sample but smaller scale structures (smaller bins in the MCMC sense, although the neural network does not explicitly perform binning). The disadvantage of MCMC is that it is time-consuming. With only 8 per cent of the sample and only 5 parameters for the model, MCMC already takes 1.5 h to converge. If we want to explore structures in the mass measurements of the entire H–R diagram, MCMC requires much more (>50) mass indices and its convergence becomes challenging within a reasonable computation time. This is the motivation for us to develop a neural network-based method in the main text.