Black Hole Coagulation: Modeling Hierarchical Mergers in Black Hole Populations

The Advanced Laser Interferometer Gravitational Wave Observatory (LIGO; the LIGO Scientific Collaboration 2015) and Virgo (Accadia et al. 2012) detectors have discovered and will continue to discover gravitational waves (GWs) from coalescing binary black holes (BBHs) and neutron stars. So far, several tens of BBH detection candidates have been reported in O3, LIGO's current observing run, and several hundreds more detections are expected over the next five years (The LIGO Scientific Collaboration & the Virgo Collaboration 2016a, 2016b). As the cosmic census these surveys provide grows more comprehensive, these observations will discriminate between formation scenarios of compact-object binaries (Mandel & O'Shaughnessy 2010; Breivik et al. 2016; Nishizawa et al. 2016; Rodriguez et al. 2016). A few formation scenarios invoke "hierarchical" growth of BBHs in which some black holes are themselves products of previous mergers. These hierarchical mergers could occur in globular clusters (Portegies Zwart & McMillan 2002; Gültekin et al. 2006), active galactic nucleus (AGN) disks (see, e.g., McKernan et al. 2012, 2019; Bartos et al. 2017; Yang et al. 2019), or nuclear star clusters (Antonini & Rasio 2016). Alternatively, the hierarchical merger components could have been produced in the early universe due to primordial density fluctuations forming primordial black holes (Clesse & García-Bellido 2015, 2017; Belotsky et al. 2019). Notably, hierarchical growth produces distinctive signatures in the mass and spin distribution (Fishbach et al. 2017; Gerosa & Berti 2017; Kimball et al. 2019; McKernan et al. 2019; Yang et al. 2019), the most generic of which is a population of spinning black holes. For some realizations of these models' parameters, several groups have made predictions about the black hole mass and spin distribution (Rodriguez et al. 2016; Belczynski et al. 2017; McKernan et al. 2019). Additional investigations have assessed whether existing observations are compatible with these models, focusing on the individual event GW 170729 (Chatziioannou et al. 2019; Kimball et al. 2019; Yang et al. 2019).

In this work, we introduce a generic, parameterized framework that accounts for BBHs that form through hierarchical mergers. The method treats black holes as particles in a box that undergo collisions based on an effective cross section. This framework can incorporate a wide range of submodels and prescriptions, enabling one to create models that are purely phenomenological or instead heavily based on detailed astrophysical investigations and simulations. We provide a concrete implementation of our framework, including astrophysically realistic initial conditions. Using existing GW observations, we perform Bayesian inference on our parameterized model.

Our paper is organized as follows. In Section 2, we describe our framework for hierarchical mergers and some parameterizations within the framework, illustrating them with simple examples. We also describe our fiducial initial conditions for BBH populations. In Section 3, we show how to constrain this parameterized model through comparison with GW observations from LIGO and Virgo's first and second observing runs. In Section 4, we discuss the results of our parameter inference on the LIGO-Virgo data, the overall efficacy of our framework, and possible extensions to the parameterizations explored herein. Finally, we summarize the results of our investigation in Section 5.

2.1. General Framework

We employ a flexible method for self-consistently generating mass and spin distributions for BBHs that include a subpopulation of hierarchical mergers. Rather than model the complex dynamics of individual stellar environments, we build a parameterized phenomenological model that describes the aggregate properties of merging binaries in the local universe using volume-averaged coupling coefficients. Our framework incorporates three generic physical processes. First, black holes coagulate when pairs of compact objects merge into single compact objects that may remain in the population. Second, we allow for depletion, where some compact objects leave dense environments and no longer have an opportunity to merge with other objects. Finally, we allow for augmentation, where some process introduces new compact objects to the hierarchical interacting environment (e.g., BHs from stellar collapse or AGN disk dynamics).

Following similar investigations (Lissauer 1993; Christian et al. 2018), we model these effects with a Monte Carlo procedure, designed to approximate a continuous-time coagulation equation (Smoluchowski 1916), which has the qualitative form

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}f(x;t)=\displaystyle \frac{1}{2}\displaystyle \int {dx}^{\prime} {dx}^{\prime\prime} f(x^{\prime} ;t)f(x^{\prime\prime} ;t){\rm{\Gamma }}(x^{\prime} ,x^{\prime\prime} ;t)\\ \times \,\delta ({x}_{\mathrm{rem}}(x^{\prime\prime} ,x^{\prime} )-x)-\displaystyle \int {dx}^{\prime} f(x;t)f(x^{\prime} ;t){\rm{\Gamma }}(x,x^{\prime} ;t)\\ +\,r(x;t)-d(x;t)\end{array}\end{eqnarray} \tag{ 1 }$$

