A PARAMETER-SPACE STUDY OF CARBON–OXYGEN WHITE DWARF MERGERS

With SPH, one uses particles as a set of interpolation points to determine continuum values of the fluid and model its dynamics. SPH is a Lagrangian method, meaning movement is automatically tracked, and regions of high density contain more particles and therefore are automatically more resolved. Moreover, SPH inherently conserves angular momentum in three dimensions, which is difficult to reproduce in grid codes except under specific coordinate systems and symmetries. SPH therefore allows one to efficiently simulate complex phenomena with a large range of length scales. It has become the method of choice for merger simulations, and so we chose it as well.

For our simulations, we use Gasoline (Wadsley et al. 2004), a modular tree-based SPH code that was designed and has been used for a wide range of astrophysical scenarios, from galaxy interactions to planet formation. It aims for tight controls on force accuracy and integration errors. Gasoline implements the Hernquist & Katz (1989) kernel—we use 100 neighbors—and uses the asymmetric energy formulation (Wadsley et al., Equation (8)) to evolve particle internal energy. In our simulations, total energy is on average conserved to 0.3%, and angular momentum to 0.006%.

By default, Gasoline uses the usual Monaghan and Gingold formulation for artificial viscosity (see Monaghan 2005), together with a Balsara switch (a standard feature of WD merger SPH simulations) to reduce viscosity in non-shocking, shearing flows. Guerrero et al. (2004) found that such a prescription did not reduce viscosity sufficiently, resulting in excess spin-up of the remnant core and associated shear heating. Yoon et al. (2007), in addition to a Balsara switch, used variable coefficients for the linear and quadratic viscosity terms in the SPH equations of motion and energy, setting these values to α = 0.05 and β = 0.1, respectively, where shocks are absent, and around unity where they are present. A similar formulation was used in Dan et al. (2011, 2012). Since Gasoline includes it as well, we have used it for our study. Excess viscosity nevertheless remains a potential problem; we investigate its effects further in Section 4.5.

We modified Gasoline to include support for degenerate gas through the Helmholtz equation of state (EOS⁴; Timmes & Swesty 2000). This code, also used in Raskin et al. (2012) and Dan et al. (2012) simulations, interpolates the Helmholtz free energy of the electron–positron plasma, along with analytical expressions for ions and photons, to determine pressure, energy, and other properties from density and temperature. It is fast, spans a large range of density and temperature, and has, by construction, perfect thermodynamic consistency. To obtain quantities as a function of density and internal energy, we utilized a Newton–Raphson inverter. To keep the energy–temperature relation positive-definite, we did not disable Coulomb corrections in cases where specific entropy became negative.

Gasoline keeps track of the internal energy of particles, using it to determine other thermodynamic properties for fluid evolution. A particle's energy will naturally fluctuate due to noise, but for nearly zero-temperature particles this could result in their energy dipping below the Fermi energy. In such situations we keep the pressure at the Fermi pressure, while letting the energy freely evolve. A consequence of the floor is that a small amount of excess energy is injected into the system through mechanical work, which eventually manifests as additional thermal energy. The accumulated energy over a simulation is typically a small fraction of the internal energy, and therefore does not significantly affect the dynamics of the merger or most properties of the remnant. In cold, degeneracy-dominated material, however, a small change in internal energy corresponds to a large temperature change, at times comparable to the physically expected values, and thus the temperatures near the centers of some of our simulations have been affected. We characterize this spurious heating in Section 4.6 and show that it does not unduly affect our work's conclusions. However, it makes it difficult to run much longer simulations.

We also place an energy floor at half the Fermi energy. This is to prevent particle energies from approaching zero (and consequently calling for tiny time steps), which under rare circumstances occurs when particles perform a great deal of mechanical work. We find this happens primarily for particles that are flung out of the system by the merger and are cooling rapidly, and therefore are confident it has only a very minor effect on our simulations.

In our work, we ignore outer hydrogen and helium layers, composition gradients, and any nuclear reactions. This is mainly because previous work has found that nuclear processing was unimportant during the merger. For instance, LIG09 found fusion released ∼10⁴¹ erg for their 0.6–0.8 M_☉ merger, orders of magnitude smaller than the ∼10⁵⁰ erg binding energy of the remnant. Only for mergers involving very massive, ≳ 0.9 M_☉ WDs might this assumption break down, with the possibility of carbon detonations arising (Pakmor et al. 2010, 2011, 2012; but see Raskin et al. 2012; Dan et al. 2012). Similarly, Raskin et al. (2012), who included standard helium envelopes of ∼1%–2% of the WD mass in their simulations, found that only for accretors with masses above ∼1 M_☉ did it make a substantial difference: a helium detonation would inject ∼10⁴⁹ erg into the merger remnant. While this led to additional heating, it was insufficient to trigger much carbon burning or unbind any portion of the remnant (helium detonations have also been found for lower-mass accretors with CO-He hybrid donors; Dan et al. 2012).

We created spherical WDs using pre-relaxed cells of particles rescaled to follow the appropriate enclosed mass–radius relation determined using the Helmholtz EOS. We assumed a composition of 50% carbon and 50% oxygen by mass, and a uniform temperature of 5 × 10⁶ K. The stars were then relaxed in Gasoline for 81 s (∼10–40 dynamical times, depending on the WD mass) with thermal energy and motion damped (to 5 × 10⁶ K and 0 cm s⁻¹, respectively) during the first 41 s, and left free during the remaining 40 s. Particle energy noise prevented cooling of ≳ 5 × 10⁶ g cm⁻³ material to below 10⁷ K. We checked that the density profile of each star after relaxation was consistent with the solution from hydrostatic equilibrium, and found this was the case—central densities, for example, agreed to within 2%. The radii of the relaxed stars, as defined by the outermost particle of a relaxed WD, on the other hand, were on average about 7% too small, reflecting our inability to model the tenuous WD outer layers.⁵

We used a constant particle mass of 10²⁸ g, so that a 0.4 M_☉ WD has 8 × 10⁴ particles, and a 1.0 M_☉ WD has 2 × 10⁵. These numbers are similar to those used by LIG09 and Yoon et al. (2007), and exceed the ∼2 × 10⁴ particles per star used by Dan et al. (2012). Raskin et al. (2012) performed a resolution test for a merger of two 0.81 M_☉ WDs, varying the number of particles per star from 10⁵ to 2 × 10⁶. They found differences of ∼2% in the mass of the core plus envelope, disk half-mass radius, and inner disk rotation frequency. The one qualitative difference they found was that at their highest particle resolution, the WDs failed to break symmetry and disrupt (note that they assumed co-rotating WDs, making such a stable contact configuration possible). We perform our own test in Section 4.7 and find similar results.

We relaxed 0.4, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.9, and 1.0 M_☉ WDs, and combined them in all possible permutations to form our parameter space of binaries. These values were chosen to represent the range of possible CO WD masses, with greater resolution near the empirical peak at ∼0.65 M_☉ of the mass distribution of (single) CO WDs (Tremblay & Bergeron 2009). We also performed additional simulations with 0.575–0.65, 0.625–0.65, and 0.64–0.65 M_☉ binaries to explore the outcomes of similar-mass mergers. We thus simulated 48 mergers in total.

We placed two relaxed, irrotational WDs in a circular orbit. We chose the initial separation a₀ such that the donor WD just fills its Roche lobe, taking the location of the donor's outermost particle as its radius and using the Roche lobe approximation (for a synchronized binary) from Eggleton (1983).

This simple initial condition is similar to that of Pakmor et al. (2010), and implies that the binary system as a whole is not equilibrated. Therefore, as the simulation begins, the two WDs react to the tides, become stretched, and strong Roche lobe overflow ensues because the donor overshoots its Roche radius (in a widely separated binary, the donor would start to pulsate). As a result, the donor disrupts after just one to two orbits. For synchronized binaries, Dan et al. (2011) showed that the onset of mass transfer is much more gentle if the WDs are relaxed in the binary potential, disruption occurring only after several dozen orbital periods. They also showed that this results in systematic changes in the merger remnants. It is not clear whether the same will hold for unsynchronized binaries, since the accretion stream hits a surface that, in its frame, counterrotates, and therefore accretion is always much less gentle than for synchronized WDs. The difference is particularly dramatic for similar-mass binaries, where, in the synchronized case, the WDs can come into gentle contact, while in the unsynchronized case, any contact is violent. Unfortunately, it is difficult to test the effect of proper equilibration for unsynchronized binaries, since one has to relax to non-trivial initial conditions. A better approximation was attempted by LIG09 and Guerrero et al. (2004), who started their WDs farther out and reduced the separation artificially until mass transfer began. In their simulations, disruption still followed very quickly. Given that, and wanting to avoid any partial synchronization, we kept our simpler setup, and tested it by running simulations with varying a₀. We will discuss these tests in Section 4.2 and compare our results with those of others in Section 5.

It is difficult to decide when a merger is "complete," since for some cases remnant properties continue to evolve long after the two WDs coalesce, with (artificial) viscosity redistributing angular momentum and heating the remnant. As a visually inspired criterion, we decided initially to use the degree of non-axisymmetry, continuing simulations until they were less than 2.5% non-axisymmetric, as measured from the ratio of zeroth to largest non-zero Fourier coefficient of particles binned in azimuth. However, this had its own issues: in dissimilar-mass mergers—where most of the particles are in the accretor, already roughly axisymmetric following the merger—our convergence criterion was achieved while the outer disk was still obviously non-axisymmetric. In equal-mass mergers, which are inherently more axisymmetric, completion also was too soon, before the densest material had reached the center of the remnant.

For the majority of our systems, however, the time required to reach 2.5% non-axisymmetry was roughly constant in units of the initial orbital period, at 6.1 ± 1.2. For about the same time, axisymmetry was also achieved (by subjective visual inspection) for both dissimilar-mass mergers (except, in extreme dissimilar-mass cases, the outermost regions of their disks) and for equal-mass mergers (where the densest material had reached the center). We therefore use six orbital periods of the initial binary as the completion time of our simulations. In Section 4.4, we discuss the effect of continuing our simulations for two further orbital periods.

As found for previous simulations, qualitatively the most important factor controlling the merger outcome is whether the WD masses are "dissimilar" or " similar." In the former case, where the donor is significantly less massive than the accretor, only the donor overflows its Roche lobe,⁷ is disrupted, and accretes onto the accretor. The accreted material is heated on impact, lifting degeneracy. Hence, the merger remnant consists of a partly non-degenerate hot envelope and small, thick sub-Keplerian disk, both surrounding a cold core containing the largely unaffected accretor.

