Dynamical masses across the Hertzsprung-Russell diagram

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 two-dimensional 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$_\odot$, 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 behaviors 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, providing independent constraints on stellar evolutionary models and the occurrence rate of compact objects.


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 et al. 2010).Therefore, to validate and improve our understanding of stellar physics and models, it is critical to have modelindependent 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 2 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 et al. 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 eclips-ing probability is proportional to the radii at a fixed semi-major 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. 2023).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(Falcon et al. , 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 doublelined 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 et al. 2016) enables the identification of millions of resolved wide binaries (El-Badry & Rix 2018;Tian et al. 2020;Hartman & Lépine 2020;El-Badry et al. 2021), critical to covering the entire H-R diagram.These resolved wide binaries have orbital periods longer than 10 3 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. 2022b,a), 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 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.Here we use both equal-mass and non-equal-mass binaries to cover the two-dimensional H-R diagram along the BP−RP colors 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 Sec. 2. In Sec. 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 Sec. 4 and conclude in Sec. 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 2 AU), regardless of whether they 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.

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 semi-major axis of the circular binary, G is the gravitational constant, and m tot = m 1 + m 2 is the total mass of the binary where the component  masses are m 1 and m 2 .The subscript '3D' indicates that they are the magnitudes of full 3-dimensional 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 tot = a 3D v 2 3D /G.In reality, full-3D magnitudes of v 3D and a 3D in Eq. 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 semi-major axis in the case of a circular orbit.Similar to Eq. 1, these projected quantities and phase-dependent separations are related to the system mass.Hence, we define u as: where v is the two-dimensional orbital velocity, s is the two-dimensional 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 −1 , and u in units of km s −1 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 disk 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 (⃗ v) and the separation vector (⃗ s).Using this method, Hwang et al. (2022b) showed that the eccentricity distribution f e (e) is close to uniform (f e (e) = 1) for wide binaries at ∼ 10 2 AU, approximately thermal (f e (e) = 2e) at 10 3 AU, and then superthermal (f e (e) ∝ e 1.3 ) at > 10 3 AU.Since our sample is dominated by wide binaries at ∼ 10 3 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).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 modeling 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 ũ are in units of km s −1 AU 1/2 M ⊙ −1/2 .Then we derive the probability distribution functions p(ũ|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 1 .p(ũ|e) is numerically normalized.Fig. 1 shows the resulting p(ũ|e), where the color represents the value of e.For a circular orbit (purple line), its p(ũ|e = 0) is nonzero up to its constant circular orbital velocity ũ = 29.8.Its strong peak at ũ = 18.4 and the tail toward ũ = 0 is caused by the projection effect of s and v.
We numerically integrate p(ũ|e) to derive the case for the thermal eccentricity distribution, p(ũ|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(ũ|f e = 2e) as where A = 4.95 × 10 −3 , B = 2.24 × 10 −3 , C = 3.85, and ũ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/ √ m tot ensures proper normalization.To simplify the notation, we use p(u|m tot ) := p(u|m tot , f e = 2e) and p(ũ) := p(ũ|f e = 2e) for the case of the thermal eccentricity distribution.One important property of u and ũ is that p(u|m tot ) and p(ũ) are independent of the semi-major axis a.For an arbitrary binary, its u = v √ s would follow p(u|m tot ) regardless of its a because the Keplerian motion is scalefree and the projection effect does not depend on a.Therefore, semi-major axes do not enter our likelihood later and there is no need to marginalize over the underlying semi-major axis distribution.

From total masses to individual masses
In the last section, we derive the likelihood p(u|m tot ) (Eq. 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, colors and absolute magnitudes), which provides the connection between the total mass (m tot ) and the component masses (m 1 and m 2 ).
For example, we can select a sample of equal-mass (twin) wide binaries where all individual components have identical colors 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% 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 A, B, and C), and there are three pairs of wide binaries, each consisting of star A+B, B+C, and A+C.Once we obtain the total masses of each wide binary from their orbital velocities and separations, m tot,1 = m A + m B , m tot,2 = m B + m C , and m tot,3 = m A + m C , the three individual masses can be solved by m A = (m tot,1 + m tot,3 − m 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 solve 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 , 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 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 color and the absolute G-band magnitude measured by Gaia as the observables in this paper.Therefore, m i,j = M((BP − RP) i,j , 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 color 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 ) can be evaluated using Eq. 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 σ ui .Here we assume a flat prior for p(u ′ ) so we can approximate p(u i |u ′ ) using p(u ′ |u i ).We discuss σ ui from the data perspective in Sec.2.4.
Our formulation of 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 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.

Neural network model for M
There are several approaches to formulate the function M(BP−RP,G).For example, one can bin the H-R diagram and assign a mass to each bin, i.e., m[k] = M(BP − RP, 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 M(BP − RP, 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% 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 color (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 2 (mean squared errors) loss function, we use the likelihood (Eq.6 with the outlier mixture model introduced in the next section) as the loss function and then adjust the model parameters using 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 minutes using a V100 GPU.

Wide binaries from Gaia
We use the wide binary catalog from El-Badry et al. (2021), where the resolved wide binaries were searched using the Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2016Collaboration et al. , 2018Collaboration et al. , 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 catalog.For each wide binary, El-Badry et al. ( 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 3 AU where the contamination from chancealignment pairs is negligible.
For each Gaia wide binary, we measure its projected physical separation s i (in AU) and projected orbital ve- 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 et al. ( 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 ∆µ The uncertainty of proper motion difference is ) 2 , and ∆µ 2 δ = (µ δ,2 − µ δ,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% precision, and parallaxes better than a few percent 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 −1 , corresponding to a proper motion difference of 0.3 mas yr −1 at parallaxes of 2.5 mas.For typical proper motion uncertainties of 0.1 mas yr −1 in Gaia EDR3, this sample selection ensures that most of the wide binaries can have signal-to-noise ratios (SNR) of proper motion difference ∆µ/σ ∆µ > 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 et al. 2022a,b).Indeed, the majority of the resolved tertiaries of inner astrometric binaries (identified by Gaia Data Release 3; Halbwachs et al. 2022) 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. 2022), 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 et al. 2022a); therefore our sample may still be subject to possible outliers due to unresolved companions.
Since p(ũ) has a sharp cutoff at ũ ∼ 40 (Fig. 1), the likelihood model from Eq. 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 et al. 2010): where p i (u i |m tot,i ) is Eq. 5 for well-behaved data, and F is the fraction of data that is better described by the outlier model p outlier (ũ).We adopt a Gaussian (truncated at 0) for the outlier model: Figure 2. Left: the number of wide binaries in our sample, including both primary and secondary stars.Right: the median u (in units of km s −1 AU 1/2 ) across the H-R diagram, where u is positively correlated with the mass in each bin.Several features of u (and therefore mass) are directly noticeable, including the higher masses in blue main-sequence stars, the lower masses in red main-sequence stars, and the higher-mass unresolved binaries (and triples) parallel to the main sequence.
where Z is the normalization factor (not exactly 2πσ 2 outlier due to the truncation at 0).After trialand-error with the Gaia data, we use ũ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 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 et al. 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 disk, 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 3 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

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% are genuine wide binaries without any unresolved companions, 40% are triples consisting of one single star and one unresolved binary, and 10% are quadruples consisting of two unresolved binaries.20% of the genuine wide binaries (i.e., 10% of all pairs) are equalmass 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 brutus2 (Speagle et al. in prep) 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 ũi from the distribution p(ũ) assuming that they follow the thermal eccentricity distribution (Eq.4), and scale it to u i = ũi √ m 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 colors, regardless of whether they have unresolved companions or not.Specifically, the mass is 2.5−0.5kM ⊙ 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 mainsequence mass prohibits some parameter space of unresolved non-equal-mass binaries.Then we apply our neural-network-based method in Sec. 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 ũ distribution using the mass measurements from the training (blue), in agreement with the likelihood model p(ũ) from Eq. 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.Astro-physically, at a fixed metallicity, the masses of stars determine the surface temperatures and therefore BP−RP colors, and the dependence of mass on absolute magnitudes at a fixed BP−RP is due to the mass ratios of unresolved binaries.The 2nd and the 3rd 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 4th and 5th panels demonstrate that the uncertainty estimates also report reasonably higher uncertainties at the high-mass end.The last panel shows that the resulting ũ using the inferred masses (blue) agrees well with the likelihood model from Eq. 4 (black line).
The tests in Fig. 3 and 4 suggest that our neuralnetwork 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.

Masses across the H-R diagram using Gaia wide binaries
Among the 99,680 pairs of Gaia wide binaries from Sec. 2.4, we randomly select 90% of them as the training sample and 10% as the validation sample.We then train our models using a learning rate of 1 × 10 −4 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 afterward.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 M in Sec.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% 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 ⊙ .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 colors.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.This discrepancy may be due to the fact that the number of white dwarfs is insufficient in the sample, only ∼ 5% of the component stars.Furthermore, the mass distribution of white dwarfs strongly peaks around 0.6 M ⊙ (Tremblay et al. 2016), and the highermass 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(ũ) of from the training result (blue).Different from the theoretical p(ũ) (orange line), the observed p(ũ) has a main component and an extended high-ũ tail up to ũ ∼ 80.The main component agrees well with the p(ũ) from the thermal eccentricity distribution (orange), supporting our observationally motivated assumption that the eccentricity  The orange line shows the theoretical p(ũ) for the thermal eccentricity distribution (Eq.4), the green line is the outlier model (Eq.10), and the black line is the adopted mixture model (Eq.9).The resulting p(ũ) from the training (blue) agrees well with the mixture model (black).
distribution follows the thermal eccentricity distribution (Hwang et al. 2022b).Overall, the behavior of the highũ tail is well described by our assumed outlier model p outlier (ũ) with an outlier fraction of F = 0.20 (Sec.2.4).The high-ũ tail up to ũ ∼ 80 cannot be explained by eccentricities because p(ũ|e) is strongly truncated above ũ = 40 for all e (Fig. 1) due to the small fraction of time spent at pericenter.The high-ũ 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-ũ tail is reduced to F ∼ 10%, suggesting that a significant fraction of the high-ũ tail is associated with fainter sources.
The high-ũ 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 et al. 2022a).Even though the nature of the high-ũ 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 4 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;Hernandez 2023;Chae 2023, but see alternative explanations in El-Badry 2019).In Fig. 7, we show that at ∼ 103 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 4 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 behaviors 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    8).The orange line is the blue line offset brighter by ∆G = 0.75 mag, the expected magnitude offset for unresolved equal-mass binaries.We investigate the mass along the green vertical line in Fig. 10.Right: the mass as a function of BP−RP.The dark blue line shows the mass along the single-star main-sequence track, and the light blue line shows the double mass of the dark blue line.The inferred masses along the ∆G = 0.75 track (orange) show that the inferred mass is consistent with equal-mass binaries at BP−RP> 2.5, and that there is a mass excess at BP−RP∼ 2 that may be due to unresolved triples.

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 3 -10 4 AU) binary fraction peaks around the solar metallicity (Hwang et al. 2021(Hwang et al. , 2022c)).Therefore, the mean metallicity of wide binaries is similar to that of the field stars in the solar neighborhood, 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 days are more common among wide pairs (i.e., wide tertiaries) (Hwang et al. 2020;Fezenko et al. 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%, 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 Sec.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-neighborhood 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 et al. 2019).In the regions of sub-giants and giants in the H-R diagram, the isochrones with different ages are overlapped, and our measurements represent their masses in an average sense.

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 (Tremblay et al. 2019;Cheng et al. 2019).In this work, we provide a different viewthe 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, 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;Tang et al. 2014;Chen et al. 2014Chen et al. , 2015) ) and from MIST (version 1.2 with rotation v initial /v critical = 0.4, Paxton et al. 2011Paxton et al. , 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 highmass (> 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 lowmass end (BP−RP> 2, or stellar masses < 0.5 M ⊙ ), our measurements agree better with the PARSEC isochrone than MIST, providing further constraints on the stellar models of low-mass stars.
Fig. 9 presents our results at the unresolved binary (and triple) track.An equal-mass binary would have the same BP−RP colors 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 PAR-SEC 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 equalmass 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 equalmass 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 lowermetallicity stars (Raghavan et al. 2010;Yuan et al. 2015;Badenes et al. 2018;Moe et al. 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.
Fig. 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(Borkovits et al. , 2020a,b;,b;Czavalinga et al. 2023;Rappaport et al. 2023;Kostov 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 behaviors.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 lowermetallicity 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(Hwang et al. , 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 ũ, and Region A has a larger fraction of wide binaries associated with the highũ tail (Sec.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-ũ tail.The fainter absolute magnitudes of Region A may also be the reason why the high-ũ tail is more related to the fainter sources, as discussed in Sec.3.2.Future investigations are needed to confirm the nature of Region A and its connection with the high-ũ 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 catalog we use for mass inference has explicitly excluded clustered regions like star-forming regions (El-Badry et al. 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' behaviors (Capistrant et al. 2022), which is a class of variability in young stellar objects that may be associated with circumstellar obscuration events, spots, or disk activities (Cody et al. 2014).These two examples have relatively low surrounding stellar densities, allowing them into our wide binary catalog.These young wide binaries are particularly important for testing the formation theory of wide binaries (Kouwenhoven et al. 2010;Moeckel & Clarke 2011;Tokovinin 2017;Xu et al. 2023;Rozner & Perets 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(Falcon et al. , 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.

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 Sec.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;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 ∼ 1% 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-ũ tail in the model, and future work is needed to investigate the nature of the high-ũ 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 2 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 the current work to perform the inference.However, it is known that the eccentricity distribution is a function of binary separations (Tokovinin 2020;Hwang et al. 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-ũ 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 like metallicities (from other spectroscopic surveys or Gaia XP spectra) to map the metallicity-dependent relation.

CONCLUSIONS
In this paper, we use wide binary stars from the Gaia survey to measure dynamical masses across the Hertzsprung-Russell diagram.Specifically, wide binaries' orbital velocities and separations measured by Gaia provide critical information for their masses (Sec.2.1 and 2.2), and we model the relation between mass and the H-R diagram using statistical inference and a neural network (Sec.2.3).Our neural network-based method enables us to model the mass along two dimensions (color 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 (Fig. 3, 4).
The resulting mass measurements provide a global view of mass structures across the H-R diagram (Fig. 5, 12).Our findings are summarized as follows: 1.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.
2. 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 mainsequence triples (Fig. 10).
3. Two regions in the H-R diagram show interesting mass behaviors (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.
4. 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.5. We discuss some future improvements, including modeling the eccentricity distribution during the mass inference and a better understanding of the high-ũ 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.

Figure 1 .
Figure 1.The simulated distributions for p(ũ|e).The colored solid lines show the p(ũ|e) for single-valued eccentricity e, where the value is represented by the color.The solid black line shows the distribution for the thermal eccentricity distribution, i.e., p(ũ|fe(e) = 2e).The dotted red line is the functional fit using Eq. 4.
shows the me-dian u in each bin of the left panel.Since u ∝ √ m tot (Eq.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 mainsequence 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.

Figure 3 .
Figure 3. Mock wide binary sample with a toy model of discrete masses.First row: true mass of the sample (left), inferred mass (middle), and the difference between the inferred mass and the true mass (right).Second row: uncertainty of the inferred mass (left), the inferred fractional error (middle), and the resulting p(ũ) distribution.

Figure 4 .
Figure 4. Mock wide binary sample with astrophysical masses.Each mock wide binary consists of two component stars, and the mass here is the mass of the component star, which can be a single star or an unresolved binary.Plot descriptions are the same as Fig. 3.These tests show that our inference can recover the underlying mass from the orbital motions of wide binaries.

Figure 5 .
Figure 5.The inferred dynamical masses across the Hertzsprung-Russell diagram.Each point is one component star in the Gaia wide binary sample.Here we use wide binaries to represent the mass model, but we emphasize that the mass model itself is a continuous function and no data binning is involved.With our neural-network approach, we measure the masses of main-sequence stars, unresolved binaries/triples of the main sequence, (sub)-giants, and white dwarfs.An annotated version is shown in Fig. 12.

Figure 6 .
Figure 6.The estimated fractional errors of the mass measurements.The mass of most regions in the H-R diagram can be measured at ∼ 10% precision.

Figure 7 .
Figure7.The resulting p(ũ) of the wide binaries from the training (blue histogram).The orange line shows the theoretical p(ũ) for the thermal eccentricity distribution (Eq.4), the green line is the outlier model (Eq.10), and the black line is the adopted mixture model (Eq.9).The resulting p(ũ) from the training (blue) agrees well with the mixture model (black).
3 .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.

Figure 8 .
Figure 8.Comparison of masses in the single-star track.Left: the H-R diagram of the field stars.The blue line shows the main-sequence fit to the field stars.The orange and the green markers are the PARSEC and MIST isochrones respectively, with a stellar age of 0.1 Gyr.Right: the mass as a function of BP−RP.The blue band is the measured masses along the main-sequence fit (blue line in the left panel).The vertical width of the line represents the 1-σ uncertainties of the mass.The orange and green markers show the masses from the isochrone models.Our measurements agree with both isochrones at BP−RP< 1.5.At BP−RP> 1.5, the measured masses are more consistent with the PARSEC model.

Figure 9 .
Figure 9.Comparison of masses in the unresolved binary/triple track.Left: the H-R diagram of the field stars.The blue line is the fit to the main-sequence stars (same as Fig.8).The orange line is the blue line offset brighter by ∆G = 0.75 mag, the expected magnitude offset for unresolved equal-mass binaries.We investigate the mass along the green vertical line in Fig.10.Right: the mass as a function of BP−RP.The dark blue line shows the mass along the single-star main-sequence track, and the light blue line shows the double mass of the dark blue line.The inferred masses along the ∆G = 0.75 track (orange) show that the inferred mass is consistent with equal-mass binaries at BP−RP> 2.5, and that there is a mass excess at BP−RP∼ 2 that may be due to unresolved triples.

Figure 10 .
Figure 10.Mass measurements (green line) with respect to the absolute G-band magnitude (brighter to the right direction) at a fixed BP−RP= 1.6 (green vertical line in Fig. 9, left panel).The single-star symbol represents the absolute magnitude and the mass of a single star.The doublestar and triple-star symbols are for unresolved equal-mass binaries and unresolved equal-mass triples.The filled star symbols consider PARSEC models with [Fe/H]= −0.35, and the open star symbols for [Fe/H]= 0. The mass peaking at 1.60 ± 0.15 M⊙ along with its brighter absolute magnitude suggests a population of unresolved triples at this location of the H-R diagram.

Figure 11 .
Figure 11.Two regions(A and B)  showing interesting behaviors in the inferred masses.The straight lines are example wide binaries associated with these regions, where the endpoints of the line are the component stars of the wide binary.The higher mass in Region A may be associated with unresolved M dwarf-white dwarf binaries.Region B is most likely due to the pre-main-sequence wide binaries.

Figure 12 .
Figure 12.Same results as Fig. 5, the measured dynamical masses across the Hertzsprung-Russell diagram.Here we highlight several interesting masses in the diagram.

Figure A1 .
Figure A1.The H-R diagram where the color code represents their mass indices k used in the MCMC method.

Figure A2 .
Figure A2.The posterior distribution of the mass parameters in the discrete mass model.

Figure A3 .
Figure A3.Left: the distribution of p(ũ) from the MCMC method.Right: comparison of mass measurements from the neural-network method (blue and orange line, same as in Fig.9), the MCMC method (black crosses), and the functional fit from Giovinazzi & Blake (2022) (red line).

Table 1 .
Table of notations.