Triage of the Gaia DR3 astrometric orbits. II. A census of white dwarfs

The third data release of Gaia was the first to include orbital solutions assuming non-single stars. Here, we apply the astrometric triage technique of Shahaf et al. to identify binary star systems with companions that are not single main-sequence stars. Gaia's synthetic photometry of these binaries is used to distinguish between systems likely to have white-dwarf companions and those that may be hierarchical triples. The study uncovered a population of nearly 3200 binaries, characterised by orbital separations on the order of an astronomical unit, in which the faint astrometric companion is probably a white dwarf. This sample increases the number of orbitally solved binary systems of this type by about two orders of magnitude. Remarkably, over 110 of these systems exhibit significant ultraviolet excess flux, confirming this classification and, in some cases, indicating their relatively young cooling ages. We show that the sample is not currently represented in synthetic binary populations, and is not easily reproduced by available binary population synthesis codes. Therefore, it challenges current binary evolution models, offering a unique opportunity to gain insights into the processes governing white-dwarf formation, binary evolution, and mass transfer.


INTRODUCTION
The detection of white dwarfs (WDs) in wide binaries goes back to the early days of modern astrometry, when in 1844 Friedrich Bessel reported on peculiar changes in the sky position of Sirius and Procyon.After considering several possibilities, Bessel concluded that Sirius and Procyon are not single stars and move under the gravitational force induced by unseen companions of unknown nature.The proximity of Sirius and Procyon to Earth eventually allowed for direct imaging of these faint companions (Bond 1862; Schaeberle 1896), and led to their identification as WDs.
The detection of WDs in binaries that are located at distances considerably greater than that of Sirius and Procyon relies on the indirect method initially proposed by Bessel-detecting the faint WD through its effect on the orbital motion of its visibly brighter, usually main-sequence (MS), primary star.However, the analysis of the binary motion alone is often insufficient to determine the nature of the faint companion.Traditionally, identifying WD companions has depended on observing their contributions to the system's brightness in short-wavelength bands (e.g.Parsons et al. 2016, and references therein).The flux ratio between the WD companions and their bright stars can sometimes be detected in those bands, as the WDs can be much hotter than their companions.
Still, in many cases, the light contribution of the WD is too small ★ These authors contributed equally.† E-mail: sahar.shahaf@weizmann.ac.il ‡ E-mail: naama.hallakoun@weizmann.ac.il to be detected, even in the short-wavelength bands.Therefore, it is likely that many WD companions to MS stars eluded detection or were incorrectly classified (Holberg et al. 2013), particularly if their MS companion is of spectral type earlier than M0.This situation has changed drastically with the release of a catalogue of nearly 170 000 astrometric binaries, many of which are presumed to have unseen WD companions, in the third data release (DR3) of Gaia (Gaia Collaboration, Arenou, et al. 2022).Compact companions in such systems can be identified by ruling out all other possibilities.These include a relatively faint single MS companion or, in case the observed binary is wide enough, a companion that is by itself a short-period binary consisting of two faint MS stars (e.g.Shahaf et al. 2019;Shenar et al. 2022;Janssens et al. 2022;Shahaf et al. 2023;El-Badry et al. 2023;Chakrabarti et al. 2023).Obtaining a large sample of binaries with WD secondaries may refine our understanding of the late stages of stellar evolution of the binary population, the WD initial-to-final mass relation, and the processes governing mass transfer and co-evolution in binary systems (Gratton et al. 2021;Venner et al. 2023;Escorza & De Rosa 2023;Zhang et al. 2023a;Sayeed et al. 2023).
The first paper of this series, Shahaf et al. (2023, Paper I henceforth), identified the astrometric binaries that are most unlikely to have single MS companions.Based on this classification, the work presented here divides those binaries into two sub-groups: those whose companions are likely to be short-period faint MS binaries by themselves (i.e., hierarchical triple systems) and those that probably have WD companions.The division between these two sub-groups is done by identifying the light contribution of the faint companion(s).Instead of searching for the short-wavelength contribution from the WD, we utilise the availability of Gaia's synthetic photometry to identify the long-wavelength contribution of faint MS close-binary companions.We assume that the binaries with no long-wavelength excess are the ones with WD companions.
The structure of this paper is as follows: In Section 2, we outline the process we used to select our sample and provide a review of its characteristics.Section 3 discusses the WD population within our sample and presents our findings.In particular, subsection 3.2 is dedicated to presenting the ultraviolet (UV) emission detected in some of our targets.In Section 4, we discuss the various biases and selection effects in our sample.Finally, in Section 5, we summarise our findings and discuss potential avenues for further research.

SAMPLE SELECTION
To identify astrometric binaries with WD companions, we compiled a sample of systems where the secondary component is less likely to be a single MS star.We then utilised the photometric measurements provided byGaia to identify excess near-infrared emission to distinguish between hierarchical triple systems and MS+WD binaries.The data underlying this analysis are obtained from table 1 of Paper I and the following Gaia DR3 tables: gaia_source, nss_two_body_orbit, binary_masses, and synthetic_photometry_gspc.

Astrometric triage
The astrometric mass-ratio function (AMRF; Shahaf et al. 2019) is defined as where  is the orbital period,  1 is the primary mass,  is the parallax, and  0 is the angular semimajor axis of the photocentric orbit.Shahaf et al. (2019) showed that based on the AMRF value, one can assign each binary with an MS primary to one of three classes: class-I, where the companion is most likely a single MS star; class-II, where the companion cannot be a single MS star, but can still be a short-period binary of two MS stars; and, class-III, where the companion cannot be a single MS star or a close MS binary; therefore the faint astrometric secondary in such systems is probably a compact object.
Paper I used the AMRF approach to analyse a sample of 101 380 astrometric binaries reported in Gaia DR3.The selected binaries have MS primary stars, orbital periods shorter than 1 000 days, and orbital solutions that passed several quality criteria (see therein), estimating the probability of each system to be a class-II or class-III astrometric binary.These probabilities, denoted Pr II and Pr III, are provided in table 1 therein.Paper I used these classification probabilities to compile a short list of 177 systems for which Pr III is close to unity.Assuming that Gaia's orbital parameters are valid, their faint astrometric secondary components are likely to be compact objects.In this work, we take a different approach, excluding class-I systems instead of looking for class-III systems.Namely, we consider astrometric orbits for which the secondary is probably not a single MS star.
The class-I probability, Pr I, was not provided in Paper I but can be calculated using the other two classification probabilities, where  = 10 5 .We include in our sample only systems that satisfy the qualifying conditions outlined in Paper I (see section 3.1, therein) and have Pr I below 10 per cent.This selection yielded 11 190 nonclass-I systems.
Binaries with a primary mass larger than ∼1.2 M ⊙ (within 1), for which WD companions are unlikely to be classified as class-II or -III systems (see figure 2 of Paper I), were excluded from our sample.This selection criterion resulted in 10 080 binaries.We then removed additional binaries for which the Gaia Synthetic Photometry Catalogue (GSPC; Gaia Collaboration et al. 2022b) did not provide estimates in the Johnson-Kron-Cousins bands (JKC; see below).Following this final step, we were left with a sample of 9 786 astrometric binaries, presented in Fig. 1.

