Wide Binaries with White Dwarf or Neutron Star Companions Discovered from Gaia DR3 and LAMOST

The Gaia Data Release 3 (DR3) mission has identified and provided about 440,000 binary systems with orbital solutions, offering a valuable resource for searching binaries including a compact component. By combining the Gaia DR3 data with radial velocities from the LAMOST spectroscopic survey, we identify three wide binaries possibly containing a compact object. For two of these sources with a main-sequence companion, no obvious excess is observed in the blue/red band of the Gaia DR3 XP spectra, and the LAMOST medium-resolution spectra exhibit clear single-lined features. The absence of an additional component from spectral disentangling analysis further suggests the presence of compact objects within these systems. On the other hand, the visible star of the third source is a stripped giant star. In contrast to most binaries including stripped stars, no emission line is detected in the optical spectra. The unseen star could potentially be a massive white dwarf or neutron star, but the possibility of an F-type dwarf star scenario cannot be ruled out. An examination of about 10 binaries containing white dwarfs or neutron stars using both kinematic and chemical methods suggests most of these systems are located in the thin disk of the Milky Way.


INTRODUCTION
The evolution path a star takes in its life depends on its initial mass, metallicity, and rotational velocity, etc.As the final product of stellar evolution, compact objects are classified into three categories: white dwarfs, neutron stars, and black holes.In the era of multi-messenger astronomy, several techniques have been developed to search and identify compact objects, including X-ray observations (Remillard & McClintock 2006), gravitational wave detections (Abbott et al. 2016(Abbott et al. , 2017)), gravitational microlensing (Lam et al. 2022;Sahu et al. 2022), radial velocity (RV) monitoring (Casares et al. 2014;Liu et al. 2019;Thompson et al. 2019), and astrometry (El-Badry et al. 2023a,b).Each method possesses its own strengths and limitations.For instance, the X-ray method can be used to detect close binaries with strong accretion, whereas the gravitational wave method is limited to binaries with two compact objects.The radial velocity method is well-suited for the search of quiescent compact objects in binaries, while the astrometric method is more suitable for identifying compact objects in longperiod binaries.Recently, a group of black holes and neutron stars have been discovered through RV monitoring with large spectroscopic surveys (e.g., Liu et al. 2019;Thompson et al. 2019;Rivinius et al. 2020;Jayasinghe et al. 2021;Yi et al. 2022;Yuan et al. 2022), although some systems are still in debate.
Gaia spacecraft is a mission designed to explore the Milky Way by observing over a billion stars, utilizing its astrometry, photometry, and spectroscopy capabilities.The spacecraft is equipped with three main instruments (Gaia Collaboration et al. 2022): an astrometric instrument, prism photometers, and a Radial Velocity Spectrometer (RVS).The RVS has a resolution of approximately R ∼ 11,500, enabling the measurement of radial velocities of bright stars (G < 12 mag).The Gaia Data Release 3 (DR3) released four non-single star (NSS) tables for 813,687 objects classified as binaries or multiple stars, including nss two body orbit, nss acceleration astro, nss non linear spectro, and nss vim f l (van Leeuwen et al. 2022).Of particular interest is the nss two body orbit, which provides orbital solutions for astrometric, spectroscopic, and eclipsing binaries, flagged as "Orbital", "SB1", "SB2", "Astro-SpectroSB1", "EclipsingBinary", etc.This comprehensive catalogue serves as a rich reservoir for studying stellar multiplicity and analyzing various aspects of binary systems.
Large Sky Area Multi-Object Fiber Spectroscopic Telescope (hereafter LAMOST), also called GuoShou-Jing Telescope, is a specialized reflecting Schmidt telescope with an effective aperture of 3.6-4.9m and a field of view of 5 • (Wang et al. 1996).The focal plane is equipped with 4000 precisely positioned fibers that are connected to 16 spectrographs (Cui et al. 2012;Zhao et al. 2012).From 2011 to 2018, LAMOST conducted its first stage focused on a low-resolution (R ∼ 1,800) spectral survey.The low-resolution spectrum (LRS) covers a wavelength range of 3650-9000 Å.Since October 2018, LAMOST started its second 5-year survey program, containing both low-and medium-resolution (R ∼ 7500) spectral surveys.The medium-resolution spectrum (MRS) spans wavelength ranges of 4950-5350 Å and 6300-6800 Å (Liu et al. 2020).The LAMOST Data Release 9 (DR9) has released more than 19.47 million spectra for stars, galaxies, and quasars, etc. Notably, the time-domain spectral survey in the second stage provides a great opportunity to achieve a breakthrough in various scientific topics, such as binary systems, stellar activity, stellar pulsation, etc (Wang et al. 2021).
By combining the optical spectra from LAMOST with orbital solutions from Gaia DR3, we searched for binaries containing a compact object.As a result, we identified three potential candidates, where the visible stars are either giants or main-sequence stars, and the unseen companions are white dwarfs or neutron stars.In Section 2, we provided an overview of the sample selection and data reduction.Detailed information about the three systems, including stellar parameters of the visible stars, Kepler orbital solutions, and masses of the unseen companions, are presented in Sections 3 to 5. Finally, a short summary is provided in Section 7.

