Hidden universality in the merger rate distribution in the primordial black hole scenario

It has been proposed that primordial black holes (PBHs) form binaries in the radiation dominated era. Once formed, some fraction of them may merge within the age of the Universe by gravitational radiation reaction. We investigate the merger rate of the PBH binaries when the PBHs have a distribution of masses around ${\cal O}(10) M_\odot$, which is a generalization of the previous studies where the PBHs are assumed to have the same mass. After deriving a formula for the merger time probability distribution in the PBH mass plane, we evaluate it under two different approximations. We identify a quantity constructed from the mass-distribution of the merger rate density per unit cosmic time and comoving volume $\mathcal{R}(m_1,m_2)$, $\alpha = -{(m_1+m_2)}^2\partial^2 \ln\mathcal{R}/\partial m_1\partial m_2 $, which universally satisfies $0.97 \lesssim \alpha \lesssim 1.05$ for all binary masses independently of the PBH mass function. This result suggests that the measurement of this quantity is useful for testing the PBH scenario.


INTRODUCTION
Recent detections of gravitational wave events (GW150914, LVT151012, GW151226, GW170104, GW170608, and GW170814) by the LIGO-Virgo collaboration (Abbott et al. 2016c(Abbott et al. ,b, 2017a revealed the existence of binary black holes (BHs) in the mass range 8-35 M . These observations clearly demonstrate that there are numerous BH-BH binaries in the Universe that have previously eluded the scrutiny by astronomers. The origin of such heavy BHs and the formation of close binary BHs which merge within the age of the Universe are widely debated. Various astrophysical scenarios for the explanations of the gravitational wave events are summarized, for instance, in Abbott et al. (2016a) and Miller (2016).
Although only five robustly identified BH-BH binary mergers with GW detections have been reported so far, merger rates are constrained to within 12-240 Gpc −3 yr −1 (Abbott et al. 2017a). With the further improvement of GW detectors, we will soon enter the era of black hole rush where a large number of BH-BH binaries are detected with their masses, spins, and locations determined. Those data will serve us important clues to clarify the origin of binary BHs as well as the formation mechanism of the binaries. Clearly, investigations of how various astrophysical scenario producing merging BH binaries can be distinguished by observations will become a fundamentally important topic.
Recently, a collaboration including three of the authors, Sasaki et al. (2016) pointed out that the GW event GW150914 could be merger events of two primordial black holes (PBHs) based on earlier studies (Nakamura et al. 1997;Ioka et al. 1998). In Nakamura et al. (1997) and Ioka et al. (1998), the formation mechanism of the PBH binaries was proposed and a connection between the PBH binaries and the gravitational wave events from the merger of binary PBHs was given 1 . 1 There are other papers in which potential detection of PBHs by LIGO was claimed (Bird et  PBHs stand for BHs that formed in the very early Universe much before the epoch of the matter radiation equality (Carr & Hawking 1974). For instance, in the well-studied scenario, PBHs form from rare high peaks of the primordial density inhomogeneities whose amplitudes are much larger than the standard deviation. In this case, the PBH mass is given by the total energy contained in the Hubble horizon at the formation time, where T is the temperature of radiation and γ = O(1) depends on the details of the BH formation. Analytic estimates give γ = 3 −3/2 ≈ 0.2 (Carr 1975). Other mechanisms of the PBH production are summarized by Carr (2005). After having formed in the very early Universe, PBHs stay on the expansion flow of the Universe. Even when PBHs are randomly distributed in space without being clustered, there is a small but non-vanishing probability that two neighboring PBHs happen to be much closer than the mean distance. Such PBHs, being initially on the cosmic expansion flow, eventually start to come closer influenced by their mutual gravity when the cosmic expansion rate becomes too low to separate them apart. As was shown by Nakamura et al. (1997), a direct collision is avoided by the tidal effect of other PBHs in their vicinity, which leads to the formation of a PBH binary with a large eccentricity. Further Ali-HaÃŕmoud et al. (2017) have recently shown that the tidal field of halos and interactions with other PBHs, as well as dynamical friction by unbound dark matter particles, do not affect PBH binaries significantly. Highly eccentric PBH binaries radiate GWs efficiently and a fraction of them can merge within 14 billion years. In Sasaki et al. (2016), under the approximation that all PBHs have the same mass of 30 M , it was shown that the expected event rate of the PBH binary mergers is consistent 2016). The binary formation path is different from that in Sasaki et al. (2016). arXiv:1709.09007v2 [astro-ph.CO] 14 Jan 2018 with the one determined by the LIGO-Virgo collaboration after the announcement of GW150914 (Abbott et al. 2016d), if the fraction of cold dark matter in PBHs is about 10 −3 . This fraction is consistent with existing observational upper limits (Gaggero et al. 2017;Horowitz 2016;Brandt 2016;Koushiappas & Loeb 2017a;Inoue & Kusenko 2017;Green 2017;Matsumoto et al. 2017;Carr et al. 2017;Poulin et al. 2017). So far, the PBH scenario proposed by Sasaki et al. (2016) is successful in explaining the LIGO event GW150914.
In the next decades, many more BH binaries will be detected, which will deliver fruitful statistical information on the merger rates in the two-dimensional BH mass plane (m 1 , m 2 ) (see Abbott et al. 2016b;O'Leary et al. 2016;Mandel et al. 2017;Zevin et al. 2017;Kovetz et al. 2017;Fishbach & Holz 2017; (Abbott et al. 2016c(Abbott et al. ,b, 2017a. In this paper, we estimate the merger rate density in the m 1 − m 2 plane predicted by the PBH scenario. We extend the formalism of previous studies (Nakamura et al. 1997;Ioka et al. 1998;Sasaki et al. 2016) to compute the merger event rate to the case in which the PBH mass function is not restricted to a single-mass but it extends over a mass range between m min and m max with m max /m min 10 2 . We assume that the PBH mass function does not extend over many orders of magnitude since in that case the dynamics may not be accurately captured by the simple physical processes adopted by Nakamura et al. (1997); Ioka et al. (1998);Sasaki et al. (2016). Quite interestingly, we find that the merger rate distribution in this case depends on the mass of the BH binary in a specific way and a quantity constructed from the massdistribution of the merger rate density per unit time and volume R(m 1 , m 2 ), is insensitive to the PBH mass function. This distinct feature is advantageous since there is no theoretically tight constraint on the shape of the PBH mass function. Identifying the information in the merger rate density which is insensitive to the BH mass function may be used to discriminate different formation channels (O'Leary et al. 2016;Kovetz et al. 2017;Zevin et al. 2017;). This information may be used to obtain the probability of mergers for given BH masses, P intr (m 1 , m 2 ) (defined by Eq. (25) below) which is essential in measuring the underlying BH mass function f (m) itself. Before closing this section, in Table 1 we list definitions of important symbols that are used in this paper.
The paper is organized as follows. We first develop a formalism to compute the event rate in the PBH scenario which can be applied to the case of a non-monochromatic 3 mass function. Then, we apply the derived formula to evaluate the mass-dependence of the merger rate in the (m 1 , m 2 ) BH mass 2 Recently, such an extension has also been done in ). Our study differs from  in that our primary purpose is to investigate the universal feature of the merger-rate distribution that is insensitive to the PBH mass function.
3 By "monochromatic mass function" we refer to a population in which all PBHs have the same mass. Physical distance to i-th outer PBH (see Fig. 1) y i Comoving distance to i-th outer PBH θ i Angle (see Fig. 1) e i Vector (see Fig. 1 and Eq. (15)) x Comoving distance between PBHs that form a binary xmax Maximum value of x to form binary (see Eq. (9) plane and show that the special quantity constructed out of the event rate density becomes almost independent of the PBH mass function.

FORMATION OF BINARY PBHS
In this section, we derive a formula of the merger rate density as a function of the masses of two BHs comprising the binary.

Formation and mass function of PBHs
There are several mechanisms to form PBHs (Carr 2005). Among them, the most natural and widely investigated mechanism is the direct gravitational collapse of the primordial density perturbation in the radiation dominated Universe. In this scenario, when an overdense region containing an extremely high density peak in which the perturbation amplitude is greater than δ th = O(1) reenters the Hubble horizon, that region directly collapses to a BH (for the estimation of δ th , see (Carr 1975;Harada et al. 2013)). Crudely speaking, all the energy inside the Hubble horizon at the time of BH formation turns into the BH. This picture enables to relate the BH mass to the comoving wavenumber k of the primordial density perturbation as (3) There are no direct observational constraints on the probability distribution of density perturbations on such small scales. Although Eq. (3) gives us a simple and approximate estimate of the PBH mass in terms of k, the relation (3) is not precisely correct since the PBH mass also depends on the amplitude of the density perturbation. Deviation of the actual PBH mass from the horizon mass becomes significant as the amplitude of the density perturbation approaches δ th (Choptuik 1993;Niemeyer & Jedamzik 1998). Thus, even if the spectrum of the primordial density perturbation is monochromatic, the resulting PBH mass function is not monochromatic (Yokoyama 1998). Furthermore, the power spectrum of the primordial density perturbations needs not be monochromatic. In the paradigm of the standard inflationary cosmology, the primordial density perturbations are produced in the inflationary era preceding the radiation dominated era. Several inflationary models have been proposed to date with different predictions for the power spectral shape of the primordial density perturbation which lead to different PBH numbers and mass functions (see (Carr et al. 2016) and references therein). To a varying degree, these models predict a nonmonochromatic power spectrum. Thus, the PBH mass function is generally not concentrated on a single mass.
The PBH mass function is determined once the inflation model is fixed and the power spectrum of the primordial density perturbation is computed 4 . Since there is no fiducial inflation model producing PBHs and different models predict different PBH mass functions, we do not restrict our analysis to any particular PBH mass function. As mentioned earlier, our only requirement is that it is confined to the mass range m max /m min 10. The case where the PBH mass function is extended over many orders of magnitude requires a separate analysis, which is beyond the scope of this paper.
In addition to the mass function, the spatial distribution of PBHs also affects the probability of binary formation. In this study, for simplicity we assume that the distribution of PBHs at their birth is statistically uniform and random in space. However, we also have to keep in mind that primordial clustering of PBHs is also possible and could be an important factor to enhance the merger event rate for a fixed mass fraction of PBHs. We define the PBH mass function f (m) such that f (m)dm is the probability that a randomly chosen PBH has mass in (m, m + dm). Thus, f (m) is normalized as We denote the comoving PBH number density as n BH . The mean comoving separation between two neighboring BHs is thus given by n −1/3 BH . Before closing this subsection, it is important to mention that we do not consider the mass growth of the PBHs following their initial formation. The mass change due to accretion is negligible when PBH is in environments similar to the cosmic average density (Carr & Hawking 1974;Custodio & Horvath 1998;Ali-Haimoud & Kamionkowski 2017). This may not be true for PBHs residing in high density regions of galaxies such as molecular clouds, accretion disks, or stellar interiors. However, since the majority of PBHs are expected to remain mostly in low density regions such as dark matter halos, we ignore the mass growth of PBHs.

Major axis and eccentricity of a binary
Just after PBHs are formed in the early Universe, they are typically separated by super-Hubble distances. Apart from a possible peculiar velocity, each PBH is attached to the flow of the cosmic expansion. Let us denote the mass of a randomly selected PBH by m 1 , and the mass of and the comoving distance to the closest PBH by m 2 and x, respectively. Denoting the physical distance between the two BHs by D (see Fig. 1), the gravitational force is given by Gm 1 m 2 /D 2 . Ignoring for the moment the subdominant effects of the other remote BHs and the initial peculiar velocity and assuming that the above gravitational force is the only dynamical effect acting on each BH 5 , the BHs attract each other and collide within the free-fall time given by In reality, the space is expanding, and the BHs will be distanced if the space expands by O(1) or more within the freefall time. Conversely, if the free-fall time is shorter than the Hubble time 1/H, then the two BHs become gravitationally bound and eventually collide. Since the free-fall time and the Hubble time respectively scale as (scale factor) 3/2 and (scale factor) 2 during the radiation dominated era, the Hubble time may eventually exceed the free-fall time in the radiation dominated era even if the BHs are initially on the cosmic expansion flow (Nakamura et al. 1997). The condition for forming the bound system can be written as where z is the cosmological redshift. Using the Friedmann equation for a flat cosmology and neglecting factors of order unity, this condition can be rewritten as where ρ(z) is the background energy density. From this expression, we can give another but equivalent physical interpretation to the criterion for forming the gravitationally bound state. The left hand side is the total mass of the two BHs, and the right hand side is the total mass of whatever matter component that dominates the background Universe. Thus, the condition for two BHs to become gravitationally bound is equivalent to the condition for the total energy m t to exceed the background energy contained in the comoving volume to the nearest PBH x 3 . In the radiation dominated era, the energy density of radiation can be written as where z eq is the redshift at the time of matter-radiation equality, ρ c,0 and Ω m respectively represent a critical density and a density parameter of the non-relativistic matter at the present, x is smaller than x max given by Eq. (7) becomes satisfied at z = z dec > z eq , where z dec is given by The physical distance of the BH pair at the time of decoupling time, which becomes the semimajor axis of the resultant binary, is given by Since the BH pair forms only for x < x max , there is an upper bound on a as a < a max = x max /(1 + z eq ).
If there is no force other than the gravitational force from the neighboring BHs, and the initial peculiar velocities vanish, such two BHs come closer by moving on the same straight line and end up with a head-on collision. However, in reality, there are other remote BHs surrounding the BHs in pair, and they exert a torque during the infall motion of the BHs in pair. As a result, the BH pair acquires an angular momentum, and the head-on collision is circumvented. The torque exerted by the i-th distant BH to the lowest order in the distance D i to the i-th BH is given by where D is the physical distance between BH1 and BH2 (see Fig.1), M i is the mass of the i-th perturber BH, and θ i is the angle between a line connecting two BHs in pair and a line connecting i-th BH and a center of mass of the BH pair (see Fig. 1). Thus, the angular momentum generated by this torque throughout the free fall becomes Taking the direction of the torque exerted by each BH into account, the total angular momentum that the BH pair acquires is given by where we have chosen the line of the major-axis to be parallel to z-axis and is the unit vector pointing to the i-th BH (see Fig. 1). For the Keplarian motion, there is a relation between the orbital angular momentum and the eccentricity e as Using this formula, we obtain where x is the comoving distance between BH1 and BH2 and y i is the comoving distance to the i-th BH. Eqs. (11) and (17) are the main results of this subsection. They are the major axis and the eccentricity of the BH binary at the time of formation. Our analysis in the next subsection is based on these formulae. Let us now estimate the value of N, namely the number of the surrounding BHs that are inside the Hubble horizon at the time of the PBH binary formation. For simplicity, only in this paragraph we assume all the PBHs have the same mass m BH and constitute a fraction f PBH of all the cold dark matter (for instance, f PBH 10 −3 is required to explain the LIGO observation Sasaki et al. 2016). First of all, we notice that N depends on the initial comoving separation of the PBHs that form a pair. For instance, if the initial comoving separation of the BHs that form a binary is sufficiently small, they form a binary at very early time. In such a case, most likely few BHs exist inside the Hubble horizon and N = 0 or N = 1 will be the typical value. Thus, what we have to estimate is the typical value of N of PBH binaries that are relevant to observations. According to Sasaki et al. (2016), the probability dP that a given BH pair forms a binary, and then undergoes a merger at short cosmic time interval (t,t + dt) is given by where T is defined by For distinction between the lifetime and merger time of binaries, see discussion around Eq. (28). The merger probability for fixed t is dominated by the binaries having eccentricity near its upper limit e upper given by Eq. (11) in (Sasaki et al. 2016), We only consider the first case t < f 37/3 PBH T which is shown to be relevant to LIGO observations (Sasaki et al. 2016). For PBH mass m PBH = 30 M , this condition becomes f PBH 10 −3 . Analysis in the second case is straightforward. PBH binaries we are interested in are those that merge on the order of the age of the Universe t = t 0 ∼ 1/H 0 . Then, when we fix the merger time and the eccentricity to t 0 and e upper , respectively, the major-axis a at the time of the binary formation is uniquely determined (see Eq. (28)). Once the typical majoraxis is determined in this way, we can convert it to the typical redshift of the PBH binary formation by using Eqs. (10) and (11), from which we can evaluate the number of PBHs inside the Hubble horizon at that redshift, namely N. The result is given by Thus, for the typical PBH binary with m PBH = O(10 M ) which we are interested in, there are in general more than ∼ 3 × 10 10 PBHs in the Hubble horizon at the time of the binary formation if t t 0 . Because of the weak dependence of the PBH number N on the merger time t, N is much bigger than unity for merger times relevant to observations. In what follows, we take N → ∞.
One may wonder if the subsequent torque exerted on the BH binary by the surrounding BHs changes significantly the orbital parameters from the ones given by Eqs. (11) and (17). Considering the contribution only from the closest BH (i = 1) for simplicity, the angular momentum that the BH pair acquires during one period T of the orbital motion is given by While D does not increase with the scale factor because the BH pair is gravitationally bound, the distance D 1 grows in proportion to the scale factor which scales as ∝ t 1/2 in the radiation dominated epoch. Then, denoting by D (0) 1 the initial value of D 1 at the time of binary formation, D 1 when the BH pair is in the n-th cycle of the orbital motion becomes n 1/2 D (0) 1 . The accumulated angular momentum becomes Thus, the subsequent change of the angular momentum of the BH binary after its formation is at most a factor of ∼ 2. This factor is not important for our main result, and we do not consider this effect in the following analysis. On the other hand, note that if a distant third black hole with mass M 1 is captured on a bound orbit around the binary in a hierarchical configuration with some orbital period T 1 T and eccentricity e 1 , it can cause significant changes in the eccentricity of the binary due to the Lidov-Kozai effect on a timescale t Kozai = [(m t +M 1 )/M 1 ](1−e 2 1 ) 3/2 T 2 1 /T (Naoz 2016). However, we neglect this possibility in this paper for simplicity.
There are also other effects that have been ignored in deriving Eqs. (11) and (17). They include peculiar velocity of the individual BH seeded in at the time of BH formation, the radiation drag, the tidal interaction with the other PBHs in the matter dominated epoch, subsequent infall of the surrounding BHs to the BH binary, tidal force from the perturbations of non-PBH dark matter, and baryon accretion onto the PBH binaries. The first three effects are investigated in (Ioka et al. 1998) and was found to be subdominant. Recent study by Ali-HaÃŕmoud et al. (2017) also confirms that the tidal forces from outer PBHs do not significantly affect the late-time evolution of PBH binaries. The subsequent infall of the surrounding BHs is also studied in (Ioka et al. 1998). Ioka et al. (1998) assumed that the dark matter consists of a single-mass PBH population. In this case, the surrounding BH that caused the angular momentum of the BH binary at early times is eventually trapped by the BH binary if the outer BHs are within the mean distance of PBHs, which can be also understood from the expression of x max given by Eq. (9). Since the dynamics of three-body problem is difficult to solve, such a case was not considered, and only the opposite case where the nearest BH is more distant than the mean distance was included in the derivation of the merger event rate in (Ioka et al. 1998). Even under this restriction, it was found that the event rate is reduced by at most by 40%. On the other hand, in the present case where PBHs constitute only a fraction f PBH of all the cold dark matter, the mean distance is enhanced by a factor f −1/3 PBH compared with the case where PBHs provide all of the dark matter. Thus the probability that the surrounding BHs are trapped by the BH binary in the latter case is smaller than the former by a factor f PBH . Because of this consideration, we make an assumption that the surrounding BHs are not gravitationally bound to the BH pair. Then, the subsequent interaction by the surrounding BH in the BH binary is not significant, and we ignore the late-time effect of the surrounding BHs in the following analysis.
The tidal force from the surrounding density perturbations of cold dark matter not in the form of PBHs, exists when PBHs constitute only a fraction of entire dark matter. This issue was addressed by Eroshenko (2016) and Ali-HaÃŕmoud et al. (2017) who showed that the tidal effect is not significant by extrapolating the primordial perturbations on CMB scales down to the PBH scales (see also Hayasaki et al. 2016). Due to the random nature of the density perturbations, they yield additional statistically independent random contribution to ζ in Eq.(14). Since the power of the dark matter perturbation on small scales is not well understood, we do not consider this effect in this paper.
Finally, baryon accretion onto PBHs was claimed to significantly affect the PBH binaries and accelerate mergers in Hayasaki et al. 2016. But, recent study by Ali-HaÃŕmoud et al. (2017), based on the simple analytic calculation, suggests that the baryon mass accumulated on PBHs in Hayasaki et al. 2016 is likely to be an overestimation and the baryonic effect is much weaker although it may still be significant with respect to angular momentum exchange. For simplicity we do not account for baryon accretion in this work.

DISTRIBUTION OF THE MERGER RATE
In the previous section, we have derived the expressions for the major axis and the eccentricity of the PBH binary in terms of the initial comoving positions and masses of PBHs. They are the basic ingredients for the evaluation of the merger rate, which is the purpose of this section.
Let us denote by R(m 1 , m 2 ,t) a merger event density per unit cosmic time t and unit comoving volume in the m 1 − m 2 plane. In other words, represents the number of merger events of PBH binaries in the mass intervals (m 1 , m 1 + dm 1 ), (m 2 , m 2 + dm 2 ) that happen during (t,t + dt) and in the comoving volume dV . Since the merger time t can be inferred from the luminosity distance (depending on the cosmological parameters), and the source frame BH masses (m 1 , m 2 ) can be also estimated from the GW waveform, R is the quantity that can be in principle determined observationally. Our strategy to derive R(m 1 , m 2 ,t) is described as follows. What we have to evaluate is the probability P intr (m 1 , m 2 ,t)dt that a given BH pair consisting of two BHs with m 1 and m 2 , respectively, forms a binary, and then undergoes a merger during the short cosmic time interval (t,t + dt). Once the quantity P intr is obtained, using the PBH mass function given by Eq. (4) and assuming that the masses of the two PBHs in the binary are independent, the merger rate density R is given by The major-axis and the eccentricity of the BH binary at the formation time are given by Eqs. (11) and (17), respectively. From these equations, we see that the initial semimajor axis is a function of the random variable x as a ≡ a(x) and the initial eccentricity is a function of the length of the random vector ζ as e ≡ e(ζ) , where ζ = | ζ|. Denoting by F the probability distribution for x and ζ, the probability that the BH binary takes the values of the parameters in the range (x, x + dx) and (ζ, ζ + dζ) is given by We can then convert this probability into the one expressed in terms of a and e as This gives the probability that the BH binary at the formation time has the major-axis and the eccentricity in the range (a, a+ da), (e, e + de). PBH binaries shrink by emitting GWs until they finally merge. The lifetime τ of the BH binary with parameters (m 1 , m 2 , a, e) until it merges due to GW emission is given by 6 (Peters 1964) Denoting by t dec the cosmic time corresponding to z dec , namely the time of binary formation, we have τ = t − t dec .
Since PBH binaries that are relevant to GW observations merge at late time t t dec (t dec < 4 × 10 5 yr), it is a good approximation to identify τ with t. Thus, in what follows, we replace τ in all of the expressions with t. Under this approximation, we can express a as a function of {t, e, m 1 , m 2 } as a = a(t, e, m 1 , m 2 ). Using this relation, Eq. (27) becomes where it should be understood that a is replaced by {t, e, m 1 , m 2 }. Initial eccentricity of the BH binary is not a quantity that can be measured directly by the GW interferometers for primordial binaries and must be integrated. There is an upper bound e m for the initial eccentricity for fixed t because of the existence of the maximum value of the major axis a max = x max /(1 + z eq ) (see Sec.2.2). It is determined by the equation Notice that in the case of the monochromatic mass function e m coincides with e upper in the second case in Eq. (20). Finally, the intrinsic probability distribution is given by Having established the general framework to compute the merger rate density, let us implement this methodology in practice. It is straightforward to derive the last three factors in the integrand of Eq. (31), and they are given by (1 − e 2 ) −7/8 .
The highly non-trivial part is the evaluation of F(x(a), ζ(e)) since ζ depends on many random variables (in fact, infinite number of variables) in a complicated manner. Formally, it can be written as where Θ(·) is the Heaviside step function and δ(·) is the Dirac's delta function. Here, we have used the parametrization Eq. (15) for e i , and introduced the notation as y 0 = x, dV i = 4πy 2 i dy i and The derivation of Eq. (34) is given in the appendix A. We evaluate F(x(a), ζ(e)) using two approximations. The first case is that only the nearest BH (i = 1) is incorporated in the calculation of ζ. This approximation was adopted in the previous studies (Nakamura et al. 1997;Ioka et al. 1998;Sasaki et al. 2016) for single-mass PBH mass functions. In that case, all the PBHs have the same mass and the nearest BH (i = 1) exerts the strongest torque on the BH binary. Given that the torque by an outer BH is suppressed by the inverse cube of the distance, the approximation of taking only the nearest BH into account is physically natural as the zero-th order approximation 7 .
On the other hand, if the mass function is multimass, a massive outer BH may exert a stronger torque than a low-mass inner one. The wider the mass function, the more likely it is that this possibility may arise. To take into account the effect of outer perturbers, in our second estimate we consider a flat mass function up to a certain BH mass m max and include the outer BHs to evaluate the torque.
In what follows, we evaluate F(x(a), ζ(e)) and the intrinsic probability distribution for these two cases, separately.
3.1. Case 1: torque only from the nearest BH In this subsection, we make an approximation that the torque is exerted only by the nearest BH. Accordingly, the 7 The cumulative torque from all objects in a logarithmic radius bin of width ∆ ln y (e.g. here we may set ∆ ln y ∼ ∆y/y ∼ n −1/3 BH /y) follows from the central limit theorem and is described by a normal distribution with zero mean and root-mean-square that corresponds to ∆N 1/2 g 1,RMS , where ∆N is the number of objects in that logarithmic radius bin and g 1,RMS = 2 −1/2 (x/y) −3 (M RMS /mt) sin(2θ) RMS . This may be estimated roughly as ∆N = 4πn BH y 3 ∆ ln y. Therefore, the relative cumulative contribution of distant objects to the torque scales with y −3/2 , and so the smallest y dominates the integral where the number of objects is ∼ 1. function g defined by Eq. (35) becomes Even after this simplification, it is hard to evaluate the integral (31) analytically. For an analytic estimate, we carry out the calculation for an arbitrary but fixed value β = sin(2θ 1 ). Our result is insensitive to the value of β as long as it is not extremely close to zero. Since the probability of realizing β 1 is suppressed (see discussion after Eq. (37) for the estimation of this probability), we think that this simplification does not lose the essential feature of the merger-rate density. The integral over y 1 simplifies to The PBH binaries at the time of their formation are highly eccentric (e ≈ 1). Since the PBH mass function is implicitly assumed to be narrow in the present case, M 1 does not differ from m t significantly, and the argument of the last Heaviside function is positive unless β is smaller than 2 3 mt M1 √ 1 − e 2 . Now, let us estimate the probability that β becomes smaller than the critical value β c for which the argument of the Heaviside function becomes zero. To this end, we again consider the monochromatic mass function and use the eccentricity given by the first case of Eq.
For β c 1, the probability that β happens to be smaller than β c is approximately given by for the fiducial values used in Eq. (38). This probability is much smaller than unity, and we replace the last Heaviside function by 1 in the following analysis. Then, the intrinsic probability distribution (31) becomes where we have introduced a dimensionless parameter K by where f PBH is the mass fraction of the PBHs to the entire cold dark matter The integration over e can be expressed in terms of the incomplete gamma function. Then, Eq (40) becomes where m c and G(x) are defined by For the monochromatic mass function, m c is given by Eq. (43) for arbitrary f (M) mass function is the final expression of the intrinsic merger probability distribution in the present case.
3.2. Case 2: torque from the outer BHs Let us next consider the case in which the PBH mass function is flat from m min = m max to m max and vanishes outside of it. As mentioned earlier, we implicitly assume that 0.1. Then the PBH mass function is given by We include not only the nearest BH but also outer BHs. It is extremely difficult to perform the integration of Eq. (34) analytically 8 . However, we can estimate the approximate behavior of F(x, ζ) in the domain n BH x 3 1 where the PBH binaries with lifetime comparable to the age of the Universe form 9 . To this end, let us first write F(x, ζ) as where P(x, ζ 0 )dζ is a probability that ζ takes value in the interval (ζ 0 , ζ 0 + dζ) for given x. For later convenience, let us defineζ by (m t /m max )ζ. Thus, we have 8 Analytic expression of the probability distribution for the eccentricity was derived for the monochromatic mass function in Ali-HaÃŕmoud et al. (2017). 9 PBH binaries with n BH x 3 ∼ 1 have larger semimajor axis and more circular orbit than those with n BH x 3 1. These two factors make the lifetime of the binaries much longer than the age of the Universe.
whereP(x,ζ 0 )dζ is the probability thatζ takes a value in the interval (ζ 0 ,ζ 0 + dζ) for given x. Looking at the definition of ζ, we expect that the typical value ofζ for given x is around n BH x 3 since y i (i = O (1)) is typically about n −1/3 BH and the contribution from y i with higher i is suppressed (see footnote 7). Noting that y i > x, the case in which ζ n BH x 3 is realized by either if y 1 n −1/3 BH or if accidental cancellation takes place among terms with different i. Since the former is suppressed exponentially as ∼ e − 4π 3 nBHy 3 1 , the latter, which is stochastic, dominates. Recalling that ζ is essentially a two-dimensional vector, the probability thatζ is in the thin ring (ζ,ζ + dζ) by the random choice is proportional to the ring area, namelỹ ζdζ. Thus, we expectP forζ n BH x 3 . On the other hand, the caseζ n BH x 3 is realized mainly when y 1 is accidentally much smaller than the typical value n −1/3 BH . The probability of such a situation is controlled by the volume element y 2 1 dy 1 , and the relatioñ ζ ∝ y −3 1 leads to y 2 1 dy 1 ∝ζ −2 dζ. Thus, we expect forζ n BH x 3 . From the definition of ζ given by Eq. (17), we have The derivation of this result is given in appendix A. One simple function that interpolates Eqs. (50) and (51) is given bỹ whereσ = (m t /m max )σ, ξ = O(1) is a fitting parameter, and the normalization condition is imposed. In order to check the validity of the approximation (53), we evaluateP(x,ζ) numerically by the Monte Carlo method. For this purpose, we first fix N and x. Then, we randomly generate a set of random variables {M i , y i , θ i , φ i } and computeζ. By repeating this process many times, we obtain the distribution ofζ for a given N and x up to the statistical uncertainty. Figure 2 shows the distribution of ten thousand realizations ofζ for N = 5 for = 0.1, 4π 3 n BH x 3 = (10 −2 , 5 × 10 −3 , 2 × 10 −3 , 10 −3 ). The red curve represents the distribution obtained by the Monte Carlo calculations, and blue one represents the analytic approximation (53) with ξ = 5.5. We find that this simple ansatz ofP(x,ζ) fairly recovers the numerically obtained probability distribution. Although we consider the flat mass function, we expect that the ansatz should work qualitatively for other mass functions since the asymptotic behaviors (50) and (51)  Here we have defined a dimensionless quantity ν by and we have changed the integration variable as w = ν −32/37 (1 − e 2 ), and w m = ν −32/37 (1 − e 2 m ). Using a relation n BH = 2ρ BH /(m max (1 + )), which is valid for a flat mass function, we have .
To estimate typical magnitude of w m , for equal mass binary (m 1 = m 2 = m BH ), w m is given by .
This shows that w m can be bigger or smaller than unity within the range of the feasible values of f PBH and m BH . Although the integration over w in Eq. (55) can be expressed in terms of the hypergeometric function, we do not write it explicitly here since it gives no useful information. Thus, Eq. (55) is the final expression of the intrinsic merger rate and the main result of this subsection.

HIDDEN UNIVERSALITY IN THE MERGER RATE DENSITY
In the previous section, we have derived the analytic expression of P intr in the m 1 − m 2 plane for the two different limiting cases corresponding to the different approximations. According to Eq. (25), the observable merger rate density is not P intr , but P intr weighted by the PBH mass function. The observable merger event density is highly dependent on the PBH mass function, and it appears at first glance that no definite prediction can be extracted for the PBH scenario without choosing the specific mass function. Contrary to this naive guess, there is a unique feature expressed as a mathematical relation for the differentiated merger rate density specific to the PBH scenario as we will show below. Such a relation could be quite useful as a powerful method for testing the PBH scenario when the sufficient number of merger events have been accumulated.
Let us first consider the case where P intr is given by Eq. (40). This expression of P intr still contains the integration over the PBH mass nearest to the BH binary. Although this integration cannot be done explicitly without choosing the specific PBH mass function, carrying out the explicit integration is not needed for our present purpose. The function G(x) appearing in the integrand is monotonically decreasing and its asymptotic behavior is given as Using this formula and noting that K, which is much smaller than unity according to Eq. (42) (47) with = 0.1. Here n BH and x is the comoving PBH number density and initial comoving distance between BHs that form binary, respectively. Blue curves represent the probability distribution given by Eq. (53) with ξ = 5.5. we find that the integrand of Eq. (40) becomes A crucial consequence of these approximate expression is that the integrand has a simple scaling property with m 1 and m 2 . Using the scalings, we find that the above integrand scales as (63) Then, the observable merger rate density R per unit time and unit volume defined by Eq. (25) can be written as 3 37 f (m) and C A , C B are quantities that are independent of m 1 and m 2 , but contain information of f (m). An interesting point of Eq. (64) is that the dependence of the merger rate density on the total mass m t is independent of the model-dependent functions h A (m) or h B (m) (namely, mass function) and is completely determined as ∝ m t 36/37 for the former case and ∝ m t 22/21 for the latter case. The mass function enters the game only through the total normalization constant (represented as C A and C B ) and the factorizable part h A (m 1 )h A (m 2 ) or h B (m 1 )h B (m 2 ). Thus, by focusing on the total mass part of merger rate density and picking it up, we can provide a definite prediction for the merger rate density which is insensitive to the shape and amplitude of the PBH mass function. Indeed, we can pick up the total mass part by taking the logarithm of R and then differentiating it for any (m 1 , m 2 ). As discussed at the beginning of Sec. 3, the merger rate density R can be determined in principle by observations if a sufficient number of BH merger events are detected and the potential detection bias can be appropriately eliminated. Thus, the quantity α on the left-hand side can be also determined observationally. In this sense, the left hand side can be determined by observations. Our PBH merger scenario predicts that this quantity is equal to 36/37 for the upper case and 22/21 for the lower case. In reality, what is realized lies between the above two cases, and the left hand side of Eq. (65) may take a value between the two values corresponding to the upper case and the lower case respectively. Given that the numerical values on the right hand side for both cases are close to 1 (within less than 5%), the left hand side of Eq. (65) in the mixture case would be also close to 1. Taking into account this possibility, we conclude that under the assumption of the uniform spatial distribution of PBHs the merger rate density satisfies the following relation This relation is robust in the sense that it is independent of the underlying mass function. Similar conclusion can be drawn to the second case where P intr is given by Eq. (55). In this case, the observable merger rate density (Equation 25) is given by (67) As we have done in the case 1, let us evaluate the integral for two limiting cases (w m 1 and w m 1), separately.
First, when w m 1, we can extend the lower limit of the integral to 0. As a result, we obtain where C 1 is a constant of order unity. Using the scaling for ν as (see Eq. (56)) Advanced VIRGO, and KAGRA (Ghosh et al. 2016). 10 We generate a 2D histogram of events and fit the value of α. Repeating this analysis 1000 times for fixed fiducial α gives an approximate posterior distribution function of the measured α. This analysis shows that a sample of 100 events is necessary to measure α to integer accuracy and 1000 events would allow to measure it with an error of 0.15 if the fiducial value of α is between 1 and 3. The current rate estimates predict R = 12-240 Gpc −3 yr −1 . Assuming a maximum detection distance of z = 0.5 for the design sensitivity of second generation instruments, a sample of ∼ 100 events (1000 events) will accumulate in between 6 and 120 days (60 days and 3.3 years).
5. SUMMARY There is a growing interest in the possibility that the merging BHs detected by LIGO are primordial. Previous study (Sasaki et al. 2016) showed that the BH binary merger event rate estimated by LIGO can be explained by the PBHs which constitute only a tiny fraction of the entire dark matter. While the estimated masses of the individual BHs show some spread 10 ∼ 30 M , it was assumed in the previous study that all the PBHs have the same mass of 30 M . Although this is a reasonable approximation when only the first event for which masses of two BHs in the binary are almost the same is observationally known, it hugely compresses the valuable information about the event rate distribution in the BH mass plane.
In this paper, we extended the formalism to compute the merger event rate to the case where the PBH mass function is not monochromatic. Our basic assumption on the mass function made throughout this paper is that it is not widely extended over many orders of magnitude in the BH mass range but is confined to the mass range ∼ 10 M . The derived formula (31) contains multiple integrations over many random variables (Eq. (34)) and is complicated enough to defeat the exact analytic computation. Based on the physical expectation that among remote BHs, the closest one gives the largest torque on average, we evaluated the simplified version of Eq. (31) in which only the closest BH is taken into account. In this case, the computation becomes much more feasible. We found that the quantity α constructed from the merger rate density R in the BH mass plane as α(m 1 , m 2 ,t) ≡ −(m 1 + m 2 ) 2 ∂ 2 ∂m 1 ∂m 2 ln R(m 1 , m 2 ,t), (78) becomes almost independent of the PBH mass function and takes a value close to unity (0.97 α 1.05). Since it is possible that several distant BHs generate the dominant torque instead of the closest one during binary formation in the early universe, we have also considered the case in which the remote BHs are taken into account for a flat PBH mass function. Even in this case, we found that the quantity α exactly coincides with the one derived for the case of the closest perturbing BH. This suggests that the determined value of α is robust to observationally test the PBH scenario once a large sample of mergers becomes available with accurately determined masses. Other astrophysical mechanisms leading to BH mergers are generally expected to yield different α values. Recently, O'Leary et al. (2016) has shown that the probability of merger is proportional to m 4 t for binary BH mergers in dense star clusters, which implies α ∼ 4 if the merger rates are nearly independent to mass ratio. PBH binaries formed in the low redshift Universe by GW emission during close encounters leads to α ≈ 1.43 (Bird et al. 2016). BH binaries formed by GW emission in mass-segregated environments such as galactic nuclei lead to α values that vary with the total binary mass m t .
The mass distribution is not the only GW observable which allows one to distinguish between different mechanisms leading to binary BH mergers. For instance, it was shown recently that PBHs are unlikely to possess large spins (Chiba P(N,V ) ………… O 1 1 radius 1 2 1 Figure 3. This figure describes a situation where the individual BHs with mass m 1 , m 2 , M 1 , · · · and M N locate at the origin, (x, x + dx), (y 1 , y 1 + dy 1 ), · · · , and (y N , y N + dy N ), respectively.
By the assumption that M i obeys the uniform distribution in the interval ( m max , m max ), we have Thus, we obtain ζ 2 = 8 45 The calculation of i=1 1/y 6 i can be done by noting that it is an expectation value of 1/y 6 where y is the distance of particles randomly distributed in the region y > x (Ioka et al. 1998 (1 + + 2 ). (A9)