Phases of Mass Transfer from Hot Subdwarfs to White Dwarf Companions and Their Photometric Properties

Binary systems of a hot subdwarf B (sdB) star + a white dwarf (WD) with orbital periods less than 2-3 hours can come into contact due to gravitational waves and transfer mass from the sdB star to the WD before the sdB star ceases nuclear burning and contracts to become a WD. Motivated by the growing class of observed systems in this category, we study the phases of mass transfer in these systems. We find that because the residual outer hydrogen envelope accounts for a large fraction of an sdB star's radius, sdB stars can spend a significant amount of time ($\sim$10s of Myr) transferring this small amount of material at low rates ($\sim 10^{-10}$-$10^{-9}\ M_\odot\,\rm yr^{-1}$) before transitioning to a phase where the bulk of their He transfers at much faster rates ($\gtrsim 10^{-8}\ M_\odot\,\rm yr^{-1}$). These systems therefore spend a surprising amount of time with Roche-filling sdB donors at orbital periods longer than the range associated with He star models without an envelope. We predict that the envelope transfer phase should be detectable by searching for ellipsoidal modulation of Roche-filling objects with $P_{\rm orb}=30$-$100$ min and $T_{\rm eff}=20{,}000$-$30{,}000$ K, and that many ($\geq$10) such systems may be found in the Galactic plane after accounting for reddening. We also argue that many of these systems may go through a phase of He transfer that matches the signatures of AM CVn systems, and that some AM CVn systems associated with young stellar populations likely descend from this channel.


