DCB & GW

# Monte Carlo population synthesis on massive star binaries: Astrophysical implications for gravitational wave sources

Iminhaji Ablimit11affiliation: Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502; iminhaji@kusastro.kyoto-u.ac.jp and Keiichi Maeda11affiliation: Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502; iminhaji@kusastro.kyoto-u.ac.jp
###### Abstract

There are important but unresolved processes in the standard formation scenarios of double compact star binaries (DCBs; BH-BH, BH-NS, NS-NS systems), such as the mass transfer and common envelope (CE). We analyze the effects of different assumptions on key physical processes and binary initial conditions on the formation of DCBs, and study the massive star binary evolution toward the GW150914-like events and GW151226-like events, including a survey of proposed prescriptions for the mass transfer and the binding energy parameter () in the CE phase, which have not been tested previously in binary population synthesis (BPS) studies. We find that different prescriptions for the mass transfer and have moderate influence on the merger rates and properties of BH-BH systems, but significantly affect those of DCBs containing a NS. The small initial separation enhances the merger rate of BH-BH systems. We confirm that a low metallicity (Population II) is a key to producing the GW150914-like BH-BH binaries, and we find that the predicted rate can be compatible to the inferred rate if an optimized combination in the initial conditions is adopted. On the other hand, the GW151226-like systems are naturally expected at the solar metallicity (Population I) within the standard CE evolution. If these two systems are to be explained this way through the CE evolution as results of Pop II and Pop I systems, respectively, then the transition in the mode of the evolution must occur at to compare to the observed relative frequency of these two types of events. On the other hand, the indirectly indicated NS-NS merger rate from observed NS-NS systems (Kim et al. 2015) can be reproduced by the BPS simulations only by the conditions inconsistent with those to explain the BH-BH merger rate through the isolated binary evolution scenario. Therefore, we predict that the NS-NS merger rate is indeed higher than so far inferred by the indirect method, or otherwise the BH-BH GW systems should be formed through a different evolutionarily pathway than NS-NS systems.

close – stars: evolution – stars: black hole – gravitational waves – X-rays: binaries

## 1 Introduction

Double compact star binaries (DCBs), such as Black hole-Black hole (BH-BH), BH-Neutron star (BH-NS) and NS-NS binaries, play an important role in binary population synthesis (BPS) studies. They provide excellent laboratories to test physical parameters in binary formation and evolution (Dominik et al. 2012). They are also expected to be potential progenitors of various astrophysical objects, such as ultra-compact X-ray binaries (van der Sluys et al. 2005) and short/long -ray bursts (Eichler et al. 1989; Narayan et al. 1992). Furthermore, the inspiral and merger of DCBs have been promising targets as Gravitational wave (GW) emitters within reach by modern GW detectors (e.g. Advanced LIGO, Advanced Virgo, and KAGRA) (e.g. Kulkarni et al. 1993; Morscher et al. 2013; Tanikawa 2013; Postnov & Yungelson 2014; Belczynski et al. 2015; O’Leary et al. 2016). Observations of NS-NS systems such as the binary pulsar PSR 1913+16 demonstrate that the merging time can be shorter than the Hubble time, making the NS-NS mergers promising detectable GW emitters (Phinney 1991). In the year 2015, gravitational wave sources GW150914 and GW151226 were finally detected by Advanced LIGO (Abbott et al. 2016a,e). More recently, another massive binary black hole merger GW170104 is also reported by Advanced LIGO (Abbott et al. 2017a). On August 14, 2017, the Advanced Virgo detector and the two Advanced LIGO detectors coherently observed GWs from the coalescence of two black holes, GW170814 (Abbott et al. 2017b). These detections gave a crucial proof for Albert Einstein’s theory of general relativity (Einstein 1918). Somewhat surprisingly, these sources are not the expected NS-NS mergers, but the BH-BH mergers. These detections confirm the existence of BHs, and challenge our understanding of the stellar evolution. Any stellar evolution model must account for the existence of the merging BH-BH binaries within the Hubble time. It opened a new window for studying DCBs.

Following this discovery, subsequent works revisited the merger rate of BH-BH binaries (Abbott et al. 2016c) and studied possible evolution scenarios especially for GW150914. Belczynski et al. (2016) proposed the standard (isolated) binary evolution scenario including the common envelope (CE) phase at a low metallicity for GW150914. After this, Woosley (2016) introduced two possible scenarios for producing GW150914. One is based on the chemically homogeneous evolution of a rapidly rotating single star of 150 and the other on the delayed merger of two black holes made by 70 and 90 stars in a binary system (initial masses and evolution path different from Belczynski et al. 2016). On the other hand, Tagawa et al. (2016) presented post-Newtonian N-body simulations of mergers of gas-accreting stellar-mass black holes (BHs) within dense stellar environment, and claimed that the BH-BH merger GW150914 was likely driven by three-body encounters accompanied by a few of gas accretion. GW151226 differs from GW150914 in the significantly lower inferred BH masses (Abbott et al. 2016e). BHs with such masses can be formed at solar metallicity (Spera et al. 2015). As we see from these possibilities, we still do not know exact evolutionary path ways to the compact star mergers (e.g., the merging BHs like GW150914).

In studying the natures of DCBs, obtaining reliable predictions of their merger rates and birth rates is still a challenge (Abadie et al. 2010). For example, the mass transfer (and/or mass loss) process, the CE phase and some other parameters in the binary evolution are critically important but the details are not well known. These can be tested by accumulating samples of GW detections from BH-BH systems as well as expected detections of BH-NS and NS-NS systems. Advanced LIGO, Advanced Virgo and KAGRA (Harry et al. 2010; Sengupta et al. 2010; Somiya 2012) will probe the universe in search for DCB signatures. To be ready for such detections from a theoretical point view, we still need to continue to present more information on the formation and evolution of DCBs. In the past decades, a series of Binary Population Synthesis (BPS) works on DCBs have been performed to investigate the effects of various physical parameters on the formation, evolution and merger rate of DCBs, but there are still debates on the treatment of uncertain physical processes.

