Gravitational wave quasinormal mode from Population III massive black hole binaries in various models of population synthesis

Focusing on the remnant black holes after merging binary black holes, we show that ringdown gravitational waves of Population III binary black holes mergers can be detected with the rate of $5.9-500~{\rm events~yr^{-1}}~({\rm SFR_p}/ (10^{-2.5}~M_\odot~{\rm yr^{-1}~Mpc^{-3}})) \cdot ({\rm [f_b/(1+f_b)]/0.33})$ for various parameters and functions. This rate is estimated for the events with SNR$>8$ for the second generation gravitational wave detectors such as KAGRA. Here, ${\rm SFR_p}$ and ${\rm f_b}$ are the peak value of the Population III star formation rate and the fraction of binaries, respectively. When we consider only the events with SNR$>35$, the event rate becomes $0.046-4.21~{\rm events~yr^{-1}}~({\rm SFR_p}/ (10^{-2.5}~M_\odot~{\rm yr^{-1}~Mpc^{-3}})) \cdot ({\rm [f_b/(1+f_b)]/0.33})$. This suggest that for remnant black hole's spin $q_f>0.95$ we have the event rate with SNR$>35$ less than $0.037~{\rm events~yr^{-1}}~({\rm SFR_p}/ (10^{-2.5}~M_\odot~{\rm yr^{-1}~Mpc^{-3}})) \cdot ({\rm [f_b/(1+f_b)]/0.33})$, while it is $3-30~{\rm events~yr^{-1}}~({\rm SFR_p}/ (10^{-2.5}~M_\odot~{\rm yr^{-1}~Mpc^{-3}})) \cdot ({\rm [f_b/(1+f_b)]/0.33})$ for the third generation detectors such as Einstein Telescope. If we detect many Population III binary black holes merger, it may be possible to constrain the Population III binary evolution paths not only by the mass distribution but also by the spin distribution.


Introduction
The final part of gravitational waves (GWs) from merging binary black holes (BBHs) is called the ringdown phase. When the remnant compact object is a black hole (BH), this phase is described by quasinormal modes (QNMs) of the BH (see, e.g., Ref. [1]). In general, the BH is expected as the Kerr spacetime [2], 2Mr dt 2 − 4Mar sin 2 θ dtdφ + dr 2 + dθ 2 + r 2 + a 2 + 2Ma 2 r sin 2 θ sin 2 θ dφ 2 , where = r 2 − 2Mr + a 2 and = r 2 + a 2 cos 2 θ , with mass M and spin a. In Eq. (1), we used the units of c = G = 1. The detection of QNM GWs not only gives a precise estimation of the BH's mass and spin, but also tests Einstein's general relativity (see the extensive review in [3]). In our previous paper [4], using the recent population synthesis results of Population III (Pop III) massive BBHs [5,6], we discussed the event rate of QNM GWs by second-generation gravitational