In the latter case of a similar-mass merger, there is a large degree of mixing between the two stars. For exactly equal masses, both stars are disrupted simultaneously, and their accretion streams impact each other near the system's barycenter. Material from the centers of both stars initially forms a thick, cold, dense torus orbiting the barycenter; this torus slowly shrinks due to viscous drag, pushing the accretion stream material above and below the equatorial plane. When the stars have slightly different mass, the lower-mass one disrupts first, forming an accretion stream (or series of streams) that mixes with accretor material down to the center of the accretor (regardless of whether or not the other also disrupts).

We show the differences between similar and dissimilar-mass merger remnants using two representative examples in Figures 1(a) and (b): a 0.4–0.8 M_☉ highly dissimilar and a 0.6–0.6 M_☉ equal-mass merger, respectively. One sees that the remnant morphologies are very different, consistent with previous work. The 0.4–0.8 M_☉ merger features a cold, nearly non-rotating and thus spherically symmetric remnant core, surrounded by a hot envelope with roughly equal degeneracy and thermal support, which itself is surrounded on the equatorial plane by a rotationally supported non-degenerate thick disk that holds most of the angular momentum. The accretor forms the core, largely undisturbed by the merger, while the envelope and disk are composed almost entirely out of donor material. The hottest points are on the interface between the core and the envelope.⁸ The 0.6–0.6 M_☉ remnant, on the other hand, has a massive, hot, partly rotationally supported and thus ellipsoidal core, and a very small but thick disk, both of which consist of material from both stars. No distinct envelope is formed. The hottest points are within the remnant core, just above and below the equatorial plane, arising from accretion stream material pushed out by the shrinking dense torus.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (48 panels)

Tarball of figure set images

A good way to visualize how mergers transition between dissimilar and similar-mass is to look at changes in the remnant properties with varying donor mass. In Figure 2, we show curves for accretors of 0.65 (left) and 1.0 M_☉ (right). One sees that remnants of highly dissimilar-mass mergers, with mass ratio q_m ≡ M_d/M_a ≲ 0.5, have properties resembling the 0.4–0.8 M_☉ merger: their donor and accretor barely mixed, their temperature curves have off-center hot plateaus, and their angular velocity profiles feature an off-center bump. The equal-mass, q_m = 1 cases resemble the 0.6–0.6 M_☉ remnant: they have flat temperature profiles and centrally peaked angular velocity profiles. Intermediate cases have intermediate profiles, with the bumps in the temperature and angular velocity profiles widening with increasing q_m. The 0.4–0.8 M_☉ and 0.6–0.6 M_☉ remnants therefore lie at the extremes of what merger remnants look like.

Zoom In Zoom Out Reset image size

The similarity between some of the curves for the 0.65 and 1.0 M_☉ accretors in Figure 2 suggests a homology. The similarity is closest for mergers with the same mass difference ΔM, as can be seen in Figure 3. For equal-mass mergers, all profiles are similar, simply scaled by a factor that depends on the total mass (except the 1.0–1.0 M_☉ merger; see below). As ΔM increases, the profiles are slightly less similar: with increasing total binary mass, the degree of mixing decreases, and the temperature and angular velocity maxima drift to slightly lower fractional enclosed mass. Nevertheless, the profiles still resemble one another far more closely than they resemble curves with other ΔM. The same holds for profiles along the rotational axis.

Zoom In Zoom Out Reset image size

It may seem surprising that the controlling parameter between these approximate homologies is the mass difference ΔM rather than the mass ratio q_m. Empirically, however, the case is clear: e.g., the 0.4–0.5 (second column, yellow) and 0.8–1.0 M_☉ (third column, black) mergers have the same q_m, but different ΔM, and their structures clearly differ from one another. The same is true for the 0.4–0.6 (third column, cyan) and 0.6–0.9 (fourth column, brown) M_☉ mergers. As we discuss below, the similarity of mergers of similar ΔM likely reflects the close relation between the ratio of central densities and mass difference.

Before discussing the homologies and trends further, we should note the one dramatic exception. The 1.0–1.0 M_☉ simulation differs fundamentally from its fellow ΔM = 0 M_☉ mergers. During the evolution of this system, unlike for all other equal-mass mergers, one WD was fully disrupted before the other, and as a result material from one star (arbitrarily designated the "donor" before the start of simulation, hence the "inverted" mixing profile in Figures 2 and 3) preferentially resides near the center of the remnant. This system also often appears as an outlier in Section 3.2. The 0.9–0.9 M_☉ merger also did not have equal mixing between the two stars, though the difference is much smaller. Raskin et al. (2012) noticed the same effect in their simulations, and concluded it reflected the fact that more massive WDs are much more concentrated and therefore harder to disrupt. This seems a likely explanation.

A major goal of our work is to establish how various global properties of the merger remnant, such as remnant core and disk mass, maximum temperature, and maximum angular velocity, vary as a function of accretor and donor mass. By quantifying these trends, we hope to help develop a parameterized model of merger remnants. Before discussing trends, however, we stress that they are necessarily approximate—second-order effects, numerical noise, and our choice of stopping time all affect the remnant properties. Moreover, while integrated values like total thermal energy do not fluctuate from time step to time step, values at specific points in the remnant do (as noted the following sections). For instance, the mass enclosed within the radius of peak equatorial temperature becomes ill-defined for similar-mass mergers because these have rather flat temperature profiles (Figure 2). To partly mitigate these fluctuations, the values presented below were determined by averaging frames from the simulation over an 8 s span, centered on the time corresponding to six orbits of the initial binary.

As might be expected from the approximate homologies described above, we found that many properties scaled well with ΔM. Of course, a scaling with a dimensional mass difference makes little sense; we believe its success reflects the fact that over the range of 0.4–1.0 M_☉, the central density ρ_c depends approximately exponentially on mass, with ρ_c ≃ 3.3 × 10⁷ g cm⁻³exp [5.64(M/M_☉ − 1)] (see Figure 4). Hence, a given mass difference ΔM corresponds to a given ratio of central densities, ρ_{c, d}/ρ_{c, a}. As argued in Section 3.2.3, ρ_{c, d}/ρ_{c, a} has a straightforward interpretation: it characterizes the degree of mixing between the donor and accretor. We therefore discuss trends as a function of q_ρ ≡ ρ_{c, d}/ρ_{c, a} from hereon. Where necessary, we refer to the mass ratio as q_m.

Zoom In Zoom Out Reset image size

3.2.1. What Constitutes Similar-mass?

As q_ρ increases from a small value toward unity, the merger remnant's morphology shifts from resembling Figure 1(a) (dissimilar-mass) to resembling Figure 1(b) (equal-mass). From Figure 2, ones sees that there is no particular q_ρ at which one transitions from "dissimilar" to "similar." Nevertheless, we can determine a rough critical value of q_ρ that separates mergers in which the core is largely unaffected from those in which it is changed significantly, a separation that likely affects the outcome of post-merger evolution.

To determine the critical value, we show in Figure 5 the ratio of central to maximum temperature, T_c/T_max, central to maximum angular velocity, Ω_c/Ω_max, and the fraction of donor to accretor material within the central core, (f_d/f_a)_cc, where we define the central core as a sphere with radius h_z. All three properties are measures of the extent to which the core has been affected: mixed regions tend to be hotter and more spun up, and contain material from both stars.

Zoom In Zoom Out Reset image size

From Figure 5, one sees that Ω_c/Ω_max approaches unity at q_ρ ≃ 0.6; at higher values, the angular velocity profile has a plateau or central peak rather than an off-center bump. Also at q_ρ ≃ 0.6, (f_d/f_a)_cc starts to deviate from zero, i.e., donor material begins to penetrate the central core. The temperature points show the transition is not abrupt: T_c/T_max starts to deviate from its downward trend (which reflects spurious heating in the most dissimilar-mass mergers; Section 4.6) around q_ρ ≃ 0.3 and continues to increase until q_ρ = 1.0; at q_ρ ≃ 0.6, T_c/T_max ≃ 0.5. Overall, this suggests that while the dependence is gradual, the morphology changes most around q_ρ ≃ 0.6. This conclusion is confirmed by looking at the two-dimensional remnant temperature structures (Figures 1(a) and (b)). At q_ρ ≪ 0.6, the remnant core has a large, spherically symmetric cold region, the nearly unperturbed accretor. This cold region shrinks with increasing q_ρ, and at q_ρ ≃ 0.6, spherical symmetry is broken. For still larger q_ρ, the cold region becomes a flat slice sandwiched between hot spots off the equatorial plane.

Given the above, we define "similar-mass" mergers as those with donor to accretor central density ratio q_ρ > 0.6, and "dissimilar-mass" mergers as those with q_ρ < 0.6. This critical density ratio corresponds to a mass difference ΔM ≃ 0.1 M_☉.

3.2.2. Structural Trends

Here and in the following subsections of Section 3.2, we describe various trends of remnant properties in detail, hoping to help attempts to interpolate between different simulations and motivate analytical and semi-analytical depictions of the merger. Readers not requiring this level of detail may wish to skip to Section 3.3. We begin our discussion of trends with size and density parameters.

The rotational axis central scale height. We define the rotational axis central scale height h_z as the characteristic width σ of a Gaussian fit to the density distribution along the z-axis at ϖ = 0. h_z is a measure of the vertical extent of the remnant. We find that the ratio h_z/h_a, where h_a is the central scale height of the accretor, is reasonably well-approximated by

$$\begin{equation} \frac{h_{\rm z}}{h_{\rm a}} = 1.03 - 0.17q_\rho ^{1/2} \quad (\pm 0.02), \end{equation} \tag{ 1 }$$

where the uncertainty listed in parentheses represents the root mean square (rms) of the residuals around the approximation (see Figure 6(a)). For highly dissimilar-mass mergers, h_z approximately equals the scale height of the accretor, while for similar-mass mergers, h_z is lower due to rotational support.

Zoom In Zoom Out Reset image size

The vertical scale height increases with increasing ϖ: the scale height at the location of maximum temperature, h(T_max), ranges from h_z to 1.21h_z, and the scale height at maximum angular velocity h(Ω_max) ranges from h_z to 1.88h_z. The prefactor for both heights increases with increasing accretor mass M_a and decreasing q_ρ.

The equatorial plane central scale height. Similar to h_z we define h_ϖ—the characteristic width of a Gaussian fit to the density distribution along the equatorial plane—as a measure of the equatorial extent of the remnant. The ratio h_ϖ/h_a can be parameterized by

