© The Authors 2025

1 Introduction

Binary star systems constitute fundamental components of the stellar population, with half of solar-type stars residing in multiple systems (Raghavan et al. 2010; Duchene & Kraus 2013; Moe et al. 2019). Binary systems provide an important context and set of constraints for most sub-fields of astrophysics (Rix et al. 2019). The statistical properties of binary populations – including the binary fraction, period distribution, mass ratio distribution, and eccentricity distribution – provide constraints on star formation theories and dynamical evolution models. Moreover, these characteristics vary systematically with stellar properties such as stellar mass, age, and metallicity, creating distinctive patterns across different stellar populations and offering valuable insights into the physical mechanisms governing both the formation and subsequent evolution of stellar systems (Duchêne & Kraus 2013; Moe et al. 2019; El-Badry et al. 2021).

The periods of binary systems span a wide range in timescales (from minutes to Myr). Current stellar surveys such as Gaia, LAMOST, SDSS-IV/APOGEE, and SDSS-V (Gaia Collaboration 2016; Liu et al. 2020; Majewski et al. 2017; Kollmeier et al. 2017) primarily identify binaries through spectroscopic (El-Badry et al. 2018; Price-Whelan et al. 2020; Zhang et al. 2022) or astrometric signatures of orbital motion (e.g., Holl et al. 2023; Penoyre et al. 2022a,b). Previous studies of unresolved binaries have generally relied on spectroscopic methods that are most effective for systems exhibiting radial velocity (RV) variations, such as identifying double-lined spectroscopic binaries (SB2s) through cross-correlation techniques (e.g., Li et al. 2021a) or disentangling component spectra from spectra taken at different orbital phases (e.g., Seeburger et al. 2024). These methods are sensitive only to systems with notable velocity separations (v_los ≳ 10−30 km s⁻¹) and high mass-ratios (q ≥ 0.7) (Zhang et al. 2022). This approach inherently limits detection to relatively short-period binaries, which represent only about one-third of the total binary population (Duchêne & Kraus 2013).

Similarly, the identification of single-lined spectroscopic binaries (SB1s) through radial velocity variability is also biased toward shorter period systems. Traditional methods are effective within specific period ranges but miss many systems outside their sensitivity windows. Spectroscopic techniques, such as identifying single-lined (SB1) and double-lined (SB2) spectroscopic binaries, are limited to systems with periods of approximately 1–1000 days, where RV variations are detectable (El-Badry 2024). Astrometric methods, on the other hand, are most sensitive to binaries with periods comparable to Gaia’s observational baseline (Castro-Ginard et al. 2024). In contrast, spectral energy distribution (SED)-based techniques remain effective across the full range of periods for unresolved binaries with distinct components, as they do not rely on RV or astrometric variability.

This methodology bridges detection gaps, complementing the short-period sensitivity of eclipsing binary searches (requiring specific orbital alignments) and long-period capabilities of direct imaging (typically exceeding 10 000 days for Gaia). SED-based approaches can detect the estimated 73% of binaries with periods greater than years (El-Badry et al. 2018; Duchêne & Kraus 2013) that remain invisible to conventional RV surveys (Gao et al. 2014; Badenes et al. 2018), while simultaneously providing independent verification of systems within the detection windows of other methods. For binaries with periods too long for dynamical detection yet too tight for spatial resolution by Gaia, our XP-based methodology may represent a key identification technique,

The SED method identifies binary stars by modeling the composite SED of two stellar components and detecting deviations from single-star models, with enhanced sensitivity when the components exhibit obvious temperature or luminosity contrasts (e.g., Aznar Cuadrado & Jeffery 2001). The SED-based methods complement these approaches by operating independently of orbital dynamics, identifying binaries through their composite spectral signatures regardless of period. These signals are most pronounced when the two components have different spectral features, as is the case with white dwarf–main-sequence binaries (Li et al. 2025).

However, the presence of an unresolved companion changes the observable spectrum in all cases where the two stars are not identical. When two main-sequence stars have similar masses, the normalized spectral features or normalized SED show minimal differences between single and binary systems, creating a potential degeneracy (El-Badry et al. 2018). However, with precise distance measurements (e.g., from Gaia parallaxes), this degeneracy can be resolved: a binary system’s total luminosity is the sum of its two components, so for stars with similar spectral features, a binary would appear twice as luminous as a single star. When forced to fit such a binary with a single-star model, the excess luminosity is misinterpreted as a hotter (bluer) temperature to compensate for the unphysically high flux at the known distance.

Together, these complementary detection techniques enable stellar surveys to deliver samples of binary stars and stellar companions orders of magnitude larger than previously known and spanning all stages of stellar evolution. Even without detectable velocity offsets, unresolved binary spectra contain distinctive signatures of unseen companions. By combining the quality of Gaia DR3 low-resolution BP/RP (XP) spectra (Gaia Collaboration 2023b) with Gaia’s precise parallax measurements (Gaia Collaboration 2016), we can independently identify binary systems through chi-squared spectral fitting (Li et al. 2025).

The characterization of binary population properties will emerge from integrated analysis across multiple stellar surveys, which collectively sample diverse stellar types, wavelengths, and observational techniques. This approach addresses gaps in binary detection, complementing both eclipsing binary searches at short periods and direct imaging at long periods, while capturing the estimated 70% of binaries with periods exceeding 10 years that traditional RV surveys miss (Duchêne & Kraus 2013; El-Badry et al. 2018) for binaries with periods too long for dynamical detection yet too close for spatial resolution.

In this paper, we focus on main-sequence main-sequence (MSMS) binary stars from Gaia XP spectra, presenting an inventory of unresolved binary systems in which both components are main-sequence stars. We employed the high quality of Gaia DR3 low-resolution XP spectra (Gaia Collaboration 2023b), combined with precise parallax measurements from Gaia, (Gaia Collaboration 2016), and we conduct a homogeneous search for binary systems throughout the Galaxy, identifying and characterizing these systems with completeness and consistency.

Our approach is based on the methodology introduced by Li et al. (2025). We utilized Gaia XP spectra, which provide low-resolution spectrophotometry for over 220 million sources in the Galaxy. Although these spectra have a lower resolution (R ~ 20−100) than traditional spectroscopic surveys, they offer sufficient signal-to-noise and self-consistent flux calibration (De Angeli et al. 2023). By modeling the composite SED of binary components, we can identify binaries regardless of their velocity separation, including long-period systems with negligible velocity offsets between components, provided they have favorable mass ratios (0.5 ≲ q ≲ 1.0).

In this work, we present an MSMS catalog derived from Gaia XP spectra, containing over 14 million binary candidates. The remainder of the paper is structured as follows: In Section 2, we describe the Gaia DR3 data and our sample selection criteria. Our methodology, described in Section 3, builds on our forward model for stellar spectra mapping of intrinsic parameters: stellar mass (M), photometric metallicity ([M/H]_phot), and extinction (E) to predicted spectra. In Section 4, we present the validation results and the method to classify the single and binary stars using parameter-dependent thresholds. Section 5 details the statistical properties of our sample, including our binary candidate selection process and the creation of a high-confidence “golden sample” of approximately one million binary systems. We discuss the implications of our findings and summarize our conclusions in Section 6. This catalog represents a valuable resource for investigating the formation and evolution of binary systems throughout the Galaxy. The size and diversity of our sample enables robust statistical analyses of binary system properties across different Galactic environments, offering new insights into the processes governing binary star formation and evolution.

2 Data

In this section, we first describe the properties and processing of the Gaia BP/RP spectra, followed by a detailed outline of our sample selection criteria and pre-processing steps.

2.1 XP: The Gaia BP/RP low-resolution spectra

Gaia, DR3 provides low-resolution spectrophotometric data (XP spectra) for about 220 million sources, obtained with the Blue Photometer (BP, 330–680 nm) and Red Photometer (RP, 640–1050 nm) instruments on board the spacecraft (Gaia Collaboration 2023b). These spectra cover the optical to near-infrared range with a wavelength-dependent resolving power of R~20–100.

