Orbital Decay in a 20 Minute Orbital Period Detached Binary with a Hydrogen-poor Low-mass White Dwarf

We discovered PTF J0533+0209 during a broad search for post-common-envelope binaries (Ivanova et al. 2013). In order to target hot young remnants of the common-envelope (CE) phase, we used the Pan-STARRS DR1 (Chambers et al. 2016) photometric survey to target blue objects, selecting all objects with a color satisfying the condition (g–r) < 0 (see Table 1—note that the values in the table are observed apparent magnitudes, and have not been de-reddened).

After imposing the photometric color cut, we cross-matched the resulting sample with the archival Palomar Transient Factory (PTF; Law et al. 2009) photometric database. The PTF (and the Intermediate PTF) was a northern-sky synoptic survey using the 48 inch Samuel Oschin Telescope at Palomar Observatory. The survey was conducted in PTF r and g bands with a typical exposure time of 60 s, resulting in limiting magnitudes of approximately 21 and 20, respectively. We restricted ourselves to a search in the r band because it is the most heavily sampled of the PTF bands. Additionally, we required a minimum of 20 epochs in the lightcurves. This cross match yielded ≈180,000 sources. We discovered PTF J0533+0209 via a period search. We detected the object with high significance because it falls in the most heavily sampled field in all of PTF, located inside the constellation of Orion, which was observed over 5000 times, primarily in the first two years of PTF. This field did not use the typical PTF exposure time of 60 s, but instead used 30 s exposures. These observations were taken in two brief intervals, one consisting of ≈3000 observations in 2009 December–2010 January, and another in 2010 December with ≈2000 observations (Table 2).

We used the high-speed photometer on the 200 inch Hale telescope at Palomar Observatory, CHIMERA (Harding et al. 2016), to obtain a well-sampled follow-up lightcurve of the object. The instrument consists of a pair of electron-multiplying charge-coupled devices (EM CCDs) and has a dichroic and dual channels, allowing us to obtain simultaneous observations in g' and i'. We used 10 s exposures for all observations operating with the conventional amplifier, and these lightcurves served as the basis for our analysis of the ellipsoidal modulation and relativistic Doppler beaming exhibited by the visible component in the system (Figure 1).

Zoom In Zoom Out Reset image size

Additionally, we used the Kitt Peak Electron Multiplying CCD demonstrator (KPED; Coughlin et al. 2019), a high-speed EM CCD photometer mounted on Kitt Peak National Observatory's 84 inch telescope, to obtain additional observations to use as timing epochs (in order to measure the orbital decay). We obtained observations in g'. Unlike our CHIMERA observations, we operated KPED with electron-multiplying gain enabled, effectively eliminating read noise. The detector acquired images at a rate of 8 Hz, which we then stacked to 10 s coadditions in order to match the exposure time used with CHIMERA. Table 2 gives a summary of these observations.

All photometric data were reduced using a custom pipeline.

On 2017 October 25, we observed the object on the Hale Telescope with four consecutive 5 minute exposures using the Double Spectrograph (Oke & Gunn 1982). While these spectra have low signal to noise, they nonetheless reveal the presence of large Doppler shifts on the order of several hundred km s⁻¹. We reduced this data using the pyraf-dbsp reduction pipeline (Bellm & Sesar 2016).

We then obtained an additional 58 two-minute exposures using the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on the 10 m W. M. Keck I Telescope on Maunakea, with the 400/8500 grism and a 2 × 2 binning on the blue arm. Eighteen of these exposures were obtained on 2017 November 15 and the remaining 40 on 2018 March 19. These exposures were sufficiently short to allow us to create a time-resolved radial velocity curve (Figure 1) and also assemble a coadded spectrum for atmospheric fitting with minimal broadening of the features from the Doppler shifts (Figure 2). In individual exposures, we averaged a signal-to-noise ratio (S/N) of 4–5, with a wavelength coverage from 3400–5600 Å, and a resolution of approximately $\tfrac{\lambda }{{\rm{\Delta }}\lambda }=700$. In order to ensure wavelength stability, we took a HeNeArCdZn arc at the telescope position of the object. For both nights, we used five consecutive dome flats and five bias frames to perform calibrations, and reduced the data using the lpipe pipeline (Perley 2019).

Zoom In Zoom Out Reset image size

