Demographics of Tidal Disruption Events with L-Galaxies

The version of L-Galaxies used in this work, L-GalaxiesBH, is based on the one presented in Henriques et al. (2015) with the modifications included in Izquierdo-Villalba et al. (2019, 2020, 2022) and Spinoso et al. (2023) to incorporate detailed massive black hole physics.

In brief, with respect to the model presented in Henriques et al. (2015), this new version adds a detailed model for the assembly (mass and spin) of nuclear MBHs, the dynamical evolution of MBH binaries, and the production of wandering MBHs (see Izquierdo-Villalba et al. 2020, 2022). Concerning the genesis of the first MBHs, L-GalaxiesBH models the spatial variations of the intergalactic metallicity and the Lyman-Werner background²²2The Lyman-Werner (LW) band is a specific energy-interval (i.e. $h\nu=[11.2\,-\,13.6]{\rm eV}$) in the UV spectrum. LW photons are responsible for the photo-dissociation of molecular hydrogen (e.g. Haiman et al. 1997). produced by star formation events to account for the formation of heavy (i.e. $M_{\rm seed}\,{\sim}\,10^{5}\,\mathrm{M_{\odot}}$) and intermediate-mass ($M_{\rm seed}\,{\sim}\,10^{3}\,{-}\,10^{4}\,\mathrm{M_{\odot}}$) MBH-seeds, respectively via the collapse of pristine massive gas clouds and runaway stellar mergers within dense high-redshift³³3not to be confused with the ones introduced later in this work, for which we follow a different treatment NSCs. The formation of light seeds ($M_{\rm seed}\,{\sim}\,10^{2}\,\mathrm{M_{\odot}}$) after the explosion of the first metal-free stars (PopIII) is accounted for by grafting/inheriting the evolved counterparts of light-seed modeled self-consistently by the GametQSO/dust model (see e.g. Valiante et al. 2016). We note that concerning the black hole seeding model presented in Spinoso et al. (2023), we adopt a slightly higher amplitude of the “grafting probability”, setting the parameter $\rm G_{p}=0.25$ (see Spinoso et al. 2023, for the definition of this parameter). This choice is motivated by the normalization of the $z\,{=}\,0$ black hole mass function in the current work.

Once the first MBH seeds are formed, galaxy mergers and disk instabilities funnel new gas to the galactic nuclei, making it available for the growth of nuclear MBHs. The gas reaching the center is progressively consumed by the MBH first in an Eddington limited phase, followed by a sub-Eddington one (Bonoli et al. 2009; Izquierdo-Villalba et al. 2020). Episodes of gas accretion, on top of triggering BH growth, also modify its spin (Izquierdo-Villalba et al. 2020). The dynamical evolution of massive black hole binaries in a post-merger galaxy is two-fold (Begelman et al. 1980). During the pairing phase, the MBH(s) from the satellite galaxy migrate towards the galactic center via dynamical friction (following Binney & Tremaine 1987), forming a hard binary upon arrival at the nucleus. Then the hardening phase, where interactions with a circumbinary disc (gas-torque model by Dotti et al. 2015 and preferential growth as in Duffell et al. 2020) or the surrounding stars (following the model by Quinlan & Hernquist 1997; Sesana & Khan 2015) assist the gravitational-wave-driven evolution, along potential triplet interactions in the cases that a third MBHs comes in during the inspiral (modeled based on the results of Bonetti et al. 2018). The latter, along with gravitational recoil after merging of MBHs (as described in Lousto et al. 2012), can result in wandering MBHs. While L-GalaxiesBH can track the evolution of wandering MBHs, we do not include that in this work for simplicity (see the discussion on TDEs from wandering MBHs in Sect. 5).

