Formation and Evolution of Compact Binaries Containing Intermediate Mass Black Holes in Dense Star Clusters

This section briefly introduces how the code handles close encounters and then describes how we implement GW mergers and TDs for IMBHs in our simulations. The IMBH embedded in a dense star cluster experiences multiple interactions with stars and BHs. A close encounter between the two particles with masses $M_{1}$ and $M_{2}$ at a distance $R$ is treated by the KS regularization scheme,¹¹1KS regularization is a treatment introduced to reduce the computational load during close encounters between particles. This should be distinguished from the merger criteria of the binaries described in the following paragraphs. provided the following conditions are met:

• •

  the two particles are closer than $r_{\rm min}$ (RMIN), and

• •

  their time-step becomes less than $\Delta t_{\rm min}$ (DTMIN).

Both $r_{\rm min}$ and $\Delta t_{\rm min}$ are user-controlled parameters in N-body units (See Appendix D). We consistently choose $r_{\rm min}=10^{-4}r_{N}$ and $\Delta t_{\rm min}=10^{-5}t_{N}$ across all simulations presented here, suppressing an auto-adjustment during the run. Here, $r_{N}\approx 1.3R_{\rm hm}$ and $t_{N}\approx 15\,{\rm Myr}(r_{N}/M_{\rm Tot})^{1/2}$ are N-body length and time unit, respectively, where $R_{\rm hm}$ is the half-mass radius of the star cluster, and $M_{\rm Tot}$ is the initial total mass in our simulations.

Our GW merger and TD models operate by exmaning the orbital properties of all objects within $10^{-4}\,{\rm pc}$ when the IMBH is in a “single” (i.e., not undergoing any close encounter that meet the criteria for KS regularization), or by analyzing the orbital properties of the companion object when the IMBH is in a KS regularized state. When the IMBH-BH binary undergoes a close encounter with other particles, their evolution is treated in the routine such as triple, quad, or chain. In this case, however, we do not apply our model, as such events are rare. See Appendix B for more justification.

If the IMBH is in a “single” state, we compute the semi-major axis $a$ and eccentricity $e$ of particles within a distance of $10^{-4}\,{\rm pc}$ at every time-step in intgrt.f. If the IMBH forms a KS pair, we check the values of $a$ and $e$ for the binary during every call of the routine ksint.f. We then decide whether to merge the particles based on the corresponding criteria regardless of KS regularization. For a SBH, we calculate the merger timescale using the equation from Peters (1964) as

where $M_{\rm IMBH}$ and $M_{\rm BH}$ refer to the masses of the IMBH and the SBH, respectively. Then the binary merges instantly if $t_{\rm GW}<1.0\,{\rm Myr}$ during the run.²²2See Appendix A for more justification This timescale is sufficiently short to assume that the influence of the third particle on the binary is negligible. Therefore, we can safely consider that the binary evolves only through GW radiation. After the merger, we increase the mass of the IMBH by $M_{\rm BH}$.

Although our primary goal is to investigate IMBH-BH collisions, implementing TD is necessary to properly consider and remove stars that are very close to the IMBH. In the case of an IMBH-star encounter, we compute the pericenter of the orbit as $r_{\rm p}=a\left(1-e\right)$. Then, we remove the star from the simulation if $r_{\rm p}$ is less than the tidal radius, $r_{\rm T}$. The tidal radius is computed by $r_{\rm T}=r_{\star}\left(M_{\rm IMBH}/M_{\star}\right)^{1/3}$, where the stellar radius computed by $r_{\star}=R_{\odot}\left(M_{\star}/M_{\odot}\right)^{1/3}$ for a star with mass $M_{\star}$. In this case, the mass of the IMBH is increased by half of the removed star’s mass, $0.5\,M_{\star}$. The choice of accreting a half of the stellar mass is based on analytical estimates (e.g., Rees, 1988; Strubbe & Quataert, 2009). The position and velocity of the IMBH are adjusted after a merger or a TD event in such a way that momentum conservation is guaranteed.

The binary of BHs with masses $M_{1}$ and $M_{2}$ shrinks and circularizes by emitting energy as gravitational waves. If external perturbations are negligible, we can analytically predict their evolution based on the 2.5PN (2.5 post-Newtonian) approximation. We compute the evolution of the IMRIs based on the orbit-averaged prescription provided by Peters (1964). Figure 1 depicts the evolution of the IMBH-BH binary with $a=10\,{\rm au}$ with three different eccentricities, as calculated from Eqs. (F1) and (F2). An IMRI with relatively modest eccentricity, $e=0.9$, requires $\sim 10^{12}$ years to merge via pure GW radiation. In contrast, it takes less than $10^{5}\,{\rm yr}$ for a highly eccentric IMRI until the merger. The merger timescale of the binary is dramatically shorter for a more eccentric binary.

We then compute the sky-averaged signal-to-noise ratio (SNR) produced by the GW source at a luminosity distance $D_{\mathrm{L}}$ by combining the results from Finn & Thorne (2000) and Barack & Cutler (2004) as