Expected colour excess
We use Gaia's spectrophotometric measurements to distinguish between hierarchical triples and MS+WD binaries and resolve this classification ambiguity.Fig. 2 shows the expected colour excess induced on the observed MS primary by an unresolved close-binary companion, using JKC's  and  bands.We define the colour excess as the difference between the triple system's and the primary star's colours.Formally, this is given by Δ ( − ) = − 2.5 log 10 10 −0.4 1 + 10 −0.4 2a + 10 −0.4 2b + 2.5 log 10 10 −0.4 1 + 10 −0.4 2a + 10 −0.4 2b − ( 1 −  1 ) , (3) Colour excess (mag) where the subscript 1 refers to the astrometric primary.The other two subscripts, 2a and 2b, refer to the two components of the astrometric secondary.
The estimates in Fig. 2 were calculated using PARSEC1 isochrones (Bressan et al. 2012;Chen et al. 2014Chen et al. , 2015;;Tang et al. 2014).The close binary is assumed to comprise two equal-mass MS stars since this configuration yields the largest mass-luminosity ratio between the astrometric components.The figure illustrates that per-cent level photometric precision should be sufficient to identify the close-binary companions for systems with primaries less massive than ∼1.2 M ⊙ .

Observed colour excess
Gaia DR3 catalogue provides flux-calibrated low-resolution spectrophotometric measurements in the wavelength range 330-1050 nm.Based on these data, the Gaia Synthetic Photometry Catalogue (GSPC; Gaia Collaboration et al. 2022b) published synthetic photometric measurements with accuracy as low as ∼10 mmag.We compared the location of each source on the Gaia colour-magnitude diagram (CMD) to its theoretically expected position, assuming it is a single star, in search of excess infrared flux.
The system's mass, age, metallicity, and extinction along the line of sight determine the theoretical CMD position.As presented below, we calculated the theoretical position assuming a fixed age of 2 Gyr.Each system's mass, metallicity and extinction are accounted for, along with their estimated uncertainties.This theoretical position is compared with the observations to identify and quantify the infrared excess flux.
We used Zhang et al. (2023b) metallicity estimates.These values were derived in a forward-modelling approach, using Gaia's fluxcalibrated low-resolution spectrophotometric measurements (i.e., Gaia's blue photometre (BP)/red photometre (RP) spectra, or "XP spectra" for short).The data-driven model Zhang et al. (2023b) developed was trained on atmospheric parameters reported by the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST) survey.Using its XP spectrum, the model provides a given target's metallicity, surface gravity, and effective temperature.The effects of binarity were not included in the model and may bias the estimated values of binaries with flux ratios of order unity Zhang et al. (2023b).However, considering the typically extreme flux ratios of the systems discussed in this work, we expect this issue to have a limited impact (but see the discussion in Section 4).
Systems flagged to have uncertain metallicity estimates by Zhang et al. (2023b, their quality_flags ≥ 8) were assigned with [M/H] = 0.00±0.25 dex.This choice is justified by the age-metallicity relation of stars in the solar neighbourhood (Rebassa-Mansergas et al. 2021).Systems with metallicity estimates laying outside of the PAR-SEC isochrones range (−2 ≤ [M/H] ≤ 0.6) were assigned with the corresponding limiting values.
We used publicly available dust maps to account for reddening and extinction.Green et al. (2019) have created Bayestar19: a threedimensional dust map that covers the northern hemisphere down to a declination of −30 deg.We used the dustmaps Python package to obtain the 16 th , 50 th , and 84 th percentiles of Bayestar19 extinction.These values were then multiplied by 0.884 to obtain the corresponding  ( − ) values.2For the southern systems, which were not included in Bayestar19, we used the three-dimensional dust map of Lallement et al. (2019). 3This map does not include error estimates.Considering systems included in both maps, the median difference between Bayestar19 and the Lallement et al. (2019) extinction is 0.02 mag.We therefore used this value as the uncertainty for the  ( − ) values of the southern hemisphere systems.Finally, we converted the  ( − ) values to  ( − ) and   , using the relations  ( − ) = (  −   ) ×  ( − ) and where the coefficients   ,   and   were taken as 3.626, 1.505 and 3.1, respectively (according to Munari & Carraro 1996;Schlafly & Finkbeiner 2011).
We estimated the colour excess using the Python package stam4 (Hallakoun & Maoz 2021): Using PARSEC main-sequence isochrones, we generated a two-dimensional interpolant grid5 representing the relation between the  −  colour index, the absolute -band magnitude, and the metallicity.We limited the maximal mass of the isochrones to 1.17 M ⊙ to avoid tracks that start to evolve off the MS, and assumed a fixed age of 2 Gyr.Stars do not significantly change their mass or appearance during their MS lifetime; therefore, this choice of isochrone age does not severely affect colour excess estimates.Our specific stellar age choice is consistent with the ages estimated for the more massive stars in our sample and justified by the observed evidence for a star formation burst in the Milky Way's disc about 2 Gyr ago (Mor et al. 2019).
The observed colour excess is defined as where ( − ) observed is the observed dereddened  −  colour index, and the ( − ) expected is the theoretically expected  −  value, assuming it is a single MS star.We calculated the theoretically expected values using the interpolant grid described above, with the system's estimated metallicity and extinction-corrected absolute -band magnitude as its independent variables.We conduct a Monte-Carlo experiment to account for the uncertainties in the observed quantities and provide confidence intervals for the estimated colour excess.For each astrometric binary, we generated 10 4 random realisations of its position on the metallicityabsolute magnitude grid.The absolute -band magnitude of the system is not expected to change by more than ∼ 0.1 mag by the presence of a WD companion (compared to a single MS star) for almost all our sample, except for those with primaries less massive than ∼ 0.3 M ⊙ and warm WD companions (see Section 4 and Fig. 11 below).We assumed that the metallicity and the absolute magnitude follow an uncorrelated bi-variate Gaussian distribution.The measurement uncertainties of the observed colour index and absolute magnitude were calculated following Hallakoun & Maoz (2021).For each instance, we calculated Δ ( − ) as described in equation ( 5).This sampled distribution's mean and standard deviation provide our estimates for the colour excess expected value and uncertainty.We provide the colour-excess estimates for all the targets in our sample in Table A1.
The colour excess median uncertainty (first percentile) value is 0.15 (0.09) mag.A systematic deviation of around 0.1 dex to our metallicity estimates will bias the colour excess by ∼0.05 mag.This value is smaller than the typical colour excess uncertainty, suggesting that our analysis is robust to metallicity systematics of this scale (also see Section 4 and Fig. 15).Similarly, our analysis is insensitive to the stellar age of a star as long as it is on the MS.For instance, the colour index of an MS star less massive than ∼1.1 M ⊙ will change by less than 0.06 mag as its age varies from 2 to 5 Gyr.However, inaccurate mass estimates may induce a more significant effect.Suppose the mass estimates have a systematic deviation of ∼5 per cent.Considering the mass range of primaries analysed in this work, such an error can induce a colour excess bias of up to around 0.25 mag.The primary masses used in this study were taken from Gaia's binary_masses table, described in Gaia Collaboration et al. (2022a).These estimates were calculated for Gaia's non-single-star catalogue, using the full orbital solution, the G-band magnitude, and the Lallement et al. (2019) extinction map.Their modelling procedure also assumed the companion is either non-luminous or a single MS star.The colour index, however, was not used in this process.
The primary mass estimation procedure can sometimes result in overestimated primary masses accompanied by significant infrared colour excess.For example, in the case of a triple system, improperly modelled flux contribution from the astrometric secondary may increase the derived primary mass due to excess G-band magnitude.This erroneous mass estimate may manifest as significant infrared colour excess.The expected deviations, on the order of ∼0.25 mag, are similar to the typical colour excess of the RCE sample (see Fig. 3).A possible fingerprint of the primary mass modelling inference in the RCE sample is discussed below.

