A seven-Earth-radius helium-burning star inside a 20.5-min detached binary

Binary evolution theory predicts that the second common envelope (CE) ejection can produce low-mass ( 0 . 32 − 0 . 36M ⊙ ) subdwarf B (sdB) stars inside ultrashort-orbital-period binary systems, as their helium cores are ignited under nondegenerate conditions. With the orbital decay driven by gravitational-wave (GW) radiation, the minimum orbital periods of detached sdB binaries could be as short as ∼ 20 minutes. However, only four sdB binaries with orbital periods below an hour have been reported so far, while none of them has an orbital period approaching the above theoretical limit. Here we report the discovery of a 20.5-minute-orbital-period ellipsoidal binary, TMTS J052610.43+593445.1, in which the visible star is being tidally deformed by an invisible carbon-oxygen white dwarf (WD) companion. The visible component is inferred to be an sdB star with a mass of ∼ 0 . 33M ⊙ , approaching that of helium-ignition limit, although a He-core WD cannot be completely ruled out. In particular, the radius of this low-mass sdB star is only 0.066 R ⊙ , about seven Earth radii, possibly representing the most compact nondegenerate star ever known. Such a system provides a key clue to map the binary evolution scheme from the second CE ejection to the formation of AM CVn stars having a helium-star donor, and it will also serve as a crucial verification binary of space-borne GW detectors in the future.

the Lomb-Scargle periodogram (LSP) computed from the TMTS light curve.The vertical dashed line indicates the frequency corresponding to maximum power ( f max ).The purple dot-dashed line represents the confidence level of 0.1%, and the red arrow shows the frequency corresponding to the orbital period (i.e., f max /2).
Since the beginning of minute-cadence observations with Tsinghua University -Ma Huateng Telescopes for Survey (TMTS) 1,2 , we have discovered a dozen unusual short-period objects 3,4 in the Galaxy.TMTS J052610.43+593445.1 (J2000 coordinates α = 81.5434,δ = +59.5792;hereafter J0526) is a newly discovered variable star with a dominant photometric period of only 10.3 min (see Fig. 1).The periodicity was cross-checked by photometric observations from the Zwicky Transient Facility (ZTF) 5,6 and the Yunnan Faint Object Spectrograph and Camera (YFOSC) mounted on the Lijiang 2.4 m telescope (LJT) 7,8 (Fig. 2).Time-resolved spectroscopic observations from the Keck I Low-Resolution Imaging Spectrometer (LRIS) 9,10 and the Gran Telescope Canarias (GTC)/Optical System for Imaging and low-Resolution Integrated Spectroscopy (OSIRIS) 11 yielded a dozen single-line spectra with various radial velocities (RVs; Fig. 3).The RV curve is modulated by a 20.5 min period and reaches its peaks and valleys at the phases of maximum light (Fig. 2), which proves that J0526 is an ultracompact ellipsoidal binary.The unequal maxima in the light curves are caused by the relativistic Doppler beaming effect 12,13 of the visible component, consistent with its large RV amplitude.This object was also recently identified as a candidate verification binary (ZTF J0526+5934) of gravitational waves (GWs) by the ZTF DR8 database and Gaia EDR3 catalog 14 .Phase-folded RV curve and double-band light curves for J0526.Upper panel: the RV curve derived from Keck/LRIS and GTC/OSIRIS observations.The dot-dashed line is the best-fit sinusoidal model.Middle and lower panels: the gand r-band phase-folded light curves provided by LJT and ZTF.The purple solid lines represent the best-fit light-curve models obtained from the ellc package 15 .Unequal maxima are caused by the relativistic Doppler beaming effect 12,13 .Orbital phase φ = 0 represents the epoch of superior conjunction when the visible star is closest to the observer.The error bars represent 68% confidence intervals throughout this paper.
Although J0526 was included in white dwarf (WD) catalogs 16,17 , the probability of being a WD (P WD ) given by the probability map 18 is only P WD = 0.0046 16 , suggesting that J0526 should have large differences from those WDs.Thus, we present a detailed analysis of J0526 in this paper.Since the nature of noneclipsing binaries is usually not well constrained from light curves alone, the physical parameters of J0526 were determined by the combination of spectroscopy, broad-band spectral energy distribution (SED), RV curve, and multicolor light curves.According to the prevailing workflow for the analysis of ellipsoidal binaries, the properties of visible stars are obtained before the orbital solutions 19,20 .Prior knowledge of the visible component helps determine the inclination of the orbital plane and the mass of the invisible component from the light curves and RV curves.We refer to the invisible component of this binary as J0526A and the visible component as J0526B.

Atmospheric parameters
As shown in the dynamical spectra (Fig. 3), the Balmer lines and He I λ 4471 show synchronous shifts against the orbital phase, favoring that both H and He features arise from the visible star in the binary system.Because there are not any significant H/He lines tracing the motion of the invisible star, the invisible component must be very faint and is assumed to be negligible in the spectral fit below.As no emission lines are visible in the spectra, mass accretion should not occur in the two components of J0526, and we thus assume that the binary system is still detached.In order to verify this assumption, we further compared the radius of the visible star with its Roche-lobe size in the following discussion.
Since the H/He absorption lines in the spectra carry key information about the atmospheric properties of the visible star, we fitted all of the GTC/OSIRIS spectra using non-local thermodynamic equilibrium (NLTE) spectral models obtained from TLUSTY and SYNSPEC software 21,22 (see Methods).The best-fit model reproduces well the main Balmer lines and He I λ 4471 seen in the observed spectra, which gives estimations of the effective temperature T eff = 25480 ± 360 K, surface gravity log g = 6.355 ± 0.068, helium abundance log y = log N He /N H = −2.305± 0.062, and projected rotational velocity ν sin i = 220 +140 −90 km s −1 for the visible component.Among the large samples obtained from current spectroscopic surveys, hydrogen-rich spectra contaminated by He I lines are very rare for WDs [23][24][25] , but common for sdB stars [26][27][28] .Observationally, the helium abundance of hot subdwarfs is overall positively correlated with their effective temperatures.However, this correlation tends to show two distinct branches for He-rich and He-weak sequences, especially at the high-temperature end 29,30 .In comparison with the existing sdB sample, J0526B locates on the low-temperature end of He-rich sequence in the T eff − log y diagram (see Fig. 4).

