Discovery of a Proto–White Dwarf with a Massive Unseen Companion

We report the discovery of SDSS J022932.28+713002.7, a nascent extremely low-mass (ELM) white dwarf (WD) orbiting a massive (>1 M ⊙ at 2σ confidence) companion with a period of 36 hr. We use a combination of spectroscopy, including data from the ongoing fifth-generation Sloan Digital Sky Survey (SDSS-V), and photometry to measure the stellar parameters of the primary pre-ELM WD. The lightcurve of the primary WD exhibits ellipsoidal variation, which we combine with radial velocity data and PHOEBE binary simulations to estimate the mass of the invisible companion. We find that the primary WD has mass M 1 = 0.18−0.02+0.02 M ⊙ and the unseen secondary has mass M 2 = 1.19−0.14+0.21 M ⊙. The mass of the companion suggests that it is most likely a near-Chandrasekhar-mass WD or a neutron star. It is likely that the system recently went through a Roche lobe overflow from the visible primary onto the invisible secondary. The dynamical configuration of the binary is consistent with the theoretical evolutionary tracks for such objects, and the primary is currently in its contraction phase. The measured orbital period puts this system on a stable evolutionary path which, within a few gigayears, will lead to a contracted ELM WD orbiting a massive compact companion.


INTRODUCTION
Corresponding author: Gautham Adamane Pallathadka gadaman1@jh.eduWhite dwarfs (WDs) are remnants left behind by main sequence stars with masses in the range 0.8-8 M ⊙ .While the typical masses of white dwarfs range from 0.5 M ⊙ up to 1.4 M ⊙ , extremely low mass WDs (ELM WDs) with masses in 0.15 − 0.35 M ⊙ have been detected in binaries (Iben & Tutukov 1985;Brown et al. 2010;Istrate et al. 2016;El-Badry et al. 2021).These ELM WDs are usually formed when a massive main sequence star loses its outer envelope to a companion.The helium core that remains with a hydrogen atmosphere is an ELM WD (Nelson et al. 2004;Li et al. 2019).They have been theorized to form only in binary systems -their mass is too low to be formed from single star evolution within a Hubble timeand this scenario is supported by observations as well (Marsh et al. 1995;Brown et al. 2010).
ELM WDs can be found in binaries with WDs (Brown et al. 2020), pulsars (Nelson et al. 2004;Antoniadis et al. 2013;Istrate et al. 2014a) or main sequence stars (Maxted et al. 2011).These binaries are interesting because they help us piece together the evolutionary pathways and test our understanding of binary evolution.Compact binaries of double WDs are particularly interesting since they are possible Type Ia supernova progenitors (Nomoto 1984;Maoz et al. 2014;Soker 2019;Jha et al. 2019).They are also dominant sources of lowfrequency gravitational waves for future observatories like the Laser Interferometer Space Antenna (LISA; Marsh 2011;Korol et al. 2022).WD mergers with neutron stars (NS) or black holes are candidate progenitors of gamma ray bursts or other types of explosive transients (Fryer et al. 1999;Margalit & Metzger 2016;Bobrick et al. 2017;Zenati et al. 2019Zenati et al. , 2020;;Bobrick et al. 2022;Kaltenborn et al. 2022).In binaries with pulsars, the mass transfer phase can lead to spin-up of the pulsar giving birth to milli-second pulsars.These systems can be used to independently verify the spin-down age of milli-second pulsars (Nelson et al. 2004).
While ELM WDs have surface gravity values in the range log g = 5 − 6.5 (where g is in cgs units), in the early stages of their evolution they can appear as bloated pre-ELM WDs with a smaller log g (Lagos et al. 2020;Wang et al. 2020;Kupfer et al. 2020a,b;El-Badry et al. 2021).Using spectroscopic and photometric data one can solve for the parameters of the visible object, hereafter called the primary.Without additional information it is not possible to solve for the mass of the invisible -secondary -object because of a degeneracy between the mass of the secondary and the system's inclination.Because of their bloated nature, the surfaces of pre-ELM WDs can be tidally distorted in the presence of a companion.This change in the shape of stellar surface leads to a periodic variation in the lightcurve.Therefore, the mass-inclination degeneracy can be further constrained with light curve data.
In this paper we report the discovery of SDSS J022932.29+713002.6 (hereafter SDSS J0229+7130), a binary consisting of a bloated pre-ELM WD and an invisible companion with a mass close to and perhaps exceeding the Chandrasekhar mass.We summarize the observations in Sec. 2. We analyze them and present the measurements of the binary parameters in Sec. 3. We discuss the results in Sec. 4. Throughout the paper g is the acceleration due to gravity in cgs units (cm s −2 ), [Fe/H] is the abundance ratio relative to the Sun on a logarithmic scale, Z is the metallicity of the star, i.e. the mass fraction of metals on a linear scale, and T eff is the stellar temperature in Kelvin.Inclination angles are measured relative to the plane of the sky, so that i = 90 deg is an edge-on system.We refer to the less massive WD as being the primary component, and the secondary is the more massive invisible companion.