INTRODUCTION
Hot subdwarf B (sdB) stars are subluminous stars of spectral type B, thought to be compact He-burning stars of mass ≈0.3-0.6 M with thin hydrogen envelopes (Heber 2009(Heber , 2016Zhang et al. 2009;Götberg et al. 2018;Geier 2020). The low masses of these He-burning objects that have not yet gone through AGB mass loss suggest that interaction with a binary companion is needed to explain the existence of sdB stars (Maxted et al. 2001;Han et al. 2002;Pelisoli et al. 2020). Indeed, in many cases sdB stars are observed in compact binary systems with orbital periods of only a few hours, implying a prior common envelope event that ejected most of the sdB star's former H envelope and caused the companion star to spiral inward toward its current orbital configuration (Han et al. 2002(Han et al. , 2003Napiwotzki et al. 2004;Nelemans 2010;Copperwheat et al. 2011;Kruckow et al. 2021).
The most compact of these systems consist of sdB (or sometimes even hotter sdO) stars with white dwarf (WD) companions with orbital periods on the order of just one hour (Vennes et al. 2012;Geier et al. 2013;Kupfer et al. 2017aKupfer et al. ,b, 2020aPelisoli et al. 2021), which can make contact and transfer mass from the sdB/sdO star to the WD via stable Roche lobe overflow within the sdB/sdO star's He-burning lifetime (Justham et al. 2009;Brooks et al. 2015). 1 Evolution of these compact sdB+WD binaries can lead to many interesting astrophysical phenomena, including thermonuclear events due to accretion on the WD and runaway subdwarf remnants in the case of a supernova that disrupts the binary (Hirsch et al. 2005;Justham et al. 2009;Wang et al. 2009;Piersanti et al. 2014;Geier et al. 2015;Brooks et al. 2015;Neunteufel et al. 2016Neunteufel et al. , 2019Bauer et al. 2017Bauer et al. , 2019Neunteufel 2020). Some observed thermonuclear tran-sients show evidence of thick He shell detonations De et al. 2019De et al. , 2020, consistent with the binary evolution predictions for mass transfer from sdB stars to WD companions.
These binaries may also be important gravitational wave sources for the Laser Interferometer Space Antenna (LISA), motivating work toward better characterization of the population of these binaries in our Galaxy (Götberg et al. 2020). For systems in which the sdB star is at or near Roche filling, high-cadence photometry can detect gravity darkening and ellipsoidal modulation of the lightcurve due to temperature variations on the distorted subdwarf surface, as Kupfer et al. (2020a,b) recently demonstrated in discovering two new Roche-filling systems with subdwarf donor stars.
Although the hydrogen envelope of an sdB star contains only a small fraction of the star's total mass (M env ∼ 10 −4 − 10 −2 M ), the envelope extends over a significant fraction of the total radius, and therefore has a large impact on the range of orbital periods for which sdB stars are Roche-filling objects that display significant ellipsoidal modulation. Previous work has emphasized the ≈20-40 minute orbital period range for compact Roche-filling sdB cores calculated neglecting the H envelope (Savonije et al. 1986;Iben & Tutukov 1991;Piersanti et al. 2014;Brooks et al. 2015). In this paper, we use the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA) to calculate models that explore the sensitivity of sdB radii to H envelope masses and map out the space of possible orbital periods where Roche-filling mass transfer can occur. We find that when accounting for the H envelope, sdB stars can fill their Roche lobes at longer orbital periods and for a much longer portion of their lifetimes. This has significant implications for the potential range of orbital periods at which these systems can be discovered by searching for ellipsoidal modulation in high-cadence photometry.
In Section 2, we calculate a grid of MESA sdB star models to demonstrate the dependence of the radius on the thin H envelope mass. In Section 3, we describe the physics of mass transfer for Roche-filling sdB stars, and we show that the mass-radius relations from our MESA models are sufficient to predict accretion rates and timescales through both the H and He mass transfer phases. Surprisingly, the H mass transfer phase can last for > 10 Myr, as long or longer than the subsequent He mass transfer phase. In Section 4, we use a grid of MESA sdB models to map out the space of orbital periods for which systems can fill their Roche lobes and begin transferring mass. We find that H envelope mass transfer can begin at orbital periods as long as ≈100 minutes, and that 40-60 minutes may be a typical period for the onset of mass transfer. In Section 5, we verify our calculations with two full binary evolution models for sdB+WD systems in which we model both stars along with the orbital evolution and mass transfer using MESA. In Section 6, we demonstrate that the ellipsoidal variation for these sdB+WD systems will be easily detectable for many systems within a few kpc, which will enable discovery of many new systems. Finally, in Section 7 we discuss the implications of our work for targeting searches for these binary systems based on ellipsoidal modulation, and predict that at least 10s of these systems may be discovered within a few kpc in the Galactic plane after accounting for reddening.

SDB STAR RADIUS AS A FUNCTION OF ENVELOPE MASS
As we will show in Section 3, the first phase of mass transfer from sdB star donors onto WD companions depends primarily on the sdB star envelope mass-radius relation, so in this section we use grids of MESA models to show this relation for use in later sections. Our MESA models all use release version 15140 . The MESA equation of state (EOS) is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000), and PC (Potekhin & Chabrier 2010). Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, and electron conduction opacities are from Cassisi et al. (2007). Nuclear reaction rates are from JINA REA-CLIB (Cyburt et al. 2010) plus additional tabulated weak reaction rates Fuller et al. (1985); Oda et al. (1994); Langanke & Martínez-Pinedo (2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al. (1996).
Our MESA sdB models employ the "predictive" mixing scheme  for convection in the Heburning core. This allows the convective core to grow over the sdB star lifetime without requiring convective overshoot (Schindler et al. 2015) while also providing a setting that prevents the core from growing too far and dividing, which can cause "breathing pulses" that cause the core size to fluctuate along with the star luminosity and radius. Breathing pulses would greatly complicate the study of mass transfer via Roche-lobe overflow from sdB stars, but they are likely numerical artifacts rather than physical phenomena , so we use the predictive mixing scheme so that our models avoid them during core He burning. For a more thorough discussion of convection treatments in MESA sdB models, see Ostrowski et al. (2021).
Our first grid of MESA sdB models represents a canonical 0.47 M sdB star with a range of envelope masses. These models are constructed from a 1 M progenitor star with metallicity Z = 0.02 that undergoes the He core flash when the core reaches a mass of 0.47 M . After the He core flash, we remove all but a thin residual H envelope, and evolve the models as sdB stars to produce the relation between sdB ra- dius R and envelope mass M env shown in Figure 1. Envelope mass M env is defined here as all the mass of material contained in the exterior layers for which the H mass fraction is X > 0.01. For this set of 0.47 M models, the envelope material has roughly solar composition matching the initial ZAMS composition of the progenitor 1 M star. Figure 1 shows three snapshots at different ages along the radius evolution in this grid of models, showing that these sdB models slowly evolve from smaller to larger radii for most of their core He burning lifetime. Figure 1 also shows radius estimates from the models of Han et al. (2002) based on the beginning of the log g-T eff tracks in figure 2 of that paper, confirming that our models generally agree with another stellar evolution code commonly used to model sdB stars. Envelope masses log(M env /M ) < −4.5 in Figure 1 give radii that converge toward the radius of a model with a bare He core and no envelope.
To illustrate the case of a less massive sdB star donor, we also run a second grid of MESA models for a 0.34 M sdB star. This smaller He-burning core is descended from a 2.8 M progenitor star that does not go through the He core flash, and the residual H-rich material that forms the sdB envelope is enriched in He due to partial nuclear processing of that region in the receding convective core on the main sequence (see Kupfer et al. 2020a for more details). Figure 2 shows the M env -R relation for this lower-mass grid of MESA sdB models. Envelope masses log(M env /M ) < −3 for this model yield radii that quickly converge to the radius of a bare He core. This is due to the smoother composition transition between the He core and H envelope, which leaves only a few percent H by mass in the residual outer envelope when log(M env /M ) < −3.

MASS TRANSFER PHYSICS
We now present the physics that describes the mass transfer rate for both the envelope and He-rich core of the sdB star. In Section 5, we will also compare these mass transfer rate predictions to some example MESA binary evolution models, and confirm that they show general agreement with the analytic calculations here.
We label the sdB as star 1 with radius R 1 and mass M 1 , and the WD as star 2 with radius R 2 and mass M 2 . The mass ratio is q ≡ M 1 /M 2 , and following Eggleton (1983), we approximate the Roche radius for the hot subdwarf as where a is the orbital separation between the two stars. Kepler's third law gives the relation between orbital separation and period: When the sdB star begins to overflow its Roche lobe, the mass transfer rate depends on the response of the star's radius to losing mass, which can be described in terms of an exponent ζ for the relation R 1 ∝ M ζ 1 (Soberman et al. 1997). If mass transfer is slow enough that the star can maintain thermal equilibrium, the relevant exponent is which simply encodes the mass-radius relation for a star in equilibrium with no structural effects from losing mass. On the other hand, if thermal transport is not efficient enough to maintain equilibrium on the timescale of mass removal, the entropy profile in the stellar interior will remain fixed, so that the relevant exponent encodes the adiabatic radius response: The overall timescale for orbital evolution is set by the rate at which gravitational-wave radiation (GWR) shrinks the orbital separation and Roche lobe: whereJ gr = − 32 5c 5 2πG P orb and In order to estimate the mass transfer rate, it is also necessary to account for the change in the Roche radius that results from transferring mass from star 1 to star 2: For conservative mass transfer (M 1 + M 2 = constant) and constant J orb , we can use Equation (7) to evaluate the first term on the right hand side as d ln a/d ln M 1 = 2(q − 1). Conservative mass transfer also yields d ln q/d ln M 1 = 1 + q. Differentiating Equation (1) gives So the full expression for the Roche radius response to conservative mass transfer can be written as (10) It will often be convenient to adopt the much simpler approximation which can be shown to agree with Equation (10) to within a few percent. Note that ζ RL < 0 for q 0.8, so in general the Roche lobe radius expands in response to mass lost from the sdB donor for the systems with M WD ≥ 0.6 M C/O WD accretors that are the focus of this paper.
Lastly, we also must account for the timescale for the donor star's radius to change due to stellar evolution independent of mass transfer We assume fully conservative mass transfer here, so thatJ gr is the only term affecting the net angular momentum of the system, and the equilibrium mass transfer rate can then be estimated as (see e.g., Ritter 1988) where ζ is determined by comparing the timescale for orbital evolution to the thermal (Kelvin-Helmholtz) timescale: We use ζ = ζ eq when τ KH < τ gr and ζ = ζ ad when τ KH > τ gr , and mass transfer is stable as long as ζ RL < ζ. Throughout this paper, we will assume that both stars in the binary maintain co-rotation with the orbit, so that we can neglect angular momentum that might be lost from the orbit e.g. due to accretion causing the WD to spin up. Direct impact accretion is thought to violate this assumption in double white dwarf systems (Nelemans et al. 2001;Marsh et al. 2004), causing dynamically unstable mass transfer in systems with q < 2/3 that would otherwise be stable. However, the sdB+WD systems that we consider here have larger orbital separations due to the larger sdB radius when it overflows its Roche lobe. Using the Nelemans et al. (2001) fit to the Lubow & Shu (1975) calculations for minimum distance of the accretion stream from the accretor, we find that sdB+WD binary systems will never come close to the direct impact regime due to the larger orbital separations required for sdB donors. These systems are therefore expected to form an accretion disk that transfers most of its angular momentum back into the orbit (Verbunt & Rappaport 1988;Priedhorsky & Verbunt 1988). Therefore, we will find that mass transfer is stable in almost all of the systems that we consider here, except the case of low mass sdB donors (M sdB 0.35 M ) transferring He in systems with q > 2/3 (see Section 3.2).

Envelope Transfer
To represent a typical subdwarf, we consider an sdB star of total mass M 1 = 0.47 M with an initial hydrogen envelope mass of 10 −3 M . Figure 1 shows that such a subdwarf will have a radius in the range R 1 = 0.16-0.20 R for most of its lifetime. For the example binary evolution scenarios presented throughout this paper, we will adopt a WD companion mass of M 2 = 0.75 M , which is representative of many of the short-period sdB+WD systems that have been observed so far (e.g., Geier et al. 2013;Kupfer et al. 2020b;Kupfer et al. 2021, in prep), likely reflecting the fact that these systems are found among young populations. For this companion mass M 2 = 0.75 M , Equations (1) and (2) together imply that the sdB star will initially fill its Roche lobe (R 1 = R RL ) at an orbital period of 50-70 min. For the sdB stars shown in Figure 1 with luminosities on the order of L ≈ 20 L , the Kelvin-Helmholtz timescale is τ KH 2 Myr. For M 1 = 0.47 M and M 2 = 0.75 M and P orb = 50 min, τ gr ≈ 150 Myr >> τ KH , and the star's radius response is therefore described by ζ eq . From the timescale for radius evolution shown in Figure 1, we can estimate τ R ≈ 400 Myr based on the change in ln R 1 according to Equation (12).
For mass transfer of the envelope, we can write the star's radius response in terms of an exponent for the envelope mass rather than the total stellar mass: which can be estimated from Figure 1, giving values of ζ eq,env ≈ 0.2-0.3. This implies that ζ eq 1 because M 1 M env . On the other hand, ζ RL is of order unity according to Equation (11) (cf. figure 4 in Soberman et al. 1997), so ζ RL can be neglected when applying Equation (13) to transferring the envelope. We can therefore approximate the instantaneous mass transfer rate for the hydrogen envelope aṡ which predicts that when the subdwarf first fills its Roche lobe at P orb = 50-70 min, mass transfer will occur at rates in the range ofṀ env ≈ 2-8 × 10 −11 M yr −1 . This mass transfer is driven primarily by the shrinking of the orbit due to GWR (τ gr term), with a smaller contribution from the slowly expanding radius of the star pushing material past its Roche lobe (τ R term). Without a hydrogen envelope, Figure 1 shows that the bare helium core has a substantially smaller equilibrium radius on the order of 0.11-0.13 R , which fills the Roche lobe at an orbital period of 28-36 min (see also Figure 3). We can write the GWR merger time of two point masses as .
This equation describes the time evolution of the system's orbital separation before the stars are in contact (M 1 and M 2 constant), and it can be shown that it also applies as long as the stars are transferring mass at a rate that is negligible for affecting the angular momentum of each star compared to GWR evolution. In this context, a rate ofṀ 10 −9 M yr −1 is negligible, so Equation (17) applies during the envelope transfer phase as well as prior to the stars coming into contact. Using this equation along with Equation (2), we find that the inspiral time to reach a period of 30 min from a period of 50 min (70 min) is 14 Myr (40 Myr). Therefore, after the hydrogen envelope begins to overflow the Roche lobe, the mass transfer must occur at an average rate on the order of 2-7 × 10 −11 M yr −1 until the 10 −3 M envelope is exhausted, which agrees with the prediction of Equation (16). A less massive sdB star (like the 0.34 M models shown in Figure 2) will not fill its Roche lobe until reaching an orbital period of P orb = 25-40 min (for a 0.75 M WD companion), and the period at which the bare He core is exposed is P orb = 18 min. The inspiral timescale according to Equation (17) is therefore only 2-10 Myr, and mass transfer of the envelope must occur at a rate on the order ofṀ env ∼ 10 −10 -10 −9 M yr −1 .
Once the hydrogen envelope is exhausted after a few to 10s of Myr, the He-rich layers underneath will be revealed and mass transfer will proceed at a much higher rate.

The Helium Core
By the time the hydrogen envelope is completely removed at an orbital period on the order of P orb ≈ 30 min for the canonical 0.47 M sdB model, the timescale for GWR to shrink the orbit is much shorter at τ gr ≈ 38 Myr. Figure 3 shows the mass-radius relation for He cores based on a grid of MESA models with no H envelope. This relation gives ζ eq ≈ 2 for He stars with no H envelope, and therefore Equations (11) and (13) giveṀ 1 ≈ 10 −8 M yr −1 when He mass transfer begins. The second panel of Figure 3 shows the periods at which sdB He cores will fill the Roche lobe with a 0.75 M companion, calculated with Equations (1) and (2) by setting R RL = R using the radius from the top panel. Figure 3 also shows the mass-luminosity relation for these cores, which allows calculating τ KH . The much lower luminosity and radius of sdB stars with M 1 0.35 M results in the orbital evolution timescale becoming shorter than the Kelvin-Helmholtz timescale, so sdB stars that lose mass beyond this point can no longer maintain thermal equilibrium and will respond to mass loss adiabatically. The index to describe the radius response for M 1 0.35 M is therefore ζ ad rather than ζ eq . We can approximate this index as ζ ad = −1/3 for an ideal gas polytrope (Soberman et al. 1997), meaning that the adiabatic response of the star is to expand slightly as it loses mass. Stable mass transfer therefore requires that ζ RL < −1/3, or q 2/3 according to Equation (11), but this only requires that the WD companion have mass M 2 > 0.525 M , so this condition will generally be satisfied. Lower mass sdB stars, or those that lose enough mass to enter this adiabatic regime, also have shorter orbital periods (2nd panel in Figure 3) corresponding to τ gr as short as 5-10 Myr (4th panel). With ζ ad −ζ RL ≈ 0.5, Equation (13)   fore predicts mass transfer rates of orderṀ 1 10 −7 M yr −1 in this lower mass regime.
The general picture for He transfer then is that a canonical 0.47 M sdB star should begin transferring He atṀ 1 ≈ 10 −8 M yr −1 , and the mass transfer rate will gradually increase up to a value approaching 10 −7 M yr −1 as the sdB star loses mass, which agrees with previous work modeling mass transfer from He stars (Savonije et al. 1986;Iben & Tutukov 1991;Yungelson 2008;Brooks et al. 2015).

POST COMMON-ENVELOPE PERIODS AND CONTACT SCENARIOS
In this section, we show the range of post-common envelope binary (PCEB) orbital periods that can lead to mass transfer from a hot subdwarf to a white dwarf companion. After a PCEB exits the common envelope at period P init , it spirals inward toward shorter orbital periods. If the GWR inspiral timescale at P init is shorter than the He-burning lifetime of the sdB star, the binary will eventually make contact and transfer mass from the sdB to the WD. On the other hand, if the orbital period is long enough that the inspiral time is longer than the He-burning lifetime, the sdB star will exhaust its nuclear fuel, and the binary will merge later as a double WD system that includes a "hybrid" He/C/O WD (Zenati et al. 2019;Schwab & Bauer 2021). Using models for sdB radius and lifetime, we can find the maximum value of P init for which a system can make contact as an sdB+WD. In order to map out this range of orbital periods, we start with a grid of MESA models run from ZAMS to construct a variety of He-burning core masses. Figure 4 shows the relation between zero-age main-sequence (ZAMS) mass M ZAMS and the resulting He core mass M core for models with no overshoot and initial metallicity Z = 0.02. The M core -M ZAMS relation is sensitive at the ≈10% level to parameters like initial metallicity and overshoot (through changes to opacities, convection, and mixing), but our grid here is sufficient for our purposes in that it produces a representative range of He-burning core masses for sdB models. A more detailed investigation of the M core -M ZAMS relation for models including a range of metallicities and overshoot treatments can be found in Ostrowski et al. (2021). After the models ignite He, we strip the outer H envelope down to a specified total H mass (either 10 −3 M or 3 × 10 −3 M ) and evolve the model as an sdB star through its core He-burning lifetime. Note that the total envelope mass can vary depending on how much He is contained in the Hrich layers above the He core, and in general M env > M H . For Orange shaded regions show orbital period ranges over which the H envelope can transfer to the WD companion, calculated using the radii from the upper panels. The blue shaded region shows the orbital period range where the orbit can become compact enough to transfer material from the He core of the sdB star. The gray shaded region to the right shows orbital periods for which the inspiral time is too long for the system to transfer any mass within the sdB star's He burning lifetime. The red shaded region shows the range of orbital periods for which the system may be able to inspiral close enough to transfer some of the H envelope, but will not have time to become compact enough for mass transfer to reach the He core. the stars descended from progenitors with M ZAMS 2.3 M that go through the He core flash, there is a sharp boundary between the He core and the H envelope, so that the envelope has roughly solar composition and its total mass is M env ≈ M H /0.7. On the other hand, sdB stars descended from progenitors with M ZAMS 2.3 M that ignite He in the center have envelope material that has been partially burned as part of the receding convective core earlier on the main sequence, so that the mass fraction of hydrogen in the region of the star that becomes the sdB envelope is X < 0.7, and it can be as low as X ∼ 0.2 (Kupfer et al. 2020a). This envelope composition affects the compactness of the envelope in our sdB models and hence the overall compactness of the sdB stars. For example, a 0.47 M sdB star can descend either from a MS star of mass M ZAMS = 1.0-1.5 M or from a higher mass MS star with M ZAMS ≈ 3.8 M . For our sdB models, the star with the higher mass progenitor has a much more He-rich envelope and a more compact overall structure, as shown in radius evolution displayed in the upper panels of Figure 5. Using these sdB models, we then calculate the maximum PCEB orbital period P init for which an sdB binary system will be able to spiral inward and transfer mass during its core He burning lifetime. For each model, we take the maximum radius during core He burning (identified by the black points in the upper panels of Figure 5), and we calculate the orbital period at which it will fill its Roche lobe using Equations (1) and (2), assuming a WD companion mass of 0.75 M . We then use the age for the maximum radius of each model in the upper panels of Figure 5 to calculate the maximum initial orbital period for which the system will have time to spiral inward and reach contact according to Equation (17). We use a similar grid of bare He models with no H envelope to identify the analogous contact and initial periods for mass transfer from the more compact He cores.
The shaded regions in the lower panels of Figure 5 map out the orbital period ranges over which H and He phases of mass transfer can occur, as well as the range of P init for which an sdB+WD system will not be able to spiral into contact with the sdB star's He core burning lifetime. For example, a 0.47 M sdB star with a sufficiently thick H-rich envelope will be able to fill the Roche lobe and transfer some hydrogen for P init as long as 175 min, but it can only spiral in far enough to transfer material from the He core for P init < 105 min. It will first fill its Roche lobe and begin transferring the H envelope at an orbital period in the range P orb = 50-150 min. It will eventually exhaust the envelope and begin transferring material from the He core at P orb < 35 min. The WD companion mass can also have a minor impact on the orbital period ranges displayed in these plots. See the Appendix for more discussion on this effect.

MESA BINARY EVOLUTION MODELS
To illustrate typical sdB+WD binary evolution scenarios, we run two example MESA binary models through the mass transfer phases up to the point of He thermonuclear runaway on the accreting WD. These models calculate the mass transfer rate according the implicit Roche lobe overflow scheme described in Paxton et al. (2015) for the Ritter (1988) mass transfer prescription, and they therefore provide independent verification for the analytic accretion rate estimates from Section 3.
The angular momentum in these binary evolution models is driven by GWR losses as described by Equations (6) and (7). We evolve both the sdB donor star and a 0.75 M WD accretor, ignoring rotation in both stars for simplicity. Our modeling for the accreting WD closely follows Bauer et al. (2017), including the 14 N(e − , ν) 14 C(α, γ) 18 O (NCO) reaction chain that may initiate the eventual 3α thermonuclear runaway depending on the accretion rate. The 14 C(α, γ) 18 O rate in our MESA models is that of Johnson et al. (2009) as in Bauer et al. (2017). Neunteufel et al. (2017Neunteufel et al. ( , 2019 have argued that rapid rotation and magnetic fields due to the angular momentum from the accreted material on the WD may significantly alter the thermal structure of the WD's accreted envelope and modify the eventual thermonuclear runaway in the accreted He. However, we interpret the calculations of Fuller & Lai (2014) as indicating that dynamical tides and dissipation due to WD pulsation modes should prevent the WD from spinning up to rotation rates near critical. Instead, tides are likely to drive the WD toward much slower rotation rates on the order of P orb , which would be far too slow to effect the envelope structure of a compact WD with dynamical time t dyn P orb . This is also consistent with the fact that observations rule out rapid rotation for the accreting WD in some AM CVn systems (Roelofs et al. 2006;Kupfer et al. 2016). Therefore, it should be safe to ignore rotation in calculating the structure of the He envelope accreted onto the WD.
As the sdB spirals inward toward its companion, tides may cause it to spin up to rotation on the order of the orbital period, though tidal dissipation may not be efficient enough to bring the donor star into fully synchronous rotation (Preece et al. 2018). Even for fully synchronous rotation, we can use Equation (2) to find that a Roche-filling donor would rotate with a frequency Ω 2 orb /Ω 2 crit = (1 + 1/q)(R RL /a) 3 , where Ω crit = GM 1 /R 3 1 is the critical rotation rate of the donor star and Ω orb = 2π/P orb . Using Equation (1) along with the above formula, we find that Ω 2 orb /Ω 2 crit ≈ 0.1 for any q 1, so distortion in the outer layers due to rotation would have only a minor impact on the donor sdB star. Since this is an upper limit on the rotation rate due to potential tidal spin-up, we feel it is justified to neglect the structural effects of rotation for the MESA sdB models in this section.
We assume conservative mass transfer (Ṁ 2 = −Ṁ 1 ) unless the WD is undergoing a nova outburst that causes it to expand and overflow its Roche lobe. In that case, we adopt the following simple prescription for removing mass from the WD to allow it to evolve through nova cycles while staying within its Roche lobe and allowing binary evolution to continue. When the WD radius expands beyond 0.85 R RL,2 (where R RL,2 is the WD Roche radius), we use a simple exponential Roche lobe overflow scheme for the WD to remove mass from system at a rate oḟ This mass is assumed to be lost from the WD as a spherical wind, and therefore carries away specific angular momentum matching the WD orbit. It should be noted that the parameters in Equation (18) are chosen for simplicity rather than based on any physical motivation. These parameters were chosen by experimenting to find values that allowed the WD to expand near its Roche lobe and then quickly ramp up mass loss to allow robustly evolving through multiple complete nova cycles, but they are not necessarily unique in accomplishing this goal. The details of this WD mass loss scheme therefore should not be construed as physical, but the small amount of mass lost from the system during these nova cycles will not have a large impact on the overall orbital evolution. For a more thorough discussion of the impact of mass loss prescriptions on MESA models of novae, see Wolf et al. (2013). For simplicity, we also do not consider any interaction between the donor and mass lost from the WD during these hydrogen novae, which may lead to some dynamical friction that enhances the rate at which the system loses angular momentum and spirals inward, though this process is uncertain (e.g., Shen 2015).

Canonical Mass sdB Donor
Our canonical sdB model is 0.47 M with an envelope mass of M env = 10 −3 M , descended from a 1.0 M model that has the H envelope removed during the He core flash. We initialize this in a binary with orbital period P init = 90 min with a 0.75 M WD companion. The left panel of Figure 6 shows that this system comes into contact after 50 Myr when P orb = 65 min, and donates the H envelope at a rate oḟ M 10 −10 M yr −1 for 30 Myr. During this time, the WD accretor model undergoes 6 H novae while accreting H-rich material, and ejects most of the 10 −3 M of accreted material from the system. The hot subdwarf donor increases its T eff to nearly 35,000 K as the Roche lobe shrinks and the hot subdwarf becomes more compact while maintaining a nearly constant luminosity before finally exposing the He core.
Once the He core is exposed and begins to overflow the Roche lobe when P orb = 34 min, the mass transfer rate in-creases toṀ ∼ 10 −8 M yr −1 for another 20 Myr, in agreement with the prediction of Section 3. The sdB donates nearly 0.2 M of He-rich material to the WD companion, and its luminosity and T eff decrease as core He burning subsides and eventually shuts off altogether in the donor star, while the remaining He core continues to grow more compact. Eventually, the WD accumulates enough He that the NCO chain initiates a dynamical thermonuclear runaway at the base of the accreted He envelope, likely leading to a detonation. The density at the ignition location in the envelope of this model is ρ = 1.7 × 10 6 g cm −3 , above the critical density threshold ρ crit = 10 6 g cm −3 for which a detonation is likely (Woosley & Weaver 1994;Woosley & Kasen 2011;Bauer et al. 2017;Neunteufel et al. 2017). This detonation may destroy the entire WD and liberate the remnant of the sdB donor to become a runaway star with a velocity close to its final orbital velocity of 780 km s −1 (Neunteufel et al. 2019;Bauer et al. 2019;Neunteufel 2020).
The fourth panel in Figure 6 gives bolometric luminosities of both the sdB donor and WD accretor as well as the estimated accretion luminosity L acc ≈ GM WDṀ /R WD . This shows that the sdB star dominates the luminosity of the system before mass transfer and during the envelope transfer phase (except during novae), but the accretion luminosity may become comparable or larger during the He transfer phase. Late in the evolution of these systems, the accretion disk may dominate observed luminosity in both the ultraviolet and optical. In this phase, these binaries would appear as AM Canum Venaticorum (AM CVn) systems, with negligible luminosity from the sdB donor once it loses enough mass to drop below ≈0.3 M . This particular model donates some carbon along with the He during this phase due to most of the He core containing some ashes from the He core flash, and the spectral signatures of this C in the accretion disk would be atypical for an AM CVn. Donors descended from higher mass progenitors that ignite He in their centers rather than through an off-center core flash do not contain any C ashes outside the convective He-burning core. These sdB stars contain a ≈0.2-0.3 M shell of unprocessed Hedominated material in their outer layers to transfer during the He overflow phase, and this would form a more typical AM CVn spectrum from the accretion disk, showing nitrogen features from the ashes of main sequence CNO burning, but not carbon features. See also Yungelson (2008) for discussion of the composition of material donated from He stars corresponding to the scenario of higher mass progenitors that did not experience an off-center He flash.

0.37 M sdB Donor
Our second binary evolution model contains a 0.37 M sdB star descended from a M ZAMS = 3.0 M progenitor (cf. Figure 4). The envelope contains a total hydrogen mass of . MESA binary evolution models of sdB stars with a 0.75 M WD companion. The left panels show the model for a canonical 0.47 M sdB star described in Section 5.1, and the right panels show the model for a 0.37 M model described in Section 5.2. The x-axis in these plots is time since He core burning started in the sdB model, which also corresponds to when the envelope was removed from the progenitor star to form a subdwarf, presumably through a common envelope that also leaves the binary systems at the initial orbital periods Pinit labeled at the top of the plots.
M H = 10 −3 M like the models in the left panel of Figure 5. However, the H mass fraction in the layers that form the eventual sdB envelope is only X ≈ 0.2 due to prior nuclear processing in the receding convective core during the MS evolution. This leads to a very compact overall structure for the sdB model because of the higher He content in the envelope, and the total mass of material in the outer layers that contain some H is M env = 7 × 10 −3 M . We initialize this sdB model in a binary with a 0.75 M WD with an orbital period of P init = 120 min. It first comes into contact after 220 Myr when the orbital period has decreased to P orb = 37 min, and it then begins donating the envelope material at a rate oḟ M ∼ 10 −9 M yr −1 . Over the next 7 Myr, the WD accretor undergoes 22 H novae and ejects much of the accreted envelope mass from the system. The downward spikes in the accretion rate in the right panel of Figure 6 are due to temporary widening of the orbit due to this small amount of mass loss during each nova cycle. These downward spikes are a direct result of our oversimplified prescription for removing mass from the WD during novae as a simple spherical wind, so the orbital evolution briefly following each nova is likely oversimplified and unphysical as well. However, this has little impact on the overall orbital evolution and mass transfer, as each nova removes only ∼ 10 −4 M from the system, and GWR quickly brings the system back into contact with a similar equilibrium mass transfer rate to that before the nova occurred.
Eventually when P orb = 26 min the underlying He core is exposed, and mass then begins to transfer to the WD at a rate on the order of 10 −8 M yr −1 for the next 9 Myr. The system transfers a total of 0.15 M before 3α burning in the degenerate accreted He leads to a thermonuclear runaway in the WD envelope. The final donor velocity in this model is 950 km s −1 , but the ignition location of the WD envelope occurs in less dense layers (ρ = 4 × 10 5 g cm −3 ) outward from the base of the accreted He due to a temperature inversion (Brooks et al. 2015;Bauer et al. 2017), so this particular model is unlikely to lead to a detonation that liberates the donor as a runaway star. If the WD accretor does not detonate, this system could continue transferring mass and evolve through a period minimum around P orb = 10 min when the donor is ≈0.2 M (Yungelson 2008; Neunteufel 2020), eventually transferring hybrid He/C/O material from the partially burned core as the system evolves back toward longer orbital periods (Nelemans 2010).

DETECTABILITY OF ELLIPSOIDAL VARIATION
Ellipsoidal modulation of the lightcurves in sdB+WD binary systems is the key observable that will enable discovery of new systems based on time-domain observations with enough epochs or short enough cadence to resolve orbital periods on the order of 20 minutes to 3 hours. In order to demonstrate the detectability of this modulation, in this section we use LCURVE 2 ) to calculate the amplitude of variation by modeling the lightcurve variability of ellipsoidally deformed sdB stars in binary systems with WD companions. We model phase folded lightcurves with 1000 data points. We use values for limb darkening and gravity darkening based on Claret (2017) and assume a typical sdB star with T eff = 25,000 K and log(g/cm s −2 ) = 5.5. We treat the WD companion as a point mass, so that there are no eclipses as we are only interested in the ellipsoidal amplitude. The geometry of the ellipsoidally deformed sdB star is the primary source of flux variation. Limb darkening and gravity darkening have small influence on the ellipsoidal amplitude in our models (order 1% level) but can increase towards hotter temperatures. Doppler beaming can also slightly enhance the overall amplitude by boosting one of the lightcurve peaks relative to the other peak (e.g., Bloemen et al. 2011), as seen in Pelisoli et al. (2021) with TESS observations of HD 265435. For actively accreting systems, an accretion disk can further alter the lightcurve through the Rossiter-McLaughlin effect or due to irradiation, as seen in Kupfer et al. (2020a). Our models in this section therefore represent a conservative estimate of the total amplitude of  Contours showing the observable amplitude of flux variability as a function of orbital inclination and Roche-filling fraction for two different mass ratios: q = 0.5 (upper panel) and q = 1.0 (lower panel). Variability amplitudes are calculated using LCURVE ) assuming a typical sdB star with Teff = 25,000 K and log(g/cm s −2 ) = 5.5. variation in an ellipsoidal system, as the amplitude can be further enhanced by eclipse features from a WD or disk at high inclination or greater limb and gravity darkening for a subdwarf with T eff > 25,000 K.
We find that the main factors influencing the amplitude of ellipsoidal modulation are the Roche-filling fraction of the sdB (R sdB /R RL ) and inclination (see equation 4 in Bloemen et al. 2012). Figure 7 therefore shows the amplitude of ellipsoidal modulation as a function of these two variables. The mass ratio can also have a minor impact on the amplitude, so we show the amplitude for two different mass ratios (q = 0.5 and q = 1.0).
As an example to illustrate the detectability of ellipsoidal variation in practice, for sdB stars within 2 kpc the variability amplitudes in Figure 7 translate into the following sensitivity to detection by the Zwicky Transient Facility (ZTF, Bellm et al. 2019;Graham et al. 2019), which was used to discover the systems of Kupfer et al. (2020a,b). Assuming an sdB star with absolute magnitude of 4 (typical in the Gaia G band) yields an apparent magnitude of 15.5 at 2 kpc without any extinction. Stars in the Galactic plane may experience an additional 0.7-1.0 mag of extinction per kpc (Krelowski & Papaj 1993;Green et al. 2019), yielding 17-17.5 mag for stars at 2 kpc. For stars at this magnitude, ZTF has a flux precision of ∆F/F ∼ 1.5 − 2% (see figure 11 in Masci et al. 2019). Therefore, a photometric amplitude of 5-6% can be easily detected with a > 3σ significance. The contours in Figure 7 therefore show that all Roche-filling sdB stars with inclination i 30 • should be easily detectable at this distance, and a significant fraction of sdB stars should be detectable even when they have only reached 70-80% of Roche-lobe filling.
In the future, the Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) of the Vera Rubin Observatory will enable pushing this detection sensitivity out to even further distances (several kpc) and smaller Roche-filling fractions for sdB+WD systems. A magnitude of 17 is the bright end for LSST, and the photometric precision goal is 5 mmag in g and r (0.5%; see Table 14 in the requirements document 3 ). An amplitude of just 1.5-2% will therefore be easily detected with > 3σ significance, enabling detection of some sdB systems that are only ≈50% of Roche-filling according to Figure 7. For a typical compact sdB+WD system that is born by exiting a common envelope at P orb ≈ 1.5 − 3 hours, the sdB could spend up to 100s of Myr with detectable ellipsoidal deformation as it spirals inward before filling its Roche lobe and beginning to transfer mass for another few 10s of Myr. Detection of the distribution of sdB+WD binaries with ellipsoidal modulation in LSST could therefore map out nearly the entire lifetime of this phase, from birth to eventual mass transfer, supernova, or merger.