Stellar radius and mass
The broad-band spectral energy distribution (SED) is a very useful tool for constraining the stellar radius and luminosity if the distance is available and reliable.With prior knowledge of the effective temperature and surface gravity of the visible star [i.e., T eff,B and log (g) B ] derived from the spectroscopic analysis above, we fit the SED (Fig. 5) using archival multicolor photometry and the distance inferred from the Gaia DR3 31 parallax (see Methods).Since this object is not included in the UV-source catalog of the Galaxy Evolution Explorer (GALEX) 32 , we also proposed UV observations with Swift/UVOT.The best-fit model suggests a radius R B = 0.0681 +0.0059 −0.0049 R ⊙ and a bolometric luminosity L bol = 1.70 +0.31 −0.24 L ⊙ for the visible star, and then updates the effective temperature and surface gravity (see Table 1).The model also yields an estimate of the line-sight-of extinction as E(B −V ) = 0.387 +0.008 −0.009 mag, well consistent with the Galactic value of E(B −V ) = 0.385 mag as queried from a three-dimensional dust extinction map 33 .With the surface gravity and radius, we computed the mass of the visible star as  Relation between He abundance and effective temperature for hot subdwarfs.All atmospheric parameters of hot subdwarfs are taken from reference 26 .The purple dashed line and green dotted-dash line represent the correlation for the He-rich sequence 29 and the He-weak sequence 30 , respectively.

5/24
M B = 0.361 +0.094 −0.073 M ⊙ using Newton's law of gravity.The radius and mass obtained from the SED suggest that J0526B is more likely a hot subdwarf with an extremely thin hydrogen-rich envelope rather than a helium-core WD (see mass-radius relation below)..Orange diamonds indicate the synthetic fluxes derived from the model spectrum and transmission curves.Lower panel: relative residuals.

Orbital dynamics
The Doppler shift of spectral lines provides key clues for us to investigate the orbital dynamics of J0526.By assuming that the orbit is circularized, the RV curve can be modeled well by a sinusoidal curve (Fig. 2), with a RV semi-amplitude of K B = 559.6 +6.4 −6.5 km s −1 inferred for the visible component.Hence, the mass function of the invisible component was computed following where M A and M B represent the masses of the invisible and visible components (respectively), i is the inclination of the orbital plane, and G is the gravitational constant.Since the RV semi-amplitude and orbital period are both well-constrained from the observations, the mass function bridges a tight relation between the masses of binary and the inclination angle, and thus aids binary parameter estimation from the light-curve fit below.

Ellipsoidal variations and compact object
Ellipsoidal modulation in the light curves is induced by tidal deformation and rotation of the visible component.Owing to the synchronization between the rotation and orbital motion, different geometric cross-sections of the visible star emerge throughout the orbit.The large amplitude of the ellipsoidal modulations suggests that the visible star almost fills its Roche lobe (e.g., f R = R B /R L,B ≳ 0.80).We modeled the gand r-band light curves of J0526 using the ellc package 15 (see Methods).The values of R B , M B , and f (M A ) obtained above were included as prior parameter distributions.The Doppler beaming effect was also included in the model for offsetting the unequal maxima.Gravity/limb-darkening coefficients and Doppler beaming factors were obtained by interpolating the grids 47 with the surface parameters derived from the spectroscopy and SED.The best-fit model  6.The characteristic strains of J0526 accompanied with dozens of verification/detectable binaries of GWs.The characteristic strains were calculated from the component masses and distances provided from reference 44 .The blue dashed line and red dotted-dash line represent the detection sensitivity curves from LISA 45 and TianQin 46 , respectively.The LISA sensitivity curve here includes the instrumental noise the foreground confusion noise, while the TianQin sensitivity curve includes only the instrumental noise.
gives i = 68.2+3.7 −5.2 deg and M A = 0.735 +0.075 −0.069 M ⊙ , and updates other physical parameters with Bayes' theorem (see Table 1).The mass suggests that the invisible star is a carbon-oxygen (CO) WD.Through the mass-radius relation of CO WDs 48 , we estimate that the radius of the invisible star is about 0.011 R ⊙ , favoring the noneclipsing scenario for J0526.With the updated radius of J0526B (0.0661 ± 0.0054 R ⊙ ), we derived the projected rotational velocity (ν B sin i) cal = 2πR B /P orb = 216 +18 −19 km s −1 , consistent with the result obtained above from the spectroscopy.
Given the ultracompact orbit and relatively high masses of the two components, J0526 is predicted to be detected by Laser Interferometer Space Antenna (LISA) 49 within the first three months of its operation (see the Methods) and thus will serve as a verification binary of GWs in the future.The GW characteristic strains of J0526 and dozens of other verification/detectable binaries of GWs are presented in Fig. 6.