where x denotes black hole parameters, $f(x;t)$ denotes the BH parameter distribution function at time t, ${\rm{\Gamma }}(x,x^{\prime} ;t)$ denotes a volume-averaged interaction rate (i.e., coagulation), and r(x; t) and d(x; t) are the augmentation and depletion rates of black holes with parameters x at time t. The first integral describes the accumulation of black holes with parameters x due to mergers of pairs of black holes with parameters $x^{\prime} ,x^{\prime\prime}$. The delta function enforces that the final parameters x are produced by a merger of BHs with parameters $x^{\prime} ,x^{\prime\prime}$. The function ${x}_{\mathrm{rem}}(x,x^{\prime} )$ computes the remnant parameters from merger component parameters x and $x^{\prime}$.⁸ The second integral accounts for the decrease of black holes with parameters x due to mergers with other black holes with parameters $x^{\prime}$, and its integrand $f(x;t)f(x^{\prime} ;t){\rm{\Gamma }}(x,x^{\prime} ;t)$ is equivalent to the merger rate as a function of parameters. In the absence of augmentation or depletion, the total number of black holes $\int {fdx}$ decreases as $-1/2\int {dxdx}^{\prime} {\rm{\Gamma }}(x,x^{\prime} ;t)f(x;t)f(x^{\prime} ;t)$, as each merger reduces the total number of black holes by one. (The factor of 1/2 is a statistical factor to avoid overcounting).

Given an initial condition f(x, t₀), an interaction rate ${\rm{\Gamma }}(x,x^{\prime} ;t)$, a map between merger components and remnants ${x}_{\mathrm{rem}}(x,x^{\prime} )$, and prescriptions for augmentation and depletion, the solution $f(x;t)$ can in principle be computed. This approach is highly modular and can incorporate complex dynamical physics via the coagulation, augmentation, and depletion functions. Additionally, existing black hole population models can be extended to include hierarchical merger effects in our framework. With this framework in hand, we first describe our method for computing these hierarchical merger distributions and then turn to astrophysically motivated choices for these functions and their application to GW data.

2.2. Monte Carlo Implementation

2.3. Model Prescriptions and Parameterizations

In this section, we describe our inputs to Equation (1), which we have chosen to be simple, computationally efficient, and astrophysically motivated. Notably, the choices we make here all assume an isotropic interaction environment, with randomly oriented spins, which may not be well suited to some environments such as AGN disks. However, we emphasize that alternative effects can be readily incorporated into this framework if desired. To limit the scope of our investigations, augmentation is not considered in this work, but future studies could include it.

2.3.1. Coagulation

For simplicity, we assume the volume- and time-averaged interaction rate Γ depends only on binary masses $(m,m^{\prime} )$, with a parametric form

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{m,m^{\prime} }\propto {\left(\displaystyle \frac{(m+m^{\prime} )}{{M}_{\mathrm{ref}}}\right)}^{a}{\left(\displaystyle \frac{\eta }{{\eta }_{\mathrm{ref}}}\right)}^{b}.\end{eqnarray} \tag{ 2 }$$

(This single interaction term is designed to capture the average effect of interactions throughout the volume, on the long timescales over which the BH mass distribution evolves appreciably through hierarchical mergers.) We include the total binary mass dependence ${\left(\tfrac{(m+m^{\prime} )}{{M}_{\mathrm{ref}}}\right)}^{a}$ for two reasons. First, bigger black holes have larger "cross-sectional areas" with which they can interact with other objects. In the limiting case of spheres in a gas with radii r, one would expect ${\rm{\Gamma }}{(r,r^{\prime} )\propto (r+r^{\prime} )}^{2}$. To account for the complex dynamics of interacting black holes, we do not fix the power to 2 and instead let it vary, and since a black hole's Schwarzchild radius is directly proportional to its mass, we replace radii with masses. The second effect this term accounts for is dynamical friction, which brings more massive black holes to dense centers of clusters where they can merge. The second term ${\left(\tfrac{\eta }{{\eta }_{\mathrm{ref}}}\right)}^{b}$ depends on the symmetric mass ratio η to account for a possible preference for mergers to choose more equal or unequal masses (see, e.g., Fishbach & Holz 2019). In globular clusters, for example, mass segregation may favor equal-mass mergers over those of unequal mass (Sollima 2008; Park et al. 2017; Rodriguez et al. 2019).

Although we assume that the black hole spins do not influence the interaction rate, we do keep track of the spin magnitudes of the black holes and calculate final black hole spins from initial component parameters. We use fits to numerical relativity simulations from Tichy & Marronetti (2008) for ${x}_{\mathrm{rem}}(x,x^{\prime} )$, the final mass and spin of a remnant black hole given the masses and spins of the individual components. To further simplify our calculations, we assume the hierarchical environment is isotropic, so only spin magnitudes χ need to be tracked since spin orientations are random. As such, we can simply write x = (m, χ) in this prescription.

2.3.2. Depletion