We discovered the PTF J0533+0209 by applying the multiharmonic analysis of variance (MHAOV) period finding routine to the sample of color-selected lightcurves (Schwarzenberg-Czerny 1996), searching periods from 7.2 minutes to 1000 days. PTF J0533+0209 exhibits strong ellipsoidal modulation in its PTF lightcurve (Figure 3), and thus exhibits significant power at half the orbital period in its power spectrum (Figure 4). We used MHAOV because of its superior sensitivity to sharp non-sinusoidal periodic features such as eclipses while still remaining sensitive to sinusoidal lightcurves such as that of PTF J0533+0209. We used the implementation of MHAOV available in the Vartools package (Hartman & Bakos 2016). In the initial search for periodic objects, we use heliocentric Julian dates as timestamps. However, all timestamps in this publication have been corrected to mid-exposure barycentric Julian dates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We measured the Doppler shift of the absorption lines in each time-resolved spectrum by fitting the centroids of the hydrogen Balmer lines of H_β and H_γ, as well as the neutral helium lines at wavelengths of 4471, 4713, and 4921 Å. All of these lines belong to PTF J0533+0209B. We derived the velocity semi-amplitude of PTF J0533+0209B, K_B, by performing a weighted least-squares fit of a sinusoid to the combined RV data, with a fixed ω corresponding to the orbital period derived from the PTF data. We derive a velocity semi-amplitude of K_B = 618.7 ± 6.9 km s⁻¹ and a systemic velocity of γ = 76.0 ± 4.3 km s⁻¹, where we have obtained the 1σ error bars from the covariance matrix of the least-squares fit to the time-resolved radial velocity measurements. The mass function of PTF J0533+0209A is given by

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{A}^{3}{\sin }^{3}(i)}{{\left({M}_{A}+{M}_{B}\right)}^{2}}=\displaystyle \frac{{{PK}}_{B}^{3}}{2\pi G}\end{eqnarray} \tag{ 1 }$$

where M_A is the mass of PTF J0533+0209A, M_B is the mass of PTF J0533+0209B, i is the orbital inclination, and G is the gravitational constant.

Thus, using Equation (1), we can constrain the relationship between the physical parameters of interest in the left-hand side of the equation with the measured P and K_B.

We calculated a grid of DBA (helium-dominated with traces of hydrogen) white dwarf atmosphere models to fit the spectra of PTF J0533+0209B (we found no evidence of any spectroscopic features associated with PTF J0533+0209A). The grid spans the effective temperature range of 11,000 ≤ T_eff ≤ 30,000 K in steps of 1000 K, the surface gravity range of 5.5 ≤ log g ≤ 9.0 in steps of 0.5 dex, and the hydrogen-to-helium number density ratio of −1.0 ≤ log H/He ≤ −10.0 in steps of 0.5 dex. A detailed description of the model atmosphere code can be found in Bergeron et al. (2011).

Analogous to the well-established spectroscopic method used for DA white dwarfs, we compared the absorption line profiles of the continuum-normalized spectrum of PTF J0533+0209B with our newly computed models. In our fit we included an E(B–V) reddening of 0.13 and assumed an extinction constant A_V = 3.1 (Green et al. 2019).

The best-fit model solution corresponds to T_eff = 20,000 ± 800 K, log g = 6.3 ± 0.1, and log H/He = −2.7 ± 0.1 (Figure 5). For limb- and gravity-darkening coefficients, an additional grid of atmosphere models was computed for log H/He = −2.56, 19,250 ≤ T_eff ≤ 20,500 K in steps of 250 K and 5.4 ≤ log g ≤ 7.2 in steps of 0.6 dex. For each model we calculated the specific intensity at 20 different angles.

Zoom In Zoom Out Reset image size

Because this binary exhibits significant ellipsoidal modulation, we can impose constraints on the system parameters by modeling its lightcurve. We see ellipsoidal variations in the lightcurve because PTF J0533+0209B is being tidally deformed, and over the course of the orbit, we see different geometric cross sections of the object. Equation (2) gives the fractional flux variation due to ellipsoidal modulation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{F}_{\mathrm{ellipsoidal}}}{F}=0.15\displaystyle \frac{(15+u)(1+\tau )}{3-u}{\left(\displaystyle \frac{R}{a}\right)}^{3}q{\sin }^{2}i,\end{eqnarray} \tag{ 2 }$$

