The merging black hole binaries detected by the LIGOVirgoKAGRA LVK gravitationalwave observatories may help us shed light on how such binaries form. In addition these detections can help us probe the hypothesized primordial black holes a candidate for the observed abundance of dark matter. In this work we study the blackhole mass distribution obtained from the LVK binary black hole merger events. We obtain that distribution by first associating a skewed normal distribution to each event detected with a signal to noise ratio SNR gt 8 and then summing all such distributions. We also simulate black hole binaries from two separate populations of merging binaries. One of these is a stellarorigin population that follows a massdistribution similar to the zeroage mass function of stars. The second population of black holes follows a Gaussian massdistribution. Such a distribution could approximate a population of black hole binaries formed from earlier black hole mergers in dense stellar environments or binaries of primordial black holes. For those populations we evaluate the number of detectable events and fit their combination to the LVK observations. In our work we rely on a wide range of stellarorigin blackhole mass distributions. We find that the observed LVK events can be fitted much better by the combination of such a stellarorigin mass distribution and a Gaussian distribution than by the stellarorigin mass distribution alone.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena', 'General Relativity and Quantum Cosmology', 'High Energy Physics - Phenomenology']
On the mass distribution of the LIGOVirgoKAGRA events
{'On the mass distribution of the LIGO-Virgo-KAGRA events': 'Mehdi El Bouhaddouti 1, ∗ and Ilias Cholis 1, †  \n1 Department of physics, Oakland University, Rochester Michigan, 48309,U.S.A  \nThe merging black hole binaries detected by the LIGO-Virgo-KAGRA (LVK) gravitational-wave observatories, may help us shed light on how such binaries form. In addition, these detections can help us probe the hypothesized primordial black holes, a candidate for the observed abundance of dark matter. In this work, we study the black-hole mass distribution obtained from the LVK binary black hole merger events. We obtain that distribution by first associating a skewed normal distribution to each event detected with a signal to noise ratio (SNR) > 8 and then summing all such distributions. We also simulate black hole binaries from two separate populations of merging binaries. One of these is a stellar-origin population that follows a mass-distribution similar to the zero-age mass function of stars. The second population of black holes follows a Gaussian massdistribution. Such a distribution could approximate a population of black hole binaries formed from earlier black hole mergers in dense stellar environments, or binaries of primordial black holes. For those populations, we evaluate the number of detectable events and fit their combination to the LVK observations. In our work, we rely on a wide range of stellar-origin black-hole mass distributions. We find that the observed LVK events can be fitted much better by the combination of such a stellarorigin mass distribution and a Gaussian distribution, than by the stellar-origin mass distribution alone.', 'I. INTRODUCTION': "The Laser Interferometer Gravitational-wave Observatory (LIGO), first detected gravitational waves (GW) from the coalescence of binary compact objects in 2015 [1]. Since then, following detector upgrades and the Virgo [2] and Kamioka Gravitational Wave Detector (KAGRA) [3] sites in Italy and Japan joining in the observations, there have been two more completed observing runs, O2 [4] and O3 [5]. As a result, the catalog of binary compact object mergers detected by the LIGO-Virgo-KAGRA collaboration (LVK) has reached 93 [6], with the O4 run currently ongoing. For each transient signal, the LVK collaboration evaluates the signal to noise ratio (SNR) and selects as merger events those with SNR>8. For every merging binary, LVK estimates the mass of the more massive compact object M 1 and of the less massive object M 2 , their ratio q = M 2 M 1 and the redshift z of the event, among other quantities. This allows for a study of the mass distribution of merging black holes (BHs). Among the 93 events in the LVK collaboration data, there are binary black holes (BBHs), binary neutron stars and neutron star-black hole binaries. In this paper we focus on the BBHs.  \nThe two LIGO detectors at Hanford Washington and Livingston Louisiana, having the better sensitivity, offer the main means by which GWs from merging binaries are detected. Virgo and KAGRA allow for an improved localization of the events. The LIGO detectors are versions of a Michelson interferometer with 4 km long arms [7], capable of detecting GWs in frequencies of 10 Hz to 1000  \nHz. Thus they are ideal for the detection of coalescing stellar-mass BBHs.  \nBBHs can have two origins, astrophysical, from collapsing stars and primordial, forming in the early universe. Astrophysical BHs can be expected to have an approximately power-law mass distribution, similar to that of their progenitor stars. For stars with a mass larger than 1 M ⊙ (including massive enough to form BHs after their collapse), their initial mass follows a power-law distribution of dN/dM ∼ M -2 . 3 ± 0 . 7 [8]. In this paper, we use a power-law model to simulate astrophysical BHs and test the resulting distribution from our simulations with the data from the LVK collaboration, by searching for the best fit power-law parameters.  \nPrimordial black holes (PBHs), differ from astrophysical ones in their origin. PBHs would form from curvature fluctuations in the primordial universe that collapse to black holes. Those fluctuations/overdensities are model dependent. They can be related to a phase transition during inflation (see for example [9]). Due to the 'no-hair' conjecture, BHs are defined by three quantities: mass, spin and charge. PBHs would be indistinguishable from astrophysical BHs. However, the dynamics of PBHs and astrophysical BHs are different, with astrophysical BH binaries closely following the baryons in galaxies while PBHs being in dark matter halos. Given their different origins and environments, astrophysical and primordial BHs constitute two separate populations. Testing PBHs is motivated given that if present, then they would compose at least a fraction of dark matter (DM) in the mass range where the current GW detectors are sensitive. Moreover, measuring their mass would shed light to the primordial curvature fluctuations in the early universe and indirectly probe an energy scale impossible to reach with high-energy particle colliders. In this work, we consider for simplicity that PBHs have a monochromatic mass-distribution. Given that the LVK detectors  \nhave a finite mass resolution, such a monochromatic mass distribution would appear in the data as a Gaussian-like component.  \nIn this paper, we present a methodology in studying the LVK BH mass-distribution data. In Section II, we discuss how we take the reported M 1 and M 2 values to construct M 1 - and q -distributions for the BBHs observed with a false alarm rate (FAR) smaller than 1 yr -1 . This is approximately similar to a SNR > 8 . Subsequently, in Section III, we simulate BBHs from the two populations and compare to the data. Our results are discussed in Section IV. We find that the LVK BBHs require the presence of a Gaussian-like component in addition to a power-law mass-distribution. Finally, in Section V, we present our conclusions and discuss how such an analysis can be used to probe for the origin of that Gaussian-like component.", 'A. Studied population': 'We use the observations from the LVK O1, O2 and O3 runs [5] obtained from the gravitational wave open science center [6]. Out of the 93 events in the catalog we keep the events that are mergers of BBHs. We use the primary mass M 1 and secondary mass M 2 > 4 M ⊙ information, for events with a FAR smaller than 1 yr -1 . These conditions restrict the number of events that we work on, to 72.', 'B. Probability density function of primary mass and binary mass-ratio': "For our study, we need to obtain, for each of the 72 events, the probability density functions (PDFs) for the primary mass M 1 and the mass-ratio q [5]. The reported error-bars from LVK on each of these quantities are asymmetric. Thus, we choose a skew-normal distribution to fit the M 1 and M 2 PDFs for each event. We use the following definition of the skew-normal distribution [10],  \nf ( x, ξ, ω, α ) = 1 √ 2 πω exp ( -( x -ξ ) 2 2 ω 2 )(1 + erf ( α ( x -ξ ) ω √ 2 ) . (1)  \nThe parameter ξ is used to describe the location along the x -axis, ω is the scale parameter and α is the shape parameter (for α = 0 we get the normal distribution). We want to fit the PDF of an event's M 1 so that the mode, Y ( ξ, ω, α ) of the distribution is located at the central value of M 1 from the LVK data. We also want the scale of the distribution to be consistent with the claimed uncertainty on M 1 . There is no analytic expression for the mode of the distribution for α = 0 . However, we can  \n̸  \nobtain the approximate expression using [10],  \nδ ( α ) = α √ 1 + α 2 (2)  \nand  \nγ 1 = 4 -π 2 ( δ ( α ) √ 2 /π ) 3 (1 -δ ( α ) 2 /π ) 3 / 2 . (3)  \nγ 1 is the skewness parameter of the distribution. Also, setting parameters,  \nµ z ( α ) = √ 2 π δ ( α ) , (4)  \nσ z ( α ) = √ (1 -µ 2 z ( α )) (5)  \nand  \nm 0 ( α ) = µ z ( α ) -γ 1 σ z ( α ) 2 -sign ( α ) 2 exp ( -2 π α ) , (6)  \nwe get the mode of the distribution expressed as,  \nY ( ξ, ω, α ) = ξ + ωm 0 ( α ) . (7)  \nFor each event, we search for the optimal combination of parameters ξ , ω and α that fits the relevant PDF of M 1 1. We also obtain the optimal combination of ξ , ω and α parameters for M 2 in the same way. In Fig. 1, we show as an example, the estimated PDF for event GW 190413 \\_ 052954 with a primary mass M 1 = 33 . 7 M ⊙ .  \nFIG. 1. Using as an example GW 190413 \\_ 052954 , we plot the PDF of primary mass M 1 = 33 . 7 M ⊙ (shown as a grey line). The primary mass, within 90% credible interval, has an upper limit M upper 1 = 33 . 7 + 10 . 4 = 44 . 1 M ⊙ and a lower mass limit M lower 1 = 33 . 7 -6 . 4 = 27 . 3 M ⊙ (light orange dashed lines).  \n<!-- image -->  \nTo obtain the PDF for q of each event, we pick a random value from the event's distribution of M 2 and divide by a random value from the event's distribution of M 1 . Since q is by definition the ratio between the less massive event and the more massive event, in the case  \nwhere the value randomly picked from the M 1 distribution ( M 1 , expl ), is smaller than the randomly picked value from the M 2 distribution( M 2 , expl ), i.e M 1 , expl < M 2 , expl , the division yields a value of q that is larger than 1. We pick in this case 1 /q instead. We repeat this process a thousand times to obtain the PDF for q for that event; and proceed in the same way to obtain the PDF for q for all events.", 'C. Mass-distribution properties of the observed binary black holes': 'With the skew-normal parameters for the quantities M 1 and M 2 , we have access to the PDFs of M 1 and q for all events. For each event, we make a histogram from the PDF of M 1 , using 15 discrete bins. We give our M 1 mass-binning in Table I. We record the number of counts per bin ( n i,j with i indexing the events and j the bins). To obtain the mass distribution of the M 1 parameter for the observed LVK BBHs, we make a histogram with bin counts n tot,j (with j the index of the bin) equal to the sum of the bin counts from the histograms of each event,  \nn tot ,j = ∑ i n i,j . (8)  \nThe bin counts of our final histogram for the M 1 distribution is then n tot ,j for bin j . In Fig. 2, we provide the derived M 1 -histogram for the 72 events with FAR of < 1 yr -1 .  \nTABLE I. The parametric binning of the histograms for M 1 and q .  \nWe proceed in the same way starting from the PDFs of q for each event to obtain the histogram for the q distribution of the 72 BBH mergers, the q -binning is also given in Table I and shown in Fig. 3.  \nSince we have a small number of events per bin from the O1, O2 and O3 runs, we use Poisson distribution to obtain the error bars for these histograms. The 1 σ -error for each bin in the mass distribution of BBHs is the 1 σ -error of the Poisson distribution drawn from the number of bin counts in the bin considered. This method is valid for the histograms of both M 1 and q .', 'III. SIMULATED BBHS AND THE χ 2 TEST': 'To reach a better understanding of the mass distribution of the BHs detected by the LVK collaboration, we simulate BHs following a power-law and a powerlaw + Gaussian peak (PP) mass-distribution. We place  \nthe BHs at a redshift following a power-law distribution for their comoving rate density,  \ndR d ( z +1) ∝ (1 + z ) κ . (9)  \nRef. [5], used the same parameterization for the redshift distribution of BBHs and found a best fit range for κ = 2 . 9 +1 . 7 -1 . 8 . We use the central value of κ = 2 . 9 . Since the observed event with the largest redshift has a redshift with a central value of 1.18, we simulate events with a redshift of up to 1.5. Then, on the simulated events, we calculate the SNR observed by the LIGO interferometers (discussed in further detail in section III B). We select the simulated events with a SNR larger than 8 and normalise the distribution of detectable events to perform our χ 2 test.', "A. The BBH properties' distributions": "We consider two mass distributions from which we will simulate BBHs. We start with the power-law distribution, in which the primary mass M 1 and q are drawn from,  \ndN dM 1 ∝ M -α 1 . (10)  \ndN dq ∝ q β . (11)  \nWhere α is the exponent of the power-law distribution from which we draw M 1 and β the exponent of the powerlaw distribution from which we draw q . We consider the power-law distribution as a simplified mass distribution of astrophysical BHs.  \nThe second distribution we consider is a PP distribution, which is a power-law distribution with an additional Gaussian component. The parameters are the exponent of the distributions for M 1 and q , i.e. α and β and the location parameter of the Gaussian peak µ . We are using the Gaussian in this distribution to simulate monochromatic PBHs or BHs from earlier mergers. Since PBHs would be indistinguishable from astrophysical black holes, we pick for the width of the Gaussian σ , to be the average error of the primary and secondary BH masses detected by LVK in an interval centered at the considered peak location. This approximation consists in averaging the errors on the detected BH masses at different redshifts. This is a simple way to estimate the LIGO detectors' sensitivity accounting for how that sensitivity changes with respect to redshift. Given the small sample of BBH detections, we consider this an adequate approximation, leaving a more sophisticated approach for future comparisons to observations with many more detections. We note that in order to have a better  \nFIG. 2. Histogram of M 1 mass distribution for BBHs, with Poisson error bars from GWTC3 events.  \n<!-- image -->  \nestimate of the mass sensitivity of the LIGO detectors in this mass range, both primary and secondary masses in the relevant mass interval are used. Our approximation allows to reduce the number of parameters probed for PP distribution, as the width of the Gaussian becomes dependant on the location of the peak.  \nUsing the PP model, we simulate events from the power-law component and events from the Gaussian component independently. If that component is associated to PBHs, this assumes that there is no term in the observed events from merging binaries consisting of a PBH and a regular BH from stellar evolution. Such binaries are expected to be very rare. PBHs occupy the entire dark matter halo. The PBH-PBH merger events occur at the smaller mass dark matter halos and at volumes within the halos that is not occupied by much regular baryonic mass and thus stars [11-14]. Only a small fraction of PBHs go through the inner part of a galaxy where BH from stellar evolution may lay. In fact, the PBHs have relatively high velocities as they cross the inner galaxies, suppressing their direct capture cross section by other BHs.", 'B. Evaluating the Signal to Noise Ratio': 'After simulating BBHs and placing them at a certain redshift, we need to calculate the SNR of the simulated events observed by LVK. Since the two LIGO interferometers are significantly more sensitive than the Virgo  \nor KAGRA ones, we evaluate the SNR from the just the LIGO observatory. We use the LIGO noise curves of Ref. [15], re-plotted in Fig. 4 for reference. Ref. [15], provides one noise curve for the Hanford site and one noise curve for the site in Livingston based on the first three months of the O3 run. These noise curves do not account for the change in noise with respect to time. We have noticed that this leads to typically less massive events being detectable in our simulations than claimed in the small sample of LVK collaboration. We leave is point to be further discussed in our Results section IV.  \nTo calculate the SNR we need the expressions of the merger and quasinormal (qnr) mode frequencies as well as the expression of the energy density spectra of the inspiral, merger and ringdown phases. All these expressions are in geometrised units ( G = c = 1 ), [16]. The merger frequency is given by,  \nf m = 0 . 02 M . (12)  \nThe qnr frequency is,  \nf qnr = 0 . 13 M . (13)  \nThe energy density spectrum of the inspiral phase for a frequency lower than f m is given by,  \ndE df i ( f ) = 1 3 π 2 / 3 µM 2 / 3 f -1 / 3 , (14)  \nFIG. 3. Histogram of the q -distribution of BBHs, with Poisson error bars from GWTC3 events.  \n<!-- image -->  \nevaluated at the source. The energy density spectrum of the merger phase for a frequency between f m and f qnr is given by,  \ndE df m ( f ) = 0 . 91 M 2 F ( µ/M ) . (15)  \nThe energy density spectrum of the ringdown phase for a frequency around f qnr is given by,  \ndE df r ( f ) = 1 8 A 2 · Q · M 2 f qnr . (16)  \nWith M the total mass: M = M 1 + M 2 , µ the reduced mass : µ = M 1 M 2 M 1 + M 2 and F ( µ/M ) = ( 4 µ M ) 2 . We pick A=0.4 and Q = 2(1 -a ) -9 20 with a=0.67 for the ringdown energy density spectrum.  \nWith the expressions on the energy spectra of the phases of BBHs coalescence, we can express the characteristic amplitude of the signal,  \nh 2 c = 2 (1 + z ) 2 ( πD L ( z )) 2 dE df ((1 + z ) f ) . (17)  \nWith z the redshift at which the BBH is placed, dE df is the energy spectrum of the coalescence, which is the sum of the energy spectra of the inspiral, merger and ringdown phases: dE df = dE df i + dE df m + dE df r . And D L ( z ) is the luminosity distance of the BBH at a regshift z,  \nD L ( z ) = (1 + z ) ∫ c H 0 √ Ω m (1 + z ) 3 +Ω Λ dz, (18)  \nwhere c is the speed of light, H 0 the Hubble constant, Ω m the matter density parameter and Ω Λ the dark energy density parameter obtained from "the best fit Euclidian Λ CDM cosmology from Plank Collaboration(2018b)" according to [17] [18]. Since h c is measured on Earth we need to account for the redshift effect on the frequency, we thus take the energy spectum at (1 + z ) f . From [15], we obtain a strain noise for the Livingston detector that we note h n,L and a strain noise for the Hanford detector that we note h n,H .  \nTo evaluate the SNR ρ we combine them as,  \n< ρ > = √ 4 5 ∫ h c / (2 f ) 2 h 2 n,L + 4 5 ∫ h c / (2 f ) 2 h 2 n,H . (19)  \nWe use the two noise curves separately in the expression of the SNR because the Livingston noise we use is consistently lower than the Hanford one (see Ref. 4). Using the two noise curves independently in the expression of the SNR instead of an average of two curves does not change the value of the SNR significantly. Nonetheless, the expression of Eq. 19 is more accurate.  \nWe simulate 10 6 BBH with masses from 4 to 150 M ⊙ and up to a redshift of 1.5. Since most generated binaries will have large luminosity distances and thus a SNR < 8 , we make a first estimate to disregard binaries for which our SNR will be well bellow 8. This allows us to optimise significantly our computing time of SNR for the remaining sample of merging binaries. In simple terms, we need to evaluate the maximum redshift at which a binary described by M 1 and q would give a SNR ≥ 8 . We  \nFIG. 4. The noise curves of LIGO Hanford and LIGO Livingston used in this work, taken from [15].  \n<!-- image -->  \nnote that we do account for the inclination angle as this is not included in Eq. 14. Discarding detactable BBHs in the study would result in potentially deforming the distributions from which we draw the quantities that make our simulated BBHs. This is why, for each range of primary masses M 1 , we pick a cutoff value for the redshift z that is slightly larger than the found maximum redshifts. This is to ensure that all detectable simulated BBHs are included in the study. In Table II, we show the largest redhsift for specific choices of M 1 . For low M 1 values, this translates to setting q = 1 . However, for very high M 1 values, the low-frequency noise curve prevents the detection of binaries with q ≃ 1 .  \nTABLE II. The largest redshift at which a binary with primary mass M 1 , would be observed with SNR > 8 . We test the range of possible values for q . We simulate 10 6 BBHs and present our result in the second column. The third column gives cutoff value of z that we pick for that M 1 value, in selecting binaries (see text for details).', 'C. Testing models with data': 'We calculate the SNR of the simulated BBHs that LIGO would detect and save the primary mass, the ratio q and the redshift of the simulated events with an SNR > 8 . These events constitute our simulated model. We make a histogram of M 1 and q distributions of the simulated model keeping the same number of bins that we use for the LVK data with the same range for M 1 and q .  \nOne of our objectives, is to find the minimal χ 2 for each choice of parameters. This minimisation is done  \nby searching for the normalisation factors that minimise the χ 2 value. For the power-law model, one normalisation factor is enough to minimise the χ 2 value. The candidates for the minimal χ 2 used to compare the M 1 distribution of the simulated events for one value of α and one value of β with the M 1 distribution of the LVK data have the following expression,  \nχ 2 M 1 = 15 ∑ i =1 ( data M 1 ,i -N · model M 1 ,i sigma M 1 ,i ) 2 . (20)  \nWhere, i labels the bins and goes from 1 to 15 for M 1 , data M 1 ,i is the number of M 1 counts from the LVK data in the bin i , model M 1 ,i is the number of M 1 counts from the simulated events in the bin i , N a normalisation factor, and sigma M 1 ,i is the 1 σ error on data M 1 ,i .  \nSimilarly, the candidates for the minimal χ 2 for the q distribution have the following expression,  \nχ 2 q = 10 ∑ i =1 ( data q,i -N · model q,i sigma q,i ) 2 . (21)  \nAgain, i labels the bins and goes from 1 to 10 for q , data q,i is the number of q counts from the LVK data in the bin i , model q,i is the number of q counts from the simulated events in the bin i , N a normalisation factor, and sigma q,i is the 1 σ error on data q,i .  \nThe total χ 2 is then,  \nχ 2 tot = χ 2 M 1 + χ 2 q . (22)  \nFor the PP model, we have two components for the distribution, a component following a power-law distribution and a component following a Gaussian distribution. The expression for χ 2 M 1 and χ 2 q that is very similar to equations 20 and 21,  \nχ 2 M 1 ,PP = 15 ∑ i =1 ( data M 1 ,i -model M 1 ,PP ,i sigma M 1 ,i ) 2 , (23)  \nwhere,  \nmodel M 1 ,PP ,i = N PL · model M 1 ,PL ,i + N Peak · model M 1 ,Peak ,i . (24)  \nSimilarly,  \nχ 2 q,PP = 10 ∑ i =1 ( data q,i -model q PP ,i sigma q,i ) 2 (25)  \nwith,  \nmodel q PP ,i = N PL · model q PL ,i + N Peak · model q Peak ,i . (26)  \nIn this case we have one normalization factor for the power-law component, indexed with "PL", and one for  \nthe Gaussian component, indexed with "Peak". Bin counts from the simulated events also come from two different components of the distribution, the same indices introduced on the normalization factors apply here. The total χ 2 result of the PP model is again given by Eq. 22 and minimised to find the best-fit combination of normalization factors N PL and N Peak .  \nNote that having two different normalization factors for the power-law component and the Gaussian components of the PP model is consistent with what the two components represent; since the BHs from stellar evolution and the PBHs have different origins and different dynamics 1 .', 'A. Power-law model': "For the single power-law model, from the χ 2 test we find the best fit values of α = 2 . 26 +0 . 21 -0 . 30 and β = -0 . 26 +0 . 62 -0 . 24 . These give a total χ 2 of χ 2 tot = 49 . 13 . The quoted confidence intervals for the best fit parameters are 90% . The χ 2 minimization was done with one parameter (the normalization factor for the power-law model), leaving us with 24 degrees of freedom for this model. The χ 2 per degree of freedom ( χ 2 d.o.f. ) in this case is then χ 2 d.o.f. = 2 . 05 . Table III, includes this result along with χ 2 values for specific parameter choices as fixing β to 0 and 1 or α to 3.5 and reporting the best-fit value for the other parameter. Compared to Ref. [5], we find significantly lower values for α , i.e. a less steep power-law distribution and lower values for β , i.e. preference for smaller values of q . This should not be surprising as the regular zero-age mass function for stars, gives α ≃ 2 . 3 for the massive end of their distribution. Given the quoted uncertainties of Ref. [5], our best fit value for α is beyond their 90 % confidence interval. That is not the case for the β parameter however. In our opinion this suggests that this method of analyzing the gravitational wave merger data is somewhat orthogonal to that of Ref. [5]. With a better understanding of the LIGO noise-curve evolution, the differences between the two methods should become less pronounced.  \nIn Fig. 5, we show the area in parameter space around the best-fit values that is within 90% confidence interval. The best-fit point is given by the red 'x' symbol. In Figs. 6 and 7, we show the M 1 and q histograms of the simulated BBHs from a power-law model with the bestfit parameters. These are shown as blue histograms with red best-fit normalization values (no errors shown). For comparison we include the histograms constructed from the LVK data (grey histogram with black Poisson errors).  \nTABLE III. Results of the χ 2 test after comparing simulated populations from the power-law model to the LVK data. The number of degrees of freedom for this test is 24. The first line shows the best fit parameters and χ 2 results, the second line the best fit results for β = 0 , the third line results for β = 1 , and last line results for the parameters reported by the LVK collaboration that yield the best fit for the merger rate of Ref. [5].  \nAs can be seen from Fig. 6, the single power-law can not explain the peak around 40 M ⊙ on the M 1 -distribution. This is also clear from the best-fit χ 2 d.o.f result for the power-law model. It is signigicantly larger than one. The power-law model, is not sufficient to describe the mass distribution of the LVK collaboration's events. It is important to note that when working with the power-law model, we are effectively considering mergers between first-generation BHs. Accounting for hierarchical mergers may require a more complex modeldistribution and lead to different results. In fact, our result suggest the need for some additional component to the simple power-law model, which at a first level we address in Section IV B.  \nThe results obtained in this section relied on calculating the SNR from simulated mergers that would be detected with by LIGO. The obtained SNR is dependent on the noise curve that we use for the two LIGO interferometers. We have access to only one noise curve for the LIGO Hanford and one for the LIGO Livingston sites [15]. LIGO's noise is time dependent and keeps changing as the detectors are undergoing improvements. While testing our code to calculate the SNRs, we compared the SNR we obtained for specific LVK's events with their reported SNR. We noticed that some events with large masses where detected by LVK with a SNR larger than 8, while we obtained a SNR smaller than 8. This can be explained by a seismic noise at the time of observation that was smaller than the seismic noise accounted for in the noise curves we have available. Having access to noise curves for the Hanford and Livingston for different times during the observation period would solve this problem.", 'B. Power-law + Peak model': "As discussed in III A, the central value µ of the Gaussian distribution has a corresponding width σ , that its value is correlated to the value of µ . To obtain σ , we average the 90% errors in the LVK catalog from all M 1 s and M 2 s, each time within an interval spanning from µ -5 M ⊙ to µ +5 M ⊙ . The derived result is then divided by 2 and  \nFIG. 5. The α -β parameter space for the power-law model probed by the LVK BBHs. The blue dots are inside the 90% credible interval around the point giving the minimum χ 2 value. (red 'x').  \n<!-- image -->  \nFIG. 6. Normalised histogram of the distribution of M 1 , of simulated BBHs detectable with a SNR > 8 . The underlying simulated distribution is that of a power-law model with α = 2 . 26 , and β = -0 . 26 (blue columns with red best-fit normalizations). Superposed to that, is the histogram of the M 1 distribution from the LVK BBHs with Poisson error-bars (grey columns with back errors).  \n<!-- image -->  \nconverted into the 1 σ error. We test 12 discrete choices for the peak value of the Gaussian distribution µ . The corresponding standard deviations are given in TableIV, for each choice of µ .  \nThe best fit values for the χ 2 test for the PP model are  \nachieved for α = 2 . 31 +0 . 31 -0 . 43 , β = -0 . 45 +0 . 71 -0 . 55 and a Gaussian peak at 40 M ⊙ . For that value of µ , the corresponding standard deviation is 6 . 38 M ⊙ . We get a χ 2 M 1 = 16 . 95 , χ 2 q = 26 . 6 and a χ 2 tot = 43 . 5 . In this case, the χ 2 per degree of freedom is χ 2 d.o.f. = 1 . 89 . This is a significant  \nFIG. 7. Normalized histogram of the distribution of q of simulated BBHs, detectable with an SNR > 8 . In our simulation the BBHs follow a power-law model with α = 2 . 26 , and β = -0 . 26 (blue columns with red normalizations). As in Fig. 6, superposed is the q distribution from the LVK BBHs (grey histogram with black Poisson errors).  \n<!-- image -->  \nTABLE IV. The location parameter values probed for Gaussian component of the PP model and their corresponding standard deviations (see text for details).  \nimprovement to the simple PL model's best fit results.  \nWith the quoted best-fit choices, about one in five of simulated events are drawn from the Gaussian component. We note that the χ 2 from the fitting of M 1 alone, gives a result very close to 1 per d.o.f. . For the PP case, it is the q -distribution that contributes the most to the χ 2 tot . This seems to suggest that modeling the BBH population as either a simple power-law for both M 1 and q , or that in combination with a Gaussian component for M 1 , is not sufficient to explain the observed q -distribution. To explain the BBH observations, we require a mechanism that favours values of around 0.6. A power-law for q , will consistently lead to underestimating  \nthe number of counts in the the q -bin centred at 0.6, or overestimate the the number of counts in the last bins.  \nTable V, includes our best-fit choices for the PP model, along with χ 2 values for specific parameter choices as fixing β to 0 and 1 or α to 3.5 and marginalizing over the other parameters. We report the χ 2 values and the best-fit values for the other free parameters.  \nFigs. 8 and 9, show the histograms of the simulated BHs from a power-law + Peak model, with the best-fit parameters compared to the histogram from the LVK data. The M 1 peak at around 40 M ⊙ seems to be insufficient in its amplitude when only looking at Fig. 8. However, we remind the reader that when adding this peak component, we use of same normalization factor for both M 1 and q (as it is the same binaries). The Peak component adds binaries almost entirely in the last q ≃ 1 bin. Any attempt to make the Gaussian peak higher to fit the M 1 data better, also raises the number of counts in the last bin for the q -histogram.  \nIn Figs. 10 and 11, we show the area in parameter space around the best-fit values that is within 90% confidence interval for the PP model. The 90% confidence interval range for the PP model overlaps with the 90% confidence interval of the power-law model. Some values of α and β give a χ 2 inside the 90% confidence interval for the PP model, regardless of the normalization of the Gaussian peak (even if very small). In these cases the PP fit ends up being effectively a power-law fit.  \nFIG. 8. Normalized histogram of the distribution of M 1 of binary black holes with an SNR > 8 that follow a PP distribution (blue columns with red normalizations). We plot here the case of α = 2 . 31 , β = -0 . 45 and µ = 40 M ⊙ . The corresponding standard deviation for a peak at 40 M ⊙ is σ = 6 . 38 M ⊙ . The simulated BBHs are superposed to the histogram of the M 1 distribution from the LVK ones with Poisson error bars (grey columns).  \n<!-- image -->  \nFIG. 9. Normalized histogram of the distribution of q of binary black holes with an SNR > 8 that follow a power-law distribution (blue columns with red normalizations). We show results for α = 2 . 31 , β = -0 . 45 and µ = 40 M ⊙ . We also show the histogram of the q distribution from the LVK BBHs with Poisson error bars (grey columns).  \n<!-- image -->", 'V. DISCUSSION AND CONCLUSION': "In this paper, we present a way of analyzing the LIGOVirgo-KAGRA observations in terms of the mass and  \nmass ratio spectra of BH binaries that it has detected. Taking each of the detected merger events of BBHs with a SNR > 8 , we construct M 1 and q distributions as we  \nFIG. 10. The α -β parameter space for the power-law model probed using the BBHs from the LVK observations for different choices of µ . The blue dots are inside the 90% credible interval around the point giving the minimum χ 2 value.  \n<!-- image -->  \nhave described in Section II. We then study how these distributions can be explained by different populations of black holes.  \nWe simulate a population of stellar-origin BHs. These BHs have a mass distribution that parametrically resembles the mass-distribution of stars at their zero age as we describe in Section III. We test a wide range of assumptions on those BHs' mass spectra. In each case we simulate a population of BBHs following specific M 1 -, q -and z -distributions and evaluate how many of these binaries would be detected by the LVK collaboration, using the reported LIGO noise curves for the Livingston and Hanford observatories. We also allow for a second component of BBHs, where their mass spectrum would be  \napproximated by a Gaussian distribution. Such a population could stem either from PBHs with a narrow massspectrum centred at the stellar mass range or from massive BHs that themselves are the product of earlier generation of BBH mergers.  \nWe find that the combination of both populations can much better explain the LVK observations (see our results in Section IV). The best-fit choice to the M 1 - and q -distributions that we study require for a stellar-origin BH population with a very similar spectrum to that expected for massive stars at their zero age, in combination with a Gaussian distribution centered at 40 M ⊙ . Our result for the need of both components to explain the observations is in agreement with the claim made by Ref. [5],  \nFIG. 11. As in Fig. 10, the α -β parameter space for the power-law model probed using the BBHs from the LVK observations for different choices of µ . The blue dots are inside the 90% credible interval around the point giving the minimum χ 2 value. (red 'x').  \n<!-- image -->  \nbut derived independently as we do not use information that could allow us to estimate rates of occurrence of individual merger events. However, we find somewhat different best fit values for the stellar-origin BH population's distribution properties.  \nOur method allows for a rapid analysis of the LVK observations as more results are soon to come from the O4 and future runs. In particular, it can be used to place limits on the merger rate of PBH binaries. We leave such an analysis for future work. It can also be used to study the conditions that would be required in dense stellar environments such that a second component of BBHs is roughly approximated by a Gaussian term in the derived  \nM 1 -distribution. We note that the models used in this paper are relatively simple, as the main aim of this work is to use the limited number of LVK BBHs as a showcase of how to address those questions (see also Ref. [19]). We leave for future work how to connect such an analysis to simulations of dense stellar environments, that would give us a prediction for the mass-properties of BHs coming from multiple/hierarchical mergers in them. This is of particular importance in relation to fitting the derived q -distribution. We noticed that neither a simple powerlaw BH mass-distribution nor a Gaussian component can optimally explain it.  \nAs a last note, we point out that for future more in-  \nTABLE V. Results of the χ 2 test after comparing simulated populations from the PP model to the LVK data. The number of degrees of freedom for this test is 23. The first line shows the best fit parameters and χ 2 results, the second line the best fit results for β = 0 , the third line results for β = 1 , and last line results for the parameters reported by the LVK collaboration that yield the best fit for the merger rate [5].  \ndepth studies, LVK-detectors noise curves that account for the change in their sensitivity during a given observing run would be necessary. In this study we had to rely only on one noise curve per LIGO detector, derived during the early O3 run, considering it valid for the entire run. This hypothesis is a rough approximation, that can be easily improved on for future runs, if similar noise curves as those given in [15], could be provided for shorter time intervals during the observation period.  \n- [1] B. P. Abbott et al. , 'Observation of gravitational waves from a binary black hole merger,' Physical Review Letters 116 (2016), 10.1103/physrevlett.116.061102.\n- [2] T. Accadia et al. , 'Virgo: a laser interferometer to detect gravitational waves,' Journal of Instrumentation 7 , P03012 (2012).\n- [3] T. Akutsu et al. , 'Overview of KAGRA: Detector design and construction history,' Progress of Theoretical and Experimental Physics 2021 , 05A101 (2020), https://academic.oup.com/ptep/articlepdf/2021/5/05A101/37974994/ptaa125.pdf.\n- [4] B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), 'Gwtc-1: A gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs,' Phys. Rev. X 9 , 031040 (2019).\n- [5] R. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration), 'Population of merging compact binaries inferred using gravitational waves through gwtc-3,' Phys. Rev. X 13 , 011048 (2023).\n- [6] R. Abbott et al. , 'Open data from the third observing run of ligo, virgo, kagra, and geo,' The Astrophysical Journal Supplement Series 267 , 29 (2023).\n- [7] J. Aasi et al. , 'Advanced ligo,' Classical and Quantum Gravity 32 , 074001 (2015).", 'ACKNOWLEDGMENTS': "IC and MEB are supported by the National Science Foundation, under grant PHY-2207912. We would like to thank J. Hodaj and J. Decoste for useful discussions.  \nThis research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation, as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.  \n- [8] P. Kroupa, 'On the variation of the initial mass function,' Monthly Notices of the Royal Astronomical Society 322 , 231-246 (2001).\n- [9] Alexandre Arbey, 'Primordial black holes, a small review,' (2024), arXiv:2405.08624 [gr-qc].\n- [10] A. Azzalini, The Skew-Normal and Related Families , 1st ed. (Cambridge University Press, 2014).\n- [11] Simeon Bird, Ilias Cholis, Julian B. Muñoz, Yacine AliHaïmoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Raccanelli, and Adam G. Riess, 'Did LIGO detect dark matter?' Phys. Rev. Lett. 116 , 201301 (2016), arXiv:1603.00464 [astro-ph.CO].\n- [12] Misao Sasaki, Teruaki Suyama, Takahiro Tanaka, and Shuichiro Yokoyama, 'Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914,' Phys. Rev. Lett. 117 , 061101 (2016), [Erratum: Phys.Rev.Lett. 121, 059901 (2018)], arXiv:1603.08338 [astro-ph.CO].\n- [13] Yacine Ali-Haïmoud, Ely D. Kovetz, and Marc Kamionkowski, 'Merger rate of primordial blackhole binaries,' Phys. Rev. D 96 , 123523 (2017), arXiv:1709.06576 [astro-ph.CO].\n- [14] Muhsin Aljaf and Ilias Cholis, 'Simulating Binary Primordial Black Hole Mergers in Dark Matter Halos,' (2024), arXiv:2408.06515 [astro-ph.GA].\n- [15] LVK, 'Noise curves used for simulations in the update of  \nthe observing scenarios paper,' (2020).  \n- [16] Éanna É. Flanagan and Scott A. Hughes, 'Measuring gravitational waves from binary black hole coalescences. i. signal to noise for inspiral, merger, and ringdown,' Phys. Rev. D 57 , 4535-4565 (1998).\n- [17] Scott Dodelson and Fabian Schmidt, Modern Cosmology (Academic Press, Elsevier Science, 2020).\n- [18] N. Aghanim et al. (Planck), 'Planck 2018 results. VI.  \n- Cosmological parameters,' Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].\n- [19] Jordan Flitter, Julian B. Muñoz, and Ely D. Kovetz, 'Outliers in the LIGO black hole mass function from coagulation in dense clusters,' Mon. Not. Roy. Astron. Soc. 507 , 743-760 (2021), arXiv:2008.10389 [astro-ph.HE]."}