Unlike conventional spectroscopic data consisting of discrete flux measurements, Gaia XP spectra implement a continuous functional representation through basis function decomposition (De Angeli et al. 2023). Each spectrum is parameterized via 55 coefficients derived from orthonormal Gauss-Hermite functions, optimized to capture both broadband morphology and higher-frequency spectral features with minimal information loss. For analytical tractability, we transform these coefficient-based representations to wavelength-flux space using GaiaXPy¹, facilitating direct identification and measurement of diagnostic spectral features as described in Appendix A. This transformation yields calibrated flux values in physical units (W nm⁻¹ m⁻²) spanning the operational wavelength range (390–1000 nm). The low-resolution XP spectrophotometry provides sufficient spectral discrimination capability for identifying key stellar atmospheric diagnostics (Zhang et al. 2023; Li et al. 2024) while maintaining the statistical power of Gaia’s extensive survey volume, rendering it suitable for our binary census and characterization program.

2.2 Sample selection criteria

To construct our parent sample, we executed the following query on the Gaia DR3 database:

SELECT *
FROM gaiadr3.gaia_source
WHERE parallax > 1
      AND phot_bp_mean_mag < 20
      AND bp_rp BETWEEN 0 AND 5
      AND parallax_over_error > 10
      AND phot_g_mean_mag + 5*log10(abs(parallax)/100) > 3
      AND has_xp_continuous = ‘true’

Our sample selection procedure implements multi-stage filtering protocols: First, we establish a volume-limited sample within a 1 kpc heliocentric radius (parallax ϖ > 1 mas) containing objects with available BP/RP (XP) continuous spectra and B-R color indices spanning 0 < B − R < 5.0, encompassing F-type to M-type stars. Second, we mitigate systematic flux calibration biases that affect faint red sources in the Gaia photometric system by implementing a magnitude threshold of B < 20 mag, where these systematic effects remain below our precision requirements. Third, we excise evolved stellar populations by imposing an absolute magnitude criterion M_G > 3.0 mag (un-dereddened), removing giants and supergiants while retaining the main sequence population. We focus on main-sequence stars to ensure homogeneity in our sample, as their spectral properties are primarily determined by mass and metallicity, simplifying the forward modeling in this study. In future work, we plan to incorporate stellar age explicitly, particularly for brighter stars, by forward-modeling spectra with mass, age, and [M/H] as parameters. Fourth, we enforce an astrometric quality constraints, requiring parallax signal-to-noise ratios ϖ/σ_ϖ > 10 to ensure accurate absolute magnitude determinations. The resulting parent sample comprises 35,831,031 unique sources with Gaia DR3 identifiers.

2.3 Training data

To build a clean training sample of main sequence stars for our neural emulator of XP spectra, we have implemented a multi-stage filtering process targeting the main sequence within ~333 parsecs (parallax> 3). Our selection methodology prioritizes purity while maintaining a representative population across the lower main sequence.

First, we applied a color-magnitude diagram (CMD) filter using PARSEC theoretical isochrones with parameters of 1 Gyr age and [M/H]= 0.2 metallicity, and a hard cut of M_G > 3 as shown in Figure 1. By excluding objects above this reference isochrone, we removed most of the main sequence binaries and pre-main sequence objects that would contaminate our sample. Second, we remove the white dwarf main sequence binary (WDMS) candidates selected from Li et al. (2025). As a final quality assurance measure, we required all spectra to have an average signal-to-noise ratio greater than 70. The resulting clean sample contains 746 548 main sequence stars with high quality XP spectra.

Fig. 1 Hertzsprung–Russell diagram of training sample stars. The color scale indicates the logarithm of the number of stars in each bin. The solid orange line represents a PARSEC isochrone with metallicity ([M/H] = 0.2, age of 1 Gyr), while the dashed orange line shows the same isochrone for equal-mass binary systems (mass ratio q = 1) as a reference. The blue dotted lines indicate iso-mass tracks, with masses ranging from 0.1 M⊙ to 1.4 M⊙ as labeled.

3 Data-driven spectral model

In this section, we present our data-driven approach for modeling the XP spectra of both single and MSMS stars using neural networks (NN). Our methodology encompasses two interconnected components: (1) a single-star spectral emulator that maps stellar parameters directly to observed spectra in Subsection 3.1, and (2) a binary star model that combines individual stellar spectra based on fundamental stellar properties in Subsection 3.2. Our NN architecture employs a simple yet streamlined design for both single and binary star modeling. For single stars, we map photometric stellar parameters directly to spectral flux using a neural network. For binary systems, we transformed physical parameters into component-specific stellar properties, and combine the resulting spectra through linear superposition.

To enhance the capability of modeling single-star XP spectra, we iteratively cleaned the training set three times, removing outliers based on fit quality metrics. This process involved evaluating the χ² statistic for each spectra in the training set, where high χ² values indicated poor agreement between observed XP spectra and the NN’s predicted spectra. Outliers, defined as spectra with χ² values exceeding 3, were removed to ensure the training data consisted of representative single-star spectra. Three iterations of this cleaning were performed to refine the dataset progressively, each time retraining the NN to improve its accuracy in mapping stellar parameters.

We adopt photometric quantities directly related to the Hertzsprung-Russell (HR) diagram as our stellar labels. Specifically, we used

$$\theta =\left(B-R,M_G\right),$$(1)

where B − R is the dust-corrected color index, and M_G is the absolute magnitude in the Gaia G band, corrected for extinction using 3D dust maps (Green 2018; Edenhofer et al. 2024).

The de-reddening process follows the method of Li et al. (2025). We utilize the dustmaps package (Green 2018) to query a 3D dust map and derive the extinction parameter E². We adopt the average extinction curve from Zhang et al. (2023), combined with the transmission curves of the Gaia passbands, to compute the extinction corrections.

This choice of M_G is motivated by the fact that M_G provides a reliable and physically meaningful representation of low-mass stars (M_G ≥ 4–5). In contrast, log g is not an ideal descriptor for low-mass stars because it becomes less sensitive to variations in luminosity in this regime. In addition, our sample spans a wide range of spectral types, from M-type to G-type stars, making it difficult to define a homogeneous set of T_eff labels. The parameters (T_eff, log g, [M/H]) of M-type stars are often parameterized differently from previous spectral types (e.g., Li et al. 2021b; Qiu et al. 2024), leading to inconsistencies when using T_eff as a universal label. Furthermore, our approach is purely data-driven, relying on observational data rather than theoretical assumptions. By focusing on photometric quantities such as B − R and M_G, we ensure that our model is both reliable and applicable to the full range of stellar types in our sample.

Fig. 2 Demonstration of our NN spectral emulators for single main-sequence stars. Each panel shows the predicted XP spectra with different stellar parameters, where the x-axis is the wavelength and the y-axis is the logarithmic flux (10−20 W m−2 nm) normalized to a distance of 10 pc. Top panels: spectra at fixed absolute G magnitude (MG) with different B − R color. Bottom panels: spectra at fixed B-R color with varying MG.

3.1 Spectral emulator

The neural network architecture is designed to be as simple as possible, consisting of three hidden layers, each containing 16 neurons activated by Rectified Linear Units (ReLU). The configuration optimizes computational efficiency and spectral modeling. While XP spectra in logarithmic absolute flux space exhibit approximately linear relationships with absolute magnitude and colors, the spectral features and detailed line profiles require nonlinear transformations to accurately represent the underlying physics of stellar atmospheres.

The neural network defines a function, f , that maps from stellar parameters to the predicted XP spectrum:

$$f\left(\lambda ;\theta ,w\right)=F\left(\lambda \right),$$(2)

where F(λ) represents the flux at wavelength λ, θ are the stellar parameters, and w denotes the trainable parameters of the network.

We trained the network using the Adam optimizer with a learning rate of 10⁻⁴ and a batch size of 16 384. The loss function is defined as the chi-squared (χ²) statistic, which measures the discrepancy between the observed and predicted spectra:

