TDEs on FIRE: Illuminating the Cosmic Evolution of Tidal Disruption Rates
Abstract
Tidal disruption events have been extensively studied in the local universe, but their prevalence at high redshifts remains largely unexplored. Using the FIRE-2 cosmological zoom-in simulations, we compute the per-galaxy tidal disruption rate (TDR) over –, covering black holes from IMBHs to SMBHs. The averaged TDR rises from the early universe, peaks at near , and declines to at . The TDR correlates strongly with host galaxy star formation rate and central stellar density at all redshifts. Qualitatively, the TDR trends with the and persist from high redshift to the local universe, suggesting similar BH-galaxy scaling across cosmic time. Satellite galaxies exhibit comparably high TDRs, with their fractional contribution increasing significantly at high redshifts, highlighting their potential for probing IMBHs and early galaxy assembly. This work demonstrates that cosmological simulations offer a promising avenue for constraining the cosmic evolution of the TDR, paving the way for future comparisons with next-generation observations.
show]rudrani.chowdhury@tifr.res.in
show]lixindai@hku.hk
I Introduction
Two-body scattering at a galactic center can send a star onto an orbit that brings it too close to a massive black hole (MBH) [41]. If the star passes within the tidal radius of an MBH, tidal forces tear it apart. Roughly half of the debris falls back onto the MBH, producing a luminous flare that peaks in weeks and fades over months [50, 16]. This type of transient is called a tidal disruption event (TDE). The properties of a TDE reveal critical information about both the MBH and the disrupted star [21, 38, 43]. Recent studies have further shown that TDEs are powerful probes of extreme accretion and outflow physics around MBHs [13], intermediate-mass black holes [10, IMBHs,], and the first stars in the universe [29]. In addition, statistical TDE samples and event rates reveal large-scale host galaxy properties [58, 48, 19, 9].
The redshift evolution of TDE rates, however, remains almost entirely unconstrained, both observationally and theoretically. While observing TDEs at would reveal black hole growth, demographics, and co-evolution with host galaxies, the current TDE sample is overwhelmingly local because surveys like ZTF are flux-limited. Only a handful of jetted TDE candidates exist at and their detection relies on large beamed luminosities [8, 5, 1]. Theoretical estimates face even deeper uncertainties. The TDE rate prediction sensitively depends on the galaxy stellar profiles, which remains highly uncertain for high-redshift galaxies. Even recent JWST observations [12] have not resolved tensions with the theory [36]. The TDE rate is also highly sensitive to the black hole mass function (BHMF), yet low-mass MBHs, which dominate the TDE population, are poorly constrained at high redshift [34, 60]. As a result, the cosmic evolution of TDE rates remains without any robust prediction or measurement.
A recent breakthrough offers a glimpse forward: [32] detected a TDE candidate at in the JWST COSMOS-Web survey [7]. This discovery opens the door to finding more early-universe TDEs. The number of such events is expected to grow substantially with upcoming high-cadence, wide-area facilities like eROSITA [53], the Einstein Probe [68], the Vera Rubin Observatory [27], and future missions such as ULTRASAT [52] and the Nancy Grace Roman Space Telescope [55]. These same missions will open the door to novel TDE classes, such as off-nuclear (ON) TDEs from IMBHs in dwarf galaxies and globular clusters (GCs), and TDEs from SMBH binaries [39, 3, 28, 67, 45], each offering new windows into galaxy assembly and merger histories [51].
It is therefore timely to investigate the TDE rates and their dependence on host galaxy properties across cosmic time, spanning black holes from the centers of massive galaxies to the small satellite scales. Previous zoom-in simulations studied TDEs at very high redshifts () [46, 47, 37], but only for small galaxy samples and at redshifts beyond the detectability of most current surveys. The second generation of the Feedback In Realistic Environments (FIRE-2) simulation [65] is ideally suited for this study, offering superior resolution and comprehensive modeling of star formation, evolution, and stellar feedback compared to other existing simulations (see Section II.2 for details). Our goal is to cover the redshift range with better-resolved stellar structures and a statistically larger sample than previous work, and to use this framework to predict TDE rates and their detectability with upcoming surveys.
In this paper, we have used the Massive Halo simulation suite of FIRE-2 available across a redshift range of . The novelty of the Massive Halo suite lies in the inclusion of the sophisticated model of the formation and growth of MBHs in the simulation [2], which is crucial for the purpose of this paper. Additionally, this simulation suite offers the unprecedented spatial resolution of FIRE-2 that reaches a few parsec scale, careful and comprehensive modeling of star formation, along with the evolution and stellar feedback physics. These are discussed in Section II.2, following the introduction of the framework for theoretical estimation of TDE rates in Section II.1. We then describe our methodology of combining the two and obtaining the TDE rate in Section II.3. Afterwards, we describe the findings from this study in Section III, including the redshift evolution of TDE rates (Section III.1), their correlation with the star formation rate (Section III.2), host galaxy and MBH masses (Section III.3) and the stellar properties (Section III.4). TDE rates in the satellite galaxies are discussed in Section III.5. Finally, in Section IV, we summarise the results, discuss future research directions based on the foundation of this paper, and certain limitations of the current work.
II Methodology
We provide a general background of calculating the TDE rates assuming two-body scattering between the stars in the following section. A basic overview of the TDE rate calculation using loss cone theory is discussed in Section II.1. We introduce FIRE-2 simulation in II.2, with a brief discussion of the important parameters of the Massive Halo suite and black hole model given in Section II.2.1 and Section II.2.2, respectively. Finally, calculation of the intrinsic TDE rates in the simulated FIRE-2 galaxies is discussed in Section II.3.
II.1 TDE Rate from Loss Cone Theory
TDE rates are calculated based on the loss cone dynamics, which we briefly describe here and refer to [57] for a more complete discussion. The loss cone can be understood as a region in phase-space through which stars diffuse into orbits close to the MBH and become disrupted. The distance to the MBH at which stars get disrupted is called the tidal radius and can be approximated as
| (1) |
where is the mass of the BH, and are the mass and radius of the disrupted star, respectively. It is common practice to define a variable, called penetration parameter (), such that
| (2) |
where is the pericenter of the stellar orbit. Partial, mild, or full TDEs are defined when , , , respectively. Correspondingly, the specific angular momentum of the stellar orbit at can be written as
| (3) |
Hence, stars with will get disrupted within . The phase space corresponding to is called the loss cone.
The number of stars per unit specific energy and scattered into the loss cone due to collisional two-body relaxation is estimated as [59, 48, 10]:
| (4) |
where
| (5) | ||||
| (6) | ||||
Hence, the total consumption rate of stars into the loss cone gives the TDE rate
| (7) |
Different terms in Equations 5 and 6 are as follows:
- •
: Stellar distribution function, which is the number density of stars with specific energy in the phase-space. This can be obtained from the density profile of the stars () and total gravitational potential () through Eddington’s formula [4]
(8) - •
: Loss cone filling factor that demarcates between the full () and empty () loss cone regimes and is given by
(9) where is the orbit-averaged angular momentum diffusion coefficient, is the orbital period of stars and
(10) being the angular momentum of a circular orbit.
- •
: The Bessel functions of zeroth and first order, respectively. is the zero of .
Differential TDE rate (Equation 4) reaches a maximum near a critical energy, , producing the highest flux into the loss cone at the energy when . This energy corresponds to a critical distance from the BH, denoted as [58]. For BHs with mass , this scale coincides with the radius of influence, , defined as the radius at which total enclosed stellar mass is equal to the mass of the BH. Hence, resolving the stellar density profile at becomes critical for a proper estimation of TDE rate in the loss-cone framework. In the following sections, we provide a description of the FIRE-2 simulation and TDE rates calculated using stellar distribution in FIRE-2 galaxies.
II.2 FIRE-2 Simulation
We briefly summarize below the components of FIRE-2 relevant for this paper and refer to [65] for a more detailed description of the simulation suites and discussion of their results.
II.2.1 FIRE-2 Massive Halo Suite Parameters
In this work we use four simulations; A1, A2, A4 and A8, collectively known as the Massive Halo suite in the public data release of FIRE-2.11 1 We note that the second data release of FIRE-2 contains more galaxies, a larger redshift range with more snapshots [66]. However, entire analysis of the paper was completed before the public release of DR2 FIRE-2 and we plan to use this bigger sample in future studies. The Massive Halo suite is re-simulated from the MassiveFIRE galaxy simulation evolved with the FIRE-2 code [17]. FIRE-2 is a cosmological zoom-in simulation, evolved with the -bodyhydrodynamics code GIZMO [25]. The simulations are developed to model galaxy formation in the cosmological simulation framework aiming to resolve multiphase interstellar medium by implementing realistic stellar evolution and feedback mechanisms. FIRE-2 includes radiative cooling, heating, low energy ionizing cosmic rays, dust, star formation, and stellar feedback in terms of stellar wind, core collapse, and Type Ia supernovae. We refer to [24] for a detailed description of the simulation models, associated parameters, star formation and stellar feedback implemented in the FIRE-2.
All four simulations in Massive Halo can achieve the mass resolution of and for the baryonic and dark matter (DM) particles, respectively. The force softening length of DM, black hole, star and gas particles are pc, pc, pc, and pc respectively. Adaptive force softening is used for gas particles with a minimum length of , whereas, and are kept fixed below . Such high resolution achieved in FIRE-2 makes it exceptional among the existing cosmological simulations. This outstanding resolution is crucial for determining the stellar distribution at the galaxy centers to accurately estimate their TDE rates. Total stellar mass at reaches , and in A1, A2, A4 and A8, respectively. We briefly discuss below the seeding criteria and evolution of black holes in the Massive Halo and refer the reader to [2] for detailed discussions.
II.2.2 Formation and Growth of Black Holes in FIRE-2
Black holes are considered to be collisionless sink particles. Seed BHs of mass are formed using the friends of friends algorithm in the DM halo by converting the most bound star particle when the stellar mass of a DM halo exceeds [14]. Such seed BHs grow afterwards through accretion and merger [56]. We note that feedback from the BH is not incorporated in the Massive Halo suite, although the accretion rate is adjusted to match the normalization of the observed relation in the local universe. The impact of the absence of active galactic nuclei (AGN) feedback on the results of this paper is discussed in Section IV.
Furthermore, the resolution of the cosmological simulation limits the ability to fully resolve the gravitational dynamics of the black holes. While in reality, BHs have much higher mass than the surrounding gas, dark matter and individual star particles, this might not always be the case in the simulation due to numerical artifacts. A seed BH could sometimes form with a lower mass compared to the other particles, which might prevent it from sinking to the gravitational potential of the galaxies or even cause the BH to be removed from the galaxy in the absence of dynamical friction. To prevent this, individual black holes are given an artificial dynamical mass of which is independent of their physical mass gained through accretion. However, if no MBH is present at the centre of the galaxy due to dynamical friction at a particular timestep, we consider the nearest black hole to the host galaxy centre as the central black hole throughout this paper.
II.3 TDE Rate Calculations From Simulated Galaxies
To compute the TDE rates of the simulated galaxies, we mainly focus on the properties and dynamics of the stars in the nuclear regions of the galaxies across the redshift range . We begin by selecting the galaxies with total stellar mass in the four simulation runs A1, A2, A4 and A8 which is the threshold halo mass above which seed BHs are formed in FIRE-2, as discussed in Section II.2.2.
We then obtain the stellar densities around these individual galaxies. Examples of stellar distribution and radial stellar density profiles at the galaxy centres are shown in Appendix A. Figure A.1 shows the individual projected stellar densities at the centres of the four primary host galaxies at for A1, A2, A4 and A8 runs, respectively. Using these stellar densities, we first obtain the radial density profiles of stars at the centre of each galaxy within kpc radius. The reason for choosing the galaxy centre instead of the central BH position lies in the occasional absence of BHs at the galaxy centre due to dynamical friction as discussed above. Furthermore, we fit each stellar density profile with a double power law function as follows:
| (11) |
where is the volume density of stars at scale radius . and are the power-law indices for the inner and outer density profiles respectively, and denotes the transition between the slopes of these two power-law profiles.
Based on these fitted parameters of the stellar density profiles, we further constrain our galaxy sample. We discard the unphysical systems with inner slopes steeper than their outer slopes and select only those with [23, 26]. A large spurious inner slope could be attributed to the numerical noise of the simulation caused by a handful of star particles present in the galaxy centers. An example of such stellar profile and the corresponding stellar distribution map are shown in Figure A.2. Moreover, we restrict our sample to those with due to numerical reasons. The lower limit on ensures a positive stellar distribution function (Equation 8), whereas the upper limit prevents the enclosed stellar mass from diverging. Furthermore, we discarded the BHs with beyond which no TDEs are produced as the stars are captured as a whole by the BHs. Figure 1 shows examples of stellar profiles at the centers of the primary host galaxies at z=1 in the A1, A2, A4, and A8 runs. It should be noted that in A2 and A8 runs. Hence, these two systems do not contribute to the TDE rates at in the subsequent analysis.
We estimate the rates of TDEs for this conservative galaxy sample using the publicly available code phaseflow 22
2
phaseflow included in the agama software library
https://github.com/GalacticDynamics-Oxford/Agama [61, 62]. The phaseflow code is used to solve for the distribution function from the density profile by Eddington inversion (Equation 8) using the fitted parameters of the double power law profile (, , , and ). phaseflow can also solve for the circular angular momentum , orbit-averaged angular momentum (), orbital period of stars and loss cone filling factor (Equation 9). With this information, we calculate the TDE rates per galaxy per unit time using Equation 7. Throughout this paper, we have assumed monochromatic stellar population, where all the stars are Sun-like with and . In the following sections we present the TDE rates across and their correlation with MBH and host galaxy stellar properties.
III Results
We study the TDE rates (TDR) in FIRE-2 galaxies and the corresponding redshift evolution in Section III.1, followed by their correlation with the star formation rate in Section III.2. Correlated properties of the TDR with their host galaxy and BH masses are studied in Section III.3. Further investigation of the impact of the host galaxy stellar properties on the TDR is done in Section III.4. Finally, we examine the TDR in the satellite galaxies and their detectability with the latest telescopes in Section III.5.
III.1 Redshift Evolution of TDE Rates
We begin by calculating the TDR for the individual galaxies at specific redshifts, and show these rates as a function of the black hole mass (Figure 2) and the galaxy mass (Figure 3). The galaxy sample, sparse at , grows over time. Small galaxies and IMBHs dominate at ; massive galaxies (and SMBHs) appears by , which gradually merge to become very massive galaxies by . Despite the large scatter in the TDRs at fixed redshift, the figures hint at an increasing trend with the and in the early universe.
To understand how the TDR evolves across the cosmic time, we plot in Figure 4 the averaged TDR at individual redshifts as a function of . TDR appears to peak around , consistent with the Figures 2 and 3. We fit the high-redshift TDR using a power-law function of , and obtain the following best-fit relation:
| (12) |
We further notice a moderately decreasing TDR at , although caution is needed when interpreting this trend, as the simulated galaxies were evolved only down to . However, we understand from the cosmic star formation history that the SFR declines in the low redshift [40]. Hence, we expect a decline in the TDR in the local universe if the TDR follows the SFR as noted in the high redshift.
Such evolution of the TDR with redshift could be attributed to the black hole model, star formation rate and their host galaxy stellar properties in FIRE-2, which we investigate in details in the following sections. The redshift evolution of the SFR is also shown in the same plot for a comparison, which will be further discussed in the next section.
III.2 Impact of Star Formation Rate on TDR
In this section, we investigate how the TDR and the star formation rates (SFR) are correlated in the host galaxies. For the simulated galaxies, we calculate the SFR within a central radius that contains of the stellar mass of host galaxy. As noted from the Figure 4, TDR closely follows SFR, where both increase as redshift decreases until , below which the SFR is appeared to get saturated. However, the comparison of the TDR with the SFR is based on the particular snapfiles at discrete redshifts. Hence, a direct comparison of these two parameters might not be straight forward. Nevertheless, we notice, at least at high redshift, the peak of TDR is following the peak of SFR. This might be indicative of some other mechanism, following a starburst phase, that is contributing to the enhanced TDR. This hints towards a similar phenomena in the local universe where an abundance of TDEs are found in the post-starburst galaxies [18, 20].
We further investigate in Figure 5 the correlation between the TDR and the SFR in individual galaxies, along with their fitted scaling relation. We notice that the averaged TDR is well correlated with the SFR across the entire redshift range of . TDR can be fitted as a simple power-law function of the SFR as follows:
| (13) |
Two galaxies in Figure 5 produce very high TDRs (). Both are located at in the A2 run, labeled ‘G1A2’ (dark green star, the primary host) and ‘G2A2’ (dark green triangle). Despite both having the highest TDRs, G2A2 has a much lower SFR (), while G1A2 shows the highest SFR among all the MassiveHalo simulations (). Therefore, it is clear that additional factors beyond the SFR must influence the TDR of host galaxies, which motivated us to investigate how the host galaxy stellar properties impacts the TDR. We will discuss this in Section III.4.
III.3 Dependence of TDE Rates on and
We further investigate how the TDR correlates with the host galaxy mass and the associated black hole mass. Since the number of galaxies at any specific redshift is small and their mass ranges are limited, we combine all the galaxies across all the redshifts for this analysis. This is shown in Figure 6, where the left and right panel represents the TDR as a function of and , respectively. It should be noted that the TDEs from the primary host galaxies are marked as stars at all the redshifts, whereas the smaller galaxies (more on this is discussed in Section III.5) are marked as points.
The left panel of Figure 6 shows that the TDR increases with the in the intermediate-mass regime, and then turns over and decreases at higher the . To quantify this trend, following [10], we equally bin the ranges on a logarithmic scale (shown as the gray shaded regions) and then calculate the averaged TDR in each mass bin. Finally, we fit the averaged TDR as a function of using a broken power-law function as follows:
| (14) |
where is the BH mass at which TDR turns around.
Intriguingly, similar trends – an increasing TDR in the IMBH regime and a turnover of the TDR around – have been found for the TDEs in the local universe in the previous studies [49, 10, 22]; however, with different normalisation, turnover BH masses and power-law slopes. Similarly, when focusing solely on SMBHs, although the overall decreasing trend of TDR agrees with the previous studies [63, 58, 48], the power-law slope and normalization show discrepancies. The differences in the fitted parameters between this work and previous studies arises from several factors: e.g., using the observed versus simulated galaxy stellar profiles, different BH mass measurement techniques, the presence or absence of dense nuclear star clusters, and – most importantly – the distinct redshift ranges probed. Most of the previous studies have used the BH samples that are limited to the local universe, while this work focuses on samples at . Hence, differences in the scaling relations are expected. However, qualitative agreement on the overall trend of the TDR with the and between this work (at high redshift) and the previous studies (in the local universe) hints at similar galaxy structures and BH-galaxy scaling across the cosmic time.
Likewise, the right panel of the Figure 6 shows that the trend between the averaged TDR and is similar to that seen for the . Specifically, the TDR rises with the for small galaxies, peaks at , and gradually declines at the higher masses. This trend can be fit to a function as follows:
| (15) |
III.4 Impact of Host Galaxy Stellar Properties on TDR
We begin with Figure 7 that shows the projected stellar density distributions of G1A2 and G2A2 galaxies, associated with the highest TDR (Figure 5). We find a dense stellar component in G1A2, which has a higher galaxy mass () and a central SMBH (). In contrast, G2A2 is a smaller galaxy () hosting an IMBH () with the lower central stellar density. This motivates us to further investigate the stellar profiles of our galaxy sample and their impact on the TDR.
Figure 8 shows the correlation of the TDR with the inner stellar density (left panel) and the inner slope (right panel) of the simulated galaxy sample used in this work, with G1A2 and G2A2 highlighted. As expected, the TDR generally increases with the higher due to a richer stellar population available for disruption. Steeper inner slopes also enhance the TDR by driving the stars into the loss cone [10]. The highest TDR in G1A2 (star symbol) results from a combination of the high and steep , whereas in G2A2 (triangle symbol), a very steep alone yields an elevated TDR despite its relatively lower .
Furthermore, a closer look at the Figure 8 reveals the distinct redshift trends in both panels. The inner stellar densities of galaxies grow with the redshift, with the most dense systems found at . The right panel further shows a bimodality in the inner power-law slope : it increases at the early times until , after which the trend reverses and the inner stellar profiles become shallower at the lower redshifts.
We further examine the black hole radius of influence (), which significantly affects the TDR (see Section II.1). We note from the left panel of Figure 9 that the IMBHs dominate in the high redshift. This could be a result of the bursty star formation activity in the FIRE-2 that depletes the available gas reservoir needed for the BHs growth at high [2]. the increase of the while keeping the almost constant (Figure 2 of [2]) explains the decreasing at in the left panel of Figure 9. Correspondingly, a rise in the TDR in this epoch due to the boosted and the rapid SFR is noted from the right panel. Star formation becomes steady below , and the BHs grows rapidly to become more massive, leading to increased , seen from Figure 9, left panel. However, much larger galaxies hosting the SMBH will lower the ratio of effective radius of the galaxy bulges, which disperse the stars, reducing the TDR in these systems [10], which is seen from the right panel (yellow and red colours). This claim is also supported by the shallower inner slope of the galaxies at low redshift in Figure 8.
Finally, in Figure 10, we show the connection of the TDR with , which is the stellar density at the influence radius. A strong correlation between the two parameters is noticed, consistent with the previous studies [48, 10]. This can be naturally explained as a consequence of the highest stellar flux at the critical radius, which is the same as for , as discussed in Section II.1. We also find a large influence density in both the G1A2 and G2A2, explaining the surge in the TDR in these galaxies.
III.5 TDEs in Satellite Galaxies
In this section, we focus on examining the TDEs in the satellite galaxies, defined as all the non-primary galaxies in our simulated galaxy sample. Because these satellite galaxies are offset from the primary (most massive) galaxy in the cluster, TDEs occurring within them may be observed as the off-nuclear TDEs, depending on the satellite’s size and the distance. Detecting such TDEs is crucial for identifying IMBHs in the globular clusters and dwarf galaxies, and for tracing the galaxy merger histories across the cosmic time.
As shown in Figure 6, satellite galaxy TDRs can be as high as those of the primary galaxies, spanning a broad range from to across all redshifts. In Figure 11, we present the fraction of the total TDR from the satellites versus the primary galaxies as a function of . It can be seen that the satellite TDR fraction increases substantially at high redshifts, highlighting their importance for studying the IMBHs and galaxy merger histories in the early cosmos.
Figure 12 shows the TDR of the individual satellite galaxies as a function of the distance from their primaries, with the transverse physical distances converted to the angular sizes in the local universe. The angular resolutions of Roman, Rubin, and ULTRASAT [54] are marked in the figure. Notably, almost all the satellite TDEs in our sample can be resolved by these missions in most cases. Nevertheless, the actual detectability depends on the additional factors, including the TDE luminosity, instrument flux limits, and survey specifications (e.g., field of view and limiting magnitude). A detailed analysis is beyond the scope of this work, and we defer a thorough investigation to a future study.
IV Summary and Discussions
In this work, we present a comprehensive analysis of how TDE rates evolve over cosmic history spanning a wide redshift range of , based on a sample of simulated galaxies drawn from the high-resolution cosmological zoom-in simulation FIRE-2. The simulation covers a diverse population of black holes and galaxies, ranging from the SMBHs in primary hosts to IMBHs in satellite galaxies, which allows us to study how TDE rates correlate with host galaxy properties and their redshift evolution. Below we summarize the key findings of this paper:
- 1.
- 2.
We observe a strong correlation of the TDR with the global SFR of their host galaxies (Figure 5) and estimate a scaling relation between TDR and SFR across the entire redshift range (Equation 13). We find that TDR closely follows the SFR, where both increase with time in the early universe before reaching a peak at , and moderately decline afterwards.
- 3.
We find that TDR correlates well with both and for the overall galaxy sample (Figure 6, Equation 14, Equation 15 ). Specifically, TDR increases with BH mass in the IMBH regime, peaks at (and ), and then declines at the high-mass end. This trend, previously seen in local galaxies, intriguingly persists at high redshifts.
- 4.
We also examine various components of the stellar distributions in galactic centers at different redshifts in connection to their associated TDE rates, shown in Figure 7, 8, 9 and 10. The combined findings from these analyses hint that TDR is large in galaxies with high inner density , steep inner slope , small influence radius and large stellar density at this scale .
- 5.
Finally, we examine the detectability of TDEs in satellite galaxies. The fraction of TDEs originating from satellite galaxies increases significantly at high redshifts, underscoring their potential as probes of IMBHs and galaxy assembly in the early universe. Encouragingly, most of these events fall within the angular resolution limits of next-generation facilities such as Roman, Rubin, and ULTRASAT, indicating that satellite TDEs are detectable with appropriate survey strategies.
Recent studies have explored the redshift evolution of the TDR, based on theoretical models [49, 42, 33]. We note that drawing firm quantitative conclusions on TDR evolution with redshift using current cosmological simulations faces several limitations. First, feedback physics such as AGN feedback, which is not included in the MassiveHalo suite, can affect host galaxy properties [15, 30, 31, 35, 6]. Its absence may artificially enhance SFRs and central stellar densities, thereby boosting TDE rates, particularly at [64, 44, 11]. On a more positive note, our sample is dominated by lower-mass BHs (), for which AGN feedback is negligible. Second, the finite resolution of FIRE-2 – although the best currently available – still prevents us from directly resolving nuclear star clusters, which are known to substantially boost local TDRs. Third, because the FIRE-2 simulation was evolved only to to mitigate over-cooling [65], our TDR constraints do not extend to . Finally, we note that being the zoom-in simulations, the halo mass range is narrow in the MassiveHalo suite. This limits the total number of galaxies available for TDE analysis. However, the combined numbers of all the galaxies across is significantly larger compared to previous studies.
Despite these limitations, we demonstrate that constraining the cosmic evolution of TDR using cosmological simulations is a promising avenue, and we encourage future work to revisit these TDR calculations with next-generation simulations. Current and next-generation telescopes, including Roman, Rubin, and ULTRASAT, are expected to detect numerous TDEs across a broad redshift range, providing an ideal testbed for comparing our predictions with observations. Such comparisons will allow us to validate and refine models of TDE rates, constrain the interplay between black hole growth, stellar dynamics, and galaxy evolution across cosmic time, and ultimately establish TDEs as a reliable probe of the co-evolution of black holes and their host galaxies from the early universe to the present day.
Appendix A Stellar Densities and Radial Profiles in FIRE-2
Figure A.1 shows the stellar densities in the most massive galaxies at in the A1 (top left), A2 (top right), A4 (bottom left) and A8 (bottom right) Massive Halo suite. Furthermore, Figure A.2 displays an example of the stellar distribution and the corresponding density profile with fitted parameters , type of systems that are discarded from our sample. It can be noticed from the left panel that a few stellar particles are present at the centre, causing a steep slope in the density profile, which could be a numerical artifact of the simulation. Finally, in Figure A.3 we show the radial stellar density profiles of all the galaxies with . We fit these profiles with double power-law functions (Equation 11) and choose only those that satisfy our selection criteria, mentioned in Section II.3 .
References
- [1] (2022) A very luminous jet from the disruption of a star by a massive black hole. Nature 612 (7940), pp. 430–434. External Links: Document, 2211.16530 Cited by: §I.
- [2] (2017) Black holes on FIRE: stellar feedback limits early feeding of galactic nuclei. MNRAS 472 (1), pp. L109–L114. External Links: Document, 1707.03832 Cited by: §I, §II.2.1, §III.4.
- [3] (2022) A fast-rising tidal disruption event from a candidate intermediate-mass black hole. Nature Astronomy 6, pp. 1452–1463. External Links: Document, 2209.00018 Cited by: §I.
- [4] (1987) Galactic dynamics. Cited by: 1st item.
- [5] (2015) Swift J1112.2-8238: a candidate relativistic tidal disruption flare. MNRAS 452 (4), pp. 4297–4306. External Links: Document, 1507.03582 Cited by: §I.
- [6] (2024) Effects of Multichannel Active Galactic Nuclei Feedback in FIRE Cosmological Simulations of Massive Galaxies. ApJ 973 (2), pp. 149. External Links: Document, 2310.16086 Cited by: §IV.
- [7] (2023) COSMOS-Web: An Overview of the JWST Cosmic Origins Survey. ApJ 954 (1), pp. 31. External Links: Document, 2211.07865 Cited by: §I.
- [8] (2012) Swift J2058.4+0516: Discovery of a Possible Second Relativistic Tidal Disruption Flare?. ApJ 753 (1), pp. 77. External Links: Document, 1107.5307 Cited by: §I.
- [9] (2026) Merger Driven or Internal Evolution? A New Morphological Study of Tidal Disruption Event Host Galaxies. arXiv e-prints, pp. arXiv:2602.06839. External Links: Document, 2602.06839 Cited by: §I.
- [10] (2025) Rates of Stellar Tidal Disruption Events around Intermediate-mass Black Holes. ApJ 980 (2), pp. L22. External Links: Document, 2407.09339 Cited by: §I, §II.1, §III.3, §III.3, §III.4, §III.4, §III.4.
- [11] (2023) The impact of AGN-driven winds on physical and observable galaxy sizes. MNRAS 523 (2), pp. 2409–2421. External Links: Document, 2303.12858 Cited by: §IV.
- [12] (2025) High-z Stellar Masses Can Be Recovered Robustly with JWST Photometry. ApJ 978 (2), pp. L42. External Links: Document, 2412.02622 Cited by: §I.
- [13] (2021) The Physics of Accretion Discs, Winds and Jets in Tidal Disruption Events. Space Sci. Rev. 217 (1), pp. 12. External Links: Document Cited by: §I.
- [14] (2008) Direct Cosmological Simulations of the Growth of Black Holes and Galaxies. ApJ 676 (1), pp. 33–53. External Links: Document, 0705.2269 Cited by: §II.2.2.
- [15] (2016) The HORIZON-AGN simulation: morphological diversity of galaxies promoted by AGN feedback. MNRAS 463 (4), pp. 3948–3964. External Links: Document, 1606.03086 Cited by: §IV.
- [16] (1989) The Tidal Disruption of a Star by a Massive Black Hole. ApJ 346, pp. L13. External Links: Document Cited by: §I.
- [17] (2017) Colours, star formation rates and environments of star-forming and quiescent galaxies at the cosmic noon. MNRAS 470 (1), pp. 1050–1072. External Links: Document, 1610.02411 Cited by: §II.2.1.
- [18] (2016) Tidal Disruption Events Prefer Unusual Host Galaxies. ApJ 818 (1), pp. L21. External Links: Document, 1601.04705 Cited by: §III.2.
- [19] (2020) The Host Galaxies of Tidal Disruption Events. Space Sci. Rev. 216 (3), pp. 32. External Links: Document, 2003.02863 Cited by: §I.
- [20] (2018) A Dependence of the Tidal Disruption Event Rate on Global Stellar Surface Mass Density and Stellar Velocity Dispersion. ApJ 853 (1), pp. 39. External Links: Document, 1707.02986 Cited by: §III.2.
- [21] (2013) Hydrodynamical Simulations to Determine the Feeding Rate of Black Holes by the Tidal Disruption of Stars: The Importance of the Impact Parameter and Stellar Structure. ApJ 767 (1), pp. 25. External Links: Document, 1206.2350 Cited by: §I.
- [22] (2024) Counting the Unseen II: Tidal Disruption Event Rates in Nearby Galaxies with REPTiDE. arXiv e-prints, pp. arXiv:2412.19935. External Links: Document, 2412.19935 Cited by: §III.3.
- [23] (2020) But what about…: cosmic rays, magnetic fields, conduction, and viscosity in galaxy formation. MNRAS 492 (3), pp. 3465–3498. External Links: Document, 1905.04321 Cited by: §II.3.
- [24] (2018) FIRE-2 simulations: physics versus numerics in galaxy formation. MNRAS 480 (1), pp. 800–863. External Links: Document, 1702.06148 Cited by: §II.2.1.
- [25] (2015) A new class of accurate, mesh-free hydrodynamic simulation methods. MNRAS 450 (1), pp. 53–110. External Links: Document, 1409.7395 Cited by: §II.2.1.
- [26] (2024) The proto-galaxy of Milky Way-mass haloes in the FIRE simulations. MNRAS 527 (4), pp. 9810–9825. External Links: Document, 2307.15741 Cited by: §II.3.
- [27] (2019) LSST: From Science Drivers to Reference Design and Anticipated Data Products. ApJ 873 (2), pp. 111. External Links: Document, 0805.2366 Cited by: §I.
- [28] (2025) An Intermediate-mass Black Hole Lurking in A Galactic Halo Caught Alive during Outburst. arXiv e-prints, pp. arXiv:2501.09580. External Links: Document, 2501.09580 Cited by: §I.
- [29] (2024) Detecting Population III Stars through Tidal Disruption Events in the Era of JWST and Roman. ApJ 966 (2), pp. L33. External Links: Document, 2401.12752 Cited by: §I.
- [30] (2020) Cosmological Simulation of Galaxy Groups and Clusters. I. Global Effect of Feedback from Active Galactic Nuclei. ApJ 889 (1), pp. 60. External Links: Document, 1911.07824 Cited by: §IV.
- [31] (2022) Cosmological Simulation of Galaxy Groups and Clusters. II. Studying Different Modes of Feedback through X-Ray Observations. ApJ 940 (1), pp. 47. External Links: Document, 2209.13349 Cited by: §IV.
- [32] (2025) JWST Discovery of a High-Redshift Tidal Disruption Event Candidate in COSMOS-Web. arXiv e-prints, pp. arXiv:2504.13248. External Links: Document, 2504.13248 Cited by: §I.
- [33] (2026) Tidal disruption event rates across cosmic time: forecasts for LSST, Roman, and JWST and their constraints on the supermassive black hole mass function. arXiv e-prints, pp. arXiv:2602.04947. External Links: Document, 2602.04947 Cited by: §IV.
- [34] (2012) Mass Functions of Supermassive Black Holes across Cosmic Time. Advances in Astronomy 2012, pp. 970858. External Links: Document, 1112.1430 Cited by: §I.
- [35] (2024) On the origin of star formation quenching in massive galaxies at z 3 in the cosmological simulations IllustrisTNG. MNRAS 534 (4), pp. 3974–3988. External Links: Document, 2310.03083 Cited by: §IV.
- [36] (2023) A population of red candidate massive galaxies 600 Myr after the Big Bang. Nature 616 (7956), pp. 266–269. External Links: Document, 2207.12446 Cited by: §I.
- [37] (2023) Growth of a Massive Black Hole in a Dense Star Cluster Via Tidal Disruption Accretion. ApJ 943 (2), pp. 77. External Links: Document, 2211.02376 Cited by: §I.
- [38] (2016) The superluminous transient ASASSN-15lh as a tidal disruption event from a Kerr black hole. Nature Astronomy 1, pp. 0002. External Links: Document, 1609.02927 Cited by: §I.
- [39] (2018) A luminous X-ray outburst from an intermediate-mass black hole in an off-centre star cluster. Nature Astronomy 2, pp. 656–661. External Links: Document, 1806.05692 Cited by: §I.
- [40] (2014) Cosmic Star-Formation History. ARA&A 52, pp. 415–486. External Links: Document, 1403.0007 Cited by: §III.1.
- [41] (1999) Rates of tidal disruption of stars by massive central black holes. MNRAS 309 (2), pp. 447–460. External Links: Document, astro-ph/9902032 Cited by: §I.
- [42] (2025) Tidal Disruption Event Demographics in Supermassive Black Hole Binaries Over Cosmic Times. arXiv e-prints, pp. arXiv:2507.08082. External Links: Document, 2507.08082 Cited by: §IV.
- [43] (2019) Weighing Black Holes Using Tidal Disruption Events. ApJ 872 (2), pp. 151. External Links: Document, 1801.08221 Cited by: §I.
- [44] (2021) Realistic mock observations of the sizes and stellar mass surface densities of massive galaxies in FIRE-2 zoom-in simulations. MNRAS 501 (2), pp. 1591–1602. External Links: Document, 2009.10161 Cited by: §IV.
- [45] (2026) JWST and Keck observations of the off-nuclear tidal disruption event TDE 2025abcr: An evolving reprocessing layer. arXiv e-prints, pp. arXiv:2604.16093. External Links: 2604.16093 Cited by: §I.
- [46] (2019) Tidal disruption event rates in galaxy merger remnants. MNRAS 488 (1), pp. L29–L34. External Links: Document, 1903.09124 Cited by: §I.
- [47] (2021) Tidal disruption events in the first billion years of a galaxy. MNRAS 500 (3), pp. 3944–3956. External Links: Document, 2006.06565 Cited by: §I.
- [48] (2020) Enhancement of the tidal disruption event rate in galaxies with a nuclear star cluster: from dwarfs to ellipticals. MNRAS 497 (2), pp. 2276–2285. External Links: Document, 2003.08133 Cited by: §I, §II.1, §III.3, §III.4.
- [49] (2024) Demographics of tidal disruption events with L-Galaxies: I. Volumetric TDE rates and the abundance of nuclear star clusters. A&A 689, pp. A204. External Links: Document, 2312.13242 Cited by: §III.3, §IV.
- [50] (1988) Tidal disruption of stars by black holes of 10-10 solar masses in nearby galaxies. Nature 333 (6173), pp. 523–528. External Links: Document Cited by: §I.
- [51] (2021) Unveiling the Population of Wandering Black Holes via Electromagnetic Signatures. ApJ 916 (2), pp. L18. External Links: Document, 2107.02132 Cited by: §I.
- [52] (2014) Science with a Wide-field UV Transient Explorer. AJ 147 (4), pp. 79. External Links: Document, 1303.6194 Cited by: §I.
- [53] (2021) First tidal disruption events discovered by SRG/eROSITA: X-ray/optical properties and X-ray luminosity function at z ¡ 0.6. MNRAS 508 (3), pp. 3820–3847. External Links: Document, 2108.02449 Cited by: §I.
- [54] (2024) ULTRASAT: A Wide-field Time-domain UV Space Telescope. ApJ 964 (1), pp. 74. External Links: Document, 2304.14482 Cited by: §III.5.
- [55] (2015) Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report. arXiv e-prints, pp. arXiv:1503.03757. External Links: Document, 1503.03757 Cited by: §I.
- [56] (2005) Modelling feedback from stars and black holes in galaxy mergers. MNRAS 361 (3), pp. 776–794. External Links: Document, astro-ph/0411108 Cited by: §II.2.2.
- [57] (2020) Rates of Stellar Tidal Disruption. Space Sci. Rev. 216 (3), pp. 35. External Links: Document, 2003.08953 Cited by: §II.1.
- [58] (2016) Rates of stellar tidal disruption as probes of the supermassive black hole mass function. MNRAS 455 (1), pp. 859–883. External Links: Document, 1410.7772 Cited by: §I, §II.1, §III.3.
- [59] (2011) Snacktime for Hungry Black Holes: Theoretical Studies of the Tidal Disruption of Stars. Ph.D. Thesis, University of California, Berkeley. Cited by: §II.1.
- [60] (2025) Broad-line AGNs at 3.5 ¡ z ¡ 6: The Black Hole Mass Function and a Connection with Little Red Dots. ApJ 986 (2), pp. 165. External Links: Document, 2409.06772 Cited by: §I.
- [61] (2017) A New Fokker-Planck Approach for the Relaxation-driven Evolution of Galactic Nuclei. ApJ 848 (1), pp. 10. External Links: Document, 1709.04467 Cited by: §II.3.
- [62] (2019) AGAMA: action-based galaxy modelling architecture. MNRAS 482 (2), pp. 1525–1544. External Links: Document, 1802.08239 Cited by: §II.3.
- [63] (2004) Revised Rates of Stellar Disruption in Galactic Nuclei. ApJ 600 (1), pp. 149–161. External Links: Document, astro-ph/0305493 Cited by: §III.3.
- [64] (2020) Measuring dynamical masses from gas kinematics in simulated high-redshift galaxies. MNRAS 497 (4), pp. 4051–4065. External Links: Document, 1908.05274 Cited by: §IV.
- [65] (2023) Public Data Release of the FIRE-2 Cosmological Zoom-in Simulations of Galaxy Formation. ApJS 265 (2), pp. 44. External Links: Document, 2202.06969 Cited by: §I, §II.2, §IV.
- [66] (2025) Second public data release of the FIRE-2 cosmological zoom-in simulations of galaxy formation. arXiv e-prints, pp. arXiv:2508.06608. External Links: Document, 2508.06608 Cited by: footnote 1.
- [67] (2025) A Massive Black Hole 0.8 kpc from the Host Nucleus Revealed by the Offset Tidal Disruption Event AT2024tvd. arXiv e-prints, pp. arXiv:2502.17661. External Links: Document, 2502.17661 Cited by: §I.
- [68] (2025) Science objectives of the Einstein Probe mission. Science China Physics, Mechanics, and Astronomy 68 (3), pp. 239501. External Links: Document, 2501.07362 Cited by: §I.