where $f_{n}=nf=n\sqrt{G\left(M_{1}+M_{2}\right)/a^{3}}$ is the frequencies of the $n$-th harmonic, $h_{c,n}=\left(\pi D_{\rm L}\right)^{-1}\sqrt{2\dot{E}_{n}/\dot{f}_{n}}$ is the characteristic amplitude of the $n$-th harmonic at distance $D_{\rm L}$, $\dot{E}_{n}$ is the energy radiated through the $n$-th harmonic, and $S_{h}(f)$ represents the noise spectral density of the detectors at frequency $f$.³³3The prefactor 5 emerges as a result of averaging the noise over all-source directions (Finn & Thorne, 2000). We explore the signals onto five detectors; aLIGO, ET, LISA, aSOGRO, and DECIGO. See Appendix H for more information. The integration domain of $f$ starts 4 years before the binary reaches the innermost stable circular orbit (ISCO).

Figure 2 shows the evolution of the first five characteristic amplitudes, $h_{c,n}$ of an IMRI. The IMRI considered consists of an IMBH with $M_{\rm IMBH}=10^{3}\,{\rm M}_{\odot}$ and a BH with $M_{\rm BH}=40\,{\rm M}_{\odot}$. The figure includes the curve of integrated GW signals and the intrinsic sensitivity curve for GW observatories. The GW frequency increases as the binary shrinks, and it remains in the LISA band for a period of time. The inspiral-plunge transition occurs approximately when the binary’s orbital frequency crosses the ISCO frequency. The frequency of its 2nd harmonic is given by

where $1/(1+z)$ reflects the effect of the redshift due to the expansion of the universe.

Due to the limited computational resources, our star clusters should have smaller particle numbers and sizes than the realistic star clusters. Therefore, we try to mimic the GCs by analyzing the IMRI rates as a function of velocity dispersion given by the following form:

to establish an empirical relation estimating the IMRI rates (see Eq.(13) and Section 4.2). Here, $\sigma_{\rm hm}$ is the 1D half-mass velocity dispersion of the cluster.

Our simulations also provide statistics on the orbital properties of IMRIs, such as the eccentricities or the masses of secondary BHs. We then estimate the detection rate of these events using the formula:

which contains the following components:

• •

  $p_{\rm IMBH}$ is the IMBH retention probability of GCs. The magnitude of the natal kick can be modeled as a Maxwellian distribution (Hobbs et al., 2005), or it can be derived from the mass fallback caused by the supernova explosion during the IMBH synthesis (Belczynski et al., 2002). About $\lesssim 30\%$ of GCs are expected to host IMBHs, and this probability may fall close to $0\%$ for GCs with masses below $10^{5}\,{\rm M}_{\odot}$ (Giersz et al., 2015; Askar et al., 2017, 2022). We utilize the retention probability of IMBH formation depending on the star cluster mass considering both high and low natal kicks referring Table 1 in Askar et al. (2022).

• •

  $f_{\rm det}$ is the fraction of IMRIs observable by a GW detector, which depends on $M_{\rm IMBH}$, redshift $z$, and the specific detector design. We assume that detectable IMRIs must have a signal-to-noise ratio $S/N>8$. This factor is modeled based on the distribution of orbital properties of IMRI events derived from our simulations.

• •

  $\mathrm{d}V_{c}/\mathrm{d}z$ is the comoving volume element at a given redshift bin, and $1/(1+z)$ accounts for the time dilation.

• •

  $n_{\rm GC}$ is the number density of GCs in each redshift bin. Figure 3 presents the “normalized” GC number density, $\nu_{\rm GC}\equiv n_{\rm GC}/n_{{\rm GC},\,z=0}$, at each redshift bin, while the absolute value of $n_{\rm GC}$ depends on the estimation of local GC density. Indeed, the GC number density at $z=0$ varies across studies — from $0.77\,{\rm Mpc}^{-3}$ (Rodriguez et al., 2015) to $8.4\,h^{3}\,{\rm Mpc}^{-3}$ (Portegies Zwart & McMillan, 1999). In our approximation, assuming GCs with masses $5\times 10^{4}$ to $8\times 10^{6}\,{\rm M}_{\odot}$ account for $0.1\%$ of the galaxy mass, we obtain $n_{{\rm GC},\,z=0}=1.45\,{\rm Mpc}^{-3}$. Note that the average mass of these GCs becomes $\left<M_{\rm GC}\right>=1.9\times 10^{5}\,{\rm M}_{\odot}$ here. We also model the halo mass function, $\phi(\bar{M},z)$, by the Schechter (1976) function (see Appendix E for detailed information), with the parameters adopted from Table.(1) of Conselice et al. (2016). More details about $n_{\rm GC}$ can be found in Appendix G.

• •

  $\mathrm{d}\nu_{\rm GC}/\mathrm{d}M_{\rm GC}$ refers to the normalized GC mass function. We assume that it follows a simple power-law, $\mathrm{d}\nu_{\rm GC}/\mathrm{d}M_{\rm GC}\propto M_{\rm GC}^{-s}$ with a slope parameter $s=2.2$ (Gieles, 2009).

