A 1 . 9 M ⊙ NEUTRON STAR CANDIDATE IN A 2-YEAR ORBIT

We report discovery and characterization of a main-sequence G star orbiting a dark object with mass 1 . 90 ± 0 . 04 M ⊙ . The system was discovered via Gaia astrometry and has an orbital period of 731 days. We obtained multi-epoch RV follow-up over a period of 639 days, allowing us to refine the Gaia orbital solution and precisely constrain the masses of both components. The luminous star is a ≳ 12 Gyr-old, low-metallicity halo star near the main-sequence turnoff ( T eff ≈ 6000 K; log (cid:0) g/ (cid:2) cm s − 2 (cid:3)(cid:1) ≈ 4 . 0; [Fe / H] ≈ − 1 . 25; M ≈ 0 . 79 M ⊙ ) with a highly enhanced lithium abundance. The RV mass function sets a minimum companion mass for an edge-on orbit of M 2 > 1 . 67 M ⊙ , well above the Chandrasekhar limit. The Gaia inclination constraint, i = 68 . 7 ± 1 . 4 deg, then implies a companion mass of M 2 = 1 . 90 ± 0 . 04 M ⊙ . The companion is most likely a massive neutron star: the only viable alternative is two massive white dwarfs in a close binary, but this scenario is disfavored on evolutionary grounds. The system’s low eccentricity ( e = 0 . 122 ± 0 . 002) disfavors dynamical formation channels and implies that the neutron star likely formed with little mass loss ( ≲ 1 M ⊙ ) and with a weak natal kick ( v kick ≲ 20 km s − 1 ). Stronger kicks with more mass loss are not fully ruled out but would imply that a larger population of similar systems with higher eccentricities should exist. The current orbit is too small to have accommodated the neutron star progenitor as a red supergiant or super-AGB star. The simplest formation scenario – isolated binary evolution – requires the system to have survived unstable mass transfer and common envelope evolution with a donor-to-accretor mass ratio > 10. The system, which we call Gaia NS1, is likely a progenitor of symbiotic X-ray binaries and long-period millisecond pulsars. Its discovery challenges binary evolution models and bodes well for Gaia’s census of compact objects in wide binaries.

1. INTRODUCTION Astrometry from the Gaia mission (Gaia Collaboration et al. 2016) has opened a new window on the population of compact stellar remnants lurking in non-interacting binaries.Precise measurements of astrometric "wobble" over the course of several years allow Gaia to detect the presence of binary companions via their gravitational effects, even when they do not emit any light (Gaia Collaboration et al. 2023b).Unlike most other binary detection methods, astrometry is most sensitive to companions that are dark and in wide orbits (P orb ∼ 1000 d for the currently published data).
Astrometric orbital solutions from the mission's third data release (DR3; Gaia Collaboration et al. 2023a) have already enabled the discovery of thousands of white Corresponding author: kelbadry@caltech.edudwarfs (WDs; Shahaf et al. 2023b;Yamaguchi et al. 2023;Shahaf et al. 2023a) and two black holes (BHs; El-Badry et al. 2023a;Chakrabarti et al. 2023;El-Badry et al. 2023b) orbited by luminous stars.The orbits of both the BH and WD systems are unexpectedly wide, falling in a period range that binary evolution models predict to be sparsely populated.BH and WD companions in au-scale orbits appear not to be rare compared to closer binaries.They were under-represented in samples before Gaia because discoveries based on photometric variability, radial velocity (RV) shifts, and signatures of accretion all favor short-period systems.
Unlike BHs and WDs, neutron star (NS) companions have not yet been unambiguously identified from the Gaia data.Several factors make astrometric identification of NSs difficult.The mass distribution of wellcharacterized NSs is peaked near 1.3 M ⊙ ( Özel et al. 2012;Kiziltan et al. 2013;Özel & Freire 2016), below the maximum WD mass (∼ 1.38 M ⊙ ;Nomoto 1987).This means that a majority of NSs cannot be distinguished unambiguously from high-mass WDs on the basis of their mass alone.Even for NSs with masses above the maximum WD mass, observational uncertainties generally make mass constraints from Gaia data alone consistent with being below the maximum WD mass at the (1−3)σ level (e.g., Shahaf et al. 2023b).Finally, most NSs are thought to form with strong birth kicks due to asymmetric supernova explosions (e.g., Hobbs et al. 2005;Janka 2012).These kicks -coupled with a sudden drop in binaries' binding energy due to mass loss -likely disrupt a large majority of would-be NS binaries.It is probably largely for this reason that fewer than 1% of all known young pulsars are in binaries (e.g., Lorimer 2008), even though most of the massive stars from which these pulsars form are in binaries and higher-order multiples (e.g., Sana et al. 2012).
After the release of Gaia DR3, we initiated an RV follow-up program targeting ∼ 50 binaries with astrometric solutions suggesting the companion might be a NS.RV follow-up of these candidates is ongoing and described in a companion paper, El-Badry et al. (2024, in prep).In this paper, we present a detailed analysis of one of the most interesting candidates from that sample, which we refer to as Gaia NS1 (Gaia DR3 source ID 6328149636482597888).Besides Gaia NS1, our followsample includes 20 other high-quality NS candidates.
Compared to other candidates we have followed-up, Gaia NS1 is unique in three respects: (a) it is the only system for which RVs alone set a minimum mass for the dark companion that is unambiguously above the Chandrasekhar mass, (b) it has the highest inferred dark companion mass of all NS candidates we have followed up, and (c) it has the lowest eccentricity among the NS candidates.These features make it unlikely that the unseen companion is anything other than a NS, and they make the system's properties particularly difficult to explain with binary evolution models.
The remainder of this paper is organized as follows.Section 2 describes how the system was identified as a compact object candidate, while Section 3 summarizes the source's basic observational properties.Section 4 describes our spectroscopic observations and measurement of RVs, while Section 5 describes fitting of the RVs and Gaia data to constrain the companion mass.In Section 6, we measure atmospheric parameters and abundances of the luminous star.Section 7 examines the system's Galactic orbit.We discuss the nature of the companion and the system's possible formation histories in Section 8 and conclude in Section 9.
2. DISCOVERY Gaia NS1 was identified as a compact object candidate soon after Gaia DR3 by several authors.Andrews et al. (2022) assumed a luminous star mass of M ⋆ = 1.21 ± 0.2 M ⊙ and calculated a dark companion mass of M 2 = 2.71 +1.50  −0.36 M ⊙ (2σ uncertainties) based on the astrometric orbit.They classified the unseen companion as a BH.Shahaf et al. (2023b) assumed M ⋆ = 1.09M ⊙ and inferred M 2 = 2.45 ± 0.2 M ⊙ , also classifying the companion as a BH.The luminous star mass estimate adopted by Shahaf et al. (2023b) was from the gaiadr3.binarymasses table (Gaia Collaboration et al. 2023b), while the estimate adopted by Andrews et al. (2022) was from the Gaia FLAME estimator (Creevey et al. 2023).Both values are based on comparison of the source's color and absolute magnitude to isochrones, with a metallicity prior that favors solar metallicity.El-Badry et al. (2023a) reported a spectroscopic metallicity of [Fe/H] ≈ −1.5 for the luminous star and a revised mass estimate of only M ⋆ = 0.78 ± 0.05 based on a metallicity-informed fit of the spectral energy distribution.This significantly lower mass estimate reflects the facts that (a) lower-metallicity stars are brighter at fixed mass, and (b) the star's CMD position implies that it is near the main-sequence turnoff, having already expanded significantly since the zero age main sequence.Given this mass estimate for the luminous star, they estimated M 2 = 2.25 +0.61  −0.26 M ⊙ based on the Gaia astrometric solution.
El- Badry et al. (2023b) subsequently reported that early radial velocity (RV) measurements were roughly consistent with the Gaia astrometric solution but cautioned that RV coverage over a majority of the orbit would be required to validate or refute the astrometric solution.We have now obtained the necessary RVs.
3. OBSERVATIONAL PROPERTIES Gaia NS1 is a bright (G = 13.3)source in the Libra constellation (α = 14:32:20.7;δ = −10:21:59) at high Galactic latitude (l = 339.5,b = 45.2).The field is sparse, with no other Gaia sources within 12 arcsec.The Gaia DR3 astrometric solution has a proper motion of 92 mas yr −1 and a parallax of 1.23 ± 0.04 mas.This implies a tangential velocity of ≈ 350 km s −1 , characteristic of a halo star.On the color-magnitude diagram (Figure 1), the source falls on top of the solar neighborhood main sequence.It has a de-reddened absolute magnitude M G,0 ≈ 3.5, about twice as bright as the Sun.This led Andrews et al. (2022) and Shahaf et al. (2023b) to assume masses of M ⋆ ≈ 1.1 − 1.2 M ⊙ .As we show below, the object is a low-metallicity turn-off star with a somewhat lower mass.

Distance
The parallax constraint from Gaia NS1's DR3 astrometric binary solution is ϖ = 1.23 ± 0.04 mas, corresponding to a distance of d ≈ 812 ± 26 pc.However, the parallax is covariant with other parameters of the astrometric model, and our joint fitting of RVs and astrometric constraints in Section 5 leads to revised constraints on those parameters, thus also updating the parallax constraint.The parallax constraint from our joint fit is ϖ = 1.36 ± 0.03 mas, corresponding to a closer distance of d ≈ 735 pc.
We adopt the parallax from the joint fit and associated closer distance of 735 pc in most of our fiducial modeling.We also explore how a larger distance of ≈ 812 pc would change our results, finding that it would change the best-fit masses of both components by less than 2%.We discuss the modest tension between the Gaia-only and Gaia+RV parallaxes in Section 5.5.

