General relativistic orbital decay in a seven-minute-orbital-period eclipsing binary system

General relativity predicts that short orbital period binaries emit significant gravitational radiation, and the upcoming Laser Interferometer Space Antenna (LISA) is expected to detect tens of thousands of such systems; however, few have been identified, and only one is eclipsing--the double white dwarf binary SDSS J065133.338+284423.37, which has an orbital period of 12.75 minutes. Here, we report the discovery of an eclipsing double white dwarf binary system with an orbital period of only 6.91 minutes, ZTF J153932.16+502738.8. This system has an orbital period close to half that of SDSS J065133.338+284423.37 and an orbit so compact that the entire binary could fit within the diameter of the planet Saturn. The system exhibits a deep eclipse, and a double-lined spectroscopic nature. We observe rapid orbital decay, consistent with that expected from general relativity. ZTF J153932.16+502738.8 is a significant source of gravitational radiation close to the peak of LISA's sensitivity, and should be detected within the first week of LISA observations.

J1539+5027, phase-folded on the 6.91 minute orbital period. At phase 0, the lightcurve exhibits a deep primary eclipse, indicating that the hot primary star is producing most of the observed light. Outside of eclipse, there is a quasi-sinusoidal modulation because the primary star heavily irradiates one side of its companion. At phases ±0.5, the secondary eclipse occurs as the hot primary transits the irradiated face of its companion. b) The phase-folded ZTF g-band lightcurve of the object. We were able to discover the object because of its periodic behavior. c) A binned g lightcurve obtained with KPED, phase-folded on the orbital period. Error bars are 1σ intervals. 3 The short orbital period means that the two components must be dense objects-white dwarfs.
Because the primary eclipse is significantly deeper than the secondary eclipse, we can infer that one white dwarf (the primary) is hotter and more luminous than its companion (the secondary), as the detected flux is almost completely attenuated when the cooler object occults the hotter. By modelling the lightcurve (Methods), we can estimate the orbital inclination, i, the radius of the primary, R 1 , and the secondary, R 2 , relative to the semi-major axis of the orbit, a (Methods).
Because of ZTF J1539+5027's extremely short orbital period, general relativity predicts that it will undergo rapid orbital decay due to the emission of gravitational radiation 10 . With CHIMERA and KPED, we can precisely measure the time of eclipse, and use these eclipse times to measure a changing orbital period. If a system has a constant orbital period derivative, we expect the deviation of eclipse times, ∆t eclipse , (compared to those of a system with constant orbital period) to grow quadratically in time. Equation 1 ∆t eclipse (t − t 0 ) = 1 2ḟ illustrates this, where t 0 is the reference epoch, P (t 0 ) is the orbital period at the reference epoch, f (t 0 ),ḟ (t 0 ), etc, are the orbital frequency and its time derivatives at the reference epoch, and t − t 0 is the time since the reference epoch.
To measure the orbital velocities of the white dwarfs in the binary, we obtained phase-resolved spectroscopy using the Low Resolution Imaging Spectrometer (LRIS) 13  Using the spectroscopic observations, lightcurve modelling, and the orbital decay, we can constrain the masses of the white dwarfs in several ways: (1): With a mass-radius relation for the hot primary and constraints from lightcurve modelling, although this depends on parameters of white dwarf models, and only weakly constrains the mass of the secondary; (2): Using the spectroscopically measured radial velocity semi-amplitudes; however, this is challenging due to the blended absorption/emission lines, the latter depending on modelling irradiation effects and a substantial center of light correction; (3): With the chirp mass inferred from the measured orbital decay; however, this approach must account for potential tidal contributions.
Because each of these methods rely on different model dependent assumptions, we chose to polynomial fit to the deviation of the measured eclipse times as a function of time, compared to a system with constant orbital period. The consistency with a quadratic deviation demonstrates that the orbital period decreases with time. The orbital decay inferred is consistent with that expected from gravitational wave emission. The initial four timing epochs come from PTF/iPTF photometry, and the remainder were obtained with CHIMERA and KPED. b) The characteristic gravitational wave strain and frequency for ZTF J1539+5027 (red star in the plot). See Table 1   immediately after the primary eclipse. This is an ideal phase to isolate the photosphere of the hot primary, because it minimizes flux contributed by the irradiated face of the secondary. The smooth blue line is a fit of a white dwarf atmospheric model to this spectrum, yielding an effective photospheric temperature of T eff,1 = 48, 900 ± 900 K and a logarithm of surface gravity log(g) 1 = 7.75 ± 0.06 log(cm s −2 ) for the hot primary. b) and c) Two phase resolved spectra of the hydrogen n = 5 to n = 2 transition at ≈ 4340Å. The smooth red line is a double gaussian fit to the absorption and emission line used to measure the Doppler shifts of these features b) illustrates a phase in which the emission line associated with the cooler secondary is redshifted, while the absorption line associated with the primary is blueshifted. c) exhibits the opposing phase. 7 estimate physical parameters by combining these constraints ( Figure 4). We present the physical parameters resulting from this analysis in Table 1, alongside those of SDSS J0651+2844 4,5,14 . We conclude that the hot primary is a ≈ 0.6 M (likely carbon-oxygen) white dwarf, while the cool secondary is a ≈ 0.2 M (likely helium-core) white dwarf.
Because of their remarkably short orbital period, the white dwarfs experience significant tidal distortion. Tidal energy dissipation may heat and spin up the white dwarfs, in addition to increasing the orbital decay rate. Based on theoretical predictions 15 , tidal torques likely cause the spin periods of the white dwarfs to synchronize with the orbital period. We expect that tidal energy dissipation may increase the orbital decay rate by ≈ 7% relative to gravitational wave emission alone (Methods), though we cannot currently measure this due to uncertainty in the white dwarf masses. Future detection of the second derivative of the orbital decay rate,f , will enable a direct measurement of the tidal contribution to the orbital decay.
When LISA 2 first begins to operate (circa 2034), we estimate that it will detect ZTF J1539+5027 within a week. At the end of LISA's 4-year mission, we estimate the SNR will be 143 +13 −13 (Methods).
Not only does the source radiate high-strain gravitational waves, but LISA's sensitivity peaks at about 5 mHz, close to the 4.8 mHz gravitational wave frequency of this source, resulting in an exceptionally large SNR. This system will serve as a crucial "verification source" for LISA 2 , because its well-constrained inclination predicts the relative amplitude of the signal in the two gravitational wave polarizations 16 , and the precisely measured orbital decay already tightly constrains its chirp mass.  together with Roche potentials. The background white bins represent the constraint imposed by the measured spectroscopic radial velocity semi-amplitudes. The red line is a 50% contour level of the constraint inferred from applying a mass-radius relation to the hot primary white dwarf 30 , and combining it with the ratio of the primary's radius to the semi-major axis, R 1 a , inferred from lightcurve modelling. The magenta dashed lines are constraints imposed by the measured chirp mass; the upper dashed line assumes orbital decay purely due to general relativity, whereas the lower dashed line includes a 10 percent tidal contribution. The blue line is a 50% contour level representing a combination of all of these constraints (Table 1). 9 ZTF J1539+5027 poses challenges for models of binary evolution and the physics of accretion.
The spectroscopically measured temperature of the hot primary is T eff,1 = 48, 900 ± 900 K. The cooling age of such a hot white dwarf is ≈ 2.5 million years 17 , significantly shorter than the > 200 million year cooling age of the secondary 18 . Thus, some recent heating must have occurred. Tidal heating could increase the surface temperature of the primary to nearly 50, 000 K, though more realistic calculations (Methods) suggest temperatures closer to half this value. A more plausible explanation is that the heating is due to recent accretion, especially since the radius of the secondary indicates it is on the brink of Roche lobe overflow. Such accretion could heat the primary to its observed temperature for accretion rates ofṀ 10 −9 M yr −1 (Methods).
However, we see no evidence for active accretion. The only other known binary systems with orbital periods shorter than 10 minutes, V407 Vul (P ≈ 9.5 minutes) 19 and HM Cancri (P ≈ 5.4 minutes) 20 , were discovered because of periodic X-ray emission, thought to arise from a hot spot formed by direct impact accretion 21 . Unlike the other two sub-10-minute binaries, ZTF J1539+5027 exhibits no detectable X-ray flux. Based on an upper limits from observations by the XRT X-ray telescope on the Neil Gehrels Swift Observatory 22 , the EPIC-pn instrument on the XMM Newton Observatory 23 , and optical constraints, we have estimated an upper limit oḟ M < 2 × 10 −8 M yr −1 , contributing 10% of its energy to a hot spot. Reaching this upper limit requires fine-tuning the accretion hot spot temperature (Methods). Active accretion could still be occurring if the accretion energy is channeled primarily into heating the optically thick atmosphere of the accretor, and is then radiated at the 50, 000 K observed.
It is also possible that the accretion proceeds intermittently. One way of temporarily halting accretion could be a recent nova eruption on the surface of the primary. We expect that a 0.6 M white dwarf accreting from a companion at a rate ofṀ ∼ 10 −9 M yr −1 should experience 24 recurrent novae on timescales of ∼ 10 5 yr. However, while mass transfer may temporarily cease after each nova, calculations suggest that it is unlikely to catch the system in this short-lived phase (Methods). We conclude that the most likely scenarios are intermittent accretion with a mechanism allowing for phases of little-to-no mass transfer, or active accretion in which the accretion energy is radiated almost entirely in the ultraviolet and optical.
The orbit of ZTF J1539+5027 will continue to decay for ≈ 130, 000 years until it reaches a period of ≈5 minutes at which point the degenerate core of the secondary will begin to expand in response to mass loss, dramatically increasing the rate of mass transfer 25 . If the mass transfer is stable, which is likely based on the mass ratio 26 of q ∼ 1/3, the binary will evolve into an AM CVn system and the orbital period will increase. Alternatively, unstable mass transfer would result in a merger that could produce an R Cor Bor star 27 , or, less likely, a detonation of accreted helium on the primary could lead to a double-detonation that disrupts the primary 28 .

