Two Dynamically Discovered Compact Object Candidate Binary Systems from LAMOST Low-resolution Survey

We report two binary systems, LAMOST J035540+381550 (hereafter J035540) and LAMOST J035916+400732 (hereafter J035916), identified through the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) low-resolution survey (LRS). Each of these two systems contains an M-type star orbiting with a invisible compact object candidate. Follow-up spectroscopic observations of Palomar 200-inch telescope (P200) enhance radial velocity measurements. We use radial velocities from LAMOST and P200, as well as light curves from Zwicky Transient Facility (ZTF) to constrain orbital parameters. The masses of the visible M-type stars are estimated by fitting the MIST isochrones and SEDs. The mass functions for the unseen companions are: $0.22\pm 0.01 M_{\odot}$ for J035540 and $0.16\pm 0.01 M_{\odot}$ for J035916. With the orbital and stellar parameters derived above and assuming the orbital inclination is 90 degree (edge-on), we find that the minimum masses of the invisible companions exceeds that of the visible stars. The single-lined feature and the dynamical evidence suggest the presence of compact objects. J035540's ZTF light curve, modeled with PHOEBE, yields a compact object mass of $0.70^{+0.12}_{-0.05} M_{\odot}$. For J035916, ellipsoidal modulation analysis constrains the light curve amplitude, yielding a compact object mass range of $0.57-0.90 M_{\odot}$. The mass estimates indicate that both are likely white dwarfs. These findings underscore the efficiency of optical time-domain surveys and dynamical methods in identifying faint, massive white dwarfs, along with other compact objects in binaries.


INTRODUCTION
White dwarfs (WDs), neutron stars (NSs), and black holes (BHs) are the products of stellar evolution in the final stages of the life cycle.Numerous observations indicate the presence of a substantial number of compact objects in binary systems (Joss & Rappaport 1984;Remillard & McClintock 2006;Rebassa-Mansergas et al. 2012).The searching and identification of these compact objects within binary systems play a crucial role in advancing our understanding of stellar formation, evolutionary history, binary interactions and evolution, as well as the properties and behavior of matter under extreme densities.Concurrently, compact object binary systems are intricately linked to a plethora of high-energy astrophysical phenomena, including novae, supernovae, X-ray bursts, and gravitational waves (Paczyński 1971;Whelan & Iben 1973;Paczynski 1976;Webbink 1984;Nomoto et al. 1984;Abbott et al. 2016Abbott et al. , 2017)).
In recent years, with the rapid development of optical time-domain surveys, the dynamical searching and identification of binary systems containing compact objects have witnessed a significant upsurge.The dynamical method relies on multi-epoch spectral observations to identify stars exhibiting significant radial velocity variations caused by compact companion stars (Trimble & Thorne 1969).Recently, several stellar-mass BHs and NSs in binary systems have been uncovered ★ E-mail: guwm@xmu.edu.cnutilizing the dynamical method based on astrometric, spectroscopic, and photometric time-domain surveys (Thompson et al. 2019;Liu et al. 2019Liu et al. , 2020;;Jayasinghe et al. 2021;El-Badry et al. 2022a;Yi et al. 2022;Tanikawa et al. 2023;El-Badry et al. 2023;Zheng et al. 2023).In the pursuit of non-interacting binary BH systems, numerous false positives have been confirmed to be actual luminous binary systems (El-Badry & Quataert 2021;El-Badry et al. 2022b;El-Badry & Rix 2022).Simultaneously, the nature of certain NS candidates identified through dynamical methods remains uncertain.This uncertainty arises from measured masses that do not surpass the Chandrasekhar limit, leaving ambiguity regarding whether they are NSs or massive WDs (Lin et al. 2023;Zhang et al. 2024).
The dynamical method offers a robust approach to search for binary systems with WDs.Rowan et al. (2024, Table 6) presents a collection of WD candidates identified using this method.Their table lists 12 WD candidates, with visible companions that include G-, K-, and M-type stars.Over the past few decades, the sample of white dwarf-main sequence (WD-MS) binaries has significantly increased.The majority of these binaries are identified through the prominent hydrogen and helium absorption lines characteristic of WDs in their blue-end spectra, which are more easily detectable when the WDs are hot (Rebassa-Mansergas et al. 2012;Li et al. 2014;Bar et al. 2017;Ren et al. 2018).Due to the intrinsic dimness of WDs, search methods relying on WD's unique spectral components are generally confined to nearby sources.Furthermore, a substantial population of young (≲ 1 Gyr) and hot (≳ 10000 K) WDs has been discovered leveraging UV excess (Parsons et al. 2016;Hernandez et al. 2021).However, the galactic population of massive WDs (≳ 0.8 ⊙ ; Hernandez et al. 2022b) requires additional supplementation.Measuring the mass distribution of WDs holds crucial implications for understanding pulsating phases on the asymptotic giant branch (Williams et al. 2009;Cummings et al. 2018), stellar population synthesis (Maraston 1998), and galaxy evolution theories (Catalán et al. 2008;Kalirai et al. 2014).Massive WDs have a more compact structure, resulting in lower luminosity, which makes them challenging to detect.Therefore, dynamical methods present an opportunity to expand the sample of massive WDs.
The Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) is a four-meter-aperture, reflective Schmidt telescope with a wide field of view (∼ 5 • ) (Cui et al. 2012).With 4000 fibers, LAMOST is capable of providing millions of stellar spectra in both medium-resolution mode ( ∼ 7500) and low-resolution mode ( ∼ 1800).As of data release 10 (DR10), LAMOST has released spectra for 22.29 million stars.Notably, LAMOST's observation strategy involves multiple-epoch observations of the same target, contributing significantly to the application of dynamical methods in the search for binary systems containing compact objects (Gu et al. 2019;Zheng et al. 2019;Yang et al. 2021;Mu et al. 2022;Li et al. 2022;Mazeh et al. 2022;Yuan et al. 2022;Qi et al. 2023).
In this work, we present two close binary systems discovered through radial velocity observations, each composed of a white dwarf candidate and an M-type companion.Section 2 describes the discovery of these systems and follow-up spectroscopic observations.In Section 3, we discuss the properties of the two sources, including the measurement of visible star's parameters and the determination of their orbital parameters.The discussion and summary are provided in Sections 4 and 5, respectively.