Spectral energy distribution
We constructed the source's broadband spectral energy distribution (SED) by combining UV photometry from GALEX (Martin et al. 2005), synthetic SDSS ugriz photometry constructed from Gaia BP/RP spectra (Gaia Collaboration et al. 2022), 2MASS JHK photometry (Skrutskie et al. 2006), and WISE W 1 W 2 W 3 photometry (Wright et al. 2010).We set an uncertainty floor of 0.03 mag in all bands to account for photometric calibration issues and imperfect models.We then fit the SED with single-star models, with a prior on the metallicity and α-abundance from spectroscopy ([Fe/H] = −1.23±0.08 and [α/Fe] = 0.20, see Section 6).
We predict bandpass-integrated magnitudes using empirically calibrated model spectral from the BaSeL library (v2.2;Lejeune et al. 1997Lejeune et al. , 1998)).We assume a Cardelli et al. (1989) extinction law with R V = 3.1 and adopt a prior on the reddening E(B − V ) = 0.08 ± 0.02 based on the Green et al. (2019) 3D dust map.We use pystellibs1 to interpolate between model SEDs, and pyphot2 to calculate synthetic photometry.
We fit the SED using emcee (Foreman-Mackey et al. 2013) to sample from the posterior, with the initial mass, metallicity, age, parallax, and reddening sampled as free parameters.We use the MINESweeper framework (Cargile et al. 2020) to interpolate on MIST isochrones (Choi et al. 2016) to predict the radius, effective temperature, and surface metallicity from the initial mass, age, and metallicity.These calculations account for atomic diffusion, which causes the present-day surface metallicity to be somewhat lower than the star's average metallicity (e.g., Dotter et al. 2017).We set a maximum age of 13.5 Gyr.

Luminous star mass constraints
The results of SED fitting are shown in Figure 1 and reported in Table 3.We find a luminous star mass of M ⋆ = 0.79 ± 0.01 M ⊙ and an age of 12.5 ± 0.5 Gyr.The age posterior samples extend to the upper limit of the prior, with a 2σ lower-limit of 11.5 Gyr.Given the star's low metallicity and halo-like orbit (Section 7), an age of 11+ Gyr is expected (e.g., Xiang & Rix 2022).
The inferred parameters place the star at the end of its main-sequence evolution and imply that it will soon become a red giant, reaching the tip of the giant branch within ≲1 Gyr.The current radius, R ≈ 1.48 R ⊙ , is already a factor of ∼ 2.1 larger than the value predicted for the same mass at the zero-age main sequence.The MIST evolutionary models predict that atomic diffusion has significantly reduced the surface metallicity: while they predict a current surface iron abundance [Fe/H] ≈ −1.23, the corresponding initial metallicity is [Fe/H] ≈ −1.08.Accounting for the star's enhanced abundance of α elements ([α/Fe] ≈ +0.20; Section 6), this implies an initial bulk metallicity of [M/H] ≈ −0.95.
We show the luminous star's inferred temperature and radius in the lower left panel of Figure 1 with a solid red point.The hollow red point shows the parameters inferred when we adopt the Gaia-only parallax: in that case, the ≈ 10% larger distance results in a 10% larger inferred radius.We compare the measured parameters to MIST evolutionary tracks with initial [M/H] = −0.95.
Because the models predict rapid evolution near the main-sequence turnoff, the inferred mass and age when we adopt the Gaia-only parallax are only slightly different: M ⋆ = 0.81 ± 0.015 M ⊙ and τ = 12.0 ± 0.7 Gyr, still corresponding to a star near the main-sequence turnoff.The tracks in Figure 1 also demonstrate that -even if our adopted distance were grossly in error -there is no plausible scenario under single-star evolution in which the mass of the luminous star is below 0.75 M ⊙ : models with M < 0.75 M ⊙ do not reach T eff ≈ 6000 K at any age below 13.5 Gyr.
There is good agreement between the inferred mass, radius, surface gravity, and distance.We discuss whether binary evolution could have produced a lower-mass luminous star with the observed parameters in Section 8.1, concluding that such a scenario is very unlikely.

SPECTROSCOPIC OBSERVATIONS
4.1.Data The promise of Gaia NS1 as a NS or BH binary prompted follow-up at several telescopes and by several groups, which are synthesized into a coherent analysis here.

FEROS
We observed Gaia NS1 19 times using 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, and P112.A-6010).The spectra have resolution R ≈ 50, 000 and cover 360-920 nm.Most of our observations used 2400 s exposures.We reduced the data using the CERES pipeline (Brahm et al. 2017), which performs bias-subtraction, flat fielding, wavelength calibration, 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.For RV measurements (Section 4.2), we used 15 orders with wavelengths between 4500 and 6700 Å.

TRES
We observed Gaia NS1 13 times 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).We measured RVs using 31 orders covering the wavelength range of 420-670 nm.

APF
We observed Gaia NS1 14 times using the Levy spectrometer on the 2.4 m Automated Planet Finder Telescope (APF; Radovan et al. 2010;Vogt et al. 2014) at Lick Observatory.We used a 2 ′′ slit, yielding spectra with resolution R ≈ 80, 000 and coverage over 373-1020 nm.Our observations did not use the iodine cell.We reduced the spectra using the California Planet Search pipeline (Howard et al. 2010;Fulton et al. 2015) and measured RVs using 35 orders orders covering the wavelength range of 440-670 nm.Lower left: comparison of the measured temperature and radius to MIST single-star evolutionary models.Solid point shows the radius inferred when adopting the parallax from our joint Gaia + RV fit; hollow point shows the parallax from Gaia alone (see Section 3.1).The observed temperature, radius, and metallicity imply a luminous star mass of 0.79 ± 0.01 M ⊙ .Lower right: dynamical constraints on the companion mass.The RVs and luminous star mass constrain the companion mass and inclination to the gray shaded region, with a minimum companion mass of 1.67 M ⊙ for an edge-on orbit.The inclination constraint from Gaia astrometry then implies M 2 = 1.90 M ⊙ .

MIKE
We observed Gaia NS1 10 times 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 for 3 observations, and the 0.7 ′′ slit for the other 6.Exposure times ranged from 300 to 2160s, yielding spectral resolution R ∼ 40, 000 on the blue side and R ∼ 55, 000 on the red side.The total wavelength coverage was 333-968 nm.We reduced the spectra with the MIKE pipeline within CarPy (Kelson et al. 2000;Kelson 2003).We measured RVs using 18 orders on the red side covering the wavelength range of 500-680 nm.

MagE
We observed Gaia NS1 twice with the Magellan Echellette spectrograph (MagE; Marshall et al. 2008) on the 6.5m Magellan Clay telescope at Las Campanas Observatory.Both observations were carried out with the 0.7 arcsec slit and a 300 s exposure.This yielded spectral resolution R ≈ 5, 400 and wavelength coverage of 350-1100 nm.We reduced the spectra using PypeIt (Prochaska et al.

2020).
4.1.6.PEPSI We obtained one spectrum 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 and the CD2 and CD5 cross-dispersers on the blue and red side, respectively, with an exposure time of 600 s.The spectrum was reduced as described in Strassmeier et al. (2018); it covers the wavelength range of 422-479 nm and 624-743 nm with spectral resolution R ≈ 50, 000.

Radial velocity measurements
We measured RVs from all the spectra by crosscorrelating the individual orders with a Kurucz template spectrum from the BOSZ library (Bohlin et al. 2017).Motivated by the spectroscopic analysis (Section 6), we used a template with T eff = 6000 K, log g = 4, and [Fe/H] = −1.25.For most instruments, we calculated uncertainties on the RVs from the standard deviation of  the RVs inferred from individual orders divided by the square root of the number of orders used in the analysis.
For the MagE spectra, where the uncertainty associated with the telluric RV correction is expected to be substantial due to lower resolution, we adopted conservative uncertainties of 3 km s −1 .
To correct for shifts in the wavelength solution due to instrumental flexure, slit centering errors, and other systematics, we cross-correlated a telluric model spectrum synthesized with TelFit (Gullikson et al. 2014) with the telluric "A" and "B" bands due to molecular oxygen in the observed spectra (Griffin 1973).We applied the thusinferred corrections to the RVs measured from the APF, MIKE, MagE, and PEPSI spectra, where they reached a maximum amplitude of a few km s −1 .We did not apply them to the FEROS and TRES spectra, where the bestfit shifts had typical amplitudes of a few tens of m s −1 , comparable to their uncertainties.We applied heliocentric corrections to all RVs after the telluric corrections.
All our measured RVs are listed in Table 1.The median uncertainty is 0.10 km s −1 .This is slightly larger than typical for our follow-up, due mainly to the star's low metallicity, but still 400+ times smaller than the star's RV variability amplitude, allowing the orbit to be tightly constrained.

ORBIT AND COMPANION MASS
To explore how the Gaia astrometry and our followup RVs separately constrain the companion mass, we present orbit fits that depend only on RVs (Section ??), on both astrometry and RVs (Section 5.2), and on astrometry alone (Section 5.3).This also allows us to assess the consistency between the astrometry and RVs and to consider the effects of underestimated astrometric uncertainties (Section 5.5).

TABLE 3
Physical parameters and 1σ uncertainties of Gaia NS1's orbit.We compare constraints on the orbit based on both Gaia and our RVs (1st block), the Gaia solution alone (2nd block), and our RVs alone (3rd and 4th block; in the 4th block, the period is fixed to the value inferred from the joint fit).Constraints with the Gaia uncertainties inflated by a factor of 3 are listed in the 5th block.We consider these the most robust constraints.

Parameters
5.1.Pure RV fit Neglecting the Gaia data entirely, we first fit the observed RVs with a Keplerian model.In addition to the standard 7 Keplerian parameters, we fit instrumental offsets for TRES, APF, and MIKE relative to FEROS.We do not fit offsets for MagE and PEPSI, where the small number of RVs would make them poorly constrained.We use emcee (Foreman-Mackey et al. 2013) to sample from the posterior, drawing 5000 samples with 64 walkers after a burn-in period of 5000 steps.We initialize the fit at the maximum-likelihood solution implied by the Gaia constraints assuming a dark companion and use broad, flat priors on all parameters.The RVs are densely sampled and cover most of a full orbit, so the posterior is unimodal.The results are reported in Table 3 under "RVs only".The measured period, P orb = 730.7 ± 2.2 d, and eccentricity, e = 0.122 ± 0.004, are fully consistent with the Gaia solution.The RV mass function is constrained to f (M 2 ) RVs = Constraints on instrumental offsets are reported at the bottom of Table 1.

