Merger rate density of binary black holes through isolated Population I, II, III and extremely metal-poor binary star evolution

We investigate the formation of merging binary black holes (BHs) through isolated binary evolution, performing binary population synthesis calculations covering an unprecedentedly wide metallicity range of Population (Pop) I, II, III, and extremely metal-poor (EMP) binary stars. We find that the predicted merger rate density and primary BH mass ($m_1$) distribution are consistent with the gravitational wave (GW) observations. Notably, Pop III and EMP ($<10^{-2}$ $Z_\odot$) binary stars yield most of the pair instability (PI) mass gap events with $m_1 = 65$--$130$ $M_\odot$. Pop III binary stars contribute more to the PI mass gap events with increasing redshift, and all the PI mass gap events have the Pop III origin at redshifts $\gtrsim 8$. Our result can be assessed by future GW observations in the following two points. First, there are no binary BHs with $m_1=100$--$130$ $M_\odot$ in our result, and thus the $m_1$ distribution should suddenly drop in the range of $m_1=100$--$130$ $M_\odot$. Second, the PI mass gap event rate should increase toward higher redshift up to $\sim 11$, since those events mainly originate from the Pop III binary stars. We find that the following three assumptions are needed to reproduce the current GW observations: a top-heavy stellar initial mass function and the presence of close binary stars for Pop III and EMP binary stars, and inefficient convective overshoot in the main-sequence phase of stellar evolution. Without any of the above, the number of PI mass gap events becomes too low to reproduce current GW observations.