Source selection
The nss two body orbit table includes 443,205 binaries or multiple star systems.For spectroscopic binaries (flagged as "SB1", "SB2", or "AstroSpectroSB1"), it provides orbital parameters such as period P , eccentricity e, and semi-amplitude K of primary, etc.For "AstroSpectroSB1" binaries, the K value can be calculated from the astrometric solutions (c thiele innes and h thiele innes).First, we selected single-lined systems ("SB1" and "AstroSpectroSB1") and crossmatched them with the LAMOST DR9 low-resolution and med-resolution general catalog1 .Second, we chose sources with more than two LAMOST spectral observations with signal-to-noise ratio (SN R) greater than 5 and clear RV variation.Third, we calculated the binary mass function of those systems using orbital parameters from the nss two body orbit, and the gravitational masses of visible stars with multi-band magnitudes and atmospheric parameters from LAMOST parameter catalog.Finally, we estimated the minimum masses of the unseen stars using the mass functions and the evolutionary masses of the visible stars.
Following these steps, three binaries with possible compact objects were selected, namely Gaia

Data reduction
We obtained all available LRS and MRS observations from the LAMOST archive for our candidate sources.These spectra have undergone various reduction steps with the LAMOST 2D pipeline, including bias and dark subtraction, flat field correction, spectrum extraction, sky background subtraction, and wavelength calibration, etc.(See Luo et al. 2015, for detalis).Figure 1 displays the LRS observations for the three targets.The cross-correlation function (CCF) was used to calculate the RV values using the red band (6300-7000 Å) of each MRS spectrum with SN R > 5. Furthermore, for each spectrum, we derived a calibration factor (∆rv = rv r1 − rv r0) with the rv r1 and rv r0 from the LAMOST MRS catalog.Here rv r1 represents the calibrated velocity of rv r0 by using RV standard stars (Huang et al. 2018).The final RV values were obtained by summing the measured value from the CCF method and the calibration factors (Table A.1).
In addition, we applied for spectral observations using the Beijing Faint Object Spectrograph and Camera (BFOSC) mounted on the 2.16 m telescope at the Xinglong Observatory and the Double Spectrograph (DBSP) mounted on the Palomar's 200 in.telescope (P200).For G8441, we carried out eight observations using the 2.16 m telescope from Feb. 20th, 2019 to Apr. 21th, 2019, with the E9/G10 grism and a 1.6 ′′ slit configuration, and four observations using the P200 telescope on Mar.14th, 2019.For G3431, we applied for two observations using the 2.16 m telescope on Oct. 4th, 2022.The observed spectra were reduced using the IRAF v2.16 software (Tody 1986(Tody , 1993) ) following standard steps, and the reduced spectra were then corrected to vacuum wavelength.The CCF method was used to calculate RV, and the barycentric corrections to the observation time and RV were made using the light travel time and radial velocity correction functions provided by the Python package Astropy (Astropy Collaboration et al. 2013, 2018, 2022).

Stellar parameters of the visible star
Each of our targets has been observed multiple times by LAMOST.Their atmospheric parameters can be estimated following (Zong et al. 2020): and where the index k is the k th epoch of the measurements of parameter P (i.e., T eff , logg, and [Fe/H]) for each star, and the weight w k is square of SNR of each spectrum corresponding to the k th epoch.(Luo et al. 2015), DD-Payne (Xiang et al. 2019;Ting et al. 2019), and CNN (Wang et al. 2020a).In the following analysis, we preferred to use the LASP LRS results followed by the LASP MRS results.The stellar parameters of G4031 are T eff = 5953±17K, logg = 4.05±0.02dex, and [Fe/H] = −0.37±0.01(Table 2).The distance from Gaia DR3 is approximately 720 pc.The Pan-STARRS DR1 3D dust map2 does not provide an extinction estimation, instead we used the SFD value with E(B − V ) = 0.02.The stellar parameters can also be obtained by fitting the spectral energy distribution (SED).We used the Python module ARIADNE3 to estimate the stellar parameters of visible star, which can automatically fit the broadband photometry by using different stellar atmosphere models, such as Phoenix4 , BTSettl5 , Castelli & Kurucz6 , and Kurucz 19937 .We constructed the SED using magnitudes from various surveys: Gaia DR3 (G, G BP , and G RP ), 2MASS (J, H, and K S ), APASS (B, V , g, r, and i) and WISE (W 1 and W 2). The Gaia DR3 parallax and extinction A V were also used as input priors.From the SED fitting (Figure 2), we obtained stellar parameters for a G-type star with an effective temperature of 5910 +36 −29 K, surface gravity of 4.05 +0.01 −0.01 dex, metallicity of −0.37 +0.01 −0.01 (Table 3), which are in good agreement with spectroscopic results.
Furthermore, we calculated stellar mass using two methods: one method involved the use of stellar evolution models, referred to as evolutionary mass; the other method utilized observed photometric and spectroscopic parameters, referred to as gravitational mass.
First, we used the isochrones, a Python module (Morton 2015), to calculate evolutionary mass of G4031 by   perature to calculate the gravitational mass following M = gR 2 /G (See details in Li et al. 2022).The final mass and corresponding uncertainty were calculated as the average values and standard deviations of the masses derived from different bands.The gravitational mass of G4031 is 1.62±0.02M⊙ .

Kepler orbit
We averaged the Barycentric Julian Dates (BJDs) and RVs taken on the same day (Appendix A) for the purpose of RV fitting.The Joker (Price-Whelan et al. 2017), a custom Markov chain Monte Carlo sampler, was employed to fit the Keplerian orbit.Figure 2 displays the phase-folded RV data along with the best-fit RV curve.Table 4 presents the orbital parameters of G4031 obtained through The Joker fitting, including orbital period P , eccentricity e, argument of periastron ω, mean anomaly at the first exposure M 0 , RV semiamplitude K, and systematic RV ν 0 .The fitting results differ slightly from those obtained from the table nss two body orbit of Gaia DR3.
The binary mass function can be calculated following where M 2 is the mass of the unseen star, q = M 1 /M 2 is the mass ratio of this system, and i is the systematic inclination angle.We obtained a mass function of ≈0.09 M ⊙ .By using the evolutionary mass of 1.46 +0.02 −0.02 M ⊙ and gravitational mass of 1.62±0.02M⊙ of the visible star, we determined the minimum mass of the unseen object as 0.77±0.03M⊙ and 0.82±0.02M⊙ , respectively.Using the orbital parameters from Gaia DR3, the mass function is ≈0.16 M ⊙ while the minimum mass of the unseen star would be 0.98 ± 0.09M ⊙ or 1.04 ± 0.09M ⊙ .These results suggest G4031 may contain a white dwarf.
G4031 has been observed by Zwicky Transient Facility (ZTF; Bellm et al. 2019) in g band and All-Sky Automated Survey for Supernovae (ASAS-SN; Kochanek et al. 2017) in V band.No clear variation caused by ellipsoidal deformation can be seen in the folded light curves (Figure 3).Therefore, no further information  such as the orbital inclination can be constrained by the light curve.