Joint fit of RVs and astrometry
We next jointly fit the RVs and Gaia astrometry assuming a Keplerian 2-body orbit.The likelihood function is a product of two terms: one comparing the measured and predicted RVs, and another comparing the predicted astrometric orbital parameters to those measured by Gaia.The fit has 14 free parameters: orbital period, P orb , eccentricity, e, luminous star mass, M ⋆ , companion mass, M 2 , inclination, i, ascending node longitude, Ω, argument of periastron, ω, periastron time, T p , centerof-mass RV, γ, parallax, ϖ, right ascension, α, declination, δ, and proper motions, µ * α and µ δ .We adopt the best-fit instrumental offsets inferred in the RV-only fit rather than re-fitting them.
Our approach closely follows El-Badry et al. (2023b) and El-Badry et al. (2023a).For each call to the likelihood function, we predict both the RVs of the luminous star at the times of our RV measurements, and the vector of 12 Gaia-constrained parameters.The latter contains both some free parameters of the fit (P orb , e, T p , α, δ, ϖ, µ * α , and µ * δ ), and some quantities that are transformations of our free parameters (the Thiele-Innes coefficients, A, B, F , and G).The astrometric term of the likelihood compares the predicted and Gaia-constrained parameters while accounting for the full covariance matrix; the RV term assumes Gaussian RV uncertainties and that RVs are not covariant with other quantities.Some fitting parameters (e.g., P orb and e) are constrained by both RVs and astrometry.Others are constrained only by astrometry (e.g.α, δ, i) or only by RVs (e.g.γ).M ⋆ is constrained only by the prior.We assume broad, flat priors on all parameters except M ⋆ , for which we assume a Gaussian prior N (0.79, 0.01).We again use emcee to sample from the posterior, drawing 5000 samples with 64 walkers after a burn-in period of 5000 steps.The results are reported in the "Gaia+RVs" block of Table 3.

Astrometry alone
We also consider constraints from the Gaia astrometry alone.In this case, we repeat the procedure described in Section 5.2 but remove the RV term from the likeli- hood function.Because the binary's center-of-mass RV, γ, is unconstrained without RVs, we fix it to the value inferred from the joint fit.The resulting constraints are essentially a re-parameterization of the Gaia constraints into Campbell elements, with the only added information coming from the prior on M ⋆ .These constraints are listed in the "Gaia only" block of Table 3.

Predicted RVs
Figure 2 compares the observed RVs to predictions from the fits described above.The top panel shows predictions from the Gaia-only fit as cyan lines.Individual lines show 50 draws from the posterior, with the spread in these predictions representing the uncertainty in the Gaia constraints.Despite some uncertainty in the pre-dicted phase at the time of our observations, there is good agreement between the observed and predicted RVs in terms of overall shape, period, and phase.The total variability amplitude of the observed RVs is a few km s −1 lower than the median value predicted by the Gaia-only constraints.
The middle panel of Figure 1 compares the observed RVs to predictions from the joint Gaia+RV fit.These are much more tightly constrained than the Gaia-only predictions and are in good agreement with the data.The bottom panel shows the residuals between the observed RVs and best-fit joint model.These are consistent with 0 in most cases.The reduced χ 2 is 1.25, suggesting that the RV uncertainties may be underestimated by ≈ 10% on average.
-Comparison of constraints from the astrometric solution alone (cyan) and a joint fit of the astrometry and RVs (black).Contours enclose 12%, 39%, 68%, and 86% probability mass (see Foreman-Mackey 2016).The constraints from astrometry+RVs are tighter for all parameters except M⋆.The two sets of constraints are consistent for most parameters, with the largest tensions being for the parallax ϖ (3σ tension) and companion mass M 2 (2σ tension).The modest tension between the solutions leads us to inflate the astrometric uncertainties in fitting (Section 5.6); the resulting constraints can be found in Figure 13.
The relative roles of the RVs and astrometry in constraining the companion mass can be seen in the lower right panel of Figure 1, which shows the companion mass implied by the RV mass function and our constraints on M ⋆ as a function of the assumed inclination.The gray shaded region shows the uncertainty in M 2 associated with the uncertainty in both M ⋆ and the RV mass function.The Gaia inclination constraint, i = 68.8± 0.5 deg, implies a unseen companion mass of M 2 = 1.902 ± 0.015 M ⊙ .If we jettison the Gaia solution completely and consider only the constraint from the RV mass function, we obtain a minimum companion mass of M 2 = 1.67 ± 0.01 M ⊙ : still ≈ 25σ above the Chandrasekhar mass.

Consistency of the astrometric and RV constraints
Figure 3 compares the posterior constraints of key physical parameters obtained from joint fitting of the RVs and Gaia astrometry (black) and astrometry alone (cyan).Constraints on most parameters are consistent across the two solutions.The parameters not shown in the corner plot due to space constraints (α, δ, µ * α , µ δ , and γ) are all consistent at the 1σ level.The best-fit joint solution also predicts Thiele-Innes coefficients that are in agreement with the Gaia constraints to 1.1σ or better.
The most significant disagreement between the two sets of solutions is in the parallax, ϖ, where the disagreement reaches 3σ.This may seem surprising, given that our RVs do not directly constrain the parallax at all!The closer distance of the joint solution can be understood as a result of several compounding factors.First, Figure 3 shows that the astrometric solution includes significant negative covariances between ϖ and the parameters ω and T p , which are constrained by RVs.The RVs favor lower values of both parameters than the astrometric solution, drawing the joint fit toward larger ϖ.Second, the Thiele-Innes coefficients constrained by Gaia are functions of the photocenter semi-major axis, a 0 .The RVs constrain the physical semi-major axis of the luminous star's orbit, a ⋆ = a 0 /ϖ, to be smaller than in the best-fit astrometry-only solution.This leads to a smaller predicted a 0 at fixed distance.The larger parallax favored by the joint fit compensates for this, resulting in a larger a 0 , which is then compatible with the Gaia constraints on the Thiele-Innes coefficients.
If we do not leave the parallax free in the joint fit, we find a higher median companion mass, M 2 = 2.03 M ⊙ , and a lower inclination, i = 66.2 deg.These covariances can be understood as arising from the fact that the RVs constrain K ⋆ , which is related to a 0 , ϖ, and i through the equation Thus, a larger ϖ or smaller i can both result in a smaller predicted K ⋆ at fixed a 0 .Figure 2 shows that the observed RVs approximately fall on top of the predictions of the Gaia-only solution: the main tension is that the predicted RV semi-amplitude is ≈ 2 km s −1 larger than observed.This suggests that the Gaia constraints are basically correct; i.e., there was no catastrophic failure in the astrometry.However, it is still possible (and indeed, likely) that the uncertainties of the Gaia-only solution are underestimated somewhat: one parameter is discrepant between the two fits at the 3σ level, and our analyses of Gaia BH1 also suggested that the uncertainties of its astrometric solutions were underestimated at the factor of ∼2 level (El-Badry et al. 2023a;Chakrabarti et al. 2023;Nagarajan et al. 2023).

Joint fit with inflated astrometric uncertainties
To assess the effects of possible unrecognized systematics and/or underestimated uncertainties in the Gaia solution, we inflated the astrometric solution by a factor of 3 while maintaining the same parameter correlations reported in Gaia DR3.This is equivalent to multiplying the covariance matrix by 9. We then repeated the joint fit as described in Section 5.2.The results are reported in the last block of Table 3.This inflation makes the astrometry-only and astrometry+RV constraints fully consistent at the 1σ level.The companion mass constraint then changes to M 2 = 1.90 ± 0.04 M ⊙ .We adopt this as our fiducial and most robust value.
The choice of a factor of 3 uncertainty inflation is likely somewhat conservative since it makes all parameters consistent within 1σ (see Appendix A).

Astrometric phase coverage
Figure 4 explores the expected orbital motion and Gaia observations of Gaia NS1.In the left panel, we show the total plane-of-the-sky motion predicted by the joint RVs+Gaia solution, separating the effects of parallax and proper motion from those of orbital motion.Note that the axis scales in RA and Dec are very different, because the source's proper motion is almost exactly aligned with the declination axis. 3In the middle and right panels, we show the orbital motion only, and compare constraints from the joint fit to those from Gaia alone.As expected, the joint RV+astrometry solution results in a tighter constraint on the orbital ellipse than the astrometry alone.
The red symbols show times when the source is predicted to have been observed by the Gaia observation scheduling tool (GOST). 4The actual epoch-level astrometry is not published in DR3, so we simply plot the coordinates predicted by the best-fit joint solution at the timestamps when GOST predicts the source would have transited across the Gaia focal plane.Due to CCD chip gaps and various issues that can cause gaps in the datastream, only about 90% of predicted focal plane transits result in a usable astrometric measurement on average.
For Gaia NS1, GOST predicts 26 transits and 15 visibility periods (i.e., 15 groups of observations separated from other groups by at least 4 days), but the gaia source table reports that only 24 transits and 14 visibility periods were used in the astrometric solution.This implies that the scan times shown in Figure 4 are a serviceable but not perfect approximation of the actual scans used in the astrometric solution.The predicted scan times result in good coverage of the astrometric orbit.No matter which two predicted observations were not used for the astrometric solution, the ellipse would still be reasonably well sampled by any combination of 24 predicted transits in 14 visibility periods.This implies better astrometric phase coverage than was available for Gaia BH1, where all the predicted astrometric observations fell on one half of the orbit (El-Badry et al. 2023a).