Beyond the strip selection effects
Fig. 7 shows that J0526B is located at an interlaced zone between hot subdwarfs 26 and (extremely) low-mass WDs 50 in the Kiel diagram (T eff − log g diagram).Since the cooling sequences of WDs are widely distributed on the blue side of the main sequence, the surface gravities and effective temperatures of hot subdwarfs are compatible with those (pre-)WDs.We can further cross-check the nature of J0526B by comparing the spectrophotometric parallax, derived from the atmospheric parameters and hypothetical nature, against the astrometric parallax obtained from Gaia DR3 31 .By assuming that J0526B is a helium-core WD, we estimated a spectrophotometric parallax ϖ spec = 1.542 ± 0.136 mas for J0526(see Methods), which shows a deviation from the Gaia DR3 parallax ϖ Gaia = 1.183 ± 0.091 mas by 2.2σ .
Following the theoretical predictions from binary evolution theory and binary population synthesis 51,52 , the second common envelope (CE) channel is responsible for the formation of sdB binaries with very short orbital periods (typically P orb < 1 hr) 52 .Because these sdB stars are produced from nondegenerate He cores, their masses are expected to be only ∼ 0.33 M ⊙ 52, 53 .However, the extreme horizontal branch (EHB) in Kiel diagram, corresponding to hot subdwarfs with canonical masses of ∼ 0.48 M ⊙ , is widely applied to confirm the natures of hot subdwarfs.This selection effect (so-called strip selection effect 52 ) would lead to a systematic absence of low-mass sdB stars with very short orbits.Additionally, for a detached sdB binary having an orbital period of only 20 min, the hydrogen-rich envelope must be extremely thin (e.g., 10 −6 M ⊙ 54, 55 ) to avoid Roche-lobe overflow (RLOF).With this prior knowledge, we ran the Modules for Experiments in Stellar Astrophysics (MESA) 56 code to reproduce evolutionary tracks for sdB stars in extremely short-orbital-period binary systems (see Methods).Although these tracks cover only a very small area relative to the regions of hot subdwarfs or WDs, as shown in Fig. 7, the atmospheric parameters of J0526B overlap exactly on the sdB tracks of 0.32 − 0.33 M ⊙ .By assuming that J0526B is an sdB star of 0.33 M ⊙ , the spectrophotometric parallax can be re-estimated as ϖ spec = 1.30 ± 0.10 mas, approximating the astrometric parallax.  26,50The blue dotted-dash lines are theoretical evolutionary tracks of helium-core WDs 35 , while the shell flash loops are clipped for the clarity of image.The purple solid lines represent the evolutionary tracks of low-mass helium-core burning stars with a very thin hydrogen envelope of 10 −6 M ⊙ (see Methods).We overplot the zero-age helium main sequence (ZAHeMS; see Methods), zero-age extreme horizontal branch (ZAEHB), and terminal-age extreme horizontal branch (TAEHB) 57 as black dashed lines. 9/24

Mass-radius relation
Since diverse equations of state (EOS) from different classes of stars lead to distinct mass-radius relations, the mass-radius diagram (MRD) [58][59][60] is a valid tool for distinguishing different types of stars.An extended version of the MRD toward the low-radius end is shown in Fig. 8, where the three dominant types of stars include CO-/He-core WDs, main-sequence (MS) stars, and hot subdwarfs.These three types of stars are located in separate regions and they hardly overlap in the MRD.Being supported by electron degeneracy pressure, WDs tend to have smaller radii at larger masses, the opposite of nondegenerate stars.The hot subdwarfs cover a wide area in the MRD owing to diverse hydrogen-envelope masses generated from different initial binaries and evolutionary channels.The observational constraints from spectroscopy, the SED, and light curves support that J0526B is located exactly at the lower tip of the hot-subdwarf domain.In other words, J0526B can be a low-mass core helium-burning (CHeB) star with an extremely thin hydrogen envelope, as also suggested by the analysis of the Kiel diagram.As a hydrogen-exhausted star inside the 20.5 min orbit, the size of J0526B is smaller than that of all previously known nondegenerate stars, even those brown dwarfs and gas planets 59,61,62 .But J0526B has an average density ∼ 1200 times greater than that of the Sun!  59 , B20 50 , and S22 63 , respectively.Six sdB stars inside ultracompact binaries having orbital periods below 100 min 44,[64][65][66][67][68][69] are shown.For comparison, 18 WD components inside 11 ultracompact binaries 64,66,[70][71][72] discovered by photometric surveys are also included.We overplot theoretical mass-radius relations for CO WDs, He WDs, MS stars, and HeMS stars.The degenerate and nondegenerate models are distinguished by solid and dashed lines, respectively.The purple area above the HeMS corresponds to the hot subdwarfs, i.e., the core-helium-burning (CHeB) stars with hydrogen-rich envelopes, in which thicker envelopes lead to larger stellar radii.A sequence for CHeB stars with an extremely thin hydrogen-rich envelope (M H,env = 10 −6 M ⊙ ) is highlighted using a pink dashed line.We present the constraints from surface gravity and SED, derived from current observations, accompanied with the mass-radius relations of an 80%/90% Roche-lobe-filling star inside a 20.5 min orbital-period binary using Paczyński's approximation 73 (see Methods).The radii of Jupiter and 2MASS J0523−1403 61 (the smallest star among known MS stars) are labeled on the left of the figure.