Population III binary population synthesis calculation
To estimate the detection rate of GWs from Pop III BBH mergers, it is necessary to know how many Pop III binaries become BBHs which merge within the Hubble time. Here, we use the binary population synthesis method of Monte Carlo simulation of binary evolutions. The Pop III binary population synthesis code [4][5][6] has been upgraded from the binary population synthesis code [12] 1 for Pop III binaries. In this paper, we calculate the same models as Ref. [6] using the same methods as Ref. [6] in order to obtain the mass ratio distribution and the spin distribution. In this section, we review the calculation method and models. Note that in this paper, we do not consider the kick models and the worst model discussed in Ref. [6] for simplicity, because in these models, BBHs have misaligned spins and the final spins after merger are too complex.
First, we need to give the initial conditions when a binary is born. The initial conditions such as primary mass M 1 , mass ratio M 2 /M 1 (where M 2 is the secondary mass), separation a, and orbital eccentricity e are decided by the Monte Carlo method with initial distribution functions such as the initial mass function (IMF), the initial mass ratio function (IMRF), the initial separation function (ISF), and the initial eccentricity function (IEF). For example, in our standard model, we use a flat IMF, a flat IMRF, a logflat ISF, and an IEF with a function ∝ e. There are no observations of Pop III binaries because they were born at the early universe. Thus, we do not know the initial distribution functions of Pop III binaries from observations. For the IMF, however, the recent simulations [13,14] may suggest a flat IMF, and therefore we adapt the flat IMF. For the other initial distribution functions, we adapt those of the Pop I case, where a Pop I star is a solar-like star. The above set of initial distribution functions is called our standard model of 140 cases with the optimistic core-merger criterion in this paper.
Second, we calculate the evolution of each star, and if the star satisfies a condition of binary interactions, we evaluate the effects of binary interactions and change M 1 , M 2 , a, and e. As the binary interactions, we treat the Roche lobe overflow (RLOF), the common envelope (CE) phase, the tidal effect, the supernova effect, and the gravitational radiation. The RLOF is stable mass transfer, while unstable mass transfer becomes the CE phase when the donor star is a giant. Here, we need some parameters for the calculation of the RLOF and CE phases.
In the case of the RLOF, we use the loss fraction β of transfered stellar matter defined aṡ PTEP 2016, 103E01 T. Kinugawa et al. whereṀ 2 is the mass accretion rate of the receiver star andṀ 1 is the mass loss rate of the donor star. In our standard model, β is determined by Hurley's function [12], which has been discussed for the Pop I case. When the receiver star is in the main sequence phase or in the He-burning phase, we assume that the accretion rate is described bẏ where τṀ is the accretion time scale defined by and the Kelvin-Helmholtz timescale τ KH,2 is defined by Here, M 2 , M c,2 , L 2 , and R 2 are the mass, the core mass, the luminosity, and the radius of the receiver star, respectively. When the receiver star is in the He-shell burning phase, we assume that the receiver star can get all the transfered matter from the donor star, i.e., Although we use the β function defined by Hurley et al. [12] in our standard model, we also treat the accretion rate of the receiver star described by the constant β parameter. This is because the accretion rate of a receiver star which is not a compact object is not understood well. Furthermore, in our previous study [6], we have shown that the Hurley fitting formula is consistent with β = 0 in the Pop III binary case. Thus, we also discuss the cases of β = 0.5 and β = 1. It is noted that the stability of the mass transfer changes if the mass transfer is non-conservative (β > 0). We use the criterion given in Ref. [15] as where M 1 and R L,1 are the mass and the Roche lobe radius of the donor star. If ζ ad = d log R ad,1 /d log M 1 < ζ L , where R ad,1 is the radius of the donor star, in the hydrostatic equilibrium of the donor star, the binary starts a dynamically unstable mass transfer such as the CE phase. When the receiver star is a compact object such as a neutron star or a BH, we always use β = 0 and the upper limit of the accretion rate is limited by the Eddington accretion rate defined bẏ   where κ T = 0.2 (1 + X ) cm 2 g −1 is the Thomson scattering opacity and X (= 0.76) is the H-mass fraction for Pop III stars. At the CE phase, the companion star plunges into the envelope of the donor star and spirals in. The orbital separation after the CE phase a f is calculated by the energy formalism [16] described by where a i , α, and λ are the orbital separation before the CE phase, the efficiency, and the binding energy parameter, respectively. In our standard model, we adopt αλ = 1. We also calculate the αλ = 0.01, 0.1. and 10 cases in this paper. Finally, if a binary becomes a BBH, we calculate the merger time from the gravitational radiation reaction, and check whether the BBH can merge within the Hubble time or not. We repeat these calculations and take the statistics of BBH mergers.
To study the dependence of Pop III BBH properties on the initial distribution functions and binary parameters, we calculate ten models with the Pop III binary population synthesis method [5,6] in this paper. Table 1 shows the initial distribution functions and the binary parameters of each model. The columns show the model name, IMF, IEF, the CE parameter αλ, and the loss fraction β of transfered stellar matter at the RLOF in each model.  show the initial mass ratio distributions and the mass ratio distributions of merging BBHs. The RLOF tends to make binaries of equal mass. Thus, the BBH mass ratio distributions depend on how many binaries evolve via the RLOF. Population III stars with mass < 50 M evolve as blue giants [5,17]. Thus, in the case of the IMF that light stars are in the majority, the binaries tend to evolve only via the RLOF, not via the CE phase. Therefore, the steeper IMFs tend to derive many equal-mass BBHs. In this calculation, since we adopt the minimum mass ratio as 10M /M 1 , the initial mass ratio distribution of models with the IMF that light stars are the majority is up to