where the fractional flux variation due to ellipsoidal variations, $\tfrac{{\rm{\Delta }}{F}_{\mathrm{ellipsoidal}}}{F}$, depends on the ratio of the radius of the tidally deformed object to the semimajor axis, R/a, the tidally deformed component's linear limb-darkening coefficient, u, its gravity-darkening coefficient, τ, the mass ratio of the binary, q, and the inclination, i (Bloemen et al. 2012).

Because this object is strongly deformed, we must also account for an additional effect known as gravity darkening, which visibly manifests itself as an asymmetry in the two flux minima of the ellipsoidal modulation. This is a result of the flux emitted on the surface of the deformed object near the first Langrange point (L1) being less than that emitted from the surface near the L2 point. Von Zeipel's theorem (von Zeipel 1924) states that the flux emitted off of the surface element of a star is proportional to the local gravitational field. Although the material at the surface of a tidally deformed star falls on a gravitational equipotential surface ϕ = const, the gradient of the gravitational potential along this surface is not necessarily constant or symmetric, i.e., $-{\rm{\nabla }}\phi =g\ne \mathrm{const}$, and this is manifested by an asymmetry between the flux emitted by surface elements at L1 and L2.

We used the ellc package to model the lightcurve (Maxted 2016). We modeled both the g' and i' lightcurves, adopting a linear limb-darkening model, and fixing the linear limb-darkening coefficient of the secondary u₂ to a value of 0.25 for the g' lightcurve, and to 0.18 for the i' lightcurve. We also fixed the passband-dependent gravity-darkening coefficient of the secondary, τ₂, to 0.25 for the g' lightcurve, and 0.19 for the i' lightcurve. These coefficients were computed using the procedure outlined in Claret (2017), adjusted for DBA atmospheres (A. Claret et al. 2019, in preparation). We computed and fixed the Doppler-beaming factors of PTF J0533+0209B to 1.94 for g' and 1.54 for i' (Loeb & Gaudi 2003) using the atmospherically derived temperature of T_eff = 20,000 K. The free parameters in the model were the ratio of the radii to the semimajor axis for the two components, ${R}_{A}/a$ and ${R}_{B}/a$, the inclination, i, the time of minimum light, t₀, the semimajor axis, a, and the mass ratio, q = M_B/M_A. Although we see no sign of luminosity from PTF J0533+0209A in the spectrum or the spectral energy distribution, as an initial test, we conducted an iteration of modeling with the surface brightness ratio, J, and albedo of PTF J0533+0209B, heat_B, as free parameters to investigate whether PTF J0533+0209B might be irradiated by its unseen companion. We found the solution did not converge to any particular value of these parameters for the g' or i' lightcurves, and that the parameters exhibited no covariance with any other parameters. We performed the same exercise considering the possibility of PTF J0533+0209A being irradiated by PTF J0533+0209B, and found a similar result. Based on this, for the final iteration of modeling, we assumed no luminosity contribution from PTF J0533+0209A, and omitted accounting for any kind of irradiation effect. The faint GALEX NUV apparent magnitude (see Table 1) also suggests it is unlikely that there is an unseen hot companion.

Because of the short orbital period and non-interacting nature of this system, we expect a possible measurable orbital decay due to the emission of gravitational radiation (Taylor et al. 1979). One can determine the orbital frequency derivative, $\dot{f}$, by measuring a change in the phase of the lightcurve over time, Δt_ellipsoidal, as illustrated in Equation (3):

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{\mathrm{ellipsoidal}}(t-{t}_{0})=\left(\displaystyle \frac{1}{2}\dot{f}({t}_{0}){\left(t-{t}_{0}\right)}^{2}+\ldots \right)P({t}_{0}),\end{eqnarray} \tag{ 3 }$$

where t₀ is the reference epoch, t − t₀ is the time since the reference epoch, $f({t}_{0}),\dot{f}({t}_{0})$, etc., are the orbital frequency and its derivatives at the reference epoch, and $P({t}_{0})=\tfrac{1}{f({t}_{0})}$ is the orbital period at the reference epoch.

In order to extract the timing epochs from PTF, KPED, and CHIMERA data, we used a least-squares fit of a sinusoid to the data and measured the phase of this sinusoid. We then used a least-squares fit of a quadratic to these timestamps as a function of epoch (Figure 6) to measure the orbital decay rate, $\dot{P}$. We report the derived ephemeris in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Table of Physical Parameters