$$\begin{equation} \frac{h_\varpi }{h_{\rm a}} = 0.96 + 0.89q_\rho ^2 \quad (\pm 0.08), \end{equation} \tag{ 2 }$$

where we excluded the 1.0–1.0 M_☉ merger remnant for our fit (see Figure 6(a)). The dependence on increasing q_ρ reflects the increased rotational support of the remnant core.

The central density of the remnant. The central density, ρ_c, is always within a factor two of the central density of the accretor, ρ_{c, a}. In Figure 6(b), one sees that for given accretor mass, ρ_c/ρ_{c, a} increases with increasing q_ρ for highly dissimilar-mass mergers due to increasing compression of the remnant core, but begins to decrease because of increasing rotational support around q_ρ ≃ 0.3. We could not find a simple parameterization for these curves. Note that for some systems, as we continue running our simulations ρ_c/ρ_{c, a} continues to increase. As discussed in Section 2.3, this is probably because artificial viscosity forces the merger remnant to undergo accelerated viscous evolution.

3.2.3. Mass Distributions

The merger mixes material between the donor and accretor. Here, we describe how this changes as a function of q_ρ, as well as how the material is distributed between the pressure-supported core and envelope and the rotationally supported disk.

The masses of the core-envelope and disk. We formally define the core-envelope as the part of the remnant inside the inner disk radius ϖ_disk, i.e., that is supported primarily by pressure (degeneracy for the core, thermal for the envelope) and not rotation. Since in every merger very little mass is ejected, a trend for either the core-envelope or the disk mass (M_ce and M_disk, respectively) suffices to determine both. The ratio of M_ce to the accretor mass M_a is well described by

$$\begin{equation} \frac{M_{\rm ce}}{M_{\rm a}} = 1 + 0.81q_\rho \quad (\pm 0.03), \end{equation} \tag{ 3 }$$

if the 1.0–1.0 M_☉ merger is neglected, and the fit's y-intercept is forced to unity. See Figure 7(a).

Zoom In Zoom Out Reset image size

The mass enclosing 50% of the donor material. The farther the donor penetrates, the smaller will be the mass enclosing half the donor's material, M_enc((1/2)M_d). For mergers with q_ρ ≲ 0.8, M_enc((1/2)M_d)/M_a ≃ 1.30, with an rms residual of 0.03 (Figure 7(c)). We present this trend mostly because we discuss similar thermodynamic and rotational enclosed masses, but it is somewhat difficult to interpret physically, since M_enc((1/2)M_d) increases with donor mass but decreases with mixing, which also depends on donor mass. The trend is easier to interpret using enclosed accretor mass rather enclosed total mass, as done below.

The accretor mass enclosing 50% of the donor material. As a different measure of the depth to which the donor penetrates, we consider just the accretor material within the mass enclosing half the donor, M_enc((1/2)M_d) − (1/2)M_d. This should equal the accretor mass if the donor is deposited above the accretor, and half the accretor mass if the two stars are completely mixed. For q_ρ ≲ 0.8, it can be approximated by (Figure 7(b))

$$\begin{equation} \frac{M_{\rm enc}\left(\frac{1}{2}M_{\rm d}\right)-\frac{1}{2}M_{\rm d}}{M_{\rm a}} = 1 - 0.190q_\rho \quad (\pm 0.009), \end{equation} \tag{ 4 }$$

where we forced the intercept to be unity. In this regime, roughly half of the donor remains outside of the accretor, though the trend discussed next indicates that the other half which does penetrate the accretor is spread across a much larger region at higher q_ρ. When q_ρ ≳ 0.8, the ratio drops sharply downward, indicative of the more thorough mixing expected for the similar-mass case. However, the existence and exact location of this drop may be a function of initial conditions (see Section 4.2).

The region over which the donor is spread. As a measure of the thickness of the region affected by the merger, we use the difference of the mass enclosing 75% of the donor material with that enclosing 25% of the donor material, i.e., ΔM_enc(M_d) = M_enc((3/4)M_d) − M_enc((1/4)M_d). Since 50% of the donor is within this range, ΔM_enc(M_d) − (1/2)M_d is a measure of the amount of accretor mixed with the donor. For q_ρ ≲ 0.8, the ratio of the latter to the total accretor mass follows

$$\begin{equation} \frac{\Delta M_{\rm enc}(M_{\rm d})-\frac{1}{2}M_{\rm d}}{M_{\rm a}} = 0.30q_\rho \quad (\pm 0.02), \end{equation} \tag{ 5 }$$

while for q_ρ ≳ 0.8 the trend curls upward until it reaches 0.5, the value expected for completely mixed remnants. See Figure 7(d).

Combining the two above trends, we can formulate a qualitative picture of mixing. For q_ρ ≲ 0.8, the donor can be thought of as being deposited onto the accretor and mixing with the accretor's outer layers, while for q_ρ ≳ 0.8, the accretor also disrupts substantially, leading to a regime where both stars mix more uniformly. The region over which the donor is spread, or thickness of the mixed layer, in both cases depends on q_ρ, which suggests that the relative densities of the donor and accretor govern mixing, i.e., the donor mixes significantly with all accretor material up to some fraction of the central density of the donor. Additional evidence of this will be seen in the thermodynamic trends below.

One might consider an alternate picture in which the donor dredges up a constant fraction of its own mass in accretor material. If this were the case, we would expect (ΔM_enc(M_d) − (1/2)M_d)/M_d to roughly be constant. Our results, however, show that for q_ρ ≲ 0.8 this quantity is nearly a straight line that is close to zero for small q_ρ ((ΔM_enc(M_d) − (1/2)M_d)/M_d = 0.35q_ρ ± 0.02). This seems more consistent with mixing being determined by density.

3.2.4. Energy Balance

The energy balance of the remnants indicate their primary means of support. Since the remnants are virialized, we consider how the ratio of degeneracy, thermal, and rotational energy to the total internal energy of the remnants varies with q_ρ.

Energy balance of the entire remnant. The support against gravity changes from being due mostly to degeneracy pressure at low q_ρ to having a substantial rotational contribution at q_ρ ≃ 1. This is because for highly dissimilar-mass mergers most of the internal energy is locked up within the accretor, which is hardly heated or spun up. For similar-mass mergers, however, donor material mixes, to some degree, with the entire accretor, causing heating and spin-up throughout the entire remnant.

The total gravitational potential energy of the merger remnant can be described adequately by a constant fraction of $GM_{\rm tot}^2/R_{\rm a}$,

$$\begin{equation} \frac{-E_{\rm pot}}{{\it GM}_{\rm tot}^2/R_{\rm a}} = 0.49 \quad (\pm 0.01). \end{equation} \tag{ 6 }$$

From the virial theorem, the internal energy should be related to the potential energy by 3(〈γ〉 − 1)E_I = −E_pot, where 〈γ〉 is an appropriately averaged equivalent to the adiabatic index. Since our remnants have cores where the electrons are becoming relativistic, one has 〈γ〉 somewhat smaller than 5/3, especially for the more massive remnants. We find that the ratio E_I/|E_pot| can be described by

$$\begin{equation} \frac{E_{\rm I}}{|E_{\rm pot}|} = 0.18\frac{M_{\rm a}}{M_\odot } + 0.42 \quad (\pm 0.01), \end{equation} \tag{ 7 }$$

which is ∼0.5 and ∼0.6 for low and high M_a, respectively.

The fraction of the internal energy carried by degeneracy and rotation is fairly well described by

$$\begin{eqnarray} \frac{E_{\rm rot}}{E_{\rm I}} &=& 0.31q_\rho ^{1/2} \quad (\pm 0.01), \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \frac{E_{\rm deg}}{E_{\rm I}} &=& 0.92 - 0.34q_\rho ^{1/2} \quad (\pm 0.02). \end{eqnarray} \tag{ 9 }$$

With these, the fraction carried by thermal energy can also be calculated; as shown in Figure 8(a), the fraction in thermal energy first increases with increasing q_ρ, but turns over at q_ρ ≃ 0.7, decreasing afterward. This reflects the competition between increased thermal energy from the two stars mixing, and increased rotational support from the spin-up of the core.

Zoom In Zoom Out Reset image size

Overall, for highly dissimilar-mass mergers, the internal energy is partitioned into degeneracy, rotational and thermal energy with a ratio of approximately 8:1:1, reflecting that, as stated above, such mergers are almost entirely supported by degeneracy pressure. Similar-mass mergers, on the other hand, partition their internal energies with the ratio 6:3:1, i.e., rotational support is significant.

Energy balance of the core-envelope. Since the variations with q_ρ seen for the remnant as a whole are almost entirely due to variations in the core and envelope rather than in the disk, the trends we find for the core-envelope are very similar to those we found above for the entire remnant,

$$\begin{eqnarray} \frac{E^{\rm ce}_{\rm rot}}{E^{\rm ce}_{\rm I}} &=& 0.28q_\rho \quad (\pm 0.02), \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \frac{E^{\rm ce}_{\rm deg}}{E^{\rm ce}_{\rm I}} &=& 0.94 - 0.32q_\rho \qquad (\pm 0.02). \end{eqnarray} \tag{ 11 }$$

Note the dependency on q_ρ, rather than on $q_\rho ^{1/2}$ as was found for the entire remnant. See Figure 8(b).

Energy balance of the disk. For the disk, we find very little dependence on q_ρ, consistent with the idea that most of the changes in the partitioning of energy have to do with increased mixing between the donor and accretor, which affects the core and envelope much more than the disk. Averaged over all mergers, we find

$$\begin{eqnarray} \frac{E^{\rm disk}_{\rm rot}}{E^{\rm disk}_{\rm I}} &=& 0.74 \quad (\pm 0.03), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \frac{E^{\rm disk}_{\rm th}}{E^{\rm disk}_{\rm I}} &=& 0.19 \quad (\pm 0.02), \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \frac{E^{\rm disk}_{\rm deg}}{E^{\rm disk}_{\rm I}} &=& 0.07 \quad (\pm 0.02). \end{eqnarray} \tag{ 14 }$$

Hence, the disk is composed of non-degenerate, primarily rotationally supported material. See Figure 8(c). (Note that we do not try to define a ratio of internal to potential energy of the disk or core-envelope, since the potential energy of either is not straightforward to determine.)

3.2.5. Temperature and Thermal Energy

Since heating of the remnant is achieved through shocks and viscous dissipation, the most heavily mixed regions should also be the hottest. We focus on equatorial thermodynamic values, but consider the rotational axis as well for similar-mass mergers.

The maximum temperature. We find that the maximum temperature on the equatorial plane, T_max, scales with the potential of the accretor (Figure 7(e)),

