The Imprint of Superradiance on Hierarchical Black Hole Mergers

The superradiant condition is simply that the boson's angular frequency, ω = μ/ℏ, must satisfy (Bekenstein 1973; Brito et al. 2020a)

$$\begin{eqnarray}&&\displaystyle \frac{\omega }{m}\lt {{\rm{\Omega }}}_{\mathrm{BH}}\equiv \displaystyle \frac{1}{2}\displaystyle \frac{{c}^{3}}{{GM}}\displaystyle \frac{\chi }{1+\sqrt{1-{\chi }^{2}}},\end{eqnarray} \tag{ 1 }$$

where m is the magnetic quantum number corresponding to a specific boson cloud mode and Ω_BH is the angular frequency of the black hole's outer horizon (Teukolsky 2015). Furthermore, the black hole's angular frequency is a function of the its mass, M, dimensionless spin, χ, and dimensionless outer radius, ${\bar{r}}_{+}\equiv {r}_{+}/{r}_{\mathrm{BH}}=1+\sqrt{1-{\chi }^{2}}$, where r_BH = GM/c² and where the characteristic black hole length scale is half the Schwarzschild radius. The energy eigenstates of the boson cloud take a similar form to those of a hydrogen atom, and are denoted by radial n, orbital l, and magnetic quantum numbers (Ternov et al. 1978; Detweiler 1980).

To highlight this comparison, we define the effective fine-structure constant for the "black hole atom":

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{{GM}\mu }{{\hslash }{c}^{3}}.\end{eqnarray} \tag{ 2 }$$

Large values of α lead to significant growth of the boson cloud. If the inequality in Equation (1) is satisfied, the boson cloud extracts rotational energy from the black hole, leading it to spin down. Unstable growth of the boson cloud ceases when the angular frequency of the bosonic wave equals mΩ_BH.

In order to incorporate the effect of superradiance, we compute the final mass and spin of a black hole according to Equation (1) (Brito et al. 2017a; Isi et al. 2019):

$$\begin{eqnarray}&&{M}_{f}\approx {M}_{i}\left(1-\displaystyle \frac{{\alpha }_{i}{\chi }_{i}}{m}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\chi }_{f}=\displaystyle \frac{4{\alpha }_{f}m}{4{\alpha }_{f}^{2}+{m}^{2}}.\end{eqnarray} \tag{ 4 }$$

Here the subscripts i and f refer to the initial and final states of the boson cloud–black hole system, respectively, if the boson cloud is given enough time to saturate a given mode. However, the cloud saturates fully only if $\bar{\lambda }\sim {r}_{\mathrm{BH}}$. Otherwise the cloud does not grow appreciably during the lifetime of the black hole.

The occupation number of a given energy state grows exponentially at the rate (Ternov et al. 1978; Detweiler 1980; Baryakhtar et al. 2017)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jlmn}}\approx 2{\alpha }^{2j+2l+5}{\bar{r}}_{+}(m{{\rm{\Omega }}}_{\mathrm{BH}}-\omega ){C}_{{jlmn}},\end{eqnarray} \tag{ 5 }$$

where C_jlmn is a dimensionless factor dependent on the quantum state inhabited by the ultralight bosons, and j is the total angular momentum quantum number. In the case of scalar ultralight bosons, j = l, and (Ternov et al. 1978; Detweiler 1980)

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{jlmn}}^{(\mathrm{scalar})} & = & \displaystyle \frac{{2}^{4l+2}(2l+n+1)!}{{\left(l+n+1\right)}^{2l+4}n!}{\left(\displaystyle \frac{l!}{(2l)!(2l+1)!}\right)}^{2}\\ & & \times \ \prod _{k=1}^{l}\left({k}^{2}(1-{\chi }^{2})+\displaystyle \frac{4{r}_{+}^{2}}{{c}^{2}}{\left(m{{\rm{\Omega }}}_{\mathrm{BH}}-\omega \right)}^{2}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

The growth rate of the occupation number can be used to approximate the timescale for the size of the cloud to grow by one e-fold, often known as the instability timescale.

Typically, the size of the boson cloud saturates at ∼180 e-folds of the instability timescale (see Equation (14) from Baryakhtar et al. 2017), which we use to compute an approximate growth timescale around a black hole (Arvanitaki et al. 2010; Arvanitaki & Dubovsky 2011; Baryakhtar et al. 2017; Ng et al. 2021a):

$$\begin{eqnarray}&&{\tau }_{{jlmn}}^{(\mathrm{scalar})}\approx \displaystyle \frac{180}{{{\rm{\Gamma }}}_{{jlmn}}}.\end{eqnarray} \tag{ 7 }$$

The growth timescale scales as

$$\begin{eqnarray}&&{\tau }_{{jlmn}}^{(\mathrm{scalar})}\propto M{\alpha }^{-4{\ell }-5},\end{eqnarray} \tag{ 8 }$$

where it is clear that when α ≪ 0.1, the instability growth rate is greatly reduced (Brito et al. 2015).

Finally, due to the greatly differing timescales between states, the mode which determines the final spin of the black hole is the mode with the lowest final spin that grows within the age of the black hole, or evolution timescale, τ_evol (e.g., see Figure 1 in Ng et al. 2021a and Figure 3 in Baryakhtar et al. 2017). Motivated by the strong fine-structure constant dependence, and the associated strong mass dependence, the quasi-adiabatic differential equations governing the evolution of the black hole–boson cloud system can be largely ignored. We consider a black hole to have undergone superradiance only if its age exceeds the growth timescale of the fastest growing mode. For the remainder of this manuscript, we focus on the first three l = m = 1, 2, 3, n = 0, states.

An important corollary of Equations (4)–(7) is that there are regions of the mass–spin parameter space, also known as the Regge plane, where a black hole cannot exist if an ultralight boson is present (Arvanitaki & Dubovsky 2011; Brito et al. 2015, 2017a). These exclusion regions are strongly dependent of the boson's mass and the evolution timescale of the black hole population. In Figure 1, we present the exclusion regions as governed by the first three l = m = 1, 2, 3, n = 0, energy eigenstates over the black hole mass range relevant for hierarchical growth from stellar-mass black holes. Each m corresponds to a particular "bump" in the excluded region, with the m = 1 mode contributing at lower black hole masses. Furthermore, increasing the evolution timescale of the systems results in the growth of boson clouds around lower mass black holes.

Zoom In Zoom Out Reset image size

To model the global properties of a cluster, we assume the heating rate from black hole binaries in the core of the cluster balances the energy flow into the whole cluster, known as Hénon's principle (Hénon 1961, 1975; Gieles et al. 2011; Breen & Heggie 2013; Kremer et al. 2019). The half-mass–radius, heating rate, and escape velocity of the cluster evolve as

$$\begin{eqnarray}&&{r}_{{\rm{h}}}(t)={r}_{{\rm{h}},0}{\left(\displaystyle \frac{3}{2}\displaystyle \frac{\zeta (t-{t}_{0})}{{\tau }_{\mathrm{rh},0}}+1\right)}^{2/3},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\dot{E}(t)={\dot{E}}_{0}{\left(\displaystyle \frac{3}{2}\displaystyle \frac{\zeta (t-{t}_{0})}{{\tau }_{\mathrm{rh},0}}+1\right)}^{-5/3},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}(t)={v}_{\mathrm{esc},0}{\left(\displaystyle \frac{3}{2}\displaystyle \frac{\zeta (t-{t}_{0})}{{\tau }_{\mathrm{rh},0}}+1\right)}^{-1/3}.\end{eqnarray} \tag{ 11 }$$

