# Mass Transfer and Stellar Evolution of the White Dwarfs in AM CVn Binaries

Tin Long Sunny Wong Department of Physics, University of California, Santa
Barbara, CA 93106, USA  Lars Bildsten Department of Physics, University of
California, Santa Barbara, CA 93106, USA  Kavli Institute for Theoretical
Physics, University of California, Santa Barbara, CA 93106, USA Tin Long Sunny
Wong<EMAIL_ADDRESS>
###### Abstract

We calculate the stellar evolution of both white dwarfs (WDs) in AM CVn
binaries with orbital periods of $P_{\mathrm{orb}}\approx 5-70$ minutes. We
focus on the cases where the donor starts as a
$M_{\mathrm{He}}<0.2\,\mathrm{M_{\odot}}$ Helium WD and the accretor is a
$M_{\mathrm{WD}}>0.6\,\mathrm{M_{\odot}}$ WD. Using Modules for Experiments in
Stellar Astrophysics (MESA), we simultaneously evolve both WDs assuming
conservative mass transfer and angular momentum loss from gravitational
radiation. This self-consistent evolution yields the important feedback of the
properties of the donor on the mass transfer rate, $\dot{M}$, as well as the
thermal evolution of the accreting WD. Consistent with earlier work, we find
that the high $\dot{M}$’s at early times forces an adiabatic evolution of the
donor for $P_{\mathrm{orb}}<30$ minutes so that its mass-radius relation
depends primarily on its initial entropy. As the donor reaches
$M_{\mathrm{He}}\approx 0.02-0.03\,\mathrm{M_{\odot}}$ at
$P_{\mathrm{orb}}\simeq 30$ minutes, it becomes fully convective and could
lose entropy and expand much less than expected under further mass loss.
However, we show that the lack of reliable opacities for the donor’s surface
inhibit a secure prediction for this possible cooling. Our calculations
capture the core heating that occurs during the first $\approx 10^{7}$ years
of accretion and continue the evolution into the phase of WD cooling that
follows. When compared to existing data for accreting WDs, as seen by Cheng
and collaborators for isolated WDs, we also find that the accreting WDs are
not as cool as we would expect given the amount of time they have had to cool.

††software: MESA (v12778; Paxton et al., 2011, 2013, 2015, 2018, 2019),
ipython/jupyter (Pérez & Granger, 2007; Kluyver et al., 2016), matplotlib
(Hunter, 2007), NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020),
Astropy (Astropy Collaboration et al., 2013, 2018), Python from python.org and
WebPlotDigitizer (Rohatgi, 2021) .

,

##  1 Introduction

AM Canum Venaticorum (AM CVn) systems are ultracompact binaries undergoing
helium mass-transfer with orbital periods, $P_{\mathrm{orb}}$, between 5 and
68 minutes (e.g., Nather et al., 1981; Nelemans et al., 2001, 2004; Ramsay et
al., 2018). Their orbital evolution is dominated by loss of angular momentum
from gravitational wave radiation (e.g., Tutukov & Yungelson, 1979, 1981), and
the gravitational wave signal from AM CVn systems should be detectable by
missions such as the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et
al., 2013, 2017; Nelemans et al., 2004; Kremer et al., 2017; Breivik et al.,
2018). The study of AM CVn systems is of interest for several reasons. Their
orbital evolution can help constrain the efficiency of tides in degenerate
stars (e.g., Piro, 2019). Unstable helium burning is expected to occur in some
AM CVn systems, leading to helium novae (e.g., Ashok & Banerjee, 2003), and in
the case of dynamical helium-burning, faint thermonuclear “.Ia” supernovae
(Bildsten et al., 2007; Shen & Bildsten, 2009). Accretion is expected to occur
via an accretion disk (except at shortest periods where direct-impact
accretion occurs; e.g., Marsh et al., 2004), which allows AM CVn systems to be
good testing grounds for the theory of nearly pure helium accretion disks
(e.g., Kotko et al., 2012).

There are three possible channels for the formation of AM CVn systems. In the
helium white dwarf (He WD) donor scenario, the initial donor is a degenerate,
low-mass ($\approx 0.1-0.2\,\mathrm{M_{\odot}}$) helium white dwarf (e.g.,
Deloye et al., 2007). In the helium star (He star) donor scenario, the donor
starts mass-transfer as a non-degenerate, helium-burning star with masses
ranging from $\approx 0.3$ to $0.7\,\mathrm{M_{\odot}}$ (e.g., Iben & Tutukov,
1991; Yungelson, 2008). Note that at long periods $P_{\mathrm{orb}}\gtrsim
40\,\mathrm{min}$, the thermal properties of the He WD and He star donors are
predicted to converge (Deloye et al., 2007; Yungelson, 2008). The two remain
distinguishable by their compositions, with the He star donors expected to
contain products of helium-burning. In the cataclysmic variable (CV) donor
channel (e.g., Podsiadlowski et al., 2003), the donor starts as a $\approx
1\,\mathrm{M_{\odot}}$ star, initiates hydrogen-rich mass transfer around the
end of core hydrogen-burning, and eventually loses its hydrogen-rich envelope.
In all three channels, the accretor is often assumed to be a
$0.6-1.0\,\mathrm{M_{\odot}}$ carbon-oxygen white dwarf (CO WD).

In this work, we revisit the calculations by Bildsten et al. (2006) on the
thermal evolution of the CO WD accretor in the context of the He WD donor
channel. In Section 2, we describe our computational setup and construction of
initial stellar models for both the donor and the accretor. We evolve donor
models of various initial central entropies under mass transfer to a point-
mass, and describe their thermal evolution in Section 3. We show that the
mass-radius relation of the donor, as well as the time evolution of
$P_{\mathrm{orb}}$, depends on their initial central entropy and whether they
can cool. In Section 4, we evolve the CO WD accretor along with the donor and
the binary orbit. Our results agree with Bildsten et al. (2006) but are more
consistent as our models track the time-dependent cooling of the donor. At
high mass transfer rates, $\dot{M}$, near period minimum, the accretor is
reheated due to accretion; as the orbit widens and $\dot{M}$ drops, the
accretor luminosity eventually becomes just that of a cooling WD. We compute
synthetic color-magnitude diagrams for our accretor models for comparisons
with observations in Section 5. We show that observed AM CVn systems appear
bluer and brighter than expected, which corresponds to a younger WD cooling
age. We conclude in Section 6.

In the following, we denote the accretor and the donor by the subscripts WD
and He. Quantities at the center of the stars are denoted by the subscript
$c$. Initial and final quantities are denoted by the superscripts $i$ and $f$.
Age since initiation of mass transfer will be used synonymously with age.

##  2 Computational Setup and Initial Stellar Model Construction

To obtain realistic, time-dependent mass transfer histories for AM CVn
binaries, we self-consistently evolve a He WD donor and a CO WD accretor using
the $\tt binary$ capability of the stellar evolution instrument Modules for
Experiments in Stellar Astrophysics (${\tt MESA}$ version 12778; Paxton et al.
2011, 2013, 2015, 2018, 2019). The equation of state (EOS) adopted for the WDs
in this study is a blend among OPAL (Rogers & Nayfonov, 2002) and SCVH (Saumon
et al., 1995) at low densities
($\log_{10}(\rho/\mathrm{g}\,\mathrm{cm}^{-3})\lesssim 2.9$), and PC at high
densities (Potekhin & Chabrier, 2010,
$\log_{10}(\rho/\mathrm{g}\,\mathrm{cm}^{-3})\gtrsim 2.9$). Radiative
opacities, where available, are taken from OPAL (Iglesias & Rogers, 1993,
1996) at $\log_{10}(T/\mathrm{K})\gtrsim 3.8$ and Ferguson et al. (2005) at
$\log_{10}(T/\mathrm{K})\lesssim 3.8$, and extrapolated otherwise (see Section
3.3). Electron conductive opacities are from Cassisi et al. (2007) with
corrections from Blouin et al. (2020) for He under moderate coupling and
moderate degeneracy. The latter corrections are implemented by modifying our
copy of the ${\tt MESA}$ source code, but will be included in future releases
of ${\tt MESA}$. Our nuclear net includes the NCO reaction chain (Bauer et
al., 2017) which may impact the helium shell thickness for unstable burning on
the accretor. The models presented in this work differ slightly from an
earlier version of models presented in van Roestel et al. (2021). The latter
adopted the HELM EOS (Timmes & Swesty, 2000) during the early phases of mass
transfer for the He WD donor, and included no electron conduction opacity
correction from Blouin et al. (2020) for both WDs. However, difference in
quantities like luminosity is insignificant (of order 10%). Our ${\tt MESA}$
input files are available at Zenodo111https://doi.org/10.5281/zenodo.5532940.

The He WD models are made using a modified version of the $\tt make\\_he\\_wd$
test suite, where a 1.5 $\mathrm{M_{\odot}}$ star with metallicity $Z=0.02$ is
evolved from the zero-age main sequence (ZAMS) to the formation of an inert
helium core of mass $0.15-0.18\,\mathrm{M_{\odot}}$. We then strip off the
envelope, and allow the bare He core to cool to the desired central entropy.
Extremely low mass (ELM) He WDs are expected to have a hydrogen envelope of
mass $\approx 10^{-3}\,\mathrm{M_{\odot}}$ (e.g., Istrate et al., 2016). When
accreted onto the CO WD the hydrogen may undergo unstable nuclear burning
(e.g., Kaplan et al., 2012), which is computationally expensive to follow and
has no thermal impact on the accretor, hence we neglect it in our study. We
leave the He WD with almost no hydrogen envelope ($M_{\mathrm{H}}\lesssim
10^{-20}\,\mathrm{M_{\odot}}$), which only impacts the initial mass transfer
phase, but has no impact on the mass transfer at late times (e.g., Kaplan et
al., 2012).

