{'fraction_non_alphanumeric': 0.03763668217241588, 'fraction_numerical': 0.03186733522624298, 'mean_word_length': 4.162482457508187, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Formation of the double white dwarf binary PTF J0533+0209 through stable mass transfer?

Double white dwarf (DWD) binaries are important for studies of common-envelope (CE) evolution, Type Ia supernova progenitors and Galactic sources of low-frequency gravitational waves (GWs). PTF J0533+0209 is a DWD system with a short orbital period of Porb ~ 20 min and thus a so-called LISA verification source. The formation of this system and other DWDs is still under debate. In this paper, we discuss the possible formation scenarios of this binary and argue that it is not likely to have formed through CE evolution. Applying a new magnetic braking prescription, we use the MESA code to model the formation of this system through stable mass transfer. We find a model which can well reproduce the observed WD masses and orbital period but not the effective temperature and hydrogen abundance of the low-mass He WD component. We discuss the possibility of using H flashes to mitigate this discrepancy. Finally, we discuss the future evolution of this system into a AM CVn binary such as those that will be detected by space-borne GW observatories like LISA, TianQin and Taiji.


INTRODUCTION
Recently, a number of double white dwarfs (DWDs) with very short orbital periods have been discovered, such as ZTF J1539+5027 (∼ 7 min orbital period, Burdge et al. 2019b and PTF J0533+0209 (∼ 20 min orbital period, Burdge et al. 2019c). In these DWDs, the primary WD is a CO WD and the companion star is an extremely low-mass (ELM) He WD with a mass smaller than ∼ 0.2 M . Because of their short orbital periods, they are important Galactic gravitational wave (GW) sources that will be detected by space-borne observatories, such as LISA (Amaro-Seoane et al. 2017), TianQin (Luo et al. 2016), Taiji (Ruan et al. 2020). However, the formation of these systems is not well understood.
The formation of DWDs has been widely investigated in the past (Iben & Tutukov 1986;Lipunov & Postnov 1988;Yungelson et al. 1994;Han et al. 1995;Iben et al. 1997;Han 1998;Sun & Arras 2018;Li et al. 2019). It is generally accepted that there are two possible formation scenarios for DWDs (see also Fig. 2 in Li et al. 2019). In the first formation scenario, a binary system consisting of a CO WD and a low-mass donor star undergoes stable mass transfer via Rochelobe overflow (RLO) and evolves into a binary system consisting of a CO WD and a He WD. The systems produced from this scenario will follow the He WD massorbital period relation (e.g. Tauris & Savonije 1999). In the second scenario, a binary system with a CO WD and a red-giant companion star is subject to unstable mass transfer and enters a common-envelope (CE) phase. If the system successfully ejects the CE, then it will evolve into a CO WD+He WD binary. Hereafter, we refer to these two scenarios as the stable RLO scenario and the CE scenario, respectively. Sun & Arras (2018) have studied the formation of ELM WDs from binary evolution using a mass range of the progenitors of the ELM WDs between 1.0−1.50 M . They found that ELM WDs can be produced from non-conservative mass transfer and are not likely to be formed via CE evolution. Li et al. (2019) have comprehensively studied the formation of ELM WDs through both stable RLO and CE evolution. They found that the ELM WDs can form from both scenarios and argue that ELM WDs with a mass smaller than ∼ 0.22 M are formed from the stable RLO scenario, whereas systems with ELM WDs > 0.22 M formed via the CE scenario.
Applying so-called convection and rotation-boosted magnetic braking, Soethe & Kepler (2021) studied the formation of ELM WDs with NS companions and found that ELM WD masses of 0.15−0.27 M can be produced via stable RLO from low-mass X-ray binaries (LMXBs) with initial orbital periods of 1 − 25 d.
For the very tight-orbit DWD system with an ELM He WD component above a certain mass limit the origin is clear. As an example, one may consider ZTF J1539+5027 with an orbital period of ∼ 7 min and an ELM He WD mass of ∼ 0.21 M . If this system were produced from stable RLO, then the orbital period at the end of mass transfer would have been around 2.4 d according to the WD mass-orbital period relation from Tauris & Savonije (1999). It would subsequently have taken the system around 3.4 × 10 12 yr to evolve to the current orbital period of 7 min via GW damping. Therefore, this system cannot have been produced from stable RLO.
Regarding the system PTF J0533+0209, which is the focus of this paper (see Table 1 for parameters), Burdge et al. (2019c) suggested that it is formed through a CE scenario. But they did not discuss whether the CE can be ejected in order to produce such a system. The results of Sun & Arras (2018) and Li et al. (2019) show that this system is not likely to be formed from a CE scenario. On the other hand, it is unclear from the binary models of Li et al. (2019); Soethe & Kepler (2021) if PTF J0533+0209 can be produced from stable RLO and subsequently decay its orbit sufficiently (once detaching from the WD mass-orbital period correlation) within a Hubble time to reproduce its present orbital period of 20.6 min. Therefore, this system is puzzling as it is apparently also hard to explain its origin with binary evolution models with stable RLO. Here we aim to investigate this system in more detail as it represents a class of DWD systems challenging current knowledge of binary stellar evolution.
Recently, Chen et al. (2021) investigated the formation of binary ELM WDs with radio millisecond pulsar (MSP) companions. Using an intermediate strength magnetic-braking prescription ("MB3") suggested by Van et al. (2019), they found that the initial parameter space of LMXBs for producing such MSP+ELM WD binaries becomes significantly larger compared to models with a standard magnetic braking, and thereby offering a potential solution to the "fine-tuning problem" of producing MSP+ELM WD binaries (Istrate et al. 2014). However, they also found that with the MB3 prescription they cannot reproduce many known wide-orbit binary MSPs. Motivated by this work, we will here study the formation of the system PTF J0533+0209 using a magnetic braking prescriptions from Van et al. (2019). The purpose of this paper is to discuss possible formation scenarios for PTF J0533+0209 and present a valid binary formation model for this system. The paper is structured as follows. In Section 2, we demonstrate that formation of PTF J0533+0209 through CE evolution is not possible from a simple energy budget argument. In Section 3, we discuss the formation of this system through stable RLO and present a viable binary evolution model for this system. We compare with observational characteristics and discuss discrepancies. In addition, we discuss the future evolution of this system as an AM CVn system. A brief discussion and a summary of our results are given in Section 4 and Section 5, respectively.

FAILURE OF A CE SCENARIO FOR PTF J0533+0209
In the CE scenario, a binary system consisting of a CO WD and a (sub)red-giant star initiates unstable mass transfer and enters a CE phase. To check for the possibility of a successful CE ejection, we apply the energy budget prescription of CE evolution (Webbink 1984;Livio & Soker 1988;Han et al. 1994Han et al. , 1995: where the left-hand side of the equation is the binding energy, E bind of the envelope of the giant star (with total mass M 2 = M core + M env ) before entering the CE. Here M (r) is the mass within radial coordinate, r of the donor star and U is the internal energy that can be used in the ejection process. The latter includes the thermal energy of a simple perfect gas, the radiation energy, as well as terms related to ionization of atoms and dissociation of molecules, as well as the Fermi energy of a degenerate electron gas (Han et al. 1994(Han et al. , 1995. The right-hand side expresses the release of orbital energy from in-spiral, ∆E orb and 0 < α CE < 1 is the efficiency of CE ejection (i.e. fraction of released orbital energy used to eject the CE). Here, M core , M env and M CO are the masses of the core, the envelope and that of the in-spiralling CO WD, respectively, and a i and a f denote the binary orbital separations before and after the CE phase, respectively. As a conservative limit, we first set a f to be the value of current binary separation. Following Burdge et al. (2019c), we assume that the progenitor mass of the ELM WD is 1.20 M . Using the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, version: 12115, Paxton et al. 2011, 2018, we computed the evolution of a single star with a mass of 1.20 M and obtain a binding energy of the envelope, E bind = −1.20 × 10 48 erg and a radius, R 2 = 3.10 R when the core mass is around 0.167 M (i.e. similar to the mass of the ELM He WD component in the PTF J0533+0209 system, Table 1). The initial binary separation before the CE can be computed from (Eggleton 1983): R rl a i = 0.49 q 2/3 0.6 q 2/3 + ln 1 + q 1/3 , where R rl is Roche-lobe radius and equal to the radius of the (sub)giant star, q = M donor /M CO is the mass ratio and M donor is the progenitor mass of the ELM WD, i.e. M donor = 1.20 M . Solving Eqs.
(1) and (2), we find a value of α CE = 1.75. However, this large value is not physical. Furthermore, we note that the binary separation of this system after a CE should be larger than the current orbital separation because of GW radiation since its formation until its current state, and thus the amount of released orbital energy from in-spiral would have been smaller. Therefore, the α CE value would be larger than 1.75, making the CE even more unlikely. It should noted that applying a larger progenitor mass of the ELM WD would only exacerbate the problem. Assuming e.g. M donor = 1.40 M increases the envelope binding energy, leading to α CE = 2.7. Assuming a smaller mass, e.g. M donor = 1.0 M leads to longer evolution time and α CE = 1.14. So far, the value of α CE is not well constrained (see Ivanova et al. 2013 for a review). Davis et al. (2010) computed the detailed structure of 1 − 8 M AGB stars and computed the binding energy of the envelopes of these stars. With these results, they found that it is possible to explain the population of post-CE binaries found by the SDSS with α CE > 0.10. Using a similar approach, Zorotovic et al. (2010) constrained the CE efficiency parameter to 0.20 < α CE < 0.30. In addition, Nandez et al. (2015) studied the formation of DWDs from CE evolution with 3D simulations and constrained the CE efficiency parameter to 0.25 < α CE < 0.50. From these studies, we expect 0.1 < α CE < 0.50 for low-mass binaries. In any case, α CE > 1.75 is not realistic. Therefore, we conclude that PTF J0533+0209 did not form through a CE scenario and we will now focus on investigating the stable RLO scenario in the rest of the paper.

Method and assumptions
Our binary evolution model is motivated by Chen et al. (2021), in which we modelled the formation of ELM WDs in binary MSP systems via RLO in LMXBs. Here we briefly summarize the assumptions and the method we adopt in our presented calculation of PTF J0533+0209 via RLO in a CV system.
We use the MESA code to model the binary evolution. In our calculation, the CO WD accretor is assumed to be a point mass, so we do not compute the detailed structure of this WD. The secondary star is assumed to initially be a zero-age main-sequence star with a chemical composition X = 0.70 (hydrogen), Y = 0.28 (helium) and Z = 0.02 (metals). We set the mixing length to be 2 times the local pressure scale height and do not consider overshooting in our calculation. The Reimer's wind prescription is used for the secondary star (Reimers 1975): (3) where η = 1.0 is our adopted scaling factor. M 2 , R 2 , L 2 are the mass, radius and luminosity of the secondary star, respectively.
In our model, we have taken different mechanisms of angular momentum loss into consideration, including angular momentum loss due to GW radiation, magnetic braking, mass loss and spin-orbital couplings. To compute the angular momentum loss due to GW radiation, we use the formula of Landau & Lifshitz (1971): where G is the gravitational constant; c is the speed of light in vacuum; M CO is the CO WD mass; M 2 is the secondary star mass; a is the orbital separation.
Regarding the angular momentum loss due to magnetic braking, the standard Skumanich law (i.e. β = ξ = α = 0 in Eq. 5) is usually applied to solar-type main sequence stars. It may be not suitable for systems with donors different than the Sun, since the stellar wind and the strength of the magnetic field can be very different which is important for the calculation of magnetic braking. By including a scaling of the strength of magnetic braking with convective turnover time and wind mass-loss rate, Van et al. (2019) modelled the evolution of LMXBs and found that the observed properties of persistent X-ray binaries are in better agreement with observations using a model with β = 0, ξ = 2 and α = 1 (MB3), compared to models with other choices of β, ξ and α. Here, in our study of DWDs we adopt their general prescription, which can be described with the following equations: (6) We adopt γ mb = 4, τ ,conv = 2.8 × 10 6 s andṀ ,w = 2.54 × 10 −14 M yr −1 . ω is the orbital angular velocity and τ conv is the turnover time of convective eddies.Ṁ w is the wind mass-loss rate of the secondary star.
A caveat of the applied model MB3 is that Chen et al. (2021) could not reproduce the many known wide-orbit binary MSPs using this magnetic braking prescription (because the drain of orbital angular momentum is simply too strong if applying MB3). Nevertheless, here we explore binary evolution using MB3 as our default model.
We adopt the "isotropic re-emission model" for the mass transfer via RLO (see e.g. Tauris & van den Heuvel 2006, for a review). In this model, the material flows conservatively from the donor star to the WD. The further fate of this material will strongly depend on the accretion rate (see Chen et al. 2019). If the accretion rate is larger than the maximum rate of steady burning, the excess material will be lost in the form of an optically thick wind (Hachisu et al. 1996). Only in a narrow range of accretion, the H or He can burn stably, and the steady burning regimes of H and He do not overlap with each other (Nomoto 1982). Therefore, even in the steady burning regime of H, not all material can be retained by the WD because of unsteady burning of He. If the accretion rate is lower than the minimum accretion rate of stable burning, the material will burning unstably and most of it will be ejected in violent shell flashes and nova outbursts and the WD may even reduce its mass if the accretion rate is too low (Hachisu et al. 1999). Keeping this and also the low mass of the CO WD of PTF J0533+0209 in mind, we simply assume that, as an average for the entire RLO phase spanning a wide range of mass-transfer rates (Ṁ 2 ), a fraction of 90% of the transferred material is lost from the CO WD accretor and takes away with it the specific orbital angular momentum of this WD.
The mass-transfer rate,Ṁ 2 from the secondary star to the WD is computed with the scheme proposed by Ritter (1988). The amount of accretion onto the WD is limited by the Eddington accretion rate (Ṁ Edd ). We set the ratio of gravitational mass to rest mass of accreted material to be 0.99 for our CO WD accretor.
In our calculation, we also take the rotation of the donor star into account. We assume that its initial rotation velocity is synchronized with the initial binary orbital period. The prescription of Paxton et al. (2015) is adopted to compute the angular momentum loss due to spin-orbital couplings. The evolution of angular velocity is computed with the following formulae: and where ω j is the angular velocity at the face of cell j and ω orb is the orbital angular velocity, τ sync is the syn-chronization time, r g is the radius of gyration, (k/T ) c is the ratio of the apsidal motion constant to the viscous dissipation timescale, R is the stellar radius and a is the binary separation. The synchronization time depends on the structure of stars. For stars with convective envelopes, turbulent viscosity in the convective region, which may operate on the equilibrium tide, is the most efficient form of viscous dissipation. For stars with radiative envelopes, the resonances of the free gravity modes are damped by radiative dissipation which operates on the dynamical tide. The value of (k/T ) c is different for these two kinds of stars and can be calculated as in Hurley et al. (2002). In our calculation, we do not include the internal heating due to the tidal fraction in the thermal evolution of the donors. Chemical mixing and transport of angular momentum induced by different kinds of instabilities resulting from rotation (Heger et al. 2000) are also taken into account, i.e.: chemical mixing induced by secular shear instability, Eddington-sweet circulation, dynamical shear instability and Goldreich-Schubert-Fricke instability. Following Heger et al. (2000), the mixing efficiency factor in our calculation is assumed to be f c = 1/30. With respect to the transport of angular momentum, we only consider the Spruit-Tayler dynamo (Spruit 2002;Heger et al. 2005). Following Heger et al. (2000), we set the parameter f µ = 0.05, which regulates the inhibiting effect of chemical gradients on the mixing process. Motivated by Istrate et al. (2016), we also include element diffusion induced by gravitational settling and chemical and thermal diffusion in our calculation.

Initial values
In order to find a model which can reproduce the properties of PTF J0533+0209, we compute a grid of binary evolution with the above assumptions. Given the low mass of the currently observed CO WD (M CO = 0.652 M ), we assume an initial CO WD mass of M CO,i = 0.55 M . The donor (secondary star) mass ranges between M 2,i = 1.0 − 1.4 M in steps of ∆M 2 = 0.05 M . The initial orbital period ranges between P orb,i = 1.0 − 100 d in steps of ∆log(P orb,i /days) = 0.1. With this grid, we evolve models resulting in final properties relatively close to the properties of PTF J0533+0209. Next, we decrease the steps of ∆P orb,i and find improved models which can better reproduce the observed properties of PTF J0533+0209. Based on these models, we also compute some models with different WD or donor masses in order to understand the influence of these parameters on our results.

Results
In this section, we present our binary evolutionary model with properties closest to those of PTF J0533+0209. In this model, the initial values are M 2,i = 1.40 M , P orb,i = 7.00 d and M CO,i = 0.55 M . Figure 1 displays the HR diagram of the secondary star (top), the evolution of mass-transfer rate (middle), and orbital and spin periods (bottom), as a function of time. In the middle panel, two phases of RLO mass transfer are identified: cataclysmic variable (CV) and AM Canum Venaticorum (AM CVn). PTF J0533+0209 is located between these two stages.
A more qualitative description of this example is as follows: when the secondary star evolves beyond the main sequence (after termination of core hydrogen burning) it enters the Hertzsprung gap while its radius increases. Meanwhile (starting already from the ZAMS) the binary separation continuously decreases because of the magnetic braking. Therefore, the secondary star fills its Roche lobe and initiates mass transfer (and the system becomes a CV system) already at an age of t = 2.95 Gyr when its radius is R 2 = 2.53 R . This result can be understood as follows. Since the He WD mass in PTF J0533+0209 is quite small, its progenitor star should start mass transfer rather early in its stellar evolution shortly after it has produced a small core mass. And the ZAMS mass of the He WD progenitor should not be too large, since in that case the RLO is no longer stable. In addition, early on the main sequence, the donor star is slightly out of synchronization. As the binary orbital period decreases significantly due to magnetic braking, the effect of tides becomes stronger and the spin of the donor star becomes synchronized with the orbit. At the very beginning of the first mass-transfer phase, because of mass loss, the donor star spins down and is out of synchronization with the orbit (see inset in bottom panel). However, due to the strong tides, it becomes synchronized again after a short time (less than 5 Myr).
At the end of the first mass-transfer phase (t = 3.23 Gyr), the secondary star detaches from its Roche lobe and slowly settles on the He WD cooling track with a mass of M 2 = 0.155 M . During the RLO, the CO WD has accreted ∼ 0.12 M and now has a mass of M CO = 0.669 M and the orbital period is P orb = 0.22 d. The system subsequently evolves as a detached DWD system, as presently observed, with decreasing orbital separation due to GW radiation. Around t = 10.661 Gyr, the He WD fills its Roche lobe and starts to transfer material onto the CO WD and the system evolves as an AM CVn system. Slightly before that stage (t = 10.660 Gyr), when the orbital period and  Table 1), while the green dot with error bars indicates the actual location of PTF J0533+0209. In the middle panel two phases of RLO is shown: CV and AM CVn. Our evolved system has an age of t = 10.660 Gyr. The middle panel shows the evolution of mass-transfer rate as a function of time. In the lower panel, the evolution of orbital (blue solid line) and spin period (green dashed line) are shown as a function of time. The two grey shaded regions indicate the two mass-transfer phases. The inset panel shows the evolution of spin and orbital periods around the onset of the first mass-transfer phase.
masses are close to those observed for PTF J0533+0209, the system has the physical parameters listed in Table 1.
During the AM CVn phase, the spin of the donor star is always out of synchronization with the orbit. This is because that the donor star spin increases due to the mass loss and the effect of tides becomes weaker as the orbital period increases. Figure 2 shows a more detailed comparison of our results with observations. Concluding from the values in Table 1 and Fig. 2, we find that our model can very well reproduce the WD masses and the orbital period of PTF J0533+0209, as well as the ELM He WD radius (and log g). However, the simulated effective temperature (T eff , and thereby also the luminosity, L) and the surface H-abundance (X surf ) of the ELM He WD represent quite a mismatch. In the grid we computed, we have checked if other models can better reproduce T eff and X surf . However, we conclude that this is not the case (see discussions in Section 4).
For our initial binary systems, if we increase the donor star mass (M 2,i ) or the orbital period (P orb,i ), we produce He WDs with too large masses and too small radii. If we decrease M 2,i or P orb,i , the results is that the timescale of nuclear evolution is too long to reproduce PTF J0533+0209 or the system never detaches from RLO (producing a so-called converging system, see e.g. Istrate et al. 2014), respectively. If we increase the initial CO WD mass (M CO,i ) of the accretor, the results are not very significantly different. To compensate for a different applied value of M CO,i we can simply adjust the accretion efficiency to reach the value of M CO = 0.652 M for PTF J0533+0209). 4. DISCUSSION 4.1. A possible way to mitigate the discrepancy As we found in the previous section, our model cannot explain the high value of T eff and the low H abundance of the ELM He WD in PTF J0533+0209. In the following, we will discuss the possible explanation for this discrepancy. However, we start by briefly summarizing the investigation by Burdge et al. (2019c).
The peculiar observed values of the above-mentioned parameters were also discussed by Burdge et al. (2019c) in their attempt to model PFT J0533+0209. Whereas we have argued that a CE origin of PTF J0533+0209 seems unlikely based on simple energy considerations (Section 2), they disregarded this aspect and instead attempted to explain PTF J0533+0209 by modelling a post-CE system in composite steps: first, they evolved a single 1.20 M star up the red-giant branch until its core mass is 0.187 M . They then stripped the outer layers until the star had a mass of 0.19 M and a total amount of hydrogen of 10 −3 M . Subsequently, they placed this model star in a P orb = 1 hr binary system with a 0.66 M companion and evolved it further using MESA. Despite introducing excessive rotational mixing, the surface abundance of hydrogen remained too high. They then computed a model with an extremely large ad-hoc mixing diffusivity, which enforces hydrogen to mix deep enough to burn and thereby reduce the value of X surf . Burdge et al. (2019c) concluded that based on the findings of Fuller & Lai (2013), the tidal heating is unlikely to heat up the ELM He WD to the observed temperature at an orbital period of 20 min. Referring to the work of Burdge et al. (2019a), they argued that the "tidal temperature" to which the ELM He WD can be heated is only ∼ 7000 K, much below the observed temperature of T eff ∼ 20 000 K. Another solution discussed by Burdge et al. (2019c) is the possibility of Hshell flashes (e.g. Driebe et al. 1998;Althaus et al. 2001;Nelson et al. 2004;Althaus et al. 2013;Istrate et al. 2016). However, due to numerical reasons, they were not able to simulate flash-induced RLO of the residual H-envelope. Such RLO could substantially reduce the amount of residual hydrogen (see below). Finally, Burdge et al. (2019c) attempted to model PTF J0533+0209 via the stable RLO channel as explored in this paper. However, due to the problems of reproducing T eff and X surf , they concluded that the CE channel is more likely.
Despite the initial issues with T eff and X surf , we are still optimistic with respect to the stable RLO channel. The residual H-envelope mass expected for ELM He WDs is typically between ∼ 10 −3 − 10 −2 M depending on WD mass and evolutionary epoch. Figure 3 displays the HR diagram of one of our grid model computations of the evolution of an ELM He WD model with a final WD mass of 0.156 M . It covers the evolution from detachment of the CV phase as a 0.163 M proto-ELM He WD (t = 3.2 Gyr and P orb = 6.92 hr) until t = 13.7 Gyr (P orb = 5.97 hr). The WD experiences several strong H-shell flashes during its evolution, leading to temporary episodes of additional RLO. As a consequence, the final WD mass is reduced to 0.156 M . From the plot it is seen that following each flash episode the surface H-abundance is considerably reduced, although in all cases X surf 0.10. Thus, we cannot explain the observed value of X surf 4.9 × 10 −4 , which may indeed suggest that extra mixing is needed, as proposed by Burdge et al. (2019c). Nevertheless, we speculate, that it is quite possible that PFT J0533+0209 is currently observed shortly after one of these shortlived H-shell flash episodes, thereby helping to explain the measured values of T eff and L. A caveat is that both X surf and L (T eff ) decay back to pre-flash values on a timescale of order 10 Myr. Hence, the time window for observing the system in an immediate post-flash evolutionary stage is rather limited.

Influence of the accretion disk
In our calculations, we did not consider the effect of an accretion disk during the AM CVn phase. In this phase, the binary separation is very small and the formation of accretion disk is not obvious. In an early paper, Nelemans et al. (2001) studied this question with an analytical method. If the minimum distance at which the accretion flow passes the accretor is comparable to the radius of the accretor, it is likely that no accretion disk will be formed and the accretion flow will directly hit the surface of the accretor. From their fig. 1, it is suggestive that no accretion disk is formed at the beginning of the AM CVn phase whereas an accretion disk will appear at later times, after the system passes through the orbital period minimum and once the donor star (He WD) mass is below ∼ 0.08 M . Without an accretion disk, loss of angular momentum from the orbit may destabilize the mass transfer unless the angular momentum lost to the accretor can be transferred back to the orbit (Marsh et al. 2004). Later, it was argued by Kaplan et al. (2012) that for AM CVn systems with ELM He WD donors, the initial contraction of the ELM He WD due to mass loss allows for more stable mass transfer than that originally found for cold He WDs studied by Marsh et al. (2004). In our MESA models, we do take into account the entropy (finite temperature) of the ELM He WD, but the detailed treatment of direct impact versus disk formation is not included.

A LISA verification source
PTF J0533+0209 is expected to be detected by LISA (i.e. a so-called LISA verification source Kupfer et al. 2018). However, during the subsequent GW damping of the system, the more extended ELM He WD component will fill its Roche lobe and initiate RLO towards the CO WD companion. During this epoch, the system will be observable as an AM CVn system (Nelemans et al. 2001;van Roestel et al. 2021). Discussions and detailed tracks of their evolution through the LISA GW frequency band during this phase is given in e.g. Tauris (2018). Figure 4 shows the evolution of characteristic strain as a function of GW frequency for our binary model shown in Fig. 1. In the calculation, we assume the distance of the system to be d = 1.5 kpc following Burdge et al. (2019c). Besides LISA, there are other low-frequency GW observatories which will be lunched in the next two decades, such as TianQin (Luo et al. 2016) and Taiji (Ruan et al. 2020). In the plot, we therefore also present the sensitive curves for the GW observatories TianQin, assuming an observational time T = 5 yr, Taiji with an observational time T = 5 yr, besides LISA with an assumed observational time T = 4 yr. We notice that a system like PTF J0533+0209 is detected both during the detached phase and the subsequent AM CVn phase. From this plot, we find a possible signal-tonoise ratio (S/N) up to 10 − 100 for simulated systems similar to, but further evolved than, PFT J0533+0209. The chirp mass of the detectable sources are between 0.10 − 0.27 M (steadily decreasing with evolutionary age, see right-hand side colour bar).  Table 1 is shown with a blue star.
In this paper, we have discussed possible formation scenarios of the DWD binary system PTF J0533+0209. We found that explaining its origin via the CE channel (as suggested by Burdge et al. 2019c) is difficult since the low mass of the ELM He WD implies CE ejection at an early evolutionary stage (Hertzsprung gap) when the binding energy of the envelope of the progenitor star is still larger (by a factor of ∼ 2) than the amount of orbital energy released during spiral-in.
We have instead presented a formation scenario based on stable RLO in a CV system. Using the MESA code in combination with a magnetic braking prescription from Van et al. (2019), we have modelled the system from the ZAMS and found a solution that matches the two measured WD component masses and the present orbital period of the system. However, our model fails to reproduce the observed T eff and, in particular, the surface H-abundance of the ELM He WD. We have discussed possible reasons for this and point to a recent H-shell flash (combined with efficient mixing) as the most plausible explanation.
The present GW signal of this system will be detected by low-frequency GW observatories such as LISA, Tian-Qin and Taiji. We calculated the signal of its detection and continued our computation for PTF J0533+0209 during its future AM CVn phase when the ELM He WD fills its Roche lobe in a bit more than 1 Myr from now.
We are strongly inclined to believe that PFT J0533+0209 formed via stable RLO. However, a complete and self-consistent model that explains all observables of the binary is still missing. We consider this system a very interesting test case for further investigations of close binary star evolution and input physics, and encourage the community to join our efforts. Chirp Mass (M ) Figure 4. Evolution of the characteristic strain as a function of GW frequency for the binary system presented in Fig. 1 for an assumed distance of d = 1.5 kpc. The red, grey and magenta solid lines represent the sensitive curves for TianQin with an observational time T = 5 yr, Taiji with an observational time T = 5 yr and LISA with an observational time T = 4 yr, respectively. The color on the evolutionary track represents the chirp mass of the binary system (see color bar; the grey and dark blue colors represents chirp masses of M chirp > 0.30 M and M chirp < 0.05 M , respectively). The black dot represents the end of first mass transfer phase and the star represents the onset of the second mass transfer phase. The cross represents the current location of PTF J0533+0209.