• •

  $\mathrm{d}M_{\rm GC}/\mathrm{d}M_{\rm IMBH}$ models the relation between $M_{\rm IMBH}$ and $M_{\rm GC}$. To simplify our calculation, here we use the linear relation from Arca-Sedda (2016),

  $$\log{\left(\frac{M_{\rm IMBH}}{\,{\rm M}_{\odot}}\right)}=\log{\left(\frac{M_{\rm GC}}{\,{\rm M}_{\odot}}\right)-2.23}.$$

  (6)

  According to this relation, IMBHs of $300-5000\,{\rm M}_{\odot}$ correspond to the GCs of masses $5\times 10^{4}-10^{6}\,{\rm M}_{\odot}$.

We generate the initial conditions of the star cluster using the code Mcluster (Küpper et al., 2011). The simulated star clusters are designed to follow the Plummer profile,

where $M_{\rm SC}$ and $a_{\rm P}$ refer to the total mass and scale size of the star cluster, respectively. The half-mass radius $R_{\rm hm}$ of the Plummer cluster is approximated by $1.3a_{\rm P}$. We initialize the clusters without including primordial mass segregation or binaries.

To construct a realistic SBH mass function, we generate the stellar population by sampling from the Kroupa (2001) initial mass function (IMF) within stellar mass range of $0.08-150\,{\rm M}_{\odot}$. We then evolve the stars using the stellar evolution code SEVN (Spera et al., 2019; Iorio et al., 2023) over $100\,{\rm Myr}$. Figure 4 presents the mass function before and after stellar evolution. It reveals that the neutron stars (NSs) and SBHs are spread over the mass range from $0.8\,{\rm M}_{\odot}$ to $60\,{\rm M}_{\odot}$. Most stars with initial masses above $10\,{\rm M}_{\odot}$ turn into the BHs, with their final mass decreasing by the effect of the pair-instability supernova (PISN) or pulsational pair-instability supernove (PPISN) (e.g., Woosley, 2017; Woosley & Heger, 2021). Meanwhile, the NSs dominate the lower-mass end of the compact object spectrum. Since our simulations run over a relatively short timescale up to $30\,{\rm Myr}$, we suppress further stellar evolution during our experiments. This is because stellar evolution processes extend over significantly longer timescales and are beyond the primary focus of our study.

The typical mass range of GCs is a few$\times 10^{3}$ to $5\times 10^{6}\,{\rm M}_{\odot}$ with a scale radius of $\sim 3\,{\rm pc}$. However, higher-mass GCs still need $\gtrsim\mathcal{O}\left(10^{5}\right)$ particles to be simulated. We therefore mimic the dense environment of GCs by reducing the size of GCs instead of increasing their mass. This approach allows us to maintain a manageable number of particles resolved in simulations, thus reducing computational costs. We note that the overall evolution of the star cluster is significantly influenced by its total mass. However, the effect of the different GC masses on the GC evolution can be marginalized, as we terminate our simulations after only $t_{\rm term}=30\,{\rm Myr}$ (more discussion in Section 4.1).

We explore the IMRI rates in various stellar system settings, by changing its initial half-mass velocity dispersion from $\sigma_{\rm hm}\sim 6\,\mathrm{km}\,\mathrm{s}^{-1}$ to $\sim 10\,\mathrm{km}\,\mathrm{s}^{-1}$ (Group A in Table 1), and by varying the initial IMBH mass from $M_{{\rm IMBH},i}=300\,{\rm M}_{\odot}$ to $5000\,{\rm M}_{\odot}$ (Group B in Table 1). Before inserting the IMBH, we gradually modify the GC centers by slowly increasing the point mass potential, in order to avoid the sudden introduction of a massive object (as described in Section 3.2). All our simulated star clusters are initialized in isolation, with no external tidal fields included.

Both observational (e.g., Schödel et al., 2009) and theoretical (e.g., Bahcall & Wolf, 1976) studies agree that star clusters with a massive central object exhibit steeply increasing profiles for density and velocity dispersion toward the center. Sigurdsson et al. (1995) suggested that an “adiabatically growing” density profile can be a reasonable way to realize an $N$-body system with an IMBH. First, before inserting the IMBH as a live particle, we generate the initial condition of star clusters following a Plummer density profile using Mcluster. Then we introduce the Plummer potential with a gradually increasing mass $M_{\rm IMBH}\left(t\right)$ as

where $\epsilon_{\rm IMBH}\sim 10^{-4}\,{\rm pc}$ is the softening parameter placed to avoid the effect of the singularity (Hong & Lee, 2015). Then, $t_{\rm IMBH}$ after the initialization, we replace the external potential with a live, moving IMBH particle. Here, the time-dependent IMBH mass is designed to grow over time as