$$\begin{array}{lll}L\hfill & =\hfill & \frac{1}{N}{\displaystyle \sum _{i=1}^N{\chi }_i^2}\hfill \\ \hfill & =\hfill & \frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(F_i^{\text{obs}}\left(\lambda \right)-f\left(\lambda ;{\theta }_i,w\right)\right)}^T{C}_i^{-1}\left(F_i^{\text{obs}}\left(\lambda \right)-f\left(\lambda ;{\theta }_i,w\right)\right),}\hfill \end{array}$$(3)

where $F_i^{\text{obs}}(\lambda )$ represents the observed spectrum of the i-th star, C_i is the covariance matrix for the i-th observation, and N is the total number of stars in the training set. Figure 2 shows the single-star spectral emulator as a function of B−R color and M_G.

3.2 Binary spectral model

3.2.1 Binary system parameterization

For modeling binary star systems, we performed stellar parameterization based on stellar properties. Each binary system is characterized by three primary parameters:

$$\varphi =\left(m_1,{[\text{M}/\text{H}]}_{\text{photo}},q\right),$$(4)

where m₁ represents the mass of the primary star, [M/H]_phot is the metallicity (assumed identical for both components), and q = m₂/m₁ is the mass ratio, constrained to 0 < q 1. This parameterization captures the key properties of binary≤ systems while minimizing the number of free parameters.

Fig. 3 Demonstration of NN spectral emulators for single main-sequence stars. Top panels: spectra at fixed stellar mass with different [M/H]phot. Bottom panels: spectra at fixed [M/H]photwith varying stellar mass.

3.2.2 Mapping stellar parameters to the HR diagram

To translate the physical parameters into observable quantities on the H-R diagram, we use the PARSEC stellar evolution models (Bressan et al. 2012; Chen et al. 2014) with stellar ages of 5 Gyr. We implement this mapping using a neural network function ℳ that captures the nonlinear transformations:

$$\theta =\mathcal{M}\left(m_1,{[\text{M}/\text{H}]}_{\text{photo}}\right),$$(5)

where θ = (B − R, M_G) represents the position in the H-R diagram expressed as dust-corrected color index and absolute magnitude.

The photometric properties of each component in a binary system are computed separately:

$${\theta }_1=\mathcal{M}\left(m_1,{[\text{M}/\text{H}]}_{\text{photo}}\right),$$(6)

$${\theta }_2=\mathcal{M}\left(q,m_1,{[\text{M}/\text{H}]}_{\text{photo}}\right)$$(7)

Here, [M/H]_photo denotes the photometric metallicity derived from broadband colors. The photometric metallicity reflects the integrated effect of metal content on the stellar atmosphere’s opacity and emergent flux, which may differ from element-specific spectroscopic determinations. Figure 3 shows the emulated XP spectra as functions of stellar mass and [M/H]_photo. Full details of this mapping implementation and its validation are provided in Appendix B.

3.2.3 Combined binary spectrum

Once the H-R diagram positions of both components are determined, we predicted the individual spectra using the same mapping function f from Section 3.1:

$$F_1\left(\lambda \right)=f\left({\theta }_1\right),$$(8)

$$F_2\left(\lambda \right)=f\left({\theta }_2\right).$$(9)

The combined spectrum of the unresolved binary system is then computed as the sum of the individual component spectra:

$$F_{\text{binary }}\left(\lambda ;m_1,{[\text{M}/\text{H}]}_{\text{photo }},q\right)=F_1\left(\lambda \right)+F_2\left(\lambda \right).$$(10)

This additive model assumes the flux contributions from the two stars are linearly superimposed, which is a valid approximation for non-interacting binary systems. The combined emulated binary spectra are shown in Figure 4.

3.2.4 Extinction modeling

We model extinction as a wavelength-dependent attenuation of stellar flux following standard extinction laws. The extinction is characterized by a parameter E, which scales the universal extinction curve R(λ). The observed flux is related to the intrinsic flux by:

$$F_{\text{obs}}\left(\lambda \right)=F_{\text{intrinsic}}\left(\lambda \right)\cdot \text{exp}\left(-E\cdot R\left(\lambda \right)\right),$$(11)

where R(λ) represents the relative extinction at each wavelength. For our analysis, we adopt the average extinction curve derived by Zhang et al. (2023), which was determined from a large sample of Gaia XP spectra spanning diverse Galactic environments.

Fig. 4 Comparison of single-star and binary-star XP spectral models for systems with varying mass ratios. Black lines represent the single-star model, while the color gradient represents binary-star models with varying mass ratios (q = m2/m1), ranging from q = 0.2 (blue) to q = 1.0 (red). The binary models incorporate consistent [M/H]phot, and distance for both components.

3.2.5 Parallax error propagation

We normalize all XP spectra to the flux at 10 pc using the parallax from Gaia DR3. Since the parallax measurement has uncertainties that cannot be neglected, we perform error propagation as follows:

$$f_{10}=f_{\text{obs}}{\left(\frac{100}{\varpi }\right)}^2,$$(12)

where f₁₀ is the flux normalized to 10 pc, f_obs is the observed flux at the star’s distance, and ϖ is the parallax in milliarc-seconds (mas). The error propagation process accounts for the uncertainty in parallax measurements when normalizing XP flux spectra to a standard distance of 10 pc.

We propagate the errors using two components:

$$C_{f_{10}}=g_{\varpi }^2{C}_{f_{\text{obs}}}+J{J}^T{\sigma }_{\varpi }^2,$$(13)

where g_ϖ is defined as the scaling factor,

$$g_{\varpi }={\left(\frac{100}{\varpi }\right)}^2.$$(14)

The first term scales the original covariance matrix $\left(C_{f_{\text{obs}}}\right)$ by the square of the scaling factor, while the second variance term accounts for parallax uncertainty. This second term is computed using the Jacobian (J) of the scaling function with respect to parallax, multiplied by the square of the parallax error (σ_ϖ):

$$J=\frac{\partial g_{\varpi }}{\partial \varpi }{f}_{\text{obs}}=-\frac{20000}{{\varpi }^3}{f}_{\text{obs}}.$$(15)

The resulting covariance matrix $(C_{f_{10}})$ incorporates both the scaled original uncertainties and the additional uncertainty introduced by the parallax measurement. This error propagation ensures that measurement uncertainties are properly accounted for in our model fitting and parameter inference.

3.3 Maximum likelihood parameter estimation

We employed a likelihood-based parameter inference of both single and binary stellar systems with a framework below.

3.3.1 Likelihood formulation

Assuming Gaussian uncertainties in the observed fluxes, we defined the likelihood function as:

$$\mathcal{L}\left(\Phi \right)\propto \text{exp}\left(-\frac{1}{2}{\chi }^2\left(\Phi \right)\right),$$(16)

where Φ represents either single-star parameters (θ, E) or binary parameters (m₁, [M/H]_photo, q, E). The χ² statistic is defined as:

$${\chi }^2\left(\Phi \right)=\text{Δ}{F}^T{C}_f^{-1}\text{Δ}F,$$(17)

where ∆F is the vector of residuals between observed and model-predicted fluxes, and C_f is the covariance matrix of flux uncertainties.

3.3.2 Modified likelihood with prior information

To incorporate prior information about extinction, we maximized a modified likelihood function:

$$\text{ln}{\mathcal{L}}_{\text{mod}}\left(\Phi \right)=-\frac{1}{2}{\chi }^2\left(\Phi \right)-\frac{1}{2}{\lambda }_p\cdot {\left(E-E_{\text{ref}}\right)}^2,$$(18)

where the second term functions as a penalty encouraging the extinction parameter E to remain close to the reference value E_ref from the extinction map (Edenhofer et al. 2024). The parameter λ_p, set to 0.05, controls the strength of the extinction constraint in the modified likelihood function. This term acts as a regularization, anchoring the fitted extinction E to the reference value E_ref from the dust map (Edenhofer et al. 2024).