The mass ratio distributions of binary black hole remnants
On the other hand, if we change the IEF, the mass ratio distribution does not change much. Thus, the dependence on the IEF is not so large (see Figs. 1 and 3).      For the CE parameter dependence (see Figs. 1 and 4), small mass ratio binaries in the model of αλ = 0.001 are much fewer than those in the other models. In the αλ = 0.001 model, all the binaries which evolve via the CE phase merge during the CE phase due to too-small αλ. Thus, the merging BBHs in this model evolve only via the RLOF, and become of equal mass by the RLOF. The change is not large between the models with CE parameters αλ = 0.1 and 10.
As for the mass loss fraction β (see Figs. 1 and 5), when β becomes large, there are three effects. First, binaries tend not to enter the CE phase. Second, the mass accretion by RLOF becomes not to be effective. Third, RLOF tends to finish early. The first effect makes binaries evolve via RLOF. However, the second and third effects have a negative impact on the tendency to become equal mass. Thus, the mass ratio distributions of the β = 0.5 and 1 models look similar to our standard model. 6

The spin distributions of binary black hole remnants
We calculate the spin evolution of binaries using the tidal friction. We use the initial spin distribution and the tidal friction calculations as in Refs. [5,12]. When the Pop III star becomes a BH, we calculate the BH's spin using the total angular momentum of the progenitor. If the estimated spin of the BH is more than the Thorne limit q Thorne = 0.998 [18], we assign the non-dimensional spin parameter q = q Thorne .
We ignore the spin up by the accretion during a mass transfer after the star became a BH for the following reason. The spin up by the accretion is calculated as where δJ , M BH , δM are the gain in angular momentum, the BH's mass, and the gain of the BH's mass, respectively. Since the accretion rate of the BH during RLOF is the Eddington rate, the gain of the BH's mass is where t life is the lifetime of the Pop III star. The Eddington accretion rate is given bẏ and the lifetime of the massive star is t life ∼ 1 Myr. As a result, we have δq ∼ 0.01, and the spin up by the accretion during RLOF is negligible. On the other hand, the accretion rate during the CE phase isṀ ∼ 10 −3 M yr −1 [19], and the timescale of the CE phase is about the thermal timescale of a red giant t KH ∼ 10 2 yr or less. As a result, we have δq ≤ 0.1, and the spin up by the accretion during the CE phase is negligible, too. Figures 6-15 show the spin distributions of merging BBHs and cross-section views of these spin distributions. The spins of merging Pop III BBH can be roughly classified into three types: group 1, in which both BHs have high spins q ∼ 0.998; group 2, in which both BHs have low spins; and group 3, in which one of the pair has high spin q ∼ 0.998 and the other has low spin.
When the BH progenitor evolves via the CE phase, the Pop III BH has low spin, and vice versa. If a Pop III star which is a giant evolves via the CE phase, the Pop III star loses the envelope and almost all the angular momentum due to the envelope evaporation. On the other hand, if the Pop III star evolves without the CE phase, the Pop III star can have a high angular momentum. Therefore, group 1 progenitors evolve without the CE phase and the envelopes of the progenitors remain. In group 2, both stars evolve via the CE phase and they lose their envelopes and almost all their angular momentum. In group 3, the primary evolves via the CE phase and the secondary evolves without the CE phase, or vice versa.
The  The distribution of q 2 for 0 < q 1 < 0.05. We can see that the q 2 distribution has bimodal peaks at 0 < q 2 < 0.15 and 0.95 < q 2 < 0.998. (b) The distribution of q 2 for 0.95 < q 1 < 0.998. We see that the large value of q 2 is the majority, so that there is a group in which both q 1 and q 2 are large.
which have high spins. In particular, in the case of the Salpeter IMF about 40% of the BBHs have spins q 1 > 0.95 and q 2 > 0.95. As for the IEF dependence, there is no tendency like the mass ratio distribution (see Figs. 6, 9, and 10).
The dependence on the CE parameter can be considered as follows (see Figs. 6, 11, 12, and 13). In the αλ = 0.01 model, almost all merging Pop III BBHs have high spins. About 60% of merging Pop III BBHs have q 1 > 0.95 and q 2 > 0.95 (i.e., group 1). The reason for this is that the progenitors which evolve via the CE phase always merge during the CE phase due to too-small αλ. Thus, the progenitors of merging Pop III BBHs in this model evolve only via RLOF and they do not lose angular momentum via the CE phase. In the case of the αλ = 0.1 model, the fraction of group 2 is lower than that of our standard model, like the αλ = 0.01 model. However, the fraction of group 1 is almost the same as that of our standard model, and the fraction of group 3 is larger than that of our standard model, unlike the αλ = 0.01 model. The reason for this is that although progenitors which enter CE phases twice merge during the CE phase due to small αλ, progenitors which enter 8      the CE phase only once do not merge during the CE phase, and the Pop III BBHs which cannot merge within the Hubble time in our standard model come to be able to merge within the Hubble time due to small αλ. In the αλ = 10 model, the shape of the spin distribution is almost the same as that of our standard model. The difference in this model from our standard model is the small increase of the fraction of group 2, because the progenitors which merge during the CE phase in our standard model come to be able to survive due to large αλ. As for the β dependence (see Figs. 6, 14, and 15), not only the stellar mass loss during RLOF but also the criterion of dynamically unstable mass transfer such as a CE phase are changed by β. In the β = 0.5 model, the fraction of group 1 is larger than that of our standard model because in this model the progenitors less frequently enter the CE phase than those of our standard model. In the β = 1 model, the fraction of group 1 is larger than that of our standard model, like a β = 0.5 model. However, the fraction of group 1 is smaller than that of the β = 0.5 model because in this model the progenitors lose a lot of angular momentum during the RLOF due to the high β. In this model, since the mass transfer cannot become dynamically unstable, the evolution passes via CE phases as follows. The progenitors of group 2 enter the CE phase when the primary and secondary become giants at the same time, and plunge into each other. On the other hand, the progenitors of group 3 enter the CE phase when the secondary plunges into the primary envelope due to the initial eccentricity.