Spectroscopic Data
We identified SDSS J0229+7130 in the firstyear data from the Milky Way Mapper, a multi-epoch Galactic spectroscopic program in the fifth-generation Sloan Digital Sky Survey (SDSS-V, Kollmeier et al. 2017;Kollmeier et al. 2019).Milky Way Mapper started in November 2020 using the Apache Point Observatory 2.5 m telescope (Gunn et al. 2006) and is now also operating at the Las Campanas Observatory 2.5 m telescope (Bowen & Vaughan 1973) using both the Baryon Oscillation Spectroscopic Survey spectrograph (BOSS; Smee et al. 2013) and the Apache Point Observatory Galactic Evolution Experiment (APOGEE) spectrographs (Wilson et al. 2019).Our target was observed for a total of 9 individual exposures between 2020-12-04 to 2020-12-06 with BOSS.Each exposure lasted 900s and covered a wavelength range from 3600 Å to 10000 Å at a resolution of R ≃ 1800.The wavelengths of the spectra are corrected to the heliocentric frame, and the absolute wavelength calibration of each exposure is accurate to <10 km s −1 .The BOSS data products in this paper were derived by using IDLspec2D v6_1_0.
We initially flagged SDSS J0229+7130 due to its significant radial velocity (RV) variation ≳ 100 km s −1 between different nights, albeit with minimal variation across successive exposures on a given night.This suggested that both the RV amplitude and orbital period are large, hallmarks of a high mass function binary.Spectroscopically, SDSS J0229+7130 looks like a pure hydrogen atmospheric DA White Dwarf, with narrow Balmer lines (Fig. 1).In addition, there is a faint sodium absorption doublet Na I D near 6000 Å and calcium absorption line Ca II K near 3933 Å, which are discussed further in Section 3.1.
Although the SDSS-V data suggested a high mass-function binary, the RV data were too sparse to better constrain the mass function and definitively measure the orbital period.To fill in the orbital phase coverage, we obtained follow-up spectroscopy with the Dual-Imaging Spectrograph (DIS) on the 3.5 m telescope at the Apache Point Observatory (APO).We used B1200/R1200 gratings with a 1.5 ′′ slit, deliver-ing resolution R ≈ 3000.We obtained twenty 900s exposures between three different nights on 2021-10-05, 2021-11-01 and 2021-11-05.The APO data were reduced by the standard IRAF pipeline.
The absolute wavelength calibration of the APO spectra showed systematic shifts, which increased with the wavelength.To account for this, we estimate the shift at 6498 Å (close to the Hα line) using the known positions of the sky lines, and measure the RVs incorporating this derived correction.We also increase the RV error by adding this correction in quadrature to the measured RV error, to get the most conservative RV values.The RVs for APO data are corrected to barycentric frame using the package barycorrpy 1 (Corrales 2015;Kanodia & Wright 2018).The low wavelength end of APO spectrum had issues and is unreliable.We only use the data with wavelength greater than 5900 Å.

Photometric Data
Our target has a Gaia parallax-inferred distance of d = 1.663±0.12kpc (Bailer-Jones et al. 2021;Gaia Collaboration et al. 2023).Looking at the Gaia color-magnitude diagram (Fig. 2), we see that the object lies well above the white dwarf track, and is only slightly below the main sequence.Thus we propose that this object is not a typical WD but rather a bloated pre-ELM WD.There are no resolved wide companions (at separation >2000 AU) near our target in Gaia, and there are no astrometric anomalies reported in Gaia DR3 for this source.So there is no evidence for additional companions to the spatially unresolved system.The object is close to the plane of the Galaxy with b ≈ 10 • , and given the inferred distance we need to carefully take into account the interstellar extinction, as explained in detail in Section 3. In addition to the spectroscopic data, this object has archival Zwicky Transient Facility (ZTF) observations (Masci et al. 2019), and we selected all observations released in DR18.The object was observed with the g filter between 2018-06-29 and 2023-02-20, and between 2018-06-14 and 2023-03-10 with the r filter.To select clean data, we choose all point sources in ZTF within 5" of target position obtained from Gaia with ZTF quality flag catflags = 0 and |sharp| < 0.25.We find that all the available data points correspond to our object and are within 5" of our target.However, quality cuts reduce the number of data points from 912 to 821, and we obtain a clean ZTF lightcurve.
This object also has archival Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) observations in the Full Frame Image (FFI) data in sectors 18 and 19.The data in sector 18 have instrumental artifacts and hence we restrict ourselves to data from sector 19.We obtain the lightcurve from FFI using eleanor2 (Feinstein et al. 2019;Brasseur et al. 2019).
We collect archival photometric data for this object using VizieR3 (Ochsenbein et al. 2000).This object has photometric data from the Sloan Digital Sky Survey in ugriz filters (Ahumada et al. 2020), from 2MASS in J, H, K s fil- ters (Skrutskie et al. 2006), from GALEX in FUV and NUV bands (Bianchi et al. 2017), from Pan-STARRS in grizy filters (Flewelling 2018) and from CatWISE in W1, W2 bands (Wright et al. 2010;Eisenhardt et al. 2020) covering the UV to infrared.These data are summarized in Table 1.We add a base error of 0.03 mags in quadrature to the existing errors to take into account systematics.