3.3.3 Optimization procedure

For parameter estimation, we employ the Trust Region Reflective (TRF) algorithm (Coleman & Li 1994), which is effective for constrained optimization problems. The TRF method constructs a local quadratic approximation of the objective function within a trust region, ensuring stability and avoiding issues such as overstepping or divergence that can occur with other optimization methods.

The parameter space exploration is bounded by physical constraints:

• Primary mass m₁ ∈ [0.1, 1.5] M_⊙

• Photometric metallicity [M/H]_phot ∈ [−2.5, 0.7] dex

• Mass ratio q strictly satisfies 0.1 < q < 1.0

The optimal parameters are obtained by solving:

$$\stackrel{‸}{\Phi }=\text{arg}\underset{\Phi }{\text{min}}\left(-\text{ln}{\mathcal{L}}_{\text{mod}}\left(\Phi \right)\right).$$(19)

In our implementation, we compute the residual vector and extinction penalty separately, then combine them into an augmented residual vector. This allows us to treat the extinction penalty as an additional “observation” in our least-squares formulation while rigorously enforcing the parameter boundaries through the TRF algorithm’s constraint-handling capabilities.

Fig. 5 Validation of parameter recovery for the MSMS mock data of S/N= 50. Top row: comparison between true and inferred values for primary mass (m1, left), [M/H]phot([M/H], middle), and mass ratio (q, right). The hard limits in the right panel result from enforcing q < 1.

4 Validation

To validate our binary parameter recovery method, we constructed a suite of mock binary systems designed to span the full range of expected stellar parameters. Since we already had the single star emulator as described in Section 3, we first generated stellar labels and then used the single star emulator to mock the single and binary XP spectra. We generated 10 000 synthetic binary systems using a uniform random distribution for primary masses (m₁) ranging from 0.1 to 1.3 M_⊙ and [M/H]_phot ranging from −2.0 to +0.5 dex. Secondary masses (m₂) were similarly drawn from a uniform distribution between 0.1 and 1.3 M_⊙. We then performed random pairing of primary and secondary stars, followed by reordering when necessary to ensure m₁ > m₂ in all cases, resulting in mass ratios (q = m₂/m₁) spanning the full physical range from 0.1 to 1.0. For each synthetic binary, we generate mock Gaia XP spectra from the XP emulator with Gaussian noise properties based on the Gaia XP error model, incorporating the appropriate signal-to-noise. These mock spectra were then processed through our parameter inference pipeline, allowing us to quantify biases and scatter in the recovery of primary mass, metallicity, and mass ratio across the entire parameter space, as shown in Figure 5.

Figure 5 demonstrates the accuracy of our parameter recovery for mock binary systems at S/N=50. The top panels show the correlation between true and inferred parameters for primary mass (m₁), [M/H]_phot, and mass ratio (q), while the bottom panels display the residuals. Primary mass recovery is excellent with minimal bias (0.001) and scatter (0.049), showing that our model accurately constrains the dominant component. [M/H]_phot is also well recovered with only slight bias (0.003) and moderate scatter (0.192), which is reasonable given the metallicity’s more subtle spectral effects. The mass ratio exhibits somewhat larger scatter (0.137) and bias (0.019), particularly for low q values where the secondary’s spectral contribution is minimal. The residuals show no systematic trends across the parameter ranges, confirming our binary model’s robustness for systems with moderate to high signal-to-noise ratios.

The best-fit masses show larger scatter for stars more massive than 1 M_⊙, as these stars evolve on the HR-diagram during their main-sequence lifetime. This highlights that we cannot neglect the effects of stellar ages for more massive stars. In future work, we will further refine our approach by explicitly modeling colors and magnitudes as functions of mass, age, and metallicity to better account for evolutionary effects.

4.1 Classification methodology

To identify binary systems, we compare the goodness-of-fit between single-star and binary models using the χ² statistic. A source is classified as a binary if

$$\text{Δ}{\chi }^2={\chi }_{\text{single}}^2-{\chi }_{\text{binary}}^2>S,$$(20)

where S represents the optimal threshold for binary classification.

The threshold S is calibrated to balance completeness and purity, with an emphasis on minimizing false positives. We optimize using the F_β-score with β = 0.5:

$$F_{\beta }=\left(1+{\beta }^2\right)\cdot \frac{\text{Purity}\cdot \text{Completeness}}{\left({\beta }^2\cdot \text{Purity }\right)+\text{Completeness}}.$$(21)

Setting β < 1 explicitly prioritizes purity over completeness, reducing contamination from single stars incorrectly classified as binaries.

4.2 Parameter-dependent thresholds

We determine optimal thresholds within bins of primary mass (m₁) and [M/H]_phot and mass ratio (q) rather than applying a uniform threshold across all stellar types. The optimal threshold S varies across these parameters as shown in Table C.1, reflecting how binary detectability depends on stellar properties.

For validation, we generate a mock dataset of 100 000 spectra (50 000 single and 50 000 binary) using our spectral emulators. We randomly sample 10 000 stars with primary masses ranging from 0.1 to 1.5 M_⊙, photometric metallicities ([M/H]_photo) from −2.0 to +0.5 dex, and for binary systems, mass ratios (q) from 0.1 to 1.0. Each spectra is then simulated at five different signal-to-noise ratios (S/N = 5, 10, 20, 50, 100) to assess classification performance across varying data quality conditions. We add Gaussian noise to the synthetic fluxes according to these S/N levels to realistically mimic observational uncertainties.

For each mock spectrum, we perform both single-star and binary model fits to obtain best-fit parameters Φ_single and Φ_binary, then calculate the χ² values and their difference (∆χ²). Within each bin of primary mass, metallicity, and mass ratio, we determine the threshold S that maximizes the F_β-score. Figure 6 shows the performance metrics as a function of S for all mock dataset. The optimal point is set to maximize the F_β value, which in this case is 2.37. This approach ensures our classification accounts for the varying sensitivity of the ∆χ² statistic across different regions of parameter space.

Binaries with nearly equal-mass components typically require lower thresholds for identification as they produce more spectral deviations from single-star models. Conversely, systems with low-mass ratios require more stringent thresholds to avoid false positives,. The resulting parameter-dependent threshold map in Table C.1 enables binary classification across the entire dataset while balancing purity and completeness.

Fig. 6 Performance metrics as a function of the difference between χ2 single and χ2 binary values. The plot shows purity (precision, purple line), completeness (recall, green line), and Fβ score (orange line) across different threshold values, where β = 0.5. The vertical dashed line is the χ2 threshold that maximum the Fβ score. Beyond this threshold, completeness decreases most rapidly, followed by Fβ score, while purity remains high.

4.3 Verification with binary catalogs

The Gaia DR3 Non-Single Star (NSS) catalog (Gaia Collaboration 2023b) identified 813 687 multiple systems, including 169 227 astrometric orbital solutions, 220 372 spectroscopic solutions, 87 073 eclipsing binaries, and 338 215 acceleration solutions (Gaia Collaboration 2023a). Spectroscopic binaries are most detectable at periods of 1–1000 days, limited by stellar rotation at short periods and the DR3 baseline at longer periods. The NSS catalog represents more than orders of magnitude increase in the number of well-characterized binary systems compared to pre-Gaia, catalogs (Gaia Collaboration 2023b).

Our verification as shown in Fig. 8 demonstrates that the vast majority of sources in the NSS catalog show preference for binary-star models when fitting their XP spectra, providing independent validation of our methodology. This correlation is particularly pronounced for double-lined spectroscopic binaries (SB2 and SB2C types), where the χ² improvement for binary models over single-star fits consistently exceeds our detection threshold. The stronger signal for SB2 systems is expected, as these binaries typically have more comparable flux ratios between components, making their composite spectral signatures more distinctive. Fig. 7 illustrates case studies of both SB2 and astrometric binaries from the NSS catalog, demonstrating how the binary model (orange dashed line) provides superior fits to the observed XP spectra compared to single-star models (blue solid line).