The mass transfer process and the CE phase are crucial for producing all kinds of compact star binaries, while the details of these processes are still unknown. Massive MS/MS close binaries (where MS means a main-sequence star, with a initial mass ) generate close compact star binaries mainly by going through a CE phase at least once. The MS/MS binary may involve the CE phase if the initial orbital separation is not large enough and the more massive one fills the Roche lobe (RLOF) due to the expansion of the star during rapid nuclear evolution. The first mass transfer may undergo on the dynamical timescale, leading to the CE evolution where the extended envelope engulfs both stars. If the stars survive the CE phase, then the binary may again evolve into the CE evolution including a BH/NS, if the second RLOF mass transfer proceeds on the dynamical timescale. The system that survived the CE phase would finally become a DCB. The MS/MS binaries in a wider orbit may also evolve into the CE evolution after the massive one evolved to a BH/NS if the second RLOF mass transfer is unstable, and the DCBs would be produced after the successful CE ejection (Iben & Livio 1993; Ivanova et al.2013). As the orbital energy and angular momentum are removed by the CE ejection (Paczyński 1976), the binary separation is decreased after the CE phase. The yet-unresolved CE phase is a too-short lived physical process to be detected, and there are some papers which address indirect observational evidence for it (e.g., Sion et al. 2012; Ivanova et al. 2013).

Theoretical studies still contain uncertainties. However, estimating the merger rates and providing various properties of merging binaries are very useful for the future detections. The theoretical predictions are crucial to understand astrophysical natures of the progenitor systems of detected GW sources (Stevenson et al. 2015). The related estimates and studies have been presented by different groups in the past decade, with different codes and different sets of physical assumptions (e.g. Lipunov et al. 1997; Nelemans et al. 2001; Voss & Tauris 2003; Dewi & Pols 2003; Yungelson et al. 2006; Belczynski et al. 2008 Mennekens & Vanbeveren 2014; Eldridge & Stanway 2016). Nearly all have drawn an attention to serious difficulties in obtaining this reconciliation, underlining our incomplete knowledge about some of these critically important evolutionary phases. It is thus always useful to add and compare the results of different approaches to the problem. Recently, de Mink & Belczynski (2015) performed a comparative BPS study on the evolutionary path of massive stars by adopting a similar approach as Dominik et al. (2012), varying some important initial inputs (such as initial mass function, distribution of mass ratio, distribution of orbital periods) based on observation results of Sana et al. (2012). They further provide the merger rates and distributions of properties for BH-BH, BH-NS and NS-NS systems by adopting the binding energy parameter from Xu & Li (2010) for the CE and assuming that half of the mass is accreted but the other half is ejected from the system for the mass transfer prescription. In this paper, we further expand the exploration on how different treatment of these processes would affect the outcome. We use three proposed expressions for the binding energy parameter in the CE model (, , and , the values derived from the MESA, see §2 for details, and Wang et al. 2016a,b; Ablimit et al. 2016), and adopt three detailed mass transfer models (three different , for more details see Section 2).

In this paper, we investigate influences that a series of binary evolution prescriptions (, , , eccentricity, initial separation, kick velocity distribution and other parameters) have on the formation of DCBs (see Table 1). We present the properties of the possible gravitational wave signals from the mergers of those DCBs in our BPS study. In §2, we describe our BPS code and our models to treat the binary physical processes including the CE evolution and mass transfer assumptions. The properties of DCBs produced by our models are presented, and the properties of resulting gravitational wave signals and the implication for GW150914 and GW151226 are discussed in §3. The paper is closed in §4 with conclusions and discussion.

## 2 Monte Carlo binary population synthesis

Our BPS code is the one developed by Hurley et al. (2002) and modified by Kiel & Hurley (2006). As compared to the version used by Kiel & Hurley (2006), we have updated and modified the code in several aspects, especially regarding the conditions for mass transfer and the treatments of CE evolution, which are briefly described below. For the initial conditions, we consider the constraints from previous observational results, as well as we use new constraints on initial conditions (e.g., high binary fraction(), preferentially short orbital periods and uniform mass ratio distribution) for massive O stars provided by Sana et al.(2012). We adopt the initial mass function (IMF) of Kroupa et al. (1993) for the primary mass distribution,

 f(M1)=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩0M1/M⊙<0.10.29056(M1/M⊙)−1.30.1≤M1/M⊙<0.50.1557(M1/M⊙)−2.20.5≤M1/M⊙<1.00.1557(M1/M⊙)−α1.0≤M1/M⊙≤150% , (1)