Remnant black holes experience recoil kicks that may eject the remnant from the environment and prevent it from merging again with another object. Here we consider two cases: 1. no depletion and 2. cluster depletion. In the first case, we assume no black holes leave the environment; in the second, we use the "V₄₅₉" fits to numerical relativity simulations from Zlochower & Lousto (2015) for recoil velocities with a prescription for the distribution of cluster escape velocities to calculate the depletion rate. We parameterize the depletion based on the magnitude of the recoil velocity v_kick and ignore the recoil direction, although future studies could incorporate the recoil direction to account for anisotropy in the merger environment. For cluster depletion, we assume that black holes are in star clusters with a variety of density profiles and hence a variety of central escape velocities. We write the escape probability as:

$$\begin{eqnarray}&&\begin{array}{l}p(\mathrm{escape}| {v}_{\mathrm{kick}},{\mu }_{M},{\sigma }_{M},{\mu }_{{r}_{0}},{\sigma }_{{r}_{0}})\\ \propto \displaystyle \int \displaystyle \int d\mathrm{log}{Md}\mathrm{log}{r}_{0}{\rm{\Theta }}\left[\displaystyle \frac{1}{2}{v}_{\mathrm{kick}}^{2}-\displaystyle \frac{{GM}}{{r}_{0}}\right]\\ \times \exp \left(-{\left(\displaystyle \frac{\mathrm{log}M-\mathrm{log}{\mu }_{M}}{{\sigma }_{M}}\right)}^{2}-{\left(\displaystyle \frac{\mathrm{log}{r}_{0}-\mathrm{log}{\mu }_{{r}_{0}}}{{\sigma }_{{r}_{0}}}\right)}^{2}\right).\end{array}\end{eqnarray} \tag{ 3 }$$

The Heavyside function enforces that remnants with kick velocities larger than the cluster escape velocity are ejected. The cluster escape potential is given by a Plummer model and the black holes are always assumed to be at the centers of clusters. The last line of terms describes the distribution of cluster masses M and effective radii r₀ in the Plummer model. We take these cluster masses and radii to be log-normally distributed and parameterized by μ_M, σ_M, ${\mu }_{{r}_{0}}$, and ${\sigma }_{{r}_{0}}$, but emphasize that other choices could be made for all of these depletion prescriptions.

2.3.3. Natal Populations

The final ingredient we need to specify in our model is the initial distribution of masses and spins $f(x;{t}_{0})$. We hereafter refer to this as the "natal" distribution, and take it to be the distribution of black hole parameters formed at black hole birth. A variety of choices could be made, but here we restrict ourselves to two cases. The first case is a simple power-law mass function in component masses with lower and upper mass cutoffs

