A population of neutron star candidates in wide orbits from Gaia astrometry

We report discovery and spectroscopic follow-up of 21 astrometric binaries containing solar-type stars and dark companions with masses near 1.4 $M_{\odot}$. The simplest interpretation is that the companions are dormant neutron stars (NSs), though ultramassive white dwarfs (WDs) and tight WD+WD binaries cannot be fully excluded. We selected targets from Gaia DR3 astrometric binary solutions in which the luminous star is on the main sequence and the dynamically-implied mass of the unseen companion is (a) more than $1.25\,M_{\odot}$ and (b) too high to be any non-degenerate star or close binary. We obtained multi-epoch radial velocities (RVs) over a period of 700 days, spanning a majority of the orbits' dynamic range in RV. The RVs broadly validate the astrometric solutions and significantly tighten constraints on companion masses. Several systems have companion masses that are unambiguously above the Chandrasekhar limit, while the rest have masses between 1.25 and 1.4 $M_{\odot}$. The orbits are significantly more eccentric at fixed period than those of typical WD + MS binaries, perhaps due to natal kicks. Metal-poor stars are overrepresented in the sample: 3 out of 21 objects (14%) have [Fe/H]$\sim-1.5$ and are on halo orbits, compared to $\sim$0.5% of the parent Gaia binary sample. The metal-poor stars are all strongly enhanced in lithium. The formation history of these objects is puzzling: it is unclear both how the binaries escaped a merger or dramatic orbital shrinkage when the NS progenitors were red supergiants, and how they remained bound when the NSs formed. Gaia has now discovered 3 black holes (BHs) in astrometric binaries with masses above 9 $M_{\odot}$, and 21 NSs with masses near $1.4\,M_{\odot}$. The lack of intermediate-mass objects in this sample is striking, supporting the existence of a BH/NS mass bimodality over 4 orders of magnitude in orbital period.

1. INTRODUCTION Most stars with initial masses ≳ 8 M ⊙ leave behind neutron stars (NSs) when they die.Several thousand NSs are known in the Milky Way, a large majority of which are radio pulsars.Most (> 99%; Lorimer 2008) young pulsars are isolated.Yet, a large majority of the massive stars from which NSs form are in binaries, triples, and higher-order multiples (e.g.Sana et al. 2012;Moe & Di Stefano 2017).The apparent mismatch in multiplicity properties of NSs and their progenitors hints that most massive binaries are destroyed during or prior to the formation of a NS.This destruction can come as a result of binary interaction leading to a stellar merger, or due to the binary becoming unbound during a supernova (SN), which can impart a kick of order 250 km s −1 on the newborn NS (e.g.Hobbs et al. 2005;Faucher-Giguère & Kaspi 2006).
All known companions to young pulsars are massive OB stars (Kaspi et al. 1996;Bassa et al. 2011;Shannon Corresponding author: kelbadry@caltech.eduet al. 2014;Lyne et al. 2015;Andersen et al. 2023;van der Wateren et al. 2023).The same is true for all detached NS + main sequence (MS) binaries detected in X-rays (e.g.Reig 2011).This likely reflects the fact that binaries containing a massive star are more likely to survive mass transfer and a SN when the companion is also a massive star.
However, NSs with low-mass stellar companions do exist, and in fact make up the majority of all known binary NSs.The known systems are all currently accreting from a binary companion ("low-mass X-ray binaries"; LMXBs) or recycled, meaning that past accretion from a companion spun up the pulsar and buried its magnetic field, slowing subsequent spin-down (e.g.Bhattacharya & van den Heuvel 1991).LMXBs and recycled pulsars are over-represented in observed samples because they can have long lifetimes (up to and exceeding the age of the Universe) compared to normal young pulsars, which are only detectable for ∼ 10 7 yrs.The companions to most recycled pulsars are white dwarfs (WDs), low-mass stars, and brown dwarfs that have transferred mass to the NS via stable Roche lobe overflow.Models predict that the initial masses of typical companions were at most (1 − 2) M ⊙ (Podsiadlowski et al. 2002).
Although they have not yet been unambiguously detected, there is little doubt that non-interacting binaries containing a low-mass MS star and a NS exist: such systems are the progenitors of LMXBs and millisecond pulsars.A few candidates for such objects have been identified, including (a) young pulsars with roughly solarmass companions and eccentric orbits (PSR B1820-11 and PSR J1954+2529 Phinney & Verbunt 1991;Parent et al. 2022), and (b) MS stars with unseen companions that may be NSs (e.g.Mazeh et al. 2022;Yuan et al. 2022;Zheng et al. 2022;Yi et al. 2022;Escorza et al. 2023;Lin et al. 2023;Zhao et al. 2023).The nature of these candidates as NS + MS binaries is quite uncertain.Among the objects in group (a), there is no doubt of the presence of a NS, but the companions -which have not been detected electromagnetically -may be massive WDs.Among those in group (b), the NS has not been detected, and the minimum dynamically implied mass is well below the Chandrasekhar limit.Many of the candidates in group (b) are likely to host massive WDs rather than NSs.
By precisely monitoring the astrometric light-centroid "wobble" of nearly two billion stars, the Gaia mission opens a new window on the Galactic binary population (see El-Badry 2024a, for a recent review).The mission's 3rd data release ("DR3") in June 2022 included orbital solutions for about 1.7 × 10 5 astrometric binaries, including 3 × 10 4 joint astrometric + radial velocity (RV) solutions (Gaia Collaboration et al. 2023b).Although DR3 employed stringent quality cuts (Halbwachs et al. 2023) and represents only a small fraction of the binary sample that will be accessible in future data releases, the DR3 binary sample was already more than an order of magnitude larger than all samples of binary orbits in the previous literature.Gaia's particular sensitivity to longperiod orbits (P orb ∼ (1 − 3) years in DR3) -and the fact that astrometric data provides constraints on binary inclinations that are not accessible with RVs alone -has already enabled the discovery of unexpected binary populations, including three stellar-mass black holes (BHs) in au-scale orbits (El-Badry et al. 2023a,b;Gaia Collaboration et al. 2024) and a population of WD + MS binaries with similar orbital separations (Shahaf et al. 2023b,a;Yamaguchi et al. 2024).The population of wide NS + MS binaries studied in this paper is closely related to these two populations.
Here we present results from a follow-up program of Gaia astrometric binaries suspected to contain NSs.One of our candidates, Gaia NS1 (J1432-1021), was already studied in detail by El-Badry et al. (2024).This object has the highest inferred dark companion mass of any of the objects in our sample, and it is the only object for which the minimum companion mass from RVs alone is well above the Chandrasekhar mass, independent of astrometric constraints on the inclination.We suspect that most of the other binaries in the sample also host NSs, but because the unseen companions have masses near the maximum WD mass, other possibilities cannot be ruled out definitively.
The remainder of this paper is organized as follows.
Section 2 describes selection of our initial sample from Gaia DR3, and Section 3 summarizes our spectroscopic follow-up and measurement of metallicities.We infer parameters of the luminous stars by fitting their spectral energy distributions in Section 4. In Section 5, we carry out joint fits of the astrometry and our follow-up RVs.In Section 6, we discuss the nature of the dark companions, the binaries' possible formation histories, their Galactic orbits and possible abundance anomalies, and the BH/NS mass distribution.We summarize our findings in Section 7. We discuss spurious astrometric solutions in Appendix A and limits on possible optical contamination from WD companions in Appendix B. Tables of RVs and orbital parameters are provided in Appendices C and D.
2. SAMPLE SELECTION We selected targets from Gaia DR3 following the general approach outlined by Shahaf et al. (2019).In brief, the astrometric "triage" algorithm seeks to identify sources whose astrometric orbits are so large -given their orbital period -that they cannot be explained by any luminous star companion, or by a companion that is a close binary containing two luminous stars.Shahaf et al. (2023b) applied this algorithm to astrometric binaries published in Gaia DR3, producing a catalog of 177 candidates in which the astrometric solution and assumed luminous star mass implies the secondary must be a WD, NS, or BH.These classifications are, however, contingent on the validity of the Gaia astrometric solutions, which are in some cases spurious.
Calculating the mass of a star's unseen companion from its astrometric solution requires an estimate of the mass of the star.In constructing their candidate sample, Shahaf et al. (2023b) used the IsocLum mass estimates calculated by Gaia Collaboration et al. (2023b).These estimates, which are available through the Gaia archive in the gaiadr3.binarymasses catalog, were inferred by comparison of the extinction-corrected colors and absolute magnitudes to a grid of PARSEC isochrones, with a prior that the metallicity is close to solar.These masses can thus be overestimated if the metallicities are subsolar, or underestimated if they are supersolar.Their validity is also contingent on the validity of the extinction estimates and on the assumption that a single star contributes to the observed photometry.We improve the mass estimates in Section 4. Shahaf et al. (2023b) noted that their compact object binary candidate sample appeared to fall within two populations in the mass-eccentricty plane: one with a mean mass close to 0.6 M ⊙ and eccentricities below 0.2, and another with an apparent mean mass of ∼ 1.3 M ⊙ and a broad eccentricity distribution.It would be natural, they noted, to identify these two populations with WDs and NSs.Our follow-up has shown that the division is probably not so simple: some systems in the high-mass, high-eccentricity sample unambiguously host white dwarfs, as revealed by strong UV excess (e.g.Ganguly et al. 2023).Other objects have high eccentricities but companion masses below 1 M ⊙ , which are implausibly low for a NS.Despite these caveats, our follow-up has shown that most of the low-eccentricity, lower-mass compact object candidates identified by Shahaf et al. (2023b) are WDs, and at least some of the higher-eccentricity, higher-mass companions are NSs.We return to this discussion in Section 6.2.
We initiated spectroscopic follow-up for a majority of the NS candidates identified by Shahaf et al. (2023b) in June 2022.We also carried out RV follow-up observations of some astrometrically-selected NS candidates not included in the Shahaf et al. (2023b) sample because their orbital periods slightly exceed 1000 days.Several of our candidates were also listed in other samples of compact object candidates from Gaia DR3, including the sample curated by Andrews et al. (2022).The observations presented here were obtained before June 2024, but our program is ongoing.

Rejection of spurious astrometric solutions
For some candidates, our RV follow-up soon showed the Gaia astrometric solution to be spurious or to have significantly underestimated uncertainties.Examples are shown in Appendix A. About a quarter of candidates with good astrometric quality flags turned out to be spurious.1While the fraction of all astrometric solutions that are spurious is small, RV follow-up over a significant fraction of an orbit is critical for vetting astrometric solutions of unusual objects: our follow-up has demonstrated that incorrect solutions do exist, even among solutions with favorable goodness of fit and other Gaia quality flags.

Completeness of the sample
Among the NS candidates identified by Shahaf et al. (2023b) and Andrews et al. (2022), our sample includes all systems that (a) are brighter than G = 15, (b) have best-fit companion masses M 2 > 1.25 M ⊙ from joint fitting of astrometry and RVs, (c) were not found to have spurious solutions or significantly underestimated astrometric uncertainties through RV follow-up, and (d) were observed over at least half an orbit.Properties of these sources are summarized in Table 1.The Shahaf et al. (2023b) and Andrews et al. (2022) samples mainly contain sources near the main sequence with inferred luminous star masses M ⋆ ≲ 1.3 M ⊙ .2While NS companions to evolved stars and more massive MS stars are likely to exist, these in most cases cannot be distinguished from luminous stars or tight luminous-star binaries based on astrometry alone.
Table 3 in Appendix A provides a summary of our RV follow-up and our current assessment of the viability of all candidates from the Shahaf et al. (2023b) and Andrews et al. (2022) samples.We defer a full description of our follow-up program -including RVs of suspected WDs and all objects that turned out to have spurious astrometric solutions -to future work.We suspect that our sample contains most of the NS-hosting binaries with M 2 ≳ 1.25 M ⊙ and astrometric solutions published in DR3.A handful of likely good candidates are not included because our RV follow-up has not yet covered enough of the orbits to confirm the astrometric solution with high confidence, and another handful were excluded because they are too faint (G > 15) to be amenable for RV follow-up with the instruments at our disposal.Our sample does not contain any low-mass NSs with M < 1.25 M ⊙ , a mass limit below which some NSs likely do exist (Ferdman et al. 2014;Martinez et al. 2015).As we discuss throughout the paper, it becomes increasingly challenging to distinguish between NSs and massive WDs at lower masses.