ANALYSIS
In this section we describe the steps followed in deriving the stellar parameters of the primary using spectroscopic and photometric data.Following this, we combine the spectroscopic analysis with the analysis of the ellipsoidal variation from the photometric data to obtain the mass of the secondary.The SDSS spectra for the target show periodic shifts in the position of absorption lines suggesting a varying radial velocity (RV).We choose to measure RV using Hα and Hβ simultaneously from SDSS-V spectra.We use just the Hα for APO spectra because the spectrum at low wavelength end is unreliable.We use corv 4 to measure the RVs for our SDSS and APO sub-exposures.corv computes the crosscorrelation between the observed stellar spectra and the template.For the template, we model each absorption line using two Voigt profiles with a common centroid.Once we obtain the RVs, we correct each spectrum for the Doppler shift associated with the measured RV and then coadd them in the target's rest-frame.One of the SDSS sub-exposures has flux calibration issues and we do not include it in coadding, but we still use it for the RV measurements since the relative flux calibration does not affect the position of the absorption lines.The coadded SDSS spectrum is shown in Fig. 1.

Stellar Parameters
Interestingly, the Na I D doublet near 5900 Å shows little sign of periodic variation, suggesting an interstellar origin rather than an association with the stellar photosphere.The measured temperature (> 8000 K, Table .2) is also far too high to observe photospheric sodium lines (Yamaguchi et al. 2023).To verify this, we coadd RV>0 sub-exposures and RV<0 subexposures separately and compare the two.We confirm that the sodium absorption is stationary and is therefore associated with the interstellar medium.We then coadd all the spectra without shifting them to the rest frame of the primary and fit a Gaussian profile to the the Na I D lines.We obtain equivalent widths for D1, centered at 5897.5 Å, and D2, centered at 5891.6 Å to be 0.79 Å and 0.82 Å respectively, and 1.62 Å for the combined D1+D2.Using the best-fit empirical relation between the Na I D equivalent widths and extinction from Poznanski et al. (2012), we derive the 2-σ lower limit on the color excess, E(B − V ), to be 0.57, 0.31, 0.68 using D1, D2 and D1+D2 respectively.The trace photospheric sodium may contribute to the absorption, and we have not incorporated the inherent error in the empirical relation, so these values are only used as indicators of high extinction, which needs to be carefully taken into account in the subsequent modeling of the system.
The Ca II K line at 3933.6 Å is from the stellar surface.This is a feature observed in many ELM WDs as seen in Gianninas et al. (2014).We find that this absorption feature shows clear RV variation.
The APO spectra of the Hα absorption line, as well as some higher order Balmer lines in SDSS sub-exposures, show weak shape features that could be interpreted as Zeeman splitting.Identifying the splitting in the coadded APO spectrum with the linear Zeeman effect (Preston 1970;Garstang 1977), we estimate that the requisite magnetic field to produce the observed splitting is around 0.2 MG.However, the Paschen series in the SDSS sub-exposures and the coadded SDSS spectrum show little sign of Zeeman splitting, and we proceed with the analysis assuming zero magnetic field.
Using 1D hydrogen-atmosphere WD model spectra5 (Koester 2010) to fit the spectrum gives a best-fitting model at the low gravity edge of the model grid, log g = 5.Therefore, to measure the stellar parameters, instead of using WD models we fit the spectrum with theoretical stellar models derived for main sequence stars with lower surface gravity, as one would expect for pre-ELM objects, and with varying metallicity.To generate theoretical spectra we use the BOSZ grid (Bohlin et al. 2017).We use the grid with grid points between temperatures of 6000 K and 12000 K, log g between 2 and 5, [Fe/H] between -2.25 and 0.5, each with spacing of 250 K, 0.5 dex and 0.25 dex respectively.
We start with theoretical spectra with resolution of R ∼ 10000, which are downgraded versions of theoretical spectra with R ∼ 30000 provided by Bohlin et al. (2017), much greater than the SDSS spectra.We convolve the theoretical spectra with a Gaussian to match the SDSS resolution of R ∼ 2000.Finally, we use a scipy function RegularGridInterpolator6 (Weiser & Zarantonello 1988;Virtanen et al. 2020) to interpolate between the stellar parameter grid points.We then fit the coadded SDSS spectrum with the interpolated and convolved theoretical spectra to obtain the best-fit log g, temperature and metallicity of the star.
For the actual fitting, we minimize the χ 2 between the rest-frame coadded SDSS spectrum and the theoretical spectra.The SDSS spectrum can have flux calibration issues and the extinction is not well determined.To take into account these long-wavelength features we multiply the theoretical spectrum by a polynomial function of wavelength.We examine the solution to make sure that the polynomial only corrects the long-wavelength features and does not over-fit the absorption lines.The results are presented for polynomial of order 6, the lowest order which gives a stable solution.The spectroscopic log-likelihood to maximize is given by where f theo (λ i ) is the theoretical spectrum evaluated at different wavelengths, while f obs,i and σ obs,i are the spectral flux densities and the associated errors at those wavelengths.We fit the extinction-corrected photometric magnitudes to obtain the radius, temperature and log g of the primary.We estimate the extinction using the Green et al. (2019) model, which is a 3D dust map derived from PanSTARRS.The extinction calculation is done using dustmaps7 (Green 2018) and pyextinction8 packages.We infer the (B −V ) color excess using E(B −V ) = 0.884×ϵ, where α is the dust reddening returned by the dustmaps package and represents the dust density in the line of sight (Green et al. 2018).We then use a Fitzpatrick (1999) extinction curve with R V = 3.1 to obtain the extinction in various filters.The Green et al. (2019) dust map is probabilistic and there is a spread associated with the predicted extinction.We use the median of all returned samples of α by dustmaps to calculate the color excess.Additionally, the conversion from the α values from Green et al. (2019) to E(B−V ) or the value of R V has uncertainty due to the varying dust size distribution and the resulting varying extinction curves.To take all of this into account we assign a conservative error of 20% to the median value provided by dustmaps.We use a Gaussian prior for E(B − V ) with mean and sigma calculated as described above and leave it as nuisance parameter in our fitting.There are other 3D extinction maps such as Capitanio et al. (2017); Lallement et al. ( 2019) which give somewhat different values for extinction.In particular, Lallement et al. (2019) predicts a much larger extinction and the disagreements between different extinction maps have been reported by them.We use the predictions from Green et al. (2019) for our fiducial solution and discuss the consequences of higher extinction in 4.
The theoretical spectra used for the photometric fitting are generated using the BaSeL library (Lejeune et al. 1997(Lejeune et al. , 1998)), which are photometrically corrected semi-empirical models.The theoretical models estimate the flux on stellar surface.We can multiply this by a factor of (R/d) 2 to obtain the apparent flux, where R is the radius of the primary and d is the distance to the object.The computation is performed using the package pystelllibs9 .We use the Gaia parallax and include the distance as a nuisance parameter using the prior provided by Bailer-Jones et al. (2021).We then use the pyphot 10 package to integrate the theoretical spectra through the appropriate filters, and fit the extinction-corrected magnitudes with the theoretical magnitudes to estimate the best-fit parameters.The actual fitting is again performed by likelihood maximization.We find that the photometric fits very weakly depend on metallicity, Z, and thus we can fix Z = 0.014 without affecting our results significantly.The likelihood, L phot , is defined similarly to the spectroscopic likelihood, replacing fluxes with appropriate magnitudes.
Finally, we perform joint spectrophotometric fitting by defining the combined likelihood as L priors is the likelihood associated with the priors for extinction and parallax as described earlier, and also the uniform priors for rest of the parameters with limits set by the BOSZ grid.For the fitting, we calculate loglikelihoods for each of these separately and add them.The posterior distribution is explored using emcee 11 (Foreman- Mackey et al. 2013).
The results showed the existence of two local minima centered around color excess of 0.41 and 0.51 respectively, and a degeneracy between E(B − V ) and other stellar parameters.In the combined fit described in Sec.3.3, we explore solutions around both these minima separately by restricting the color excess around these minima and setting a nominative error of 5%.