$$\begin{eqnarray}p(m)\propto \left\{\begin{array}{ll}{m}^{-\alpha }, & \mathrm{if}\ {m}_{\min }\leqslant m\leqslant {m}_{\max }\\ 0, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

In all cases we use the fiducial value m_min = 5 M_⊙ for the sake of simplicity. For the other parameters, we either fix them to fiducial values of α = 2.35 (from a Salpeter initial mass function (IMF)) and ${m}_{\max }=20\,{M}_{\odot }$ (from early stellar evolution modeling), or we allow the data to tune them, assuming uniform priors. The blue curves in Figure 1 show two examples of black hole natal mass distributions with the Salpeter prescription. The fiducial Salpeter mass distribution for black holes is a basic model that assumes that the fraction of mass retained from stellar birth to black hole formation, m/m_ZAMS, is constant across all masses. This is unlikely to be true in reality, as the processes undergone by a star depend strongly on its mass. To take things a step further in complexity, we still assume the mass distribution follows a power law, but with an index α that differs from the IMF's value. This is still fairly unrealistic, as the black hole natal mass spectrum is not expected to be this simple, but this at least lets the data determine the general trend of the spectrum.

Zoom In Zoom Out Reset image size

Our second model has a better footing in physical principles, but loses some flexibility. We assume a pure Salpeter IMF for the zero-age main-sequence (ZAMS) masses in the range $[5,\infty )\,{M}_{\odot }$ and evolve them to black holes using the Fryer et al. (2012) Rapid model. (Our calculations implicitly adopt the same wind mass-loss model as employed in that study.) This introduces an additional hidden variable, the stellar metallicity Z_metal for each progenitor star. The red and green curves in the bottom panel of Figure 1 show our inferred progenitor distributions for two choices of Z. In principle, this should be a random variable, obeying some distribution that may correlate with the IMF. For simplicity, however, and motivated by the approximate similarity between these two distributions, we fix this to a constant ${Z}_{\mathrm{metal}}^{* }$, assumed to be the same for every progenitor.

Now we turn to to the black hole natal spins. Black hole natal spins remain a matter of considerable observational and theoretical debate. Motivated by LIGO's observations and recent modeling (Belczynski et al. 2017; Farr et al. 2018; The LIGO Scientific Collaboration & the Virgo Collaboration 2018; Fuller & Ma 2019), we adopt a simple fiducial choice: all BHs in our original population have small characteristic spin magnitudes, drawn from a Beta distribution with mean(χ) = 0.047 and Var(χ) = 0.002. We also assume that the spin directions in the natal population are randomly oriented, but again we emphasize that other choices could be made.

2.3.4. Merger Rates

As described in Section 2.2, our Monte Carlo procedure works with a finite set of black holes. We take these black holes to be a proxy for the entire population and assume that the overall merger rate of black holes is simply a scaled population of those generated in our Monte Carlo simulations. We also stipulate that the merger rate density is constant in comoving volume. Future studies could certainly incorporate more detailed effects, but here we opt for simplicity. In the following section, we show normalized distributions of the masses and spins of black holes, but in Section 3 we present inference results that allow the merger rate density to be inferred by the data.

2.4. Characterizing the Parameters

To elucidate the effect of each of the parameters described in the previous section, we take the reader through a sequence of examples. The examples we present here are primarily for illustration and do not necessarily represent parameters that describe the observed population of black hole mergers to date. Note that the histograms and kernel density estimate curves shown here are not explicitly used in our analysis; they are simply representations of the samples from our Monte Carlo procedure.

2.4.1. Time Evolution

As hierarchical mergers occur, a secondary population of high-mass, high-spin black holes begins to form alongside the natal population. In our Monte Carlo procedure, the time evolution of the population is reduced to individual time steps, as described in Section 2.2. Figure 2 illustrates how the population changes with each time step. Starting with a Salpeter IMF with ${m}_{\max }=20{M}_{\odot }$ and Beta-distribution spin magnitudes as the natal population (which we take as our fiducial natal population) we evolve the population forward for three iterations, allowing 5% of the black holes to merge at each step and setting the coupling strength to a = 2 and b = 0. The red, blue, and black lines show the distributions of the total masses of mergers for time steps 0, 1, and 2, respectively.

Zoom In Zoom Out Reset image size

Since the remnant black holes inherit angular momentum from their parents and from their orbit, hierarchical mergers also produce a strong evolution of BH spins (Fishbach et al. 2017; Gerosa & Berti 2017). With successive mergers, the total mass distribution tends toward higher masses, and an island of high-mass, high-spin black holes begins to grow. Figure 3 shows 90% and 99.9% confidence intervals for the joint mass–spin distributions at T = [0, 1, 2]. Notably, hierarchical mergers of comparable-mass binaries introduce a characteristic peak near $\chi \simeq 0.7$, which is why the top panel of Figure 3 shows a surplus of black holes near that spin magnitude. Generically, hierarchical mergers should produce a similar subpopulation of high-mass, high-spin black holes, since general relativity predicts that a postmerger remnant black hole is always more massive than either of its premerger components and its final spin is away from zero. The χ_eff versus chirp mass distribution in the lower panel shows that while χ₁ tends to be large for the hierarchically produced mergers, the χ_eff distribution is smoothed out around 0, since the black hole spin directions are isotropically distributed.

Zoom In Zoom Out Reset image size

2.4.2. Coupling Strength

The overall mass and spin distributions are sensitive to the average coupling strength of black holes. Figure 4 shows the total mass distribution (top panel) and the primary mass distribution (bottom panel) after four time steps for different values of a and b. Increasing the total mass coupling parameter a drives the most massive mergers to occur, causing the total mass distribution to quickly expand to higher masses, while increasing the symmetric mass-ratio coupling b simply forces most mergers to be of equal-mass components. Cranking up a and b simultaneously gives particularly interesting behavior. In those cases, the heaviest black holes take place in mergers, and the products of those mergers are likely to merge again, which can create multiple distinct peaks in the mass distributions. As a result, in the Salpeter natal distribution example, our procedure produces a characteristic "smoothed staircase" mass distribution, with "steps" in the mass distribution appearing at multiples of the primordial maximum mass m_max,0. At very high mass, these "step" features become smoothed out.

Zoom In Zoom Out Reset image size

The mass ratio and spin distributions also have characteristic features. When a is large but b is small, a population of highly unequal mass mergers can be produced, as seen in the purple curves of Figure 5. A near-flat mass ratio distribution (shown in green) is found for a = 2 and b = 0 in this case, because the natal mass distribution power-law slope (α = 2.35) is nearly matched to the total mass coupling, so the dearth of higher mass black holes is exactly counteracted by their higher likelihood of participating in mergers. As b is increased, the distribution begins to favor equal-mass mergers, as shown in the black curve.

Zoom In Zoom Out Reset image size

Figure 6 shows contours of the joint primary mass and χ₁ distribution for different coupling strengths after four time steps. While a high-mass, high-spin subpopulation is present in all the cases considered here, they are notably affected by the coupling strength parameters. When b is large, the subpopulation is more concentrated at χ₁ ∼ 0.7, because the mergers tend to be equal mass and therefore have a similar final remnant spin.

Zoom In Zoom Out Reset image size

2.4.3. Depletion

The most widely proposed hierarchical scenario involves hierarchical formation in globular clusters. Merging black holes will be very frequently ejected from these low binding energy environments, strongly suppressing the prospects for hierarchical mergers through multiple generations (Favata et al. 2004; Merritt et al. 2004; Gerosa & Berti 2019; Rodriguez et al. 2019). To illustrate how depletion impacts the observed merger distributions, we incorporate the cluster depletion model from Section 2.3.2 into a hierarchical merger population. Figure 7 plots three total mass distributions, one without depletion effects, one with "light" clusters (${\mu }_{M}=5\times {10}^{4}{M}_{\odot }$, ${\mu }_{{r}_{0}}=10$ pc, ${\sigma }_{M}={\sigma }_{{r}_{0}}=1$), and one with "heavy" clusters (${\mu }_{M}=5\times {10}^{5}{M}_{\odot }$, ${\mu }_{{r}_{0}}=5$ pc, ${\sigma }_{M}={\sigma }_{{r}_{0}}=1$). These hierarchical distributions are evolved forward four time steps from a fiducial natal distribution under these three depletion prescriptions and with a = 2, b = 0. This figure shows that as the confining potentials become shallower, remnant black holes are kicked from the environment so that hierarchical mergers are strongly suppressed, as known from previous work.

Zoom In Zoom Out Reset image size

2.4.4. Natal Distributions

As we have seen in the previous examples, the hierarchical distributions produced in our framework contain imprints of the natal populations. Figure 8 plots three hierarchical merger total mass distributions after three time steps assuming the strong coupling parameters a = 4, b = 20. Unsurprisingly, the natal distributions with support at higher masses quickly evolve to have high-mass mergers. Additionally, the more complex structure in the Fryer natal mass distributions is imprinted in the evolved hierarchical distributions, while the Salpeter-based mass distributions are more smoothed out. In sum, the natal distribution is crucially important to the evolution of the mass distribution when hierarchical mergers can take place.

Zoom In Zoom Out Reset image size

We use an updated version of the PopModels population inference code (Wysocki et al. 2019) to compare our hierarchical formation model to real GW observations from GWTC-1 (the LIGO Scientific Collaboration & the Virgo Collaboration 2019). For each collection of observations ${ \mathcal D }$, this code evaluates the inhomogeneous Poisson likelihood

$$\begin{eqnarray}&&{ \mathcal L }({ \mathcal R },{\rm{\Lambda }})\propto {e}^{-\mu ({ \mathcal R },{\rm{\Lambda }})}\prod _{n=1}^{N}\int d\lambda \,{{\ell }}_{n}(\lambda )\,{ \mathcal R }\,p(\lambda | {\rm{\Lambda }}),\end{eqnarray} \tag{ 5 }$$

where ${{\ell }}_{n}(\lambda )=p({d}_{n}| \lambda )$ is the likelihood of data d_n given binary parameters λ, $\mu ({ \mathcal R },{\rm{\Lambda }})$ is the expected number of detections, ${ \mathcal R }$ is the merger rate, and Λ refers to any relevant model parameters: all parameters needed to characterize our hierarchical evolution equations, along with the choice of metallicity and initial conditions. Unlike Wysocki et al. (2019), we evaluate the integrals $\int d\lambda {{\ell }}_{n}(\lambda )p(\lambda | {\rm{\Lambda }})$ by using Monte Carlo integration via samples drawn from our hierarchical model $p(\lambda | {\rm{\Lambda }})$, combined with an analytic likelihood ${{\ell }}_{n}(\lambda )$.

We perform inference with four models, which are described in Table 1. The right-hand side describes prescriptions for fixed values and priors in each model. The posterior distributions on our inference parameters are shown in Figures 9–11. Model 1 is our most basic phenomenological model, the natal distribution having a power law in component mass and zero spin. The number of iterations and mass coupling parameters are inferred from the data. The blue curves in the right panel of Figure 9 show that the data have a slight preference for $a\sim 2$ and a strong preference for large b values around b ∼ 30. The overall rate density of mergers in our hierarchical model ${{ \mathcal R }}_{{\rm{h}}}$, shown in the left panel, is consistent with the rates inferred in the LIGO Scientific Collaboration & the Virgo Collaboration (2018). Additionally, the natal distribution power-law index and maximum mass are constrained to similar values to those found in the LIGO Scientific Collaboration & the Virgo Collaboration (2018), as seen in the blue curves of Figure 10. Given that our hierarchical model reduces to a nonhierarchical model in the low–time step limit and the data favors fewer time steps, it is not surprising that our natal population parameters match the overall population parameters in the LIGO Scientific Collaboration & the Virgo Collaboration (2018). The inference on the number of time steps is shown in Figure 9 in terms of the variable N_gen = T + 1, which is the highest allowed generation of black holes in the population.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Hierarchical Merger Models Fit to O1/O2 Data

Model	Description
Model 1	Natal population: power law in component mass
	${m}_{\min }=5{M}_{\odot }$
	$\alpha \in [-3,5]$, uniform
	${m}_{\max }\in [15,50]$, uniform
	${\mathbb{E}}[\chi ]=0.047$
	$\mathrm{Var}[\chi ]=0.002$
	Coagulation parameters:
	a ∈ [1, 6], uniform
	b ∈ [1, 100], log uniform
	T ∈ [0, 9], uniform
	w = 0.05
Model 2	Same as Model 1, except T ∈ [0, 5] and
	Beta-distribution natal spins
	${\mathbb{E}}[\chi ]\in [0,1]$, uniform
	Var[χ] ∈ [0.25], uniform
	Depletion
	${\mu }_{M}=5\times {10}^{5}{M}_{\odot }$
	${\mu }_{{r}_{0}}=5$ pc
	${\sigma }_{M}={\sigma }_{{r}_{0}}=1$
Model 3	Same as Model 2 except
	${\mu }_{M}\in [{10}^{5},{10}^{10}]{M}_{\odot }$, uniform
	${\mu }_{{r}_{0}}\in [5,{5}^{5}]$ pc, uniform
Model 4	Mixture of Model 2 and Model A of the LIGO Scientific Collaboration & the Virgo Collaboration (2018)
Model 5	Same as Model 1 except
	Fryer rapid SN natal population
	Mixture with Model A of the LIGO Scientific Collaboration & the Virgo Collaboration (2018)
	T = 2

Download table as:  ASCIITypeset image

The most widely proposed hierarchical scenario, however, involves hierarchical formation in globular clusters. Merging black holes will be very frequently ejected from these low binding energy environments, strongly suppressing the prospects for hierarchical merger (Gerosa & Berti 2019; Rodriguez et al. 2019). Model 2 adds a depletion prescription to Model 1 with fixed cluster mass and radius distribution parameters. In this case, similar coupling and natal distribution parameters to Model 1 are inferred, which is shown in orange in Figures 9 and 10. Notably there is a slight preference for higher total mass couplings a for Model 2 compared to Model 1, because the depletion effects strongly suppress hierarchical mergers, and therefore higher masses from the natal population are favored. The strong depletion in this case also results in no preference on the number of time steps.

We then allow the cluster mass and radius distribution parameters to vary in Model 3. The results of inferring cluster sizes are shown in Figure 11. Interestingly, the cluster radii and masses are pushed to large values, far greater than those of real star clusters. This is partially an artifact of the parameterization chosen here. The gravitational potential in the Plummer profile is sensitive only to the ratio of cluster mass to cluster radius, so if we consider the ratios of μ_M to ${\mu }_{{r}_{0}}$, the inferred values are roughly similar in gravitational potential to the fixed values used in Model 2. In other words, the data prefer somewhat shallow potentials wherein hierarchical mergers are suppressed. A future parameterization may instead opt for a distribution of gravitational potentials rather than cluster parameters.

Next we consider two mixture models. In the first (Model 4), we fit a mixture of our Model 2 with Model A from the LIGO Scientific Collaboration & the Virgo Collaboration (2018). Then in Model 5 we create a mixture of Model A and our hierarchical model applied to the Fryer rapid SN natal population with no depletion and exactly 3 time steps of evolution. In these mixture analyses, we simultaneously fit the parameters of Model A (power-law index, maximum mass cutoff, overall rate) and the parameters of the hierarchical model. Figure 10 shows the distributions of population parameters for the "field" (Model A) and "hierarchical" mixtures. The mixtures complicate the picture significantly. Model 4 (shown in red) in particular has little discerning power on its underlying population parameters due to the additional model freedom. Model 5 (purple) on the other hand, for which the natal distribution and number of time steps are fixed, shows some interesting behavior. In particular, the "field" (i.e., Model A) component parameters are driven to a near-flat distribution in component masses with a slightly lower mass cutoff than for the other models' considered inferences. Meanwhile, the coagulation parameters a and b are pushed to lower values. These shifts in the inferred parameters are likely due to fixing the number of time steps to 3 with no depletion. Fixing the number of time steps to 3 favors the existence of some hierarchical mergers that tend to be higher mass. To counteract the buildup of too many high-mass black holes compared with the data, the mass distribution of the field population is cut off at a lower m_max and the coagulation mass coupling is decreased. Also, the inferred metallicity Z_m of the natal population also slightly favors higher values, which pushes the natal mass distribution to lower masses, alleviating some of the unwarranted buildup of high-mass black holes. Lastly, the contribution of the hierarchical population is subdominant to the field population, as seen in Figure 11.

A hierarchical formation scenario provides an efficient way to produce binaries that would otherwise be challenging to generate: high masses, exceptional mass ratios, and characteristically high spins. The identification of binaries with characteristically extreme properties could provide a clear indication of hierarchical formation. In this section we explore our posterior predictive constraints on these scenarios, within the framework of the constrained fiducial model described above. We also discuss further extensions of the models presented here and the overall effectiveness of this framework.

4.1. Posterior Predictive Distributions

Figure 12 shows our inferred posterior mass distributions, both intrinsic and detection-weighted, which resemble the conclusions in the LIGO Scientific Collaboration & the Virgo Collaboration (2018). Specifically, we infer a mass distribution for the more massive component in merging black holes (m₁) that is approximately a power law between 10 M_⊙ and 30 M_⊙, followed by a rapid decrease at higher mass. Notably, this figure shows characteristic decay and "echo" features at about 30 M_⊙, inherited by our formation model; these features could be probed by future observations and used to better constrain hierarchical formation.

Zoom In Zoom Out Reset image size

Figure 13 shows our inferred mass ratio distribution. Because GWTC-1 does not include a significant component of asymmetric binaries, our posterior necessarily strongly favors hierarchical models that preferentially produce binaries with q ≃ 1. Constraints on binary mass ratios will very strongly constrain prospects for hierarchical formation, particularly insofar as some hierarchical scenarios produce significant numbers of highly asymmetric mergers (Yang et al. 2019).

Zoom In Zoom Out Reset image size

4.2. Possible Extensions

4.3. Efficacy of this Hierarchical Merger Population Framework

Observing hierarchically formed black holes is an exciting prospect for GW detectors. We have presented here a self-consistent framework for generating black hole merger populations that includes hierarchical formation. This framework evolves arbitrary natal black hole populations, enabling any existing black hole distributions to be extended to include hierarchical mergers. With the cases we explore here, our fits suggest that scenarios with many hierarchical mergers are disfavored.

In this work, we simulate coagulation and depletion effects while assuming the BBH population is not continuously repopulated from another reservoir of black holes. We will explore self-consistent repopulation in later work. We also perform simplified averaging, not allowing for a distribution of initial conditions like metallicity or for trends versus redshift. Our scheme ignores higher order correlations and multibody effects, thus averaging everything into an effective cross section that is constant. Our approximation is reasonable in the limit of weak hierarchical reprocessing dominated by a low-mass seed population; we defer more sophisticated averaging to future work. Additionally, we adopt a simple dependence of Γ on total mass, allowing it to increase without bound according to a single power law as the binary mass increases. More detailed investigations will produce more complex dependence of Γ on mass. The results of our inference on GWTC-1 with partially constrained versions of our model framework show consistency with the LIGO Scientific Collaboration & the Virgo Collaboration (2018), but more detections are required to make more definitive statements about whether hierarchical formation of black holes is at work. Our fits to GWTC-1 hint that hierarchical merger scenarios are not required to fit the population, but are also not ruled out by the population. However, the conclusions drawn herein are subject to the simplifying assumptions made for this preliminary analysis. Incorporation of more astrophysically motivated inputs to the framework will be necessary to put the tightest constraints on the formation environments and scenarios. With the wealth of black hole merger detections we expect to see in the coming years, the prospects for uncovering hierarchically formed black holes are promising, and the framework we have presented herein is well suited for such investigations.

The general purpose hierarchical population code will be released for public use in the near future.

Z.D., R.O.S., and D.E.H. initiated this work at the Kavli Institute for Theoretical Physics, supported in part by NSF PHY-1748958. Z.D. was supported by the NSF Graduate Research Fellowship Program, grant DGE-1144082. R.O.S. and D.W. are supported by NSF AST-1707965 and PHY-1909534. Z.D. and D.E.H. are supported by NSF grant PHY-1708081 and the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation. D.E.H. also gratefully acknowledges support from the Marion and Stuart Rice Award. B.F. is supported by NSF grant PHY-1807046. The authors thank the LIGO Scientific Collaboration for access to computing resources and public data and gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, and the Max Planck Society (MPS) for support of the construction of Advanced LIGO. Additional support for Advanced LIGO was provided by the Australian Research Council. For their use in implementing our hierarchical inference code, we would like to acknowledge Numpy and Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007), and h5py (Collete 2013).