Summary of the sample
Basic properties of our NS candidates are listed in Tables 1 and 4. Figure 1 compares our candidates to the full sample of astrometric binaries (solution types Orbital and AstroSpectroSB1) published in DR3.Three candidates have AstroSpectroSB1 solutions: J0152-2049, J2145+2837, and J1150-2203; the rest have Orbital solutions.The upper left panel of Figure 1 shows the sources on the extinction-corrected color-magnitude diagram.We estimate the extinction for sources in the north (δ > −30 deg) using the 3D dust map from Green et al. (2019); we use the map from Lallement et al. (2022) for sources farther south.All our targets are solar-type stars on the main sequence, with absolute magnitudes and colors suggesting luminous star masses of (0.7 − 1.3) M ⊙ .Several candidates are near the blue edge of the main sequence.A potential concern is that this could be due to blue excess from a hot WD companion.However, our analysis of the sources' full spectral energy distributions -in particular, the lack of UV excess -speaks against this possibility (Appendix B).
The upper right panel of Figure 1 shows orbital periods and distances.Most of the binaries in our sample are within 1 kpc of the Sun and have periods between 100 and 1000 days.There are no binaries with orbital periods close to 1 year owing to the degeneracy between such orbits and parallactic motion.The period distribution peaks near 600 days, likely because short-period orbits are smaller and can only be resolved at close distances, while significantly longer orbits would not have been well-sampled during the ∼1000-day observing window for DR3 solutions.The period distribution of the NS candidate sample is fairly similar to that of all astrometric binaries.
The lower left panel of Figure 1 shows the sources' distribution on the sky in Galactic coordinates, with the Galactic center at the center.Both our candidates and the full astrometric binary sample are distributed all across the sky.Some evidence of the Galactic disk is evident in the distribution of all binaries, but the distribution is heavily affected by the Gaia scanning law, and most of the NS candidates are at high latitude.
Finally, the lower right panel of Figure 1 shows the masses of the luminous stars in our sample and the minimum mass ratio inferred from their astrometric mass ratio function (AMRF; Shahaf et al. 2019).The luminous star masses plotted here are taken from the gaiadr3.binarymasses table following Shahaf et al. (2023b); more accurate masses for our candidates are measured in Section 4. The astrometric mass ratio functions are also calculated based on Gaia data alone following Shahaf et al. (2023b), without accounting for our follow-up RVs.By virtue of our selection, the objects in our sample have among the largest minimum mass ratios of binaries with solutions published in DR3.
That the unseen companions to objects in our sample are massive can be appreciated intuitively from Figure 2, which compares their periods and the physical size of their photocenter orbits to other binaries with astrometric solutions in DR3, including Gaia BH1 and BH2.According to Kepler's 3rd law, more massive and darker companions are found above and to the left of the population of luminous binaries and triples in this parameter space (see El-Badry et al. 2023a).NS companions are expected to be found between the populations of luminous binaries and BH companions, and this is indeed where our candidates are clustered.Close inspection of Figure 2 will reveal that a few additional binaries fall above the candidates studied here, suggestive of higher masses.These are primarily sources where our RV follow-up showed the astrometric solution to be spurious, and red giants, for which a more massive MS companion could not be ruled out.

SPECTROSCOPIC FOLLOW-UP
We measured multi-epoch RVs for all targets, primarily using high-resolution spectrographs on two 2m-class telescopes.Our targets are bright, with most having apparent magnitudes G = 13 − 14.We additionally obtained single-epoch higher-SNR spectra for most targets with 8-10m class telescopes.We describe all the spectroscopic observations below.

FEROS
We obtained 129 spectra with the Fiberfed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) on the 2.2m ESO/MPG telescope at La Silla Observatory (programs P109.A-9001, P110.A-9014, P111.A-9003, P112.A-6010, and P113.26XB).Some observations used 2 × 2 binning to reduce readout noise at the expense of spectral resolution; the remainder used 1 × 1 binning.The resulting spectra have resolution R ≈ 40, 000 (2 × 2 binning) and R ≈ 50, 000 (1 × 1 binning).Exposure times ranged from 1200 to 3600 seconds, depending on source brightness.We reduced the data using the CERES pipeline (Brahm et al. 2017), which performs bias-subtraction, flat fielding, wavelength cali-  -Comparison of NS candidates presented here (red stars) to the rest of the Gaia DR3 astrometric binary sample.At fixed period, dark and massive companions produce larger photocenter orbits than normal stellar companions.For typical solar-type primaries, NS companions produce photocenter orbits that are smaller than those of BH binaries, but larger than those of luminous binaries and triples.
bration, and optimal extraction.The pipeline measures and corrects for small shifts in the wavelength solution during the course a night via simultaneous observations of a ThAr lamp with a second fiber.

TRES
We obtained 155 spectra using the Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész 2008) mounted on the 1.5 m Tillinghast Reflector telescope at the Fred Lawrence Whipple Observatory (FLWO) on Mount Hopkins, Arizona.TRES is a fibrefed echelle spectrograph with a wavelength range of 390-910 nm and a resolving power of R ∼ 44, 000.The spectra were extracted as described in Buchhave et al. (2010).

MIKE
We obtained 7 spectra with the Magellan Inamori Kyocera Echelle (MIKE) spectrograph on the Magellan Clay telescope at Las Campanas Observatory (Bernstein et al. 2003).We used the 0.5" slit and exposure times ranging from 600 to 1200 seconds, yielding spectral resolution R ∼ 40, 000 on the blue side and R ∼ 55, 000 on the red side.The typical SNR at 5800 Å was ∼ 30 and the total wavelength coverage was ∼ 3330 − 9680 Å.We reduced the spectra with the MIKE Pipeline within CarPy (Kelson et al. 2000;Kelson 2003) and subsequently flux-calibrated them using observations of a standard star.We merged the orders into a single spectrum, weighting by inverse variance in the overlap regions.

HIRES
We obtained 6 spectra using the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the 10m Keck I telescope on Maunakea.The data were obtained and reduced using the standard California Planet Survey setup (CPS; Howard et al. 2010), including use of the C2 decker (0.86 arcseconds × 14 arcseconds), which yields spectra with R ≈ 55, 000 and wavelength coverage over most of 3700-8000 Å.We used 600 second exposures, yielding a typical SNR of 40 per pixel at 6000 Å.The CPS reduction includes sky-subtraction using the long C2 decker.We merged spectra from individual orders using the same procedure as with the MIKE data.

PEPSI
We obtained 10 spectra using the Potsdam Echelle Polarimetric and Spectroscopic Instrument (PEPSI; Strassmeier et al. 2015) spectrograph on the Large Binocular Telescope in binocular mode.We used the 300 µm fiber, which has a diameter of 2.3 arcsec on sky, and the CD2 and CD5 cross-dispersers on the blue and red side, respectively.Exposure times ranged from 300 to 1200s.The spectra were reduced and orders were merged as described in Strassmeier et al. (2018); they cover the wavelength ranges of 4222-4792 Å and 6236-7433 Å with spectral resolution R ≈ 50, 000.

ESI
We obtained 1 spectrum with the Echellette Spectrograph and Imager (ESI; Sheinis et al. 2002) on the 10m Keck II telescope on Maunakea.We used a 300 second exposure with the 0.3 arcsec slit, yielding a resolution R ≈ 12000 and SNR ≈ 50, with useful wavelength coverage of 3900-10000 Å.We reduced the data using the MAuna Kea Echelle Extraction (MAKEE) pipeline, which performs bias-subtraction, flat fielding, wavelength calibration, and sky subtraction, and we refined the wavelength solution using telluric absorption lines.

Other spectroscopy
Early in our follow-up program, we measured RVs for several candidates using spectra from lower-resolution spectrographs, including Keck/DEIMOS (Faber et al. 2003), Magellan/MagE (Marshall et al. 2008), and Palomar P200/DBSP (Oke & Gunn 1982).These RVs were used to rule out some candidates with spurious solutions but were not included in our final analysis due to their larger uncertainties (≳ 3 km s −1 ).A few objects in our sample also have archival spectra from the LAMOST survey (Cui et al. 2012).However, these typically have RV uncertainties of at least a few km s −1 -about 100 times larger than our observations -so we opted not to include them in our analysis.

RVs and offsets
We measure RVs for each echelle order by crosscorrelating a synthetic spectral template with the normalized spectrum.We report the median across orders as the measured RV for each epoch and calculate the uncertainty as the standard deviation across orders divided by the square root of the number of orders.We use a Kurucz spectral template from the BOSZ library (Bohlin et al. 2017) with surface gravity log(g/cm s −2 ) = 4.0 and effective temperature and metallicity matched to each target (Section 3.9).We use non-rotating templates for all targets except J0230+5950, which is rotating with v sin i ≈ 15 km s −1 .
For the FEROS spectra, we use 15 orders covering wavelengths between 450 and 670 nm when calculating RVs.For the TRES spectra, we use 31 orders covering 420-670 nm.The median uncertainty of the FEROS RVs is ≈ 0.06 km s −1 , while the median uncertainty of the TRES RVs is ≈ 0.05 km s −1 .
As described by El-Badry et al. (2024), we fit for a single global RV offset between TRES and FEROS.We found an offset of 0.16 km s −1 , in the sense that RV FEROS = RV TRES − 0.16 km s −1 .This offset has already been applied to our reported RVs.Because we generally only obtained one spectrum per target for MIKE, HIRES, PEPSI, and ESI, we did not fit for offsets for these instruments, but instead adopted a conservative 0.5 km s −1 minimum uncertainty for all the RVs measured with these instruments.All the RVs are listed in Appendix D and provided in machine-readable form as supplementary material.

Metallicities
The radius and temperature of a MS star of a given mass depend on its metallicity, so metallicity estimates are important for reliable mass estimates of the luminous stars, and in turn, mass estimates of their companions.We collate metallicity measurements from several sources: 3.9.1.BACCHUS For the MIKE, HIRES, and FEROS spectra, we measured metallicities and atmospheric parameters using the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS; Masseron et al. 2016;Hayes et al. 2022).The code performs 1D LTE spectral synthesis to determine stellar parameters from Fe excitation/ionization balance; i.e., the requirement that lines with different excitation potentials all imply the same Fe abundance.The reported metallicity [Fe/H] is the mean Fe abundance calculated over lines in the VALD atomic linelist with a wavelength coverage of 4200 to 6700 Å.Here we assume that the detailed abundance pattern traces solar values; detailed abundances of these stars will be investigated in future work.The errors reported by BACCHUS represent the scatter in the implied abundances between the different lines and methods of abundance calculations but do not take into account other systematic uncertainties.
3.9.2.SPC We fit the TRES spectra using the Stellar Parameter Classification (SPC) tool (Buchhave et al. 2012), which cross-correlates the normalized TRES spectra with a grid of synthetic spectra in the wavelength range of 5050 to 5360 Å, centered on the Mg I b triplet.SPC infers the metallicity [M/H], effective temperature, T eff , and surface gravity log g by fitting the peaks of the crosscorrelation functions with a three-dimensional polynomial in stellar parameters.Given systematic uncertainties in the synthetic stellar spectra, error floors on the derived [M/H] and T eff values are ∼ 0.08 dex and ∼ 50 K, respectively (Buchhave et al. 2012;Furlan et al. 2018).