Remnant mass and spin
Based on Ref. [11] (see also Refs. [20,21]), we calculate the remnant mass M f and spin q f from given BH binary parameters, M 1 , M 2 , q 1 , and q 2 (see Ref. [4] for a detailed discussion). The remnant mass and spin for each case is shown in Figs. 16-25. Here, we have normalized the distribution, and used binning with M f = 10 M for M f and q f = 0.1 (thick, red) and 0.02 (thin, blue) for q f .
The IMF dependence shown in Figs. 16, 17, and 18 is described below for the remnant mass and spin. When we treat the steeper IMF, we have a lower number of high-mass remnants. On the other hand, the number of high-spin remnants increases slightly in the steeper IMF cases. This is because in the steeper IMF models we have a large number of progenitors with mass smaller than 50 M .
As for the IEF dependence, we find that in Figs. 16,19, and 20 there is no strong tendency. Next, from Figs. 16,21,22, and 23 the CE parameter dependence can be described. In the αλ = 0.01 model, the maximum of the remnant mass becomes much smaller than that of our standard model. This is because the high-mass progenitors merge during a CE phase due to too-small αλ. As for the remnant spin, we do not have remnant spins which are smaller than 0.55 since BBHs tend to be of equal mass. If a light BH falls into a non-spinning massive BH, the remnant BH can have a small spin (q f < 0.6). However, in the above model many BBHs are equal mass. In the αλ = 0.1 model, the maximum remnant mass is smaller than that of our standard model again. In this model, the fraction of remnant spins with 0.7 < q f < 0.8 is larger than that for our standard model because the fraction of group 3 in this model is larger than in our standard model.    As for the β dependence, we find from Figs. 16, 24, and 25 that the maximum remnant mass becomes lower for higher β, due to the mass loss during RLOF.