The nature of the unseen object
We calculated the size of Roche-lobe by using the standard Roche-lobe approximation (Eggleton 1983), where a = (1 + q)a 1 = (1 + q)P (1 is the separation of binary system.The filling factor is ≈ 3.2% assuming i = 90 • , indicating this system is a detached binary system.For the visible star, the absolute magnitude in the G band is G = 3.04 mag; while for a main-sequence star with a mass of 0.77M ⊙ or 0.98M ⊙ , the absolute magnitude in the G band would be around 6.20 mag or 4.80 mag.By comparing the flux-calibrated Gaia XP spectrum (Huang et al. 2023, in preparation) with the Phoenix template, no flux excess is observed (Figure 2), indicating the companion of G4031 is a compact object if the secondary has a mass of 0.98M ⊙ .However, if the secondary has a mass of 0.77M ⊙ , it cannot be distinguished based on the SED or spectral continuum.Given this uncertainty, we decided to investigate the binary nature of G4031 with the LAMOST MRS observations by using the spectral disentangling method.
We employed the algorithm of spectral disentangling proposed by Simon & Sturm (1994).Before applying this algorithm, we performed preliminary tests to determine the detection limits for our targets.We derived the effective temperature and surface gravity of the secondary (assuming a main-sequence star) using the mini-mum mass8 following different mass ratio (Table 5).The corresponding Phoenix9 model was used as the theoretical spectra of the secondary.Then, we combined the observed spectra and theoretical spectra, with a reduced resolution of R = 7500 and different rotational broadening (vsini = 10 km/s, 50 km/s, 100 km/s and 150 km/s), to produce synthetic binary spectra.For the simulated binary spectra, the wavelength range of 6400 Å to 6600 Å was used for spectral disentangling.The results by visual check are presented in Table 5.It can be seen that for G4031, when vsini is less than 100 km/s, this method can successfully separate the spectra of primary and secondary components.
The MRS observations within the wavelength range of 6400 Å to 6600 Å were utilized for spectral disentangling.As an example, Figure 4 shows the results for a mass ratio of 1.4 and 1.8.No feature of absorption spectrum is apparent for an additional component (red lines in Figure 4), implying that the component star is a compact object.
We also used single-and binary-star spectral models (Liu et al. 2023, in preparation) to perform fitting for the LRS observations.In brief, the spectra from LAMOST LRS and the stellar parameters from the Apache Point Observatory Galactic Evolution Experiment (APOGEE) DR16 were used as the training set to develop a model using the neural network version of Stellar LAbel Machine (SLAM) (Zhang et al. 2020a,b) for single-star spectra.The binary-star model was then created by combining two single-star models with appropriate radial velocities.Spectral observations were fitted with both the single-and binary-star models, and their corresponding χ 2 values were compared to determine the best fit.The results of the semi-empirical spectroscopic experiments show that one system is a single star if the difference in χ 2 (i.e., ∆χ 2 ) between the single-star and binary-star fitting is less than 0.4.For G4031, the ∆χ 2 value is about −0.04, further evidencing that it includes one compact object.

Stellar parameters of the visible star
The visible star of G3431 is a G-type main-sequence star.The stellar parameters from the LASP (LRS) method are used in the following analysis, with an effective temperature of T eff = 6402±13K, surface gravity of logg = 4.23±0.02dex, and metallicity of [Fe/H] = −0.06±0.01(Table 2).According to Gaia DR3, the  −0.01 (Table 3 and Figure 5).The isochrones fitting based on evolutionary models yields an estimated evolutionary mass and radius for G3431 of 1.42 +0.03 −0.03 M ⊙ and 1.83 +0.03 −0.03 R ⊙ , respectively.On the other hand, the gravitational mass of G3431, obtained by averaging six gravitational mass estimations, is determined to be 2.06±0.13M⊙ , which differs from the evolutionary mass estimation.It is important to note that even within the gravitational mass measurements, there are discrepancies between the results derived from Gaia bands (≈ 1.88 M ⊙ ) and those from 2MASS bands (≈ 2.24 M ⊙ ).By reviewing the calculation process, we found that the bolometric magnitudes from Gaia bands (BP , G, RP ) are 3.13, 3.09, 3.04, while the bolometric magnitudes from 2MASS bands (J, H, K S ) are 2.93, 2.89, 2.86.This gradual variation with wavelength can be attributed to the measurement accuracy of stellar parameters, such as T eff and E(B − V ).We noticed that the difference in temperature between different methods is about 250 K (Table 1), and the temperature from SED fitting is lower than spectral results.In this case, the mass measurement from stellar model fitting is more reliable than the mass calculated with atmospheric parameters.Additionally, the evolutionary mass of ≈1.42 M ⊙ is more reasonable for a late-F main-sequence star.
Taking all these factors into consideration, we adopted the evolutionary mass measurement for further analysis.

Kepler orbit
The Barycentric Julian Dates (BJDs) and RV values taken on the same day (Appendix A) were averaged to improve the accuracy of the RV fitting.The phasefolded RV data and the best-fit RV curve from The Joker are shown in Figure 5. Based on the RV fitting results, we obtained a mass function of ≈0.32 M ⊙ and a minimum mass of 1.36 ± 0.09M ⊙ for the unseen star, which are slightly lower than the values calculated from Gaia DR3 solution.As G4031, no clear variation can be seen in the light curve of G3431 (Figure 6).Using the Roche-lobe approximation (Eq.4), we estimated that the filling factor of G3431 is approximately 3.4% (assuming i = 90 • ), suggesting this system is a detached binary system.The possibility that G3431 contains two main-sequence stars can be ruled out according to the luminosity ratio.For instance, the absolute G-band magnitude of the visible star is G = 3.06 mag, while the absolute magnitude of a main-sequence star with a mass of 1.36M ⊙ is G = 3.10 mag.

The nature of the unseen object
We also applied the spectral disentangling technique to the observed spectra of G3431.The results of the spectral disentangling tests (Table 5) indicate that for binary systems with two normal stars having a mass set similar to G3431, this method can successfully separate the binary components including two normal stars with high significance when vsini is less than 150 km/s.No additional component with an optical absorption spectrum can be detected (Figure 7), excluding the possibility of a normal binary system.Furthermore, the ∆χ 2 values of G3431 from single-star and binary fitting are about −0.18 and −0.24 for its two LRS observations, indicating that these spectra are more likely from a single star.

Stellar parameters of the visible star
The stellar parameters of G8441 obtained from LASP (LRS) are as follows: T eff = 4194±2 K, logg = 1.83±0.03dex, and [Fe/H] = −0.74±0.01(Table 2).According to Gaia DR3, the estimated distance to G8441 is approximately 857 pc, and the Pan-STARRS DR1 3D dust map suggests an extinction of E(B − V ) ≈ 0.01 at this distance.SED fitting using the ARIADNE package (Figure 8a) yields an effective temperature of 4144 +23 −20 K, surface gravity of 1.86 +0.11 −0.09 dex, and metallicity of −0.77 +0.10 −0.10 , which are in good agreement with the spectroscopic estimation.
Unlike G4031 and G3431, the light curves of G8441 exhibit significant ellipsoidal variations (Section 5.2).The large amplitude of variation in V -band light curve (≳ 0.4 mag) implies that the visible star of G8441 is (close to) Roche lobe-filling.As a result, the single-star evolution model is unsuitable for assessing the mass of the visible star.For a stripped star with Roche lobe overflow, its density can be calculated following (Frank et al. 2002), ρ = 110P −2 hr g cm −3 . (5) Using the radius from SED fitting (R = 16.55R ⊙ ; Table 3), we calculated the mass of the visible star to be M = 4π ρR 3 /3 ≈ 0.28±0.03M ⊙ .

Kepler orbit
The phase-folded RV data and the best-fit RV curve from The Joker are shown in Figure 8.The RV fitting returns a mass function of ≈0.84 M ⊙ and a minimum mass of 1.25 ± 0.09 M ⊙ for the unseen star, while using the orbital parameters of Gaia DR3 (Table 4), the mass function is ≈0.67 M ⊙ and the minimum mass of the unseen star is 1.07 ± 0.10 M ⊙ .
The Lomb-Scargle algorithm (Lomb 1976;Scargle 1982) was used to calculate the orbital period (≈46.885days) with the ASAS-SN g-band light curve (Figure 9).To obtain a more precise period estimation, we further folded the light curves with periods ranging from 46.5 to 47 days, in a step of 0.001 day.Through visual inspection of the folded light curves, we determined a period of ≈46.895 days.Figure 10 displays the folded light curves in multiple bands.
We employed the software Physics of eclipsing binaries (Prša et al. 2016;Horvat et al. 2018;Conroy et al. 2020) (PHOEBE) to fit the ASAS-SN V -and g-band light curves.The ellipsoidal modulation (Figure 10) and the substantial variation amplitude indicate the visible star is a stripped star.During the LC fitting, we used the eclipse method = only horizon as the eclipse model and specified a semi-detached system for the binary configuration.Two scenarios were considered: one with a stripped star and a compact object, and the other with a stripped star and a standard main-sequence star.In the first case, a cold (T eff = 300 K) and small (R = 3×10 −6 R ⊙ ) blackbody was used as the secondary.
We used two different sets of asini values: one from The Joker fitting and the other from the Gaia DR3 solution (Table 4).The bolometric gravity darkening coefficients was adopted as β = 0.58, which was derived from Claret & Bloemen (2011) in the V band for stars with similar atmospheric parameters.Table 6 presents the parameter estimates from the PHOEBE fitting.For the scenario of a stripped star and a compact object, the fitting results show that the unseen star is a neutron star or a massive white dwarf, while for the scenario of a stripped star and a main-sequence star, the best-fit models indicate the unseen dwarf is an A-or F-type star.

The nature of the unseen object
G8441 has been reported as a possible compact object candidate in previous studies (Gu et al. 2019;Zheng et al. 2019).In this section, we will try to investigate the nature of the unseen object by examining the optical spectra, and the UV and X-ray emission.Due to the limited number of spectral observations, we did not perform the spectral disentangling.
No emission line is detectable in the optical spectra (Figure 1).It is evident that the H α lines in all spectral observations are absorption lines (Figure 11), indicating the absence of an accretion disk.This contrasts with most binaries containing a stripped giant star, such as V723 Mon and 2M04123153+6738486 (Jayasinghe et al. 2021(Jayasinghe et al. , 2022)), which display clear H α emission lines.G8441 is likely similar to a few binaries without emission lines, explained by temporary halted or slowed accretion (El-Badry & Rix 2022).
G8441 has been observed in FUV and NUV bands by GALEX.By using the UV magnitudes of f FUV = 23.28 µJy and f NUV = 1142.51µJy, we calculated the luminosities as L FUV = 1.16 × 10 31 erg s −1 and L NUV = 6.19 × 10 32 erg s −1 .This luminosity can be attributed to the giant star, given that K-type giants typically exhibit FUV and NUV luminosities within the ranges of ≈ 5 × 10 29 -5 × 10 32 erg s −1 and ≈ 3 × 10 29 -7 × 10 33 erg s −1 , respectively (Wang et al. 2020b).However, it is uncertain whether a stripped star has the same magnetic dynamo mechanism as a normal giant.If the the unseen Figure 10.Folded light curves (from top to bottom: ASAS-SN g, ASAS-SN V , TESS) of G8441 by using the orbital period of 46.895 day.The left panels show the fitting for the scenario of a stripped star and a compact object, while the right panels show the fitting for the scenario of a stripped star and a main-sequence star.The red and green lines represent the best-fitting results from PHOEBE, using the orbital solution from The Joker and Gaia DR3, respectively.
star is an A-type star (T eff ≈ 8200 K from PHOEBE fitting), its FUV and NUV emission would exceed the observed levels significantly (Figure 12).Therefore, the A-type star scenario can be ruled out.If the unseen star is an F-type star (T eff ≈ 7200 K from PHOEBE fitting), its UV emission are compatible with the observation.On the other hand, if the unseen object is a white dwarf, the upper limit of its surface temperature can be constrained with the FUV emission.White dwarfs with effective temperatures greater than 20,000 K can be ruled out since models predict higher fluxes than the GALEX FUV observation (Figure 12).Thus, the UV luminosities (especially the high NUV luminosity) can be either totally from the stripped star, or from a combination of the stripped star and the unseen object.
G8441 was observed by XMM-Newton on May 4, 2020.No X-ray counterpart was detected from the EPIC instruments, including M1, M2, and PN.Using a radius of 10 ′′ , an upper limit of the net count rate in the 0.3-8 keV range was calculated with the PN instrument, yielding approximately 0.0004 counts/sec.We utilized the PIMMS tool10 to convert the count rate into flux.Applying the standard linear relation between the hydrogen column density N H and the reddening N H = 5.8 × 10 21 /E(B − V ) (Schlegel et al. 1998), the hydrogen column density N H was calculated to be 5.8 × 10 19 cm −2 with the optical reddening of E(B − V ) = 0.01 mag.By using a power-law photon index of 1.5 or 2, the upper limit of the unabsorbed flux in the 0.3-8 keV band was derived to be 1.1 × 10 −15 erg cm −2 s −1 or 7.9 × 10 −16 erg cm −2 s −1 , respectively.At the adopted distance of 857 pc, the upper limit of the emitted luminosity was estimated to be L X ≲ 9.7 × 10 28 erg s −1 or 6.9 × 10 28 erg s −1 .
To sum up, the substantial variation amplitude of light curves indicates the visible star is most likely a stripped star, although there is no emission feature in the optical spectra.The unseen object could be a compact object (white dwarf or neutron star) or a main-sequence F star.Unfortunately, present observational data (e.g., UV and X-ray data, SED, and optical spectra) do not provide sufficient evidence to distinguish between these scenarios.

Evolutionary track
The Binary Population And Spectral Synthesis code (BPASS) (Eldridge et al. 2017;Stanway & Eldridge 2018) is a tool to calculate the evolutionary tracks of single stars or binary systems.We used the hoki package (Stevance et al. 2020) to search for binary models consistent with the observed properties of G8441 based on the results from BPASS.The selection criteria included and M 2 in the range of 1.1-1.5 M ⊙ , where M 1 and M 2 represent the masses of the giant and the invisible star, respectively.We identified four systems including an F-or G-type star and one system with a compact object (Table 7).The parameters of these systems are not entirely consistent with the analysis in Section 5.2 and 5.3.For example, the masses of the stripped star in these models are higher than our calculation (≈0.28 M ⊙ ), and for the four models containing normal stars, the temperatures of these normal Zhao et al. star are lower than the PHOEBE fitting results (Table 6).
One interesting finding is that the mass of the stripped star (≈0.28 M ⊙ ) follows the relationship between the mass and orbital period (M WD /M ⊙ = ( P orb /day 1.1×10 5 ) 1/4.75 + 0.115) of white dwarfs or the core mass of low-mass giants (Tauris & Savonije 1999).This suggests that the mass transfer in G8441 might have been close to termination and most of the hydrogen envelope of the stripped star has likely been stripped away.This donor star could eventually evolve into an extremely low mass white dwarf.

SPATIAL DISTRIBUTION IN THE MILKY WAY
Until now, about ten binaries including white dwarfs or neutron stars have been discovered by RV.We collected these binary systems (Mazeh et al. 2022;Yuan et al. 2022;Li et al. 2022;Yi et al. 2022;Zhang et al. 2022;Zheng et al. 2022a,b) to investigate their positions in the Milky Way.
To identify whether one source belongs to the thin disk, thick disk, or halo, first, we defined one probability assuming that the velocities (U LSR , V LSR , W LSR ) follow a 3-D Gaussian distribution (Guo et al. 2016 where c = (2π) −3/2 (σ U σ V σ W ) −1 normalizes the expression.Here V asym is the asymmetric drift, and σ U , σ V , and σ W are the velocity dispersions in three dimensions, all of which vary for different components 11 .Second, we defined the probability ratio of belonging to different components as One star was classified as a candidate in the thin disk, thick disk, or halo when the corresponding ratio exceeded 50%.8).

SUMMARY
By combining the Gaia DR3 astrometric solution and LAMOST DR9 LRS and MRS data, we discover three binaries with possible compact components (i.e., G4031, G3431 and G8441).
G4031 is a binary system with a period of P ≈ 140 day consisting of a G-type main-sequence star.The evolutionary and gravitational masses of the visible star are estimated to be 1.46 +0.02  −0.02 M ⊙ and 1.62 ± 0.02M ⊙ , respectively.The Kepler solution from The Joker matches better with the RV data than the orbit from Gaia DR3.By using The Joker results, the mass function is calculated to be f (M ) ≈ 0.092 ± 0.006M ⊙ .These lead to a minimum mass of the unseen star of 0.77 ± 0.03M ⊙ (calculated with the evolutionary mass of the visible star) or 0.82 ± 0.02M ⊙ (calculated with the gravitational mass of the visible star), suggesting the unseen star is likely a white dwarf.
G3431 has a orbital period of P ≈ 120 day containing a G-type main-sequence star.Using stellar evolution models, the mass of the visible star is determined to be M = 1.42 +0.03 −0.03 M ⊙ .The orbit solutions obtained from The Joker and Gaia DR3 are in good agreement, indicating a mass function of f (M ) ≈ 0.324 ± 0.038M ⊙ or 0.392 ± 0.023M ⊙ , respectively.Therefore, an invisible star with a minimum mass of 1.36 ± 0.09M ⊙ or 1.49 ± 0.05M ⊙ is obtained, suggesting the presence of a white dwarf or neutron star.
G8441 is a binary with a period of P ≈ 46.895 day.Different with G4031 and G3431, the light curves of G8441 strongly suggest that the visible star is a stripped star.Using the period derived from Gaia DR3, we calculated the mass of the stripped star to be 0.28±0.03M ⊙ .The RV fitting indicates a mass function of f (M ) ≈ 0.84 ± 0.12M ⊙ and a minimum mass of 1.25 ± 0.17M ⊙ for the unseen star.By fitting the multi-band light curves, we determined the orbital inclination angle to be around 70 degrees.The mass range of the invisible star is 1.1-1.5 M ⊙ , implying it is possibly a compact object (white dwarf or neutron star) or an F-type mainsequence star.Although the visible star is recognized as a stripped star, there are no detectable emission lines in the optical spectra.This intriguing feature may be attributed to a temporary interruption or reduction in the accretion process.
For the two sources with a main-sequence companion, by examining the Gaia XP spectra and applying spectral disentangling to the LAMOST MRS observations, we have not identified any additional optical component in addition to the visible star, supporting the scenario that these sources are binary systems with compact components.An examination of about ten binaries containing  white dwarfs or neutron stars using both the kinematic and chemical methods suggest most of these systems are located in the Galactic thin disk.
RV and astrometry have opened up new possibilities for searching for binaries with compact components.In addition to LAMOST, the RV data from other large spectroscopic surveys, including APOGEE, Radial Velocity Experiment (RAVE), GALactic Archaeology with HERMES (GALAH), and Dark Energy Spectroscopic Instrument (DESI), should also be collected and utilized.By combining the epoch RV data and astrometric data from Gaia, particularly with epoch astrometric measurements from Gaia DR4 in the future, numerous compact objects hidden in wide binaries can be discovered.A comprehensive sample of compact objects with accurate dynamical mass measurements will significantly contributes to building the mass function of compact objects.Especially, a sample of massive neutron stars and low-mass black holes can help confirm the maximum mass of neutron stars and minimum mass of black holes, and help solve the mass gap problemthe scarcity of black holes with masses ranging from 2.5 to 5 M ⊙ (Bailyn et al. 1998).These discoveries, along with the compact objects found by other methods (e.g., X-ray, gravitational wave, and gravitational microlensing), are expected to provide valuable insights into the late evolution of massive stars (in binary systems) and help us better understand the supernova mechanism and the black hole formation process.We thank the anonymous referee for helpful comments and suggestions that have improved the paper.We thank Dr. Kareem El-Badry for very useful discussion about binaries including stripped stars.The Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences.Funding for the project has been provided by the National Development and Reform Commission.LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.This work uses data obtained through the Telescope Access Program (TAP), which has been funded by the TAP member institutes.This work presents results from the European Space Agency (ESA) space mission Gaia.Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC).Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Mul-tiLateral Agreement (MLA).The Gaia mission website is https://www.cosmos.esa.int/gaia.The Gaia archive website is https://archives.esac.esa.int/gaia.We acknowledge use of the VizieR catalog access tool, operated at CDS, Strasbourg, France, and of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013).This research made use of Photutils (Bradley et al. 2020)