MA	${0.652}_{-0.040}^{+0.037}\,{M}_{\odot }$
MB	${0.167}_{-0.030}^{+0.030}\,{M}_{\odot }$
i	${72.8}_{-1.4}^{+0.8}\,\mathrm{degrees}$
RB	${0.057}_{-0.004}^{+0.004}\,{R}_{\odot }$
TB	${20000}_{-800}^{+800}\,{\rm{K}}$
$\mathrm{log}{({\rm{g}})}_{{\text{}}B}$	${6.3}_{-0.1}^{+0.1}$
$\mathrm{log}\,{\rm{H}}/{\mathrm{He}}_{{\text{}}B}$	$-{2.7}_{-0.1}^{+0.1}$
T0	${2458145.096848}_{-0.000046}^{+0.000046}\,{\mathrm{BJD}}_{\mathrm{TDB}}$
P(T0)	${1233.97298}_{-0.00017}^{+0.00017}\,{\rm{s}}$
$\dot{P}({T}_{0})$	$(-3.94\pm 0.80)\times {10}^{-12}\,{\rm{s}}\,{{\rm{s}}}^{-1}$

Download table as:  ASCIITypeset image

If purely due to general relativity, the decay should be governed by the expression given in Equation (4):

$$\begin{eqnarray}&&{\dot{f}}_{\mathrm{GW}}=\displaystyle \frac{96}{5}{\pi }^{\tfrac{8}{3}}{\left(\displaystyle \frac{{{GM}}_{c}}{{c}^{3}}\right)}^{\tfrac{5}{3}}{f}_{\mathrm{GW}}^{\displaystyle \frac{11}{3}},\end{eqnarray} \tag{ 4 }$$

where ${M}_{c}=\tfrac{{({M}_{1}{M}_{2})}^{\tfrac{3}{5}}}{{({M}_{1}+{M}_{2})}^{\tfrac{1}{5}}}$ is the chirp mass, c is the speed of light, and ${f}_{\mathrm{GW}}=\tfrac{2}{P}$ is the gravitational-wave frequency. Thus, we can use the measured decay rate to constrain the chirp mass of the system. Such decay has been observed in several known detached double degenerate systems, including ZTF J1539+5027 (Burdge et al. 2019) and SDSS J0651+2824 (Hermes et al. 2012), as well as some accreting systems such as HM Cancri (Strohmayer 2005).

Equation (5) gives an estimate for tidal contribution to the measured orbital decay:

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{P}}_{\mathrm{tide}}}{{\dot{P}}_{\mathrm{GW}}}\simeq \displaystyle \frac{3({M}_{A}+{M}_{B})}{{M}_{A}{M}_{B}}\left[{\kappa }_{A}{M}_{A}{\left(\displaystyle \frac{{R}_{A}}{a}\right)}^{2}+{\kappa }_{A}{M}_{2}{\left(\displaystyle \frac{{R}_{B}}{a}\right)}^{2}\right],\end{eqnarray} \tag{ 5 }$$

where ${\dot{P}}_{\mathrm{tide}}$ is the tidal contribution to the orbital decay, expressed as a fraction of the contribution to orbital decay from the emission of gravitational radiation, ${\dot{P}}_{\mathrm{GW}}$, which is a function of the ratio of the radii to the semimajor axis of the two components, the component masses, and also κ_a and κ_b, dimensionless constants reflecting the internal structure of the white dwarf (Burdge et al. 2019).

Using the masses derived from the combined analysis, as well as κ_a = 0.14 and κ_b = 0.066 based on white dwarf models, we estimate a tidal contribution of $\tfrac{{\dot{P}}_{\mathrm{tide}}}{{\dot{P}}_{\mathrm{GW}}}\simeq 0.02$. We correct the chirp mass inferred from $\dot{P}$ by accounting for this tidal contribution, though it does not significantly alter the solution, as it is less than one-sixth of the measurement uncertainty on $\dot{P}$. As the measurement of $\dot{P}$ improves with time, the tidal contribution will become large compared to the measurement uncertainty, and will become the dominant source of uncertainty on the chirp mass.