Methods 1 Summary of Observations
Extended Data Table 1 provides a summary of all observations used in our analysis.

Period Finding
We identified ZTF J1539+5027 by using the conditional entropy 31

Lightcurve Modelling
To model the lightcurve, we used data from three nights of CHIMERA g observations (July 5-7 2018), with a total of 12, 999 individual 3-s exposures.
To fit the CHIMERA data, we used the ellc package 32 to model the lightcurve and fit for the ratio of the radii to the semi-major axis, surface brightness ratio, J, of the unheated face of the secondary compared to the hot primary, and the mid-eclipse time of the primary eclipse, t 0 . We adopted a linear limb darkening model, with limb darkening coefficients for the primary (ldc 1 ) and secondary (ldc 2 ). We treated the primary as a spherical object, but invoked a Roche approximation for the secondary. We also included a gravity darkening coefficient for the secondary (gdc 2 ) and a single free heating parameter (heat 2 ) to attempt to fit the sinusoidal variation due to the irradiation of the secondary, which acts as an albedo of the secondary, but in this system must be larger than 1 to achieve a good fit because a significant amount of incident ultraviolet light is reprocessed into optical wavelengths.
We allowed the limb and gravity darkening coefficients to vary as free parameters in the fit, using uniform priors for each with values of ldc 1 = 0.15 ± 0.15, ldc 2 = 0.4 ± 0.2, and gdc 2 = 0.6 ± 0.1, based on extrapolations of existing models 33,34 .
We performed the final fit using the period derived from the quadratic fit to the timing epochs (Table 1). We left r 1 , r 2 , i, J, q, t 0 , ldc 1 , ldc 2 , gdc 2 , and heat 2 as free parameters. Extended Data Figure 1 illustrates the corner plots from this fit, but excludes t 0 , which did not exhibit significant covariance with any other parameter. We fixed the eccentricity to 0, because measuring this quantity depends on the shape of the secondary eclipse, which in turn depends on the poorly understood irradiation of the secondary. We ruled out the possibility of significant eccentricity, because we failed to detect any sign of apsidal precession in the eclipse time measurements, and furthermore, we tried fitting for this using the lightcurve modelling, and found a value consistent with 0. Because we lack a good physical model for the irradiation of the secondary, we also do not account for Doppler beaming, which we expect to alter the shape of the lightcurve at the few percent level, peaking at phases between the eclipses, where the irradiation dominates the behavior of the lightcurve.

