https://orcid.org/0000-0001-9458-821X
ABSTRACT
A tidal disruption event (TDE) occurs when a star is destroyed by the strong tidal shear of a massive black hole (MBH). The accumulation of TDE observations over the last years has revealed that post-starburst galaxies are significantly overrepresented in the sample of TDE hosts. Here we address the post-starburst preference by investigating the decline of TDE rates in a Milky-Way like nuclear stellar cluster featuring either a monochromatic (1 |$\, \mathrm{M}{\odot {}}$|) or a complete, evolved stellar mass function. In the former case, the decline of TDE rates with time is very mild, and generally up to a factor of a few in 10 Gyr. Conversely, if a complete mass function is considered, a strong TDE burst over the first 0.1–1 Gyr is followed by a considerable rate drop, by at least an order of magnitude over 10 Gyr. The decline starts after a mass segregation time-scale, and it is more pronounced assuming a more top-heavy initial mass function and/or an initially denser nucleus. Our results thus suggest that the post-starburst preference can be accounted for in realistic systems featuring a complete stellar mass function, even in moderately dense galactic nuclei. Overall, our findings support the idea that starbursting galactic nuclei are characterized by a top-heavy initial mass function; we speculate that accounting for this can reconcile the discrepancy between observed and theoretically predicted TDE rates even in quiescent galaxies.
1 INTRODUCTION
The disruption of a star operated by a massive black hole (MBH) is known as tidal disruption event (TDE). These phenomena result in a multiband electromagnetic emission that can exceed the MBH Eddington luminosity by orders of magnitude; in addition, TDEs feature a distinctive decline in their bolometric luminosity, that is expected to drop roughly as t−5/3 over time-scales of months to years (e.g. Rees 1988; Lodato, King & Pringle 2009; Lodato & Rossi 2011). This makes TDEs ideal sources to unveil an otherwise quiescent population of MBHs, and their detection can provide crucial insights on the low-mass end of the MBH mass function, since the TDE rates seem to get larger if the central MBH is less massive (Wang & Merritt 2004; Stone & Metzger 2016).
The sample of observed TDEs has quickly grown to nearly a hundred in the last decades, thanks to the advent of wide-field transient surveys in the optical (as the Palomar Transient Factory, Rau et al. 2009; Pan-STARRS, Chambers et al. 2016; ASASSN survey, Shappee et al. 2014; the Zwicky Transient Facility, Bellm et al. 2019) X-ray (as ROSAT, XMM, Jansen et al. 2001), and gamma-ray band (e.g. Swift, Chandra observations). The number of detections is expected to get larger by orders of magnitude in the next years, thanks to the new observations from the presently operating eROSITA (Merloni et al. 2012) and the coming LSST/Vera Rubin (Ivezić et al. 2019) and Athena (Barcons et al. 2012) observatories.
The opening of this golden era for TDEs has come with its own surprises. Together with the accumulation of the first dozens of detections, population studies of TDE hosts have revealed that these events are significantly overrepresented in post-starburst galaxies (Arcavi et al. 2014; French, Arcavi & Zabludoff 2016, 2017; Law-Smith et al. 2017; Graur et al. 2018; French et al. 2020; Hammerstein et al. 2021). The spectra of these galaxies lack strong emission lines, suggesting low current star formation rates, but feature strong Balmer absorption lines, a signature that reveals the presence of many young stars, and it is associated with a starburst over the last ≈1 Gyr (see e.g. French 2021 for a recent review). Post-starburst galaxies are typically subdivided into the more bursty, rarer E+A galaxies (that make up |$\approx 0.2\rm \, per\, cent$| of low-redshift galaxies, and whose recent starburst created |$\gt 3\rm \, per\, cent$| of their current stellar mass over 25–200 Myr) and quiescent Balmer-strong galaxies (making up |$\approx 2\rm \, per\, cent$| of local galaxies, and which formed |$\gt 0.1\rm \, per\, cent$| of their current stellar mass over 25 Myr−1 Gyr, French et al. 2017, 2018). The overrepresentation of TDEs in post-starbursts was initially estimated to be as large as a factor ∼100 for E+A and ∼30 for quiescent Balmer strong hosts (French et al. 2016; Graur et al. 2018). Law-Smith et al. (2017) showed that, once accounting for selection effects, the E+A post-starburst preference should be reduced by a factor ≈4, bringing the overall post-starburst preference to nearly 15–30, as later confirmed by Hammerstein et al. (2021).
In order to investigate the origin of the post-starburst preference, it is important to properly model the dynamical processes that trigger the generation of TDEs. The most robustly understood and ubiquitous TDE production mechanism is stellar two-body relaxation, according to which the granularity of the stellar system, and more specifically the two body interactions between pairs of stars, slowly deflect stellar orbits, some of which may eventually reach the MBH vicinity and get destroyed. Stone et al. (2018) explored a set of possible dynamical drivers for the post-starburst overrepresentation, and find that the most compelling explanation is that the host may be initially exceptionally cuspy near its centre (with density profile as steep as ρ∝−γ, γ ≳ 2.5, see red lines in Fig. 6), so that its initially extremely short relaxation time – or equivalently, extremely efficient two body relaxation – can guarantee a substantial drop of TDEs in time. However, this picture requires stellar densities to reach impressively large values about the MBH. Even if post-starburst galaxies appear to harbour steep density profiles (see e.g. Stone & van Velzen 2016; French et al. 2020 and references therein), currently available observational campaigns can hardly resolve the MBH influence radius, and more resolved observations are required to probe the presence of the extreme cusps suggested by Stone et al. (2018). Furthermore, it is unclear whether star formation can proceed very efficiently so close to the MBH; molecular clouds should in fact get disrupted by its extreme gravitational field prior to forming stars (Sanders 1998). Alternative theories to explain the post-starburst overrepresentation are the fact that a large number of stars can be initially set on radial orbits at the starburst stage (Stone et al. 2018); otherwise, the formation of eccentric nuclear discs about the MBH – resulting from the gas infall accompanying the starburst – can initially ferry a large amount of stars in the vicinity of the MBH (Madigan et al. 2018; this theory requires a large fraction of post-starburs galaxies to host eccentric discs; those could develop in galaxy mergers that then trigger the starburst, as suggested by Hopkins & Quataert 2010).
Most theoretical studies estimating dynamically induced TDE rates only account for idealized, monochromatic stellar systems composed by |$1\, \mathrm{M}{\odot {}}$| stars. However, real stellar systems are composed by an extended stellar population that evolves with time, and accounting for this aspect could sensibly impact the dynamical modelling of TDE rates. In fact, more massive (but generally rarer) stellar objects are expected to segregate towards the centre of the galactic nucleus, leaving lighter objects farther from the MBH, a process known as mass segregation (Bahcall & Wolf 1977; Alexander & Hopman 2009; Preto & Amaro-Seoane 2010; Alexander 2017). Mass segregation generally occurs over a fraction of the system two-body relaxation time (Merritt 2013), and its characteristic time-scale can be much shorter than a Hubble time in galaxies harbouring nuclear star clusters. This process may be very relevant for the nuclei of post-starburst galaxies, as those appear to be concentrated, and their typical galaxy and MBH masses suggest they likely host nuclear stellar clusters (Neumayer, Seth & Böker 2020; Hoyer et al. 2021), implying they may have relatively short mass-segregation time-scales. Furthermore, recent studies (e.g. Zhang et al. 2018; Toyouchi et al. 2021) found that the stellar initial mass function (IMF) characterizing starbursting systems is generally top-heavy (note that the Mikly Way nucleus, which is by far the best known galactic nucleus, is also characterized by a top-heavy IMF, see e.g. Bartko et al. 2010; Lu et al. 2013), meaning that more massive stars are formed more efficiently in comparison with the commonly adopted IMFs (e.g. Kroupa 2001). A top-heavy IMF can enhance the effects of mass segregation and relaxation, that ultimately set the TDE rates, and accounting for this aspect is crucial to properly track the time evolution of TDE rates.
Motivated by this, in this study we evaluate the importance of accounting for a complete, evolved initial stellar mass function in setting the time evolution of TDE rates; we also explore how varying the IMF and galaxy density profile affects the evolution of TDE rates. Those interested in a cursory reading can find a summary of the methodology at the beginning of Section 2 and a summary of the results in Section 4. The paper is structured as follows: in Section 2 we detail the methodology adopted in this study; in Section 3 detail our results; in Section 4 we summarize and discuss our findings and present the final remarks.
2 METHODS
In this section, we detail our implementation for addressing TDE rates. Below we briefly summarize the key aspects of our methodology, which are described in detail in the following subsections.
We adopt a Milky-Way like MBH and galaxy profile. We perform a series of runs accounting for different components: (i) the stellar bulge only; (ii) the stellar bulge and the observed nuclear star cluster (forming instantaneously at t = 0), exploring different values for the compactness of the latter (iii) the stellar bulge and the nuclear cluster, 30 per cent of which is formed during the simulation over a time-scale of 200 Myr.
We account for a complete mass function for the system, assuming either a Kroupa IMF, or a top-heavy IMF with different power-law slopes (∝m−α, α = {1.5, 1.7, 1.9}) for masses above 0.5 |$\, \mathrm{M}{\odot {}}$|.
The IMF is converted into a mass function at a given time, and large mass main sequence stars are mapped into white dwarfs (WDs), neutron stars (NSs), stellar black holes (BHs) depending on their initial mass and on their age. Stars are associated a stellar radius that sets their value for the tidal disruption radius rt.
Only main sequence stars can undergo TDEs in our implementation (compact objects can directly enter the MBH horizon and increase its mass, but we do not account for these events in the reported rates).
Even if we account for an extended mass function, this cannot be evolved at run time in the simulation. For this, we implement a series of post-processing strategies that correct for the fact that stars age with time; note that these corrections change the event rates typically by a few per cent, and up to 25 per cent in the most extreme configurations.
The system is evolved in time for 14 Gyr using the one-dimensional Fokker-Plank code Phaseflow (Vasiliev 2017).
2.1 Definition of tidal disruption event
2.2 Background galaxy density profile and star formation treatment
Our model for the background galaxy is a Milky Way model. We implement it accounting for the stellar bulge and the central nuclear star cluster, both modelled as Sersic (1968) profiles: following Pfister et al. (2020), we adopt a bulge stellar mass of |$9.1\times 10^9 \, \mathrm{M}{\odot {}}$| (Licquia & Newman 2015), effective radius 1.04 kpc, and Sersic index 1.3 (Davis, Graham & Cameron 2019). If the nuclear star cluster is present, we model it following Pfister et al. (2020), whose analysis is based on the observational data by Schödel et al. (2018): the cluster effective radius is set to 6 pc, its Sersic index to 2, and its total mass is |$4\times 10^7 \, \mathrm{M}{\odot {}}$|. In all configurations, we add a central MBH of |$4\times 10^6 \, \mathrm{M}{\odot {}}$| (Gravity Collaboration 2020) that increases its mass by accreting stellar objects. We stress that the Milky Way properties and the MBH mass adopted for this study are similar to those inferred for many TDE hosts (e.g. French et al. 2020). The systems featuring a central nuclear star cluster expand their size over time (e.g. Merritt 2013; Vasiliev 2017); it is thus likely that the nuclear star cluster we observe today was initially more concentrated. For this reason, we also explore the evolution of analogous galaxies with the same properties as above but with the nuclear star cluster Sersic index set to 3 and to 5.
In most runs we assume the nuclear star cluster (that sets TDE rates) to be formed instantaneously at t = 0. This is a reasonable assumption as the starburst duration is much shorter than the age of the system (the same is assumed in Stone et al. 2018); in addition, the fraction of starbursting mass in post-starburst galaxies is compatible with the possibility of the entire nuclear cluster to be formed in the starburst. For comparison, we also run two cases for which the nuclear star cluster is present and relaxed, but it undergoes extended star formation after 5 Gyr; the stellar cluster grows its mass by 30 per cent over 200 Myr (|$6\times 10^{-2}\, \mathrm{M}{\odot {}}$| yr−1), after which star formation is quenched and the system is further evolved in time. Star formation is implemented by adding a source term in the Fokker-Plank treatment, it is constant in time and occurs within r ≈ 6 pc, with the form |$\dot{\rho }\propto r^{-1/2}$| (see Vasiliev 2017 for more details).
2.3 The IMFs
2.4 Stellar evolution and stellar parameters
We now describe a series of prescriptions adopted to map the initial mass of stars into their mass and radius at any given time. Note that the age of the stellar population needs to be set at the beginning of the integration, but it remains constant over the course of the simulation, meaning that stellar masses and radii are not evolved at run time. However, we adopt a series of post-processing routines to correct for this simplification, as described in Section 2.4.4.
2.4.1 Main sequence stars
2.4.2 Compact stellar remnants
Stars whose mass is initially above |$8.7\, \mathrm{M}{\odot {}}$| are mapped into the mass of the associated compact remnant (NS, BH). We map their initial mass into their final mass adopting the PARSEC evolutionary tracks at Solar metallicity (Bressan et al. 2012; Tang et al. 2014); the mass of the compact remnant is obtained according to Spera, Mapelli & Bressan (2015), assuming the delayed supernova model; the initial to final mass of such massive stars can be seen in fig. 1 of Bortolas, Mapelli & Spera (2017).
The panels show the time evolution of TDE rates for different choices of the background system and IMF. The top left-hand, top right-hand, bottom left-hand, bottom right-hand panels respectively show the evolution of the TDE rates associated with a top heavy IMF with α = 1.5, 1.7, 1.9 and a Kroupa IMF. In all panels, the thin grey lines show the TDE rates associated with a monochromatic IMF (1|$\, \mathrm{M}{\odot {}}$| stars). The solid, dashed, dotted lines in each plot refer to nuclear star clusters whose Sersic index is equal respectively to 2, 3, 5. The TDE rates significantly decline at late times for all the chosen IMFs, except for the monochromatic one. The ‘x’ on the top of each coloured line marks the associated value of the mass segregation time-scale, computed via equation (10). The event rates start dropping significantly in all runs after this time-scale.
WDs, NSs, and BHs are lost into the MBH (and all their mass is added to the MBH) if they get closer than r• = 4GM•/c2 (here G is the gravitational constant, and c is the speed of light in vacuum). These events may give rise to detectable gravitational wave signals (the so-called EMRIs, see e.g. Zwick et al. 2020, 2021; Vázquez-Aceves et al. 2022) but they are not TDEs.1
2.4.3 Mass families in the Fokker–Plank implementation
The Fokker Plank code adopted for the evolution of the system does not allow us to include a continuous mass function. For this, we divide the initial stellar mass range into 100 logarithmic bins,2 each characterized by its initial and evolved stellar mass, number of objects, stellar radius and so forth. Each mass bin is called a mass family.
2.4.4 Post-processing corrections for stellar evolution
Our choice for mbreak is |$2 \, \mathrm{M}{\odot {}}$| (i.e. the population is ≈1.8 Gyr old; see the beginning of Section 2.4); all mass bins characterized by an initial stellar mass exceeding this value are treated as compact objects during the simulation. We checked that our specific choice of mbreak does not significantly impact the results of this work. We implement two post-processing corrections that account for the fact that the stellar population starts out younger and eventually gets older than the time associated with mbreak.
First of all, for t > t⋆(mbreak), we neglect all TDEs that are associated with a mass family that already became a compact stellar remnant at that given t. Secondly, for t < t⋆(mbreak), we estimate the TDE rate that would be associated with each k-th mass family with initial mass larger than mbreak as |$\dot{N}_k = \dot{N}_{\rm ref} (n_k/n_{\rm ref})$|; this allows us to correct for the fact that in the code we cannot destroy very massive stars that quickly become compact remnants. Here |$\dot{N}_k, n_k$| respectively represent the TDE rate and number of objects in the system associated with the k-th mass family, while |$\dot{N}_{\rm ref}, n_{\rm ref}$| represent the same quantities associated with a reference mass family that can be used as a proxy for the initially heavier k-th family. Here we choose the reference mass family to be the largest stellar mass still in the main sequence (2|$\, \mathrm{M}{\odot {}}$|). Both the implemented corrections are small, and together they never change the event rates by more than ≈25 per cent (the maximum correction is associated with the α = 1.5 IMF).
2.4.5 Correction for the total mass of the system
Even neglecting TDEs and focusing only on the stellar population, it is important to caution that its total mass declines in time, as main sequence stars are turned into compact objects of lower mass; for a given a system, the conserved quantity is the total number of objects. In order to consistently compute the TDE rates within a given system and normalize its total mass consistently for all the chosen IMFs, here we opt for defining the total mass of the system (Section 2.2) when the stellar population has a specific age. In what follows, we always define the total mass at t⋆=10 Gyr, or equivalently at |$m_{\rm break, ref} =1\, \mathrm{M}{\odot {}}$|, meaning that the total mass is set for a 10 Gyr old population, and the real total mass evolved in the simulation is (13, 7, 10, 13) per cent larger for a (Kroupa, α = 1.5, 1.7, 1.9) IMF. We do not apply any correction when we treat a monochromatic mass function.
2.4.6 Further approximations in the treatment of the mass function
In this paper, we do not account for the presence of binaries and multiplets. The larger effective mass of these systems compared to normal stars make them subject to mass segregation, and neglecting them may result in slightly underestimating the effectiveness of mass segregation. However, we expect this effect to be relatively small, and we believe the inclusion of multiplets can potentially enhance even further the early TDE rates.
Concerning giants, albeit their very large radii could in principle make them ideal candidates for TDEs, their low density stellar envelope makes them prone to less violent and thus more elusive events (e.g. Rossi et al. 2020), as the so-called ‘spoon-feeding’ (MacLeod et al. 2013). In addition, the giant stages are typically short compared to the other stages considered here. For this, we do not account for giants in the present implementation, and we assume stars whose mass is above mbreak to directly turn into compact remnants.
2.5 Fokker-Plank integration
The evolution of the system is performed by means of the Phaseflow Fokker-Plank integrator (Vasiliev 2017), which is part of the AGAMA toolkit (Vasiliev 2019). The code allows us to evolve in time an isotropic, spherically symmetric distribution of stars by solving the coupled system of Poisson and orbit-averaged Fokker-Planck 1D equations for the gravitational potential, the density, and the distribution function. The evolution of the system is driven by two-body relaxation and mass segregation. The system can account for the presence of a sink term representing the loss of objects into the MBH (whose mass is allowed to grow, and rt, r• are consistently adjusted accordingly) and possibly a source term that mimics star formation. The code also has the capability of evolving multiple stellar families with their own properties (described above). The code has been successfully tested in the context of e.g. nuclear star clusters evolution (Generozov et al. 2018; Emami & Loeb 2020), TDEs in massive high redshift clumps (Pestoni et al. 2021) and primordial BHs (Zhu et al. 2018; Stegmann et al. 2020).
In this work, we set the simulation parameters as follows: the phase-volume grid is chosen to have 400 logarithmically spaced points; we select an integration accuracy of 10−3 and the quadratic discretization method. The Coulomb logarithm for the evolution of the system is set to ln Λ = ln (0.4M•/〈m〉) (Spitzer & Hart 1971), with 〈m〉 defined in equation (6).
2.6 Static estimates of TDE enhancements
The TDE rate variation associated with adopting a complete, instead of a mono- or bi- chromatic stellar mass function, has been first explored by Magorrian & Tremaine (1999), and subsequently addressed by Kochanek (2016), Stone & Metzger (2016), Pfister et al. (2022);3 in particular, Stone & Metzger (2016) presented a simple back of the envelope estimate that we repropose and generalize here, based on the stellar evolution framework described above. The rate change in TDEs when considering a complete mass function, rather than a monochromatic mass function of 1 |$\, \mathrm{M}{\odot {}}$| stars, can be estimated by accounting for a series of effects.
In order to properly normalize the total mass of the system at any given time, we estimate the ratio in equation (8) accounting for the correction described in Section 2.4.5. This practically means that 〈m〉 appearing at the denominator (which in fact tunes the number of objects in the system) should be evaluated considering a population whose mf(m) is evaluated for a population of t⋆ = 10 Gyr, while all the other quantities are computed at the chosen time of the stellar population. The rate enhancements as computed from equation (8) are equal to respectively {7.3, 8.7, 7.9, 4.3} for IMFs with α = 1.5, 1.7, 1.9 and the Kroupa IMF; those values are evaluated for t⋆ =1.8 Gyr, corresponding to |$m_{\rm break}=2\, \mathrm{M}{\odot {}}$|.4It is important to note that the estimate neglects the time evolution of the system, the fact that stars have a different rt depending on their mass and the fact that mass segregation affects the evolution of the system. The same values computed from the PhaseFlow integrator (whose rates are presented in Section 3) at t = 0 are in the range |${\dot{N}_{\rm MF}}/{\dot{N}_{\rm mc}}=3-6$| if the nuclear cluster Sersic index is set to 2, and in the range |${\dot{N}_{\rm MF}}/{\dot{N}_{\rm mc}}=4-7$| if it is set to 5. The presented rough estimate is thus not too bad (even if necessarily approximate) at the beginning of the evolution, but it becomes completely unreliable as time passes, as shown below.
3 RESULTS
3.1 TDE rates from the stellar bulge
First of all, we address the TDE rates for a galaxy featuring the stellar bulge only, without a nuclear stellar cluster. We evolve the system, modelled with a Kroupa IMF, for 14 Gyr. The density profile of the system at the end of the integration is virtually indistinguishable from the initial conditions, suggesting relaxation is too weak to impact the evolution of the system. Accordingly, the TDE rate remains very low, nearly ≈2.3 × 10−7 yr−1 at all times. This implies the bulge contribution to the TDE rates is negligible compared to that of the nuclear star cluster, reported in the next section (in agreement with Pfister et al. 2020). This also means that TDEs can be assumed to come from stars in the nuclear star cluster only to a very good approximation, while the bulge virtually acts only as a background potential.
3.2 TDE rates accounting for a nuclear star cluster
Fig. 1 shows the time evolution of the TDE rates for a set of models that account for both the bulge and the nuclear star cluster. The results are shown for different IMFs and for nuclear cluster Sersic index equal to 2, 3, 5; in all panels, we also show the TDE rates associated with a monochromatic mass function. Remarkably, in all cases, the TDE rates significantly drop over time if one adopts an extended mass function, while the drop is very mild in the monochromatic case. Table 1 reports the peak to valley ratio for the event rate in each run; it also shows the ratio between the average event rate at t = 250–750 Myr and the same value at t = 10–12 Gyr. Except for the monochromatic case, the early event rates are significantly larger than the late time ones (∼×10–100, except for the less concentrated Kroupa case). For this specific choice of galaxy and MBH mass, the early to late TDE rates are thus compatible with the post-starburst preference (which is by nearly a factor 15–30) especially if considering top-heavy IMFs and/or initially more concentrated clusters.5
The table displays the ratio between the initially large event rates and the lower, late event rates for the runs shown in Fig. 1; each run is defined by its IMF (columns) and nuclear cluster Sersic index (rows). In particular, numbers out of the parenthesis show the ratio between the maximum and minimum value of TDE rates for each given run. The numbers in parenthesis, instead, show the ratio between the early TDE rate computed between 250 and 750 Myr, and the late event rate between 10 and 12 Gyr. The inclusion of a complete, evolved IMF induces a significant decline of the TDE rate with time; the same is not observed adopting a monochromatic mass function.
| Sersic | Top-heavy | ||||
|---|---|---|---|---|---|
| index | α = 1.5 | α = 1.7 | α = 1.9 | Kroupa | Monochr. |
| 2 | 29 (14.5) | 19 (11.3) | 11 (8.1) | 5 (3.7) | 2 (1.3) |
| 3 | 77 (19.7) | 47 (16.3) | 26 (12.8) | 9 (6.9) | 3 (2.2) |
| 5 | 312 (24.9) | 178 (21) | 91 (16.6) | 27 (9.4) | 7 (4.5) |
| Sersic | Top-heavy | ||||
|---|---|---|---|---|---|
| index | α = 1.5 | α = 1.7 | α = 1.9 | Kroupa | Monochr. |
| 2 | 29 (14.5) | 19 (11.3) | 11 (8.1) | 5 (3.7) | 2 (1.3) |
| 3 | 77 (19.7) | 47 (16.3) | 26 (12.8) | 9 (6.9) | 3 (2.2) |
| 5 | 312 (24.9) | 178 (21) | 91 (16.6) | 27 (9.4) | 7 (4.5) |
The table displays the ratio between the initially large event rates and the lower, late event rates for the runs shown in Fig. 1; each run is defined by its IMF (columns) and nuclear cluster Sersic index (rows). In particular, numbers out of the parenthesis show the ratio between the maximum and minimum value of TDE rates for each given run. The numbers in parenthesis, instead, show the ratio between the early TDE rate computed between 250 and 750 Myr, and the late event rate between 10 and 12 Gyr. The inclusion of a complete, evolved IMF induces a significant decline of the TDE rate with time; the same is not observed adopting a monochromatic mass function.
| Sersic | Top-heavy | ||||
|---|---|---|---|---|---|
| index | α = 1.5 | α = 1.7 | α = 1.9 | Kroupa | Monochr. |
| 2 | 29 (14.5) | 19 (11.3) | 11 (8.1) | 5 (3.7) | 2 (1.3) |
| 3 | 77 (19.7) | 47 (16.3) | 26 (12.8) | 9 (6.9) | 3 (2.2) |
| 5 | 312 (24.9) | 178 (21) | 91 (16.6) | 27 (9.4) | 7 (4.5) |
| Sersic | Top-heavy | ||||
|---|---|---|---|---|---|
| index | α = 1.5 | α = 1.7 | α = 1.9 | Kroupa | Monochr. |
| 2 | 29 (14.5) | 19 (11.3) | 11 (8.1) | 5 (3.7) | 2 (1.3) |
| 3 | 77 (19.7) | 47 (16.3) | 26 (12.8) | 9 (6.9) | 3 (2.2) |
| 5 | 312 (24.9) | 178 (21) | 91 (16.6) | 27 (9.4) | 7 (4.5) |
We stress that the different TDE rates for varying IMFs here are obtained for systems with the same mass (at t = 10 Gyr). In order to properly compare these results to observations, one should in principle normalize different systems depending on their luminosity in a given band, that necessarily depends on the chosen IMF and age of the stellar population. This is beyond the scope of this work, but accounting for this would necessarily partially affect the TDE rates normalization and time evolution.
3.3 Mass segregation as the driver of the evolution
Before tMS, mass segregation is not yet completed, implying that the cusp is growing about the MBH both for heavier and lighter stellar objects; as a reference, Fig. 2 shows the density profile (for BHs, main sequence stars, and total) at the beginning, at the moment the cusp is fully grown (t ≈ 1.8 Gyr) and at the end of the evolution for the run with a Kroupa IMF and Sersic index set to 2. After the cusp is fully grown, the system further relaxes, expanding (thus lowering the normalization of the density) in time as a result of the central MBH acting as an heat source (Vasiliev 2017). The system expansion and density drop is more efficient in systems featuring a top-heavy IMF, as relaxation is more efficient. Fig. 2 also marks the radius from which most TDEs are generated. Once segregation is completed, this radius tends to move outwards, where the density is lower. The event rates thus drop after a mass segregation time-scale as the radius from which the TDE flux is maximum gets larger, and the density at all radii gets lower. The presented picture remains nearly the same for all runs.
Density profile of the system initialized with a Kroupa IMF and with nuclear cluster Sersic index set to 2 (the density evolution in all the other runs is analogue to this). The green dash-dotted line shows the total initial density profile. The blue and red curves respectively refer to the profiles at the moment the density cusp is fully developed, and at the end of the run (14 Gyr). The solid line refers to the density profile of main sequence stars, the dotted line refers to the profile of compact stellar remnants with mass above 2|$\, \mathrm{M}{\odot {}}$| (thus mainly BHs), and the upper edge of the shaded areas show the total density profile. The vertical, thin solid lines mark the radii at which the TDE flux is maximum. The black segments show as a reference the power laws ρ∝r−1.75 (top) and ρ∝r−1.35 (bottom). The density profile for main sequence stars and BHs respectively follow those power law trends in the inner region to a good approximation.
In all explored cases, the late time event rates are very similar if one assumes the same IMF, regardless of the initial cluster profile. This is due to the fact that the system density profile relaxes to a very similar density profile, regardless of the initial condition, and whose normalization slowly drops as time passes; if the system is initially more concentrated, its initial TDE rates are larger, but the system expands more efficiently, so that the two effects compensate each other and the final rates are similar.
3.4 Preference for the disruption of more or less massive stars
The complete mass function evolved in our runs allows us to track which stars are more likely to be disrupted as a function of their mass. An educated guess would be that stars of a given (initial, or equivalently here, main sequence) mass are destroyed proportionally to their IMF, i.e. |$\dot{N}(m)/\dot{N}\propto \chi (m)$|, where |$\dot{N}(m)$| is the TDE rate for main sequence stars with mass m, |$\dot{N}$| is the total TDE rate [|$\dot{N}=\int \dot{N}(m)\mathrm{ d}m$|], and again χ represents the stellar mass function. Based on this, Fig. 3 (top panel) shows the ratio |$\tfrac{\dot{N}(m)}{\dot{N}}/\chi (m)$|;6 we show this ratio for the early (t < 1 Gyr) and late (t = 14 Gyr) evolution.7 The plotted ratio remains close to 1 for all values of m, and it is interesting to notice that the shape of the residual curve is very similar to the trend of the tidal disruption radius rt(m), also plotted in the figure with a dotted black line. This suggests that the value of rt associated to each mass m also has an impact on the TDE rates.
The plots shows the probability of disrupting a star of mass m, normalized to the number of stars with that given m in the IMF. We show this value for TDEs occurring over the first Gyr (solid line) and in the final snapshot (14 Gyr, dashed line); in the former case, we show an extrapolation for large masses with a dotted line – see footnotes 7,8 for more details; in the latter, we set the event rates for stars with |$m\gt 1\, \mathrm{M}{\odot {}}$| to 0 as those would have become WDs. For reference, we show the tidal disruption radius rt arbitrarily normalized to its value for a star of |$0.5 \, \mathrm{M}{\odot {}}$| with a dotted black line. The top plot refers to the run with TDE radius defined via equation (2); the bottom one assumes η = 1. Both plots show the runs with a Kroupa IMF and Sersic index equal to 2.
In order to better explore this aspect, we also perform an additional run in which we fix η = 1 (the value that defines the TDE radius, see equation 1) instead of allowing it to vary with m (as in the standard runs, equation 2). The case with η = 1 is shown in the bottom panel of Fig. 3. Again, the residual curve resembles the trend of rt, confirming that the choice of the TDE radius has an impact on the number of events to be expected for each given m.8 A proper definition of rt is thus crucial to properly compare observations and theoretical work. We checked that the results presented in this subsection remain very similar if one varies the IMF or the Sersic index of the nuclear stellar cluster.
In the presented runs, the value of q at the energy at which the TDE flux is maximum is always greater than 1, implying that most TDEs are produced in the full loss cone regime. Depending on the choice of the IMF, the value of q is ≈2–10 at the beginning of the evolution and it reaches ≈10–20 towards the end, since the maximum flux of TDEs progressively gets to lower binding energies (thus larger radii). It can be shown (see. e.g. Stone et al. 2020 and references therein) that the flux in the loss cone scales nearly as ∝1/ln (llc/l0) if q ≪ 1, and it is thus very little dependent on the size of the loss cone; instead, it scales nearly as ∝(llc/l0)2 if q ≫ 1, so that the loss cone size can have a non-negligible impact on TDE rates. Since in the studied system most TDEs came from the full loss cone, the residual curve shown in Fig. 3 resembles the shape of rt(m).
3.5 Preference for the disruption of stars with |$m\gt 1.3 \, \mathrm{M}{\odot {}}$|
Recently, Mockler et al. (2022) were able to constrain the nitrogen-to-carbon abundance ratios for a subset of the observed TDEs. Since only stars initially more massive than |$\approx 1.3 \, \mathrm{M}{\odot {}}$| undergo the CNO cycle, these moderately massive stars are nitrogen-rich and carbon-deficient, and they can be spotted via spectral analysis. Considering the galaxy sample presented by French et al. (2020), they find that 2 of 5 TDE candidates (40 per cent) in E+A galaxies and 3 of 13 (23 per cent) in quiescent Balmer-strong galaxies have enhanced nitrogen abundances; this is a strict lower limit, since part of the TDEs in the sample do not have C/N measurements.
In order to compare our results with their findings, Fig. 4 shows the fraction of TDEs coming from stars more massive than 1.3|$\, \mathrm{M}{\odot {}}$|, obtained from an analytical estimate based only on the mass function and age (equation 4); in addition, we show the same fraction as extracted from the simulations – with a further large mass correction. Specifically, on the top of the early time mass corrections described in Section 2.4.4, we further correct for the fact that massive stars are progressively more likely to be destroyed due to their larger radii, according to the |$m\gt 2\, \mathrm{M}{\odot {}}$| extrapolations presented in Section 3.4 (computed for all IMFs).
The plots show the probability of disrupting a star with initial mass above 1.3 |$\, \mathrm{M}{\odot {}}$| (i.e. a star that is expected to exibit a C/N overabundance in the TDE spectrum) as a function of time for different choices of the IMF, as colour-coded in the legend. The dotted lines are obtained by a simple estimate based on the IMF and on the time needed by each star to become a WD. The solid thin lines show the results as extracted from the Fokker-Plank runs, while the solid thick lines show the same curves corrected based on the results of Section 3.4, i.e. accounting for the fact that heavy stars – which are not self-consistently simulated – have larger stellar radii and are thus more likely to be disrupted (see e.g. footnotes 7,8); this thick line should be considered the reference one to make comparisons with the observational data. Both panels show runs with nuclear cluster Sersic index set to 2. The top plot shows results from the default runs with TDE radius defined via equation (2); the bottom one assumes η = 1. The step-like shape of the lines is due to the fact that the simulations feature discrete mass families instead of a continuous stellar mass function.
The fraction of relatively massive disrupted stars is very close to the analytical estimate (top panel, dotted) if rt is defined via equation (2), i.e. η = η(m); this is due to the fact that in this case massive stars are disrupted nearly as predicted by the IMF, and the ratio in the top panel of Fig. 3 is very close to 1. However, if one fixes η = 1 (bottom panel, Fig. 4) the fraction of disrupted stars above |$1.3\, \mathrm{M}{\odot {}}$| gets significantly larger compared with the analytical estimate; this is due to the fact that if η = 1, stars above 1.3 |$\, \mathrm{M}{\odot {}}$| are significantly more likely to be disrupted compared with what would be estimated based on the IMF only (bottom panel of Fig. 3).
In general, we can conclude that more top-heavy IMFs can result in up to ≈15 (25) per cent of TDEs to come from moderately massive stars over the first hundreds of Myr if η = η(m) (η = 1). The results by Mockler et al. (2022) suggest that the observed fraction is likely larger; however, it is important to keep in mind that currently their sample is very small and thus subject to large statistical errors. Notably, the results also imply that a different prescription for rt can significantly affect the probability of disrupting moderately massive stars. 9
3.6 TDE rates in case of an extended central star formation
In order to better constrain the effect of star formation, instead of assuming the entire stellar cluster to be formed instantaneously at t = 0, here we test the effect of an extended star formation prescription. We take the snapshot at t = 5 Gyr of the nuclear cluster simulations described above (whose rates are shown in Fig. 1), and in particular we refer to the cases with a Kroupa and a top-heavy, α = 1.7 IMF, both with Sersic index equal to 3. We switch on star formation in the inner ∼6 pc for 200 Myr, growing the central cluster mass by 30 per cent, and then we quench star formation (see Section 2.2 for more technical details). The obtained rates are shown in Fig. 5. Even if only 23 per cent of the new mass of the cluster comes from the recent star formation, the TDE rates grow by a factor of a few (× 2–4) compared to the values right before it started; the impact is more dramatic in the case of a top-heavy IMF. Not only the event rate is boosted by TDEs coming from newly formed stars, but also the former population gets its rate enhanced as star formation occurs, as shown in Fig. 5. The fraction of TDEs coming from the newly formed stars with respect to the total rate for the top-heavy (Kroupa) case reaches a maximum of about 60 (45) per cent near the peak after star formation, and declines to nearly 39 (41) per cent by the end of the evolution.
TDE rates for the Kroupa (blue) and top-heavy, α = 1.7 case (green) if star formation is switched on at t = 5 Gyr, grows the mass of the central cluster by 30 per cent and is quenched after 200 Myr. The TDE rates before 5 Gyr (dotted lines) are the same shown in Fig. 1; the dashed lines show the event rate during star formation, and the solid lines show it after star formation is quenched. In all cases, the shown rate comes from the full stellar population; after star formation, the lower edge of the shaded area shows the event rate coming from the former population as a function of time: this alone also grows as a result of star formation.
The TDE rate boost necessarily depends on the choice of the fraction of newly formed stars, on where they form and so on. Still, the presented results clearly show that even assuming only a fraction of the cluster to be formed in the starburst leads to a significant growth in TDE rates, as the perturbed former population has its rates enhanced while on its way to a new equilibrium state.
4 SUMMARY AND DISCUSSION
This paper explores the effect of a complete, evolved stellar mass function on the rate of TDEs in a nucleated galaxy analogous to observed TDE hosts. We study the event rates in systems (i) featuring a Kroupa and a set of top-heavy IMFs, (ii) featuring different central cluster concentrations (iii) assuming an instantaneous and a more extended star formation. Below we summarize our key results.
The bulge virtually does not contribute to TDE rates.
Assuming that the nuclear cluster is present and it instantaneously forms at t = 0, we find that the TDE rate drops significantly (by a factor ∼10–100, Fig. 1) over 14 Gyr if one accounts for a complete mass function, while the drop is only a factor of a few if we consider a population of 1 |$\, \mathrm{M}{\odot {}}$| stars. The rate drop is more pronounced if one starts with an initially more concentrated stellar cluster, and if the IMF is more top-heavy. The ratio between the early and the late event rates (Table 1) can explain the overabundance of post-starburst galaxies among TDE hosts, especially if a top-heavy IMF is assumed and/or the initial density profile of the central cluster is more concentrated.
The initial peak in the event rates starts declining after a mass segregation time-scale. Within this time-frame, both mass segregation occurs and a Bahcall & Wolf cusp develops near the MBH (Fig. 2). Afterwards, the event rate declines as a result of the stellar system expansion induced by two-body relaxation, that brings to a drop in the density normalization, while the radius from which most TDEs are produced grows larger.
The fraction of TDEs coming from stars of different mass can be described by the stellar mass function to the 0-th order. However, most of the destroyed stars come from the full loss cone regime, for which the event rate scales nearly linearly with the tidal disruption radius rt(m), which grows with the stellar mass m; it follows that there is a residual modulation in the fraction of destroyed stars as a function of m, that resembles the trend of rt(m) (Fig. 3). Depending on the definition of rt(m), the number of destroyed stars as a function of their mass can vary significantly.
We find that, a few hundreds of Myr after the starburst, up to ≈15–25 per cent of destroyed stars have |$m\gt 1.3\, \mathrm{M}{\odot {}}$| (Fig. 4). This fraction is lower compared with what has been observationally constrained by Mockler et al. (2022), but still relatively close to it; varying the definition of rt(m) significantly affects this estimate.
If we assume only 30 per cent of the cluster mass to be formed in a 200-Myr star formation episode (happening when the older cluster stellar population is already relaxed), the TDE rate grows by a factor of a few after the starburst. A sensible part of the enhancement comes from the older stellar population, whose equilibrium state is perturbed by the freshly injected mass.
One of our key results is that accounting for a complete stellar mass function can explain the post-starburst preference of TDE hosts. The alternative, currently most widespread explanation is that post-starburst galaxies are initially born with an extremely steep density profile (Stone et al. 2018). Fig. 6 compares the initial conditions adopted in our runs with the ultrasteep initial conditions adopted by Stone et al. (2018) to explain the post-starburst preference, rescaled to our MBH mass according to their footnote 12 (see also the caption of our Fig. 6). The central density in our and their profile differs by several orders of magnitude, and in general it is not straightforward to imagine how star formation could generate such an extremely steep cusp. In fact, Sanders (1998) points out that, too close to the MBH, molecular clouds should get disrupted by its gravitational field preventing star formation; even if the Milky Way nucleus features a population of young stars very close to the MBH (Morris 1993; Ghez et al. 2003; Lu et al. 2013), supporting instead the idea of possible central star formation, it seems unlikely it can develop such steep profiles. For instance, in these systems the time-scale for stellar collisions and mergers can be shorter than the relaxation time-scale at small radii (<10−2 pc, Stone et al. 2018). Note that if the density profile is initially very steep, relaxation processes render it milder in time, as the equilibrium configuration is the Bahcall & Wolf (1976) solution (ρ∝r−1.75): this means that the ultrasteep configuration should be in place since the starburst stage, i.e. stars should be basically born with ρ∝r−γ, γ = 2.5–2.75. In order to compare our results with the ones by Stone et al. (2018), we evolve their (monochromatic) profiles shown in Fig. 6 with ρ∝r−2.5, r−2.75 and we compare the obtained rates with our complete mass function models in Fig. 7. Even if their rates are overall larger, the early to late event rate is compatible with what we find for much milder profiles featuring a complete mass function; in fact, the ratio between the average TDE rate in the range 250–750 Myr (i.e. roughly the age of post-starburst systems) and 10–12 Gyr is equal to 8–14 for the Stone et al. (2018) cases, i.e. they are compatible with or smaller than the same ratios computed in our runs (see numbers in parenthesis in Table 1).
Comparison between our initial conditions based on a Milky-Way like model (black lines) with those employed by Stone et al. (2018). They explore a range of cusps ρ∝r−γ, where γ = 1.75–2.75, and they conclude that the post-starburst overrepresentation in their monochromatic runs can be explained if γ > 2.5 (red area); their density normalization has been set via ρ(1 pc|$)=0.617^\gamma \times 10^{6.32}\, \mathrm{M}{\odot {}}\, {\rm pc}^{-3} (M_\bullet /10^6\, \mathrm{M}{\odot {}})^{-1.05+0.425\gamma }$| (see their footnote 12). The density profiles that justify the post-starburst preference in Stone et al. (2018) are several orders of magnitude denser in the inner regions compared with our models.
Comparison between the TDE rates obtained by Stone et al. (2018) with the ultrasteep cusp (pink, thick lines featuring the highest rates; solid line: ρ∝r−2.5; dotted line: ρ∝r−2.75) and our TDE rates (our rates are also shown in Fig. 1, where it is easier to distinguish each different line). The different IMFs adopted in our runs are colour-coded as in the legend. The solid, dashed, and dotted lines refer to our models with Sersic index respectively equal to 2, 3, 5. The shaded regions mark where the average early (250–750 Myr) and late (10–12 Gyr) TDE rates are computed in Table 1 – numbers in parenthesis. For comparison, in the same regions, the TDE enhancements for the Stone et al. (2018) models are respectively equal to 8, 14 if γ = 2.5, 2.75.
It is also worth stressing that the two-body relaxation-driven expansion that results in a more or less prompt rate drop acts in both the complete mass function and monochromatic case; in fact, it ultimately represents the driver of the rate decline in the ultrasteep cusps of Stone et al. (2018). Expansion starts later and is less efficient in the monochromatic cases, but Fig. 1 clearly shows that the TDE rates start declining in a nearly power-law fashion at late times even in the single-mass runs. Since lower mass MBHs are likely to be embedded in systems with shorter relaxation times, it is reasonable to expect that these systems would relax over a time-scale significantly shorter than a Hubble time even in the monochromatic case, and thus start expanding and initiate their rate decline much earlier. As a consequence, in systems featuring low mass MBHs, the ratio between the early and late TDE rate could be larger than what found here for the single mass scenario. None the less, a complete mass function would likely enhance the rate drop even about lower mass MBHs; in addition, we stress that the post-starburst galaxies that host observed TDEs harbour an MBH whose mass estimate is compatible with the one chosen in this work (French et al. 2020).
Our results support the idea that starbursting systems are likely characterized by at least a moderately top-heavy IMF, as already suggested in the literature (Zhang et al. 2018). In fact, our computed early to late event rates are similar to the observations-based post-starburst overrepresentation (15–30) if α ≲ 1.7. Furthermore, a top-heavy IMF can help explain the large fraction of moderately massive disrupted stars in post-starburst galaxies (Mockler et al. 2022, Fig. 4). Interestingly, this implies that in principle TDE observations can be adopted to constrain the IMF of starbursting galaxies as more TDE observations become available. Another interesting point concerns the choice of the prescription for computing the tidal disruption radius rt = η(m/M•)1/3⋆. Changing the prescription of η, the probability of disrupting stars with different masses changes significantly. This means that more work is needed to best constrain the functional form of η, as our results suggest that e.g. accounting for partial TDEs can affect the relative probability of disrupting more or less massive stars.
Concerning the overall estimate of TDE rates, theoretical studies (typically assuming a monochromatic mass function) tend to predict values that are generally larger than the observed ones (respectively a few × 10−4 versus 0.3 − 1 × 10−4 events per galaxy per year, see e.g. Stone & Metzger 2016; French et al. 2020); this discrepancy is worsened by the post-starburst overrepresentation, that makes observed rates for quiescent galaxies lower by a factor ∼3. Interestingly enough, the ratio between the late event rates of a monochromatic and a complete mass function |$\dot{N}_{\rm mc}/\dot{N}_{\rm MF}$| varies significantly depending the chosen IMF. In the Kroupa case, this ratio is very close to 1 (|$\dot{N}_{\rm mc}/\dot{N}_{\rm MF}=0.75-0.86$|); however, choosing a top-heavier IMF makes the late time event rates much smaller than their corresponding value for a monochromatic mass function. In fact, |$\dot{N}_{\rm mc}/\dot{N}_{\rm MF}\approx 6,3,1.5$| if α = 1.5, 1.7, 1.9 implying that assuming, in general, a top heavy IMF to be in place in galactic nuclei would help us to reconcile the observed and detected TDE rates. This obviously would work only if all galaxies dominating TDE rates (e.g. more concentrated ones, featuring a nuclear star cluster and a relatively low MBH mass) are characterized by a top-heavy IMF. To better constrain this aspect, it would be important to adopt an observationally motivated galaxy and MBH sample, in order to assess in more detail the statistics of TDE rates when assuming different complete mass functions. This will be the subject of a future work.
To conclude, it is worth mentioning that one important assumption introduced in most of our runs is that the whole galactic nucleus is formed at t = 0. It may seem unrealistic the whole nuclear cluster is formed instantaneously; however, this assumption would be compatible with the mass fraction generated by the starburst in a post-starburst galaxy (Section 1, also note that the same assumption of instantaneous formation is made in Stone et al. 2018). The main process that leads to the formation of nuclear star clusters is still a debated topic, and there are two main proposed formation channels: in the in-situ (wet) scenario, the nuclear cluster is generated via efficient star formation occurring at the galaxy centre (e.g. Loose, Kruegel & Tutukov 1982); in the migration scenario, stellar clusters formed elsewhere migrate via dynamical friction in the galaxy centre, to form a nuclear overdensity (e.g. Tremaine, Ostriker & Spitzer 1975; Capuzzo-Dolcetta 1993).10 In the in-situ framweork, the instantaneous approximation for the cluster formation is reasonable as the star formation time-scale in post-starburst systems can be as low as a few tens of Myr (French et al. 2018), much shorter than the typical galaxy evolution time-scales. Even in the migration scenario, the sudden perturbation in the potential caused by a cluster inspiral would result in a strong dynamical perturbation for the central orbits, so that it is reasonable to assume the system instantaneously forms at that moment. More in general, since mass segregation is in fact the main trigger of the burst in TDE rates, any event violent enough to strongly perturb the central galaxy potential well and quash the previously established, mass-segregated cusp would likely result in a strong rate enhancement. Phenomena like a significant inflow of gas possibly resulting in a strong burst of star formation, a galaxy merger (Gualandris & Merritt 2012), the inspiral of a massive cluster near the nucleus (Bortolas, Mapelli & Spera 2018; Arca-Sedda et al. 2019), could thus induce a TDE burst. Addressing in detail these aspects would require different strategies than those adopted in this work, and we plan to explore this in a forthcoming study.
ACKNOWLEDGEMENTS
We warmly thank the anonymous referee for their very useful comments and suggestions. EB warmly thanks Matteo Bonetti, Luca Broggi, Alberto Sesana, and Nicholas Stone for fruitful discussion, and Eugene Vasiliev for the support in the use of Phaseflow. EB acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program ERC-2018-COG under grant agreement N. 818691 (B Massive). EB also acknowledges support from the University of Milano-Bicocca through the ‘Giovani Talenti’ prize.
DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonable request to the corresponding author.
Footnotes
The supernova natal kick experienced by NSs and BHs can also funnel them within r• (Bortolas & Mapelli 2019), but we did not account for this possibility in this study.
We performed a convergence study varying the number of mass bins, and we found our choice to differ by only 0.1 per cent in TDE rates compared with a case with a sampling of 500 bins.
Foote, Generozov & Madigan (2020) instead focus on the effect of mass segregation and TDE rate enhancement in eccentric nuclear discs with a bi-chromatic mass function, but they do not explore the more general spherical and isotropic scenario and the effect of a complete mass function.
The introduced static estimate is valid in the empty loss cone regime, and it neglects the fact that stars have a different rt depending on their mass. Kochanek (2016) introduce a more elaborate static correction for the presence of a complete mass function, based on the previous works by Wang & Merritt (2004) and Magorrian & Tremaine (1999); in particular, according to their equation (3), the rate enhancement compared with a monochromatic mass function would be
The rate enhancements computed via the equation above within our stellar evolution framework are very modest, in the range 1–1.6, and they are in reasonable agreement with what found by Kochanek (2016); however, they fail to reproduce the rate enhancements obtained with Phaseflow even at the very beginning of the integration.
During the 14 Gyr of evolution, the MBH reaches a mass of |$4.4-5.1\times 10^6\, \mathrm{M}{\odot {}}$| depending on the chosen initial conditions, with a more efficient growth in case of higher initial Sersic index and lower value of α.
Here χ(m) has been normalized to 1 between the minimum allowed mass and |$m_{\rm break} = 2\, \mathrm{M}{\odot {}}$|, in order to avoid accounting for compact objects that cannot undergo a TDE.
We also extrapolate the ratio for the large stellar masses that would be present at early times by performing a linear fit of the shown data between |$m=1.3\, \mathrm{M}{\odot {}}$| and |$m=2\, \mathrm{M}{\odot {}}$|, obtaining Ratio|$(m)=0.85\, \mathrm{M}{\odot {}}+0.11m$| (this will be used in the next subsection).
In this case, the large mass linear regression between 1.3 and 2 |$\, \mathrm{M}{\odot {}}$| results in Ratio|$(m)=1.2\, \mathrm{M}{\odot {}}+0.14m$|.
This also implies that assuming stellar radii to be fixed in time may impact the relative rates as a function of the stellar mass, since stellar radii can vary non-negligibly along the main sequence (e.g. Bressan et al. 2012).
More recently, a combined, wet migration scenario has also been proposed (Guillard, Emsellem & Renaud 2016).