with in this study. The secondary mass distribution is determined by the distribution of the initial mass ratio,

 n(q)={0q>1μqν0≤q<1, (2)

where , is the normalization factor for the assumed power law distribution with the index . We consider a flat distribution ( and constant) for the initial mass ratio distribution (IMRD). The distribution of the initial orbital separation, , is assumed to be given by the following formalism (Davis et al. 2008),

 n(ai)={0ai/R⊙<3 or ai/R⊙>1060.078636(ai/R⊙)−13≤ai/R⊙≤106 . (3)

We test two possibilities for the orbit eccentricity. In one case, the uniform (flat) initial eccentricity distribution is assumed in a range between 0 to 1. We also test the possibility that binaries are in circle orbits or have negligible eccentricities (Dominik et al. 2015; Nishizawa et al. 2016) with tighter orbital distribution (Hurley et al. 2002). This is motivated by observational results of Sana et al.(2012). We assume the binary fraction of (). We fix the metallicity to be for the parameter study for the key binary evolution processes, while we consider the low metallicity case of in discussing the origins of GW150914 and GW151226 (§3.3).

Regarding the key physical processes, our simulations have three main tunable parameters, i.e., , and . One main aim of this paper is to investigate the effects that these prescriptions have on the evolution toward compact star binaries. The CE evolution is an important but unsolved phase in the binary evolution process. The widely used -model considers energy conservation (Webbink 1984),

 Ebind=αCEΔEorb , (4)

where , , and are the binding energy of the envelope, the efficiency parameter and the change in the orbital energy during the CE phase, respectively.

In the -formalism, the efficiency parameter () and the binding energy parameter () are vital to the fate of the CE evolution, but they have been treated as constants in most previous works. However, in reality, they should change with the mass and evolutionary stage of the primary star. The binding energy of the envelope is expressed by the following:

 Ebind=−GM1MenλR1, (5)

where , and are the total mass, envelope mass and radius of the primary star, respectively. The two stars coalesce if the radius of either of the stars satisfies Equation(75) of Hurley et al.(2002) and if the CE evolution is longer than the dynamical timescale (for details of the criterion for surviving or merging during the CE phase, see Hurley et al.2002).

The values of for evolving stars can be calculated by considering gravitational energy only (hereafter ), adding inner energy (), or adding entropy of the envelope (). The binding energy can be described as follows if it is only due to the gravitational binding energy,

 Ebind=∫M1Mc−GM(r)rdm , (6)

where is the gravitational constant, is the primary mass and is its core mass. If both of the thermal energy and the recombination energy contribute to the binding energy, then

 Ebind=∫M1Mc[−GM(r)r+U]dm , (7)

where is the internal energy (Han et al. 1994). If the enthalpy model (Ivanova & Chaichenets 2011) is used, then

 Ebind=−∫M1Mc[−GM(r)r+U+Pρ]dm , (8)

where and are the pressure and the density of the gas, respectively. For more details, see Ivanova & Chaichenets(2011) and Wang et al.(2016a). We test all the three prescriptions independently in our BPS study, adopting the calculation results with the MESA code by Wang et al. (2016a,b). This is a new treatment, comparing to previous BPS studies.

Most BPS studies take merely as a constant (typically 1 or 3, see Hurley et al. 2002 and Dominik et al. 2012). For relatively low mass stars, Davis et al. (2012) derived a relation of with from analysis of observed post-CE white dwarf (WD) binaries. On the other hand, a relation between and for massive stars leading to a NS/BH DCBs is not clarified yet. As a demonstration, in this paper we test the description by Davis et al. (2012) for . While it would not provide a physically motivated relation for the massive stars, it would give an estimate about how the results are dependent on variations of . In typical situation, the most values of calculated from this relation are lower than 1 (see Davis et al. 2012 and Ablimit et al. 2016), similar to a constant value adopted for massive binary evolution. We adopt the following formula from Davis et al. (2012):

 log10αCE=ϵ0+ϵ1log10(q) , (9)

where and are constants taken from Davis et al. (2012). We simulate our BPS by treating the CE phase with this equation or with for the comparison with other theoretical assumptions.

The critical mass ratio () is another physical key parameter that determines the stability of the mass transfer. Shao & Li (2012) computed the critical mass ratio considering the possible response of the accreting star (i.e., spin-up and rejuvenation) under three different assumptions (including conservative and non-conservative cases): (1) Half of the transferred mass is accreted by the secondary, and the other half is lost from the system, also removing the specific orbital angular momentum of the accretor (see also de Mink et al. 2007). (2) The transferred mass is assumed to be accreted by the secondary unless its thermal timescale () is much shorter than the mass transfer timescale (). The accretion rate is limited by– (Hurley et al. 2002). Rapid mass accretion may drive the accretor out of thermal equilibrium, which will expand and become over-luminous. Shao & Li (2012) found the values of are usually much lower than that of the same star in thermal equilibrium. Therefore, it is always as , and the mass transfer is generally conservative. (3) The accretion rate onto a rotating star is reduced by a factor of (), where is the angular velocity of the star and is its critical value. In this prescription, a star cannot accrete mass when it rotates at . The remaining material is ejected out of the binary by the wind, and it takes away the specific orbital angular momentum from the accretor. The critical mass ratios corresponding to these three cases are denoted as , and , respectively. For each model as described above, we run the BPS simulations with these different treatment of the critical mass ratio (i.e., adopting , , or ). The three mass transfer models are also a major update in this BPS code.

For , we track (the angular rotational velocity of a star) in the same way as given in Hurley et al. (2000). We adopt a uniform distribution of the rotational velocity within a star (van den Heuvel 1968). The spin evolution of each component in a close binary system could be influenced by tidal synchronization, which is taken into account in our calculations. For the evolution of angular velocity of each star, including the tidal evolution of a binary, we follow default prescriptions in the BPS code of Hurley et al.(2002).

Another main issue is the evolution of ‘single’ massive stars and the birth of the compact stars. For massive stars, stellar winds have non-negligible influence on final outcomes (Belczynski et al. 2010), but their loss rates have not been accurately determined. For example, Hurley et al. (2000) (wind 1) and Vink et al. (2001) (wind 2) gave different prescriptions for O and B stars in different stages. Another main difference between the default model wind 1 and the wind 2 model exists for treatment of luminous blue variable stars (LBVs). In the wind 1 for LBVs with and :

 ˙M=0.1×[10−5(R/R⊙)(L/L⊙)0.5−3]3((L/L⊙)6×105−1)M⊙yr−1, (10)

where and are the radius and luminosity of the star, respectively. In the wind 2, it is given as,

 ˙M=flbv×10−4M⊙yr−1, (11)

with the standard factor (for more details see Humphreys & Davidson 1994). Here, we adopt two prescriptions for the wind loss rate: the wind 1 model based on Hurley et al. (2000,2002) and wind 2 model based on Vink et al. (2001) and Humphreys & Davidson (1994). For low mass H-rich stars, it is taken from Hurley et al. (2000).

Supernova (SN) explosions play an important role on the birth of BH and NS. We consider the prescription which successfully reproduces the observed mass distribution of BH and NS star. This is based on the neutrino-driven convection-enhanced SN mechanism (the so-called rapid SN mechanism). This prescription results in a ‘failed’ SN for a massive star with the initial mass giving the birth of BH. The gravitational remnant mass ( in ) under the rapid SN mechanism is described by Fryer et al. (2012) as follows,

 Mrem=0.9(Mproto+Mfb) . (12)

Here is the mass of the proto-compact object and

 Mfb=⎧⎪⎨⎪⎩0.2M⊙MCO<2.5M⊙0.286MCO−0.514M⊙2.5M⊙

with

 ffb=⎧⎪⎨⎪⎩1.06M⊙≤MCO<7M⊙a1(MCO/M⊙)+b17M⊙≤MCO<11M⊙1.0MCO≥11M⊙, (14)

where is the CO core mass, here is the mass of the progenitor star mass just before the core-collapse, (see Wang et al. 2016a) and .

The formation of the BH is connected to the compactness of the stellar core at the time of collapse in this prescription. Recent study on the SN explosion mechanism suggests that the compactness parameter is an important parameter, which is defined as a mass within a given radius in the core of the progenitor star (O’Connor & Ott 2011). The stars with a small compactness parameter are more likely to explode as SNe and produce NSs, while stars with a large compactness parameter tend to evolve to failed SNe that produce BHs. The He core masses or CO core masses prior to core collapse determine the compactness parameter, thus the newly born BH masses (e.g. Kochanek 2014, 2015; Sukhbold & Woosley 2014; Clausen et al. 2015). In our BPS simulation, the He and CO core masses are tracked following Hurley et al. (2000). We require that . We do not consider the formation of BHs by an accretion-induced collapse (AIC) of a NS in X-Ray Binaries while the formation of NS by the AIC of the WD is taking into account (Nomoto & Kondo 1991; Ablimit & Li 2015). Default values of Hurley et al.(2002) are used for the upper limit masses of NS and WD.

We also assume that a natal kick is imparted on the newly born BHs, similar to the case of the NS formation. The kick velocity is set to be inversely proportional to the remnant mass, i.e., , where is the kick velocity for NSs, which follows the Maxwellian distribution as

 P(υk)=√2πυk2σk3e−υk2/2σk2 . (15)

We use the velocity dispersion of (Hansen & Phinney 1997) or (Hobbs et al. 2005).

In Table 1, we summarize our models, where the ‘case’ distinguishes different combinations of and , while ‘model’ variants are for the other parameters (see Table 1). In our standard model, we assume , for the CE phase, , a flat IMRD for the initial conditions. For we adopt the default prescription in the BSE code by Hurley et al. (2002) in the standard run. In our standard run, the Belczynski et al.(2002)’s prescription for the SN mechanism is used. For the other model parameters in the standard run, we adopt the values used in model 3. Note that the tidal evolution is on in all the simulations run. Combining models 1–3 with cases 1–9, there are thus 27 different models. Including models 4–6 with cases 1–3 and the standard model, we cover 37 different models in our numerical calculations. The default values from Hurley et al.(2002) are adopted for other physical parameters which are not mentioned in this study.

We present our merger rates as follows. We assume a constant star formation rate (SFR) of over the past 10 Gyr. This is a simplified but frequently adopted treatment in the field of BPS study given that the star formation history of our galaxy is not well known (Wyse 2009). To allow straightforward comparisons to previous studies, we have decided to adopt this simple model for the SFR (see Dominik et al. 2012). We however take into account the evolution of SFR in §3.3, in the discussing the natures of GW150914 and GW151226-like sources.

## 3 Results

We perform the binary population synthesis simulations to see the effect of different combinations of the CE parameters, the critical mass ratio models and other parameters on the property distributions of the final BH-BH, BH-NS and NS-NS systems. We use the BPS code to generate initial MS/MS binaries under different models with different cases as described in Table 1 for obtaining these DCBs. The metallicity is fixed to be in §3.1-3.2, while a low metallicity is also considered in §3.3-3.4. Our Monte Carlo simulations have statistical fluctuations due to their finite size. The relative statistical error (, where is the number of simulated DCBs who can merge in 10 Gyr) is estimated to be 5% - 35% based on the results of model 1–6.

### 3.1 Merger rates

In Table 2 we list the double compact binary galactic merger rates calculated for a fiducial Milky Way-like galaxy for a single metallicity () and different variations of models. This is visualized in Figure 1 (which also includes the low metallicity case to be discussed in §3.3). The merger rates of BH-BH, BH-NS and NS-NS systems in our standard run are 15.1 , 1.21 and 40.3 , respectively. It is seen that influences obviously the merger rates while other parameter variations on models 1–3 do not have any clear effect. With our model 4/case 1, we have the highest merger rate of BH-BH systems (27.5 ), which means that combining the assumptions of the tighter initial separation, eccentricity with 0, varying , and is an optimistic setup for producing the larger number of BH-BH systems in a solar metallicity condition. We also find that the other initial conditions have effects within a factor of a few. The highest rates of BH-NS and NS-NS systems (13.8 and 153 ) are produced at by model 4 with constant , and . Namely, the prescriptions of and affect the rates of all the DCB merger rates in a similar way, and the treatment of affects the merger rates of the BH-BH and the other systems (BH-NS and NS-NS) in the opposite way.

The initial eccentricity distribution does not have clear effect on merger rates. If the distribution of orbital angular momentum is used to determine the initial state of each binary, rather than one of semi-major axis or period, the results do not depend on the form of any chosen eccentricity distribution. In fact, eccentricity need not be a free parameter. This becomes evident when considering that tidal interaction conserves angular momentum and that almost all systems circularize before RLOF.

It is seen that the CE, wind mass loss, SN mechanism and mass transfer are very crucial for the massive star binary evolution. Especially, a tight initial orbital period distribution and increase the merger rates of all kinds of DCBs, and the mass transfer prescription ( ) is crucial to determine the relative rates of BH-BH and NS-NS mergers. By taking the spin of the accretor into account, the mass accretion rate is reduced, enhancing the NS survival/formation. However, this may let the system evolve into a CE easily due to the angular momentum loss by the spin-enhanced mass loss from the system. At the same time, the prescription helps to enhance the envelope ejection and thus increases the chance that DCBs survive the CE. The expected rate of merging DCBs is maximized for the initial orbital separations of instead of , because the tighter initial orbital period distribution could bring more systems closer and making them more likely to interact. Indeed, our BPS models tend to predict high NS-NS merger rates, which will be discussed further in §3.4.

### 3.2 Distributions of final properties

Figures 2 shows the orbital period distributions of BH-BH, BH-NS and NS-NS systems (upper, middle and lower panels of the Figures 2), respectively. From the results of models 1, 3 and 6, we see that the orbital period is mostly affected by the kick velocity. The model with produces a larger number of very short orbital period DCBs than the model with . A larger number of BH-NS and NS-NS systems in the short orbits are produced with and also with . From the results, it is seen that the eccentricity distribution influences the orbital period distributions. The two prescriptions of and other initial conditions do not have a significant effect.

The chirp mass distributions of BH-BH, BH-NS and NS-NS systems are given in Figures 3, 5 and 6, respectively. The chirp mass is given as follows:

 Mchirp=μ3/5M2/5 , (16)

where and (with and the masses of the primary and secondary (in ), respectively). The wind mass loss is the most important function to determine the chirp mass of BH-BH systems. With the wind 1 model, the chirp mass of BH-BH systems ranges only from 5 to 15 , while it ranges from 4 to 20 with wind 2 model (see Figure 3).

The upper limits of the chirp mass in BH-BH systems in our calculation (15 and 20 with wind 1 and 2 models) are larger than the results of Belczynski et al.(2002, 2010) ( 11 and 15 with wind 1 and 2 models). In order to clarify the reason which causes this difference in the chirp mass distribution, we additionally test MS-MS binaries and single stars evolution (SSE). For the SSE simulation, we use the SSE code developed by Hurley et al.(2000). We compare the two prescriptions for the SN mechanism (the fate of the core collapse), one by Fryer et al.(2012) (F12, adopted in this paper) and the other one used by Belczynski et al.(2002, hereafter B02). We adopt and for the CE phase, the wind 2 model and solar metallicity in this exercise to compute the BH mass distributions. The other initial conditions are same as default ones in the code (for more details see Hurley et al. 2000, 2002). The calculated remnant mass from the single star evolution shows that the star which has a initial mass between 22 and 38 could leave a more massive remnant under the rapid SN explosion mechanism of F12 than that of B02 (see the lower panel of Figure 4). This mass range dominates the formation of the BH-BH systems following the IMF. The mass distributions of primary BHs in the BH-BH binaries under two SN mechanisms are given in the upper and middle panels of Figure 4. We find that the rapid SN explosion mechanism of F12 could produce more massive BHs than the results under B02. The consideration here demonstrates that the BH mass distribution is sensitive to the treatment of the SN mechanism (which determines the range of the BH formation as a function of the core mass) and the wind mass loss (which determines the mass of the pre-collapse progenitor star).As such, an increasing sample of the GWs from DCBs merging systems may tell this important information to understand the debated fate of massive stars and the unresolved mechanism of SN explosions.

Recipes for the wind mass loss, SN explosion mechanism, initial circular orbit and eccentric orbit affect final BH-BH populations. Other parameters have little influence on the chirp mass of BH-BH systems. Comparing with the standard model, the chirp mass distribution of BH-NS systems tends to peak in a smaller value if changes with the stellar evolution rather than it is treated merely as a constant (see Figure 5). Most parameters affect the BH-NS population, but there is not a clear trend in the dependence of the chirp mass distribution for different parameters. In Figure 6, it is shown that the , wind mass loss and eccentricity have clear effect on the chirp mass distribution of NS-NS systems while others have little influence.

### 3.3 Implications for the GW150914-like and GW151226-like binary black hole mergers

GW150914 is identified as a merging BH-BH binary. The estimated pre-merging primary and secondary BH masses are in the ranges of – 41 and – 33 (with a chirp mass around ), respectively (Abbott et al. 2016b). The detection of GW150914 improves our understanding of the BH-BH binary and BH mass which previously came only from X-ray BH binary studies. However, formation and evolution pathway of the progenitor of the GW150914-like binary black merger are still under debate. A number of recent works with different assumptions have explored the binary BH formation and merger rate: Kinugawa et al. (2014) investigated the formation of binary BHs through the population III stars (zero metallicity stars), and predicted high mass BHs (high chirp mass) as a typical outcome as consistent with the nature of GW150914. However, with different methods and arguments, Hartwig et al. (2016) and Dvorkin et al. (2016) claimed that GW150914 is unlikely formed from the population III stars. Furthermore, Marchant et al. (2016) and Mandel & de Mink (2016) proposed the so-called chemically homogeneous evolution model for the formation of massive BHs. Their models also have possible shortcomings – for example, the chemically homogeneous evolution needs a high rotation of stars and would only occur in a binary with a very short orbital period (1.5 - 2.5 days) having a massive companion star ( 40 ). Finally, Belczynski et al. (2016) argued that GW150914-like source is a natural outcome of the standard isolated binary evolution through the CE phase.

In sum, the formation and evolution of the progenitor system to realize GW150914 are still not fully known. Abbott et al. (2016d) summarized appropriate conditions of the isolated binary evolution scenario for the GW150914 as follows: (1) the massive star wind cannot be strong, which means a low metallicity environment is needed (see also Belczynski et al. 2016), (2) if the stars do not have high enough rotation and get small natal BH kicks, a successful evolution through the CE must be possible.

In this part, we use population II () stars under our three optimal models (other suitable conditions mentioned above are included in these models) to perform the BPS calculations in order to investigate a possible formation and evolution scenario of GW150914. As shown in the dotted region of Figure 1, we obtain the binary BH merger rates of for models 4–6 with 111The predicted NS-NS merger rates are 150–240 at the low metallicity (Pop II).. The common envelop parameters and eccentricity thus have a factor of few effect on the binary BH merger rate. On the other hand, the low metallicity leads to large enhancement of the BH-BH merger rate.

In our calculations at the solar metallicity (), the BH-BH system having the primary and secondary masses in the range observed for GW150914 is very hard to realize. Figure 7, showing the chirp mass distribution for the low metallicity model, shows that a larger number of massive BH-BH binaries can be formed from Pop II stars than from Pop I stars. For Models 4–6 under the low metallicity condition, a number of GW150914-like BH-BH systems are formed. Namely, a number of GW150914-like BH-BH systems can be produced from the low metallicity, Pop II stars. Our results thus support the possible origin of GW150914-like progenitors as proposed by Belczynski et al.(2016) and Abbott et al. (2016d), and we show that more realistic treatments on the binding parameter () indeed lead to the enhancement of GW150914-like BH-BH systems.

The nature of GW151226 is very different from that of GW150914. GW151226 has a chirp mass around . These BH masses are similar to a BH component found in the X-ray binaries (e.g. Ozel et al. 2010). A question is if these two sources are evolved through the same route, and what is then is the origin of the difference. In Figure 8, we investigate mass distribution of the merging BH-BH binaries for four different models and two different metallicities. By comparing the models with different metallicities in Figure 8, it is indeed seen that GW151226-like systems are typically expected at the solar metallicity (Pop I) while it is more likely to have GW150914-like systems at the lower metallicity (Pop II). We show, especially, that this different behavior at different metallicities is not obvious in our standard model, while it is the case for our optimized models with the realistic treatment in the CE evolution adopting the stellar wind model 2. We also note that for the low metallicity case the BH-BH binary systems typically have a more massive ‘secondary’ than the primary. 222In our definition, the primary is a more massive star at the formation of the binary system. This is likely caused by the less significant wind mass loss and more significant mass transfer for the lower metallicity environment. Finally, the delay time distributions of the BH-BH merger as shown in the right panel of Figure 9 also support the previous conclusions.

If the origins of GW150914 and GW151226 are both the isolated binary systems of massive stars through the CE evolution but with the difference in the metallicity (Pop II vs. Pop I), then an interesting question is if the combined constraint could be consistently satisfied within the same framework. Especially, the relative rate of the GW150914-like events and GW151226-like events will provide the information about the cosmic metallicity (thus the cosmic time) when the transition might have happened on the formation of the BH with to that with . Alternatively, this will provide a sanity check if the standard CE evolution toward the BH-BH binary formation is compatible to the frequencies of these GW sources.

Under the assumption that BH-BH binaries from Pop II stars represent GW150914-like events and those from Pop I stars represent GW151226-like events, we could estimate the relative rate of these two types of events by a single parameter that roughly discriminates the massive star binary evolution toward a GW150914-like event and that toward a GW151226-like event. Note that in reality we expect that this transition from the GW150914-like events to GW151226-like events as a function of the metallicity would not be extremely sharp, since the expected chirp mass both for Pop I and Pop II cases has a wide distribution and there is an overlap between the two cases. Still, deriving the characteristic metallicity as a dividing line between the two cases can provide us an important insight. Indeed, once a larger number of detected events would give good statistics, one may investigate further the distribution of the chirp mass as a function of the metallicity.

Irrespective of the value of required for particular binary evolution models (i.e., in our model must be between ), we vary as a parameter and compute the resulting relative ratio of the two types of the GW sources. This way, we can provide a strong test for the origins of the two GW sources as Pop II and Pop I systems, or alternatively the insight into the still-unknown treatment of the binary evolution and the mass loss rate.

Figure 10 shows the relative ratio of the Pop I to the Pop II BH-BH merging systems as a function of . We adopt the cosmic star formation history (as a function of redshift) from Madau & Dickinson (2014),

 SFR(z)=0.015(1+z)2.71+[(1+z)/2.9]5.6M⊙Mpc−3yr−1, (17)

where z is the redshift. We take the fitted values from the Figure 17 of Savaglio et al. (2005) for the relation of the cosmic metallicity evolution with the cosmic time and the redshift. Convolving the formation rates and the delay times of the BH-BH systems based on our BPS study with our three ideal models (models 4–6), we obtain the relative merger rate at the present Universe. While the observed statistics is still very poor, an interesting picture is already emerging. The detections of GW151226 and GW150914 indicate that the relative merger ratio of these different types of BH-BH systems should be (the thick orange line in the Figure 10), or (the gray area in the Figure) with the still poor statistics (Abbott et al. 2016f). This requires . This is indeed consistent with the binary evolution model we used, adding another support for the origins of GW150914 and GW151226 as the massive star binaries though the CE evolution at different metallicities (Pop II and Pop I). As such, further observationally constraining the relative rate of the similar systems will provide a powerful test on the natures of the progenitors of these GW emitters.

### 3.4 Implications for the NS-NS merger rate

The NS-NS merger rates calculated in our BPS study (i.e., 21 and 23.2 ) could be consistent with the NS-NS merger rate estimated indirectly from the observed NS-NS systems (median value of 21 ) indicated by Kim et al. (2015), if the wider orbital period distribution, binding energy parameter and (or ) are adopted. However, the same conditions in our BPS simulations, lead to the lower BH-BH merger rate than the BH-BH merger rate obtained through the statistics of the detected GW events from the BH-BH systems. For the conditions required to explain the observed BH-BH GW systems through the isolated binary evolution models, we thus predict that the NS-NS merger rate should be around the upper limit (79 Myr-1) so far obtained (Kim et al. 2015) or even higher than the previous estimate.

The mass transfer process deserves further discussion. For the most plausible (psychically motivated) prescription adapting qcr=qcr3, the predicted NS-NS merger rates become very high (see §3.1). We note that the predicted NS-NS merger rates, under the condition to explain the BH-BH merger rate in the same context, is much higher than that of Chruslinska et al. (2017), who did not consider the BH-BH merger rate and the NS-NS merger rate within the same context nor consider the spin evolution.

Our study highlights the importance of studying the BH-BH merger and NS-NS merger self-consistently in the same scheme. We show that (1) If the BH-BH mergers are the outcome of the isolated binary evolution, the NS-NS rate will be near the upper limit of the currently available estimate from the observed NS-NS systems or likely even larger, or (2) if the real NS-NS merger rate is consistent with the indirect estimate, then the BH-BH mergers require to come from the evolution beyond the isolated binary scenario. The future GW detection from NS-NS merger would provide a wealth of information on the binary evolution, thus the origin(s) of different DCB systems and the GW sources.

## 4 Conclusions and Discussion

With the Monte Carlo BPS calculations, we have studied binary evolutionary paths toward detached double compact star binaries (BH-BH, BH-NS, NS-NS binaries). In doing this, we have investigated effects of different recipes to simulate key binary evolution processes. Main conclusions of this paper are summarized as follows:

1. A larger number of BH-BH mergers can be produced by taking the gravitational energy, internal energy and entropy of the envelope () into account for the CE ejection, and considering one non-conservative mass transfer (). More BH-NS and NS-NS binaries are generated with and the mass transfer model 3 taking into account the spin of the star (). We have shown that two recipes of and other initial conditions have small effect on the parameter distributions and rate of the resulting double BH systems, while metallicity and wind mass loss models substantially affect natures of the BH-BH systems. The different and (including the conservative and non-conservative mass transfer prescriptions) affect the final features and rate of NS-NS binaries.

2. More BH-BH sources having the higher chirp masses are produced with the parameter combination of , , wind 2 model and SN mechanism of Fryer et al.(2012). The metallicity is another key parameter, and the low metallicity results in a larger chirp mass.

3. The origins of GW150914 and GW151226 are consistent with the isolated binary evolution model through the CE evolution, where the main difference in the nature (masses) is the metallicity at the formation of the binary systems. In this interpretation, GW150914 can be explained by a Pop II system while GW151226 is from a Pop I system.

4. The relative ratio of Pop I BH-BH mergers and Pop II BH-BH mergers as inferred by detections of GW151226 and GW150914-like events indicates that the transition in the binaries evolution toward the GW150914-like events and GW151226-like events would happen at a metallicity . This is consistent with the binary evolution model we adopt, adding another support for this interpretation. The promising future detections of a number of the similar events will provide a very strong constraint on the binary evolution and/or the origin of these GW sources.

5. Our calculated merger rates can be consistent with the NS-NS merger rate indicated by the indirect observational result (Kim et al. 2015). However, the same model would not explain the BH-BH merger rate constrained by the direct GW defections. We find that our different assumptions could change the NS-NS merger more than a factor of 7. We conclude that if the detected BH-BH mergers mergers are the outcome of isolated binary evolution, then the NS-NS merger rate must be at the upper limit of the presently available indirect estimate or even higher. Otherwise, a different mechanism toward the BH-BH mergers should be invoked.

This work was funded by JSPS Postdoctoral fellowship of Japan (P17022). The work was also supported by JSPS KAKENHI (No. 26800100, 17H02864) from MEXT and by WPI Initiative, MEXT, Japan. The authors thank the Kyoto University Foundation for their support. REFERENCES Abadie, J., et al. 2010b, Classical and Quantum Gravity, 27, 173001 Abbott, B. P., et al. 2016a, arXiv, arXiv:1602.03837 Abbott B. P., et al., 2016b, arXiv, arXiv:1602.03840 Abbott, B. P., et al. 2016c, arXiv, arXiv:1602.03842 Abbott B. P., et al., 2016d, arXiv, arXiv:1602.03846 Abbott B. P., et al., 2016e, arXiv, arXiv:1606.04855 Abbott B. P., et al., 2016f, arXiv, arXiv:1606.04856v2 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, PRL, 118, 221101 Abbott, B. P., et al. 2017b, arXiv, arXiv:1709.09660v1 Ablimit, I., & Li, X.-D. 2015, ApJ, 800, 98 Ablimit, I., Maeda, K. & Li, X.-D. 2016, ApJ, 826, 53 Belczynski, K., Kalogera, V. & Bulik, T. 2002, ApJ, 572, 407 Belczynski, K., Kalogera, V., Rasio, F. A., Taam, R. E., Zezas, A., Bulik, T., Maccarone, T. J., & Ivanova, N. 2008, ApJS, 174, 223 Belczynski, K., Bulik, T., Fryer, C., et al. 2010, ApJ, 714, 1217 Belczynski, K., Repetto, S., Holz, D., O’Shaughnessy, R., Bulik, T., Berti, E., Fryer, C., & Dominik, M. 2015, arXiv, arXiv:1510.04615 Belczynski, K., Holz, D., Bulik, T., & O’Shaughnessy, R. 2016, arXiv:1602.04531v1 Chruslinska, M., Belczynski, K., Klencki, J. & Benacquista, M. 2017, arXiv:1708.07885 Clausen, D., Piro, A. L. & Ott C. D. 2015, ApJ, 799, 190 Davis, P. J., Kolb, U., Willems, B. & Gänsicke, B. T., 2008, MNRAS, 389, 1563 Davis, P. J., Kolb, U., & Knigge, C. 2012, MNRAS, 419, 287 Dewi, J. D. M., & Pols, O. R. 2003, MNRAS, 344, 629 de Mink, S. E., Pols, O. R., & Hilditch, R. W. 2007, A&A, 467, 1181 de Mink, S. E. & Belczynski, K. 2015, ApJ, 814, 58 Dominik, M., Belczynski, K., Fryer, C., Holz, D. E., Berti, E., Bulik, T., Mandel, I. & O’Shaughnessy, R. 2012, ApJ, 759, 52 Dominik, M., Berti, E., O’Shaughnessy, R. et al. 2015, ApJ, 806, 263 Dvorkin, I., Vangioni, E., Silk, J., Uzan, J., Olive, K., 2016, ArXiv: 1604.04288v1 Eldridge, J. J., & Stanway, E. R. 2016, MNRAS, 462, 3302 Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126 Einstein, A. 1918, Sitzungsberichte der Koniglich Preubischen Akademie der Wissenschaften (Berlin), Seite 154-167., 154 Fryer, C. L., Belczynski, K., Wiktorowicz, G. et al. 2012, ApJ, 749, 91 Han, Z., Podsiadlowski, P. & Eggleton, P. P. 1994, MNRAS, 270, 121 Hansen, B. M. S., & Phinney, E. S. 1997, MNRAS, 291, 569 Harry, G. M. and LIGO Scientiﬁc Collaboration, 2010, Classical and Quantum Gravity 27, 084006 Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E., Magg M. & Stacy A. 2016, ArXiv: 1603.05655 Heggie D.C., 1975, MNRAS, 173, 729 Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 Humphreys, R. M. & Davidson, K. 1994, PASP, 106, 1025 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 89 Iben, Jr.,I. & Livio, M. 1993, PASP, 105, 1373 Ivanova, N. & Chaichenets S. 2011, ApJ, 731, L36 Ivanova, N., Justham, S., Chen, X., et al. 2013, ARA&A, 21, 59 Kiel, P. D., & Hurley, J. R. 2006, MNRAS, 369, 1152 Kim C., Perera B. B. P. & McLaughlin M. A. 2015, MNRAS, 448, 928 Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, MNRAS, 442, 2963 Kochanek, C. S. 2014, ApJ, 785, 28 Kochanek C. S. 2015, MNRAS, 446, 1213 Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545 Kulkarni, S. R., Hut, P., & McMillan, S. 1993, Nature, 364, 421 Lipunov, V. M., Postnov, K. A., & Prokhorov, M. E. 1997, MNRAS, 288, 245 Madau, P. & Dickinson, M., 2014, ARA&A, 52, 415 Mandel, I., de Mink, S. E., 2016, MNRAS, 458, 2634 Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., Moriya, T. J., 2016, A&A, 588, A50 Mennekens, N., & Vanbeveren, D. 2014, A&A, 564, A134 Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJ, 763, L15 Narayan, R., Paczynski, B., & Piran, T. 1992, ApJ, 395, L83 Nelemans, G., Yungelson, L. R., & Portegies Zwart, S. F. 2001, A&A, 375, 890 Nishizawa, A., Berti, E., Klein, A. & Sesana, A. 2016, arXiv:1605.01341v1 Nomoto, K., & Kondo, Y. 1991, ApJL, 367, L19 O¡¯Connor, E. & Ott, C. D. 2011, ApJ, 730, 70 O’Leary, R. M., Meiron, Y., & Kocsis, B. 2016, arXiv:1602.02809 Ozel, F., Psaltis, D., Narayan, R. & McClintock, J. E. 2010, ApJ, 725, 1918 Paczyński, B. 1976, in Structure and Evolution of Close Binary Systems, eds. P.Eggleton, S. Mitton, & J. Whelan (Dordrecht: Kluwer), IAU Symp., 73, 75 Phinney E. S., 1991, ApJL 380, L17 Postnov, K. A., Yungelson, L. R. 2014, LRR, 17 Sana, H., et al. 2012, Science, 337, 444 Savaglio, S., Glazebrook, K., Le Borgne, D. et al., 2005, ApJ, 635, 260 Sion, E. M., Bond, H. E., Lindler, D., et al. 2012, ApJ, 751, 66 Shao, Y., & Li, X.-D. 2014, ApJ, 796, 37 Sengupta, A. S., LIGO Scientific Collaboration and Virgo Collaboration, 2010, Journal of Physics Conference Series, 228, 012002. Spera, M., Mapelli, M. and Bressan, A., 2015, MNRAS, 451, 4086 Somiya, K. 2012, Classical and Quantum Gravity, 29, 124007 Stevenson, S., Ohme, F., & Fairhurst, S. 2015, ApJ, 810, 58 Sukhbold, T. & Woosley, S. E. 2014, ApJ, 783, 10 Tagawa, H., Umemura, M., & Gouda, N. 2016, arXiv:1602.08767 Tanikawa, A. 2013, MNRAS, 435, 1358 van den Heuvel, E. P. J., 1968, BAN, 19, 449V van der Sluys, M. V., Verbunt, F., & Pols, O. R. 2005, A&A, 431, 647 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 Voss, R., & Tauris, T. M. 2003, MNRAS, 342, 1169 Wang, C., Jia, K. & Li, X.-D. 2016, MNRAS, 457, 1015 Wang, C., Jia, K. & Li, X.-D. 2016b, RAA, 16h, 9 Webbink, R. F. 1984, ApJ, 277, 355 Woosley, S. E. 2016, arXiv, arXiv:1603.00511v1 Wyse, R. F. G. 2009, in IAU Symposium, Vol. 258, IAU Symposium, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, 11-22 Xu, X.-J., & Li, X.-D. 2010, ApJ, 716, 114 Yungelson, L. R., Lasota, J.-P., Nelemans, G., Dubus, G., van den Heuvel, E. P. J., Dewi, J., & Portegies Zwart, S. 2006, A&A, 454, 559
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters