Extreme close approaches in hierarchical triple systems with comparable masses
Niv Haim Boaz Katz


firstpage–lastpage
[
    [
    
=====





We study close approaches in hierarchical triple systems with comparable masses using full N-body simulations, motivated by a recent model for type Ia supernovae involving direct collisions of white dwarfs (WDs). For stable hierarchical systems where the inner binary components have equal masses, we show that the ability of the inner binary to achieve very close approaches, where the separation between the components of the inner binary reaches values which are orders of magnitude smaller than the semi-major axis, can be analytically predicted from initial conditions. The rate of close approaches is found to be roughly linear with the mass of the tertiary. The rate increases in systems with unequal inner binaries by a marginal factor of ≲ 2 for mass ratios 0.5≤ m_1/m_2 ≤1 relevant for the inner white-dwarf binaries. For an average tertiary mass of ∼ 0.3 M_⊙ which is representative of typical M-dwarfs, the chance for clean collisions is ∼ 1% setting challenging constraints on the collisional model for type Ia's. 





gravitation – celestial mechanics – white dwarfs – supernovae: general – planets and satellites: dynamical evolution and stability








§ INTRODUCTION


The progenitor problem of type Ia supernovae (SNe Ia) - what triggers the thermonuclear explosion of some white-dwarfs (WDs)? remains open despite decades of observational and theoretical research <cit.>. It was recently suggested that a primary channel for SNe Ia may be the direct collisions of WDs <cit.>. A successful explosion in such collisions has been established and several simple and robust properties have been shown to agree with SNe Ia observations <cit.>. 
 
The main challenge of this model is whether the collision rate is sufficient to account for the type Ia rate, which requires about 1% of WDs to explode <cit.>. The collision rate of free floating WDs in the field is off by many orders of magnitudes. Collisions in the dense cores of globular clusters are insufficient given the small 10^-4 fraction of stars that reside in these systems <cit.>. Recent advances in the study of the long term evolution of few-body systems <cit.> have lead to the realization that extremely close approaches often occur in such systems. Close approaches between WDs in multiple stellar systems may result in mergers due to the significant gravitational wave emission <cit.> or direct collisions <cit.>. If a sufficient amount of WDs are in relevant multiple systems, the rate of direct collision may be as high as the type Ia rate <cit.>. 

Is the collision rate sufficient to account for the the type Ia rate? An estimate of the collision rate requires the knowledge of the multiplicity properties of WDs and an understanding of the dynamics of multiple systems. Both of these aspects are not sufficiently known to provide a reliable answer at the moment <cit.>. WDs in multiple systems are hard to observe since their main-sequence companions greatly outshine them <cit.>. In particular, the closest known WDs - Sirius B and Procyon B, would probably not have been discovered had they been a few times farther away, and it is likely that the majority of such systems have not been discovered yet even within the local solar neighborhood <cit.>. Note that the question of how many WDs have an additional WD companion is much more well constrained and it is estimated that about 10% of WDs have a lighter WD companion on a relevant wide-orbit (0.5< period <5000yr)<cit.>. A significant improvement in the census of WDs in multiples is expected in the near future by the now operating European astrometric mission Gaia. 

On the theoretical side, the WD collision model faces the following challenge - if a given system is sufficiently active dynamically as to lead to a collision of the WDs, why didn't it lead to the collision of the main sequence progenitors of the WDs which are much bigger <cit.>? Possible ways around this can be related to changes in the system's configuration occurring on time scales which are longer than the stellar evolution time due to passing stars <cit.> or multiplicity which is higher than 3 <cit.> or to small changes in the configuration in the last stages of mass loss leading to the formation of the WDs.
Indeed, population synthesis calculations of isolated triple systems, assuming isotropic, adiabatic mass loss have obtained very low collision rates <cit.>. While this is a serious challenge, our poor knowledge of the primordial systems and mass loss processes do not allow clear conclusions. 

In this paper we focus on the simplest dynamical aspect of the collision model which has only been partly explored - the conditions for close approaches in triple systems after the inner binary stars became WDs. We use direct numerical n-body simulations to map the conditions for close approaches for the relevant parameter space of configurations and mass values of the 2 WDs and the third star. This work extends the results presented in <cit.> in three important ways: First, the dependence on the mass is studied down to low tertiary masses of 0.1  M_⊙ which are relevant to the abundant M-dwarfs. Second, we derive analytic criteria that allow the collision probability to be estimated with good confidence for the case of nearly equal mass WDs. Third, we study the dependence on the mass ratio of the inner-binary. For the purpose of this work, a new code was written in Python by N. Haim which is now publicly available.  

The structure of this paper is as follows. In section <ref> we describe the numerical integrations that are performed and show the main numerical results. In section <ref> we show the evolution in a few numerical examples  and emphasize the main characteristics of triple systems leading to extreme close approaches. In section <ref> we derive an analytic criterion for extreme close approaches in systems with equal-mass inner binaries, and in section <ref> we summarize the results and discuss the implications for the collision model of type Ia's.





§ DIRECT NUMERICAL INTEGRATIONS SHOW THAT THE CLOSE-APPROACHES RATE IS ROUGHLY LINEAR WITH M 3



We consider hierarchical triple systems of three gravitationally interacting bodies of masses m_1,m_2,m_3, where the distance from m_3 to either m_1 or m_2 is much larger than the distance between m_1 and m_2. On short time scales the hierarchical triple system behaves as two 2-body (Keplerian) systems: the smaller inner orbit of m_1 and m_2 and the larger outer orbit of m_3 and the center of mass of m_1 and m_2. There is very little exchange of energy between the two orbits and therefore both inner and outer semi-major axes (a and a_out) are almost constant. Accordingly the hierarchical arrangement of the system stays stable on very long timescales. In contrast, there is an exchange of angular momentum. The orbital parameters other than the semi-major axis, and in particular the eccentricity, change slowly on time-scales much larger than both orbital periods P and P_out.

It is useful to quantify the hierarchy of the system as the initial ratio between the pericenter of the outer orbit r_p,out and a, namely the minimal separation between m_3 and the center of mass of the inner binary, expressed in units of a. Systems with small hierarchy experience substantially more exchange of energy and they may be quickly disrupted, with one of the masses ejected to infinity. On the other hand, systems with higher hierarchy stay stable for longer periods. For example, systems of three comparable masses with hierarchies larger than 3-4 are stable for millions of orbits <cit.>.

In this paper we focus on hierarchical triples with moderate hierarchy 1.5 < r_p,out / a < 10. The inner binary includes two White Dwarfs (WD) of mass 0.6 < m_1, m_2 < 1.2 M_ with semi-major axis 1 < a < 1000 AU. The perturber is a stellar object of mass 0.1 < m_3 < 1.2 M_. Most of our runs include an equal mass binary m_1=m_2=0.6M_ (where there is a peak in the WD mass distribution, e.g. <cit.>) but we also perform simulations with inner binaries having unequal masses.

We study the conditions for extreme close approaches in the inner binary, where the inner binary separation (r=r_1-r_2) becomes smaller by orders of magnitude compared to the semi-major axis a. A system is said to have experienced a close approach if r / a becomes smaller than a pre-defined threshold. For example, in order that two white dwarfs in an orbit with a=10AU experience a collision, the separation has to be smaller than twice the white-dwarf radius or r/a < 2R_WD/10 AU∼ 10^-5.



 §.§ Numerical simulations



Our goal is to estimate the probability of achieving close approaches. This is a statistical property which can be estimated based on large ensembles of integrations with varying initial conditions. 

We run ∼120,000 full N-body simulations of three bodies using a second-order symplectic integrator with an adaptive time-step <cit.> (PTMT). Details of the integrator are provided in appendix <ref> and in <cit.>.  

Ensembles Multiple runs are prefomred for ensembles of ∼ 500-2000 randomly chosen initial conditions (see Table <ref> for exact details). All the simulations in a given ensemble have the same mass values and initial .  The initial orbital parameters of each system in an ensemble are chosen randomly from the following distribution:

Initial conditions Both the inner and the outer eccentricities are chosen uniformly at random in the range 0 < e, e_out < 0.9. The z-axis is chosen along the initial direction of the outer angular momentum and the x-axis along the outer orbit's eccentricity vector (Runge-Lenz vector). The orbital orientation of the inner orbit is chosen randomly from an isotropic distribution, with uniform distributions of 0 < Ω,ω < 2π and -1 < cosi < 1, where Ω is the longitude of ascending node, ω is the argument of periapsis and i is the (mutual) inclination. The mean anomalies are chosen uniformly at random in [0,2π]. 

Stopping conditions The evolution of the system is stopped if either

  
  * t / P_0 > 2 × 10^6, where P_0 is the (initial) inner period.
  
  * The system is disrupted. The disruption condition is that the outer orbit has a positive energy and that its separation is more than 50 times larger than the separation of the inner orbit.
  
  * In section <ref> we add a condition and stop if 5 Gyr have passed. Note that each run can be scaled in distance and time, so that this condition is enforced in post-analysis.


For each ensemble we estimate the close approach probability as the fraction of systems that experienced a close approach out of the total number of systems in the ensemble. In section <ref> we discuss clean collisions, in which only systems that did not have close approaches prior to the collision are counted <cit.>.  
Convergence is demonstrated in Appendix <ref>.

We provide the simulation code, written in Python. The performance is enhanced using the numba package for generating just-in-time machine instructions. The code can be found here: <https://github.com/nivha/three_body_integration>




 §.§ Direct estimates of the close approach probability


Figure <ref> shows the close approaches probability as a function of m_3. The systems shown have masses m_1=m_2=0.6M_ and 0.1 < m_3 < 1.2 M_ and initial hierarchies of =5 (black) and =8 (blue). We show the results for close approach thresholds of r/a < 10^-4,10^-5,10^-6. The dashed line is the analytic expected asymptotic rate of r/a→0 which is derived in section <ref>, and applied to initial conditions of 500,000 systems sampled from the same distributions as the numerical simulations.





As can be seen in Figure <ref>, the close approaches probability increases with smaller initial hierarchy. The rate is insensitive to the exact value of the close approach threshold for the relevant extreme values r/a≲ 10^-5 and is close to the analytic asymptotic value. The rate decreases with the mass of the perturber m_3 roughly linearly. 

Figure <ref> shows the close approach probability as a function of the hierarchy  for a threshold of r/a<10^-5. Results are shown for different mass choices including three inner-binary mass ratios m_2/m_1=1,7/6, 3/2, 2, 3, and two third-body mass ratios, quantified as m_3/(m_1+m_2) = 1/3,1/6 (upper and middle panels resp.). 
 
As can be seen in Figure <ref>, for the mass ranges that we consider, the dependence of the close approaches probability on the inner-binary mass ratio is not very large; For a rare case of WD-WD binary of mass ratio m_2/m_1=2 the probability is about twice as large as the rate for equal mass binary with similar perturber. For the more probable masses of WD-WD binaries, this difference will be even smaller.

The close approach probability is suppressed at small hierarchies due to the quick disruption of the systems and at large hierarchies due to the smaller range of initial conditions leading to high eccentricities.
We next focus on the dependence of the close approach probability at higher hierarchies. For a study of the effects of disruptions on the close approach rates see <cit.>.  





§ NUMERICAL EXAMPLES OF A TRIPLE SYSTEMS



In Figure <ref> are shown examples for the evolution of the orbital parameters in hierarchical triples. All three systems (shown in black, green and red) have masses m_1=m_2=m_3=0.6M_ and initial hierarchy =7.

The evolution of the pericenter of the inner binary r_p relative to the (almost constant) inner semi-major axis a is shown in the upper panel. As can be seen it is possible for r_p/a to reach very small values of <10^-2 and that the behavior is approximately periodic. The long term periodic evolution of hierarchical triple systems is well explained by the Lidov-Kozai (LK) mechanism (<cit.>) and can be approximated using the following basic assumptions:

  
  * There is no exchange of energy between the inner and outer orbit (a,a_out=const.).
  
  * The Hamiltonian is separated to the two Keplerian Hamiltonians of the inner and outer orbits and a small perturbing potential which is expanded up to quadratic terms as a function of the small ratio between the inner and outer separations.
  
  * The small changes in the Keplerian orbital parameters (e.g. Δ j, and Δ e) due to the perturbation are averaged twice: once over one outer period and once over one inner period (this is called double averaging).


Under these assumptions, the double averaged quadratic potential turns out to be axisymmetric with respect to the direction of the outer angular momentum Ĵ_out. Therefore no torques exist along this axis, and the magnitude of the outer angular momentum is conserved:


    | J_out| = const.


One consequence is that the number of constants turns out to be one less than the number of degrees of freedom and the resulting long-term evolution is periodic. Another consequence is that J has a lower bound set by the constants of motion. In fact, given that the total angular momentum | J_tot| is conserved it follows that


    | J | > || J_tot| - | J_out||.


Equation (<ref>) can be used to put a lower bound on the separation of the inner binary, which is shown in the upper panel of figure <ref> (as dashed lines). As can be seen, this lower bound works well for the systems shown in black and green but it fails for the integration shown in red. The reason for the failure is that in reality J_out is fluctuating on short time-scales and these fluctuations are not taken into account by the double-averaging approximation <cit.>. 

In order to derive a condition for close approaches, the fluctuations in J_out need to be quantified. Given that even small relative fluctuations in J_out may be important it is more convenient to quantify the corresponding fluctuations in the inner-binary's angular momentum J. In particular in the test-particle limit, where one of the inner masses is negligible, the angular momentum vector of the outer orbit is fixed. In this case, the axisymmetry of the averaged perturbing potential implies that the component of the inner angular momentum along the axis of symmetry is constant. Here it is customary to choose the z direction along the direction of the constant outer angular momentum vector. It is therefore useful (in the test-particle limit) to study the short-term fluctuations in J_z, which is fixed under the double-averaging approximation.

In the non test-particle case Ĵ_out is not conserved and it is useful to study the fluctuations in the following (double-averaging) constant which reduces to J_z in the test particle limit (e.g. <cit.>): 

    = | J_tot|^2 - | J_out|^2/2| J_out| J_circ = j·Ĵ_out + j^2 J_circ/2| J_out|,

where

    = μ√(G(m_1 + m_2 )a)

is the value of the angular momentum of the inner orbit if it were circular (the maximal possible value),

    j = J/

is the normalized inner angular momentum and j=|j|. Note that j is directly related to the eccentricity through j^2 + e^2 = 1, and that at high eccentricities the pericenter can be expressed as:

    r_peri = 1/2 j^2.



In orbit with close approaches j≪1 and ≈j·Ĵ_out<j (eq. <ref>) implying a lower limit for j:

    j ≥||.

The lower limit Eq. (<ref>) is equivalent to the lower limit (<ref>) for close approaches where |J|≪|J_tot|. To see this note that in this case, |J_tot|≈ |J_out| so that |J_tot|^2-|J_out|^2≈ (|J_tot|-|J_out|)(2|J_out|). Eq. (<ref>) then follows from Eqs. (<ref>) and (<ref>). By plugging  into Eq. (<ref>) we obtain a lower limit for the separation of the inner binary:

    r_min = 1/2^2.

Note that the exact value of the minimal separation can be calculated using all the initial parameters and the conservation of the averaged perturbing hamiltonian <cit.> and will be equal or higher than equation (<ref>).    

In the middle panel of Figure <ref> we show the evolution of  for all three systems. First, note that  is not strictly constant but rather fluctuating around a constant mean  (shown in dashed blue). Second, note that for the system shown in red, the fluctuations in  enable it to cross zero. In this case, j is not bounded from below and the inner binary achieves very small separations, as is clearly seen in the upper panel. This is the major difference between the system shown in red and the other systems.

In Figure <ref> we show an example of the evolution of the pericenter and  for the case of unequal binary masses (m_1=0.6,m_2=1.2,m_3=0.6 M_). As can be seen  is no longer constant. In this case, the quadratic expansion of the potential (assumption (b) of the LK mechanism) is not sufficient and the observed oscillations are caused by the next term, the octupole <cit.>. Due to these oscillations  can cross zero and lead to extreme close approaches <cit.>. The octupole depends on the masses through the combination (m_2-m_1) / (m_2+m_1) <cit.>. As can be seen in the bottom panel of Fig <ref> the close approaches probability increases linearly with the octupole term. Note that the short-term (outer period timescale) fluctuations have a long-term effect on the evolution <cit.>. This is not taken into account in existing secular codes <cit.>, thus, calculating the evolution of  requires more work.

In Figure <ref> the initial values of || and the minimal separations obtained are shown for all runs performed in this work. Each point in the figure represents the result of a triple system integrated to 2·10^6 inner periods. The y-axis shows the minimal separation r_min/a achieved by the inner binary throughout its evolution, and the x-axis is the initial value of || (the absolute value is taken for convenience, since  can be also negative). Red color means that  of that system crossed zero at some point during its evolution; Blue color means  did not cross zero[During simulations we compute  only at apocenters. This means that in rare cases some 'crossings' of  might have been missed (blue points that should be red)]. The three example systems from Figure <ref> are shown here as black dots. The black line is the lower limit for the separation of the inner binary from the Lidov-Kozai theory (eq. <ref>).

As can be seen in Fig. <ref>, systems whose  did not cross zero (blue) roughly follow the predictions of the LK model (black line) with a tail of systems reaching smaller separations (due to the lower values of || obtained throughout the integration). On the other hand, most systems whose  crossed zero (red) reach much closer separations. In such cases, the value for the minimal separation of the inner binary is a stochastic parameter. Assuming that each system scans its phase space uniformly (<cit.>, section II.B), we expect that after 2·10^6 inner periods, systems would reach separations of r/a∼1.5·10^-7. As can be seen in the figure, most systems where  crossed zero reached such small separations. It is therefore likely that for longer integrations, the minimal separation obtained by these systems will decrease linearly with time. Note that the median is somewhat lower than the expected value for 2·10^6 orbits, indicating that the distribution of pericenters is not strictly uniform. Further analysis is beyond the scope of this paper.


















§ A CRITERION FOR EXTREME CLOSE APPROACHES IN HIERARCHICAL TRIPLE SYSTEMS WITH EQUAL MASS BINARIES


In this section we derive an approximate analytic criterion for  to cross zero (allowing extreme close approaches) for stable triple systems in which the inner binary has equal masses (m_1=m_2).



 §.§ Estimating the maximal fluctuations in j z eff



The fluctuations in  occur within each outer orbit. The angular momentum can be written as a sum of a slowly varying mean j̅ which is assumed constant within each outer orbit and the fluctuating part Δj:

 
    j   = j̅ + Δj
    J_out   = J̅_out + ΔJ_out


The fluctuations have been calculated to linear order by <cit.>. Note that while the long term corrections in <cit.> are restricted to the test particle limit, the short term corrections are not.
The mean values of the vectors j̅ and e̅ can be calculated from the instantaneous values of j, e and the outer true anomaly f_out using equation (31)  from <cit.> (see also appendix B there). 

The fluctuations of the outer orbit's angular momentum can be obtained from those of the inner orbit using the conservation of the total angular momentum, ΔJ_out = - ΔJ, or in normalized form (see eq. <ref>):


    ΔJ_out = - Δj


Expanding  (eq. <ref>) to first order in Δj yields the fluctuating term in :

 
    =  + Δ     where, 
       = j̅·Ĵ̅̂_out + j̅^2/2J̅_out
    Δ   = (1+/J̅_out) Ĵ̅̂_out·Δj


The mean values hardly change in one outer orbit, so we can set a coordinate system of ẑ along the direction of the mean outer angular momentum Ĵ̅̂_out and x̂ along the direction of the outer eccentricity vector ê̅̂_out. In that case eq. <ref> becomes:
   
 
    Δ = (1+/J̅_out) Δ j_z


The maximal fluctuation in j_z during a close approach can be obtained using eqs. (35,33,20) from <cit.> in the limit j̅→ 0 and e̅→ 1 yielding:


    Δ j_z,max = ϵ_SA(1 - e̅_z^2) 15/8(1+2√(2)/3e_out) ,

where ϵ_SA is a dimensionless parameter that sets the scale of the fluctuations and is given by

    = (a/a_out)^321/(1-e_out^2)^32m_3/√((m_1+m_2)(m_1+m_2+m_3)) ,

and e̅_z^2 is the projection of e̅ onto the direction of Ĵ̅̂_out. The value of e̅_z^2 can be obtained from the initial conditions using the following conserved quantity of the double averaged equations which is approximately constant (a consequence of the conservation of the double averaged potential, <cit.>):


    C = -e̅^2+5/2e̅_z^2-1/2j̅_z^2 .


The conserved quantity C is computed from the initial conditions using Eq. (<ref>).

In the vicinity of close approaches, we can assume the limits j⃗→0 and e⃗→1 to obtain 

    C = -1 + 5/2e̅_z^2 ,

which allow e̅_z to be inferred from C. 

By calculating C from the initial conditions using Eq.  <ref> and using equations <ref>, <ref> , <ref> we obtain an expression for the maximal value of the fluctuations in  as a function of the initial conditions:
 
 
    =   
       3/8ϵ_SA(1+/J̅_out) (3+2e̅^2-5e̅_z^2+j̅_z^2)(1+2√(2)/3e̅_out)|_t=0

All parameters in eq. <ref> are obtained from the initial conditions of the system. A code for calculating the maximal fluctuations using this equation is provided and described in section  <ref>.




 §.§ The criterion


As shown in section <ref>, extreme close approaches occur when  crosses zero. Using the approximations described in section <ref> the condition that  crosses zero is that 

    >    (at t=0)

where  and  are calculated from the initial conditions (an implementation in code is provided in section <ref>). A comparison of the analytic criterion (<ref>) is compared to the results of numerical simulations in Figure <ref>. Each dot in the figure represents the results of a simulation of a triple system of masses m_1=m_2=0.6M_ and 0.1<m_3<1.2M_ and hierarchies =5,6,7,8,10. The other initial conditions are sampled from the distributions described in section <ref>. Red (blue) dots correspond to systems in which  crossed (did not cross) 0. As can be seen, the criterion eq. (<ref>) succeeds at predicting whether or not  would cross zero during its evolution to an excellent approximation. Systems whose  crossed zero would achieve extreme close approach of the inner binary, as confirmed in Figure <ref>.




In Figs  <ref> and  <ref> the analytic fractions, based on the criterion  <ref>, are shown (labeled analytical). The analytic fractions are obtained by sampling initial conditions for 500,000 triple systems (using the same distribution as the simulations, described in section <ref>), and checking for each system whether or not the initial conditions meet the criterion. As seen in the figures, the criterion allows a good estimate to the close approaches fraction for triples with equal mass inner binaries at above-moderate hierarchies ≳ 5.

For systems with small hierarchies, the system may be quickly disrupted before a close approach in the binary is attained. In these cases the analytic criterion, which is derived based on the assumption that the systems stay intact forever, predicts a close approach fraction which is too high, as can be seen in Fig. <ref> at low hierarchies. For systems with unequal mass binaries, the octupole term adds another oscillation to  (see figure <ref> and discussion in the text) thus increasing its chance to cross zero. In these cases the fraction predicted based on the analytic criterion is too low. Estimating the fluctuations in  due to the octupole term is beyond the scope of this paper. However, for this case, criterion <ref> can be used to compute a lower bound to the close approaches fraction.




 §.§ Implementation of the criterion in computer code



Criterion <ref> can be computed from the initial conditions of a triple system by computing:  from equation <ref>,  from equation <ref>, and j̅ from equations (20), (31) and appendix B in <cit.>. All the parameters and the criterion can be computed using a Python file which is attached to this paper.









§ SUMMARY AND DISCUSSION




 §.§ Summary

It was shown by <cit.> that triple systems consisting of an equal-mass inner white-dwarf binary and a comparable mass tertiary at intermediate hierarchies 3≤≤10 have a few percent chance of resulting in extreme close approaches r/a≲ 10^-5 and the collision of the white-dwarfs. In this paper we reproduced and significantly extended the study of the conditions for such extreme close approaches, both analytically and numerically using a new code written by N. Haim which is now publicly available (see <ref>). As described in section <ref> we extended the numerical experiments for a wide range of inner and tertiary masses. We showed that extreme close approaches are possible for low mass M-dwarf tertiaries (down to  at least 0.1  M_⊙) but with a declining probability that is roughly linear with the mass (see figures <ref>,<ref>). For equal mass inner binaries in stable systems (usually ≥5), we derived in <ref> an analytic criterion (Equations <ref>,<ref>, with an implementation provided in <ref>) that allows extreme close approaches to be predicted from the initial conditions with good confidence (see Figs <ref>,<ref> for individual runs and Figs <ref>,<ref> for statistical comparisons). To achieve this, the short-term (outer period) fluctuations of the orbital parameters were calculated using the results of <cit.>. Finally, we explored a wide range of unequal mass inner binaries and demonstrated that the collision probability increases linearly as a function of the octupole mass coefficient (m_2-m_1)/(m_2+m_1) (see figure <ref>).



 §.§ Implications for the collision model of type Ia supernovae


In principle, any close approach of two white dwarfs may lead to interaction. However, in cases of a slowly decreasing inner binary pericenter, energy dissipation or a grazing encounter may affect the orbits and suppress the chances of experiencing a direct collision. It is therefore useful to consider clean collisions <cit.>, in which all close passages prior to the collision itself are sufficiently far to ignore such effects. In other words, in all early close passages, the separation of the inner binary is always larger than a factor r>R_dissip. Only in its last passage the separation becomes small enough to allow for a collision r<R_col. 

In figures <ref> and <ref> the clean collisions probabilities are shown for different masses and a wide range of inner-semi major axis 1 < a < 1000AU, assuming a log-uniform distribution of the initial hierarchy in the range 3≤≤10 <cit.>. In our implementation for clean collisions we use R_col=2R_WD and R_dissip=4R_WD, where R_WD=10^9cm. 

These results tighten the constraints on the collision model as a primary channel for type Ia supernovae. About 0.01 of WDs need to explode within a Hubble time to account for the SN Ia rate <cit.>. Given that only 10% of WDs are likely to have a (lighter, wide-orbit) companion WD <cit.>, collision probabilities in the relevant systems need to be at least ∼ 10%. As can be seen in figures <ref> and <ref>, the probability seems to be too low by a factor of few for ∼ M_⊙ tertiaries and an order of magnitude too low for low mass tertiaries M≲ 0.2 M_⊙. The problem may be much worse if only a small fraction of double WDs have tertiaries and if the relevant high inclination systems lead to interaction before the stars become WDs <cit.>. 

In order to estimate the collision rate, the multiplicity of white-dwarfs needs to be reliably measured. Massive main sequence tertiaries M≳ 0.5 M_, for which the collision probabilities are higher, are bright M_V<10, M_K≲ 6 <cit.> and can be detected with high completeness by adaptive optics surveys of intermediate mass stars <cit.> or as common-proper motion companions to A-stars and WDs <cit.>. Direct observational constraints can thus be obtained on such triple systems but this is beyond the scope of this paper. Finaly, we note that while the results of this paper suggest that the collision probabilities in triple systems with low-mass tertiaries are low, the role of low-mass tertiaries in enhancing collisions in higher multiplicity systems  <cit.> remains to be explored. 














§ ACKNOWLEDGEMENTS


We thank Liantong Luo, Doron Kushnir, Subo Dong, Andrew Gould, Silvia Toonen, and Hagai Perets for useful discussions. This research was supported by the ICORE Program (1829/12) and the Beracha Foundation. The computations were performed with a high performance computing facility which is partly supported by the Israel Atomic Energy Commission - The Council for Higher Education - Pazi Foundation and partly by a research grant from The Abramson Family Center for Young Scientists.







mnras
 










§ RELATION OF J Z EFF TO THE MUTUAL INCLINATION


In the test particle approximation, the condition =0 reduces to j_z=0 or a mutual inclination between the inner and outer orbit of 90 degrees.
It is therefore useful to express the relation between  and the mutual inclination i in the non-test particle case. Eq. <ref> can be expressed as:

 
    = j ( cosi + j /2| J_out|) .


The condition =0, where close approaches are expected, requires initial retrograde configurations i>90 as can be inferred also by the equivalent condition |J_tot|=|J_out| <cit.>.




§ INTEGRATION METHOD



 Similar to <cit.> we use the PTMT, second-order, symplectic integrator with adaptive time-step <cit.>. The propagation of the positions and velocities of the bodies x and v to the next time step is performed in a two stage leapfrog scheme:
 
 
    v_next   =v + a(x) · dt_0 ·( U(x)/U_0)^-3/2
    x_next   =x + v_next· dt_0 ·( E_0-K(v_next)/U_0)^-3/2

Where a are the accelerations, E_0 is the initial total energy of the system, K and U are the  kinetic and potential energies respectively and dt_0 and U_0 are constants. U_0 is chosen as U_0 = -Gm_1m_2 / a_0 where a_0 is the initial semi major axis. The time-step amplitude is parametrized as dt_0=P_in,0/N_s, where P_in,0 is the initial period of the inner binary and for most of our runs we use N_s=1000. For the convergence check (section <ref>) we also use N_s=250,50.




§ CONVERGENCE



In Figure <ref> we show the close approaches probability (the fraction of systems, out of the total systems for the same parameters, whose inner binary separation became smaller than the semi major axis by 10^-5) for triple systems with masses m_1=m_2=0.6M_ and 0.1 < m_3 < 1.2 M_. We show the results for runs with different time-step amplitudes with N_s=1000,250,50 shown in solid, dashed and dotted lines respectively, see appendix <ref>), and for two initial hierarchies r_p,out=5,8. The solid lines are the same lines as in Figure <ref> for the relevant parameters. As can be seen fast convergence is achieved as a function of N_s.

It is important to note that most of the simulations are not individually converged. In figure <ref> the results of two integrations are shown which are identical in all parameters except for position of m_1 along the x-axis, which was modified by 10^-8a. The figure shows the evolution of the difference between the positions of m_1 along the x-axis in the two integrations. As can be seen, after a few thousand inner periods, the difference grows to order unity. The Lyapunov timescale for this system is thus of order a few hundreds of inner orbits on average. In our simulations we integrate each system to millions of inner periods, and it is therefore futile to aim for the convergence of an individual integration. Having said that, As can be seen in Fig. <ref>, The statistical properties of the systems, such as the close approaches probabilities, are indeed converged to within statistical error. 

Another issue with the integrator is that for some systems with very close approaches the time-step became so small that the integration was not completed. For the small-time step amplitudes that we used with N_s=1000, this happened for only a handful of systems (out of ∼100,000) and we did not consider them in the results. For cruder time resolutions (smaller N_s), this happens more frequently. For our convergence tests, out of a total of ∼10,000 runs for each N_s, the number of such cases were 70 for N_s=250 and 910 for N_s=50. However, even for N_s=50, where 10% of the runs were stuck, the overall convergence statistics that we show is not significantly affected.








§ SUMMARY OF THE PROPERTIES OF THE NUMERICAL ENSEMBLES


A summary of all numerical ensembles that were calculated in this work is given in table <ref>.
For low hierarchies 1.5≤≤4 all N_batch=500 sampled initial conditions are integrated. For higher hierarchies 4.5≤≤10, N_batch=2000 initial conditions are sampled but systems with large initial || are assumed to not lead to close approaches and only a subsample satisfying a prescribed upper limit ||<||_max are simulated (the values of ||_max are given in table <ref> and the number of simulations actually performed is given in table <ref>). The range of || is empirically chosen so that the vast majority of the systems that experience close approaches lie deep within the region of || < ||_max. In Figure <ref> the initial values of || for the systems for different initial hierarchies  are shown. Systems that  experienced (did not experience) a close approach of r/a<10^-5 are shown in red (blue). As can be seen, most systems that experience a close approach have small initial =0. Using such plots we verify that the chosen values ||_max for each ensemble captures the majority of systems that experience a close approach. We emphasize that this is a conservative approach since all systems that are not simulated are assumed to not lead to a close approach.
















