The double helium-white dwarf channel for the formation of AM CVn binaries

To investigate the stable mass transformation of two HeWDs, 16 pairs of HeWD models have been calculated on a grid with increments of 0.05 M_⊙ in mass, as shown in Figure 1.

Zoom In Zoom Out Reset image size

In our calculations, we treat the evolution of both WDs in detail, using the binary module of the stellar evolution code MESA Modules for Experiments in Stellar Astrophysics v8118; Paxton et al. 2011, 2013, 2015). This avoids the approximation of treating the accretor as a point mass.

To construct an HeWD, starting with a 1.5 M_⊙ zero-age main sequence star with metallicity of Z = 0.02, evolution is computed until the He core reaches the required mass, e.g. 0.30 M_⊙. Nucleosynthesis is switched off and a high mass-loss rate is applied to remove the hydrogen envelope. Low-mass HeWDs are considered to retain a small hydrogen envelope (Althaus et al. 2001; Steinfadt et al. 2010; Hall et al. 2013). Mass loss is therefore switched off when the hydrogen envelope is reduced to 10⁻⁴ M_⊙. The model evolves straight to the WD cooling track and does not ignite helium at this stage, i.e. log (L/L_⊙) = −2. This procedure produces the initial HeWD models which are used as input to the binary module; the evolution of both stars is then computed simultaneously, with mass transfer rates and orbital parameters computed self-consistently. The mass transfer is computed assuming optically thick overflow (Kolb & Ritter 1990), as also used for AM CVn calculations by Deloye & Taam (2006). In the optically thick overflow model, the mass transfer rate is sensitive to the difference between the stellar and Roche lobe radii, depending particularly on the detailed structure of the donor's outer layers. The pressure scale height H_p in the outer layer of the donor is much smaller than its radius R, i.e. H_p ≪ R.

Other details of the MESA physics are as follows. The ratio of mixing length to local pressure scale height is set to α = l/H_p = 1.9179, as found by the solar calibration of Paxton et al. (2011). The opacity tables are from Iglesias & Rogers (1996) and Ferguson et al. (2005). Because the abundances of carbon and oxygen in the interior change after an He flash, we use the OPAL Type 2 opacity tables. The outer boundary condition is chosen to be an Eddington gray photosphere. Orbital-decay by GW radiation has been switched on. We do not consider rotational, semiconvective or thermohaline mixing, or mass loss due to a stellar wind.

The birth rate of AM CVn systems from the double WD path is very sensitive to the synchronization timescale (τ_s) of the WD's spin with the orbital motion (Nelemans et al. 2001b). Using τ_s < 1000 yr, the birth rate can increase by more than a factor of 2 (Marsh et al. 2004). This is related to whether accretion is direct or mediated through a disk (Gokhale et al. 2007). A strong coupling, i.e. τ_s → 0, is used in our calculation in order to produce more AM CVn systems. We assume that mass transfer is conservative, so that the total mass and orbital angular momentum (apart from that lost by GW radiation) remains constant during evolution.

As the binary orbit decays due to GW radiation, 2×HeWD systems reach a semidetached state (in which the lower-mass component fills its Roche lobe), with minimum period of about 5 min. After mass transfer begins, evolution may proceed in one of two ways. For a high mass ratio system, the accretor will expand and fill its Roche lobe soon after accretion has begun. The two components will form a CE and then merge. For a low mass ratio system, the accretor will remain detached and stable mass transfer through the L₁ point will continue; such systems are identified as AM CVn binaries. From the MESA calculations, we make the assumption that a binary will merge when both stars completely fill their Roche lobes; the calculation is halted at this point. These are shown as solid triangles in Figure 1. For the models which do not merge, we compute the detailed evolution including stable Roche lobe overflow (RLOF: squares in Fig. 1) The calculation is stopped when the surface luminosity of the donor drops to log (L/L_⊙) = −5.