SPECTROSCOPIC ANALYSIS
We analyzed the highest-SNR MIKE spectrum of Gaia NS1 to obtain a spectroscopic estimate of the luminous star's atmospheric parameters and abundance pattern.The spectrum was obtained on JD 2460037.6527 with a 2160 s exposure and has SNR = 102 at 6500 Å.Our analysis follows the hybrid spectroscopic/isochrone fitting method described in Reggiani et al. (2022).In brief, we first infer the effective temperature, T eff , surface gravity log g, metallicity, [Fe/H], and microturbulent velocity, ξ, through an iterative fit of the broadband SED and the equivalent widths (EWs) of Fe I and Fe II atomic absorption lines.We use an iterative approach that minimizes the dependence of the iron abundances inferred from individual lines on reduced equivalent width, while requiring that the atmospheric parameters, parallax, and SED are consistent with MIST isochrones.After the atmospheric parameters have been constrained, we infer individual elemental abundances while fixing the atmospheric parameters to the values measured in the first step.
For the isochrone fitting, we include multiwavelength photometry from the ultraviolet to the near-infrared: GALEX NUV, Gaia DR2 G, 2MASS J, H, and Ks, and WISE W1 and W2.We also include our updated parallax (ϖ = 1.36 ± 0.03; Table 3) 4 our LTE abundance inferences and uncertainties.Using spectral synthesis we also derived the abundances of lithium, cobalt, and europium.The lithium abundance was derived through spectral synthesis of the 6707 Å transition, for cobalt we used the 3894, 3995, 4020, 4110, and 4118 Å transitions, and for europium the abundance was derived via the 4129, 4205, and 6645 Å transitions.All linelists used for the synthesis analyzes are available upon request.For lithium we also performed corrections for 3D, and non-LTE effects using the breidablik code (Wang et al. 2021).
Our fit yields atmospheric parameters in good agreement with those inferred from the SED fit (Section 3.3).In most respects, the abundance pattern is unremarkable for an halo star.The α elements O, Ca, Mg, and Ti are enhanced ([α/Fe] ≈ +0.2) and there is no evidence for s-process enhancement ([Ba/Fe] = −0.17±0.13).The most unusual feature of the abundance pattern is a strong enhancement in lithium, A(Li) = 2.92±0.09after the 3D NLTE abundance correction.This value is almost 1 dex higher than typical for stars with temperatures similar to Gaia NS1.We discuss the lithium enhancement further  in Section 6.2.
6.1.Comparison of the spectrum to similar stars To search for possible anomalous spectral features, we compared the spectrum of Gaia NS1 to spectra of stars with similar stellar parameters and abundances observed by the GALAH survey (De Silva et al. 2015;Buder et al. 2021).To this end, we degraded the highest-SNR MIKE spectrum to R = 28, 000, shifted it to rest frame, and identified its nearest neighbor in pixel space among stars observed by GALAH with SNR > 50.The closest match was Gaia DR3 source 6251344007644295680 (GALAH DR3 object ID 170220004601155), which has T eff = 5980 ± 78 K, log g = 4.01 ± 0.18, [Fe/H] = −1.38 ± 0.08, and [α/Fe] = 0.32 ± 0.02.These parameters are all similar to those we find for Gaia NS1, providing independent validation of our inferred parameters.In Figure 5, we compare the normalized and resolution-matched spectra of the two sources in the 1st and 3rd GALAH wavelength windows, which respectively contain Hβ and Hα.The spectra are very similar.The most obvious differences are in the wings of the Hβ line -which could reflect either differences in continuum normalization or slight differences in log g -and in small telluric absorption lines near the Hα line.The two sources also have identical dereddened G BP − G RP colors within 0.01 mag, and identical absolute magnitudes within 0.03 mag when we assume our fiducial parallax of ϖ = 1.36 mas for Gaia NS1.Like Gaia NS1, the star is on a halo orbit, with a plane-of-thesky tangential velocity v ⊥ ≈ 370 km s −1 .The similarity of the two observed spectra and the sources' absolute magnitudes both speaks against the presence of a luminous secondary and validates the revised parallax from our joint Gaia+RV fit.