In the text, we consistently employ a concrete Monte Carlo implementation of hierarchical mergers. This powerful method allows us to efficiently incorporate the best merger physics, but uses discrete generations. In this appendix, for pedagogical purposes we provide a toy model implementation of a true continuous-time coagulation equation. In this approach, we consider only the evolution of binary mass, assuming no mass is lost during mergers; we neglect spins and depletion. After these simplifications, our model is essentially analytically tractable, and can be understood by both perturbation theory and direct numerical simulation. In this appendix, we present a few supplementary illustrations of these hierarchical calculations to further illuminate our model's behavior at very high mass and in the absence of depletion.

As our first example, to illustrate the parameters ζ and a, Figure A1 shows the results of evolving Equation (1) starting with an initial power-law mass distribution through different ranges of interaction parameter ζ ≪ 1, for two choices of a and for b = 0. The dotted curves show the results of a direct numerical time integration; the solid curves show our Monte Carlo procedure; and different colors indicate different choices for x and a, respectively.

Zoom In Zoom Out Reset image size

This example first shows how the parameter ζ controls the effective number of generations at the reference parameters, absorbing factors present in the overall interaction time T and in the interaction cross section Γ. As expected based on perturbative arguments, higher-order generations increase in significance in proportion to x^g for g the number of generations. At very high mass, the hierarchical mass distribution approaches an exponentially decaying function of m, which increases exponentially with x².¹⁰