where $M_{{\rm IMBH},i}$ is the initial mass of the IMBH when it is finally inserted as a live particle, and $t_{\rm IMBH}$ is the mass growth timescale. Any choice of $t_{\rm IMBH}$ is suitable as long as it is not too short (Sigurdsson et al., 1995), thus we set it to be $t_{\rm IMBH}=50t_{\rm N}$. This is a relatively short timescale in our simulations, $\sim$ a few Myrs. The definition of the N-body units can be found in Appendix D. Figure 5 shows an example of “adiabatic growth” of a star cluster to introduce an embedded IMBH. At $t=0$, the cluster core has a flat density profile. As the analytic IMBH potential gradually grows, the inner profile evolves into a cuspy profile of $\rho\propto r^{-\gamma}$ with a slope index $\gamma\approx 1.4$. This is slightly shallower than the classical prediction by Bahcall & Wolf (1976), $\gamma=1.75$. This difference can be attributed to the nature of the star clusters consisting of stars with different masses (Baumgardt et al., 2004a).

We begin by examining the temporal evolution of the star clusters. The characteristic timescale for the evolution can be estimated by the relaxation timescale, $t_{\rm rlx}$, which is defined as

where $\sigma\left(r\right)$ is the local velocity dispersion, $\left<M_{\star}\right>$ is the average stellar mass, $\rho\left(r\right)$ is the mass density, and $\ln\lambda\approx 10$ is the Coulomb logarithm. This expression suggests that denser star clusters evolve more rapidly compared to less dense clusters.

The evolution of the half-mass radius, $R_{\rm hm}$, and the Lagrangian radius enclosing 5% of the initial cluster mass, $R_{5\%}$ of runs in Group B is depicted in Figure 6. These radii serve as key indicators of the structural evolution of star clusters over time. In general, the star clusters undergo expansion due to the relaxation processes. The rate of this expansion notably depends on the cluster’s relaxation timescale. For example, the GCs with initial properties $M_{\rm GC}=5\times 10^{4}\,{\rm M}_{\odot}$ and $R_{\rm hm}=0.7\,{\rm pc}$ undergoes an expansion of $R_{\rm hm}$ from $0.7\,{\rm pc}$ to $1.3\,{\rm pc}$. In contrast, the GCs with initial properties $M_{\rm GC}=10^{5}\,{\rm M}_{\odot}$ and $R_{\rm hm}=0.7\,{\rm pc}$ expands from $R_{\rm hm}=0.7\,{\rm pc}$ to only $R_{\rm hm}=1.0$, reflecting the relaxation timescale determined by the GC mass.

The evolution of the simulations in Group B is illustrated in Figure 7. All simulations are initialized with the same $R_{\rm 5\%}$ and $R_{\rm hm}$ at $t=0$ except for the small discrepancies caused by the adiabatic insertion of IMBH before initialization (See 3.2). As the GC evolves, the rate of expansion of the cluster size varies across different $M_{\rm IMBH,i}$. For a lower-mass IMBH ($M_{\rm IMBH,i}=300\,{\rm M}_{\odot}$), the half-mass radius increases from $R_{\rm hm}=0.7\,{\rm pc}$ to $R_{\rm hm}=1.1\,{\rm pc}$. For a more massive IMBH ($M_{\rm IMBH,i}=5000\,{\rm M}_{\odot}$), the half-mass radius grows significantly more, reaching $R_{\rm hm}=2.1\,{\rm pc}$. This trend highlights the crucial role of IMBH mass in driving the dynamical expansion of the host GC as discussed in previous studies (e.g., Baumgardt et al., 2004b; Freitag et al., 2006a; Hong & Lee, 2015).

The detailed evolution of the cluster core is depicted in Figure 8. The two left panels display the evolution of the central density, $\rho_{\rm c}$, for the populations within $r<R_{5\%}$ from the cluster’s center. In all simulations, the central density decreases monotonically as the core expands over time. The evolution of the central density is significantly affected by both the initial star cluster properties and the mass of the central IMBHs. This trend is already expected from the evolution of $R_{5\%}$ (See Figure 6 & 7).

The right panel of Figure 8 presents the evolution of the central velocity dispersion, $\sigma_{\rm c}$. Unlike the central density, the overall evolution of $\sigma_{\rm c}$ exhibits a more gradual and steady decline over time. In Group A, the evolution of $\sigma_{\rm c}$ closely follows the trend of central density decline, supporting the interpretation that the decrease in velocity dispersion is primarily driven by the expansion of the cluster core. As in the central density, the velocity dispersion of the clusters with longer $t_{\rm rlx}$ experiences more gradual dynamical evolution. For example, the GC with a total mass $M_{\rm GC}=5\times 10^{4}\,{\rm M}_{\odot}$ and an initial half-mass radius of $R_{\rm hm}=0.5\,{\rm pc}$ undergoes a more than $30\%$ reduction in central velocity dispersion, decreasing from $20\,{\rm km}\,{\rm s}^{-1}$ to $13\,{\rm km}\,{\rm s}^{-1}$. In contrast, the GCs with the same mass but a larger initial half-mass radius of $R_{\rm hm}=1\,{\rm pc}$ experiences a relatively smaller ($\sim 20\%$) decline, with $\sigma_{\rm c}$ decreasing from $14\,{\rm km}\,{\rm s}^{-1}$ to $11\,{\rm km}\,{\rm s}^{-1}$.

In Group B, the presence of an IMBH introduces additional dynamical effects. The strong gravitational potential of a massive IMBH significantly enhances the velocity dispersion of stars near it, leading to distinct differences in the initial values of $\sigma_{\rm c}$ depending on the IMBH mass. For an IMBH with $M_{\rm IMBH}=5000\,{\rm M}_{\odot}$, the initial central velocity dispersion is $\sim 24\,{\rm km}\,{\rm s}^{-1}$ at $t=0\,{\rm Myr}$, gradually decreasing to $18\,{\rm km}\,{\rm s}^{-1}$ by $t=30\,{\rm Myr}$ as cluster expands. In contrast, for an IMBH with $M_{\rm IMBH,i}=300\,{\rm M}_{\odot}$, the initial velocity dispersion is significantly lower, starting at $\sim 14\,{\rm km}\,{\rm s}^{-1}$ at $t=0\,{\rm Myr}$ and declining more modestly to $11\,{\rm km}\,{\rm s}^{-1}$ by $t=30\,{\rm Myr}$. In summary, our result suggests that the evolution of the cluster center is highly governed by the properties of GCs, and the mass of the IMBH exerts a significant gravitational influence on the stars in its vicinity.

Figure 9 displays the evolution of the average mass of stars and SBHs within $R_{5\%}$ from the center. The mean stellar mass exhibits a rapid increase early on and stabilizes by $10\,{\rm Myr}$ in all simulations. In the top panel, the initial mean stellar mass differs between clusters of different total masses, being higher in the $5\times 10^{4}\,{\rm M}_{\odot}$ GCs compared to the $10^{5}\,{\rm M}_{\odot}$ cluster. This discrepancy reflects a more rapid relaxation process during the adiabatic growth of IMBH prior to the initialization of simulations (See Section 3.2). Across all simulations in Group A, the central mean stellar mass stabilizes commonly at $\left<M_{\star}\right>\sim 0.9\,{\rm M}_{\odot}$. However, the subsequent evolution of the mean stellar mass depends on the initial properties of the star cluster, reaching values between $\sim 0.7\,{\rm M}_{\odot}$ and $\sim 1\,{\rm M}_{\odot}$ by $t=30\,{\rm Myr}$. Meanwhile, the lower panel highlights a clear decreasing trend in the average stellar mass at the cluster center with increasing IMBH mass. This demonstrates a significant influence of IMBH mass on the mean stellar mass at the cluster center.

We now explore the implications of the initial conditions of a star cluster on the motion of an IMBH. One key parameter is the IMBH’s distance from the center, or the wandering radius, representing the typical distance the IMBH moves from the cluster center due to gravitational interactions with surrounding stars. The wandering radius of the IMBH can be estimated using the following equation (Stone et al., 2017):

While the equation originally uses $R_{\rm c}$ as the core size of the star cluster, there is ambiguity in how the core is defined for a star cluster in our simulations. Therefore, we instead adopt $R_{5\%}$ as a practical proxy. The top-left panel of Figure 10 confirms that the wandering radius increases with the size of the star clusters, consistent with our expectation that larger cores provide more room for the IMBH to move around. The bottom-left panel, however, does not show a clear dependence of $d_{\rm IMBH}$ on $M_{\rm IMBH}$, even though Equation 11 predicts an inverse-square-root scaling.

The principle of energy equipartition provides an estimate for the speed of the IMBH, $|\mathbf{v}_{\rm IMBH}|$ as follows:

where $\left<M_{\star}\right>$ is the average stellar mass, and $\sigma_{\rm c,3D}$ is the 3D velocity dispersion of the stellar system given by $\sigma_{\rm c,3D}=\sqrt{3}\sigma_{\rm c}$. The two left panels in Figure 10 illustrate the evolution of the velocity dispersion in Group A and Group B, respectively. In Group A, $|\mathbf{v}_{\rm IMBH}|$ exhibits a mild decreasing trend with increasing GC size for fixed $M_{\rm GC}$. Meanwhile, Group B shows a clear inverse dependence on $M_{\rm IMBH}$, consistent as expected by Equation 12. However, the scales of IMBH’s wandering radius and the measured speeds in both groups are found to be an order of magnitude larger than the estimates provided by Equations 11 & 12. For example, inserting $M_{\rm IMBH,i}=1000\,{\rm M}_{\odot}$, $\left<M_{\star}\right>=1.0\,{\rm M}_{\odot}$, and $\sigma_{\rm c}=14\,{\rm km}\,{\rm s}^{-1}$ (as seen for $M_{\rm IMBH,i}=1000\,{\rm M}_{\odot}$ in Group B) into Equation 12, we obtain $|v_{\rm IMBH}|\approx 1\,{\rm km}\,{\rm s}^{-1}$. This theoretical estimate is slightly lower than the measured value of $|v_{\rm IMBH}|\approx 3\,{\rm km}\,{\rm s}^{-1}$ presented in the bottom right panel of Figure 10. This suggests that additional effects such as close encounters with SBHs might be influencing the IMBH’s motion beyond the simple equipartition model.

The evolution of the GC cores in the real universe might not be as rapid as seen in these simulations, because the real GCs have longer relaxation timescales due to their greater masses. The shortened relaxation time is an artifact of our simulation design as discussed in Section 3.1. However, as we focus only on the central region surrounding the embedded IMBHs, and the statistical properties of IMRIs during the short runtime, we expect that our chosen GC properties would not significantly affect the conclusion from our analyses.

In this section, we explore the rates and statistical properties of the IMRIs derived from our suite of $N$-body simulations. Figure 11 depicts the IMRI event rates, $\Gamma_{\rm IMRI}$, in the GCs with initial 1D half-mass velocity dispersions ranging $\sigma_{\rm hm}\sim 6\,\mathrm{km}\,\mathrm{s}^{-1}$ to $\sim 10\,\mathrm{km}\,\mathrm{s}^{-1}$ (Group B in Table 1). They are calculated from the stellar and SBH populations within the half-mass radius in each simulation. Each data point represents a value averaged over 10 different random initializations. The initial mass of the IMBH is $M_{\rm IMBH,i}=1000\,{\rm M}_{\odot}$. Here we observe that $\Gamma_{\rm IMRI}$ increases with $\sigma_{\rm hm}$. The IMRI rate is approximately $\sim 0.06\,{\rm Myr}^{-1}$ in a cluster system with $\sigma_{\rm hm}\sim 6\,\mathrm{km}\,\mathrm{s}^{-1}$, but it rises to $\sim 0.2\,{\rm Myr}^{-1}$ in a cluster with $\sigma_{\rm hm}\sim 10\,\mathrm{km}\,\mathrm{s}^{-1}$. This trend suggests that denser clusters, characterized by higher velocity dispersions, are more conducive to the IMRI formation.

Figure 12 highlights the relation between $M_{\rm IMBH}$ and $\Gamma_{\rm IMRI}$ for fixed cluster properties. The IMRI event rate increases from $\sim 0.08\,{\rm Myr}^{-1}$ for a IMBH with $M_{\rm IMBH,i}=300\,{\rm M}_{\odot}$ to $\sim 0.2\,{\rm Myr}^{-1}$ for $M_{{\rm IMBH},i}=5000\,{\rm M}_{\odot}$. By integrating all our simulation results and fitting them to Eq.(4), we obtain the best linear fit that estimates the IMRI event rate:

Note that this relation is built by employing a limited IMBH mass range from $300\,{\rm M}_{\odot}$ to $5000\,{\rm M}_{\odot}$. Although it remains uncertain whether this relation can be extrapolated to the IMBH masses beyond this range, our linear-fit is sufficient to predict the lower mass IMBHs that would occupy the majority of the IMBH population, which is more sensitive to the nature of natal kick. See Section 5.2 for more information.

We now study the orbital properties of the IMRI events from our suite of $N$-body simulations. The left panel of Figure 13 displays a histogram of the orbital eccentricities found in all of our simulations of Group B (see Table 1). It shows that the majority of the IMRIs have high eccentricities, $\left(1-e\right)\sim 10^{-3}-10^{-4}$, at $1\,{\rm Myr}$ before the final merger. Binaries with even higher eccentricities, $\left(1-e\right)\lesssim 10^{-5}$, are less common because it is expected that such high values are more difficult to achieve in the first place. Similarly, encounters with lower eccentricities, $\left(1-e\right)\gtrsim 10^{-2}$, are uncommon because they need to achieve an extremely tiny semi-major axis to satisfy our merger criteria (i.e., only the binaries with $t_{\rm GW}<1.0\,{\rm Myr}$ are merged in the simulation; see Section 2.1). To validate this, in the middle panel of Figure 13 we show a relationship between the semi-major axes $a$ and the eccentricities $e$ from all Group B simulations. Requiring $t_{\rm GW}=1.0\,{\rm Myr}$ for a fixed combination of $M_{\rm IMBH}$ and $M_{\rm BH}$ provide us with a relation $a^{4}\left(1-e\right)^{3.5}=\mathrm{constant}$ (see Eq.(1) in Section 2.1; here we assume $\left(1-e\right)\ll 1$). Indeed, the black dashed line marks such a relation for $M_{\rm IMBH}=10^{3}\,{\rm M}_{\odot}$ and $M_{\rm BH}=60\,{\rm M}_{\odot}$ with $t_{\rm GW}=1.0\,{\rm Myr}$. As expected, we observe that most data points are located below this line. Lastly, the right panel of Figure 13 presents the mass distribution of the secondary BHs involved in the IMRI event. We find that an IMBH tends to form a binary with the most massive BHs in the star cluster.

Figure 14 presents the fraction of observable IMRI events for different GW telescopes, assuming events are detectable if its SNR is larger than 8 (Eq.(2) in Section 2.2). In the first panel, aLIGO is unable to detect IMBHs with $M_{\rm IMBH}\gtrsim 500\,{\rm M}_{\odot}$ because its detection band is limited to frequencies $\gtrsim 10\,\mathrm{Hz}$. Even for $M_{\rm IMBH}\lesssim 300\,{\rm M}_{\odot}$, it can only detect very close IMRIs in the local universe ($z\lesssim 0.06$). Consequently, we can reasonably conclude that aLIGO would not be able to discover a merger event involving an IMBH, unless there is one with mass $\lesssim 300\,{\rm M}_{\odot}$ in the vicinity of Milky Way ($\lesssim 260\,{\rm Mpc}$).

ET has a better capability to detect IMRIs compared to aLIGO, especially in the lower-mass range of $M_{\rm IMBH}$. For example, the second panel of Figure 14 shows that it has the potential to detect almost all $300\,{\rm M}_{\odot}$ IMBHs up to $z\sim 1.0$, and some up to $z\sim 1.5$. ET also offers a relatively wider detection mass spectrum up to $M_{\rm IMBH}\sim 2000\,{\rm M}_{\odot}$, although the horizon redshift decreases to $z\lesssim 0.1$. Note that the maximum redshift of detection — or the “horizon redshift” — decreases as $M_{\rm IMBH}$ increases.

On the other hand, we find that LISA, aSOGRO, and DECIGO can observe all IMRIs involving IMBHs across the entire mass range explored. However, their horizon redshifts vary significantly from each other. Figure 15 illustrates the horizon redshifts of each GW observatory across the IMBH masses $M_{\rm IMBH}=300\,{\rm M}_{\odot}$ to $10^{4}\,{\rm M}_{\odot}$. LISA has a growing horizon redshift for increasing $M_{\rm IMBH}$ from $z\sim 0.2$ to $\sim 0.5$. Especially, it can better cover the detection of IMRIs including $M_{\rm IMBH}\gtrsim 2000\,{\rm M}_{\odot}$ than ET. Meanwhile, aSOGRO’s horizon redshift is confined within the local universe $z\lesssim 0.05$. Although its detection range diminishes for higher mass IMBHs, DECIGO can detect IMRI events occurring at redshifts as high as $z\sim 10$.

In conclusion, both aLIGO and DECIGO may serve as a local probe for low-mass IMBHs, but aSOGRO covers a broader mass range. Meanwhile, ET and LISA can extend the cosmic volume under our detection up to $z\sim 0.5-1.0$, complementing each other in their detection capabilities across different mass ranges. Lastly, DECIGO can push the observational limits beyond $z\sim 10$ across all mass ranges.

We now estimate the IMRI detection rates by three different GW observatories, employing the GC number density as discussed in Section 2.3. In this work, the distribution of half-mass velocity dispersions, $\sigma_{\rm hm}$, of GCs is derived from

Here, we use the relation $\sigma_{\rm hm}=\sqrt{GM_{\rm GC}/6R_{\rm hm}}$ to estimate velocity dispersion, assuming a typical GC size of $R_{\rm hm}\sim 2.5\,{\rm pc}$, consistent with observations of GCs in the Milky Way (Portegies Zwart et al., 2010; Shin et al., 2013).

Combining and integrating Eqs.(6), (13) and (14) yields the expected IMRI detection rates per unit comoving volume at each redshift. Figure 16 shows the resulting detection rates for two natal kick scenarios (see Section 2.3). In both scenarios, all observatories except for aLIGO have the rate $\sim$ a few $\times 10^{-8}\,\mathrm{cMpc}^{-3}\,{\rm yr}^{-1}$ at $z\lesssim 0.3$. First of all, its detection range encompasses the entire history of IMRI mergers, since DECIGO may cover all events up to a redshift $\lesssim 10$. LISA can detect some IMRIs up to the $z\lesssim 0.5$, while ET’s detection range extends slightly farther, up to $z\lesssim 1.0$. Despite their similar detection rates at $z<0.3$, it is important to note that they are complementary in terms of the mass spectrum they can observe (see Figure 15 and Section 5.1).

Meanwhile, aLIGO and aSOGRO have much narrower detection ranges ($z\lesssim 0.06$); However, aLIGO exhibits a remarkably lower detection rate per unit comoving volume of $\lesssim 2\times 10^{-9}\,\mathrm{cMpc}^{-3}\,{\rm yr}^{-1}$ compared to aSOGRO. This gap can be attributed to aLIGO’s limited detectable mass range even in the optimistic scenario with low natal kicks. In the pessimistic-case scenario in which many lower-mass IMBHs are ejected due to strong natal kicks, aLIGO will barely detect any IMRIs, even in very close proximity ($z<0.01$). Figure 16 illustrates that, for some GW telescopes, the IMRI detection rate can be susceptible to the intensity of natal kicks assumed in our model.

By integrating the above results up to the horizon redshift, we derive the final conclusions of our study — that is the IMRI detection rates per year for different GW telescopes under two scenarios of BH natal kicks. The detection rate for aLIGO is found as