Here t₀ is the time at which the first black holes begin to heat the cluster, and ζ ≃ 0.1 (Gieles et al. 2011; Alexander & Gieles 2012) is a dimensionless factor. The initial half-mass–radius, heating rate, and escape velocity are dependent only on the cluster mass, M_cl, and initial density, ρ_h,0,

$$\begin{eqnarray}&&{\dot{E}}_{0}\simeq 2.3\times {10}^{5}\,{M}_{\odot }{\left({\mathrm{kms}}^{-1}\right)}^{2}\,{\mathrm{Myr}}^{-1}{M}_{5}^{2/3}{\rho }_{5,0}^{5/6},\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{rh},0}\simeq 7.5\,\mathrm{Myr}\,{M}_{5}{\rho }_{5,0}^{-1/2},\end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray}&&{v}_{\mathrm{esc},0}\simeq 50\,{\mathrm{kms}}^{-1}\,{M}_{5}^{1/3}{\rho }_{5,0}^{1/6},\end{eqnarray} \tag{ 12c }$$

where M₅ ≡ M_cl/10⁵ M_⊙ and ρ_5,0 ≡ ρ_h,0/10⁵ M_⊙ pc⁻³, and ${\rho }_{{\rm{h}}}=3{M}_{\mathrm{cl}}/8\pi {r}_{{\rm{h}}}^{3}$. All the quantities are dependent on only the initial density of the cluster and the cluster mass.

Within our cluster model, we must first create a 1G black hole population. 1G black holes' masses are assumed to follow the inferred Power Law + Peak model (Talbot & Thrane 2018) from the second gravitational-wave transient catalog from the LIGO/Virgo Collaboration (Abbott et al. 2021a), with an additional strict mass cutoff at 45 M_⊙. This model is a combination of a power-law describing low-mass black holes, and a Gaussian peak at higher masses motivated by the prediction of a pair-instability supernova upper-mass gap precluding the formation of black holes with masses between ∼45 and ∼130 M_⊙ (Barkat et al. 1967; Fryer et al. 2001; Heger & Woosley 2002; Belczynski et al. 2016; Spera & Mapelli 2017; Stevenson et al. 2019).

All 1G black holes are considered to be initially nonspinning, motivated by studies which indicate that isolated black holes are likely to form with small spins (χ ≲ 0.01; Fuller & Ma 2019). The total mass of the 1G black holes within the cluster is fixed to 2% of the cluster mass (Löckmann et al. 2010; Antonini et al. 2019; Kremer et al. 2020b). Additionally, each black hole is initialized with a natal supernova kick velocity (Fryer & Kalogera 2001; Hobbs et al. 2005), though this has a minimal effect on the initial population. Finally, the black holes are deposited within the cluster at the initial half-mass–radius, r_h,0. Each black hole is expected to settle within the core of the cluster on the dynamical friction timescale, τ_df(r_h,0, t₀) (Binney & Tremaine 2011),

$$\begin{eqnarray}&&{\tau }_{\mathrm{df}}(r,t)\simeq 0.346\displaystyle \frac{{r}^{2}{v}_{\mathrm{esc}}(t)}{{GM}\mathrm{ln}\,{\rm{\Lambda }}},\,\end{eqnarray} \tag{ 13 }$$

where $\mathrm{ln}\,{\rm{\Lambda }}\simeq 10$ and assumes the King (1966) cluster model to relate the cluster's velocity dispersion to escape velocity.

Once the first black holes settle within the cluster's core, a black hole binary might form through dynamical three-body interactions (Ivanova et al. 2005; Morscher et al. 2015; Park et al. 2017). We let only one black hole binary exist at any one time such that it dominates the heating of the cluster. The binary is formed according to a mass distribution given by $\propto {\left({M}_{1}+{M}_{2}\right)}^{4}$ (O'Leary et al. 2016). The required heating rate of the cluster is balanced with the loss of energy from the binary in the core of the cluster (Antonini et al. 2019; Kremer et al. 2019). The timescale during which dynamical interactions dominate the energy flow of the cluster is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{dyn}}\simeq \displaystyle \frac{{{GM}}_{1}{M}_{2}}{2{a}_{m}}{\dot{E}}^{-1}(t),\,\end{eqnarray} \tag{ 14 }$$

where ${a}_{m}=\max ({a}_{\mathrm{ej}},{a}_{\mathrm{GW}})$. Here a_ej is the binary separation at which the binary is ejected, and a_GW is the separation at which gravitational-wave radiation dominates. Since we are focusing on denser stellar environments, the majority of binary systems are likely to merge within the cluster rather than be ejected. The ejection of the interloper black holes is also incorporated (see Appendix A.3). We calculate the separation at which gravitational-wave radiation dominates by equating the separation evolution due to dynamical interactions and gravitational-wave emission (Peters 1964). We compute the timescale for the binary to merge due to gravitational-wave radiation as ${\tau }_{\mathrm{GW}}={a}_{m}/| {\dot{a}}_{\mathrm{GW}}|$.