In order to estimate the physical parameters of the system, we combined the measurements of orbital kinematics and orbital period decay with the lightcurve models fit to the CHIMERA g' and i' data taken on 2017 December 14–15. The atmospheric fit parameters only entered this modeling as a basis for deriving the limb- and gravity-darkening coefficients used in the lightcurve modeling. The free parameters of this analysis were the masses M_A and M_B, the inclination i, the radius R_B, and the time of minimum light t₀. We used the multinest algorithm (Feroz et al. 2009) to sample over these free parameters, and present the corner plots of the resulting analysis in Figure 7. The inferred parameters from this analysis, as well as the atmospheric analysis and timing analysis, are listed in Table 3.

Zoom In Zoom Out Reset image size

PTF J0533+0209 has a corresponding entry in the Gaia DR2 catalog (Gaia Collaboration et al. 2018), with a measured parallax of $\bar{\omega }=0.4741$ and uncertainty ${\sigma }_{\bar{\omega }}=0.4786$. Since we are in a regime where ${\sigma }_{\bar{\omega }}\gt \bar{\omega }$, we cannot simply infer the distance as the reciprocal of the parallax. We adopt an exponentially decreasing space density prior (Bailer-Jones et al. 2018), with an assumed characteristic length scale of 400 pc (Kupfer et al. 2018). Using this technique, we infer a distance of $D={1.5}_{-0.5}^{+0.7}\,\mathrm{kpc}$, and we adopt this distance for our estimate of the LISA gravitational-wave strain of the system. If we assume a length scale of 200 pc, the solution becomes $D={1.05}_{-0.3}^{+0.4}\,\mathrm{kpc}$, and for a length scale of 800 pc, $D={2.3}_{-0.9}^{+1.4}\,\mathrm{kpc}$. We would like to note that using our inferred radius and temperature, we estimate a distance of $D\,={2.7}_{-0.2}^{+0.2}\,\mathrm{kpc}$, farther than any of these estimates.

To calculate the characteristic strain (Korol et al. 2017), S_c, we use the expression in Equation (6), where c is the speed of light, and T_obs is the operation time of the LISA mission:

$$\begin{eqnarray}&&{S}_{c}=\displaystyle \frac{2{\left({{GM}}_{c}\right)}^{5/3}{\left(\pi f\right)}^{2/3}}{{c}^{4}D}\sqrt{{{fT}}_{\mathrm{obs}}}.\end{eqnarray} \tag{ 6 }$$

When computing the characteristic strain, we marginalized over the posterior distributions for the masses derived from our combined analysis, and the distance estimated using the parallax measurement and the assumption of a 400 pc length scale. The point and uncertainties corresponding to PTF J0533+0209 in Figure 8 reflect the posterior distribution of the characteristic strain of this object.

Zoom In Zoom Out Reset image size

However, in order to compute the estimated S/N of the source in LISA, we must account for instrument response, as well as the source position in the sky, and the source inclination. Thus, to compute the S/N, we calculate the amplitude of the signal at the detector, $A=\sqrt{| {F}_{+}{| }^{2}| {h}_{+}{| }^{2}+| {F}_{\times }{| }^{2}| {h}_{\times }{| }^{2}}$, which depends on the two gravitational-wave polarizations, h₊ and h_×, and LISA's response patterns in coupling to these polarizations, F₊ and F_× (Korol et al. 2017).

Based on its component masses, inclination, and distance of $D={1.5}_{-0.5}^{+0.7}\,\mathrm{kpc}$ based on the parallax measurement and 400 pc length scale, we estimate that LISA will detect PTF J0533+0209 with an S/N of ${8.4}_{-3.0}^{+4.2}$ at the end of the nominal four-year mission lifetime (Figure 8). If we instead use the distance of $D={2.7}_{-0.2}^{+0.2}\,\mathrm{kpc}$ inferred from the radius, temperature, and apparent magnitude of the object, the estimated signal to noise after four years is ${4.7}_{-0.9}^{+0.9}$. By the launch of LISA in 2034, we estimate that the orbital decay of PTF J0533+0209 will have caused the ellipsoidal modulation to have drifted by >900 s in phase since the initial PTF epochs, allowing for a precise constraint on the chirp mass with regular monitoring over this time period. With the chirp mass constrained, the strain measured by LISA will enable a direct probe of the distance to the object, helping verify the estimates of radius and surface temperature in this analysis. Perhaps the best example of LISA's utility is that it will be able to provide an independent measurement of the inclination angle of the system (Shah et al. 2012). This is because the gravitational-wave signal can be decomposed into two polarization components, h₊ and h_×, and the strain amplitudes of these components include factors of $\left(1+{\cos }^{2}(i)\right)$ and 2 cos(i), respectively (Thorne 1987). Combined with the measurements of K_B, and $\dot{P}$, this will provide a constraint on both component masses independent of the model-dependent fit of the ellipsoidal modulation.