Gaia XP metallicities
All of the stars in our sample have metallicity estimates (along with T eff and log g) calculated as described by Andrae et al. (2023) using Gaia XP low-resolution spectra.For bright stars within the temperature range of our sample, the expected precision of these metallicities is ∼ 0.1 dex.However, the metallicities published by Andrae et al. (2023) are expected to be less reliable for astrometric binaries -even those with dark companions -because they are calculated based on the gaiadr3.gaiasource parallax, which comes from a 5-parameter astrometric solution that neglects orbital motion.These parallaxes can be significantly biased for astrometric binaries with large photocenter wobbles.We therefore recalculate metallicities for our targets using the parallaxes in the gaiadr3.nsstwo body orbit catalog, otherwise using 2.0 1.5 1.0 0.5 0.0 0.5 [Fe/H], high resolution spectra Fig. 3.-Metallicities measured from high-resolution spectra (yaxis), and from low-resolution Gaia XP spectra (x-axis).The XP metallicities are inferred with a modified version of the Andrae et al. (2023) model that uses the parallax from the astrometric binary solution rather than the one from the 5-parameter solution reported in the gaia source catalog.The two independent metallicity measurements are in good agreement.
the same empirically trained machine-learning models as Andrae et al. (2023).

Metallicity comparisons
For each source, we report both the Gaia XP metallicity and a metallicity measured from a high-resolution spectrum in Table 2. Several sources were observed with more than one high-resolution spectrograph.In these cases, we adopt the metallicity from the highest-SNR spectrum.The results are reported in Table 2 and shown in Figure 3.The agreement between the Gaia XP and high-resolution metallicities is good, implying that the XP metallicities are robust for this sample and that highresolution spectra will not be essential for bulk metallicity measurements in future follow-up efforts, though they are required for RV measurements.
We also investigated the consistency of metallicities measured from high-resolution spectra by different pipelines.For sources with more than one high resolution spectrum, we found a median absolute deviation of 0.08 dex between measurements, which is comparable to the reported uncertainties.

SPECTRAL ENERGY DISTRIBUTIONS AND
LUMINOUS STAR MASSES We inferred the masses and evolutionary states of the luminous stars by fitting their broadband spectral energy distributions (SEDs) with single-star models.In the optical, we use the synthetic photometry in the SDSS ugriz bands constructed from Gaia BP/RP spectra and empirically calibrated as described by Gaia Collaboration et al. (2023a).We supplement this with near-infrared photometry from the 2MASS (Skrutskie et al. 2006) and WISE (Wright et al. 2010) surveys.To search for possible light contributions from hot WD companions, we also retrieved UV photometry from GALEX (Martin et al. 2005), but we did not include it in our SED fits.
The SEDs of all sources are shown in Figure 4.While optical and near-infrared photometry is available for all sources, several sources are outside the published footprints of GALEX in one or both of its UV bands.In cases where a source was within the footprint of a GALEX observation but was was not detected, we plot the 3σ upper limit.Sources with no upper limits or detection shown in the NUV or FUV were not observed in that band.
We fit the SEDs using MINESweeper (Cargile et al. 2020), a code designed for joint modeling of stellar photometry and spectra.We only use the photometric modeling capabilities of the code but place a prior on the present-day metallicity from spectroscopy.The free parameters of the fit are the initial mass and metallicity of each star, its age (as parameterized by the "equivalent evolutionary point"), the parallax, and the foreground extinction.For each call to the likelihood function, the mass, metallicity, and equivalent evolutionary point are converted to a radius and effective temperature using MIST isochrones (Choi et al. 2016), and are then used to predict a model SED that is compared to the data.We place priors on the parallax from the Gaia 12-parameter binary solution and set an age upper limit of 13 Gyr.We also use a metallicity prior from the high-resolution spectra (Table 2), and an extinction prior from 3D dust maps.We again use the Green et al. (2019) map for sources with declination δ > −30 deg, and the Lallement et al. (2022) map for sources farther south.
The SEDs and best-fit models of all candidates are shown in Figure 4, and best-fit parameters for each target are listed in the upper right corner of each panel.Most of the luminous stars have masses near 1 M ⊙ .This is also true for the full astrometric binary catalog published in Gaia DR3 and mostly reflects the distance and flux limits of the astrometric binary sample.Several of the MS stars are somewhat evolved (i.e., their radii are up to 70% larger than stars of the same mass at the zero-age main sequence).
We removed one initial target, Gaia DR3 ID 747174436620510976, from the sample because it was not possible to obtain a good single-star fit to the SED without changing the Gaia parallax.The star has near solarmetallicity star and an SED that suggests T eff ≈ 5450 K and R ≈ 0.60 R ⊙ .However, this radius is too small for any main sequence star of the observed effective temperature.A possible explanation is that the Gaia parallax is underestimated and the true radius is ≈ 0.9R ⊙ , as expected for a main-sequence star.This would also cast doubt on the other parameters of the astrometric solution, and our follow-up RVs are in modest tension with the predictions of the astrometric solution (Appendix A), so we removed the source from the sample.
The best-fit parameters from SED fitting are listed in Table 2.We set a minimum uncertainty of 0.03 M ⊙ on M ⋆ to account for systematic uncertainties in the stellar models.This uncertainty propagates to the uncertainties on M 2 calculated in the next section.
Our mass constraints are based on single-star evolutionary models.This is reasonable because the expected lifetimes of the companions' progenitors are much shorter than the lifetimes of the luminous stars.This means than any mass transfer would have occurred before the  2. Parameters of the luminous stars.We list both metallicities measured from high-resolution spectra and those measured from low-resolution Gaia XP spectra following Andrae et al. (2023).The high-resolution metallicities are used as a prior when fitting the broadband SEDs.
The metallicities are compared in Figure 3, and the SED fits are shown in Figure 4.
luminous stars experienced any significant evolution.As discussed by El-Badry et al. (2024), there is also no plausible evolutionary scenario in which the luminous stars have significantly lower masses than implied by singlestar models.

JOINT FITTING OF ASTROMETRY AND RVS
Our joint modeling of the RVs and Gaia astrometry assumes a Keplerian 2-body orbit.The luminous star has mass M ⋆ , while the companion has mass M 2 .The orbit is specified by its period, P orb , eccentricity, e, periastron time, T p , and three angles describing the orientation: the inclination, i, longitude of the ascending node, Ω, and argument of periastron, ω.The center-of-mass motion of the binary is described by a systemic velocity, γ, and five astrometric parameters: the right ascension α and declination δ, the parallax ϖ, and the proper motions µ * α and µ δ .
The semimajor axis a is set by Kepler's 3rd law: while the semimajor axis of the luminous star's orbit is where q = M 2 /M ⋆ is the mass ratio.The RV semiamplitude of the star's orbit is Finally, the angular photocenter semimajor axis is given by å0 Here d is the distance to the binary and ϵ is the flux ratio in the G-band.We define ϵ = F G,2 /F G,⋆ , where F G,2 and F G,⋆ represent the G-band flux of the companion and the luminous star, respectively.For a dark companion (ϵ = 0), the photocenter simply traces the primary (i.e., the star whose RVs are being measured).A nonzero value of ϵ increases the primary semimajor axis, a 1 , that corresponds to a given photocenter semimajor axis.
The Gaia astrometric solutions are expressed as constraints on 4 or 6 Thiele-Innes parameters, for Orbital and AstroSpectroSB1 solutions, respectively: In the convention used here and in the Gaia archive, the astrometric parameters A, B, F , and G have angular units (mas), while the spectroscopic parameters C and H have physical units (au).In the Gaia data processing, C and H are constrained by the measured RVs, while A, B, F , and G are constrained by astrometry (see Pourbaix et al. 2022).Our standard fit has 14 free parameters: P orb , e, M ⋆ , M 2 , i, Ω, ω, ϖ, α, δ, µ * α , µ δ , T p , and γ.In the fiducial modeling, we assume a dark companion (ϵ = 0).We also experimented with a 15-parameter fit in which ϵ is left free (Section 5.4).
For each call to the likelihood function, we construct the predicted vector of Gaia-constrained parameters, θ Gaia .For objects with Orbital solutions, this is given by ) while for those with AstroSpectroSB1 solutions, We then calculate a likelihood term that quantifies the difference between these quantities and the Gaia constraints: (13) Here µ Gaia and Σ Gaia represent the vector of best-fit parameters constrained by Gaia and their covariance matrix, which we construct from the corr vec parameter reported in the Gaia archive.
We additionally predict the RVs of the luminous star at the array of times t i at which we obtained spectra.This requires calculating K ⋆ from Equation 3 and iteratively solving Kepler's equation (e.g.Murray & Dermott 1999).When then calculate a RV term in the likelihood, where RV i and RV pred (t i ) are the measured and predicted RVs, with their uncertainties σ RV,i .The full likelihood is then given by To assess the relative importance of the Gaia data and our follow-up RVs in constraining the orbits, we also carry out a fit in which we omit the ln L RVs term from Equation 15.This essentially samples from the Gaia covariance matrix, but it incorporates our constraints on M ⋆ , transforms from the Thiele-Innes coefficients to a more physically interpretable set of parameters, and allows a direct constraint on M 2 .
We use flat priors on all parameters except M ⋆ , for which we use a normal distribution informed by our SED fit (Section 4).We sample from the posterior using emcee (Foreman-Mackey et al. 2013) with 64 walkers, taking 3000 steps per walker after a burn-in period of 3000 steps.We choose the initialization point for the sampler by transforming the best-fit Gaia parameters to Campbell elements using nsstools (Halbwachs et al. 2023).

Inclination sign degeneracy
Astrometric orbits describe the plane-of-the-sky motion of a binary's photocenter.With astrometry alone, the data are always equally consistent with an orbit that has inclination +i and one that has inclination −i.Mathematically, this is evident from the fact that A, B, F , and G depend only on cos i, an even function.RVs break this degeneracy.However, because orbit-fitting posteriors are usually multi-modal, our MCMC sampling often fails to converge if initialized far from the best-fit solution.For this reason, we separately initialize two sets of MCMC chains near the positive and negative astrometry-only solutions, and retain only the chain that reaches a higher maximum likelihood.

Results
The results of all our joint fits are reported in Table 4 in Appendix C. Constraints on a few parameters are also listed in Table 1.Figures 5 and 6 compare our measured RVs to predictions of the Gaia-only solution (cyan) and the Gaia+RVs solution (black).We plot predicted RV curves for 50 random samples from the posterior in order to show the uncertainty in these predictions.For the Gaia-only fits, we fix the center-of-mass RV (which is not constrained by astrometry) to the maximum-likelihood value from the joint fit.For all the targets shown here, the RVs are consistent with predictions of the Gaia-only solution, meaning that at least some samples predict RV curves that overlap with the observed RVs.
Not surprisingly, the predicted RV curves from the joint RVs+astrometry fits are much more tightly constrained than those from the Gaia-only fits.The residuals in Figure 6 show that for most targets, the RVs are generally consistent with the joint solution within (1 − 2)σ.Some systems have χ 2 /N RVs < 1.0, suggesting that their RV uncertainties are overestimated.This could occur if a small number of outlier orders dominate the standard deviation in the RV uncertainty calculation.A few candidates have χ 2 /N RVs = 1.0 − 2.0, implying modestly underestimated RV uncertainties.