At this point in the binary's evolution, we modify the black holes' masses and spins according to Equations (4) and (3) to incorporate the effects of superradiance, ensuring each black hole's age exceeds to the growth timescale of the boson cloud. Our treatment of superradiance does not consider the impacts of dense cluster environments or binary systems on the spin-down process. It might be expected that close, hyperbolic encounters within the cluster could tidally disrupt the boson cloud. However, we can estimate the typical close-encounter timescale (Malmberg et al. 2007; Binney & Tremaine 2011) for clusters considered within this manuscript. To estimate the possibility of tidal disruption of the boson clouds, we note that the boson cloud radius is r_cloud ≳ α⁻² r_BH (Arvanitaki et al. 2010; Arvanitaki & Dubovsky 2011). Furthermore, since the boson cloud growth follows Equation (8), we are only interested in systems where α ≳ 10⁻² so that the black hole can sufficiently spin down. Therefore, we can consider the limiting case where α ∼ 10⁻² and r_cloud ∼ 10⁴ r_BH. We assume typical values of a half-mass–radius of 1 pc, and a cluster mass of 10⁷ M_⊙. Under these assumptions, a 50 M_⊙ black hole would undergo a close encounter on the timescale of ∼400 Myr. This is longer than both the typical timescale required for the majority of black holes to be spun down via superradiance (see Figure 1; Arvanitaki et al. 2010; Arvanitaki & Dubovsky 2011; Baryakhtar et al. 2017; Ng et al. 2021a) and the typical timescale of the black hole forming a binary (Binney & Tremaine 2011; Antonini et al. 2019). Furthermore, systems with higher values of α will have less extended boson clouds, thereby reducing the possibility of tidal disruption. Therefore, close encounters within the dense stellar environment are very unlikely to impact the spin-down process.

There is also the possibility that interactions between the boson clouds and the binary system can influence the evolution of the black hole–boson cloud system (Arvanitaki & Dubovsky 2011; Baumann et al. 2019; Berti et al. 2019; Zhang & Yang 2019; Baumann et al. 2020; Ng et al. 2021a). Of particular interest are the phenomena of "floating orbits" (Zhang & Yang 2019; Baumann et al. 2020) and angular momentum transfer due to level mixing (Baumann et al. 2019; Berti et al. 2019; Zhang & Yang 2020). Floating orbits arise from a driven transition between a growing and decaying mode of the boson cloud, which can increase the angular momentum of the orbiting black holes (Zhang & Yang 2019). However, it has been shown that such orbits exist for up to ∼100 yr for the intermediate mass ratio systems which arise during hierarchical black hole growth (see Figures 16 and 17 in Baumann et al. 2020). This is a negligible correction to the dynamical merger timescales considered above. In addition, the impact of cloud disruption during the binary black hole merger event might spoil the statement that black holes possess low spins prior to merging. However, studies have shown that cloud disruption might actually lead to reduced black hole spin (Baumann et al. 2019; Berti et al. 2019; Zhang & Yang 2019; Baumann et al. 2020; Zhang & Yang 2020). Due to the selection rules governing the behavior of a black hole–boson cloud system, typically boson cloud states with negative angular momentum fall in toward the black hole (Zhang & Yang 2020; Ng et al. 2021a), which leads to the further reduced spin. Ultimately, there is rich phenomenology associated with the behavior of binary systems in the presence of ultralight bosons. However, their impact on hierarchical growth appears minimal in the regime considered and are therefore ignored for the remainder of the study.

A vital component to the semi-analytic model is the computation of the black hole merger remnant properties. In particular, due to the conservation of linear momentum, the final remnant black hole receives a kick from the anisotropic emission of gravitational waves (Merritt et al. 2004; Campanelli et al. 2006, 2007a; González et al. 2007a; Campanelli et al. 2007b; González et al. 2007b; Lousto et al. 2010; Lousto & Zlochower 2013; Varma et al. 2019, 2020). This recoil velocity, v_k , determines whether the remnant remains in the core, is ejected from the core and has to settle through dynamical friction, or is ejected from the cluster entirely. We utilize the Precession code (Gerosa & Kesden 2016) to determine the final mass (Barausse et al. 2012), spin (Barausse & Rezzolla 2009), and recoil velocity of the remnant black hole. The details of the calculation of the recoil velocity are outlined in Appendix B.

The salient features of the recoil velocity calculation for the study of hierarchical mergers arise from black hole binary spins and mass ratios. It has been found that the largest kicks are typically a result of special spin configurations known as superkicks (Campanelli et al. 2007a; González et al. 2007a) and hang-up kicks (Lousto & Zlochower 2011, 2013). These spin configurations can lead to recoil velocities up to ∼5000 km s⁻¹. Conversely, nonspinning binaries typically lead to the smallest recoil velocities. For example, the largest recoil velocity a nonspinning binary can produce is only ∼170 km s⁻¹ for q ≃ 1/3 (González et al. 2007b). In order to quantify the binary's total spin and explore its contribution to hierarchical mergers, we introduce a mass-weighted spin magnitude:

$$\begin{eqnarray}&&\langle \chi \rangle \equiv \displaystyle \frac{{\chi }_{1}+{q}^{2}{\chi }_{2}}{{\left(1+q\right)}^{2}}=\displaystyle \frac{{M}_{1}^{2}{\chi }_{1}+{M}_{2}^{2}{\chi }_{2}}{{\left({M}_{1}+{M}_{2}\right)}^{2}}.\end{eqnarray} \tag{ 15 }$$

After calculating the properties of the merger remnant, we determine if or where it should be reintroduced into the black hole population. If the recoil velocity exceeds the escape velocity of the cluster at the time of the merger, v_k ≥ v_esc(t_merge), then the remnant is ejected from the cluster. Alternatively, if v_k < v_esc(t_merge), the remnant is retained and placed at a radius (Antonini et al. 2019)

$$\begin{eqnarray}&&{r}_{\mathrm{in}}\simeq {r}_{{\rm{h}}}({t}_{\mathrm{merge}})\sqrt{\tfrac{{v}_{\mathrm{esc}}^{4}({t}_{\mathrm{merge}})}{{\left({v}_{\mathrm{esc}}^{2}({t}_{\mathrm{merge}})-{v}_{k}^{2}\right)}^{2}}-1}.\end{eqnarray} \tag{ 16 }$$

