Black Hole Mass Function of Coalescing Binary Black Hole Systems: Is There a Pulsational Pair Instability Mass Cutoff?

We analyze the LIGO/Virgo GWTC-2 catalog to study the primary mass distribution of the merging black holes. We perform hierarchical Bayesian analysis, and examine whether the mass distribution has a sharp cutoff for primary black hole masses below $65 M_\odot$, as predicted in pulsational pair instability supernova model. We construct two empirical mass functions. One is a piece-wise function with two power-law segments jointed by a sudden drop. The other consists of a main truncated power-law component, a Gaussian component, and a third very massive component. Both models can reasonably fit the data and a sharp drop of the mass distribution is found at $\sim 50M_\odot$, suggesting that the majority of the observed black holes can be explained by the stellar evolution scenarios in which the pulsational pair-instability process takes place. On the other hand, the very massive sub-population, which accounts for at most several percents of the total, may be formed through hierarchical mergers or other processes.


INTRODUCTION
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. 2016) on September 14, 2015, opened a brand-new window into observing the Universe. So far, three observing runs (O1-O3) have been finished by LIGO and Virgo, and the data of several tens of events are released to the public, including two confirmed binary neutron star (BNS) events, 44 confident BBH events, and one neutron star−black hole (NSBH) candidate (Abbott et al. 2020c;The LIGO Scientific Collaboration et al. 2020, see however Han et al. (2020) for a dedicated investigation on the possible NSBH nature of GW190425 (Abbott et al. 2020a)).
In population studies, the black hole mass functions (BHMFs) in different binary systems are one of the key subjects since such information can help us reveal the stellar evolution physics and the origin of these systems. The BHMF of coalescing binary black hole systems can not be tightly constrained with the events detected in the O1 run of the Advanced LIGO. Assuming a power-law BHMF with an exponential cutoff at the mass of ∼ 40M , Liang et al. (2017) predicted that the birth of the lightest intermediate mass black hole (LIMBH, which has a final mass of ≥ 100 M ) is very likely to be caught by the Advanced LIGO/Virgo detectors in their O3 run. The data of O1 and O2 observing runs (see Abbott et al. (2019b) for the first Gravitational-Wave Transient Catalog (GWTC-1)), however, strongly favor an abrupt cutoff of the BHMF much sharper than the exponential one. Such a sharp cutoff in the mass spectrum is anticipated in the pulsational pair-instability supernovae (PPISN) scenarios (Fowler & Hoyle 1964;Belczynski et al. 2016;Woosley et al. 2020). In this case, the pre-merger BH can not have a mass substantially more massive than ∼ 40M and no LIMBHs formed in the BBH mergers are expected in the O3 run. However, some BBHs observed in O3a have masses much larger than those in GWTC-1. In particular, even the secondary mass of GW190521 is more massive than 45M at 90% credibility and an LIMBH was formed in this event (Abbott et al. 2020b). The LIGO Scientific Collaboration et al. (2020) further showed that the simple truncated power-law model is disfavored compared with other more complected models, and the mass spectrum must extend to masses much higher than 40M . The collaboration also claimed the detections of non-zero χ p and anti-aligned spins in the population. These new discoveries might indicate the presence of multiple formation channels, and some studies were carried out in order to investigate the origins and branch ratios of the sub-populations of BBHs (Hütsi et al. 2020;Kimball et al. 2020;De Luca et al. 2021).
In stellar evolution theories, the exact position for the lower edge of the mass gap produced by PPISN depends on the details of the evolution models. Previous studies typically placed the lower edge of the gap at ≤ 65M (Belczynski et al. 2016;Spera & Mapelli 2017;Woosley 2017;Stevenson et al. 2019;Mapelli et al. 2020). After the discovery of GW190521, some recent works addressed that the gap is most sensitive to the 12 C(α, γ) 16 O reaction rate, and it could be also affected by other factors like the evolution of H-rich envelope (Farmer et al. 2019;Costa et al. 2021). By adopting low 12 C(α, γ) 16 O rates, the maximum black hole mass below the gap can reach ∼ 90M (Farmer et al. 2019), or the mass gap is completely removed (Costa et al. 2021).
PPISNe and pair instability supernovae (PISNe) could affect all formation channels in which the BHs are of stellar origin, including but not limited to classic isolated binary evolution, dynamical capture in different environments, and formation in triple/quadruple systems (see Stevenson et al. (2019) and the references therein), so the presence of mass gaps are expected in these channels. BBHs formed from dense stellar environments, hierarchical mergers or primordial black holes could populate the mass gap, however, the branch ratios of these channels are generally considered to be small. Under such a recognition, in this work we assume that: (1) The majority of the observed BBHs are originated from stellar evolution, and their distribution can be described by a truncated power law; (2) The position of the hard cutoff of the power law is limited by the traditional prediction on the PPISN/PISN gap, i.e., ≤ 65M . Our main goal is to construct simple models that can be directly compared to the most preferred models in The LIGO Scientific Collaboration et al. (2020) under the framework of Bayesian inference (i.e., to clarify whether the presence of a sharp cutoff in the mass function at the mass ≤ 65M is still consistent with the data). The other purpose is to examine whether the data of GWTC-1 events and those obtained in O3a run can be reasonably interpreted within the same mass function model.
The rest of the paper is structured as follows: in Sec.2, we introduce the models for parameterizing the mass distributions, the likelihood of hierarchical inference, and the selection effects. We report our results and make comparison between models in Sec.3, and discuss about constraints on the sub-populations in Sec.4. Sec.5 is our Conclusion and Discussion.