Consistency of the astrometric and RV constraints
Figure 7 compares our constraints on several fitting parameters, as well as parameters that are transformations of them, from the Gaia-only and Gaia+RVs fits.For most systems, the uncertainties in all parameters are much smaller in the Gaia+RVs fits than in the Gaia-only fits.This is true even for parameters such as the inclination that are not constrained directly by RVs, because the Gaia-only constraints include significant covariances between parameters.
The two sets of parameters are generally consistent within (1−2)σ, suggesting that the astrometric solutions and their uncertainties are generally reliable.There is little evidence of systematic biases in most of the Gaiaonly parameters, though the inferred M 2 values are on average lower in the joint fits than in the Gaia-only fits.As we discuss in Section 5.3, this is likely a consequence of how the candidates were selected.
El-Badry et al. ( 2024) found that the astrometry-only and astrometry+RVs constraints for Gaia NS1 (J1432-1021) were inconsistent at the 3σ level, suggesting that the astrometric uncertainties for that source were underestimated.Similarly, Chakrabarti et al. (2023) and Nagarajan et al. (2023) found evidence for underestimated uncertainties in the astrometric orbit of Gaia BH1. Figure 7 implies that most of our candidates have more robust astrometric uncertainties: indeed, Gaia NS1 has the most discrepant solution among all the objects in the sample.This could arouse suspicion that the companion's mass is simply overestimated.However, El-Badry et al. (2024) showed that RVs alone require the companion to be the most massive in the sample, even if the astrometric orbit constraints are disregarded.
Considering the 4 Thiele-Innes parameters A, B, F , and G, only 22 out of 84 parameters across 21 objects differ by more than 1σ, and only 5 differ by more than 2σ.The Gaia-only and Gaia+RVs parallaxes are consistent within 1σ for 17 candidates, and within 2σ for 20 candidates.

Mass functions
Figure 8 compares several limits on the unseen companion masses between the Gaia-only and Gaia+RVs fits.The left panel shows the RV mass function, which represent a lower limit on the companion mass that corresponds to the limiting case where the orbit is edgeon and the luminous star is massless.K ⋆ is the RV semiamplitude as measured from the joint fit (Equation 3).The middle panel shows the astrometric mass function, For a dark companion, That is, f (M 2 ) ast represents a lower limit on the companion mass that incorporates the known inclination from astrometry but treats the luminous star as massless.Finally, the right panel shows the minimum companion mass implied by the RV mass function and our assumed luminous star mass, for an edge-on orbit.We obtain this by setting i = 90 deg and solving Equation 17for M 2 , setting M ⋆ to the value inferred from the SED (Table 2).We propagate uncertainties on all parameters by calculating these mass functions directly from the MCMC samples.
Although the Gaia-only and Gaia+RVs mass functions are generally consistent, the Gaia+RVs constraints on average favor lower mass functions.This means that, on average, the observed RV semi-amplitudes are somewhat lower than predicted by the astrometric solutions.Since our targets are selected from the tails of the companion mass function distribution (e.g. Figure 1), it is expected that more objects in the sample will scatter toward higher masses than toward lower masses.That is, binaries with noise scattering their best-fit masses toward lower values would not have entered the sample in the first place.

Flux ratio constraints
Our modeling thus far has assumed that the companion is dark.In this case, the semimajor axis of the photocenter orbit, å0 , will be identical to the semimajor axis of Individual lines show random draws from the posterior.For AstroSpectroSB1 solutions, the systemic RV is set to the value measured by Gaia; for Orbital solutions, it is set to the best-fit value from our joint astrometry+RV fits (Figure 6).Systems with long periods (P orb ≳ 700 days) have significant uncertainties in their Gaia-only periods and phases at the time of our observations.The RVs are broadly consistent with the Gaia predictions.
Gaia only Gaia + RVs Fig. 7.-Orbital and astrometric parameters inferred from joint fits of the Gaia solution and our RVs (vertical axis) and the Gaia solution alone (horizontal axis).The photocenter semi-major axis, å0 , and the Thiele-Innes coefficients, A TI , B TI , F TI , and G TI , are derived quantities, while the other parameters are free parameters in our fits.Uncertainties on most parameters are much smaller in the joint fits, but the two sets of constraints are generally consistent at the (1 − 2)σ level.Inclinations are mirrored about 90 deg (i.e., 100 deg is shown as 80 deg) for visualization only.
the star whose RVs are being measured, å1 .To explore the possibility of a luminous companion, we tried repeating the fit while adding the flux ratio ϵ (Equation 4) as a free parameter.
We first tried fitting with no bounds on ϵ.In this case, we found negative best-fit values of ϵ for a ma-0.0 0.5 1.0 Gaia only Gaia + RVs Fig. 8.-Comparison of mass functions for the unseen companions inferred from the Gaia solution alone (horizontal axis) and joint fitting of that solution and our RVs (vertical axis).Left and middle panels show the spectroscopic and astrometric mass functions.Right panel shows the minimum companion mass implied by the RV mass function if the orbit were edge-on.The mass functions inferred from Gaia+RVs are in most cases consistent with those from Gaia alone, but they are systematically lower.Because our candidates were selected on the basis of high astrometric mass functions, and massive companions are intrinsically rare, their companion masses are overestimated on average.That is, systems with underestimated masses did not enter the sample.
jority of the sources, with a median and ±1σ range of ϵ = −0.01 ± 0.03 across the full sample.Most sources have ϵ constrained with a 1σ uncertainty of σ ϵ ≲ 0.02.Negative values of ϵ are unphysical, but mathematically allowed (e.g.Equation 4).A preference for negative ϵ simply reflects observed RV variability with a lower amplitude than predicted by the best-fit astrometric solution with ϵ = 0; i.e., it is a consequence of the same selection effect that causes the mass functions in Figure 8 to be overestimated.
Next, we imposed a constraint that ϵ must be positive.In this case, the posterior constraints on ϵ pile up against 0 for most sources.The median upper limit on ϵ is ϵ < 0.012 (1σ), or ϵ < 0.032 (2σ).That is, the combination of RVs and astrometry strongly disfavor luminous companions, because luminous companions would imply a larger a 1 for fixed å0 , and this would imply a larger RV variability amplitude than is observed.

Astrometric phase coverage
One possible concern is that the Gaia observations could have resulted in poor orbital phase coverage.This can occur, for example, for orbital periods that are close to a year, half a year, or several periods related to the Gaia scanning law (El-Badry et al. 2023a;Holl et al. 2023).While actual epoch-level astrometric data are not published in DR3, it is possible to predict when a source should have been observed using the Gaia observation scheduling tool (GOST) 3 .Given a source's coordinates, GOST returns a list of observation times when the scanning law predicts that the source will transit across the Gaia focal plane.
Figure 9 shows sky-projected photocenter orbits for all candidates.Individual gray lines show predictions for 50 random draws from the posterior of the joint fit.Red points show interpolations of the scan times predicted by GOST onto the orbit corresponding to the marginal-3 https://gaia.esac.esa.int/gost/ized posterior median.Overall, most of the sources have good astrometric phase coverage, with Gaia observations sampling most of the predicted orbital ellipse.
It is not guaranteed that Gaia will actually obtain data on each source at the predicted times, as GOST does not account for gaps between CCDs and issues such as micrometeoroids impacts that cause temporary gaps in the datastream.In each panel of Figure 9, we list both the number of visibility periods predicted by GOST and the number actually used in calculating the astrometric solution (visibility periods used; a visibility period is a group of observations separated from other observations by at least 4 days).The number of visibility periods used ranges from 80 − 100% of the number predicted by GOST, suggesting that the scan times shown in Figure 9 are a reasonable but not perfect approximation of those used in producing the astrometric solutions.The median number of visibility periods used is 21, significantly higher than the minimum of 12 that is required to constrain a 12-parameter astrometric binary solution.

Astrometric solution quality diagnostics
As an additional check on the reliability of our candidates' astrometric solutions, Figure 10 compares several quality diagnostics of these candidates to the larger sample of all astrometric binaries from Gaia DR3.Cyan points show the Gaia astrometric solutions alone, while red points show joint Gaia + RV fits.Hollow points show three sources we found to have spurious solutions after RV follow-up (Appendix A).
First considering the Gaia-only solutions, the NS candidates have goodness of fit values and parallax uncertainties that are broadly representative of the larger astrometric binary sample.They have higher-significance photocenter ellipse measurements (i.e., larger å0 /σ å0 ) than most other astrometric binaries, reflecting the fact that NS companions produce large photocenter orbits at fixed period (Figure 2).The sources with spurious solutions have similar goodness of fit values to those with good solutions, though their parallax uncertainties are on average larger at fixed apparent magnitude.
The uncertainties from the joint Gaia+RV fits are often significantly smaller than those from Gaia alone.
Even the parallax, which is not directly constrained by RVs, has smaller uncertainties in the joint Gaia + RVs fit, because the Gaia solutions include significant covariances between parallax and other orbital parameters that are directly constrained by RVs.We emphasize, however, that the uncertainties shown in Figure 10 are purely statistical.The parallax zeropoint and uncertainties of the astrometric binary solutions have yet to be investigated in detail.Nevertheless, we conclude from Figure 10 and the good agreement between measured and predicted RVs that the astrometric solutions of the candidates in our sample are robust.We also conclude that spurious solutions cannot always be easily identified on the basis of their astrometric uncertainties and quality flags.
6. DISCUSSION 6.1.Nature of the unseen companions The orbits of the binaries in our sample are wellcharacterized, and the companion masses are measured with typical uncertainties of a few percent.However, we have not detected light from the companions, leaving their astrophysical nature uncertain.Our joint fitting of astrometry and RVs places upper limits of a few percent on the G−band flux ratio, ϵ = F G,2 /F G,⋆ (Section 5.4).As we show below, this rules out all plausible non-degenerate companions.
We begin by considering the astrometric mass function (Equation 19).For a luminous companion with flux ratio ϵ and mass ratio q = M 2 /M ⋆ , this quantity can be expressed as (e.g.Shahaf et al. 2019): Given an observationally constrained M ⋆ and f (M 2 ) ast , we can solve Equation 20 for q, and thus for M 2 , for any possible ϵ.In Figure 11, the black line shows this constraint for a representative source in our sample, J2244-2236, which has M ⋆ = 1.00 ± 0.03 M ⊙ and f (M 2 ) ast = 0.503 ± 0.007 M ⊙ .The dynamically-implied mass of the companion increases with its assumed luminosity.Physically, this reflects the fact that a brighter companion dilutes the light from the primary more, such that a higher companion mass is required to explain the same photocenter orbit.
We can then consider possible companion types: 1.A single MS star: The dotted cyan line in Figure 11 shows the expected flux contribution versus mass relation for a MS companion.We take this from a 5 Gyr-old, solar-metallicity MIST isochrone.Because the mass supplied by such a companion is much less than the astrometric constraint for all flux ratios, a single MS companion is ruled out.
2. An inner binary containing two MS stars: The yellow dashed line in Figure 11 shows the mass and flux expected if the companion is a binary containing two MS stars, so the full system is a triple.We assume the two stars have equal mass, since this yields the highest mass-to-light ratio.We use the same isochrone as for a single star but assume twice the total mass and light.This still falls short of the astrometric mass constraint, so the companion cannot be a close MS+MS binary.
3. An inner triple containing three MS stars: The red dot-dashed line in Figure 11 shows the case where the companion is a triple containing three MS stars (i.e., the full system is a quadruple).We assume all three inner components have the same mass.Even this falls below the astrometric constraint for all flux ratios, so the companion cannot be an inner triple.In addition, a quadruple system with four stars in such a tight orbit would be unlikely to be stable over long timescales.
4. An inner binary containing a WD and a MS star: The blue dashed line in Figure 11 corresponds to an inner binary containing a 1 M ⊙ WD and a MS star.We assume that the WD is sufficiently cold that it does not contribute any light in the optical.Such an inner binary could reproduce the observed astrometric mass function if the MS star has a mass of ≈ 0.6 M ⊙ , in which case it would contribute ≈ 4% of the total light.However, as discussed in Section 5.4, the RVs rule out a flux ratio above 2% for this source, because a 1.6 M ⊙ companion would produce larger-amplitude RV variability.Only a tight binary containing a WD with M ≳ 1.1 M ⊙ and a low-mass MS star could match the data.
5. An inner binary containing two WDs: All of the objects in our sample could in principle be explained by a tight inner binary containing two cold and relatively massive WDs.The total mass of the inner binaries would need to be near the Chandrasekhar mass.Chandrasekhar-mass close WD+WD binaries are relatively rare: despite decades of dedicated searches, none have been conclusively identified.Empirical limits suggest a space density comparable to the space density of neutron stars (e.g.Badenes & Maoz 2012).As discussed by El-Badry et al. ( 2024), there are significant evolutionary challenges to forming a massive WD+WD binary within the orbits of the observed luminous stars: such triples are likely to become dynamically unstable during their evolution, particularly given that the initial orbits of the luminous stars would have to have been significantly tighter than observed today when mass loss from the inner binary is accounted for.Nevertheless, tight WD+WD binaries cannot be ruled out based on the observed data alone.
6.A single ultramassive WD: Several of our candidates have best-fit masses above the maximum WD mass mass, which is near 1.37 M ⊙ (Althaus et al. 2022).However, about half have masses below this limit, and given the astrometric uncertainties, a majority of our candidates are consistent at the few-σ level with having M 2 ≲ 1.37 M ⊙ .If the companions are single WDs, all of them would be among the highest-mass WDs known (e.g.Cognard et al. 2017;Caiazzo et al. 2021;Miller et al. 2023;Yamaguchi et al. 2024).
Given that our candidates are selected from the upper tail of the inferred companion mass distri-   bution, the possibility that companion masses are overestimated systematically should be taken seriously.However, at least one object in the sample -Gaia NS1 (J1432-1021) -has a companion mass high enough that it cannot be a single WD: this object has M 2 = 1.90 ± 0.03, about 17σ above the maximum WD mass.
As we discuss below, almost all the binaries in our sample have higher eccentricities than typical WD+MS binaries at similar periods.A simple interpretation is that the companions are NSs and the eccentricities are a result of natal kicks, but the data also suggest that systems containing massive WDs have higher typical eccentricities than those containing ∼ 0.6 M ⊙ WDs, so we cannot completely rule out a scenario in which many of the unseen companions are ultramassive WDs.