where we include the dependency on $n_{{\rm GC},\,z=0}$ in the formula. The aLIGO mainly targets low-mass IMBHs which are vulnerable to the recoil kick in general (Heggie, 1975; Antonini & Rasio, 2016; Arca Sedda et al., 2023). This makes it challenging to constrain the detection frequencies for IMBHs with $M_{\rm IMBH}<10^{3}\,{\rm M}_{\odot}$. In extreme cases, these IMBHs may be completely ejected from their host GCs during evolution. This uncertainty is reflected in Table 1 of Askar et al. (2022), which shows the highest variance in the retention of low-mass IMBHs. In summary, our results indicate that aLIGO will likely detect very few IMBHs, even under the most optimistic assumptions.

Meanwhile, the detection rate of upcoming GW observatories will significantly improve the chance of IMRI detection. The detection rate of the rest of GW observatories are as follows. The rate of ET is estimated as

offering ample chances of IMRI detection for $M_{\rm IMBH}\lesssim 10^{3}\,{\rm M}_{\odot}$. The detection rate of LISA is

LISA appears to be less affected by the natal kicks, as more massive IMBHs dominate its detectable IMRI populations (see Figure 15 and Section 5.1). The detection rate of aSOGRO is given as

Lastly, the detection rate of DECIGO is estimated as

making it the most promising GW observatory for detecting IMRIs among those considered in this study.

According to our analyses, ET, LISA, and DECIGO are expected to provide robust detections of IMRIs, thereby greatly improving our ability to study the IMBH population over a wide range of masses and redshifts. Although aSOGRO exhibits a lower detection rate compared to other future observatories, it remains capable of detecting IMRI events within a practical observational timescale ($\sim 5-10\,{\rm yr}$).

In this section, we discuss the limitations of our study, focusing on the factors that are not considered in our analyses.

• •

  Our simulations assume that the IMBH-BH binary tightens only via GW radiation. However, a binary may undergo multiple encounters with other stars or SBHs nearby; and, these “third” objects passing close to the center of mass of the hard binary can deform the binary configuration (Saslaw et al., 1974; Merritt & Milosavljević, 2005). As a result, on average, the actual merger timescale could be shorter than our simulation results. More importantly, these interactions can significantly alter the characteristics of the observed GW signals. Also, the general relativistic effect such as the capture by gravitational capture (Quinlan & Shapiro, 1989) is not included in our simulation. Given the rare frequency of the gravitational capture events, this would not affect our result significantly (See Appendix C).

• •

  Our study also primarily focuses on the GW signals during the inspiral stage since it is the dominant source. However, despite its short duration, the merger and ringdown phases can also produce loud signals in some cases. Consequently, aLIGO and ET may have the potential to observe a broader mass spectrum than our predictions.

• •

  In a realistic GW observation, multiple aspects can affect such as the Earth rotation (e.g., Isi et al., 2015). However, the SNR calculation utilizing Eq.(2) will not be significantly different from the precise measurement as an order of magnitude sense. Therefore, the actual IMRI detection rate will not be greatly far from our estimation.

• •

  The time evolution of IMRI rates may depend on the formation mechanism of the IMBH. IMBHs can form through two distinct scenarios: the slow and fast scenarios (Giersz et al., 2015). In the slow scenario, an IMBH grows gradually over time by accreting stars and BHs, allowing IMRI events persists up to the present day. In contrast, in the fast scenario, where the IMBH forms via the runaway merger of massive stars, the majority of IMRI events are expected to occur within the first $5\,{\rm Gyr}$, depleting the sources of IMRIs by low redshifts (Hong et al., 2020).

• •

  The mass range of IMBHs is constrained to $M_{\rm IMBH}<10^{4}\,{\rm M}_{\odot}$ in our analysis. However, observatories such as LISA and DECIGO may be also sensitive to $M_{\rm IMBH}>10^{4}\,{\rm M}_{\odot}$. Consequently, the actual detection rates for these two observatories are likely to exceed the prediction of our analysis.

• •

  The low-mass IMBHs are vulnerable to the GW recoil kicks from the IMBH-BH merger (González Prieto et al., 2022). The GW kick may eject most of the low-mass IMBH from the dense star clusters, and this would result in the decrease of the IMRI rate for IMRIs with $\left(M_{\rm IMBH}<10^{3}\,{\rm M}_{\odot}\right)$. The impact of GW kick on the IMRI detection rate will be rather limited since detectors such as LISA and aSOGRO are sensitive to higher mass regimes. Meanwhile, ET and aLIGO may detect less frequently than our estimations although we already concluded that aLIGO would be less likely to detect any IMRIs in a practical timescale.

• •

  The detailed evolution of the IMRI rate across cosmic time may exhibit significant variations depending on the assumptions on the GC formation and evolution history. For example, a semi-analytic model work of Fragione et al. (2018) accounts for the mass loss of GCs due to the galactic tidal field. Their study predicts that at low redshift, the IMRI event rate falls within the range of $\sim 0.5-3\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$. However, during the peak GC formation epoch, this rate increases as high as $\sim 100\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$. In contrast, our simulations, which focus on short-term GC evolution, predict an order of magnitude higher IMRI rates than those found in (Fragione et al., 2018), while it declines at higher redshifts because of the decreasing GC number density.