If r_in < r_h(t_merger)/10, a conservative estimate of the core radius (Hurley 2007; Madrid et al. 2012), then the remnant black hole remains in the core. Otherwise, the binary must first settle in the core due to dynamical friction (see Equation (13)). After the merger remnant has either been ejected (over a period τ_dyn) or retained (over a period of τ_dyn + τ_GW) either in the core or in the cluster, we repeat the process outlined above by generating a new black hole binary to support the heating rate condition. The simulation is concluded when either only one black hole remains, too few black holes remain for dynamical hardening, or 13 Gyr have passed.

In Figure 2 we plot the total mass of the binary black hole system as a function of the coalescence time. These merger evolution tracks are plotted for four different ultralight boson masses (and for the case of no ultralight boson) using two different initial cluster densities. The shape of the points contained to retained (blue squares) or ejected (red crosses) binary systems. Different initial densities lead to different escape velocities and dynamical friction timescales. The denser cluster (left), with an initial density ρ_h,0 = 1.9 × 10⁶ M_⊙ pc⁻³, has an initial escape velocity v_esc,0 = 391.8 km s⁻¹, whereas the less dense cluster (right) only has an initial escape velocity v_esc,0 = 224.2 km s⁻¹.

Zoom In Zoom Out Reset image size

There are two distinct epochs during a black hole population's evolution in a cluster. Initially, the black hole population undergoes random, dynamical interactions which lead to mergers and formation of 2G or third-generation (3G) black holes. Since the black hole binaries formed early within the cluster have similar total masses, the chance any two black holes become a binary is small, though heavier black holes are slightly more likely to form binaries (O'Leary et al. 2016). This epoch is categorized by remnant masses typically less than ∼200 M_⊙ in the early cluster evolution. Eventually, however, one black hole tends to become sufficiently massive to dominate, leading to the second evolutionary phase. Since the heavy black hole will primarily form binaries with much smaller black holes, the resulting small mass ratio can lead to significantly reduced gravitational-wave recoil velocities (González et al. 2007b). This black hole therefore seeds hierarchical growth in the cluster. This period of evolution is seen by a track of increasing total binary mass with the merger coalescence time in Figure 2.

In the left column of Figure 2, we show the binary black hole mergers which occur in a cluster simulation with M_cl = 1.1 × 10⁷ M_⊙ and ρ_h,0 = 1.9 × 10⁶ M_⊙ pc⁻³. We see that all simulations are able to support hierarchical black hole growth. However, the existence of an ultralight boson at some masses can change the time at which hierarchical growth starts to occur. This is due to the spin down of 2G+1G and 2G+2G black hole mergers (see Figure 1) when bosons of masses ∼2 × 10⁻¹⁴–2 × 10⁻¹³ eV exist, helping retain more merger remnants. By retaining more massive black holes, the formation of a hierarchical black hole seed is more likely to occur earlier in the simulation. The simulation result in the top panel (μ = 7.6 × 10⁻¹⁵ eV) is similar to the case when no ultralight boson exists because the boson cloud instability growth rate is significantly reduced for stellar-mass black holes. However, there is some indication of spin-down from superradiance in the high-mass, hierarchical-growth regime.

On the right of Figure 2, we show the results from a cluster with M_cl = 1.1 × 10⁷ M_⊙ and ρ_h,0 = 6.6 × 10⁴ M_⊙ pc⁻³. Only one particular ultralight boson, with a mass μ = 4.1 × 10⁻¹⁴ eV, can facilitate hierarchical black hole growth. The clear difference between this simulation and the remaining four is that the binary systems with a total mass ∼80–200 M_⊙ have mass-weighted spins 〈χ〉 ≲ 0.15. The important result from these individual simulations is that the presence of ultralight bosons can impact what cluster properties support hierarchical growth. For the remainder of the manuscript, we use the scenario where μ = 4.1 × 10⁻¹⁴ eV as our best-case scenario example.

In Figure 3, we plot the distributions of the binary black hole merger properties for all 2G+2G and 2G+1G mergers for cluster simulations presented in Figure 2. We use the full set of 30 simulations for each set of initial cluster parameters, and show the distributions for the μ = 4.1 × 10⁻¹⁴ eV, 2.2 × 10⁻¹³ eV and no ultralight boson simulations. The left plot shows the distributions in the case where all clusters lead to hierarchical growth. The bulk of the kick-velocity distributions lay below the initial escape velocity, v_esc,0, of the cluster. Crucially, in the μ = 4.1 × 10⁻¹⁴ eV results (orange), we see two spin populations corresponding to systems which have undergone substantial superradiance (〈χ〉 ≲ 0.15), and a population which has not (〈χ〉 ∼0.3–0.45). The low-spin population is entirely retained within the cluster.

Zoom In Zoom Out Reset image size

In the right panel, only the μ = 4.1 × 10⁻¹⁴ eV (orange) distribution corresponds to a system capable of hierarchical growth. In these results, the majority of black holes are spun down to the low-spin population, and are still retained. Conversely, in the cases of μ = 2.2 × 10⁻¹³ eV and no ultralight boson, the black holes are not spun down enough to reduce the velocity distribution significantly. Therefore the majority of 2G black holes are ejected upon merging with another black hole.

For each set of simulations with a given ultralight boson mass and cluster properties, we compute the fraction of the repeated 30 simulations which result in the formation of a black hole with a mass ≥1000 M_⊙. ⁸ We denote this fraction f_hier (Antonini et al. 2019). We use f_hier as a proxy for the fraction of simulations leading to hierarchical growth within the cluster. We compute the region within which more than 50% of the simulations for a given initial effective radius, R_eff, and cluster mass undergo hierarchical growth, i.e., f_hier > 0.5. We compare the region of the effective radius–cluster mass (R_eff − M_cl) parameter space where clusters support hierarchical growth with a population of globular clusters (Baumgardt et al. 2019) and nuclear star clusters (Georgiev et al. 2016). For the nuclear star cluster population, we plot the 152 nuclear star clusters with cluster mass uncertainty less than an order of magnitude, and effective radius estimates with only positive radius support. However, we use the full population of 228 nuclear star clusters from Georgiev et al. (2016) in future calculations. These results are presented in Figure 4 for μ = 4.1 × 10⁻¹⁴ eV (orange), as well as simulations where no ultralight boson is present (hatched blue).