For example, the 0.40+0.25 M_⊙ model will merge and form a hot subdwarf (cf. Zhang & Jeffery 2012). However, the 0.40+0.15 M_⊙ model will survive and continue to transfer mass. The orbital period will increase from 5 min to ∼ 40 min over a timescale of 1 Gyr; increasing to ∼ 50 min after a further 2–3 Gyr. Although the orbital periods of semi-detached 2×HeWD binaries are small, the orbital expansion due to mass transfer is faster than contraction due to GW radiation. Thus, the orbit becomes wider during a long-lived evolution which can survive for a few Gyr. In summary, there are two different evolutionary directions during mass transfer in a 2×HeWD, i.e. to merge or to form an AM CVn binary. Those two directions depend on the initial HeWD masses and are shown in Figure 1.

There are some differences with previous calculations. A mass transfer rate = 10⁻⁵ M_⊙ yr⁻¹ was assumed by Nelemans et al. (2001b) as the boundary between stable and unstable mass transfer. Binaries with > 10⁻⁵ M_⊙ yr⁻¹ would have unstable mass transfer and then merge. For example Nelemans et al. (2001b) find that a 0.40 + 0.25 M_⊙ model does not merge but reaches stable RLOF because < 10⁻⁵ M_⊙ yr⁻¹. However, from detailed calculations for the same model (0.40 + 0.25 M_⊙) we find that, although the maximum mass transfer rate _max = 6.3 × 10⁻⁶ M_⊙ yr⁻¹, helium-burning in the accretor can be ignited. In this model, the accretor radius increases and the binary components merge. This result emphasizes the need for evolution of both accretor and donor to be included in the calculation.

To investigate AM CVn binaries formed through the 2×HeWD channel, we need to know the mass distributions of both HeWDs. We therefore compute the properties of a population of 10⁷ primordial binary systems using a Monte Carlo approach. Using the rapid binary evolution code (BSE, Hurley et al. 2000, 2002) and the technique of population synthesis, we identify those binaries which evolve to become 2×HeWDs with q < 2/3.

The BSE input parameters are the same as those used for modeling the Galactic rate of double WDs by Han (1998) and Zhang et al. (2014). We consider a single population of 10⁷ coeval binaries, and only count the pairs of 2×HeWDs which will merge or become semidetached. In the Monte Carlo simulation, all stars are assumed to be members of binaries and have circular orbits. The masses of the primaries are generated according to the formula of Eggleton et al. (1989), i.e. they follow the initial mass function of Miller & Scalo (1979) and are in the mass range 0.08 to 100 M_⊙. The secondary mass, also with a lower limit of 0.08 M_⊙, is subsequently obtained assuming a constant mass-ratio distribution. The distribution of orbital separations, p(a), is that of Han (1998)