Second CE ejection and AM CVn stars
Following the theoretical predictions from the second CE ejection channel of sdB stars 51,52 , some sdB stars are born from a pair of WD and (evolved) MS stars, which is the stellar remnant that survived the first CE ejection 74 or stable RLOF 75 .In this channel, the WD companion has a very small radius and can penetrate deeply into the CE before CE ejection, which allows the formation of a sdB binary with a very short orbital period.In particular, if the MS component has an initial mass larger than the critical mass for a star to experience the helium flash at the end of its first giant branch (FGB, also red-giant branch), e.g.M MS,0 ≳ 2.0 M ⊙ 51, 53 , its primary fusion reactions are through the CNO cycle instead of the proton-proton chain during its MS stage.Consequently, the MS component can retain a higher central temperature and thus produces a more massive, nondegenerate helium core (∼ 0.2 M ⊙ ) when leaving the MS, compared to those with lower initial masses.Given the higher central temperature, the helium core can be potentially ignited under nondegenerate conditions 51,76 even if the envelope is lost during passage through the Hertzsprung gap.SdB stars formed through this subchannel would have very low masses (M sd = 0.32 − 0.36 M ⊙ ) 52,53 .Envelopes of the Hertzsprung-gap stars are generally more tightly bound than those of stars near the tip of the FGB.Consequently, the inner binary system must release more orbital energy to counterbalance the binding energy, resulting in an extremely short final orbital period (typically of order tens of minutes) after the second CE ejection.J0526 is such an ultracompact binary.Its extremely short orbital period, WD companion, and low-mass sdB component exactly follow the theoretical predictions from the second CE ejection channel of sdB stars 51,52 .By assuming that the mass of J0526B is equal to 0.33 M ⊙ , its Roche-lobe radius is estimated as about 0.079 R ⊙ , which is properly wider than the size of J0526B, consistent with the above inference that this binary system is currently detached.With orbital contraction driven by gravitational-wave radiation (GWR), after about 1.5 million years, J0526B will overflow its Roche lobe and transfer mass toward the WD at an orbital period of around 14 min (see Fig. 9), leading to the formation of an AM CVn star through the helium-star channel 77,78 .Owing to the nondegenerate nature of J0526B, its radius will shrink in response to mass loss induced by RLOF, supporting further orbital contraction driven by GWR.Since the mass transfer quenches the helium burning, J0526B will begin a transition to a degenerate state, e.g.becoming a He-core WD.When the electron-degeneracy pressure becomes dominant, J0526B will reach the minimum orbital period of ∼ 9 min and it will start to expand with its mass loss, leading to an increasing orbital period as predicted by binary evolution theory.Ultimately, the donor star in such an AM CVn system 11/24 either evolves into a planet orbiting the WD companion 78,79 , or is tidally disrupted by the WD accretor when its mass becomes smaller than ∼ 0.002M ⊙ 80 .In summary, J0526 could be the shortest-orbital-period single-degenerate detached binary, which provides crucial observational evidence supporting the complete evolutionary scheme ranging from initial binary MS stars, to MS+WD binary, to sdB+WD binary, to AM CVn star 77,78 .With the operations of the Large Synoptic Survey Telescope (LSST) 81 , Wide Field Survey Telescope (WFST) 82 , and space-borne gravitational wave observatory 49,83 , more previously unknown extremely-shortorbital-period sdB binaries will be discovered and thus aid our understanding of the formation of sdB stars and AM CVn stars.

Photometric observations and orbital period
In the first two-year survey, TMTS has discovered more than 1100 variable stars with periods shorter than 2 hr 4 .J0526 is one of the shortest-period variable stars in the catalog.Its 10.3 min periodicity was first revealed by the Lomb-Scargle periodogram (LSP) 84 derived from the 12 hr minute-cadence observations on 18 December 2020 (UTC dates are used throughout this paper; Fig. 1).The TMTS Light-curve Analysis Pipeline automatically estimated the dereddened color (B p − R p ) 0 = −0.41± 0.02 mag and absolute magnitude M G = 6.86 ± 0.21 for J0526 using its embedded Gaia DR2 database 85 and DUSTMAPS Python package 86 .A color-magnitude diagram for J0526 with some sdB stars, low-mass WDs is presented in Fig. 10 , using the Gaia DR3 database 31 .The ultrashort period and extraordinary color drove us to trigger further photometric and spectroscopic observations of this object.  26,50 and are then cross-matched with the Gaia DR3 database 31 .The Gaia magnitudes and colors of models are computed using the bolometric correction tables from MESA Isochrones & Stellar Tracks (MIST) 87,88 .
The gand r-band photometric data were obtained from ZTF Public Data Release 14 (DR14) 5,6 and LJT/YFOSC observations 7,8 conducted on 19 December 2022.The LJT observations lasted for 60 min in r and 44 min in g, with a common exposure time of 30 s and a readout time of ∼ 3 s.All LJT photometric data were reduced according to standard procedures, including bias subtraction, flat-field correction, and cosmic-ray removal.For ZTF data, the measurements with cat f lag = 32768 were excluded.All Modified Julian Days (MJDs) in both ZTF and LJT data were converted into BJD TDB .All observed fluxes were normalized by average fluxes in each band.We computed the LSPs from light curves of ZTF and LJT/YFOSC observations, and thus confirmed the periodicity of J0526.Thanks to the long-term observational coverage from ZTF, we obtained a precise photometric period P ph = 10.2531213 ± 0.0000026 min from the g-band light curve, and thus the orbital period is P orb = 20.5062426± 0.0000053 min.Since the uncertainty in the orbital period is tiny, P orb was fixed to 20.5062426 min for the analysis below.