As an example of the capability of L-GalaxiesBH to produce a realistic population of MBHs, in Fig. 1 we show the MBH mass function $\phi(M_{\bullet})$, where $M_{\bullet}$ is the black hole mass (throughout this work), and the MBH median spin distribution $\tilde{\chi}_{\bullet}(M_{\bullet})$ at $z\,{=}\,0$ and $z\,{=}\,5.5$ for the version of L-GalaxiesBH adopted in the current work. We compare our results with available data. Regarding the black hole mass function, as noted by Shankar et al. (2019), all observationally-derived values seem to converge at the high mass-end ($M_{\bullet}\,{>}\,10^{7.5}\,\mathrm{M_{\odot}}$ Merloni & Heinz 2008; Cao 2010; Gallo & Sesana 2019; Shankar et al. 2009, 2013; Mutlu-Pakdil et al. 2016; Vika et al. 2009; Aversa et al. 2015). However, constraints in the intermediate-mass range $M_{\bullet}\,{\sim}\,10^{5}-10^{6}\,\mathrm{M_{\odot}}$ are much less stringent. L-GalaxiesBH agrees with the most optimistic estimates from Shankar et al. (2009) and over-estimates with respect to all other available constraints. Regarding the spin distribution, the model agrees with the constraints from Reynolds (2021) at high masses, yet it predicts high-to-maximal spin for MBHs in the intermediate mass range, where there are no observational constraints. TDEs can potentially offer as a new probe for both the MBH mass and spin distribution in this range (see discussion in Sect. 4.3).

As mentioned above, the novelty of this work is the inclusion of TDEs within a full galaxy evolution context. To encompass this ambitious task it is necessary to model the stellar environment in which MBHs are embedded, on top of their formation and evolution. In this section, we describe how disks, bulges, and NSCs are included in L-GalaxiesBH. Together with black hole masses, the properties of the nuclear stellar environment will be the input for predicting time-evolving TDE rates with PhaseFlow.

Bulges in L-GalaxiesBH grow after galaxy interactions (major and minor mergers) and disk instabilities occurring in massive stellar disks. The specific properties of these events fully determine the final mass and extension (i.e. effective size) of the bulge. We stress that the scale length for disks $R_{\rm eff,d}$ and the effective size of bulges $R_{\rm eff,b}$ (equivalent to the scale radius of a Sérsic profile) are computed self-consistently inside L-GalaxiesBH by tracing the spin evolution of the galaxy components and applying energy conservation during mergers and disk instabilities. The profile of each bulge is assumed to follow a Sérsic model, whose steepness (i.e. Sérsic index) is associated with each bulge by using the observational results of Gadotti (2009), approximately a Gaussian distribution peaking at Sérsic index $n_{s}\,=\,3-4$, as implemented in Izquierdo-Villalba et al. (2019). As already mentioned, galactic disks arise as a consequence of gas cooling and star formation events occurring in the center or dark matter halos. Together with galaxy encounters, these events determine the extent of the disk radial profile. Taking this into account, in this work we assume that pure disks with no bulge are characterized by the Sérsic index $n_{s}\,{=}\,1$, regardless of redshift, and a scale radius of $R_{\rm gal}\,{=}\,1.68\,R_{\rm eff,d}$. For the rest of this work, we refer to the sum of the disk and bulge mass as galaxy stellar mass $M_{\ast}$ (the halo and intracluster stellar mass are neglected during our analysis), while $M_{\rm gal}$ is reserved for the integrated mass of a Sérsic profile (either a bulge or a disk) with radius $R_{\rm gal}$. As we will see later, TDE events due to encounters between the nuclear MBHs and stars belonging to the bulge or disk component will be calculated assuming that these density profiles extend all the way to the center of the galaxy.

NSCs observed in the centers of a great fraction of dwarf and massive galaxies in the local Universe are the densest stellar structures known (see Neumayer et al. 2020, for a review). This inevitably suggests that they might be an ideal nursery for TDEs. In the following paragraphs, we describe a simple phenomenological model that we are introducing in L-GalaxiesBH to incorporate NSCs in galaxies (“nucleation”). An extensive and self-consistent model of the birth and evolution of NSCs will be presented in a future paper (Hoyer et al. in prep).

In our implementation, we assume that the stellar system surrounding the MBH is composed by a bi-chromatic population of stars made up of main-sequence stars of 0.38 M_⊙ (encompassing $\approx 93$% of the total stellar mass) and $16\,\mathrm{M_{\odot}}$ stellar black holes⁵⁵5The choice for these values for the mass of main sequence stars and stellar black holes comes from the fact that those are close to the average masses for those objects assuming an evolved Kroupa (2001) stellar mass function; in addition, the second moment of the mass function, which sets the relaxation rate of the system and thus the TDE rates, attains a value which is very close to the actual one if we were to assume a complete and evolved stellar mass function (see e.g. Kochanek 2016; Stone & Metzger 2016; Pfister et al. 2022; Bortolas 2022). (encompassing the remaining $\approx$7%). Stars are considered to be destroyed if their separation to the MBH gets below $r_{\rm t}$ (Eq. 3), where we used $r_{\star}=0.44R_{\odot}$, which is the expected radius of a 0.38 $\rm M_{\odot}$ star (see e.g. Bortolas 2022, Eq. 3, for more details). Once accretion occurs, 30% of the stellar mass⁶⁶6At a full disruption, the MBH captures 50% of the stellar mass. We selected a subsequent mass loss due to radiation and winds removing 40% of the bound material based on results from observations (Mockler & Ramirez-Ruiz 2021) and simulations (Bu et al. 2023). Simulations including radiative transfer from Steinberg & Stone (2024) yield a lower percentage of 15%. is added to the MBH, while the remainder is assumed to be lost in radiation and the interstellar medium. Stellar black holes are instead captured by the MBH if they get closer than $8GM_{\bullet}/c^{2}$ from the MBH, where $c$ is the speed of light in vacuum and $G$ the gravitational constant. During the evolution of the system, the stellar populations undergo the traditional dynamical phenomena expected in the vicinity to an MBH, i.e. they develop a Bahcall & Wolf (1976) cusp, stellar black holes segregate in the center dominating relaxation in the closest vicinity to the MBH and finally, once the system reaches a dynamical equilibrium, expands and lowers its TDE rates as a result of dynamical heating. All these phenomena are captured by PhaseFlow.

Given the fast performance of the code, we have generated multi-dimensional tables spanning a range of initial central MBH masses and a range of host-environment properties (mass, scale radius, compactness for disks, bulges, and NSCs), encompassing all values predicted by L-GalaxiesBH, as described in what follows. In particular, all runs in the multi-dimensional grid were initially evolved for a Hubble time using 300 bins in the phase-volume, the variable used in the code to parameterize the distribution function ⁷⁷7The phase-volume $h(E)$ is the volume in phase space spanned by all the orbits with energy $E$; unlike the orbital energy it is invariant under adiabatic changes of the potential, like in the case of a central MBH accreting via TDEs, and this makes $h(E)$ a very good choice for this problem..

We model rates depending on their reservoir origin, which could be either a bulge/disk or an NSC. Rates are saved in a large multi-dimensional grid, which includes the parameter space of the MBH, the stellar environment properties, and the time dimension. Rates are not static but instead vary with time (not necessarily a steady state is reached). This adds substantial realism to our computation, as many systems dominating the overall TDE rate are often characterized by a very large initial rate which drops by orders of magnitude with time. Assuming the static rate (which would coincide with the $t=0$ rate) thus results in a non-negligible overestimation of TDE rates, especially for the case of NSCs.

Once the runs with PhaseFlow are completed and the grid in the parameter space has been fully mapped, we store the resulting event rates, $\Gamma$, which are a function of the aforementioned parameters and time $t$. For the bulge/disk contribution, the rates are encapsulated in the function $\Gamma_{\rm gal}(M_{\rm gal},R_{\rm gal},n_{s};M_{\bullet},t)$, while for the NSCs, the rates are given by the function $\Gamma_{\rm NSC}(M_{\rm NSC};M_{\bullet},t)$. We store the rates differently for late times:

and for early times

with 50 and 60 evenly spaced logarithmic bins, respectively for the two time ranges. We consider separately TDE rates at times below $10^{7}{\rm yr}$ given the time resolution of L-GalaxiesBH which spans between $dt_{\rm step}\rm\sim 5\,{-}\,20\,Myr$, depending on redshift. Specifically, for MBH mass $M_{\bullet}<10^{6}\,\mathrm{M_{\odot}}$ the peak TDE rate due to NSCs falls below the time resolution of L-GalaxiesBH $dt_{\rm step}\sim 10$Myr, therefore the initial phase of high TDE rates and growth (dubbed here as prompt phase) would not be treated properly by our model (see Fig. 3). To resolve that, when the TDE phase due to NSCs is initialized for these small MBHs, we draw a random time between $1$ and $dt_{\rm step}$ and assign rates from the prompt phase tables according to this time. At the same time, we add the integrated mass during the unresolved prompt phase to the MBH.

In the following section, we outline how we merge L-GalaxiesBH with the $\Gamma_{\rm gal}$ and $\Gamma_{\rm NSC}$ rates constructed from the multi-dimensional grid described above. To guide the reader, Fig. 3 shows several examples of the time-evolving TDE rates for a variety of black hole masses and for a fixed set of parameters. In the figure, for each black hole mass, we assume that the corresponding galaxy stellar mass is given from the relation $\log_{10}\,(M_{\bullet}/10^{7.43}{\rm M}_{\odot})=0.62\log_{10}(M_{\rm gal}/10^{10.5}{\rm M}_{\odot})$, which is the best-fit relation we get in L-GalaxiesBH at $z\,{=}\,0$ for MBHs up to ${\sim}\,10^{8}\,\mathrm{M_{\odot}}$. The scale radius assigned to the Sérsic profile is here assumed to be an average value of $R_{\rm gal}\,{=}\,R_{\rm scl}$. The disk corresponds to $n_{s}\,{=}\,1$ and the bulge to $n_{s}\,{=}\,4$. The matching NSC profiles are created by using the scaling relations presented in Eq. 1 and Eq. 8.

As shown from Fig. 3, the main general feature of our model is that the typical, average rates increase with increasing MBH mass, as the stellar reservoirs assigned are increasingly more massive for more massive MBHs. However, there is a strong differentiation between the galaxy component and the NSC component; while bulges and disks profiles often obtain a constant event rate given the large relaxation timescale of these systems, NSCs matching to $M_{\bullet}\,{<}\,10^{6}\,\mathrm{M_{\odot}}$ show decay in their initial TDE rate due to their fast relaxation dynamics. Typically, the decay happens earlier for smaller MBHs and smaller reservoirs. Notably, NSCs with $M_{\bullet}$ from $10^{4}\,\mathrm{M_{\odot}}$ to $10^{6}\,\mathrm{M_{\odot}}$ start with similar TDE rates at $10^{7}\,\mathrm{yrs}$, but maintain it shorter times compared to the Hubble time, as the systems quickly reach their equilibrium in the form of a Bahcall & Wolf (1976) cusp and subsequently expand due to the dynamical heating resulting from efficient relaxation. That proves the importance of including the analysis of the initial rates that fall below L-GalaxiesBH time resolution.

As soon as an MBH forms inside a galaxy in L-GalaxiesBH, we set a TDE rate depending on the existence of a stellar bulge/disk and/or an NSC. Effectively, we look up the multi-dimensional grid generated with PhaseFlow for the rate $\Gamma_{i}$, with $i\,{=}\,{\rm gal}$ for the galaxy component and $i\,{=}\,{\rm NSC}$ for the NSC. From this, we derive the time-evolving rates for each MBH and stellar environment. For galaxy and NSC properties we look up for the closest values in the grid, while we interpolate in the MBH and time dimension.⁸⁸8In only a few instances values predicted by L-GalaxiesBH are outside the parameter space covered by the grid. In these rare cases, we use the closest valid value.

When TDEs start being produced within a galaxy, either from the relaxation of the bulge/disk or the NSC, we start a clock, $\Delta t_{\rm TD,i}$. In the successive steps of the L-GalaxiesBH run, the PhaseFlow grid is checked at the corresponding subsequent times. As described in the next section, we reset this clock to zero only when the MBH and/or the galaxy have substantially changed, which can happen, for example, after mergers.

Provided the TDE rate and the accretion fraction $f_{\ast,\,\mathrm{TDE}}$ per event, the mass accreted by the central MBH from each stellar component¹¹¹¹11Note that we consider an NSC as a decomposed part of the central bulge/disk, that we subtract separately from. during each time-step of L-GalaxiesBH is given by:

and the mass subtracted ($i=\rm gal$ or $i=\rm NSC$) is:

Notice that, if both a bulge/disk and an NSC are present, we sum the two TDE rate contributions. The results of the mass growth of MBHs through this mechanism are discussed briefly in Sect. 4.1.

As mentioned above, PhaseFlow is able to capture the time evolution of the stellar density profile as a result of relaxation and MBH growth due to star accretion. Thus, once TDEs start taking place within a galaxy in L-GalaxiesBH, we follow their time evolution based on the grid provided by PhaseFlow.

However, dramatic events, such as mergers, disk instabilities, and starbursts, can lead to major changes in the stellar environment surrounding the MBH. In these cases, we reset the clock for TDE rates ($\Delta t_{\rm TD,i}$), and the new, unrelaxed stellar system starts again evolving towards a new relaxed state until the next violent dynamical event.

The conditions under which we consider it relevant to check the state of the TDE process taking place around an MBH are the following ones:

Changes in MBH mass:

If the mass accreted by the MBH via gas accretion is three times larger than the mass accumulated from TDEs, then the $\Delta t_{\rm TD,NSC}$ clock is reset again and the NSC mass and properties are adapted to the galaxy stellar mass at that time¹²¹²12This condition guarantees the MBH does not move by more than three mass bins away from its initial mass on the PhaseFlow grid, at which point the assigned TDE rates we expect to be significantly different. For bulges and disks, rates are less time-dependent (see discussion in Sec. 2.2), so we omit this condition.

Changes in the galaxy component mass:

If during one L-GalaxiesBH time-step the bulge or disk stellar mass increases by more than 20%, the clock $\Delta t_{\rm TD,gal}$ is reset. Such variations in stellar mass in a short amount of time is a typical change in post-starburst galaxies (Kaviraj et al. 2007; van Velzen 2018; Wild et al. 2016, 2020). Note that events that do not modify the total stellar mass can also lead to a reset of the rates (e.g., the transfer of mass from the disk to the bulge during disk instabilities).

Changes in NSC mass:

We do not expect all changes of the galaxy component to affect the TDE rate evolution in NSCs, since these can continue relaxing at their own pace unaffected by the changes at larger galactic scales. However, as described in Sect. 2.1, an NSC can undergo significant mass changes following the host galaxy changes. Following the same approach described above, we select that NSC that grow more than ${>}\,$20% within one L-GalaxiesBH time-step will satisfy the condition of resetting the clock $\Delta t_{\rm TD,NSC}$. This threshold is selected to be equal to that of the galaxy component in order to minimize the free parameters of the model.

Changes in NSC to MBH mass ratio:

While there are no clear constraints on the limits of the mass ratio between an NSC and the MBH at its center, namely $\mathcal{D}\,=\,M_{\rm NSC}/M_{\bullet}$, we decided to put a lower limit on this ratio so that an NSC gets to be (re)generated only if $\mathcal{D}\,{>}\,\mathcal{D}_{\rm nuc}$. Only few physically-motivated arguments can provide some insight for the allowed value of $\mathcal{D}$: a) MBHs cannot be arbitrarily massive if they form in the cluster, e.g. simulations of runaway stellar collisions by Kritos et al. (2023) indicate that $\mathcal{D}\,>\,1$ and b) the NSC destruction from binary MBHs e.g. total and partial disruption of clusters from intermediate-mass MBH binaries takes place for values $\mathcal{D}\,{<}\,5\;{\rm and}\;<15$ according to Khan & Holley-Bockelmann (2021). Given the theoretical uncertainties of the NSC-MBH symbiosis and provided local NSCs are observed with even lower values ($\mathcal{D}\,\sim\,0.01$ Neumayer & Walcher 2012), we select conservatively $\mathcal{D}_{\rm nuc}\,=\,0.1$. We stress that this parameter (along with $M_{\ast,\rm NSCcut-off}$) can be later substituted from a physically motivated formulation when the NSC model is implemented in L-GalaxiesBH by Hoyer et al. (in prep).