Partitioning by colour excess
We define Pr(red) as the fraction of realisations where the dereddened observed colour index was larger than the expected colour excess.We have set the threshold above which a system is considered to have infrared excess to be the maximal Pr(red) value for which the measured colour excess is smaller than 10 −2 , Pr(red) ⩾ 56 per cent.
(6) This threshold ensures that we are sensitive to the expected colour excess from all triple systems in our sample (see Fig. 2).We provide the colour-excess probability for all the targets in our sample in Table A1.
Our non-class-I sample includes 6 620 systems with red colour excess (RCE).The remaining 3 166 targets do not exceed the abovementioned threshold and will be referred to as no-colour-excess (NCE) systems.Fig. 3 shows the colour excess distribution of the entire population and the RCE and NCE subsamples.
Fig. 4 presents our sample, with red points indicating RCE systems and grey points representing NCE systems.A black line shows a 2 Gyr isochrone with a metallicity of [M/H] = −0.2, for demonstration (the colour-excess determination was performed individually for each system as described above).Using the same isochrone, we depict the positions of MS binaries, hierarchical triples, and MS+WD binaries.The black dashed line shows the location of binaries with a mass ratio of 0.8; the dashed-dotted line represents triple systems with an equalmass close binary as the faint astrometric companion, assuming a mass ratio of 0.8 of the astrometric components; the dotted lines represent MS+WD binaries at various WD cooling ages, derived assuming a hydrogen-dominated WD with a mass of 0.6 M ⊙ (using the synthetic colours of Bergeron et al. 1995;Holberg & Bergeron 2006;Kowalski & Saumon 2006;Tremblay et al. 2011;Blouin et al. 2018;Bédard et al. 2020). 6he colour excess estimate was obtained using synthetic photometry in the JKC system.Other photometric systems, such as the SDSS system or Gaia's  BP and  RP bands, can also be used.We found that the JKC's  −  colours are more sensitive to the colour excess than the broader Gaia  BP and  RP bands, and are available for almost all targets in our sample, as opposed to the SDSS synthetic photometry in the  band.

A CENSUS OF WHITE DWARFS
In the previous section, we identified a sample of astrometric binaries where the companion is unlikely to be a single MS star.Within this sample, we used the  −  colour index to identify targets exhibiting red colour excess, allowing us to identify the contribution of light from close binary companions in hierarchical triple systems.Thus, we now have a sample (referred to as NCE) in which the secondaries are unlikely to be (i) single MS stars, or (ii) close binaries of two MS stars.
The following section shows that the NCE sample is consistent with a predominantly MS+WD population.

