Generation of gravitational waves and tidal disruptions in clumpy galaxies

Obtaining a better understanding of intermediate-mass black holes (IMBHs) is crucial, as their properties could shed light on the origin and growth of their supermassive counterparts. Massive star-forming clumps, which are present in a large fraction of massive galaxies at $z \sim$ 1-3, are amongst the venues wherein IMBHs could reside. We perform a series of Fokker-Planck simulations to explore the occurrence of tidal disruption (TD) and gravitational wave (GW) events about an IMBH in a massive star-forming clump, modelling the latter so that its mass ($10^8 \,{\rm M}_{\odot}$) and effective radius ($100$ pc) are consistent with the properties of both observed and simulated clumps. We find that the TD and GW event rates are in the ranges $10^{-6}$-$10^{-5}$ and $10^{-8}$-$10^{-7}$ yr$^{-1}$, respectively, depending on the assumptions for the initial inner density profile of the system ($\rho \propto r^{-2}$ or $\propto r^{-1}$) and the initial mass of the central IMBH ($10^5$ or $10^3\,{\rm M}_{\odot}$). By integrating the GW event rate over $z$ = 1-3, we expect that the Laser Interferometer Space Antenna will be able to detect $\sim$2 GW events per yr coming from these massive clumps; the intrinsic rate of TD events from these systems amounts instead to a few $10^3$ per yr, a fraction of which will be observable by, e.g. the Square Kilometre Array and the Advanced Telescope for High Energy Astrophysics. In conclusion, our results support the idea that the forthcoming GW and electromagnetic facilities may have the unprecedented opportunity of unveiling the lurking population of IMBHs.

In the last decades, observational astronomers have been pointing to an apparent dearth of BHs with masses between the heaviest SBHs (∼100 M ) and the lightest SMBHs (∼10 6 M ) (for a review, see Mezcua 2017): BHs inhabiting this loosely populated mass range have been named intermediate-mass BHs (IMBHs). While the lack of SBHs with masses above 100 M can be justified within the current stellar evolutionary framework (Mapelli 2018), 1 it is 2007; Maccarone & Servillat 2008;Noyola et al. 2010;Lützgendorf et al. 2012Lützgendorf et al. , 2013Lanzoni et al. 2013;Haggard et al. 2013;Reines et al. 2013;Casares & Jonker 2014;Baldassare et al. 2015;Kızıltan et al. 2017;Perera et al. 2017;Lin et al. 2018;Mezcua et al. 2018) and many of them remain debated. It is unclear whether the paucity of IMBH detections has to be interpreted as an intrinsic lack of these objects, or rather observational biases render their detection inherently more challenging (e.g. Mezcua 2017).
The idea that IMBHs could inhabit the centre of globular clusters dates back to the late 1960s (Wheeler 1968;Wyller 1970); from that moment on, a number of theoretical studies focussed on the formation path of IMBHs within their hosts. For instance, Miller & Hamilton (2002) proposed the formation of an IMBH starting from a massive SBH that increases its mass thanks to the coalescence with smaller SBHs. Alternatively, a series of collisions among massive stars, enhanced by the segregation of these bodies at the centre of their hosts system and the resulting increase in the central density, could lead to the growth of a supermassive object in a runaway fashion (e.g. Portegies Zwart & McMillan 2002a;Portegies Zwart et al. 2004a), possibly aided by the presence of gas (Reinoso et al. 2020). Additionally, Giersz et al. (2015) suggested that an IMBH could result from dynamical interactions of hard binaries containing a SBH with other bodies in the globular cluster. A significant number of the proposed formation scenarios relies on dynamical processes that bring to a series of collisions within the host system. For this, very dense stellar environments subject to mass segregation, core collapse, and binary interactions seem to be ideal environments for the formation of IMBHs (e.g. Gürkan et al. 2004;Gültekin et al. 2004;Gürkan et al. 2006;Freitag et al. 2006).
In the near future, thanks to new observatories such as the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017; Barack et al. 2019), TianQin (Luo et al. 2016), the Square Kilometre Array (SKA; Dewdney et al. 2009), the extended ROentgen Survey with an Imaging Telescope Array (eROSITA; Merloni et al. 2012), and the Advanced Telescope for High Energy Astrophysics (Athena; Barcons et al. 2012), we will have the capability to observe IMBHs via new strategies, detecting the signal that these objects may generate as a consequence of the interaction with their surroundings up to the relatively high-redshift Universe. Such events could manifest themselves both in the form of gravitational-wave (GW) signals -in particular from extreme and intermediate mass-ratio inspirals (EMRIs and IMRIs; e.g. Amaro-Seoane et al. 2007) -and tidal disruption events (TDEs; e.g. Kashiyama & Inayoshi 2016). The search for this kind of events about IMBHs is crucial, considering that these objects could constitute the high-redshift seeds for the SMBHs we observe at the present day, thus the aforementioned detections could provide new crucial insights on the cosmic origin and growth of SMBHs (Volonteri 2010).
Recent observations have shown that a large fraction ( 50 per cent; Shibuya et al. 2016) of massive galaxies in the redshift range z ∼ 1-3 present massive star-forming clumps, with typical masses of ∼10 7 -10 8 M and, in some cases, even ∼10 9 M (e.g. Elmegreen et al. 2009). These clumps are likely to constitute a privileged setting for the formation of IMBHs and, consequently, for the occurrence of IMRIs and TDEs that could be potentially detected by the forth-coming facilities. The existence of such numerous massive clumps in many galaxies at the peak of star formation and galaxy merger rates is corroborated both by spatially resolved observations (e.g. Huertas-Company et al. 2020) and numerical simulations (e.g. Tamburello et al. 2015).
For this reason, in this paper we adopt a numerical approach to estimate the number of absorption events 2 occurring in a typical massive stellar clump, whose properties are extrapolated from the isolated simulation of a star-forming galaxy by Tamburello et al. (2015), in the assumption that an IMBH is found within the host system.
The paper is organised as follows: in Section 2, we describe the numerous aspects of our model and the numerical setup; in Section 3, we present our results; in Section 4, we discuss our findings, and we conclude in Section 5.