As shown in Table 1, for Gaia NSS catalog, our methodology correctly classifies 67% of astrometric/orbital binaries, 81% of SB1, and 93% of SB2 systems as binaries based on our parameter-dependent thresholds in Subsection 4.2. The APOGEE binary sample (El-Badry et al. 2018), derived from ~20 000 main-sequence star spectra in APOGEE DR13, includes ~2500 binaries and ~200 triples identified through spectral fitting with composite models, alongside ~700 velocity-variable systems. For APOGEE, the correct classification rates are 52% for SB1 and 85% for SB2 systems, confirming the robustness of our spectral fitting approach. Additionally, only 34% of APOGEE’s single-star sample is classified as binaries in our XP fitting, further validating our results, as some of these “single” stars may be long-period q≈1 binaries indistinguishable from single stars in earlier models due to limited parallax precision.

The astrometric-only binaries (Orbital) show a more varied response in our analysis, with detection rates correlating strongly with the magnitude difference between components. These patterns across different NSS solution types provide insight into the sensitivity limits and complementary nature of spectrophotometric binary detection compared to spectroscopic and dynamic approaches. To improve usability, we also provide the cross-matched catalog between our binary fits and the NSS catalog in Appendix C.3.

We further validate the mass estimates derived from our XP fitting method by comparing them with those provided in the NSS catalog, specifically the gaiadr3.binary_masses table (Gaia Collaboration 2023a). This table contains pre-computed masses for the primary (m₁) and secondary (m₂) components of binary systems. The NSS mass estimates are derived using a variety of observational techniques, primarily combining high-precision astrometric, spectroscopic, and, in some cases, photometric (e.g., eclipsing binary) data. We adopt the binary mass values from all available solutions, including the solution type, such as ‘Orbital+SB2’ (astrometric orbit combined with double-lined spectroscopy), ‘EclipsingSpectro(SB2)’ (eclipsing binary photometry with SB2 spectroscopy), or ‘SB1+M1’ (single-lined spectroscopy with primary mass from isochrone fitting). These methods leverage Kepler’s Third Law and, where applicable, radial velocity semi-amplitudes and orbital inclinations to resolve individual masses. As shown in Fig. 9, the m₁ and m₂ values derived from our XP fitting method exhibit strong consistency with the NSS catalog’s mass estimates. This agreement validates the robustness of our approach, particularly given the NSS catalog’s reliance on diverse and complementary observational data.

4.4 Validation with dark companion systems

The identification of single stars represents the inverse challenge of binary detection. Achieving high purity in binary star identification inevitably comes at the cost of both binary sample completeness and single star sample purity. The chi-squared comparison method requires careful calibration, with threshold settings that vary as demonstrated in Section 5.2. To validate our classification approach, we sought benchmark “single-star-like” systems to verify whether our method correctly categorizes them. Such validation would be demonstrated by these systems showing minimal improvement in fit quality when modeled as binaries rather than single stars. For this purpose, we utilized two specialized samples.

First, we examined the neutron star companion sample from El-Badry et al. (2024), which contains 21 astrometric binaries with solar-type stars orbiting dark companions with masses near 1.4 M_⊙. Since these companions contribute virtually no light to the systems, their spectra should closely resemble those of single stars despite their binary nature. Second, we analyzed the white dwarf-main sequence (WDMS) binary sample from Li et al. (2025), which identifies approximately 30 000 WDMS systems from XP spectra. In most of these systems, the main sequence primary is at least 10 times brighter than the white dwarf secondary (Li et al. 2025), though some were flagged as MSMS binaries using the same chi-squared method.

The right panel in Figure 8 compares the chi-squared distribution of these samples against a control sample consisting of all Gaia sources within 1 kpc. The results demonstrate that for the neutron star sample, the binary model provides minimal improvement over the single-star model $\left({\text{log}}_{10}\left({\chi }_{\text{MS}}^2\right)-{\text{log}}_{10}\left({\chi }_{\text{MSMS}}^2\right)<0.2\right)$. Similarly, 84% of the WDMS sample shows no notable χ² improvement with the binary model, although some systems in this sample could potentially be MSMS binaries with white dwarfs as tertiary companions. This validation confirms that our method appropriately classifies systems with dark or faint companions as effectively single stars from a spectroscopic perspective, despite their astrometric binary nature.

We also validate the χ² difference for Gaia Black Hole 1 (BH1) as shown in Figure 10, we see that the total ${\chi }_{\text{single }}^2$ is 47.3, ${\chi }_{\text{binary }}^2$ is 53.3, which also proves our method can successfully identify it as a single star/dark companion. This validation confirms that our method appropriately classifies systems with dark or faint companions as effectively single stars from a spectroscopic perspective.

Fig. 7 Gaia XP spectra and best-fitting models for four MSMS binaries. In each panel, the black line represents the observed flux. The blue solid line shows the best-fit single-star model, while the orange dashed line represents the best-fit binary model. The pink dotted line corresponds to the primary component of the binary model, and the dark blue dotdash line shows the secondary component. For each system shown, the binary model is preferred as it achieves a lower χ2 value than the single-star model. Top row: example NSS SB2 binary. Bottom row: example NSS astrometric binaries.

Fig. 8 Distribution of the logarithmic χ2 improvement, $\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)$, showing the statistical evidence favoring binary over single-star spectral models. Left panel: distributions for Gaia DR3 NSS sources: systems with orbital solutions (blue), single-lined spectroscopic binaries (SB1/SB1C; orange), and double-lined spectroscopic binaries (SB2; green). Right panel: distributions for specialized populations with likely dark companions: neutron star binaries (El-Badry et al. 2024) (red) and WDMS (Li et al. 2025) (blue). The Gaia DR3 NSS SB2 sample (green dotted line) is shown for reference. Common elements in each panel: the distribution for all stars within 1 kpc (black dashed line) acts as a control representing the general population. A vertical dashed line marks 0.2 for reference. Higher values of the χ2 improvement indicate stronger evidence for binarity; known binaries typically show significant positive values.

Fig. 9 Comparison of primary (left) and secondary (right) masses derived from the XP binary model (m1, m2) against those from the Gaia non-single star catalog (m1, m2). The dashed line represents the one-to-one relation.

Fig. 10 Gaia XP spectra and best-fitting models for Gaia BH1. The black line represents the observed flux. The blue solid line shows the best fit single-star model, while the orange dashed line represents the best fit binary model. Bottom inset panel: residuals (∆ log flux) for both models.

5 The MSMS binary catalog

5.1 Data products

We have compiled and published the MSMS binary catalog based on Gaia DR3 XP spectra. This catalog contains 14 million binary candidates identified from our parent sample of 35 831 031 stars within 1 kpc of the Sun. The candidates were selected through a combination of spectroscopic modeling techniques that compare the χ² values between single-star and binary-star model fits.

For each entry, we include the best-fit results from both single-star and binary spectral models, allowing users to evaluate the binary classification (12). The binary model parameters include primary mass, system metallicity, mass ratio, and derived photometric properties such as colors and absolute magnitudes for both components.

The full MSMS binary catalog dataset, along with NN weights enabling users to forward model XP single/binary spectra according to their requirements, is now available for public access Li (2025)³.

5.2 Binary candidate selection

To identify binary systems from our stellar sample, we employ a statistical approach based on the χ² difference between single-star and binary-star model fits. Our method leverages the model structure described in Section 3, extending it to accommodate binary systems by combining the fluxes of two stars with different masses but sharing the same [M/H]_phot, and distance. This approach is particularly sensitive to binaries with high mass-ratios, which were identified as outliers during the self-cleaning process described in Section 3.

We begin with a parent sample of 35 million stars within 1 kpc. For each star, we perform two separate model fits: one assuming that the object is a single star and another assuming it is a binary system. The chi-square difference (∆χ²) between these two models provides a statistical measure of how much the data favors the binary hypothesis over the single-star hypothesis.