Sample characteristics
Fig. 5 shows the relationship between the AMRF and the mass of the primary star for the NCE and RCE populations.The dotted line separates between class-I and -II systems, and the dashed line separates class-II from -III.The grey stripes represent the expected position of non-luminous companions in the mass range of 0.45 − 0.75 M ⊙ .Notably, the two populations are differently distributed: The NCE sample is concentrated around the grey stripe, which corresponds to the expected mass range for WD companions; in contrast, the RCE sample is more spread out and populates a larger area on the  1 − A plane.The distribution of the NCE sample indicates that the population is mostly comprised of MS+WD binaries.However, the RCE sample is not a pure sample of hierarchical triples, as suggested by the RCE class-III binaries found on the grey stripe.These are likely MS+WD binaries that were found to have an apparent colour excess (see below).
Fig. 6 compares the distances, primary masses, -band extinction, and Galactic latitudes of both populations.This figure shows that RCE systems are located at greater distances and closer to the Galactic plane than the NCE sample.The extinction histogram indicates that this may be due to our colour classification, as systems closer to the Galactic plane suffer from greater extinction and appear redder, suggesting that some of our extinction or metallicity estimates may be inaccurate.We further discuss this possibility in the following.An alternative explanation is that triple systems, which may be more luminous than MS+WD binaries, are observed at greater distances and show larger extinction values.The RCE primary-mass histogram in the top right panel of Fig. 6 shows an excess of ∼1 M ⊙ targets, accompanied by a deficiency of ∼0.75 M ⊙ primaries (which is also apparent in the right panel Fig. 5).These features do not appear in the NCE primary mass estimates.
We speculate these artefacts stem from the primary mass estimation procedure, as they suggest that infrared excess is associated with its products.As mentioned above, the masses were derived using the absolute magnitudes of the targets without considering their colour index.Additionally, while the process did account for possible flux contribution from a single MS companion in a binary system, hierarchical triple systems were not modelled.It is, therefore, possible that triple systems were analysed, assuming their companion is nonluminous, in a process that affected their primary mass estimates, leaving a significant infrared colour excess signature.
To test whether the differences between the populations are caused by vulnerabilities in the orbital fitting scheme, we checked for systematic differences in the number of measurements (visibility periods used; VPU), the goodness of fit (GoF), and the significance of the angular semi-major axis ( 0 /Δ 0 ) and parallax (/Δ).The corresponding histograms are plotted in Fig. 7.While the two populations follow similar trends in these four parameters, the NCE population seems to have slightly better significance and GoF values.This is expected, considering the larger distances and lower Galactic latitudes of the RCE sample.
Fig. 8 presents the eccentricity versus the orbital period of the stars in the sample.Most of the eccentricities of the RCE sample (shown in the right panel) are below ∼0.8.This is probably an observational bias caused by the significant amount of time highly eccentric systems spend close to apastron, where the motion is not easily detected.Furthermore, most of the RCE binaries are found below the dotted line that indicates orbits that reach separations of ∼100 R ⊙ at periastron passage, assuming a total mass of 1 M ⊙ .On the other hand, the eccentricities of the NCE sample (shown in the left panel) are mostly below ∼0.3, as can be seen in the left panels of Fig. 9.
The secondary mass distribution of the NCE targets, shown in Fig. 9, peaks sharply around ∼0.6 M ⊙ , consistent with them being MS+WD binaries.In contrast, the RCE sample displays a broader secondary mass distribution peaking at ∼0.7 M ⊙ .However, even if we assume that the RCE sample is purely comprised of triple systems, inferring the secondary mass or mass-ratio distribution of hierarchical triple systems based on it is not straightforward.The secondary mass estimates assume negligible light contribution from the astrometric secondary (see table 1 in Paper I).This assumption is invalid for hierarchical triple systems, making these estimates merely lower limits of the secondary mass (Shahaf et al. 2019).
Fig. 9 also presents histograms of the orbital periods and the cosine of the inclination angle.Compared to the NCE population, the occurrence of systems with orbital periods below one year is smaller in the RCE population.One simple explanation is that RCE systems are found at greater distances, where relatively close separations are harder to detect.However, this difference may also reflect the interaction between the present-day primary MS and the WD progenitor when it was a red giant or, alternatively, the stability of hierarchical triple systems.In terms of the inclination angle, the two populations exhibit similar distributions.

UV excess
Newly formed WDs can significantly contribute to the binary system's brightness in the visible bands (e.g., Fig. 4).This contribution becomes more pronounced in the UV, where MS primaries below ∼1.2 M ⊙ are relatively faint.To identify these WD contributions, we searched the Galaxy Evolution Explorer (GALEX; Morrissey et al. 2007) database for UV counterparts to the NCE binaries in our sample.A similar analysis of the class-III binaries from Paper I was recently done by Ganguly et al. (2023).
We cross-matched the NCE targets in our sample with the GALEX all-sky imaging survey (AIS; Bianchi et al. 2011Bianchi et al. , 2017)).The Gaia coordinates were propagated back to the coordinates at the AIS midpoint using the measured proper motion of each system.The fiducial search cone radius around each coordinate was 3 arcsec (Bianchi & Shiao 2020).This cone was enlarged to account for the proper motion of each system during the GALEX mission duration.
We have identified 3 673 matches between GALEX and Gaia targets, out of which 1 436 belong to the NCE population.Table A2 provides a list of the matched sources, including the Gaia and GALEX identifiers, apparent UV magnitude, and the angular separation between the matched sources.The median separation between the two positions of the matched sources is ∼0.6 arcsec.Only 7 sources are separated by 3 arcsec or more, with a maximal separation of ∼3.6 arcsec.There are 9 targets for which we identified a second source within the search cone; we chose the closest one in these cases.
We use a simple single-valued relation to estimate the extinction in the GALEX near-UV (NUV) band. 7The relation is given by where  NUV and  G represent the extinction in the NUV and G bands, respectively.Their corresponding extinction coefficients,  NUV and  G , were obtained from Zhang & Yuan (2023).Fig. 10 presents the NUV absolute magnitudes of the NCE sample.The solid black line depicts the corresponding MS mass-luminosity relation obtained using our default isochrone (see above).The PAR-SEC NUV magnitudes are provided in the Vega photometric system; therefore, the derived relation was shifted by 1.699 mag to match the AB system of the AIS (Bianchi 2011).To illustrate the effect of extinction on the UV emission from the system, we added dashed and dotted black lines representing an extinction of 0.75 and 1.5 mag, respectively.These values roughly correspond to the 80 th and 95 th percentiles of the NCE sample.
Most points in Fig. 10 fall between the solid and dotted black lines, as expected.To illustrate the expected contribution of the WD, we plotted as dotted blue lines the expected absolute NUV magnitudes of 0.6 M ⊙ WDs with hydrogen-dominated photospheres at different cooling ages.For this purpose we used the tables 6 of WD synthetic colours by Bergeron et al. (1995); Holberg & Bergeron (2006); Kowalski & Saumon (2006); Tremblay et al. (2011); Blouin et al. (2018); and Bédard et al. (2020).The diagram suggests that WDs with cooling ages below ∼0.5 Gyr should be easily detectable if their primary MS companions are of spectral types later than K.  Around 8 per cent of the NCE sample, 117 targets, show significant excess UV emission.As with the infrared colour excesses, we somewhat arbitrarily defined excess UV as a system brighter than expected MS single-star curve by more than 1 (see Table A2).We accounted for the NUV magnitude, parallax, extinction and primary mass uncertainties.The RCE sample, for comparison, contains 118 binaries (∼5 per cent) showing NUV excess detected by the same criterion.However, nearly 40 of these systems have primaries with masses of ∼1 M ⊙ , which is unlikely, given the young WD cool-ing ages implied.Considering the RCE primary mass distribution in Figure 6, it is likely that the expected NUV value of some of these systems is erroneous.The fraction of RCE systems with NUV excess is probably closer to ∼3.5 per cent.
As discussed in the following section, the colour-based sample segmentation aims to ensure that the NCE sample is less contaminated by triple systems, not vice versa.Therefore, we expect the RCE sample to contain triple systems and MS+WD binaries.Fig. 5 shows that the NUV-bright systems are localised around the WD grey stripe on the A− 1 diagram, both for the NCE and RCE samples, as expected.The higher occurrence of NUV-bright binaries in the NCE sample and the consistency of their secondary masses with the typical mass of a WD corroborates our classification scheme.The ratio between the two samples suggests that the contamination of the RCE sample can reach, at worst, about 60 per cent.Still, a more detailed analysis of the spectral energy distribution of each system is required (e.g., Ganguly et al. 2023).