Zoom In Zoom Out Reset image size

The lower-right region, corresponding to denser and heavier clusters, unsurprisingly supports hierarchical growth. Our bounds are similar to those presented in Antonini et al. (2019) for the scenario with no ultralight boson, though here we empirically compute the region of parameter space supporting hierarchical growth. The inclusion of superradiance extends the region within which hierarchical growth is supported further into the astrophysical population of clusters, though it still does not permit hierarchical growth in globular clusters. Importantly, the orange contour includes a higher percentage of the nuclear star cluster population. Furthermore, the highlighted contour region only corresponds to the initial conditions, whereas it is inevitable that the clusters have evolved since their formation. To illustrate this, we evolve clusters with initial conditions along the contour to an age of 10 Gyr, shown by the dashed curves in Figure 4. A larger proportion of the nuclear star cluster population is contained under this contour. However, given that the age of present-day clusters is not well known, it is difficult to interpret whether hierarchical growth might have previously occurred in observed clusters.

In order to interpret the significance of the difference between the hierarchical-growth regions for simulations with and without superradiance, we compute the fraction of nuclear star clusters from Georgiev et al. (2016), which presently would support hierarchical growth, f_NSC, i.e., those nuclear star clusters that lie within the contours in Figure 4. See Appendix C for details about the calculation and uncertainty estimation. As an example, for the contours presented in Figure 4, the fraction of present-day nuclear star clusters which can sustain hierarchical growth is ≈4.5% in the absence of ultralight bosons and ≈7% in the presence of a boson with a mass of μ = 4.1 × 10⁻¹⁴ eV. The fraction in the absence of ultralight bosons is less than the value presented in Antonini et al. (2019 ; f_NSC ≃ 10%) as we empirically generate the contour from simulations, as well as explicitly integrate under it to calculate f_NSC. The values of f_NSC for different boson masses are presented in Figure 5. The fraction peaks at μ ≃ 4.1 × 10⁻¹⁴ eV, corresponding to a ≈60% increase in the number of clusters capable of supporting hierarchical growth presently. There is also an increased number of clusters currently capable of supporting hierarchical growth if bosons with masses μ ∼ 2 × 10⁻¹⁴–2 × 10⁻¹³ eV exist.

Zoom In Zoom Out Reset image size

To understand the distribution of binary black hole mergers in present-day clusters, we generate a synthetic population, which we evolve to a similar state as the observed population. We generate cluster masses from ${\mathrm{log}}_{10}({M}_{\mathrm{cl}}/{M}_{\odot })\sim { \mathcal N }(\mu =6.3,\sigma =0.8)$ and initial densities from ${\mathrm{log}}_{10}({\rho }_{{\rm{h}},0}/{M}_{\odot }\,{\mathrm{pc}}^{-3})\sim { \mathcal N }(\mu =4.5,\sigma =1.5)$. Each cluster is evolved to an age drawn at random from ${\mathrm{log}}_{10}({t}_{\mathrm{age}}/\mathrm{Gyr})\sim { \mathcal N }(\mu =0.31,\sigma =1)$. These distributions are chosen such that the evolved population visually appears similar to the observed nuclear star cluster properties from Georgiev et al. (2016). We evolve 4.8 × 10⁴ clusters with properties drawn randomly from these distributions, for scenarios where μ = 4.1 × 10⁻¹⁴ eV and no ultralight boson exists. The final merger is then included in the merger population. The individual black hole properties from both scenarios, along with the spin-down limits set by μ = 4.1 × 10⁻¹⁴ eV, are shown in Figure 6. The densities of the merger fraction, f_m , in mass and spin are also presented.

Zoom In Zoom Out Reset image size

From the one-dimensional marginal distributions, the inclusion of superradiance from bosons with μ = 4.1 × 10⁻¹⁴ eV clearly leads to an enhancement of the number of black holes with masses above 100 M_⊙ in the population of merging binary systems. Furthermore, a distinct 2G population with spins ∼0.7 is observed. This is present regardless of the presence of the ultralight boson, as the majority of black holes in this population are not old enough to undergo any significant superradiant boson cloud growth. Additionally, the superradiant instability leads to a strong correlation between the black hole spin and mass in the mass range from ∼80–300 M_⊙, as the black holes are maximally spun down through the formation of a boson cloud. Finally, there is a trend of reduced spins as the black holes become heavier. This is due to the fact that these black holes are formed through high mass ratio mergers, which typically reduce the spin of the merger remnant. This is present in both scenarios.

In Figure 7, we plot the mass–spin distribution while taking into account gravitational-wave detector selection biases. The primary black hole mass–spin 1σ and 2σ two-dimensional credible intervals for GW190412, GW190517, and GW190521, as well as the secondary black hole intervals for GW190521, are also shown (shown as dashed contours; Abbott et al. 2020a, 2020b, 2020d).

Zoom In Zoom Out Reset image size

By comparing the gravitational-wave detection-weighted binary black hole populations, we see a number of distinct features. First, the observed population is restricted to black holes with individual masses ∼5–200 M_⊙. These selection biases are well known and consistent between both scenarios presented here (Messenger & Veitch 2013; Farr 2019). Additionally, the black hole mass–spin correlation is still observable in the population in the presence of ultralight bosons. Finally, there is an increase in the number of heavier black hole detections in the ultralight boson scenario. To quantify the increased number of heavier black hole merger populations, we compute the merger rates of both the total population, ${ \mathcal R }$, as well as the intermediate-mass (100 <M/M_⊙ < 1000) population; please refer to Appendix D for details. We summarize the merger rates calculated and present the observed rates (Abbott et al. 2020b, 2021a, 2021e) in Table 1. The values calculated in the absence of ultralight bosons are consistent with the results from Antonini et al. (2019).