PTF J0533+0209 provides interesting constraints on compact binary formation physics. In particular, the combination of the extremely low-mass white dwarf (ELM WD), high surface temperature, and low surface hydrogen (H) abundance is difficult to understand, since most ELM WDs are H-rich (DA) WDs due to helium (He) sedimentation. H-deficient CO WDs, such as DB WDs, are thought to originate from AGB stars undergoing a late He-burning shell flash during the early WD cooling phase. The shell flash induces convection that subsumes the H envelope, mixing H deep into the star where it is quickly burnt (Althaus et al. 2005). However, the He-core ELM WD in PTF J0533+0209 could not have undergone such a He-burning pulse, so the H could not be destroyed in the same way that it is for DB WDs. Nor is it likely the H has been totally stripped due to mass transfer like in AM CVn systems, because there is no reason to think that mass transfer has yet begun in this system. Several H-poor proto-ELM WDs have been discovered (see, e.g., Kaplan et al. 2013; Gianninas et al. 2016), though these have higher H-abundances and low surface gravities, and can likely be explained by weaker gravitational settling opposed by rotational mixing (Istrate et al. 2016).

To understand the formation of PTF J0533+0209, we constructed models of binary stars using the MESA stellar evolution code (Paxton et al. 2011, 2013, 2015). First, we evolve a 1.2 M_⊙ star up the red giant branch until its core helium mass is 0.187 M_⊙. We then strip the outer layers until the star has a mass of 0.19 M_⊙, leaving it with 10⁻³ M_⊙ of hydrogen, and place it in a binary with a 0.66 M_⊙ companion in a 1 hour orbit. This mimics a CE event that births a compact binary system with masses consistent with our measurements for PTF J0533+0209. We evolve this system using MESA's binary module, allowing for orbital decay via gravitational radiation, while simultaneously evolving the structure of the ELM WD. We enforce tidal synchronization at all times, as short-period binary WDs are expected to be nearly tidally synchronized (Fuller & Lai 2012a). Rotational mixing is included via MESA's implementation of Eddington-Sweet circulation with am_D_mix_factor = 1 and D_ES_factor = 1. We note that this amount of rotational mixing is very optimistic, as it is a factor of 30 larger than expected (Heger et al. 2000), yet we shall see below that even this excessive estimate of rotational mixing cannot explain the observed hydrogen depletion.

Figure 9 shows the evolution of the system described above. As the system ages, the WD contracts, and the model radius is $R\simeq 0.05\,{R}_{\odot }$ at a 20 minute orbital period. Because of extra mixing in our model (see below), the temperature initially increases as H is mixed into the interior and burned, but the temperature then decreases as H is depleted. The mass, temperature, and radius of our model at 20 minutes are all approximately consistent with PTF J0533+0209. If the post-CE period is substantially longer than one hour, the ELM WD cools to temperatures lower than observed. So, the system can only be explained by the CE formation channel if the ELM WD is born at a fairly short orbital period.

Zoom In Zoom Out Reset image size

While the match with observations above is encouraging, it is difficult to explain the low H abundance from this model, which initially has a roughly equal surface mass fraction of H and He after the CE event. The Eddington-Sweet mixing keeps the outermost layers of the WD well mixed due to the tidally enforced rapid rotation of the ELM WD, but it does not mix H into the deep interior. In the limited number of models we have explored, rotational mixing is not enough to mix H down into the He core where it can be burned, allowing the surface H to be depleted. In order to greatly deplete the surface H abundance, we have added an ad hoc mixing diffusivity of ${D}_{\mathrm{mix}}=2.5\times {10}^{3}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, which causes H to mix deep enough to burn. Figure 9 shows that the surface mass fraction of H is depleted by a factor of 50 down to ∼10⁻² in this model, though it is still much larger than our measurements indicate. Without the extra mixing, the surface mass fraction remains nearly constant. Our models do not include gravitational settling and diffusion, which would cause the surface mass fraction of H to increase as He gravitationally settles. The burning of H in our model keeps it warmer compared to a model without mixing, so our model cools slower than a normal WD of the same mass.