Event rates for ringdown gravitational waves
To estimate the event rate for ringdown gravitational waves, it is necessary to have the merger rate density of Pop III BBHs. The merger rate density R m [Myr −1 Mpc −3 ] has been derived for various models in Ref. [6], and can be approximated by a fitting formula for low redshift. This is summarized in Table 2. In practice, we have considered the fitting for R m in terms of redshift z up to z = 2, but 15    the above R m is derived by using z ∝ D where D denotes the (luminosity) distance, because we use it only up to z ∼ 0.2 in this paper.
Using Ref. [22], we calculate the angle-averaged signal-to-noise ratio (SNR) as where we assume r = 3% of the total mass energy radiated in the ringdown phase. Note that for simplicity, any effect of the cosmological distance is ignored here. The symmetric mass ratio 16 PTEP 2016, 103E01 T. Kinugawa et al. are obtained from the remnant BH's mass and spin (see Ref. [23]). We evaluate the above SNR of the QNM GWs in the expected KAGRA noise curve S n (f ) [9,10] [bKAGRA, VRSE(D) configuration] (see Ref. [4] for the detailed calculation). This noise curve is presented in Ref. [26], and we use the fitting noise curve obtained in Ref. [24], based on Ref. [26]. Then, the event rate for a given SNR is derived by using the merger rate density in Table 2. In the right column of Table 2, we present the event rate with SNR > 8. Here, the event rates [yr −1 ] have been divided by dependence on the star formation rate SFR p and the fraction of the binary f b . In Fig. 26, based on Ref. [24], we show the parameter estimation in a case with SNR = 35 for the typical case [5,6] (with M rem = 57.0904 M and α rem = 0.686710). The (black) thick line shows the Schwarzschild limit and the ellipses are the contours of 1 σ , 2 σ , 3 σ , 4 σ , and 5 σ . In general relativity, the top-left side of the thick black line is prohibited. Thus, using the event with SNR > 35, we summarize various event rates in Tables 3 and 4. Table 3 shows the total event rates [yr −1 ] divided by dependence on the star formation rate SFR p and the fraction of the binary f b for ten models, and those for a remnant BH with q f > 0.7, 0.9, and 0.95.
In Table 4, we present the detection rate [yr −1 ] divided by dependence on the star formation rate SFR p and the fraction of the binary f b as a function of the lower limit of the solid angle of a sphere 4π C by which we can estimate the contribution of the ergoregion. The relation between this C and the spin parameter q was obtained in Ref. [27] (see also the recent studies in [25,[27][28][29]). It is noted that q f > 0.9 corresponds to C 0.97.

Discussions
In this paper, we extended our previous work [4] (the standard model in this paper) by looking at the dependence on various parameters of the Pop III binary population synthesis calculation. As shown in the right column of Table 2, the detection rate with SNR > 8 for the second-generation GW detectors such as KAGRA was obtained as 5.9-500 events yr −1 (SFR p /(10 −2.  by the mass distribution but also by the spin distribution. In particular, as described in Sect. 4, the spin of a black hole depends strongly on whether the progenitor of black hole enters the CE phase or not. Thus, we can check whether a BBH progenitor evolved via the CE phase or not by the spins of the BBH. One of the interesting outputs from the QNM GWs is whether we can confirm the ergoregion of the Kerr BH. From Table 4, the event rate for the confirmation of > 50% of the ergoregion is 0.040-3.1 events yr −1 (SFR p /(10 −2.5 M yr −1 Mpc −3 )) · ([f b /(1 + f b )]/0.33) with SNR > 35.
When we consider extracting the rotational energy of BHs using the Penrose process [31] or the Blanford-Znajek process [32], for example, we want to observe highly spinning remnant BHs. For remnant BHs with spin q f > 0.95, the event rate with SNR > 35 is 0.0027-0.037 events yr −1 (SFR p /(10 −2.5 M yr −1 Mpc −3 )) · ([f b /(1 + f b )]/0.33) in the KAGRA detector from Table 3. A third-generation GW observatory, the Einstein Telescope [33] will have an improvement in sensitivity of about a factor of ten over second-generation detectors. This means that we have roughly 1000 times higher expected event rates, and, for example, the ringdown event rate with 19 Here, we have introduced r as the fraction of the BH mass radiated in the ringdown phase, and assumed r = 3% to calculate the SNR and the event rates in this paper. If r = 0.3%, we will still have the possibility of detecting QNM GWs from highly spinning remnant BHs. Finally, Pop III BBH mergers can be a target for space-based GW detectors such as eLISA [34] and DECIGO [35]. Study in this direction is one of our future works.