METHODS
In order to estimate the rates of absorption events occurring in stellar clumps, we first obtain the characteristic parameters of a typical clump from an isolated simulation of a starforming galaxy (Section 2.1). We then use such quantities to build analytical models with different density profiles (Section 2.1.1), that are used as the initial conditions adopted in the Fokker-Planck simulations (see Section 2.2) performed to keep track of the event rates.

Initial Conditions
In order to build a realistic model for the star-forming clumps in which IMBHs may reside, we rely on the numerical work of Tamburello et al. (2015), where the properties of massive star-forming clumps in isolated galaxies were studied in detail. They performed a large suite of galaxy simulations, to assess the dependence of gas fragmentation properties in a galaxy on sub-grid physics, galaxy mass, structural parameters, and resolution, using the smoothed particle hydrodynamic code Gasoline2 , an updated version of Gasoline (Stadel 2001;Wadsley et al. 2004). In this paper, we specifically analyse the clumps forming in their run 27 (see tables 1 and 2 of Tamburello et al. 2015), in which they simulated a dark matter halo with an embedded stellar and gaseous disc, with a virial and baryonic mass of 2.5 × 10 12 and 7.8 × 10 10 M , respectively, a gas fraction of 50 per cent, and a virial concentration of 6. This choice of parameters was particularly conducive to the formation of long-lasting massive clumps. Star formation and stellar (blast-wave) feedback were modelled following the prescriptions of Stinson et al. (2006), with a star formation efficiency, threshold temperature, and threshold density of 0.01, 3 × 10 4 K, and 10 mH g cm −3 , respectively, and an energy injection of 4×10 50 erg per supernova explosion. The initial numbers of particles were 1.2 × 10 6 , 10 5 , and 10 5 , for 2 In this work, we define absorption an event in which a stellarmass object (star or compact object) is "absorbed" by the central IMBH. We then distinguish between TDEs and GW events, which can be subsequently divided into EMRIs, IMRIs, or plunges, depending on the mass ratio and on the number of orbits of the stellar-mass compact object around the IMBH. dark matter, gas, and stars, respectively, whereas the softening, identical for all particles, was set equal to 100 pc.
Using Skid (Stadel 2001), a group finder for N -body simulations, we extracted the mass and half-mass radius of the largest clumps formed in the simulation, finding a typical stellar clump mass compatible with what was observationally found by Elmegreen et al. (2009) in their Hubble Ultra-Deep Field analysis of clump-cluster and chain galaxies, and a characteristic half-mass radius r hm ≈ 100 pc.
We caution that the mass and length resolution in the simulation by Tamburello et al. (2015) are too coarse to perform a self-consistent modelling of the clumps starting from the density profile obtained from the simulation; for this reason, we create analytical models using Equations (1-2) as constraints, as described in Section 2.1.1.

Composition and structure of the cluster
For simplicity, we assume that the entire gaseous mass residing in the clusters described in Tamburello et al. (2015), which is of the order of the stellar mass, collapses in stars and/or is evaporated from the clump by means of stellar/BH feedback. 3 Such assumption allows us to model the cluster as a dissipationless system composed only by stellar objects. Specifically, we assume the cluster to initially follow a spherical density profile belonging to the Dehnen (1993) family (see also Tremaine et al. 1994): where a is the scale radius of the system, M clus its total mass, and γ ∈ [0; 3) its logarithmic inner density slope. The half-mass radius of such systems is related to the scale radius via the relation In this work, we only model two cases: a Jaffe (1983) profile (γ = 2, default case) and a Hernquist (1990) profile (γ = 1). The choice of the default configuration (γ = 2) is motivated by the fact that violent relaxation and dissipational collapse could produce a system with a density slope γ ≈ 2, as hinted by past studies (see e.g. Hjorth & Madsen 1991;Treu & Koopmans 2002); we additionally modelled the cluster with a Hernquist density profile, to assess the dependence of our results on the choice of the initial profile (see Section 3.2.1).
We assume that the clump is composed of one central IMBH (see below) and four stellar families, each following the same density profile at the beginning of the integration: main sequence stars (MSs), white dwarfs (WDs), neutron stars (NSs), and stellar-mass BHs (SBHs). The Fokker-Planck treatment, and in particular the PhaseFlow code described below (see Section 2.2) can only handle a finite number of components, each having a fixed mass. Hence, we could not adopt a smooth and non-discretized mass function, but assigned instead a typical mass to each of the four components, whose values are listed in Table 1. Such values are derived assuming that stars were born with a Kroupa (2001) initial mass function and the stellar population is 1 Gyr old (i.e. the cluster is assumed to have a turnoff mass of 2.512 M ), which is the typical time-scale within which massive clumps are expected to dissolve in the tidal field of their host galaxy (Tamburello et al. 2015).
In particular, we adopt the initial-to-final mass relation proposed by Merritt (2013). The mass of each object in a given species is shown in Table 1: for each component, we selected a default mass value (in bold in the table) and performed additional runs varying it (as described in Section 3.2.3). Combining the object masses with Equation (1), we obtain the numbers of objects in the default configuration, as shown in Table 1.
We assume that an IMBH of mass m• exists at the centre of the massive clump. There are multiple ways to achieve this in Nature. In particular, the IMBH could be either formed in-situ via runaway collapse of relatively massive stars, or be captured by the massive stellar clump. We speculate on the possible channels of in-situ formation or capture in Section 4. Here we assume that an IMBH with an initial mass m•0 = 10 5 M lies at equilibrium at the centre of the system; this mass choice implies that the cluster and its IMBH sit on the extrapolated bulge mass-BH mass relation observed in more massive systems in the local Universe (Kormendy & Ho 2013). We additionally performed a test with m•0 = 10 3 M , to assess how our results change when we vary the initial IMBH mass (see Section 3.2.2).