We have checked that variations within an order of magnitude of the thresholds mentioned above do not significantly affect the results on the volumetric rates presented in this work. However, we will further discuss some of these assumptions in Sect. 4.4.

When MBHs experience phases of gas accretion, they shine as active galactic nuclei (AGN). These phases can be coeval with events of stars disruptions. The two processes may not be completely unrelated to each other. For instance, it has been proposed that TDE rates can be enhanced due to the alignment of retrograde orbits of NSC stars with the MBH accretion disk (Generozov & Perets 2023; Nasim et al. 2023) or due to the turn-off of the disk (Wang et al. 2024) or due to the quadrupole moment of the disk modifying the loss cone(Kaur & Stone 2024). However, by construction, our model does not account for such possible rate enhancement since gas physics are not included in PhaseFlow. Regarding the impact of TDEs on the AGN disk, the few TDE interpretations of fast flares in AGN light curves suggest a (temporal) post-flare increase (Blanchard et al. 2017; Liu et al. 2020) and others a decrease (Cannizzaro et al. 2022; Cao et al. 2023) of the AGN luminosity. Given the high uncertainty, we assume that the disk recovers quickly (compared to the time to the next TDE) or is not interrupted at all by TDEs so that both the TDE and AGN luminosity are kept as independent properties of each MBH. We do not expect this assumption to affect the global rates, due to the small AGN fraction at $z\,{\sim}\,0$ (${<}$10% of MBHs are active) unless the TDE rates are significantly boosted in the presence of an accretion disk (by a factor of 10 or greater).

Once TDEs are carefully linked to the MBHs and their galaxies evolving within L-GalaxiesBH, we need to transform the individual rates to observable events to finally compare our predictions with the events picked up by time-domain surveys. Here we keep the following minimum number of assumptions to translate simulated TDE rates to observed ones:

• $\bullet$

  We assume that all TDE events are full disruptions. In reality, for massive stars, only those with a pericenter that is a fraction of the tidal radius are fully disrupted (otherwise they are partial), and the energy output of the events may depend on the penetration factor of the star. Moreover, provided that a handful of partial disruption events have been identified (Payne et al. 2021; Malyali et al. 2023), we discuss the impact of such a choice in Sect. 5.

• $\bullet$

  We assume that all TDEs give a luminosity that falls within the sensitivity of current surveys, i.e., all events are detectable, if there is no AGN contribution. If there is no AGN contribution. In the presence of an AGN, we ignore TDEs in MBHs with AGN bolometric luminosity $L_{\rm bol}\,{>}\,10^{42}$erg/s (fiducial model) unless noted otherwise. Note also that we are not including any obscuration and line-of-sight effects.