Table 1. Total and Intermediate Mass Black Hole Merger Rates from Observations of Gravitational Waves from Stellar-mass Black Hole Mergers (Abbott et al. 2020b, 2021a, 2021e) Compared to the Calculated Merger Rates From Our Simulated Population (with 68% Credible Intervals)

	Total Merger Rate, ${ \mathcal R }$	IMBH Merger Rate, ${{ \mathcal R }}_{\mathrm{IMBH}}$
	(Gpc−3 yr−1)	(Gpc−3 yr−1)
Observed	${23.9}_{-8.6}^{+14.3}$	${0.08}_{-0.07}^{+0.19}$
No ultralight boson	≈1.2 fnMBH	≈0.04 fnMBH
μ = 4.1 × 10−14 eV	≈0.9 fnMBH	≈0.08 fnMBH

Note. The observed intermediate-mass black hole (IMBH) merger rate is a single-event rate determined from GW190521 (Abbott et al. 2020b, 2021e). While the total merger rate cannot be justified by hierarchical mergers in nuclear star clusters, the observed rate of intermediate-mass mergers might be explained by such events. The inclusion of superradiance due to the existence of bosons with μ = 4.1 × 10⁻¹⁴ eV leads to a doubling in the IMBH merger rate for our simulated population. Here, f_nMBH corresponds to the fraction of galaxies without a central massive black hole.

Download table as:  ASCIITypeset image

From these calculations, nuclear star clusters, regardless of whether superradiance occurs, cannot explain the total observed merger rate from gravitational-wave detectors (Abbott et al. 2020d). This result is expected, and it is anticipated that field binaries and/or dynamical mergers in globular clusters can explain the bulk of the observed binary black hole mergers (Dominik et al. 2013; Eldridge et al. 2019; Neijssel et al. 2019; Kremer et al. 2020b; Mapelli et al. 2021; Santoliquido et al. 2020; Rodriguez et al. 2021; Wong et al. 2021; Zevin et al. 2021). However, the inclusion of superradiance leads to a doubling of the merger rate of intermediate-mass black hole (IMBH) mergers to ≈0.08 f_nMBH Gpc⁻³ yr⁻¹ for μ = 4.1 × 10⁻¹⁴ eV bosons. Both inferred values of ${{ \mathcal R }}_{\mathrm{IMBH}}$ are consistent with the single-event merger rate determined for GW190521 (Abbott et al. 2020b, 2021e). Future gravitational-wave observations of heavy-mass binary black hole systems will continue to reduce the uncertainty in the merger rate of these heavier systems.

Finally, within our simulated population, the interpretations of gravitational-wave observations remain unchanged in the presence of ultralight bosons. For GW190521 (Abbott et al. 2020a, 2020b), we conclude that the system is likely a 2G+2G binary black hole merger (discussed in Gayathri et al. 2020, Kimball et al. 2021, and Romero-Shaw et al. 2020). Recently, Ng et al. (2021b) confirmed that the primary black holes from GW190412 (Abbott et al. 2020c) and GW190517 (Abbott et al. 2020d) exclude non-self-interacting ultralight scalar bosons with masses 1.3–2.7 × 10⁻¹³ eV. From our simulations, we suspect GW190517 might be a 2G+1G merger too light to be spun down through superradiance. This is tentatively supported by the results from Kimball et al. (2021). GW190412 is inconsistent with our simulated population irrespective of the presence of ultralight bosons—likely a direct result of our mass cutoff or spin distribution used. We can nevertheless interpret GW190412 as a direct result of hierarchical and/or dynamical formation (as in Abbott et al. 2020c, Gerosa et al. 2020, Rodriguez et al. 2020, Safarzadeh & Hotokezaka 2020, and Zevin et al. 2020). For the most impactful boson masses from our analysis, GW190412 and GW190517 do not provide any information about their presence.

Currently, the existence of ultralight bosons is purely hypothetical. In the mass range we are focused on in this manuscript, μ ∼ 2 × 10⁻¹⁴–2 × 10⁻¹³ eV, their detection would likely be made by gravitational-wave detectors via either the direct observation of continuous gravitational waves, as a stochastic background, or as a feature of the population (Brito et al. 2017b; Isi et al. 2019; Tsukada et al. 2019; Sun et al. 2020; Ng et al. 2021a, 2021b; Tsukada et al. 2021). Other ultralight boson searches would likely cover a different boson-mass range, which impacts black holes of different masses. The detection of ultralight bosons in this energy range is likely only possible with future generations of gravitational-wave detectors (Isi et al. 2019; Sun et al. 2020). Isi et al. (2019) found that the horizon distance for continuous gravitational-wave emission from boson clouds with masses 10⁻¹⁴–10⁻¹³ eV is ≲100 Mpc for aLIGO detectors (Aasi et al. 2015). 3G gravitational-wave detectors such as the Cosmic Explorer (Reitze et al. 2019) or Einstein Telescope (Maggiore et al. 2020) would be able to detect continuous gravitational-wave emission from the desired ultralight bosons out to a horizon of ∼1 Gpc. Analyzing the mass–spin distributions of the stellar-mass black holes may be a more successful method for indirectly determining the existence of ultralight bosons (Ng et al. 2021a, 2021b). Ng et al. (2021a) found that ${25}_{-15}^{+95}$ or ${80}_{-70}^{+210}$ gravitational-wave detections with signal-to-noise ratios exceeding 30 were required to determine the presence of ultralight bosons with μ = 10⁻¹³ eV, for two fiducial 1G spin distributions. It is unlikely that these detection numbers can be applied to a ∼10⁻¹⁴ eV boson detection prediction. Finally, even though we find interesting features beyond the spin-down exclusion regions, they will likely be masked by current uncertainties in the modeling of binary black hole cluster dynamics (Rodriguez et al. 2018; Antonini et al. 2019; Rodriguez et al. 2019; Mapelli et al. 2021). Therefore, in the near future, the detection of ultralight bosons will still rely on a direct gravitational-wave measurement or on spin-population studies of gravitational-wave events (Brito et al. 2017b; Dergachev & Papa 2019; Isi et al. 2019; Palomba et al. 2019; Tsukada et al. 2019; Sun et al. 2020; Zhu et al. 2020; Abbott et al. 2021d; Ng et al. 2021a, 2021b; Tsukada et al. 2021). Our study has explored the possible differences superradiance might bring to the observed population of binary black hole mergers.