7.
A single neutron star: The inferred masses of the unseen companions in our sample are typical of neutron stars in pulsar binaries ( Özel & Freire 2016).We consider this the simplest and most plausible scenario.
6.2.Period-eccentricity relation Figure 12 (left panel) shows the period-eccentricity relation for our candidates.We compare them to the larger sample of ∼ 3000 WD+MS binary candidates selected from Gaia astrometry by Shahaf et al. (2023a) on the basis of their AMRF.Our candidates -which differ from the Shahaf et al. (2023a) sample only in that they have inferred companion masses above 1.25 M ⊙ -have significantly higher eccentricities than those systems at fixed period.The low eccentricities of the WD+MS binaries are likely a result of tidal circularization when the progenitors of the WDs were giants.If the unseen companions in our sample are NSs, their higher eccentricities could could be naturally understood as a result of natal kicks.
Although the eccentricities of the WD+MS binary candidates in Figure 12 are much lower on average than those of ordinary MS+MS binaries at the same periods, the large majority of them are not consistent with zero, and are much higher than expected for binaries that have gone through long periods of stable mass transfer (Phinney 1992;Lorimer 2008).This may be a consequence of these systems having first initiated mass transfer on the asymptotic giant branch rather than on the first giant branch and/or having gone through common envelope evolution rather than stable mass transfer.Similar eccentricities are observed in other populations of WD+MS binaries at these periods, such as blue stragglers and barium stars (Mathieu & Geller 2009;Jorissen et al. 2019;Escorza et al. 2019).
The right panel of Figure 12 shows eccentricities and dark companion masses.A majority of the objects from the Shahaf et al. (2023a) sample with M 2 > 1.25 M ⊙ are included in our sample: those that are not are either faint (G > 15.0) or were excluded because our follow-up showed them to have spurious astrometric solutions or lower dark companion masses.It is evident that while most WD+MS candidates with M 2 ≲ 0.8 M ⊙ have loweccentricity orbits (e < 0.2), systems with M 2 ≳ 1.0 M ⊙ tend to be more eccentric.It is tempting to simply attribute this dichotomy to WD vs. NS companions (e.g.Shahaf et al. 2023b).However, the eccentric population seems to extend to lower masses than expected for NSs.We have also identified a handful of systems in the course of our follow-up that have M 2 > 1 M ⊙ , e > 0.5, and clear UV excess, which points to a WD companion rather unambiguously.These eccentricities could be a result of faster orbital inspiral and/or more asymmetric mass loss in the super-AGB progenitors of mas-sive WDs (Izzard et al. 2010;El-Badry & Rix 2018), or eccentricity-pumping due to massive circumbinary disks formed through the mass transfer process (Dermine et al. 2013).Kozai-Lidov oscillations in systems containing a wide tertiary are another possibility.

Galactic orbits
To explore the Galactic stellar populations our NS candidates are members of and imprints of possible natal kicks, we show their locations in the Toomre diagram in Figure 13.For each binary, we use the center-of-mass velocity from the joint fit and the parallax and proper motion from the Gaia binary solution to calculate the current 3D motion of the binary's center of mass in Galactocentric cylindrical coordinates.We assume that the Sun is 20 pc above the Galactic midplane and 8.12 kpc from the Galactic center and has a 3D velocity vector (V R,⊙ , V ϕ,⊙ , V Z,⊙ ) = (−12.9,245.6, 7.78) km s −1 (Drimmel & Poggio 2018).We perform the same calculation for the other 33,467 binaries with AstroSpectroSB1 solutions.We do not consider Orbital solutions because their center-of-mass RVs are not known.
The results are shown in Figure 13.We compare the NS candidates in our sample to all binaries with AstroSpectroSB1 solutions (left panels) and to those with d = 0.4 − 1.0 kpc (right panels; these have a similar distance distribution to the NS candidates).Dashed lines centered on V ϕ = 225 km s −1 show approximate boundaries of the Galactic thin disk (total velocity < 70 km s −1 ), thick disk (70 − 180 km s −1 ), and halo (> 180 km s −1 ) (e.g.Bensby et al. 2014).Three objects in our sample -J1739+4502, J1432-1021, and J0152-2049are on halo orbits, with total velocities of 370, 350, and 290 km s −1 with respect to the local standard of rest.These objects all have metallicities [Fe/H] < −1.2, implying that their high velocities are mainly the result of membership to an old stellar population, not natal kicks.On the other hand, the 6 systems in the "thick disk" region of the Toomre diagram have metallicities close to solar.This suggests that natal kicks may have played a role in kinematically heating their orbits.The median and middle 68% range of midplane distances, |z|, are respectively 0.37 kpc and (0.11−0.55) kpc.This makes the objects in our sample kinematically cold compared to NS LMXBs, which typically have |z| ≳ 1 kpc (van Paradijs & White 1995; Jonker & Nelemans 2004).This is not surprising since only NSs formed with weak kicks will remain bound in the wide orbits to which astrometry is sensitive.
Metal-poor stars on halo orbits seem to be significantly overrepresented in the NS candidate sample.In particular, 3 of the 21 objects in the sample are unambiguously on halo orbits.In comparison, 0 of the 319 WD + MS binary candidates selected by Shahaf et al. (2023b) with AstroSpectroSB1 solutions have halo-like orbits!Among 11420 AstroSpectroSB1 orbits with d = 0.4−1.0kpc, 61 (i.e., 0.5%) have halo-like orbits, 1636 (14.3%) have thick disk-like orbits, and 9723 (85.1%) have thin disk-like orbits.In contrast, in the NS sample, 3/21 (14%) have halo-like orbits, 6/21 (29%) have thick disklike orbits, and 12/21 (57%) have thin disk-like orbits.Since the halo orbits of the NS candidates cannot be the result of natal kicks, the high fraction of halo or-  2023b) from Gaia DR3 on the basis of their astrometrically-constrained mass ratios.For NS+MS candidates, uncertainties are smaller than symbols.Most of the WD candidates have low eccentricities, probably resulting from (incomplete) tidal circularization when the WD progenitor was a red giant.Our candidates have significantly higher eccentricities than the WD candidates at all periods, which may be the result of kicks during NS formation.Right: mass-eccentricity diagram.Our candidates are all selected to have M 2 > 1.25 M ⊙ .Solid and dashed lines show median and middle 68% across both red and black points.Binaries with M 2 = (1.0 − 1.2) M ⊙ , most of which likely contain high-mass WDs, also have higher eccentricities than systems containing lower-mass WDs.
bits suggest that NSs formed from low-metallicity stars are more likely to survive in Gaia-detectable binaries.As discussed by El-Badry (2024b), BH companions also seem to be overrepresented in the Gaia sample at low metallicity.