The isolated binary scenario, however, has been challenged by the discovery of GW190521 with a ∼ 90 M BH with the 90 % credible interval (Abbott et al. 2020a), which is prohibited to form by the pulsational pair instability (PPI) and pair instability (PI) supernovae (Heger et al. 2003;Belczynski et al. 2016b;Woosley 2017Woosley , 2019), in contrast to the hierarchical merger scenario in star clusters (Rodriguez et al. 2019;Di Carlo et al. 2020a;Liu & Bromm 2020b;Anagnostou et al. 2020b;Tagawa et al. 2021b;Kimball et al. 2021;Mapelli et al. 2021b;Liu & Lai 2021;Gayathri et al. 2021;Gerosa & Fishbach 2021;Arca Sedda et al. 2021). There are, however, still some way outs to accommodate the existence of GW190521 in the framework of the isolated binary scenario, which can be classified into three categories. The first possibility is that GW190521 actually straddles the PI mass gap; the primary and secondary BHs have masses beyond and below the PI mass gap, respectively (Fishbach & Holz 2020;Nitz & Capano 2021;Estellés et al. 2021). The second is that the PI mass gap may shift upward because of microphysics in single star evolution different from what has been thought: e.g., 12 C(α, γ) 16 O reaction rate three times smaller than the standard one (Farmer et al. 2019(Farmer et al. , 2020Costa et al. 2021), or dark matter effects (Croon et al. 2020;Ziegler & Freese 2021). In particular, Belczynski (2020) have shown that isolated binary stars can form the PI mass gap events consistent with GW190521-like event rate, taking into account the reduced 12 C(α, γ) 16 O reaction rate. The third is that the lower PI mass gap, ∼ 65-100 M , may disappear with more detailed consideration of macrophysics in single star evolution, such as metal-poor star evolution (Kinugawa et al. 2021a;Farrell et al. 2021;Vink et al. 2021;Tanikawa et al. 2021a), strength of stellar wind mass loss (Belczynski et al. 2020b), and stellar rotation effect . In this paper, we focus on the third possibility, to be more precise metal-poor star evolution, to accommodate the formation of GW190521-like events in the isolated binary scenario. Tanikawa et al. (2021a) have shown that Pop III stars can produce sufficient PI mass gap events to be reconciled with the GW190521-like event rate inferred by GW observations, if convective overshoot between their convective core and radiative envelope is inefficient during the main-sequence (MS) phase. Pop III stars, however, produce fewer binary BH mergers in the lower mass range of 50 M than the expectation from GW observations. Thus, we need to take into account not only Pop III stars but also Pop I/II stars to explain all the observed BH mergers in the framework of the isolated binary scenario.
For this purpose, we perform binary population synthesis (BPS) calculations for isolated Pop I/II/III and extremely metal-poor (EMP) binary stars in an unprecedentedly wide metallicity range in this paper. The reason for remarking the "unprecedentedly wide metallicity range" is that we cover not only Pop I/II stars but also EMP and Pop III stars, which are strictly defined at the end of this section. Previous studies of binary BHs have focused on only Pop I/II stars as most of BPS codes support only Pop I/II stars like SeBa (Portegies Zwart & Verbunt 1996;Toonen et al. 2012), BSE (Hurley et al. 2002), binary c (Izzard et al. 2009), MOBSE , COSMIC (Breivik et al. 2020), and COM-PAS (Team COMPAS et al. 2021). Some studies have examined Pop III stars (Kinugawa et al. 2014;Hartwig et al. 2016), and StarTrack (Belczynski et al. 2002 can cover Pop I/II/III stars. However, no studies have included EMP stars. We investigate isolated Pop I/II binary evolution by our own calculations rather than making use of the results of previous studies on isolated Pop I/II binary stars as done by Kinugawa et al. (2021c). We do this for consistency. We construct our own Pop III star evolution model as in Tanikawa et al. (2021a), decreasing the metallicity of our own Pop I/II star evolution model, rather than adopting the widely-used SSE Pop I/II star evolution model (Hurley et al. 2000). It is unclear if our own Pop I/II star evolution models yield similar binary BHs to the SSE Pop I/II star evolution models, or if our own model Pop I/II star evolution models can produce binary BHs in the mass range of 50 M consistently with the GW observations, and do not overproduce the PI mass gap events compared to the GW observations. We expect that the SSE and our own Pop I star evolution models would yield similar binary BHs, since these models are similar because of similar calibration to nearby Pop I stars including the Sun (Pols et al. 1998;Hurley et al. 2000;Yoshida et al. 2019). On the other hand, this may not be the case for the SSE and our own Pop II evolution models. Therefore, we use the SSE model for Pop I stars, and our own model for Pop II stars by our own calculations. We could have borrowed Pop I binary BHs from previous studies. However, we do not so. This enables us to examine Pop I binary BHs as well as Pop II, III and EMP binary BHs in more detail, since we obtain their properties by ourselves. For example, we can investigate the redshift evolution of binary BH mergers formed from Pop I/II/III and EMP binary stars.
In this paper, we examine the isolated binary scenario by the BPS calculations, and demonstrate that from the isolated Pop I/II/III and EMP binary stars the observed BH merger rate density and its primary BH mass distribution over ∼ 5-100 M are successfully reproduced under the assumption of inefficient convective overshoot. This means that isolated Pop III (and EMP) binaries can form GW190521-like events, and that isolated Pop I/II star evolution models can produce binary BHs in the mass range of 50 M consistently with GW observations without overproducing PI mass gap events regardless of the choice of Pop I/II star evolution models. We also predict the redshift evolution of merging binary BHs, which should help assessing the validity of our model by future GW observations, e.g., by Cosmic Explorer (Reitze et al. 2019) and Einstein Telescope (Punturo et al. 2010). We also point out the key assumptions allowing the predicted BH properties to be consistent with the observations, considering large uncertainties in the initial conditions and evolution of Pop III stars due to its non detection so far (Frebel & Norris 2015;Magg et al. 2018Magg et al. , 2019. The structure of this paper is as follows. We present our numerical method to calculate the BH merger rate density at a given redshift in section 2. We present and discuss the numerical results in section 3 and 4, respectively. We summarize this paper in section 5. We adopt the cosmological parameters based on the latest observation by Planck Collaboration et al. (2020). We define metallicities of Pop I, Pop II, EMP, and Pop III stars as Z/Z > 0.16, 10 −3 < Z/Z ≤ 0.16, 0 < Z/Z ≤ 10 −3 , and Z/Z = 0, respectively, where Z is the metallicity, Z is the solar metallicity set to Z = 0.02. We treat Z/Z = 10 −8 stars as Pop III stars, Z/Z = 10 −6 and 10 −4 as EMP stars, Z/Z = 10 −2 , 0.025, 0.05 and 0.1 as Pop II stars, and Z/Z = 0.25, 0.5 and 1 as Pop I stars.

METHOD
We describe initial conditions for our BPS calculations in section 2.1. We overview our single and binary star models for our BPS calculations in section 2.2. As we take into account several different initial conditions and single/binary star models, we describe our parameter sets for them in section 2.3. In section 2.4, we formulate how to obtain the BH merger rate density and its derivatives, using our BPS calculation results.

Initial conditions
We account for the star formation rate (SFR) density with two components as where z is the redshift, and φ I/II (z) and φ III (z) are the Pop I/II and Pop III SFR densities, respectively. We adopt the Pop I/II SFR density of Madau & Fragos (2017), expressed as We assume the metallicity distribution of Pop I/II components at a given z as a logarithmic normal distribution similarly to Belczynski et al. (2020a) and Santoliquido et al. (2021), given by where Z I/II (z) is the average metallicity of the Pop I/II component at a given redshift z, and σ Z is the metallicity dispersion set to σ Z = 0.35. The average metallicity evolution in the Pop I/II component Z I/II (z) is considered in the same way as Madau & Fragos (2017): We use the following Pop III SFR density, obtained by simplifying the simulation result of Skinner & Wise (2020) as shown in their figure 4: where t indicates the cosmic time since the Big Bang. The three formulae in Eq. (5) are not continuous, since the numbers in the three formulae are written to 3 significant digits. We regard that Z = 10 −8 Z stars, which are the most metal-poor in our single star model, are in our simulations equivalent to Pop III (i.e. zero-metal) stars. This is because their evolution is almost identical. Thus, we suppose the metallicity distribution of the Pop III component as where the δ function has the usual property as seen above. This means that Z = 10 −8 Z (i.e. Pop III or zero-metal) stars have the SFR density described in Eq.
(5) throughout all the times (or redshifts). Then, we can get the SFR density differentiated by the metallicity: This equation is used in section 2.4. We assume that the binary number fraction is 0.5 for all metallicities, where the fraction is close to an intrinsic binary fraction 0.69 +0.09 −0.09 found by Sana et al. (2012). The initial mass function (IMF) of the single stars and primary stars of the binary stars is given by a combination of the following two mass functions. The first is the Kroupa's IMF (Kroupa 2001) written as where m is the stellar mass. The second is a function flat in the logarithmic scale (hereafter, the logarithmically flat IMF), such that We add the above two mass functions with a weight such that more stars follow the logarithmically flat IMF at lower metallicities. We adopt three different methods to weight the two functions. The first is to adopt just the Kroupa's IMF for all the metallicities. We call this IMF model the "top-light IMF model". The second is to use the Kroupa's IMF for Z/Z ≥ 10 −5 and otherwise the logarithmically flat IMF. This is because the stellar IMF is theoretically predicted to transition from toplight to top-heavy at Z/Z ∼ 10 −5 (Bromm & Larson 2004;Omukai et al. 2005;Schneider et al. 2006;Maio et al. 2010). We call this IMF model the "stepwise IMF model". In the third method, the IMF shifts gradually from Kroupa's IMF to the logarithmically flat one with decreasing metallicity. All the stars follow the Kroupa's IMF for Z/Z > 10 −2 . At Z/Z = 10 −2 , 60 % of them obey the logarithmically flat IMFs in terms of mass, while this fraction increases to 90 % at Z/Z = 10 −4 . At metallicity as low as Z/Z ≤ 10 −6 , all the stars follow the logarithmically flat distribution. This behavior is inspired from the simulation results of Chon et al. (2021). We call this IMF model the "transitional IMF model".
For the initial distributions of binary parameters we refer to Sana et al. (2012), who proposed the following expressions for the mass ratios (q), periods (P ), and orbital eccentricities (e) of binary stars: f (log P ) ∝ (log P ) −0.55 (0.15 < log(P/day) < 5.5) f (e) ∝ e −0.42 (0 ≤ e ≤ 1), where q is defined as the ratio of the secondary mass to the primary mass. We call this binary star initial condition the "Sana's model". The period range in Eq. (12) is the same as one adopted by de Mink & Belczynski (2015). On the other hand, our binary number fraction (0.5) is different from that of de Mink & Belczynski (2015) who have adopted 1 for the binary number fraction as default by careful consideration of the difference between their and Sana's period ranges.
In addition to this, we introduce another binary star initial condition, where we exclude binary stars with pericenter distances of less than 200 R when Z ≤ 10 −2 Z . This condition is inspired by the following consideration. Z/Z 10 −2 stars tend to expand to ∼ 100 R in radius during their protostellar evolution thanks to the high mass accretion rate (Hosokawa & Omukai 2009a,b;Omukai et al. 2010). Even if a metalpoor star initially has a close companion star, they may merge during their protostar phase, and may not make a close binary system. We call this binary initial condition the "wide binary model". Nonetheless, we still adopt the Sana's model as the fiducial one. This is because the radial expansion depends on the geometry of protostellar accretion, which is uncertain. If the accretion shifts from spherical to disk accretion relatively early on, the maximum expansion radius would be less than 100 R ). In addition, even with Z < 10 −2 Z close binaries can be formed through Nbody interaction after their maximal expansion phase (e.g. ).

Single and binary star models
We employ the SSE model for single stars with 5 × 10 −3 ≤ Z/Z ≤ 1.5 (Hurley et al. 2000), and M and L models 1 for single stars with Z/Z = 10 −8 , 10 −6 , 10 −4 , and 10 −2 ≤ Z/Z ≤ 0.1 (Tanikawa et al. 2020(Tanikawa et al. , 2021a. An L-model star experiences more efficient convective overshoot than an M-model star does, for which the convective overshoot happens between the convective core and radiative envelope in the MS phase. The L-model star has a larger He core than the M-model star at the end of their MS phase if these stars have the same initial mass. The L-model star expands much more in radius than the M-model star at their post-MS phase. We summarize the details of the M/L models, and compare them with the SSE models in appendix A. We describe later in this section how to adopt the SSE, M, and L models. We take into account stellar wind mass-loss in the same way as Tanikawa et al. (2020), based on the wind model of Belczynski et al. (2010) with modification of metallicity-dependent luminous blue variable winds . For simplicity, we do not include rotational enhancement of stellar winds (unlike as in Tanikawa et al. 2020).
Our supernova model is the rapid model  with the PI mass loss. The PI mass loss modifies a remnant mass given by the rapid model (m rapid ) as 1 The two models are named after the Milky way and the Large Magellanic Cloud, whose stars were used to calibrate the overshoot parameter of the M and L models, respectively ). follows: where m rem is the final remnant mass, m c,He is the He core mass, and m c,He,PPI , m c,He,PISN , and m c,He,DC are the lower He core mass limits of pulsational PI (PPI: Heger & Woosley 2002;Woosley et al. 2007;Yoshida et al. 2016;Leung et al. 2019), PI supernova (PISN: Barkat et al. 1967;Fraley 1968;Bond et al. 1984;El Eid & Langer 1986;Fryer et al. 2001;Heger & Woosley 2002;Umeda & Nomoto 2002), and direct collapse (DC), respectively. Note that, although m c,He,PPI is fixed to a constant value in the above equation, in reality it should depend on m c,He in some complex way (see figure 1 in Stevenson et al. 2019). We consider two sets of the mass limits, since the mass limits strongly depend on the uncertain 12 C(α, γ) 16 O reaction rate (Farmer et al. 2020;Costa et al. 2021). The first and second sets are (m c,He,PPI , m c,He,PISN , m c,He,DC ) = (45M , 65M , 135M ) and (90M , 90M , 180M ), which are similar to Belczynski et al. (2016b) and Belczynski (2020), respectively. The first set is based on the standard 12 C(α, γ) 16 O reaction rate with respect to the STARLIB (Sallaska et al. 2013). Thus, we call the PI model "the standard PI model". On the other hand, the second set can be obtained when the 12 C(α, γ) 16 O reaction rate is lower than the standard one by 3σ. Thus, we call the PI model "the 3σ PI model".
We have four different combinations of single star models and PI models as seen in  (Pols et al. 1998;Hurley et al. 2000). Rather, there is a possibility that our results are affected by the choice of the SSE and M/L models for more metal-poor stars, i.e. Pop II stars. However, we can see that there is little difference yielded by the choice of the SSE and M/L models for Pop II stars as seen in appendix B. Despite that the SSE model supports metallicity down to 5 × 10 −3 Z , we do not use the SSE model with less than 10 −2 Z . The practical reason for this is that we would like to adopt the same metallicities for Pop II stars among the SSE and M/L models, and that the M/L models do not support 5 × 10 −3 ≤ Z/Z < 10 −2 star evolution. However, we can see that the not using the SSE models for this range has no impact (see details in appendix C). We adopt the M/L models for Pop III+EMP stars, since the SSE model does not fully support Pop III+EMP star evolution. Especially, we regard the M/L models with Z/Z = 10 −8 as Pop III stars. We model the BH spin angular momentum (J rem ) as follows. We know spin angular momenta of He core (J c,0 ) and H envelope (J e,0 ) for each star at the presupernova phase from the stellar evolution model. We set where m c,0 and m e,0 are the masses of He core and H envelope at the pre-supernova phase, and m c,1 and m e,1 are the masses of He core and H envelope contributing to the BH mass, i.e., the sum of m c,1 and m e,1 is the BH mass. We obtain the above relation, assuming that the He core and H envelope are in rigid rotation, and they each have uniform densities. In this prescription, we may underestimate the BH spin angular momentum, since stellar rotation and density become slower and smaller outward, respectively. We can calculate the dimensionless BH spin as where G is the gravitational constant, and c is the light speed.
We include natal kicks of neutron stars and BHs because of asymmetric supernova explosions. The velocity distribution is given by a single Maxwellian with 265 km s −1 (Hobbs et al. 2005) if the remnant acquires no fallback mass. The natal kick velocity is reduced by 1 − f b , where f b is the fraction of the fallback mass . The direction of the natal kick is assumed to be isotropic. Our binary star model is based on the BSE model (Hurley et al. 2002). We overview the parameters of our choice, and the different points between the BSE and our models. The BSE model includes wind accretion, tidal evolution, stable mass transfer, common envelope evolution, magnetic braking, and orbital decay through GW radiation. The wind accretion in our model is limited by the Eddington limit expressed by Cameron & Mock (1967), while the wind accretion in the BSE model is not. We adopt this limitation, since wind accretion may exceed the Eddington limit when a star has a massive companion. The BSE and our models have the same prescription for tidal evolution of stars with convective envelopes. On the other hand, we adopt two different prescriptions for tidal evolution of stars with radiative envelopes, one of which is the same as the BSE model (described in detail later). The BSE and our models have different criteria for whether a star has radiative or convective envelope. In the BSE model, a massive star has radiative (convective) envelope at its core (shell) He burning phases, respectively. On the other hand, in our model, a massive star has radiative or convective envelope when its effective temperature (log(T eff /K)) is above or below 3.65, respectively. The same criteria are applied also for the stability of mass transfer. We need the different criteria, since the M-and L-model stars can have radiative envelopes in their shell He burning phases, and convective envelopes in their core He burning phases unlike the BSE model. In a stable mass transfer phase, the maximum fraction of transferred mass is 0.5 in our model, while it is unity in the BSE model. Such non-conservative mass transfer is also adopted by the latest BPS calculations (e.g. Belczynski et al. 2020a;Kinugawa et al. 2020). Our model as well as the BSE model apply the α formalism for the common envelope evolution (Webbink 1984). We set α = 1, and adopt the formulae of Claeys et al. (2014) for λ. We assume that when a star in the Hertzsprung gap phase fills its Roche lobe and its mass transfer is unstable, the star merges with its companion star. This is because a star in its Hertzsprung gap phase does not have steep density gradient between the He core and H envelope (Ivanova & Taam 2004). This treatment is similar to other binary population synthesis models (Dominik et al. 2012;. Note that such a star experiences common envelope evolution in the BSE model. We briefly describe the prescription for tidal evolution of stars with radiative envelope, which is similar to Kinugawa et al. (2020). In this case, the tidal evolution is subject to the radiative damping of the dynamical tide (Zahn 1975). The tidal coupling parameter k/T is given by (eq. (42) in Hurley et al. 2002, see also Zahn 1977Hut 1981), where R is the stellar radius and a is the binary semi-major axis. In the BSE model, the variable E depends on the stellar mass as Recently, Yoon et al. (2010, see also Qin et al. 2018 proposed an alternative fitting formula for the E based on the ratio of the convective core radius (R conv ) to the stellar radius, such that E = 10 −0.42 (R conv /R) 7.5 for H-rich stars 10 −0.93 (R conv /R) 6.7 for He-rich stars . (20) We set the convective core radius as We obtain the convective core radius for H-rich stars from our L-model star's data. We adopt the same convective core radius for He-rich stars as Kinugawa et al. (2020), although this value is fitted to extremely metalpoor stars. Hereafter, the prescriptions for tidal evolution expressed as Eqs. (19) and (20) are called "originaltide" and "new-tide" models, respectively.

Parameter sets
We investigate eight parameter sets with different IMFs, binary star initial conditions, and single and binary star models. We summarize the parameter sets in Table 2.
We regard as the fiducial one the parameter set in which its IMF, binary star initial condition, single star model, and binary star model are the transitional IMF, the Sana's model, the M and SSE-Pop I model with the standard PI model, and new-tide model, respectively. We present the results of the SSE-std set only in appendices B and C.
For parameter sets with the transitional IMF (e.g. the fiducial, wide-binary, original-tide, L-std, L-3σ, and SSE-std parameter sets), we prepare 12 groups of 10 6 binary stars, and perform BPS calculations for these groups. The 12 groups include 10 different metallicities: Z/Z = 10 −8 , 10 −6 , 10 −4 , 10 −2 , 0.25, 0.05, 0.1, 0.25, 0.5, and 1. The two most metal-poor groups follow the logarithmically flat IMF, and the six most metalrich groups follow the Kroupa's IMF. The two groups with Z/Z = 10 −4 and 10 −2 have a mixture of the Kroupa's and logarithmically flat IMFs. For the toplight and stepwise parameter sets, we prepare 10 groups of 10 6 binary stars with the 10 different metallicities. In the top-light parameter set, all the groups follow the Kroupa's IMF. For the stepwise parameter set, the two most metal-poor groups follow the logarithmically flat IMF, and the other groups follow the Kroupa's IMF.
We generate 10 6 binary stars without any restriction on the upper and lower mass limits in the case of the logarithmically flat IMF. The total mass of these binary stars is 8.2 × 10 7 M . Since we assume the binary number fraction to 0.5 as described above, the 10 6 binary star formation should be accompanied by 10 6 single star formation. The total mass of the 10 6 single stars and the 10 6 binary stars is 1.3 × 10 8 M . On the other hand, we generate 10 6 binary stars whose primary and secondary masses are more than 10 M for the Kroupa's IMF. This is because the number fraction of stars leaving behind BHs is small for the Kroupa's IMF. The total mass of these binary stars is 4.2 × 10 7 M . These binary stars should be accompanied by a large number of single stars and binary stars whose primary and secondary masses are less than 10 M . The total mass of the single and binary stars (including the 10 6 binary stars) is 5.2 × 10 8 M .

BH merger rate density
We use the results of BPS calculations, and obtain the BH merger rate density (R) and its derivatives by binary BH properties (P), such as primary mass, mass ratio, and their combination in the source frame. First, we express the differential BH merger rate density at a given look-back time as where t h is the Hubble time, t lb is the look-back time, and t d is the delay time of a BH merger from star formation. We can regard N (Z) as the total number of BH mergers per stellar mass formed at a given metallicity Z.
We call N (Z) "production efficiency of BH mergers". 2 We discretize the integral of Z in Eq. (22), and obtain the following equation: N Z is the number of different metallicities (i.e. 10), and Z i /Z = 10 −8 , 10 −6 , · · · , 0.5, 1 for i = 1, 2, · · · , 9, 10, respectively. ∆φ i (z) is the SFR density at a given redshift z in a metallicity range around metallicity Z i . The lower and upper boundaries of the metallicity range are the geometric means with its lower and upper neighbor metallicities: (Z i−1 Z i ) 1/2 and (Z i Z i+1 ) 1/2 , respectively. We define Z 0 = 0 and Z 11 = ∞. Thus, we can calculate ∆φ i (z) as where dφ(z)/dZ can be expressed as Eq. (8). Note that we use the Z = Z result for this metallicity range, We can also discretize the integral in Eq. (23) as follows. The derivative in Eq. (23) can be divided into two parts: binary stars subject to the Kroupa's and log-arithmically flat IMFs, such that where N IMF is the number of IMFs (i.e. N IMF = 1 or 2), f i,j is the mass fraction of stars formed in each IMF, and N i,j is production efficiency of BH mergers for each Z i and each IMF. Note that f i,1 + f i,2 = 1. Eq. (25) can be written as where N bps,i,j is the number of simulated binary stars for each Z i and each IMF (i.e. N bps,i,j = 10 6 ), M bps,i,j is the total stellar mass required for forming the N bps,i,j binary stars (i.e. M bps,i,1 = 5.2×10 8 M and M bps,i,2 = 1.3 × 10 8 M as described in section 2.3), and t i,j,k d and P i,j,k are the delay time and binary BH properties of each binary star. Using Eqs. (25) and (26), we can discretize Eq. (23) as The BH merger rate density can also be obtained as 3.

RESULTS
We present the results for the fiducial parameter set in section 3.1. We see how the results change in the cases with other parameter sets in section 3.2. Figure 1 shows the redshift evolution of the BH merger rate density for the fiducial set. We find R ∼ 20 yr −1 Gpc −3 at z = 0, comparable to the observed rate, 19.3 +15.1 −9.0 yr −1 Gpc −3 according to the LIGO-Virgo Gravitational-Wave Transient Catalog 2 (GWTC-2) (Abbott et al. 2021a) 3 . Note that our prediction is also consistent with the GWTC-2. merger rate density up to z = 1 also falls within the 90 % credible interval inferred by GWTC-2. The BH merger rate density monotonically increases with increasing redshift up to z 2.5, and then decreases monotonically for higher redshift. The BH merger rate density is dominated Pop I/II stars at redshifts of z 7. Pop III and EMP stars contribute only by ∼ 1 (10) % to the BH merger rate density at the redshift of z = 0 (z = 5, respectively). Their contribution becomes dominant at redshifts higher than z ∼ 7.

Fiducial parameter set
We find that the redshift at which the BH merger rate reaches the maximum (z ∼ 2.5) is slightly higher than the redshift of the maximum SFR (z ∼ 2). In naive expectations, the BH merger rate should achieve the maximum at some epoch after the peak SFR as binary BHs merge after their progenitor formation. However, this is not the case for the following reason. At z < 3, the SFR for Z ∼ 0.25 Z is lower than that for Z ∼ 0.5 (1) Z by a factor of about 3 (10). On the other hand, the production efficiency of BH mergers over the Hubble time (η(Z) = t h 0 dt d dN (Z)/dt d ) 4 for Z ∼ 0.25 Z is higher than those for Z ∼ 0.5 and 1 Z by a factor of about 30 and 1000, respectively, such that η(Z) = 9.1 × 10 −6 M −1 for Z = 0.25 Z , 3.0 × 10 −7 M −1 for Z = 0.5 Z , and 7.6 × 10 −9 M −1 for Z = 1 Z . Such sharp dependence on metallicity has also been pointed out by  and Klencki et al. (2018).
Merging binary BHs are thus mainly yielded by stars of Z 0.25 Z at z < 3. Moreover, about half of them have much shorter inspiral times than the Hubble time ( 100 Myr). The redshift evolution of the BH merger rate density follows the SFR of Z 0.25 Z stars, which achieves the maximum at z 2.5. We compare our obtaining BH masses with BH masses inferred by GWTC-2 (Abbott et al. 2021b). Figure 2 shows the distributions of the primary BH masses (top) and mass ratios (bottom) at z = 0, which are defined as the BH merger rate density differentiated by the primary BH mass and mass ratio, respectively. The primary BH mass distribution is roughly consistent with the observations in GWTC-2 in three respects. First, our primary BH mass distribution has the global maximum at m 1 ∼ 10 M . Second, the mass distribution suddenly decreases at m 1 ∼ 45 M due to the PPI effect. It should be noted that merging binary BHs with m 1 ∼ 45 M are overproduced while those with m 1 ∼ 35 and 50 M are underproduced, relative to GWTC-2, because of the standard PI modeling in which PPI produces BHs with a fixed mass (45 M ). If we had adopted a PPI model whose remnant mass has a broad distribution (e.g. Yoshida et al. 2016;Leung et al. 2019;Umeda et al. 2020), this discrepancy would have disappeared. Third, we find BHs within the PI mass gap, m 1 ∼ 65-90 M . As shown in Tanikawa et al. (2021a), Pop III stars of 65-90 M keep their small radii ( 40 R ), and do not lose their H envelopes through binary interactions. They have He cores with 45 M , and thus leave behind 65-90 M BHs without either PPI or PISN effect. As seen in comparison between the blue and cyan curves in Figure 2, not only Pop III binary stars but also EMP binary stars contribute to PI mass gap mergers. The results in Figure 2 also eliminate our concerns described in section 1. In other words, Pop I/II stars can produce binary BHs in the mass range of 50 M , and do not overproduce the PI mass gap events.
We can see that only Pop III and EMP stars can produce merging binary BHs within the PI mass gap. In the mass range 65-90 M , their radii remain small, and thus they can keep their H envelope until the collapse to BHs. They thus leave behind the PI mass gap BHs. On the other hand, Pop I/II stars in the mass range 65-90 M expand greatly, and lose their H envelope through binary interactions. They leave behind at most ∼ 45 M BHs. The radius evolution can be seen in the Appendix A. In advance, we should remark that the formation of merging binary BHs within the PI mass gap strongly depends on the choice of the single star model for Pop III and EMP stars. If we choose the L model, we do not have such binary BHs as seen in section 3.2.
The bottom panel of Figure 2 shows the mass ratio distribution of merging binary BHs. The distribution roughly follows the power law ∝ q 3 (dashed curve). According to Abbott et al. (2021b), the power-law index of the GWTC-2 events is less than or equal to 3. Our mass ratio distribution is thus marginally consistent with it.
We find that Pop I/II binary stars preferentially result in equal mass ratio merger events of q ∼ 1. Because of the large radii in the post-MS phase, they tend to experience binary interactions, such as common envelope and stable mass transfer. After the binary interactions, which work to equalize the binary component masses, they have similar BH masses. On the other hand, Pop III and EMP stars hardly experience such binary interactions, and keep their initial mass ratios and thus the initial flat q distribution (see section 2.1). In fact, their mass ratios follow nearly flat distribution as seen in the lower panel in Figure 2 (blue). Figure 3 shows the cumulative distribution of effective spins for the fiducial set of the parameters. The effective spin of a binary BH system is defined as where χ 1 and χ 2 are the dimensionless spin of the primary and secondary BHs, respectively, and θ 1 and θ 2 are the inclination angle between the primary and secondary spin vectors and the orbital angular momentum vector of the binary BH, respectively. All of merging binary BHs have zero or positive χ eff . This means that they have either no spin or spins aligned to their orbital angular momentum vectors, since they have θ 1 ∼ 0 and θ 2 ∼ 0. Such spinning BHs acquire their spin angular momenta in the progenitor phase through tidal interactions with the companion stars. About 10 % of merging BHs have χ eff ∼ 0.4, which are formed from Pop I binary stars in the following pathway. After the primary star in a binary system collapses to a BH with little spin, the secondary star becomes a giant, and fills its Roche lobe. The binary system then experiences the common envelope evolution, and shrinks its separation down to 10 R . The secondary star loses its H envelope and eventually becomes a naked He star. Because of the small separation, it is efficiently spun up by the tidal interaction (see the formulae in Eqs. (18), (20), and (21)), and collapses to a BH conserving its spin angular momentum (see Eqs. (15) and (16)). Since the resultant binary BHs have small separation, they merge after a short delay time of ∼ 100 Myr from their progenitor star formation. This has been argued by Kushnir et al. (2016) and Hotokezaka & Piran (2017).
Pop II binary stars leave behind merging BHs with a wide range of positive χ eff ∼ 0-1. These binary BHs especially dominate the population of merging BHs with χ eff 0.4. The reason for the wide range of χ eff is because their progenitor stars are spun up by tidal interactions, and lose their H envelope, i.e., large reservoir of spin angular momenta, to some extent through stable mass transfer. Note that their progenitors lose little mass through stellar winds because of their low metallicity, and tend not to experience the common envelope evolution since they become blue supergiant stars (see appendix A). The spin magnitudes of these binary BHs depend on whether their progenitors keep the H envelope.
Most of Pop III and EMP binary stars also leave behind merging BHs with higher χ eff than Pop II binary stars. Pop III and EMP stars with smaller radii, tend to keep more massive H envelope than more metal-enriched stars (see appendix A). Merging BHs resulting from Pop III and EMP binary stars however, do not much affect the cumulative χ eff distribution of all the merging binary BHs. This is because they account for only 1 % of the total BH merger rate density at z = 0.
There are, however, discrepancies between our cumulative χ eff -distribution and observations in three respects. First, no binary BH with negative χ eff is found in our sample. Second, our binary BHs are fewer than observed in the range χ eff ∼ 0-0.1. Third, there are more binaries BHs with χ eff ∼ 0.4 in our sample than what is observed. We nevertheless expect that we can mitigate these deviations once we account for observational uncertainties as in Bavera et al. (2020). Figure 3 also shows the cumulative distribution of binary BHs with m 1 /M ≥ 65, which are within the PI mass gap, at z = 0 and 1. More than 90 and 70 % of binary BHs within the PI mass gap have χ eff > 0.5 at the redshifts of z = 0 and 1, respectively. The events at higher redshift have smaller χ eff since the progenitor binary stars have smaller separations at higher redshift, and lose their H envelope more promptly.
We have presented the cumulative distribution of binary BHs with PI mass gap BHs both for z = 0 and 1, since some of the observed PI mass gap BHs are located closer to z = 1 rather than z = 0. Here, we regard the observed binary BHs with m 1 /M ≥ 65 as PI mass gap events even if they are not conclusive. For example, GW190521 is conclusive, since its m 1 is in the PI mass gap with the 90 % credible intervals (Abbott et al. 2020a,b with the 90 % credible intervals). This suggests that our model may still be consistent with the observation even considering the high probability of PI mass gap BHs with high spins in our sample.
We do not try to compare our model with the observations in terms of the so-called spin precession, i.e., the BH spin component perpendicular to the binary orbital angular momentum vector. This is because the estimates of spin precessions have larger uncertainties than those of effective spins, even though meaningful constraints may be imposed on spin precessions of a few events in GWTC-2 data (Abbott et al. 2021a), such as GW190412 (Abbott et al. 2020c Figure 4. BH merger rate density as a function of primary BH mass at the redshifts of 0 (black), 2 (blue), 4 (green), and 8 (red) for the fiducial set. The gray, cyan, light-green, and pink curves indicate Pop III contributions at the redshifts of 0, 2, 4, and 8, respectively. The shaded gray regions show the 90% credible interval inferred by the LIGO-Virgo Gravitational-Wave Transient Catalog 2 (GWTC-2), which is relevant for the local redshift, z < 2. The PI mass gap region is on the right-hand side of the dashed line. Figure 4 shows the redshift evolution of the primary BH mass distribution for the fiducial set of the parameters. The global maximum at m 1 ∼ 10 M increases from z = 0 to z = 2 − 4 by an order of magnitude, and then again decreases from z = 4 to z = 8 by an order of magnitude. This redshift evolution reflects that of the Pop I BH merger rate. Note that merging binary BHs around the global maximum results from the Pop I progenitors. We find that the BH merger rate density at m 1 ∼ 10 M grows much more than at m 1 ∼ 15 M from z = 0 to z = 2 and 4, despite that Pop I binary stars can yield binary BHs with both m 1 ∼ 10 and 15 M at the redshift of z = 0 (see Figure 2). This makes a spike in the BH merger rate density at m 1 ∼ 10 M at z = 2 and 4. This spike comes from different formation pathways between Pop I binary BHs with m 1 ∼ 10 and 15 M . The former binary BHs are formed through common envelope evolution, and typically have short delay time of ∼ 100 Myr from their progenitor star formation. On the other hand, the latter binary BHs are formed through stable mass transfer, and their typical merger time is ∼ 10 Gyr. Thus, the former merger rate density grow much more than the latter from z = 0 to z = 2 and 4. In advance, we note that this feature can be seen in all the sets (see Figure 6).
Interestingly, the BH merger rate density in the mass range m 1 ∼ 20-30 M does not change over the time interval z = 0-8. This is because they are originated mainly from Pop II binary stars (see Figure 2), and their merger rate density does not change over those redshifts (see Figure 1).
The merger rate density of binary BHs within the PI mass gap (∼ 65-90 M ) monotonically increases up to z = 8. They are formed from the Pop III and EMP progenitors (see Figure 2), whose merger rate increases with redshift up to z ∼ 8 (see Figure 1). This is because the Pop III SFR increases toward higher redshifts up to z ∼ 20.
We can see from Figure 4 that Pop III contributions in the PI mass gap become larger with increasing redshift. At z = 0, EMP binary stars contribute to the PI mass gap equally to Pop III binary stars. However, Pop III binary stars start dominating the PI mass gap from above. At z = 2, 4, and 8, Pop III binary stars have great roles on forming binary BHs with m 1 90, 80, and 70 M in the mass gap, respectively. We can expect that a mass gap event has its origin in a Pop III binary star if it is detected at a redshift z 8.  Figure 5 shows the redshift evolution of the cumulative effective spin distribution for the fiducial set. The contribution of binary BHs with χ eff ∼ 0.4 increases from z = 0 to z = 4. As described above, they are originated from Pop I binary stars, and have short delay time since their progenitor formation. Thus, their merger events occur more frequently at redshifts when star formation (especially of Pop I stars) is more active (i.e. z = 2 and 4). At z = 8, the fraction of binary BHs with χ eff 0.4 is low, while the fraction with higher χ eff is high. This is because at this redshift binary BHs are originated dominantly from the Pop II, Pop III and EMP progenitors, especially Pop III and EMP progenitors. Figure 6 shows the distributions of primary BH masses at the redshifts of z = 0, 2, 4, and 8 for parameter sets other than the fiducial one. We can see that in the case of the top-light set (top-left) the BH merger rate for m 1 /M 50 at z = 0 is definitely smaller than observed, although the rate for m 1 /M 50 at z = 0 is comparable to the observation. For the stepwise set (top-right), the BH merger rate of m 1 /M 50 at z = 0 is consistent with the observation, although the PI mass gap rate is smaller than that in the fiducial set. We also find that there are few binary BHs within the PI mass gap in the case of the wide-binary set (middleleft). These results demonstrate that the initial conditions of metal-poor binary stars have great impacts on the formation of binary BHs within the PI mass gap. Metal-poor binary stars should follow a top-heavy IMF and have small initial pericenter distance in order to produce a sufficient number of the PI mass gap events.

Other parameter sets
The stepwise set (top-right) shows an interesting result as for Pop III binary stars. Pop III binary stars dominate the PI mass gap events at any redshifts from z = 0 to 8. In other words, EMP binary stars form little PI mass gap event. This is because all the Z/Z = 10 −4 stars have a top-light IMF, and the star formation rate of Z/Z = 10 −6 stars is small. If IMFs depend on metallicity stepwisely like the stepwise set, the PI mass gap events immediately indicate the presence of Pop III stars.
We next see the results of parameter sets with the L-model stars, where the L-model stars experience efficient convective overshoot. The result for the L-std set (bottom-left) shows that the primary BH mass distribution for m 1 /M 50 is consistent with the observed one, while that for m 1 /M 50 is not. Even Pop III and EMP binary stars fail to form binary BHs within the PI mass gap. This result is consistent with Tanikawa et al. (2021b,a). In the L model, Pop III stars with ∼ 80 M expand in radius up to a few 10 3 R . Their H envelope is stripped through stable mass transfer or common envelope evolution. Thus, they cannot form PI mass gap BHs through binary evolution.
In contrast to the L-std set, the L-3σ set (bottomright) can produce binary BHs within the PI mass gap because of the different criteria for PPI and PISN effects. The primary BH mass distribution at z = 0 is consistent with the observation, as is already pointed out by Belczynski (2020). The primary BH mass distribution in the L-3σ set has no sudden drop at m 1 /M ∼ 50 unlike in the fiducial set. Recall that the sudden drop is caused by PPI effects for the fiducial set. From this point of view, the fiducial set might be slightly preferred over the L-3σ set.
The striking difference between the fiducial and L-3σ sets is seen in the redshift evolution of the BH merger rate density within the PI mass gap (see Figure 4 and the bottom-left panel of Figure 6). The merger rate density increases monotonically from z = 0 to z = 8 for the fiducial set, while it stops increasing at z = 4 for the L-3σ set. This can be confirmed in Figure 7. The merger rate density has a peak at a redshift of ∼ 11 in the fiducial set, and at a redshift of ∼ 6 in the L-3σ set. This is because Pop III stars dominantly form the PI mass gap events in the fiducial set, and do not in the L-3σ set. Note that Pop III SFR has a peak at a redshift of ∼ 20. There are two reasons why the Pop III stars have smaller contribution to PI mass gap events in the L-3σ set than in the fiducial set. First, Pop II stars can also form PI mass gap events. Second, the number of Pop III stars forming PI mass gap BHs becomes smaller in the L-3σ set than in the fiducial set, since more massive Pop III stars form PI mass gap BHs in the L-3σ set. In the fiducial set, Pop III stars with 65-90 M form the PI mass gap BHs, since they do not lose their H envelopes through binary interactions, and collapse to BHs without much mass loss. On the other hand, in the L-3σ set, Pop III stars with 130 M form the PI mass gap BHs. They expand to a few 10 3 R , lose their H envelopes through binary interactions, such as mass transfer or common envelope evolution, and become totally naked He stars with 65 M . Finally, they leave behind the PI mass gap BHs. Note that even stars with 65-90 M do not cause PPI nor PISN in the L-3σ set.
Surprisingly, we can find that Pop III binary stars have a significant impact on the formation of the PI mass gap events at high redshift even in the L-3σ set. As seen in the bottom panel of Figure 6, 1 of 2-3 PI mass gap events originates from Pop III binary stars at z 3. Even if stars experience efficient convective overshoot, Pop III binary stars can produce sufficient the PI mass gap events, depending on the PI model.
The primary BH mass distribution in the original-tide set is similar to that in the fiducial one. Pop III contributions become more important with increasing redshift, similarly to those in in the fiducial set. The prescription for tidal interaction has no impact on the primary BH mass distribution. Figure 8 shows the cumulative distribution of effective spins for parameter sets other than the fiducial one. Except for the original-tide set, the redshift evolution of the cumulative χ eff distributions is the same as for the fiducial set. This is because merging binary BHs are dominated by those with m 1 50 M for all the sets. These binary BH progenitors acquire their spin angular momenta in the same process for all the sets: Pop I binary stars acquire spins at the phase of naked He stars while all but Pop I binary stars are tidally spun up keeping their H envelope.
The cumulative distribution in the original-tide set is different from that in the fiducial set. In the original tidal prescription, naked He stars are less spun up. L-3 GWTC-2 z=8 z=4 z=2 z=0 Figure 6. BH merger rate density differentiated by primary BH mass at the redshifts of 0 (black), 2 (blue), 4 (green), and 8 (red) for parameter sets other than the fiducial one. The gray, cyan, light-green, and pink curves indicate Pop III contributions at the redshifts of 0, 2, 4, and 8, respectively. The shaded gray regions show the 90% credible interval inferred by the LIGO-Virgo Gravitational-Wave Transient Catalog 2 (GWTC-2).
Thus, binary BHs with χ eff ∼ 0.4 in the fiducial set are shifted to χ eff ∼ 0.2-0.4 in the original-tide set.
In the original tidal prescription, blue supergiant stars are more strongly spun up. Thus, binary BHs with χ eff 0.5 in the original-tide set contribute more than in the fiducial one. About 50 % of binary BHs have χ eff ∼ 0 in the fiducial set, while few binary BHs have χ eff ∼ 0 in the original-tide set. This is because blue supergiant stars are spun up more as described above.
In the original-tide set, the contribution of binary BHs with χ eff 0.5 is too large to be reconciled with the observation. The new tidal prescription seems to have advantages in reproducing the GW observations, which show the dominance of binary BHs with small χ eff in GWTC-2 (∼ 70-90 %) (Abbott et al. 2021a,b;Galaudage et al. 2021;Roulet et al. 2021).

DISCUSSION
As described in section 3.1, the BH merger rate density and its derivatives for the fiducial set are roughly consistent with the results of GWTC-2 and GWTC-2.1.
In this section, we compare the results of the fiducial set with those of previous studies.
First, we compare our results with two other isolated binary scenarios, Belczynski (2020) and Kinugawa et al. (2021c), who have claimed that isolated binary stars can reproduce the BH mass distribution in the range of 5-100 M inferred by the GW observations. Belczynski (2020) have reported that Pop I/II binary stars with metallicities 0.05-1.5 Z can make merging binary BHs consistent with the GW observations with regard to the BH mass and mass ratio distributions. Since we construct the 3σ PI model with reference to the model of Belczynski (2020), the results of the L-3σ set are naturally similar to those of Belczynski (2020). Here, we thus use the results of the L-3σ set as a substitute for those of Belczynski (2020) in discussing the differences between Belczynski (2020)'s and our model predictions. As described in section 3.2, the BH merger rate density with PI mass gap BHs evolve differently with redshift between the models: it stops increasing at z ∼ 6 in the L-3σ set, while it is still increasing up to z ∼ 11 in the fiducial set. This is because Pop III stars dominate the PI mass gap BHs in the fiducial set, and partially form them in the L-3σ set. With future GW observatories, such as Cosmic Explorer (Reitze et al. 2019) and Einstein Telescope (Punturo et al. 2010), which might be able to detect merging binary BHs in this mass range up to z ∼ 10, we can pick out the right scenario. Kinugawa et al. (2021c) have argued that merging binary BHs with m 1 30M are dominated by Pop III binary stars. On the other hand, not only Pop III but also Pop II and EMP binaries dominantly yield merging binary BHs with m 1 50M in our model. The future GW observatories can also distinguish these two models because of their capability to detect binary BH mergers at high redshifts. In the result of Kinugawa et al. (2021c), Pop III binary stars typically make binary BHs with m 1 ∼ 30M . Thus, the BH merger rate density with m 1 ∼ 30M should increase up to at least z ∼ 8. On the other hand, the BH merger rate density does not change from z = 0 to 8 in our model (see Figure 4).
In section 3.2, we show that Pop III and EMP binary BHs strongly depend on their IMFs, period distributions, and star evolution models. Here, we mention that Pop I/II binary BHs are also sensitive to several uncertainties. For example, the following two factors would make the local BH merger rate density much lower than the observed one. First, the actual common envelope parameters largely deviate from our choices, α = 1 and λ ∼ 0.1 (Dominik et al. 2012). Second, the metallicity distribution (σ Z ) is much smaller than our choice, σ Z = 0.35 . In those cases, we need to consider dynamical capture scenarios described below. In addition, the PI model (see Eq. (14)) can affect the BH mass distribution. Even if the 12 C(α, γ) 16 O reaction rate is fixed to the standard one, there are several possibilities of the PPI effects . The PPI effects can change the BH mass distribution around m 1 ∼ 40 M .
The observed primary BH mass distribution can also be explained by dynamical capture scenarios in globular clusters (e.g. Rodriguez et al. 2019), open clusters (e.g. Kumamoto et al. 2020;Di Carlo et al. 2020a;Santoliquido et al. 2020), or active galactic nuclei (AGN) disks (e.g Tagawa et al. 2020b). The difference between those models and ours is as follows. While binary BHs with m 1 ∼ 100 − 130M are in principle allowed to form in the dynamical capture scenarios, they are not in our scenario. Note that binary BHs with m 1 130M can form if 300 M stars are present (Mangiagli et al. 2019;Tanikawa et al. 2021b;Hijikawa et al. 2021), since they overcome PPI and PISNe. Thus, the primary BH mass distribution should smoothly extend beyond m 1 ∼ 100M if the dynamical capture is the dominant formation process of merging binary BHs. If our scenario is correct, the BH merger rate density should suddenly decrease at m 1 ∼ 100 − 130M . The fourth and fifth observing runs of the current GW observatories possibly assess the validity of these two scenarios if these observing runs detect binary BHs with m 1 ∼ 100 − 130M , or put some constraints on the upper limit of a BH merger rate in this m 1 range.

SUMMARY
We have derived the BH merger rate density and the distributions of merger properties by means of BPS calculations for Z/Z = 0-1 binary stars equivalent to Pop I/II/III and EMP binary stars, which span an unprecedentedly wide metallicity range. We adopted various initial conditions of binary stars and single and binary star models. With the fiducial set of the parameters, we successfully obtained the BH merger rate density, the primary BH mass in the range 5-100 M , and effective spin distributions consistent with the current GW observations. In our results, Pop III and EMP binary stars are responsible for the PI mass gap events like GW190521. Moreover, Pop III binary stars become more responsible for the PI mass gap events with increasing redshift, and contribute most of them at z 8. We have found that Pop III and EMP binary stars should have a topheavy IMF, small pericenter distance at the initial time, and inefficient convective overshoot in order to produce a sufficient number of the PI mass gap events to explain the PI mass gap event rate inferred by GW observations. We examine the dependence of binary BH properties on Pop III (and EMP) IMFs, period distributions, and star evolution models with different convective overshoot. We note that Pop III stars also contain other uncertain parameters, such as Pop III SFR and binary fraction, since Pop III stars have not yet been discovered. These uncertainties would directly affect our predicted BH merger rate density, especially the PI mass gap event rate mainly yielded by Pop III (and EMP) binary stars in our model. Apart from Pop III and EMP stars, Pop I/II binary BHs can also depend on initial conditions (SFR and metallicity distribution) and binary evolution parameters (common envelope parameters and mass transfer parameters). If they largely deviate from our choice, the local BH merger rate density can be much smaller (or larger) than the observed one. However, examining the dependence is beyond of this paper.
Finally, we emphasize the Pop III contributions to the formation of the PI mass gap events. In the parameter sets which reproduce the observed primary BH mass distribution (the fiducial, stepwise, original-tide, and L-3σ sets), Pop III binary stars have important roles in forming the PI mass gap events. In the fiducial, stepwise, and original-tide sets, Pop III binary stars dominate the PI mass gap at z 8. Even in the L-3σ set, 1 of 2-3 PI mass gap events originates from Pop III binary stars at z 3. This means that future GW observations of the PI mass gap events will be very useful for Pop III studies, if one of the above isolated binary scenarios (Belczynski 2020;Kinugawa et al. 2021c, and this paper) is correct.
The authors appreciate the anonymous referee for the thorough reading and many fruitful suggestions. The authors wish to express their cordial gratitude to Prof. Takahiro Tanaka, general PI of Innovative Area Grants-in-Aid for Scientific Research "Gravitational wave physics and astronomy: Genesis" for his continuous interest and encouragement. This research could not been accomplished without the support by Grantsin-Aid for Scientific Research (17H01101, 17H01102, 17H01130, 17H02869, 17H06360, 17K05380, 19K03907, 19H01934, 20H00158, 20H05249, 21K13914) from the Japan Society for the Promotion of Science. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org/), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. Software: Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo (Abbott et al. 2021c); HOSHI (Takahashi et al. 2016(Takahashi et al. , 2018(Takahashi et al. , 2019Yoshida et al. 2019); BSE (Hurley et al. 2000(Hurley et al. , 2002Tanikawa et al. 2020Tanikawa et al. , 2021a; Matplotlib (Hunter 2007); NumPy (van der Walt et al. 2011).

A. SINGLE STAR MODELS
In this section, we show the evolutions of the M-and L-model stars with various metallicities, and compare them with those of the SSE-model stars as a reference. We remark again that the M-and L-model stars have inefficient and efficient convective overshoot, respectively. Before showing their evolutions, we overview how to construct evolution tracks of the M-and L-model stars, which is described in detail in Tanikawa et al. (2020). We perform 1D hydrodynamics simulations of stars with various metallicities and masses. Based on the simulation results, we make fitting formulae of stellar evolution tracks, and incorporate the fitting formulae into the BSE code. The main components of the 1D hydrodynamics simulations are as follows. We use a 1D stellar evolution code, HOSHI code (Takahashi et al. 2016(Takahashi et al. , 2018(Takahashi et al. , 2019Yoshida et al. 2019). We take into account convection, semiconvection, and convective overshoot for chemical mixing, and do not account for rotation and rotational mixing. We input stellar physics with nuclear reaction network with 49 species of nuclei (Takahashi et al. 2018), stellar equation of state (Blinnikov et al. 1996;Vardya 1960;Iben 1963), OPAL, molecular, and conductive opacities (Iglesias & Rogers 1996;Ferguson et al. 2005;Cassisi et al. 2007, respectively), and neutrino energy loss (Itoh et al. 1996).
The M-and L-model stars are modeled as stars with the less and more efficient convective overshoot, respectively. In the HOSHI code, we treat convective overshoot as a diffusive process above convective regions. The diffusion coefficient of the convective overshoot exponentially decreases with the distance from the convective boundary as where D cv,0 and H P0 are the diffusion coefficient and the pressure scale height at the convective boundary, respectively, and ∆r is the distance from the convective boundary. The overshoot parameter f ov is set to be 0   Z/Z = 10 2 Figure 11. The same as Figure 9, except for the SSE model. Note that the SSE model does not support stars with Z/Z = 10 −4 , 10 −6 , and 10 −8 . Figure 11 shows the HR diagrams of the SSE-model stars as a reference. Note that the SSE model does not prepare stars with Z/Z < 5 × 10 −3 . They expand to more than 10 3 R for all masses and Z. These features are similar to the L-model stars.

B. COMPARISON BETWEEN THE SSE-AND M-MODEL STARS
In this section, we compare between the BH merger rate densities in the fiducial and SSE-std sets as a reference. Additionally, we also make comparison between the L-std and SSE-std sets. The difference between these sets is single star models with Pop II metallicities. Thus, we focus only on their results. The top panels of Figure 12 show the mass and spin distributions of merging binary BHs at the redshift of z = 0 in the fiducial and SSE-std sets. Roughly speaking, the mass and spin distributions in the SSE-std set are in good agreement with those in the fiducial set.
Especially, their spin distributions are quite similar. So, we focus only on their mass distributions. The primary BH masses in the SSE-std set are cut off at m 1 /M = 45 because of PPI and PISN effects, similarly to those in the fiducial set. The SSE-model stars expand to more than 10 3 R and evolve to red supergiant stars. They lose their H envelopes through stable mass transfer or common envelope evolution. Thus, they cannot form PI mass gap BHs.
The BH merger rate density in the SSE-std set is a few times smaller than in the fiducial set in the range of m 1 /M ∼ 30-40. The reason is as follows. The SSE-model stars have larger radii in their MS phases than the M-model stars. Thus, the SSE-model stars lose more masses in their MS phases through stable mass transfer. Then, they leave smaller He cores. In their post-MS phases, they further lose their masses through stable mass transfer, and finally leave naked He stars. Note that they have companion stars close enough to merge within the Hubble time. Thus, the resulting BHs in the SSE-std set have smaller mass than in the fiducial set, and the BH merger rate density with m 1 /M ∼ 30-40 is smaller.
In summary, the M-model and SSE-std stars generate consistent binary BHs, although there are several small differences. Moreover, when it comes to Pop II results, the M/L-model and SSE-std stars yield similar binary BHs.

C. EFFECTS OF THE METALLICITY CHOICE IN SSE-MODEL STARS
For the SSE-std set, we choose SSE-model stars with metallicities 0.25, 0.5, and 1Z for Pop I stars, SSE-model stars with metallicities 10 −2 , 0.025, 0.05, and 0.1Z for Pop II stars, and M-model stars with metallicities 10 −8 , 10 −6 , and 10 −4 Z for Pop III and EMP stars as seen in Tables 1 and 2. We do not use SSE-model stars with 5 × 10 −3 Z , despite that the SSE model supports metallicities down to 5 × 10 −3 Z . In order to examine this, we prepare the SSE-plus set. In the SSE-plus set, we add SSE-model stars with 5 × 10 −3 Z for Pop II stars in addition to the metallicity list of the SSE-std set. Figure 13 compares the mass and spin distribution of binary BHs in the SSE-std and SSE-plus sets. There is little difference between the two sets. Thus, we can conclude that the not using SSE-model stars with 5 × 10 −3 Z has little effect on our results.