Lightcurve and Orbital Parameters
We can derive the orbital parameters by analyzing the radial velocity and the lightcurve data.We first use corv to measure RVs from both APO and the SDSS data, as decribed in Sec.3.1.Then we use the Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) from astropy on these RVs to obtain the orbital period, but find it unreliable because the sparse phase coverage results in a broad and unconvincing power spectrum.
11 https://emcee.readthedocs.io/en/stable/W can also use the lightcurve data to measure the period.Stars in close binary systems show periodic changes in flux due to eclipses, heating of the companion due to the primary (reflection binary ;Vaz 1985;Schaffenroth et al. 2022), or distortions of their surfaces induced by the tidal forces of the companion (ellipsoidal modulation; Kopal 1959;Green et al. 2023).We combine the data from both the g and r filters of ZTF and normalize the data to capture only the fractional variations about the median flux.We then use the Lomb-Scargle periodogram to derive the period of the orbit, and phase fold both the lightcurve as well as the RV data to this period.The periodogram is shown in Fig. 3. Amplitude of flux variation can differ depending on the filters used.However, for the purpose of period determination, we choose to collate data from both filters and treat them as single dataset.This should result in a more robust period determination while having negligible downsides.
Depending on the physical reasons that cause the flux variability, the orbital period of the system can correspond to the peak of the periodogram or be twice as long.The dominant effect of the tidal distortion is to produce a prolate ellipsoid pointed at the companion and this produces a quasi-sinusoidal variation in the emitted flux due to varying surface area.As discussed in El-Badry et al. ( 2021), the quasi-sinusoidal variation would have different minima at two different end-on configurations (when the longest axis of the ellipsoid is directed towards the observer) due to different gravity darkening.However, this is a second order effect which could be virtually invisible in a noisy lightcurve.If the minima are not distinguishable, then the periodogram peaks at half the orbital period.We find that phase folding to the period P = 35.8703hours, twice the period corresponding to the maximum power of the periodogram, fits both the lightcurve and the RV data well, and the resulting fit is shown in Fig. 4.
The variation of both the datasets behave as we would theoretically expect for ellipsoidal variation -with one of the maxima of lightcurve aligning with minima of RV curve and other maximum coinciding with the RV maxima.This strongly supports our argument that the origin of the photometric variation is orbital in nature and is due to the ellipsoidal variation.We do not observe any signs of eclipses, which would be signaled by a difference in the depth of the minima.When we fold the data to the period corresponding to the maximum power in the periodogram, we find that the RV data are not well fit.This rules out the possibility of a reflection binary, where the photometric variability is due to the differences in the temperature of the primary and the companion (Wilson 1990).For systems where only one of the binary stars is visible, reflection effect is seen for very hot primaries (> 30000K) (Hilditch et al. 1996).It is reasonable that this should not be the dominant effect in our source.
To verify the obtained period, we repeat the process with TESS data.There are two objects at 12" and 19" from our target, while each TESS pixel covers a region of 21".Hence, the observed flux from our source is blended with that of the contaminants and the absolute flux variability measured by TESS does not accurately represent the variability of our target.Therefore, we analyze the TESS lightcurve separately from the ZTF lightcurve to account for the potential blending with nearby sources.Indeed, we find that the fractional variability of the TESS lightcurve is about 0.4%, much smaller than that of the ZTF lightcurve (2%).Thus, we restrict ourselves to using the ZTF data for further analysis while using TESS to verify the period.We obtain periods of 35.8703 ± 0.0006 hours and 35.56 ±0.07 hours from ZTF and TESS data respectively.
From the Lomb-Scargle periodogram, we compute the un-normalized power spectral density -which is proportional to χ 2 defined between the lightcurve data and the periodogram model (Scargle 1982;VanderPlas 2018).To estimate the errors, we calculate the frequencies where this χ 2 increases by magnitude of one from the minimum.While there is a mis-match in the derived periods, at 5-σ level they are not inconsistent.Using TESS we also verify that the period of 36 hours is robust and not associated with aliasing because of the different cadences associated with each dataset.Once we obtain the period, to measure the amplitude of variations and to verify that the periodicity is indeed due to ellipsoidal variation, we re-fit the lightcurve and RV data to sinusoid (i.e assuming eccentricity e = 0) -first individually and then both simultaneously.We fix the period and use the phase-folded data while leaving semi-amplitude of flux variation, RV semiamplitude (K), systemic velocity (v γ ) and phase Figure 4. Radial velocity and ZTF lightcurve phase-folded at the best period.In the top plot, solid data points are the binned lightcurve while the raw data is plotted in the background.We plot the PHOEBE simulated RV and lightcurve data for the best-fit parameters.In the top two panels, for comparison, we also plot the best-fit sinusoidal models as black dashed lines and also plot 100 random samples of sinusoidal models, in orange, as of the errors involved.
(ϕ) as free parameters.We find that both variations behave exactly as we expect and are exactly in phase.We also find that the flux of our target varies by about (2.02 ± 0.14)% from the minimum to the maximum of the light curve and the RV semi-amplitude to be K = (169±4) km s −1 .Once we know the RV semi-amplitude, orbital period and mass of one of the stars we can use the binary mass function to express the inclination in terms of the mass of the second star: (2) Fig. 5 illustrates the parameter space that we need to explore.We see that based on the relatively long period combined with a relatively large RV variation, the secondary must be quite massive.Normally, one of the constraints on masses would be that the star should not overflow its Roche lobe radius (Eggleton 1983).However, in our case, the period of the binary is quite large and hence the orbital separation for reasonable masses (primary mass M 1 > 0.1 M ⊙ and secondary mass M 2 > 1 M ⊙ ) is large enough that the Roche lobe radius is much larger than the radius of the star and this constraint is never important.
Fig. 5 gives a lower limit on the secondary mass M 2 of about 0.8 M ⊙ .Given that we only see the primary, we suspect that the companion is a compact object.For any star, flux near the Rayleigh-Jeans tail of the stellar spectrum is proportional to R 2 T .Thus, the ratio between the flux of two stars at large wavelengths can be written as (R 1 /R 2 ) 2 (T 1 /T 2 ).We can use this to our advantage and rule out main sequence companions.If the companion is a 0.8 M ⊙ main sequence star -which is the lowest possible companion mass allowed by the radial velocity measurements -the radius and temperature is approximately equal to 0.7 R ⊙ and 0.84 T ⊙ , respectively.Using the derived mass and temperature of the primary, we get the flux ratio of the main sequence companion to the primary to be 0.83.So the main sequence companion contributes almost as much as the primary at large wavelengths and the total flux of the system would be twice the expected flux from the primary.We do not see such a behavior in 2MASS or WISE photometry.For higher companion masses, this effect would be more pronounced.Thus, we rule out a main sequence companion for this system.Depending on the composition, the maximum mass of a white dwarf ranges between 1.35−1.40M ⊙ (Caiazzo et al. 2021).Fig. 5 indicates that for nearly edge-on orbits there is a possibility of the companion being a massive WD.For all moderate inclinations we expect the companion to be more massive than any reasonable limit on the WD mass, and therefore it would need to be a neutron star.
Along with the RVs we also use the lightcurve data to constrain the mass of companion by simulating the lightcurves and RV variation.To simulate the periodic ellipsoidal variation we use PHOEBE12 (Prša 2018).We use the stellar and orbital parameters that we derived in earlier sections as inputs to PHOEBE and simulate the lightcurve and the RV variation.To make this process computationally efficient, after phasefolding the lightcurve we bin it into two-hour time bins using lightkurve, reducing the total number of data points from 821 to 36.
If the secondary is a compact object, then for M 2 ∼ 0.9 M ⊙ , the upper limit on the radius of the secondary of about 0.008 R ⊙ can be obtained from the mass-radius relation of WDs (Chandra et al. 2020), and the flux change due to the eclipse of primary due to secondary would be about 0.02% for edge-on orbits, which is much smaller than the sensitivity of our data.Thus, as mentioned earlier, we ignore eclipses.We can also ignore effects of microlensing because their effects are of the same order as eclipses and hence negligible (Marsh 2001).Limb darkening and gravity darkening are second order effects and given the small amplitude of variation and noisiness of our data their effects are nearly negligible as well.Thus we model the surface using the default CK2004 atmospheric models (Castelli & Kurucz 2003) provided in PHOEBE, which should work well enough to simulate gravity and limb darkening.Our target has temperature T > 8000 K and thus we model it as having purely radiative atmosphere with gravity darkening coefficient β = 1 (von Zeipel 1924;Claret 1999).We simulate the lightcurve in the ZTF g and r bands and fit them to observed lightcurve.
Since the effects of limb darkening and gravity darkening are small, the dependence of lightcurve on the temperature should also be minimal, and the PHOEBE-generated lightcurves showed this behavior.Therefore, we fix the temperature of the primary to the spectrophotometric value previously obtained for ease of computation with little downside.We leave the surface gravity, radius of the primary, mass of the secondary and the inclination as free parameters.In addition to these, we also leave the systemic velocity, v γ , and initial phase of the orbit, t 0,sup−conj , as free parameters.