Data Selection
In this population study, we include 34 BBH (m 1 > m 2 > 3M ) events observed in O3a with FAR < 1yr −1 , as well as the 10 BBHs presented in GWTC-1. This choice of events is in the same way as for The LIGO Scientific Collaboration et al. (2020). To test the robustness of our results, the analysis is performed using different subsets of the data: i.e., all of the selected events in GWTC-2, the 34 events detected in O3a alone, and the 10 events in GWTC-1 alone.
The samples from parameter estimation (PE) of individual event are taken from LIGO Public Document Database (see Sec.2.3 for details), and for each event, we use 1000 random draws of the primary and secondary mass pairs (m 1 , m 2 ) from the PE samples in the analysis. The detector frame component masses in the PE samples are transferred to the masses in source frame using "Planck15" cosmology in Astropy.

Parameterized Mass Spectra
In The LIGO Scientific Collaboration et al. (2020), the authors concluded that the primary mass distribution is more consistent with a broken power law, or a power law with a Gaussian feature. We also include these two models in our work to compare their preferences by the data with other models. For the BROKEN POWER LAW model (hereafter Model I), we keep its definition and parameter names identical to the descriptions in The LIGO Scientific Collaboration et al. (2020). For the POWER LAW + PEAK model (hereafter Model II), we modify the formula of primary mass distribution to π(m 1 | λ, α, m min , δ m , m max , µ m , σ m ) = (1 − λ)P (m 1 | α, m min , δ m , m max ) + λG (m 1 | m min , δ m , µ m , σ m ), with P (m 1 | α, m min , δ m , m max ) = A 1 P(m 1 | α, m min , m max )S(m 1 | m min , δ m ) (2) being a truncated power-law distribution P (with spectral index α, minimum mass m min , and maximum mass m max ) modulated by a smooth function S (see The LIGO Scientific Collaboration et al. (2020) for details), and being a Gaussian distribution (with mean µ m and standard deviation σ m ) modulated by the smooth function. A 1 and A 2 are constants to normalize the distributions, and their values are calculated numerically according to the model parameters. The difference between our expression of Model II and Eq.(B5) of The LIGO Scientific Collaboration et al. (2020) is that both P and G are firstly multiplied by the smooth term S before the normalization. The purpose of making this modification is to ensure that we reproduce Fig.1  As mentioned in Sec.1, the stellar evolution scenarios typically predict hard cutoff < 65M . If such a cutoff exists, other evolution channels are needed to explain events with higher masses, and the shape/magnitude of the overall primary mass distribution might change significantly after the cutoff. We first consider a relatively simple case, in which the mass distributions before and after the cutoff are shaped by the power-laws, but with different spectra indices and magnitudes. The corresponding formula for this case (hereafter Model III) is a piece-wise function with two segments, where α 1 and α 2 are the power law indices for the the segments before and after the cutoff mass m max at the lower edge of mass gap, respectively. F represents the ratio between the possibility densities of the two segments at m max . Motivated by the typical predictions about the mass gap, we restrict the prior on m max to be uniform with a maximum of 65M . We expect the second segment to have probability density much smaller than the first segment, so we adopt a log-uniform prior for F . To illustrate Model III, we present the representative distribution of the primary mass with arbitrary choice of parameter values in Fig.1. This model is motivated by some astrophysical theories in which the merging black holes can have different origins. For instance, some BBH systems, in particular those residing within the accretion disks of the Active Galactic Nuclei (AGN) (Bartos et al. 2017;Stone et al. 2017;McKernan et al. 2018) may be able to accrete material from the surrounding and have higher masses. It has also been proposed that heavy black holes can be dynamically formed by lighter objects through hierarchical mergers or through runaway collisions (Rodriguez et al. 2015;Antonini & Rasio 2016;Mapelli 2016;Fishbach et al. 2017). For example, the BBH system of GW170729 may be formed through hierarchical mergers in the migration traps that developed in the accretion disks of AGN (Yang et al. 2019). So it is reasonable to expect an extended tail of the mass distribution of the BHs or a population of "high" mass objects following a sudden drop of the BHMF at m max .
Both Abbott et al. (2019a) and The LIGO Scientific Collaboration et al. (2020) showed that an extra Gaussian component peaking at ∼ 30 − 40M might exist. We introduce Model IV, a modified version of the MULTI PEAK model in The LIGO Scientific Collaboration et al. (2020), to partially study the existence of different components and their influence to the overall shape of the spectrum. Model IV is expressed as where P and G are described in Eq.
(2) and Eq. (3), (1-λ) is the fraction of binaries in the main truncated power-law component, and λλ 1 is the fraction of binaries in the first modulated Gaussian, respectively. For Model IV, we also adopt the astrophysically motivated prior in which m max is uniformly distributed below 65M . In addition, comparing with the original MULTI PEAK model, in Model IV the lower bound for the prior on µ 2 is set to 30M , while the upper bound for the prior on σ 2 is set to 50M . These changes on the prior for the secondary Gaussian component increase its flexibility on describing the spectrum at high masses. The right panel of Fig.1 shows two extreme cases  The priors on the parameters listed above are all uniform of Model IV: when µ 2 m max and σ 2 10, the second modulated Gaussian actually represents a very shallow component that extends to high masses, while if µ 2 m max and σ 2 10, it represents a clear and narrow Gaussian component after m max .
For all of the four models described above, we use a conditional mass ratio (q) distribution that is consistent with Eq.(B8) of The LIGO Scientific Collaboration et al. (2020), leading to the inclusion of an additional free parameter β q in our inference. While different components in the models may have diverse mass ratio distributions, our work mainly focuses on the primary mass distribution, so we leave this issue to future studies. We summarize the parameters of Model I-IV, as well as their priors in Tab.1.