To establish appropriate thresholds for binary classification, we divide our sample into bins based on primary mass ([M/H]_phot) and mass ratio (q), as shown in Table C.1. These thresholds are calibrated to balance completeness and contamination in our binary identification, accounting for the varying sensitivity of our method across different stellar parameters.

A star is classified as a binary candidate if the χ² difference exceeds the threshold corresponding to its mass bin:

$$\text{Δ}{\chi }^2={\chi }_{\text{single}}^2-{\chi }_{\text{binary}}^2>\text{threshold}\left(m_1,{[\text{M}/\text{H}]}_{\text{phot}},q\right).$$(22)

Fig. 11 shows the distribution of ${\log }_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)$ values for our parent sample. Binary candidates systematically exhibit larger chi-squared differences, with the distribution showing a clear bimodality.

Using the thresholds defined in each m₁-[M/H]_phot-q bin, we identify 14 183 917 binary candidates from our parent sample, corresponding to an overall binary fraction of approximately 39%. This is consistent with expected binary fractions from previous studies (e.g. Gao et al. 2014; Liu 2019; Moe et al. 2019), although the exact fraction varies with stellar mass and metallicity environment.

Fig. 11 Distribution of ${\log }_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)$ for a subset of Gaia DR3 XP spectra, illustrating the statistical separation between single-star and binary models. The histogram shows the frequency of sources as a function of ∆χ2, with the classifiation based on an optimal threshold S (calibrated to maximize the Fβ-score with β = 0.5) used to classify binary systems.

Fig. 12 Venn diagram showing the hierarchical relationship between our sample selections. The parent sample consists of 35 million stars within 1 kpc from Gaia DR3. From this population, we identify 14 million binary candidates based on spectral model fitting criteria. The golden sample is identified through selection criteria in Table 2.

5.3 Golden binary sample

To facilitate studies requiring higher confidence binary identifications, we defined a “golden sample” of 959 594 binary systems, representing the most reliable binary candidates. These golden sample members satisfy stringent selection criteria and quality cut detailed in Table 2.

The “golden sample” of 959,594 systems, selected from an initial catalog of 35 831 031 systems (Table C.2), represents the highest confidence binary candidates identified in this work. These systems were chosen to ensure robust binary identification and focus on well-characterized main-sequence pairs suitable for detailed analysis.

These golden sample members satisfy stringent criteria including a significant chi-square improvement with ${\log }_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)>0.2$ , good fit quality for both single or binary models with χ² < 3 × 61 (3 per degree of freedom), successful convergence of model fits, primary masses restricted to M₁ < 1.0 M_⊙, [M/H]_phot from the binary fit > −1, absolute G-band magnitude M_G ≥ 4.5 to focus primarily on main sequence stars, and extinction E_binary < 0.2 to ensure the extinction is not too severe. Additionally, systems must have reliable mass ratios between 0.4 ≤ q ≤ 1.0, and flux ratios between 0.2 and 1.0. These cuts are based on practical considerations from mock tests and are not overly stringent, though they may exclude some binaries.

Fig. 12 illustrates the relationship between our parent sample, binary candidates, and golden binary sample using a Venn diagram, highlighting the nested nature of these increasingly selective subsets. Fig. 13 shows HR diagrams for our binary systems across different mass ratio ranges. The distinct positions of equal-mass binaries (right panel) compared to lower mass-ratio systems (left panel) in the color-magnitude diagram support the effectiveness of our classification approach, with nearly equal-mass binaries appearing systematically brighter than their single-star counterparts of the same color. Our catalog includes the stellar parameters of both the primary and the inferred secondary companion for all binary candidates, along with the ∆χ² value and classification confidence metrics. The MSMS binary catalog represents an extensieve uniform samples of binary systems to date, enabling population studies across different Galactic environments and investigating the dependence of binary fraction on stellar and environmental parameters.

6 Discussion and conclusions

We discuss binary star systems identified from Gaia DR3 XP spectra in the following sections. In Section 6.1, we first validate our spectral fitting results with enhanced parallax error treatment, demonstrating the robustness of our methodology. We then examine potential caveats and limitations of our binary analysis in Section 6.2.

6.1 Parallax errors

Two factors contribute to the underestimation of uncertainties for true binaries in our analysis. First, treating binaries as single stars introduces errors in parallax estimates (Gaia Collaboration 2023a). We validate the fitting results using parallax solutions derived from binary orbit assumptions. Second, binaries exhibit poor astrometric fits due to orbital motion, leading to underestimated uncertainties in Gaia’s astrometric solutions, including parallax. We apply a simulation-driven method (El-Badry 2025) to correct the parallax error and compare the results with our findings.

Fig. 13 HR diagrams for golden binary sample of mass-ratio >0.9. The background blue shows the density distribution of selected single stars in, with the scale bar indicating the logarithmic number of stars per bin. The pink solid and dashed lines show PARSEC isochrones (Bressan et al. 2012) for [M/H]=0, and q = 1 binary configurations.

6.1.1 Comparison with parallax from binary solutions

The NSS catalog fits complex models to describe the behavior of non-single stars, accounting for orbital motion in binary systems (Gaia Collaboration 2023a). It includes solutions for astrometric, spectroscopic, and eclipsing binaries, providing parallax estimates derived from these models. The comparison between DR3 source and NSS parallax solutions (Fig. 14) reveals a mean difference of 0.004 mas with an RMS scatter of 0.3 mas, indicating that source parallaxes slightly underestimate the true values of the NSS binary systems. We further compare the χ² improvement using different parallax sources in Figure 15, where the right panel demonstrates that the difference between NSS parallax and source parallax is negligible. This suggests that for sources without NSS orbital solutions, using the DR3 source parallax for XP fitting may be a reasonable approach.

Fig. 14 Comparison of parallax from Gaia DR3 source catalog and NSS calalog. The dashed red line represents the one-to-one relation, and the color scale denotes the number density.

6.1.2 Inflating parallax uncertainties

To address underestimation, we adopt the inflation model from El-Badry (2025) to adjust uncertainties based on the Renormalized Unit Weight Error (RUWE), a metric of astrometric fit quality where higher values indicate poorer fits. Many binaries with close companions have elevated RUWE values due to photocenter oscillations from orbital motion, which single-star models cannot accommodate (Gaia Collaboration 2023a). For binaries with orbital periods of a few years, photocentric acceleration further increases RUWE. The model of El-Badry (2025) adjusts reported uncertainties as a function of RUWE, parallax, and apparent magnitude, transforming error distributions into Gaussian forms.

We compare the χ² improvement $\left({\text{log}}_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)\right)$ in Figure 15. For astrometric (orbital) binaries, the difference before and after the parallax inflation is negligible. However, for SB1 systems, we observe modest differences: For binaries with ${\text{log}}_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)<0.5$, the XP fit incorporating parallax error inflation demonstrates better performance than without inflation. Specifically, binaries previously showing ${\log }_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)$ between 0 and 0.2 shift to the 0.2–0.5 bin after applying parallax error inflation, indicating improved classification capability. For binaries with ${\text{log}}_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)>0.5$, parallax error inflation produces no significant difference in results.

Fig. 15 Distribution of the logarithmic χ2 improvement for binary models. In each panel, blue lines denote astrometric binary stars classified by the NSS catalog, while orange lines represent SB1 binaries. In the left panel, the dashed lines are derived from the original DR3 source catalog parallax errors, while solid lines show results after error inflation following the model of El-Badry (2025). In the right panel, the dashed lines use parallaxes from the DR3 source catalog, while solid lines use parallaxes from the NSS catalog, which incorporate binary orbital solutions in their parallax estimates.

6.2 Caveats and limitations

6.2.1 Extinction modeling

Our extinction modeling relies on an average extinction curve derived from Gaia XP spectra (Zhang et al. 2023). This approach, while effective for population-level analysis, may not fully capture the variation in extinction properties across different Galactic environments. Recent studies by Zhang & Green (2025) demonstrate that R(V), which characterizes the slope of the extinction curve in optical wavelengths and is related to the typical size of dust grains, can vary depending on the line of sight. This variation could introduce systematic uncertainties in our derived stellar parameters, particularly in regions with anomalous dust properties or high extinction values. Future work could incorporate spatially dependent extinction curves to further refine the binary parameter estimates.