$$\begin{equation} \frac{kT_{\rm max}}{{\it GM}_{\rm a}m_{\rm p}/R_{\rm a}} = 0.20 \quad (\pm 0.03). \end{equation} \tag{ 15 }$$

This scaling is natural in the limit of highly dissimilar-mass merger—for each nucleon of the order of GM_am_p/R_a is liberated and converted into thermal energy. The temperature and thermal energy profiles in Figure 2 show that with increasing q_ρ, additional thermal energy is deposited into the remnant, but this energy is spread over a larger region, such that the maximum temperature remains roughly the same even as q_ρ approaches unity.

For a dissimilar-mass merger, the highest temperature along the rotational axis, $T^{\rm z}_{\rm max}$, is found at the tenuous outer edge of the hot envelope. It is slightly higher than the maximum temperature found in the equatorial plane. With increasing q_ρ, however, the difference increases noticeably due to the two off-center hot spots found along the rotational axis in similar-mass mergers. Fitting $T^{\rm z}_{\rm max}$, we find $kT^{\rm z}_{\rm max}/(GM_{\rm a}m_{\rm p}/R_{\rm a}) = 0.24\pm 0.03$, though this does not capture the upturn for similar masses well.

All remnants with q_ρ ≳ 0.8 have convectively unstable cores along the equatorial plane. Artificially mixing these cores to make them isentropic decreases their maximum temperatures by 10%–50% (not shown in Figure 7, but see the left panel of Figure 16). All remnants are stable against convection along the rotational axis.

The density at the point of maximum temperature. For dissimilar-mass mergers, the density at the hottest equatorial point, ρ(T_max), depends mostly on the donor (see Figure 7(f)). For q_ρ ≲ 0.5, we find

$$\begin{equation} \frac{\rho (T_{\rm max})}{\rho _{\rm c,d}} = 0.34 \quad (\pm 0.02). \end{equation} \tag{ 16 }$$

This proportionality again suggests that, at least for dissimilar-mass mergers, the donor mixes with the accretor up to a fraction of the central density of the donor, as alluded to earlier.

At q_ρ ≳ 0.6, the dependence becomes less obvious, with ρ(T_max) varying from ∼25%–90% of ρ_{c, d}. Since for these density ratios the donor material starts to penetrate the central core of the accretor—and the accretor starts to disrupt as well—the simple picture of the donor mixing up to a fraction of its own central density may be breaking down. Furthermore, part of the spread in density reflects that for high q_ρ the equatorial temperature profiles become nearly flat (see Figure 2, left column), thus increasing the sensitivity to noise in the determination of the location (but not the value) of maximum temperature. This also affects our results for the enclosed mass, M_enc(T_max), below.

Since for dissimilar-mass mergers $T^{\rm z}_{\rm max}$ is located near the tenuous outermost regions of the hot envelope, where particle noise is high, the density at the point of maximum rotational axis temperature, $\rho (T^{\rm z}_{\rm max})$ (plus symbols in Figure 7(f)), varies wildly between about ρ(T_max) and one order of magnitude below ρ(T_max). For similar-mass mergers, $\rho (T^{\rm z}_{\rm max})$ appears to be ∼20%–50% of ρ(T_max).

The mass enclosed within the radius of maximum temperature. For dissimilar-mass mergers, the radius of maximum temperature occurs at an enclosed mass of M_enc(T_max) ≃ M_a, while for mergers with q_ρ ≳ 0.8, maximum temperature occurs at the center and M_enc(T_max) ≃ 0 (see Figure 7(g)). For q_ρ ≲ 0.5, we find

$$\begin{equation} \frac{M_{\rm enc}(T_{\rm max})}{M_{\rm a}} = 1-0.28q_\rho \quad (\pm 0.01), \end{equation} \tag{ 17 }$$

where the fit's y-intercept is forced to unity. Note that the maximum temperature occurs near the bottom of the mixed zone, which is why M_enc(T_max) is substantially smaller than M_enc((1/2)M_d). The reasons it starts to deviate from a tight trend at q_ρ ≃ 0.6 are the same as those for ρ(T_max): the breakdown of the simple mixing picture and the difficulty in determining the location of peak temperature for a broader plateau.

Since M_enc(T_max) < M_a, it may be surprising that ρ(T_max) is not higher than ρ_{c, d}. This is because the additional thermal and rotational support against gravity reduces the density gradient that would be required if degeneracy pressure were the only source of support.

For dissimilar-mass mergers, $M_{\rm enc}(T^{\rm z}_{\rm max})$ is only slightly higher than M_enc(T_max), since the remnant core is nearly spherically symmetric, apart from the fact that the hot envelope is slightly more extended in the vertical direction. For similar-mass mergers, the difference increases, reflecting the development of the off-center hot spots, until $M_{\rm enc}(T^{\rm z}_{\rm max})/M_{\rm a}\simeq 0.25$ for equal-mass mergers.

The mass enclosing half the remnant thermal energy. As a more robust measure of where thermal energy is deposited during the merger, we consider the mass enclosing half the remnant thermal energy, M_enc((1/2)E_th) (see Figure 7(h)). We find this is very close to the mass of the accretor,

$$\begin{equation} \frac{M_{\rm enc}\left(\frac{1}{2}E_{\rm th}\right)}{M_{\rm a}}= 1.06 \quad (\pm 0.04), \end{equation} \tag{ 18 }$$

if the 1.0–1.0 M_☉ merger is neglected. One sees turnovers at the extremes, for q_ρ ≲ 0.2 and q_ρ ≳ 0.8. The former likely is because thermal energy is deposited into a narrow strip right on the surface of the accretor, while the latter is probably due to the disruption of the accretor.

While smaller than M_enc((1/2)M_d), M_enc((1/2)E_th) is always larger than M_enc(T_max). This reflects that high density degenerate material has lower specific heat, so that for the same energy per unit mass the temperature is higher (see Figure 2).

The width of the remnant thermal energy. The mass enclosed between the 25th and 75th percentiles of thermal energy, ΔM_enc(E_th) = M_enc((3/4)E_th) − M_enc((1/4)E_th), is a measure of the extent of the remnant that has been heated (see Figure 7(i)). Ignoring the 1.0–1.0 M_☉ merger, it can be fit by

$$\begin{equation} \frac{\Delta M_{\rm enc}(E_{\rm th})}{M_{\rm a}} = 0.11+0.94q_\rho \quad (\pm 0.03). \end{equation} \tag{ 19 }$$

Here, we did not force the y-intercept to go to zero, which is expected physically but gives a substantially poorer trend.

3.2.6. Angular Velocity and Rotational Energy

For a dissimilar-mass merger, the donor carries most of the angular momentum. As a result, the hot envelope contains more angular momentum and features higher angular velocities than the core, since the envelope is where most of the accreted donor material resides. Spin-up of the accretor is accomplished through shocks, PdV work and shearing forces. For a similar-mass merger, the two stars carry similar amounts of angular momentum and thoroughly mix. Conservation of angular momentum then implies that the entire remnant rotates rapidly.

The maximum angular velocity. On the equatorial plane, the highest angular velocity, Ω_max, scales linearly with the orbital angular velocity of the pre-merger binary, Ω_orb = 2π/P_orb (see Figure 7(j)),

$$\begin{equation} \frac{\Omega _{\rm max}}{\Omega _{\rm orb}} = 3.8 \quad (\pm 0.6). \end{equation} \tag{ 20 }$$

The ratio of core-envelope to total angular momentum. For more similar-mass mergers, more angular momentum is deposited in the accretor and ends up in the core and envelope (see Figure 7(k)). The ratio of core-envelope to total angular momentum, L_ce/L_tot, is approximately

$$\begin{equation} \frac{L_{\rm ce}}{L_{\rm tot}} = 0.70q_\rho \quad (\pm 0.03), \end{equation} \tag{ 21 }$$

where we fit only for q_ρ > 0.25 and ignore the 1.0–1.0 M_☉ merger. For q_ρ ≲ 0.25, the trend becomes shallower, resulting in a non-zero intercept. This suggests that even in cases where the donor has negligible mass some angular momentum is transferred to the accretor.

The mass enclosed inside the radius of maximum angular velocity. For dissimilar-mass mergers, both M_enc(Ω_max) and M_enc(T_max) are about equal to M_a, with M_enc(Ω_max) slightly larger than M_enc(T_max): for q_ρ ≲ 0.55, M_enc(Ω_max)/M_a = 1.05 ± 0.06. This is consistent with the idea that the hottest and most spun-up regions are those where the donor mixed most strongly with the accretor. For q_ρ ≃ 0.6, the off-center angular velocity peak is replaced by a plateau, and M_enc(Ω_max) becomes ill-defined; for even larger q_ρ, the highest velocities occur in the center, and M_enc(Ω_max) ≃ 0. See Figure 7(l).

The mass enclosing half the remnant rotational energy. Like for the thermal energy, for very dissimilar masses, the mass enclosing half the rotational energy, M_enc((1/2)E_rot), is similar to the accretor mass (see Figure 7(m)). For q_ρ ≲ 0.8, we find

$$\begin{equation} \frac{M_{\rm enc}\left(\frac{1}{2}E_{\rm rot}\right)}{M_{\rm a}} = 1.12+0.27q_\rho ^{1/2} \quad (\pm 0.01). \end{equation} \tag{ 22 }$$

Note that unlike M_enc(Ω_max), M_enc((1/2)E_rot) continues to increase with q_ρ (except for exactly equal-mass mergers), a consequence of particles with lower angular velocity but large lever arm that carry substantial rotational energy (see Figure 2). Near q_ρ ≃ 0.8, the trend breaks as both stars are significantly disrupted. However, even exactly equal-mass mergers have more of their rotational energy stored in the outskirts (otherwise one would have M_enc((1/2)E_rot)/M_a ≃ 1).

The width of the remnant rotational energy. We measure the extent to which the remnant is affected by spin-up through the difference between the masses enclosing 25% and 75% of the rotational energy, ΔM_enc(E_rot) = M_enc((3/4)E_rot) − M_enc((1/4)E_rot) (see Figure 7(n)). Ignoring the 1.0–1.0 M_☉ merger, it is well described by

$$\begin{equation} \frac{\Delta M_{\rm enc}(E_{\rm rot})}{M_{\rm a}} = 0.12+0.77q_\rho \quad (\pm 0.03). \end{equation} \tag{ 23 }$$

Like for the thermal energy, one sees that for more similar-mass mergers, rotational energy is spread more widely throughout the remnant.