• $\bullet$

  Regarding the event horizon suppression, the predicted TDE rates for each MBH are considered either visible or not depending on the Hills mass of the black hole itself. The Hills mass depends on (i) the mass $M_{\bullet}$ and spin $\chi_{\bullet}$ of the MBH, (ii) the age of the loss-cone since the cluster/galaxy component was last reset $\Delta t_{\rm TD,i}$ and (iii) the mass of the infalling star. For this last point, we draw the stellar mass $m_{\ast}$ of the disrupted star from a truncated Kroupa mass function ($m_{\ast}\,{<}\,1.5\,\mathrm{M_{\odot}}$). Given the values above, each MBH is assigned a Hills mass $M_{Hills}(m_{\ast},\Delta t_{\rm TD,i},\chi_{\bullet})$ based on the tabulated results of Huang & Lu (2024, Table D1) who have performed a series of simulations of solar-metallicity stars disrupted by MBHs¹³¹³13These authors may have used slightly different environmental set-ups, but we expect the values to be quantitatively close.. If $M_{\bullet}\,{>}\,M_{Hills}$, we omit this MBH from the sum of the bulk rates. Otherwise, we assume that all events are full TDEs and have the rate assigned to the MBH.¹⁴¹⁴14From theory, e.g. MacLeod et al. (2012) and Kochanek (2016), for a given $M_{\bullet}$, TDE rates scale with the initial mass function $dN\,{/}dm_{\ast}$ and power-law $m_{\ast}^{1/6}$. This will matter mostly beyond the maximum Hills mass for $m_{\ast}\,{=}\,0.38\,\mathrm{M_{\odot}}$ ($M_{\bullet}\,{>}\,10^{8}\,\mathrm{M_{\odot}}$), a mass range that TDE rates are shown only indicatively. Although the dependence on the initial mass function is captured with the sampling of stars, we do not rescale rates with this power-law dependence for massive monochromatic stars when adding them to the bulk rates.

Note also that, given the mass resolution of the dark matter simulation used and of the multidimensional grid of TDE rates, we omit from the analysis the TDE rates from galaxies with masses $M_{\ast}\,{<}\,10^{5.5}\,\mathrm{M_{\odot}}$ and/or MBHs $M_{\bullet}\,{<}\,10^{2.5}\,\mathrm{M_{\odot}}$. Finally, note that given all the assumptions above, the volumetric rates presented below can be considered optimistic.

Having captured the $z\,{=}\,0$ global rates, we now discuss their redshift evolution and how that compares with the evolution of the underlying MBH properties and their environment.

In Fig. 7 we display the volumetric TDE rates for redshifts $z\,=\,0,1,2$. Regarding the case of TDE rates as a function of MBH mass, there does not seem to be a significant change in the rates between $z\,{=}\,0$ and $z\,{=}\,1$. This is due to the very mild evolution evolution predicted by the model for $z\,{<}\,1$. This result is in line with the general assumptions that the volumetric TDE rates do not evolve significantly with redshift (see e.g.van Velzen & Farrar 2014; van Velzen 2018; Yao et al. 2023 see, however, the approach of Kochanek 2016). Despite this, there is a slight decrease of half an order of magnitude when comparing resolved TDEs at $z\,{=}\,0$, $1$ and $z\,{=}\,2$, which will be an important feature to study with deep sky surveys such as LSST (Hambleton et al. 2023; Bricman & Gomboc 2020; Bučar Bricman et al. 2023) and UV transient surveys like QUVIK (Zajaček et al. 2024). However, it is worth mentioning that our predictions are model-dependent. We assume that MBHs are always associated to an NSC at all redshifts (see Sect. 2.1). If this was not true at early times, high-redshift TDE rates would be lower.

The volumetric rates shown in Fig. 7 can be divided into contributions from the galaxy component (“gal”) and the NSC. The first, although it steadily increases from $z\,{=}\,2$ to $z\,{=}\,0$ with the increase of bulges in galaxies, still provides a negligible contribution to the volumetric TDE rates at all masses and all redshifts. In practice, the volumetric rates are entirely due to galaxies hosting NSCs (this contribution is not shown in the figure as it essentially coincides with the total “fiducial” volumetric TDE rates). While the per-galaxy rates for NSCs remain the same at high masses towards higher redshifts, there is a slight enhancement (factor less than three) compared to $z\,{=}\,0$ for $M_{\bullet}\,<\,10^{6}\,\mathrm{M_{\odot}}$ MBH and $M_{\ast}\,<\,10^{9}\,\mathrm{M_{\odot}}$ host-galaxy stellar masses, because of continuous reset of TDE rates during frequent mergers at cosmic noon. This somewhat counteracts the decrease of the number density of inactive MBHs towards higher redshifts and results in a moderate evolution of the volumetric TDE rates.