# Supermassive black holes from runaway mergers and accretion in nuclear star clusters

Konstantinos Kritos, Emanuele Berti, Joseph Silk

<figure><img src="image_1.png" /></figure>

## To cite this version:

Konstantinos Kritos, Emanuele Berti, Joseph Silk. Supermassive black holes from runaway mergers and accretion in nuclear star clusters. Mon.Not.Roy.Astron.Soc., 2024, 531 (1), pp.133-136. 10.1093/mnras/stae1145. hal-04576593

HAL Id: hal-04576593

https://hal.science/hal-04576593

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Supermassive black holes from runaway mergers and accretion in nuclear star clusters

Konstantinos Kritos,<sup>1*</sup> Emanuele Berti<sup>1</sup> and Joseph Silk<sup>1,2,3</sup>

Accepted 2024 April 26. Received 2024 April 22; in original form 2024 January 18

## ABSTRACT

Rapid formation of supermassive black holes occurs in dense nuclear star clusters that are initially gas-dominated. Stellar-mass black hole remnants of the most massive cluster stars sink into the core, where a massive runaway black hole forms as a consequence of combined effects of repeated mergers and Eddington-limited gas accretion. The associated gravitational wave signals of high-redshift extreme mass-ratio inspirals are a unique signature of the nuclear star cluster scenario.

**Key words:** galaxies: nuclei – quasars: supermassive black holes – black hole mergers.

## 1 INTRODUCTION

The emergence of supermassive black holes (SMBHs) in the early Universe is an outstanding problem in modern cosmology and astrophysics. Many active galactic nuclei (AGN) powered by SMBHs have been identified at redshift $z \gtrsim 6$, with spectroscopically inferred black hole (BH) masses in the range $10^6–10^8 M_\odot$ (Greene et al. 2024; Juodžbalis et al. 2024; Maiolino et al. 2024). As suggested by recent spectroscopic surveys with the JWST, the abundance of red-selected AGN is a hundred times that expected in this redshift range (Maiolino et al. 2023; Matthee et al. 2024). The consensus is that SMBHs have grown out of massive seed BHs (Volonteri 2012), and the problem is then reduced to explaining the formation of those intermediate-mass BH seeds (Greene, Strader & Ho 2020). One possible formation channel is the direct collapse of pristine gas clouds at $z > 10$ (Begelman, Volonteri & Rees 2006; Latif et al. 2022). Alternative seeding mechanisms involve repeated BH mergers in dense star clusters (Quinlan & Shapiro 1987, 1989; Davies, Miller & Bellovary 2011; Lupi et al. 2014; Antonini, Gieles & Gualandris 2019; Fragione & Silk 2020; Kroupa et al. 2020; Atallah et al. 2023; Kritos, Berti & Silk 2023) and feeding by debris of successive tidal disruption events (TDEs; Stone, Küpper & Ostriker 2017; Inayoshi, Visbal & Haiman 2020; Rizzuto et al. 2023). Moreover, several works have studied the rapid formation of supermassive stellar objects with masses in the range $10^3–10^5 M_\odot$ within metal-poor massive star clusters (Boekholt et al. 2018; Tagawa, Haiman & Kocsis 2020; Das et al. 2021; Schleicher et al. 2022; Reinoso et al. 2023) and nuclear star clusters (NSCs; Escala 2021; Vergara et al. 2023) through runaway collisions of stars and gas accretion. These supermassive stellar objects can then collapse into intermediate-mass BH seeds.

NSCs are extremely dense and massive stellar systems that have been observed to reside in the centres of most galaxies in the local

Universe (Neumayer, Seth & Böker 2020). In the nearby Universe, NSCs span the mass range from $10^5$ to $10^9 M_\odot$, with evidence for in situ star formation in the more massive systems (Fahrion et al. 2022). At high redshift, a number of lensed star clusters with effective radii as small as a parsec (pc) and with masses of order $10^6 M_\odot$ have been identified (Vanzella et al. 2023; Adamo et al. 2024). There is some theoretical support for the fragmentation of the central region of protogalaxies into millions of metal-poor stars, forming clusters with half-mass radii $\lesssim 0.5$ pc (Devecchi et al. 2010).

In this paper, we use a semi-analytic NSC model and demonstrate the possibility of rapidly forming SMBHs through the combined processes of repeated BH mergers and Eddington-limited gas accretion, assuming that a substantial amount of gas can be retained within the cluster for a few tens of millions of years (Myr).

## 2 NUCLEAR STAR CLUSTER MODEL

Our numerical model for the NSC implements a two-mass system consisting of $N_*$ stars, each of mass $m_*$, and $N_{BH}$ stellar BHs, each of mass $m_{BH}$ and non-spinning, with an initial mass $M_g$ of residual gas not turned into stars. We find that typically one of the centrally located and initially stellar mass BHs grows dramatically via gas, star, and BH accretion, and we denote its mass and spin by $M_{BH}$ and $\chi_{BH}$, respectively. The total cluster mass is therefore $M_{cl} = m_*N_* + m_{BH}N_{BH} + M_g + M_{BH}$. Since the growing BH is one of the stellar BHs, initially $M_{BH,0} = m_{BH}$ and $\chi_{BH,0} = 0$. This ‘oligarchic’ growth scenario can occur in dynamical environments (see e.g. Matsubayashi, Shinkai & Ebisuzaki 2004; Kovetz et al. 2018), where the formation of multiple massive BHs within the same system is significantly suppressed by their preferential accretion by the most massive object (Kritos et al. 2023).

The gas, composed primarily of ionized hydrogen with sound speed $c_s \simeq 10 \text{ km s}^{-1}$ (Krumholz & Matzner 2009), is removed exponentially over a characteristic gas expulsion time-scale $\tau_{ge}$ (Baumgardt & Kroupa 2007). The growing BH may accrete some of

<footnote n="1"><sup>1</sup>William H. Miller III Department of Physics and Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA</footnote>