DISCUSSION AND CONCLUSIONS
We have calculated characteristic timescales and orbital periods for mass transfer from sdB stars to WD companions. Our results imply that many undiscovered compact sdB+WD binary systems with Roche-filling or near Rochefilling sdB donors have yet to be discovered. We can estimate a lower limit for the number of observable systems based on the following argument. Kupfer et al. (2020a,b) found two short-period (39 and 56 min) subdwarf+WD binaries within 2 kpc in which the subdwarf is currently overflowing its Roche lobe and transferring mass to the WD. These systems were discovered because the subdwarf donors are very hot (T eff > 33,000 K), which allowed discovery due to ellipsoidal modulation in a sample with strict color cuts even after reddening due to dust in the Galactic plane where these systems reside (Kupfer et al. 2020b).
Modeling in Kupfer et al. (2020a,b) indicated that these observed systems are currently evolving through a short-lived ≈1 Myr phase in which the subdwarf becomes hotter just after core He burning has ended. Similar systems that exited a common envelope at slightly shorter orbital periods would come into contact Myr earlier while the donor is still a cooler He burning sdB star. Our calculations in this paper indicate that the mass transfer phases for both the outer H envelope and the He core last for ∼10 Myr when the binary comes into contact while the sdB star is still burning He in its core. Section 6 shows that the variability of these systems will be easily detectable by ZTF for all Roche-filling systems within 2 kpc except those with small inclinations i 30 • . Therefore, at least 10s of binary systems with Roche-filling donors in both of these phases should be observable, and there are likely 100s more that are close enough to Roche-filling to exhibit detectable ellipsoidal modulation according to Figure 7. Clearly, many more sdB+WD systems await discovery among the millions of variable stars observed in the Galactic plane (Kupfer et al. 2021).
Most of these systems have not yet been discovered because they likely reside in the Galactic plane, where the strict color cuts of previous searches for sdB stars have resulted in excluding these systems due to reddening. The demographics of the few systems that we do know of so far (Vennes et al. 2012;Geier et al. 2013;Kupfer et al. 2017aKupfer et al. ,b, 2020a indicate that they come from a young population in the Galactic plane. Indeed, measured WD companion masses are often significantly larger than 0.6 M , indicating short MS progenitor lifetimes prior to the common envelope. The population synthesis models of Han et al. (2003) also indicate that sdB+WD binaries with P orb < 3 hr primarily descend from stars with M ZAMS > 2 M that do not experience a degenerate core He flash, another indicator that these systems must descend from a young population.
Additionally, Geier et al. (2019) present a catalog of ≈40,000 sdB candidates from Gaia DR2, but their selection criteria are likely to exclude many sdB stars in the Galactic plane associated with a young population in dense stellar regions. Their selection criteria include quality metrics from Gaia DR2 which remove blended sources and sources with bad astrometry. This affects in particular dense stellar regions, which makes the selection based on the Geier et al. (2019) catalog very incomplete at low Galactic latitudes with high stellar densities. Therefore, searches based on this catalog will miss a significant number of sdBs from a young population. sdB stars descended from such a young population would have the more compact envelope configurations of stars that do not experience the He core flash in Figure 5, and therefore they should typically experience the onset of H envelope mass transfer at orbital periods of 40-60 minutes according to the darker orange shaded region in the lower panels of Figure 5. We therefore predict that this period range should be the most fruitful in searching to discover new Roche-filling sdB+WD binary systems.
Discovering and characterizing this population of sdB+WD binary systems will provide useful constraints on common envelope outcomes. The lower panels of Figure 5 show that systems which spiral close enough to transfer mass from the sdB star must exit the common envelope at orbital periods shorter than the shaded regions on the right side of those plots. Modeling of well-characterized systems will therefore provide a direct link to post-common-envelope outcomes for this class of binary systems.
Extending our estimates to the total number of systems in the Galaxy based on the volume in which observed systems have been detected so far, we estimate that our Galaxy contains at least ∼ 10 3 -10 4 sdB+WD binaries with orbital periods shorter than 1-2 hours. This lower limit is based on the two systems of Kupfer et al. (2020a,b) that were discovered within a distance of d ≈ 2 kpc using ZTF in the Northern Hemisphere. We can then scale the 10s to 100s of observable systems that we estimate should exist within a few kpc by a factor of M /[(π/2)d 2 Σ d ] ∼ 100 to arrive at the total number of systems in the Galaxy, where Σ d is the local stellar surface density and M is the total Galactic stellar mass (McGaugh 2016). These sdB+WD systems will therefore be a significant source of Galactic foreground gravitational waves for LISA, and some may be resolvable LISA sources. Our estimates of the stars in this population are closely related to the work of Götberg et al. (2020), who focused on He stars in binaries with slightly longer orbital periods P orb ≥ 60 min just outside the range of where their He star models may begin to fill their Roche lobes and transfer mass. Our models therefore represent a subset of the population of ≈ 10 5 such systems that they calculate, with the evolution extended toward shorter orbital periods to continue through mass transfer.
Many of the sdB+WD systems with P orb 2 hr are likely to produce an eventual thermonuclear explosion of a thick (∼0.1 M ) He shell on the accreting WD. With typical inspiral times on the order of 10 8 yr, our estimate of 10 3 -10 4 systems therefore predicts that the rate of such explosions in our Galaxy could be as high as one per 10 4 yr, comparable to the rate of ".Ia" supernovae predicted from double WD AM CVn systems (Bildsten et al. 2007;Shen & Bildsten 2009).
We have also found that during the He transfer phase, the accretion luminosity is often the dominant source of luminosity in the system, and therefore many of these systems may spend several Myr in a state that appears very similar to the properties of AM CVn systems (Nelemans et al. 2001Marsh et al. 2004). For sdB donors descended from M ZAMS 2.3 M progenitors that do not experience the He core flash, the outer ≈0.2 M of He-dominated material that transfers onto the WD has not undergone any He-burning to process the CNO abundances, and so the composition of accretion disks in these systems would be indistinguishable from the He WD donor channel for AM CVns ). Therefore we argue that it is possible that a significant number of relatively short-period (10 min < P orb < 30 min) AM CVn-like systems associated with a young stellar population in the Galactic plane may be descended from sdB+WD binaries, even though the donor composition may point towards CNO cycle material that has traditionally been interpreted as indicating a He WD donor.
In fact, one such AM CVn-like system may already be known. SDSS J190817.07+394036.4 is an AM CVn system (Fontaine et al. 2011) with an orbital period of 18 minutes, and Kupfer et al. (2015) found from the relation between the superhump and the orbital period that this system likely has a mass ratio of q = 0.33, which is problematic for a double WD interpretation of this system. According to Marsh et al. (2004), a double WD system with this mass ratio would likely experience unstable mass transfer and merge, while we have found that q < 2/3 sdB+WD systems will experience stable mass transfer. Our binary model in Section 5.2 evolves through stable He mass transfer with a mass ratio and orbital period very similar to those observed for SDSS J1908.
In principle, these systems would also be distinguishable from the classic He WD AM CVn scenario by the direction of their orbital period evolution. The sdB+WD systems we discuss in this paper haveṖ orb < 0, evolving toward shorter orbital periods, while AM CVn systems with He WD donors spend most of their lifetime evolving in the opposite direction toward longer orbital periods. A detection of negativeṖ orb could therefore serve as another indicator that an AM CVn system comes from the sdB+WD binary channel. We expect orbital period changes on the order ofṖ orb = −(3/8)P orb /τ merge when negligible mass is being transferred so that we can use Equations (2) and (17) directly to calculate a purely GWRdriven value forṖ orb . This equation gives orbital period evolution rates of |Ṗ orb | 10 −12 -10 −11 s s −1 for sdB+WD binary systems. He mass transfer from the less massive sdB to the more massive WD slows orbital evolution to make |Ṗ orb | smaller relative to purely GWR-driven binary evolution (cf. Equation (7) and the top panels in Figure 6). In practice, the magnitude of this orbital period evolution may be detectable after a few years of regular monitoring, at least for some eclipsing systems with precise timing.
We thank the anonymous referee for an insightful report that led to many valuable improvements to this manuscript, including much of the material in Section 6 in particular. We thank Tom Marsh for valuable discussions that originated ideas about sdB hydrogen envelope transfer that grew into this work. We thank Josiah Schwab and the other organizers of The Beginnings and Ends of Double White Dwarfs hosted at the Niels Bohr Institute in Copenhagen during July 2019, where those discussions were stimulated. This work also benefited from the subsequent workshop held at DARK in July 2019 that was funded by the Danish National Research Foundation (DNRF132). We thank Warren Brown, Lars Bildsten, Tom Marsh, and Patrick Neunteufel for comments and suggestions that improved this manuscript as it was in preparation. This research benefited from interactions that were funded by the Gordon and Betty Moore Foundation through grant GBMF5076. This work was supported by the National Science Foundation through grants PHY-1748958 and ACI-1663688.