The Likelihood and Selection Effects
The likelihood for hierarchical Bayesian inference is constructed based on Poisson process. For a series of measurements of N obs events d, assuming a non-evolving merger rate R 0 , the likelihood for the hyper-parameters Λ (including where N = R 0 V 0 T obs is the expected number of mergers during the observation period T obs and within the astrophysical volume V 0 . Here we take V 0 T obs = 167.6 Gpc −3 yr −1 for O3 and 154.5 Gpc −3 yr −1 for O1-O2. In Eq.(6), the n i posterior samples for the i-th event, the evidence Z ∅ (d i ) as well as the default prior π(θ k | ∅) are available for both of GWTC-1 1 and GWTC-2 2 events. η(λ) is the detection efficiency for a particular λ, and we follow the procedures described in the SensitivityTutorial in LIGO Public Document Database to compute this quantity 3 . We use the same criteria that define the detectable events as The LIGO Scientific Collaboration et al. (2020), i.e., SNR > 8 for O1-O2 and FAR < 1/yr for O3. Finally, we use the python package Bilby and PyMultinest sampler to obtain the Bayesian evidence and posteriors of the hyper-parameters for each models.

MODELS COMPARISON
We compute Bayes factors between models to quantify their relative preferences by data. Tab.2 shows the Bayes factors B for each mass model relative to Model II, and in the context we interpret B of < 1/3 as moderate, < 1/30 as strong, and < 1/100 as decisive evidence for the first model is less favorable by the data compared with the second model (Model II) (Jeffreys 1998). For the analysis including all data, the Bayes factor between Model I and Model II is 0.13, which is consistent with the reported value of 0.12 between the BROKEN POWER LAW model and the POWER LAW+PEAK model in The LIGO Scientific Collaboration et al. (2020). Both Model III and Model IV has larger B comparing with Model I in the GWTC-2 analysis and O3a-only analysis, while they are less preferred by the data in the GWTC-1 analysis. This result is understandable, since Model III and Model IV are more flexible on describing the high mass spectrum and have larger prior volumes. GWTC-1 contains relatively lighter BHs on average, so Model III and Model IV suffer from Occam factor penalty in Bayes factor; the fraction of BHs with large masses increase significantly in O3a, hence the introducing of larger prior volumes allows the models to fit the data better, giving higher likelihoods; the outcome of GWTC-2 analysis can be regarded as the average over the GWTC-1 and O3a-only analyses. On the other hand, there is no strong evidence (B < 1/30) that one of the four models is better supported by the data against others. In general, one conclusion that can be made at this stage is that Model III and Model IV give comparable goodness of fit to the data compared with Model I and Model II respectively, indicating they are also acceptable approximations to the observed population. Since Model III and Model IV are not strongly disfavored by the GWTC-1 data, we conclude that there is no strong tension between the data collected in different runs. The table shows the median and 90% credible intervals of posterior distributions, inferred from the analysis using all 44 samples. The analysis performed with only O3a events gives consistent and slightly looser constraints on the parameters.

CONSTRAINTS OF THE POTENTIAL COMPONENTS
In this section, we mainly focus on the constraints obtained using all data (the GWTC-2 analysis). To present the constraints on each model, we summarize the median and 90 percent credible intervals of the hyper-parameters for the GWTC-2 analysis in Tab.3. The posterior distribution for the parameters of Model III and Model IV are also shown in Fig.2 and Fig.3 respectively. The inferred parameters for Model I and Model II in our work are consistent with the results in The LIGO Scientific Collaboration et al. (2020), despite the peak of the posterior distribution for λ in Model II is shifted to a higher value due to our modification on the formula of POWER LAW + PEAK model (see Eq.1 for details).
For Model III, the inferred parameters that describe the first power-law segment are in agreement with the ones in Model I. The maximum mass (m max ) of this segment is constrained to 47.13 +10.75 −8.75 M , above which the probability density function (PDF) of the primary mass spectrum falls by a factor of 0.01 − 0.72 (according to the credible interval of log 10 F ) and the power-law index α 2 changes to 2.54 +5.64 −5.63 . Interestingly, the position of m max , as well as the first power-law index α 1 of Model III is also in agreement with the inferred α and m max for the truncated power law model in previous study (see the results for Model B in Abbott et al. (2019a)) using GWTC-1 data. We also find that the constraint on m max is insensitive to the PPISNe motivated prior of m max ≤ 65M . By changing the prior for m max to 30M ≤ m max ≤ 95M , and fixing M edge = 100M , the resulting constraint is m max = 49.15 +22.26 −10.59 M . Together with other constrained parameters, we can infer that 0.1% − 5.8% of the primary BHs have masses larger than m max . We show the credible region of the mass spectrum for Model III in the left panel of Fig.4. Since α 2 is poorly constrained, we further fix it to 1, which represents a very shallow tail after m max , and the analysis using all data also gives an acceptable Bayes factor of B = 0.28 for Model III. We can see that by introducing an abrupt drop on the mass spectrum at m max , Model III allows a much shallower segment at high masses compared with Model I.
For Model IV, m max is less constrained. As shown in Fig.3 the posterior of m max is broadly distributed across the range of prior, and there is significant posterior support around the median of m max inferred from Model III. The result also shows a Gaussian component peaking at ∼ 33M , regardless of whether we use all events or O3a-only events in the analysis. Both the study in Abbott et al. (2019a) and the O3a-only analysis in our work have recovered consistent peaks for the Gaussian component, which enhances the evidence about the presence of this sub-population of black holes. The second Gaussian component, which contains 0.4% − 4.9% of all primary BHs, is poorly constrained. From the posterior distributions, we can only exclude small values (< 10M ) of σ 2 . The µ 2 for this component can  be either smaller than µ 1 of the first Gaussian component, or larger than m max . The credible region for Model III is shown in the right panel of Fig4, and due to the large uncertainties on the second Gaussian component, it could alternatively represent a weak and shallow component rising below m max and extending to higher masses rather than a real Gaussian one (which is similar to the case marked with dashed lines in Fig.1). To quantify the preference for the shape of this component by data, we further fix (µ 2 = 83, σ 2 = 17) and (µ 2 = 40, σ 2 = 44) (chosen according to the 68% upper and lower bounds of the posterior distributions for µ 2 and σ 2 ) to reanalyze the data. Posterior distributions for the parameters and the credible regions of the primary mass spectra of these two cases are presented in Fig.3 and Fig.5. The resulting Bayes factor between the two cases is B = 0.17, which indicates a modest support for the component being a shallow one peaks below m max . Nevertheless, it is still lack of strong evidence to exclude the case in which the component has a relatively narrower (σ 2 20) Gaussian shape and peaks after m max . Another possible effect of PPISNe in addition to the formation of mass gap is leaving an excess of BHs near the lower edge of the gap (Talbot & Thrane 2018). If this excess fully accounts for the first Gaussian component in Model IV, the number of black holes in such component should be no more than the number of black holes that would have been formed by the power-law component continued to the upper limit of the mass gap m PI (Talbot & Thrane 2018), i.e., where we take m PI = 150M . To test if this hypothesis is supported by our posterior, we compute the fraction of black-holes, pertaining to the power-law component of the mixture, with a mass above m max . We do this using the posteriors for the parameters of Model IV. We find that only ∼ 13% of the posteriors predict such a fraction larger than λ λ 1 . We further apply the restriction in Eq.(7) to the Bayesian inference of Model IV, and find that the resulting Bayes factor (with respect to Model II) is 0.23. Since the B for Model IV without this restriction is 1.31, the PPISNe origin of the first Gaussian component is less preferred (but not excluded).

CONCLUSION AND DISCUSSION
In this work, we study the primary mass distribution of the merging BBHs and focus on probing the presence of a sudden drop of the mass function at the black hole mass of ≤ 65M , as predicted in pulsational pair instability supernova model. We construct two empirical mass functions, and by performing Bayesian inference, we find that these two models are still comparable with the most preferred empirical models (i.e., Model I and Model II in this work) found in previous studies and a cutoff of the black hole mass function at ∼ 50M is indeed consistent with the data (see Fig.4). The very massive sub-population, which accounts for at most several percents of the total merging black holes, may be from hierarchical mergers or other processes. Note that in this work we concentrate on the BBH systems. In the future a reasonably large neutron star−black hole merger event sample is expected to be available, with which the mass function of these black holes can be reconstructed ) and it would be quite interesting to see whether the black hole mass functions are significantly different among different binary systems.
As already discussed in Sec.4, there is a consistency between the inferred parameters for the first segment in Model III and the ones for Model B in Abbott et al. (2019a). We therefore suspect that Abbott et al. (2019a) has already recovered spectral shapes that roughly match the actual primary mass distribution below the cutoff, while missed the sub-population above the cutoff due to the small number of the GWTC-1 events. Based on our analysis about Model IV, we find moderate support for the low mass Gaussian component being not originated from PPISNe. On the other hand, if PPISNe truly account for this excess of BHs, the first Gaussian component should be relatively weak, and as shown by the right panel of Fig.5, the credible region for Model IV in this case is very similar to the credible region for Model III (the right panel of Fig.4). It is difficult to place good constrains on all of the parameters in Model IV. As shown in Fig.3 and Fig.4, the three superimposed components may cover up the cutoff feature of the power-law component. Nevertheless, since their origins are unclear, the two sub-dominant components in model IV might essentially belong to one population peaking at ∼ 33M and having an extended tail. For the purpose of this work, we do not attempt to recover the exact shape of these components by introducing more assumptions (making more complicated models) with current data.
If the primary mass distribution indeed consists of different sub-populations, evidences may be found elsewhere in addition to the mass spectrum. For example, the inclusion of spin data (although additional considerations are needed to construct the spin model) would enhance the ability for model comparison in the inference. The expected spins for BHs with primary mass above the mass cutoff could be larger and more isotropic if we assume these BBHs are formed by dynamical capture. It is worthy of noting that GW190521 has χ p ∼ 0.6 (Abbott et al. 2020b), which is larger than the χ p of other observed events with lower masses. On the other hand, the evidence for anti-aligned spin and the non-zero χ p of the whole population found in The LIGO Scientific Collaboration et al. (2020) may suggest the presence of different formation channels. However, The LIGO Scientific Collaboration et al. (2020) also pointed out that there is no strong evidence for variation of the spin distribution with mass. More events are needed if we want to take all these possibilities into account to make solid conclusions.