$$\begin{eqnarray}p(a)=\left\{\begin{array}{lr}0.070{(a/{a}_{0})}^{1.2} & a\le {a}_{0},\\ 0.070 & {a}_{0}\le a\le {a}_{1},\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where a₀ = 10 R_⊙ and a₁ = 5.75 × 10⁶ R_⊙ = 0.13 pc. According to this, approximately 50 per cent of stellar systems have orbital periods greater than 100 rm yr (Han 1998) and are here considered to be single stars.

We find that these models occupy a region outlined in Figure 1. These results give a probability distribution of all possible mergers or semidetached HeWDs. The majority has initial masses in the range 0.35+0.2 – 0.47+0.28 M_⊙. Another concentration is found with masses in the range 0.3+0.14 – 0.4+0.16 M_⊙. The majority of WDs in those two regions is formed after two phases of CE ejection. However, binaries in the less massive zone require critical initial conditions. The two stars must be very close to each other. The more massive star evolves faster and forms an HeWD after the first CE ejection (red giant branch (RGB)+main sequence). The second CE phase (HeWD+RGB) must occur at the base of the RGB phase while the red giant has a very low-mass helium core, which then forms a second low-mass HeWD following CE ejection. Moreover, the first HeWD that forms must be massive enough to eject the high mass envelope of the RGB star. Hence, with a narrow region of initial mass ratios and separations, the number of 2×HeWDs in this subgroup is very small.

Eleven of the 16 original 2×HeWD models reach a stable RLOF phase. Of these, five are in the region where BSE population synthesis predicts 2×HeWDs to occur, namely the models with initial HeWD masses of 0.30+0.15, 0.30+0.20, 0.35+0.15, 0.35+0.20 and 0.40+0.15 M_⊙. For each theoretical system, we record the mass ratio, the donor and accretor masses, the orbital radial velocity and the GW strain amplitude as a function of orbital period, which is a proxy for elapsed time after the components come into contact. These will be compared with observed binary systems.

Figure 2 shows the evolution of the 0.35+0.15 M_⊙ model. The mass transfer rate increases very sharply at the beginning of RLOF, and then decreases exponentially as the orbital period increases. Meanwhile, the accretor mass increases as the orbital period increases.

Zoom In Zoom Out Reset image size

Mass is an important parameter for AM CVn systems. However, AM CVn star masses are difficult to measure because the systems are mostly faint and the orbits are short, making the components difficult to resolve. Mass determinations for most AM CVn stars depend on assuming the formation channel, and have been estimated by comparing observations with evolutionary models. We are particularly interested in those AM CVn systems which may have formed through the 2×HeWD channel. In order to compare our models with observations, data for eight AM CVn stars have been taken from Solheim (2010); Campbell et al. (2015) and Kalomeni et al. (2016). These are: AM CVn itself, HM Cnc, J0926, V406 Hya, J1240, GP Com, Gaia14aae and V396 Hya. All have relatively precise mass ratios from either Doppler tomography or radial velocity and eclipse measurements, but there are still no good measurements for their component masses.

Mass ratios have also been published for other AM CVn systems including CP Eri, HP Lib, CR Boo, KL Dra and V803 Cen. In these cases, the mass ratios were estimated from an empirical superhump period–mass ratio relation and have much larger systematic uncertainties.

For stable mass transfer in AM CVn systems, the donor's stellar radius R₂ is assumed to be equal to its Roche lobe radius, R_L, which can be approximated as (Paczyński 1967)

$$\begin{eqnarray}&&{R}_{{\rm{L}}}\approx 0.46a{\left(\frac{{M}_{2}}{{M}_{1}+{M}_{2}}\right)}^{1/3},\end{eqnarray} \tag{ 4 }$$

when q ≡ M₂/M₁ < 0.8 and the orbit is circular. Here a is the orbital separation, and M₁ and M₂ are the masses of the accretor and donor respectively. From Kepler's Third Law,

$$\begin{eqnarray}&&\frac{{P}^{2}}{{a}^{3}}=\frac{4{\pi }^{2}}{G({M}_{1}+{M}_{2})},\end{eqnarray} \tag{ 5 }$$

where P is the orbital period. Hence, equating R₂ = R_L,

$$\begin{eqnarray}&&P\propto {R}_{2}^{3/2}{M}_{2}^{-1/2}.\end{eqnarray} \tag{ 6 }$$

Since the radii of WDs can be approximately related to their masses, the donor mass M₂ is a function of orbital period.

Figure 3 shows that direct measurements of M₂ provide a strong constraint on the evolution models, since the latter require that the donor mass M₂ must be a function of orbital period. HeWDs must all follow approximately the same mass–radius relation, regardless of age. Thus, a single M₂–P relation correctly describes all five models computed, as seen in the figure.

Zoom In Zoom Out Reset image size

We obtain a polynomial fit to the theoretical M₂–P relation for the 2×HeWD channel

$$\begin{eqnarray}\begin{array}{lll}{\rm{log}}\left(\frac{{M}_{2}}{{M}_{\odot }}\right) & = & -1.7505\times {\rm{log}}\left(\frac{P}{{\rm{\min }}}\right)\\ & & +0.6659\times {\rm{log}}{\left(\frac{P}{{\rm{\min }}}\right)}^{2}\\ & & -0.2325\times {\rm{log}}{\left(\frac{P}{{\rm{\min }}}\right)}^{3}\\ & & +\mathrm{0.1758.}\\ & & \end{array}\end{eqnarray} \tag{ 7 }$$

For comparison, estimates of minimum donor masses have been taken from Deloye et al. (2005). The latter assumes a fully degenerate donor that fills its Roche lobe. As Deloye et al. (2005) note, this assumption may not always be satisfied since the thermal history (entropy) of the HeWD contributes strongly to the mass–radius relation used to derive the mass. Nevertheless, these masses were also used by Bildsten et al. (2006) to calculate the heating of WDs by accretion, demonstrating that the latter significantly contributes to the optical and ultraviolet emission of AM CVn stars.

Figure 3 demonstrates that the Deloye et al. (2005) minimum donor mass relation for AM CVn stars matches our 2×HeWD relation well. There remains a small difference between the relation of Deloye et al. (2005) and our models. This difference is related to the mass–radius relations adopted. The minimum donor masses from Deloye et al. (2005) are based on a mass–radius relation of HeWDs with central temperature log(T_c/K) = 4. In our models, the central temperature is higher than log(T_c/K) = 4 and varies with time during evolution, see the example shown in Figure 4. The higher internal temperatures will lead to larger radii and hence longer periods for a given mass.

Zoom In Zoom Out Reset image size

It is therefore suggested that our M₂–P relation and the mass ratio q could be used to estimate both minimum donor and accretor masses with good accuracy. Based on the M₂–P relation, the donor masses can be obtained from the orbital periods. Hence, the accretor masses may be calculated from the observed mass ratios. By this means we obtain masses for

HM Cnc (${0.289}_{-0.060}^{+0.101}$ + 0.144 M_⊙),

AM CVn (${0.218}_{-0.011}^{+0.013}$ + 0.039 0.039 M_⊙),

J0926 (${0.493}_{-0.042}^{+0.051}$ + 0.021 M_⊙),

V406 Hya (0.463 + 0.017 M_⊙),

J1240 (0.375 + 0.015 M_⊙),

GP Com (0.596 + 0.011 M_⊙),

Gaia 14aae (0.511 + 0.010 M_⊙) and

V396 Hya (0.490 + 0.006 M_⊙).

The major source of error is the observed mass ratio; only for HM Cnc, AM CVn and J0926 are errors given by Solheim (2010).

The donor's mass can be written as a function of mass ratio and total mass of the binary, i.e. M₂ = q/(1+q) × M_total. Hence, the mass ratio is a function of orbital period and total mass. The evolution of mass ratio q with orbital period P for each of the five selected models is shown in Figure 5. The model pairs of 0.30+0.20 and 0.35+0.15 M_⊙, and 0.35+0.20 and 0.40+0.15 M_⊙ both share the same relations, because of the same total masses. The mass ratios for seven of the eight AM CVn stars lie very close to the predicted q–P relation for stable RLOF in 2×HeWD stars, i.e. HM Cnc, J0926, V406 Hya, J1240, GP Com, Gaia14aae and V396 Hya. The AM CVn star was suggested to have a helium star donor (Solheim 2010).

Zoom In Zoom Out Reset image size

The stars J1240, GP Com and V396 Hya are classified as low-state AM CVn stars, which show almost no outbursts, except for J1240 (Levitan et al. 2015). They have very low mass transfer rates, in the range 10⁻¹³–10⁻¹² M_⊙ yr⁻¹. The principal difference between these three stars and other AM CVn stars is that their helium dominated spectra show lines with a triple-peaked profile. This is interpreted as a double-peaked accretion disk profile plus a sharp, low-velocity central component (Marsh 1999; Ruiz et al. 2001; Morales-Rueda et al. 2003; Roelofs et al. 2005; Kupfer et al. 2016). The origin of the sharp central component, or "central spike," is unclear. They have only ever been observed in some AM CVn systems and He-rich dwarf novae but never in hydrogen-dominated cataclysmic variables. The central spikes have been suggested to originate from the WD accretor or at least from near the accretor. Moreover, the central spike shows significant radial velocity shift as a function of orbital phase, which also indicates a likely origin close to the accretor, so that the semi-amplitude of the velocity shift should be equivalent to the orbital velocity of the accretor (assuming an orbital inclination close to 90^°).

During mass transfer in a binary system, stars exchange mass and angular momentum, which are here assumed to be conserved. The total angular momentum is

$$\begin{eqnarray}&&J={M}_{1}{M}_{2}{\left(\frac{Ga}{{M}_{1}+{M}_{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 8 }$$

The change of separation due to mass transfer is given by

$$\begin{eqnarray}&&\frac{\delta a}{a}=2\frac{\delta {M}_{2}}{{M}_{2}}\left(\frac{{M}_{2}}{{M}_{1}}-1\right).\end{eqnarray} \tag{ 9 }$$

In MESA, the orbital separation is calculated from the mass exchanged in each timestep of the evolution. Hence, the orbital periods are obtained from Kepler's laws; the orbital velocities follow directly.

Figure 6 shows the orbital velocities v of the accretor stars in our 2×HeWD models as a function of P. Observations generally provide the projected semi-amplitude, K = v sin i, where i is the inclination of the orbit with respect to the plane of the sky. The minimum radial velocities of the central spikes of three AM CVn systems are calculated with the estimated maximum inclination angles, i.e. 53°, 78° and 79° for J1240, GP Com and V396 Hya respectively. These agree well with the radial velocity of our accretormodels. However, to confirm whether central spikes originate from the WD accretor, more accurate observations are still required.

Zoom In Zoom Out Reset image size

3.3. Abundance

3.4. Gravitational Waves

AM CVn stars are an important source of GW emission and are expected to contribute strongly to the Galactic foreground detected by the Laser Interferometer Space Antenna (LISA) project (Nelemans et al. 2004). The amplitude of a GW signal emitted by a binary system in a circular orbit is given by

$$\begin{eqnarray}\begin{array}{lll}h & = & 0.5\times {10}^{-21}\times {(\frac{ {\mathcal M} }{{M}_{\odot }})}^{5/3}\\ & & \times {(\frac{P}{{\rm{h}}})}^{-2/3}\times {(\frac{d}{{\rm{kpc}}})}^{-1},\end{array}\end{eqnarray} \tag{ 10 }$$

where ${\mathcal M} ={({M}_{1}{M}_{2})}^{3/5}/{({M}_{1}+{M}_{2})}^{1/5}$ is the chirp mass, P is the orbital period and d is the distance to the system (Evans et al. 1987; Yu & Jeffery 2010). The GW frequency f = 2/P. To investigate the GW properties of our models, we take the 0.40+0.15 M_⊙ model as an example and calculate the GW strain amplitude for the model as it evolves from high frequency and high amplitude to low frequency and low amplitude (Fig. 7). Three representative distances are considered, i.e. 75 pc (the distance of GP Com), 500 pc (most observed AM CVn stars are closer than this distance) and 3 kpc (a few more distant AM CVn stars). The predicted negative shift in GWfrequency (chirp) due to orbital expansion in accreting systems is well known (Marsh et al. 2004; Kremer et al. 2017).

Zoom In Zoom Out Reset image size

We have recalculated the strain amplitude using the new masses of HM Cnc, AM CVn, J0926, V406 Hya, J1240, GP Com, Gaia14aae and V396 Hya predicted by our model (Section 3.1) and assuming lower limits on the distances of 1000, 513, 460, 100, 350, 75, 215 and 92 pc respectively (Solheim 2010; Campbell et al. 2015). The strain amplitudes are shown in Figure 7. With a one-year mission and significance threshold of 1σ, up to five AM CVns could be detected by LISA. Furthermore, more precise masses for the WD components could be obtained from extended LISA measurements. This is most likely for HM Cnc; the other AM CVn stars discussed here will be harder to resolve from the foreground of unresolved double WD (DWD) binaries (Nelemans et al. 2004; Yu & Jeffery 2010; Liu et al. 2010). For instance, the dot-dashed line shown in Figure 7 indicates the average foreground noise due to unresolved close WD binaries (Nelemans et al. 2004). AM CVn stars close to or beneath such a line will be difficult to resolve from such a foreground. Others, such as HM Cnc, AM CVn and V406 Hya, are potentially detectable above the LISA noise threshold and the DWD foreground.

In a recent study on the implications of accreting DWD binaries for LISA, Kremer et al. (2017) note that there are no observed DWDs with negative chirps. The best AM CVn candidate, HM Cnc, has a positive chirp suggesting that GW radiation dominates the mass transfer. Although Deloye & Taam (2006) argue for a low turn-on timescale for mass transfer which would allow the gravitational radiation to dominate, Kremer et al. (2017) find that the turn-on timescale should be very short (∼ 100 yr) and hence that the 2×HeWD model is implausible for HM Cnc.

Figure 8 shows the evolution of the chirp (rate of frequency change) and mass-transfer rate with period for the model of a 0.30+0.15 M_⊙ 2×HeWD, with masses believed to be similar to HMCnc. Combined with the observed orbital period of 5.4 min, we infer ≈ 10⁻⁷ M_⊙ yr⁻¹, which is similar to that inferred in previous work, e.g. Bildsten et al. (2006). However, the chirp is about half the value of 3.63±0.06 × 10.16Hz s⁻¹ reported by Strohmayer (2005). Figure 9 shows the chirp and frequency evolutionwith time. The turn-on timescale is about 2.5 × 10⁴ yr. A possible region for HMCnc is in between two vertical dashed lines, within which the period (5.4 min) is constant and the chirp is positive. With such characteristics, LISA observations will be crucial in resolving questions about the formation channel.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have computed the evolution of 2×HeWD binaries covering a wide mass range with mass ratios q < 2/3, and including the detailed evolution of both stars. Two different evolutionary pathways were considered, i.e. for low-mass ratio pairs to show stable RLOF over several Gyr, and for high-mass ratio pairs to form mergers. The response of the accretor to its increasing mass has a modest effect on the boundary between RLOF pairs and mergers, but otherwise we have found negligible effect on the overall evolution of the system.

The value q_crit = 2/3 was calculated assuming a simple mass–radius relation for the HeWDs. Generally speaking, such relations are derived for coldWDs. There are likely to be instances when the HeWD structure departs significantly from this condition. Hence, some pairs of 2×HeWDs with q > 2/3 may evolve through stable RLOF and do not merge. However, the most recent version of MESA does not provide a stable calculation for such models. This will be examined carefully in future work. Furthermore, we have not yet excluded He star+CO WD or He+CO WD binaries as progenitors for these systems; this will require additional models.

From the models, an orbital period–donormass relation has been derived and used to estimate the component masses for eight AM CVn stars. From radial velocity arguments, the models provide evidence that central spikes in the triple-lined emission spectra of J1240, GP Com and V396 Hya originate from or close to the accretor in 2×HeWD systems. In terms of surface composition, our 2×HeWD models show surface N/C and N/O ratios very similar to GP Com, the only AM CVn in which such measurements have been made.

As AM CVn stars are among the most likely resolvable sources of Galactic GWradiation, we recalculate the predicted signal using our new mass estimates for the eight suspected 2×HeWD AM CVn systems. We confirm that the stars HM Cnc, AM CVn and V406 Hya are potentially detectable by LISA; the GWstrain-amplitude would confirm the masses and hence possible origin for these systems.

Future work will examine the boundary between merger and stable RLOF channels in more detail, address its impact on the predicted number densities for AMCVn stars in the solar neighborhood, and compare with models for stable RLOF in He+CO WD binaries.