PRELUDE TO A DOUBLE DEGENERATE MERGER: THE ONSET OF MASS TRANSFER AND ITS IMPACT ON GRAVITATIONAL WAVES AND SURFACE DETONATIONS

White dwarfs (WDs) are the end point of the evolution of the vast majority of the stars in the universe. Our Galaxy harbors of order 10¹⁰ WDs (Napiwotzki 2009) and about 2.5 × 10⁸ close binary WDs (Nelemans et al. 2001b), often referred to as double degenerates (DDs). For about half of them, a Hubble time is long enough so that gravitational wave emission can drive them to a phase of mass transfer. Such interacting DDs are suspected to produce observable objects such as sdb stars (Webbink 1984; Iben & Tutukov 1986; Saio & Jeffery 2000; Han et al. 2003), R Corona Borealis stars (Webbink 1984; Clayton et al. 2007), AM CVn stars (Paczyński 1967; Nelemans et al. 2001a; Solheim 2010; Kilic et al. 2011), and maybe most spectacularly, Type Ia supernovae (Iben & Tutukov 1984; Webbink 1984). While for many years there was a strong inclination toward the single degenerate scenario, the recent discovery of "super-Chandrasekhar" events (Howell et al. 2006; Hicken et al. 2007; Yuan et al. 2007; Silverman et al. 2011) has re-strengthened the interest in DD mergers as possible Type Ia progenitors. The lead-up to a merger includes a phase of mass transfer that is crucial for understanding the evolution and fates of double WD systems. This phase will affect the orbital evolution, the angular momentum and thermodynamical state of the final merger product, the associated gravitational wave signal, the environment into which any potential ejecta produced by an thermonuclear explosion will expand into, and whether or not a particular binary will merge at all. The numerical study of the onset of this mass transfer phase is the main topic of this investigation.