Combined Fit
With all the datasets assembled, we finally perform joint fitting where we add all the loglikelihoods and the associated log-priors, and maximize the resulting likelihood using emcee.Each dataset constrains different stellar and orbital parameters in a distinct way.The spectroscopic fit constrains T eff , log g, [Fe/H] independent of extinction.The lightcurve and RV fits also constrain log g and radius of the primary independent of extinction, and also constrain M 2 , inclination.The photometric fit is most sensitive to extinction and constrains T eff , log g, [Fe/H], R, E(B − V ).
We explore two solutions centered around the two local minima we found from the initial spectrophotometric fitting.The result from the combined fit is shown in Fig. 1 and Fig. 4. Looking at the final fits to the lightcurve and based on the theoretical expectations and empirical data on such binary systems (e.g., Fig. 8) we present the solution with the lower extinction value as the fiducial solution.The results are tabulated in Table .2. We explain this choice further in Sec.4.The posterior distributions for stellar parameters are shown in Fig. 6.The mass of the primary is found to be M 1 = 0.18 M ⊙ and the mass of the secondary to be M 2 = 1.19 M ⊙ , with 1σ upper limit of 1.41 M ⊙ and 3σ lower limit around 0.9 M ⊙ as shown in Fig. 7, suggesting a massive WD or a neutron star companion.The 1-sigma errors reported are calculated at 16th and 84th percentile of the posterior distributions.