Second, this example shows how the coagulation equation successively reprocesses each generation, potentially with different interaction scales. In this example and generally in the usual case that a, b > 0, binaries with smaller masses or more asymmetric mass ratios by construction interact even less frequently. Conversely, within our framework high-mass binaries rapidly undergo multiple generations of mergers. As a result, in this power-law example, our procedure produces a characteristic "smoothed staircase" mass distribution, with "steps" in the mass distribution appearing at multiples of the primordial maximum mass ${m}_{\max ,0}$. At very high mass, these "step" features become smoothed out.

Third, this example shows the importance of the interaction cross section: because we adopt b = 0 (no preference to any mass ratio) and because low-mass BHs are dramatically more prevalent than high-mass binaries, the overall merger rate $f(x)f(x^{\prime} ){{\rm{\Gamma }}}_{x,x^{\prime} }$ for binaries with one massive component x is overwhelmingly dominated by mergers where $x^{\prime}$ is drawn from this scenario's ubiquitous low-mass black holes. As a result of these frequent minor mergers, features in the mass spectrum proportional to the primordial maximum mass are rapidly smoothed out, both in x and as we go to higher multiples of the maximum mass, except for the first feature.

In sum, our coalgulation model naturally "builds up" self-consistent hierarchical populations, producing mass distributions that can (but need not) possess clear features reflecting the number of generations and any sharp cutoffs present in the seed distribution.