Spectroscopic observations
We observed two series of spectra for J0526 within two independent observations runs, one with the 10 m Keck-I telescope equipped with LRIS (blue grism 600/4000, R∼1000; red grating 1200/7500, R∼2000) 9,10 , and the other with the 10.4 m GTC plus OSIRIS instrument (grism R1000B, 2 × 2 binning, R∼1000 ) 11 .The Keck/LRIS spectra were observed at six sequential orbital phases on 23 September 2022, with an exposure time of 180 s for the first spectrum and 240 s for the others.A total of seven GTC/OSIRIS spectra were observed on 26 January 2023, with each having an exposure time of 180 s and a readout time of ∼ 25 s.Because of terrible weather conditions during the Keck I observations, the signal-to-noise ratio (SNR) of the Keck/LRIS spectra is significantly lower than that of the GTC/OSIRIS spectra.
The GTC/OSIRIS spectra were reduced following standard tasks in IRAF via the graphical user interface FOSCGUI, which designed to extract SN spectra and obtained with FOSC-like instruments.It was developed by E. Cappellaro, and a package description can be found at http://sngroup.oapd.inaf.it/foscgui.html.The raw images were first corrected for bias, overscan, trimming, and flat fielding, and subsequently one-dimensional (1D) spectra were optimally extracted from the 2D images.Wavelength calibration was performed using spectra of comparison lamps that were produced two days earlier than the observation night, while flux calibration was done via observations of spectrophotometric standard stars.These calibration images were taken with the same instrumental configuration and on the same night as the spectra of J0526.Finally, the J0526 spectra were fine-tuned with the coeval broad-band photometry data, and the broad absorption bands (e.g., H 2 O, O 2 ) due to Earth's atmosphere were removed using the spectrum of the standard star.The Keck/LRIS spectra were reduced through a dedicated pipeline LPipe 89 following similar procedures.

Bayesian inference
In order to constrain the physical parameters of J0526 from spectra, broad-band SED, RV curve, and light curves together, we linked the model parameters derived from each observational clue using the Bayes' theorem 90 , where p(θ |D) is the posterior distribution of model parameters θ with given data D, L (D|θ ) is the likelihood function (also sampling distribution) for D with given θ , and p(θ ) is the prior density distribution of model parameters θ .Following prevailing methods for resolving the nature of ellipsoidal binaries 19,20 , we first determined the physical parameters of the visible component, which provides better parameter constraints to obtain correct orbital solutions.The GTC spectra were fitted with non-local thermodynamic equilibrium (non-LTE) model atmospheres for obtaining atmospheric parameters (including log g, T eff ) of the visible component independently.Then these atmospheric parameters were taken as prior density distributions [p(θ )] to derive the radius and mass of J0526B from the broad-band SED.We also fitted the RV curve to obtain the RV semi-amplitude (K B ) and epoch of superior conjunction (T 0 ) when the visible star was closest to the observer.The model parameters [log (g) B , R B , K B , and T 0 ] were used to construct prior density distributions for light-curve modeling, and were refined by final posterior distribution.

Atmospheric model
Because the SNR of the LRIS spectra is too low (owing to poor weather) to yield correct atmospheric parameters for J0526B, we adopted only the OSIRIS spectra to derive atmospheric parameters.To reduce smearing effects and potential contributions from the invisible component, we fit all GTC spectra simultaneously (see Fig. 11), using metal-free non-LTE TLUSTY (v207) model atmospheres and synthetic spectra with SYNSPEC (v53) 21,22 .The best-fit model was obtained from the iterative spectral analysis procedure (XTGRID) 91 which applies a steepest-descent χ 2 minimization to simultaneously optimise all free parameters, including effective temperature, surface gravity, chemical abundances, and the projected rotational velocity.All comparisons were run globally using the spectral range 3780-6850 Å, and the model was piecewise normalized to the observations.In parallel, the RVs were also determined by shifting each observation to the model.The procedure converges once the relative changes of all model parameters and χ 2 drop below 0.5% in three consecutive iterations.Finally, parameter uncertainties were calculated by mapping the parameter space around the best solution.Notice that, when modeling the spectra of sdB stars, the systematic differences caused by replacing the metal-free non-LTE model atmospheres with the metal-line blanketed LTE model atmospheres (i.e., with [Fe/H]=−2 to +1) could be in the levels of ∆T eff ≈ ±500 K and ∆ log g ≈ ±0.05 [92][93][94] , which are comparable to the uncertainties quoted in our results.All atmospheric parameters of J0526B are summarised in Table 1.

Spectrophotometric parallax
A common method to test the hypothetical nature of sources is to compare their spectrophotometric parallaxes with available astrometric parallaxes.Following the approach of calculating the spectrophotometric parallax 50,95 , we first assumed J0526B to be an extremely low-mass (ELM) WD or a sdB star.For the case of an ELM WD, we can obtain its mass M ′ B = 0.237±0.018M ⊙ by interpolating the grids from evolutionary tracks of He-core WDs 35 using the T eff , log g obtained from spectra.The grids have included the evolutionary tracks of shell flash loops, and thus lead to multiple solutions at same region of T eff − log g diagram.The mass uncertainty here included the errors of surface parameters and multiple solutions of model tacks.For the case of an sdB star, its mass is assumed to be 0.33 M ⊙ based on the position of J0526B in the T eff − log g diagram (Fig. 7).With Newton's law of gravity and the Stefan-Boltzmann law, we estimated its stellar radius and luminosity from the atmospheric parameters and the assumed masses.The absolute magnitude in the Gaia G filter can be obtained from where M bol,⊙ = 4.74, and BC G is the bolometric correction for the G band.Hence, the spectrophotometric parallax was calculated by following where G = 17.563 ± 0.003 mag, and A G is the interstellar dust extinction in the G band.Both BC G and A G were obtained by interpolating the bolometric correction tables from MESA Isochrones & Stellar Tracks (MIST) 87,88 with the atmospheric parameters and solar metallicity.

Radial-velocity curve
As introduced above, the RVs were determined by shifting each observed spectrum to the TLUSTY/SYNSPEC synthetic spectrum.All RV measurements were corrected to the barycentric rest frame.For such a compact orbit, it is reasonable to assume that the 14/24 orbit is highly circularized (i.e., eccentricity e = 0).Therefore, the RV curve can be fitted by a sinusoidal function, and the likelihood function for RV measurements was calculated by where V r,model, j = K B sin [ 2π P orb × (t ν,obs, j − T 0 )] − γ is the model RV at time t ν,obs,j , and γ is the RV of barycenter of the binary system.Then t ν,obs,j , V r,obs,j , and σ ν,obs,j are the BJD, observed RV, and uncertainty obtained from the j-th spectroscopic observation, respectively.