Absorption events
Our main aim is to obtain the rates of TDEs and GW events occurring in the modelled stellar clump hosting an IMBH.
A TDE occurs when a MS or a WD (when we assume that WDs can be tidally disrupted; see Section 3.2.3) approaches the close vicinity of the IMBH and is tidally disrupted as the IMBH gravitational attraction overcomes the self gravity of the stellar object. The threshold separation for a TDE to occur is named the tidal radius, and depends on the stellar mass, m , radius, r , and internal structure, scaling as (Hills 1975;Rees 1988) Table 1. Parameters associated to the four stellar-remnant families -MSs, WDs, NSs, and SBHs -included in our runs: (i) initial mass range, (ii) name of the stellar-remnant species, (iii) final number fraction associated to each species (assuming a 1 Gyr old stellar population), (iv) mass of each single object in the family (we explored three values; the default is in bold), (v) respective stellar radius (for MSs and WDs), and (vi) number of objects. The total number of objects is 1.6 × 10 8 for all cases (meaning that the total stellar mass can vary by up to ∼30 per cent when we change the object mass; see Section 3.2.3).
set equal to 1. The radii for MSs and WDs used in our runs are listed in Table 1: we extrapolated the MS radii from the empirical mass-radius relation in Chen et al. (2014), whereas we used the relationship in Magano et al. (2017) for the WD radii.
If the stellar object is too compact, it is expected to enter the IMBH horizon prior to undergoing a TDE. In fact, a gravitational capture, and the possibly associated EMRI or IMRI (depending on the mass ratio), 4 occurs when a stellarmass compact object (such as a SBH, a NS, or a WD -if the latter is not tidally disrupted) finds itself on an orbit that brings it inside the IMBH event horizon. The capture radius can be approximated as where G is the gravitational constant, c is the speed of light in vacuum, and α depends on the orbital properties of the inspiralling object and on the IMBH spin. 5 If we consider a compact object inspiralling on a nearly circular orbit about a non-spinning IMBH, then α = 6, so that rcap equals the radius of the IMBH innermost stable circular orbit (ISCO). However, the majority of compact objects are expected to reach the IMBH vicinity via two-body scatterings with other stars, that settle them on nearly radial orbits (e.g. Hils & Bender 1995;Bortolas & Mapelli 2019); if such orbits do not circularise due to dynamical friction or GW emission, the critical α for a capture to occur is α = 8, in the assumption of a non-spinning central IMBH . The aforementioned parameter takes different values if the central IMBH is spinning: in particular, α = 2 (11.6) for a compact object inspiralling on a non-circular, prograde (retrograde) orbit on to a maximally spinning IMBH . While all compact objects crossing rcap are expected to be accreted on to the central BH, not all of them can produce a detectable GW event. In fact, it is common use to distinguish between direct plunges and inspirals. In the former case, the compact object is directly scattered on to 4 In this work, we assume the inspiral's mass ratio to be extreme (EMRI), intermediate (IMRI), or moderate, when the mass ratio of the two objects is q < 10 −4 , 10 −2 > q ≥ 10 −4 , and q ≤ 10 −2 , respectively. The choice of these mass-ratio thresholds does not affect the results. 5 We do not consider the effect of the inspiralling compact object's spin, as it is assumed to be negligible (e.g. Barack & Cutler 2004; however, see Han 2010).
rcap; in the latter case, the compact object gradually inspirals on to the IMBH via emission of GWs, possibly producing a long-lived signal that can be detected by forthcoming observatories such as LISA. The Fokker-Planck approach does not allow us to distinguish between EMRIs/IMRIs and direct plunges; for simplicity, we will refer to both as GW events, and we will comment upon the rates of both events in Section 4.
In what follows, we refer to both the tidal or capture radius (depending on the species) as the absorption radius, r abs .