Lithium enhancement
Although most of the spectrum of Gaia NS1 is very similar to the metal-poor turnoff stars observed by GALAH, there is one respect in which the object is unique: it has an unusually strong Li I λ6708 line, indicative of a high lithium abundance.This is illustrated in Figure 6.The equivalent width of the line is 114 m Å, more than 3 times larger than the median value for ≈ 35 m Å for the top 20 closest spectral doppelgangers.Our inferred lithium abundance, A(Li) = 2.92 ± 0.08, is higher than that of any star observed by GALAH with similar parameters and a high-SNR spectrum (right panel of Figure 6).Within a temperature range of ±200 K and a metallicity range [Fe/H] < −0.7, there are 1097 stars observed by GALAH with SNR > 30, and none has lithium enhancement comparable to Gaia NS1.We inspected the spectra of the most lithium-rich GALAH sources with temperatures comparable to Gaia NS1 and found most to have spectral artifacts near the Li I λ6708 line, pointing to spurious A(Li) measurements.We did not find any GALAH sources with similar parameters to Gaia NS1 and comparably strong Li I λ6708 lines.The lithium enhancement of Gaia NS1 is thus very significant.
Lithium enhancement is commonly associated with young stars, since lithium is destroyed in stellar interiors.Gaia NS1 is obviously not young given its low metallicity and CMD position near the main-sequence turnoff.Strong lithium enhancement has also been found in the donors of several BH and NS X-ray binaries (e.g., Martin et al. 1992;Marsh et al. 1994;Martin et al. 1994;Filippenko et al. 1995;Martín et al. 1996).A clear reason for this enhancement has never been identified, with proposed explanations including lithium production in supernovae (e.g., Dearborn et al. 1989), spallation on the surface of the companions due to the neutron flux from the compact object's accretion flow (e.g., Guessoum & Kazanas 1999;Fujimoto et al. 2008), spallation within the accretion flow and subsequent pollution of the companion by winds (Yi & Narayan 1997), and slowing of lithium destruction in the companions due to their rapid rotation (Maccarone et al. 2005).In Gaia NS1, the wide separation of the star and companion, negligible expected accretion rate of the compact object, and slow rotation of the companion make all these proposed explanations seem unlikely.
Lithium is also produced in AGB and super-AGB stars (e.g., Cameron & Fowler 1971;Sackmann & Boothroyd 1992;Ventura & D'Antona 2010), raising the question of whether the companion could be a WD.Our measured companion mass is far too high for it to be a single WD, but WD+WD binary models are in principle possible (Section 8.2).Lithium enhancement is not, however, commonly observed in stars with WD companions and other abundance anomalies associated with pollution from AGB companions (Pinsonneault et al. 1984;Allen & Barbuy 2006).Moreover, most lithium-enhanced stars do not exhibit significant RV variability (Castro-Tapia et al. 2023;Sayeed et al. 2023).We have not found evidence of lithium enhancement in our follow-up of astrometric binaries from Gaia with inferred masses significantly below the Chandrasekhar mass.
The other two metal-poor NS candidates we have identified (El-Badry et al. 2024, in prep)  -Spectral cutouts of Gaia NS1 (red) compared to a reference star (black), a halo star with similar stellar parameters and abundances observed by the GALAH survey.The two spectra are very similar, ruling out significant light contributions from a companion and validating the stellar parameters and abundances inferred in Section 6.
eccentric orbits (e = 0.66 and e = 0.68).These masses are typical of normal NSs, and the higher eccentricities are consistent with what is expected from natal kicks.We have found little evidence for lithium enhancement in any of the candidates with metallicities near solar.
Assuming the enhancement is due to some form of pollution from a companion, it is expected that it should be stronger at low metallicity: low-metallicity stars have thinner convective envelopes, meaning that the enriched material is mixed with less material from the luminous star.This is a fairly significant effect: from MESA models of 0.8 M ⊙ stars, we find that in a model with initial metallicity Z = 0.0014, the convective envelope contains only ≈ 0.2% of the star's mass at the observed evolutionary phase of Gaia NS1, compared to ≈ 10% in a model with Z = 0.014.That is, accreted material is diluted ∼ 50 times less on the surface of the luminous star in Gaia NS1 than it would be in a solar-metallicity companion.
7. GALACTIC ORBIT To investigate the Galactic orbit of Gaia NS1, we used the parallax, proper motion, and center-of-mass RV from the joint Gaia+RV fit as starting points to compute its orbit backward in time for 1 Gyr using galpy (Bovy 2015).We used the Milky Way potential from McMillan (2017).The system's current LSR-corrected Galactic space velocity is UVW = (218, −250, −88) km s −1 .
The orbit is shown in Figure 7; for comparison, we also show the orbit of the Sun.The orbit is typical of a halo star: it is highly eccentric and reaches several kpc above the disk midplane.This conclusion is insensitive to variations of the source's astrometric parameters within their uncertainties and/or the adopted Galactic potential.
We note that the wide orbit and low eccentricity of Gaia NS1 rule out a large natal kick to the NS.The system's large space velocity must thus be a consequence of its being very old.
8. DISCUSSION 8.1.Could the luminous star have an unusually low mass?Several dormant BH candidates published in recent years turned out to be mass transfer binaries in which the luminous star had a much lower mass than assumed when its temperature and radius were compared to single-star models (e.g., Shenar et al. 2020;Bodensteiner et al. 2020; El-Badry & Quataert 2021; El-Badry & Burdge 2022).It is worth considering whether the luminous star in Gaia NS1 could similarly be a product of mass transfer, in which case the companion's mass could be overestimated as well.
To test whether a pre-WD could conceivable masquerade as a main-sequence star similar to the luminous star in Gaia NS1, we calculated a grid of MESA binary evolution models (Paxton et al. 2011(Paxton et al. , 2015)).We considered donor stars with initial masses of (1 − 2) M ⊙ and initial periods of 50-150 d, modeling the companion as a (1 − 2) M ⊙ point mass.Several of these models detach from their Roche lobes at periods near 730 d, contract, and become helium WDs.However, their temperatures and radii never simultaneously approach those of mainsequence stars similar to Gaia NS1: when they have R ∼ 1 R ⊙ , these models have T eff = 30, 000 − 40, 000 K. By the time they cool to T eff ∼ 6000 K, they are on the  WD cooling track, with R ∼ 0.01 R ⊙ .We thus conclude that stripped stars are unlikely to masquerade as solarmass main sequence stars.
Moreover, for a pre-WD formed by stable mass transfer, the expected mass at a period of 730 d is ≈ 0.45 M ⊙ (e.g., Rappaport et al. 1995).If the mass of the luminous star were somehow this low, the dynamically-required mass of the companion would still be M 2 ≈ 1.55 M ⊙ .This is significantly above the Chandrasekhar limit and would thus most simply imply a NS companion.
At short periods, "semi-stripped" donors with thick hydrogen envelopes can have radii and temperatures similar to solar-mass main sequence stars (e.g., El-Badry et al. 2021a,b), but the evolutionary scenario that gives rise to such donors requires stable mass transfer to a close companion (P orb < 1 d), which is clearly not applicable here.

Nature of the companion
All our evidence of the companion comes from its gravitational effects on the luminous star, whose motion we can probe with both astrometry and RV measurements.While we can constrain the total mass of the companion dynamically, this provides only indirect information about its nature.We discuss several possibilities below.
1.A single neutron star: This is the most straightforward explanation and our default assumption in most of the text.was ≈ (2.5 − 2.7) M ⊙ .Given our lower mass constraint of ≈ 1.90 M ⊙ , a BH seems unlikely: there are no well-understood astrophysical channels to produce BHs of this mass, and several objects that are definitely NSs have fairly well-constrained masses above 1.9 M ⊙ (e.g., Antoniadis et al. 2013;Cromartie et al. 2020).
3. One, two, or three main-sequence stars, or other kinds of dim but luminous companions: A luminous companion that is a single star, a tight mainsequence binary, or even a tight main-sequence triple, is ruled out.This is because any companion that is not dark would require the luminous star's orbit, a ⋆ , to be larger than the photocenter orbit, a 0 , and thus would require an even larger companion mass.This is illustrated in Figure 8, which compares the mass-to-light ratios of various types of companions (dashed and dotted colored lines) to the dynamical constraint from the astrometric mass function (black line).For example, if a hypothetical dim companion contributed 15% of the G-band light, then the companion mass implied by the astrometric mass function would rise to 3.0 M ⊙ .Figure 8 shows that a WD+MS inner binary is ruled out for the same reason.Joint fitting of astrometry and RVs directly constrains the fraction of the light coming from the companion, because a luminous companion implies larger RV variability amplitude at fixed a 0 .To quantify this limit, we tried fitting ϵ, the G-band flux ratio of the companion to the primary, as a free parameter (see El-Badry et al. 2023b).We conservatively inflated the astrometric uncertainties by a factor of 3 (see Section 5.6) and found 1 and 2σ upper limits of ϵ < 0.02 and ϵ < 0.03, respectively.These limits are quite stringent and rule out almost all scenarios involving non-degenerate companions.
4. A neutron star + white dwarf binary: The unseen object could also be a close binary containing a neutron star and a white dwarf, with masses of, e.g., 1.4 M ⊙ + 0.5 M ⊙ or 1.7 M ⊙ + 0.2 M ⊙ .Such an inner binary could be effectively invisible in the optical.How such a system could have formed is quite uncertain.In the 1.7 M ⊙ + 0.2 M ⊙ scenario, the inner binary could be a result of stable mass transfer and the NS should be a recycled millisecond pulsar.If the inner binary was very tight, it may subsequently have evolved through an ultra-compact X-ray binary phase, possibly leaving behind an isolated NS (Tauris & van den Heuvel 2023).In the 1.4 M ⊙ + 0.5 M ⊙ scenario, it would have to be a product of common envelope evolution to fit inside the luminous star's orbit today (e.g., Rappaport et al. 1995).Either scenario requires a common envelope episode prior to the formation of the NS.This seems difficult to achieve, because a lowmass companion is only expected to survive such an interaction with a red supergiant progenitor if it begins in an orbit of several au, beyond the current orbit of the luminous star (e.g., Dewi & Tauris 2000).
There is one known NS in a triple: the millisecond pulsar PSR J0337+1715 (Ransom et al. 2014;Kaplan et al. 2014).That system contains a 1.43 M ⊙ NS, in a 1.6 d orbit with a 0.20 M ⊙ WD, all orbited by a 0.41 M ⊙ WD in a 327 d orbit.How that system formed is not well understood: two intriguing but yet uncertain scenarios involving a double common envelope have been proposed (Tauris & van den Heuvel 2014;Sabach & Soker 2015).In any case, it is plausible that Gaia NS1 could have formed through a similar pathway to that system.This could involve, for example, the luminous star being born in a wider orbit than it inhabits today and spiraling inward following interaction with the NS progenitor's unbound envelope after its ejection through a common envelope interaction with an inner star (e.g., Tauris & van den Heuvel 2014).

5.
A white dwarf + white dwarf binary: Our dynamically-inferred mass is also consistent with two WDs.Both would have to be rather massive, e.g., 1.0 M ⊙ + 0.9 M ⊙ or 1.3 M ⊙ + 0.6 M ⊙ .
No close WD+WD binaries are known with total mass exceeding ∼ 1.4 M ⊙ (e.g., Napiwotzki et al. 2020), so such a binary would be the most massive WD+WD binary known by a large margin.Given that Gaia NS1 is ≳ 12 Gyr old and the WD progenitors would have had MS lifetimes of at most a few Gyr, both WDs would be ancient, colder than the luminous star, and essentially impossible to detect.As with a NS+WD binary, the main weakness of the WD+WD binary scenario is the difficulty of forming such a massive inner binary while retaining the luminous star as a bound tertiary in a stable orbit.Population synthesis simulations of triples with inner WD+WD binaries predict the tightest triples to have outer semimajor axes a outer ≈ 50 au, much wider than Gaia NS1 (Shariat et al. 2023).
We discuss some possible evolutionary scenarios below.
The simplest WD+WD binary model is one in which the binary was born tight and the luminous star never interacted with the progenitor of either WD.Such a scenario, is, however, difficult to make work in practice.Most formation channels for massive WD+WD binaries require the binary to expand to dimensions of order 1 au (e.g., Ruiter et al. 2013), since tighter orbits lead to early envelope stripping and formation of lower-mass WDs.The total initial mass of the inner binary would have to have been significantly greater than the mass of the dark companion today, since (0.9 − 1.0) M ⊙ WDs are formed from stars with initial masses of (4−6) M ⊙ (e.g., El-Badry et al. 2018;Cunningham et al. 2024).For a conservative total initial inner binary mass of 8 M ⊙ , the initial orbital period of the luminous star would have been ≈ 69 d (Jeans 1924).
In order to remain stable, the inner binary period then could not have exceeded ≈ 69/5 d during its evolution (e.g., Mardling & Aarseth 2001).We experimented with a variety of MESA binary models in search of an evolutionary scenario in which a binary remains tight and produces two WDs with total mass ∼ 1.9 M ⊙ but were unable to find one: mass transfer widens the orbit too much, unless the initial mass ratio is highly unequal.Binaries with highly unequal initial mass ratios can remain tight, but in these cases, mass transfer becomes dynamically unstable and likely leads to a stellar merger.
A more speculative scenario that might lead to a WD+WD inner binary places the initial orbit of the luminous star wider than its current orbit, e.g. with an initial period of ∼ 10 yr.In this case, the (hypothetical) third star could have been born with an initial period of order 1 yr, gone through a common envelope when the primary became an AGB star, and gone through another phase of interaction when it terminated its main-sequence evolution.The tertiary could, perhaps, spiral inward due to interaction with the AGB star's wind or common envelope ejected preferentially in the orbital plane (e.g.Bodenheimer & Taam 1984;Passy et al. 2012).An attractive feature of such a scenario -which is similar to the one explored by Tauris & van den Heuvel (2014) for PSR J0337+1715 -is that the tertiary could accrete Li-enriched material from the AGB star's wind.However, it is quite uncertain whether this scenario can work in practice, because the presence of the inner companion will likely suppress mass transfer to the outer tertiary, preventing its inspiral.
We conclude that a single massive NS is the simplest explanation for the data.A single WD is untenable.Scenarios involving an inner binary (WD+WD or WD+NS) cannot be ruled out based on our data alone, but are challenging to explain with evolutionary models.Such scenarios could be tested with high-precision RV observations (e.g., Nagarajan et al. 2023), which would be sensitive to inner binaries with periods longer than about a week.The NS+WD scenario can potentially be tested with radio observations, since the NS would likely be a recycled pulsar.However, the beaming fraction even of a millisecond pulsar is only ∼ 0.3 (Lorimer & Kramer 2012), so it may not beam in our direction in any case.In the following discussion, we focus on the scenario in which the companion is a single NS.

Difficulty of explaining the current wide orbit
The current separation between the luminous star and dark companion is ≈ 2.2 au.This separation is uncomfortable because it is both too tight for the progenitor of the NS to have fit inside as a red supergiant or AGB star, and too wide to be a result of common envelope evolution.Essentially the same problem applies to the orbits of Gaia BH1 and BH2, to the 20 other NS candidates we have followed-up (El-Badry et al. 2024, in prep.), and perhaps also to most of the ∼3,000 WD+MS binaries with astrometric orbits in Gaia DR3 (Yamaguchi et al. 2023;Shahaf et al. 2023a).It is possible that these systems can form through common envelope evolution if the compact object progenitor's envelope is very weakly bound and/or if significant energy is released through recombination (Webbink 2008;Izzard & Jermyn 2018;Yamaguchi et al. 2023;Belloni et al. 2024).Alternatively, mass transfer may remain stable even under extreme mass ratios (e.g., Nelemans et al. 2000).The fact that a rather large population of such wide post-interaction binaries exists suggests that their formation does not require exotic conditions or fine-tuned scenarios.The fact that most of the wide WD binaries have near-circular orbits (as does Gaia NS1) suggests that these systems are not primarily assembled dynamically, as has been proposed for the BH systems (Rastello et al. 2023;Di Carlo et al. 2023;Tanikawa et al. 2024).

Effects of low metallicity
We next consider whether the low metallicity of the Gaia NS1 system could have resulted in a compact NS progenitor that never overflowed its Roche lobe.At the current orbital separation of a ∼ 2.2 au, a (10 − 20) M ⊙ red supergiant would overflow its Roche lobe when its radius reached (280−310) R ⊙ .Observed red supergiants reach maximum radii of ∼ 1000 R ⊙ (Yang et al. 2020) at SMC metallicity, which is about twice the metallicity of Gaia NS1.Observations of massive stars with Z ≲ 0.1 Z ⊙ are scarce (e.g., Gull et al. 2022), making their late-stage evolution uncertain.
A realistic constraint on the maximum radius of the progenitor if there was no Roche lobe overflow is more restrictive than the current orbit, because the progenitor would have been more massive than the NS and the orbit would have expanded in response to mass loss.An initial mass of (10 − 20) M ⊙ implies an initial separation of (0.3 − 0.6) au, which would lead to Roche lobe overflow of the progenitor if its radius exceeded (40 − 75) R ⊙ .-Constraints on kicks and the pre-SN orbit of Gaia NS1.We assume that the binary had a circular orbit before the SN and that the NS received a kick with velocity v kick when it formed.We simulate 10 8 pre-SN orbits with a uniform period distribution P orb /d ∼ U (0, 1000), a pre-SN mass distribution M He star /M ⊙ ∼ U (1.9, 5.0), and a kick velocity distribution v kick /(km s −1 ) ∼ U (0, 100).We then show properties of the systems whose post-SN orbits have periods and eccentricities similar to Gaia NS1 (upper left).Most simulated systems that can produce the observed orbit have weak kicks (v kick ≲ 20 km s −1 ) and formed from low-mass He stars (M He star ≲ 3 M ⊙ ), though there is a long tail of higher kick velocities and more massive He stars that could produce the observed orbit for fine-tuned kick orientations.
Models have been proposed for low-metallicity massive stars that never become red supergiants and indeed reach maximum radii of ≲ 40 R ⊙ .These include, for example, some early models for the progenitor of SN 1987A, which exploded at R ∼ 40 R ⊙ (e.g., Weiss 1989).These models do not, however, match the observed distribution of evolved stars in the HR diagram at LMC metallicity: more red supergiants are observed there than at solar metallicity, while the models predict the opposite (Pod-siadlowski 1992).This is likely a result of weaker winds at low metallicity, which were not accounted for in early calculations.
More recent models predict that (10 − 20) M ⊙ stars at low metallicity do expand to become red supergiants, but only do so during their late-stage evolution, during core He burning rather than in the Hertzsprung gap.At Z = 0.01 Z ⊙ (lower than Gaia NS1), a 14 M ⊙ model calculated by Klencki et al. (2020)  -Predicted eccentricity of a 1.9 M ⊙ NS + 0.8 M ⊙ MS binary following a supernova.We assume the orbit is initially circular, that the NS forms from a helium star of mass M He star , and that it receives a natal kick of velocity v kick in a random direction at formation.Solid lines and shaded regions show median and middle 68% eccentricity regions of the surviving orbits.In the absence of a natal kick, loss of more than 0.3 M ⊙ leads to a higher eccentricity than observed.Kicks generally increase the eccentricity.The low observed eccentricity of Gaia NS1 thus suggests that the NS formed from a low-mass He star with a small natal kick.
its He-burning lifetime with R < 50 R ⊙ , only expanding to R > 100 R ⊙ in the final 10 5 yr before supernova, when the envelope is tenuous and might be easier to unbind.However, similar models at Z = 0.1 Z ⊙ expand to R = 300 R ⊙ soon after beginning He burning.We conclude that Gaia NS1's low metallicity is unlikely to have prevented the NS progenitor from overflowing its Roche lobe.-Probability of remaining in a bound orbit with eccentricity e < 0.13, the value observed for Gaia NS1.As in Figure 10, we assume that the orbit is initially circular, that the NS forms from a helium star of mass M He star , and that it receives a natal kick of velocity v kick in a random direction at formation.For very low v kick values, a low He star mass is required, because mass loss will otherwise result in an eccentric orbit.For higher v kick , a bound and near-circular orbit can be maintained if the kick is fortuitously aligned with the companion's instantaneous velocity vector, but the probability of this occurring is low.

Expected helium star progenitor
Models predict that mass transfer from red supergiants proceeds until all but a few tenths of a solar mass of hydrogen have been stripped, and that most of the residual hydrogen will be removed by winds (e.g., Yoon et al. 2017;Gilkis et al. 2019).Given the presence of the luminous star, it is likely that the progenitor of the NS was stripped by binary mass transfer and terminated its evolution as a hydrogen-poor supernova.
The mapping between He star mass and compact remnant mass is expected to be complex and non-monotonic, and is complicated further by mass transfer in a binary (e.g., Laplace et al. 2021).In one recent set of models, Vartanyan et al. (2021) predict that NSs with masses near 1.9 M ⊙ are produced from stripped He stars with masses of (5 − 7) M ⊙ .These form from stars with initial masses of (17 − 20) M ⊙ in their solar-metallicity models, and likely from somewhat lower initial masses at lower metallicity.Lower-mass He stars are predicted to produce lower-mass NSs, with M < 1.6 M ⊙ .Other models in the literature predict minimum He star masses for massive NSs (M ≳ 1.8 M ⊙ ) of ∼ 3 to 7 M ⊙ (Mandel & Müller 2020;Burrows et al. 2020;Antoniadis et al. 2022).

Constraints on mass loss and kicks
We use Monte Carlo simulations to study the combinations of pre-supernova orbits and kick velocities and orientations that could have produced Gaia NS1.This approach follows Tauris et al. (2017) and Larsen et al. (2024); for a general review on the effects of SNe in binaries, see Tauris & van den Heuvel (2023).In brief, we simulate a large number (N = 10 8 ) of orbits and kicks and use the formalism from Brandt & Podsiadlowski (1995) to predict the post-SN parameters.We assume that prior to the formation of the NS, the orbit was circular and the NS progenitor was a He star of mass M He .This object's mass instantaneously drops to M NS = 1.90 M ⊙ when the NS forms.In addition, the NS receives a kick of velocity v kick in a random direction.We fix the luminous star mass to 0.79 M ⊙ and the NS mass to 1.9 M ⊙ .Finally, we consider the initial parameters of the simulations for which the final period and eccentricity are close to the observed values for Gaia NS1: P orb,final = 700 − 760 d, and e final = 0.11 − 0.13.Narrowing these ranges reduces the number of surviving Monte Carlo samples but has no significant effect on their distribution.
We begin with a uniform distribution of pre-SN orbital periods between 0 and 1000 d.We simulate kicks assuming velocities distributed uniformly between 0 and 100 km s −1 , and take He star masses distributed uniformly between 1.9 and 5.0 M ⊙ .This is not a population synthesis simulation -pre-SN orbits are unlikely to be uniformly distributed in all parameters -but should be viewed as an exploration of what combinations of orbits, kicks, and mass loss could have produced the observed orbit.
The results of this experiment are shown in Figure 9. Pre-SN periods ranging from ∼ 400 to ∼ 900 d can produce the observed orbit for suitable combinations of mass loss and kicks.Most of the accepted kicks have ϕ kick ≈ 180 deg, meaning that the kick points in the opposite direction of the He star's motion at the time of the SN.The median v kick for SNe that match the observed orbit is 7.9 km s −1 , and 90% have v kick < 32 km s −1 .The median He star mass is 2.5 M ⊙ , and 90% of accepted simulations have M He star < 4.1 M ⊙ .
Although the parameters of Gaia NS1 can be produced with a significant kick (up to ∼ 80 km s −1 ) for fine-tuned orientations and He star masses, most NS + MS star binaries are expected to end up in more eccentric orbits in this case.This is illustrated in Figure 10, which shows the predicted median and middle 68% range of the eccentricities as a function of He star mass and kick speed for a pre-SN period of 730 d.Only binaries that remain bound are considered.In the case of no kick (top panel), the eccentricity depends only on the mass lost during the NS's formation: e = (M He − M NS ) / (M NS + M ⋆ ) (Blaauw 1961;Hills 1983).In this case, the He star progenitor could not have had a mass larger than 2.2 M ⊙ , corresponding to a total mass loss of 0.3 M ⊙ during the SNe.Given that neutrinos alone are expected to remove 0.2 − 0.3 M ⊙ during the formation of a 1.90 M ⊙ NS (Lattimer & Prakash 2001), this leaves little room for any other mass loss in the absence of a kick.
The lower panels of Figure 10 show predictions of simulations in which the kick velocity is increased to 5, 10, and 20 km s −1 .This generally increases the predicted eccentricity, such that in simulations with v kick > 5 km s −1 , most binaries end up with a higher eccentricity than is observed for Gaia NS1, regardless of the progenitor He star's mass.
Although larger kicks tend to produce higher eccentricities and unbound orbits, there is always a small probability that a fortuitously aimed kick with velocity comparable to the luminous companion's pre-supernova orbital velocity can result in a bound, low-eccentricity orbit.This essentially requires the kick to be oriented such that the NS chases the fleeing companion, and so such kicks result in a net systemic velocity to the binary.
Figure 11 explores the probability of such a fine-tuned kick for a variety of He star masses and kick velocities.For each kick velocity, there is a He star mass that maximizes the probability of a bound, low-eccentricity or-bit; this is approximately the mass for which the Blaauw (1961) kick velocity is equal to the NS's natal kick.The probabilities in question are all rather low except when the kick is small and the He star has low mass.For v kick = 10 km s −1 , the maximum probability of producing a system like Gaia NS1 occurs for M He star ≈ 3.5 M ⊙ and is about 8%.For v kick = 15 km s −1 , the maximum probability falls to ≈ 4% at M He star ≈ 5 M ⊙ .
Light curve modeling of stripped-envelope supernovae suggests typical ejecta masses of (1−3) M ⊙ (Lyman et al. 2016;Taddia et al. 2018), which for Gaia NS1 would correspond to He star masses of ∼ (3 − 5) M ⊙ .This is within the range of plausible progenitors we infer for kick velocities of v kick = 5−20 km s −1 , though only a minority of kick orientations would result in an orbit as circular as Gaia NS1 under these conditions.It remains to be seen whether other binaries like Gaia NS1 exist with higher eccentricities.

Comparison to other NSs
Figure 12 compares our constraints on the mass and orbital period of Gaia NS1 to known NSs in binaries with well-measured masses, the parameters of which we take from Özel & Freire (2016) and Fonseca et al. (2021).Cyan points show non-recycled neutron stars; i.e., systems with NS spin periods of order 1 s and magnetic fields of order 10 12 G, which have (presumably) not been spun up by accretion.Most of these systems are NS+NS binaries with eccentric orbits and high-precision mass constraints from post-Keplerian effects.Red points show recycled pulsars, most of which have spin periods of a few milliseconds and low-mass WD companions.NSs in spider binaries (e.g Romani et al. 2022) are not included.
The recycled NSs generally have larger mass uncertainties because the orbits are circular (leading to fewer measurable post-Keplerian parameters) and the constraints in some cases depend on spectroscopic estimates of the mass of the WD companions.Most of the highest-mass NSs are recycled, meaning that they have accreted significant material and their present-day mass may be significantly higher than their birth mass.Most of the nonrecycled NSs have masses near 1.3 M ⊙ , with only two having masses above the Chandrasekhar mass.
Gaia NS1 is the fourth most massive object in Figure 12, with only the recycled pulsars J0740+6620, J0348+0432, and J1614-2230 having higher masses.If the dark object is indeed a single NS, there is no plausible scenario in which it gained mass since its formation, so it would be one of the strongest known cases for a NS being born massive.That being said, J0348+0432 is only mildly recycled, implying that it may have been born with a mass near 2.0 M ⊙ (Antoniadis et al. 2013), andTauris et al. (2011) have argued on evolutionary grounds that J1614-2230 was born with a mass of at least 1.6 M ⊙ , and probably higher.There are also some reports of young NSs with relatively high masses in high-mass Xray binaries (e.g., Barziv et al. 2001;Falanga et al. 2015); these systems are not included in Figure 12 because their masses generally carry higher uncertainties.
Gaia NS1 is also the longest-period NS in a binary with a precise mass estimate.A few pulsars in longerperiod orbits are known (e.g., Wang et al. 2004;van der Wateren et al. 2023), but without precise mass mea-surements.When the luminous star ascends the giant branch, it will begin transferring mass to the companion, first by winds, and then by stable Roche lobe overflow.This likely makes Gaia NS1 the first known progenitor of a symbiotic X-ray binary (a NS accreting from a red giant; e.g.Yungelson et al. 2019).

CONCLUSION
We have presented discovery of Gaia NS1, a relatively nearby (d ≈ 735 pc) binary containing a low-metallicity turn-off star (M ⋆ ≈ 0.79 M ⊙ ) and an unseen massive companion that we suspect is a neutron star (NS) in a 2yr orbit.Joint modeling of the system's Gaia astrometry and our follow-up RVs spanning most of an orbital period allows us to constrain the system's parameters with fewpercent precision (Tables 3).Our main conclusions are as follows: 1. Properties of the luminous companion: The luminous companion is a halo star near the main sequence turnoff (Figure 1).Its low metallicity and α−enhanced abundance pattern ([Fe/H] = −1.23;[α/Fe] = +0.20),high space velocity (v pec ≈ 300 km s −1 ) and ancient age (τ ≳ 12 Gyr) all suggest that it is a halo star, probably formed in a since-accreted satellite of the Milky Way.The star's abundance pattern (Table 4) is in most respects unremarkable.However, the star is highly enriched in lithium, with a surface abundance higher than any of the 1000+ stars with similar parameters observed by the GALAH survey (Figure 6).This lithium enhancement is almost certainly related to the companion, but exactly how is unclear.Unexplained lithium enhancement has been found in the companions to several BH and NS X-ray binaries.Lithium is also produced in massive AGB stars, suggesting possible pollution from the progenitor of a WD or a NS formed through an electron-capture supernova.However, these scenarios are expected to produce a compact object below the Chandrasekhar mass and are likely only viable for Gaia NS1 if the companion is a tight binary.The star's parameters and evolutionary state imply a mass of 0.79 ± 0.01 M ⊙ when interpreted with single-star evolutionary models.We expect the mass inferred in this way to be reliable because any mass transfer would have occurred many Gyr ago in all the evolutionary scenarios we consider.The star's broadband spectral energy distribution (Figure 1) and high-resolution spectra (Figure 5) are well-matched by single-star models.Even if we consider the (remote) possibility that the distance could be significantly different from our fiducial assumptions, the star's temperature and metallicity imply a minimum mass of 0.75 M ⊙ .We consider and reject evolutionary scenarios in which the star is a low-mass stripped object, as these do not produce stars with temperatures and surface gravities near the observed values (Section 8.1).
2. Orbit and companion mass: Our RV follow-up covers the full dynamic range of the orbit and thus tightly constrains its parameters.Joint fitting of the astrometry and RVs constrains the companion mass to M 2 = 1.90 ± 0.01 M ⊙ when the Gaia astrometric uncertainties are taken at face value (Table 3).The RVs evolve in a way that is basically consistent with the predictions of the pureastrometry solution (Figure 2) and there are no significant deviations from a Keplerian orbit.However, the observed RV variability amplitude is ∼ 3 km s −1 smaller than the median value predicted by the astrometric solution, implying a smaller physical scale for the orbit.For this reason, the best-fit companion mass inferred from our joint RV+astrometry fit is ∼ 2σ lower than the value inferred from astrometry alone.
There is moderate tension between the parameter estimates from the astrometric orbit alone and our joint astrometry + RV fit (Figure 3).The most discrepant parameter between the two sets of constraints is the parallax, for which the Gaia-only constraint is ϖ = 1.24±0.04mas but the constraint from joint fitting of the astrometry and RVs is ϖ = 1.36±0.01mas.We ascribe this tension to underestimated uncertainties in the astrometric solution, as was also found for Gaia BH1 (Chakrabarti et al. 2023;Nagarajan et al. 2023).We therefore inflated all the astrometric uncertainties by a factor of 3 and repeated the joint RV + astrometry fit.
In this case, constraints between the astrometryonly and astrometry+RVs fit are consistent at the 1σ level, and the companion mass is constrained to M 2 = 1.90 ± 0.04 M ⊙ .The orbit is nearly circular, with an eccentricity e = 0.122 ± 0.003.The inferred companion mass is only weakly dependent on the Gaia astrometry, because the orbit is relatively edge-on (i = 69 deg; Figure 1).If the orbit were edge-on (i = 90 deg), the RVs and luminous star mass would imply a companion mass of M 2 = 1.67 M ⊙ .
3. Nature of the companion: We have not detected any electromagnetic radiation from the companion, so our speculation about its nature is informed only by our constraints on its mass and considerations of possible formation histories.Its mass is clearly too high to be a white dwarf, and likely too low for it to be an astrophysical BH.Luminous companions are ruled out because they would imply even higher companion masses (Figure 8) and RV amplitudes that are inconsistent with our observations; indeed, joint fitting of the astrometry and RVs sets an upper limit on the flux ratio of ϵ < 0.02, even with the astrometric uncertainties inflated.A single, relatively massive NS thus appears to be the simplest possible explanation.
The only plausible alternative is that the dark companion is a tight binary containing two massive WDs, or a WD and a NS.This scenario is consistent with the data but difficult to explain with evolutionary models (see Section 8.2 for discussion).Models in which the inner binary remained tight over its whole evolution fail to form sufficiently massive WDs, and to explain the observed  -Masses and periods of NSs with well-characterized orbits and mass measurements compared to Gaia NS1.Cyan points show radio pulsars that are non-recycled neutron stars (mainly in double NS binaries and a few NS + massive WD binaries), while red points show recycled NSs (millisecond pulsars with low-mass WD companions).Compared to these systems, Gaia NS1 hosts one of the most massive NSs in the widest orbit.Its parameters are particularly unusual given that the NS presumably has not accreted from the luminous star.This suggests that the NS was born massive (it is not recycled, and is most similar to the cyan points).lithium pollution of the luminous star.Models in which both inner and outer companions to a massive star went through simultaneous common envelope events provide a potential solution but are very uncertain.
4. Formation history: The system's formation history is also difficult to explain if the dark object is a single high-mass NS.The low eccentricity makes dynamical formation channels (i.e., exchange interactions in a dense cluster) unlikely.Assuming the system formed from a primordial binary, the progenitor of the NS likely had an initial mass of 10 − 20 M ⊙ .A star of this mass is expected to become a red supergiant that would have overflowed its Roche lobe in the current orbit.Because the donor-to-accretor mass ratio would be highly unequal in this scenario, mass transfer is expected to be unstable.However, unstable mass transfer (i.e., common envelope evolution) is predicted to often result in a much tighter orbit than observed.These considerations suggest that (a) mass transfer was stable despite a highly equal mass ratio, (b) common envelope ejection was very efficient, or (c) the system did not form via isolated binary evolution.This is essentially the same "problem" that has been pointed out for isolated binary evolution channels forming Gaia BH1 and BH2.

Constraints on natal kicks:
The wide and nearcircular orbit of Gaia NS1 today places strong constraints on mass loss and kicks during the NS's formation (Figures 9,10 and 11).The simplest explanation for the observed orbit is a weak natal kick (v kick ≲ 20 km s −1 ) and little ejected mass (M He star ≲ 3 M ⊙ , corresponding to an ejecta mass of ≲ 1 M ⊙ .)Stronger kicks are possible if coupled with more mass loss, but only for fine-tuned kick orientations.
Like Gaia BH1 and BH2, Gaia NS1 is an unexpected discovery.The system's wide orbit and low-mass luminous star make it difficult to form through the standard evolutionary channels invoked for the formation of lowand high-mass X-ray binaries and millisecond pulsars.On the other hand, a few low-mass red giant + NS binaries are known in symbiotic binaries with orbital periods comparable to Gaia NS1 (e.g., Hinkle et al. 2006).Such systems must have formed from wide MS + NS binaries like Gaia NS1 and are X-ray bright only for the short period in which the companion is a giant.There is thus little doubt that wide low-mass MS + NS binaries exist, and we suspect that Gaia NS1 is the first clear example of one to be identified.
search Council for the ERC Advanced Grant [101054731].This research benefited from discussions in the ZTF Theory Network, funded in part by the Gordon and Betty Moore Foundation through Grant GBMF5076, and from collaboration at the "Renaissance of Stellar Black-Hole Detections in The Local Group" workshop hosted at the Lorentz Center in June, 2023.
This research made use of pystrometry, an open source Python package for astrometry timeseries analysis (Sahlmann 2019).This work made use of Astropy, 7 a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration et al. 2022).
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.