BIASES AND SELECTION EFFECTS
The first obvious selection effect in our NCE sample stems from the pre-selection of primary MS stars since low-mass MS stars will appear as photometric primaries in the Gaia  band only if their WD companion is cool enough.However, this will affect only the smallest MS primaries, as demonstrated in Fig. 11, which shows the expected photometric excess in the Gaia  band of an MS+WD binary compared to the expected emission from the primary MS alone, as a function of primary mass and WD effective temperature.
From this figure it is clear that only binaries with a ∼ 0.3 M ⊙ MS star and a WD warmer than ∼ 16 000 K, or a ∼ 0.2 M ⊙ MS star and a WD warmer than ∼ 12 000 K, would be a priori excluded from our sample, since the WD would dominate the Gaia  band.This explains the sparsity of systems with smaller primary masses (see Fig. 5).The expected photometric excess was calculated assuming black-body spectra for both components.The radius of the MS component was estimated using a 2 Gyr, [M/H] = 0, PARSEC isochrone, while the WD radius was estimated using the WD_models8 Python package with the models of Bédard et al. (2020), assuming a 0.6 M ⊙ hydrogendominated DA WD.
Our NCE sample, presumably primarily consisting of systems with WD companions, suggests that many MS+WD binaries exhibit  a modest eccentricity of ∼0.15.This eccentricity might be linked to processes during the final stages of stellar or binary evolution (see Section 5.1.1 below).Thus, it is noteworthy that the eccentricity distribution of our sample is significantly biased.Fig. 12 presents the orbital eccentricity versus the system's distance, indicating a tendency towards higher eccentricities with increasing distances.One explanation for this bias could be that the eccentricity is derived more accurately at closer distances; if the inference procedure is sensitive to correlated noise, this situation can result in biased eccentricity estimation (Lucy & Sweeney 1971;Bashi et al. 2022; also see, for example, a similar discussion by Hara et al. 2019 for the case of radial-velocity modelling of exoplanetary orbits).However, a comprehensive analysis to determine the origin of this bias requires a detailed study of Gaia's selection function, which falls outside the scope of this study.
Another selection bias is related to the masses of the primary stars.The triage classification limits vary as a function of the primary star's mass (see Fig. 5, for example).As a result, the WD mass distribution is biased since the minimal WD mass grows with the mass of the primary star, as Fig. 13 illustrates.One immediate consequence of this selection effect is a bias towards high WD masses and high mass ratios in the NCE sample.However, other biases may also stem from this selection effect.For example, since all the primary stars in our sample are on the MS, the correlation between the secondary and primary masses can be translated into a correlation between the WD's mass and the binary system's absolute magnitude.As a result, the maximal distance-and hence the orbital separation and eccentricity-of the detected systems might also be affected.
Finally, the assumption that the light contribution from the WD is negligible compared to that of the MS primary star is not necessarily accurate.Despite their compact nature, their contribution to the overall luminosity of the binary system may be non-negligible.If unaccounted for, this contribution may bias the WD mass estimates.The relative intensity of the WD compared to its MS companion depends on several factors, including the mass of the MS primary, the mass of the WD itself, and its cooling age.Since WDs cool with time over billions of years, younger WDs are brighter and bluer.In Fig. 14, we plot the minimal cooling age below which WDs contribute more than 5 per cent of the system's light in Gaia's bandpass.WDs younger than this critical age can introduce bias more significant than ∼10 per cent to the mass estimation process.For systems with ≳ 0.5 M ⊙ MS primaries and WDs with cooling ages older than ∼0.5 Gyr, the contribution of light from the WD cannot distort their mass estimates significantly.Each line represents a system with a different primary MS mass.
It is important to mention that the RCE sample does include some MS+WD binaries.For example, some of the RCE binaries are located above the class-III limit and within the region of the WD mass stripe of Fig. 5. Another indication of the WD contamination in the RCE sample is the NUV-bright systems.The position of these systems on the diagram is generally consistent with the WD stripe.Hence, these are likely misclassified MS+WD binaries identified as having excess infrared emission.This misclassification is probably the result of erroneous metallicity estimates for systems with primary masses smaller than ∼ 0.6 M ⊙ (see Fig. 15) that lead to an unreliable colourexcess assignment.This is consistent with fig.18 of Zhang et al. (2023b), demonstrating the low fraction of reliable stellar parameter estimates for M dwarfs that were not included in the LAMOST training set used in that study.Luckily, almost the entire expected range of M-dwarf primaries with WD secondaries resides within class III on the AMRF plot (Fig. 5).Thus, all class-III systems with primary masses ≲ 0.6 M ⊙ are very likely to be MS+WD binaries, regardless of their assigned colour excess.Since there are very few class-II systems with M-dwarf primaries, the impact of this issue on our final MS+WD catalogue is small and can affect only systems with extremely low-mass WD companions.Better estimates of Mdwarf metallicities will enable more accurate classifications of these systems in the future.
Class-III systems with primaries more massive than ∼ 0.6 M ⊙ that are included in the RCE sample are more likely to have erroneous primary mass estimates (especially those with masses around ∼ 1 M ⊙ , see discussion above).Alternatively, the assumption of 2 Gyr-old isochrones used for the colour-excess estimation (see Section 2.2) might be inaccurate for the systems with the most massive primaries, which could be either younger or already in the process of leaving the MS.Finally, we mention that the RCE sample tends towards lower Galactic latitudes.In more crowded fields, excess flux from foreground or background sources could potentially impact the mass estimates, orbital fitting, or colour excess measurements.It is possible that crowding contributed to the artefacts detected in the primary mass histogram and affected the GoF distributions of the orbital solutions (see Fig. 6 and 7, respectively).