The He WD models have initial specific central entropies of
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=2.23,2.61,3.07,3.52$,
where $N_{\mathrm{A}}$ is Avogadro’s number. In Deloye et al. (2007), the
donor models are labeled by their initial central degeneracy parameter,
$\psi^{i}_{c,\mathrm{He}}=E_{\mathrm{F,c}}/k_{\mathrm{B}}T_{\rm
c}\approx\rho_{\rm c}/(1.2\times 10^{-8}T_{\rm
c}^{3/2})\,\mathrm{K}^{3/2}\mathrm{cm}^{3}\mathrm{g}^{-1}$, where $T_{\rm c}$,
$\rho_{\rm c}$ and $E_{\mathrm{F,c}}$ are the initial temperature, density and
electron Fermi energy at the center of the He WD. Our models correspond to
$\log_{10}(\psi^{i}_{c,\mathrm{He}})=2.89,2.39,1.82,1.34$ respectively. Our
models are chosen to be similar in range of $\psi^{i}_{c,\mathrm{He}}$ to
those in Deloye et al. (2007). Deloye et al. 2007 determined the distribution
of $\psi^{i}_{c,\mathrm{He}}$ at contact (see their Fig. 3) based on post-
common envelope conditions of the surviving double WD systems in the
population synthesis study of Nelemans et al. (2001), and modeled the mass
transfer of He WDs with $\log_{10}(\psi^{i}_{c,\mathrm{He}})\approx 1.1-3.5$.
Our coldest model is less degenerate than theirs due to inadequate EOS
coverage in ${\tt MESA}$ (see below).

The CO WD models are similarly made using a modified version of the $\tt
make\\_co\\_wd$ test suite. We evolve a ZAMS model of mass
$3-6\,\mathrm{M_{\odot}}$ with metallicity $Z=0.02$ up to helium shell burning
and the formation of a CO core of the desired mass. We then remove the
envelope and allow the CO core to cool to the desired starting core
temperature.

The He WD and CO WD models are then taken as initial models in ${\tt MESA}$
$\tt binary$ which self-consistently evolves the stellar structures of both
stars and the orbital parameters. We assume fully conservative mass transfer
at a rate $\dot{M}$ and model it using the Ritter scheme (Ritter, 1988). We
assume that orbital angular momentum loss is only due to the emission of
gravitational waves. While the accretor may be significantly spun up by the
accreted material at the expense of orbital angular momentum (Marsh et al.,
2004; Gokhale et al., 2007), we choose to model both stars as non-rotating as
a first approximation. Tides may play an important role in synchronizing both
components, which helps stabilize the orbit, and in impacting their thermal
evolution through tidal dissipation (e.g., Marsh et al., 2004; Fuller & Lai,
2012, 2014). Observations show that at long periods
($P_{\mathrm{orb}}=46\,\mathrm{min}$ for GP Com and $65\,\mathrm{min}$ for
V396 Hya) the accretor is rotating far slower than critical (Kupfer et al.,
2016).

The binary simulations begin with an initial orbital period of $\approx 20$
minutes with both components detached. The orbit decays due to gravitational
wave radiation and the He WD eventually fills its Roche lobe. We follow the
evolution of the system until the donor has cooled down sufficiently that
$\log Q=\log\rho-2\log T+12\geqslant 5$ in any of the cells. This choice was
made because for $\log Q\geqslant 5$, ${\tt MESA}$ switches to an ideal gas
equation of state which affects the donor mass-radius relation and
subsequently the orbital evolution. This typically occurs only when
$P_{\mathrm{orb}}\gtrsim 60\,\mathrm{min}$, allowing us to compare our
predictions with all known AM CVn systems. Future work with ${\tt MESA}$
incorporating the EOS’s of Chabrier et al. (2019) and Jermyn et al. (2021) can
probe the further evolution of the donor at even lower temperatures with
better numerical accuracy.

For the accretor, we set $\tt tau\\_factor=30$ for numerical convenience,
which places the outermost cell at an optical depth of $\tau=2/3\times 30=20$.
We also opted to adopt a grey Eddington $T-\tau$ relation instead of
interpolating from the DB tables provided by Odette Toloza and Detlev Koester.
The level of discrepancy in the accretor effective temperature for these two
approaches is $\lesssim 10\%$. We allow the latent heat of crystallization to
be released at a Coulomb coupling parameter $\Gamma=\langle Z^{5/3}_{i}\rangle
e^{2}/a_{e}k_{\mathrm{B}}T$ between 225 and 235 (see equation 1 of Bauer et
al., 2020), where $\langle Z^{5/3}_{i}\rangle$ is an average of $Z^{5/3}_{i}$
over all ion species with charge $Z_{i}$, and $a_{e}=(3/4\pi n_{e})^{1/3}$ is
the electron separation.

As we show in Section 3, the initial thermal properties and subsequent thermal
evolution of the donor are important in setting the binary orbital evolution
and $\dot{M}$ history. Therefore, for each
$M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ model of fixed
$S^{i}_{c,\mathrm{He}}$, we generate an adiabatic mass-radius relation, by
stripping the $0.15\,\mathrm{M_{\odot}}$ He WD model to the desired mass via
the $\tt relax\\_mass$ control in $\tt MESA$. The rate of mass change for this
construction is chosen such that the core of the model evolves adiabatically.
We obtain the cold WD (fully degenerate) limit similarly by stripping the
$0.15\,\mathrm{M_{\odot}}$ He WD models to the desired mass, and subsequently
allowing them to cool to a central temperature of
$T_{\mathrm{c,He}}=10^{5}\,\mathrm{K}$. To avoid the $\log Q\geqslant 5$
stopping condition, we irradiate the surface layer of the He WD model. This is
done in part for numerical convenience, but is also expected for realistic AM
CVn systems where, on geometrical grounds, the accreting WD should strongly
irradiate the donor.

For comparison, we also evolved a
$M^{i}_{\mathrm{He}}=0.35\,\mathrm{M_{\odot}}$ He star model taken from Brooks
et al. (2015). We chose $M^{i}_{\mathrm{WD}}=0.5\,\mathrm{M_{\odot}}$ and
$P^{i}_{\mathrm{orb}}=20\,\mathrm{min}$, which corresponds to the first model
shown in Table 1 of their paper and of Yungelson (2008). We also obtained an
adiabatic mass-radius relation for the He star using the $\tt relax\\_mass$
control described above. The initial model for this construction is taken from
the binary simulation when core He burning is quenched in the He star and its
mass decreases to $0.2\,\mathrm{M_{\odot}}$.

##  3 Donor Evolution

Figure 1:  Figure 2:  Donor mass-radius relation (top panel) and the
associated power-law index $n=\mathrm{d}\ln R_{\mathrm{He}}/\mathrm{d}\ln
M_{\mathrm{He}}$ (bottom panel), color-coded by their initial central specific
entropy, $S^{i}_{c,\mathrm{He}}$ (labeled by white lines in the colorbar). We
take $M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ (solid lines) or
$0.18\,\mathrm{M_{\odot}}$ (dashed lines), both with
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$. Square symbols indicate where
the donor thermal timescale, $\tau_{\mathrm{th}}=\int
c_{\mathrm{p}}T\mathrm{d}m/L_{\mathrm{He}}$, first becomes smaller than the
mass transfer timescale, $\tau_{\mathrm{m}}=M_{\mathrm{He}}/\dot{M}$. Circle
symbols indicate where the donor first becomes fully convective. For the
$M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ models, we show the adiabatic
(dot-dashed lines) and fully degenerate (taken as
$T_{\mathrm{c,He}}=10^{5}\,\mathrm{K}$; dotted grey line) mass-radius
relations. We also show a He star model (grey solid line) taken from Brooks et
al. (2015), and the corresponding adiabatic mass-radius relation starting from
$0.2\,\mathrm{M_{\odot}}$. For comparison, we show lines of constant
$P_{\mathrm{orb}}$ (light blue lines) using the mean density - period relation
for a Roche-lobe filling object.  Figure 3:  Figure 4:  Time evolution of the
mass transfer rate, $\dot{M}$ (top panel), orbital period, $P_{\mathrm{orb}}$
(middle panel), and donor mass, $M_{\mathrm{He}}$ (bottom panel), for various
initial donor central specific entropies, $S^{i}_{c,\mathrm{He}}$. The colored
solid (grey dot-dashed) lines are for
$M^{i}_{\mathrm{WD}}=0.75\,(1.05)\,\mathrm{M_{\odot}}$ and
$M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$, the colored dashed lines are
for $M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ and
$M^{i}_{\mathrm{He}}=0.18\,\mathrm{M_{\odot}}$, and the solid grey line is for
the He star model. For comparison, we show $\dot{M}\propto t^{-1.3}$ (dashed
grey line; Bildsten et al., 2006) as analytically expected.