Fig. 1 .
Fig. 1.-Observational properties and mass constraints of Gaia NS1.Upper left: extinction-corrected color-magnitude diagram.Gaia NS1 falls on the main sequence and is slightly bluer and brighter than the Sun.Upper right: broadband spectral-energy distribution.Red points show observed photometry; open black squares show the integrated fluxes predicted for the best-fit model spectrum (black line).Lower left: comparison of the measured temperature and radius to MIST single-star evolutionary models.Solid point shows the radius inferred when adopting the parallax from our joint Gaia + RV fit; hollow point shows the parallax from Gaia alone (see Section 3.1).The observed temperature, radius, and metallicity imply a luminous star mass of 0.79 ± 0.01 M ⊙ .Lower right: dynamical constraints on the companion mass.The RVs and luminous star mass constrain the companion mass and inclination to the gray shaded region, with a minimum companion mass of 1.67 M ⊙ for an edge-on orbit.The inclination constraint from Gaia astrometry then implies M 2 = 1.90 M ⊙ .
Fig.2.-Observed and predicted RVs.In the top panel, cyan lines show predictions of the Gaia astrometric solution.The center-of-mass RV is set to the value inferred from the joint-fit, but the predictions are otherwise uninformed by the observed RVs.The observed RVs broadly match the behavior predicted by the astrometric solution.In the middle panel, black lines show predictions of a joint fit of the RVs and the astrometric constraints.The combination of RVs and astrometry leads to tight constraints on the orbit.Bottom panel shows residuals of the data compared to the best-fit joint model.These are generally consistent with 0, indicating that the RVs are well-described by a Keplerian orbit.
Fig. 4.-Observation times of Gaia NS1 predicted by the Gaia observation scheduling tool (GOST).Black line shows the best-fit orbit from our combined fit.Red points show the predicted photocenter positions at the times when Gaia observed the source according to GOST.The actual measured ∆RA and ∆Dec are not published, only the predicted scan times.The expected per-epoch astrometric uncertainties are about 0.2 mas in the along-scan direction.The predicted scan times show that the orbit has been sampled well.
, and the extinction from theGreen et al. (2019) 3D dust map.We use the isochrones package 5(Morton 2015) to fit the MESAIsochrones and Stellar Tracks (MIST;Dotter 2016;Choi et al. 2016) library to the photospheric stellar parameters and multiwavelength photometry, parallax, and extinction data using MultiNest 6(Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019) ) via PyMultinest(Buchner et al. 2014).We measured EWs using the Spectroscopy Made Harder (smh) python package(Casey 2014), fitting unblended atomic transitions with Gaussian profiles.The atomic data for each line as well as the measured EWs are available upon request.We assume Asplund et al. (2021) photospheric solar abundances.The atmospheric parameters from the hybrid spectroscopic/isochrone fit are listed in Table3.We next inferred elemental abundances of O I, Na I, Mg I, K I, Ca I, Ti I, Ti II, Cr I, Cr II, Mn I, Fe I, Fe II, Ni I, Zn I, Y II, Zr II, and Ba II from their EWs, including isotopic/hyperfine splitting details where needed.We again measure the EWs by fitting Gaussian profiles using smh.We assume local thermodynamic equilibrium (LTE) and use the 1D plane-parallel, α-enhancedCastelli & Kurucz (2003) model atmospheres and the 2023 version of the LTE radiative transfer code MOOG(Sneden 1973) to infer elemental abundances based on our EWs.We report in Table Fig.5.-Spectral cutouts of Gaia NS1 (red) compared to a reference star (black), a halo star with similar stellar parameters and abundances observed by the GALAH survey.The two spectra are very similar, ruling out significant light contributions from a companion and validating the stellar parameters and abundances inferred in Section 6.

