Monte Carlo population synthesis of post-common-envelope white dwarf binaries and type Ia supernova rate
MMonte Carlo population synthesis of post-common-envelope whitedwarf binaries and type Ia supernova rate
Iminhaji Ablimit , Keiichi Maeda , and Xiang-Dong Li , ABSTRACT
Binary population synthesis (BPS) study provides a comprehensive way to un-derstand evolutions of binaries and their end products. Close white dwarf (WD)binaries have crucial characteristics in examining influence of yet-unresolvedphysical parameters on the binary evolution. In this paper, we perform MonteCarlo BPS simulations, investigating the population of WD/main sequence(WD/MS) binaries and double WD binaries, with a publicly available binarystar evolution code under 37 different assumptions on key physical processes andbinary initial conditions. We considered different combinations of the bindingenergy parameter ( λ g :considering gravitational energy only, λ b : considering bothgravitational energy and internal energy, and λ e :considering gravitational energy,internal energy, and entropy of the envelope, the values of them derived with theMESA code), CE efficiency, critical mass ratio, initial primary mass functionand metallicity. We find that a larger number of post-CE WD/MS binaries intight orbits are formed when the binding energy parameters is set by λ e than thecases adopting the other prescriptions. We also find effects of the other inputparameters on orbital period and mass distributions of post-CE WD/MS binariesas well. Containing at least one CO WD, the double WD system evolved fromWD/MS binaries may explode as type Ia supernovae (SNe Ia) by merging. Inthis work, we also investigate a frequency of two WD mergers and compare it tothe SNe Ia rate. The calculated Galactic SNe Ia rate with λ = λ e is comparablewith observed SNe Ia rate, ∼ . × − yr − – ∼ × − yr − depending on the Key Laboratory for Optical Astronomy, National Astronomical Observatories, Chinese Academy ofSciences, Beijing 100012, China Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo, 5-1-5Kashiwanoha, Kashiwa, Chiba 277-8583, Japan Department of Astronomy, Nanjing University, Nanjing 210046, China Key Laboratory of Modern Astronomy and Astrophysics, Ministry of Education, Nanjing 210046, China a r X i v : . [ a s t r o - ph . H E ] M a y ∼ . ≥ M (cid:12) andmass ratio > .
8, results in a much lower SNe Ia rate than observed.
Subject headings: binaries: close – stars: evolution – white dwarfs – stars:supernovae – X-rays: binaries
1. Introduction
The mass transfer process and the common envelope (CE) phase are crucial for pro-ducing all kinds of compact star binaries. Among the compact star binaries, white dwarf +main sequence (WD/MS) binaries or double WD systems serve as excellent laboratory tounderstand yet-unclarified physics of binary evolution during a common envelope phase, andof type Ia supernovae (SNe Ia). MS/MS close binaries generate close compact star binariesmainly by going through a CE phase at least once. In the MS/MS binary (the orbital sep-aration between ∼
10 and 1000 R (cid:12) ), once the more massive one (primary) evolves into thefirst giant branch (FGB) or the asymptotic giant branch (AGB), it fills its Roche lobe andthe mass transfer may be dynamically unstable, thus the envelope of the primary engulfs theless massive one (secondary) (Iben & Livio 1993; Webbink 2008). Paczy´nski (1976) proposedthat the orbital energy and orbital angular momentum are removed by the CE ejection. Thestellar parameters of detached WD/MS binaries are easily observed among all compact starbinaries, thus post-CE WD/MS binaries (PCEBs) are ideal systems to test yet-unresolvedphysical processes involved in the CE phase.Binary population synthesis (BPS) is a very useful tool to study the physical parametersand the evolution processes involved in the formation of various types of close binaries.After de Kool & Ritter (1993), a bunch of BPS works have been performed on PCEBs withdifferent CE phase models. The CE evolution is crucial for producing compact star closebinaries. However, the CE evolution is still not well understood. The α -formalism and γ -formalism have been widely used to study the CE evolution, the former one considers energyconservation (Webbink 1984; Dewi & Tauris 2000) while the later one considers angularmomentum conservation (Nelemans & Tout 2005). In the α -formalism, the outcome of theCE phase is related with the efficiency parameter α CE , E bind = α CE ∆ E orb (1)where E bind and ∆ E orb are the binding energy of the envelope and the change in the orbitalenergy during the CE phase, respectively. The two stars coalesce if the system satisfies the 3 –equation 75 of Hurley et al. (2002) and the CE evolution is longer than dynamical timescale(for details of the criterion for surviving or merging during the CE phase see Hurley et al.(2002)). The binding energy of the envelope is expressed as the following: E bind = − GM M en λ R (2)where M , M en and R are the total mass, envelope mass and radius of the primary star,respectively. It is commonly believed that there is no mass accretion during the CE phase.The binding energy parameter, λ , depends on the mass and evolutionary stage of the primarystar. In previous BPS studies of PCEBs (de Kool & Ritter 1993; Politano & Weiler 2006,2007; Wang et al. 2015), the binding energy parameter has been treated merely as a constant(0.5 or 1.0), or it is fixed with α CE (Toonen & Nelemans 2013). However, a number of papers(Han et al. 1994; Dewi & Tauris 2000; Webbink 2008) claimed that assuming a constant valuefor λ is not a promising way to address all types of primary stars and their evolution phases.They found that the value of λ changes as the star evolves. According to Davis et al. (2010)and Xu & Li (2010), the values of λ for evolving stars can be calculated by consideringgravitational energy only (hereafter λ g ), adding internal energy ( λ b ), or adding entropy ofthe envelope ( λ e ). For a 6 M (cid:12) star, the values of λ g , λ b and λ e range from 0.11 and 0.64, 0.4and 1.5, 0.6 and 9.98, respectively. We emphasize that all these parameters (e.g. α CE , λ , thecritical mass ratios q cr ) (see § λ depends on thestructure of the primary star, α CE may change with the secondary mass (de Marco et al.2011). Therefore, we need further simulations of PCEBs and additional observational datato clarify the nature of the CE evolution.PCEBs with a relatively massive MS secondary may produce type Ia supernovae (SNeIa) in the two most popular scenarios proposed so far: one is the single degenerate model 4 –(SD), and the other is the double degenerate model (DD) (Wang & Han 2012, for a re-view). Most of PCEBs can evolve (either with an unstable or a stable mass transfer) intodouble WD binaries. The DD scenario considers that a part of the double WD systemsmerge via gravitational wave radiation on timescales shorter than the age of the universeas a potentially significant population leading to SNe Ia (Webbink 1984). The explosionmechanism involves complicated interaction between hydrodynamics evolution and nuclearreactions. Several models exist that vary in the physical processes leading to the explosion:for example, the thermonuclear flame can be either detonation or deflagration, the mass ofthe immediate progenitor WD can either be the Chandrasekhar mass or sub-Chandrasekharmass (Hillebrandt& Niemeyer 2000). Moreover, there is a serious concern about the capa-bility of the merging product to explode as an SN Ia, since the remnant could lead to aformation of a neutron star rather than an SN Ia (Nomoto& Iben 1985; Piersanti et al. 2003;Pakmor et al. 2013; Dan et al. 2014; Ablimit & Li 2015). In sum, both observationally andtheoretically, the exact nature of the SNe Ia progenitors remains unclear. Studying the birthrate and delay time distribution can help understand the binary evolution channel towardSNe Ia. Indeed, the birth rate and delay time distribution of the observed SNe Ia are notredily satisfied by either of the SD or DD model (Maoz et al. 2014, for a recent review). Inmost theoretical works including CE evolution, the standard merger model for SNe Ia, i.e.,merging two CO WDs, is not able to reproduce the observed SN Ia rate (Ruiter et al. 2009;Mennekens et al. 2010; Toonen et al. 2012; Chen et al. 2012; Bours et al. 2013).The discrepancy between the theoretical and observational SN Ia rates is an open ques-tion. Claeys et al. (2014) considered effects of various physical parameters on the binaryevolution on the rate of SNe Ia in the SD and DD models. Their BPS results show thateven for their optimistic models the predicted SN Ia rate is a factor of three less than thegalaxy-cluster SN Ia rate. Chen et al. (2012) used different mass transfer models to computethe SNe Ia rate, with a constraint from a sample of double CO WDs which can merge withinHubble time. Even though they arbitrarily relaxed the critical mass ratio leading to an SN Iain a violent merger scenario from the value obtained by hydrodynamic simulations (Pakmoret al. 2013), their BPS models also could not get the SN Ia rate as large as what is observed.Recently, He WD donors have been considered as a possible important channel for SNe Ia,given that the number of He WDs can be much higher than that of CO WDs. If mergersof a CO WD with a He WD are hypothesized as a SN Ia progenitor in the DD scenario, itcould be a major contributor to the SN Ia population (Napiwotzki et al. 2007; Badenes &Maoz 2012; Ruiter et al. 2011; Ruiz-Lapuente 2014).In this paper, we investigate the influence that a series of binary evolution prescrip-tions (especially λ and q cr ) has on the formation of WD binaries (WD/MS and double WDsystems). We consider both mergers of CO WDs with CO WDs and those of CO WDs 5 –with He WDs to compute the SNe Ia rate in our BPS study. In §
2, we describe our BPScode and our models to treat the binary physical processes including the CE evolution. Theproperties of WD/MS binaries produced by our models are presented in §
3. The resultingdouble WD systems and SNe Ia rate are discussed in §
4. Finally, the paper is closed in §
2. Monte Carlo BPS simulations
We use the BPS code developed by Hurley et al. (2002) and modified by Kiel & Hurley(2006) to generate 10 initial MS/MS binaries for each model. For the distribution of themasses of the primary stars, we adopted the initial mass function (IMF) of Kroupa et al.(1993), f ( M ) = M /M (cid:12) < . . M − . . ≤ M /M (cid:12) < . . M − . . ≤ M /M (cid:12) < . . M − . . ≤ M /M (cid:12) , (3)The masses of the secondary stars are determined by the distribution of the initial massratio, n ( q ) = (cid:26) q > µq ν ≤ q < , (4)where q = M /M , µ is the normalization factor for the assumed power law distributionwith the index ν . We consider two cases for the initial mass ratio distribution (IMRD): aflat IMRD ( ν = 0 and n ( q ) =constant) and an IMRD proportional to q ( ν = 1). For thedistribution of the initial orbital separation, a i , we adopt the following formalism (Davis etal. 2008): n ( a ) = (cid:26) a i /R (cid:12) < a i /R (cid:12) > . a i − ≤ a i /R (cid:12) ≤ (5)We assume that all binaries are in circle orbits. We calculate Pop. I and II binarieswith the metallicity given as Z = 0 .
02 and 0.001, respectively.Regarding the key physical processes, our simulations have three main tunable param-eters, i.e., α CE , λ and q cr . The main aim of this paper is to investigate the effects that 6 –these prescriptions have on the evolution toward WD/MS binaries and double WD binaries.As introduced in §
1, the CE evolution is an important but unsolved phase in the binaryevolutionary process. In most works treating the CE evolution, the two main parametersdescribing the CE evolution ( α CE and λ ) have been set as constants. However, in reality,they should change with the evolutionary process. de Marco et al. (2011) and Davis et al.(2012) proposed that α CE may change with the WD mass, the secondary mass, the massratio, or the orbital period. Here we consider the following formula,log α CE = (cid:15) + (cid:15) log q (6)in Davis et al. (2012) with the values of (cid:15) and (cid:15) taken from their Table 6. We calculatethe CE evolution either by adopting this equation or by fixing α CE = 1.There have been three prescriptions proposed so far to describe the binding energyparameters ( λ g , λ b , and λ e : see § q cr ) is another physical key parameter that determines whetherthe mass transfer is stable or not. Shao & Li (2014) computed the critical mass ratioconsidering the possible response of the accreting star (i.e., spin-up and rejuvenation) underthree different assumptions: (1) Half of the transferred mass is accreted by the secondary,and the other half is lost from the system, also taking the specific orbital angular momentumof the accretor (Also see de Mink et al. 2007). (2) The transferred mass is assumed to beaccreted by the secondary unless its thermal timescale ( τ KH ) becomes much shorter thanthe mass transfer timescale ( τ ˙ M ). The accretion rate is limited by–[min(10( τ ˙M /τ KH ) , (Hurley et al. 2002). Rapid mass accretion may drive the accretor out of thermal equilibrium,which will expand and become overluminous. Shao & Li (2014) found the values of τ KH are usually much lower than that of the same star in thermal equilibrium. Therefore, it isalways as τ KH < τ ˙ M , and the mass transfer is generally conservative. (3) The accretionrate onto a rotating star is reduced by a factor of (1 − Ω / Ω cr ), where Ω is the angularvelocity of the star and Ω cr is its critical value. In this prescription, a star cannot accretemass when it rotating at Ω cr . The remaining material is ejected out of the binary by theisotropic wind, and it takes away the specific orbital angular momentum from the accretor.The critical mass ratios are denoted as q cr1 , q cr2 and q cr3 , respectively. For each model asdescribed above, we run the BPS simulations with these different treatment of the criticalmass ratio (i.e., adopting q cr = q cr1 , q cr2 , or q cr3 ). For models 2–13, there are thus 36 differentsimulations. Including the standard model (model 1), we cover 37 different models in ournumerical calculations. In the Table 1, we summarize our models. In our standard model,we assume α CE = 1, λ = 1 for the CE phase, Z = 0 .
02, the default q cr prescription in theBSE code and a flat IMRD for the initial conditions. For other initial physical inputs, we 7 –adopt the default values in the BPS code given by Hurley et al. (2002).
3. Results3.1. Post common envelope WD/MS binaries (PCEBs)
Figure 1 shows the orbital period distributions of PCEBs in our standard model (model1) and other twelve models (from model 2 to 13) with three different prescriptions for q cr .The dashed black line, shown in all panels of Figure 1, shows the result of our standardmodel. In the upper left panel, the solid, dotted, dashed colored lines show results of model2, 3, and 4, respectively. In each model, there are three lines corresponding to three differentprescriptions for q cr . Hereafter, different line-styles are used to describe results adoptingdifferent prescriptions for λ , and different colors are used to describe results adopting differentprescriptions for q cr . It is hard to construct reliable observational distributions, but it is stillinteresting to give some observational information. We use the solid black line to demonstratethe observed PCEBs sample from Tables 1 and 2 of Zorotovic et al. (2011).From the upper left panel we see that the different prescriptions for q cr do not have cleareffects on the orbital period distribution of WD binaries. However, the orbital distributionis sensitive to the prescription for λ . The orbital separation is the smallest for λ = λ e ,while the largest for λ = λ g , and the difference in the typical orbital period is more than anorder of magnitude. There are a larger number of PCEBs with short orbital periods whenadopting λ = λ e , where the shortest orbital period is ∼ .
008 day. This means that the morePCEBs can survive CE evolution when λ = λ e (or more precisely α CE ∗ λ is higher. Mostof observed PCEB samples have short orbital periods, and their orbital period range seemsto be covered by the calculated results with λ = λ e and model 1. The upper right panelshows the results in the models 5–7, where α CE changes with q , while other parameters aresame as the models in the upper left panel. The prescription of α CE affects the distributionclearly. If λ = λ g or λ = λ b , as seen by comparing the right panel with the left panel, thedistribution becomes narrower and more sharply peaked. The distribution is however notsensitive to the prescription for α CE if λ = λ e . In the lower left panel of Figure 1, the resultsin models 8-10 are shown where the initial mass ratio distribution is given as ∝ q , while theother parameters are same with models 2–4. The orbital period distribution is similar to thecorresponding models 2–4, showing that it is not sensitive to the initial mass distribution.In the lower right panel, the results in models 11–13 with low metallicity (Z = 0 . λ e leads 8 –to the large number of short orbital period systems.The distributions of the secondary masses are given in Figure 2. The results from models2, 3, and 4 are similar to those from the standard model (in the upper left panel of the Figure2). Nearly all PCEBs (in the solid black line) have low mass secondaries. We can also seein the upper right panel that the different prescriptions for α CE (models 5–7) do not lead toclear difference, as compared to the upper left one (models 2–4 and the standard model). Formodels 8–10 (in the lower left panel), the distribution of the secondary mass moves towarda larger value. That is, more massive secondary stars are produced and can dispel the CEwhen the IMRD is given by n ( q ) ∝ q rather than adopting the flat IMRD distribution.For models with low metallicity (models 11–13, in the lower right panel), the low masssecondaries become more abundant than in the standard model. While the prescriptionsfor q cr hardly affect the secondary mass for the solar metallicity, it has some influence onthe distribution when Z = 0 . λ = λ g , λ b and λ e respectively, the peak in the WDmass distribution shifts to a lower value and low-mass WDs become more abundant. Fromthe distributions of used PCEB samples and our results, it is seen that most PCEB samplescontain low mass WDs as well. The prescription for q cr has no clear effect. From the resultsof model 2-10 and our standard model we can see that massive WDs are more abundantwhen λ = λ g and λ = λ b (especially in model 7), while low mass WDs are more abundantwhen λ = λ e (a larger number of WD binaries are produced when λ = λ e is adopted ratherthan λ b or λ g ). For models with low metallicity (models 11–13) the distribution moves tothe right (toward more massive WD) than the solar metallicity models (models 1–3). Thisshift in the WD mass is likely caused by different evolutionary age and the RLOF momentof the different metallicity star.Figure 4 shows the probability distributions in the orbital period–WD mass (left), theorbital period–secondary mass (middle), and the secondary mass–WD mass (right) planes.Shown here are for the standard model and models 2–4. As we take into account threedifferent prescriptions for λ ( λ g , b , e ) and three different prescriptions for q qr ( q cr1 , , ), there 9 –are 10 models shown in Figure 4. From the left panel we see that systems with short orbitalperiods are most abundant when λ = λ e . Also, low-mass WDs become most abundantwhen λ = λ e . From the middle panel, we see that the relation between the period and thesecondary mass is quite universal, which is not affected substantially by the treatment ofthese unresolved physical processes. In the right panel, the relations of the masses of thetwo binary members are different for different parameters, and more binaries survive CEevolution with λ = λ e . We let all WD/MS binaries continue their evolution to double WD binaries. The closedouble WD binaries which can merge within the Hubble time are produced by going throughthe CE evolution at least once. They may evolve into the CE phase during the first or secondmass transfer or may have CE evolution twice. There is possibility that some systems mayindeed explode as an SN Ia within the SD scenario e.g., (Li & van den Heuvel 1997; Wang& Han 2012; Ablimit et al. 2014), but we only consider double degenerate systems. Figure 5shows the orbital distribution of binaries containing a CO WD primary and either a CO orHe WD secondary. These binaries have experienced the CE phase at least once. Therefore,the orbital periods of these systems are less than 40 days as seen in Figure 5. It is seenthat a larger number of double WD systems have a short orbital period when the value of λ is larger (i.e., when λ = λ e is adopted). The number of WD binaries with a short orbitalperiod also increases for low metallicity. We donot know the full observational distributionof double WD systems yet, however, we have a number of confirmed samples. The solidblack lines (in Figs.5, 6, 7) show the distributions of the observed double WD binaries whichcan merge within the Hubble time (See tables of Marsh 2011 & Kaplan 2010), and they havesimilar orbital period, first WD mass and mass ratio distribution ranges as our results (seeFig.6 and 7 for WD mass and mass ratio distributions).Figure 6 shows the distribution of the primary CO WD masses. For all the models,most of the primary CO WD masses are in the range of 0 . − . M (cid:12) . As described above,the double WD may go through first stable or unstable mass transfer (for the second masstransfer we follow Hurley et al. (2002), and we use the same prescription for alpha as inthe first CE phase, if it is unstable), thus the first stable mass transfer contributes to thedistribution of the primary CO WD masses, and the distribution in this Figure more or lessdiffers from that of Figure 3. The various parameters influence the primary mass distributionto some extent, but generally a global pattern in the distribution of the primary CO WD 10 –masses is not sensitive to these assumptions in the BPS study. The distribution of the massratios between the secondary WD and the primary WD is given in Figure 7, the distributionhas double peaks for almost all the models.Santander-Garc´ıa et al. (2015) analyzed the observed data of the central object in thePlanetary star Henize 2-428, and claimed that Henize 2-428 has a double WD system witha mass ratio of nearly unity. They also reported its combined mass (1.76 M (cid:12) , which is wellabove the Chandrasekhar limit mass) and its short orbital period (4.2 hours). Based onthese values, they suggested that the system should merge within 700 million years, beingthe first candidate progenitor of the super-Chandrasekhar-mass channel in the context of theDD model of SNe Ia. On the other hand, Garc´ıa-Berro et al. (2015) reanalyzed the results ofSantander-Garc´ıa et al. (2015) and suggested that the central object of Henize 2-428 is notlikely a double degenerate system. In particular, Garc´ıa-Berro et al. (2015) argued that thepossibility of forming double WD systems with mass ratio of ∼ In the DD scenario, the orbit of a double WD binary shrinks through the gravitationalwave emission and eventually the two WDs merge. The timescale of this process is given asfollows (Landau & Lifshitz 1971). t GW = 8 × × ( M + M ) / M M P / year , (7)where P is the orbital period (in units of hour), M and M are the masses of the primaryWD and the secondary WD (in units of M (cid:12) ), respectively. Recently, it has been proposed,both theoretically and observationally, that both the SD and DD channels might have own(non-negligible) contributions to produce SNe Ia. As for the SNe Ia birth rate, it has beensuggested from BPS studies that the prediction from the DD scenario is closer to the observedrate, but there is still some difference between theory and obervation (Claeys et al. 2014).Hereafter, we focus on the DD scenario in this paper, and a consistent modeling of the SDand DD scenarios is beyond the scope of this paper. 11 –The criteria for a DD system to explode as an SN Ia have not been clarified yet, exceptfor the requirement that the system should merger within the Hubble time. Adding to thetotal mass of the system, a mass ratio of the two WDs ( q = M /M ) could also be animportant factor, but it depends on yet-unclarified explosion mechanism. Different criteriaare adopted in different BPS studies. In this section, we investigate how the BPS parametersaffect the predicted SN Ia rate, by adopting the same criteria as Chen et al. (2012). Weconsider three conditions as follow: (1) the combined mass ≥ the Chandrasekhar limit mass(1 . M (cid:12) in this work), (2) q = M /M ≥ /
3, (3) the two WDs must merge within the Hubbletime. Note that while Chen et al. (2012) obtained the SN Ia rate in their conservative modelas consistent with the observation, the adopted criterion on q might indeed be optimistic.Also, Chen et al. (2012) found that the evolution of the SN Ia rate with time did not fit theobservation. The aim of this section is to clarify how these conclusions would be affected bythe different choice of the BPS parameters.We compute the expected SN Ia rate under the DD scenario, including those frommergers of two CO WDs and mergers of a CO WDs with a He WD. Galaxies have complicatedstar formation histories, but in this paper we adopt two simple models for demonstration; (1)a constant star formation rate (SFR) over the past 13.7 Gyr, and (2) a single star burst (i.e.,a delta function). For the case of the single star burst, we assume that the burst producethe stellar mass of 10 M (cid:12) . In the case of the constant SFR, we assume it to be 5 M (cid:12) yr − (Willems & Kolb 2004).Figure 8 displays the evolution of the SN Ia birthrates with time, for the two modelsof the star formation history. The colors and styles of lines are the same as those used inthe previous figures. In the left panels for the single star burst, we also show a fit with theuncertainty of the delay-time distribution (DTD) inferred from observations (the three blackthin-solid lines, Marsh 2011; Maoz et al. 2012, see). Specifically, the middle line representsthe one obtained with the formula for the DTD, φ ( t ) = 4 × − SNyr − M (cid:12)− ( t ) − (Maoz & Mannucci 2012), and the upper and lower lines are for the ±
50% uncertainties.As a general prediction from our BPS simulations, we note that the DD model could resultis a power-law like distribution, but it is not described as simple as a single power law asfrequently attributed to the DD scenario. The peak in the DTD results from the peak inthe orbital periods in the DD systems (Figure 5). Hereafter in comparing the observed DTD(for the single burst case), we mainly focus on the peak in the predicted rate. In Figure 8,the upper two panels show the results for the models with α CE = 1 and solar metallicity(models 2-4). For a constant SFR case, the predicted birthrates of SNe Ia are lower thanthat of the standard model ( ∼ . × − yr − ). The lowest birthrate among these models( ∼ . × − yr − ) corresponds to the case with λ = λ g and q cr = q cr2 . The birthrate isthe highest ( ∼ . × − yr − ) when λ = λ e and q cr = q cr1 . As shown in the left panel for 12 –the single star burst case, the evolution of the delay time marginally fits the observationallower limit. The second two panels of Figure 8 (models 5-7, in which α CE changes with q )show that there is no obvious change in both the evolution of the birthrate and the DTD,except that the peaks are just slightly lower and a little change happens in the evolution ofbirthrate when λ = λ b .The third two panels of Figure 8 (models 8-10) show the results with IMRD n ( q ) ∝ q . Inthe constant SFR case, the birthrates are higher than the standard model when λ = λ e , butcomparable or lower for λ = λ g and λ b . The range of the predicted birthrate for these modelsis between ∼ . × − yr − and ∼ . × − yr − . For the single star burst case, the delaytime evolution is marginally consistent with the observationally derived rate, considering anuncertainty of ∼ Z = 0 . λ = λ e . Toonen et al. (2012) estimated that the low metallicity has no effecton the delay-time evolution. However, our results show that the low metallicity changes thedelay-time evolution to some extent.In summary, our DD models generally under-predict the SN Ia rate, except for themodel with λ = λ e and q cr = q cr1 . The criteria of q = M /M ≥ / > . M (cid:12) implies that the mass of secondary WDs must be at least 0.55 M (cid:12) , so thesecondary He WDs donot significantly contribute to the SN Ia rate. This motivates us tofurther investigate a condition to increase the SN Ia rate under the DD scenario. Anotherkey parameter to describe the nature and outcome of the DD systems is the distribution ofthe initial orbital separation. Figure 9 shows the predicted SN Ia birthrate when we adoptthe range of the initial orbital separations as 3 ≤ a i /R (cid:12) ≤ (Hurley et al. 2002), insteadof 3 ≤ a i /R (cid:12) ≤ in our reference models. The predicted SN Ia rate increases, rangingbetween ∼ . × − and ∼ × − yr − depending on the BPS parameters. This brings alarge fraction of our BPS moldes to the values as high as observationally derived, under theparticular criteria we assumed for a DD system to explode as an SN Ia. We note that the results shown in the previous section does not necessarily provide afair comparison to the absolute values of the predicted SN Ia rates to the observed rate. Therate is highly dependent on the criteria for a DD system to become an SN Ia, where differentscenarios have different criteria (and frequently the criteria are not accurately determinedfrom the first principle). In the previous section, we adopted the same criteria as those 13 –adopted by Chen et al. (2012), so that we can focus on dependence of the SN Ia rate ondifferent BPS parameters by calibrating our BPS models with previous study. In this section,we study the SN Ia rate based on physically-motivated criteria, taking into account recentdevelopment in the explosion simulations.Given a lack of the detailed knowledge on the real explosion mechanisms, it is notpossible to cover all the possible models. However, investigating the following two scenariosprovides useful insight, as these models likely represent the lower and upper limits on SN Iarate within the DD scenario.1.
Carbon-Ignited Violent Merger Model:
When two CO WD merge, the high ac-cretion rate would create hot spots on the primary WD’s surface where the temperatureis so high that carbon burning would be ignited explosively and produce a detonationwave (Pakmor et al. 2013). While there is still a numerical convergence issue (see, e.g.,Tanikawa et al. 2015), the numerical simulations agree that this mode is likely a resultfor a merging DD system if a total mass well exceeds the Chandrasekhar limiting massand the mass ratio is nearly unity. Since this prediction is robust, this will give a lowerlimit for the SN Ia rate from the DD system. We adopt the following criteria: (1)0 . M (cid:12) < M < M (Sato et al. 2015) (the combined mass of two CO WDs is at least1 . M (cid:12) ), (2) q > . Chandrasekhar Mass Model:
If a DD system avoid a prompt detonation at themerging, the system is then represented by a massive CO WD that accretes materialsfrom a thick accretion torus and a hot envelope. Given the high accretion rate, theWD is suggested to become a ONeMg WD by a carbon deflagration and then thesystem would not explode as an SN Ia (Saio & Nomoto 1985). However, it could stillavoid the deflagration (Yoon et al. 2007), and in this case the primary WD is expectedto evolve into a Chandrasekhar mass WD. Thus, adopting the following criteria, weshould obtain an upper limit for the SN rate under the DD scenario: (1) the combinedmass of the two WD ≥ . M (cid:12) , and (2) the merging of two WDs within the Hubbletime.The upper two panels of Figure 10 show the case of the ‘Chandrasekhar mass model’.Depending on the BPS parameters, the range of SNe Ia birthrate is between ∼ . × − and ∼ . × − yr − . From these results, we can see that the mass ratio criterion tobecome SNe Ia has some effects on SNe rates but not that much in this context (i.e., eitherthe criterion is set at q ∼ / ∼ . M (cid:12) ,the systems should in any case have the mass ratio of q > .
5. 14 –However, this modest dependence on the critical mass ratio is not necessarily the caseif we consider the criterion much tighter than q ∼ .
5. The simulated SNe rates with models1-4 for the ‘Carbon-Ignited Violent Merger Model’ are shown in the lower two panels ofFigure 10. Compared to the results with the modest criterion ( q > / ∼ yeas,unlike what is observationally derived. This also results in a rapid rise of SN Ia rate inthe case of the constant SFR. With this criterion, a large fraction of systems having thecombined mass exceeding 1 . M (cid:12) are now rejected as SN Ia progenitors. The peak in theDTD corresponds to the peak in the distribution of the mass ratio at q > . ∼ . × − yr − , and those in models 2-4fall in the range between ∼ . × − and ∼ . × − yr − . The model with λ = λ e (model 2) and q cr = q cr1 gives the highest birthrate, and model 4 gives the lowest rate underthis tight criterion.
4. Discussion and conclusions
With the Monte Carlo method, we have performed detailed BPS simulations of theevolution of WD+MS binaries and double degenerate systems for 37 models with differentrecipes for the key binary evolution processes. The final systems in our simulations come frombinaries that have experienced CE evolution at least once, and the effects of the followingkey processes/conditions have been investigated: the CE efficiency ( α CE ), the binding energyparameter ( λ ), the critical mass ratio ( q cr ), and the initial mass ratio distribution, andmetallicity.de Kool & Ritter (1993) performed pioneering simulations for the formation of PCEBsby adopting λ = 0 .
5. Later, Dewi & Tauris (2000) demonstrated that the binding energyparameter λ changes with the evolutionary phases. The works on the PCEB simulationswere further updated (Willems & Kolb 2004). However, these previous results were differentfrom observations in some degree. Recently, a number of researchers performed BPS forPCEBs with different assumptions, and altogether these studies would lead to comprehensiveinvestigation of the formation and evolution of PCEBs. For example, Politano & Weiler(2007) used very low values for the CE efficiency that changes with the mass of the secondarystar. Davis et al. (2010, 2012) relaxed the assumption that the binding energy parameter isa constant, and considered both λ g and λ b to describe λ . Similar comprehensive studies havealso been made by Zorotovic et al. (2010) and Toonen & Nelemans (2013). It seems thateither the modeling assumptions did not treat the parameters in the realistic ways or themodel results were not comparable with observations. More recently, Camacho et al. (2014) 15 –presented detailed analysis of the selection effects that affect the sample of observed PCEBsobtained through SDSS. They also provided a thorough comparison between their BPSresults and the observed sample of these systems. However, their sample was very limited.Zorotovic et al. (2014) presented a systematic investigation that includes the contributionfrom the recombination energy to the energy budget of the CE evolution. They found thatthe recombination energy leads to a large number of PCEBs with long orbital period (longerthan 10 days), considering only three different input parameters in BPS simulations. In ourpresent work, we have further updated and extended the previous works by systematicallyconsidering different recipes for the key binary evolution processes (e.g., by using threedifferent prescriptions for λ , three different recipes for the critical mass ratios).In our BPS simulations for PCEBs, the main features that characterize the distributionsof resulting binary parameters for the different models can be summarized as follows:1. The three different prescriptions for the binding energy parameter ( λ ) influence thedistributions clearly. Binaries can eject the CE more easily and result in shorter orbits(for both PCEBs and double WD systems), for a higher value of λ (i.e., λ = λ e ). Theorbital period of PCEBs can be as short as 0.008 day when λ = λ e . If α CE is treatedto change with q , a larger number of systems survive at the CE phase than in the casewhere α CE = 1.2. The effect of different initial mass ratio distributions are mainly reflected by the re-sulting mass distributions.3. For binaries with a low metallicity, a larger number of systems can dispel the CE andtend to have shorter orbits. A choice of the prescription for the critical mass ratios turnout to have only little effect on the distributions of binary parameters of the resultingPCEBs.Double WD systems considered to be possible progenitors of SNe Ia are descendants ofWD/MS binaries. There have been a number of BPS works on the SD and DD scenarios andobservational constraints (Maoz et al. 2014; Ruiz-Lapuente 2014, for recent reviews). TheSNe Ia rate is one of the important constraints to understand the natures and progenitors ofSNe Ia. In our BPS study of the SN Ia rate, we tested different prescriptions for five maininput physical parameters in the binary evolution. By considering several conditions for thesystems to explode as an SN Ia, our findings are summarized as follows:1. A larger number of double WD systems can survive the CE phase when λ changeswith the evolution than a case where it is set as a constant. A number of double WD 16 –systems have mass ratios around 1. We have shown that α CE , λ , n ( q ) and Z affect thedistributions of the resulting double WD parameters.2. Considering a merging of CO WD with a CO WD or a He WD, the simulated SNe Iabirthrate ranges between ∼ . × − to ∼ . × − yr − with different parameters.With our first criteria, He WDs donot have significant contribution to the SN Ia rate.This range applies when the criterion for the mass ratio to become an SNe Ia is notlarger than q ∼ .
8. The simulated birthrates are marginally comparable with theobserved SNe Ia rate if λ = λ e and q = q cr1 .3. If we adopt the initial orbital separations of 3 ≤ a i /R (cid:12) ≤ instead of 3 ≤ a i /R (cid:12) ≤ , then we obtain the higher birthrate ranging between ∼ . × − and ∼ × − yr − (see the Figure 9). The observed SNe Ia rate is well explained by this modelif λ = λ e .To summarize, we conclude that the CE parameters, the metallicity and other parame-ters affect the evolution of SNe Ia rate with time. While even our optimistic models still showa discrepancy in the evolution of the SN Ia rate as compared to the observational inferredpower law distribution (i.e., the ratio of the ‘young’ population to the ‘old population’ isquite high in our BPS models), the simulated SNe Ia rate can be comparable to the obser-vations, depending on the treatment of the parameters. Especially, a combination of λ = λ e and q cr = q cr1 results in the highest SNe Ia rate as being compatible to the observations.While the above argument should be correct for the dependence of the resulting SNe Iarate to the BPS parameters, the absolute rate should depend on the particular criteria for aDD system to explode as an SNe Ia. The criteria are different for different explosion models.To see the effect and to link the present BPS study to the existing explosion models, we alsocalculate the simulated SNe Ia rate under different assumptions on the criterion of the massratio for the SN Ia progenitors.1. Even if we totally remove the criterion of the mass ratio as compared to our referencemodel (where q > / ∼ . × − and ∼ . × − yr − .This prescription corresponds to the classical Chandrasekhar mass scenario.2. The high value of the mass ratio is an essential ingredient in the so-called violentmerger scenario. Setting the critical mass ratio higher than ∼ .
8, the resulting SNeIa rate is decreased substantially. By adopting 0 . M (cid:12) < M < M (Sato et al. 2015)and q > .
8, the SN Ia birthrates are ∼ . × − yr − (model 1), and ∼ . × − – ∼ . × − yr − (models 2–4 with different prescriptions for q cr ). Even for the 17 –optimistic model with λ = λ e (model 2) and q cr = q cr1 , the simulated SNe Ia rate is farbelow the observationally derived SNe Ia rate.This work was funded by China Postdoctoral Science Foundation, the Natural ScienceFoundation of China under grant numbers 11390371, 11133001 & 11333004, the StrategicPriority Research Program of CAS under grant No. XDB09000000. The work also waspartly supported by JSPS KAKENHI (No. 26800100) from MEXT and by WPI Initiative,MEXT, Japan. REFERENCES
Ablimit, I., & Li, X.-D. 2015, ApJ, 800, 98Ablimit, I., Xu, X.-J., & Li, X.-D. 2014, ApJ, 780, 80Badenes, C., & Maoz, D. 2012, ApJ, 749, L11Bours, M. C. P., Toonen, S., & Nelemans, G. 2013, A&A, 552, A24Camacho, J., Torres, S., Garc´ıa-Berro, E., et al. 2014, A&A, 566, A86Chen, X., Jeffery, C. S., Zhang, X., & Han, Z. 2012, ApJ, 755, L9Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J. & Verbunt, F.W. M. 2014, A&A, 563,A83Dan, M., Rosswog, S., Brggen, M., & Podsiadlowski, P. 2014, MNRAS, 438, 14Davis, P. J., Kolb, U., Willems, B. & G¨ansicke, B. T., 2008, MNRAS, 389, 1563Davis, P. J., Kolb, U., & Willems, B. 2010, MNRAS, 403, 179Davis, P. J., Kolb, U., & Knigge, C. 2012, MNRAS, 419, 287de Kool, M. & Ritter, H. 1993, A&A, 267, 397de Marco, O., Passy, J., Moe, M., et al. 2011, MNRAS, 411, 2277de Mink, S. E., Pols, O. R., & Hilditch, R. W. 2007, A&A, 467, 1181Dewi, J. D. M. & Tauris, T. M. 2000, A&A, 360, 1043 18 –Garc´ıa-Berro, E., Soker, N. & Althaus, L.G. 2015, arXiv:1503.01739v1Han, Z., Podsiadlowski, P., & Eggleton, P. P. 1994, MNRAS, 270, 121Hillebrandt, W., & Niemeyer, J. C. 2000, ARA&A, 38, 191Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897Iben, Jr.,I. & Livio, M. 1993, PASP, 105, 1373Kaplan, L.D. 2010, ApJ, 717, L108Kiel, P. D., & Hurley, J. R. 2006, MNRAS, 369, 1152Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545Landau, L. D., & Lifshitz, E. M. 1971, Classical Theory of Fields (Oxford: Pergamon)Li, X.-D., & van den Heuvel, E. P. J. 1997, A&A, 322, L9Maoz, D., Mannucci, F., Li, W., et al. 2011, MNRAS, 412, 1508Maoz, D., Mannucci, F., & Brandt, T. D. 2012, MNRAS, 426, 3282Maoz, D., & Mannucci, F. 2012, PASA, 29, 477Maoz, D., Mannucci, F., & Nelemans, G. 2014, ARA&A, 52, 107Marsh, T. R. 2011, Class. Quant. Grav., 28, 094019Mennekens, N., Vanbeveren, D., De Greve, J. P., & De Donder, E. 2010, A&A, 515, A89Napiwotzki, R., Karl, C.A., Nelemans, G., et al. 2007, in 15th European Workshop on WhiteDwarfs, ASP Conference Series, 372, ed. R. Napiwotzki & M.R. Burleigh (Astron. Soc.Pacific, San Francisco), p. 387Nebot Gmez-Morn, A., A., G¨asicke, B. T., Schreiber, M. R., et al. 2011, A&A, 536, A43Nelemans, G. & Tout, C. A. 2005, MNRAS, 356, 753Nomoto, K., & Iben, Jr., I. 1985,ApJ, 297, 531Paczy´nski, B. 1976, in Structure and Evolution of Close Binary Systems, eds. P.Eggleton, S.Mitton, & J. Whelan (Dordrecht: Kluwer), IAU Symp., 73, 75Pakmor, R., Kromer, M., Taubenberger, S., & Springel, V. 2013, ApJ, 770, L8 19 –Parsons, S. G., G¨asicke, B. T., Marsh, T. R., et al. 2013, MNRAS, 429, 256Piersanti, L., Gagliardi, S., Iben, Jr., I., & Tornamb´e, A. 2003, ApJ, 598, 1229Politano, M. & Weiler, K. P. 2006,ApJ, 641, L137Politano, M. & Weiler, K. P. 2007, ApJ, 665, 663Rebassa-Mansergas, A., Nebot G´omez-Mor´an, A., Schreiber, M. R., et al. 2012, MNRAS,419, 806Ruiter, A. J., Belczynski, K., & Fryer, C. 2009, ApJ, 699, 2026Ruiter, A.J., Belczynski, K., Sim, S.A., Hillebrandt, W., Fryer, C.L., Fink, M., & Kromer,M. 2011, MNRAS, 417, 408Ruiz-Lapuente., P. 2014, New Astron. Rev., 62, 15Saio, H. & Nomoto, K. 1985, A&A, 150, L21Sato, Y., Nakasato, N., Tanikawa, A., et al. 2015, ApJ, 807, 105Shao, Y., & Li, X.-D. 2014, ApJ, 796, 37Santander-Garc´ıa, M., Rodriguez-Gil, P., Corradi, R. L. M., Jones, D., Miszalski, B., Boffin,H. M. J., Rubio-Diez M. M., & Kotze, M. M. 2015, Nature, 519, 63Tanikawa, A., Nakasato, N., Sato, Y., Nomoto, K., Maeda, K. & Hachisu, I. 2015, ApJ, 807,40Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70Toonen, S. & Nelemans, G. 2013, A&A, 557, A87Wang, B., & Han, Z. 2012, New Astron. Rev., 56, 122Wang, B., Ma, X., Liu, D.-D. et al., 2015, A&A, 576, A86Wang, C., Jia, K. & Li, X.-D. 2016a, MNRAS, 457, 1015Wang, C., Jia, K. & Li, X.-D. 2016b, RAA, submittedWebbink, R. F. 1984, ApJ, 277, 355Webbink R. F., 2008, in Milone E. F., Leahy D. A., Hobill D., eds, Short-Period BinaryStars: Observations, Analyses, and Results. Springer,Berlin, p. 233 20 –Willems, B., & Kolb, U. 2004, A&A, 419, 1057Xu, X.-J., & Li, X.-D. 2010, ApJ, 716, 114Yoon, S.-C., Podsiadlowski, Ph. & Rosswog, S. 2007, MNRAS, 380, 933Zorotovic, M., Schreiber, M. R., G¨asicke, B. T., & Nebot G´omez-Mor´an, A. 2010, A&A,520, A86Zorotovic, M., Schreiber, M. R. & G¨asicke, B. T. 2011, A&A, 536, A42Zorotovic, M., Schreiber, M. R., Garc´ıa-Berro, E., et al. 2014, A&A, 568, A68
This preprint was prepared with the AAS L A TEX macros v5.2.
21 –Table 1: Different models used in our calculationModel α CE λ n ( q ) Z1 1 1 1 0.022 1 λ e λ b λ g λ e λ b λ g λ e ∝ q λ b ∝ q λ g ∝ q λ e λ b λ g λ , which specifies a model) and the critical mass ratios ( q cr for which the three cases areconsidered for each model), there are three lines (except for the standard model) shown ineach panel. The models are indicated as follows: The solid red, green and blue lines for λ = λ e with q cr1 , q cr2 and q cr3 , respectively. The dotted cyan, magenta and yellow lines arefor λ = λ b with q cr1 , q cr2 and q cr3 , respectively. The dash dotted orange, green+yellow andgreen+cyan lines for λ = λ g with q cr1 , q cr2 , and q cr3 , respectively. The sold black line is forthe selected observed samples. 23 –Fig. 2.— Distribution of the secondary MS masses. See the caption of Figure 1 for the modeldescription. 24 –Fig. 3.— Distribution of the WD masses. See the caption of Figure 1 for the model descrip-tion. 25 –Fig. 4.— The left(1st) panel shows the distribution of the orbital period and the WD masse.The middle(2nd) panel shows the distribution of the orbital period and secondary mass.The right (3rd) panel shows the distribution of the secondary mass and the WD mass all atthe PCEB phase. In each panel, 10 different models are shown, including model 1 (the firstpanel), model 4 with q cr1 , , (the 2nd–4th panels), model 3 with q cr1 , , (the 5th–7th panels)and model 2 with q cr1 , , (the 8th–10th panels), respectively. Here, the panel numbers arecoordinated as follows: 1st (the first row, left), 2nd (the first row, right), 3rd (the secondrow, left), 4th (the second row, right), and so forth. 26 –Fig. 5.— Distributoin of the orbital periods of double WD systems. See the caption ofFigure 1 for the model description. he sold black line is for the selected observed samples(for the details see the text). 27 –Fig. 6.— Distribution of the primary CO WD masses. See the caption of Figure 1 for themodel description. 28 –Fig. 7.— Distribution of the mass ratios between secondary WDs and primary CO WDs.See the caption of Figure 1 for the model description. 29 –Fig. 8.— The SN Ia birth rates calculated for various models, for a single starburst of10 M (cid:12) (left panels) and for a constant star formation rate of 5 M (cid:12) yr − over the past 13.7Gyr(right panels). The letter represents the Milky Way. For the panels(models) & lines(sameas Fig. 1), and for the three sold black lines of observational results, see the text. 30 –Fig. 9.— In this figure, the range of the initial orbital separation is set to be 3 ≤ a i /R (cid:12) ≤ ,models are the same as in Figure 8. 31 –Fig. 10.— This Figure shows the results of model 1 and models 2–4 with three prescriptionsfor q cr , under our different criteria for a DD system to explode as an SN Ia. The upperpanels show the results from the ‘Chandrasekhar Mass Scenario’ and the lower panels showthe results from the ‘Carbon-Ignited Violent Merger Scenario’ (see §§