The mass-transfer rate, $\dot{M}$, determines whether the accretor is heating
or cooling, and sets the evolution of binary parameters (e.g., orbital period,
$P_{\mathrm{orb}}$; Bildsten et al., 2006). It depends on the donor’s mass-
radius relation, which in turn is set by its thermal evolution (Deloye et al.,
2007). Thus we start by exploring the effect of varying the donor’s initial
central specific entropy, $S^{i}_{c,\mathrm{He}}$, on $\dot{M}$ and the age-
period relation. We run ${\tt MESA}$ models for initial donor masses of
$M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ or $0.18\,\mathrm{M_{\odot}}$
with $S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=2.23,2.61,3.07,3.52$.
As we are first exploring the $\dot{M}(t)$ and $P_{\mathrm{orb}}(t)$
evolution, we model the accretor as a point-mass with an initial
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$. We show the effect of an
initially more massive accretor by another model with
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.07$ and
$M^{i}_{\mathrm{WD}}=1.05\,\mathrm{M_{\odot}}$.

###  3.1 Donor Mass-Radius Relation

The resulting donor’s mass-radius relations are shown in the top panel of
Figure 2, and the associated power-law index $n=\mathrm{d}\ln
R_{\mathrm{He}}/\mathrm{d}\ln M_{\mathrm{He}}$ in the bottom panel. Initially,
the donor contracts ($n>0$) as the outermost radiative layer is stripped off
(see Deloye et al. 2007 and Kaplan et al. 2012 for in-depth discussions).
Eventually, the underlying layers drive an expansion in radius ($n<0$), with
the lowest entropy model closest to $n=-1/3$ as expected for a fully
degenerate object. In general, the less degenerate, higher entropy models show
a more positive $n$. As $M_{\mathrm{He}}$ decreases, the power-law slope
becomes more positive, in part because Coulomb effects become more important
(Deloye & Bildsten, 2003).

For less degenerate donors (especially the
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.52$ donor), $n$ becomes
more positive when they start to cool. Initially, the donors evolve
adiabatically as the mass change timescale,
$\tau_{\mathrm{m}}=M_{\mathrm{He}}/\dot{M}$, is much shorter than the thermal
timescale, $\tau_{\mathrm{th}}=\int
c_{\mathrm{p}}T\mathrm{d}m/L_{\mathrm{He}}$. This is confirmed by the
agreement of the $R(M)$ relations from the binary calculations (solid lines)
and the adiabatic models (dot-dashed lines) in Figure 2, below
$P_{\mathrm{orb}}\approx 30\,\mathrm{min}$ when
$\tau_{\mathrm{m}}\ll\tau_{\mathrm{th}}$. As $\dot{M}$ drops with time,
eventually $\tau_{\mathrm{th}}\lesssim\tau_{\mathrm{m}}$
($\tau_{\mathrm{th}}=\tau_{\mathrm{m}}$ is labeled by square symbols). Then
the donor becomes fully convective (labeled by circle symbols) and starts to
cool (Deloye et al., 2007), as we discuss in the end of this section.
Therefore, less degenerate donors start to converge to the mass-radius
relation of a fully degenerate He WD, as seen in Figure 2 (see also Fig. 8 of
Deloye et al., 2007).

The mass-radius relation depends only slightly on $M^{i}_{\mathrm{He}}$, as
seen in Figure 2. At a given $S^{i}_{c,\mathrm{He}}$, the difference between
the $M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ (solid lines) and
$0.18\,\mathrm{M_{\odot}}$ (dashed lines) tracks is only evident at larger
$M_{\mathrm{He}}$ ($P_{\mathrm{orb}}\lesssim 20\,\mathrm{min}$), and is
negligible at small $M_{\mathrm{He}}$ especially when donor cooling happens.

###  3.2 Mass Accretion Rate Histories

We show the time evolution of $\dot{M}$ in the top panel of Figure 4. As the
He WD donor fills its Roche lobe, the mass transfer rate rapidly rises to
$\sim 10^{-7}-10^{-8}$ $\mathrm{M_{\odot}}\,{\rm yr}^{-1}$. All $\dot{M}$’s
peak within $10^{6}$ yrs after mass transfer initiates, with more degenerate
donors having higher $\dot{M}$ due to their smaller radii (e.g., Kaplan et
al., 2012). Afterwards, $\dot{M}\propto t^{-1.3}$ as analytically expected
(Bildsten et al., 2006). At a fixed age, $\dot{M}$ is lower for more
degenerate donors. The ordering of $\dot{M}$ with donor degeneracy both at
peak $\dot{M}$ and in the power-law phase agrees with Deloye et al. 2007
(their Fig. 6). The top panel of Figure 4 also shows that $\dot{M}$ does not
vary significantly with $M^{i}_{\mathrm{WD}}$, as shown by the models with
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ and $1.05\,\mathrm{M_{\odot}}$
at $S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.07$.

The time evolution of $P_{\mathrm{orb}}$ is tied intimately to the donor’s
mass-radius relation and its thermal evolution, as shown in the middle panel
of Figure 4. Less degenerate donors have a larger $P_{\mathrm{orb}}$ at period
minimum, but $P_{\mathrm{orb}}$ evolution with time slows down after an age of
$\approx 10^{8}$ yrs, so that at the end, more degenerate donors catch up and
have larger $P_{\mathrm{orb}}$ at a fixed age. For the less degenerate donors,
the reason for their slower time evolution is that they eventually cool when
$\tau_{\mathrm{th}}\lesssim\tau_{\mathrm{m}}$.

Due to increased orbital angular momentum loss from gravitational wave
emission, a higher total mass also speeds up the $P_{\mathrm{orb}}$ evolution
with time. This is more clearly shown by the pair of
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.07$ models with
different $M^{i}_{\mathrm{WD}}$, and barely observable in the pairs of models
with different $M^{i}_{\mathrm{He}}$. However, a higher total mass increases
the likelihood of a helium flash being triggered on the accretor, since less
accumulated mass on the accretor is required given a higher $\dot{M}$ (a
higher $M^{i}_{\mathrm{He}}$ leads to a higher peak $\dot{M}$; see top panel
of Fig 4) and a higher $M^{i}_{\mathrm{WD}}$ (e.g., Bauer et al., 2017).
Indeed, our calculations that simultaneously evolve the stellar structure of
the accretor show that the $M_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ accretor
undergoes a helium flash for the
$M^{i}_{\mathrm{He}}=0.18\,\mathrm{M_{\odot}}$,
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.07$ model at
$P_{\mathrm{orb}}\approx 6\,\mathrm{min}$, while the
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.52$ model avoids a
helium flash due to the countering effect of a higher central entropy.

We show the evolution of $\dot{M}$ with $P_{\mathrm{orb}}$ in Figure 6. For
$P_{\mathrm{orb}}<40$ min, an initially hotter donor gives a higher $\dot{M}$
at fixed $P_{\mathrm{orb}}$. But starting at $P_{\mathrm{orb}}\approx 30$ min,
the $\dot{M}$ for the initially hotter donors start to converge to that for
cold donors, reflecting the eventual cooling of the donors to a fully
degenerate configuration. During this cooling, $\dot{M}$ drops more sharply
with $P_{\mathrm{orb}}$. As analytically derived by Cannizzo & Nelemans
(2015), $\dot{M}\propto P_{\mathrm{orb}}^{\xi}$, where $\xi=4(5-6n)/3(3n-1)$.
Initially as $n\approx-1/3$, $\xi\approx-4.67$, but during cooling, $n\approx
0$ and $\xi\approx-6.67$. Thus, the sharp drop in $\dot{M}$ with
$P_{\mathrm{orb}}$ is consistent with analytical expectations. Even though
hotter donors are at longer $P_{\mathrm{orb}}$ at an age of $10^{8}$ yrs (due
to a larger minimum $P_{\mathrm{orb}}$), this ordering reverses before an age
of $10^{9}$ yrs.

For comparison, we semi-analytically derive $\dot{M}-P_{\mathrm{orb}}$
relations for an adiabatic donor evolution (dot-dashed lines in Figure 6).
Given an adiabatic mass-radius relation from ${\tt MESA}$, we numerically
integrate the binary orbital parameters ($M_{\mathrm{He}},M_{\mathrm{WD}}$ and
binary separation $a$) assuming conservative mass transfer and orbital angular
momentum loss solely due to gravitational waves (see eqns 1 & 2 in Brooks et
al., 2015). The initial values are taken to be when $M_{\mathrm{He}}$
decreases to $0.14\,\mathrm{M_{\odot}}$ in the ${\tt MESA}$ binary run with
$M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ and
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$, which correspond roughly to
the moment of peak $\dot{M}$ in Figure 6. The excellent overlap between the
solid and dot-dashed lines in Figure 6 at $P_{\mathrm{orb}}\lesssim
30\,\mathrm{min}$ illustrates the initially adiabatic donor evolution. The
resulting $S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.52$ adiabatic
line corresponds approximately to the “hot” donor lines in Fig 2 of Bildsten
et al. (2006) who assumed adiabatic donor evolution as well. Since our models
eventually cool, they track closely “cool” models regardless of
$S^{i}_{c,\mathrm{He}}$ once $P_{\mathrm{orb}}\gtrsim 40-50\,\mathrm{min}$.

Figure 5:  Figure 6:  Mass transfer rates, $\dot{M}$, as a function of orbital
period, $P_{\mathrm{orb}}$. The models shown are the same as in Figure 4, with
the addition of semi-analytical lines for adiabatic donor evolution (dot-
dashed lines). The open circles denote fixed ages since initiation of mass
transfer. For comparison we show the theoretical critical mass transfer rates
from the disk instability model (Kotko et al., 2012),
$\dot{M}^{+}_{\mathrm{cr}}$ ($\dot{M}^{-}_{\mathrm{cr}}$) with an orange
(green) dashed line. The helium equivalent of dwarf novae are expected to
occur between these two limits (see Appendix A).