Using Hénon's principle, the heating rate from binary black hole systems in the core of the cluster is balanced against the energy flow into the cluster as a whole (Hénon 1961, 1975; Gieles et al. 2011; Breen & Heggie 2013). The rate of energy generation, $\dot{E}$, is therefore directly related to the properties of the cluster,

$$\begin{eqnarray}&&\dot{E}=\zeta \displaystyle \frac{| E| }{{\tau }_{\mathrm{rh}}},\end{eqnarray} \tag{ A1 }$$

where $E\simeq -0.2{{GM}}_{\mathrm{cl}}^{2}/{r}_{{\rm{h}}}$ is the total energy of the cluster, M_cl is the cluster mass, r_h is the half-mass–radius, ζ ≃ 0.1 (Gieles et al. 2011; Alexander & Gieles 2012) is a dimensionless factor, and τ_rh is the relaxation timescale of the cluster, which can be approximated as (Spitzer & Hart 1971; Antonini et al. 2019)

$$\begin{eqnarray}&&{\tau }_{\mathrm{rh}}=0.138\sqrt{\displaystyle \frac{{M}_{\mathrm{cl}}{r}_{{\rm{h}}}^{3}}{G}}\displaystyle \frac{1}{\langle {m}_{* }\rangle \psi \mathrm{ln}{\rm{\Lambda }}}.\end{eqnarray} \tag{ A2 }$$

Here, 〈m_*〉 is the mean stellar and compact object mass within the cluster and $\mathrm{ln}{\rm{\Lambda }}\simeq 10$ is the Coulomb logarithm. The factor ψ is related to the distribution of masses within the half-mass–radius and approximated in Spitzer & Hart (1971) as $\psi =\langle {m}_{* }^{2.5}\rangle /\langle {m}_{* }{\rangle }^{2.5}$. For the remainder of the manuscript, we follow Antonini et al. (2019) and assume low-mass stars dominate the mass distribution of the cluster, such that 〈m_*〉 = 0.6 M_⊙, and ψ ≃ 5.

Finally, by defining the half-mass density, ${\rho }_{{\rm{h}}}=3{M}_{\mathrm{cl}}/8\pi {r}_{{\rm{h}}}^{3}$, and escape velocity, ${v}_{\mathrm{esc}}\propto \sqrt{{M}_{\mathrm{cl}}/{r}_{{\rm{h}}}}$, we can express the initial heating rate, relaxation time, and escape velocity at t = t₀ uniquely in terms of the initial cluster density and its mass. From Equations (A6) and (A7) (Antonini et al. 2019),

$$\begin{eqnarray}&&{\dot{E}}_{0}\simeq 2.3\times {10}^{5}\,{M}_{\odot }{({\mathrm{kms}}^{-1})}^{2}\,{\mathrm{Myr}}^{-1}{M}_{5}^{2/3}{\rho }_{5,0}^{5/6},\end{eqnarray} \tag{ A3a }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{rh},0}\simeq 7.5\,\mathrm{Myr}\,{M}_{5}{\rho }_{5,0}^{-1/2},\end{eqnarray} \tag{ A3b }$$

$$\begin{eqnarray}&&{v}_{\mathrm{esc},0}\simeq 50\,{\mathrm{kms}}^{-1}\,{M}_{5}^{1/3}{\rho }_{5,0}^{1/6},\end{eqnarray} \tag{ A3c }$$

where M₅ ≡ M_cl/10⁵ M_⊙ and ρ_5,0 ≡ ρ_h,0/10⁵ M_⊙ pc⁻³. Here, t₀ corresponds to the time at which the first black hole binaries begin to heat the cluster (i.e., when they have entered the cluster's core; Breen & Heggie 2013). Furthermore, we have implicitly assumed that mass loss of the cluster is a negligible contribution to the evolution of the cluster. Below, we explicitly require a constant M_cl in the calculation of the cluster model's evolution. The proportionality constant for the initial escape velocity is derived from the King (1966) cluster model with W₀ = 7, where W₀ is the central value of the dimensionless form of the cluster's gravitational potential, uniquely describing the potential's shape. This model is used throughout for the calculation of the dynamical features of the clusters.

To derive the time dependence of the cluster properties, we assume no mass loss (constant M_cl) and that the cluster is always in virial equilibrium. Under these assumptions, the rate of expansion of the cluster is related to the heating rate as (Hénon 1965)

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{r}}_{{\rm{h}}}}{{r}_{{\rm{h}}}}=\displaystyle \frac{\dot{E}}{| E| }=\displaystyle \frac{\zeta }{{\tau }_{\mathrm{rh}}},\,\end{eqnarray} \tag{ A4 }$$

following Equation (A6). Fixing the cluster mass, we have ${\tau }_{\mathrm{rh}}\propto {r}_{h}^{3/2}$, and solve Equation (A4) for the time evolution of the half-mass–radius from t₀,

$$\begin{eqnarray}&&{r}_{{\rm{h}}}(t)={r}_{{\rm{h}},0}{\left(\displaystyle \frac{3}{2}\displaystyle \frac{\zeta (t-{t}_{0})}{{\tau }_{\mathrm{rh},0}}+1\right)}^{2/3},\end{eqnarray} \tag{ A5 }$$

when t > t₀. When t < t₀, the cluster has not yet extracted energy from the binary population and is assumed to have the initial half-mass–radius. Equation (A5) can be directly substituted into the expressions for the heating rate and escape velocity to find

$$\begin{eqnarray}&&\dot{E}(t)={\dot{E}}_{0}{\left(\displaystyle \frac{3}{2}\displaystyle \frac{\zeta (t-{t}_{0})}{{\tau }_{\mathrm{rh},0}}+1\right)}^{-5/3},\end{eqnarray} \tag{ A6 }$$

and

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}(t)={v}_{\mathrm{esc},0}{\left(\displaystyle \frac{3}{2}\displaystyle \frac{\zeta (t-{t}_{0})}{{\tau }_{\mathrm{rh},0}}+1\right)}^{-1/3}.\end{eqnarray} \tag{ A7 }$$

After a binary black hole system has formed, the system undergoes dynamical hardening (Leigh et al. 2018). This process tightens the binary orbit while heating the cluster via Hénon's principle (Hénon 1975). This implies that the rate of energy loss from the binary is approximately equal to the heating rate of the cluster in Equation (A6):