DISCUSSION
In this paper we present the detection of a pre-ELM WD in a binary with a massive compact companion -likely a massive white dwarf or a neutron star.We perform a joint analysis of the radial velocity curve and of the photo- † The reported errors are numerical and could be underestimated.A more conservative error estimate would be 100 K and 0.1 dex -approximately half the BOSZ model grid spacing metric lightcurve of the visible pre-ELM WD to determine the allowed range of the masses of companion and find it to be 1.19 +0.21 −0.14 M ⊙ .The formation of such systems has been theoretically studied in detail as ELM WDs were found in binaries with pulsars and white dwarfs (Liu & Chen 2011;Shao & Li 2012;Istrate et al. 2014a;Li et al. 2019).We follow the various evolutionary tracks derived in these works and piece together the evolutionary pathway of our target.Both ELM WD-WD double degenerate binaries and ELM WD-NS binaries can form via two channels -stable Roche lobe over- In white we show the 1-sigma contour.
studied, especially with neutron star companions (Istrate et al. 2016), common envelope ejection channel is unlikely to produce the orbital period and primary mass combination that we observe.Thus, regardless of whether the unseen companion in our target is a massive white dwarf or a neutron star, we suggest that the system formed via a stable Roche lobe overflow from the main sequence star (which later became the primary pre-ELM WD) to the compact companion.Based on final orbital period, masses and evolutionary tracks, we deduce that the system must have arisen from a binary with initial orbital period of the same order as final orbital period, with the initial mass of primary between 1.2-2 M ⊙ .So far we have guessed the initial properties of the binary regardless of the nature of the companion.The initial companion mass is sensitive to the nature of the companion.Following Li et al. (2019), we find that if the companion is a WD, then its initial mass pre-Roche lobe interaction is around 0.6-0.8M ⊙ .WDs less massive than this range cannot accrete enough mass to grow to the observed mass of > 1M ⊙ because the mass accumulation is limited by the hydrogen burning rate on the surface of the white dwarf and by the resulting stellar winds.WDs which are more massive than this range cannot appreciably gain mass due to hydrogen-shell flashes causing inefficient mass accumulation.If the companion is a neutron star, results from Shao & Li (2012) suggest that a 1.1-1.2M ⊙ neutron star in a binary system with properties similar to those of our target can accrete about 0.3 M ⊙ worth of mass over the course of the evolution of ELM WD and lead to a system like ours.While the exact parameter values depend on various factors such as magnetic braking index, metallicity and, more importantly, on the exact nature of the companion, these are the representative values characteristic of systems such as our target and gives ideas about plausible scenarios such systems could have evolved from.
If the companion is a neutron star, the accretion of 0.3 M ⊙ mass would spin up the neutron star and lead to a milli-second pulsar with a spin-period of few milli-seconds (Liu & Chen 2011).The system is not detected in the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) at 1.4 GHz down to a ∼ 2.5 mJy flux level.X-ray signals can be expected from the surface of pulsars, especially from the magnetosphere of their polar caps in the case of recycled pulsars (Zavlin et al. 2002;Zavlin 2006;Kilic et al. 2011).Looking at data from X-ray survey ROSAT (Truemper 1982) we find no X-ray signal near our target.However, these are all-sky surveys with short exposure times and consequently high flux limits and thus only the brightest pulsars would be visible (Danner et al. 1994;Becker et al. 1996).Using ROSAT flux limit of 10 −13 erg cm −2 s −1 and assuming zero extinction in X-ray wavelengths we obtain upper limit on X-ray luminosity to be L X < 3.5 × 10 31 ergs −1 .Typical milli-second pulsars have X-ray luminosities be-low this limit (Bogdanov et al. 2006;Lee et al. 2018) and therefore would remain undetected.The presence of extinction would make their detection in ROSAT even more difficult.
The final stage of the close binary evolution through Roche lobe overflow results in a relationship between the primary masses and the orbital periods.These values are shown in Fig. 8, where we compare our target with ELM WD-WD binaries from the ELM survey (Brown et al. 2020) and and ELM WD-NS binaries from Gao & Li (2023).The short-period binaries (P < 10 hrs) with a broad distribution of the primary masses are from the common envelope channel, whereas the Roche lobe channel forms a relatively narrow locus of objects with a strong mass-period dependence.Our target with its period of P = 35.87hours and M 1 = 0.18 M ⊙ falls on the long-period end of the expected Roche lobe channel track.For comparison, we show the theoretically expected mass-period relations from binary simulations consisting of a neutron star and a WD primary.
The gravitational wave merger timescale for our target is of order a few hundred billion years.First, the primary will cool down, contract, and settle onto the WD cooling track in a time of the order of 0.2-2 Gyr (Istrate et al. 2014b).Comparing the two time scales we conclude that this system will end up as a stable binary of an ELM WD and a massive compact object whose orbital period will remain at its current value for much longer than the Hubble time.This justifies the comparisons of our pre-ELM WD with ELM WDs from the literature.
While we have assumed that the primary is a pre-ELM WD based on the spectrum and the surface gravity, we verify it using the derived parameters.We use the derived parameters to determine extinction-corrected absolute magnitudes or, equivalently, the luminosity of the primary which we find to be 1.5L ⊙ .Using this value and the mass-luminosity relation for main-sequence stars (Duric 2003;Eker et al. 2018) we estimate that a main-sequence star of equivalent brightness should have mass around 1.1 M ⊙ , which is much greater than the derived mass and therefore is inconsistent with the observed log g value.In addtion to this, the primary is too cold and much less massive than an sdB star (T sdB > 20000 K and M sdB ∼ 0.5 M ⊙ ) (Heber 2009).We conclude that despite its position on the color-magnitude diagram (Fig. 2) which is very close to the main-sequence track, the primary is indeed a pre-ELM WD.We also see that the primary has an anomalously large luminosity (or, equivalently, a larger radius) compared to pre-ELM sample from El-Badry et al. ( 2021), including those which are still mass transferring.This can be explained in the context of the Roche lobe channel as well: our target has a much greater orbital period com-pared to their sample, and consequently during the binary evolution the Roche-lobe mass transfer is expected to terminate earlier leading to a larger radius.One of the major sources of error in our analysis is the estimation of extinction.As we have discussed, we tried to take into account the various uncertainties associated with it.The uncertainty in extinction along line of sight and the degeneracy between the color excess and other stellar parameters gives two local χ 2 minima in the multi-dimensional parameter space of the fit.While high values of extinction are not ruled out based on the value of the χ 2 or of the likelihood, looking at Fig. 8, we find that these leads to solutions which are physically disfavored.If we allow the high-extinction solutions, we expect a higher primary mass than for the lowextinction solutions.Based on the Fig. 5, we deduce that this would strengthen the case for a neutron star companion.
In addition to the extinction, the other uncertain part is the polynomial fitting of the continuum.While we have presented results for the polynomial of order 6, we refit the spectrum by increasing the order of the polynomial to make sure this is a stable solution.We find that the final primary mass does not change appreciably and remains within 1-σ of our quoted value.We therefore present the results from the fitting with the lowest order polynomial which gives a stable solution.We also performed our analysis by fitting continuum normalized absorption lines alone instead of the full spectral fitting.This technique was sensitive to the choice of absorption lines and to the the continuum normalization.We did manage to recover our solution -albeit with a higher uncertainty, making it less reliable than our fiducial solution.
As seen in El-Badry et al. ( 2021), there is a possibility of emission lines being present in the spectrum due to the nature of the such systems which involves mass transfer.We see a faint sign of such emission at higher wavelengths near Paschen lines.However, because they are so faint, and the primary is not expected to be mass transferring at the moment, this could be due to noisy data.
As has been previously observed, the eccentricity for systems similar to ours is very small (Edwards & Bailes 2001;Istrate et al. 2014a).The binary must have gone through a common envelope phase before the more massive companion formed a compact remnant (Li et al. 2019) resulting in a circular orbit.While one would expect that the formation of a neutron star would give a kick to the primary leading to an eccentric orbit, the subsequent Roche-lobe overflow would circularize the orbit, as has been observed.The phase-folded RV curve and the lightcurve we obtain are consistent with a circular orbit.In principle, the eccentricity could be left as a free parameter in the analysis.However, given the strong theoretical footing for a nearly zero eccentricity, we do not expect the results to change in a significant way and leave it at zero.thors and do not necessarily reflect the views of the National Science Foundation.
Based on observations obtained with the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project.ZTF is supported by the National Science Foundation under Grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories.Operations are conducted by COO, IPAC, and UW.
This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.
Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, the National Science Foundation, and the Participating Institutions.SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah.The SDSS web site is www.sdss.org.
SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Carnegie Institution for Science, Chilean National Time Allocation Committee (CN-TAC) ratified researchers, the Gotham Participation Group, Harvard University, Heidelberg University, The Johns Hopkins University, L'Ecole polytechnique fédérale de Lausanne (EPFL), Leibniz-Institut für Astrophysik