OBSERVATION AND DATA ANALYSIS
LAMOST J035540+381550 (hereafter J035540) and LAMOST J035916+400732 (hereafter J035916) have been identified as singlelined spectroscopic binaries each consisting of a K/M type companion and an unseen compact object candidate by Mu et al. (2022) based on the LAMOST low-resolution survey (LRS).J035540 and J035916 were selected due to their significant radial velocity variations.Since the photometric data from TESS of J035540 and J035916 do not show significant photometric periodicity, Mu et al. (2022) estimated a theoretical lower limit of the orbital period based on the Roche-lobe radius assumption and the mass-radius relation for low-mass main sequence star.

LAMOST spectra
The LAMOST LRS covers a wavelength range of 3690 to 9100 Å with a spectral resolution of approximately 1800.An LAMOST LRS observation typically involves three consecutive exposures, each lasting 10-30 minutes.Subsequently, the spectra from these exposures are merged to improve the signal-to-noise ratio (SNR).The singleexposure spectra are archived in the first data release of LAMOST LRS single-epoch spectra (Bai et al. 2021).These spectra consist of two arms: the blue arm covers 3690-6000 Å, and the red arm covers 5700 to 9100 Å.
The combined LAMOST spectra for J035540 and J035916 are displayed in Figure 1 (gray).For J035540, LAMOST LRS conducted six observations on December 30, 2015, and December 2, 2016, each with three consecutive 1800-second exposures.J035540 exhibited distinct radial velocity variations on both observation nights.On December 30, 2015, a rapid change of approximately 88 km s −1 was observed within one hour, while on December 2, 2016, a similar velocity change of about 84 km s −1 occurred within a similar timeframe.
In the case of J035916, LAMOST LRS made a total of nine observations.These observations were scheduled on December 30, 2015, December 2, 2016, and December 6, 2021, with three consecutive exposures.Each exposure lasted 1800 seconds.Notably, the most significant change in radial velocity occurred during the observation on December 2, 2016, transitioning rapidly from −127 km s −1 to −14.6 km s −1 within one hour.Detailed data are presented in Table 1.