Spectral Energy Distribution
The broad-band SED was constructed from Swift 37 Target of Opportunity (ToO) observations (UV-W2/M2/W1 bands) and archival photometry from Gaia DR3 31,42 (BP, G, and RP bands), Pan-STARRS 38,39 (grizy bands), ZTF 5,6 (gr bands), and AllWISE 40,41 (W1-W4 bands).The photometric measurements in three Swift-UV bands were obtained by running the HEASOFT (ver.6.31) command UVOTPRODUCT.Since J0526 was accidentally located in the bad area of the detector during the observation on 15 March 2023 (ID: 00015916001), we adopted the measurements from the observations on 22 March 2023 (ID: 00015916002).The optical and infrared photometric fluxes were directly obtained from the Virtual Observatory SED Analyzer (VOSA) 96 online tool.We did not use the upper-limit measurements in the WISE W2-W4 bands in the SED fitting.
We constructed the grid of synthetic photometry by sequentially integrating extinction factors and filter transmission curves over the TMAP 43 synthetic spectra.The extinction factors were derived from the Fitzpatrick extinction curve 97,98 with reddening law R V = 3.1, and the grid of reddening E(B −V ) spans 0.00-1.00mag with a step of 0.05 mag.The transmission curves were obtained from the Filter Profile Service of Spanish Virtual Observatory (SVO) 99 .The TMAP spectrum grid covers 20, 000 K ≤ T eff ≤ 66, 000 K with a step of 2000 K, and 4.50 ≤ log g ≤ 6.50 with a step of 0.25.We applied a 3D linear interpolation to approximate the synthetic flux at arbitrary [T eff , log g, E(B −V )] coordinates within the grid.
In order to consider additional flux uncertainties caused by ellipsoidal variations of J0526, a free parameter σ F,sys was included in the likelihood function, where F obs,j and σ F,obs,j represent the observed flux and uncertainty in the j-th photometry, respectively.F syn,j , the synthetic flux in the j-th band, is a function of T eff , log g and E(B −V ).

Light curves
In order to obtain the inclination angle ( i) of the orbital plane and thus calculate the mass of the primary M A through the mass function (Eq.1), we modeled all gand r-band light curves from both ZTF and LJT observations simultaneously, using the ellc (ver 1.8.7) package 15 .For an ellipsoidal binary, the flux variation can be approximated as where ∆F ell /F represents the fractional semi-amplitude of the ellipsoidal modulation, and a is the orbital separation; µ B and τ B are the Limb-darkening coefficient (LDC) and gravity-darkening coefficient (GDC) of the visible star, respectively.Grids of LDC and GDC were generated from the tables of reference 47 .We applied 2D cubic interpolations to approximate the g and r coefficients at arbitrary (T eff , log g) coordinates inside the grids.With the T eff and log g determined above, we obtained τ B = 0.45 and µ B = 0.31 for the g light curves, and τ B = 0.41 and µ B = 0.26 for the r band.The Doppler beaming effect leads to higher flux (φ ≈ −0.25) when the visible component is approaching the observer, than that measured in the case of moving away (φ ≈ 0.25).We set the beaming factor as b = 1.30 for the g band and b = 1.47 for the r band 47 .
In order to respect the RV semi-amplitude K B and epoch of superior conjunction T 0 obtained from the RV curves, we introduced a prior density distribution of the mass function f (M B ) (i.e., Eq. 1) to bridge a tight relation among M A , M B , and inclination i.In addition, the T 0 derived from the RV curve was also included as a prior distribution for the epoch of φ = 0.Because the mass of the visible star (M B ) in an ellipsoidal binary is poorly constrained by both orbital dynamics and ellipsoidal variations, we adopted the prior density distribution of M B derived from spectroscopy and the broad-band SED.
Since the photometric data were obtained with various instruments and filters, we introduced four free parameters σ f,sys,k (k = 1, 2, 3, 4) for offsetting the systematic errors in ZTF g, ZTF r, LJT g, and LJT r light curves, respectively.Hence, the likelihood function was expressed as where N k is the number of photometry points in the k-th light curve, and f ellc, j,k represent the ellc model flux for the k-th band at time t j,k .Then t j,k , f obs, j,k , and σ f,obs, j,k represent the BJD, observed flux, and uncertainty of j-th photometry point in the k-th light curve, respectively.

Mass-radius relation within Roche lobe
Given the orbital period and fillout factor (the ratio between stellar radius to its Roche-lobe radius, f R = R B /R L,B ), a mass-radius relation for the star 60,100 can be obtained from Kepler's third law, and the Paczyński approximation (preferred for q = M B /M A ≲ 0.8) 73 for the Roche-lobe radius (in units of orbital separation), where a is the orbital separation.With R B = f R R L,B = f R f (q) × a to bridge the relation between Eq. 9 and Eq. 10, the terms q in the two equations cancel each other out, leading to a q-independent expression (if q ≲ 0.8 ) for the mass-radius relation,

Galactic orbit
We computed the components of the Galactic velocity for J0526.We set the Sun at a distance of R 0 = 8.27 ± 0.29 kpc from the Galactic center 101 and its peculiar motion in relation to the Local Standard of Rest 102 at (U ⊙ ,V ⊙ ,W ⊙ ) = (11.1,12.24, 7.25) km s −1 .The rotational speed of the Milky Way at the Solar circle 101 was set to be V c = 238 ± 9 km s −1 .The computed Galactic velocity components resulted in (U = 43 ± 4 km s −1 , V = 233 ± 9 km s −1 , W = 5 ± 3 km s −1 ), suggesting a membership in the thin-disk category 103 .