Fokker-Planck Simulation
The number of GW events and TDEs occurring in the stellar clump is obtained via the PhaseFlow code (Vasiliev 2017), which is part of the publicly available Agama library (Vasiliev 2019). PhaseFlow is a Fokker-Planck code that computes the dynamical evolution of a spherical isotropic system by solving the coupled system of Poisson and orbitaveraged Fokker-Planck one-dimensional equations for the gravitational potential, the density, and the distribution function; note that PhaseFlow is the first code adopting the phase volume instead of the energy as the argument for the distribution function, thus rendering the code extremely efficient. PhaseFlow can simultaneously handle several species of stars (components), each with its own distribution function. In addition, it can account for the presence of a central BH, and it incorporates a sink term that can mimic the loss of stars on to the central BH, thus delivering the rates of TDEs and GW events. 6 PhaseFlow has been successfully tested in the context of, e.g. nuclear star clusters (Generozov et al. 2018;Emami & Loeb 2020), Bahcall & Wolf (1976) cusps (Vasiliev 2017), primordial BH mass functions (Zhu et al. 2018;Stegmann et al. 2020), and TDE rates in galaxy merger remnants and nuclear star clusters (Pfister et al. 2019(Pfister et al. , 2020.
In the present study, the IMBH is incorporated in the simulation assuming it is at equilibrium with the surrounding system since the beginning of the integration, and the initial distribution function is constructed in the combined potential of the IMBH and the stellar-mass objects. While this is of course an approximation, violent relaxation that could have acted in the system at the stages of violent in-situ star formation and IMBH formation may have brought the system to dynamical relaxation even if its lifespan is shorter than its relaxation time.

Simulation parameters
In this section, we describe the parameters of our simulations, focusing first on the default configuration.
We account for the presence of four species of stellarmass objects -MSs, WDs, NSs, and SBHs -each distributed according to the same density profile: a Jaffe profile (see the grey-shaded region in the top panel of Figure 1) with total mass of 10 8 M and half-mass radius r hm,D = 100 pc. As mentioned in Section 2.1.1, the number fraction assigned to each stellar species and the mass of each single object are displayed in Table 1.
The central IMBH is given an initial mass m•0 = 10 5 M , and its mass increases via TDEs and GW events as the system evolves, the 'absorption' boundary for such events being set according to Equations (5) and (6) (with α = 8), respectively; in both cases, this radius increases during the simulation as the IMBH mass grows. We also assume that WDs get tidally disrupted rather than undergoing an IMRI event, in agreement with, e.g. Haas et al. (2012) and Anninos et al. (2018).
When a stellar object (MS or WD) gets tidally disrupted, 30 per cent 7 of its mass is accreted on to the IMBH, whereas the whole compact object's mass is added to the IMBH when a GW event occurs.
We set the Coulomb logarithm, ln Λ, which tunes the significance of two-body relaxation in the Fokker-Planck code, equal to 16, since ln Λ ≈ ln (r hm /b90), where b90 ∼ 10 −5 pc is the 90 • deflection radius for stellar encounters.
Finally, we initialized the Fokker-Planck solver so that the distribution function is computed on a grid in phase volume with 400 logarithmically spaced points, and we chose the quadratic discretization method (method 2 in Phase-Flow; Vasiliev 2017) and an accuracy parameter for the time integration of 10 −2 . We checked that our results do not depend on the choice of such parameters. We evolved each realisation for 1 Gyr, the typical massive-clump lifetime (Tamburello et al. 2015).

RESULTS
In this section, we first present the rate of absorption events when considering the default scenario (Section 3.1) and then show how the results change when we vary the physical parameters of our simulation (Section 3.2). Figure 2 shows the rate of absorption events (TDEs for MSs and WDs; GW events for NSs and SBHs) as a function of Figure 2. Rate of absorption events (TDEs for MSs and WDs; GW events for NSs and SBHs) as a function of time occurring in the default configuration (black solid line; see Section 3.1), in the case adopting a Hernquist instead of a Jaffe initial density profile (red dashed line; see Section 3.2.1), and in the case assuming an IMBH with initial mass of 10 3 M instead of 10 5 M (blue dotted line; see Section 3.2.2). The event rates associated to the default case are typically larger than those obtained in the other two scenarios by at least one order of magnitude, except for the case of SBHs.

Default scenario
time for the default case (together with other scenarios discussed below), whereas Table 2 lists the average rates of events per yr and the relative number fraction of objects belonging to each species captured or disrupted by the IMBH. Overall, the highest event rates in the default scenario are associated to MSs (3.7 × 10 4 TDEs per Gyr), followed by SBHs (3.7 × 10 2 IMRIs per Gyr), WDs (49 TDEs per Gyr), and NSs (42 EMRIs per Gyr). This reflects the fact that MSs are the most abundant species in the system ( Table 1). The absorption event rate of SBHs is larger than that of WDs and NSs, despite the former species having a lower number fraction than the latter two. This can be understood in terms of mass segregation: SBHs have the largest mass amongst the four components by a large margin, hence they sink towards the centre of the system while pushing the lighter objects at larger distances from the centre. This fact can be appreciated by looking at the inner region in the top panel of Figure 1, where the SBH profile has initially the lowest normalisation and, at the end of the integration, ends up with the largest.
The SBH mass segregation can occur over ∼0.2 relaxation times either in the strong or in the weak (standard) regime (Preto & Amaro-Seoane 2010;Merritt 2013). In the weak (standard) regime, SBHs are numerous enough at all radii within the IMBH influence sphere (defined as the region enclosing twice the IMBH mass) to dominate the relaxation rate, and they develop a so-called Bahcall & Wolf cusp, with a logarithmic density slope ρ ∝ r −7/4 . In the strong regime, a range of radii within the IMBH influence sphere exists where SBHs do not dominate the density profile, even after segregation occurred; in that outer region, they settle in a steeper cusp, ρ ∝ r where ρH is the mass density of heavy bodies (i.e. SBHs), ρL that of light bodies (i.e. MSs; we neglect here WDs and NSs, as they are relatively rare), and, in this particular case, m H = 20 M and m L = 0.6 M are the masses of a single body of the two populations; if ∆ 1 ( 1), the system is in the strong (weak) segregation regime. In our system, ∆ 1 within the innermost parsec, implying that the strong regime is never in place. The same applies to the other systems described in the next sections.
The rate of TDEs and GW events is not constant along the simulation, as displayed in Figure 2. In the first 0.1 Myr, the rate of events can be orders of magnitude larger than that at later times. In fact, within the IMBH influence sphere, the system almost immediately develops a Bahcall & Wolf cusp in the SBHs and a profile that scales as ρ ∝ r −1.45 in the other components, resulting in a drop of the density profile of all species. After the cusp has developed, the system further slowly evolves, expanding in a self-similar fashion and producing a slow drop of the inner density normalization. This expansion is powered by the energy source of the IMBH, which absorbs stars with large negative energy, thus acting as a heat source, as discussed in detail in Vasiliev (2017). The system expansion and associated drop in the central density produces the slow decrease in all the event rates at late times. If one compares the late event rate for any species (from Figure 2) with the average rate shown in Table 2, the values differ at most by factors of a few. This means that, even neglecting the initial burst of events, our results remain valid at least at the order-of-magnitude level.
Finally, it is worth mentioning that the fractional mass accreted by the central IMBH via TDEs and GW events by the end of the integration is ∼0.14 (approximately half of which from MSs and the other half from SBHs), i.e. the IMBH does not grow by a significant amount, suggesting that the IMBH mass growth has little or no impact on the event rates.

Dependence on physical parameters
In order to assess the dependence of the results on the physical parameters of our system, we ran a suite of control simulations in which we varied several quantities. In particular, we explored the effects of varying (i) the initial stellar density profile (Section 3.2.1) and (ii) the initial mass of the central IMBH (Section 3.2.2). In Section 3.2.3, we further appraised the impact of having a different mass (and, in the case of MSs and WDs, corresponding radius) of the stellarmass objects, a different size of the capture radius for GW events (i.e. varying α), a different absorption process for WDs (EMRIs instead of TDEs), and a different stellar population age.

The initial stellar density profile
The initial stellar density profile likely plays an important role in the determination of the absorption event rates. In this section, we study the evolution of the system when we replace the Jaffe model (γ = 2 in Equation 3) with a Hernquist model (γ = 1), keeping the initial cluster mass and half-mass radius equal to 10 8 M and 100 pc, respectively, thus setting up a model with an initial scale radius a = 41.4 pc (as opposed to 100 pc).
The average absorption event rates in this alternative scenario are listed in Table 2, whereas Figure 2 shows the rate of events as a function of time. A Bahcall & Wolf cusp develops within the IMBH influence radius (≈ 2 pc at the beginning of the integration, ≈ 2.3 pc at the end), as expected (Figure 1). The event rate for each species drops by a factor ∼20-30 with respect to the default case. This can be explained in terms of the different density at small scales (compare the two top panels in Figure 1): as an example, the average initial (final) total density within a radius of 1 pc is 2.4 × 10 5 M pc −3 (7.9 × 10 4 M pc −3 ) for the Jaffe case, and 1.4 × 10 4 M pc −3 (7.8 × 10 3 M pc −3 ) in the Hernquist case, i.e. a factor of ∼20 (∼10) smaller; the proportional dependence of the event rates on the inner density is expected from the loss-cone theory (e.g. Merritt 2013). If the original Hernquist density profile had a larger central density in the beginning, it is reasonable to expect that the event rates would converge to the values obtained in the default scenario.
Note that, in this case, the event rates do not always decrease with time: there is typically an initial rise that culminates at ≈0.2 Gyr, followed by a relatively slow decay. This behaviour is a result of a sequence of events: first, the cusp develops within the IMBH influence radius, resulting in the growth of the event rates (in this scenario, the initial profile is less steep than the Bahcall & Wolf cusp, as opposed to the default case); after it is fully developed, the whole system again slowly expands in a self-similar fashion, thus the central density drops, and so does the event rate of all species. It is interesting that the late-time evolution of the event rates is similar to that of the default scenario, with just a different normalization, owing to the difference in the central mass density of the two systems.
Here the IMBH only grows by ≈700 M by the end of the integration. Most of the growth is due to accretion of MSs and SBHs (in almost equal part), as in the default case.

The initial mass of the intermediate-mass black hole
Given the uncertainty on the mass of the IMBH hosted in the system, we here explore the impact on the event rates when assuming that the IMBH has an initial mass m•0 = 10 3 M , i.e. 100 times smaller than in the default scenario. The IMBH mass is a critical quantity for the estimation of the TDE and GW event rates, as both the tidal and capture radii grow with the IMBH mass, even if in a different fashion (see Equations 5 and 6). Table 2 and Figure 2 report the rate of events in this configuration, showing that lowering the IMBH mass brings to a drop -compared to the default scenario -in the event rates for all species, although by vastly different amounts. The drop is by a factor of ∼30 for MSs and WDs, the tidally disrupted objects, and by a factor of ∼90 for NSs. On the contrary, GW events associated to SBHs drop only by a factor of ∼3. In addition, in this situation the IMBH grows by a substantial amount, en- . Radius at which the flux of events peaks, for three different scenarios. Such radius depends on the considered species, and in particular here we show its value for MSs (empty symbols) and SBHs (full symbols). Most of the flux comes from a radius that is roughly one order of magnitude smaller than the IMBH influence radius (evaluated at late times) for MSs, and from a radius which is nearly 5 (15) times smaller than the MS maximum flux radius if we consider SBHs, in the cases for which the IMBH has m •0 = 10 5 M (10 3 M ).
hancing its mass by a factor of ∼4 in 1 Gyr. Nearly 91 per cent of such growth is to be attributed to SBHs, whereas MSs grow the IMBH mass by about 8 per cent. The crucial contribution of SBHs is likely due to the following facts: (i) they start to dominate the density profile at larger radii ( Figure 1) and (ii) the radius from which most of the SBH flux comes from, shown in Figure 3, is much smaller in this situation (almost by a factor of 100), compared to the other two discussed scenarios. It follows that SBHs contribute more to the total capture rate.
We also note that the IMBH influence radius increases from an initial value of ∼2×10 −3 pc to ∼0.12 pc by the end of the integration, in response to the mass growth of the IMBH itself and to the evolution of the inner density profile.

Robustness of the results
In the previous sections, we computed the absorption rates of different types of stellar remnants after fixing the values of several quantities related to both the stellar-mass objects and the IMBH. Here we show that our choice of parameters does not significantly affect the results.
We first varied the mass of MSs (between 0.4 and 0.8 M ), WDs (between 0.7 and 1.3 M ), NSs (between 1.6 and 3 M ), and SBHs (between 10 and 30 M ). When varying the mass of MSs and WDs, we altered the stellar radius accordingly (see Table 1). To better assess the effect of the change in each stellar population, we modified the stellar masses in a separate run for each component.   (5) and (6), where for MSs and WDs the default values listed in Table 1 were used.
Given the uncertainty in the size of the capture radius (which depends on quantities that we cannot follow, such as the orbital properties of the inspiralling body and the IMBH spin), we then assessed the dependence of our results on the choice of the parameter α (see Equation 6), using the values listed in Table 3.
We also investigated the impact of having a different absorption process for WDs, by assuming that, instead of undergoing a TDE (the default case), they produce an EMRI (in practice, employing Equation 6 instead of Equation 5).
Finally, we computed how our results would change if we considered an older stellar population, by assuming an age of 5 Gyr (instead of the default value of 1 Gyr), i.e. considering a cluster turnoff mass of 1.320 M .
In Table 4, we report the relative variation in the total event rates with respect to the default scenario, for all these cases. The event rates as a function of time (not shown) vary by similar amounts.
Altering the mass of WDs and NSs is inconsequential, except for a very mild increase (decrease) of ∼20-30 per cent in the event rate of the NSs when we increase (decrease) the NS particle mass, and a comparable decrease (increase) of ∼20-30 per cent in the event rate of the WDs when we increase (decrease) the WD particle mass. Increasing (decreasing) the mass of MSs slightly enhances (reduces) the event rates for all species, with changes varying from a few per cent in the case of WDs, to ∼10 per cent for NSs, to ∼30 per cent for MSs and SBHs. Varying the mass of SBHs has the largest effect of all, with positive (negative) changes in the event rates -ranging between ∼15 per cent for SBHs and ∼150 per cent for NSs -when the particle mass decreases (increases).
Varying the definition of the capture radius has the expected result of reducing (enhancing) the event rates when decreasing (increasing) α. We note, however, that the change is at most ∼55 per cent.
Changing the absorption process for WDs has virtually no effect on the other species and only a mild increase in the event rates (∼130 per cent), due to the fact that the absorption radius for TDEs is smaller than that of EMRIs (see Table 3).
number of objects, the mass of the stellar cluster has to vary, although by a little amount (∼30 per cent when varying the MS mass, ∼1 per cent in the WD case, and 1 per cent in the NS and SBH cases).  Table 4. Percentage differences in the total number of events at 1 Gyr, with respect to the default case (whose numbers are listed in the third column of Table 2 Increasing the stellar population age from 1 to 5 Gyr has almost no effect on the event rates, except for the TDEs from WDs, which increase by ∼160 per cent, due to the fact that an older stellar population contains a larger number of WDs (in this comparison, four times larger in number fraction).
Note that the variation in the event rates discussed here is difficult to interpret, owing to the fact that each variation impacts (i) the relaxation rate, (ii) the degree of mass segregation, and (iii) the structure of the models. For this reason, we do not attempt here to enter in detail on this aspect.
What we can clearly see is that, overall, absorption event rates are only mildly sensible to variations in the mass of the stellar-mass objects, to the definition of the capture radius, to the type of absorption process of WDs, and to the assumed stellar population age, as they change at most by a factor of 3.

DISCUSSION
In this paper, we studied the rates of TDEs and GW events about an IMBH that could inhabit the centre of a massive, star-forming clump as the ones that are typically both observed and predicted by simulations at z 3. However, we did not address how the IMBH could have formed or found itself within this system. Although many possible scenarios are possible (see Section 1), here we discuss the most favoured.
A particularly interesting scenario for the formation of IMBHs in the explored systems is via runaway collisions (e.g. Portegies Zwart et al. 2004b): in particular, Portegies Zwart & McMillan (2002b) showed that, if the system has a halfmass relaxation time smaller than ∼25 Myr (i.e. roughly the lifespan of the most massive stars in the cluster), mass segregation could induce stellar collisions that can produce a star whose mass can reach ∼10 −3 cluster masses. This supermassive star is then expected to turn into an IMBH. Note that, in our system, runaway collisions could not occur on the scale of the entire cluster, whose half-mass relaxation time is longer than a Hubble time, but rather in the innermost region (∼0.1 pc), where the relaxation time-scale is 25 Myr in the assumption of an initial Jaffe profile. An alternative possibility is that an IMBH is brought in the clumpy host galaxy via a minor merger; if the main galaxy hosts a large number of massive clumps, as in Tamburello et al. (2015), one of them could capture the wandering IMBH which would then sink in the cluster centre via dynamical friction. Several close interactions between a secondary BH and massive clumps are indeed observed in Tamburello et al. (2017), wherein a central and a secondary BH were added to the galaxies of Tamburello et al. (2015) to study BH pairing time-scales.

Fraction of direct plunges
In Section 3, we reported the rate of gravitational captures due to the central IMBH. However, when a compact object enters the IMBH ISCO, it can do so either almost radially (i.e. via a so-called direct plunge), or via the gradual inspiral induced by GWs (via an EMRI or IMRI). Only in the latter case, the event could be detected by future facilities such as the forthcoming LISA. Although a careful analysis of the EMRI/IMRI and direct-plunge orbits is beyond the scope of this paper, and the GW inspiral is not followed in the current framework, here we propose a simple analysis to obtain at least a rough estimate of the fraction of plunges.
Following Amaro-Seoane (2018) and Bortolas & Mapelli (2019), we can estimate the critical semimajor axis ap above which only plunges can occur as where C = 0.5 (Merritt et al. 2011;Brem et al. 2014) is a constant that defines how much shorter the GW emission time-scale should be, compared to the two-body relaxation time, in order for a potential EMRI/IMRI not to be deflected by two-body relaxation; t rel is the relaxation timescale at ap; and m is the mass of the objects that drive two-body relaxation (i.e. SBHs at the small scales considered for this process). In order to obtain ap, we estimate t rel (E) ∼ D −1 (E), where D(E) is the sum of the angular momentum diffusion coefficients associated to each species in the Fokker-Planck treatment; D −1 (E) is evaluated at ap by recursive iteration. Finally, we compute the fraction of direct plunges as a function of time as the fraction of GW events associated to SBHs or NSs coming from a semimajor axis larger than ap (more precisely, an energy E larger than the energy associated to ap).
The fraction of plunges as a function of time for the three main scenarios explored in this work is shown in Fig the central IMBH has an initial mass of 10 3 M . In all cases, the plunge fraction for NSs is larger than that for SBHs, as the lower mass associated to NSs renders their inspiral timescale longer, making them more susceptible to deflections by two-body relaxation.

Cosmological event rates
Clumpy galaxies represent a significant fraction of the population of UV-bright star-forming galaxies, and their abundance peaks at the same epoch at which the cosmic star formation history reaches its maximum, at around z ∼ 2 (Shibuya et al. 2016). We set out to compute the expected cosmological IMRI rate resulting from SBHs accreting on to a 10 5 M IMBH at the centre of star-forming clumps, assuming that this is present in at least the most massive clumps (∼10 8 M ). Clumpy galaxies such as those simulated by Tamburello et al. (2015) and Tamburello et al. (2017) contain 1-5 of such massive clumps, which is actually a conservative result compared to other published simulations (see, e.g. Mandelker et al. 2014). Conservatively, we will assume just one such massive clump per clumpy galaxy. To compute the cosmological IMRI rate, we count events in the (comoving) cosmological volume between z = 1 and 3. This is because at lower and higher redshift the number of clumpy galaxies drops drastically, whereas it is estimated to be between 40 and 60 per cent of the star-forming galaxy population within such redshift range (Shibuya et al. 2016). Note that clumpy galaxies are present also at lower and higher redshift, hence our estimate will be conservative also on this side. We note, however, that at z > 3 IMRIs occurring on an IMBH of the mass considered here would be hardly detectable by LISA due to low signal-to-noise ratio (Babak et al. 2017). With the chosen volume, our estimate should thus yield a lower limit on the number of such IMRIs that should be detectable by LISA.
In order to compute that, we take into account that, when the IMBH mass is 10 5 M , ∼30 per cent of such events are expected to be accessible to the LISA observatory, the exact number depending on a range of assumptions (see Babak et al. 2017, and in particular their figure 10, for more details on this). Moreover, as discussed in Section 4.1, we need to subtract the direct plunges which, in the default scenario, account for ∼80 per cent of the GW events.
The cosmological LISA-detectable IMRI event rate in clumpy galaxies is thus given by RIMRI,cosmo = nSFf clumpy rIMRIV∆zf det , where nSF is the number density of star-forming galaxies, which we assume to be constant in the chosen redshift range (a reasonable approximation within a factor of two; see, e.g. Brammer et al. 2011), f clumpy is the fraction of clumpy galaxies (which we also assume to be constant: f clumpy = 0.5), V∆z is the comoving cosmological volume between z = 3 and 1 (∼1000 Gpc 3 ; Wright 2006;Planck Collaboration et al. 2020), rIMRI is the IMRI rate per clump computed in the previous sections (4 × 10 −7 yr −1 , see Figure 2), and f det is a suppression factor, which we estimate to be roughly 0.05, that accounts for (i) the detectability of such events by LISA and (ii) the fact that some of these events would be direct plunges. Note that, according to Brammer et al. (2011), the number density of star-forming galaxies varies slightly as a function of luminosity/stellar mass considered, but by less than a factor of two in the mass range in which a disc is massive enough to fragment (4-10 × 10 11 M ; Tamburello et al. 2015), being around ∼2 × 10 −4 Mpc −3 . For the conservative assumption of one star-forming clump per galaxy that hosts a central IMBH, we obtain RIMRI,cosmo ∼ 2 events per yr.
A similar calculation can be performed for TDEs. Since our computations show that the TDE rate is ∼100 times higher than the IMRI rate, we expect an intrinsic rate RTDE,cosmo of a few thousand events per yr in the cosmological volume between z = 1 and 3. Even though current observations of TDEs are all within z 1, future observatories should be able to detect several events also at high redshift. For example, SKA may be able to detect several hundred TDEs per year below z ∼ 2.5 (with a peak rate redshift ∼ 0.5). Using follow-up observations by X-ray facilities such as Athena, the number of identified TDEs could be ∼400 per year, almost up to z = 2 (Donnarumma & Rossi 2015).

Model limitations and caveats
The approach adopted in this paper to address the rates of TDEs and GW events is necessarily approximated in many ways. Here we discuss the main limitations and caveats.
First of all, we opted for a Fokker-Planck approach, which is very fast and computationally cheap. However, its advantages come at a price: in particular, the Fokker-Planck method adopted in PhaseFlow handles a one-dimensional distribution function that only depends on energy, and not on the angular momentum; the dependence of the Fokker-Planck solution on angular momentum in the context of collisional stellar systems has been shown to be weak (Cohn 1980) and, for the IMBH masses we are addressing, the TDE rates have been shown not to get crucially suppressed by, e.g. an initial tangential anisotropy in the velocity space (Lezhnin & Vasiliev 2015). More in general, PhaseFlow can only handle spherically symmetric systems with isotropic velocity distributions. Neglecting a possible net rotation, radial anisotropy in the velocities, and deviations from spherical symmetry likely results in an underestimate of the losscone flux (Holley-Bockelmann & Sigurdsson 2006;Holley-Bockelmann & Khan 2015;Lezhnin & Vasiliev 2016;Stone et al. 2020); in this sense, our estimates can be considered lower limits to the actual event rates. However, one should also consider that the regions near the IMBH are expected to turn into spherical and isotropic owing to two-body relaxation effects, and possibly also to the orbital scattering of stars owing to the presence of the IMBH; both effects occur over a time-scale of the order of the relaxation time (Lake & Norman 1983;Gerhard & Binney 1985;Merritt & Poon 2004;Bortolas et al. 2018).
Furthermore, our treatment assumed an initial lack of mass segregation in the system (i.e. all components shared the same initial density profile); in fact, the central regions of our system relax very quickly in our simulations, and they only retain memory of the initial central density normalization, thus we do not expect this aspect to significantly impact our results. The adoption of four distinct stellar families instead of a smooth and non-discretized mass function is also an approximation, but it is already a significant improvement over the more common approximation of a singlemass family.
An additional caveat concerns the IMBH at the centre of the clump. First of all, we do not address the formation of the IMBH; the occupation fraction of massive clumps by IMBHs is very uncertain; however, our estimate of the event rates detailed in the previous section is conservative on the former aspect (we assumed only one clump hosting an IMBH per galaxy). Clearly, the dynamics of the IMBH cannot be captured in our idealized simulations: in fact, the IMBH Brownian wandering resulting from two-body interactions with the surrounding stars (Lin & Tremaine 1980;Merritt 2005) could influence the event rates in real physical systems, but it was not captured in the Fokker-Planck simulations. Brownian motion is expected to be significant in the system for which we assumed a low IMBH mass (m•0 = 10 3 M ), while it has been shown to be negligible for central objects with masses ∼10 5 M in the context of massive BH binaries 9 (Bortolas et al. 2016).
Connected to this, one should also account for the GW recoil that the IMBH could undergo in response to the anisotropic emission of gravitational radiation following an EMRI or IMRI. This phenomenon is well studied and known to scale with the spins and the mass ratio between the two merging compact objects. The magnitude of the kick associated to small mass ratios is typically very low (Morawski et al. 2018;Baker et al. 2008), much lower than the escape speed of the system considered here.
In our treatment, we could not account for the presence of gas, and we assumed that this component either collapses into stars or is evaporated from the clump by means of stellar feedback. This is likely a good approximation after the first supernovae in the system occurred, thus after a few tens of Myr from the beginning of our integration. 9 Massive black hole binaries experience super-elastic scatterings with the surrounding stars, thus the Brownian wandering in Bortolas et al. (2016) is enhanced compared to single massive BHs (Merritt 2001).
Finally, our treatment of gravity was purely Newtonian, and did not account for general relativistic effects. In addition, we did not take into account resonant relaxation (Rauch & Tremaine 1996), but its effect has been shown to be virtually washed out by its coupling with relativistic precession (Merritt 2015;Alexander 2017), thus this process is not expected to significantly affect the event rates.

CONCLUSION
In this paper, we explored the rates of TDEs and GW events that could occur about IMBHs hosted in massive starforming clumps, typical structures that are found in many galaxies in the high-redshift Universe. Our reference massive clump was modelled so that its proprieties are in good agreement with what found in simulations of clumpy galaxies at z ≈ 1-3 (Tamburello et al. 2015;Tamburello et al. 2017). In fact, clumpy galaxies constitute an abundant fraction (about 50 per cent) of the star-forming galactic population in this range of redshifts (Shibuya et al. 2016).
We found that the TDE rate for MSs is in the range 10 −6 -10 −5 yr −1 per massive clump, depending on the mass of the central IMBH and on the density profile of the clump. We estimated the rate of GW events involving SBHs to be around 10 −7 yr −1 for IMBHs with masses of ∼10 3 -10 5 M ; such rate decreases by an order of magnitude if the initial density profile in the inner regions is assumed to scale as r −1 instead of r −2 . We estimated a relevant fraction (at least ≈20 per cent for SBHs) of GW events to be EMRIs and IMRIs, thus to be in principle detectable by the forthcoming LISA mission.
Finally, we estimated that the forthcoming LISA observatory will be able to detect nearly 2 GW events per yr originating from such star-forming clumps; we also estimated the intrinsic rate of TDEs of MSs in these clumps, obtaining a few 10 3 events per yr; a fraction of these can be observed by forthcoming instruments, and in particular by SKA in combination with Athena.
We highlight that, in the future, the observed event rates in the redshift range z ≈ 1-3 could be adopted to constrain the occurrence of massive clumps hosting IMBHs.
In conclusion, both TDE and EMRI/IMRI detections will be crucial tools to explore the occupation fraction of IMBHs in massive star-forming clumps in the early Universe, and will provide new details on the origin and evolution of massive BHs across the cosmic time.

DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonable request to the corresponding author.