The P200 spectra
We conducted follow-up observations on J035916 and J035540 using the Double Spectrograph (DBSP) mounted on the Palomar Observatory's 200-inch telescope (P200).The observations were taken during a dark night on September 11, 2023.The DBSP was configured with a wavelength range of 3900-5400 Å, achieving a resolving power R∼3400 for the blue arm, and a wavelength range of 6300-7900 Å, with R∼4500 for the red arm.Wavelength calibration for the blue and red arms involved the utilization of HeNeAr and FeAr lamps, respectively.
For J035916, we took a total of 6 spectra.Except for the initial spectrum with a 1200-second exposure, the remaining spectra were exposed for 1500 seconds each.As for J035540, we obtained a total of 3 spectra, each with an exposure time of 1200 seconds.Observation details are listed in Table 1.
We used the Python package Pypeit 1 (Prochaska et al. 2020) to reduce the spectra in the blue arm.Since Pypeit does not support the wavelength calibration of the red arm spectrum of the DBSP 1200/7100 D55, we turned to the IRAF to process the spectra in the red arm.The heliocentric correction was performed using the pyasl package PyAstronomy2 (Czesla et al. 2019).The resulting reduced spectra are depicted in Figure 1.

The Xinglong 2.16m Telescope Spectra
We applied for spectral observations of J035916 and J035540 using the 2.16 m telescope situated at the Xinglong Observatory on October 12, 2023.We chose the G7 grism paired with a 1.8 arcsec slit configuration.We exposed both targets to two consecutive observations, each lasting 1800 seconds.We used the IRAF to reduce the spectra using standard procedures.The heliocentric correction was performed using the PyAstronomy.
Due to the extreme faintness of both target sources and adverse weather conditions on the nights of observation, the obtained spectral quality was exceptionally poor.The spectra of J035540 exhibited an extremely low signal-to-noise ratio, rendering the measurement of radial velocity unfeasible.However, two spectra of J035540 did reveal a clear H emission line (see Figure 1).In contrast, the spectra of J035916 did not exhibit any distinct spectral structures due to the low quality, thus are not presented in Table 1.
Table 1.Statistics of the observed spectra of J035540 and J035916.Notes: ID: the target id; Telescope: the telescopes used for the observations; HMJD: the modified heliocentric Julian date at the middle of an exposure; Obs.Data: the UTC time for the middle of an exposure; Exptime: duration of each observation; RV(blue): the radial velocity measured from the blue side of the spectrum; SNR(blue): the signal-to-noise ratio of the blue side of the spectrum; RV(red): the radial velocity measured from the red side of the spectrum; SNR(red): the signal-to-noise ratio of the red side of the spectrum.The spectra from Xinglong 2.16m, due to poor quality, did not yield measurable radial velocities.* The difference between RV(red) and RV(blue) is attributed to atmospheric refraction causing greater displacement in the blue arm of the spectrum at lower zenith angles.During the exposure, the airmass for J035916 was 1.722, corresponding to a zenith angle of 35.5 • .According to Law et al. (2015, Figure 8), we estimate a deviation of ∼ 30 km s −1 in the blue arm RV measurement due to atmospheric refraction, which is consistent with the deviation between blue and red arms.

The Gaia information
We cross-match the two sources with Gaia data release 3 (Gaia DR3; Gaia Collaboration et al. 2023) using a matching radius of 2 ′′ .J035540's Gaia DR3 ID is 220660614918463232, and the phot g mean mag is recorded at 16.37 mag.Gaia DR3 provides a parallax of  = 4.64 ± 0.06 mas, corresponding to a distance of  = 215.5 ± 2.8 pc.For J035916, its Gaia DR3 ID is 229858064046205312.with phot g mean mag of 16.75 mag.Gaia DR3 provides a parallax measurement of  = 2.90 ± 0.06 mas, inverting to a distance of  = 344.8±7.0 pc.Astrometric measurements by Gaia are listed in Table 2.

Radial Velocity Measurements
To measure radial velocities, we employ the MARCS template, which is a set of theoretical stellar model atmospheres and flux sample files based on the Uppsala software by the same name (Gustafsson et al. 1975(Gustafsson et al. , 2008)).The optimal matching spectral template is selected by interpolating the MARCS model grid, minimizing the chi-square between the observation and the template.Figure 2 illustrates the observed spectrum and the best-matched template spectrum.Both the blue and red sides of LAMOST LRS and P200 spectra are used to measure the radial velocities, using wavelength ranges of 4500-5600 Å (blue) and 6000-8700 Å (red).Regions containing emission lines and strong telluric absorption lines are masked.We use the Cross-Correlation Function (CCF) method to measure the radial velocity by cross-matching the template spectrum with the observed spectrum.Simultaneously, we verify that two systems are single-lined spectroscopic binaries by confirming the single-peaked nature of the CCF profile (Merle et al. 2017;Li et al. 2021).The radial velocity uncertainties are estimated using the flux randomization/random subset sampling (FR/RSS) method, employing a Monte Carlo approach.This involves randomly generating a new flux point within the error range for each spectrum flux point to create a new spectrum with random removal of some data points.This Monte Carlo process is repeated 1000 times, and the standard deviation of the resulting mock radial velocity distribution is adopted as the velocity uncertainty.Table 1 summarizes the measured radial velocities and their uncertainties.