As our second example, we consider how features of the high-mass mass distribution are inherited from the low-mass mass spectrum, using a broad, featureless power-law distribution initially $f(m)\propto {m}^{-\alpha }$ at low mass. For simplicity and unlike in the example used above, we consider interactions with $b\gg 1$, insuring that almost all mergers occur between comparable-mass binaries.

Qualitatively speaking, coagulation requires the formation rate of BHs with mass $2m$ must be ${\rm{\Gamma }}f{(m)}^{2}\propto {m}^{a-2\alpha }$, a slope that can be shallower or steeper than the low-mass slope (− α) depending on the sign of a − α. Evidently, as corroborated by Figure A2, larger a favors higher mass black holes and a more extended tail in the mass distribution. As in the previous example, at very high masses the distribution decays exponentially, with a coefficient that depends on a.

Zoom In Zoom Out Reset image size

Rather than use a continuous-time coagulation equation, which implicitly allows BH remnants to participate in subsequent hierarchical mergers immediately, we employ discrete time steps with a specific fraction w of BHs that participate in mergers at each iteration. As $w\to 0$, our algorithm converges to continuous coagulation, because the population changes slowly over many iterations, ensuring that ${\rm{\Delta }}f(x,t)/{\rm{\Delta }}t$ is small at each step and that postmerger remnant black holes are immediately available to merge again. As w increases, our iterative process increasingly differs from continuous evolution. In effect, w encodes a "recycling delay time," i.e., the time for a postmerger remnant black hole to be reintegrated into the population. We emphasize that while larger w loses fidelity to the continuous coagulation equation, high w can still model a real compact-object population. For example, it is conceivable that all natal black holes merged at an early time (i.e., w = 1), and then all participate in second-generation mergers at later times.

As a concrete example, we apply our algorithm to our fiducial power-law natal population using different w values while holding constant the total number of mergers. In other words, we ensure wT is a constant, where T is the number of iterations of our Monte Carlo method. The coupling constants a and b are set to 0. Figure B1 shows the total mass distribution of mergers after $T=2,4,8,16$ iterations with $w=0.08,0.04,0.02,0.01$, respectively. At lower masses, the distributions roughly agree, but only the small w cases have tails extending to higher masses, since those cases are closer to continuous coagulation wherein there is no recycling delay time and postmerger remnants are free to remerge immediately.

Zoom In Zoom Out Reset image size