Astrophysical implications
The astrometric binaries presented here provide a uniquely large sample of nearly 3 200 probable MS+WD binaries with orbital separations of ∼1 au.These systems populate a region of parameter space largely unexplored by other observational techniques-only a few MS+WD binary systems with orbital periods of hundreds of days are currently known (see Anguiano et al. 2022;Parsons et al. 2023, and references therein), enlarging the sample size by about two orders of magnitude.This sample has already yielded some insights into the WD population and binary evolution (Hallakoun et al. 2023).

Binary evolution
On the one hand, the present orbits of many of the binaries are too small to allow for uninterrupted evolution of the WD progenitors through their red-giant phase; the eccentricity-period diagram displays some binaries with large eccentricities, such that the binary separation at the periastron is smaller than, say, 100  ⊙ .On the other hand, their separation is likely too large to assume that they went through a common-envelope phase, as post-common envelope binaries with a WD and M-dwarf components have been found to exist only at the shortest periods (less than about two weeks, peaking at ∼8 hr; Nebot Gómez-Morán et al. ( 2011 2021), but see also Yamaguchi et al. (2024) for some systems with periods of up to seven weeks).Therefore, some of these binaries, after the WD nature of the unseen companion is confirmed, might necessitate special binary evolutionary tracks (e.g.Perets & Fabrycky 2009).
Stable mass transfer is a known mechanism to form long-period binaries with WD components.In fact, theory predicts a relation between the period and the WD mass for systems that have evolved this way (Rappaport et al. 1995;Chen et al. 2013).For WD masses between 0.45 − 0.5 M ⊙ and donor stars on the Red Giant Branch (RGB), Chen et al. (2013) predicts periods between 600 − 935 days and 275 − 428 days at metallicities of  = 0.01 and  = 0.0001, which are in good agreement with the periods of our sample.However, most of our WD masses are above 0.55 M ⊙ , so their direct progenitors cannot be RGB stars but must have already been on the Asymptotic Giant Branch (AGB) at the onset of the mass transfer.While mass transfer with an AGB donor is typically expected to lead to a common-envelope phase (e.g.Hjellming & Webbink 1987;Hurley et al. 2002, but see Mohamed & Podsiadlowski 2007), our systems imply stable mass transfer is possible, likely aided by mass loss through non-conservative mass transfer (Soberman et al. 1997).In addition, the eccentricities of our systems cannot easily be explained by stable mass transfer alone since, in classical binary evolution theory, tides would circularise the orbit before the onset of the mass transfer.Possible eccentricity pumping mechanisms are enhanced wind mass loss at periastron (Van Winckel et al. 1995;Bonačić Marinović et al. 2008), Roche-lobe overflow at periastron (Soker 2000;Vos et al. 2015), formation of a circumbinary disc (Waelkens et al. 1996;Dermine et al. 2013;Vos et al. 2015), and WD natal kicks (Izzard et al. 2010).
If some of these binaries underwent mild mass-transfer interaction, we might find evidence for heavy s-process element pollution in the luminous MS star atmosphere.Very few known barium stars are in our sample, and more generally in the Gaia non-single star catalogue.Interestingly, the period-eccentricity distribution in the left panel of Fig. 8 resembles that seen in barium stars at a similar range of orbital periods (e.g., Jorissen et al. 2016Jorissen et al. , 2019;;Escorza & De Rosa 2023).The orbital elements and amount of enrichment may uncover the properties of the WD progenitors and the underlying mass-transfer process.Sayeed et al. (2023), for example, recently proposed binary evolution as a possible formation pathway of fast-rotating chemically enriched red giants.Studying the element abundances of the systems in our catalogue may uncover the properties and significance of these mechanisms in a different range of separations and binary evolution stages.
Fig. 16 compares the period-mass distribution of the observed NCE sample with the predicted distribution of various binary population synthesis (BPS) models.The synthetic populations shown here were generated using the SeBa BPS code (Portegies Zwart & Verbunt 1996;Nelemans et al. 2001;Toonen et al. 2012).The models differ in their treatment of the common-envelope phase.The process is parameterised using scalar factors governing the efficiency with which energy or angular momentum transfer unbind the envelope.These parameters are often denoted "" (Paczynski 1976;Webbink 1984;Livio & Soker 1988) and "" (Nelemans et al. 2000;Nelemans & Tout 2005;van der Sluys et al. 2006), respectively.The expulsion of the envelope shrinks the orbit of the binary in a process that can lead to the formation of a close double-WD system.These systems presumably experienced at least two phases of mass transfer, of which at least one occurred during a common-envelope phase.
The BPS samples used here rely on two model families,  and  (Toonen et al. 2012), defined based on the formalism used for each common-envelope phase.The  prescription, which is the standard model in most BPS codes (i.e. = 2, where  is a parameter that depends on the structure of the donor star), has two additional variants: The 2 model, capable of producing a more significant orbital shrinkage during the common-envelope phase (i.e. = 0.25), was developed to match known populations of close MS+WD binaries (Toonen & Nelemans 2013;Zorotovic et al. 2010;Camacho et al. 2014;Scherbak & Fuller 2023); The 2 model with isotropic reemission, in which mass is allowed to leave the system with the  specific angular momentum of the accretor (instead of 2.5 times that of the orbit), is also capable of producing wider orbits.The resulting period versus WD mass distributions of all four prescriptions are presented in Fig. 16, derived during the MS+WD evolutionary stage of the synthetic population.Note that in the context of this study only the treatment of the first mass transfer phase resulting in an MS+WD system is relevant, and the two-phase model names are kept to ease the comparison with previous BPS studies.
The observed NCE population is plotted as red dots in Fig. 16.Although this population is affected by the selection criteria and observational biases, it is useful to consider whether current BPS codes are at all capable of producing it.We note that we only consider simulated MS+WD systems where the mass of the MS star is smaller than 1.2 M ⊙ , to reflect the selection in our observed sample.The  model and its variants have gaps around orbital periods of a few hundred days, where our observed NCE sample resides.Therefore, while some  variants might further improve by adjusting their parameters, this prescription seems less suitable for describing the observed NCE population.In contrast, the  model, which employs the  formalism for the first mass transfer event, succeeds in generating some MS+WD systems in part of the observed parameter range.This model was developed in order to match known populations of close double-WD binaries (Nelemans et al. 2000;Nelemans & Tout 2005;van der Sluys et al. 2006).We note, however, that all these models fail to reproduce MS+WD binaries with eccentricities larger than zero, in contrast with our observed sample.
If indeed the NCE sample represents a homogeneous MS+WD population, Fig. 16 suggests that some BPS prescriptions do not fully describe the binary population at the ∼1 au separation range.This unprecedentedly large sample of MS+WD systems with intermediate orbital separations should help constrain and adjust these models to better fit the observations.