The merger process of two WDs has been modeled by a number of groups (Benz et al. 1990; Rasio & Shapiro 1995; Segretain et al. 1997; Guerrero et al. 2004; Yoon et al. 2007; Pakmor et al. 2010), all of which use the smoothed particle hydrodynamics (SPH) method. In all of these simulations the mass transfer phase lasted for very few orbital timescales. Recently, grid-based simulations (Motl et al. 2002; D'Souza et al. 2006; Motl et al. 2007) were performed in which the authors carefully tried to reduce the angular momentum non-conservation due to advection errors that often plague grid-based codes (New & Tohline 1997). The latter results differed from the previous SPH calculations in that they showed long-lived mass transfer for many orbital periods and this has sparked some discussion about the stability of mass transfer in DD systems during these late stages. Based on our experience in modeling the dynamics of neutron star black hole systems (Rosswog et al. 2004; Rosswog 2005), we had suspected that double WD binaries may be very sensitive to the exact initial conditions (ICs). As we will demonstrate below, the ICs used have a decisive impact and are one of the major reasons for the disagreement between various mass transfer calculations.

This is the companion paper of Guillochon et al. (2010), in which we focused on the fate of accreted helium in unstable, direct impact systems. The work showed that, contrary to prior belief, the mass transfer phase between two WDs may be eventful and can be accompanied by Kelvin–Helmholtz-triggered thermonuclear surface explosions. These explosions may lead to the ignition of the accreting CO WD in a double detonation scenario, which could be the long sought-after mechanism to initiate a Type Ia supernova from WD binary merger. For the described work, we had chosen a hybrid approach by combining results using the SPH method for the orbital evolution and the adaptive mesh refinement (AMR) code FLASH (Fryxell et al. 2000) to focus on the evolution of the low-density helium material. We consider this a "best-of-both-worlds" approach because it combines the major strengths of both Lagrangian and Eulerian methods. In this paper we present the details of the mass transfer/orbital evolution calculations that have been obtained with SPH, as well as a detailed comparison with the previously published FLASH calculations.

The main focus of this paper is the onset of the mass transfer and its impact on the orbital evolution; the structure of the resulting remnants will be discussed elsewhere. In Section 2 we briefly summarize conditions for the stability of mass transfer, and in Section 3 we focus on the numerical aspects of the mass transfer modeling. We describe in detail how we construct our ICs and demonstrate their impact on the outcome of essentially every aspect of DD merger simulations. We also briefly summarize the essential features of the SPH code that we are using. In a set of test calculations, we explore to which extent the numerical resolution, the form of artificial viscosity, and the exact numerical ICs affect the results. In Section 4, we highlight results from two systems representative of a set of binary evolution calculations. In Sections 5 and 6, we discuss the implications of the extended mass transfer phase on the gravitational wave signal and the triggering of surface detonations prior to merger, respectively. We conclude the paper in Section 7 with a summary.

In what follows we briefly collect the standard arguments on the stability of mass transfer (see, e.g., Frank et al. 2002). Once stars are close enough that mass transfer via Roche lobe overflow can set in, the further evolution of the binary system depends on how both the orbit/Roche lobe size and the star react to the re-arrangement of mass and angular momentum. Due to their particular mass–radius relation, WDs react to mass loss by expansion which tends to enhance the mass transfer rate. If mass is transferred to the heavier star, the orbit has to widen to ensure the conservation of the center of mass and, therefore, mass transfer is tendentially quenched. The interplay of both effects determines the detailed orbital evolution. The total angular momentum of a point mass binary is given by

$$\begin{equation} J= M_1 M_2 \left(\frac{G a}{M_{\rm tot}}\right)^{1/2}, \end{equation} \tag{ 1 }$$

where the M_i are the masses, G is the gravitational constant, M_tot is the total mass, and a is the separation of the binary system. In the remainder of this work, we will always denote the donor by label 2, the accretor by 1. On logarithmic differentiation one finds

$$\begin{equation} \frac{\dot{a}}{2 a}= \frac{\dot{J}}{J} - \frac{\dot{M_2}}{M_2} \left(1 - q \right), \end{equation} \tag{ 2 }$$

where q = M₂/M₁. As expected, angular momentum loss tends ($\dot{J}<0$) to shrink the orbit, while mass transfer ($\dot{M}_2<0$) tends to increase the orbital separation. The change of the donor Roche lobe size can be related to the change in orbital separation by logarithmically differentiating a simple Roche lobe size estimate (Paczyński 1971)

$$\begin{equation} R_{\rm L,2}= 0.462 \; a \; \left(\frac{q}{1+q}\right)^{1/3}, \end{equation} \tag{ 3 }$$

which yields

$$\begin{equation} \frac{\dot{R}_{\rm L,2}}{R_{\rm L,2}}= \frac{\dot{a}}{a} + \frac{1}{3} \frac{\dot{M}_2}{M_2}. \end{equation} \tag{ 4 }$$

Using this to eliminate $\dot{a}$ from Equation (2) then yields

$$\begin{equation} \frac{\dot{R}_{\rm L,2}}{R_{\rm L,2}}= \frac{2 \dot{J}}{J} + \frac{2 \dot{M}_2}{M_2} \left(q-\frac{5}{6}\right), \end{equation} \tag{ 5 }$$

and therefore for a mass ratio q > 5/6 a non-zero mass loss will trigger an expansion of the donor WD and at the same time shrink the Roche lobe size R_{L, 2}. This will be further accelerated by any additional angular momentum loss, so with the used assumptions the mass transfer is clearly unstable in such cases.

Of course, the case discussed so far is highly idealized. A point mass binary was assumed where the mass lost by the donor was completely incorporated into the accretor, while in reality mass can be lost from the system. Depending on whether or not the circularization radius exceeds the accretor radius, the transferred mass can either form a disk or it can impact directly onto the accretor surface. If a disk is formed, it would be tidally perturbed by the companion star and can thus feed back some fraction of the angular momentum into the orbit. In the direct impact case, most of the angular momentum is invested into spinning up the accretor, which removes orbital angular momentum from and is therefore expected to destabilize the orbit. It should be noted, however, that a deformed donor can also feed back some angular momentum into the orbit. The orbit and the different binary components can exchange angular momentum both by advection of mass and by tides. Marsh et al. (2004) have generalized the simple treatment leading to Equation (2) for direct impact and tidal spin evolution,

$$\begin{equation} \frac{\dot{a}}{2 a}= \frac{\dot{J}}{J_{\rm orb}} - \left(\frac{\dot{J}_1 + \dot{J}_2}{J_{\rm orb}}\right)_{\rm tid} - \frac{\dot{M_2}}{M_2} \big[ 1 - q - \sqrt{(1+q) r_{\rm c}}\big], \end{equation} \tag{ 6 }$$

where r_c is the circularization radius in units of the orbital separation. The second term accounts for the exchange of angular momentum due to tides, and the third term accounts for the spin of the accretor. If tides and accretor spin are ignored, Equation (6) reduces to Equation (2). In Figure 1, we show the stability criteria for binaries of different donor and accretor masses, superimposed with the initial masses of our simulations (Table 1). Marsh et al. (2004) find the upper left of the diagram to be guaranteed unstable, the lower right guaranteed stable, and an uncertain region along the diagonal of the diagram. In this work we mainly investigate the "uncertain" region in the diagram.

Zoom In Zoom Out Reset image size

The main tool of this investigation is a three-dimensional SPH code. For recent reviews of the method see, e.g., Monaghan (2005) and Rosswog (2009). The orbital dynamics of a binary system is very sensitive to the (re-)distribution of angular momentum, and therefore a numerical simulation has to ensure very accurate angular momentum conservation. Its purely Lagrangian nature and the built-in conservation of mass, energy, and linear and angular momenta (e.g., Section 2.4 in Rosswog 2009) make SPH a very powerful tool for this type of study. The conservation is exact up to effects that result from finite accuracy in gravitational forces due to the use of a tree (see below) and in the numerical integration of the ordinary differential equation (ODE) set. Both issues, however, are fully under control and can be improved to a desired level of accuracy by choosing the corresponding parameters for the tree and the ODE integration in an appropriate way, at the expense of a larger computational effort, though. For our parameters, even in the 84 orbit P1, corresponding to as many as 16,930 dynamical timescales of the 0.8 M_☉ star, energy and angular momentum are conserved to better than 1%. It has to be stated, however, that the advantage of exact mass conservation via constant SPH particle masses turns into a disadvantage if we are interested in the resolution of low-density portions of the flow. That is why we turned to the AMR code FLASH in Guillochon et al. (2010) to investigate the settling of the transferred material onto the accretor.

Details of our SPH code can be found in Rosswog et al. (2008). It uses an artificial viscosity scheme (Morris & Monaghan 1997) that reduces the dissipation terms to a very low level away from discontinuities, together with a switch (Balsara 1995) to suppress the spurious viscous forces in pure shear flows. The system of fluid equations is closed by the Helmholtz equation of state (EOS; Timmes & Swesty 2000). It accepts an externally calculated nuclear composition and allows a convenient coupling to a nuclear reaction network. We use a minimal nuclear reaction network (Hix et al. 1998) to determine the evolution of the nuclear composition and to include the energetic feedback onto the gas from the nuclear reactions. A set of only seven abundance groups greatly reduces the computational burden, but still reproduces the energy generation of all burning stages from He burning to nuclear statistical equilibrium accurately. We use a binary tree (Benz et al. 1990) to search for the neighbor particles and to calculate the gravitational forces.

3.1. Initial Conditions

In what follows we will summarize how ICs are constructed for double WD binary systems, both in this and in previous work.

3.1.1. Previous Work

In their pioneering work, Benz et al. (1990) followed the dynamical evolution of a 0.9 + 1.2 M_☉ WD binary system. They constructed cold, isolated WDs in hydrostatic equilibrium and placed them without spin-rotation at an initial separation such that the secondary's Roche radius was smaller than its unperturbed initial radius by 8%. One of the major conclusions of Benz et al. was that within about two orbital periods the secondary is completely destroyed and transformed into a thick accretion disk around the accretor. As we will show below, this short time to merger is an artifact of the initial separation being too small. We have followed the recipe of Benz et al. for the ICs and find exactly the same result as presented in their paper. We have further compared the orbital evolution of a 50,000 particle binary system, once according to the original Benz et al. description and according to our scheme, as described in Section 3.1.2. Again, with their ICs the merger set in within about one orbital period. With our ICs and everything else being the same, the binary survived for 14 orbits of continued mass transfer.

In Segretain et al. (1997) and Yoon et al. (2007), the stars were placed at a mutual separation so that the secondary fills its Roche lobe. In Segretain et al. the WDs had no initial spin, while in Yoon et al. the WDs begin tidally locked. Guerrero et al. (2004) and Lorén-Aguilar et al. (2009) initially placed the stars in a detached configuration and then slowly decreased the separation to the point where mass transfer began. However, their relaxation procedure did not allow the stars to reach complete hydrostatic equilibrium. If the stars are brought together too quickly, this can cause undesired (and unrealistic) oscillations in the WDs and may again lead to an abbreviated mass transfer phase. For example, the mass transfer phase in their 0.4 + 1.2 M_☉ system (Figure 10 in Guerrero et al. 2004) only lasts for a single orbit. This is also the case for the simulation of a 0.6 + 0.8 M_☉ system (see Figures 1 and 2 in Lorén-Aguilar et al. 2009).

Zoom In Zoom Out Reset image size

Our ICs are constructed in a way that is similar to Rasio & Shapiro (1995; Section 3.1.2). Rasio & Shapiro constructed a synchronized binary system with mass ratio of q = 0.5 and also followed the dynamical evolution once a few particles had crossed the L₁ point. They found a merger after five orbital periods and also questioned the earlier use of non-equilibrium ICs. We have tried to reproduce their calculations as closely as possible by adopting the same particle number, EOS, and the same prescription for evolving the entropic function K(s) (Equation (12)), but still we only found a merger after 15 orbital periods. Fryer & Diehl (2008) place the two stars close enough such that the donor is already overfilling its Roche radius, which leads to a premature merger. The recent work of Pakmor et al. (2010) focuses primarily on the coalescence of the binary, but their ICs are not appropriate for studies of long-term stability.

3.1.2. Our Approach

We primarily focus on synchronized binary systems, but we also perform some tests with initially non-spinning WDs. Our procedure is similar to the one described in Rosswog et al. (2004). We begin our initialization procedure with stars that are set up according to one-dimensional, cold (T = 10⁵ K) equilibrium models. These stars are then mapped into three dimensions and are then driven to an accurate numerical equilibrium by applying a velocity-dependent acceleration $\vec{f}_{\rm damp}$ while keeping the temperature constant. The stars are then placed at an initial separation that is large enough to avoid any immediate mass transfer. The full SPH code is then used to relax the ICs to a tidally deformed binary system at the brink of numerically resolvable mass transfer. This phase of the initialization is performed in a corotating frame, i.e., centrifugal and Coriolis forces must be accounted for.³ We calculate the orbital frequency ω in a way that the total forces on the individual stellar centers of mass vanish exactly (Rosswog et al. 2004) and thus account for finite size deviations from the Kepler frequency. The acceleration applied to a particle i then reads

$$\begin{eqnarray} {\bm f}_i&=& {\bm f}_{{\rm grav}, i} + {\bm f}_{{\rm hyd}, i} - {\bm \omega } \times \left({\bm \omega } \times {\bm r}_i\right) - 2 {\bm \omega } \times {\bm v}_i\nonumber\\ && - \frac{{\bm v}_i - {\bm \Omega }_{1,2} \times ({\bm r}_{1,2} - {\bm r}_i)}{\tau _i}. \end{eqnarray} \tag{ 7 }$$

Here, ${\bm f}_{{\rm grav}, i}$ and ${\bm f}_{{\rm hyd}, i}$ are the gravitational and hydrodynamic acceleration, ${\bm r}_i$ and ${\bm v}_i$ are the position and velocity in the corotating frame, ${\bm \Omega }_{1,2}$ is either the primary or the secondary star's angular velocity, and τ_i is a damping timescale. In this frame, ${\bm \Omega }_{1,2} = 0$ for synchronized systems and ${\bm \Omega }_{1,2} = -{\bm \omega }$ for non-spinning stars. During this relaxation process, we adiabatically decrease the orbital separation a in a way that the orbital shrinkage timescale, $\tau _{\rm shrink}=a/\dot{a}$, is substantially longer than the dynamical timescale of the secondary τ_{dyn, 2},

$$\begin{equation} \tau _{\rm shrink} = \frac{a}{\dot{a}} = \frac{1}{\epsilon \sqrt{G \bar{\rho }_2}} \gg \tau _{\rm dyn, 2}, \end{equation} \tag{ 8 }$$

where $\bar{\rho }_2$ is the average density of the secondary. We found that spurious oscillations are minimized for values of ≲ 0.05. As soon as the first SPH particle crosses the Lagrange point L₁, we stop the relaxation process, transform to a space-fixed frame with the same coordinate origin, and start the dynamical simulation. This moment is adopted as the time origin (t = 0) of the simulation, shown in Figure 2 for a 0.9 + 1.2 M_☉ system. In Figure 3, we show the potential of the ICs, where Φ_pm is the usual Roche potential as calculated for a point mass binary,

$$\begin{equation} \Phi _{\rm pm}({\bm r})= - \frac{{\it GM}_1}{|{\bm r} - {\bm r}_1|} - \frac{{\it GM}_2}{\left|{\bm r} - {\bm r}_2\right|} - \frac{1}{2}\left({\bm \omega }_{\rm k} \times {\bm r}\right)^2, \end{equation} \tag{ 9 }$$

with M₁/M₂ and ${\bm r}_1/{\bm r}_2$ the stellar masses and center of mass positions and ${\bm \omega }_{\rm k}$ is the point mass Kepler frequency. The quantity Ψ is the corresponding quantity as calculated for the actual fluid configuration,

$$\begin{equation} \Psi ({\bm r})= \Phi \left({\bm r}\right) - \frac{1}{2}\left({\bm \omega } \times {\bm r}\right)^2, \end{equation} \tag{ 10 }$$

where ${\bm \omega }$ is the orbital frequency of our (fluid) binary and Φ is the gravitational potential due to the fluid.

Zoom In Zoom Out Reset image size

As we will show below in detail, the accuracy of the ICs has a strong impact on all major aspects of the simulations, prolonging the mass transfer stage by dozens of orbits in all investigated cases.

3.2. Numerical Treatment

In this section we explore the sensitivity of the binary evolution on physical and numerical factors, including numerical resolution, shock heating treatment, the artificial viscosity prescription, and the ICs.

3.2.1. Dependence on Numerical Resolution

The duration of mass transfer and the separation at which it begins depend on the extent to which we can resolve the outer layers of the WD donor. In SPH, where usually the particle masses are not evolved in time, this means that even in a decent numerical effort with a million SPH particles the initial resolvable mass transfer rate is already a substantial fraction of the total system mass. The numerically resolvable mass transfer $\dot{M}_{\rm lim}$ can be estimated as

$$\begin{eqnarray} \dot{M}_{\rm lim}\approx \frac{1\ {\rm particle \; mass}}{{\rm orbital \; period}} &\approx& 2 \times 10^{-8} \left(\frac{N_{\rm p}}{10^6}\right)^{-1} \left(\frac{M_{\rm tot}}{1\, M_\odot }\right)^{3/2}\nonumber\\ &&\times \left(\frac{a_0}{2 \times 10^9\, {\rm cm}}\right)^{3/2}\frac{M_\odot }{\rm s}, \end{eqnarray} \tag{ 11 }$$

where a₀ is the separation between the stars at the onset of mass transfer. In order to better resolve regions of lower density, SPH particles of different masses could be used, but we refrain from such a possibility as it can introduce a substantial level of numerical noise.

3.2.2. Comparison with the StarCrash Code

As a validation of our code, we have also performed a comparison with the results obtained with the StarCrash code. This code was developed originally by Rasio and Shapiro to study merging binaries (Rasio & Shapiro 1992, 1994, 1995). The major conceptual difference in their code is the use of a fast Fourier transform instead of a tree to calculate the gravitational forces.

To enable a valid comparison, for this test we also use a polytropic EOS, P = Kρ^5/3, with Γ = 5/3, we evolve the entropic function according to

$$\begin{equation} \frac{{\it dK}_i}{dt}=\frac{\Gamma -1}{2\rho ^{\Gamma -1}_i}\sum _j m_j\Pi _{ij} {\bm v}_{ij}\cdot \nabla _i W_{ij}, \end{equation} \tag{ 12 }$$

where m_j is the mass of particle j, Π_ij is the artificial viscosity tensor depending on parameters α and β, ${\bm v}_{ij} \equiv {\bm v}_i-{\bm v}_j$ is the velocity difference between particles i and j, and W is the smoothing kernel. We fix our artificial viscosity parameters in this comparison to the values recommended in their StarCrash documentation (α = 0.5, β = 1). We use the subroutines of the StarCrash code to set up a synchronized binary system of mass ratio q = 0.5 with 40,000 SPH particles; both simulations start from these ICs. The codes yield nearly indistinguishable results, and after about 15 orbital periods the donor WD becomes tidally disrupted. Although Rasio & Shapiro (1995) also utilize the StarCrash code, they start with slightly different ICs, with the two stars being closer to one another at t = 0. As a result, they find that the donor only survives five orbits, despite using the same EOS and SPH particle number. We suspect that they started from a slightly too small initial separation, but the very good agreement of the two codes for the same ICs is encouraging. In passing we note that these results are very similar to those D'souza et al. (2006), which will be discussed in more detail in Section 4.4.

3.2.3. Effect of Shock Heating Treatment

As mentioned previously, recent grid-based simulations (Motl et al. 2002; D'souza et al. 2006; Motl et al. 2007) showed long-lived mass transfer phases, in contrast to the earlier SPH simulations. In these calculations artificial viscosity was only used in the momentum equation, but not in the energy equation, which means that entropy production in shocks was neglected. It was suspected that at least part of the difference in the mass transfer duration may be due to the suppression or emergence of shocks/ignoring the change of entropy (Fryer & Diehl 2008).

We perform two test calculations at the example of a q = 0.5 binary, similar to D'souza et al. (2006). We follow their approach and use a Γ = 5/3 polytrope, but in one case we evolve the entropy function according to Equation (12) and in the other we keep K constant. A mass ratio of 0.5 results in a direct impact of the accretion stream onto the accretor surface. In the setup that does not include shocks, the accreted matter is smoothly incorporated into the accretor (Figure 4, lower row), while the setup that includes shocks produces a high-entropy halo around the accretor (Figure 4, upper row). This high-entropy halo applies a greater torque on the binary as it has a larger effective lever arm than the low-entropy halo that is produced when shocks are neglected. As a result, the matter within the high-entropy halo is able to remove more angular momentum from the orbit over the same number of orbital periods (Figure 5), leading to a more rapid evolution of the binary separation. This supports the findings of Fryer & Diehl (2008) that the large increase in orbital separation reported by D'souza et al. (2006) is in part an artifact of shock suppression.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.4. Impact of the Artificial Viscosity Prescription

As we have shown that the entropy generation associated with shocks can affect the stability of DD systems, we need to consider the importance of other entropy-generating mechanisms. To this end we compare different implementations of artificial viscosity. All have the "standard" form originally suggested by Monaghan & Gingold (1983), but differ in terms of constant or non-constant parameter values and additional "limiter switches:"

• 1.  
  AV₁: constant parameters α = 1 and β = 2;
• 2.  
  AV₂: constant parameters α = 2 and β = 4;
• 3.  
  AV₃: time-dependent parameters (Morris & Monaghan 1997), controlled as described in Rosswog et al. (2008), but without "Balsara switch" (Balsara 1995); and
• 4.  
  AV₄: time-dependent parameters plus Balsara switch as described in Rosswog et al. (2008).

As a test case we run simulations with a 0.3 He WD and a 0.6 CO WD, each one modeled by 20,000 particles (T4–T7, Table 1). Synchronized ICs are constructed as described in Section 3.1.2 and all tests used the Helmholtz EOS. Artificial viscosity does indeed have an impact: in the lowest viscosity case, AV₄, the system merges after 19 orbital periods; in the highest viscosity case, AV₂, after 22 orbital periods. To quantify the dissipation we calculate the thermal energy produced by artificial viscosity before the merger,

$$\begin{equation} E_{\rm diss}= \sum _i m_i \int _0^{t} \left(\frac{du_i}{dt}\right)_{\rm diss} dt^{\prime }, \end{equation} \tag{ 13 }$$

where

$$\begin{equation} \left(\frac{du_i}{dt}\right)_{\rm diss}=\frac{1}{2}\sum _j m_j\Pi _{ij}{\bm v}_{ij}\nabla _i W_{ij} \end{equation} \tag{ 14 }$$

is the artificial viscosity contribution to the energy equation of particle i. The quantity E_diss as a function of time (for t < t_disr, where t_disr is the time when the donor star is tidally disrupted) is shown in Figure 5 (right panel), for all four artificial viscosity prescriptions. Interestingly, the prescription with the lowest average artificial dissipation parameters, AV₄, dissipates the most energy. The reason for this is that the combination of low-viscosity parameters and low resolution allow the SPH particles to build up a large level of velocity noise (i.e., fluctuations around the local average value). This leads to an increase in the typical value of ${\bm v}_{ij}$, and thus more dissipation, despite the lower viscosity parameters and suppression of dissipation in shear flows.

The increased rate of dissipation results in a more efficient spin-up of the accretor, with the angular momentum for this spin-up coming at the orbit's expense. This leads to a reduction in the number of orbital periods the binary can survive prior to merger. We perform our production runs with prescription AV₄ which yields a lower limit on the duration of the mass transfer phase.

3.2.5. Dependence on Initial Conditions

As we have discussed, most of the earlier WD merger simulations were constructed using approximate ICs where the separation was estimated by comparison of the stellar radius with the Roche lobe size. To illustrate the importance of carefully constructed ICs, here we compare the evolution of a binary initialized using both "approximate" and "accurate" ICs.

To construct the approximate ICs, we equate Eggleton's formula (Eggleton 1983) for the Roche radius with the radius R₂ of the secondary as found from the SPH particle positions. The initial separation a₀ is thus given by

$$\begin{equation} a_0^{\rm Egg}= R_2 \; \frac{0.6 q^{2/3}+\ln (1+q^{1/3})}{0.49q^{2/3}}. \end{equation} \tag{ 15 }$$

The stars, previously relaxed in isolation (i.e., without the companion's tidal field), are subsequently placed at this mutual separation with a spin that has the same period as the orbit of the binary. We refer to these ICs as "approximate" and compare them to the carefully constructed ICs described in Section 3.1.2. The described procedure is plausible, as we will see below, and its approximate nature introduces a slew of artifacts.

Two representative examples of binaries are compared, a 0.2 M_☉ He + 0.8 M_☉ CO binary and a 0.6 M_☉ CO + 0.9 M_☉ CO binary. Each system is simulated using both the approximate (corresponding to runs T1 and T2 in Table 1 and Figure 1) and accurate (P1 and P5) ICs. Some marked differences are seen between the simulations initialized using these two different prescriptions. The approximate prescription systematically underestimates the separation at which mass transfer sets in by about 14% in both cases. This is mainly because it neglects the tidal deformation of the donor star. As a result, binaries initialized using approximate ICs carry a deficit of angular momentum when compared to those evolved using the more accurate prescription. The difference in angular momentum between the two IC prescriptions is 11% and 13% for the 0.2 + 0.8 M_☉ and 0.6 + 0.9 M_☉ cases, respectively. This slight discrepancy severely impacts their subsequent evolution, with the duration of the mass transfer phase being underestimated by about two orders of magnitude. In the 0.2 + 0.8 M_☉ case, the approximate ICs yield two orbits (478 s) of mass transfer while the accurate ICs show a mass transfer for as long as 84 orbital periods (20,328 s). The difference in the 0.6 + 0.9 M_☉ case is similarly dramatic (1 versus 29 orbital periods).

The effect of varying ICs on the orbital dynamics is most clearly seen when comparing the resulting gravitational wave signals (Figures 6 and 7). The approximate ICs give rise to a rapidly fading signal, while in the more accurate case the amplitude remains essentially constant up to the point of tidal disruption. The 0.6 + 0.9 M_☉ case, which is less stable than the 0.2 + 0.8 M_☉ system (Figure 1), shows only a slightly noticeable increase in the gravitational wave amplitude up to disruption. In both cases, the approximate ICs overestimate the peak amplitude and frequency. The approximate case gives an orbital frequency of 46 mHz in the 0.6 + 0.9 M_☉ case and 10 mHz in the 0.2 + 0.8 M_☉ case, while the accurate ICs yield 38 and 8.4 mHz, respectively (a difference of about 20%). This is mainly due to the difference in separation at final plunging (40% difference for the 0.2 + 0.8 M_☉ case and 30% for the 0.6 + 0.9 M_☉ case).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a result of the more abrupt artificial decrease in separation, the approximate ICs also give rise to large spurious oscillations of the central remnant core. Such oscillations are not only imprinted on the gravitational wave signal (see Figure 7) but are also present on the thermodynamical trajectories of the merger products. As shown in Figures 8 and 9, both peak densities and temperatures are overestimated when using the approximate ICs. The use of accurate ICs helps to better characterize the time evolution of the temperature, density, and angular momentum during mass transfer as well as the final morphology of the merger product, as seen in Figure 10. Due to the larger initial separation for the more accurate ICs, the matter possess more total angular momentum, which results in more mass in the trailing arm and less mass in the disk than when starting from the approximate IC. The material ejected on eccentric orbits will fall back later onto the remnant (similar to the case of neutron star mergers; see, e.g., Rosswog 2007; Lee & Ramirez-Ruiz 2007), the timescale of this fallback matter being seriously underestimated when using the approximate ICs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The use of accurate ICs and a correct treatment of shock heating is necessary for a reliable temporal evolution of the temperature, density, and angular momentum during the mass transfer phase. The evolution of two representative binary systems is described here in detail: the merger of a synchronized binary system with a 0.2 He WD and a 0.8 CO WD (P1, Table 1) and a tidally locked binary with 0.6 + 0.9 M_☉ CO WDs (P5). Their evolution is then subsequently compared with the wide range of binary systems that we have simulated as well as with the results of earlier work.

4.1. 0.2 + 0.8 M_☉: A Borderline Unstable System Near the Disk Stability Limit

According to the stability analysis of Marsh et al. (2004), a 0.2 He + 0.8 CO M_☉ binary system (marked as "P1" in Figure 1) should still be in the "direct impact" regime where the pericenter separation of the transferred mass is smaller than the accretor radius, but near the stable, disk-forming regime. This is clearly seen in simulation snapshots of the matter distribution plotted in Figure 11. In this case, we find that the rotation of the donor star remains synchronized with the orbital motion throughout the mass transfer phase up to the point of tidal disruption. Due to a not completely accurate numerical relaxation process, there is a slight deviation from the angular frequency that exactly produces a circular orbit. This oscillation, as expected, also modulates the mass transfer rate (Figure 12).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The average orbital separation in this system is observed to secularly increase and, as a result, the mass transfer is stabilized, similar to what is expected for AM CVn stars. The transfer rate stays nearly constant at a level of ≈10^6.2 $\dot{M}_{\rm Edd}$ for about 50 orbital periods and continues until the donor star is completely disrupted. This occurs after as many as 84 orbital or 16,930 dynamical timescales of the accretor. Before being disrupted, the donor star looses about 0.08 M_☉. Our results are in agreement with the Marsh et al. analysis, which suggests the mass transfer phase to be long lived but still unstable.

The binary system during the mass transfer phase sheds matter through the outer Lagrange points, resulting in the loss of approximately 4% of the total angular momentum. As we have argued, the orbital dynamics is very sensitive to the loss of angular momentum which, in turn, prevents efficient mass transfer stabilization. The orbital angular momentum that is carried away is shared between the accreting star and the disk (Figure 13), and, due to the constant change in orbital separation, is observed to oscillate (Figure 13).

Zoom In Zoom Out Reset image size

The densities and temperatures of the final remnant core are shown in Figure 14 for the 0.2 + 0.8 M_☉ case. The highest temperatures are reached in the shearing region between the central core and the surrounding thick disk. The peak temperatures of around 4 × 10⁸ K are reached at densities of only about 10⁵ g cm⁻³, so that no efficient carbon burning is expected to occur.

Zoom In Zoom Out Reset image size

4.2. 0.6 + 0.9 M_☉: An Unstable Binary System with Direct Impact Mass Transfer

A 0.6 + 0.9 M_☉ binary system is predicted to be in the unstable, direct impact regime (P5 in Figure 1). This is seen in the three-dimensional density renderings plotted in Figure 15, which show the binary system at three different stages of its merger evolution. Although clearly unstable, the mass transfer continues for about 29 orbital periods before disruption. In contrast to the 0.2 + 0.8 M_☉ case, the orbit systematically shrinks from the onset of the mass transfer phase (Figure 16). The depletion of the donor star in the 0.6 + 0.9 M_☉ and 0.2 + 0.8 M_☉ binary systems is shown in Figure 17. In the 0.2 + 0.8M_☉ case, the donor mass decreases very regularly to about 90%, followed by a phase where the mass loss speeds up dramatically and the donor becomes subsequently disrupted during the remaining 20 orbits. Once mass transfer begins in the 0.6 + 0.9 M_☉ case, the rate of mass loss from the donor star increases exponentially with time, leading to a swift depletion of the donor.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The motion of the SPH particles within the Roche-type potential Ψ (see Equation (10)) are shown in Figure 18. As expected from a direct impact configuration, the angular momentum carried across is efficiently converted into the spin of the accreting star (Figure 19), with a small fraction being stored in the accretion torus. Similar to the 0.2 + 0.8 M_☉ binary system, the donor star is observed to remain close to synchronization during the entire duration of the mass transfer phase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The thermodynamical properties of the remnant core are displayed in Figure 20. The peak temperatures now reach about 5 × 10⁸ K at densities of about 10⁶ g cm⁻³. Under these conditions, carbon burning is observed to occur, but the rate of energy injection is not large enough to be dynamically dominant. The structure and further evolution of our DD merger remnant products will be discussed in more detail elsewhere.

Zoom In Zoom Out Reset image size

4.3. Summary of All Other Simulations

4.4. Comparison with Other Mass Transfer Calculations

The Laser Interferometer Space Antenna (LISA) will provide the largest observational sample of (interacting) double WD binaries (Farmer & Phinney 2003; Ruiter et al. 2010). Such systems enter the LISA observational window (0.1–100 mHz) when they reach a period ≈5 hr. They subsequently secularly evolve through radiation reaction across the LISA band and are so numerous as to create a stochastic foreground. When a binary reaches a gravitational wave frequency of a few mHz and is close enough to be resolved individually, the recorded signal shows an intrinsic frequency evolution over a typical one year observation time. Figure 21 shows the gravitational wave amplitudes (for an assumed distance of 10 kpc) for all our production runs together with the LISA sensitivity curve. Since none of the orbital frequencies change by more than 1% in the last year prior to merger, we chose this as a typical integration time. Provided they are in the Milky Way, all of the investigated systems would be detectable by LISA.

Zoom In Zoom Out Reset image size

We show that, depending primarily on the ICs, the frequency evolution of the binary can be calculated incorrectly. Simulations that use the approximate ICs, as discussed in Section 3.2.5, severely underestimate the initial separation when mass transfer sets in. As a result, the binary only survives for a few orbits, resulting in a rapidly fading gravitational wave signal (see, e.g., Figure 7). The accurate ICs, on the other hand, produce binary systems in which the mass transfer phase is dramatically extended by almost two orders of magnitude in time, resulting in a gravitational wave signal with amplitude and frequency that remain essentially constant up until merger. Population synthesis models have shown that these binaries are more numerous than dynamically stable systems for a given mass ratio (Nelemans et al. 2001b; Ruiter et al. 2010). As we show that these binaries can survive at small separation for hundreds of orbital periods (keep in mind our mass transfer duration estimates are—due to finite numerical resolution—strict lower limits on the durations realized in nature), their associated gravitational wave signal should be included when calculating the LISA gravitational wave foreground (although expected to below LISA's sensitivity at these high frequencies).

With the use of accurate ICs as well as realistic shock heating treatment, the frequency evolution of the binary can thus be explored in detail. Although there is a slight decrease in the orbital separation over the evolution of the 0.6 + 0.9 M_☉binary system, the gravitational wave amplitude remains close to the point mass prediction up to almost complete disruption, as seen in Figure 23. For the marginally unstable 0.2 + 0.8 M_☉ system, the increase in the orbital separation is more clearly imprinted on the gravitational wave signal, as seen in Figure 22.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

During the final phase of mass transfer prior to merger, the rate of mass exchanged approaches the mass of the donor divided by the orbital timescale,

$$\begin{equation} \dot{M}_2 \approx \left({M_2 \over R_2}\right)^{3/2} \sqrt{G \over \mu _1} \approx 0.1\, M_{\odot }\, {\rm s}^{-1}, \end{equation} \tag{ 16 }$$

where μ₁ is the reduced mass of the accretor. As shown in Guillochon et al. (2010), Kelvin–Helmholtz instabilities develop within the accretion stream once the accretor has built up a torus of material acquired from the donor, and this leads to the formation of dense knots of material that can periodically strike the accretor's surface. This striking action can lead to the initiation of a detonation within the torus about the accretor, assuming that matter acquired from the donor mostly consists of helium.

The question of whether knots are capable of detonating the accretor's helium layer depends on the ram pressure applied by individual knots, and by the ICs of the helium torus at the time of compression. From the SPH simulations presented in this paper, we have an accurate record of the time evolution of the accretion stream's density at L1 and the conditions at the base of the helium torus. As shown in our FLASH simulations (Figure 24), a standing shock develops between the stream and the torus, equalizing the pressure at the stream–torus interface, and reducing the velocity of the fluid flow in the torus relative to the accretion stream to subsonic velocities. As the interaction is subsonic and no heat transport processes can operate effectively on a dynamical timescale, the accretion stream smoothly evolves relative to the background pressure, and thus maintains constant entropy. The exterior pressure experienced by the accretion stream is given by the pressure within the torus.

Zoom In Zoom Out Reset image size

Given the initial density ρ_L1 and temperature T_L1 of the stream material as it leaves L1, the pressure at the base of the helium torus P_tor, and an EOS that parameterizes the entropy s and pressure P as functions of density and temperature, we can determine the density of the stream when it reaches the base of the helium torus by solving the following system of equations for ρ_knot and T_knot,

$$\begin{equation} s(\rho _{\rm L1}, T_{\rm L1}) = s(\rho _{\rm knot}, T_{\rm knot}) \end{equation} \tag{ 17 }$$

$$\begin{equation} P(\rho _{\rm tor}, T_{\rm tor}) = P(\rho _{\rm knot}, T_{\rm knot}), \end{equation} \tag{ 18 }$$

where "tor" refers to the base of the helium torus and "knot" refers to the stream post-compression. For all of the SPH simulations presented in this paper, the maximum increase in density achieved within the stream as it falls from L1 to the accretor's surface is ≳ 10. The small initial entropy of the stream (T_L1 ∼ 10⁵ K) and minimal shock heating means that the final temperature reached within the stream is substantially smaller than the temperature required for dynamical burning, and thus burning is not expected to be important within the stream itself.

With our estimate for ρ_knot, we can estimate the increase in density and temperature of the torus after it is compressed by the stream. As in Guillochon et al. (2010), we solve for the trajectory of a test particle released at L1 with initial velocity given by the sound speed of the donor to estimate v_⊥, the component of the velocity perpendicular to the accretor's surface. Given the density of the stream ρ_knot, we can determine the ram pressure of the stream,

$$\begin{equation} P_{\rm ram} = \rho _{\rm knot} v_\perp ^2. \end{equation} \tag{ 19 }$$

We can then use P_ram and solve a set equations similar to Equations (17) and (18) to determine the post-compression state of the torus,

$$\begin{equation} s(\rho _{\rm tor}, T_{\rm tor}) = s(\rho _{\rm comp}, T_{\rm comp}) \end{equation} \tag{ 20 }$$

$$\begin{equation} P_{\rm ram} = P(\rho _{\rm comp}, T_{\rm comp}), \end{equation} \tag{ 21 }$$

where "comp" refers to the conditions in the helium torus after compression. This enables us to calculate the triple-alpha burning timescale τ_3α and compare it to the dynamical timescale τ_dyn at the accretor's surface to determine if runaway burning can be achieved, leading in turn to a detonation.

In Figure 25, we show the curves for which τ_dyn = τ_3α and 10τ_dyn = τ_3α within the convex hull of the SPH simulations presented in this paper. For the two simulations which exhibited surface detonations in Guillochon et al. (2010) the two timescales are comparable, whereas t_dyn is far shorter than τ_3α for the simulation that showed no surface detonation. The SPH simulations show that the conditions for dynamical burning being triggered on the surface of the accretor are achieved for accretor masses larger than ∼0.8 M_☉, with a cutoff at low donor masses as the circularization radius of the accretion stream begins to exceed the radius of the accretor. Such systems could be possible candidates for sub-Chandrasekhar double detonation scenarios, as argued in Guillochon et al. (2010).

Zoom In Zoom Out Reset image size

The onset of (numerically resolvable) mass transfer in DD binary systems and its impact on the orbital evolution forms the main theme of our paper. We find for all investigated cases that mass transfer ensues for at least dozens of orbits, even for systems which according to the analysis of Marsh et al. (2004) are guaranteed to be unstable. Due to the finite numerical resolution, the duration of the mass transfer phase calculated here should be considered as a strict lower limit.

Part of our motivation has been to settle a previous debate concerning the duration of mass transfer in such systems. Almost all of the earlier SPH simulations had found quick disruption of the donor star within a few orbital periods after the onset of mass transfer, while recent grid-based simulations (D'souza et al. 2006) found a longer lived mass transfer phase. Using the same assumptions and ICs, we have reproduced essentially all previous results. In addition, we have compared our results to those obtained by using the StarCrash code developed by Faber and Rasio and found excellent agreement when using the same ICs and EOS.

Concerning the mass transfer duration debate, we find that the reported discrepancies can be accounted for by the following two effects. The first is the self-consistent inclusion of entropy production in shocks, which we find to have a major impact on the orbital dynamics. In direct-impact systems, we find the realistic inclusion of shock heating to produce an extended, high-entropy halo around the accretor. When ignored, a smooth incorporation of the transferred material into the stellar surface is observed. Since substantial amount of angular momentum can be stored in the shocked material (which is provided at the expense of the orbital angular momentum), the realistic inclusion of shocks yields shorter lived mass transfer phases. These results confirm the earlier concerns raised by Fryer & Diehl (2008).

The second and even more important reason is attributed to the initialization procedure of the binary in the simulations. The results of simulations constructed using approximate corotating ICs have been compared with those using carefully constructed binaries and are found to be markedly different. Simulations that use the approximate ICs underestimate the initial separation at which mass transfer sets in by about 15%, which yields a binary that only survives for a few orbits and thus a rapidly fading gravitational wave signal. This is because once mass transfer sets in, the orbital evolution is no longer driven by angular momentum losses due to gravitational wave emission but instead is dominated by the redistribution of angular momentum in the system. The use of approximate ICs systematically underestimates the number of orbits and overestimates the gravitational wave peak frequencies and amplitudes by about 20%. In contrast, the accurate ICs produce a binary system in which the mass transfer phase is extended by almost two orders of magnitude in time, resulting in a gravitational wave signal with amplitude and frequency that remain essentially constant up until merger.

The use of accurate ICs and a correct treatment of shock heating allows us to better characterize the time evolution of the temperature, density, and angular momentum, which are important when considering thermonuclear events that may occur during the mass transfer phase and/or after merger. Approximate ICs, for example, tend to significantly overestimate the densities and temperatures of the final merger remnant. This fact may have adverse consequences for possible Type Ia candidate systems. Our treatment allows us to also accurately identify when surface detonations may occur in the lead-up to the merger. The results of this study have been used as boundary conditions for FLASH simulations that focus on the fate of the low-density accreted material and that are described in a companion paper (Guillochon et al. 2010). This combined study has shown that contrary to earlier beliefs, the mass transfer can be very eventful. We have shown that for the case of unstable, direct impact He accretion onto the accretor a hot (>5 × 10⁸ K) helium torus builds up in which the incoming accretion stream triggers Kelvin–Helmholtz instabilities. Such Kelvin–Helmholtz "knots" repeatedly impact the accretor's helium surface and, once the triple-alpha timescale is shorter than the local dynamical timescale, can trigger violent surface explosions. This explosion launches shocks into the accretor's interior that can subsequently ignite the carbon–oxygen core. This could possibly be the long sought-after explosion mechanism on how to trigger supernova explosion in a DD scenario which, together with the promising rates, could make DD mergers a viable Type Ia progenitor system.

We thank Lars Bildsten, Marcus Brüggen, Philipp Podsiadlowski, John Fregeau, Dan Kasen, Gijs Nelemans, and Ken Shen for very useful discussions. We acknowledge support from DFG grant RO 3399 (M.D. and S.R.), the David and Lucille Packard Foundation (J.G. and E.R.-R.), NSF grant AST-0847563 (E.R.-R.), and the NASA Earth and Space Science Fellowship (J.G.).