Verification binary of gravitational waves
For ultracompact binaries, GWs are thought to dominate the orbital angular momentum loss (AML) and thus lead to the orbital contraction.Thanks to the compact orbit, the GWs generated from J0526 are expected to be detected by space-borne GW observatories, such as LISA 49 and Tianqin 83 .With the component masses, orbital period, and distance of J0526 (Table 1), we can estimate the detection SNR of its GW signal after three-month 104 and four-year observations, respectively.For a binary system, the GW strain amplitude can be calculated 105 by where c is the speed of light in a vacuum, d is the source distance, f GW = 2/P orb represents the GW frequency, and is the chirp mass.The characteristic strain h c is thus given by h c = √ f GW T obs A , where T obs is the integration observation time, typically 4 yr for LISA and TianQin.We present the characteristic strains for dozens of verification/detectable binaries with detection sensitivity curves from LISA 45 and TianQin 46 in Fig. 6.The chirp mass of J0526 is calculated using the binary masses given by the light-curve analysis, as shown in Table 1.Benefiting from its close distance, J0526 is one of the ultracompact binaries that can generate the strongest GW characteristic strains.For the first-three-month GW observation of LISA, the SNR of J0526 can reach ∼ 3 (see Table 1).
Owing to the absence of densely-sampled observations, the orbital contraction rate of J0526 is not available yet.By assuming that the AML is driven only by GW radiation, we can obtain a theoretical decay rate of the orbital period using Ṗorb Hence, the theoretical period derivative of J0526 is Ṗorb /P orb = −1.72

Evolutionary models
J0526 is an excellent object in developing gravitational-wave astronomy and investigating key process in binary evolution such as second CE ejection.In order to figure out the nature of J0526, we employed the stellar evolution code MESA 56, 106-109 to investigate its origin and final fate.

Single-star evolution models
According to the observational clues introduced above, J0526 probably consists of a CO WD and a low-mass sdB star.To construct models for the sdB star, we first create a series of He-pre-main-sequence stars, whose core masses range from 0.32 M ⊙ to 0.36 M ⊙ in a step of 0.01 M ⊙ , with He-mass fraction of 0.98 and metallicity of 0.02.The nuclear reaction network adopted in our simulations is "approx21.net"-a 21-isotope α-chain nuclear network.The "OPAL type 2" radiative opacities tables for a carbon-/oxygen-rich mixture with a base metallicity Z = 0.02 is employed in our simulations 110,111 .The adopted metallicity is consistent with that of Galactic thin disk membership.Some physical processes were not included in our models such as overshooting, rotation, diffusion, etc.
To generate the HeMS models (i.e., the CHeB models without H envelope), we create a series of pre-MS models with initial He abundance 0.98 and metallicity 0.02 (we refer to these models as "He-pre-MS" models) and evolve them until central-He ignition.
In order to generate the sdB models (i.e., the CHeB models with a very thin H envelope), we evolve the "He-pre-MS models" until their central temperatures reach logT c = 7.8, and then turn off all the nuclear reactions and implant different masses of a pure H shell onto the surface of the He core using a mass accretion rate of 10 −7 M ⊙ yr −1 .The mass of the H shell is in the range of 10 −6 to 10 −2 M ⊙ .After the formation of the H shell, we restore the nuclear reactions to evolve the sdB models forward with time.When the nuclear reactions are restored, all the sdB models undergo Kelvin-Helmholtz contraction and thus ignite the central-He burning.At the onset of central-He burning, their central temperatures reach logT c = 8.0, the central densities range from logρ c = 4.6 to 4.8, and the central He mass fractions are about 0.96.The sdB stars with core masses of 0.32 − 0.33 M ⊙ and H-envelope mass of 10 −6 M ⊙ are consistent with J0526B (see Fig. 7).
For the CHeB models with 0.32 − 0.33 M ⊙ , the stars will evolve toward the WD cooling sequence if the central He is exhausted, while other models (0.34 − 0.36 M ⊙ ) will experience unstable shell-He burning and finally become a WD when the shell He abundance cannot maintain further He burning.