###  3.3 Donor Thermal Evolution Uncertainties Under Mass Transfer

The impact of the thermal evolution of the donor is substantial for the
observed $\dot{M}(P_{\mathrm{orb}})$ and $P_{\mathrm{orb}}(t)$ relations and
so deserves some additional scrutiny. This evolution is set by both the
initial entropy and whether entropy can be lost as mass is transferred, which
we now discuss.

The donor can only cool when its thermal timescale, $\tau_{\mathrm{th}}$, is
comparable to or shorter than the mass transfer timescale $\tau_{\mathrm{m}}$.
The former is inversely proportional to the donor luminosity,
$L_{\mathrm{He}}$, so lowering $L_{\mathrm{He}}$ may delay the cooling of the
donor, and speed up the evolution of $P_{\mathrm{orb}}$ with age.

Figure 7:  Figure 8:  Density-temperature profiles of the donor (left panel),
and entropy profiles of the donor (right panel), both color-coded by
$P_{\mathrm{orb}}$. Convection zones are labeled by thick blue lines. All
profiles range from $P_{\mathrm{orb}}=10$ to 60 min in 10 min intervals, and
we also include the initial model ($P_{\mathrm{orb}}\approx 20$ min).

The thermal evolution of the donor with
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ and
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.07$ is shown in Figure
8. The left panel shows a sequence of temperature-density profiles color-coded
by $P_{\mathrm{orb}}$, and the right panel shows the entropy profiles of the
same models. The black dashed line in the left panel shows the evolution of
the core in $T-\rho$ space. At $P_{\mathrm{orb}}\lesssim 30\,\mathrm{min}$,
the core evolves adiabatically as
$\tau_{\mathrm{m}}\lesssim\tau_{\mathrm{th}}$, and starts to cool for
$P_{\mathrm{orb}}\gtrsim 30\,\mathrm{min}$. This is corroborated by the
entropy profiles, which show that the central entropy is roughly constant for
$P_{\mathrm{orb}}\lesssim 30\,\mathrm{min}$ and only drops due to cooling at
longer periods. This decrease of central entropy happens even before the donor
becomes fully convective. Due to radiative diffusion, heat is always lost from
the center, i.e., $\epsilon_{\mathrm{grav},c}\equiv-
T_{c}(\mathrm{d}S_{c}/\mathrm{d}t)>0$. The change in central entropy becomes
apparent when $P_{\mathrm{orb}}\gtrsim 30\,\mathrm{min}$ as the timescale
increases from $\sim 10^{7}\,\mathrm{yr}$ to $\sim 10^{8}-10^{9}\,\mathrm{yr}$
and $\tau_{\mathrm{th}}$ eventually becomes comparable to, or shorter than,
$\tau_{\mathrm{m}}$. However, becoming fully convective accelerates the loss
of entropy from the center, the occurrence of which we now describe.

Before the donor becomes fully convective and cools, it has an outer
convective zone (blue thick lines in Figure 8), and a radiative core. The
outgoing luminosity in the thin radiative surface is set by
$L/M_{\mathrm{He}}=(64\pi G/3)(\sigma_{\mathrm{sb}}T^{4}/\kappa
P)\nabla_{\mathrm{ad}}\propto 1/\kappa$, at the top of the convective boundary
(Arras & Bildsten, 2006). Thus, the cooling luminosity of the donor is
inversely proportional to the opacity at the top radiative-convective boundary
for a fully convective star. At $P_{\mathrm{orb}}\gtrsim 20\,\mathrm{min}$,
this is at $\log_{10}(\rho/\mathrm{g}\,\mathrm{cm}^{-3})\approx 0$ and
$\log_{10}(T/\mathrm{K})\lesssim 4$. However, there are no radiative opacity
tables currently at this range, and instead ${\tt MESA}$ extrapolates from the
Ferguson et al. (2005) low-temperature opacity tables at the same
$\log_{10}T$. We defer investigation of the opacity to future studies, and
note that a higher opacity at the top radiative-convective boundary can delay
cooling of the donor and affect the age-period relation. We also note that
irradiation can play an important role in keeping the top radiative-convective
boundary at a higher $T$ (and potentially $\kappa$), but this is subject to
the same opacity uncertainties described above.

The same uncertainty in the donor cooling physics exists for the He star
donors at late times. From table 1 of Yungelson (2008), the typical age and
$P_{\mathrm{orb}}$ for $M^{i}_{\mathrm{He}}=0.35-0.4\,\mathrm{M_{\odot}}$ are
a few hundred Myrs and $\approx 40$ min at the end of their calculations. If
we take the corresponding binary parameters, and integrate forward taking
$n=-0.16$, we get an age since mass transfer between $\approx 1.5$ and $3$
Gyrs at $P_{\mathrm{orb}}=65$ min. These are lower limits nonetheless, since
$n=-0.16$ is fitted for $P_{\mathrm{orb}}\approx 10-35$ min. Eventually the He
star donors should start to cool (Deloye et al., 2007), as hinted by $n\gtrsim
0$ at low masses for the mass-radius relation in Fig. 4 of Yungelson (2008).
Indeed, this is explicitly shown by our He star model in the middle panel of
Figure 4, where the age increases significantly around
$P_{\mathrm{orb}}\approx 40\,\mathrm{min}$ as the He star starts to cool (see
also Figure 2).

To summarize, we have quantified the scatter in the age-$P_{\mathrm{orb}}$
relation considering all He WD models with donor cooling. In the middle panel
of Figure 4, the range of ages at which $P_{\mathrm{orb}}\approx 65$ min is
$4-5$ Gyrs. The importance of the period-age relation, itself tied closely to
the thermal properties of the donor, will become clear once we discuss the
heating and cooling of the accreting WD in Section 4.

##  4 Heating and Cooling of the Accreting WD

Figure 9:  Density-temperature profiles of the WD accretor during core-heating
(left panel, color-coded by $M_{\mathrm{WD}}$) and core-cooling (right panel,
color-coded by $P_{\mathrm{orb}}$). The core temperature reaches a maximum at
$P_{\mathrm{orb}}\approx 20$ min and $M_{\mathrm{WD}}\approx
0.87\,\mathrm{M_{\odot}}$, after which the WD cools and crystallizes (at a
Coulomb coupling parameter $\Gamma_{\mathrm{C/O}}=225$ for the core
composition of the ${\tt MESA}$ model, as shown by the grey dashed line in the
right panel). For each profile, we indicate the transition from the He
envelope to the CO core by an open circle, and convection zone by overplotting
a thick line. Core evolution is shown by the black dashed line and black
arrow. We indicate the transition to degeneracy by where the electron pressure
is twice the ideal gas electron pressure (dotted grey line), the start of
helium burning by where the helium burning rate is balanced by thermal
neutrino losses (dot-dashed black line), and two lines of constant Coulomb
coupling parameter $\Gamma_{\mathrm{He}}$ for pure He (dashed grey lines). The
two profiles with the labels * have $T_{\mathrm{eff}}$ falling within the DBV
instability strip ($\num{22400}\leqslant
T_{\mathrm{eff}}/\mathrm{K}\leqslant\num{32000}$; Córsico et al., 2019).

We now focus on the effects of mass transfer on the accretor properties,
particularly $T_{c,\mathrm{WD}}$. Given the small scatter in period-age
relation over a range of $T^{i}_{c,\mathrm{He}}$ and $M^{i}_{\mathrm{WD}}$, we
fix $M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$,
$S^{i}_{c,\mathrm{He}}/N_{\mathrm{A}}k_{\mathrm{B}}=3.07$, and study
$M^{i}_{\mathrm{WD}}=0.75$ and $0.85\,\mathrm{M_{\odot}}$.

Once settled onto the more massive CO WD, the accreted He releases heat while
it is compressed deep in the envelope. At the initially high $\dot{M}\sim
10^{-8}-10^{-7}\,\mathrm{M_{\odot}}\,{\rm yr}^{-1}$, the rate of compression
is higher than that of heat transport, resulting in a local temperature
increase (Nomoto, 1982; Townsley & Bildsten, 2004). This “compressional
heating” creates a temperature inversion in the envelope, clearly seen in
Panel (a) of Figure 9, as the sequence of temperature-density profiles of the
accretor at various times (color-coded by accretor mass). Heat transport from
the temperature peak occurs both towards the surface and into the colder core.
The latter can be seen as a thermal front propagating to the center in Figure
9, which reaches the center in a few $10^{7}$ years after the start of mass
transfer. Before then, the core evolves adiabatically. The
$(1\mathrm{-}3)\times 10^{7}$ years for the core to become heated is
consistent with the heat conduction time, $t_{\mathrm{cond}}$, as derived in
Townsley & Bildsten 2004 (their equation A2). All models presented in this
section avoided the occurrence of a helium flash, as the temperature of the
helium layer never reached a value adequate to trigger a thermonuclear
instability (see Figure 9).

When the age since mass transfer is comparable to
$t_{\mathrm{cond}}\approx(1\mathrm{-}3)\times 10^{7}$ yr, $T_{c,\mathrm{WD}}$
reaches a maximum and the WD cools. By then, $P_{\mathrm{orb}}\gtrsim 20$
minutes and the donor has reached $0.03\,\mathrm{M_{\odot}}$. The mass
accretion rate, $\dot{M}$, has dropped from peak by 1-2 orders-of-magnitude
(see Figure 4), diminishing the importance of compressional heating and
leading to a transition from core heating to core cooling (Bildsten et al.,
2006). Panel (a) of Figure 9 shows that near peak $T_{c,\mathrm{WD}}$, the
temperature inversion gradually disappears as the core is heated due to the
inward heat flow. We show the accretor profiles in the subsequent cooling
phase in Panel (b) of Figure 9, color-coded by orbital period. As the WD
further cools, the surface convection zone deepens and eventually
crystallization starts when $\Gamma_{c}=225$, at $P_{\mathrm{orb}}\approx 50$
min.