Figure 1 .
Figure 1.The flux-corrected coadded SDSS spectrum (green line) and archival photometry (orange points) are shown along with the best-fit theoretical spectrum (blue line).The bottom panel shows the spectral region covered by the SDSS spectroscopic data.

Figure 2 .
Figure 2. Color-magnitude diagram.Our target is marked in red.We also show the extinctioncorrected position.The blue points are ELM WDs from Brown et al. (2020) and the green points are the pre-ELM WDs from El-Badry et al. (2021).The sequence in the lower left, below the ELM candidates, is the normal WD track.Our target is just below the main sequence.

Figure 3 .
Figure 3.The normalized periodograms from various datasets: blue for TESS photometry, red for ZTF photometry, grey for spectroscopy.

Figure 5 .
Figure5.The parameter space of masses of the binaries allowed by the binary mass function.M 1 is the mass of the visible star, M 2 is the mass of the invisible companion, and i is the orbital inclination counted from the plane of the sky.The colored regions are allowed while the white region in the bottom right would require an unphysical sin i > 1.We also show the Chandrashekhar limit (1.4 M ⊙ ) for M 2 .

Figure 6 .Figure 7 .
Figure 6.Posterior distribution from the joint parameter fitting of spectroscopic, photometric, RV and lightcurve datasets flow channel or common envelope ejection channel -both of which involve the massive compact companion accreting mass from a main sequence star (Nelson et al. 2004; Istrate et al. 2014a; Li et al. 2019).Following results from Li

Figure 8 .
Figure 8. Mass vs. orbital periods for ELM WDs binaries is plotted along with our target.Binaries with WD companions are selected from Brown et al. (2020), while binaries with NS companions are selected fromGao & Li (2023).The solid line is the theoretical expectation for the Roche lobe channel fromTauris & Savonije (1999) while the dashed line is fromLin et al. (2011).
Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania State University, Smithsonian Astrophysical Observatory, Space Telescope Science Institute (STScI), the Stellar Astrophysics Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Illinois at Urbana-Champaign, University of Toronto, University of Utah, University of Virginia, Yale University, and Yunnan University.Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, the National Science Foundation, and the Participating Institutions.SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah.The SDSS web site is www.sdss.org.SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Carnegie Institution for Science, Chilean National Time Allocation Committee (CN-TAC) ratified researchers, the Gotham Participation Group, Harvard University, Heidelberg University, The Johns Hopkins University, L'Ecole polytechnique fédérale de Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania State University, Smithsonian Astrophysical Observatory, Space Telescope Science Institute (STScI), the Stellar

Table 2 .
Observed and derived parameters for SDSS J0229+7130