This model is not meant to be a "fit" to the observed properties, and some other combination of mass, initial orbital period, extra mixing, etc., can likely provide a better match to the data. Our main claim is that some sort of extra mixing or other method of removing H from the ELM WD must be at play. Since PTF J0533+0209 is in a very short-period binary compared to other H-rich ELMs, we suspect the extra mixing may be related to rotational or tidal effects that are much stronger in tight binaries. Also, if PTF J0533+0209 evolved from the CE channel, the system was likely born at an orbital period less than a couple hours in order for the ELM WD to be so hot. Based on the findings of Fuller & Lai (2013), we find that tidal heating is unlikely to heat up the ELM WD to the observed temperature at an orbital period of 20 minutes. From Equation (10) of Burdge et al. (2019), the "tidal" temperature to which the ELM WD can be heated is only ∼7000 K, much below the observed temperature. A tidally induced nova (Fuller & Lai 2012b) could heat up the ELM WD, but these events are not expected to occur at such long orbital periods.

It may be possible that H-burning flashes contribute to the depletion of H near the surface, as such flashes in a compact binary will cause the ELM to overflow its Roche lobe and lose much of its expanding H envelope. The model described above had a small enough initial H mass that it did not encounter such flashes, but we have attempted to simulate flashing models with a larger initial H mass. For numerical reasons, we have not been able to simulate the flash-induced mass loss in a compact binary; however, the models indicate that the H at the base of the envelope would not overflow the Roche lobe during these events and would likely be retained. Hence, we believe H flashes may have contributed to the depletion of H, but some extra mixing (e.g., rotational mixing, or mixing induced via tidally excited gravity waves) is likely required to further deplete the remaining H to the observed level. One possibility is that the H-burning flashes induce more mixing below the convective H-burning shell than predicted by our models, a phenomenon also thought to occur during classical novae (e.g., Alexakis et al. 2004; Denissenkov et al. 2013). This could mix H downward where it is burnt and dilute the H shell with dredged-up He, potentially reducing the observed H abundance to the level observed.

We have also attempted to create a system like PTF J0533+0209 from a stable mass transfer scenario like that described in Sun & Arras (2018), Tauris (2018), and Li et al. (2019). We evolve systems near the orbital period bifurcation (Podsiadlowski et al. 2003) such that the donor is evolved (and hence H-poor) and the binary shrinks to sub-hour orbital periods via magnetic braking and gravitational radiation. This channel typically forms low-mass ELM WDs with masses in the range $0.15\,{M}_{\odot }\,\lesssim M\lesssim 0.2\,{M}_{\odot }$ that detach at very short orbital periods, so it can certainly produce systems similar to PTF J0533+0209. However, the donor in these models usually has lower surface temperature and higher H abundance than observed in PTF J0533+0209. Hence, while it might be possible to form a system like PTF J0533+0209 through stable mass transfer, we consider the CE channel to be more likely.

Based on its measured orbital decay rate, PTF J0533+0209 has a characteristic decay timescale of ${\tau }_{c}=\tfrac{3}{8}\tfrac{P}{| \dot{P}| }\approx {\rm{3,700,000}}$ yr. However, in a fraction of this time, PTF J0533+0209B will overflow its Roche lobe and initiate mass transfer onto PTF J0533+0209A, which occurs at an orbital period of ≈14 minutes in Figure 9, though precisely when this occurs depends on the rate at which PTF J0533+0209B cools and shrinks compared to its merger timescale. The binary will continue to evolve to shorter orbital periods due to the emission of gravitational radiation, and eventually mass transfer will begin to remove mass from the degenerate helium core of PTF J0533+0209B (Kaplan et al. 2012), which will expand in response to this mass loss, consequently increasing the mass transfer rate and producing a feedback loop. Given the mass ratio of the binary, this means that PTF J0533+0209 will likely evolve into an AM CVn system at longer orbital periods (Marsh et al. 2004) or, alternatively, could result in a merger producing an R Coronae Borealis–type star (Paczyński 1971).