Fig. 6 .
Fig.6.-Left:observed spectrum of Gaia NS1 (red) compared to the reference star from the GALAH survey whose spectrum most closely matches Gaia NS1 (Figure5).Gaia NS1 is strongly enhanced in lithium compared to the reference star.Right: lithium abundance of Gaia NS1 compared to stars observed by GALAH with [Fe/H] < −0.7 and SNR > 30.Gaia NS1 is clearly an outlier in lithium abundance: it has higher A(Li) than any of the 1000+ stars with T eff within ±200 K.

Fig. 7 .
Fig. 7.-Galactic orbit of Gaia NS1, integrated backwards for 1 Gyr.The orbit of the Sun is shown for comparison.Gaia NS1 has an eccentric orbit that travels several kpc from the Galactic disk, characteristic of a halo star.

F
Fig.8.-Constraints on the mass of the unseen companion in Gaia NS1 as a function of the G−band flux ratio.Solid black line shows the constraint from the astrometric mass function, assuming M⋆ = 0.79 M ⊙ ; a completely dark companion would imply M companion ≈ 1.90 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 main-sequence companion.Because this is always below the black line, no single MS companion can explain the orbit.The same is true for an equalmass inner binary (yellow dashed) or an equal-mass inner triple (red dot-dashed).Even an inner binary containing a 1.3 M ⊙ WD and a MS star (blue dashed) line has a mass-to-light ratio that is too low to match the observed mass function.This leaves a single NS or a WD-WD binary as the only viable companions.
Fig.9.-Constraints on kicks and the pre-SN orbit of Gaia NS1.We assume that the binary had a circular orbit before the SN and that the NS received a kick with velocity v kick when it formed.We simulate 10 8 pre-SN orbits with a uniform period distribution P orb /d ∼ U (0, 1000), a pre-SN mass distribution M He star /M ⊙ ∼ U (1.9, 5.0), and a kick velocity distribution v kick /(km s −1 ) ∼ U (0, 100).We then show properties of the systems whose post-SN orbits have periods and eccentricities similar to Gaia NS1 (upper left).Most simulated systems that can produce the observed orbit have weak kicks (v kick ≲ 20 km s −1 ) and formed from low-mass He stars (M He star ≲ 3 M ⊙ ), though there is a long tail of higher kick velocities and more massive He stars that could produce the observed orbit for fine-tuned kick orientations.
Fig.10.-Predictedeccentricity of a 1.9 M ⊙ NS + 0.8 M ⊙ MS binary following a supernova.We assume the orbit is initially circular, that the NS forms from a helium star of mass M He star , and that it receives a natal kick of velocity v kick in a random direction at formation.Solid lines and shaded regions show median and middle 68% eccentricity regions of the surviving orbits.In the absence of a natal kick, loss of more than 0.3 M ⊙ leads to a higher eccentricity than observed.Kicks generally increase the eccentricity.The low observed eccentricity of Gaia NS1 thus suggests that the NS formed from a low-mass He star with a small natal kick.