Orbital parameters
We apply the Lomb-Scargle algorithm (Lomb 1976;Scargle 1981) to search for orbital periods ( orb ) in the radial velocity data points of J035540 and J035916.However, the sparse phase coverage of the radial velocity data makes these results less reliable.To obtain 4000 4500 5000 5500 6000 6500 7000 7500 8000 more robust period information, we turn to their photometric data.We collect the Zwicky Transient Facility (ZTF; Bellm et al. 2019) DR19 data for both targets.For J035540, its ZTF g-band photometric dataset has 406 points during the observation from 2018 July 2 to 2023 February 16.The r-band has 851 points during the observation from 2018 July 13 to 2023 July 3. Similarly, J035916's ZTF gband and r-band photometric datasets consist of 386 and 810 data points, spanning from 2018 March 27 to 2023 April 2, and July 21, 2018, to July 3, 2023, respectively.The Lomb-Scargle method is then employed on the photometric data points to search for their respective photometric periods.
In the case of J035540, the corresponding peak period does not fit well with the radial velocity curve.Subsequently, we employ the period corresponding to the second-highest peak in the periodogram to fit both the radial velocity curve and the light curve.Remarkably, twice this period effectively fits the variations in both curves.This phenomenon is prevalent in ellipsoidal variables, where the light curve exhibits two peaks and valleys due to tidal deformation of the visible companion under gravitational forces.This results in a quasi-sinusoidal flux variation caused by changes in surface area.Consequently, we deduce the orbital period as 0.4768791 days.For J035916, fitting the highest power peak in the periodogram to both the light curve and radial velocity data yields an orbital period of 0.3904818 days.
We assume a circular orbit for the binary system and utilize the sinusoidal curve  R = − 1 sin[2( −  0 )/ orb ] +  to fit the radial velocities, where  1 is the radial velocity semi-amplitude, t is the mid-time of each exposure,  0 is the epoch of superior conjunction, and  is the systemic velocity.We use the radial velocities measured from the red arm of the spectrum for fitting, since the red arm spectra exhibits higher SNR and resolution than the blue one.For J035540, the fitting results indicate that  1 = 164.8± 2.1 km s −1 and  = −8.3± 1.0 km s −1 , while for J035916,  1 = 159.4± 1.7 km s −1 and  = −9.3± 1.2 km s −1 .The fitted radial velocity curves are shown in Figure 3 and the fitted orbital parameters are listed in Table 2.