From our empirical results above, a qualitative picture of a merger emerges. A dissimilar-mass merger has the donor overflowing its Roche lobe and forming an accretion stream. This stream mixes with the accretor up to approximately the central density of the donor pre-merger, ρ_{c, d}. Those layers of the accretor that are denser than ρ_{c, d} are hardly affected, and form the cold core of the merger remnant, while the mixed material will form a partly thermally supported outer envelope, which somewhat compresses the core, as well as a rotationally supported disk.

At q_ρ ≳ 0.6, the above picture begins to break down, as portions of the donor start to penetrate to the center of the accretor. This results in substantial heating and spin-up of the central core: the merger becomes a similar-mass merger. As the masses become more similar, the distinction between donor and accretor is lost and both stars disrupt and form accretion streams. For all q_ρ ≳ 0.6, the remnants are similar: a large, ellipsoidal and partly rotationally supported hot core with two hot spots off the equatorial plane, surrounded by a small, hot disk.

For all mergers, the maximum temperature reached by dissipation of orbital energy is proportional to the accretor's gravitational potential energy. For increasing q_ρ, the maximum temperature remains similar, but the region over which the thermal energy is deposited widens. The density at maximum temperature is of the same order of magnitude as the central density of the donor, consistent with the mixing picture discussed above. The latter no longer holds for q_ρ ≳ 0.6, when the entire remnant is mixed and heated.

For a dissimilar-mass merger, the angular momentum remains in the outer regions, since most of it was originally carried by the donor. Angular momentum can be transferred between regions through shocks, PdV work or shearing forces, all of which becomes increasingly important as the donor penetrates deeper. As a result, with increasing q_ρ, the remnant core is spun up further. Where both WDs disrupt, leading to colliding accretion streams, even the densest regions of the remnant have high rotational velocities.

Ignoring fusion, WD mergers should be insensitive to changes in composition, since the dominant electron degeneracy pressure only depends on the mean molecular weight per electron, which is close to μ_e ≃ 2 for all likely compositions. To confirm this, we ran simulations assuming pure helium and pure magnesium for an equal-mass case (0.4–0.4 M_☉) and an dissimilar one (0.4–0.8 M_☉). The results are shown in Figure 9. One sees that most quantities indeed have very similar profiles.

Zoom In Zoom Out Reset image size

The set of profiles showing most variation are those of the temperature. These have similar shape, but different normalization. Since the thermal energy curves are very similar, it is clear that this reflects differences in heat capacity, which does depend on composition: for He composition, there are more (non-degenerate) ions than for our standard CO mixture, boosting the heat capacity (and thus lowering the temperature for given thermal energy), while for Mg composition, there are fewer, lowering the heat capacity (and increasing the temperature). As a result, the maximum equatorial plane temperatures for the 0.4–0.4 M_☉ simulations are 0.95, 1.48 and 1.68 × 10⁸ K for He, CO, and Mg, respectively, while for the 0.4–0.8 M_☉ simulations, they are 2.92, 4.05, and 4.47 × 10⁸ K.

The smaller differences seen for the other profiles reflect small differences in initial conditions. All WDs are constructed assuming T = 5 × 10⁶ K throughout, which implies more thermal energy for higher heat capacity. As a result, the relaxed He and Mg WDs are slightly larger and smaller, respectively, than the CO WD. These slight differences in radius translate into differences in initial separation, which in turn cause small differences in the angular velocity and rotational energy curves.

For our simulations, we chose an initial orbital separation a₀ for which a co-rotating donor would fill its Roche lobe. Since our (non-rotating) WDs are equilibrated in isolation, once the simulation starts they immediately begin to adjust to the tides and hence disrupt quickly. Ideally, one would allow them to adjust to the binary potential and start mass transfer properly. For non-synchronous rotation, however, this is not straightforward (see Section 5.2). Nevertheless, we try to get a sense of the influence of this by running simulations for two cases—0.6–0.6 M_☉ and 0.6–0.8 M_☉—with a₀ increased and decreased by 10% (see Figure 10).

Zoom In Zoom Out Reset image size

Our default simulations were considered complete at six orbits of the initial binary. For runs where a₀ was changed, we used 2.5% non-axisymmetry, and a requirement for the density to be highest at the remnant's center, as the completion criteria.⁹ Not surprisingly, runs with larger a₀ needed longer to achieve these criteria: with a 10% increase, the 0.6–0.8 M_☉ merger required 861 s, or ∼15 orbits, to complete, while the one with a 10% decrease required 230 s, or ∼5.5 orbits. Similarly, the 0.6–0.6 M_☉ merger with a 10% increase in a₀ required ∼13 orbits (725 s) to complete, compared to the seven orbits of a 10% decrease (283 s). In both cases, the increase simply reflects that it takes longer for the donor to be disrupted fully if a₀ is increased.¹⁰ For instance, for the 0.6–0.8 M_☉ binary with increased separation, it took almost a dozen orbits before full disruption, while at the standard separation disruption occurred after just 1.5 orbits. This feature is of particular interest because for mergers of synchronously rotating WDs Dan et al. (2011, 2012) and Raskin et al. (2012) all note almost immediate disruption of the donor when approximate initial conditions are used, and much delayed disruption for more accurate initial conditions (up to ∼30 orbits; see Section 5.2).

We find that the density profiles of the merger remnants are remarkably insensitive to varying a₀ (ρ_c changing by ≲ 2% for the 0.6–0.8 M_☉ merger, and ∼20% for the 0.6–0.6 M_☉ merger), and show substantial systematic changes only in the outer regions. The latter can be understood from the mixing and rotational profiles, where one sees that with increasing a₀, donor material is mixed less deeply into the accretor, and rotational energy is shifted outward, causing the rotational frequency to peak at lower values and larger radii. This reflects the increase in angular momentum with increasing a₀, which creates a more rotationally supported remnant (for both our systems, a 10% increase (decrease) in a₀ results in a 5% increase (decrease) in angular momentum). In the 0.6–0.8 M_☉ merger, the decreased mixing causes the accretor to be spun up less, thus lowering the rotational energy of the core, and narrowing the thermal energy plateau. These effects are also seen in the similar-mass case, where the center of the remnant receives less rotational support and becomes denser with increasing separation, and the mixing becomes less uniform.

Qualitatively, with increasing a₀, the properties change in a way that is similar to the changes seen with decreasing q_ρ, i.e., similar to mergers with more dissimilar mass: reduced mixing, larger disks and less core rotational support, and shifts in the thermal and rotational energies toward larger radii. The converse is also true, decreasing a₀ has similar effects as increasing q_ρ, i.e., the mergers become similar to those with more equal masses. The changes are substantial at times: e.g., with a 10% increase in a₀ for the 0.6–0.6 M_☉ remnant, the maximum equatorial temperature is reduced by 40%, while the corresponding density increases by 25% (for the rotational axis hot spots, the values are a 13% and 30% reduction, respectively), and the mass of the disk increases by 65%. Similar, though far less extreme, changes are seen for the properties of the 0.6–0.8 M_☉ remnant. All this makes a₀ one of the parameters to which our mergers are most sensitive.

In our simulations, the WDs have zero spin, i.e., we assume that tidal dissipation is too weak to synchronize their rotation. Whether or not this is correct is currently unknown, but to see what the effect could be, we ran simulations assuming synchronized rotation for 0.6–0.6 M_☉ and 0.6–0.8 M_☉ binaries (see Figure 11). We used approximate initial conditions for these systems, identical to those of Section 2.2 except that the stars rotate at the orbital angular frequency.

Zoom In Zoom Out Reset image size

As in Section 4.2, our unsynchronized runs are from our parameter space study, and use six orbits as their completion criterion, while our synchronized runs use 2.5% non-axisymmetry, and a requirement for the density to be highest at the remnant's center. Completion occurred at 407 s (∼8.5 orbits) for the synchronized 0.6–0.6 M_☉ simulation, and 314 s (∼6.5) orbits for the synchronized 0.6–0.8 M_☉.

The asynchronous and synchronous mergers differ mostly in the amount of heating and spin-up. This has two causes. First, for the synchronized binary, the total amount of angular momentum is about 10% larger, and high angular momentum material has more difficulty penetrating the accretor, as is evident in the 0.6–0.8 M_☉ mixing profile (because of the larger amount of angular momentum, the mergers also take about 1.5 orbits longer to achieve 2.5% non-axisymmetry). Second, in a synchronized binary, the donor and accretor have much less differential rotation with respect to each other, leading to much less spin-up and heating.

Both effects are largest for the equal-mass case. In particular, in a synchronized, equal-mass binary contact can occur without any friction, while in an unsynchronized one it involves shocks at the full orbital velocity. In consequence, for the synchronized case, rotational support is weaker in the center and stronger in the outskirts, causing the central density of our 0.6–0.6 M_☉ remnant to increase by ∼70% and the disk mass to increase by a factor of two. Furthermore, while for the non-synchronized case, the maximum temperature along the equatorial plane was found in the center, for the synchronized case it is found in the outskirts, and is more than a factor two lower (1.3 × 10⁸ K instead of 2.9 × 10⁸ K). The hot spots on the rotational axis also have much reduced temperature, 2.3 × 10⁸ K instead of 3.6 × 10⁸ K.

For the dissimilar-mass merger, the effects of synchronization are less dramatic: the accretor still spins up substantially, and rotational and thermal energy are deposited in roughly the same way. The main difference is that the synchronized case has slightly less mixing, causing a drop in total thermal energy and maximum temperature (from 4.1 × 10⁸ K to 3.5 × 10⁸ K on the equatorial plane, and from 4.6 × 10⁸ K to 4.3 × 10⁸ K on the rotational axis).

We considered our mergers completed after six orbits, since in that time they on average had reached our convergence criterion of 2.5% non-axisymmetry (see Section 2.3). To test the robustness of our results, we also determined properties attained after eight orbits. In Figure 12, we compare the six and eight orbit results for 0.6–0.6 M_☉ and 0.6–0.8 M_☉ binaries.

Zoom In Zoom Out Reset image size

We find our mergers show little evolution between six and eight orbits, with the largest changes seen for the angular velocity profiles. For the 0.6–0.8 M_☉ merger, the rigidly rotating core sped up and the off-center peak decreased in height and moved out, while for the 0.6–0.6 M_☉ merger, the center spun down and the rotational profile became flatter. In the dissimilar-mass merger, T_max and ρ(T_max) changed by ≲ 5% and the density and temperature structures nearly overlap, while in the equal-mass merger T_max and ρ(T_max) changed by ∼20% (2.9 to 2.3 × 10⁸ K and 1.7 to 2.0 × 10⁶ g cm⁻³), reflecting an increase in central density and a shifting of the temperature profile, with temperature decreasing in the center but increasing elsewhere. The evolution of all properties is consistent with viscous evolution—expected to follow the merger proper—with the core driven into rigid rotation, and angular momentum transferred outward to the disk. In the dissimilar-mass merger, the net effect is spin-up of the core and spin-down of the envelope, while in the equal-mass case it is the reverse. Of course, in the process, rotational energy is turned into thermal energy, heating the remnants.

