Tidal disruption events in the first billion years of a galaxy
Hugo Pfister, Jane Dai, Marta Volonteri, Katie Auchettl, Maxime Trebitsch, Enrico Ramirez-Ruiz
MMNRAS , 1–13 (2019) Preprint 12 June 2020 Compiled using MNRAS L A TEX style file v3.0
Tidal disruption events in the first billion years of a galaxy
Hugo Pfister, , (cid:63) Jane Dai, , Marta Volonteri, Katie Auchettl, , , , Maxime Trebitsch , and Enrico Ramirez-Ruiz , DARK, Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, DK-2100 Copenhagen, Denmark Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong, China Institut d’Astrophysique de Paris, Sorbonne Universit´e, CNRS, UMR7095, 98bis boulevard Arago, F-75014, Paris, France Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA School of Physics, The University of Melbourne, Parkville, VIC 3010, Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) Max-Planck-Institut f¨ur Astronomie, K¨onigstuhl 17, 69117 Heidelberg, Germany Zentrum f¨ur Astronomie der Universit¨at Heidelberg, Institut f¨ur Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Accretion of stars on massive black holes (MBHs) can feed MBHs and generate tidaldisruption events (TDEs). We introduce a new physically motivated model to self-consistently treat TDEs in cosmological simulations, and apply it to the assembly ofa galaxy with final mass × M (cid:12) at z = . This galaxy exhibits a TDE rate of ∼ − yr − , consistent with local observations but already in place when the Universewas one billion year old. A fraction of the disrupted stars participate in the growth ofMBHs, dominating it until the MBH reaches mass ∼ × M (cid:12) , but their contributionthen becomes negligible compared to gas. TDEs could be a viable mechanism to growlight MBH seeds, but fewer TDEs are expected when the MBH becomes sufficientlymassive to reach the luminosity of, and be detected as, an active galactic nucleus.Galaxy mergers bring multiple MBHs in the galaxy, resulting in an enhancement ofthe global TDE rate in the galaxy by ∼ order of magnitude during
100 Myr aroundmergers. This enhancement is not on the central MBH, but caused by the presenceof MBHs in the infalling galaxies. This is the first self consistent study of TDEs ina cosmological environment and highlights that accretion of stars and TDEs are anatural process occurring in a Milky Way-mass galaxy at early cosmic times.
Key words: transients: tidal disruption events – quasars: supermassive black holes– galaxies: evolution – galaxies: high-redshift
Gas falling onto black holes, increasing their mass and re-leasing gravitational energy, can explain the growth of mas-sive black holes (MBHs) with masses (cid:38) M (cid:12) in the centerof most massive galaxies (Kormendy & Ho 2013). However,the energy released from gas falling onto MBHs as well asfrom nearby supernovae (also known as “feedback”) can heatand eject gas, preventing MBH growth in low-mass galax-ies ( e.g. Dubois et al. 2015; Habouzit et al. 2017; Trebitschet al. 2018). As a result, it is challenging to explain, with gasaccretion only, the presence of MBHs with masses (cid:38) M (cid:12) at z > when the Universe was only 1 Gyr (Tenneti et al.2017; Ba˜nados et al. 2018).However, the vicinity of MBHs is not only composedof gas, but also of stars (Sch¨odel et al. 2018) that can also (cid:63) Sophie and Tycho Brahe Fellow; hugo.pfi[email protected] be accreted by MBHs and increase their mass. Furthermore,contrary to gas, stars are not subject to feedback, thereforethey could provide a continuous source of material for MBHsto accrete and grow (Alexander & Bar-Or 2017). In order toknow the contribution of stars to the growth of MBHs, oneneeds to know the rate at which stars get close enough toa MBH to be swallowed, whole or in part, and increase itsmass. For a Solar-like star and MBHs with mass (cid:46) M (cid:12) ,stars are not swallowed whole but tidally disrupted, produc-ing a unique signature known as a tidal disruption event(TDE; Lacy et al. 1982; Rees 1988). This allows us to ob-servationally measure the TDE rate, thus providing an esti-mate of the rate at which stars get close enough to MBHsto increase their masses.With a handful of observed TDEs, for central massiveMBHs in quiescent galaxies at z = , the typical rate is ( . − . ) × − yr − (Donley et al. 2002; Gezari et al. 2008;van Velzen & Farrar 2014; Holoien et al. 2016; Blagorod- © a r X i v : . [ a s t r o - ph . GA ] J un H. Pfister et al. nova et al. 2017; Auchettl et al. 2018; van Velzen 2018).This average rate can be well understood theoretically withthe analytical loss-cone theory (Lightman & Shapiro 1977;Magorrian & Tremaine 1999; Wang & Merritt 2004): a MBHwith mass M • is embedded in a stellar density profile andfed by stars scattered toward its direction through 2-bodyinteractions. This formalism also predicts that there is a neg-ative correlation between the TDE rate and the mass of theMBH (the TDE rate scales as M − (cid:15) • with < (cid:15) < . , Wang& Merritt 2004; Stone & Metzger 2016; Pfister et al. 2020),which is confirmed by observations for MBHs with a mass (cid:38) M (cid:12) (van Velzen 2018).This justifies the use of the the loss-cone theory to esti-mate the growth of MBHs through TDEs, and Alexander &Bar-Or (2017) showed that MBHs can reach masses above ∼ × M (cid:12) regardless the initial MBH seed mass and red-shift formation only through accretion of stars. While thissuggests that TDEs are efficient in growing light MBH seeds,this result was derived under idealized assumptions: MBHsare embedded in a singular isothermal sphere (Binney &Tremaine 1987) with an inner Bahcall-Wolf cusp (Bahcall &Wolf 1976), the M • − σ relation (Ferrarese & Merritt 2002;Kormendy & Ho 2013) is always verified, MBHs are fixed inthe center of the stellar distribution and the only relevantprocess is 2-body interactions. To relax these idealized as-sumptions, N-body simulations in which stars getting closeenough to the MBH can be directly counted have been car-ried out (Baumgardt et al. 2004a,b; Brockamp et al. 2011;Zhong et al. 2014). Baumgardt et al. (2004a) and Brock-amp et al. (2011) find that MBHs with masses (cid:38) M (cid:12) can only double their mass within a Hubble time. To sum-marize, the contribution of TDEs to the growth of MBHs isstill uncertain.Furthermore, these previous studies assume that MBHsand their host galaxies are isolated, while most galaxies un-dergo several mergers during their life ( e.g. Fakhouri et al.2010). These mergers drastically affect galaxies, triggeringstar formation and substantially changing the stellar densityprofile near the central MBHs ( e.g.
Van Wassenhove et al.2014; Capelo et al. 2015). As a consequence, it is natural tobelieve that galaxy mergers affect the TDE rate. In addition,this enhancement of the TDE rate during galaxy mergersis somewhat motivated by observations, as E+A galaxieswhich are post-mergers galaxies, are found to have an en-hanced TDE rate of − yr − (French et al. 2016; Stone &van Velzen 2016). To test this, N-body simulations of galaxymergers have been performed ( e.g. Li et al. 2017; Sakuraiet al. 2018), and they indeed find that mergers enhance theTDE rate. There are two reasons to this, and both of themresult in an enhancement of the loss cone feeding: (i) thestellar distribution is overall more triaxal due to the merger;and (ii) when the MBHs get close to each other, stars boundto one MBH see their dynamics greatly perturbed by thecompanion MBH.However, neither the loss-cone formalism nor N-bodysimulations include gas, which can cool and turn into newstars that can then be disrupted. Therefore these frame-works cannot provide a fully consistent picture. To partiallyovercome this issue, Pfister et al. (2019a) adopted a trade-off between the ability to estimate exactly the TDE rateby resolving stars getting close enough to MBHs and in-cluding the physics of galaxies (star formation, supernovae etc...): they used an isolated hydrodynamical simulation ofa galaxy merger starting from idealized initial conditions,and post-processed the TDE rate applying loss-cone the-ory (Vasiliev 2017, 2019) onto the self-consistently evolvingdensity profiles. They found that, indeed, during mergers,nuclear starbursts around MBHs enhance the central stellardensity, naturally resulting in an enhancement of the TDErate (Stone & van Velzen 2016).Furthermore, as the TDE rate results from a combi-nation of the properties of the MBH and its surroundingstellar density profile, it is natural that it varies from galaxyto galaxy (French et al. 2020b). We mentioned the enhance-ment in post-mergers E+A galaxies, but ultra-luminous in-frared galaxies could have a TDE rate as high as − yr − (Tadhunter et al. 2017; Kool et al. 2020), and high redshiftgalaxies which are more star forming and compact (e.g.,Madau & Dickinson 2014; Allen et al. 2017) as well activegalactic nuclei (AGNs) could also exhibit a different TDErate.The current status of the field hints to a diversity ofTDE rates based on galaxy properties, environment andcosmic epoch. The role and importance of stellar accretionon MBH growth and the evolution of the TDE rate mustbe investigated in a fully cosmological context, in whichgalaxies grow over time by accretion of cosmic filaments,where galaxy mergers are numerous, especially at early cos-mic times, and galaxies are more “messy” than in an idealset-up (compare Fig. 3 of Capelo et al. (2015) with Fig. 1).In this paper, we introduce in § × M (cid:12) at z = described in §
3. Thisallows us not only to study the contribution of TDEs tothe growth of MBHs, but also the evolution of the TDErate during mergers and AGN phases, as this galaxy suffersseveral mergers and sometime has an AGN. We discuss ourresults in § § We present how we estimate the TDE rate given the prop-erties around MBHs. We first recall analytical estimates in § § Ramses (Teyssier 2002) and finish by the caveats of our model in § It is customary to express TDEs as sourced by two differentregions (Syer & Ulmer 1999; Merritt 2013): the empty losscone (Wang & Merritt 2004), close to the MBH ( r < r c , r c is defined in the following paragraph), where the diffusiontimescale T r is longer than the radial period; and from thefull loss cone (Pfister et al. 2019a), farther away ( r > r c ),where the diffusion timescale is shorter than the radial pe-riod.We assume a MBH with a mass M • , embedded in astellar density and stellar velocity dispersion profiles ρ and σ , all stars having a mass m (cid:63) and radius r (cid:63) . In this situation, r c is the radius at which the contributions of the full and MNRAS , 1–13 (2019)
DEs under the grid empty loss cone to the flux of stars match, meaning that r c is solution to (Pfister et al. 2019a): G ρ ( r c ) r c σ ( r c ) = r (cid:63) q / (1) ⇔ ρ ( r c ) r c M ( r c ) + M • = r (cid:63) q / , (2)where q = M • / m (cid:63) ; and we have assumed the velocity disper-sion to be σ ( r ) ∼ G ( M ( r ) + M • )/ r , where M ( r ) is the enclosedstellar mass within r .Following Wang & Merritt (2004), we estimate the TDErate coming from the empty loss cone as: Γ empty = M ( r c ) m (cid:63) T r ( r c ) , (3)where T r is the (Spitzer & Harm 1958): T r ( r ) = √ σ ( r ) π G m (cid:63) ρ ( r ) ln ( . M • / m (cid:63) ) . (4)Following Pfister et al. (2019a), we estimate the TDErate coming from the full loss cone as: Γ full = π G q / r (cid:63) ρ ( r c ) σ ( r c ) . (5)The total TDE rate, Γ , can be expressed as the sum ofthe two: Γ = Γ empty + Γ full . (6)To get a step further, we assume the stellar density pro-file to be a power law, with logarithmic slope − < γ ≤ : ρ ( r ) = ρ − γ (cid:18) rr (cid:19) − γ (7) M ( r ) = π r ρ (cid:18) rr (cid:19) − γ , (8)where ρ corresponds to the mean density within r . In thissituation, we can rewrite Eq. (2) as: + ˜ ρ ˜ r − γ c = ρ ˜ r − γ c (9)where we introduce ˜ ρ = ρ / ρ u and ˜ r c = r c / r u , with: ρ u = m (cid:63) / π r (cid:63) (cid:18) r (cid:63) r (cid:19) γ (cid:18) − γ π (cid:19) − γ q ( γ − )/ (10) r u = π − γ q / r (cid:63) . (11)Unfortunately, even for this simple density profile, ingeneral, no explicit expression of r c can be written. However,in some limiting regimes we have: r c ∼ r u (cid:18) ρ ρ u (cid:19) /(− + γ ) if ρ (cid:28) ρ u r u if ρ = ρ u r u if ρ (cid:29) ρ u . (12)allowing us to estimate the TDE rate. Here we detail how we go from the theoretical analysis de-rived in § Ramses .At each timestep of the simulation, between t and t + ∆ t ,we estimate the mean stellar density, ρ , S k ( k is 4 or 8), inthe sphere S k , of radius k ∆ x centered on the MBH ( ∆ x beingthe minimum cell size in the simulation) as: ρ , S k = / π ( k ∆ x ) (cid:213) i ∈S k m i (13) = / π ( k ∆ x ) M k , (14)where m i is the mass of the stellar particle i ; and M k is thetotal stellar mass enclosed within k ∆ x around the MBH. Ifwe assume the density around the MBH to be expressed asgiven in Eq. (7), we obtain: ρ = ρ , S (15) r = ∆ x (16) γ = − ln (cid:18) M M (cid:19) . (17)The mass of the MBH, M • , is measured directly in thesimulation, and we assume stars to be all solar–like, that is r (cid:63) = R (cid:12) and m (cid:63) = M (cid:12) . Finally, we estimate r c as: r c = r u (cid:18) ρ ρ u (cid:19) /(− + γ ) − tanh (cid:16) log (cid:16) ρ ρ u (cid:17)(cid:17) + +
12 tanh (cid:16) log (cid:16) ρ ρ u (cid:17)(cid:17) + + (18) (cid:18) − − + /(− + γ ) (cid:19) exp (cid:18) − log (cid:18) ρ ρ u (cid:19)(cid:19)(cid:21) , which approximates the true value of r c , solution to Eq. (2),within less than 30% error (see Appendix A).With all this, we can estimate the TDE rate onto theMBH, Γ (Eq. (6)), as shown in § Γ m (cid:63) ∆ t is then removed from surrounding starswithin ∆ x and the three following steps are done:(i) A mass (cid:219) M • , star ∆ t = f a Γ m (cid:63) ( − (cid:15) r ) ∆ t is added to theMBH, where f a = . is the fraction of mass which falls ontothe MBH ; and (cid:15) r is the radiative efficiency, which dependson the spin of the MBH (4% for a non rotating MBH andup to 42% for a highly spinning MBH; in this paper we usefixed value of (cid:15) r = , see § (cid:219) M d ∆ t = ( − f a ) Γ m (cid:63) ∆ t does not fall onto theMBH and returns into cells containing disrupted stars as gas(see Eq. (20)). Note that in the paper we clearly make the difference betweenthe TDE rate in yr − corresponding to the number of stars beingdisrupted, and the “stellar accretion rate” (stars are not accreted per se , gas falling back from the disrupted stars is) in M (cid:12) yr − corresponding to the total mass of disrupted stars falling onto theMBH. This difference is mainly “syntactic” as we assumed thatall stars are solar like and f a = . , therefore the stellar accretionrate and the TDE rate differ by a factor of two in their respectiveunits.MNRAS000
12 tanh (cid:16) log (cid:16) ρ ρ u (cid:17)(cid:17) + + (18) (cid:18) − − + /(− + γ ) (cid:19) exp (cid:18) − log (cid:18) ρ ρ u (cid:19)(cid:19)(cid:21) , which approximates the true value of r c , solution to Eq. (2),within less than 30% error (see Appendix A).With all this, we can estimate the TDE rate onto theMBH, Γ (Eq. (6)), as shown in § Γ m (cid:63) ∆ t is then removed from surrounding starswithin ∆ x and the three following steps are done:(i) A mass (cid:219) M • , star ∆ t = f a Γ m (cid:63) ( − (cid:15) r ) ∆ t is added to theMBH, where f a = . is the fraction of mass which falls ontothe MBH ; and (cid:15) r is the radiative efficiency, which dependson the spin of the MBH (4% for a non rotating MBH andup to 42% for a highly spinning MBH; in this paper we usefixed value of (cid:15) r = , see § (cid:219) M d ∆ t = ( − f a ) Γ m (cid:63) ∆ t does not fall onto theMBH and returns into cells containing disrupted stars as gas(see Eq. (20)). Note that in the paper we clearly make the difference betweenthe TDE rate in yr − corresponding to the number of stars beingdisrupted, and the “stellar accretion rate” (stars are not accreted per se , gas falling back from the disrupted stars is) in M (cid:12) yr − corresponding to the total mass of disrupted stars falling onto theMBH. This difference is mainly “syntactic” as we assumed thatall stars are solar like and f a = . , therefore the stellar accretionrate and the TDE rate differ by a factor of two in their respectiveunits.MNRAS000 , 1–13 (2019) H. Pfister et al. (iii) An energy f a Γ m (cid:63) (cid:15) r ∆ tc is emitted by the MBH. Atthe moment, we consider that the feedback is similar for ac-creted stars than for accreted gas: either thermal or kineticdepending on the Eddington ratio (see § Ramses ). Assuming a similar expres-sion for the feedback following gas or stellar accretion is notabsurd, indeed, once the accretion disk is formed followingthe disruption of a star, whether the material was originatedfrom a star or a gaseous clump should not change the behav-ior. Note that this should be a upper limit of the radiativefeedback, since the radiative efficiency likely has a smallervalue than the thin disk one, and a fraction of the boundstellar debris can become outflows.In addition, to conserve total momentum, we updatethe velocity of the MBH, gas and stars accordingly. In theend, we have:stars (cid:40) m i ( t + ∆ t ) = m i ( t ) − Γ m (cid:63) ∆ t f i v i ( t + ∆ t ) = v i ( t ) (19)gas ρ g , i ( t + ∆ t ) = ρ g , i ( t ) + (cid:219) M d ∆ t f i ∆ x u i ( t + ∆ t ) = ρ g , i ( t ) ∆ x u i ( t ) + (cid:219) M d ∆ t f i v i ( t ) ρ g , i ( t ) ∆ x + (cid:219) M d ∆ t f i e i ( t + ∆ t ) = e i ( t ) + (cid:219) M d ∆ t f i ∆ x v i ( t ) (20)MBH v • ( t + ∆ t ) = v • ( t ) M • ( t ) + (cid:104) v (cid:63) (cid:105) (cid:219) M • , star ∆ tM M • ( t ) + (cid:219) M • , star ∆ tM • ( t + ∆ t ) = M • ( t ) + (cid:219) M • , star ∆ t , (21)where (cid:104) v (cid:63) (cid:105) is the mass-weighted velocity of stars, with veloc-ities v i , within ∆ x from the MBH; v • is the velocity of theMBH; f i = m i ( t )/ M is the contribution of the stellar parti-cle i to the TDE rate ( M = M k for k = is the enclosedstellar mass within ∆ x ); and ρ g , i , u i and e i are respectivelythe density, velocity and total energy density of the cell con-taining the stellar particle i . We discuss here a few numerical and physical caveats of theimplementation: • If the MBH is off-center from its host galaxy, and istherefore not in a spherical density profile, Eq. (17) couldgive negative γ . When γ < , we set Γ = . • If the available mass of stars ( M ) is lower than the dis-rupted mass Γ∆ t , then there are not enough stars. In thissituation, we set Γ = M / ∆ t and remove all available stars(note that in practice this did not happen in our simula-tions). • Even if the density profile is spherical around the MBH,it is possible that it does not follow a simple power law.Our “bet” is that, if the resolution of the simulation is highenough, then the estimate of the inner slope γ is enoughfor an estimate of the TDE rate. In practice as our sim-ulation reaches a resolution ∆ x ∼ (see Table. 1), thistranslates into assuming a constant slope within ∼
60 pc forour estimate of the TDE rate. Note that observed galaxiesat z (cid:28) are usually well fitted with fixed inner slope within ∼
100 pc ( e.g. Lauer et al. 2007), and there seem to be a correlation between density at these scales and the TDE rate(French et al. 2020a). However, this excludes the presenceof a nuclear star cluster around MBHs (Pechetti et al. 2019;S´anchez-Janssen et al. 2019) which could enhance the TDErate by orders of magnitude (Pfister et al. 2020). • It is currently not known what is the fraction of dis-rupted material which falls back onto the MBH ( f a ), norhow long it takes. If a star comes with a highly parabolicorbit, i.e. with total energy “close to zero”, then we expecthalf of the debris to remain bound and half to be unbound.We assume here that all bound debris immediately falls backonto the MBH ( f a = . ). • Pfister et al. (2019a) and Wang & Merritt (2004)only give an approximate TDE rate in the full and emptyloss cone regime. More detailed analytical framework ex-ist (Stone & Metzger 2016; Vasiliev 2017), but it would benumerically inefficient (it involves computing numerous “in-tegrals”) and meaningless (we assume a spherical densityprofile and all stars a solar like which are “larger” approxi-mations than the full/empty loss cone) to use them. • Assuming that all stars are all solar–like is clearly sim-plistic, however, Stone & Metzger (2016) have shown thatusing a stellar mass distribution function varies the TDErate by only ∼ with respect to the monochromatic Solarpopulation we consider. • Although stellar accretion can be super–Eddington, westill use the feedback thermal mode (see § (cid:44) adowskiet al. 2016; Dai et al. 2018) which find that at high accre-tion rate, the feedback is more likely to be mechanical andpossibly jetted if the conditions are optimal. We leave thisdevelopment of super–Eddington accretion to a future study. In order to study the evolution of TDE rate in a galaxyevolving in a realistic context, we run a cosmological zoomon a halo whose properties are described in § Ramses (Teyssier 2002).
Ramses follows the evolution of the gas using thesecond-order MUSCL-Hancock scheme for the Euler equa-tions; and the approximate Harten-Lax-Van Leer ContactRiemann solver, with a MinMod total variation diminishingscheme to reconstruct the interpolated variables from theircell-centered values, is used to compute the unsplit Godunovfluxes at cell interfaces (Toro 1997). An equation of state ofperfect gas composed of monoatomic particles with adiabaticindex / is assumed to close the full set of fluid equations.The Courant factor is set to 0.8 to define the timestep.Collisionless particles (dark matter, stars and MBHs)are evolved using a particle-mesh solver with a cloud-in-cell(CIC) interpolation. The size of the CIC is that of the localcell for MBHs and stars. As dark matter (DM) particles arelarger in mass, we smooth their distribution to reduce theircontribution to shot noise, and they can only project theirmass on the grid down to a minimum cell size of ∆ x DM ,corresponding to the highest level unlocked when runningthe DM only simulation with the same mass resolution.The AMR grid is refined using a quasi-Lagrangian crite- MNRAS , 1–13 (2019)
DEs under the grid Name Value Comments L box
59 Mpc Size of the box at z = M vir × at z = . l max
23 Maximum level of refinement of the AMR grid ∆ x ∆ x DM pc Spatial resolution of dark matter m partDM M (cid:12) Mass of high resolution dark matter particles m part (cid:63) × M (cid:12) Mass of stellar particles M • , seed M (cid:12) Seed mass of MBHs
Table 1.
Simulation parameters rion: a cell is refined if M cellDM + ( Ω m / Ω b − ) M cell b ≥ × m partDM ,where M DM and M cell b are respectively the mass of darkmatter and baryons in the cell; Ω m and Ω b are the totalmatter and baryon density and m partDM is the mass of high-resolution dark matter particles. The minimum cell size, ∆ x ,is kept roughly constant in proper physical size with red-shift: an additional level of refinement is added every timethe expansion factor a exp increases by a factor of two, suchthat the maximum level, l max , is reached at a exp = . . Forsimplicity, we further assume that ∆ x = L box / l max , where L box is the size of the box at z = .The subgrid physics is described below in § § The initial conditions are produced with
Music (Hahn &Abel 2013) and are the same as in Trebitsch et al. (2019). Weassume a Λ CDM cosmology with total matter density Ω m = . , baryon density Ω b = . , dark energy density Ω Λ = . , amplitude of the matter power spectrum σ = . , n s = . spectral index and Hubble constant H = .
74 km s − Mpc − consistent with the Planck data (PlanckCollaboration et al. 2016).Low resolution dark matter particles with mass m coarseDM = × M (cid:12) are placed onto the box with an effectiveresolution of elements. Additional high-resolution darkmatter particles, with an effective resolution of ele-ments corresponding to a mass m partDM = M (cid:12) , are placedaround a halo of mass M vir = × M (cid:12) at z = . . Gas is allowed to cool by hydrogen and helium with a contri-bution from metals using cooling curves from Sutherland &Dopita (1993) for temperatures above K . For gas below K and down to our minimum temperature of
10 K , weuse the fitting functions of Rosen & Bregman (1995).The effect of reionization is modelled with a uni-form heating from the UVB background from Haardt &Madau (1996) below z = . . In addition, to take intoaccount that the center of dense regions can be shieldedby neutral hydrogen, the UV photo-heating is reduced by exp (− ρ g / ρ shield ) , where ρ g is the gas density of the cell and ρ shield = .
01 amu cm − . During each timestep ∆ t , in leaf cells with gas den-sity ρ g > − , N stellar particles with mass m part (cid:63) = × M (cid:12) are drawn from a Poisson distribu-tion with parameter λ = M SF / m part (cid:63) , where M SF is the massof newly formed stars (Rasera & Teyssier 2006). M SF is com-puted so that the star formation rate follows a Kennicutt–Schmidt Law (Schmidt 1959; Kennicutt 1998), that is M SF = (cid:15) ρ g ∆ x ∆ t / t ff , where (cid:15) is the star formation efficiency and t ff = (cid:112) π /(
32 G ρ g ) is the free fall time. (cid:15) depends on the local properties of gas and is estimatedusing the multi-ff PN model from Federrath & Klessen(2012).
21% of the mass of each stellar particles is re-emitted inthe medium in supernovae 5 Myr after their formation, re-leasing a (kinetic) energy of × erg M − (cid:12) . The amountof momentum depositted depends on the local density andmetallicity of each neighbouring cell, and depends on thestages of the Sedov-Taylor blast wave (see Kimm & Cen2014). In addition, modifications from Kimm et al. (2017),using the results of Geen et al. (2015), to take into accountpre-heating of the interstellar medium by radiation beforethe supernovae explosion, are used. Our model for MBHs follows closely from Dubois et al.(2012).
MBHs are represented with sink particles, with an initialmass M • , seed = M (cid:12) . They are formed in Jeans unsta-ble cells containing enough gas to form the MBH, and with min ( ρ (cid:63) , ρ gas ) >
100 amu cm − , where ρ (cid:63) ( ρ gas ) corresponds tothe stellar (gas) density in the cell. As this criterion forma-tion is local, i.e. we do use any halo finder to enforce MBHseeding in the exact center of halos/galaxies ( e.g. Vogels-berger et al. 2013), this could result in multiple MBHs pergalaxy. In order to avoid this, an exclusion radius of 50 kpcis used.
Each MBHs are surrounded by massless cloud particlesequally spaced by ∆ x / on a regular grid lattices within asphere of radius ∆ x around the MBH. These cloud particlesare used to measure the averaged gas quantities around theMBH. For instance, the mean gas density is obtained as: ˜ ρ g = (cid:213) i ∈ cloud particles ρ g , i exp (cid:32) − r i r • (cid:33) , (22) This corresponds to the mass fraction of stars more massivethan (cid:12) assuming a Kroupa initial mass function (Kroupa 2001)with stars having a mass in between 0.08 and
100 M (cid:12) .MNRAS000
100 M (cid:12) .MNRAS000 , 1–13 (2019)
H. Pfister et al.
Figure 1. Left:
Stellar density projection.
Right: gas–density–weighted gas density projection. In all cases, the images are centered onthe main galaxy.
Top:
Moment at which the main galaxy undergoes a minor 1:10 merger, the satellite galaxy is in the bottom right ofthe main galaxy and contains the “minor” MBH (red dot).
Bottom:
Moment at which the main galaxy undergoes a major 1:4 merger,the satellite galaxy is on the top right of the main galaxy and contains the “major” MBH (blue dot). In all cases, we show all the MBHs(dots): the central MBH of the main galaxy (black), the minor MBH (red), the major MBH (blue) as well as all the other MBHs in thefield of view (green). Finally, we indicate MBHs which have, at the time of the snapshot, a TDE rate larger than − yr − with a yellowring. The colors used for the 3 “special” MBHs are the same as in Fig. 5. where ρ g , i is the gas density of the cell the cloud particle liesin and r i is the distance of the cloud particle to the MBH. r • is defined as: r • = ∆ x if r B < ∆ x r B if ∆ x < r B < ∆ x ∆ x if ∆ x < r B , (23) where r B = G M • /( c s + v • , g ) is the Bondi radius (Bondi 1952); c s and v • , g are respectively the sound speed and relativevelocity of the MBH with respect to the gas, in the cell theMBH lies in.From these averaged quantities we can estimate the gasaccretion rate (cid:219) M • , gas , using the minimum between the Bondi MNRAS , 1–13 (2019)
DEs under the grid and the Eddington (Eddington 1916) accretion rate: (cid:219) M B = π G M • ˜ ρ g ( ˜ c s + ˜ v • , g ) / (24) (cid:219) M Edd = π G M • m p (cid:15) r σ T c , (25)where m p is the proton mass; c is the speed of light; σ T isthe Thompson cross-section and (cid:15) r is classically fixed to 10%as the spin is not followed in the simulation.In addition, stellar accretion onto MBHs through TDEs,as described in § (cid:219) M • . Following accretion, between t and t + ∆ t , the energy releasedin the medium is E AGN = (cid:15) r (cid:15) f (cid:219) M • c ∆ t , (26)where (cid:15) f is the coupling efficiency, indicating how does theenergy released couples with the gas and depends on themode the MBH is in.At high accretion rate (Eddington ratio, χ = (cid:219) M • / (cid:219) M Edd > ), the (thermal) energy is uniformly distributed in all cellswithin r AGN = ∆ x from the MBH: this is the thermal mode.In this situation we set (cid:15) f = . , lower than the value fromDubois et al. (2012) or Trebitsch et al. (2019), but largerthan Lupi et al. (2019) and similar to Capelo et al. (2015).At low accretion rate ( χ < ) the (kinetic) energy isreleased through a cylindrical bipolar jet centered on theMBH, with radius/height r AGN and direction parallel to theangular momentum of surrounding gas: L g = (cid:213) i ∈ cloud particles ρ g , i r i × u i , (27)where r i and u i are respectively the distance and velocityrelative to the MBH of the gas cell hosting the cloud particle i . The rate at which momentum is deposited depends on theradial distance r to the axis of the cylinder: (cid:219) p Jet ( r ) = ψ ( r ) η (cid:219) M • × (cid:115) (cid:15) r (cid:15) f η c , (28)where (cid:15) f = as in Dubois et al. (2012); η = is themass loading factor, corresponding the the enhancement ofthe mass due to swept up gas , and: ψ ( r ) ∝ exp (cid:32) − r r (cid:33) (29)sums up to 1 over the whole cylinder. Contrary to many simulations where MBHs are anchored tothe center of galaxies ( e.g.
Vogelsberger et al. 2013), we Note that the speed of the jet is km s − with the parameterchosen, whereas in reality jets are relativistic. The difference isdue to our lack of resolution (7 pc) and the jet should instead beconsidered as a wind. allow MBHs to freely move in the potential. Being massive,they suffer dynamical friction (Chandrasekhar 1943; Binney& Tremaine 1987; Tremmel et al. 2015), some of which isunresolved due to lack of resolution (Pfister et al. 2017). Forthis reason, additional forces, in the opposite direction of thevelocity of the MBH, are added to correct the dynamics.Dynamical friction from stars/dark matter is detailedin Pfister et al. (2019b) (using analytical work from Chan-drasekhar 1943), and dynamical friction from gas is detailedin Dubois et al. (2014) (using analytical work from Ostriker1999). To our knowledge, Ramses is currently the only codewhich physically treats both collisional and collisionless un-resolved dynamical friction.Finally, we stress that we have chosen a relatively mas-sive MBH seed ( M • , seed = M (cid:12) > m part (cid:63) ), as such, theseMBHs are not subject to spurious 2-body interactions andno additional correction is needed for the dynamics (Pfisteret al. 2019b). When two MBHs get closer than ∆ x , and if the gravita-tional energy of the binary is larger than the kinetic en-ergy, i.e. the binary would be bound in vacuum, MBHs arenumerically merged. Note that this could lead to spuriousmergers (Volonteri et al. 2020), which we do not explore inthis paper. We use
AdaptaHOP (Aubert et al. 2004) on dark matter(stellar) particles to detect gravitationally bound structures, i.e. halos (galaxies), containing at least 50 particles. We thenconstruct the history of halos (galaxies) using
TreeMaker (Tweed et al. 2009), which match halos (galaxies) from oneoutput to the other using the IDs of particles forming thestructures.We then match galaxies to halos, selecting the closestgalaxy in position. As the zoom has been made on a par-ticular halo, which is the most massive one unpolluted, i.e. containing only high resolution dark matter particles, thegalaxy of this halo is the “main” galaxy. Galaxies which areidentified and are matched to other unpolluted halos arecalled “satellite” galaxies.Finally, we match MBHs to galaxies. A MBH is assumedto belong to a galaxy if it is within the effective radius of thegalaxy (see definition in Appendix B), and the closest to thecenter is the central MBH of this galaxy. If a MBH can beassociated to many galaxies, we assign the MBH the mostmassive galaxy. In what follows, we refer to the “central”MBH as the central MBH of the main galaxy at the end ofour simulation (at z ∼ ).In Fig. 1 we show the stellar (gas) density projection ofthe main galaxy during a minor 1:10 and a major 1:4 merger.We indicate MBHs with dots: the central MBH (black), thecentral MBH of the satellite galaxy of the minor merger (the“minor”MBH in red), the central MBH of the satellite galaxyof the major merger (the “major” MBH in blue) as well asall the other MBHs in the field of view (green). Finally, weindicate MBHs which have, at the time of the snapshot, aTDE rate larger than − yr − with a yellow ring. MNRAS000
TreeMaker (Tweed et al. 2009), which match halos (galaxies) from oneoutput to the other using the IDs of particles forming thestructures.We then match galaxies to halos, selecting the closestgalaxy in position. As the zoom has been made on a par-ticular halo, which is the most massive one unpolluted, i.e. containing only high resolution dark matter particles, thegalaxy of this halo is the “main” galaxy. Galaxies which areidentified and are matched to other unpolluted halos arecalled “satellite” galaxies.Finally, we match MBHs to galaxies. A MBH is assumedto belong to a galaxy if it is within the effective radius of thegalaxy (see definition in Appendix B), and the closest to thecenter is the central MBH of this galaxy. If a MBH can beassociated to many galaxies, we assign the MBH the mostmassive galaxy. In what follows, we refer to the “central”MBH as the central MBH of the main galaxy at the end ofour simulation (at z ∼ ).In Fig. 1 we show the stellar (gas) density projection ofthe main galaxy during a minor 1:10 and a major 1:4 merger.We indicate MBHs with dots: the central MBH (black), thecentral MBH of the satellite galaxy of the minor merger (the“minor”MBH in red), the central MBH of the satellite galaxyof the major merger (the “major” MBH in blue) as well asall the other MBHs in the field of view (green). Finally, weindicate MBHs which have, at the time of the snapshot, aTDE rate larger than − yr − with a yellow ring. MNRAS000 , 1–13 (2019)
H. Pfister et al. M M ,gas M ,star MM /100067891011 z Figure 2.
Masses of the central MBH (green) and the maingalaxy (orange) as a function of time in our simulation (solidlines) and in the simulation of Trebitsch et al. (2019) (dashedlines). We also show the accreted mass of gas (blue) and of starsfollowing TDEs (black). MBH mergers are indicated with dots(this work) or triangles (Trebitsch et al. 2019). In the end, the to-tal contribution of accretion following TDEs is negligible, exceptat early time, where accretion from TDEs and gas is similar.
As we have used the exact same initial conditions asTrebitsch et al. (2019), we can make a fair comparison be-tween the global properties of the two simulations, keepingin mind that details may vary, as some parameters are notexactly the same (particle stellar mass, seed mass of MBHs,use of boost for gas dynamical friction, absence of TDEsetc... see § M (cid:12) , its mass isindependent of the detailed parameters of the simulation.On the same Figure, we show the mass of the centralMBH (green) in the two simulations, as well as the momentsat which the central MBH undergoes a MBH merger (mark-ers). The final masses, which differ by a factor of 3, matchremarkably well considering that (in unranked order) (i) theinitial MBH seed masses are different; (ii) Trebitsch et al.(2019) uses a boost for gas dynamical friction, “encourag-ing” the MBH to remain in gas dense regions and reducingits relative velocity to surrounding gas, enhancing the accre-tion rate (which scales as the density and the inverse cubicof the relative velocity, see Eq.(24)), sometime by orders ofmagnitude; (iii) the number of mergers, and the total “ac-creted” mass through mergers greatly differ: 3 mergers inour simulation corresponding to 6% of the final mass, and20 mergers in Trebitsch et al. (2019) corresponding to 24%of the final mass (this is likely to be related to (i) and (ii) ,but we leave this for future investigations, as we are inter-ested in the TDE rate in this paper); (iv)
Trebitsch et al.(2019) do not include MBH growth through TDEs; and (v) the AGN feedback coupling efficiency in the thermal modediffers by a factor of 10 in the two simulations. M , s t a r M , s t a r + M , g a s z Figure 3.
Ratio of the mass accreted through TDEs ( M • , star ) withthe total mass accreted from gas and stars ( M • , star + M • , gas ), as afunction of time, for the central MBH. We also indicate the lasttime at which fraction is larger than 50% (solid black line), 10%(dashed black line) and 1% (dotted black line). In the end, about . of the mass is gained from TDEs, and their contribution isnegligible for massive MBHs. However, during the first 300 Myr,their contribution is larger than 10%. On the same Figure, we show the mass accreted throughgas (blue), and through stars following TDEs (black). Asthe total contribution of TDEs is only M (cid:12) out of the × M (cid:12) of the MBH final mass, this suggests that thedifference between Trebitsch et al. (2019) and our simulationis not due to (iv) , and including TDEs is not mandatory toproperly estimate the final mass of the MBH. However, atearly time, when the MBH is lighter than ∼ × M (cid:12) , thecontribution of stars appear to be similar to that of gas.In Fig. 3 we show the fraction of mass accreted throughstar as a function of time. At early time, TDEs and theirsubsequent stellar accretion have a significant contributionto the growth of the MBH. Indeed, more than 10% of theaccreted mass of the central MBH is coming from stars dur-ing the first 300 Myr of its life, until its mass is larger than ∼ × M (cid:12) and the MBH is massive enough to accrete atabout the Eddington rate and mostly grow through gas ac-cretion. Unfortunately, for numerical reasons (see § §
50 kyr ; dark color, 10 Myr).Although they are shown with the same frequency, the stel-lar accretion rate is smoother than gas accretion rate. Thereason is twofold: (i) the stellar density is spatially smootherthan the gas density (see Fig. 1 for projections maps), there-fore changes in the MBH position will change the gas den-sity (and MBH gas accretion), leaving the stellar density(and the MBH stellar accretion) unchanged; and (ii) starsare not subject to feedback while gas is, so at a given spatialposition, the stellar density is temporally smoother than thegas density (Prieto et al. 2017). More quantitatively, we sim-ply estimate smoothness of a quantity u as the time average MNRAS , 1–13 (2019)
DEs under the grid M / y r M ,gas M ,star M ,Edd z Figure 4.
Gas (blue), stellar (black) and Eddington (red) accre-tion rate of the central MBH. Light colors are direct outputs ofthe simulation (every
50 kyr ) and dark colors are averaged with a10 Myr window. The stellar accretion rate is strikingly smootherthan the gas accretion rate, although the final contribution of thelatter is larger (see Fig. 2 and Fig. 3). of the relative variation throughout the simulation: u = (cid:28)(cid:12)(cid:12)(cid:12)(cid:12) ∆ uu (cid:12)(cid:12)(cid:12)(cid:12)(cid:29) , (30)where ∆ u is the variation of the quantity u between twoconsecutive timesteps (about 50 kyr); u is the mean valueof the quantity u on two consecutive timesteps and (cid:104) . (cid:105) indi-cates an average over the duration of the simulation. We findthat ( (cid:219) M • , star , (cid:219) M • , gas , ρ , S , ˜ ρ g ) = ( , , , ) . The rela-tive variation in the gas accretion ( (cid:219) M • , gas ) reproduces wellthe relative variations of the the gas density in the vicinityof the MBH ( ˜ ρ g ). While the relative variations of the stellaraccretion ( (cid:219) M • , star ) does not reproduce as well the relativevariations of the stellar density around the MBH ( ρ , S ), werecall that, contrary to gas accretion, stellar accretion doesnot scales directly linearly with the stellar density aroundthe MBH.This confirms however that the rapidly (slowly) varyinggas (stellar) density around the MBH results in a rapidly(slowly) varying gas (stellar) accretion. We note that, be-cause it is much smoother, at any times (hence MBHmasses), the accretion rate following TDEs can be ordersof magnitude larger than the gas accretion rate. This sug-gests that, at any time, it is possible that the properties ofthe emitting MBH are those of a MBH accreting stars only.If the composition of stars differ from the composition ofsurrounding gas ( e.g. stars have a higher nitrogen to carbonabundance), this confirms that, at any time, nitrogen richquasar could be due to TDEs (Kochanek 2016; Liu et al.2018).To summarize, stellar accretion due to TDE is smootherthan gas accretion, simply due to that the stellar density inthe vicinity of the MBH is smoother than the gas density,and stellar accretion can be much larger than gas accretionat all time. However, in the end, growth through TDEs isefficient only for MBHs with a mass lower than × M (cid:12) ,more massive MBHs mostly grow through gas accretion andthe final TDEs contribution is negligible. / y r centralminormajorgal z Figure 5.
TDE rate as a function of time of the central MBHof the main galaxy (black), of the central MBH of the secondarygalaxy during the 1:10 minor merger (red), and of the centralMBH of the secondary galaxy during the 1:4 major merger (blue).These MBHs are shown with the same colors as in Fig. 1. We alsoshow the total TDE rate of the main galaxy (orange). All TDErates are averaged with a 10 Myr window. The two thick blackvertical areas indicate two moments at which the main galaxyundergoes a merger, and for which we show stellar/gas densityprojection maps in Fig. 1. When MBHs are not in the main galaxy,we indicate their evolution with dashed lines, their subsequentevolution following the galaxy merger is marked with solid lines.We find a clear enhancement of ∼ order of magnitude of thetotal TDE rate of the main galaxy during mergers, however, theenhancement does not occur on the central MBH. Our simulation allows us to estimate the TDE rate of everyMBHs as a function of time. Since we also know which MBHsbelong to the main galaxy, we can estimate the total TDErate of the galaxy as: Γ gal = (cid:213) i ∈ BHs in the main galaxy (cid:104) Γ i (cid:105)
10 Myr , (31)where Γ i is the TDE rate of MBH i and (cid:104) . (cid:105)
10 Myr indicates anaverage over a 10 Myr window (our results are unchangedwith a 5 or 50 Myr window).In what follows we will focus on the three “special”MBHs presented in § We show in Fig. 5 the TDE rates of the 3 MBHs Γ central / Γ minor / Γ major in black/red/blue (colors are the sameas the dots representing these MBHs in Fig. 1) as well asthe total TDE rate of the galaxy Γ gal (orange) as a func-tion of time. The minor and major galaxy mergers shown inFig. 1 are indicated with thick vertical black areas. Whenthe MBHs of the satellite galaxies are not in the main galaxy(they are brought by the galaxy merger), we indicate theirevolution with a dashed line.The total TDE rate of the galaxy (orange) is few MNRAS000
10 Myr indicates anaverage over a 10 Myr window (our results are unchangedwith a 5 or 50 Myr window).In what follows we will focus on the three “special”MBHs presented in § We show in Fig. 5 the TDE rates of the 3 MBHs Γ central / Γ minor / Γ major in black/red/blue (colors are the sameas the dots representing these MBHs in Fig. 1) as well asthe total TDE rate of the galaxy Γ gal (orange) as a func-tion of time. The minor and major galaxy mergers shown inFig. 1 are indicated with thick vertical black areas. Whenthe MBHs of the satellite galaxies are not in the main galaxy(they are brought by the galaxy merger), we indicate theirevolution with a dashed line.The total TDE rate of the galaxy (orange) is few MNRAS000 , 1–13 (2019) H. Pfister et al. − yr − . This value is in good agreement with local es-timates (Donley et al. 2002; Gezari et al. 2008; van Velzen& Farrar 2014; Holoien et al. 2016; Blagorodnova et al. 2017;Auchettl et al. 2018; van Velzen 2018) but already in placeat z (cid:38) . We recall that, by construction only one fairlymassive galaxy is studied here (this is a zoom-in simula-tion), and a more statistical analysis should be performed,but this suggests that some galaxies could already have awell established TDE rate of few − yr − at z ∼ when theuniverse is 1 Gyr.Initially, the total TDE rate of the galaxy (orange) issimilar to that of the central MBH (black), i.e. the TDE rateof the galaxy is dominated by TDEs occuring on the centralMBH. However, MBHs brought by successive mergers (allthe dots but the black one in Fig. 1), which can take verylong time to sink toward the center of the galaxy throughdynamical friction (Pfister et al. 2019b), also contribute tothe total TDE rate of the galaxy, sometime dominating it.For instance, during the first minor merger we consider(at t = .
73 Gyr ), the MBH of the satellite galaxy (the minorMBH in red), which has a high TDE rate ( × − yr − )penetrates the main galaxy, resulting in an enhancement thetotal TDE rate. This high TDE rate around the minor MBHis due to a merger induced nuclear starburst at t = .
70 Gyr ,time at which the star formation rate within ∆ x =
28 pc from the minor MBH is enhanced by 30. This picture is inagreement with previous theoretical results who find thatmergers trigger nuclear starbursts, enhancing the TDE rate(Pfister et al. 2019a).During the second major merger we study (at t = .
90 Gyr ) the major MBH (blue) penetrates the main galaxyand completely dominates the rate. The picture here is how-ever different than that of the first merger, as the TDE ratearound this major MBH is not enhanced per se : it was of × − yr − since t = .
60 Gyr . Instead, the major MBHpenetrates the main galaxy while being surrounded by analready dense stellar cusp (see bottom left panel of Fig. 1),therefore its already high TDE rate is not affected.Overall, we find that during the two mergers we dis-cussed, the TDE rate is enhanced by 1 order of magnitudeduring about 100 Myr. This enhancement is due to a nuclearstarburst for the first minor merger, and to that a MBH witha well established stellar cusp enters the main galaxy for thesecond major merger. Other processes resulting in an en-hancement of the TDE rate could happen during mergers:dynamical effects in dry mergers (Li et al. 2017); or simplya MBH on an eccentric orbit periodically crossing the densecenter of the main galaxy. We did not find such examples inour simulation.Finally, we note that the TDE rate of the central galaxyis dominated by off-centered MBHs during about 200 Myrout of the 1 Gyr our simulation lasts, suggesting that duringup to 20% of the time, the TDE rate could be dominatedby off-centered TDEs. While surveys designed to find TDEs( e.g. van Velzen et al. 2020) usually look for central TDEsto exclude most supernovae, blind surveys may already haveobserved off-centered TDEs (Lin et al. 2018; Margutti et al.2019).To summarize, we find that, for some galaxies at least,the TDE rate at z > could already be similar to the oneat z = . We also confirm that the TDE rate is globallyenhanced by about 1 order of magnitude during 100 Myr L X (erg/s)10 ( / y r ) all masses M < 4 × 10 M M > 1 × 10 M 10 ( , L X ) dd L X Figure 6.
Fraction of time spent at a a given TDE rate ( Γ ) andX-ray luminosity ( L X ). We show the mean TDE rate at fixed L X for all MBH masses (red), for light MBHs (solid black line) andfor more massive MBHs (dashed black line). For all MBH masses,there is a decrease in the TDE rate in AGNs. This is due to thatAGNs are usually powered by more massive MBHs, with a lowerTDE rate. At fixed mass, the TDE rate is independent of theX-ray luminosity. around mergers, but not necessarily for the central MBHof the main galaxy. MBHs brought by successive mergerscould see their TDE rate larger than the one of the cen-tral MBH, and actually dominate the total TDE rate of thegalaxy, resulting in fairly frequent ( ∼ of the time in oursimulation) off-centered TDEs. As AGNs and TDEs share the properties of having strongvariability and being quite luminous, it is challenging to de-tect TDEs in AGNs using standard methods and, in general,AGNs are excluded from searches of TDEs ( e.g. van Velzenet al. 2020). For these reasons, few candidates of TDEs inAGNs have been suggested ( e.g.
Blanchard et al. 2017), andit is currently difficult to constrain the TDE rate in AGNsfrom observations. Nonetheless, several groups suggest thatup to 10% of AGNs are powered by TDEs (Milosavljevi´cet al. 2006; Merloni et al. 2012). With our simulation, wecan directly test what is the TDE rate when the galaxy hasan AGN.First, we have to define when the main has an AGN.We follow Brightman & Nandra (2011) (see § −
10 keV band of the central MBH, L X , is larger than erg s − . To this purpose, we use thefollowing bolometric correction (Hopkins et al. 2007; Shenet al. 2020): L X = L bol k (32) L bol = (cid:15) r − (cid:15) r (cid:219) M • , gas c (33) k = . (cid:18) L bol L (cid:12) (cid:19) . + . (cid:18) L bol L (cid:12) (cid:19) − . . (34)We exclude here the stellar accretion when computing L bol .The reason is that including stellar accretion would result MNRAS , 1–13 (2019)
DEs under the grid in an X-ray background: the central MBH is constantly ac-creting stars at about − M (cid:12) yr − , and taking into accountstellar accretion would result in constant minimum X-rayluminosity of ∼ erg s − . This artifact comes from ourpoor ( ∼
50 kyr ) temporal resolution: in reality, TDEs occurson ∼ yr timescale with much brighter luminosity (Auchettlet al. 2017), therefore do not produce this unphysical X-raybackground. In other words, because the number of TDEsduring one timestep is small ( Γ∆ t (cid:46) ), if we were to observethe galaxy during one timestep, the fraction of time duringwhich the luminosity would be the one of a TDE would bevery small ( ∼ Γ × ∼ − for a typical duration of 1yr), and at much brighter L X . We stress here that we donot pretend to capture the details of the luminosity curveto differentiate between AGNs and TDEs: both our spatialand temporal resolution are far too poor. Our goal here isto know what would be the typical TDE rate in an AGN.Note that we also exclude all the wandering MBHs of themain galaxy, which could also produce X-rays. We did sobecause none of them has an accretion rate similar to thatof the central MBH: the second most massive MBH is only × M (cid:12) (it is the major blue MBH from § Γ and L X at all times, we can compute the jointdistribution P( Γ , L X ) , such that P d Γ d L X corresponds to thefraction of time spent a L X and Γ , as: P d Γ d L X = (cid:205) ∆ t i τ • , (35)where i corresponds to timesteps during which the X-rayluminosity and TDE rate are respectively in [ L X , L X + d L X ] and [ Γ , Γ + d Γ ] ; ∆ t i is the duration of these timestep and τ • ∼ .
76 Gyr is the time during which the MBH is followedin the simulation.We show P d Γ d L X in Fig. 6. We find a large scatter,suggesting no clear relations between TDE rate and X-rayluminosity. We compute the mean TDE rate at fixed L X (redline): ˜ Γ ( L X ) = (cid:18)∫ Γ P Γ d Γ (cid:19) (cid:30) (cid:18)∫ Γ P d Γ (cid:19) . (36)On average, the TDE rate increases with L X until erg s − where it plateaus at few − yr − . When L X reaches erg s − and the MBH is classified as an AGN (Brightman& Nandra 2011), the TDE rate starts decreasing, suggest-ing that the TDE rate is lower in AGNs. However, we recallthat the TDE rate is lower for more massive MBHs (Wang &Merritt 2004), and that more massive MBHs can shine more(assuming their luminosity is a fraction of the Eddington lu-minosity). Therefore it could be that this lower TDE ratein AGNs is simply due to that MBHs in AGNs are usuallymore massive.To test this, we split the simulation in two sub-samples:when the MBH is less massive than × M (cid:12) ( t < .
69 Gyr ; τ • ∼ .
44 Gyr ), and when it is more massive than M (cid:12) ( t > .
79 Gyr ; τ • ∼ .
22 Gyr ) . We then recompute ˜ Γ for thesetwo sub-samples (black lines). Regardless of the X-ray lumi-nosity, the TDE rate is larger for lighter MBHs, in agree- The third part, when the mass of the MBH is in between × M (cid:12) and M (cid:12) is excluded to avoid spurious results due toarbitrary transition. ment with Wang & Merritt (2004). Regarding the enhance-ment, or not, of the rate in AGNs, we find that, as long as L X > erg s − , the TDE rate is fairly constant at all L X ,confirming that the lower TDE rate in AGNs is due to moremassive MBHs.To summarize, our simulation suggests that, at fixedMBH mass, there is no enhancement of the TDE rate inAGNs. However, in general, the TDE rate should be lower inAGNs simply because AGNs are powered by massive MBHs,for which the TDE rate is lower. We have developed a physically motivated subgrid model toinclude stellar accretion on MBHs and TDEs in cosmologicalsimulations, and we have performed a cosmological zoomsimulation of a × M (cid:12) galaxy at z ∼ . Our main findingsare the following:(i) Overall, TDEs and stellar accretion do not contributemuch to the growth of MBHs, in our particular case only0.2% of the final mass comes from stars. However, TDEsare particularly efficient in growing MBHs in their early life,when they are lighter than ∼ × M (cid:12) , with more than 10%of the total accreted mass coming from stars during the first300 Myr. We stress that this value could be underestimatedas the minimum MBH mass allowed in our simulation is M (cid:12) , and that the TDE rate increases with decreasingmass. All this suggests that accretion following TDEs is apromissing channel to rapidly grow light MBHs.(ii) Stellar accretion is much smoother than gas accre-tion, this results naturally from the stellar density beingtemporally and spatially smoother than the gas density. Atany time, the gas accretion rate can be orders of magnitudelower or higher than the stellar accretion rate.(iii) When a galaxy merger occurs, the global TDE ratein a galaxy can be enhanced by up to 1 order of magnitudeduring 100 Myr. This enhancement occurs on the centralMBH of the satellite galaxy and it is caused by a nuclearstarburst or a MBH entering the main galaxy with a densestellar cusp (hence with a high TDE rate).(iv) As galaxy mergers bring many MBHs which may takea long time to sink toward the center of the main galaxy,the amount of off-centre TDEs could be fairly high. In oursimulation, the TDE rate of the main galaxy is dominatedby off-centre TDEs during 20% of the time.(v) Some galaxies with mass comparable to that of theMilky Way today could already have a well established TDErate of − − − yr − , comparable with local estimates, at z > .(vi) At fixed MBH mass, the TDE rate is independent ofthe X-ray luminosity of the central MBH, and no enhance-ment is expected in AGNs. However, since luminous AGNare powered by MBHs with mass > M (cid:12) and the TDErate decreases as M • increases, for a population of AGNsthe TDE rate is expected to be < − yr − .This is the first study of TDEs and their evolution overcosmic time using cosmological hydrodynamic simulations.While only one galaxy has been studied in this analysis,we are planning to run a cosmological volume in order to in-crease the statistical validity of our investigation and explore MNRAS , 1–13 (2019) H. Pfister et al. how stellar accretion and TDEs depend on the environmentand properties of their galaxies.
ACKNOWLEDGMENTS
HP is indebted to the Danish National Research Foundation(DNRF132) and the Hong Kong government (GRF grantHKU27305119) for support. KAA and ERR are supportedby the Danish National Research Foundation (DNRF132).Parts of this research were supported by the Australian Re-search Council Centre of Excellence for All Sky Astrophysicsin 3 Dimensions (ASTRO 3D), through project numberCE170100013. MT is supported by Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) underGermany’s Excellence Strategy EXC-2181/1 - 390900948(the Heidelberg STRUCTURES Cluster of Excellence). Theauthors thank the Yukawa Institute for Theoretical Physicsat Kyoto University. Discussions during the YITP work-shop YITP-T-19-07 on International Molecule-type Work-shop ”Tidal Disruption Events: General Relativistic Tran-sients” were useful to complete this work. This work wasmade possible with the access to the HPC resources ofCINES under allocations DARK n ° A0060406955 made byGENCI. This work has made use of the Horizon Clusterhosted by Institut d’Astrophysique de Paris; the authorsthank St´ephane Rouberol for running smoothly this cluster.
REFERENCES
Alexander T., Bar-Or B., 2017, Nature Astronomy, 1, 0147Allen R. J., et al., 2017, ApJ, 834, L11Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376Auchettl K., Guillochon J., Ramirez-Ruiz E., 2017, ApJ, 838, 149Auchettl K., Ramirez-Ruiz E., Guillochon J., 2018, ApJ, 852, 37Ba˜nados E., et al., 2018, Nature, 553, 473Bahcall J. N., Wolf R. A., 1976, ApJ, 209, 214Baumgardt H., Makino J., Ebisuzaki T., 2004a, ApJ, 613, 1133Baumgardt H., Makino J., Ebisuzaki T., 2004b, ApJ, 613, 1143Binney J., Tremaine S., 1987, Galactic Dynamics, first edn.Princeton Series in Astrophysics, Princeton University PressBlagorodnova N., et al., 2017, ApJ, 844, 46Blanchard P. K., et al., 2017, ApJ, 843, 106Bondi H., 1952, MNRAS, 112, 195Brightman M., Nandra K., 2011, MNRAS, 413, 1206Brockamp M., Baumgardt H., Kroupa P., 2011, MNRAS, 418,1308Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L.,Governato F., 2015, MNRAS, 447Chandrasekhar S., 1943, ApJ, 97, 255Dai L., McKinney J. C., Roth N., Ramirez-Ruiz E., Miller M. C.,2018, ApJ, 859, L20Donley J. L., Brandt W. N., Eracleous M., Boller T., 2002, TheAstronomical Journal, 124, 1308Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420,2662Dubois Y., et al., 2014, MNRAS, 444, 1453Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., TeyssierR., 2015, MNRAS, 452, 1502Eddington A. S., 1916, MNRAS, 77, 16Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406,2267Federrath C., Klessen R. S., 2012, ApJ, 761, 156Ferrarese L., Merritt D., 2002, ArXiv Astrophysics e-prints, French K. D., Arcavi I., Zabludoff A., 2016, ApJ, 818, L21French K. D., Wevers T., Law-Smith J., Graur O., Zabludoff A. I.,2020a, arXiv e-prints, p. arXiv:2003.02863French K. D., Arcavi I., Zabludoff A. I., Stone N., Hiramatsu D.,van Velzen S., McCully C., Jiang N., 2020b, ApJ, 891, 93Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MN-RAS, 448, 3248Gezari S., et al., 2008, ApJ, 676, 944Haardt F., Madau P., 1996, ApJ, 461, 20Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935Hahn O., Abel T., 2013, MUSIC: MUlti-Scale Initial Conditions,Astrophysics Source Code Library (ascl:1311.011)Holoien T. W. S., et al., 2016, MNRAS, 455, 2918Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731Kennicutt Jr. R. C., 1998, ApJ, 498, 541Kimm T., Cen R., 2014, ApJ, 788, 121Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., SlyzA., 2017, MNRAS, 466, 4826Kochanek C. S., 2016, MNRAS, 458, 127Kool E. C., et al., 2020, arXiv e-prints, p. arXiv:2006.01518Kormendy J., Ho L. C., 2013, ARAA, 51, 511Kroupa P., 2001, MNRAS, 322, 231Lacy J. H., Townes C. H., Hollenbach D. J., 1982, ApJ, 262, 120Lauer T. R., et al., 2007, ApJ, 664, 226Li S., Liu F. K., Berczik P., Spurzem R., 2017, ApJ, 834, 195Lightman A. P., Shapiro S. L., 1977, ApJ, 211, 244Lin D., et al., 2018, Nature Astronomy, 2, 656Liu X., Dittmann A., Shen Y., Jiang L., 2018, ApJ, 859, 8Lupi A., Volonteri M., Decarli R., Bovino S., Silk J., Bergeron J.,2019, MNRAS, 488, 4004Madau P., Dickinson M., 2014, ARAA, 52, 415Magorrian J., Tremaine S., 1999, MNRAS, 309, 447Margutti R., et al., 2019, ApJ, 872, 18Merloni A., et al., 2012, arXiv e-prints, p. arXiv:1209.3114Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei.Princeton University PressMilosavljevi´c M., Merritt D., Ho L. C., 2006, ApJ, 652, 120Ostriker E. C., 1999, ApJ, 513, 252Pechetti R., Seth A., Neumayer N., Georgiev I., Kacharov N., denBrok M., 2019, arXiv e-prints, p. arXiv:1911.09686Pfister H., Lupi A., Capelo P. R., Volonteri M., Bellovary J. M.,Dotti M., 2017, MNRAS, 471, 3646Pfister H., Bar-Or B., Volonteri M., Dubois Y., Capelo P. R.,2019a, MNRAS, p. L87Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019b,MNRAS, 486, 101Pfister H., Volonteri M., Lixin Dai J., Colpi M., 2020, arXiv e-prints, p. arXiv:2003.08133Planck Collaboration et al., 2016, AAP, 594, A13Prieto J., Escala A., Volonteri M., Dubois Y., 2017, ApJ, 836, 216Rasera Y., Teyssier R., 2006, AAP, 445, 1Rees M. J., 1988, Nature, 333, 523Rosen A., Bregman J. N., 1995, ApJ, 440, 634Sakurai Y., Yoshida N., Fujii M. S., 2018, preprint,( arXiv:1810.01985 )S´anchez-Janssen R., et al., 2019, ApJ, 878, 18Schmidt M., 1959, ApJ, 129, 243Sch¨odel R., Gallego-Cano E., Dong H., Nogueras-Lara F.,Gallego-Calvente A. T., Amaro-Seoane P., Baumgardt H.,2018, AAP, 609, A27Shen X., Hopkins P. F., Faucher-Gigu`ere C.-A., Alexander D. M.,Richards G. T., Ross N. P., Hickox R. C., 2020, arXiv e-prints,p. arXiv:2001.02696S (cid:44) adowski A., Lasota J.-P., Abramowicz M. A., Narayan R., 2016,MNRAS, 456, 3915Spitzer Jr. L., Harm R., 1958, ApJ, 127, 544Stone N. C., Metzger B. D., 2016, MNRAS, 455, 859Stone N. C., van Velzen S., 2016, ApJ, 825, L14MNRAS , 1–13 (2019)
DEs under the grid Sutherland R. S., Dopita M. A., 1993, ApJ, 88, 253Syer D., Ulmer A., 1999, MNRAS, 306, 35Tadhunter C., Spence R., Rose M., Mullaney J., Crowther P.,2017, Nature Astronomy, 1, 0061Tenneti A., Di Matteo T., Croft R., Garcia T., Feng Y., 2017,preprint, ( arXiv:1708.03373 )Teyssier R., 2002, AAP, 385, 337Toro E. F., 1997, Riemann Solvers and Numerical Methods forFluid Dynamics. Springer-Verlag Berlin HeidelbergTrebitsch M., Volonteri M., Dubois Y., Madau P., 2018, MNRAS,478, 5607Trebitsch M., Volonteri M., Dubois Y., 2019, MNRAS, 487, 819Tremmel M., Governato F., Volonteri M., Quinn T. R., 2015,MNRAS, 451, 1868Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009,AAP, 506, 647Van Wassenhove S., Capelo P. R., Volonteri M., Dotti M.,Bellovary J. M., Mayer L., Governato F., 2014, MNRAS, 439Vasiliev E., 2017, ApJ, 848, 10Vasiliev E., 2019, MNRAS, 482, 1525Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V.,Hernquist L., 2013, MNRAS, 436, 3031Volonteri M., et al., 2020, arXiv e-prints, p. arXiv:2005.04902Wang J., Merritt D., 2004, ApJ, 600, 149Zhong S., Berczik P., Spurzem R., 2014, ApJ, 792, 137van Velzen S., 2018, ApJ, 852, 72van Velzen S., Farrar G. R., 2014, ApJ, 792, 53van Velzen S., et al., 2020, arXiv e-prints, p. arXiv:2001.01409
APPENDIX A: AN ESTIMATE OF THECRITICAL RADIUS
In general, there exists no simple solution to Eq. (2). Thisstill holds when the density profile is very simple such as apower-law ( Eq. (7)) for which Eq. (2) reduces to Eq. (9).However, in some situations ( γ =
0; 1; 2 ), Eq. (9) is a poly-nomial with simple solutions (which we do not report here)and r c can be expressed.In Fig. A1, we show the in the top panel the exact so-lution solving the polynomial (thick lines) and our approxi-mate solution given by Eq. (18) (thin lines), and in the bot-tom panel relative difference between solutions. For γ span-ning between 0 and 2, i.e. almost all the value allowed inour subgrid model, and for ρ / ρ u spanning 6 orders of mag-nitude, the relative diffence peaks at 30%, which we consideras “reasonable” given the assumptions of the model. APPENDIX B: EFFECTIVE RADIUS OFGALAXIES
Contrary to halos, for which the virial radius can be definedto obtain the “size” of the structure, there are no clear def-inition for the size of a galaxy. In this Appendix we definethe effective radius R eff which we use for the “size” of thegalaxy.Once gravitationnally bound structures have been de-tected with AdaptaHOP , we compute the pseudo–inertiatensor: (cid:101) I ij = (cid:213) k m k x i , k x j , k , (B1)where the sum is made on stellar particles k belonging to / u r c / r u = 0= 1= 210 / u | r c , r e a l r c , r c p a p e r |/ r c , r e a l Figure A1. Top:
Exact solution of Eq. (9) (thick lines) and ap-proximate solution from Eq. (18).
Bottom:
Relative differencebetween r c , real , the real solution of Eq. (2), and r c , fit , the approx-imate value given by Eq. (18). the galaxy, with masses m k and positions ( x , k , x , k , x , k ) = ( x k , y k , z k ) from the center of the galaxy.From (cid:101) I we can obtain the principal ellipsoid of thegalaxy. The eigenvectors are the principal directions, andthe eigenvalues ( I , I , I ) are related to the principal axis ( a , a , a ) by: I i = Ma i , (B2)where M is the total mass of the galaxy, and the 1/5 factoris added so that the equation is correct for a homogeneousellipsoid.Once the principal ellipsoid is known, we compute themass in concentric ellipsoid and find the one which contains90% of the total mass of the galaxy. The principal axis ofthis ellipsoid are ( a , eff , a , eff , a , eff ) = ( α eff a , α eff a , α eff a ) , α eff > , so that the effective radius is given by: R eff = ( a , eff a , eff a , eff ) / . (B3) This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000