Stellar Parameters
For J035540 and J035916, we conduct broadband spectral energy distribution (SED) fitting to constrain the stellar parameters.We use the python package astroARIADNE3 , which employs the Nested Sampling algorithm to automatically fit the SED of target stars (Vines & Jenkins 2022).
The SED fitting procedure involves the following steps.Firstly, astroARIADNE uses the python package astroquery to automatically collect multi-band photometric data, including GALEX (Mar-  (Hauschildt et al. 1999;Allard et al. 2012), and Kurucz (Kurucz 1993) to fit the SED and derive stellar parameters such as effective temperature, surface gravity, metallicity, distance, radius, and V-band extinction.The resulting best-fitting stellar parameters are summarized in Table 2, and the corresponding models are illustrated in Figure 4.
After the SED fitting, we estimate the masses of the visible stars in J035540 and J035916 using stellar evolution models.The SED best-fitting parameters and photometry are served as inputs for deriving isochrone masses with the MESA Isochrones & Stellar Tracks (MIST; Dotter 2016) isochrones.This procedure is implemented in the python package isochrones 4 .The isochrone masses are presented in Table 2.

Mass constraints
Using the the radial velocity semi-amplitude  1 for the visible stars and  orb , the mass function for the invisible companion can be calculated as: 4 https://isochrones.readthedocs.io/en/latest/ where  1 is the mass of the visible star and  2 is the mass of the invisible companion,  is the orbital inclination angle, and G is the gravitational constant.The mass function provides a robust lower limit on the mass of the unseen companion when  = 90 • and  1 = 0 ⊙ .The mass functions for J035540 and J035916 are 0.22 ± 0.01 ⊙ and 0.16 ± 0.01 ⊙ , respectively.We adopt the isochrone mass as the mass estimate for the visible stars.Utilizing the previously determined mass functions for the unseen companions, we calculate the lower mass limit for these invisible stars, assuming a binary orbital inclination of 90 degrees.The  min 2 of the unseen companion in J035540 is 0.64 ± 0.04  ⊙ , while for J035916, it is 0.57 ± 0.05  ⊙ .
The results indicate that the lower mass limit of the unseen companions in both J035540 and J035916 exceeds the mass of the visible stars.In typical binary systems with similar masses, the companion stars contribute significantly to the overall flux and exhibit distinct double-lined spectroscopic features, which are not observed in the obtained spectra.Therefore, it can be concluded that both systems harbor hidden compact stars.

Light curve and filling factor
Due to the relatively faint nature of J035540 and J035916, the photometric data from ZTF exhibit significant uncertainties.Nevertheless, distinctive variations in the photometric flux are visually apparent in the noisy phased light curves.The folded light curves reveal the characteristic signature of ellipsoidal modulated variability, as discussed in Section 3.3.The standard ellipsoidal modulation light curve features a double-peak-double-valley shape, a consequence of tidal distortion of the star by a companion.The star is elongated to an ellipsoidal shape.Furthermore, gravity-darkening induces different minima at the two distinct peaks when the longest axis of the ellipsoidal body is oriented towards the observer (El-Badry et al. 2021).However, this second-order effect is nearly imperceptible in the noisy light curve.
To enhance the clarity of the light curve profile, we use a threeharmonic model for fitting the normalized light curve.This model can be expressed as (Morris & Naftilan 1993): where  1 ,  2 , and  3 are the fitted parameters.The fitted model reveals that the light curve of J035540 exhibits a standard ellipsoidal modulation, while the light curve of J035916 displays an asymmetric double-peaked variation.The light curves and the fitting model for two targets are shown in Figure 5 After determining their binary orbital parameters and stellar parameters of the visible star, we calculate their Roche lobe filling factor.The Roche lobe filling factor is defined as  ≡  1 / L1 , where  L1 is the Roche lobe radius of the visible star.We calculate  L1 using the equation (4) of Qi et al. (2023): where  day is the orbital period in the unit of days.We use Equation (3) to estimate  L1 , while  1 is measured from SED fitting.Following this, the filling factors for J035540 and J035916 are calculated as 0.62 ± 0.06 and 0.65 ± 0.05, respectively, aligning with the observed significant variations in the light curve.

Constraining the invisible object's mass
Thus far, we have established the lower bounds for the masses of the compact objects in J035540 and J035916 using the mass function, assuming an orbital inclination of  = 90 • .However, to determine the actual masses, we need to establish the orbital inclination of the binary systems.J035540 displays a characteristic double-peaked ellipsoidal modulation, with key influencing factors being the binary inclination , mass ratio  =  2 / 1 , Roche lobe filling factor  , limb-darkening  During the light curve fitting, we set distortion method = none, treating the invisible companion as an object with no flux or eclipse contribution.Simultaneously, we set the orbital period and effective temperature of the visible star to the values listed in Table 2. Constraints on the projection of the orbital semi-major axis are applied using the orbital period and radial velocity semi-amplitude.The limb-darkening coefficient is manually set, following a logarithmic limb-darkening law, with values of [0.8, 0.3] (see Claret & Bloemen 2011).The gravity-darkening coefficient is treated as a free parameter, with a prior value of  1 ∼ N (0.5, 0.2) (see Claret & Bloemen 2011).We employ the emcee package for Markov Chain Monte Carlo (MCMC) fitting, using the provided values for radius, mass, and distance in Table 2 as priors.Multiple parallel chains are run, each consisting of 10,000 steps.
In Figure 6, the posterior samples for , , and  1 are presented.The fitting results indicate  = 71.30+13.55  −12.61 degrees for J035540, along with an associated mass for the unseen companion of  2 = 0.70 +0.12 −0.05  ⊙ .The optimal PHOEBE model outcomes are represented as the black lines in Figure 5, indicating that the ellipsoidal model can fit the observed light curve aptly.
The light curve of J035619 exhibits an asymmetric double-peaked structure, possibly attributed to spots on the visible star's surface (Luger et al. 2021;Rowan et al. 2024).Modeling this asymmetric light curve with PHOEBE is highly dependent on the spot model.In our analysis, we use the work proposed by Gomel et al. (2021) to fit the light curve of J035916 for providing certain constraints on the inclination.This approach is derived from the work of Morris & Naftilan 1993 (MN93) (Morris & Naftilan 1993) and offers an analytical approximation for ellipsoidal modulation under the assumption of a circular orbit.The leading corrected approximation  can be written as where  is the mass ratio  ≡  2 / 1 ,  is the average luminosity,  0 is the stellar brightness without secondary,  is the orbital phase and the  2 is a function of the linear limb-and gravity-darkening coefficients of the detected star, the  () is the Eggleton (1983) approximation for the volume-averaged Roche lobe radius and  (,  ) is the correction factor related to  and  .The correction factor begins at 1 for  = 0 (no correction) and rises monotonically to a maximum value of ∼ 1.5 at  ≳ 0.90 (Gomel et al. 2021).The modified formula demonstrates that the ellipsoidal amplitude depends on the filling factor, the inclination, and the mass ratio.The mass function is invoked to place an additional constraint between the inclination angle and mass ratio for a given visible star's mass.Based on the known stellar properties and orbital parameters, the observed ellipsoidal variability can provide a combined constraint on the mass of the invisible component of the binary system.
Figure 7 illustrates the constraints on the orbital inclination and the mass of the invisible star using the aforementioned analysis model.The graph illustrates the correlation between the half-amplitude of ellipsoidal variations and the orbital inclination.The results reveal that the orbital inclination of J035916 falls within 50 • − 90 • , corresponding to an inferred mass of the compact object ranging from 0.57 − 0.90  ⊙ .
We use PHOEBE to confine the inclination of J035540 and subsequently utilize the analysis model to restrict the inclination of J035916.However, due to the noisy nature of the photometric data used, our constraints on the true orbital inclination remain relatively broad.To further tighten the constraints on the inclination, one can enhance measurements by measuring the rotational broadening of the visible stars from high-resolution spectra (Marsh et al. 1994).

The 𝑣 rot sin 𝑖 measurement
As discussed in Section 4.2, the ZTF light curve of J035916 exhibits asymmetrical features, preventing us from accurately fitting the orbital inclination using an ellipsoidal model.An alternative approach is to independently constrain the inclination by measuring the projected rotational broadening velocity ( rot sin ).We selected the P200 red-side spectrum at HMJD 60198.416755 to measure the  rot sin .This spectrum was taken at the quadrature orbital phase ( ∼ 0.25), therefore the radial velocity smearing effect is minimized.We estimated that the velocity smearing is only around 6.4 km s −1 due to the star's orbital motion during the 1500s on-target exposure.
Major stellar spectral-line broadening mechanisms include rotational broadening, macro-turbulence broadening, and instrumental broadening.We determined the instrumental broadening by measuring the FWHM (full width at half maximum) of the lamp spectrum taken during the P200 observation.Specifically, a Gaussian kernel was used to fit the lamp spectrum line profile, yielding the FWHM = 1.4 Å at 6402 Å, which corresponds to a spectral resolution of R ∼ 4400, or equivalently, a broadening of Δ ins ∼ 68 km s −1 in the velocity space.The theoretical rotational broadening kernel (Gray 2005) is as follows: where  = / rot sin , and  is the linear limb-darkening coefficient.
We set the macro-turbulence broadening to be 2 km s −1 .We select the spectral range of 6400-6500 Å for fitting, where a few relatively strong and distinct absorption lines are present.Additionally, it is free of telluric absorption lines.We generate the spectral template from MARCS model atmospheres using the stellar parameters listed in Table 2.The template is first convolved with the rotational broadening kernel (Equation 5), then convolved with Gaussian instrumental broadening kernel.We obtain the  rot sin  by minimizing the chi-square  2 = (( − )/) 2 between the model and observed spectra, where  and  are the fluxes of the observation and model, respectively;  is the flux uncertainty, and the summation sign is taken over all the flux points with the spectral range.
The fitting result is  rot sin  = 59 +6 −5 km s −1 .Figure 8 presents the optimal fitting outcomes.The short orbital period of J035916 implies that the system is tidally synchronized.Consequently, we adopt the assumption that the rotation period ( rot ) of the visible star is equivalent to  orb .Therefore, we can derive the rotational velocity as  rot = 2 vis / orb = 70±3 km s −1 , where  vis is the radius of the visible star.Thus, we estimate that J035916's inclination  = 57 +18 −13 degrees, consistent with that estimated from the analytical ellipsoidal model (Section 4.2).
In Figure 9, we depict the locations of two targets in the Toomre diagram, both falling within the region associated with the thin disk.Following this, we compute the thick-to-thin disk probability ratio (TD/D) using the equation ( 7) from Li & Zhao (2017).The TD/D values are 0.08 for J035540 and 0.02 for J035916, indicating that they are most likely thin disk stars.

The nature of two targets
As described above, we constrain the orbital inclinations of J035540 and J035916 using the light curve, thereby restricting the mass of the unseen stars in these systems.The inferred mass of the unseen stars suggests that these two systems are binary systems composed of a WD candidate and an M-type star.The most reliable method to confirm the presence of WDs and to measure their parameters is through their UV spectra (Parsons et al. 2016;Hernandez et al. 2022b).Since a WD has radiation peaked at near UV to UV wavelengths, we tried crossmatching J035540 and J035916 with GALEX using a cross-matching radius of 20 ′′ , to search for potential UV counterparts.However, J035540 and J035916 are not covered in GALEX's footprints.
The future Chinese Space Station Telescope (CSST) will conduct a wide-field imaging in the near-ultraviolet-optical-nearinfrared (NUV-OP-NIR) bands (Zhan 2021).In the context of white dwarfs, we utilize WD cooling models (Bédard et al. 2020) to assess their detectability by CSST based on assumed temperatures.The tables from http://www.astro.umontreal.ca/~bergeron/CoolingModels provide bolometric and absolute magnitudes for WDs on various photometric systems including NUV colors.We use the cooling model for a typical DA-type WD with mass of  2 = 0.60 ⊙ to calculate the absolute NUV magnitudes across a grid of WD temperatures.The visual NUV magnitude from CSST is  NUV =  NUV + 5log() +  NUV − 5, where  NUV is the absolute NUV magnitude obtained from the cooling model,  is the distance from Gaia DR3, and  NUV is the NUV extinction value.The  NUV = 0.25 is derived from the Milky Way Average Extinction Curve (Gordon et al. 2009).CSST has a detection limiting magnitude for NUV at approximately 25.4 mag.We conclude that the compact objects in J035540 and J035916, if they are DA-type WDs, will be detected by CSST when their temperatures exceed 6000K and 6500K, respectively.
We use the dynamical method for the identification of two binary systems that potentially harbor unseen WDs.This method distinguishes itself from both the spectroscopic decomposition method (Li et al. 2014;Rebassa-Mansergas et al. 2016) and the UV-excess selection method (Parsons et al. 2016).More than ten sources identified through these methods have also undergone the dynamical measurement to further investigate their binary properties (Rebassa-Mansergas et al. 2016;Hernandez et al. 2021Hernandez et al. , 2022a;;Parsons et al. 2023).The dynamical approach is particularly advantageous in this context due to the inherent faintness of WDs compared to brighter main sequence stars.
Spectroscopic decomposition methods are typically limited to detecting WD-M dwarf binaries, as these systems exhibit a prominent white dwarf signature in their spectra.Notably, the temperature of WDs identified through spectroscopic decomposition tends to be sufficiently high to show a discernible white dwarf component in the spectrum.On the other hand, the UV-excess method faces challenges in identifying sources lacking UV data points, as observed in objects like J035540 and J035916.Additionally, the UV radiation from stellar chromospheres can significantly contaminate the UV flux of candidates (Linsky 2017).In contrast, the dynamical method proves to be a more robust approach.By utilizing radial velocity curves to derive orbital parameters and employing mass functions, this method allows us to confidently determine the mass of unseen objects.

SUMMARY
J035540 and J035916 are two binary systems containing compact object candidates identified through meticulous screening of LAMOST LRS, with the visible components being M-type stars.Follow-up spectroscopic observations from P200 augment the radial velocity measurements of two sources.We obtain their orbital periods and radial velocity semi-amplitudes from LAMOST LRS and P200 radial velocities, in conjunction with the ZTF light curves.The mass functions of the unseen companions are calculated, with  ( 2 ) being 0.22 ± 0.01 ⊙ for J035540 and 0.16 ± 0.01 ⊙ for J035916.Mass estimates for the visible M-type stars are derived through SED fitting, which in turn constrains the minimum mass of the invisible companions via mass functions.Given that the invisible components' minimum masses exceed those of the visible stars and that the spectra are consistent with single-lined spectroscopic binaries, we infer that these unseen stars are likely compact objects.J035540's ZTF light curve is modeled using PHOEBE with a pure ellipsoidal modulation, providing constraints on the system's inclination.Consequently, we obtain a mass of  2 = 0.70 +0.12 −0.05  ⊙ of the compact object.Due to the asymmetric double-peaked nature of J035916's light curve, which exhibits strong model dependence when simulated with a spot model, an analysis based on ellipsoidal modulation is employed.This analysis constrains the light curve amplitude, thereby restricting the orbital inclination and yielding a mass range of 0.57 − 0.90 ⊙ for the compact object.We measure the  rot sin  of J035916 as 59 +6 −5 km s −1 from the red-side P200 spectrum at the quadrature orbital phase, thereby constraining the orbital inclination as 57 +18 −13 degrees.This is consistent with the result from the analytical ellipsoidal model.
The mass estimates suggest that both compact objects are likely WDs.These results highlight the efficacy of optical time-domain survey data and dynamical methods in the search for WDs.This approach provides a promising avenue for detecting WDs, even for those extremely faint WDs.Moreover, this methodology holds the potential for uncovering other compact objects within binaries, including stellar-mass black holes and neutron stars, even during their quiescent states.

DATA AVAILABILITY
This paper uses the data from the LAMOST low-resolution survey.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.ZTF is a public-private partnership, with equal support from the ZTF Partnership and from the U.S. National Science Foundation through the Mid-Scale Innovations Program (MSIP).

WavelengthFigure 1 .
Figure 1.Normalized spectra for J035540 and J035916.Observational information is annotated next to each spectrum.LAMOST LRS spectra observed on the same night have been combined while P200 and Xinglong 216 spectra are single exposures.

Figure 2 .
Figure2.The orange curves represent the best-matching model spectra used for radial velocity measurements.The grey curves represent the P200 spectrum for the blue and red arms.The grey shading indicates the wavelength of the observed spectrum contaminated by emission or telluric lines, which is masked during the fitting process.

)Figure 3 .
Figure3.The fitted radial velocity curves of J035540 and J035916.The dots represent the radial velocities from LAMOST and P200, while dashed lines depict the best-fitting radial velocity curves.