v
Fig. 11.-Probability of remaining in a bound orbit with eccentricity e < 0.13, the value observed for Gaia NS1.As in Figure10, we assume that the orbit is initially circular, that the NS forms from a helium star of mass M He star , and that it receives a natal kick of velocity v kick in a random direction at formation.For very low v kick values, a low He star mass is required, because mass loss will otherwise result in an eccentric orbit.For higher v kick , a bound and near-circular orbit can be maintained if the kick is fortuitously aligned with the companion's instantaneous velocity vector, but the probability of this occurring is low.

Fig
Fig.12.-Masses and periods of NSs with well-characterized orbits and mass measurements compared to Gaia NS1.Cyan points show radio pulsars that are non-recycled neutron stars (mainly in double NS binaries and a few NS + massive WD binaries), while red points show recycled NSs (millisecond pulsars with low-mass WD companions).Compared to these systems, Gaia NS1 hosts one of the most massive NSs in the widest orbit.Its parameters are particularly unusual given that the NS presumably has not accreted from the luminous star.This suggests that the NS was born massive (it is not recycled, and is most similar to the cyan points).

Fig
Fig. 13.-Similar to Figure3, but for the case with the Gaia uncertainties inflated by a factor of 3. Constraints on all parameters are then consistent within 1σ between the astrometry only and astrometry+RVs fit.

TABLE 1
RVs.The best-fit instrumental offsets (lower block) have been subtracted so that all RVs are reported on the same scale. +0.01