Binary evolutionary models
The strong GWR from J0526 leads to orbital contraction and RLOF of J0526B in the future.In order to predict its final fate, we evolve binary models with an initial orbital period of 20.5 min.Among these models, the donor stars are the zero-age CHeB stars obtained above from single-star evolution.The accretor in these binaries is a 0.735 M ⊙ CO WD, which is treated as a particle in the models.To calculate the change rate of orbital angular momentum, we took into account contributions from both GWR and mass loss.The AML driven by GWR is introduced by using the standard formula 112 , where M WD and M sd are the masses of the WD accretor and sdB donor, respectively.Mass transfer begins after the donor fills its Roche Lobe.However, the WD cannot accumulate all the material transferred from the donor, while some material is lost from the WD through the stellar wind, leading to the mass loss and AML of the system.The AML caused by the mass loss can be calculated following the expression JML,WD = −(1 − η He ) Ṁacc ( where Ṁacc is the mass-accretion rate of the WD, which equals the mass-transfer rate from the donor to the WD.Referring to previous work [113][114][115][116] , the mass-accumulation efficiency of WD for helium (η He ) can be approximated as where Ṁlow,He and Ṁcr,He are the lowest and highest critical accretion rate (respectively), and Ṁst,He is the minimum accretion rate for stable He-shell burning.If Ṁacc ≤ Ṁlow,He , all of the accreted material will be blown away by the strong stellar wind, which means that the mass-accumulation rate of the WD is almost negligible; as Ṁlow,He < Ṁacc < Ṁst,He , the massaccumulation efficiency under weak He-shell flashes (η ′ He ) is adopted from the simulation results 117 ; for Ṁst,He ≤ Ṁacc ≤ Ṁcr,He ,

17/24
the He-shell burning is stable and thus all accreted material is accumulated on the surface of the WD (i.e., η He = 1); when Ṁacc > Ṁcr,He , the WD accumulates material with its extreme rate of Ṁcr,He , and excess material is blown away by the optically thick wind.The sdB stars in our models begin to transfer material to the WD after their RLOF.We evolved these models until the luminosities of the donors were below 10 −6 L ⊙ .All these models experience AM CVn phases and then the mass-transfer rates begin to decrease with increasing orbital periods, as can be seen from Fig. 9.The final outcome is probably a WD+planet system 78,79 or a single WD when the donor is tidally disrupted by the accretor 80 .

Figure 1 .
Figure 1.TMTS light curve and Lomb-Scargle periodogram for J0526.Upper panel: the TMTS L-band light curve over a 12 hr night on 18 December 2020.The magnitudes are presented as mean values ± 1σ .Middle panel: a 3000 s subset of the TMTS L-band light curve.The solid red line represents the best-fit sinusoidal model with a period of 10.3 min.Lower panel: the Lomb-Scargle periodogram (LSP) computed from the TMTS light curve.The vertical dashed line indicates the frequency corresponding to maximum power ( f max ).The purple dot-dashed line represents the confidence level of 0.1%, and the red arrow shows the frequency corresponding to the orbital period (i.e., f max /2).

Figure 2 .
Figure 2.Phase-folded RV curve and double-band light curves for J0526.Upper panel: the RV curve derived from Keck/LRIS and GTC/OSIRIS observations.The dot-dashed line is the best-fit sinusoidal model.Middle and lower panels: the gand r-band phase-folded light curves provided by LJT and ZTF.The purple solid lines represent the best-fit light-curve models obtained from the ellc package15 .Unequal maxima are caused by the relativistic Doppler beaming effect12,13 .Orbital phase φ = 0 represents the epoch of superior conjunction when the visible star is closest to the observer.The error bars represent 68% confidence intervals throughout this paper.

Figure 3 .
Figure 3. Dynamical spectra of J0526 from GTC/OSIRIS and Keck/LRIS observations.Upper panels: line profiles of Hα, Hβ , Hγ, and He I λ 4471 at all epochs of observations.Lower panels: spectral lines phased with the orbital period of 20.5 min.Color scales indicate the continuum-normalized flux and red solid lines represent the best-fit RV curve of the visible star of J0526.
-rich sequence Correlation for He-weak sequence

Figure 4 .
Figure 4. Relation between He abundance and effective temperature for hot subdwarfs.All atmospheric parameters of hot subdwarfs are taken from reference26 .The purple dashed line and green dotted-dash line represent the correlation for the He-rich sequence29 and the He-weak sequence30 , respectively.

Figure 7 .
Figure 7. Kiel diagram for hot subdwarfs, low-mass WDs, and J0526B.The atmospheric parameters of hot subdwarfs and (extremely) low-mass WDs are taken from the references26,50 .The blue dotted-dash lines are theoretical evolutionary tracks of helium-core WDs35 , while the shell flash loops are clipped for the clarity of image.The purple solid lines represent the evolutionary tracks of low-mass helium-core burning stars with a very thin hydrogen envelope of 10 −6 M ⊙ (see Methods).We overplot the zero-age helium main sequence (ZAHeMS; see Methods), zero-age extreme horizontal branch (ZAEHB), and terminal-age extreme horizontal branch (TAEHB)57 as black dashed lines.

HeFigure 8 .
Figure 8. Mass-radius diagram for main-sequence stars, WDs, hot subdwarfs, and J0526B.Masses, radii of MS stars, WDs, and hot subdwarfs are collected from C17 59 , B2050 , and S2263 , respectively.Six sdB stars inside ultracompact binaries having orbital periods below 100 min44,[64][65][66][67][68][69] are shown.For comparison, 18 WD components inside 11 ultracompact binaries64,66,[70][71][72] discovered by photometric surveys are also included.We overplot theoretical mass-radius relations for CO WDs, He WDs, MS stars, and HeMS stars.The degenerate and nondegenerate models are distinguished by solid and dashed lines, respectively.The purple area above the HeMS corresponds to the hot subdwarfs, i.e., the core-helium-burning (CHeB) stars with hydrogen-rich envelopes, in which thicker envelopes lead to larger stellar radii.A sequence for CHeB stars with an extremely thin hydrogen-rich envelope (M H,env = 10 −6 M ⊙ ) is highlighted using a pink dashed line.We present the constraints from surface gravity and SED, derived from current observations, accompanied with the mass-radius relations of an 80%/90% Roche-lobe-filling star inside a 20.5 min orbital-period binary using Paczyński's approximation73 (see Methods).The radii of Jupiter and 2MASS J0523−140361 (the smallest star among known MS stars) are labeled on the left of the figure.

Figure 9 .
Figure 9. Binary evolution models for extremely-short-orbital-period sdB binaries.Two models are differentiated owing to the different core masses of sdB stars.Mass transfers are expected to begin at around 14 and 17 min for M sd = 0.33 M ⊙ and M sd = 0.36 M ⊙ , respectively.The red arrows denote the direction of evolution.

Figure 10 .
Figure10.Color-magnitude diagram for hot subdwarfs, low-mass WDs, and J0526.The samples of hot subdwarfs and (extremely) low-mass WDs are taken from the references26,50 , and are then cross-matched with the Gaia DR3 database31 .The Gaia magnitudes and colors of models are computed using the bolometric correction tables from MESA Isochrones & Stellar Tracks (MIST)87,88 .

TFigure 11 .
Figure 11.GTC/OSIRIS spectra of J0526 with their best-fit TLUSTY/SYNSPEC synthetic spectra.Residuals between observations and models are listed in the lower panel.The main H/He absorption features are indicated by grey lines.