Figure 4 .
Figure 4. SED fitting results for J035540 (top panel) and J035916 (bottom panel).Cyan points denote observed fluxes, purple diamonds indicate synthetic fluxes, and the grey curve represents the best-fitting SED model.

Figure 5 .Figure 6 .
Figure 5. ZTF light curves (grey dots) of J035540 and J035916 folded by their respective orbital periods.The green curves depict the three-harmonic model based on Equation 2. In the top panel, the black dotted line represents the best-fitting PHOEBE model for J035540.

Figure 7 .
Figure 7.The panel depicts the semi-amplitude of ellipsoidal variability versus orbital inclination.The horizontal shaded region represents the ellipsoidal modulation amplitudes observed in J035916's light curve, while the vertical shaded region indicates the corresponding range of orbital inclinations.

Figure 8 .
Figure8.Spectral wavelength region used to measure the rotational broadening of J035916.The grey curve represents the observed spectrum and the orange curve represents the rotationally broadened best-fit template.

Figure 9 .
Figure 9. Toomre diagram of our sources.Dashed lines indicate constant peculiar space velocities in steps of 50 km s −1 .Blue dots represent the thick disk stars and red dots represents the thin disk stars adopted from Bensby et al. (2003).

Table 2 .
Summary information for J035540 and J035916, including astronomical parameters from Gaia DR3, orbital parameters derived from ZTF light curves and radial velocities measured from LAMOST LRS and P200 spectroscopic data, and stellar parameters determined through SED fitting.