The time evolution of $T_{c,\mathrm{WD}}$ and $L_{\mathrm{WD}}$ is shown in
Figure 11. Consistent with Figure 9, the core first evolves adiabatically,
then becomes significantly heated once the conduction front reaches the center
at $10^{7}\,\mathrm{yrs}$, and cools thereafter. An initially cooler WD is
heated to a lower peak $T_{c,\mathrm{WD}}$ and hence shows a lower
$L_{\mathrm{WD}}$ around peak $T_{c,\mathrm{WD}}$. We note that, if a steadily
accreting disk is present, its bolometric accretion luminosity
$L_{\mathrm{acc}}\equiv GM_{\mathrm{WD}}\dot{M}/2R_{\mathrm{WD}}$ would always
dominate $L_{\mathrm{WD}}$. We discuss the challenges of differentiating
between luminosity from accretion and from the cooling WD in the Appendix.

Recent calculations by Blouin et al. (2020) show a lower conductive opacity in
WD envelopes, than the Cassisi et al. (2007) conductive opacities which are
used in typical WD cooling models (e.g., Bédard et al., 2020). We incorporated
this correction in all AM CVn models presented in this work. Compared to a
model without this correction, the reheated accretor reaches a slightly lower
peak central temperature, due to faster transport of heat out of the envelope.
In the cooling phase, the accretor has a higher luminosity at the same central
temperature (by $\approx 30\%$) and cools faster. Equivalently, in well
agreement with Blouin et al. (2020), this leads to a lower luminosity at a
fixed cooling age. However, the difference in inferred cooling ages is much
less than 1 Gyr, for a $0.9\,\mathrm{M_{\odot}}$ DB WD with
$T_{\mathrm{eff}}\approx 10^{4}\,\mathrm{K}$ (their Fig. 4), so its effect on
our calculations is small.

We compare our results with Bildsten et al. (2006), who modeled the thermal
state of the accretor using quasistatic envelope methods. For early times,
before the core gets heated, they found that the WD luminosity varies only
with $\dot{M}$ and $M_{\mathrm{WD}}$ as
$L_{\mathrm{WD}}\propto\dot{M}^{1.4}M_{\mathrm{WD}}^{-0.3}$ for
$\dot{M}>10^{-9}\mathrm{M_{\odot}}\,{\rm yr}^{-1}$. We compare this expression
with our ${\tt MESA}$ results in Figure 13. For $\dot{M}>3\times
10^{-9}\,\mathrm{M_{\odot}}\,{\rm yr}^{-1}$, the ${\tt MESA}$ models show no
dependence on $T^{i}_{c,\mathrm{WD}}$ as expected (Bildsten et al., 2006).
Compressional heating dominates $L_{\mathrm{WD}}$ at early times. The power-
law dependence of $L_{\mathrm{WD}}$ on $\dot{M}$ is slightly shallower than
predicted. For $\dot{M}\lesssim 3\times 10^{-9}\,\mathrm{M_{\odot}}\,{\rm
yr}^{-1}$, the WD cooling luminosity starts to dominate, and an initially
cooler WD is heated to a smaller peak $T_{c,\mathrm{WD}}$ (see Figure 11).
Thus between $\dot{M}=10^{-10}\,\mathrm{M_{\odot}}\,{\rm yr}^{-1}$ and
$3\times 10^{-9}\,\mathrm{M_{\odot}}\,{\rm yr}^{-1}$, a dependence of
$L_{\mathrm{WD}}$ on $T^{i}_{c,\mathrm{WD}}$ is seen, with a lower
$T^{i}_{c,\mathrm{WD}}$ giving a lower $L_{\mathrm{WD}}$ at fixed $\dot{M}$.

Figure 10:  Figure 11:  Time evolution of the accretor core temperature,
$T_{c,\mathrm{WD}}$ (top panel), and accretor luminosity (bottom panel), since
the start of mass transfer. The solid (dashed) lines show models with an
initial accretor mass $M^{i}_{\mathrm{WD}}=0.75\,(0.85)\,\mathrm{M_{\odot}}$,
each set with a range of initial core temperatures for the accretor,
$T^{i}_{c,\mathrm{WD}}$ ($5\times 10^{6}$ for deep-blue to $2\times 10^{7}$ K
for light green). We show where the models start to crystallize
($\Gamma_{c}=225$) by a circle (square) symbol for the
$M^{i}_{\mathrm{WD}}=0.75\,(0.85)\,\mathrm{M_{\odot}}$ models. We also show
the accretion luminosity, $L_{\mathrm{acc}}\equiv
GM_{\mathrm{WD}}\dot{M}/2R_{\mathrm{WD}}$ (dot-dashed line).  Figure 12:
Figure 13:  Luminosity of the accretor, $L_{\mathrm{WD}}$, as a function of
the mass transfer rate, $\dot{M}$, during the rapid accretion phase for the
same set of models as in Figure 11. We show
$L_{\mathrm{WD}}\propto\dot{M}^{1.4}$ (dot-dashed grey line) as expected from
Bildsten et al. (2006) for comparison.  Figure 14:  Figure 15:  Luminosity of
the accretor, $L_{\mathrm{WD}}$, as a function of the core temperature of the
accretor, $T_{c,\mathrm{WD}}$, during the core-cooling phase. The models are
the same as in Figure 11. At late times, $L_{\mathrm{WD}}\propto T_{\rm
c}^{2.5}$ as expected for cooling WDs, and is independent of
$T^{i}_{c,\mathrm{WD}}$ (Bildsten et al., 2006).

Bildsten et al. (2006) also predicted that at late times
($P_{\mathrm{orb}}\gtrsim 30\,\mathrm{min}$), the luminosity of the accretor
should follow that of a cooling WD,
$L_{\mathrm{WD}}\approx\mathrm{L_{\odot}}(T_{\rm
c}/10^{8}\,\mathrm{K})^{2.5}(M_{\mathrm{WD}}/0.6\,\mathrm{M_{\odot}})$. Figure
15 shows excellent agreement between our models and the predicted relation.
The $L_{\mathrm{WD}}$ tracks for the
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ models are roughly 0.89 times
those for the $M^{i}_{\mathrm{WD}}=0.85\,\mathrm{M_{\odot}}$ models,
consistent with the linear $M_{\mathrm{WD}}$ dependence of $L_{\mathrm{WD}}$
and reflects the ratio $0.9/1.0=0.9$ in the final WD mass ($M_{\mathrm{WD}}$
is nearly constant beyond $P_{\mathrm{orb}}\gtrsim 20\,\mathrm{min}$, see
bottom panel of Fig 4). More importantly, at late times
($P_{\mathrm{orb}}\gtrsim 30\,\mathrm{min}$), the accretor evolves as a
cooling WD independent of $\dot{M}$ and $T^{i}_{c,\mathrm{WD}}$. Then, the
binary age since mass transfer can be taken as a cooling age for the accretor,
since the few $10^{7}$ yrs required to heat the WD is insignificant compared
to the expected 4-5 Gyrs required for the binary to evolve to, say,
$P_{\mathrm{orb}}\approx 65$ min (see Section 3).

However, inferring the cooling ages of AM CVn accretors (see Section 5) is
subject to the caveat that, compared to normal CO WDs, AM CVn accretors may
have a significant fraction of their total mass in a thick He shell, if
unstable He burning can be avoided (e.g., our
$M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$,
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ model would have
$q_{\mathrm{He}}=M_{\mathrm{He}}/M_{\mathrm{WD}}=0.15/0.9\approx 0.17$ in the
end). Since He has a different specific heat capacity than a C/O mixture, the
cooling luminosity of an AM CVn accretor would be different than that of a CO
WD of the same total mass, at a given cooling age.

If we combine $L_{\mathrm{WD}}=-C_{\mathrm{V}}(\mathrm{d}T_{\rm
c}/\mathrm{d}t)$, where $C_{\mathrm{V}}$ is the total heat capacity, with
$L_{\mathrm{WD}}\propto T_{\rm c}^{5/2}$ (e.g., Townsley & Bildsten, 2004) and
solve for the cooling time $\tau_{\mathrm{cool}}$, then
$\tau_{\mathrm{cool}}\propto C_{\mathrm{V}}L_{\mathrm{WD}}^{-3/5}$. As seen in
Figure 9, at the ages $1-5$ Gyrs that we are interested in, the C/O is near
the crystallization limit and the He is approximately in the liquid state.
Then we can approximate, for an AM CVn accretor that has a C/O core mass
$M_{\mathrm{CO}}$ and He shell mass $M_{\mathrm{He}}$, $C_{\mathrm{V}}\approx
3(k_{\mathrm{B}}M_{\mathrm{CO}}/\mu_{\mathrm{CO}}m_{\mathrm{p}})+2(k_{\mathrm{B}}M_{\mathrm{He}}/\mu_{\mathrm{He}}m_{\mathrm{p}})$,
where $\mu_{\mathrm{He}}=4$, $\mu_{\mathrm{CO}}\approx 14$, and the specific
heat capacities for C/O and He are motivated by our ${\tt MESA}$ models (see
also, Baiko & Yakovlev, 2019). Similarly, for a CO WD of the same total mass
$M_{\mathrm{tot}}=M_{\mathrm{CO}}+M_{\mathrm{He}}$, $C_{\mathrm{V}}\approx
3(k_{\mathrm{B}}M_{\mathrm{tot}}/\mu_{\mathrm{CO}}m_{\mathrm{p}})$. Hence, for
a CO WD and an AM CVn accretor of the same total mass $M_{\mathrm{tot}}$ and
cooling luminosity $L_{\mathrm{WD}}$, the ratio of their cooling ages is
$\tau_{\mathrm{cool,AM\,CVn}}/\tau_{\mathrm{cool,C/O}}\approx(M_{\mathrm{CO}}/M_{\mathrm{tot}})+(2/3)(\mu_{\mathrm{CO}}/\mu_{\mathrm{He}})(M_{\mathrm{He}}/M_{\mathrm{tot}})$.
For our fiducial case, $M^{i}_{\mathrm{He}}=0.15\,\mathrm{M_{\odot}}$ and
$M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$, this ratio is about $1.22$.