<footnote n="2"><sup>2</sup>Insitut d'Astrophysique de Paris, Sorbonne Université, UMR 7095 CNRS and UPMC, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Department of Physics, University of Oxford, Oxford OX1 3RH, United Kingdom</footnote>

<footnote>*E-mail: kkritos1@jhu.edu</footnote>
# Supermassive black holes from runaway mergers and accretion in nuclear star clusters

Konstantinos Kritos,<sup>1*</sup> Emanuele Berti<sup>1</sup> and Joseph Silk<sup>1,2,3</sup>

Accepted 2024 April 26. Received 2024 April 22; in original form 2024 January 18

## ABSTRACT

Rapid formation of supermassive black holes occurs in dense nuclear star clusters that are initially gas-dominated. Stellar-mass black hole remnants of the most massive cluster stars sink into the core, where a massive runaway black hole forms as a consequence of combined effects of repeated mergers and Eddington-limited gas accretion. The associated gravitational wave signals of high-redshift extreme mass-ratio inspirals are a unique signature of the nuclear star cluster scenario.

**Key words:** galaxies: nuclei – quasars: supermassive black holes – black hole mergers.

## 1 INTRODUCTION

The emergence of supermassive black holes (SMBHs) in the early Universe is an outstanding problem in modern cosmology and astrophysics. Many active galactic nuclei (AGN) powered by SMBHs have been identified at redshift $z \gtrsim 6$, with spectroscopically inferred black hole (BH) masses in the range $10^6-10^8 M_{\odot}$ (Greene et al. 2024; Juodžbalis et al. 2024; Maiolino et al. 2024). As suggested by recent spectroscopic surveys with the JWST, the abundance of red-selected AGN is a hundred times that expected in this redshift range (Maiolino et al. 2023; Matthee et al. 2024). The consensus is that SMBHs have grown out of massive seed BHs (Volonteri 2012), and the problem is then reduced to explaining the formation of those intermediate-mass BH seeds (Greene, Strader & Ho 2020). One possible formation channel is the direct collapse of pristine gas clouds at $z > 10$ (Begelman, Volonteri & Rees 2006; Latif et al. 2022). Alternative seeding mechanisms involve repeated BH mergers in dense star clusters (Quinlan & Shapiro 1987, 1989; Davies, Miller & Bellovary 2011; Lupi et al. 2014; Antonini, Gieles & Gualandris 2019; Fragione & Silk 2020; Kroupa et al. 2020; Atallah et al. 2023; Kritos, Berti & Silk 2023) and feeding by debris of successive tidal disruption events (TDEs; Stone, Küpper & Ostriker 2017; Inayoshi, Visbal & Haiman 2020; Rizzuto et al. 2023). Moreover, several works have studied the rapid formation of supermassive stellar objects with masses in the range $10^3-10^5 M_{\odot}$ within metal-poor massive star clusters (Boekholt et al. 2018; Tagawa, Haiman & Kocsis 2020; Das et al. 2021; Schleicher et al. 2022; Reinoso et al. 2023) and nuclear star clusters (NSCs; Escala 2021; Vergara et al. 2023) through runaway collisions of stars and gas accretion. These supermassive stellar objects can then collapse into intermediate-mass BH seeds.

NSCs are extremely dense and massive stellar systems that have been observed to reside in the centres of most galaxies in the local

Universe (Neumayer, Seth & Böker 2020). In the nearby Universe, NSCs span the mass range from $10^5$ to $10^9 M_{\odot}$, with evidence for in situ star formation in the more massive systems (Fahrion et al. 2022). At high redshift, a number of lensed star clusters with effective radii as small as a parsec (pc) and with masses of order $10^6 M_{\odot}$ have been identified (Vanzella et al. 2023; Adamo et al. 2024). There is some theoretical support for the fragmentation of the central region of protogalaxies into millions of metal-poor stars, forming clusters with half-mass radii $\lesssim 0.5$ pc (Devecchi et al. 2010).

In this paper, we use a semi-analytic NSC model and demonstrate the possibility of rapidly forming SMBHs through the combined processes of repeated BH mergers and Eddington-limited gas accretion, assuming that a substantial amount of gas can be retained within the cluster for a few tens of millions of years (Myr).

## 2 NUCLEAR STAR CLUSTER MODEL

Our numerical model for the NSC implements a two-mass system consisting of $N_{\star}$ stars, each of mass $m_{\star}$, and $N_{BH}$ stellar BHs, each of mass $m_{BH}$ and non-spinning, with an initial mass $M_{\text{g}}$ of residual gas not turned into stars. We find that typically one of the centrally located and initially stellar mass BHs grows dramatically via gas, star, and BH accretion, and we denote its mass and spin by $M_{BH}$ and $\chi_{BH}$, respectively. The total cluster mass is therefore $M_{\text{cl}} = m_{\star}N_{\star} + m_{BH}N_{BH} + M_{\text{g}} + M_{BH}$. Since the growing BH is one of the stellar BHs, initially $M_{BH,0} = m_{BH}$ and $\chi_{BH,0} = 0$. This ‘oligarchic’ growth scenario can occur in dynamical environments (see e.g. Matsubayashi, Shinkai & Ebisuzaki 2004; Kovetz et al. 2018), where the formation of multiple massive BHs within the same system is significantly suppressed by their preferential accretion by the most massive object (Kritos et al. 2023).

The gas, composed primarily of ionized hydrogen with sound speed $c_s \simeq 10 \text{ km s}^{-1}$ (Krumholz & Matzner 2009), is removed exponentially over a characteristic gas expulsion time-scale $\tau_{\text{ge}}$ (Baumgardt & Kroupa 2007). The growing BH may accrete some of

<footnote n="1"><sup>1</sup>William H. Miller III Department of Physics and Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA</footnote>

<footnote n="2"><sup>2</sup>Insitut d'Astrophysique de Paris, Sorbonne Université, UMR 7095 CNRS and UPMC, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Department of Physics, University of Oxford, Oxford OX1 3RH, United Kingdom</footnote>