6.2.2 Low-mass secondary extrapolation

For binary systems with secondary masses (m₁ × q) below 0.1 M_⊙, our parameter estimates should be treated with caution. The stellar structure models we employ (Bressan et al. 2012) have a lower-mass limit of 0.1 M_⊙, meaning that binary fits with low-mass secondaries rely on extrapolation from our NN emulators. We provide a quality flag flag_m2_extrap in our catalog to identify these systems.

6.2.3 Evolutionary effects for higher-mass stars

For stars more massive than approximately 0.8 M_⊙, which corresponds to the highest mass that remains on the main sequence after 12.5 Gyr at [Fe/H] = –2 (Paxton et al. 2011; Choi et al. 2016), evolutionary effects become increasingly important. These stars likely have evolved from their zero-age main sequence positions in the color-magnitude diagram, making our modeling from mass and metallicity to color-absolute magnitude less accurate without accounting for age variations. Current parameter estimates for these higher-mass systems do not incorporate age as a variable, potentially introducing systematic biases in the inferred mass ratios and metallicities. According to mock tests, stars of M > 0.8 M_⊙ have a systematic error < 0.05 M_⊙ due to the age change from 1 Gyr to 10 Gyr. A future refinement of this work should consider stellar ages and initial masses of binary systems when modeling their spectral and photometric properties.

6.3 Summary

This study presents the MSMS Catalog, identifying ~14 million binary star candidates from Gaia DR3 BP/RP spectra. The analysis uses spectrophotometric data for 35 million stars within 1 kpc of the Sun, selected from 220 million Gaia DR3 sources.

Binaries are identified when the binary model yields a lower χ² than the single-star model, with a parameter-dependent χ² improvement threshold for mass ratios 0.5 ≤ q ≤ 1.0. A subset of approximately one million systems (“golden sample”) is defined using a stricter chi-squared improvement, achieving 93% completeness of SB2s validated using the Gaia non-single star catalog.

This approach overcomes limitations of spectroscopic methods, which rely on radial velocity shifts and favor short-period binaries (periods < 1000 days). By comparing spectral fits, the method detects binaries across mass ratios 0.5 ≤ q ≤ 1.0. Validation with synthetic spectra confirms recovery of primary mass and[M/H]_photwith errors < 10%. The 14 million candidates reveals a binary fraction of 39% among main-sequence stars, with mass- and metallicity-dependent trends in the mass-ratio distribution.

Acknowledgements

We thank the anonymous referee for their constructive comments that improved the manuscript. J.L. thanks Coryn Bailer-Jones for constructive discussion and insights. We acknowledge support from the European Research Council through ERC Advanced Grant No. 101054731. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Y.S.T. is supported by the National Science Foundation under Grant AST-2406729.

References

Appendix A Sampled XP spectra

To accurately process XP spectra, we must account for covariances between observed fluxes across wavelengths. Our sampling range of 392–992 nm and step size of 10 nm were chosen to trim noisy spectral edges and ensure positive-definite covariance matrices (Zhang et al. 2023).

The χ² statistic between predicted and observed fluxes is calculated as:

$${\chi }^2=\Delta f^T{C}_{f_{\text{obs}}}^{-1}\Delta f=|L\Delta f{|}^2,$$(A.1)

where ∆f ≡ f_pred − f_obs represents the residuals between predicted and observed fluxes.

To ensure numerical stability in this calculation, we diagonalize each covariance matrix using its eigendecomposition:

$$C_{f_{\text{obs}}}=UD{U}^T,$$(A.2)

where U is orthonormal and D is diagonal. We impose a minimum value of 10⁻⁹ (in units of 10⁻³⁶ W² m⁻⁴ nm⁻²) on the diagonal elements of D (yielding $\widehat{D}$ ) and define $L\equiv {\widehat{D}}^{-1/2}{U}^T$.

The covariance matrix $C_{f_{\text{obs}}}$ incorporates several uncertainty components:

$$\begin{array}{r}\hfill C_{f_{\text{obs}}}=C_{f_{\text{gaia}}}+\text{diag}{\left(0.005f_{\text{gaia}}\right)}^2\\ \hfill +{0.005}^2{f}_{\text{gaia}}{f}_{\text{gaia}}^T+{0.001}^2{f}_{\text{BP}}{f}_{\text{BP}}^T+{0.001}^2{f}_{\text{RP}}{f}_{\text{RP}}^T,\end{array}$$(A.3)

where f_gaia represents the sampled flux, and f_BP and f_RP are zero-padded vectors containing the sampled BP and RP spectra, respectively.

The terms in this expression account for the following:

1. The sample-space covariance matrix from Gaia DR3, denoted as $C_{f_{\text{gaia}}}$ .

2. A 0.5% independent flux uncertainty at each wavelength.

3. A 0.5% uncertainty in the zero-point calibration of the overall spectral flux.

4. A 0.1% uncertainty in the zero-point calibration of the BP spectrum.

5. A 0.1% uncertainty in the zero-point calibration of the RP spectrum.

Following Zhang et al. (2023), we inflate the original covariance matrix to avoid dominance by small uncertainties reported in Gaia DR3. The zero-point correction involves three components: the overall spectral zero-point and the relative zero-points for BP and RP spectra.

Appendix B From physical parameters to CMD

The mapping function ℳ that translates stellar physical parameters to H-R diagram coordinates is implemented as a neural network with three hidden layers, maintaining consistency with the spectral emulator architecture described in Section 3.1.

Figure B.1 validates our stellar parameter-to-observable mapping function ℳ, showing notable accuracy in translating physical parameters to H-R diagram positions. The correlation between true and predicted color indices (B-R) shows high precision with minimal bias (0.001) and scatter (0.021), demonstrating that our NN implementation accurately captures the nonlinear relations between mass, [M/H]_phot, and color. Similarly, absolute magnitude predictions (M_G) exhibit agreement with true values, with negligible bias (0.003) and low scatter (0.026). The residual plots (bottom panels) show no notable systematic trends across the parameter space, confirming that our implementation of the PARSEC evolutionary tracks through the NN mapping function successfully reproduces the photometric properties of stars with diverse masses and metallicities. This approach ensures that the physical parameters inferred from our spectral analysis can be converted to observable quantities for comparison with photometric data and for visualization on standard H-R diagrams.

Fig. B.1 Validation of the neural network-based mapping function ℳ that translates stellar physical parameters to H-R diagram coordinates. Left panels: Comparison between true and predicted B-R color, showing good agreement with minimal bias (0.001) and scatter (0.021). Right panels: Validation of absolute G-magnitude (MG) predictions, displaying negligible bias (0.003) and low scatter (0.026). The bottom panels show residuals (predicted minus true values) across the parameter range, demonstrating no systematic trends.

Appendix C The catalog

Appendix C.1 Binary candidate selection

We provide both the single and binary fitting results for ~35 million sources in Gaia database. For each source, we include the optimal parameters from both single-star and binary-star models, along with quality flags and derived quantities. The catalog contains 14,183,917 binary candidates, with a subset of 959,594 sources meeting our stricter criteria for the "golden sample." For each source, we provide the primary and secondary stellar parameters (masses, metallicities), the mass ratio, flux ratios in different bands, and quality metrics including χ² values and fit convergence flags. The catalog also includes Gaia parallax and extinction estimates from both our model fits and external dust maps. Users can filter the catalog based on the provided quality flags to select samples optimized for either completeness or purity depending on their scientific objectives.

Appendix C.2 Quality flags