Therefore, when comparing the contours of constant age for CO WDs with thin He
envelopes ($q_{\mathrm{He}}=10^{-2}$) in Figure 19 to observed data points,
the AM CVn accretors may have older actual cooling ages than inferred, by tens
of percent. Nonetheless, as we will show in Section 5, there exists a cooling
age discrepancy of $\approx 3-4$ Gyrs at $P_{\mathrm{orb}}\approx
65\,\mathrm{min}$ between theory and observations, so the corrections to the
cooling age presented here are insignificant.

Finally, the accretor effective temperature, $T_{\mathrm{eff,WD}}$, is shown
as a function of $P_{\mathrm{orb}}$ in Figure 17. Models with
$M^{i}_{\mathrm{WD}}=0.65,0.75$ and $0.85\,\mathrm{M_{\odot}}$, all with
$T^{i}_{c,\mathrm{WD}}=2\times 10^{7}$ K and identical initial donor models,
are shown. The WD starts cooling at $\approx 20$ min, so subsequent evolution
of $T_{\mathrm{eff,WD}}$ with $P_{\mathrm{orb}}$ is set by the period-age
relation of the binary. For comparison, we show the $M^{i}_{\mathrm{WD}}=0.65$
models with a “hot”/“cold” donor from Bildsten et al. (2006). Our
$M^{i}_{\mathrm{WD}}=0.65$ accretor model initially tracks closely the “hot”
donor model and later at long periods ($P_{\mathrm{orb}}\gtrsim
60\,\mathrm{min}$) the “cold” donor model. This is because our donor model
starts with high-entropy and we allow for its subsequent cooling. The period-
age relation then starts to deviate from the “hot” donor model which cools
adiabatically, and to converge towards the “cool” donor model. We also show
the DBV strip ($\num{22400}\leqslant
T_{\mathrm{eff}}/\mathrm{K}\leqslant\num{32000}$; e.g., Córsico et al., 2019),
though no AM CVn accretors are presently known to be pulsating. Finding one
would certainly enable a new probe of both the WD thermal state and thickness
of the He layer.

Figure 16:  Figure 17:  Effective temperatures of the accretor,
$T_{\mathrm{eff,WD}}$, as a function of orbital period, $P_{\mathrm{orb}}$.
Here we show our models with $T^{i}_{c,\mathrm{WD}}=2\times 10^{7}$ K and
$M^{i}_{\mathrm{WD}}=0.65,0.75$ or $0.85\,\mathrm{M_{\odot}}$, the
$M^{i}_{\mathrm{WD}}=0.65\,\mathrm{M_{\odot}}$ with “hot” or “cold” donor
models (dotted orange or light-blue lines) from Bildsten et al. (2006), and
the DBV instability strip (blue stripe; e.g., Córsico et al., 2019) for
comparison.

##  5 Comparison to Observations

We showed that at $\dot{M}\lesssim 3\times 10^{-9}\,\mathrm{M_{\odot}}\,{\rm
yr}^{-1}$, the accretor luminosity is just that of a cooling WD. If in
addition we assume that, active accretion is not contaminating the data
(Bildsten et al., 2006), then combining a cooling WD model with the
age-$P_{\mathrm{orb}}$ relation yields a theoretical prediction for the
luminosity of an AM CVn system at a given $P_{\mathrm{orb}}$. We compare this
prediction for the accretor to observed systems, and show that the observed
systems appear $\approx 3-4\,\mathrm{Gyrs}$ younger than predicted at
$P_{\mathrm{orb}}\approx 65\,\mathrm{min}$.