<footnote>*E-mail: kkritos1@jhu.edu</footnote>
this gas. We assume an Eddington-limited Bondi accretion rate on to the BH seed $M_{\text{acc}}$ with Salpeter time-scale of 50 Myr, and a time evolving Plummer gas density $\rho_g = 3M_g/[4\pi(r_h/1.3)^3]$, where $r_h$ is the half-mass radius. As a consequence of gradual gas removal, $r_h$ expands adiabatically (Hills 1980).

The half-mass relaxation time $\tau_{rh}$, modified for two-mass systems following equation (10) from Antonini & Gieles (2020), is the characteristic time over which a fraction $\zeta \simeq 0.1$ of the total cluster energy $E_{\text{cl}} = -\kappa GM_{\text{cl}}^2/(2r_h)$ (with $\kappa \simeq 0.4$, see Spitzer & Hart 1971) can be transferred throughout the cluster via two-body relaxation (Breen & Heggie 2013). The core is dominated by BHs, which partially decouple from stars due to the impossibility of energy equipartition ('Spitzer instability'). The core collapses on a time-scale $\tau_{cc}$ given by equation (5) from Portegies Zwart & McMillan (2002).

At this point, three-body hard binaries (3bb) can form and heat the system. The rate of energy production $\dot{E}_{\text{cl}} = -\zeta E_{\text{cl}}/\tau_{rh}$ depends on $r_h$, and is controlled by the efficiency of two-body relaxation ('Hénon's principle'). We compute the energy generated by a single 3bb $Q_{3bb}$, and use this condition to calculate the 3bb rate $R_{3bb} = \dot{E}_{\text{cl}}/Q_{3bb}$ (as in Heggie & Hut 2003, page 244). If all BHs but one have been ejected from the system, we assume that energy generation proceeds via tidal disruption of stars accreting on to the runaway BH. The TDE rate is $R_{\text{TDE}} = \dot{E}_{\text{cl}}/Q_{\text{TDE}}$, where $Q_{\text{TDE}} = GM_{\text{BH}}\bar{m}_{\star}/(2r_{\text{T}})$ is the heat released from a single TDE and $r_{\text{T}}$ is the tidal radius (Rees 1988). The main advantage of our method, as in Antonini & Gieles (2020), is that it is not necessary to calculate the properties of the core, while we can incorporate the effect of a cusp around the SMBH.

The differential equations that govern the evolution of the system are given by

$$ \frac{d\bar{m}_{\star}}{dt} = -\nu \frac{\bar{m}_{\star}}{t} \Theta(t - t_{\text{se}}), \qquad (1a) $$

$$ \frac{dN_{\star}}{dt} = -\xi_{e} \frac{N_{\star}}{\tau_{rh}} - R_{\text{TDE}}, \qquad (1b) $$

$$ \frac{dN_{\text{BH}}}{dt} = -k_{\text{ejs}} R_{3bb}, \qquad (1c) $$

$$ \frac{dM_{\text{g}}}{dt} = -\frac{M_{\text{g}}}{\tau_{\text{ge}}} - \dot{M}_{\text{acc}}, \qquad (1d) $$

$$ \frac{dM_{\text{BH}}}{dt} = \dot{M}_{\text{acc}} + f_{\text{TDE}} \bar{m}_{\star} R_{\text{TDE}}, \qquad (1e) $$

$$ \frac{dr_{\text{h}}}{dt} = \left( \frac{\zeta r_{\text{h}}}{\tau_{\text{rh}}} + \frac{2r_{\text{h}}}{M_{\text{cl}}} \frac{dM_{\text{cl}}}{dt} \right) \Theta(t - \tau_{\text{cc}}) - \frac{r_{\text{h}}}{M_{\text{cl}}} \left( N_{\star} \frac{d\bar{m}_{\star}}{dt} - \frac{M_{\text{g}}}{\tau_{\text{ge}}} \right). \qquad (1f) $$

Equations (1a) and (1b) can be compared with Alexander et al. (2014). The main difference is a loss term due to TDEs, $R_{\text{TDE}}$, on the right-hand side of equation (1b). The first three terms on the right-hand side of equation (1f) follow Antonini & Gieles (2020), but we also include an adiabatic expansion term due to gas expulsion. The quantity $k_{\text{ejs}}$ in equation (1c) is the number of single BHs ejected per unit binary. We set $\nu = 0.07$, $t_{\text{se}} = 2$ Myr (Antonini & Gieles 2020), and $\xi_{e} = 0.0074$ (Alexander et al. 2014).

At each iteration, we sample $k_{3bb}$, the number of 3bbs formed, from a Poisson distribution with parameter $\lambda = R_{3bb}dt$, where $dt$ is the time-step of the simulation. Every 3bb is formed with components $M_{\text{BH}}$ and $m_{\text{BH}}$, with its initial semimajor axis $a_0$ determined by sampling the initial hardness parameter $\eta_0 \ge 1$ (Morscher et al. 2015). Binaries harden until they reach their ejection ($a_{\text{ejb}}$) or gravitational wave (GW) radius ($a_{\text{GW}}$), assuming thermal eccentricity. To compute the

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Time evolution of the growing BH mass $M_{\text{BH}}$ for different initial amounts of gas in the system. The gas expulsion time-scale is set to $\tau_{\text{ge}} = 100$ Myr, the initial stellar number to $N_{\star,0} = 10^7$, and the initial half-mass radius to $r_{h,0} = 0.1$ pc. Thick (thin) lines correspond to simulations with $m_{\text{BH}} = 10 M_{\odot}$ ($m_{\text{BH}} = 50 M_{\odot}$). For comparison, we also show how a $10 M_{\odot}$ seed (thick grey) and a $50 M_{\odot}$ seed (thin grey) would grow under the sole effect of Eddington accretion.</figcaption></figure>

merger remnant mass $m_r$, spin $\chi_r$, and relativistic recoil – or GW kick $v_k$ – we implement the equations in Sec. V.A-C of Gerosa & Kesden (2016) and references therein, which are based on numerical relativity fits and assume isotropic spin directions. If $a_{\text{GW}} > a_{\text{ejb}}$ the binary merges within the cluster through the emission of GWs, and if $v_k < v_{\text{esc}}$, the merger remnant is retained in the cluster. In this case, we discretely update the runaway BH mass, spin, and number by the substitutions $M_{\text{BH}} \rightarrow m_r$, $\chi_{\text{BH}} \rightarrow \chi_r$, and $N_{\text{BH}} \rightarrow N_{\text{BH}} - 1$, while in any case of ejection, $M_{\text{BH}} \rightarrow m_{\text{BH}}$, $\chi_{\text{BH}} \rightarrow 0$, and $N_{\text{BH}} \rightarrow N_{\text{BH}} - 2$.

### 3 BLACK HOLE GROWTH

We implement a fourth-order Runge-Kutta scheme to numerically solve the system of equations (1) using a time-step of $dt = 0.01$ Myr (chosen to be much smaller than $\tau_{rh}$). We set $N_{\star,0} = 10^7$, $r_{h,0} = 0.1$ pc (where the '0' refers to initial values) and $\tau_{\text{ge}} = 100$ Myr. Using a Kroupa initial mass function in the range [0.08, 150] $M_{\odot}$, the average initial stellar mass is $\bar{m}_{\star,0} = 0.59 M_{\odot}$, while the fraction of BH progenitors with mass >20 $M_{\odot}$ is $\simeq 0.18$ per cent. We assume $m_{\text{BH}} = 10 M_{\odot}$, but to quantify uncertainties, we also repeat the analysis with $m_{\text{BH}} = 50 M_{\odot}$. In Fig. 1, we show the growth of $M_{\text{BH}}$ for four selected values of $M_{\text{g,0}} = 0, 10^6, 10^7, 10^8 M_{\odot}$.

All simulated systems have $v_{\text{esc,0}} > 300 \text{ km s}^{-1}$, allowing repeated mergers to form a runaway BH (Antonini et al. 2019). The BH rapidly grows by hundreds of repeated BH mergers (Fig. 2, left panel) with a merger rate that peaks at $\sim 100 \text{ Myr}^{-1}$ right after core collapse (see inset in the left panel of Fig. 2) and steadily drops as the system expands.

Significant amounts (20 per cent–30 per cent) of gas can be accreted within 200 Myr, rapidly bringing the BH mass into the SMBH regime. The residual gas does not fragment because of heating by UV radiation from cluster stars. Our accretion channel is similar to the direct collapse scenario but occurs on the Salpeter time-scale, rather than on the dynamical time-scale. The combination of mergers and accretion helps accelerate the growth because Eddington-limited accretion cannot itself grow a stellar-mass BH seed above $1000 M_{\odot}$ on a similar time-scale (compare the grey lines in Fig. 1, which ignore BH mergers).
this gas. We assume an Eddington-limited Bondi accretion rate on to the BH seed $M_{\text{acc}}$ with Salpeter time-scale of 50 Myr, and a time evolving Plummer gas density $\rho_g = 3M_g/[4\pi(r_h/1.3)^3]$, where $r_h$ is the half-mass radius. As a consequence of gradual gas removal, $r_h$ expands adiabatically (Hills 1980).

The half-mass relaxation time $\tau_{rh}$, modified for two-mass systems following equation (10) from Antonini & Gieles (2020), is the characteristic time over which a fraction $\zeta \simeq 0.1$ of the total cluster energy $E_{\text{cl}} = -\kappa GM_{\text{cl}}^2/(2r_h)$ (with $\kappa \simeq 0.4$, see Spitzer & Hart 1971) can be transferred throughout the cluster via two-body relaxation (Breen & Heggie 2013). The core is dominated by BHs, which partially decouple from stars due to the impossibility of energy equipartition ('Spitzer instability'). The core collapses on a time-scale $\tau_{cc}$ given by equation (5) from Portegies Zwart & McMillan (2002).

At this point, three-body hard binaries (3bb) can form and heat the system. The rate of energy production $\dot{E}_{\text{cl}} = -\zeta E_{\text{cl}}/\tau_{rh}$ depends on $r_h$, and is controlled by the efficiency of two-body relaxation ('Hénon's principle'). We compute the energy generated by a single 3bb $Q_{3bb}$, and use this condition to calculate the 3bb rate $R_{3bb} = \dot{E}_{\text{cl}}/Q_{3bb}$ (as in Heggie & Hut 2003, page 244). If all BHs but one have been ejected from the system, we assume that energy generation proceeds via tidal disruption of stars accreting on to the runaway BH. The TDE rate is $R_{\text{TDE}} = \dot{E}_{\text{cl}}/Q_{\text{TDE}}$, where $Q_{\text{TDE}} = GM_{\text{BH}}\bar{m}_{\star}/(2r_{\text{T}})$ is the heat released from a single TDE and $r_{\text{T}}$ is the tidal radius (Rees 1988). The main advantage of our method, as in Antonini & Gieles (2020), is that it is not necessary to calculate the properties of the core, while we can incorporate the effect of a cusp around the SMBH.

The differential equations that govern the evolution of the system are given by

$$ \frac{d\bar{m}_{\star}}{dt} = -\nu \frac{\bar{m}_{\star}}{t} \Theta(t - t_{\text{se}}), \quad (1a) $$

$$ \frac{dN_{\star}}{dt} = -\xi_{e} \frac{N_{\star}}{\tau_{rh}} - R_{\text{TDE}}, \quad (1b) $$

$$ \frac{dN_{\text{BH}}}{dt} = -k_{\text{ejs}} R_{3bb}, \quad (1c) $$

$$ \frac{dM_{\text{g}}}{dt} = -\frac{M_{\text{g}}}{\tau_{\text{ge}}} - \dot{M}_{\text{acc}}, \quad (1d) $$

$$ \frac{dM_{\text{BH}}}{dt} = \dot{M}_{\text{acc}} + f_{\text{TDE}} \bar{m}_{\star} R_{\text{TDE}}, \quad (1e) $$

$$ \frac{d\bar{r}_{\text{h}}}{dt} = \left( \frac{\zeta r_{\text{h}}}{\tau_{\text{rh}}} + \frac{2r_{\text{h}}}{M_{\text{cl}}} \frac{dM_{\text{cl}}}{dt} \right) \Theta(t - \tau_{\text{cc}}) - \frac{r_{\text{h}}}{M_{\text{cl}}} \left( N_{\star} \frac{d\bar{m}_{\star}}{dt} - \frac{M_{\text{g}}}{\tau_{\text{ge}}} \right). \quad (1f) $$

Equations (1a) and (1b) can be compared with Alexander et al. (2014). The main difference is a loss term due to TDEs, $R_{\text{TDE}}$, on the right-hand side of equation (1b). The first three terms on the right-hand side of equation (1f) follow Antonini & Gieles (2020), but we also include an adiabatic expansion term due to gas expulsion. The quantity $k_{\text{ejs}}$ in equation (1c) is the number of single BHs ejected per unit binary. We set $\nu = 0.07$, $t_{\text{se}} = 2$ Myr (Antonini & Gieles 2020), and $\xi_{e} = 0.0074$ (Alexander et al. 2014).

At each iteration, we sample $k_{3bb}$, the number of 3bbs formed, from a Poisson distribution with parameter $\lambda = R_{3bb}dt$, where $dt$ is the time-step of the simulation. Every 3bb is formed with components $M_{\text{BH}}$ and $m_{\text{BH}}$, with its initial semimajor axis $a_0$ determined by sampling the initial hardness parameter $\eta_0 \ge 1$ (Morscher et al. 2015). Binaries harden until they reach their ejection ($a_{\text{ejb}}$) or gravitational wave (GW) radius ($a_{\text{GW}}$), assuming thermal eccentricity. To compute the

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Time evolution of the growing BH mass $M_{\text{BH}}$ for different initial amounts of gas in the system. The gas expulsion time-scale is set to $\tau_{\text{ge}} = 100$ Myr, the initial stellar number to $N_{\star,0} = 10^7$, and the initial half-mass radius to $r_{h,0} = 0.1$ pc. Thick (thin) lines correspond to simulations with $m_{\text{BH}} = 10 M_{\odot}$ ($m_{\text{BH}} = 50 M_{\odot}$). For comparison, we also show how a $10 M_{\odot}$ seed (thick grey) and a $50 M_{\odot}$ seed (thin grey) would grow under the sole effect of Eddington accretion.</figcaption></figure>

merger remnant mass $m_r$, spin $\chi_r$, and relativistic recoil – or GW kick $v_k$ – we implement the equations in Sec. V.A-C of Gerosa & Kesden (2016) and references therein, which are based on numerical relativity fits and assume isotropic spin directions. If $a_{\text{GW}} > a_{\text{ejb}}$ the binary merges within the cluster through the emission of GWs, and if $v_k < v_{\text{esc}}$, the merger remnant is retained in the cluster. In this case, we discretely update the runaway BH mass, spin, and number by the substitutions $M_{\text{BH}} \rightarrow m_r$, $\chi_{\text{BH}} \rightarrow \chi_r$, and $N_{\text{BH}} \rightarrow N_{\text{BH}} - 1$, while in any case of ejection, $M_{\text{BH}} \rightarrow m_{\text{BH}}$, $\chi_{\text{BH}} \rightarrow 0$, and $N_{\text{BH}} \rightarrow N_{\text{BH}} - 2$.

### 3 BLACK HOLE GROWTH

We implement a fourth-order Runge-Kutta scheme to numerically solve the system of equations (1) using a time-step of $dt = 0.01$ Myr (chosen to be much smaller than $\tau_{rh}$). We set $N_{\star,0} = 10^7$, $r_{h,0} = 0.1$ pc (where the '0' refers to initial values) and $\tau_{\text{ge}} = 100$ Myr. Using a Kroupa initial mass function in the range [0.08, 150] $M_{\odot}$, the average initial stellar mass is $\bar{m}_{\star,0} = 0.59 M_{\odot}$, while the fraction of BH progenitors with mass $> 20 M_{\odot}$ is $\simeq 0.18$ per cent. We assume $m_{\text{BH}} = 10 M_{\odot}$, but to quantify uncertainties, we also repeat the analysis with $m_{\text{BH}} = 50 M_{\odot}$. In Fig. 1, we show the growth of $M_{\text{BH}}$ for four selected values of $M_{\text{g,0}} = 0, 10^6, 10^7, 10^8 M_{\odot}$.

All simulated systems have $v_{\text{esc,0}} > 300 \text{ km s}^{-1}$, allowing repeated mergers to form a runaway BH (Antonini et al. 2019). The BH rapidly grows by hundreds of repeated BH mergers (Fig. 2, left panel) with a merger rate that peaks at $\sim 100 \text{ Myr}^{-1}$ right after core collapse (see inset in the left panel of Fig. 2) and steadily drops as the system expands.

Significant amounts (20 per cent–30 per cent) of gas can be accreted within 200 Myr, rapidly bringing the BH mass into the SMBH regime. The residual gas does not fragment because of heating by UV radiation from cluster stars. Our accretion channel is similar to the direct collapse scenario but occurs on the Salpeter time-scale, rather than on the dynamical time-scale. The combination of mergers and accretion helps accelerate the growth because Eddington-limited accretion cannot itself grow a stellar-mass BH seed above $1000 M_{\odot}$ on a similar time-scale (compare the grey lines in Fig. 1, which ignore BH mergers).
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Time evolution of the cumulative number of BH mergers $N_{\text{me}}$ (left panel), merger rate $R_{\text{me}}$ (left panel, inset), half-mass radius $r_{\text{h}}$ (right panel), and cluster mass $M_{\text{cl}}$ (right panel, inset) under the same initial conditions used in Fig. 1.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Characteristic GW strain amplitude $h_c$ (computed as in Bonetti & Sesana 2020) as a function of detector-frame GW frequency for five extreme mass-ratio inspirals (EMRIs) at redshift $z = 10$ (coloured lines). We use cosmological parameters from Planck 18 (Aghanim et al. 2020). The initial orbital frequency and eccentricity are indicated in the legend; labels on the coloured curves refer to the mass of the primary. The bands indicate the $h_c$ range as the secondary mass is varied between $10 M_{\odot}$ (lower bound) and $50 M_{\odot}$ (upper bound). The black solid curves correspond to the noise amplitudes of present and planned GW observatories: aLIGO (https://dcc.ligo.org/LIGO-T1800044/public), aVIRGO (https://git.ligo.org/sensitivity-curves/observing-scenario-paper-2019-2020-update), ET (https://www.et-gw.eu/index.php/etsensitivities), CE (https://cosmicexplorer.org/sensitivity.html), DECIGO (Yagi & Seto 2011), LISA (Robson, Cornish & Liu 2019), and $\mu$Ares (Sesana et al. 2021). Finally, the shaded regions correspond to lower redshift ($z = 1–2$) binary sources from stellar-mass compact astrophysical binaries (CAB), EMRIs, and SMBH binaries.</figcaption></figure>

Starting from the initial radius $r_{\text{h,0}} = 0.1$ pc, the cluster inflates by a factor of 10–100 (depending on $M_{\text{g,0}}$), so that $r_{\text{h,0}}$ is on the parsec scale after a Hubble time (see Fig. 2, right panel). When $m_{\text{BH}} = 50 M_{\odot}$, a larger amount of energy is produced, causing $r_{\text{h}}$ to expand to above 100 pc (a feature that resembles the bulges of present day dwarf galaxies). The rapid mass-loss in the first few hundred Myr is due to gas expulsion (see the inset of the right panel of Fig. 2). After core collapse, the gas-richest system (dotted red line) undergoes a phase of contraction due to balanced evolution that starts at $\simeq 22$ Myr. Once the gas has been expelled, the energy production term in equation (1f) dominates, and the system starts to expand again at $t \simeq 200$ Myr.

We have demonstrated with an example that a combination of gas accretion and BH mergers is capable of growing stellar-mass BH

seeds into the SMBH regime within only 200 Myr. Most notably, if such a cluster forms at $z = 15$, then (depending on the initial gas content) an SMBH can assemble by $z \simeq 10$. These SMBH seeds can then further grow by subsequent gas accretion episodes during gas inflows into the NSC following a galaxy merger. This hierarchical scenario can be studied by combining our formalism with galaxy merger trees.

We recall that very massive NSCs are observed locally in many galaxies. Many of them are inferred to have gas because we observe young stars. We find that NSC evolution inevitably results in expansion by an order of magnitude or more (see the right panel of Fig. 2). In contrast, at formation, many NSCs were massive and compact. This supports our contention that NSCs are a promising environment for generating
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time evolution of the cumulative number of BH mergers $N_{\text{me}}$ (left panel), merger rate $R_{\text{me}}$ (left panel, inset), half-mass radius $r_{\text{h}}$ (right panel), and cluster mass $M_{\text{cl}}$ (right panel, inset) under the same initial conditions used in Fig. 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Characteristic GW strain amplitude $h_c$ (computed as in Bonetti & Sesana 2020) as a function of detector-frame GW frequency for five extreme mass-ratio inspirals (EMRIs) at redshift $z = 10$ (coloured lines). We use cosmological parameters from Planck 18 (Aghanim et al. 2020). The initial orbital frequency and eccentricity are indicated in the legend; labels on the coloured curves refer to the mass of the primary. The bands indicate the $h_c$ range as the secondary mass is varied between $10 M_{\odot}$ (lower bound) and $50 M_{\odot}$ (upper bound). The black solid curves correspond to the noise amplitudes of present and planned GW observatories: aLIGO (https://dcc.ligo.org/LIGO-T1800044/public), aVIRGO (https://git.ligo.org/sensitivity-curves/observing-scenario-paper-2019-2020-update), ET (https://www.et-gw.eu/index.php/etsensitivities), CE (https://cosmicexplorer.org/sensitivity.html), DECIGO (Yagi & Seto 2011), LISA (Robson, Cornish & Liu 2019), and $\mu$Ares (Sesana et al. 2021). Finally, the shaded regions correspond to lower redshift ($z = 1–2$) binary sources from stellar-mass compact astrophysical binaries (CAB), EMRIs, and SMBH binaries.</figcaption></figure>

Starting from the initial radius $r_{\text{h,0}} = 0.1$ pc, the cluster inflates by a factor of 10–100 (depending on $M_{\text{g,0}}$), so that $r_{\text{h,0}}$ is on the parsec scale after a Hubble time (see Fig. 2, right panel). When $m_{\text{BH}} = 50 M_{\odot}$, a larger amount of energy is produced, causing $r_{\text{h}}$ to expand to above 100 pc (a feature that resembles the bulges of present day dwarf galaxies). The rapid mass-loss in the first few hundred Myr is due to gas expulsion (see the inset of the right panel of Fig. 2). After core collapse, the gas-richest system (dotted red line) undergoes a phase of contraction due to balanced evolution that starts at $\simeq 22$ Myr. Once the gas has been expelled, the energy production term in equation (1f) dominates, and the system starts to expand again at $t \simeq 200$ Myr.

We have demonstrated with an example that a combination of gas accretion and BH mergers is capable of growing stellar-mass BH

seeds into the SMBH regime within only 200 Myr. Most notably, if such a cluster forms at $z = 15$, then (depending on the initial gas content) an SMBH can assemble by $z \simeq 10$. These SMBH seeds can then further grow by subsequent gas accretion episodes during gas inflows into the NSC following a galaxy merger. This hierarchical scenario can be studied by combining our formalism with galaxy merger trees.

We recall that very massive NSCs are observed locally in many galaxies. Many of them are inferred to have gas because we observe young stars. We find that NSC evolution inevitably results in expansion by an order of magnitude or more (see the right panel of Fig. 2). In contrast, at formation, many NSCs were massive and compact. This supports our contention that NSCs are a promising environment for generating
SMBHs via an enhanced merger rate density in the high-redshift Universe.

The ultimate test of our model would be the observation of individual GW signals from SMBH-BH mergers in NSC cores, which may be possible with future GW observatories (Fig. 3). The GW background produced by these sources is a distinctive signature of the NSC scenario that deserves further study. Our predictions motivate the development of the science case for next-generation GW observatories spanning the µHz frequency range between LISA and Pulsar Timing Arrays. Our proposed GW background should be taken into account when estimating astrophysical confusion noise sources that can affect the sensitivity of proposed µHz detectors, such as µAres. For comparison, the black-dotted line in Fig. 3 shows an estimate of the confusion noise from SMBH and white dwarf binaries in the µAres band (Sesana et al. 2021).

## ACKNOWLEDGEMENTS

We thank Fabio Antonini, Mark Ho-Yeuk Cheung, and Giada Caneva Santoro for discussions. This paper was supported by the Onassis Foundation - scholarship ID: F ZT 041-1/2023-2024. KK and EB are supported by NSF grant nos AST-2006538, PHY-2207502, PHY-090003, and PHY-20043, by NASA grant nos 20-LPS20-0011 and 21-ATP21-0010, by the John Templeton Foundation Grant 62840, by the Simons Foundation, and by the Italian Ministry of Foreign Affairs and International Cooperation grant no. PGR01167. This work was carried out at the Advanced Research Computing at Hopkins (ARCH) core facility (rockfish.jhu.edu), which is supported by the NSF grant no. OAC-1920103.

## DATA AVAILABILITY

The output of the simulations is available upon a reasonable request to the corresponding author. A PYTHON implementation of our model through the code NUCE (for 'NUclear Cluster Evolution') is publicly available on GITHUB at the URL https://github.com/Kkritos/Nuce.

## REFERENCES

* Adamo A. et al., 2024, preprint (arXiv:2401.03224)

* Aghanim N. et al., 2020, A&A, 641, A6

* Alexander P. E. R., Gieles M., Lamers H. J. G. L. M., Baumgardt H., 2014, MNRAS, 442, 1265

* Antonini F., Gieles M., 2020, MNRAS, 492, 2936

* Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008

* Atallah D., Trani A. A., Kremer K., Weatherford N. C., Fragione G., Spera M., Rasio F. A., 2023, MNRAS, 523, 4227

* Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589

* Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289

* Boekholt T. C. N., Schleicher D. R. G., Fellhauer M., Klessen R. S., Reinoso B., Stutz A. M., Haemmerlé L., 2018, MNRAS, 476, 366

* Bonetti M., Sesana A., 2020, Phys. Rev. D, 102, 103023

* Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779

* Das A., Schleicher D. R. G., Leigh N. W. C., Boekholt T. C. N., 2021, MNRAS, 503, 1051

* Davies M. B., Miller M. C., Bellovary J. M., 2011, ApJ, 740, L42

* Devecchi B., Volonteri M., Colpi M., Haardt F., 2010, MNRAS, 409, 1057

* Escala A., 2021, ApJ, 908, 57

* Fahrion K., Leaman R., Lyubenova M., van de Ven G., 2022, A&A, 658, A172

* Fragione G., Silk J., 2020, MNRAS, 498, 4591

* Gerosa D., Kesden M., 2016, Phys. Rev. D, 93, 124066

* Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257

* Greene J. E. et al., 2024, ApJ, 964, 39

* Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: a Multi-disciplinary Approach to Star Cluster Dynamics. Cambridge Univ. Press, Cambridge, p. 372

* Hills J. G., 1980, ApJ, 235, 986

* Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27

* Juodžbalis I. et al., 2024, preprint (arXiv:2403.03872)

* Kovetz E. D., Cholis I., Kamionkowski M., Silk J., 2018, Phys. Rev. D, 97, 123003

* Kritos K., Berti E., Silk J., 2023, Phys. Rev. D, 108, 083012

* Kroupa P., Subr L., Jerabkova T., Wang L., 2020, MNRAS, 498, 5652

* Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352

* Latif M. A., Whalen D. J., Khochfar S., Herrington N. P., Woods T. E., 2022, Nature, 607, 48 Nature, 607, 48

* Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MNRAS, 442, 3616

* Maiolino R. et al., 2023, preprint (arXiv:2306.01230)

* Maiolino R. et al., 2024, Nature, 627, 59

* Matsubayashi T., Shinkai H., Ebisuzaki T., 2004, ApJ, 614, 864

* Matthee J. et al., 2024, ApJ, 963, 129

* Morscher M., Pattabiraman B., Rodriguez C., Rasio F. A., Umbreit S., 2015, ApJ, 800, 9

* Neumayer N., Seth A., Böker T., 2020, A&AR, 28, 4

* Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576, 899

* Quinlan G. D., Shapiro S. L., 1987, ApJ, 321, 199

* Quinlan G. D., Shapiro S. L., 1989, ApJ, 343, 725

* Rees M. J., 1988, Nature, 333, 523

* Reinoso B., Klessen R. S., Schleicher D., Glover S. C. O., Solar P., 2023, MNRAS, 521, 3553

* Rizzuto F. P., Naab T., Rantala A., Johansson P. H., Ostriker J. P., Stone N. C., Liao S., Irodotou D., 2023, MNRAS, 521, 2930

* Robson T., Cornish N. J., Liu C., 2019, Class. Quantum Gravity, 36, 105011

* Schleicher D. R. G. et al., 2022, MNRAS, 512, 6192

* Sesana A. et al., 2021, Exp. Astron., 51, 1333

* Spitzer Lyman J., Hart M. H., 1971, ApJ, 164, 399

* Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS, 467, 4180

* Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 892, 36

* Vanzella E. et al., 2023, ApJ, 945, 53

* Vergara M. C., Escala A., Schleicher D. R. G., Reinoso B., 2023, MNRAS, 522, 4224

* Volonteri M., 2012, Science, 337, 544

* Yagi K., Seto N., 2011, Phys. Rev. D, 83, 044011

This paper has been typeset from a TEX/LATEX file prepared by the author.
SMBHs via an enhanced merger rate density in the high-redshift Universe.

The ultimate test of our model would be the observation of individual GW signals from SMBH-BH mergers in NSC cores, which may be possible with future GW observatories (Fig. 3). The GW background produced by these sources is a distinctive signature of the NSC scenario that deserves further study. Our predictions motivate the development of the science case for next-generation GW observatories spanning the µHz frequency range between LISA and Pulsar Timing Arrays. Our proposed GW background should be taken into account when estimating astrophysical confusion noise sources that can affect the sensitivity of proposed µHz detectors, such as µAres. For comparison, the black-dotted line in Fig. 3 shows an estimate of the confusion noise from SMBH and white dwarf binaries in the µAres band (Sesana et al. 2021).

## ACKNOWLEDGEMENTS

We thank Fabio Antonini, Mark Ho-Yeuk Cheung, and Giada Caneva Santoro for discussions. This paper was supported by the Onassis Foundation - scholarship ID: F ZT 041-1/2023-2024. KK and EB are supported by NSF grant nos AST-2006538, PHY-2207502, PHY-090003, and PHY-20043, by NASA grant nos 20-LPS20-0011 and 21-ATP21-0010, by the John Templeton Foundation Grant 62840, by the Simons Foundation, and by the Italian Ministry of Foreign Affairs and International Cooperation grant no. PGR01167. This work was carried out at the Advanced Research Computing at Hopkins (ARCH) core facility (rockfish.jhu.edu), which is supported by the NSF grant no. OAC-1920103.

## DATA AVAILABILITY

The output of the simulations is available upon a reasonable request to the corresponding author. A PYTHON implementation of our model through the code NUCE (for 'NUclear Cluster Evolution') is publicly available on GITHUB at the URL https://github.com/Kkritos/Nuce.

## REFERENCES

* Adamo A. et al., 2024, preprint (arXiv:2401.03224)

* Aghanim N. et al., 2020, A&A, 641, A6

* Alexander P. E. R., Gieles M., Lamers H. J. G. L. M., Baumgardt H., 2014, MNRAS, 442, 1265

* Antonini F., Gieles M., 2020, MNRAS, 492, 2936

* Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008

* Atallah D., Trani A. A., Kremer K., Weatherford N. C., Fragione G., Spera M., Rasio F. A., 2023, MNRAS, 523, 4227

* Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589

* Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289

* Boekholt T. C. N., Schleicher D. R. G., Fellhauer M., Klessen R. S., Reinoso B., Stutz A. M., Haemmerlé L., 2018, MNRAS, 476, 366

* Bonetti M., Sesana A., 2020, Phys. Rev. D, 102, 103023

* Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779

* Das A., Schleicher D. R. G., Leigh N. W. C., Boekholt T. C. N., 2021, MNRAS, 503, 1051

* Davies M. B., Miller M. C., Bellovary J. M., 2011, ApJ, 740, L42

* Devecchi B., Volonteri M., Colpi M., Haardt F., 2010, MNRAS, 409, 1057

* Escala A., 2021, ApJ, 908, 57

* Fahrion K., Leaman R., Lyubenova M., van de Ven G., 2022, A&A, 658, A172

* Fragione G., Silk J., 2020, MNRAS, 498, 4591

* Gerosa D., Kesden M., 2016, Phys. Rev. D, 93, 124066

* Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257

* Greene J. E. et al., 2024, ApJ, 964, 39

* Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: a Multi-disciplinary Approach to Star Cluster Dynamics. Cambridge Univ. Press, Cambridge, p. 372

* Hills J. G., 1980, ApJ, 235, 986

* Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27

* Juodžbalis I. et al., 2024, preprint (arXiv:2403.03872)

* Kovetz E. D., Cholis I., Kamionkowski M., Silk J., 2018, Phys. Rev. D, 97, 123003

* Kritos K., Berti E., Silk J., 2023, Phys. Rev. D, 108, 083012

* Kroupa P., Subr L., Jerabkova T., Wang L., 2020, MNRAS, 498, 5652

* Krumholz M. R., Matzner C. D., 2009, ApJ, 703, 1352

* Latif M. A., Whalen D. J., Khochfar S., Herrington N. P., Woods T. E., 2022, Nature, 607, 48 Nature, 607, 48

* Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MNRAS, 442, 3616

* Maiolino R. et al., 2023, preprint (arXiv:2306.01230)

* Maiolino R. et al., 2024, Nature, 627, 59

* Matsubayashi T., Shinkai H., Ebisuzaki T., 2004, ApJ, 614, 864

* Matthee J. et al., 2024, ApJ, 963, 129

* Morscher M., Pattabiraman B., Rodriguez C., Rasio F. A., Umbreit S., 2015, ApJ, 800, 9

* Neumayer N., Seth A., Böker T., 2020, A&AR, 28, 4

* Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576, 899

* Quinlan G. D., Shapiro S. L., 1987, ApJ, 321, 199

* Quinlan G. D., Shapiro S. L., 1989, ApJ, 343, 725

* Rees M. J., 1988, Nature, 333, 523

* Reinoso B., Klessen R. S., Schleicher D., Glover S. C. O., Solar P., 2023, MNRAS, 521, 3553

* Rizzuto F. P., Naab T., Rantala A., Johansson P. H., Ostriker J. P., Stone N. C., Liao S., Irodotou D., 2023, MNRAS, 521, 2930

* Robson T., Cornish N. J., Liu C., 2019, Class. Quantum Gravity, 36, 105011

* Schleicher D. R. G. et al., 2022, MNRAS, 512, 6192

* Sesana A. et al., 2021, Exp. Astron., 51, 1333

* Spitzer Lyman J., Hart M. H., 1971, ApJ, 164, 399

* Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS, 467, 4180

* Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 892, 36

* Vanzella E. et al., 2023, ApJ, 945, 53

* Vergara M. C., Escala A., Schleicher D. R. G., Reinoso B., 2023, MNRAS, 522, 4224

* Volonteri M., 2012, Science, 337, 544

* Yagi K., Seto N., 2011, Phys. Rev. D, 83, 044011

This paper has been typeset from a TEX/LATEX file prepared by the author.