Lithium enhancement
To search for spectroscopic anomalies that could result from mass transfer from the unseen companions' progenitors, we compared their spectra to spectra of similar stars observed by the GALAH survey (Buder et al. 2021).We followed the same procedure described by El-Badry et al. (2024): in brief, we compared the rest-frame, resolution matched, and continuum-normalized spectra of our candidates pixel-by-pixel to all high-SNR spectra in GALAH DR3 and visually inspected the closest matches.Most of our candidates have unremarkable spectra and abundances, and we defer a full abundance analysis to future work.However, we highlight one feature that is striking: all three of the metal-poor halo stars in the sample are strongly enhanced in lithium.
Figure 14 compares the spectra of these three stars to those of their closest spectral doppelgängers.The latter are identified using the wavelength range 5650 − 5850 Å, which does not contain any lithium lines.The three NS candidates have Li I λ6708 lines that are clearly significantly stronger than those of the comparison stars shown in Figure 14, and indeed, stronger than those of any stars observed by GALAH with similar atmospheric parameters.
The origin of the excess lithium is not yet understood.Possibilities include pollution by products of hot bottom burning in super-AGB stars, pollution by supernova ejecta, and spallation by high-energy particles (see El-Badry et al. 2024, for further discussion).Because metalpoor stars have thin convective envelopes, only their surface layers are expected to have been polluted.This is likely the reason we only find strong lithium enhancement in the metal-poor stars in our sample: the thickness of the convective envelope is more than 50× greater at solar metallicity than at [Fe/H] = −1.5, meaning that any accreted material will be diluted 50 times more at solar metallicity, and abundance anomalies will be much more subtle.Assuming lithium is mixed uniformly into the outer 0.2% of the stars by mass, (1 − 3) × 10 −11 M ⊙ of lithium must have been accreted in order to explain the observed enhancement.
Our inferred A(Li) are higher than the predicted yields of super-AGB stars with M ≲ 8 M ⊙ , but compatible with predictions for 8 M ⊙ models (see Ventura & D'Antona 2010, who note that these yields depend sensitively on the assumed winds).This could support a scenario where the companions are NSs formed through electron-capture SNe, or ultramassive WDs.However, this scenario is difficult to reconcile with the fact that one of the binaries in question has M 2 = 1.90 ± 0.03 M ⊙ , whereas electron capture SNe are expected to produce low-mass NSs.

Constraints on kicks and mass loss
Models predict that dynamical channels are less efficient for forming wide NS binaries than for BH binaries (Tanikawa et al. 2024).Therefore, we assume the systems in our sample formed from primordial binaries and use their current orbits to place limits on their pre-SN masses and natal kicks.We consider the evolutionary state of the binary that likely immediately preceded the current state: a circular orbit containing a stripped helium star of mass M He star and a 1 M ⊙ MS companion.Given that the orbits in our sample are too tight to accommodate red supergiants at their maximum radii of ∼ 1000 R ⊙ , the binaries would likely have gone through a common envelope event prior to reaching this stage.Survival of such a common envelope event and a wide final orbit is in fact nontrivial (e.g.Kotko et al. 2024), but here we consider only the effects of the SN.
If there is no natal kick, the binary will be unbound if its total mass decreases by more than a factor of 2 during the SN (e.g.Blaauw 1961).The expected post-SN eccentricity in this case is  Where M NS is the mass of the NS, ∆M = M He star −M NS is the mass lost during the explosion, and M ⋆ is the mass of the luminous star.For typical targets in our sample, e ≈ 0.4 and M ⋆ = 1 M ⊙ .This would imply that only of order 1 M ⊙ was lost during the SN if natal kicks were weak.
Kicks complicate this analysis: it is possible for a wellaimed natal kick to keep a binary bound -and even maintain a low eccentricity orbit -with arbitrarily large mass loss during the SN.This, however, requires finetuning, and most orbits will simply be unbound if kicks are strong and mass loss is significant.Inferring the full distribution of kick velocities from observations of binaries like those in our sample is difficult, since NSs born with strong kicks will not make it into the sample in the first place.We can, however, still constrain the kicks and mass loss that likely occurred in the binaries that did survive.
We model the combined effects of natal kicks and instantaneous mass loss using Monte Carlo simulations.-Survival probability (top) and median predicted eccentricity (bottom) for wide helium star + MS binaries when the He star explodes and leaves behind a NS.We assume the initial orbit is circular and consider initial helium star masses of 3, 4, and 5 M ⊙ , and initial periods of 100, 300, and 1000 days.In all cases, we model the companion as a 1 M ⊙ MS star and adopt a final NS mass of 1.4 M ⊙ .The NS is born with a velocity v kick in a random direction.For a 3 M ⊙ He star, orbits are most likely to remain bound if the kick is slow.For M He star > 3.8 M ⊙ , the binary will be unbound by mass loss unless a fortuitously aligned kick allows the NS to catch the escaping companion.
Following Brandt & Podsiadlowski (1995), we predict post-SN orbits for a range of v kick , M He star , and pre-SN orbital period.For each choice of these values, we generate N = 10 6 kick directions distributed uniformly on a sphere and calculate the post-SN orbit via energy and angular momentum conservation.In all cases, we assume a 1 M ⊙ companion, a 1.4 M ⊙ NS, and an initially circular orbit.
The results are shown in Figure 16.The top panels show the fraction of all orbits that remain bound, and the bottom panels show the median eccentricity of those that do.The probability of remaining bound is highest for low M He star and low v kick .However, a significant fraction of orbits remain bound even for M He star of 5 or 10 M ⊙ when paired with a suitable kick velocity.Higher values of M He star require stronger kicks to remain bound: for a given M He star , binaries are most likely to survive when v kick is about half the pre-SN orbital velocity of the MS star.The median predicted post-SN eccentricity of surviving binaries is typically 0.4 − 0.8.These values are slightly higher than the median eccentricity of 0.4 within our sample, but given that Gaia is less sensitive to higheccentricity orbits, the observed eccentricity distribution appears consistent with originating mainly from kicks.
Light curve modeling of stripped-envelope SNe suggests typical ejecta masses of (1 − 3) M ⊙ (e.g.Lyman et al. 2016), which would correspond to He star masses of ∼ (3−5) M ⊙ for the NSs in our sample.For kicks with v kick ≲ 50 km s −1 , about 30% of binaries forming from such He stars at periods typical of our sample would be expected to survive.Such kick velocities are quite low compared to typical values inferred from observations of young pulsars (e.g.Hobbs et al. 2005;Faucher-Giguère & Kaspi 2006), but within the range required for NSs to be retained within globular clusters (Ivanova et al. 2008).Although only a minority of NSs formed in binaries with the periods and companions typical of our sample are likely to survive, the relatively low space density of wide NS+MS binaries suggested by our sample compared to the predicted space density of all dormant NSs 4 allows for a scenario in which only the low-v kick tail of the NS population survives in binaries.
Curiously, the most massive NS in the sample, Gaia NS1 (J1432-1021), has the lowest-eccentricity orbit, requiring the lowest kicks.This appears inconsistent with SNe simulations, which predict high-mass NSs to receive the strongest kicks (Müller et al. 2019;Burrows et al. 2023).It is possible that Gaia NS1 formed through a different channel than the other objects in the sample, but larger sample is required to investigate this possibility quantitatively.6.6.Future evolution to symbiotic X-ray binaries When the MS stars in our sample leave the main sequence, they will begin transferring mass to the dark companions: first by winds, and eventually by Roche lobe overflow.Rodriguez et al. (2024) have calculated 4 An accurate estimate of the space density of wide NS+MS binaries is presently difficult to calculate due to uncertainties in the Gaia selection function.However, the nearest NS candidate in our sample is at d ≈ 250 pc.The can be compared to d ≈ 120 pc, the distance to the nearest known young NS (Walter et al. 2010), and d ≈ 20 pc, a predicted distance to the nearest NS based on the Galactic model of Sweeney et al. (2022).These distances suggest that there are ∼ 10 3 dormant NSs for every one in a binary like the objects in our sample.models representative of this evolution.Their calculations suggest that the binaries will be X-ray bright for at least ∼ 5 Myr during Roche lobe overflow, and most likely for a few tens of Myr due to wind accretion at earlier times.Since the X-ray bright phase is predicted to be 100-1000 times shorter-lived than the current X-ray faint phase, one expects X-ray bright systems to be rare, but they are discoverable to large distances.
At least a few neutron stars with solar-mass red giant companions are known to exist in symbiotic X-ray binaries (e.g.Hinkle et al. 2006Hinkle et al. , 2019)).One system that is likely representative of our candidates' future evolution is GX 1+4 (Davidsen et al. 1977;Hinkle et al. 2006), which contains a ≲ 1 M ⊙ giant in a 1161-day orbit at a distance of ≈ 4 kpc.Another similar system is the symbiotic Xray binary IGR J16194-2810 (Masetti et al. 2007), which contains a ∼ 1 M ⊙ red giant at a distance of ≈ 2.1 kpc in a 193-day orbit with a NS companion (Nagarajan et al. 2024;Hinkle et al. 2024).Given their X-ray properties, these systems undoubtedly contain NSs.This provides some support for the interpretation that the dark objects in our binaries are also NSs, despite the uncertainties associated with forming NS+MS binaries in such wide orbits.6.7.Is there a BH/NS mass gap? Figure 17 shows the mass distribution of dark companions identified thus far from Gaia astrometry: the 21 NS candidates presented here, and three BHs studied previously (El-Badry et al. 2023a,b;Nagarajan et al. 2023;Gaia Collaboration et al. 2024).All the NS candidates have best-fit masses between 1.25 and 2.0 M ⊙ , while the BHs have masses above 9 M ⊙ .We have not detected any dark companions with intermediate masses, and our follow-up has been complete for sources published in DR3 with astrometric solutions implying M 2 ≳ 2 M ⊙ .The fact that we detected 21 dark companions with M 2 < 2 M ⊙ -which produce smaller astrometric wobbles than would low-mass BHs -strongly suggests that our search would have detected 3 − 5 M ⊙ BHs if they existed in our search sample.
Figure 18 compares the masses and orbital periods of BH and NS binaries identified with Gaia astrometry to BHs and NSs found with other methods.Yellow points show NSs in radio pulsars binaries ( Özel & Freire 2016;Fonseca et al. 2021).The hollow yellow symbol shows the companion to PSR J0514-4002E, whose nature is uncertain (Barr et al. 2024).Red points show highand low-mass BH X-ray binaries (Remillard & McClintock 2006;Corral-Santana et al. 2016;Miller-Jones et al. 2021).Cyan lower limits show spectroscopic BH binary candidates (Giesers et al. 2019;Shenar et al. 2022;Mahy et al. 2022).The BH and NS samples are both far from complete.Nevertheless, Figure 18 shows rather unambiguously that the BH/NS mass distribution is bimodal over at least 4 orders of magnitude in orbital period.
Our results thus support the presence of mass gap between NSs and BHs, similar to the gap reported for X-ray binaries (Bailyn et al. 1998;Kreidberg et al. 2012) and gravitational wave sources (Abbott et al. 2023).The existence of a mass bimodality does not depend much on the nature and exact mass of objects proposed to be in the mass gap (e.g.Casares et al. 2022;Barr et al. 2024), as such objects are rare, at least in binaries.Of course, it is possible that the mass distribution of BHs in binaries is different from the mass distribution of all BHs, since (for example) low-mass BHs could experience stronger kicks and more frequently be unbound (Burrows et al. 2023).

CONCLUSIONS
The first set of orbital solutions from the Gaia mission led to the identification of more than 50 candidate neutron star (NS) + main sequence (MS) binaries in auscale orbits.We are carrying out a spectroscopic followup program that yields radial velocities (RVs) and stellar parameters for these binaries.We have now covered a majority of an orbital period for most of these candidates with RVs.This allows us to tighten constraints on orbits through joint fitting of RVs and astrometry to and root out spurious astrometric orbits when the predictions of the Gaia solution do not match RVs (Appendix A).Our main results are as follows: 1. Summary of the sample: We have constructed a sample of 21 candidate NS + MS binaries containing solar-type MS stars and dark companions with best-fit masses of 1.25 − 1.90 M ⊙ (Figure 1).Most systems have periods of 400 to 1000 days and distances of 0.4 to 1 kpc.We select candidates as those with photocenter wobbles too large to be explained by a normal luminous companion or an unresolved inner binary containing two luminous stars.At fixed period and luminous star mass M ⋆ , NS companions produce photocenter orbits that are larger than WDs and non-degenerate companions can produce, but smaller than those produced by BH companions (Figure 2).
2. Quality of the Gaia orbital solutions: RVs for all the objects in our final sample are in good agreement with predictions of the Gaia orbital solutions (Figure 5), and joint fitting of RVs and astrometry leads to well-constrained orbits (Figure 6).For most objects, constraints from the Gaia-only and Gaia+RVs solutions are consistent within 1σ (Figure 7).However, some initial candidates were removed from the sample after RV follow-up revealed tensions with the Gaia solution (Figure 19 and Table 3).Most but not all of these objects have large goodness of fit, indicative of a problematic astrometric solution.On average, companion masses constrained by joint RV+astrometry fits are lower than the pure-astrometry estimates (Figure 8); this is an expected consequence of selecting candidates from the upper tail of the companion mass distribution.The objects all have unproblematic Gaia quality flags and astrometric uncertainties typical of their apparent magnitudes (Figure 10).Most of the orbits are predicted to have been sampled uniformly in phase by Gaia (Figure 9).