One curious aspect for equal-mass mergers is the evolution of the off-center hot spots (Figure 1(b)). Over time, these broaden parallel to the equatorial plane, yet become narrower along the rotational axis. As a result, the hourglass shape is lost, and the center of the remnant stops being one of the hottest point in the system. After eight orbits, the system resembles more closely what we find for typical similar-mass mergers, which have more pancake-shaped hot spots flanking a colder, denser region on the equatorial plane.

Comparing more broadly the six- and eight-orbit results, we find changes at the ∼5% level. The trends presented in Section 3.2 continue to hold to within ∼10% except for h_ϖ (which becomes $h_\varpi /h_{\rm a} = 0.98 + 0.74q_\rho ^2 (\pm 0.07)$), the fraction of disk energy in degeneracy energy (E_{deg, disk}/E_{I, disk} = 0.07(± 0.02)), mass enclosed within the radius of maximum temperature (M_enc(T_max)/M_a = 1 − 0.21q_ρ(± 0.01)), maximum rotational frequency (Ω_max/Ω_orb = 3.4(± 0.5)) and the widths of the regions in which thermal and rotational energy are deposited (ΔM_enc(E_th)/M_a = 0.13 + 0.87q_ρ(± 0.03); ΔM_enc(E_rot)/M_a = 0.15 + 0.70q_ρ(± 0.02)). Changes to these trends are consistent with the viscous evolution described above: the remnant is beginning to spin down, lose its rotational support, and energy is being redistributed.

The addition of artificial viscosity is required in SPH to accurately capture shocks, but no consensus exists on how best to implement it. We ran two additional simulations of a 0.6–0.8 M_☉ merger to check the robustness of our results with respect to changes in the viscosity, one with small and one with large artificial viscosity (fixed (α, β) = (0.05, 0.1) and (1, 2), respectively). As before, these additional runs use 2.5% non-axisymmetry and a requirement for the density to be highest at the remnant's center as their completion criteria, and the low viscosity run completed at 198 s (4.1 orbits) while the high completed at 242 s (5.0 orbits). Here, we expect that low values of α will lead to large particle noise and inaccurate shock capturing, while high values result in large viscous heating and rapid loss of differential rotation. Our results confirm this (Figure 13): the simulation with low artificial viscosity leads to a remnant with stronger differential rotation, with the disk carrying 34% more rotational energy (and the remnant 33% less) than in the standard variable α simulation. Lower viscosity also leads to greater mixing of donor and accretor material, reflecting the stronger diffusion associated with the larger particle noise inherent to low viscosity.

Zoom In Zoom Out Reset image size

Aside from the mixing and spin-up, the results for the three different viscosity prescriptions do not differ greatly. While one might have expected greater dissipation of rotational into thermal energy for higher viscosity, the maximum temperatures and rotation rates vary by ≲ 10%, and the thermal and rotational profiles are quite similar. The density profiles are virtually identical except near the outer parts, where the low viscosity simulation leaves matter with greater rotational support.

As discussed in Section 2.1 and throughout Section 3, noise combined with a pressure floor in the EOS lead to small increases in internal energy. While this energy has a negligible effect on most remnant properties, in the most degenerate regions of the remnant it can cause significant temperature increases. Here, we discuss the extent to which spurious heating affects our results.