WD initial to final mass relation
Fig. 10 suggests we can estimate the cooling ages of many WDs in our sample based on their UV excess (but see the caveats discussion below).Since the WD masses were dynamically measured, estimating the binary's age would allow us to derive an empirical initial-to-final mass relation (IFMR) for these WDs.In most cases, these relations are obtained for single WDs in clusters (e.g.Cummings et al. 2018).
This sample allows for empirical estimation of the IFMR in binary systems with orbital separations potentially small enough to induce interaction during the late stages of the WD progenitor's evolution.However, determining the age of MS stars presents a challenge.Techniques like Gyrochronology (Bouma et al. 2023) or abundance analysis may be helpful, yet their reliability, considering the configuration and evolutionary history of the binaries in our sample, is unclear (Silva-Beyer et al. 2023;Gruner et al. 2023).An alternative approach would be associating MS+WD binaries with open clusters.Even though the numbers are relatively small, these systems provide a reliable age estimate for the WD and the binary as a whole.

Caveats and limitations
As in the first paper, we rely heavily on the orbital solution's validity, the astrometric primary's mass estimate, and its classification as an MS star, as provided by Gaia.Spurious astrometric measurements may lead to false detection of binary motion (see appendix E of El-Badry et al. 2023).Inaccurate estimates of the primary might affect the estimates regarding the nature of the faint companion, even if the orbit is valid (e.g., El-Badry et al. 2024).However, the follow-up campaign for the relatively small class-III sample in Paper I required dozens of observations over many months.Scaling up an initiated spectroscopic campaign to dynamically validate the orbits of thousands of systems might be impractical.Alternative approaches, such as non-dynamical validation via detailed modelling of the spectral energy distribution or archival velocity measurements, may be useful.
The infrared-excess classification scheme aims to identify hierarchal triples, not MS+WD binaries.As a result, the RCE sample is not purely comprised of triple systems.The classification probabilities and estimated B-I colour excess values are provided in Table A1 to enable other studies of this population that may require different confidence levels.Furthermore, as discussed in Section 4, the samples presented in this work are biased; astrophysical conclusions drawn based on them must take these biases into account.
Astrometric binaries with an MS primary more massive than ∼1.2M ⊙ will appear as class-I binaries even if their companion is a WD.This is not a trivial limitation, as detecting A-type stars with WD companions is beyond the reach of the triage method.A complementary approach, using interferometric measurements, was recently presented in a series of publications by Waisberg et al. (2023m,n).Their analysis enabled the detection of companions to A-type MS stars at orbital periods on the order of a few tens of years which, in some cases, were identified as WD candidates or close binaries (also see Waisberg et al. 2023a,b,c,d,e,f,g,h,i,j,k,l).
An astrophysical scenario not considered in this work is a triple system in which the secondary companion is a binary comprised of a close binary of a WD and a low-mass MS star.While the occurrence and evolutionary path leading to such systems in this separation range is unclear, empirical constraints on their occurrence might be of interest (Shariat et al. 2023).In the context of this work, such a system could, in principle, demonstrate excess emission both in short and long wavelengths.
We also note that UV excess by itself is not a definitive proof of a WD companion, as a low-mass chromospherically active companion can produce UV excess (e.g., 2MASS J06281844-7621467; see Lagos et al. 2022b).We emphasise that the classification scheme in this work relied on the astrometric orbit.UV emission, if detected, was used to corroborate the astrometric conclusion.Nevertheless, the possibility of activity-induced UV emission should be considered when studying individual systems from this sample, for example, by examining their spectral properties, persistence, and variability.This is particularly important when inferring the IFMR, as the WD cooling ages may be biased by activity-induced excess UV emission.

Prospects for future work
The samples presented in this work offer an opportunity to probe the biases and selection effects in the Gaia catalogue.The NCE sample represents a homogeneous population of binary systems: The light contribution of both components in these binaries can be reliably constrained, their orbits are expected to be fairly circularised, and they are found at various Galactic latitudes.One can, therefore, use it to derive Gaia's selection function for binary orbits (e.g., Rix et al. 2021;Cantat-Gaudin et al. 2023;Castro-Ginard et al. 2023).
For some science cases, validation of the orbits and classification validity may be required.Medium-resolution broadband spectroscopy offers a straightforward path to determine the nature of systems in our NCE sample.The spectra will allow us to characterise the primary star, leading to better constraints on its mass and metallicity.Excess emission at wavelengths ≤ 450 nm (≥ 750 nm) might reveal the presence of a WD (late-type M-dwarf) companion.Fig. 17 shows the signal-to-noise ratio (SNR) required to detect a WD at a given effective temperature (i.e.cooling age) versus wavelength assuming an MS primary of 1 M ⊙ .The component radii were estimated similarly to those of Fig. 11 (see Section 4).
As the median magnitude of a target in our sample is relatively bright ( = 14.33) and the measurements are not time-critical, a flexible program for telescopes operating in queue mode will be an effective follow-up strategy.The upcoming Son of X-Shooter (SoXS) spectrograph, expected to see first light in early 2024 on the European Southern Observatory (ESO) 3.6 m New Technology Telescope (NTT), is our case study for such a program.The spectrograph offers a resolution of ∼4 000 from 0.35 to 2.1μm in a single exposure (Schipani et al. 2020;Rubin et al. 2020).Using the SoXS Exposure Time Calculator,9 we find that a ∼150 sec exposure is sufficient to achieve an SNR of 10 at the short wavelength range of the spectrograph.At ∼1 μm a 25 sec exposure is required to reach the same SNR, similar to the exposure time required at ∼0.4 μm.
Some studies, such as determining the WD IFMR, may require additional spaceborne UV measurements.For example, detailed studies of individual systems containing young WDs can utilise the Cosmic Origins Spectrograph (COS; Osterman et al. 2011) or the Advanced Camera for Surveys (ACS) on board the Hubble Space Telescope (HST) to determine their effective temperature, surface gravity and composition.On the other hand, the Wide-Field Camera 3 (WFC3; Teplitz et al. 2013) on the HST can be used to determine open clusters and WD cooling ages efficiently.
The samples in this work can provide insights into binary and stellar evolution and identify biases in the astrometric binary population.With the forthcoming fourth data release of Gaia, our analysis can be extended to longer orbital periods.This will clarify the processes at the onset of mass transfer and their impact on binary evolution and the Galatic WD population.

