key: cord-0439688-4jrnsl1f
authors: Payne, Ethan; Sun, Ling; Kremer, Kyle; Lasky, Paul D.; Thrane, Eric
title: The imprint of superradiance on hierarchical black hole mergers
date: 2021-07-25
journal: nan
DOI: nan
sha: 5832e828bd49a0b40a8776079f29fb3b9bb850da
doc_id: 439688
cord_uid: 4jrnsl1f

Ultralight bosons are a proposed solution to outstanding problems in cosmology and particle physics: they provide a dark-matter candidate while potentially explaining the strong charge-parity problem. If they exist, ultralight bosons can interact with black holes through the superradiant instability. In this work we explore the consequences of this instability on the evolution of hierarchical black holes within dense stellar clusters. By reducing the spin of individual black holes, superradiance reduce the recoil velocity of merging binary black holes, which, in turn, increases the retention fraction of hierarchical merger remnants. We show that the existence of ultralight bosons with mass $ 2times10^{-14}lesssim mu/textrm{eV} lesssim2times10^{-13}$ would lead to an increased rate of hierarchical black hole mergers in nuclear star clusters. An ultralight boson in this energy range would result in up to $approx60%$ more present-day nuclear star clusters supporting hierarchical growth. The presence of an ultralight boson can also double the rate of intermediate mass black hole mergers to $approx0.08$,Gpc$^{-3}$,yr$^{-1}$ in the local Universe. These results imply that a select range of ultralight boson mass can have far-reaching consequences for the population of black holes in dense stellar environments. Future studies into black hole cluster populations and the spin distribution of hierarchically formed black holes will test this scenario.

A number of extensions to the Standard Model of particle physics propose the existence of a theoretical ultralight boson with a mass between µ ∼ 10 −33 -10 −10 eV (Arvanitaki et al. 2010; Arvanitaki & Dubovsky 2011; Ringwald 2013) . These particles can provide solutions to outstanding problems in cosmology and fundamental particle physics such as by being viable candidates for dark matter (Jaeckel & Ringwald 2010; Hui et al. 2017; Hu et al. 2000) or solving the strong charge-parity problem (Peccei & Quinn 1977a,b; Weinberg 1978) . While ultralight bosons are not expected to strongly interact with particles from the Standard Model (Dine et al. 1981; Shifman et al. 1980; Kim 1979) , the weak equivalence principle requires them to gravitate in a similar manner to visible matter (Detweiler 1980) .