As a comparison for the spurious heating seen in some of the simulations, we relaxed a 0.8 M_☉ isolated WD for 489 s longer than the standard 81 s we used for relaxing single stars. While the total energy of the WD (potential, degeneracy, and thermal energy combined) increased by ∼1% of the original total energy over the additional period of time, the change in thermal energy was enough to raise the central WD temperature from 1.2 × 10⁷ K (increased from 5 × 10⁶ K due to particle noise) to 1.2 × 10⁸ K. In Figure 14, we compare the thermal energy profile of this isolated 0.8 M_☉ WD with those of a 0.4–0.8 M_☉ and a 0.7–0.8 M_☉ merger, with the times for the isolated WD taken at 489 and 224 s (the mergers' respective completion times) longer than the standard 81 s. The total thermal energy generated in the WD at an additional 489 s is ∼10% of the thermal energy generated in a 0.4–0.8 M_☉ merger.¹¹ Indeed, given how well the specific thermal energy profiles match in the interior, it is clear that spurious heating dominates there.

Zoom In Zoom Out Reset image size

Spurious heating is much less important for the 0.7–0.8 M_☉ merger, since its core has mixed to a much greater extent and less time was needed for the merger to complete. The total thermal energy generated in the isolated WD at 224 s is only ∼3% of the thermal energy generated during merger, and even at the very center spurious heating contributes ∼35% (rather than nearly all) of the thermal energy. As a result, the central temperature of the merger, 1.4 × 10⁸ K, is substantially higher than that of the isolated WD, 7.8 × 10⁷ K. It is important to note this comparison overestimates spurious heating, since the isolated WD will have had many more particles that dip below the Fermi energy than the much hotter merger remnant core.

Overall, we conclude that spurious heating is present, but is recognized fairly easily and does not influence our conclusions. In particular, its effects on remnants should be small in both high-density regions with T ≳ 3 × 10⁸ K and in lower-density ≲ 10⁶ g cm⁻³ regions. Other simulations may suffer from spurious heating as well. In this respect, it is intriguing that our equatorial temperature curves for a 0.6–0.8 M_☉ merger are a good match those of LIG09, even in the center (Figure 15; Section 5.1).

Zoom In Zoom Out Reset image size

To determine whether or not the numerical resolution matters for our results, we ran three additional simulations of a 0.625–0.65 M_☉ merger, with roughly a quarter, half, and double the number of particles (63,736, 127,525, and 510,047, respectively), corresponding to 0.63, 0.79, and 1.26 times the SPH smoothing length (resolution) we normally use. We also ran a second simulation with the same number of particles (255,035 particles). Here, we chose 0.625–0.65 M_☉ to see if numerical resolution has any effect on whether a merger is "similar-mass." All simulations were considered complete at six orbits of the initial binary, though we also checked the 2.5% non-axisymmetry convergence times.

From Figure 13, one sees that the two runs using the same number of particles—and identical initial conditions and the same version of Gasoline—still give slightly different results. This is due to the inherent nonlinear nature of fluid dynamics, coupled with small, random perturbations, e.g., from differences in the order of force addition in parallel processing, round-off errors, and slight inconsistencies in converting thermal energy to temperature. Overall, merger remnant properties change by ∼3.5% between the two runs. The most prominent differences are seen in properties determined from low numbers of particles, such as T_c (varies by ∼10%), and those involving finding maxima of temperature plateaus, such as $\rho (T^{\rm z}_{\rm max})$ (40%) and $\rho (T^{\rm cv}_{\rm max})$ (a factor of three—in one case convection shifts the temperature maximum off-center).

The differences for different resolutions are larger. While to first order the equatorial density and mixing profiles are very similar, there is a systematic ∼20% drop in the equatorial density—∼25% in the rotational axis density—near the center of the remnant with increasing numerical resolution (top right panel of Figure 13). The angular velocity and rotational energy profiles are again very similar, except in the central regions, where there is a ∼20% increase in Ω_max. For the temperature the effects are larger: with higher resolutions, most of the equatorial plane is colder, with a ∼20% drop in the value of the temperature plateau near M/M_tot = 0.5. The temperatures along the rotational axis, however, increase with increasing resolution, by ∼10% across the range of resolutions, as does the upturn in equatorial temperature near the center of the remnant, by ∼50% (∼30% if we do not include the lowest resolution run). The latter effects are due to increasing prominence of the off-center hot spots at higher resolutions, which also tend to look more hourglass-shaped. Indeed, for our lowest resolution, the densest material in the two stars remains relatively cold throughout the entire merger, resembling the synchronized systems described in Section 4.3. Finally, we find that if we do not consider the lowest resolution run, the disk half-mass radius varies by 4%, angular velocity at the half-mass radius varies by 4%, and the core-envelope mass changes by 3%. This is similar to the results of the resolution tests of Raskin et al. (2012).

The 2.5% non-axisymmetry convergence times for the half and double-particle number runs are within 14 s of the 275 s non-axisymmetry convergence time of the default run, a small difference that implies a negligible amount of post-merger evolution. Only the quarter-particle number run deviated substantially, converging 57 s earlier. This may simply reflect the smaller number of particles in the disk, where the system is most asymmetric.

We stress that even though the order-of-magnitude change in particle number (factor of two change in resolution) generates 10%–30% variations in some remnant properties, the overall shapes of the profiles in Figure 13 are very similar. In particular, the merger remnant does not look more or less "similar-mass" (except, arguably, the temperature curve at the lowest resolution). Exact values of properties, therefore, will vary depending on resolution (and will vary on similar or larger levels if initial conditions like a₀ are changed), but the overall picture of the merger and trends should be more robust.

LIG09 simulated a number of WD mergers, and gave detailed temperature, surface density, and rotational frequency curves for three. In Figure 15, we compare their results (from their Figures 3 and 4) with ours for two of these, 0.6–0.6 M_☉ and 0.6–0.8 M_☉ (the third was a 0.4–0.8 M_☉ He-CO WD merger, whose temperature profile cannot be compared directly). We note that they used different initial conditions, starting their systems with an orbital separation too large for mass transfer to begin, and then slowly reducing the separation until it does. This point defines their t = 0 and the start of the merger simulation proper. Given this different setup, their merger completion times cannot be compared directly to ours. In their simulations, however, coalescence (the final consolidation of the two WDs into one) also occurs after just about one orbit, so the differences should not be too large. To give a sense of the effect of different completion criteria, we compare their results both with our standard results, taken after six orbits, and our results taken at their merger completion times.

For both mergers, the surface density curves are similar, although in their 0.6–0.6 M_☉ merger, the central peak is ∼30% higher (∼10% if we use their completion time of 514 s). For the 0.6–0.8 M_☉ merger, the temperature profiles are also very similar, with maxima¹² differing by only ∼10%–15%, and having nearly identical shapes. Larger differences are seen for the 0.6–0.8 M_☉ rotational frequency profile, where the angular velocity peaks further out and at lower value (∼0.3 s⁻¹ compared to our 0.45 s⁻¹—or 0.50 s⁻¹ using their completion time of 164 s). Indeed, our entire remnant is more spun-up than theirs.

For the 0.6–0.6 M_☉ merger, LIG09 have a plateau in their angular frequency profile, with Ω_max ≃ 0.25 s⁻¹, while our profile is much more peaked and reaches a much higher frequency of 0.60 s⁻¹. By their completion time, our rotation curve is not as sharply peaked, but still reaches 0.44 s⁻¹. The temperature profiles are also much less similar: our maximum temperature in the equatorial plane is a factor of two lower than theirs (factor of 3.3 at their completion time), and even our maximum temperature along the rotational axis is a factor 1.6 lower (factor 1.9 at their completion time).

Finally, we can compare how mass is distributed. In both our simulations and those of LIG09, negligible mass is lost, so only the distribution between disk and core-envelope matters. For our 0.4–0.8 M_☉, 0.6–0.6 M_☉, and 0.6–0.8 M_☉ simulations, we infer disk masses of 0.31, 0.10, and 0.40 M_☉, respectively, which are reasonably close to the 0.28, 0.10, and 0.30 M_☉, respectively, listed by LIG09 (their Table 1), especially considering that we likely use a different definition of what is a "disk."

Overall, the primary differences between our simulations appear to be the amount of spin-up and heating of the equal-mass merger. We believe it is unlikely that this reflects differences in initial conditions: we found much smaller changes in the angular velocity profile with increasing a₀ (see Figure 10), and in the simulation of LIG09 the stars still seem to be quite close to spherically symmetric at the start and disrupt quickly (their Figure 1), even though they were more properly relaxed. Instead, we believe the more likely explanation is that the viscosity prescription of LIG09, based on Riemann solvers, yields larger effective viscosity. This would explain both the reduction in angular velocity and increase in temperature (since viscous evolution converts rotational into thermal energy), as well as the fact that similar-mass mergers are affected more (they mix more, and LIG09 ran their equal-mass merger for a very long time). If we ran our simulations longer and thus included further viscous evolution, the similarity with their simulations would likely be closer.

Simulations of WD mergers have also been presented by Yoon et al. (2007), Pakmor et al. (2010, 2011, 2012), Dan et al. (2011, 2012), and Raskin et al. (2012). Unfortunately, comparison with those results is difficult, since, unlike LIG09, all these authors are sparse with quantitative details about their results.

An exception is the 0.81–0.9 M_☉ merger simulated by Dan et al. (2012), shown in their Figure 1. While that simulation is for synchronized WDs, it is still particularly useful to compare with, since Dan et al. show results for both approximate and accurate initial conditions. We find that their spherically enclosed mass profile is very similar to ours, with, e.g., M = 0.9 M_☉ at 4.5 × 10⁸ cm in both (though since spherically enclosed mass is a cumulative quantity, significant structural differences can remain hidden). Our spherically averaged density profile looks most similar to the profile they found using approximate initial conditions. Our central density, 1.9 × 10⁷ g cm⁻³, is within ∼10% of theirs, and the density profile remains similar up to r ≃ 5 × 10⁸ cm (ρ ≃ 10⁶ g cm⁻³). Beyond, their profile becomes shallower while ours continues to decline; at r = 10⁹ cm, they find ρ ≃ 3 × 10⁵ g cm⁻³, while we find ρ ≃ 10⁵ g cm⁻³. This may be a consequence of the additional angular momentum associated with synchronized rotation. With accurate initial conditions, a difference with our results is that the density profile becomes flat beyond 10⁹ cm.

Comparing temperature profiles, we roughly reproduce their spherically averaged one for approximate initial conditions, including the off-center peak—their T_max is ∼20% lower (to be expected since their binary is synchronized; see Section 4.3), but is also located at 4 × 10⁸ cm (or M_r ≃ 0.9 M_☉). However, our central temperature (2.2 × 10⁸ K) is an order of magnitude higher than theirs (2 × 10⁷ K) and at r ≳ 10⁹ cm our temperatures are systematically hotter, perhaps a result of the much larger dissipation expected for non-rotating WDs. With accurate initial conditions they found an even narrower temperature peak than the one with approximate conditions, which thus deviates even more from our curve. This trend is similar to what we see when increasing a₀ (Section 4.2), so it seems likely we would reproduce their simulations more closely if we used the same initial conditions.

Many of the recent simulations (Dan et al. 2011, 2012; Raskin et al. 2012) assume co-rotating WDs. This assumption is numerically convenient, in that it is relatively straightforward to start the simulation in the physically correct state: since in the co-rotating frame there are no flow velocities, one can easily relax a simulated binary within an appropriate potential in the co-rotating frame, damping out any velocities resulting from an initial mismatch.

As a result, it has been possible to study the onset of mass transfer in detail. As first pointed out by D'souza et al. (2006) from simulations using a grid code, the disruption of the donor is preceded by a rather long—dozens of orbits—phase of mass transfer. Further simulations by Dan et al. (2011, 2012) showed that in this initial phase a significant fraction, ∼10% of the donor mass, is transferred. As a result, e.g., the disk is substantially colder and more extended. The remnant core seems more subtly affected, in that its appearance becomes "more dissimilar," reflecting that coalescence is between two WDs whose masses have become more disparate than they were initially. As a consequence, e.g., even for similar-mass binaries, the hottest point of the merger is found to be well outside the center. Indeed, Raskin et al. (2012) find that even for equal-mass binaries, the final outcome for more massive mergers is one where the core of one of the WDs is virtually undisturbed.

At present, it is not clear how important accurate initial conditions would be for asynchronous mergers. Qualitatively, we expect the effects to be smaller than for synchronous mergers for three reasons. First, from the analytic study of Lai et al. (1994), in which tidal and rotational distortion are approximated by ellipsoids, co-rotating binaries always reach contact or Roche lobe overflow before becoming dynamically unstable, while irrotational binaries become dynamically unstable first. While an exact treatment of the irrotational case found that, in fact, Roche contact preceded dynamical instability (Uryū & Eriguchi 1998), it suggests that WDs in irrotational binaries will disrupt much sooner. Second, the simulations of LIG09 use initial conditions that should be quite close to correct, yet their WDs disrupt quickly (see Section 5.1). Third, the two components are counterrotating in the rotating frame. Hence, any mass transferred will hit the accretor with a larger relative velocity than would be the case for co-rotating WDs. Indeed, in the limit of equal-mass WDs, very little would happen for co-rotating WDs when one reaches contact, while a strong shock would be expected for the irrotational case. In general, one expects part of the shocked material to enter a high-entropy halo around the accretor. For co-rotating WDs, Dan et al. (2011) found that this halo helps remove angular momentum from the orbit, leading to a shorter start-up phase. For the irrotational case, given the stronger expected shocks, the start-up phase would likely be reduced even further. On the other hand, we saw in Section 4.2 that remnant properties are sensitive to changes in angular momentum content through changes in a₀. Simulating more realistically the onset of mass transfer through accurate initial conditions will likely change a₀.

Ideally, one would still simulate the initial mass transfer phase accurately. Unfortunately, even though the equilibrium solution is known (Uryū & Eriguchi 1998), it is not straightforward to set up the initial binary properly, since it is difficult to relax to a state that includes substantial fluid motion, and to slowly evolve such a state to contact, while ensuring viscosity remains low enough that there is no artificial tidal dissipation. Such dissipation is seen in our tests with varying initial distance a₀ in Section 4.2 (and may affect the simulations of LIG09 as well). Prior to coalescence, strong dissipation of tidal bulges heats the outer envelope of the donor (both stars for similar-mass mergers), and spin-orbit coupling due to both tides and the direct-impact accretion stream result in both donor and accretor becoming 25%–50% synchronized by coalescence.

Since it significantly affects the merger and merger outcome, whether or not tidal dissipation causes real CO WD binaries to synchronize before the merger remains a major source of uncertainty. For the radiative stellar envelopes appropriate for WDs, tidal dissipation is expected to be inefficient, with a timescale 10¹² to 10¹⁵ yr, suggesting that WDs do not synchronize (Marsh et al. 2004 and references therein). However, coupling of the tides to pulsations may dramatically increase dissipation (Fuller & Lai 2012). Fortunately, it may be possible to determine the rate of synchronization observationally. For instance, Piro (2011) suggested that tidal dissipation is responsible for the relatively high temperature of the primary WD in the 13 minute eclipsing binary SDSS J065133.33+284423.3, predicting that it would be about halfway to being synchronized. This could be tested by either measuring the rotational broadening of the narrow cores of the hydrogen lines, or looking for velocity deviations through the transit of the more massive secondary.

Following the merger, processes that happen on timescales slower than the dynamical time can become important. These include viscous evolution, neutrino emission, radiative or convective thermal adjustment, and magnetic dipole radiation spin-down. Out of these, convection acts on the fastest timescale, and we already included its effect on the core in Figure 16 (left). Next fastest would almost certainly be viscous evolution. The merger remnant is unstable to both the magneto-rotational instability (Balbus & Hawley 1991) and Tayler–Spruit dynamo (Spruit 2002). Radiative adjustment is expected to be much slower, except at the surface, where radiative losses may also lead to convection in some systems (Shen et al. 2012; Schwab et al. 2012; Raskin et al. 2012). Using the standard Shakura & Sunyaev (1973) α-prescription for the viscosity ν = αc_sH, where c_s is the local sound speed and H is the scale height of the system, the viscous evolution timescale for the remnant disk is

$$\begin{equation} t_{\rm visc} = \frac{R_{\rm disk}^2}{\nu } = \frac{1}{\alpha }\left(\frac{R_{\rm disk}}{H}\right)^2t_{\rm dyn} \sim \frac{10}{\alpha }t_{\rm dyn}, \end{equation} \tag{ 24 }$$

implying a timescale t_visc ∼ 10³–10⁵ s for α ∼ 10⁻³ to 10⁻¹ and t_dyn ∼ 10 s. This is orders of magnitude smaller than both the neutrino loss timescale (≳ 10³ yr; see Figure 16 (left)) and thermal adjustment timescale (≳ 10⁴ yr; Shen et al. 2012).

It is possible that the strong differential rotation during a merger results in substantial amplification of magnetic fields. The one known probable WD merger remnant, RE J0317−853, has a surface magnetic field of 3.4 × 10⁸ G (Barstow et al. 1995; Külebi et al. 2010). If mergers lead to strongly magnetized WDs, and these WDs additionally drive an ionized outflow, the magnetic coupling between the outflow and the WD could serve to transport angular momentum out of the system, spinning down the WD. The timescale for such a spin-down is roughly given by

$$\begin{equation} t_{\rm msd} \sim \frac{L}{\dot{M}R_A^2\Omega } \sim \frac{L}{(\dot{M}\Omega)^{3/5}(B^2R^6)^{2/5}}, \end{equation} \tag{ 25 }$$

where L is the angular momentum of the remnant, $\dot{M}$ is the mass-loss rate, R_A is the Alfvén radius, and Ω is the angular spin frequency. For the second approximation, we used that $R_A\sim (B^2R^6/\dot{M}\Omega)^{1/5}$, with B the surface magnetic field and R the remnant radius. Scaling to B = B₈10⁸ G and $\dot{M}={\dot{M}}_{-7}10^{-7}\,M_{\odot }{\rm \,yr^{-1}}$ (similar to what is observed for RE J0317−853—see above—and [WR] cores of planetary nebulae; Hamann 1997), and using the properties inferred for a 0.6–0.6 M_☉ remnant (L ≃ L_tot ≃ 10^50.5 g cm² s⁻¹, R ≃ R_disk ≃ 10⁹ cm, Ω ≃ Ω_max ≃ 10^−0.3 s⁻¹, we find $t_{\rm msd}\simeq 8\times 10^3\,B_8^{-4/5}{\dot{M}}_{-7}^{-3/5}{\rm \,yr}$, which is of the same order as the neutrino cooling timescale of ∼10⁴ yr at the ignition line (for the whole range of remnants, 2 × 10³ ≲ t_msd ≲ 5 × 10⁴ yr).

Accretion from the disk, loss of rotational support, and possible cooling of the hot envelope could all compress and heat the remnant core. A detailed study of this is beyond the scope of this paper, but we can make first-order estimates of the effects on our merger remnants, and compare these with the more detailed analysis of Schwab et al. (2012) in one specific case.

For our estimates, we make four assumptions: (1) spin-down and accretion are much faster than thermal processes, and do not lead to local dissipation (i.e., particles entropies are constant in time); (2) all angular momentum is carried away to large distances; and (3) corresponding matter ends up with zero total energy (i.e., is at large distances and has negligible kinetic and internal energy). From energy conservation, the last assumption implies that the remaining object will have the same total energy as our merger remnant (but a lower mass), the first that it will have the same entropy structure, and the second that it has no rotational support. To determine the properties, we first determine the entropy profile of the merger remnant by averaging entropy over isopotential surfaces. We then use this entropy profile and an estimated central density to construct a spherically symmetric (non-spinning) hydrostatic model, iterating on the central density until it has the correct total energy (inside of the zero-pressure surface). This automatically gives the mass contained in this object, which will be lower than our remnant mass, the remainder representing material that, due to dissipation of rotational energy, has expanded out to large distances and therefore provides negligible weight. To determine the evolution of hot spots, we order remnant particles by potential, and map them to their new positions in the final object, calculating new temperatures from the new densities, again assuming their entropy did not change (entropy is not constant over isopotential surfaces, so these temperatures are not strictly consistent with the hydrostatic model).

In Figure 17, we show the results of our evolutionary estimate for our fiducial 0.4–0.8 and 0.6–0.6 M_☉ systems. For the former case, the core-envelope, originally 0.90 M_☉, accretes 0.06 M_☉ from the disk, the remaining 0.25 M_☉ going to large distances. The central core is not significantly heated, while the lower-density hot envelope is, with the outer hot envelope along the rotational axis passing the ignition line. Since this material is almost non-degenerate, the resulting nuclear burning will likely be stable, or be extinguished by expansion. Thus, not unexpectedly, the hot envelopes of dissimilar-mass mergers are not good candidates for a nuclear runaway.

Zoom In Zoom Out Reset image size

For the 0.6–0.6 M_☉ system, the center of the final, spun-down object is at much higher density and temperature than the remnant, while much of the outer regions have become less dense and cool. The latter happens because similar-mass mergers have strong rotational support, and if this is removed their binding energy increases significantly. To compensate for this, a large amount of mass has to expand to large distances, causing the core-envelope mass to decrease from 1.11 M_☉ to 0.91 M_☉. In the final object, the hottest point on the equatorial plane does not reach the ignition line, but the significantly hotter points above and below the equatorial plane do, at densities under which degeneracy pressure still dominates. Hence, if the hot spots indeed compress with the rest of the remnant, a nuclear runaway could be triggered. (Of course, a nuclear runaway would start as soon as the heating timescale becomes shorter than the compression timescale, which may happen closer to the ignition line.)

In the right panel of Figure 16, we show the results of applying our estimates to all our merger remnants. One sees that all compress and heat, and almost every remnant whose accretor mass is above 0.8 M_☉ will reach ignition somewhere on the equatorial plane, in many cases under degenerate conditions. We also chart the evolution of the off-center hot spots in similar-mass mergers (triangles), and while they are at lower density, they remain degenerate and are all pushed substantially farther above the ignition line than their counterparts on the equatorial plane. Almost all similar-mass mergers with an accretor mass above 0.5 M_☉ could therefore experience nuclear runaways due to their hot spots, though at least some of them will become non-degenerate before an explosion can occur.

The above suggests it is at least plausible that many of our mergers would eventually ignite in degenerate conditions, and that it thus is worthwhile to simulate their evolution in detail. Suitable simulations have recently been pioneered by Shen et al. (2012) and Schwab et al. (2012). Shen et al. started with a one-dimensional simulation, where they ported the remnant of a 0.6–0.9 M_☉ merger (from Dan et al. 2011), and evolved it assuming a γ = 5/3 polytropic EOS and an α = 10⁻² viscosity. They find the system spins down completely due to outward angular momentum transport, and the rotationally supported thick disk is transformed into a tenuous, thermally supported envelope that hardly affects the core. Over longer, thermal evolution timescales (simulated using MESA; Paxton et al. 2011), this tenuous hot envelope cools, compresses the core, and lights off-center convective carbon burning, eventually turning the remnant into an ONe WD (that may end its life in an accretion induced collapse).

Schwab et al. (2012) went a step further, porting the same 0.6–0.9 M_☉ simulation, as well as seven other systems, into two-dimensional ZEUS-MP2 simulations (Hayes et al. 2006), using the Helmholtz EOS and an α = 3 × 10⁻² viscosity. They confirm the one-dimensional results, finding complete spin-down and transformation of the rotationally supported disk into a tenuous, spherically symmetric, hot envelope. They find a 50% increase in the temperature of the hottest point, and a factor of three increase in the corresponding density. They also find entropy to roughly be constant in the remnant, except in the outer regions and at the very center, where dissipation of rotational energy leads to heating.

It is encouraging that the results of the above detailed simulations are similar to what we find using our first-order estimates. For our 0.6–0.9 M_☉ remnant, our estimate gives increases for the hottest equatorial point of a factor of 2.5 in density and 1.6 in temperature, reasonably close to what is found by Schwab et al.¹⁴ Thus, our simplifying assumptions appear to be appropriate at least for dissimilar-mass mergers, where most of the rotational dissipation will be in the disk and envelope, and the structure of the remainder is roughly spherically symmetric (both in density and temperature). It is not clear if our estimates would be equally good for similar-mass mergers, where rotational dissipation should occur throughout the star, heating the entire remnant, and where there are substantial differences between the remnant's density and temperature structures. In particular, the rotational axis hot spots may dissipate, which would potentially make it more difficult for a similar-mass system to achieve a runaway. It will thus be particularly interesting to simulate the further evolution of those remnants in more detail. Unfortunately, no such remnants were included by Shen et al. and Schwab et al.

From our estimates, it seems that, as suggested by vKCJ10, many merger remnants will ignite carbon fusion. If a detonation is triggered, the resulting explosion may well resemble an SN Ia. Indeed, if the remnants spun down before ignition, their structures are sufficiently close to that of a cold WD that the calculations of Sim et al. (2010) should apply. From our estimates, for mergers that have a total mass between 1.2 and 1.4 M_☉ (which should be the most common ones), the final objects have masses between ∼0.9 and ∼1.1 M_☉ which matches fairly nicely the range of ∼1 to ∼1.2 M_☉ required to reproduce the observed range of SN Ia luminosities (Sim et al. 2010).

Of course, it is far from clear whether ignition leads to a detonation, since we do not currently understand how detonations are triggered (Seitenzahl et al. 2009; Woosley et al. 2011 and references therein). Generally, it should help that ignition in our remnants is at much lower density (a few 10⁷  g cm⁻³) than is the case for near-Chandrasekhar models (∼10⁹ g cm⁻³), because complete burning leads to much larger relative overpressures (e.g., Mazurek et al. 1977; Seitenzahl et al. 2009). Also, if a deflagration is started, plausible mechanisms to transition to a detonation all seem to require densities around 10⁷  g cm⁻³, where the conductive flame speeds are slower and the separation between the various burning fronts increases (e.g., Woosley et al. 2009, 2011).

An interesting aspect of our results is that for all cases ignition likely happens off-center: in shells for dissimilar-mass mergers and in hot spots along the rotational axis for similar-mass ones. Previous one-dimensional simulations suggested off-center ignition would lead to a slow deflagration flame that turns the CO WD into a ONe WD (e.g., Saio & Nomoto 1985). However, these calculations assumed a hot spot many pressure scale heights above the center. For ignition closer to the center, a deflagration plume is produced (e.g., Aspden et al. 2011), which may transition to a detonation (Seitenzahl et al. 2011) and unbind the star.

Given our findings, it seems likely that, if a detonation occurs, it will be triggered off-center. It would be interesting to simulate the resulting explosion, and see whether one could reproduce the observational evidence for asymmetries, which have been interpreted in terms of off-center ignition (though so far only in the context of near-Chandrasekhar models; Maeda et al. 2010a, 2010b).

Finally, while we simulated only mergers of CO WDs, we can extrapolate our results to more massive ONe WDs. For these, the temperatures would be at least as high as for our 1 M_☉ accretors, and, after further viscous evolution, the mergers should become hot enough to ignite Ne burning. If this also leads to a detonation, the lower fusion energy released would likely lead to a less energetic explosion than expected for a CO WD merger, but it would produce far more ⁵⁶Ni and have a very large mass. Plausibly, it would resemble an SN Ia like SN 2009dc, which had unusually low ejecta velocities, produced ∼1.8 M_☉ of ⁵⁶Ni and had a total ejecta mass of ∼2.8 M_☉ (Taubenberger et al. 2011).