Eccentricities:
Most of the binaries in our sample have fairly eccentric orbits, with a median eccentricity of 0.4.Their eccentricities are significantly higher than those of typical WD companions at similar periods (Figure 12), and are consistent  et al. 2023a,b).A third BH has significantly higher mass, M dark ≈ 33 M ⊙ (Gaia Collaboration et al. 2024).No dark companions with mass between 2 and 8 M ⊙ have been detected, even though the effective search volume for such sources is larger than for the NSs.-Masses and orbital period of Galactic BHs and NSs in binaries with well-constrained masses.The candidates presented here are shown in magenta and are compared to NSs in pulsar binaries (yellow), BHs in astrometric binaries (black), BHs in X-ray binaries (red), and BHs in spectroscopic binaries (cyan).A mass bimodality exists over at least 4 orders of magnitude in orbital period.
with eccentricities expected due to kicks during NS formation (Figure 16).However, high-mass WDs (M ≳ 1 M ⊙ ) also have higher eccentricities than ∼ 0.6 M ⊙ WDs, so the eccentricities of our candidates do not guarantee that they are NSs.

Metallicities and chemical abundances:
We measured metallicities for the MS stars from highresolution spectra (Figure 3).A majority of them have metallicities near solar (−0.5 < [Fe/H] < 0.5), but low metallicity halo stars are over-represented.Three out of 21 targets (14%) have [Fe/H] ∼ −1.5 and space velocities of ∼ 300 km s −1 (Figure 13).Only 0.5% of all binaries with Gaia astrometric solutions at comparable distance are on halo orbits.This suggests that low-metallicity massive stars are more likely to form NSs in wide orbits with lowmass companions.
All three low-metallicity stars are strongly enhanced in lithium (Figure 14 and 15).Lowmetallicity stars have much thinner convective envelopes than solar-metallicities stars of the same mass and evolutionary state, so accreted material will be diluted less and abundance anomalies will be more detectable at low metallicity.The origin of the lithium is unclear: one possibility is accretion of Li-rich winds from super-AGB stars (Cameron & Fowler 1971;Ventura & D'Antona 2010).However, at least one the Li-enhanced objects has a companion with M 2 = 1.90 ± 0.03 M ⊙ , which is considerably higher than the expected mass of high-mass WDs or NSs formed from electron capture SNe.
5. Nature of the companions: We have not detected radiation from the companions, but constraints on their masses and orbits can narrow down the possibilities.Joint fitting of astrometry and RVs places tight upper limits on the flux ratio, with a median 2σ upper limit of ϵ < 0.03.This rules out all plausible nondegenerate companions (Figure 11).Several different possibilities remain, including single NSs, single massive WDs, and tight WD+WD, WD+NS, or WD+MS binaries.For most of the companions to be single massive WDs, their masses would have to be overestimated at the few-σ level.Scenarios involving inner binaries are difficult to explain with evolutionary models -triple systems tend to become unstable during their evolution when the outer orbit is as tight as our candidates -but should not be dismissed entirely for want of imagination.It is quite possible that the sample contains more than one kind of dark companion.

BH/NS mass distribution:
The Gaia mission has now enabled astrometric discovery of three BHs, 21 candidate NSs, and thousands of WDs in astrometric binaries.The mass distribution of BHs and NS candidates is shown in Figure 17.There is a conspicuous gap in the mass distribution between 2 and 8 M ⊙ .This may be in part a result of small number statistics, but the distribution does seem to disfavor a population of lower-mass BHs that significantly outnumber higher-mass BHs.Since the astrometric search volume is larger for higher-mass objects, our sample suggests that low-mass BHs are significantly less common in astrometric binaries than are NSs.When combined with samples of BHs and NSs discovered with other methods, our sample reveals a clear mass bimodality that extends over 4 orders of magnitude in orbital period (Figure 18).
It will likely remain difficult to conclusively establish the nature of the companions for some time.Radio detection of the companions as pulsars could prove that they are NSs, and efforts are underway to search for radio pulsations from most of our candidates.However, radio detection of any individual candidate seems a priori unlikely: young NSs are only detectable as radio pulsars for ≲ 10 Myr -only ∼ 0.1% of the expected lifetime of these binaries -and in the simplest evolutionary scenarios for their formation, there is no reason to expect the NSs to be recycled.Given the wide separations of the binaries and weak winds of the main-sequence companions, an X-ray detection due to accretion is also not expected (e.g.Rodriguez et al. 2024).WD companions are not expected to be detectable at any wavelength unless they are fortuitously young.Models involving inner binaries can be tested with high-precision RV observations (Nagarajan et al. 2023), but only if the inner binary has a period longer than a few days.
We briefly discuss two possible avenues for determining the nature of the companions in the absence of a direct detection.First, higher-precision astrometric orbits from future Gaia data releases will tighten mass constraints.The companion mass uncertainties are in most cases limited by the astrometric inclination uncertainties, and more precise measurements could more firmly rule out the possibility that the companions are ultramassive WDs.Future Gaia releases will also include epoch astrometry for all these sources, which will allow for further checks of the consistency between the astrometry and RV data.
Second, ongoing RV follow-up of astrometric binaries containing WDs will provide a more complete mapping of the period -eccentricity relation for post-interaction binaries containing WDs.This relation is particularly uncertain for massive WDs in wide orbits, making it difficult to distinguish between WDs and NSs on the basis of their eccentricities.Features in the relation near the WD-NS transition mass may make it possible to statistically differentiate between WD and NS scenarios, even if the nature of companions in individual systems remain uncertain.tained at Keck Observatory, which is a private 501(c)3 non-profit organization operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration.The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the Native Hawaiian community.We are most fortunate to have the opportunity to conduct observations from this mountain.

REJECTION OF SPURIOUS SOLUTIONS AND
SUMMARY OF ALL FOLLOW-UP For some candidates, our RV follow-up showed the Gaia astrometric solution to be spurious or to have significantly underestimated uncertainties.Measured and predicted RVs for three such sources are shown in Figure 19.The source shown in the top panel, with DR3 source id 6593763230249162112, was characterized as a high-quality NS candidate by both Andrews et al. (2022) and Shahaf et al. (2023b).The Gaia orbital solution has goodness of fit = 1.67, indicative of an unproblematic astrometric solution.Its DR3 astrometric solution is based on 20 visibility periods, which is typical for solutions published in DR3.There does not appear to be any reason to mistrust the solution based on the quality flags published in DR3.Yet, comparison of the observed and predicted RVs leaves little doubt that the astrometric solution is seriously in error.
The same is true for source 3869650535947137920, which is shown in the middle panel of Figure 19 and is included in the sample from Shahaf et al. (2023b).The source's goodness of fit is 5.46, which is quite normal for a source with G < 13 (Figure 10).The first several months of our follow-up showed the RVs to be evolving in a manner consistent with the astrometric solution's predictions.However, the RVs began to diverge from those predictions in the second season of our follow-up.Although the qualitatively similar shape of the predicted and observed RV curves suggests the astrometric solution is not completely spurious, the observed degree of disagreement suggests that the astrometric uncertainties are significantly underestimated.
Finally, source 747174436620510976, shown in the bottom panel of Figure 19, is an example of a case where the shape of the observed RV curve is consistent with the astrometric predictions, but there is a few-σ offset in phase.Our analysis of the source's SED also suggested that the parallax is overestimated, leading us to exclude the source (Section 4).The Gaia solution has goodness of fit = −0.96,indicative of a good solution, and was included in both the Shahaf et al. (2023b) and Andrews et al. (2022)

candidate lists.
These examples highlight the importance of RV followup: while the sources shown in Figure 19 clearly are binaries, and least two of the three likely have orbital parameters not too far from those inferred from astrometry alone, their astrometric uncertainties are unlikely to be reliable.Since most NSs have masses near the Chandrasekhar limit, small problems with the astrometric solution can seriously change our conclusions about the nature of a candidate.Epoch astrometry from Gaia DR4 may allow us to obtain more useful astrometric constraints from sources like those shown in Figure 19.
In Table 3, we provide a summary of our follow-up of all NS candidates from Shahaf et al. (2023b) and Andrews et al. (2022).For the Shahaf et al. (2023b) sample, we list candidates for which they inferred M 2 > 1.25 M ⊙ .The Andrews et al. (2022) sample only includes targets for which they infer M 2 > 1.4 M ⊙ .A majority of the candidates ruled out by RV follow-up have significantly higher goodness of fit than typical sources of the same apparent magnitude (Figure 10).Several candidates with incomplete RV follow-up, particularly those with G < 15 that are inaccessible with FEROS and TRES, remain to be studied.

BLUE EXCESS DUE TO WD COMPANIONS
As can be seen in Figure 1, several of our candidates are at the blue edge of the main sequence.Compared to stars near the middle of the main sequence, this amounts to a blue color excess of 0.05 to 0.1 mag.Most of our candidates have also been observed by GALEX, with NUV detections that rule out significant (≳ 0.2 mag) UV excess.Here we investigate whether the optical blue excess could be due to a WD companion.
We first consider a 1 M ⊙ MS star in a binary with a single 1.3 M ⊙ WD companion.We construct the combined spectral energy distribution of the pair, modeling the WD with a Koester DA model with log(g/cm s −2 ) = 9.0   -Examples of NS+MS candidates we rejected after RV follow-up.In the top panel, the RVs are grossly inconsistent with the predictions of the Gaia solution, suggesting that the solution is spurious, or at least has a large accumulated phase error not captured in the Gaia uncertainties.Middle panel shows an example in which the RVs evolve in a manner somewhat consistent with the Gaia solution prediction, but do not match it quantitatively.This indicates that the uncertainties of the astrometric parameters are likely underestimated, casting doubt on the companion mass constraint.Bottom panel shows a marginal case: the RVs qualitatively track the astrometric solutions predictions, but a phase offset suggests that the period uncertainty is underestimated.We exclude all these objects -even though they may indeed host NSs or massive WDs -retaining only sources with good agreement between astrometry and RVs.(Tremblay & Bergeron 2009;Koester 2010) and the MS star with the same models used in Section 4. We then calculate the combined total magnitude of both objects using pyphot.The solid lines in Figure 20 show the resulting color excess as a function of WD effective temperature.We define the color excess in each band as the color of the combined MS + WD pair minus the color of the MS star alone.
Because a 1.3 M ⊙ WD is very small, the predicted blue excess in the optical is negligible except at very high temperatures.On the other hand, effective temperatures above ≈ 25, 000 K would lead to a detectable FUV source (here we assume a distance of 500 pc), which represents a strong UV excess easily distinguishable from a single MS star.A weaker, but still significant, NUV excess is also predicted.On the basis of these calculations, we conclude that -given their high dynamical masses -there is no plausible way for single WD companions to significantly change the optical colors of our candidates without simultaneously causing a large UV excess, which is not observed.
Next, we consider the possibility of a close WD+WD binary of total mass 1.3 M ⊙ .In this case, a broader range of combined SEDs are possible, depending on the mass ratio and cooling ages of both WDs.For simplicity, we consider the case where both WDs have the same mass and effective temperature.In this case, the blue excess is significantly stronger at fixed T eff , reflecting the fact that two 0.65 M ⊙ WDs have ≈ 18 times the total surface area of a single 1.3 M ⊙ WD.
With two 0.65 M ⊙ WDs, a ≈ −0.04 mag blue excess is possible in the optical if both WDs have T eff ∼ 50, 000 K. However, this would be accompanied by an excess of several magnitudes in both the NUV and FUV bands, and is thus ruled out.The same is true for a 1.0 + 0.3 M ⊙ WD binary.In this case, the flux is dominated by the 0.3 M ⊙ WD, which could cause a significant optical blue excess for T eff ≳ 30, 000 K.Here too, the excess would be accompanied by a much larger excess in the UV.-Predicted color excess due to a WD companion.In all cases, we assume a 1 M ⊙ solar-type primary with a 1.3 M ⊙ companion, representative of the objects in our sample.Solid blue lines show the case where the companion is a single WD, while dashed black line shows an equal-mass WD+WD binary and dotted cyan lines show a (1.0 + 0.3) M ⊙ WD binary.Top panel shows the optical blue excess in the Gaia G BP −G RP bands, while the middle and bottom panels show the UV excess in the GALEX NUV and FUV bands.Lighter lines in the bottom panel show cases where the system would be fainter than 20.0 mag in the FUV band at a distance of 500 pc and would likely not be detected.WD+WD binaries lead to much more blue and UV excess than single WD companions of the same mass.Even WD companions that are quite hot produce negligible blue excess in the optical.

J2102+3703;
Gaia DR3 ID 1871419337958702720  This paper was built using the Open Journal of Astrophysics L A T E X template.The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike.Learn more at http://astro.theoj.org.

RADIAL VELOCITIES
Fig. 1.-Black points show all binaries with astrometric orbital solutions published in Gaia DR3.Red points show the 21 objects presented in this work, which have astrometrically-inferred companion masses M 2 > 1.25 M ⊙ .Upper left: dereddened color-magnitude diagram.Most candidates are solar-type stars near the main-sequence.Upper right: orbital period and distance.Candidates have orbital periods of ∼ (100 − 1000) d and are within ∼ 1 kpc of the Sun.Lower left: Galactic coordinates, with the Galactic center in the middle of the plot.Imprints of the Gaia scanning law are visible; most NS candidates are at high latitude.Bottom right: luminous star mass and minimum mass ratio, M 2 /M⋆.The NS candidates have some of the highest estimated mass ratios in the astrometric binary sample.
Fig. 2.-Comparison of NS candidates presented here (red stars)to the rest of the Gaia DR3 astrometric binary sample.At fixed period, dark and massive companions produce larger photocenter orbits than normal stellar companions.For typical solar-type primaries, NS companions produce photocenter orbits that are smaller than those of BH binaries, but larger than those of luminous binaries and triples.

Fig. 4 .
Fig. 4.-Spectral energy distributions and best-fit single-star models.The parameters of the best-fit model are listed in each panel.UV photometry is plotted where available, but not included in the fit.All targets' SEDs are reasonably well-fit by a single-star model.Parameters from SED fitting are listed in Table 2.

Fig
Fig.5.-Comparison of our measured RVs (red) to predictions of the Gaia astrometric solutions (cyan).Individual lines show random draws from the posterior.For AstroSpectroSB1 solutions, the systemic RV is set to the value measured by Gaia; for Orbital solutions, it is set to the best-fit value from our joint astrometry+RV fits (Figure6).Systems with long periods (P orb ≳ 700 days) have significant uncertainties in their Gaia-only periods and phases at the time of our observations.The RVs are broadly consistent with the Gaia predictions.
Fig.6.-Jointfits of our RVs (red) and the Gaia constraints.Lower sub-panels show residuals compared to the best-fit solution.Inclusion of RVs in the fit yields much tighter constraints than those from Gaia alone (Figure5).The RV residuals are generally consistent with 0 within their (1 − 2)σ uncertainties.

N
Fig.9.-Sky-projected orbits of the luminous stars.The angular orbits have been multiplied by distance, so that the displayed orbital sizes correspond to relative physical sizes.Red points mark Gaia observations predicted by GOST, showing that all of the targets in our sample are expected to have good astrometric phase coverage.Note that we do not have access to the actual measured ∆RA and ∆Dec values; only to the predicted scan times.Red text in each panel indicates the number of visibility periods that are predicted by GOST, as well as the number actually used in calculating the astrometric solution.
Fig. 10.-Comparison of our NS candidates (cyan/red) to all sources with astrometric orbital solutions published in DR3 (black).Hollow symbols show three initial candidates we found to have spurious solutions (Appendix A).Upper left: astrometric goodness of fit.Large values indicate a poor formal fit.There is a discontinuity at G = 13 owing to a quirk of the Gaia data processing: brighter sources typically have larger goodness of fit values because the uncertainties in their individual-epoch astrometric data are underestimated.Lower left: parallax uncertainty.Cyan points show the uncertainty reported in Gaia DR3, while red points show the uncertainty from our joint RV + astrometry fit.Right panels: parallax over error (top) and photocenter semimajor axis over error (bottom) as a function of orbital period.The diagonal "cliff" is a result of quality cuts imposed on the published solutions.The NS candidates have typical goodness of fit values and astrometric uncertainties for their apparent magnitudes and periods, indicative of unproblematic astrometric solutions.The sources with spurious solutions have slightly larger parallax uncertainties than average but otherwise do not stand out from the sources with reliable solutions.
Fig. 11.-Constraints on the mass of the unseen companion in J2244-2236, a typical system, as a function of the G−band flux ratio.Solid black line shows the constraint from the astrometric mass function (Equations 20), assuming M⋆ = 1.00 M ⊙ ; a completely dark companion would imply M 2 ≈ 1.44 M ⊙ .If the companion contributes some light, its astrometrically-implied mass increases.Dotted cyan line shows the expected flux ratio and mass for a single MS companion.Because this is always below the black line, no single MS companion can explain the orbit.The same is true for an equal-mass inner binary (yellow dashed) or an equal-mass inner triple (red dot-dashed).An inner binary containing a 1 M ⊙ WD and a ≈ 0.6 M ⊙ MS star (blue dashed) could explain the astrometric mass function with the WD contributing 4% of the light in the G−band, but this violates the flux ratio constraint from RVs (see text).A tight WD+WD binary (not shown, because it could have F 2 /Ftot ≈ 0), could explain the orbit, but it is unclear how such a system could form.
Fig. 12.-Left: Period-eccentricity diagram for NS+MS candidates from this work (black) compared to WD+MS binary candidates selected byShahaf et al. (2023b) from Gaia DR3 on the basis of their astrometrically-constrained mass ratios.For NS+MS candidates, uncertainties are smaller than symbols.Most of the WD candidates have low eccentricities, probably resulting from (incomplete) tidal circularization when the WD progenitor was a red giant.Our candidates have significantly higher eccentricities than the WD candidates at all periods, which may be the result of kicks during NS formation.Right: mass-eccentricity diagram.Our candidates are all selected to have M 2 > 1.25 M ⊙ .Solid and dashed lines show median and middle 68% across both red and black points.Binaries with M 2 = (1.0 − 1.2) M ⊙ , most of which likely contain high-mass WDs, also have higher eccentricities than systems containing lower-mass WDs.

Fig. 13 .
Fig. 13.-Toomre diagram of NS candidate binaries compared to all binaries with AstroSpectroSB1 solutions published in Gaia DR3.In the left panels, black points show the full DR3 catalog.In the right panels, they show only binaries with distances 0.4-1 kpc, similar to the NS candidates.Points are colored by metallicity.Dashed lines show total velocities of 70 and 180 km s −1 with respect to the local standard of rest, approximately separating the thin disk, thick disk, and halo.The three candidates with halo orbits (labeled) all have low metallicities.Halo and thick-disk orbits are significantly over-represented in the NS candidate sample.Natal kicks may be responsible for kinematically heating the orbits of these objects somewhat, but our results imply that NS candidates are over-represented in low-metallicity populations.
Fig. 14.-Evidence for lithium enhancement in the three halo stars.Red lines show the spectra of the three NS candidates.Black lines show spectra of three stars observed by the GALAH survey with similar physical parameters to each target and spectra that closely match those of the target in most lines.All three of our targets have unusually strong Li I λ6708 lines, indicative of lithium enhancement.
Fig. 15.-Lithium abundances of the three metal-poor halo stars in our sample (Figure 14), compared to metal-poor stars observed by the GALAH survey.All three stars are strongly enhanced in lithium, with A(Li) higher than any of the ∼ 10 3 stars with similar stellar parameters observed by GALAH.
and median predicted eccentricity (bottom) for wide helium star + MS binaries when the He star explodes and leaves behind a NS.We assume the initial orbit is circular and consider initial helium star masses of 3, 4, and 5 M ⊙ , and initial periods of 100, 300, and 1000 days.In all cases, we model the companion as a 1 M ⊙ MS star and adopt a final NS mass of 1.4 M ⊙ .The NS is born with a velocity v kick in a random direction.For a 3 M ⊙ He star, orbits are most likely to remain bound if the kick is slow.For M He star > 3.8 M ⊙ , the binary will be unbound by mass loss unless a fortuitously aligned kick allows the NS to catch the escaping companion.

Fig
Fig. 17.-Mass distribution of dark companions revealed by follow-up of Gaia astrometric binaries.The 21 sources with masses M dark = 1.25 − 1.9 M ⊙ are presumed NSs introduced in this paper.The two objects with masses near 9 M ⊙ are presumed BHs (El-Badryet al. 2023a,b).A third BH has significantly higher mass, M dark ≈ 33 M ⊙ (GaiaCollaboration et al. 2024).No dark companions with mass between 2 and 8 M ⊙ have been detected, even though the effective search volume for such sources is larger than for the NSs.
Fig.18.-Masses and orbital period of Galactic BHs and NSs in binaries with well-constrained masses.The candidates presented here are shown in magenta and are compared to NSs in pulsar binaries (yellow), BHs in astrometric binaries (black), BHs in X-ray binaries (red), and BHs in spectroscopic binaries (cyan).A mass bimodality exists over at least 4 orders of magnitude in orbital period.

P
Fig. 19.-Examples of NS+MS candidates we rejected after RV follow-up.In the top panel, the RVs are grossly inconsistent with the predictions of the Gaia solution, suggesting that the solution is spurious, or at least has a large accumulated phase error not captured in the Gaia uncertainties.Middle panel shows an example in which the RVs evolve in a manner somewhat consistent with the Gaia solution prediction, but do not match it quantitatively.This indicates that the uncertainties of the astrometric parameters are likely underestimated, casting doubt on the companion mass constraint.Bottom panel shows a marginal case: the RVs qualitatively track the astrometric solutions predictions, but a phase offset suggests that the period uncertainty is underestimated.We exclude all these objects -even though they may indeed host NSs or massive WDs -retaining only sources with good agreement between astrometry and RVs.

Fig
Fig. 20.-Predicted color excess due to a WD companion.In all cases, we assume a 1 M ⊙ solar-type primary with a 1.3 M ⊙ companion, representative of the objects in our sample.Solid blue lines show the case where the companion is a single WD, while dashed black line shows an equal-mass WD+WD binary and dotted cyan lines show a (1.0 + 0.3) M ⊙ WD binary.Top panel shows the optical blue excess in the Gaia G BP −G RP bands, while the middle and bottom panels show the UV excess in the GALEX NUV and FUV bands.Lighter lines in the bottom panel show cases where the system would be fainter than 20.0 mag in the FUV band at a distance of 500 pc and would likely not be detected.WD+WD binaries lead to much more blue and UV excess than single WD companions of the same mass.Even WD companions that are quite hot produce negligible blue excess in the optical.

TABLE 3
All the NS candidates fromAndrews et al. (2022, A22)andShahaf et al. (2023b, S23).We list the orbital period according to the Gaia astrometric solution, the apparent magnitude and goodness of fit, and the status of our follow-up.Candidates for which the status is colored in green are modeled in detail in this work.Those for which status is colored in red are disfavored by RV follow-up.Blue text indicates that detailed follow-up is presented elsewhere, and black text indicates that insufficient RVs have been obtained for a verdict to be reached.
Table 4 lists constraints from joint fitting of RVs and Gaia astrometry for all sources.

TABLE 4
Orbit fitting results for all candidates.

TABLE 5
Table 5 lists all the RVs used in our analysis.TRES Radial velocities for all targets.A machine-readable version of the table is included as supplemental material.