If an ultralight bosonic wave exists in the vicinity of a spinning black hole, the wave can become gravitationally bound to the black hole. The bosonic field may be amplified by the extraction of rotational energy from the black hole through a process known as superradiance (Klein 1929; Dicke 1954; Zel'Dovich 1971 , 1972 Press & Teukolsky 1972; Bekenstein 1973; Bekenstein & Schiffer 1998; Brito et al. 2015 Brito et al. , 2020 . Superradiant amplification of the bosonic field can be unstable, leading to the exponential growth of the bosonic wave (Press & Teukolsky 1972) , forming a macroscopic quantum object, or boson cloud (e.g., Balakumar et al. 2020) . This superradiant instability occurs most rapidly when the Compton wavelength of the particle is comparable to the outer horizon radius of a spinning black hole, and only ceases when the angular frequency of the black hole's rotation equals the frequency of the bosonic wave (Arvanitaki et al. 2010; Bekenstein 1973; Brito et al. 2015; Bekenstein & Schiffer 1998) . The final astrophysical system, following unstable black hole superradiance, is a lighter black hole with a reduced spin, arXiv:2107.11730v1 [gr-qc] 25 Jul 2021 surrounded by a macroscopic boson cloud (Brito et al. 2017a ). This matter configuration can also lead to longlasting, nearly monochromatic gravitational-wave radiation from the rotation of the boson cloud-black hole system (Yoshino & Kodama 2014) , often referred to as a continuous gravitational-wave signal.

Directed searches for continuous waves from ultralight boson clouds around black hole merger remnants have not yet been undertaken because the signal is likely too weak to observe with current observatories (Isi et al. 2019; Sun et al. 2020) . Nonetheless, non-detections of an incoherent stochastic background from the continuouswave emission have placed tentative constraints on boson masses of (2-3.8) × 10 −13 eV given particular assumptions about the black hole formation spin distribution (Tsukada et al. 2019; Brito et al. 2017b; Tsukada et al. 2021) .

Furthermore, by studying the spins of the population of resolved binary black hole (BBH) mergers, constraints exclude boson masses between (1.3 -2.7) × 10 −13 eV assuming negligible boson self-interactions (Ng et al. 2021a,b) . Other observations of black holes such as the recent images of M87 * from the Event Horizon Telescope (Akiyama et al. 2019; Davoudiasl & Denton 2019) , and radial velocity and photometric data from Cygnus X-1 (Iorio 2008; Orosz et al. 2011; Reid et al. 2011; Middleton 2016; Isi et al. 2019) have excluded the presence of ultralight bosons in different regions of the boson mass parameter space. Ultimately, there is no strong evidence yet for either the existence or absence of ultralight bosons.

A population of low-spin black holes produced by superradiance would have wide-ranging astrophysical consequences. In dense stellar environments, such as nuclear star clusters or globular clusters, low-spin, firstgeneration (1G) black holes can merge hierarchically to form second-generation (2G) black holes (Zwart & McMillan 2000; Miller & Lauburg 2009; Downing et al. 2010; Rodriguez et al. 2015 Rodriguez et al. , 2016 Rodriguez et al. , 2018 Kremer et al. 2020a) . By decreasing blackhole spins through superradiance, the gravitationalwave recoil velocities are reduced (Campanelli et al. 2006 (Campanelli et al. , 2007a González et al. 2007b,a; Lousto et al. 2010; Lousto & Zlochower 2013; Varma et al. 2019) . Lower recoil velocities lead to a higher retention fraction of binary black hole merger remnants (Merritt et al. 2004; Varma et al. 2020) , and subsequent enhancement of hierarchical black hole growth as a result (Antonini et al. 2019; Rodriguez et al. 2018 Rodriguez et al. , 2019 .

The remainder of the manuscript is structured as follows. In Sec. 2, we discuss the theory of black hole superradiance in the context of scalar bosons and its impact on individual black holes. 1 We summarize our semi-analytic model for studying black hole mergers in dense star clusters (based on Antonini et al. 2019) , in Sec. 3. We present results from cluster simulations in relation to both individual clusters and the population as a whole in Sec. 4. Finally, we outline the implications of superradiance for recently observed black hole mergers, and the possibility for the detection of ultralight bosons in this boson mass regime in Sec. 5.

In this section we summarize the key expressions from Brito et al. (2020) to provide the relevant theory for including the effects of ultralight bosons in cluster simulations. We work with the analytic approximations for the evolution of boson clouds around spinning black holes (Ternov et al. 1978; Detweiler 1980; Baryakhtar et al. 2017) , as opposed to the coupled differential equations governing the black hole-boson cloud system assuming quasi-adiabatic evolution (Brito et al. 2015) . Though these approximations are accurate in the limit α 0.1, the errors are insignificant when extrapolated (Pani et al. 2012 ).

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

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,r + ≡ r + /r BH = 1+ 1 − χ 2 , where r BH = GM/c 2 where the characteristic black hole length scale which 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) . 1 We focus on scalar (spin-0) bosons (Damour et al. 1976; Ternov et al. 1978; Detweiler 1980; Yoshino & Kodama 2014) . However, the general conclusions can be applied to spin-1 (East 2017; Frolov et al. 2018; Siemonsen & East 2020 ) and spin-2 (de Rham 2014; Hinterbichler 2012) bosons if efficient black hole spindown is possible.

To highlight this comparison, we define the effective fine-structure constant for the "black hole atom" as the ratio of the characteristic length scale to the boson's Compton wavelength, λ = c/µ,

Large values of α lead to significant growth of the boson cloud. If the inequality in Eq. (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 Eq. (1) (Isi et al. 2019; Brito et al. 2017a ),

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 λ ∼ r 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) ,

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 )

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 (cf. Eq. 14 from Baryakhtar et al. 2017) , which we use to compute an approximate growth timescale around a black hole (Arvanitaki & Dubovsky 2011; Arvanitaki et al. 2010; Baryakhtar et al. 2017; Ng et al. 2021a) ,

The growth timescale scales as

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 Fig. 1 in Ng et al. (2021a) , Fig. 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 Eqs. (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 (Brito et al. 2015 (Brito et al. , 2017a Arvanitaki & Dubovsky 2011) . These exclusion regions are strongly dependent of the boson's mass and the evolution timescale of the black hole population. In Fig. 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.

To model the evolution of dense stellar clusters such as nuclear star clusters, we follow Antonini et al. (2019) which applies Hénon's principle (Hénon 1975) to simulate the evolution of dense stellar environments. In this section, we briefly summarize the key components of the cluster model and the evolution of the cluster's black hole population. Please refer to App. A for the full details of the model. Brito et al. 2017a ). The darker-coloured regions bounded by the solid curves correspond to black hole mass-spin parameter space within which any black hole would be spun down via superradiance after τ evol = 10 Myr. Since binary black hole merger delay time distribution has little support for systems merging within 10 Myr of formation (e.g. Britt et al. 2021) , we expect all black holes to allow for superradiant cloud growth for at least 10 Myr. Therefore, the observation of a black hole within the darker region of the parameter space can rule out some range of boson mass. The lighter region bounded by the dashed curves corresponds to when the black holes are given 10 Gyr to evolve.

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 (Hénon , 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

Here t 0 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 , 

where M 5 ≡ M cl /10 5 M and ρ 5,0 ≡ ρ h,0 /10 5 M pc −3 , and ρ h = 3M cl /8πr 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 first generation black hole population. First generation 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. 2020a) , with an additional strict mass cut-off 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 M and ∼ 130 M ( Barkat et al. 1967; Heger & Woosley 2002; Belczyn-ski 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 first generation 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 (Hobbs et al. 2005; Fryer & Kalogera 2001) , 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,

where ln Λ 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 cluster. The binary is formed according to a mass distribution given by ∝ (M 1 + M 2 ) 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

where a m = max(a ej , a 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 can be ejected. The ejection of the interloper black holes is also incorporated (see App. 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 τ GW = a m /|ȧ GW |. At this point in the binary's evolution, we modify the black holes' masses and spins according to Eqs. (4) and (3) to incorporate the effects of superradiance, ensuring each black hole's age exceeds to the growth timescale of the boson cloud.

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 (Campanelli et al. 2006 (Campanelli et al. , 2007a González et al. 2007b,a; Lousto et al. 2010; Lousto & Zlochower 2013; Varma et al. 2019; Merritt et al. 2004; Varma et al. 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 App. 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 −1 . Conversely, non-spinning binaries typically lead to the smallest recoil velocities. For example, the largest recoil velocity a nonspinning binary can produce is only ∼ 170 km s −1 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,

After calculating the properties of the merger remnant, we determine if or where it should be re-introduced 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) 

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 (cf. Eq. (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 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.

We create a grid of ≈ 1800 simulations in the space of cluster mass M cl , initial density ρ h,0 , and boson mass µ. We select log-uniformly spaced cluster masses between 10 6 and 2 × 10 8 M , initial densities between 10 3 and 10 7 M pc − 3, and boson masses between 5 × 10 −15 and 5 × 10 −13 eV. We run the simulations with and without ultralight bosons. 2

In Fig. 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 square) or ejected (red cross) binary systems. Different initial densities lead to different escape velocities and dynamical friction timescales. The denser cluster (left), with an initial density of ρ h,0 = 1.9 × 10 6 M pc −3 , has an initial escape velocity of v esc,0 = 391.8 km s −1 . Whereas, the less dense cluster (right) only has an initial escape velocity of v esc,0 = 224.2 km s −1 .

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 thirdgeneration (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 suffi-ciently 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 Fig. 2 .

In the left column of Fig. 2 , we show the binary black hole mergers which occurr in a cluster simulation with M cl = 1.1 × 10 7 M and ρ h,0 = 1.9 × 10 6 M pc −3 . We see that all simulations are able to support hierarchical black hole growth. However, the existence of 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 generation black hole mergers (cf. Fig. 1 ) when bosons of masses ∼ 2×10 −14 -2 × 10 −13 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 −15 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 spindown from superradiance in the high mass, hierarchical growth regime.

On the right of Fig. 2 , we show the results from a cluster with M cl = 1.1 × 10 7 M and ρ h,0 = 6.6 × 10 4 M pc −3 . Only one particular ultralight boson, with a mass of µ = 4.1 × 10 −14 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 −14 eV as our best-case scenario example.

In Fig. 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 Fig. 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 −14 eV, 2.2 × 10 −13 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 −14 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 lowspin population is entirely retained within the cluster.

In the right panel, only the µ = 4.1 × 10 −14 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 −13 eV and no ultralight boson, the black holes are not spun-down enough to reduce the velocity distribution significantly. Therefore the majority of second-generation 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 3 . 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 Fig. 4 for µ = 4.1 × 10 −14 eV (orange), as well as simulations where no ultralight boson is present (hatched blue).

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 Fig. 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 Fig. 4 . See App. C for details about the calculation and uncertainty estimation. As an example, for the contours presented in Fig. 4 , the fraction of present-day NSCs 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 −14 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 Fig. 5 . The fraction peaks at µ 4.1 × 10 −14 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 −14 -2 × 10 −13 eV exist.

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 log 10 (M cl /M ) ∼ N (µ = 6.3, σ = 0.8) and initial densities from log 10 (ρ h,0 /M pc −3 ) ∼ N (µ = 4.5, σ = 1.5). Each cluster is evolved to an age drawn at random from log 10 (t age /Gyr) ∼ N (µ = 0.31, σ = 1). These distributions are chosen such that the evolved population visually appears similar to the observed nuclear star cluster for all binaries forming a third generation (3G) black hole, two second generation black holes (2G+2G) and second generation and first generation (2G+1G) merger events. The black lines corresponds to the initial escape velocity of the cluster, vesc,0. The distributions correspond to six panels in Fig. 2 , when µ = 4.1 × 10 −14 , 2.2 × 10 −13 eV and no ultralight boson is present. Of these ultralight boson masses, only µ = 4.1 × 10 −14 eV (orange) facilitates hierarchical black hole growth in both example clusters considered. Therefore, black holes spinning down due to superradiance directly leads to a higher retention fraction of black holes and consequently a higher chance of hierarchical growth.

properties from Georgiev et al. (2016) . We evolve 4.8 × 10 4 clusters with properties drawn randomly from these distributions, for scenarios where µ = 4.1 × 10 −14 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 spindown limits set by µ = 4.1 × 10 −14 eV, are shown in Fig. 6 . The densities of the merger fraction, f m , in mass and spin are also presented.

From the one-dimensional marginal distributions, the inclusion of superradiance from bosons with µ = 4.1 × 10 −14 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.

The process of black-hole superradiance can increase the number of nuclear star clusters where hierarchical black hole growth occurs while impacting the black hole merger population. The signature of this distinct population may be detectable by gravitational-wave detectors such as Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2014; Abbott et al. 2020b,c,d) .

In Fig. 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. 2020d,b,c) . . Effective radius (R eff ) -cluster mass (M cl ) parameter space that can support hierarchical growth, and present-day populations for globular clusters (Baumgardt et al. 2019 ) and nuclear star clusters (Georgiev et al. 2016 ). The solid orange and hatched blue contours correspond to the regions of the parameter-space where more than 50% of cluster simulations at a given point can support hierarchical growth if µ = 4.1 × 10 −14 eV or if no ultralight boson exists, respectively. We assume the effective radius is approximately given by the half-mass radius as R eff 0.75r h . Any observed cluster within this contour is capable of undergoing hierarchical black hole growth now or in the future. The dashed contours indicate the evolved state of the bounded region after 10 Gyr. A cluster under the dashed contours might have supported hierarchical growth in the two scenarios.

By comparing the gravitational-wave detectionweighted 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 mergers populations, we compute the merger rates of both the total population, R, as well as the intermediate mass (100 < M/M < 1000) population; please refer to App. D for the details. We summarize the merger rates calculated and present the observed rates (Abbott et al. 2020c (Abbott et al. ,a, 2021 in Tab Fig. 4 . In the absence of ultralight bosons, we find only ≈ 4.5% of observed nuclear star clusters from (Georgiev et al. 2016) can support hierarchical growth (black interval shown), whereas bosons with masses of µ ∼ 2×10 −14 -2×10 −13 eV lead to an increased fraction (shown in orange).

in the absence of ultralight bosons are consistent with the results from Antonini et al. (2019) . 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; Neijssel et al. 2019; Eldridge et al. 2019; Mapelli et al. 2020; Santoliquido et al. 2020; Kremer et al. 2020b; Zevin et al. 2021; Wong et al. 2021; Rodriguez et al. 2021 ). However, the inclusion of superradiance leads to a doubling of the merger rate of intermediate mass black-hole mergers to ≈ 0.08f nMBH Gpc −3 yr −1 for µ = 4.1 × 10 −14 eV bosons. Both inferred values of R IMBH are consistent with the single-event merger rate determined for GW190521 (Abbott et al. 2020c (Abbott et al. , 2021 . 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 Table 1 . Total and intermediate mass black hole merger rates from observations of gravitational waves from stellar mass black hole mergers (Abbott et al. 2020c (Abbott et al. ,a, 2021 compared to the calculated merger rates from our simulated population (with 68% credible intervals). The observed intermediate mass black hole merger rate is a single-event rate determined from GW190521 (Abbott et al. 2020c (Abbott et al. , 2021 . 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 −14 eV leads to a doubling in the intermediate mass black hole merger rate for our simulated population. Here, fnMBH corresponds to the fraction of galaxies without a central massive black hole. GW190521 (Abbott et al. 2020b,c) , we conclude that the system is likely a 2G+2G binary black hole merger (discussed in Kimball et al. 2020; Romero-Shaw et al. 2020; Gayathri et al. 2020) .

Recently, Ng et al. (2021b) confirmed that the primary black holes from GW190412 (Abbott et al. 2020e) and GW190517 (Abbott et al. 2020d) exclude non-self-interacting ultralight scalar bosons with masses 1.3 -2.7 × 10 −13 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. (2020) . GW190412 is inconsistent with our simulated population irrespective of the presence of ultralight bosons -likely a direct result of our mass cut-off or spin distribution used. We can nevertheless interpret GW190412 as a direct result of hierarchical and/or dynamical formation (as in Abbott et al. 2020e; Safarzadeh & Hotokezaka 2020; Rodriguez et al. 2020; Zevin et al. 2020; Gerosa 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 −14 -2 × 10 −13 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 (Isi et al. 2019; Ng et al. 2021a,b; Tsukada et al. 2019; Sun et al. 2020; Brito et al. 2017b; Tsukada et al. 2021) . Other ultralight boson searches would likely cover a different boson mass range, which impact black holes of different masses. The detection of ultralight bosons in this energy range is likely only possible with future generations of gravitationalwave 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 −14 -10 −13 eV is 100 Mpc for aLIGO detectors (Aasi et al. 2015) . Third-generation gravitationalwave detectors such as 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. Analysing the mass-spin distributions of the stellar mass black holes may be a more successful (Aasi et al. 2015) for the µ = 4.1 × 10 −14 eV or no ultralight boson scenarios. The 1σ and 2σ credible intervals inferred for the primary black hole properties from GW190412 (Abbott et al. 2020e) , GW190517 (Abbott et al. 2020d) , and GW190521 (Abbott et al. 2020b,c) are shown as solid contours. The properties of GW190521's secondary black hole is given by the dashed contours. The existence of a µ = 4.1 × 10 −14 eV ultralight boson would lead to an enhancement of heavier binary systems with lower spins. Currently observations from gravitational waves cannot rule out the existence of ultralight bosons in the mass range 10 −14 -10 −13 eV. method for indirectly determining the existence of ultralight bosons (Ng et al. 2021a,b) . Ng et al. (2021a) found that 25 +95 −15 or 80 +210 −70 gravitational-wave detections with signal-to-noise ratios exceeding 30 were required to determine the presence of ultralight bosons with µ = 10 −13 eV, for two fiducial first-generation spin distributions. It is unlikely that these detection numbers can be applied to a ∼ 10 −14 eV boson detection prediction. While current ground-based observatories provide ultralight boson detection opportunities from nearby known black holes, future gravitational-wave detectors provide a significantly greater chance.

A.1. Cluster evolution 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 (Hénon , 1975 Gieles et al. 2011; Breen & Heggie 2013) . The rate of energy generation,Ė, is therefore directly related to the properties of the cluster,

where E −0.2GM 2 cl /r 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 )

Here, m * is the mean stellar and compact object mass within the cluster and ln Λ 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 ψ = m 2.5 * / m * 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, ρ h = 3M cl /8πr h 3 , and escape velocity, v esc ∝ M cl /r h , we can express the initial heating rate, relaxation time, and escape velocity at t = t 0 uniquely in terms of the initial cluster density and its mass. From Eqs. (A1) and (A2) (Antonini et al. 2019) , 

where M 5 ≡ M cl /10 5 M and ρ 5,0 ≡ ρ h,0 /10 5 M pc −3 . Here, t 0 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 0 = 7, where W 0 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 

following Eq. (A1). Fixing the cluster mass, we have τ rh ∝ r 3/2 h , and solve Eq. (A4) for the time evolution of the half-mass radius from t 0 ,

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

and v esc (t) = v esc,0

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 Eq. (A1),

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

where a is the semimajor axis of the binary's orbit. Therefore, differentiating both sides of Eq. (A9) and applying Eq. (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 

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

where we have used the fact that a h a m , as seen in Eq. (14). The lower bound semimajor axis, a m , is given by

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) 

where M 123 = M 1 +M 2 + M core is the average total mass of the binary-single interaction, and q 3 = M 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), M 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 

where, e is the eccentricity of the binary, and 

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 τ GW = a m /|ȧ 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 )

where the interloper mass ratio, q 3 , 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.

We utilize the Precession code (Gerosa & Kesden 2016) to determine the final recoil velocity of the remnant black hole. The kick velocity fitting formula are constructed from mass weighted spin vector combinations (Gerosa & Kesden 2016) ,

where χ 1,2 are the spin magnitudes of the individual black holes,Ŝ 1,2 are their associated unit vectors, and q = M 2 /M 1 is the mass ratio. These vectors can also be decomposed into the components parallel and orthogonal to the orbital angular momentum, L, given asχ =χ ·L,χ ⊥ = |χ ×L|, ∆ = ∆ ·L, and ∆ ⊥ = |∆ ×L|. The formulation of the final magnitude of the kick velocity of the black hole, v k , is constructed from the recoil velocities from mass, v m , and spin, v s and v s⊥ , asymmetries in the binary system. The velocity component from mass asymmetry lies perpendicular to the orbital angular momentum vector, whereas the spin component has some contribution both parallel to the vector, and in the orbital plane of the binary. This leads to the simple expression for the magnitude of the kick velocity (Campanelli et al. 2007a) ,

where β is the angle between the components from mass and spin asymmetry orthogonal to the orbital angular momentum. Each of the individual terms in Eq. (B20) has been determined through analytical fits to numerical relativity simulations. The velocity components v m , v s , and v s⊥ are calculated as v m = Aη 2 1 − q 1 + q (1 + Bη),

where η = M 1 M 2 /(M 1 + M 2 ) 2 , and the coefficients have been found through fitting to numerical relativity simulations: A = 1.2 × 10 4 km s −1 , B = −0.93 (González et al. 2007b ), H = 6.9 × 10 3 km s −1 (Lousto & Zlochower 2008) , V 11 = 3677.76 km s −1 , V A = 2481.21 km s −1 , V B = 1792.45 km s −1 , V C = 1506.52 km s −1 (Lousto et al. 2012) , C 2 = 1140 km s −1 , C 3 = 2481 km s −1 (Lousto & Zlochower 2013) , and β = 145 • (Lousto & Zlochower 2008) . Finally, Θ is the angle between ∆ ×L and the infall direction of the black holes at merger, with an additional offset of ∼ 200 • (Brügmann et al. 2008; Lousto & Zlochower 2009 ). This angle is assumed to be uniformly distributed from 0 to π (Gerosa & Kesden 2016) , and is therefore drawn at random for each merger calculation.

Galactic dynamics

Annales d'Astrophysique

Black Hole Spin: Theory and Observation

25th Rencontres de Blois on Particle Physics and Cosmology

Pis'ma Zh

We thank Richard Brito for insightful discussions about black hole superradiance. We also thank Susan Scott and Karl Wette for helpful discussions about ul-tralight boson cloud detection prospects. This work is supported through Australian Research Council (ARC) Centre of Excellence CE170100004. EP acknowledges the support of ARC CE170100004's COVID-19 support fund. PDL is supported through ARC Future Fellowship FT160100112, and ARC Discovery Project DP180103155. KK is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2001751. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. Computing was performed on LIGO Laboratory computing cluster at California Institute of Technology.APPENDIX

We use Monte Carlo integration to determine the fraction of the individual nuclear star cluster property distributions, p i (R eff , M cl ), contained within the contour (which itself has uncertainty), where i indexes the cluster. To estimate the uncertainty on the contours, the uncertainty on the fraction of simulations that demonstrate hierarchical growth must first be calculated. By assuming the fraction follows a binomial distribution, the uncertainty on f hier is estimated with a Wilson score interval (Wilson 1927) . The median and confidence interval are fit with a generalized extreme value (GEV) distribution following Possolo et al. (2019) . This allows us to draw 1000 possible fractions at each point in the parameter space, from which we create 1000 possible hierarchical growth contours. We fit the N NSC = 228 nuclear star cluster properties (Georgiev et al. 2016) with GEV distributions (Possolo et al. 2019 ) as an approximation for p i (R eff , M cl ) to generate N p = 4 × 10 3 plausible (R eff , M cl ) values for each cluster. We compute the integral aswhere Θ j (R i,k eff , M i,k cl ) evaluates to one for a point within the j th , f hier > 50% contour, and zero otherwise. The value of f NSC is taken as the median from the distribution f NSC,j , with an uncertainty set by the distribution. where dN gx / dV c ≈ 0.01 Mpc −3 (Conselice et al. 2005 ) is the number density of nucleated galaxies, and f nMBH is the fraction of galaxies without a central massive black hole. For each nuclear star cluster, we compute the reciprocal of the binary formation timescale, τ dyn , as an approximate for the rate of binary mergers within the cluster, since a single binary dominates the cluster evolution at any one time. Here, f nMBH is kept as a parameter in the expression as its value is somewhat unconstrained, though expected to be ∼ 0.2 -1 (Seth et al. 2008; Antonini et al. 2015; Nguyen et al. 2018) . Furthermore, the merger rate for which at least one member of the binary is within the black hole mass range 10 2 M < M < 10 3 M is calculated aswhere Θ(10 2 < M 1,i /M < 10 3 ) evaluates to one when the condition is satisfied and zero otherwise. This mass-range corresponds to the lower end of the intermediate mass black hole (IMBH) range (e.g., Greene et al. 2020 ). The merger remnant (and possibly also the primary black hole) for GW190521 lies conclusively within this range (Abbott et al. 2020b,c) . The calculation in Eq. D26 can be compared to the observed IMBH merger rates.