GW190521 and the GWTC-1 Events: Implication on the Black Hole Mass Function of Coalescing Binary Black Hole Systems
Yuan-Zhu Wang, Shao-Peng Tang, Yun-Feng Liang, Ming-Zhe Han, Xiang Li, Zhi-Ping Jin, Yi-Zhong Fan, Da-Ming Wei
DDraft version September 9, 2020
Typeset using L A TEX default style in AASTeX63
GW190521 and the GWTC-1 Events: Implication on the Black Hole Mass Function of CoalescingBinary Black Hole Systems
Yuan-Zhu Wang, Shao-Peng Tang,
1, 2
Yun-Feng Liang, Ming-Zhe Han,
1, 2
Xiang Li,
1, 2
Zhi-Ping Jin,
1, 2
Yi-Zhong Fan,
1, 2 and Da-Ming Wei
1, 2 Key Laboratory of dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing, 210033,People’s Republic of China. School of Astronomy and Space Science, University of Science and Technology of China, Hefei, Anhui 230026, People’s Republic ofChina. Guangxi Key Laboratory for the Relativistic Astrophysics, Department of Physics, Guangxi University, Nanning 530004, People’sRepublic of China.
ABSTRACTWith the black hole mass function (BHMF; assuming an exponential cutoff at a mass of ∼ M (cid:12) ) ofcoalescing binary black hole systems constructed with the events detected in the O1 run of the advancedLIGO/Virgo network, Liang et al. (2017) predicted that the birth of the lightest intermediate massblack holes (LIMBHs; with a final mass of (cid:38) M (cid:12) ) is very likely to be caught by the advancedLIGO/Virgo detectors in their O3 run. The O1 and O2 observation run data, however, strongly favora cutoff of the BHMF much sharper than the exponential one. In this work we show that a power-lawfunction followed by a sudden drop at ∼ M (cid:12) by a factor of ∼ a few tens and then a new power-lawcomponent extending to ≥ M (cid:12) are consistent with the O1 and O2 observation run data. Withthis new BHMF, quite a few LIMBH events can be detected in the O3 observation run of advancedLIGO/Virgo. The first LIMBH born in GW190521, an event detected in the early stage of the O3 runof advanced LIGO/Virgo network, provides additional motivation for our hypothesis. INTRODUCTIONThe stellar-mass (i.e., ≤ M (cid:12) ) black holes (BHs) are expected to form when very massive stars collapse at the endof their life cycle (Woosley et al. 2002) or merger/accretion from lighter compact objects. Traditionally, the massesof the BHs can be measured if they are within the binary systems. Indeed, a few dozens of BHs of stellar mass havealready been detected in X-ray binaries within the Milky Way and some nearby galaxies (Farr et al. 2011). On 2015September 14, the successful detection of a gravitational wave (GW) signal from the merger of a binary black hole(BBH) by Advanced LIGO (aLIGO; Abbott et al. 2016a) opens a brand-new window into observing the universe. Sofar, about two dozens of BBH GW events have been reported and even more events are expected to be released soon.In the population studies, the black hole mass functions (BHMFs) in different binary systems are one of the keysubjects since such information can help us reveal the stellar evolution physics and/or the origin of these systems.For instance, there could be the suppression on the number of the low mass BHs and the absence of very massivestellar-mass BHs with masses above ∼ − M (cid:12) , as expected in the modern supernova explosion models (e.g.Belczynski et al. 2016; Woosley et al. 2020). In principle, the BHMFs for different binary systems (for instance, theX-ray binaries, the BBH merger events and the neutron star-BH merger events) may show different characters (Tanget al. 2020). It is therefore essential to construct them with the observation data and some interesting features maybe revealed. Before the recent report of a BH candidate with a mass of 3 . +2 . − . M (cid:12) identified in the binary system2MASS J05215658+4359220 (Thompson et al. 2019), the BHs identified in the X-ray binaries have the lightest massof ∼ M (cid:12) , which is much heavier than the upper limit of neutron stars ( ¨Ozel et al. 2010) and suggests the presenceof a mass gap between the neutron stars and the BHs. Since the discovery of the BBH merger event, the GW datahave been extensively adopted to reconstruct the BHMF in such a specific group of objects (e.g., Abbott et al. 2016b;Fishbach & Holz 2017; Kovetz et al. 2017; Bai et al. 2018; Abbott et al. 2019a). While for the O1 events, the BHMF Corresponding authors: [email protected] (YZF) and [email protected] (DMW) a r X i v : . [ a s t r o - ph . H E ] S e p was only loosely constrained due to the rather small sample. Interestingly, assuming an exponential cutoff of the massfunction at ∼ M (cid:12) , the predicted birth rate of the lightest intermediate mass BH (i.e., M f ≥ M (cid:12) , where M f isthe final mass of the new BH formed in the merger) is high and the detection prospect is found to be quite promisingin the O3 run of the advanced LIGO/Virgo (Liang et al. 2017). Such a result is indeed encouraging. However, theanalyses of both O1 and O2 BBH events find a “termination” of the BHMF at ∼ M (cid:12) (Fishbach & Holz 2017; Baiet al. 2018; Roulet & Zaldarriaga 2019; Talbot & Thrane 2018; Abbott et al. 2019a; Wysocki et al. 2019), which hasbeen taken as the evidence for the (pulsational) pair instability processes of the supernova explosions. Consequently,the detection prospect of IMBH is unpromising. However, heavier stellar-mass BHs certainly exist in the Universe (forinstance in GWTC-1 there were BHs formed with M f ∼ − M (cid:12) (Abbott et al. 2019b); see Sec.2.1 for extendeddiscussion of this issue) and they may also be involved in the BBH mergers. We hence speculate the presence ofa new component in the “high” mass range but there is a sudden drop in the BHMF at the mass of m max due tothe supernova explosion mechanism limit. Note that in the final stage of preparing for this work , the LIGO/Virgocollaboration reported the first IMBH formed in GW190521, a merger of two BHs with the masses of ∼ − M (cid:12) (The LIGO Scientific Collaboration et al. 2020a). These authors also pointed out that the detection of an IMBH in theearly stage of O3 run of advanced LIGO/Virgo is in tension with the O1-O2 data (The LIGO Scientific Collaborationet al. 2020b). As we will demonstrate in this work, such a tension can be solved in our new BHMF model.This work is organized as the follows: in Sec.2, we introduce the models for parameterizing the mass distributions, thelikelihood of hierarchical inference, and the selection effects. We report our results in Sec.3 and discuss the implicationof predicting the IMBH in O3 run in Sec.4. And Sec.5 is our Conclusion and Discussion. DATA, MODELS AND SELECTION EFFECTS2.1.
Parameterized Mass Spectrums
We first introduce the BHMF models adopted in Abbott et al. (2019a, hereafter LVC19), where three models (namelyModel A, B, and C) are discussed. Model A and B share the same formula that are constructed by a power-law (PL)distribution with a hard cutoff at the high mass end, p ( m , m | Λ) ∝ N ( m ) m − α q β q if m min ≤ m ≤ m ≤ m max , N ( m ) is chosen so that the marginal distribution is a power-law in m . In Model A, m min and β q are fixed to5 and 0 respectively, while in Model B all of the parameters Λ = ( α, m max , m min , β q ) are allowed to vary. Since ModelA of LVC19 is a simplified version of Model B, and there is only mild evidence that the data is better fitted by thecomplex Model C (Abbott et al. 2019a), we only take their Model B as the fiducial model and its name is unchangedto keep the “consistency” among the literature. To further examine the properties of the mass function around thepossible high mass gap, we propose two new BHMFs, i.e., Model Bcut and Model 2PL. Model Bcut is analogue toModel B, but ends smoothly, which is described by p ( m , m | Λ) ∝ m − α exp[ − ( m /m cut ) k ] q β q if m min ≤ m ≤ m ≤ m max , m cut is governed by the “cut-off index” k . For k ≤
1, there isa good fraction of BHs with masses above m cut . While for k (cid:29)
1, our case reduces to Model B of LVC19. And theModel 2PL consists of two power-law segments differing greatly in their magnitudes, which reads p ( m , m | Λ) = (1 − lg F ) A m − α q β q , if m min ≤ m ≤ m ≤ m max , lg F A m − α q β q , if m max ≤ m ≤ m edge and m min ≤ m ≤ m , (3)where A and A represent the normalization factor of the first and second segment respectively, and F is the fractionof BHs between the second segment and the whole population. This model is motivated by some astrophysical theoriesin which the merging black holes can have different origins. For instance, some BBH systems, in particular those residewithin the accretion disks of the Active Galactic Nuclei (AGN) (Bartos et al. 2017; Stone et al. 2017; McKernan et
10 100 M / M P D F
10 100 M / M P D F Figure 1.
Representative mass distributions of the models. Left panel: the Bcut model, the blue solid line has parametersof ( m min , α, m cut , k ) = (8 , . , , k = 5; right panel: the 2PLmodel, the blue solid line has parameters of ( m min , m max , m edge , α , α , lg F ) = (8 , , , . , , − , , , . , , − . Table 1.
Priors of the Parameters for Different ModelsParameters/Models Model B Model Bcut Model 2PL α or α [-4, 12] [-4, 12] [-4, 12] m max [ M (cid:12) ] [20, 100] N/A [20, 100] m min [ M (cid:12) ] [5, 10] [5, 10] [5, 10] β q or β q , [-4, 12] [-4, 12] [-4, 12] m cut [ M (cid:12) ] N/A [20, 100] N/A k N/A [0, 12] N/A α N/A N/A 1 and 2 β q , N/A N/A 0 m edge [ M (cid:12) ] N/A N/A 100 and 150lg F N/A N/A fixed values1. All of the parameters are uniformly distributed in their domains. We have also used a narrower prior volume to calculatethe Bayes factors discussed in Sec.3, which is m max ∈ U(30 , m min ∈ U(6 , α ∈ U(0 ,
12) and β ∈ U(0 , F in Model 2PL are fixed to − , − . , − , − . , − . , − . − al. 2018) may be able to accrete material from the surrounding and have higher masses. It has also been proposedthat heavy black holes can be dynamically formed by lighter objects through hierarchical merger or through runawaycollisions (Rodriguez et al. 2015; Antonini & Rasio 2016; Mapelli 2016; Fishbach et al. 2017). For example, theBBH system of GW170729 may be formed through hierarchical mergers in the migration traps that developed in theaccretion disks of AGN (Yang et al. 2019). So it is reasonable to expect an extended tail of the mass distribution of theBHs or a population of “high” mass objects following a sudden drop of the BHMF at m max . Since only GW170729 hasa primary mass with inferred median value exceeding 40 M (cid:12) during the O1 and O2 runs, we expect that the secondsegment cannot be well constrained from current data. Thus we only take the fixed values of β = 0, α = 1 (or 2), m edge = 100 (or 150), and a series of lg F ranging from − − The Likelihood and Selection Effects
The likelihood for hierarchical inference is constructed based on Thrane & Talbot (2019). For a series of measurementsof N events (cid:126)d , the likelihood for the hyperparameter Λ can be inferred via L ( (cid:126)d | Λ) = N (cid:89) i Z ∅ ( d i ) n i n i (cid:88) k π ( θ ki | Λ) π ( θ ki | ∅ ) . (4)In Eq.(4), the n i posterior samples for the i -th event, the evidence Z ∅ ( d i ) as well as the default prior π ( θ k | ∅ ) areavailable for the released O1 and O2 events . Taking the selection effect into account (Thrane & Talbot 2019), thelikelihood can be revised to L ( d, N | Λ , det ) = p det (Λ | N ) L ( d, N | Λ , R ) if ρ ≥ ρ th , p det (Λ | N ) is the probability of detecting N events, which is dependent of population hyperparameter Λ. If weassume that the Poisson-distributed rate R has a uniform-in-log prior, the marginalized p det (Λ | N ) (Thrane & Talbot2019) can be expressed as p det (Λ | N ) ∝ (cid:18) V (Λ) V tot (cid:19) N = f (Λ) N , (6)where V tot is the total spacetime volume and V (Λ) is the visible volume of the population described by Λ, and f (Λ)is the fraction of detectable sources.The f (Λ) for specific hyperparameters Λ can be obtained by carrying out injection campaign into the searchingpipeline, however it is very computational expensive especially for models with high dimensional parameter space.Abbott et al. (2016b) points out that the searching process can be approximated by a semi-analytic method, whichassumes that a source is detectable when the signal-to-noise ratio (SNR) produced in a single detector ≥ ρ ( θ ) = w (R . A ., dec ., ι, ψ ) ρ oo ( m , m , z ) , (7)where ρ oo is the SNR of an optimally oriented (i.e., faced-on/directly overhead) source with component masses m , at redshift z , and w is the angular factor (Finn & Chernoff 1993) accounting for different sky position and orientation.The distribution of w is available in (Finn & Chernoff 1993), thus f (Λ) can be evaluated (Vitale 2020) by f (Λ) = (cid:90) d m d m d z CCDF w (cid:18) ρ thr ρ oo (cid:19) π ( m , m , z | Λ) , (8)where CCDF w denotes the complementary cumulative distribution function of w , and ρ thr is chosen to be 8. Toimplement fast generating of ρ oo ( m , m , z ), we first calculate the optimal SNRs for a series of sources with detectorframe masses ( m det1 , m det2 ) (yielded from a regular grid) and a reference luminosity distance d refL . The IMRPhenomPv2 waveform model (Hannam et al. 2014; Khan et al. 2019) and the aLIGOEarlyHighSensitivity
PSD taken from
PyCBC (a good approximation to the PSD during the O1 and O2 runs (Fishbach & Holz 2017)) are used in the calculation.Then ρ oo ( m , m , z ) can be derived by interpolating the results of the above calculation and rescaling to a specificdistance, ρ oo ( m , m , z ) = ρ ref ( m det1 , m det2 , d refL ) × d refL d L ( z ) . (9)Therefore, the f (Λ) in Eq.(8) can be obtained by using Monte-Carlo integration. Finally, we use the python package Bilby and
PyMultinest sampler to obtain the Bayesian evidence and posteriors of the parameters for each models. RESULTSIn this section, we show the results of hierarchical inference described above, and the credible intervals are reportedwith 90% uncertainties unless otherwise stated. The constraints on the parameters of each models are reported in https://dcc.ligo.org/LIGO-P2000193/public Table 2.
Summary of Constraints on the Parameters Considered in Tab.1Parameters/Models Model B Model Bcut Model 2PL α or α . +1 . − . . +1 . − . ( − . , . m max [ M (cid:12) ] 42 . +16 . − . N/A (31 . , . m min [ M (cid:12) ] 7 . +1 . − . . +1 . − . (5 . , . β q or β q , . +4 . − . . +4 . − . (0 . , . m cut [ M (cid:12) ] N/A 41 . +22 . − . N/A k N/A 7 . +3 . − . N/A α N/A N/A fixed β q , N/A N/A 0 m edge [ M (cid:12) ] N/A N/A fixedlg F N/A N/A fixedThe table shows the 90% credible intervals of posterior distributions. The intervals for Model 2PL are obtained by combiningthe results of all cases.
Tab.2. The result of Model B in our inference consist with that of Abbott et al. (2019a). The common parametersbetween Model B and Model Bcut are consistent with each other. The “cut-off index” k in Model Bcut has a largemedian value of 7.48, which indicates that the BHMF, if its main component has a power-law shape, does have a veryhard cut-off after a certain limit of mass. In Fig.2, we show the dependencies of the inferred α , m max , m min and β on the choice of fixed values for Model 2PL. One can find from Fig.2 that the constraints on the parameters of thefirst segment are insensitive to the choice of α and m edge , while the constraints on α and m max are sensitive to thechoice of lg F .Since Model 2PL can be regarded as Model B followed by an additional weak component, we investigate the possibilityof the existence of such a weak component, by computing the Bayes factors BF B2PL between the two models. Weinterpret Bayes factors of < −
10 as moderate, 10 −
30 as strong, 30 −
100 as very strong, and > B2PL is most sensitive to the choice of lg F . The BFs in all groups are still lowerthan 10 even the lg F has reached − . R = 0 . F reaching − . R = 0 . the evidence of excluding an extra segmentof BHMF after m max that extends to much higher mass is not strong enough, and this component is unlikely to occupy > of all primary BHs, but the fraction can still be as high as . The resulting BF is affected by the choice of prior. In the analysis above, we have used relatively broad and un-informative priors for the parameters (see Tab.1 for details), and we expect that the information used in our modelcomparison mainly comes from the data. To study the influence from the choice of prior, we narrow down the volumesof the common priors for the four models (also shown in Table.1) to recalculate the BFs above, and find that theinfluence is relatively small: for example, the BF
B2PL with ( α , m edge ) = (1 , B2PL with ( α , m edge , lg F ) = (1 , , − .
5) is 4.88 for the default and 3.61 forthe narrower prior volume (note that the BF derived from the nest sampling also has its uncertainty, which is ∼ . ∼ .
05 for the default and narrower prior volume respectively). As a consequence, our conclusion will still holdfor un-informative priors with reasonably broader/narrower volumes. THE CHANCE OF HEARING THE BIRTH OF IMBH IN O3In the above analysis, we showed that the fraction of high mass ( (cid:38) M (cid:12) ) BHs can be as high as a few percents.Since the ground-based GW detectors is most sensitive to the BBHs with total mass ∼ − M (cid:12) , the chance ofdetecting IMBHs born from mergers can be promising. Previously, we have already studied this topic in Liang et al.(2017). However, since there were only a few BBH events being reported at that time, the BHMF took there needsto be modified now. Based on the O1-O2 data, we derive the expected detections and compare it with the current lg F lg F m m a x lg F m m i n lg F Figure 2.
Dependency of the inferred parameters of Model 2PL on the choice of fixed parameters. The medians of theparameters are marked with filled squares, and the dashed lines show the 90% credible intervals. Different choices of α and m edge are represented by different colors. Blue: α = 1 , m edge = 100; green: α = 2 , m edge = 100; red: α = 1 , m edge = 150;orange: α = 2 , m edge = 150. lg F B F B P L = 1, m edge = 100 = 2, m edge = 100 = 1, m edge = 150 = 2, m edge = 150 Figure 3.
The Bayes factor between Model B and Model 2PL with different fix parameters of the second segment. The twogray dashed lines mark the BF
B2PL of 10 and 3, above which the evidence of “Model B is more preferred by the data” is strongand moderate respectively.
20 40 60 80 100 M /M ⊙ m min = 7⊙61m max = 37⊙45α = 0⊙89β = 6⊙85
20 40 60 80 100 120 140 M /M ⊙ m min = 7⊙72m max = 36⊙85α = 0⊙89β = 7⊙04
20 40 60 80 100 M /M ⊙ m min = 7⊙94m max = 38⊙65α = 1⊙22β = 6⊙63
20 40 60 80 100 120 140 M /M ⊙ m min = 8⊙00m max = 40⊙34α = 1⊙48β = 6⊙92 Figure 4.
The observed primary (blue bars) and secondary (orange bars) mass distribution in LIGO/Virgo’s O3 run predictedby Model 2PL with fix parameters. The parameters of the 2nd segment are selected by the Bayes factor criterion discussedin Fig.3, and represent the most optimistic cases in our analysis for observing IMBH. The cases in the upper row are selectedby BF
B2PL <
10, and ∼
16% of their remnants are IMBHs; the cases in the lower row are selected by BF
B2PL <
3, and ∼ m edge = 100 , α = 1 and lg F = − .
2; top right: the case with m edge = 150 , α = 2 and lg F = − .
2; bottom left: the case with m edge = 100 , α = 2 and lg F = − .
5; bottom right: the casewith m edge = 150 , α = 1 and lg F = − .
8. The parameters of the first segment is the one with the highest likelihood in theBayes inference. released information of O3. The Gravitational-Wave Candidate Event Database (GraceDB ) reports a total numberof 56 detection candidates, among which 38 events have the probability (cid:38)
90% to be BBHs (and 4 additional eventshave probability ≥
95% lie within the lower mass gap). GW190521 is confirmed to be a merger yielding IMBH, whilethe identification and properties of other BBH candidates are yet to be published. In the following we assume these38 candidates are all BBHs.Among all the parameter sets presented in Fig.3, lg F can be as large as − . α , m edge ) = (2 , α , m edge ) = (1 , α , m edge ) =(2 , m , m ) are drawn frompopulation described by the parameters with the highest likelihood in the inference. Then we count the fraction ofmergers with SNR ≥ aLIGOLateLowSensitivity PSD in
PyCBC ) and M tot > / . M (cid:12) (consideringapproximately 5% of the total gravitational mass is radiated by GW), and find that the fraction of mergers formingIMBH can be as high as ∼
16% for both groups, which means that 6 out of the “actual” 38 detections have IMBHremnants. A more conservative case may be accessed by choosing the parameter sets with BF <
3. By applying thiscriteria, the most optimistic expectation on the fraction of mergers that form IMBH is ∼
7% (for the BH populationwith ( α , m edge , lg F ) = (2 , , − .
5) and ( α , m edge , lg F ) = (1 , , − . https://gracedb.ligo.org/superevents/public Fig.4. Note that these are the expected mean values of IMBH detections, and due to poisson fluctuation the observednumbers have variances equal to the mean values. It is worthy of pointing out that although the above expectationsfor IMBH detections are similar for m edge = 100 and m edge = 150, the population with α = 1, m edge = 150 andlg F = − . <
3) predicts 3% (1 out of 38) of its mergers having m > M (cid:12) . This implies that we havethe chance to clarify whether a single pre-merger star can be the lightest IMBH (with a mass ≥ M (cid:12) ) in O3. CONCLUSION AND DISCUSSIONIn this work, we have studied the BHMFs based on the 10 BBHs observed in LIGO’s O1-O2 runs, and consideredthree models presented in Tab.1. By performing Bayesian hierarchical inference, we find that the data do prefer avery sharp cut-off at the mass ∼ M (cid:12) in the BHMF. However, there is still a good fraction of parameter space inwhich the model with an extra power-law segment can not be ruled our with strong evidence by the single power-lawmodel. The fraction of the extra segment to the overall distribution can be high up to 3%. Such a new segment cangive rise to rather exciting results in LIGO/Virgo’s O3 and future observation runs. The fraction of mergers that formIMBHs can be as large as ∼
16% in the most optimistic case and can be ∼
7% in a more conservative case. Thus, thediscovery of GW190521 is within our expectation.There are space of improvement in the future works. In our analysis, we only include the information from themass distributions inferred from each event. The inclusion of spin data (although additional thoughts are needed toconstruct the spin model) would enhance the ability for model comparison in the inference. For example, the expectedspins in the second segment of Model 2PL could be larger and more isotropic if we assume the BBHs in the firstsegment is formed by field binary stars and the second segment is formed by dynamical capture. It is worthy to notethat the most heavy BBH GW170729 (with source frame primary mass of ∼ M (cid:12) ) has χ eff ∼ .
37 (Abbott et al.2019b), while the GW190521 has χ p ∼ . ≥ M (cid:12) , the BHMF may show another interesting structure. Onthe one hand, the supernova explosions can directly produce such massive black holes (Woosley et al. 2002). On theother hand, the mergers of some black hole binaries can also produce IMBHs with the masses ≥ M (cid:12) , as alreadyobserved in GW190521 (The LIGO Scientific Collaboration et al. 2020a). All these objects could be involved in theBBH mergers and might be detected by the advanced LIGO/Virgo in the next decade. Software:
Bilby (Ashton et al. 2019, version 0.6.9, ascl:1901.011, https://git.ligo.org/lscsoft/bilby/), PyCBC (ThePyCBC Team 2018, version 1.13.6, ascl:1805.030, https://github.com/gwastro/pycbc), PyMultiNest (Buchner 2016,version 2.6, ascl:1606.005, https://github.com/JohannesBuchner/PyMultiNest)REFERENCES