To compare with observations, we compute absolute magnitudes from our ${\tt
MESA}$ models. Given an effective temperature, $T_{\mathrm{eff}}$, and surface
gravity, $\log g$, from ${\tt MESA}$, we interpolate from the pure-helium
model atmosphere table from the Montreal group (Bergeron et al., 2011; Blouin
et al., 2018; Bédard et al., 2020)222https://www.astro.umontreal.ca/
bergeron/CoolingModels/ (catalog
https://www.astro.umontreal.ca/$∼$bergeron/CoolingModels/) and obtain a
bolometric correction for various photometric systems.

Our V-band absolute magnitude, $M_{\mathrm{V}}$, when plotted against
$P_{\mathrm{orb}}$, follows closely the “cool” donor track of Bildsten et al.
(2006). However, we discovered that their $M_{\mathrm{V}}$ is always lower by
$\approx 0.7\,\mathrm{mag}$ (i.e., brighter) than the $M_{\mathrm{V}}$ we
obtain using the Montreal bolometric corrections with the corresponding
$T_{\mathrm{eff}}$ and an assumed $\log g=8.5$. For example, at
$P_{\mathrm{orb}}\approx 63\,\mathrm{min}$, their “cool” donor model has
$T_{\mathrm{eff,WD}}\approx 6000\,\mathrm{K}$ which should give
$M_{\mathrm{V}}\approx 14.949$, but their $M_{\mathrm{V}}$ reads $\approx
14.3$ instead. This explains why our $M_{\mathrm{V}}-P_{\mathrm{orb}}$ track
agrees better with their “cool” donor track than with their “hot” donor track,
despite better agreement in $T_{\mathrm{eff}}-P_{\mathrm{orb}}$ with the
latter.

Regardless, as shown by Ramsay et al. (2018) (their Fig. 2), our $V$ or
$g$-band absolute magnitudes, $M_{\mathrm{V}}$ or $M_{\mathrm{g}}$, are lower
than observed by 2 magnitudes at $P_{\mathrm{orb}}\approx 60$ min. We show
this discrepancy instead by comparing our ${\tt MESA}$ tracks on a color-
magnitude diagram ($M_{\mathrm{g}}$ vs $g-r$ for SDSS and Pan-STARRS-1
filters; Fig 19), to observed AM CVn systems with $P_{\mathrm{orb}}>50$ min,
both color-coded by $P_{\mathrm{orb}}$. By $P_{\mathrm{orb}}\approx 20$ min,
our $M^{i}_{\mathrm{WD}}=0.75\,\mathrm{M_{\odot}}$ (or
$0.85\,\mathrm{M_{\odot}}$) accretor has accumulated most of the donor’s mass,
and starts to cool as a $M^{f}_{\mathrm{WD}}=0.9\,\mathrm{M_{\odot}}$ (or
$1.0\,\mathrm{M_{\odot}}$) WD. This explains the good agreement of our ${\tt
MESA}$ tracks even before $P_{\mathrm{orb}}\approx 40$ min with the
corresponding $0.9\,\mathrm{M_{\odot}}$ (or $1.0\,\mathrm{M_{\odot}}$) cooling
tracks from Bédard et al. (2020).

Nevertheless, even though the observed long-period systems lie close to the WD
cooling tracks, they appear brighter and bluer than expected from our ${\tt
MESA}$ tracks for the corresponding $P_{\mathrm{orb}}$. If we interpret the
observed source as dominated by the luminosity of the cooling WD, then their
photometric cooling ages, as inferred from the CMD, are between 0.5 and 2 Gyr.
In contrast, the ${\tt MESA}$ age-period relation (Sec 3) would give a cooling
age of $4-5$ Gyr at $P_{\mathrm{orb}}\approx 65$ min.

Figure 18:  Figure 19:  Color-magnitude diagram (Mg vs $g-r$) of the ${\tt
MESA}$ models for SDSS (Pan-STARRS 1) bands in the top (bottom) panel. The red
(orange) dashed line shows the ${\tt MESA}$ models with
$T^{i}_{c,\mathrm{WD}}=0.75\,(0.85)\,\mathrm{M_{\odot}}$, where we apply
bolometric corrections from the DB atmosphere models from the Montreal group
(Bergeron et al., 2011; Blouin et al., 2018; Bédard et al., 2020). We label
where the ${\tt MESA}$ models reach a certain $P_{\mathrm{orb}}$ with open
circle symbols. For comparison we show their DB cooling models at
$M_{\mathrm{WD}}=0.6-1.0\,\mathrm{M_{\odot}}$ (blue solid lines) and contours
of constant cooling age (grey dot-dashed lines), and a contour of when ${\tt
MESA}$ CO WD cooling models start crystallizing ($\Gamma_{c}>225$; dotted blue
line). The error bars show observed AM CVn systems (Ramsay et al., 2018) with
$P_{\mathrm{orb}}>50$ min, color-coded by $P_{\mathrm{orb}}$.

##  6 Conclusion

In this work, we perform binary calculations with ${\tt MESA}$ for AM CVn
binaries with a He WD donor. In Section 3, we evolved He WD donor models of
various initial central specific entropies, $S^{i}_{c,\mathrm{He}}$, and
showed that the initial entropy and the subsequent thermal evolution of the
donor dictates the mass transfer history and hence the age-period evolution of
the binary. We find that at $P_{\mathrm{orb}}\approx 40\,\mathrm{min}$, the
donor starts to cool as its thermal timescale, $\tau_{\mathrm{th}}$, is
comparable to or shorter than its mass transfer timescale,
$\tau_{\mathrm{m}}$, although this is subject to uncertainties in the surface
opacity as described in Section 3.3. We then evolve the accretor along with
the donor and the orbit in Section 4. At the initially high $\dot{M}$ near
period minimum, the accretor is heated due to compressional heating. Later as
$\dot{M}$ drops while the orbit widens, the accretor behaves simply as a
cooling WD (at $P_{\mathrm{orb}}\gtrsim 30\,\mathrm{min}$). Our calculations
show well agreement with the semi-analytic predictions of Bildsten et al.
(2006), but are more accurate since we self-consistently consistent both
binary components.

Given that we theoretically expect that accretor is simply a cooling WD at
$P_{\mathrm{orb}}\gtrsim 30\,\mathrm{min}$, we compute synthetic color-
magnitude diagrams for our accretor models and compare with observations
(Ramsay et al., 2018) in Section 5. We show that the observed systems, if
interpreted as the cooling WD, appear much younger (in terms of cooling age)
than expected at the corresponding $P_{\mathrm{orb}}$ (by $\approx 3-4$ Gyrs
at $P_{\mathrm{orb}}\approx 65$ min). This may be attributed in part to
uncertainties in the initial entropy distribution and subsequent cooling of
the donor, which affects the age-period relation of our models. Better
understanding of the donor cooling can be obtained by incorporating opacities
of warm dense helium in ${\tt MESA}$. Another possibility is that the observed
luminosities are contaminated by the accretion disk or boundary layer
(discussed in the Appendix), although this appears unlikely in the case of
Gaia14aae as revealed by eclipse modeling (Green et al., 2018). Finally,
isolated WDs show evidence of a cooling delay (e.g., Cheng et al., 2019; Kilic
et al., 2020), and we speculate that a similar phenomenon in the accretor may
well explain the observed discrepancy. The recent discovery of 5 new eclipsing
AM CVn systems with $P_{\mathrm{orb}}=35-62\,\mathrm{min}$ by van Roestel et
al. (2021), which constrains the luminosity contributions from different
components and the degree of cooling of the donor, will help disentangle these
various scenarios.

We thank the referee for their constructive suggestions that have greatly
improved our manuscript. We thank Evan Bauer and Josiah Schwab for helpful
conversations about WD cooling and guidance on using ${\tt MESA}$, and the
latter in addition for implementing the Blouin et al. (2020) electron
conductive opacity corrections in ${\tt MESA}$. We are grateful to Matthew
Green, Thomas Kupfer, Tom Marsh, and Jan van Roestel, for valuable discussions
about observations of AM CVn systems. We thank Sihao Cheng for making the
python package WD_models publicly available. This research was supported in
part by the National Science Foundation through grant PHY-1748958, and the
Gordon and Betty Moore Foundation through grant GBMF5076. Use was made of
computational facilities purchased with funds from the National Science
Foundation (CNS-1725797) and administered by the Center for Scientific
Computing (CSC). The CSC is supported by the California NanoSystems Institute
and the Materials Research Science and Engineering Center (MRSEC; NSF DMR
1720256) at UC Santa Barbara.

##  Appendix A Accretion Luminosity

We now address the possibility that observed AM CVn systems have contributions
from the accretion disk or boundary layer. As the bottom panel of Figure 11
shows, if accretion is steady, the accretion luminosity,
$L_{\mathrm{acc}}\equiv GM_{\mathrm{WD}}\dot{M}/2R_{\mathrm{WD}}$, would be
more than 5 times greater than the WD luminosity, $L_{\mathrm{WD}}$, at all
times.

If we assume steady-state accretion, and model the disk as a multicolor
blackbody, we would expect the disk to dominate over the cooling WD in the
V-band. We note that Bildsten et al. (2006) proposed the opposite, because
they assumed a uniform, averaged disk temperature defined by
$L_{\mathrm{acc}}=S\sigma T^{4}$, where $S$ is the disk surface area. Their
approach underestimated the temperature in the inner disk and hence the disk
contribution to the $V$-band (again, assuming steady state). This implies we
cannot expect a priori that the cooling WD dominates over the disk
contribution in the V-band.

However, the disk may not be in a steady state between 20 and 60 minutes
(Ramsay et al., 2018; Rivera Sandoval et al., 2020, 2021), as it is expected
to be thermally and viscously unstable if the mass transfer rate, $\dot{M}$,
is between the limits $\dot{M}^{+}_{\mathrm{cr}}$ and
$\dot{M}^{-}_{\mathrm{cr}}$ (Kotko et al., 2012, their eqn. A2 for $Z=0.02$),
shown in Figure 6. The upper limit, $\dot{M}^{+}_{\mathrm{cr}}$, is evaluated
at the tidal radius $R_{\mathrm{t}}=0.6a/(1+q)$, where $a$ is the binary
separation and $q=M_{\mathrm{He}}/M_{\mathrm{WD}}$ is the mass ratio. The
lower limit, $\dot{M}^{-}_{\mathrm{cr}}$, is evaluated at $R_{\mathrm{WD}}$,
the accretor radius at the end of the simulation (essentially the radius of a
$0.75+0.15=0.9\,\mathrm{M_{\odot}}$ WD).

As noted by Cannizzo & Nelemans (2015), Kotko et al. (2012) did not account
for the zero-torque inner boundary condition in converting from
$T_{\mathrm{eff}}(R)$ to $\dot{M}(R)$ locally in the disk. This affects the
lower stability condition $\dot{M}^{-}_{\mathrm{cr}}$, since the global
condition for the disk to be cold and stable is
$T_{\mathrm{eff}}(R)<T^{-}_{\mathrm{eff,cr}}(R)$ for every radius $R$ within
the disc. Because the disk effective temperature peaks at
$(49/36)R_{\mathrm{WD}}$, with peak value given by $0.488T_{*}$ where
$T_{*}=3GM_{\mathrm{WD}}\dot{M}/(8\pi
R_{\mathrm{WD}}^{3}\sigma_{\mathrm{sb}})$, the lower stability condition
should actually be
$0.488T_{*}<T^{-}_{\mathrm{eff,cr}}(R=(49/36)R_{\mathrm{WD}})$333The Kotko et
al. 2012 approach essentially takes
$T_{*}<T^{-}_{\mathrm{eff,cr}}(R=R_{\mathrm{WD}})$ and solves for $\dot{M}$..
Assuming $T^{-}_{\mathrm{eff,cr}}\propto R^{\beta}$, where
$\beta\approx-0.09$, and solving for $\dot{M}$, the actual
$\dot{M}^{-}_{\mathrm{cr}}$ should be a factor
$(49/36)^{4\beta}/0.488^{4}\approx 15.8$ times
$\dot{M}^{-}_{\mathrm{cr}}(R=R_{\mathrm{WD}})$. Thus, the green dashed line in
the top panel of Figure 6 should be 15.8 times higher if we account for the
zero-torque inner boundary condition, and we would expect no disk outbursts
for $P_{\mathrm{orb}}\gtrsim 50$ min. However, this is to be taken with a
grain of salt since AM CVn disk outbursts are yet to be fully understood
(e.g., Kotko et al., 2012; Rivera Sandoval et al., 2020).

Moreover, observations suggest an insignificant disk contribution for long-
period AM CVn systems. Eclipse modeling of Gaia14aae suggests that the WD
contribution is $\approx 80\%$ (Green et al., 2018, though a He WD donor is
unlikely for this system). In addition, optical spectra and photometry of
long-period systems can be fitted with a single blackbody of temperature
$T_{\mathrm{eff}}\approx 10^{4}\,\mathrm{K}$ (e.g., SDSS J1137+4054,
$P_{\mathrm{orb}}=59.6$ min, and SDSS J1505+0659, $P_{\mathrm{orb}}=67.8$ min;
Carter et al., 2014), whereas an accretion disk spectrum is expected to be
flatter (e.g., Nagel et al., 2009).

These observational constraints can be met by an optically thick boundary
layer. A boundary layer can be consistent with eclipse modeling by virtue of
its compact emitting area. Coincidentally, the “accretion temperature”
$T_{\mathrm{acc}}$, defined by $L_{\mathrm{acc}}=4\pi
R_{\mathrm{WD}}^{2}T_{\mathrm{acc}}^{4}$, gives the inferred blackbody
temperature $10^{4}\,\mathrm{K}$. Note that, however, this is difficult from a
theoretical point-of-view, because the boundary layer may be optically thin
for $\dot{M}\lesssim 10^{-10}\,\mathrm{M_{\odot}}\,{\rm yr}^{-1}$, and it may
not spread over the entire WD surface which implies an effective temperature
much higher than $T_{\mathrm{acc}}$ (e.g., Piro & Bildsten, 2004).

## References

  * Amaro-Seoane et al. (2013) Amaro-Seoane, P., Aoudia, S., Babak, S., et al. 2013, GW Notes, 6, 4.  https://arxiv.org/abs/1201.3621
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786.  https://arxiv.org/abs/1702.00786
  * Arras & Bildsten (2006) Arras, P., & Bildsten, L. 2006, ApJ, 650, 394, doi: 10.1086/506011
  * Ashok & Banerjee (2003) Ashok, N. M., & Banerjee, D. P. K. 2003, A&A, 409, 1007, doi: 10.1051/0004-6361:20031160
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Baiko & Yakovlev (2019) Baiko, D. A., & Yakovlev, D. G. 2019, MNRAS, 490, 5839, doi: 10.1093/mnras/stz3029
  * Bauer et al. (2017) Bauer, E. B., Schwab, J., & Bildsten, L. 2017, ApJ, 845, 97, doi: 10.3847/1538-4357/aa7ffa
  * Bauer et al. (2020) Bauer, E. B., Schwab, J., Bildsten, L., & Cheng, S. 2020, ApJ, 902, 93, doi: 10.3847/1538-4357/abb5a5
  * Bédard et al. (2020) Bédard, A., Bergeron, P., Brassard, P., & Fontaine, G. 2020, ApJ, 901, 93, doi: 10.3847/1538-4357/abafbe
  * Bergeron et al. (2011) Bergeron, P., Wesemael, F., Dufour, P., et al. 2011, ApJ, 737, 28, doi: 10.1088/0004-637X/737/1/28
  * Bildsten et al. (2007) Bildsten, L., Shen, K. J., Weinberg, N. N., & Nelemans, G. 2007, ApJ, 662, L95, doi: 10.1086/519489
  * Bildsten et al. (2006) Bildsten, L., Townsley, D. M., Deloye, C. J., & Nelemans, G. 2006, ApJ, 640, 466, doi: 10.1086/500080
  * Blouin et al. (2018) Blouin, S., Dufour, P., & Allard, N. F. 2018, ApJ, 863, 184, doi: 10.3847/1538-4357/aad4a9
  * Blouin et al. (2020) Blouin, S., Shaffer, N. R., Saumon, D., & Starrett, C. E. 2020, ApJ, 899, 46, doi: 10.3847/1538-4357/ab9e75
  * Breivik et al. (2018) Breivik, K., Kremer, K., Bueno, M., et al. 2018, ApJ, 854, L1, doi: 10.3847/2041-8213/aaaa23
  * Brooks et al. (2015) Brooks, J., Bildsten, L., Marchant, P., & Paxton, B. 2015, ApJ, 807, 74, doi: 10.1088/0004-637X/807/1/74
  * Cannizzo & Nelemans (2015) Cannizzo, J. K., & Nelemans, G. 2015, ApJ, 803, 19, doi: 10.1088/0004-637X/803/1/19
  * Carter et al. (2014) Carter, P. J., Gänsicke, B. T., Steeghs, D., et al. 2014, MNRAS, 439, 2848, doi: 10.1093/mnras/stu142
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Chabrier et al. (2019) Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, 51, doi: 10.3847/1538-4357/aaf99f
  * Cheng et al. (2019) Cheng, S., Cummings, J. D., & Ménard, B. 2019, ApJ, 886, 100, doi: 10.3847/1538-4357/ab4989
  * Córsico et al. (2019) Córsico, A. H., Althaus, L. G., Miller Bertolami, M. M., & Kepler, S. O. 2019, A&A Rev., 27, 7, doi: 10.1007/s00159-019-0118-4
  * Deloye & Bildsten (2003) Deloye, C. J., & Bildsten, L. 2003, ApJ, 598, 1217, doi: 10.1086/379063
  * Deloye et al. (2007) Deloye, C. J., Taam, R. E., Winisdoerffer, C., & Chabrier, G. 2007, MNRAS, 381, 525, doi: 10.1111/j.1365-2966.2007.12262.x
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fuller & Lai (2012) Fuller, J., & Lai, D. 2012, MNRAS, 421, 426, doi: 10.1111/j.1365-2966.2011.20320.x
  * Fuller & Lai (2014) —. 2014, MNRAS, 444, 3488, doi: 10.1093/mnras/stu1698
  * Gokhale et al. (2007) Gokhale, V., Peng, X. M., & Frank, J. 2007, ApJ, 655, 1010, doi: 10.1086/510119
  * Green et al. (2018) Green, M. J., Marsh, T. R., Steeghs, D. T. H., et al. 2018, MNRAS, 476, 1663, doi: 10.1093/mnras/sty299
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Iben & Tutukov (1991) Iben, Icko, J., & Tutukov, A. V. 1991, ApJ, 370, 615, doi: 10.1086/169848
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752, doi: 10.1086/172958
  * Iglesias & Rogers (1996) —. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Istrate et al. (2016) Istrate, A. G., Marchant, P., Tauris, T. M., et al. 2016, A&A, 595, A35, doi: 10.1051/0004-6361/201628874
  * Jermyn et al. (2021) Jermyn, A. S., Schwab, J., Bauer, E., Timmes, F. X., & Potekhin, A. Y. 2021, ApJ, 913, 72, doi: 10.3847/1538-4357/abf48e
  * Kaplan et al. (2012) Kaplan, D. L., Bildsten, L., & Steinfadt, J. D. R. 2012, ApJ, 758, 64, doi: 10.1088/0004-637X/758/1/64
  * Kilic et al. (2020) Kilic, M., Bergeron, P., Kosakowski, A., et al. 2020, ApJ, 898, 84, doi: 10.3847/1538-4357/ab9b8d
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas: Proceedings of the 20th International Conference on Electronic Publishing, IOS Press, 87
  * Kotko et al. (2012) Kotko, I., Lasota, J. P., Dubus, G., & Hameury, J. M. 2012, A&A, 544, A13, doi: 10.1051/0004-6361/201219156
  * Kremer et al. (2017) Kremer, K., Breivik, K., Larson, S. L., & Kalogera, V. 2017, ApJ, 846, 95, doi: 10.3847/1538-4357/aa8557
  * Kupfer et al. (2016) Kupfer, T., Steeghs, D., Groot, P. J., et al. 2016, MNRAS, 457, 1828, doi: 10.1093/mnras/stw126
  * Marsh et al. (2004) Marsh, T. R., Nelemans, G., & Steeghs, D. 2004, MNRAS, 350, 113, doi: 10.1111/j.1365-2966.2004.07564.x
  * Nagel et al. (2009) Nagel, T., Rauch, T., & Werner, K. 2009, A&A, 499, 773, doi: 10.1051/0004-6361/200811562
  * Nather et al. (1981) Nather, R. E., Robinson, E. L., & Stover, R. J. 1981, ApJ, 244, 269, doi: 10.1086/158704
  * Nelemans et al. (2001) Nelemans, G., Portegies Zwart, S. F., Verbunt, F., & Yungelson, L. R. 2001, A&A, 368, 939, doi: 10.1051/0004-6361:20010049
  * Nelemans et al. (2004) Nelemans, G., Yungelson, L. R., & Portegies Zwart, S. F. 2004, MNRAS, 349, 181, doi: 10.1111/j.1365-2966.2004.07479.x
  * Nomoto (1982) Nomoto, K. 1982, ApJ, 253, 798, doi: 10.1086/159682
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9, 21
  * Piro (2019) Piro, A. L. 2019, ApJ, 885, L2, doi: 10.3847/2041-8213/ab44c4
  * Piro & Bildsten (2004) Piro, A. L., & Bildsten, L. 2004, ApJ, 610, 977, doi: 10.1086/421763
  * Podsiadlowski et al. (2003) Podsiadlowski, P., Han, Z., & Rappaport, S. 2003, MNRAS, 340, 1214, doi: 10.1046/j.1365-8711.2003.06380.x
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82, doi: 10.1002/ctpp.201010017
  * Ramsay et al. (2018) Ramsay, G., Green, M. J., Marsh, T. R., et al. 2018, A&A, 620, A141, doi: 10.1051/0004-6361/201834261
  * Ritter (1988) Ritter, H. 1988, A&A, 202, 93
  * Rivera Sandoval et al. (2021) Rivera Sandoval, L. E., Maccarone, T. J., Cavecchi, Y., Britt, C., & Zurek, D. 2021, MNRAS, 505, 215, doi: 10.1093/mnras/stab1246
  * Rivera Sandoval et al. (2020) Rivera Sandoval, L. E., Maccarone, T. J., & Pichardo Marcano, M. 2020, ApJ, 900, L37, doi: 10.3847/2041-8213/abb130
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064, doi: 10.1086/341894
  * Rohatgi (2021) Rohatgi, A. 2021, Webplotdigitizer: Version 4.5.  https://automeris.io/WebPlotDigitizer
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Shen & Bildsten (2009) Shen, K. J., & Bildsten, L. 2009, ApJ, 699, 1365, doi: 10.1088/0004-637X/699/2/1365
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Townsley & Bildsten (2004) Townsley, D. M., & Bildsten, L. 2004, ApJ, 600, 390, doi: 10.1086/379647
  * Tutukov & Yungelson (1979) Tutukov, A. V., & Yungelson, L. R. 1979, Acta Astron., 29, 665
  * Tutukov & Yungelson (1981) —. 1981, Nauchnye Informatsii, 49, 3
  * van Roestel et al. (2021) van Roestel, J., Kupfer, T., Green, M. J., et al. 2021, arXiv e-prints, arXiv:2107.07573.  https://arxiv.org/abs/2107.07573
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Yungelson (2008) Yungelson, L. R. 2008, Astronomy Letters, 34, 620, doi: 10.1134/S1063773708090053