Figure 1 .
Figure 1.Left panel: Position of visible stars on the Hertzsprung-Russell diagram for G8441 (blue star), G3431 (red triangle) and G4031 (blue triangle).The gray points are from the Gaia DR3 with distance d < 200 pc, Gmag between 4-20 mag, and galactic latitude |b| > 10.No extinction correction for those stars.Right panel: LAMOST LRS observations of G4031, G3431 and G8441.

Figure 2 .
Figure 2. Panel a: SED fitting of G4031.Panel b: Folded LAMOST RV data and the RV curve from The Joker.The blue dots are the RV data from LAMOST MRS spectra, and the green dots are the RV data from LAMOST LRS spectra.Panel c: Folded LAMOST RV data and the RV curve from Gaia nss two body orbit solution.Panel d: Comparison of flux-calibrated Gaia XP spectrum and the Phoenix template (T eff = 6000K, logg = 4 dex, [Fe/H]= −0.5).

Figure 3 .
Figure 3. Folded ASAS-SN V -band and ZTF g-band light curves of G4031 with the orbital period of ≈137.98 day.

Figure 4 .
Figure 4. Left panel: spectral disentangling of G4031 with q = 1.4.Right Panel: spectral disentangling with q = 1.8.The vertical panels show spectra in different phases (close to the minimum or maximum RV phase of the visible star).The blue lines mark the reconstructed spectra of the visible star, while the red lines represent the second component in each spectra.The green lines are the sum of the two components, and the black lines represent the observed spectra.estimated distance to G3431 is approximately 338 pc.The Pan-STARRS DR1 3D dust map indicates an extinction with E(B − V ) = 0.03 at this distance.SED fitting returns stellar parameter estimations for the visible star, with T eff = 6156 +69 −39 K, logg = 4.22 +0.04 −0.03 dex, [Fe/H] = −0.06+0.01 −0.01 (Table3and Figure5).The isochrones fitting based on evolutionary models yields an estimated evolutionary mass and radius for