The catalog includes quality flags to assess the reliability and convergence status of single- and binary-model fits. The fit_status_single and fit_status_binary codes indicate termination reasons: -1 for invalid inputs, 0 for exceeding function evaluations, and 1–4 for convergence via gtol, ftol, xtol, or multiple tolerances. The Boolean fit_success_single and fit_success_binary flags denote overall success (typically True for status codes ≥ 1, subject to additional criteria like physically plausible parameters). The flag_extrap_m2 highlights cases where the binary fit yields a secondary mass < 0.1 M_⊙, signaling potential extrapolation beyond the model’s calibrated range. These flags collectively help users gauge the robustness of fitted parameters and identify fits requiring scrutiny.

The least_squares optimizer in Scipy uses three key tolerance parameters to determine convergence:

gtol Gradient tolerance

• Optimization stops when the maximum gradient magnitude falls below this threshold.

• Corresponds to quality flag: fit_status_* = 1

ftol Function tolerance

• Stops when relative improvements in the cost function become negligible.

• Corresponds to quality flag: fit_status_* = 2

xtol Parameter tolerance

• Terminates when parameter adjustments become sufficiently small.

• Corresponds to quality flag: fit_status_* = 3

both ftol + xtol

• Stopping condition met for both function and parameter tolerances.

• Corresponds to quality flag: fit_status_* = 4

Typical default values are 10⁻⁸ for all tolerances. Lower values enforce stricter convergence.

Appendix C.3 Cross-match with Gaia non-single star catalog

In addition to our primary binary star catalog, we provide a supplementary dataset cross-matching our identified binary systems with the Gaia DR3 Non-Single Star (NSS) catalog. Table C.3 summarizes the key NSS parameters included in this cross-matched catalog. This cross-matched catalog is particularly valuable for systems with both spectrophotometric signatures (identified through our neural network approach) and dynamical evidence (from the Gaia NSS solutions), allowing for more robust mass ratio determinations and complete orbital characterization.

All Tables

All Figures

Fig. 1 Hertzsprung–Russell diagram of training sample stars. The color scale indicates the logarithm of the number of stars in each bin. The solid orange line represents a PARSEC isochrone with metallicity ([M/H] = 0.2, age of 1 Gyr), while the dashed orange line shows the same isochrone for equal-mass binary systems (mass ratio q = 1) as a reference. The blue dotted lines indicate iso-mass tracks, with masses ranging from 0.1 M⊙ to 1.4 M⊙ as labeled.

In the text

Fig. 2 Demonstration of our NN spectral emulators for single main-sequence stars. Each panel shows the predicted XP spectra with different stellar parameters, where the x-axis is the wavelength and the y-axis is the logarithmic flux (10−20 W m−2 nm) normalized to a distance of 10 pc. Top panels: spectra at fixed absolute G magnitude (MG) with different B − R color. Bottom panels: spectra at fixed B-R color with varying MG.

In the text

	Fig. 3 Demonstration of NN spectral emulators for single main-sequence stars. Top panels: spectra at fixed stellar mass with different [M/H]phot. Bottom panels: spectra at fixed [M/H]photwith varying stellar mass.
In the text

	Fig. 4 Comparison of single-star and binary-star XP spectral models for systems with varying mass ratios. Black lines represent the single-star model, while the color gradient represents binary-star models with varying mass ratios (q = m2/m1), ranging from q = 0.2 (blue) to q = 1.0 (red). The binary models incorporate consistent [M/H]phot, and distance for both components.
In the text

	Fig. 5 Validation of parameter recovery for the MSMS mock data of S/N= 50. Top row: comparison between true and inferred values for primary mass (m1, left), [M/H]phot([M/H], middle), and mass ratio (q, right). The hard limits in the right panel result from enforcing q < 1.
In the text

Fig. 6 Performance metrics as a function of the difference between χ2 single and χ2 binary values. The plot shows purity (precision, purple line), completeness (recall, green line), and Fβ score (orange line) across different threshold values, where β = 0.5. The vertical dashed line is the χ2 threshold that maximum the Fβ score. Beyond this threshold, completeness decreases most rapidly, followed by Fβ score, while purity remains high.

In the text

Fig. 7 Gaia XP spectra and best-fitting models for four MSMS binaries. In each panel, the black line represents the observed flux. The blue solid line shows the best-fit single-star model, while the orange dashed line represents the best-fit binary model. The pink dotted line corresponds to the primary component of the binary model, and the dark blue dotdash line shows the secondary component. For each system shown, the binary model is preferred as it achieves a lower χ2 value than the single-star model. Top row: example NSS SB2 binary. Bottom row: example NSS astrometric binaries.

In the text

Fig. 8 Distribution of the logarithmic χ2 improvement, $\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)$, showing the statistical evidence favoring binary over single-star spectral models. Left panel: distributions for Gaia DR3 NSS sources: systems with orbital solutions (blue), single-lined spectroscopic binaries (SB1/SB1C; orange), and double-lined spectroscopic binaries (SB2; green). Right panel: distributions for specialized populations with likely dark companions: neutron star binaries (El-Badry et al. 2024) (red) and WDMS (Li et al. 2025) (blue). The Gaia DR3 NSS SB2 sample (green dotted line) is shown for reference. Common elements in each panel: the distribution for all stars within 1 kpc (black dashed line) acts as a control representing the general population. A vertical dashed line marks 0.2 for reference. Higher values of the χ2 improvement indicate stronger evidence for binarity; known binaries typically show significant positive values.

In the text

	Fig. 9 Comparison of primary (left) and secondary (right) masses derived from the XP binary model (m1, m2) against those from the Gaia non-single star catalog (m1, m2). The dashed line represents the one-to-one relation.
In the text

	Fig. 10 Gaia XP spectra and best-fitting models for Gaia BH1. The black line represents the observed flux. The blue solid line shows the best fit single-star model, while the orange dashed line represents the best fit binary model. Bottom inset panel: residuals (∆ log flux) for both models.
In the text

Fig. 11 Distribution of ${\log }_{10}\left({\chi }_{\text{single }}^2/{\chi }_{\text{binary }}^2\right)$ for a subset of Gaia DR3 XP spectra, illustrating the statistical separation between single-star and binary models. The histogram shows the frequency of sources as a function of ∆χ2, with the classifiation based on an optimal threshold S (calibrated to maximize the Fβ-score with β = 0.5) used to classify binary systems.

In the text

	Fig. 12 Venn diagram showing the hierarchical relationship between our sample selections. The parent sample consists of 35 million stars within 1 kpc from Gaia DR3. From this population, we identify 14 million binary candidates based on spectral model fitting criteria. The golden sample is identified through selection criteria in Table 2.
In the text

	Fig. 13 HR diagrams for golden binary sample of mass-ratio >0.9. The background blue shows the density distribution of selected single stars in, with the scale bar indicating the logarithmic number of stars per bin. The pink solid and dashed lines show PARSEC isochrones (Bressan et al. 2012) for [M/H]=0, and q = 1 binary configurations.
In the text

	Fig. 14 Comparison of parallax from Gaia DR3 source catalog and NSS calalog. The dashed red line represents the one-to-one relation, and the color scale denotes the number density.
In the text

Fig. 15 Distribution of the logarithmic χ2 improvement for binary models. In each panel, blue lines denote astrometric binary stars classified by the NSS catalog, while orange lines represent SB1 binaries. In the left panel, the dashed lines are derived from the original DR3 source catalog parallax errors, while solid lines show results after error inflation following the model of El-Badry (2025). In the right panel, the dashed lines use parallaxes from the DR3 source catalog, while solid lines use parallaxes from the NSS catalog, which incorporate binary orbital solutions in their parallax estimates.

In the text

Fig. B.1 Validation of the neural network-based mapping function ℳ that translates stellar physical parameters to H-R diagram coordinates. Left panels: Comparison between true and predicted B-R color, showing good agreement with minimal bias (0.001) and scatter (0.021). Right panels: Validation of absolute G-magnitude (MG) predictions, displaying negligible bias (0.003) and low scatter (0.026). The bottom panels show residuals (predicted minus true values) across the parameter range, demonstrating no systematic trends.

In the text