Contribution of Population III Stars to Merging Binary Black Holes

A large number of mergers of binary black holes (BHs) have been discovered by gravitational wave observations since the first detection of gravitational waves 2015. Binary BH mergers are the loudest events in the universe, however their origin(s) have been under debate. There have been many suggestions for merging binary BHs. Isolated binary stars are one of the most promising origins. We have investigated the evolution of isolated binary stars ranging from zero metallicity (Population III stars or Pop III stars) to the solar metallicity by means of so-called rapid binary population synthesis simulation. We have found that binary BHs formed from isolated binary stars reproduce the redshift evolution of the merger rate density and the distribution of primary BH masses and mass ratios inferred by Gravitational-Wave Transient Catalog 3 (GWTC-3). Pop III stars have a crucial role in forming merging binary BHs in so-called the pair instability mass gap. Note that we choose the conventional prescription of pair instability mass loss, based on the standard $^{12}$C($\alpha$,$\gamma$)$^{16}$O reaction rate. Finally, we have shown the redshift evolution of the rate density of pair instability supernovae, and have predicted that a few pair instability supernovae would be discovered in the next few years. The discoveries would validate our results of merging binary BHs.


INTRODUCTION
The origin of merging binary black holes (BHs) is one of the biggest problems to solve in astronomy and astrophysics.In this paper, we suggest that isolated binary stars can be their origin when we consider metal-free or Population (Pop) III stars.In this section, we introduce three ingredients in this paper.In section 1.1, we overview gravitational wave observations of merging binary black BHs, and focus on binary BHs with ∼ 100 ⊙ in the so-called pair instability (PI) mass gap.In section 1.2, we describe pair instability supernovae (PISNe) forming PI mass gap.In section 1.3, we introduce Pop III stars, and briefly describe how Pop III stars produce binary BHs in the PI mass gap.Finally, we shortly summarize this section.