Figure 5 .
Figure 5. Panel a: SED fitting of G3431.Panel b: Folded LAMOST RV data and best-fit RV curve from The Joker.The blue and green dots are the RV data from LAMOST MRS and LRS observations, and the black star marks the RV data from 2.16 m observation.Panel c: Folded RV data and the RV curve from Gaia nss two body orbit solution.Panel d: Comparison of the flux-calibrated Gaia XP spectrum and the Phoenix template spectrum (T eff = 6200K, logg = 4.5 dex, [Fe/H]=0.0)for G3431.

Figure 6 .
Figure6.Folded ASAS-SN V -band light curve of G3431 by using the orbital period of ≈120.56 day.

Figure 8 .Figure 9 .
Figure 8. Panel a: SED fitting of G8441.Panel b: Folded LAMOST RV data and the RV curve from Gaia nss two body orbit solution.The green and blue dots are the RV data from LAMOST LRS and MRS, respectively.The purple square marks the RV data from P200, while the black stars represent the RV data from 2.16 m telescope.

Figure 11 .
Figure 11.Profiles of the RV-corrected Hα lines of G8441.

Figure 12 .
Figure12.SED of G8441.Top panel: The black points are the multi-band observations of G8441.The red line represents the spectrum of the stripped star, while the blue line marks the spectrum of an A-type main-sequence (T eff ≈ 8200 K from PHOEBE fitting).Middle panel: The blue line represents the spectrum of an F-type dwarf (T eff ≈ 7200 K from PHOEBE fitting).Bottom panel: The blue line represents the spectrum of a white dwarf with a mass of 1.2M⊙ and a temperature of 20,000 K.
Figure 13 (left panel) shows the Toorme diagram of those binary systems.Furthermore, we collected the [α/Fe] and [Fe/H] values from Wang et al. (2022) and segregated these systems into the thin disk and thick disk based on the [α/Fe]-[Fe/H] diagram (Figure 13).Both the kinematic and chemical methods suggest most of these systems belong to the Galactic thin disk (Table

Figure 13 .
Figure 13.Spatial distribution of binaries including white dwarfs or neutron stars.Left panel: The Toorme diagram of the binary systems.The circles and triangles represents stars classified as in the thin disk or thick disk using the kinematic method.Right panel: The [Fe/H] -[α/Fe] diagram of the binary systems.The green dashed line indicates the approximate border of thin disk and thick disk from metallicity (Masseron & Gilmore 2015).

Table 1 .
Estimation of atmospheric parameters with different methods (LASP, DD-Payne, CNN) for our systems.The LASP results are from the LAMOST DR9 parameter catalog.

Table 2 .
Stellar parameters for our systems, including atmospheric parameters, distance, extinction and multi-band photometric magnitudes.BP , G RP , J, H, and K S ) and the effective tem- Zhao et al.

Table 4 .
Keplerian orbit fitting results from the The Joker and Gaia DR3 table nss two body orbit.The minimum mass of secondary (M2,min) is calculated by using the evolutionary mass of the visible star.

Table 5 .
Spectral disentangling test results by visual check."✔✔" means the spectra are well disentangled, "✔" means the spectra can be disentangled but in low significance, and "✕" means the spectra can not be disentangled.

Table 6 .
Parameter estimates from PHOEBE for G8441.

Table 8 .
Spatial distribution parameters of binaries including white dwarfs or neutron stars.
a The classification is from the kinematic and chemical methods.