$$\begin{eqnarray}&&\dot{E}(t)\simeq -{\dot{E}}_{\mathrm{bin}}.\end{eqnarray} \tag{ A8 }$$

The total energy of the binary during the process of dynamical hardening is

$$\begin{eqnarray}&&{E}_{\mathrm{bin}}=-\displaystyle \frac{{{GM}}_{1}{M}_{2}}{2a},\end{eqnarray} \tag{ A9 }$$

where a is the semimajor axis of the binary's orbit. Therefore, differentiating both sides of Equation (A9) and applying Equation (A8) to express the semimajor axis evolution of the binary as a function of the heating rate of the cluster leads to (Antonini et al. 2019)

$$\begin{eqnarray}&&\dot{a}=-\displaystyle \frac{2{a}^{2}}{{{GM}}_{1}{M}_{2}}\dot{E}(t).\end{eqnarray} \tag{ A10 }$$

We integrate this expression with respect to the semimajor axis from the initial semimajor axis of formation, a_h , to the minimum semimajor axis for which dynamical hardening dominates, a_m . The heating rate is assumed to be approximately constant over the binary's evolution. This defines the timescale of dynamical hardening for a binary to be

$$\begin{eqnarray}&&{\tau }_{\mathrm{dyn}}=-{\int }_{{a}_{m}}^{{a}_{h}}\displaystyle \frac{{{GM}}_{1}{M}_{2}}{2{a}^{2}}{\dot{E}}^{-1}(t){\rm{d}}a\simeq \displaystyle \frac{{{GM}}_{1}{M}_{2}}{2{a}_{m}}{\dot{E}}^{-1}(t),\end{eqnarray} \tag{ A11 }$$

where we have used the fact that a_h ≫ a_m , as seen in Equation (14).

The lower bound semimajor axis, a_m , is given by

$$\begin{eqnarray}&&{a}_{m}=\max ({a}_{\mathrm{ej}},{a}_{\mathrm{GW}}),\end{eqnarray} \tag{ A12 }$$

where a_ej is the semimajor axis for which the binary system is ejected from the cluster through binary–single interactions, and a_GW is the semimajor axis at which the energy loss from dynamical hardening is equal to the energy loss from gravitational-wave emission. If the semimajor axis of the binary decreases below a_ej, then the binary is ejected from the cluster prior to merging (Quinlan 1996; Miller & Hamilton 2002; Miller & Lauburg 2009). The value of a_ej is given by (Antonini & Rasio 2016)

$$\begin{eqnarray}&&{a}_{\mathrm{ej}}\simeq 0.2\displaystyle \frac{{M}_{1}{M}_{2}}{{M}_{123}}\displaystyle \frac{{q}_{3}}{{v}_{\mathrm{esc}}^{2}(t)},\end{eqnarray} \tag{ A13 }$$

where ${M}_{123}={M}_{1}+{M}_{2}+\langle M{\rangle }_{\mathrm{core}}$ is the average total mass of the binary–single interaction, and ${q}_{3}=\langle M{\rangle }_{\mathrm{core}}/({M}_{1}+{M}_{2})$ is the mass ratio of the interaction. Since we do not track the individual trajectories and interactions within the cluster, we use the average black hole mass in the core (excluding the binary), $\langle M{\rangle }_{\mathrm{core}}$, to calculate a_ej.

If a_GW > a_ej, the binary will merge within the cluster. Due to the massive and dense clusters considered in this manuscript, this occurs for almost all binaries. The evolution of the semimajor axis through gravitational-wave emission is given by (Peters 1964)

$$\begin{eqnarray}&&{\dot{a}}_{\mathrm{GW}}=-\displaystyle \frac{65}{4}\displaystyle \frac{{G}^{3}{M}_{1}{M}_{2}\left({M}_{1}+{M}_{2}\right)}{{c}^{5}{a}^{3}{\left(1-{e}^{2}\right)}^{7/2}}g(e),\end{eqnarray} \tag{ A14 }$$

where e is the eccentricity of the binary, and

$$\begin{eqnarray}&&g(e)=1+\displaystyle \frac{73}{24}{e}^{2}+\displaystyle \frac{37}{96}{e}^{4}.\end{eqnarray} \tag{ A15 }$$

By equating Equations (A10) and (A14), solving for a leads to the expression (Antonini et al. 2019)

$$\begin{eqnarray}&&{a}_{\mathrm{GW}}^{5}=\displaystyle \frac{32}{5}\displaystyle \frac{{G}^{4}{M}_{1}^{2}{M}_{2}^{2}({M}_{1}+{M}_{2})g(e)}{{c}^{5}{\left(1-{e}^{2}\right)}^{7/2}}{\dot{E}}^{-1}(t).\end{eqnarray} \tag{ A16 }$$

For every binary formed, we generate an eccentricity from the thermal eccentricity distribution ∝ e (Jeans 1919). The timescale for gravitational-wave emission from the binary is then given by ${\tau }_{\mathrm{GW}}={a}_{m}/| {\dot{a}}_{\mathrm{GW}}|$. Just prior to merger, the effects of superradiance on the individual black holes are incorporated.

Although it is very unlikely that the binary system is ejected from binary–single interactions, this is not the case for interlopers. To account for this, the expected number of ejected interlopers as (Antonini et al. 2019)

$$\begin{eqnarray}&&{N}_{\mathrm{ej}}\approx 6\,\mathrm{ln}\left(\displaystyle \frac{{a}_{\mathrm{ej}}}{{q}_{3}^{2}{a}_{m}}\right),\end{eqnarray} \tag{ A17 }$$

where the interloper mass ratio, q₃, is calculated from the mean mass of black holes in the core. With the approximate number of ejected black hole interlopers calculated, we discard N_ej black holes randomly. Generally, this process would favor the ejection of low-mass black holes. However, given that the majority of black holes in the cluster are from the first generation, this averaged procedure does not change the expected results. The key difference with the inclusion of interloper ejection is a restriction on the upper mass of the black hole formed through hierarchical black hole mergers. Since this process removes a fraction of black holes from the cluster, in some simulations this limits the available black holes for hierarchical growth.