Figure 1 .
Figure 1.The astrometric mass-ratio function (AMRF) versus the primary mass of the 9 786 selected non-class-I systems.The dotted line separates between class-I and class-II regions, and the dashed line separates between class-II and -III.

Figure 2 .Figure 3 .
Figure2.The expected  −  colour excess induced by an unresolved binary system versus the mass of the astrometric primary star.Each line in the diagram corresponds to a different mass ratio, , between the astrometric primary and secondary.The close binary is assumed to consist of two equalmass MS stars.

Figure 4 .
Figure 4.A colour-magnitude diagram for the 9 786 objects in our sample, not corrected for extinction.The red-colour excess (RCE) and no-colour excess (NCE) populations are plotted as red and grey points, respectively.The number of objects in each population is given in the figure legend.A 2 Gyr isochrone of [M/H] = −0.2 is plotted as a solid black line.The dashed and dash-dotted black lines illustrate the position of binary and hierarchical triple systems, respectively, and the dotted lines represent the position of MS+WD binaries at different WD cooling ages (see text).

Figure 5 .
Figure 5.The AMRF versus the mass of the primary star, presented separately for the NCE (left panel) and the RCE (right panel) populations.The dotted line separates between class-I and class-II systems, and the dashed line separates class-II from -III.A grey strip highlights the expected location of systems with non-luminous companions in the mass range of 0.45−0.75M ⊙ , where typical MS+WD binaries are expected to reside.The figure shows that, as a population, NCE systems tend to populate this area in parameter space, as opposed to the RCE sample.The white points depict systems in which excess emission in the near-UV band of GALEX was detected, as discussed in Section 3.2, which we ascribe to the contribution of the WD companion.See Section 4 for a discussion of MS+WD contamination of the RCE sample.

Figure 6 .Figure 7 .
Figure 6.Distance, primary mass, -band extinction, and Galactic latitude histograms (clockwise from top left) for the RCE sample (red line) and the NCE sample (grey bars).The median and the 36 − 84 percentile range for each distribution are written on the top left corner of each panel, where the upper and lower rows describe the NCE and RCE samples, respectively.

Figure 8 .
Figure8.A period-eccentricity diagram for the NCE and RCE populations is plotted on the left and right panels, respectively.The gap corresponds to orbital periods of about one year.The dotted line represents an upper envelope that corresponds to periastron distance of ∼100 R ⊙ , assuming a total mass of 1 M ⊙ .The white points depict systems in which excess emission was detected (see Figure5and Section 3.2).

Figure 9 .
Figure 9. Histograms of the eccentricity, orbital period, cosine of the inclination angle, and secondary mass (clockwise from top left).The colour coding is the same as in Fig. 6.The two left panels, showing the eccentricity and secondary mass, also show the median and 36 − 84 percentile range for the NCE (top line) and RCE (bottom line) populations.

Figure 10 .
Figure 10.The GALEX NUV absolute magnitude versus the primary star's mass.The grey circles represent the NCE binaries we have identified near-UV counterparts in GALEX (see text).The solid, dashed and dotted black lines represent the expected MS relation, based on a PARSEC isochrone of 2 Gyr, for extinction values of 0, 0.75 and 1.5 mag, respectively.The horizontal lines show the expected NUV absolute magnitude of a 0.6 M ⊙ WD at various cooling ages.

Figure 11 .
Figure 11.The expected photometric excess in the Gaia  band of an MS+WD binary compared to the expected emission from the primary MS alone.This calculation assumes black-body spectra for both components and a 0.6 M ⊙ hydrogen-dominated DA WD.

Figure 12 .
Figure 12.The orbital eccentricity versus the parallax distance for the systems in the NCE sample.

Figure 13 .Figure 14 .
Figure 13.The secondary versus primary masses of the systems in the NCE sample.Systems with secondary masses larger than 0.8 M ⊙ are not included in this diagram.

Figure 15 .
Figure 15.The metallicity error as a function of primary mass of all non-class-I systems.Reliable [M/H] estimates (quality_flags < 8 in Zhang et al. 2023b) are shown in the greyscale two-dimensional histogram.Unreliable [M/H] estimates (quality_flags ≥ 8) are plotted as violet circles.These systems were assigned a metallicity of 0 ± 0.25 dex in our analysis (see Section 2.2).

Figure 16 .
Figure 16.A comparison to a synthetic population of MS+WD binaries with  MS ≤ 1.2 M ⊙ .Each panel shows the predicted period vs. WD mass distribution (coloured in logarithmic scale) assuming a different SeBa model (see text for details; clockwise from top left): standard   model, standard   model with severe shrinkage in the common-envelope phase, the same as the previous but with isotropic re-emission, and the standard   model.Note that in the context of this study, only the treatment of the first mass transfer phase is relevant, and the two-phase model names are kept to ease the comparison with previous BPS studies.The Gaia NCE sample is plotted as red dots.

Figure 17 .
Figure17.Minimal SNR required to detect a 0.6 M ⊙ WD companion to a 1 M ⊙ MS star.Each line corresponds to a different WD effective temperature (or cooling age).This calculation assumes black-body spectra for both components.The calculation is performed at a spectral resolution of 1Å.