Orbital Decay
To measure the orbital decay, we independently fit each night of KPED and CHIMERA data for mid-eclipse times using the lightcurve modelling performed with ellc. We convert all timestamps to Barycentric Julian Dates (BJD) in Barycentric Dynamical Time (TDB) to achieve the required timing precision. Because the eclipse time is not strongly covariant with model dependent parameters such as gravity and limb darkening, we omitted these from the fits to reduce complexity (but still fitting for all other parameters described in the previous section). We also extracted photometry from archival PTF/iPTF data, and because of significant smearing due to the exposure time of 60-s (particularly for the primary eclipse), we extracted timing epochs by performing a least squares fit of a sinusoid to this data (Extended Data Figure 2).
After obtaining the timing epochs, we measured the deviation of each eclipse time since the start of observations relative to a model with a constant orbital period ( Figure 2). This was non-trivial, because the eclipse time had drifted by multiple orbits since the PTF epochs. However, upon allowing for an integer number of orbits to have passed, only one solution yielded a significant fit when using a quadratic model (corresponding to 5 orbits since the initial epoch). times to estimateṖ (this fit is the red curve down in (Figure 3a), which resulted in an adjusted R 2 = 0.999995. This quadratic fit is independent of any model assumptions about decay due to tidal contribution, general relativity, or any other mechanism; it only indicates that the orbital period is decreasing with with an approximately constant rate as a function of time over the course of the observations.
In addition to extractingṖ , we verified our measurement by fitting a quadratic to only mid-eclipse times extracted from CHIMERA and KPED data (Extended Data Figure 3 WithṖ measured, we assessed how much contribution there was due to "secular acceleration" 35 and differential Galactic acceleration 36 . For the former we used the proper motion as measured from the second data release of Gaia 37 , µ = 5.1 ± 2.2 mas yr −1 . This caused an apparent excesṡ P Shk of (8.3 ± 3.6) × 10 −7 (d/1 kpc) ms yr −1 where d is the distance. This term must be subtracted from the measured value to obtain the intrinsicṖ . Similarly, differential Galactic acceleration leads to an apparentṖ DGR of −3.5 × 10 −7 ms yr −1 at a distance of 2.4 kpc, computed using the All error bars on the timing epochs are 1 sigma.

5 Spectroscopic Analysis
To perform the spectroscopic analysis, first we coadded 317 individual spectra (all taken with a 52s exposure) into 12 phase-bins, using the emphemeris for the mid-eclipse time of the primary eclipse to define phase 0. We then fit stellar atmospheric models 39 and obtained measurements of the logarithm of the surface gravity of the primary, log(g) 1 , and its effective photospheric temperature, T eff,1 , using the spectrum immediately after the phase of the primary eclipse (phase 0.0833), to minimize the flux contributed by the irradiated face of the secondary (Figure 3a).
In order to measure radial-velocity semi-amplitudes of the objects, K 1 and K 2 , we used a double Gaussian model to fit the line associated with the hydrogen n = 5 to n = 2 transition in each phase-binned spectrum to extract the Doppler shifts of the emission and absorption lines.
To derive overall radial velocities, we adopted an out-of-phase sine wave model for the absorption and emission features v 1 (j) = K 1measured sin 2πj 12 + A, v 2 (j) = K 2measured sin 2πj 12 where j encodes the index of the phase-bin and A and B are the systemic velocities of each white dwarf, K 1measured and K 2measured are the observed velocity semi-amplitudes, and v 1 (j) and v 2 (j) are the Doppler shifted velocities associated with each phase-bin (Extended Data Figure 4). We required the two systemic velocities to satisfy A > B when sampling, because the gravitational redshift of the primary should exceed that of the secondary because of the primary's larger surface gravity. We used a Kernel Density Estimator (KDE) applied to the posterior distributions for the measured velocity semi-amplitudes K 1measured and K 2measured to assign the probabilities when 24 sampling.
In order to derive masses from these spectroscopic fits, we applied corrections to the measured velocity semi-amplitudes. First, we applied a smearing correction of 20% to both K 1measured and K 2measured , due to our phase binned spectra each being co-additions of spectra taken over a third of the orbital phase (though each individual spectrum was taken over only an eight of the orbital phase, broader coadditions were necessary to reach sufficient SNR for measuring radial velocities).
We used a prior of C = 0.68 ± 0.04 for the center of light correction term, based on the mass ratio estimate 40 . This modified the velocity amplitude of the secondary by: where R 2 /a was taken from the lightcurve fitting, and the mass ratio was defined as q = K 1 /K 2 (thus requiring us to solve this expression to isolate K 2 ). We applied the center-of-light correction to the secondary, because one face of it is heavily irradiated (meaning that the emission lines arise from a location different than the center of mass).

Mass and Radius Analysis
There are three main constraints that contribute to the overall mass and radius estimates. Because of the low SNR of the spectra, these fits have large uncertainties (especially in the case of the primary, with its shallow and broad absorption lines). This is reflected in the broad distribution of possible masses associated with the spectroscopic constraint illustrated in Figure 4. All error bars are 68% CIs.
where G is the gravitational constant. Because we have measured both K 1 and K 2 , as well as P and i, we can write two such equations, and derive constraints on the two masses ( Figure 4).

White Dwarf Model Constraint:
We used a mass-radius relation from models for a carbon-oxygen white dwarf 30 , which depends on mass, metallicity, hydrogen fractional mass in the atmosphere, and radius, together with the spectroscopic measurements of T eff,1 and log(g) 1 , to derive constraints on these properties. In these fits, we marginalized over metalicity and hydrogen mass fraction. We used the measured ratio of the radius of the primary to the semi-major axis R 1 a (from the lightcurve modelling), and combined this with Kepler's law, and the white dwarf model constraints 30 , to constrain the system parameters. This also weakly constrains the mass of the secondary, which enters into Kepler's law as contributing to the total mass of the system. , which is be related to the orbital decay rate by Equation 5: where c is the speed of light. Because the chirp mass inferred fromḟ assumes the decay is caused purely by general relativity, we also computed a chirp mass assuming a 10 percent tidal contribution toḟ as a lower bound on M c .

Combined Fit:
To combine all of the constraints described above, we created a KDE based on the M 1 and M 2 estimates from the spectroscopic, white dwarf model, and chirp mass constraints.
The posterior distribution of this combined analysis yielded the values for the masses reported in Table 1. The radii were computed using r 1 and r 2 from the lightcurve modelling, and the semi-major axis a determined with the masses and the orbital period.

Distance Estimate
Because no reliable parallax measurement exists for ZTF J1539+5027, we instead used its bolometric luminosity to estimate the distance. First, we measured the apparent magnitude of the hot primary without contribution from the irradiated face of the secondary by performing a least squares fit of a sinusoid to the ZTF-g lightcurve, omitting data from the eclipse, and measuring the minimum of this sinusoid. This yields an apparent ZTF-g band apparent magnitude of 20.38 ± 0.05.
Next, we used the atmospherically determined temperature of the primary, T eff,1 ≈ 48, 900 K, to infer the absolute g luminosity of the 0.6 M DA white dwarf 17 . The spectroscopic temperature puts ZTF J1539+5027 approximately between two temperatures in the models-one corresponding to T eff,1 = 45, 000 K and the other to T eff,1 = 55, 000 K. These corresponded to g absolute magnitudes of 8.71 and 8.35, respectively, and we adopted these values as the lower and upper bounds on the object's absolute luminosity, and assume a uniform distribution of possible absolute luminosities between these values. Using this to solve for the distance of the object, we estimated a distance of d = 2.34±0.14 kpc, though we emphasize that the error bars derived using atmospherically determined quantities tend to be underestimated. We also incorporated a uniform distribution of 28 E(B − V ) in the range of 0 to 0.04 to account for the effects of extinction at these coordinates 41 .

Gravitational Wave Strain
The expression for the characteristic strain 42 used in Figure 2b (including the value plotted for ZTF J1539+5027) is: where d is the distance to the object, T obs is the integrated observation time of the LISA mission.
Though it is the conventional quantity used to construct such diagrams, the characteristic strain does not capture any information about source inclination, detector response, etc.
In order to compute the signal to noise 42 for LISA, we directly invoked the signal amplitude at the detector A = |F + | 2 |h + | 2 + |F × | 2 |h × | 2 , where h + and h × are the two gravitational wave polarization amplitudes, and F + and F × are the LISA detector response patterns corresponding to these polarizations. The h + polarization amplitude includes a factor of (1 + cos 2 (i)) and the h × polarization a factor of 2 cos (i), meaning that systems like ZTF J1539+5027 with an inclination close to 90 degrees exhibit a gravitational wave signal up to a factor of √ 8 smaller than an equivalent face-on system with an inclination close to 0 degrees in situations where F + ≈ F × .
In Table 1, we include an estimate for both the SNR of ZTF J1539+5027, and SDSS J0651+2844 computed using the same technique.

Tidal Effects
Tidal Contribution to Orbital Decay Tidal dissipation can transfer orbital energy into rotational and thermal energy within the stars. The former effect can cause the orbit to decay slightly faster than gravitational radiation alone, while the latter effect can increase the surface temperatures of the stars. Studies of tidal synchronization and heating predict that tidal energy dissipation scales more strongly with orbital period than gravitational radiation 15,43 . As the orbit decays and the white dwarfs draw nearer to each other, tides begin to act on a shorter timescale than orbital decay, spinning up the stars toward synchronous rotation. The "critical" orbital period, P c , below which tidal spin-up can occur faster than orbital decay, lies in the range P c ∼ 45 − 130 minutes 43 , depending on the white dwarf mass and age. In any case, at an orbital period of only 6.91 minutes, we expect the stars in ZTF J1539+5027 to be spinning synchronously with the orbit. In this regime, the rate of tidal energy transfer from the orbit to the stellar interiors iṡ where I = I 1 +I 2 is the combined moment of inertia of the two stars, andṖ GW is the orbital period decay caused by gravitational waves.
Comparing Equation 7 with the orbital energy lost to gravitational wave emission, the tidal contributionṖ tide to the total orbital decay rateṖ is given bẏ The moment of inertia of star 1 is I 1 = κ 1 M 1 R 2 1 (and similarly for star 2), where κ 1 is a dimensionless 30 constant determined by the internal structure of the white dwarf. Equation 8 can also be written aṡ From the white dwarf models (Section 10), we find κ 1 0.14 and κ 2 0.11. Using M 1 0.61M , M 2 0.21M , R 1 /a 0.14, R 2 /a 0.28, we findṖ tide /Ṗ GW 0.067. Thus, we expect the orbit to decay several percent faster than gravitational radiation alone, provided that the stars are tidally locked. There is an additional effect of orbital decay caused by orbital energy used to raise a tidal bulge in each star 16, 44 , but we find it to be more than an order of magnitude smaller than Equation 9 .
Tidal Heating The tidal energy dissipation within the white dwarfs is partitioned between kinetic energy used to spin up the white dwarfs and heat that can diffuse out and increase their surface temperatures. For an aligned, sub-synchronous and rigidly rotating white dwarf, the ratio of tidal heating to tidal energy dissipation isĖ heat /Ė tide = 1 − P/P spin , where P spin is the spin period of the white dwarf 43 . At orbital periods below the critical value P c , the spin period, P spin , decreases to approach the orbital period. Calculations suggest that the tidal heating rate in this regime is expected to beĖ heat ∼Ė tide (P/P c ) 15 , such that the tidal heating rate is much smaller than the energy dissipation rate, because most of the energy is converted to rotational energy.
To estimate an upper limit on the surface temperature of each white dwarf that can be obtained from tidal heating alone, we assumeĖ heat ∼Ė tide ≈ 6π 2 I/(P 2 t GW ), where t GW = 3 2 P |Ṗ | is the gravitational wave timescale. If the tidal heat is instantaneously reradiated, this corresponds 31 to a surface temperature where σ B is the Stefan-Boltzmann constant, and M is the mass of the white dwarf. Note that this temperature depends only weakly on the white dwarf mass and moment of inertia, and is independent of the white dwarf radius. Using the same values as above and the measured value t GW ≈ 830, 000 yr, we find an upper limit of T tide ≈ 44, 000 K for the primary white dwarf.
However, accounting for the suppression factor P/P c ∼ 6.9/60, such thatĖ heat ∼Ė tide (P/P c ), yields a more realistic temperature T tide ≈ 25, 000 K. Hence, while tidal heating may be able to heat the primary to temperatures near that observed, our best estimate suggests substantially lower temperatures. Applying Equation 10 to the secondary white dwarf predicts an upper limit due to tidal heating of T tide ≈ 33, 000 K, but our best estimate is T tide ≈ 19, 000 K. A measured nightside temperature near this value would be consistent with tidal heat powering the nightside flux of the secondary. However, the value of the surface brightness ratio in g suggests a secondary temperature of T eff,2 < 10, 000 K.

Binary Models
To understand the evolutionary state of ZTF J1539+5027, we have constructed several binary models using the MESA stellar evolution code 45 Figure 5. We initialize these runs at orbital periods of 1 hour to mimic the end of a common envelope event that formed the tight binary. The helium white dwarf is evolved simultaneously with the orbit, with angular momentum losses due to gravitational waves and fully non-conservative mass transfer.
Extended Data Figure 5 shows a plot of the mass loss rate from the secondary as function of orbital period. The secondary overflows its Roche lobe and begins mass transfer at orbital periods ranging from P ∼ 6.5 − 13 minutes, and the expected mass loss rates at P = 7 minutes are typicallyṀ ∼ 3 × 10 −9 M yr −1 (ranging from 0 up to 10 −8 M yr −1 , depending on the white dwarf masses and hydrogen envelope mass). During the initial phase of slow mass transfer, the secondary loses its non-degenerate hydrogen envelope as the Roche lobe contracts inward. The mass accreting onto the primary can greatly heat it as gravitational potential energy is converted to heat. The energy released by accretion isĖ accrete ∼ GM 1Ṁ /R 1 , with order unity corrections due to its non-zero kinetic and gravitational energy when it is lost from the secondary. We do not evolve the primary (accretor), but we may crudely estimate its temperature by assuming that the accretion energy is uniformly radiated as a blackbody over its surface. The corresponding accretion temperature (assuming 100% efficiency) is The bottom panel of Extended Data Figure 5 shows that once mass transfer begins, it can easily increase the primary's temperature to T 1 50, 000 K. Systems begin at large orbital period and move towards smaller period due to gravitational radiation, and in some cases they move back out due to stable mass transfer. Except for high-mass donors with thin hydrogen envelopes, mass transfer is expected to begin at orbital periods longer than 7 minutes. Bottom: Corresponding accretion temperature from Equation 11.
Although accretion could explain the high temperature of the primary, we have not detected any evidence of ongoing mass-transfer. We have constrained the possibility of active accretion contributing luminosity to an accretion hot spot using both optical and X-ray data. The upper limit we inferred isṀ < 2 × 10 −8 M yr −1 , and this is only possible in a very narrow hot spot temperature range (Extended Data Figure 6). For both the optical and X-ray constraint we have assumed that only 10% of accretion energy is being converted to the luminosity of the hot spot. This is based on a model where we assume 90% of the accretion energy is being deposited into heating the optically thick photosphere of the white dwarf, while only 10% is contributing to luminosity immediatly re-radiated in the form of a hot spot. We estimated the upper limit by computing the X-ray flux using NASA's WebPIMMS mission count rate simulator, using a 3-sigma upper limit on the background count rate determined from the X-ray images. We assumed a hydrogen column density 46 of nH = 1.5 × 10 20 cm −2 . We used the distance of d = 2.4 kpc to convert the upper limit on the unabsorbed X-ray flux to an X-ray luminosity. For the optical constraint, we computed an upper limit on the optical luminosity of the hot spot as 10% the luminosity of the photosphere of the hot primary, based on an absence of emission lines in our coadded spectra, which have an SNR of approximately 10. We calculated the upper limit onṀ at various temperatures by integrating a Planck function at the corresponding temperature over the instrument passbands, and then computed the maximum bolometric luminosity, L accretion , an emitting region at this temperature could have and still be consistent with the non-detections in these passbands.
We then determined the corresponding upper limit onṀ by equating L accretion = 0.1GM 1Ṁ /R 1 .
We obtained the XMM EPIC-pn data from the public XMM-Newton science archive (observation ID: 0800971501 PI: Pratt, Gabriel). The SWIFT observation was obtained with our own program (observation ID: 00010787001, 00010787002 PI: Kulkarni).

Novae
The high temperature of the primary may plausibly be explained if the white dwarf is cooling after having recently undergone a nova outburst, caused by accretion of hydrogen from the secondary or a tidally induced nova 47 . The nova likely ejects an amount of mass comparable to the hydrogen shell mass which must be accreted, which for a 0.6M white dwarf is M H ∼ 10 −4 M 24 .
Following the nova, the orbit widens slightly, and the system is brought out of contact such that mass transfer from the secondary ceases. The change in semi-major axis following the loss of the nova shell is ∆a/a ∼ M H /(M 1 + M 2 ). The length of time the binary is detached before gravitational wave emission brings the system back into contact is Using M H = 10 −4 M , we estimate t detach ∼ 100 yr. This can be compared to the time spent accreting mass between subsequent novae.
As the non-degenerate hydrogen envelope of the low-mass secondary is stripped off (see discussion in 25 ), the approximate mass transfer rate is expected to beṀ 10 −8 M yr −1 (Extended Data in terms of the mass accretion rate contributing to the accretion luminosity of a hypothetical hot spot. The solid red curve illustrates the constraint imposed by the XMM EPIC-pn X-ray non-detection, which rules out significant mass transfer contributing to a hot spot with temperatures greater than ≈ 150, 000 K, while the green dotted line illustrates a weaker upper limit imposed by the non-detection in a SWIFT XRT observation. We constructed the dashed blue curve, which represents the optical constraint, by requiring that any accretion luminosity originating from a hot spot should contribute < 10% to the luminosity in the band ranging from 320 to 540 nm, as we know from the optical spectrum ( Figure 3) that this light is dominated by the ≈ 50, 000 K photosphere of the hot primary, and also see no signature of a hot spot in the CHIMERA lightcurve ( Figure 1). We chose the threshold of < 10%, because given the SNR of the spectra, we expect we should be able to detect optically thin emission with an amplitude at the 10% level. Other white dwarfs with such a hot spot (such as HM Cancri) exhibit such emission, particularly in lines associated with ionized helium. spent in a detached state relative to an accreting state is t detach /t accrete < 10 −2 . Hence, while it is possible that the system is in a detached state following a nova caused by mass transfer, the chances of catching the system in this state are small. To help rule out the possibility, we used the WASP instrument on the Hale telescope to obtain a deep H-α image of the field and found no evidence for a remnant nova shell; however, this analysis was limited by the lack of an off-band image.

Population Implications
From 48 , the merger rate of He+CO WDs in the Milky Way is roughly 0.003 yr −1 . This number is reached from both observational and population synthesis arguments. The number of systems with decay time equal or less than the ∼210 kyr decay time of ZTF J1539+5027 is thus ∼630.
Out to the distance of 2.3 kpc, given a local surface density of 68 M pc −2 from 49 , the stellar mass is ∼ 10 9 M , roughly 2% of the total disk mass of ∼ 5 × 10 9 M . We thus expect to find ∼13 binaries with a similar distance and merging timescale as ZTF J1539+5027. The fraction of eclipsing systems is roughly R/a ∼ 0.25 for our measured parameters, hence we may expect ∼3 eclipsing systems like ZTF J1539+5027. ZTF can detect such systems in most of the volume out to its distance, as long as they are as bright as this system. We may be missing slightly longer period systems that are dimmer because they have not yet started mass transfer. We comment that the estimate from 48 found that many double WDs must be born at short orbital periods in order to explain the abundance of short period systems relative to longer period systems, and ZTF J1539+5027 may support that conclusion. 38 Extended Data Upon request, the first author will provide reduced photometric and spectroscopic data, and available ZTF data for the object. We have included the eclipse time data used to construct the orbital decay diagram in Figures 2a, Extended Data Figure 2, and Extended Data Figure 3. The X-ray observations are already in the public domain, and their observation IDs have been supplied in the text. The proprietary period for the spectroscopic data will expire at the start of 2020, at which point this data will also be public and readily accessible.

Code Availability
Upon request, the first author will provide code (primarily in python) used to analyze the observations and data such as posterior distributions used to produce the figures in the text (MATLAB was used to generate most of the figures).