Merging binary black holes
Gravitational waves are ripples in spacetime generated by accelerated masses, and natural consequences of general relativity.Long after the presence of gravitational waves was predicted 1916 by Albert Einstein, it was indirectly confirmed by the orbital decay of the Hulse-Taylor binary pulsar (Hulse & Taylor 1975).About 100 years after the Einstein's prediction, gravitational waves have been directly detected by Laser Interferometer Gravitational-wave Observatory (LIGO), and the gravitational waves are generated by merging binary BHs (Abbott et al. 2016).This is also the first discovery of binary BHs and their mergers.Gravitational wave signals have been detected about 90 times until 2021 (Abbott et al. 2019(Abbott et al. , 2021a(Abbott et al. , 2023)), and most of them have been generated by merging binary BHs except for a few events (Abbott et al. 2017(Abbott et al. , 2020b(Abbott et al. , 2021b)).This means that merging binary BHs are the loudest objects in the LIGO band at a gravitational wave frequency of 10-100 Hz.
A large number of merging binary BHs enables us to assess the origin of merging binary BHs statistically, using BH masses, BH spins, positions, and redshifts.In particular, the mass distribution primary (or heavier) BHs is used most frequently, because primary BH masses are determined with an accuracy of ∼ 10 % (Abbott et al. 2019(Abbott et al. , 2021a(Abbott et al. , 2023)).Note that effective spins of merging binary BHs (or combinations of two BHs' spin components aligned to binary BHs' spin angular momenta) contain much larger relative errors than primary BH masses, since their typical values are zero (Abbott et al. 2019(Abbott et al. , 2021a(Abbott et al. , 2023)).Based on Gravitational-Wave Transient Catalog 3 (GWTC-3; The LIGO Scientific Collaboration et al. 2021), the primary BH mass distribution ranges from ∼ 5 ⊙ to ∼ 85 ⊙ with the global maximums of ∼ 10 ⊙ .The distribution does not decrease monotonically from ∼ 10 ⊙ to ∼ 85 ⊙ , and has an overdensity at ∼ 35 ⊙ and possible another overdensity at ∼ 17 ⊙ .
In this paper, we mainly focus on merging binary BHs with primary BHs of 65 ⊙ .Such merging binary BHs were thought not to be formed through isolated binary evolution because of PISNe as described in section 1.2.Thus, the discovery of GW190521 with ∼ 85 ⊙ primary BH (Abbott et al. 2020a,c) seemed to favor the formation of merging binary BHs in dense stellar clusters (Rodriguez et al. 2019;Di Carlo et al. 2020a;Tagawa et al. 2021b).After the discovery, several suggestions have been made that isolated binaries can also form such binary BHs (Belczynski 2020;Kinugawa et al. 2021a).We have also proposed that Pop III stars (see in section 1.3) can form GW190521-like binary BHs (Tanikawa et al. 2021a).In this paper, we review these results and related works.

Pair instability supernovae
Stellar evolution theory predicts that very massive stars (say 100-300 ⊙ at the zero-age main-sequence times, ZAMS times for short) cause thermonuclear explosions driven by PI, called "PISNe" (Fowler & Hoyle 1964;Barkat et al. 1967;Fraley 1968;Ober et al. 1983;Bond et al. 1984;El Eid & Langer 1986;Fryer et al. 2001;Heger & Woosley 2002;Umeda & Nomoto 2002;Langer 2012).A PISN proceeds as follows.A very massive star contains a carbon-oxygen (CO) core embedded in helium core and hydrogen envelope several Myr after its formation.The CO core has high temperature enough to cause pair production in which a photon is converted into an electropositron pair.The pair production decreases pressure in the CO core radiation-dominated, and then the CO core implodes.This implosion suddenly increases temperature in the CO core, and triggers explosive carbon and oxygen burning.This explosive burning finally overcomes the implosion, and conversely explodes the very massive star.This explosion leaves behind no stellar remnant.PISNe have not been conclusively discovered so far, although there are several promising candidates (Terreran et al. 2017;Schulze et al. 2023).
PISNe shape a mass range without BHs, so-called PI mass gap.This is because PISNe leave behind no stellar remnant, and stars with ZAMS masses of 100 ⊙ and 300 ⊙ leave behind BHs with quite different masses (Heger & Woosley 2002;Spera & Mapelli 2017).A 300 ⊙ star experiences PI, however fails to explode.The carbon and oxygen burning products are photodissociated into hydrogen and helium before explosion.Since this is endothermic reaction, the CO core pressure is decreased, and the whole star collapses to a BH.The BH mass is at least ∼ 130 ⊙ , which corresponds to helium core mass when PI sets in.A 100 ⊙ star experiences PI, and partially explodes because of weak carbon and oxygen burning.This explosion is so-called pulsational pair instability supernovae (PPISNe) (Heger & Woosley 2002;Woosley et al. 2007;Langer 2012;Chen et al. 2014;Yoshida et al. 2016;Woosley 2017;Marchant et al. 2019;Farmer et al. 2019;Leung et al. 2019;Renzo et al. 2020a).PPISNe leave behind ∼ 40 ⊙ BHs.Further lower-mass stars experience core-collapse supernovae, and leave behind 60 ⊙ BHs for Pop I/II stars.Thus, PI mass gap ranges from ∼ 60 ⊙ to ∼ 130 ⊙ for Pop I/II stars.
The PI mass gap can deviate from the above range, and is especially sensitive to the 12 C( , ) 16 O reaction rate (Takahashi 2018;Farmer et al. 2020;Costa et al. 2022).The 12 C( , ) 16 O reaction rate controls the ratio of carbon to oxygen at the ending time of the core helium burning phase; the carbon abundance becomes larger with the 12 C( , ) 16 O reaction rate decreasing.After the core helium burning phase, the higher carbon abundance makes carbon-carbon reactions more active.The more active carbon-carbon reactions effectively reduce helium core mass.Thus, the PI mass gap is shifted upward with smaller 12 C( , ) 16 O reaction rate.If the 12 C( , ) 16 O reaction rate is 3 (or 3 times) smaller than the standard value (Sallaska et al. 2013), the PI mass gap is shifted from ∼ 60-130 ⊙ to ∼ 90-180 ⊙ .In that case, GW190521 can be formed through isolated binary evolution.However, in this paper, we just adopt the standard value unless stated.

Population III stars
Pop III stars, also known as first stars or metal-free stars, are made of pristine gas at the high-redshift universe.Pop III stars have been theoretically studied.They are typically formed as massive stars from 10 ⊙ to 1000 ⊙ (Omukai & Nishi 1998;Abel et al. 2002;Bromm & Larson 2004;Yoshida et al. 2008;Hosokawa et al. 2011;Stacy et al. 2011Stacy et al. , 2012;;Bromm 2013;Susa 2013;Susa et al. 2014;Hirano et al. 2015;Sugimura et al. 2020).They can have low-mass, ∼ 1 ⊙ , as minor components (Nakamura & Umemura 2001;Machida et al. 2008;Clark et al. 2011b,a;Greif et al. 2011Greif et al. , 2012;;Machida & Doi 2013;Susa et al. 2014;Chiaki et al. 2016).Massive Pop III stars have not been directly observed so far (Rydberg et al. 2013), because they are located at the high-redshift universe due to their short-lived nature.Low-mass Pop III stars have been discovered neither (Frebel & Norris 2015;Magg et al. 2018Magg et al. , 2019)).Extremely metal-poor stars may be Pop III stars polluted by interstellar medium (Komiya et al. 2015).However, several studies have pointed out that Pop III stars cannot be polluted up to the level of extremely metal-poor stars (Tanaka et al. 2017;Suzuki 2018;Tanikawa et al. 2018;Kirihara et al. 2019).
It is under debate if Pop III binary stars contribute to merging binary BHs observed.Several studies have suggested that Pop III binary stars little dominate merging binary BHs (Hartwig et al. 2016;Belczynski et al. 2017).On the other hand, it has been argued that Pop III binary stars shape the overdensity at ∼ 35 ⊙ in the primary BH mass distribution (Kinugawa et al. 2021b).Pop III binary stars may dominate merging binary BHs in the PI mass gap (Tanikawa et al. 2021a(Tanikawa et al. , 2022b)), which is the main topic of this paper.
We briefly explain how Pop III binary stars fill the PI mass gap.Pop III stars are expected to lose little mass through stellar winds, since stellar winds becomes weaker with metallicity smaller (Vink & de Koter 2005).A Pop III star with 60-100 ⊙ can keep its mass until it ends its life.It forms a CO core not large enough to cause PI, and then produces a iron core.The iron core collapses at some point.The Pop III star leaves behind 60-100 ⊙ BHs without mass ejection.Then, we can find that Pop III stars can fill the lower side of the PI mass gap.In section 3.1, we will assess if such BHs can be genuinely members of merging binary BHs.

Short summary
In the above sections, we overview merging binary BHs, in particular, those with BHs in the PI mass gap.Such merging binary BHs are suggested to be possibly formed from Pop III binary stars.In sections 2 and 3, we describe the method to confirm the suggestion, and its results, respectively.In section 4, we summarize this paper.

BINARY POPULATION SYNTHESIS SIMULATION
We describe the method to investigate isolated binary evolution.In section 2.1, we show initial conditions we adopt.In section 2.2, we overview single and binary evolution models.

Initial conditions
We take into account the redshift evolutions of the star formation rate (SFR) densities for Pop I/II and Pop III stars separately.Figure 1 shows Pop I/II and Pop III SFR density evolution as a function of redshift.We adopt the formula of Madau & Fragos (2017) for the Pop I/II SFR density evolution.Note that the Pop I/II SFR density is updated by recent observations, especially for high redshift (Harikane et al. 2022(Harikane et al. , 2023)).On the other hand, we simplify the simulation result of Skinner & Wise (2020) for the Pop III SFR density evolution.
Here, we note Pop III SFR density evolution we adopt.The simulation of Skinner & Wise (2020) has considered Pop III star formation and feedback, and followed the ending of Pop III star formation due to metal enrichment.We construct our Pop III SFR density evolution to approximate the Pop III SFR density evolution in fig. 4 of Skinner & Wise (2020) by connecting three line segments shown in eq. ( 5) in Tanikawa et al. (2022b).We use a Pop III SFR density evolution based on a theoretical study, because Pop III stars have not yet been observed unlike Pop I/II stars.We note that other studies have also suggested Pop III SFR density evolution (Jaacks et al. 2019;Liu & Bromm 2020a;Visbal et al. 2020;Hartwig et al. 2022), and some studies have constrained the upper limit of Pop III SFR density evolution from observations (de Souza et al. 2011;Inayoshi et al. 2016).Our model may stop Pop III star formation earlier than other models, for example, suggested by Visbal et al. (2020) who have shown that reionization makes Pop III star formation long-lived.
Figure 2 shows the average metallicity of Pop I/II stars formed at each redshift.For each redshift, stellar metallicities are distributed as a logarithmic normal distribution centered on the average metallicity with dispersion of 0.35 dex in logarithmic scale.We base the average metallicity evolution on the analysis of observational data by Madau & Fragos (2017), and simplify the metallicity dispersion referring to the simulation results of Chruslinska & Nelemans (2019).Figure 2 also indicates metallicities lower and higher than the average one by 1 for each redshift.For Pop III stars, we assume that all of them have zero metallicity.
The initial conditions of single and binary stars are as follows.We assume that the binary fraction is 0.5 in number for all metallicities.This binary fraction is similar to an intrinsic binary fraction obtained by Sana et al. (2012).The initial mass function of single stars and primary stars is set as a combination of two mass functions.The first mass function is Kroupa's IMF (Kroupa 2001).The minimum and maximum masses are 0.08 and 150 ⊙ , respectively.The second one is a function flat in the logarithmic scale (hereafter, logarithmically flat IMF) from 10 ⊙ to 150 ⊙ , which is modeled as Pop III IMF (Hirano et al. 2014;Susa et al. 2014;Hirano et al. 2015).We change weights of the two IMFs depending on metallicities, which is motivated by star formation simulations with various metallicities (Chon et al. 2021(Chon et al. , 2022)).We show the weights in Figure 3.The weights of Kroupa's IMF are unity and zero for metallicity of ≥ 0.1 ⊙ and ≤ 10 −6 ⊙ , respectively, where we assume ⊙ = 0.02.Between these metallicities, the weights decreases gradually.We  adopt Sana's model (Sana et al. 2012) for binary initial conditions: the distributions of mass ratios, periods, and orbital eccentricities.
Similarly to our Pop III SFR density evolution, we need to rely on theoretical studies for determining Pop III IMF.We construct our Pop III IMF, referring to fig. 9 in Hirano et al. (2014) and fig.6 in Hirano et al. (2015).Our logarithmically flat IMF does not deviate from their results by an order of magnitude.

Single and binary evolution model
We investigate the formation of merging binary BHs through isolated binary evolution, employing rapid binary population synthesis simulation with the BSEEMP code (Tanikawa et al. 2020(Tanikawa et al. , 2022b)), which is based on the BSE code (Hurley et al. 2000(Hurley et al. , 2002)).Note that the BSE code simultaneously solves single star evolution from ZAMS stars to stellar remnants (e.g.white dwarfs, neutron stars, and BHs) with stellar wind mass loss, and binary star evolution, such as wind accretion, tidal interaction, stellar mass transfer, common envelope evolution, and orbital decay due to gravitational wave radiation.The BSEEMP code adopts different single and binary evolution model from the BSE code.Here, we overview the single and binary evolution model adopted by the BSEEMP code.In particular, we concentrate on different parts between the BSE and BSEEMP codes.
In the BSEEMP code, single stars evolve along with fitting formulae of their total masses, helium core masses, radii, and luminocities.Note that the BSE-based code like the BSEEMP code can rapidly follow single and binary star evolution, since such fitting formulae are prepared in advance.In the BSEEMP code, single stars with metallicities of > 0.1 ⊙ evolve the same as in the BSE code.Note that the single star evolution model has been constructed by Pols et al. (1998).On the other hand, single stars with metallicities of ≤ 0.1 ⊙ evolve along with different fitting formulae from the BSE code.Additionally, as one of the important features of the BSEEMP code, the BSEEMP code supports ≤ 5 × 10 −3 ⊙ , while the BSE code does not.The BSEEMP code prepares two types of fitting formulae called "M model" and "L model" constructed in Yoshida et al. (2019) with the HOSHI code (Takahashi et al. 2016(Takahashi et al. , 2018(Takahashi et al. , 2019)).The "M" and "L" stand for "M"ilky Way and "L"arge Magellanic Cloud, respectively.Yoshida et al. (2019) have referred stars in these galaxies when they construct these models.The difference between the M and L models is the efficiency of overshoot at the convective boundary, and the overshoot in the M model is smaller than in the L model as summarized in Appendix A of Tanikawa et al. (2022b).For extreme metal-poor cases, radius evolutions of M-and L-model stars are quite different.Figure 4 shows the Hertzsprung-Russell diagram of M-and L-model stars with = 10 −8 ⊙ equivalent to Pop III stars.M-model stars expand up to ∼ 100 ⊙ even for a 160 ⊙ star, while L-model stars exceed 10 3 ⊙ for 80 and 160 ⊙ stars.If we switch off the overshooting effect, the maximum radii should be smaller than those of M-model stars.
The reason for the difference of the maximum radii between M-and L-models is as follows.A L-model star's core gains fresh hydrogen from its envelope due to more efficient overshoot during its mainsequence (MS) phase.Eventually, a L-model star has a larger helium core mass than a M-model star at the ending of their MS phases by 20 %.As a result of this mass difference, a L-model star has much larger core luminosity than a M-model star at their post-MS (PMS) phases.Since larger core luminosity makes stellar radius larger, a L-model star has a much larger maximum radius than a M-model star.
Here, we compare the M-and L-models with other models.The radius evolution of a L-model star is similar to the Pop III model by Marigo et al. (2001).Stars with 80 ⊙ finally evolve to red supergiant stars in both the Marigo's and L models.On the other hand, M-model stars evolve like the Pop III model by Farrell et al. (2021).Even stars with ∼ 80 ⊙ keep blue supergiant stars throughout their lives in both the Farrell's and M models.
Here, we describe how to separate four stellar phases: MS, Hertzsprung-gap, core helium burning, and shell helium burning phases.When stars run out of hydrogen at their core, the whole stars shrink, and their effective temperature once become larger.We define the ending of the MS phase as the time when the effective temperature gets the local maximum, which can be seen at the effective temperature of ∼ 10 4.7 K in Figure 4.The ending of Hertzsprunggap, core helium burning, and shell helium burning phases can not be seen in Hertzsprung-Russell diagram like Figure 4. We read required information from the results of the HOSHI code when we construct our fitting formulae.We define the ending of the Hertzsprung-gap phase as the time when helium is ignited at the stellar core, or the central temperature exceeds ∼ 2 × 10 8 K.We define the ending of the core helium burning phase as the time when the central helium fraction becomes sufficiently small.We set the fraction to 0.1 %.We define the ending of the helium shell burning phase as the time when carbon is ignited at the stellar core, or the central temperature exceeds ∼ 10 9 K.Note that Pop III stars have no Hertzsprung phase unlike Pop I/II stars.They have enough high central temperature to ignite helium at the ending of their MS phases.Thus, Pop III stars skip Hertzsprung-gap phases in the BSEEMP code.
Massive stars ( 8 ⊙ ) finally experience supernovae, and leave behind neutron stars, BHs, or no stellar remnant.Supernova model in the BSEEMP code is the Fryer's rapid model (Fryer et al. 2012) with the modification of PI mass loss modeled as the moderate PPISN model (Leung et al. 2019;Belczynski & Banerjee 2020).Note that several PPISN models have been suggested (Stevenson et al. 2019;Renzo et al. 2020b).
We evolve 20-160 ⊙ stars from their ZAMS phases to their remnant phases, taking into account stellar wind mass loss and supernova mass loss.We can see the effects of stellar winds and supernovae in Figure 5. Remnant masses increase with metallicities decreasing because of weakening stellar winds.Pop III remnant masses are suddenly reduced at ZAMS masses of ∼ 80-85 ⊙ and ∼ 120-135 ⊙ due to PPISN and PISN effects, respectively.Note that helium core masses just before core collapse exceed ∼ 40 ⊙ and ∼ 60 ⊙ at ZAMS masses of ∼ 80-85 ⊙ and ∼ 120-135 ⊙ , respectively.PPISN and PISN can equalize remnant masses of = 10 −2 ⊙ and Pop III stars.Especially for the M model, Pop III stars with ZAMS masses of 60 -85 ⊙ can fill the lower side of the PI mass gap.This is because they keep their initial mass before core collapse, and do not experience PPISN nor PISN.
The ZAMS-remnant mass relations are similar between the M and L models, however slightly different.ZAMS mass ranges for PPISNe and PISNe are shifted downward for the L model.This is because L-model stars create larger helium core masses due to more efficient overshoot than M-model stars during their MS phases.This effect is noticeable in the difference of the maximum BH mass between the L and M models (∼ 60 ⊙ and ∼ 85 ⊙ , respectively).
Neutron stars and BHs should receive natal kicks due to asymmet-ric supernova explosions following core collapses.We assume this natal kicks, such that the velocity distribution is given by a single Maxwellian with 265 km s −1 (Hobbs et al. 2005) if there is no fallback mass.We reduce the natal kick velocity by 1 − b where b is the fraction of the fallback mass (Fryer et al. 2012).We suppose that the directions of the natal kicks are isotropic.We construct our binary evolution model, the BSEEMP model, based on the BSE model (Hurley et al. 2002) with several modifications.The BSE model considers wind accretion, tidal evolution, stable mass transfer, common envelope evolution, magnetic braking, and orbital decay through gravitation wave radiation.We modify prescriptions of tidal evolution, stable mass transfer, and common envelope evolution.Here, we focus on different points of these processes from the BSE model.
A star of a binary star evolves to a giant star at some point, and fills its Roche lobe unless the orbit of the binary star is too a wide.The star (hereafter, donor) transfers its mass toward its companion star (hereafter accretor).If the mass transfer is stable, the binary star experiences stable mass transfer.Otherwise, the binary star experiences stellar coalescence, or common envelope evolution, depending on their stellar types as described in detail below.A donor can sustain stable mass transfer if it has a radiative envelope, whereas a donor with a convective envelope can cause unstable transfer.This is because a star with a radiative envelope more sensitively decreases its radius with decreasing its mass.We modify the boundary criteria between radiative and convective envelopes.In the BSE model, a massive star changes its envelope from radiative to convective when it finishes core helium burning, and starts shell helium burning.On the other hand, in the BSEEMP model, a massive star has a radiative (convective) envelope when the effective temperature of its envelope is above (below) 10 3.65 K.In general, whether an envelope is radiative or convective depends on the effective temperature.Moreover, in the M and L models, a massive star can start shell helium burning at an effective temperature much higher than 10 3.65 K, especially for metal-poor stars, such as Pop III stars.Thus, this criteria are very reasonable.Similar criteria have been suggested (Kinugawa et al. 2014;Olejak et al. 2021).
We also modify the prescription of stable mass transfer.In the BSE model, all the mass lost by a donor are transferred onto an accretor unless the mass accretion rate exceeds Eddington accretion limit.In the BSEEMP model, we can control the fraction of transferred mass.Here, we fiducially set the fraction to 0.5.This choice can be seen in several rapid binary population synthesis codes (Belczynski et al. 2020;Kinugawa et al. 2020).Since the fraction can be not wellestablished, and range from zero to unity (Meurs & van den Heuvel 1989), we take the middle.This is consistent with a recent work that the fraction should be 0.5 (Stevenson et al. 2017).
In the BSE model, if post main-sequence stars including Hertzsprung-gap, core helium burning, and shell helium burning stars fills its Roche lobe, and its mass transfer is unstable, a binary star experiences common envelope evolution.However, in the BSEEMP model, this is true only for core helium burning and shell helium burning stars.If a Hertzsprung-gap star fills its Roche lobe, and its mass transfer is unstable, the binary star experience stellar coalescence.This is because a Hertzsprung-gap star does not yet develop steep density gradient between the helium core and hydrogen envelope (Ivanova & Taam 2004).Such a prescription is also equipped with several rapid binary population synthesis codes (Dominik et al. 2012;Giacobbo et al. 2018).
Pop III stars, even at the ZAMS time, the core region contracts sufficiently so that the central temperature reaches ∼ 10 8 K, where carbon is produced by helium burning.Thus, hydrogen burning stably proceeds via the CNO cycle.Nevertheless, Pop III stars do not have steep density gradient between their core and hydrogen envelope in their MS phases.This is because Pop III stars contracts not only core regions but also overall regions at the ZAMS time.
In common envelope evolution, a binary orbit shrinks as a back reaction of common envelope ejection.In rapid binary population synthesis codes, such a common envelope evolution is modeled as so-called the formalism (Webbink 1984).In this formalism, two parameters and have to be set.A parameter controls the efficiency at which a binary orbital binding energy is converted to energy driving common envelope ejection through gas drag.We set = 1.A parameter is a dimensionless parameter of the binding energy of a donor's envelope or common envelope.We adopt the formulae of Claeys et al. (2014) for .
We adopt a new prescription for tidal evolution of stars with radiative envelopes.The tidal evolution follows the radiative damping of the dynamical tide (Zahn 1975).In the BSE model, the tidal coupling parameter is given by Zahn (1977) and Hut (1981).We replace this tidal coupling parameter with a new formula derived by Yoon et al. (2010) and Qin et al. (2018).This formula is also used by Kinugawa et al. (2020).

Realization
We prepare two parameter sets which adopt the M and L models.We consider a parameter set with the M model as the fiducial model, and call it "fid" set.We call the other one the "L" set.
We generate primary stars of binary stars, such that the minimum and maximum masses are 10 and 150 ⊙ .Moreover, we set the lower limit of the secondary stars to be 10 ⊙ .This is because binary stars with < 20 ⊙ in total can form no BH even if the binary components merge.

RESULTS
In this section, we show several features of our fid set.In sections 3.1 and 3.2, we present properties of merging binary BH and PISNe, respectively.

Population III binary black holes
In this section, we compare our theoretical results with observed binary BH population.First of all, we can see the redshift evolution of binary BH merger rate densities for the fid and L sets in Figure 6.For both the fid and L sets, the merger rate densities are ∼ 20 yr −1 Gpc −3 at the present day.The merger rate densities increase up to redshifts of 2 or 3, and decrease beyond redshifts of 2 or 3.The features follow the redshift evolution of SFR densities (see Figure 1).Note that the merger rate densities are less sensitive to redshifts than SFR densities.This is because binary BH mergers can be delayed by several Gyr from star formation.
For the fid and L sets, the merger rate densities are quite consistent with the GWTC-3 result.The merger rate densities are within 90 % credible bounds inferred by observations.Note that we obtain this consistency because of prohibition of common envelope evolution for Hertzsprung-gap stars (see section 2.2).Unless we prohibit the common envelope evolution of Hertzsprung-gap stars, we would obtain ∼ 10 times larger BH merger rate density than obtained in Figure 6.The merger rate densities are very sensitive to binary evolution model.
We explain why the criteria of common envelope evolution affects the merger rate density in more detail.It may be unclear why prohibition of common envelope evolution for Hertzsprung-gap stars largely affects the merger rate density, because Hertzsprung-gap stars cannot enter into common envelope evolution due to their radiative envelopes.Certainly, Hertzsprung-gap stars are not easy to enter into common envelope evolution, but not impossible.If Hertzsprunggap stars have sufficiently larger masses than their companion stars, they can enter into common envelop evolution.Furthermore, it may look strange why M-model stars enter into common envelope evolution, because Pop III stars always have radiative envelopes.Actually, Pop II stars in the M model have convective envelopes just before core collapse.Pop III stars in the M model can always have radiative envelopes, not only because of their inefficient overshoot, but also the small opacities of their metal-free envelopes.We can see Hertzsprung-Russell diagrams of Pop II stars in the M model from fig. 9 of Tanikawa et al. (2022b).
At the present day, Pop III binary BHs contribute to only 1 and 3 % of the total merger rate densities for the fid and L sets, respectively.They are always minor for the L set, while they become dominant beyond a redshift of 10 for the fid set.The current observatories detect binary BHs with non Pop III origins, although the future observatories, such as Einstein Telescope (Punturo et al. 2010) and Cosmic Explorer (Reitze et al. 2019), will discover many binary BHs originating from Pop III stars.
Figure 7 shows the primary BH mass distributions of merging binary BHs at the present day, where a primary BH mass (or 1 ) is the heavier BH mass of two BHs.For the fid set, the primary BH mass ranges from 5 ⊙ to 85 ⊙ , similarly to the GWTC-3 result.Moreover, Pop III binary BHs dominate binary BHs with 1 50 ⊙ .These BHs come from Pop III stars with their ZAMS mass of 50-90 ⊙ (see Figure 5).In other words, Pop III binary BHs fill the lower half of the PI mass gap (60 -130 ⊙ ).This feature is quite different from that of the L set.For the L set, there is no binary BHs with 1 55 ⊙ .Even Pop III binary BHs cannot fill the PI mass gap at all.As seen in Figure 7, there is overproduction at 1 ∼ 20 ⊙ for both the fid and L sets.This comes from binary BHs with Pop II origins.We may solve this overproduction, constructing Pop I/II SFR density evolution more carefully (van Son et al. 2023).
In Figure 8, we draw the picture of a typical formation channel of a merging binary BH in the PI mass gap for the fid set.We prepare a Pop III binary star with the primary mass of ∼ 85 ⊙ , the secondary mass of ∼ 31 ⊙ , and the separation of ∼ 30 ⊙ .The primary star evolves with little mass loss.Pop III stars experience weak stellar wind mass loss.Moreover, for the fid set, a ∼ 85 ⊙ star expands up to a few 10 ⊙ (see Figure 4), and loses little mass through stable mass transfer.Eventually, the primary star leaves behind a BH with ∼ 80 ⊙ in the PI mass gap.The secondary star also experiences little mass loss through stellar winds and stable mass transfer.Because it experiences supernova mass loss, it leaves behind a BH with ∼ 20 ⊙ .Finally, they evolve to a binary BH with a wide but eccentric orbit.Thus, they can merge within the Hubble time.
As seen in Figure 8, progenitors of binary BHs in the PI mass gap experience stable mass transfer.Nevertheless, their formation channel does not depend on the choice of formulation of stable mass transfer, such as the mass-transfer termination conditions (Inayoshi et al. 2017).This is because Pop III stars in the M model experience little stable mass transfer due to their small maximum radii.
In Figure 8, we can see that the primary star loses 0.07 ⊙ from = 0 Myr to 2.37 Myr, while the secondary star gains only 0.01 ⊙ during this time, despite that the fraction of transferred mass is 0.5.This is because the primary star its mass mainly due to stellar wind mass loss, although it also experiences mass transfer.This is the reason why the secondary star increases only 0.01 ⊙ .
For the L set, the channel in Figure 8 is prohibited.Since the primary star expands beyond 10 3 ⊙ , it loses a large amount of mass in its hydrogen envelope through stable mass transfer or common envelope.Then, they leave behind ∼ 40 ⊙ , not BHs in the PI mass gap.Instead of the above initial condition, let's consider a binary star with the separation of 10 3 ⊙ .In that case, they do not experience mass loss through binary interactions.Thus, they can leave ∼ 80 ⊙ BHs.However, their separation is still 10 3 ⊙ .They cannot merge within the Hubble time.This is the reason why there are no merging BHs in 1 55 ⊙ .
From the above, we can conclude that Pop III stars can explain the presence of merging binary BHs in the PI mass gap if their convective overshoot is as inefficient as the M model.In the future, we need whether the convective overshoot is efficient or not.We mention other three features of the primary BH mass distribution.First, for both the fid and L sets, the binary BH population suddenly decreases below 1 ∼ 8 ⊙ .This is because the supernova model we adopt tend not to form BHs with 8 ⊙ .We note that this sudden decrease can be explained not by the supernova model but by binary interaction (van Son et al. 2022b).Second, there is the overdensity at ∼ 35 ⊙ in the GWTC-3 result.However, our results cannot form the overdensity.Although PPISNe are expected to form the overdensity they fail to form (Renzo et 2022;van Son et al. 2022a).We may not fully explain the primary BH mass distribution, although we can form merging binary BHs in the PI mass gap.Third, there may be another overdensity at ∼ 17 ⊙ although the overdensity cannot be seen in the GWTC-3 result of Figure 7.It seems that the fid set has an overdensity at ∼ 17 ⊙ .We have no intention of claiming that the overdensity in the fid set is consistent with the overdensity at ∼ 17 ⊙ in the GWTC-3 result.
We can see the mass ratio distribution of merging binary BHs at the present day in Figure 9.The mass ratio is defined as the ratio of the secondary BH mass ( 2 ) to the primary BH mass, where the secondary BH mass is the lighter BH mass of two BHs.As seen in Figure 9, the mass ratio distributions in our results are within 90 % credible bounds of the GWTC-3.Especially for the fid set, binary BHs with low mass ratios (∼ 0.2) are formed from Pop III binary stars with low mass ratios at the initial time.As described above, Pop III binary stars leave behind merging binary BHs with little mass loss through binary interaction.They can keep the initial mass ratio.If Pop III binary stars have a mass ratio distribution similar to Pop I/II binary stars, they can form merging binary BHs with such low mass ratios.
Figure 10 shows the cumulative distribution of effective spins of merging binary BHs at the present day for the fid and L sets.An effective spin can be expressed as where − → 1 and − → 2 are, respectively, dimensionless spin vectors of the primary and secondary BHs, and − → is the orbital angular momentum of a binary BH.Roughly speaking, an effective spin indicates BH spin component parallel to the orbital angular momentum of a binary BH.For both the fid and L sets, most of merging binary BHs have positive eff , about 80 % of them have eff 0.1, and the rest of them have eff 0.1.
Merging binary BHs with eff 0.1 have highly spinning secondary BHs.Such BHs are mainly formed as follows (Kushnir et al. 2016;Hotokezaka & Piran 2017).Let's consider a binary with one BH and one main-sequence star, where the BH mass is smaller than the main-sequence star.Such a binary can appear after the primary star collapses to the BH.When the main-sequence star evolves off to a giant star, it fills its Roche lobe, and then starts mass transfer.Regardless of the mass transfer stability, the binary separation becomes smaller with the giant star losing its mass.The giant star completely loses its hydrogen envelope, and evolves to a naked helium star.The binary interaction shrinks the binary separation down to ∼ 10 ⊙ .Then, the naked helium star is tidally spun-up.Finally, the naked helium star collapses to a highly spinning BH.When the binary BH merges, we can observe it as a binary BH with high eff .
As seen in Figure 10, about 80 % of merging binary BHs have eff 0.1.They have large separation (say 10 ⊙ ) when the giant star evolves to the naked helium star, losing its hydrogen envelope completely.The naked helium star cannot be tidally spun-up when the binary separation is large.Eventually, binary BHs have larger effective spins when their progenitors have smaller separation after common envelope evolution.Initial separations of binary BH progenitors indirectly affect effective spins of binary BHs, since binary BH progenitors have smaller separations after common envelope evolution when their initial separations are smaller.
As described above, tidal interaction between naked a helium star and a BH is the primary source of high spins.This is the case for Pop I/II binary BHs.In the case of Pop III binary BHs, the primary source of high spins is stable mass transfer.As seen in Figure 10, a large fraction of Pop III binary BHs have eff 0.4 due to spin up in stable mass transfer.This can be also seen in Figure 8.However, we emphasize that this spin-up process is not dominant for overall binary BHs, since Pop III binary BHs are minor compared to Pop I/II binary BHs.
As seen in Figure 6, Pop III binary BHs have small contribution.Thus, their effective spin distribution does not affect the overall distribution.However, we explain how to form their effective spin distribution especially for the fid set, because their distribution is completely different from the overall distribution.About 40 % of Pop III binary BHs have eff 0.1.The reason is as follows.They can keep their hydrogen envelope without stable mass transfer nor common envelope evolution because of their small radii (see Figure 4).Thus, they can keep their spin angular momenta throughout their evolution.Finally, they can collapse to highly spinning BHs, and merge as binary BHs with high eff .
For the L set, Pop III binary BHs have an effective spin distribution similar to the overall distribution.Pop III stars expand similarly to Pop I/II stars.Thus Pop III BHs obtain spin angular momenta in the same process as Pop I/II BHs: tidal spin-up during naked helium stars as described above.
For both the fid and L sets, the cumulative distributions appear different from the GWTC-3 result.In particular, our results contain no merging binary BHs with negative eff in contrast to the GWTC-3 result.However, this may not be much of a problem.Several studies have suggested that binary BH population with positive eff is real, but BH population with negative eff is spurious (Galaudage et al. 2021;Roulet et al. 2021).If this is true, our results are consistent with the GWTC-3 result.

Pair instability supernova rate
We investigate features of PISNe for the fid set in order to forecast current and future PISN surveys (Pan et al. 2012;Moriya et al. 2019;Wong et al. 2019;Regős et al. 2020;Moriya et al. 2022b,a).As reference, we prepare a different parameter set which can produce merging binary BHs in the PI mass gap despite of choice for the L model.In this section, we first introduce the parameter set, called the "L-3 " set.
For the L-3 set, we adopt the L model for the single star evolution model, and the Fryer's rapid model for the supernova model, which is the same as the L set.The L-3 set has a different PI mass loss from the L set.In the PI mass loss model, no PPISNe happen, and PISNe are generated from stars with helium cores of 90 -180 ⊙ .PISN progenitor masses are shifted upward because of reducing the 12 C( , ) 16 O reaction rate by 3 or 3 times.This is based on the results of Farmer et al. (2020) and Costa et al. (2021), and similar to the PI mass loss model adopted by Belczynski (2020).
Figure 11 shows the relation between ZAMS and remnant masses for the L-3 set.Pop III and = 10 −2 ⊙ stars cause PISNe for ZAMS masses of 150 ⊙ and 250 ⊙ , respectively.
= 10 −2 ⊙ stars need to have larger ZAMS masses than Pop III stars, because they largely lose their masses through their stellar winds.= ⊙ stars can not generate PISNe.They cannot keep their masses large enough to cause PISNe just before their core collapses.
Pop III stars with ZAMS masses of ∼ 150 ⊙ leave behind ∼ 125 ⊙ BHs.This is because they have their hydrogen envelopes until their core collapses.For Pop III stars, PISN progenitors can have hydrogen envelopes.On the other hand, = 10 −2 ⊙ stars with ZAMS masses of ∼ 250 ⊙ leave behind ∼ 90 ⊙ BHs.They completely lose their hydrogen envelopes until their core collapses.
= 10 −2 ⊙ PISN progenitors are naked helium stars.Actually, this is also true for the fid set.Observational features of PISNe depend on whether hydrogen envelopes are present or not (Kasen et al. 2011).
Figure 12 shows the primary BH mass distribution of merging binary BHs at the present day for the L-3 set.We can see that the L-3 set also produces merging binary BHs in the PI mass gap.Such merging binary BHs can be yielded by both of Pop III and non-Pop III binary stars.This reason is as follows.A Pop III star with 150 ⊙ in a binary star loses its hydrogen envelope through binary interaction and leaves behind a naked helium star with 90 ⊙ , because its maximum radius exceeds 10 3 ⊙ (see Figure 4).This naked helium star can collapse to a BH with 90 ⊙ .Note that PISNe only happen for stars with helium core masses of 90 ⊙ .
We investigate different observational features between the fid and L-3 sets both of which can produce merging binary BHs in the PI mass gap (Tanikawa et al. 2023).We expect different PISN features, because these sets have different PI mass loss models.
Figure 13 shows the redshift evolution of PISN rate densities for the fid and L-3 sets.Note that we include not only binary stars but also single stars.For both the sets, PISN rate densities increase with redshifts increasing up to a redshift of ∼ 12.For the higher-redshift universe, more metal-poor stars are formed, and metal-poor stars cause PISNe more easily because of weaker stellar winds.Up to a redshift of ∼ 6, PISN progenitors are dominated by naked helium stars.This is because metal-rich stars (say 10 −2 ⊙ ) lose their hydrogen envelopes until their core collapses.Beyond a redshift of ∼ 6, PISN progenitors can have hydrogen envelopes, and a part of PISN progenitors are Pop III stars.
For the L-3 set, PISNe happen in spite of the fact that even Pop III PISNe happen only when their ZAMS masses are larger than 150 ⊙ (see Figure 11).Note that we set the maximum ZAMS mass to 150 ⊙ as described in section 2.1.For the L-3 set, merged stars and stars accreting mass from their companion stars cause PISNe.This is partly because the PISN rate density is small.However, we make clear that the PISN rate density increases only slightly even if we increase the maximum ZAMS mass to 300 ⊙ (Tanikawa et al. 2023).Star cluster observations find at most ∼ 300 ⊙ stars (Doran et al. 2013;Schneider et al. 2018b,a).300 ⊙ ZAMS stars would be large enough to be regarded as the maximum ZAMS mass.
We focus on the detectability of Euclid space telescope (Laureĳs et al. 2011), because Euclid space telescope has a PISN survey program (Moriya et al. 2022a).Euclid space telescope can detect PISNe up to a redshift of ∼ 4. The PISN rate density of the fid set is larger than that of the L-3 set by 2 orders of magnitude.Actually, this difference has a large impact on the PISN detectability by Euclid space telescope.For the fid and L-3 sets, Euclid space telescope can detect two and zero PISNe, respectively, for 5-year operation (Tanikawa et al. 2023).Thus, if Euclid space telescope will detect PISNe, the fid set may be correct.
There was one concerning point for the above argument.This is because there was no reference which describes observational features of PISNe with the 3 -smaller 12 C( , ) 16 O reaction rate.Recently, the amount of Nickel-56 mass yielded from PISNe with various 12 C( , ) 16 O reaction rates have been calculated (Kawashimo et al. 2023).Nickel-56 plays an important role in observational features of supernovae, because its radioactive decay powers supernova luminosity.According to Kawashimo et al. (2023), Nickel-56 mass decreases with the 12 C( , ) 16 O reaction rate decreasing if ZAMS mass is fixed.This means that PISNe in the L-3 set have not only a smaller rate density but also smaller luminocities than those in the fid set.Thus, we can confirm the above argument that no PISNe can be detected by Euclid space telescope.
We have to make caveats on the above arguments.First, stellar wind mass loss contains uncertainties, especially for massive stars causing PISNe (Gräfener 2021;Vink et al. 2021).This can change ZAMS masses generating PISNe.Thus, the detection number of PISNe for the fid set is changeable, although it is robust that the detection number of PISNe for the fid set is much larger than for the L-3 set.Second, as described above, a part of PISNe are yielded from merged stars and stars accreting masses from their companion stars.In the BSEEMP model, merged stars and stars accreting masses from their companion stars evolve in the same way as stars with the same mass which do not experience mergers nor mass accretion.In fact, merger and accretion history can alter stellar evolution.This is because these effects cause efficient mixing of stellar interior, providing fresh fuel (Hellings 1983(Hellings , 1984;;Renzo & Götberg 2021).These effects can enlarge helium core masses, and finally affect PISN feasibility.Actually, our calculations generate a large fraction of PISNe which do not happen without mergers nor mass accretion.The fractions are ∼ 20 and ∼ 50 % for Pop III and = 10 −2 ⊙ stars, respectively.We may need to treat merger and accretion history in more detail in order to obtain more accurate PISN rate density.

SUMMARY
We have investigated features of merging binary BHs formed from isolated binary stars over all metallicities from Pop III stars to the solar-metal stars by means of rapid binary population synthesis simulation.We have found that the redshift evolution of the binary BH merger rate density in the fid set is consistent with the GWTC-3 result.The fid set can reproduce the primary BH mass and mass ratio distributions inferred by the GWTC-3 result.Notably, Pop III binary BHs dominantly contribute to merging binary BHs in the PI mass gap.The binary BH population does not contain binary BHs with negative effective spins.However, several studies have suggested that binary BHs with negative effective spins in the GWTC-3 result might be spurious.Thus, the absence of binary BHs with negative effective spins in our result may not matter.Overall, our result successfully reproduces the GWTC-3 result.
Our fid set can generate a large number of PISNe at the lowerredshift universe, compared to the L-3 in which binary BHs in the PI mass gap can also merge within the Hubble time.If Euclid space telescope will detect a few PISNe, we will conclude that the fid set is likely to be correct, and that Pop III binary stars dominantly produce binary BHs in the PI mass gap.Conversely, this means that we can use binary BHs in the PI mass gap as a probe for Pop III stars.In any case, ongoing and future observations will shed light on the origin(s) of merging binary BHs and Pop III studies.

Figure 5 .
Figure 5. Relation between ZAMS and remnant masses (solid) and between ZAMS and helium core mass just before core collapse (dotted) for = ⊙ , = 10 −2 ⊙ and Pop III stars.The top and bottom panels show the M and L models, respectively.Light-gray and dark-gray regions indicate those in which PPISN and PISN are effective, respectively, for = 10 −2 ⊙ and Pop III stars.= ⊙ stars do not experience PPISN nor PISN because of their small helium cores just before core collapse.

Figure 6 .
Figure 6.Redshift evolution of binary BH merger rate densities for the fid and L sets, indicated by solid and dashed curves, respectively.Black and blue curves show the merger rate densities of all binary BHs and Pop III binary BHs, respectively.The gray shaded region indicates the corresponding merger rate density inferred by GWTC-3 within 90 % credible bounds.

Figure 7 .
Figure 7.Primary BH mass distribution of merging binary BHs at a redshift of 0 for the fid (top) and L (bottom) sets.Black and blue curves indicate all the binary BHs and Pop III binary BHs, respectively.The gray shaded region indicates the corresponding primary BH mass distribution inferred by GWTC-3 within 90 % credible bounds.

Figure 8 .
Figure 8.Typical formation channel of a merging binary BH in the PI mass gap for the fid set.MS, PMS, and CCSN are short for main-sequence star, post main-sequence star, and core collapse supernova, respectively.The smaller spheres in PMSs indicate helium cores in PMSs.A dimensionless spin parameter ( ) is attached to each star.

Figure 9 .
Figure 9. Mass ratio distribution of merging binary BHs at the present day for the fid (top) and L (bottom) sets.Black and blue curves indicate all the binary BHs and Pop III binary BHs, respectively.The gray shaded region indicates the corresponding mass ratio distribution inferred by GWTC-3 within 90 % credible bounds.

Figure 10 .
Figure 10.Cumulative distribution of effective spins of merging binary BHs at the present day for the fid (top) and L (bottom) sets.Black and blue curves indicate all the binary BHs and Pop III binary BHs, respectively.The gray shaded region indicates the corresponding the cumulative distribution of effective spins inferred by GWTC-3 within 90 % credible bounds.

Figure 11 .Figure 12 .
Figure 11.Relation between ZAMS and remnant masses for = ⊙ , = 10 −2 ⊙ and Pop III stars in the L-3 set.Dark-gray and light-gray regions indicate that in which PISN is effective for Pop III and = 10 −2 ⊙ stars.=⊙ stars do not experience PISN because of their small helium cores just before core collapse.

Figure 13 .
Figure 13.Redshift evolution of intrinsic rate densities of PISNe for fid (solid) and L-3 (dashed) sets.Black and red curves indicate the rate densities of all PISNe and PISNe with naked helium stars as progenitors.