NIHAO XXII: Introducing black hole formation, accretion and feedback into the NIHAO simulation suite

TL;DR

This paper integrates black hole physics into the NIHAO galaxy simulation suite, enabling detailed study of black hole effects on galaxy evolution, especially in high-mass elliptical galaxies, and calibrates the models against observed galaxy relations.

Contribution

It introduces algorithms for black hole formation, accretion, and feedback into NIHAO, expanding the simulation suite to include a large, diverse sample of galaxy masses.

Findings

Successfully calibrates black hole models against observed stellar and black hole mass relations.

Demonstrates that black hole feedback quenches star formation in massive galaxies.

Findings on the scatter of galaxy relations decreasing over time, aligning with observations.

Abstract

We introduce algorithms for black hole physics, i.e., black hole formation, accretion and feedback, into the NIHAO (Numerical Investigation of a Hundred Astrophysical Objects) project of galaxy simulations. This enables us to study high mass, elliptical galaxies, where feedback from the central black hole is generally thought to have a significant effect on their evolution. We furthermore extend the NIHAO suite by 45 simulations that encompass $z=0$ halo masses from $1 \times 10^{12}$ to $4 \times 10^{13}\,\mathrm{M}_{\odot}$, and resimulate five galaxies from the original NIHAO sample with black hole physics, which have $z=0$ halo masses from $8 \times 10^{11}$ to $3 \times 10^{12}\,\mathrm{M}_{\odot}$. Now NIHAO contains 144 different galaxies and thus has the largest sample of zoom-in simulations of galaxies, spanning $z=0$ halo masses from $9 \times 10^{8}$ to $4 \times 10^{13}\,\mathrm{M}_{\odot}$. In this paper we focus on testing the algorithms and calibrating their free parameters against the stellar mass versus halo mass relation and the black hole mass versus stellar mass relation. We also investigate the scatter of these relations, which we find is a decreasing function with time and thus in agreement with observations. For our fiducial choice of parameters we successfully quench star formation in objects above a $z=0$ halo mass of $10^{12}\,\mathrm{M}_{\odot}$, thus transforming them into red and dead galaxies.

Tables

Table 1. Table 1: Galaxies and their properties: particle number N 𝑁 N , dark matter particle number N DM subscript 𝑁 DM N_{\mathrm{DM}} , star particle number N ⋆ subscript 𝑁 ⋆ N_{\star} , halo mass M 200 subscript 𝑀 200 M_{200} , stellar mass M ⋆ subscript 𝑀 ⋆ M_{\star} and black hole mass M BH subscript 𝑀 BH M_{\mathrm{BH}} . First section: galaxies from the original NIHAO sample from the 60 Mpc box, second section: new galaxies from the 60 Mpc box, third section: new galaxies from the 90 Mpc box. Galaxies without black hole mass are re-simulations without black holes.

galaxy	$N$	$N_{\mathrm{DM}}$	$N_{\star }$	$\log \left(M_{200}\left[M_{☉}\right]\right)$	$\log \left(M_{\star }\left[M_{☉}\right]\right)$	$\log \left(M_{\mathrm{BH}}\left[M_{☉}\right]\right)$
g7.55e11	932961	431968	327418	11.92	10.31	7.41
g8.26e11	1070666	496667	424753	11.97	10.42	7.80
g1.12e12	1062175	516394	430205	11.98	10.42	7.78
g1.92e12	1960845	1035887	736102	12.28	10.66	8.25
g2.79e12	3383907	1663073	1395558	12.49	10.93	8.39
g7.55e11	1185149	455930	483752	11.95	10.49	-
g8.26e11	1513265	518236	739749	12.01	10.67	-
g1.12e12	1977112	564785	1222359	12.05	10.90	-
g1.92e12	4018048	1200667	2467821	12.37	11.20	-
g2.79e12	5598386	1800095	3099962	12.55	11.29	-
g1.26e12	1253780	563159	462669	12.03	10.46	7.55
g1.27e12	1150753	471966	605247	11.94	10.57	7.84
g1.55e12	1469159	600775	679269	12.06	10.61	7.86
g1.62e12	1063633	533023	327877	12.00	10.31	7.61
g2.37e12	2159202	1098881	881920	12.30	10.73	8.30
g2.71e12	2164500	1195318	596766	12.35	10.55	8.19
g3.74e12	3397933	1715201	1182008	12.51	10.83	8.47
g4.41e12	463537	250487	204647	12.55	10.83	8.84
g4.55e12	503192	270807	164411	12.60	10.88	8.69
g4.81e12	594745	306880	211336	12.66	11.01	8.70
g4.84e12	523621	267346	185147	12.60	10.93	8.66
g5.22e12	592556	313138	187130	12.67	10.96	8.58
g5.41e12	644421	290846	352249	12.62	11.23	9.14
g5.53e12	547812	306764	183274	12.65	10.83	8.69
g6.53e12	591850	276231	266355	12.61	11.10	8.73
g6.57e12	763439	356974	303556	12.73	11.15	8.75
g6.86e12	838063	430471	303160	12.80	11.03	8.90
g6.70e12	1290678	431794	693967	12.83	11.54	8.83
g7.50e12	648889	387253	259728	12.74	11.08	9.32
g7.55e12	860302	433196	262410	12.81	10.95	8.79
g7.71e12	823620	415601	312032	12.79	11.17	8.94
g7.92e12	1323285	446337	731566	12.84	11.55	9.06
g8.08e12	1003594	486702	364796	12.86	11.22	9.00
g8.45e12	779820	368069	308291	12.74	11.17	8.80
g1.05e13	1146928	573974	489758	12.92	11.34	9.18
g1.14e13	1354041	646068	532474	12.98	11.41	9.16
g1.17e13	2085075	889400	782176	13.14	11.32	8.92
g1.25e13	1217182	673068	367429	13.00	11.18	9.18
g1.33e13	1918776	809672	751999	13.09	11.58	9.24
g1.44e13	2836084	876310	1474963	13.14	11.86	8.94
g1.54e13	1456656	722348	500173	13.03	11.36	9.17
g1.57e13	889955	436443	310259	12.81	11.16	8.85
g1.63e13	2089314	953300	684168	13.16	11.43	9.25
g1.87e13	2278215	1101553	636663	13.23	11.39	9.29
g2.02e13	1947099	1056092	840843	13.18	11.59	9.59
g2.07e13	3177787	1497578	1020594	13.36	11.47	9.45
g2.10e13	2972176	1249734	1072224	13.29	11.68	9.33
g2.11e13	2434943	1250655	955862	13.26	11.63	9.45
g2.20e13	1437878	847511	578127	13.08	11.41	9.46
g2.37e13	2872994	1497621	1030754	13.34	11.59	9.54
g2.58e13	3951084	1782390	1347000	13.44	11.66	9.47
g3.26e13	4333704	2100752	1118595	13.51	11.45	9.52
g3.42e12	2620934	1068134	1313357	12.30	10.91	8.12
g3.78e13	3462973	1607283	1074439	13.39	11.51	9.45
g3.89e13	5426667	2281160	1751607	13.55	11.81	9.70
g5.41e12	1265524	337582	778618	12.73	11.59	-
g6.86e12	1656812	453120	975692	12.86	11.69	-
g7.50e12	1478432	443431	834490	12.85	11.62	-
g7.92e12	1546499	455160	925496	12.85	11.65	-
g1.05e13	2286298	658752	1301680	13.02	11.75	-
g1.44e13	2738604	878724	1374316	13.14	11.83	-

Equations

$${\dot{M}}_{\mathrm{BHL}}=\frac{4\text{uppi}\alpha G^2{M}_{\mathrm{BH}}^2\rho }{(c_{\mathrm{s}}^2+v^2{)}^{3/2}},$$

$${\dot{M}}_{\mathrm{Edd}}=\frac{M_{\mathrm{BH}}}{{ϵ}_{\mathrm{r}}{\tau }_{\mathrm{S}}}$$

$${\dot{M}}_{\mathrm{BH}}=\min ({\dot{M}}_{\mathrm{BHL}},{\dot{M}}_{\mathrm{Edd}}).$$

$$L={ϵ}_{\mathrm{r}}{\dot{M}}_{\mathrm{BH}}{c}^2$$

$$\dot{E}={ϵ}_{\mathrm{f}}{ϵ}_{\mathrm{r}}{\dot{M}}_{\mathrm{BH}}{c}^2$$

$${\Delta }_{\mathrm{IMF}}=\min \left(0.28+\left(0.14m-11\right),0.3\right)$$

Full text

††pubyear: 2019††pagerange: NIHAO XXII: Introducing black hole formation, accretion and feedback into the NIHAO simulation suite–1

1 Introduction

It is now well established that black holes exist in the centres of almost all galaxies (e.g., Kormendy & Richstone, 1995; Magorrian et al., 1998), and that they are connected to their hosts. The general idea outlined by Di Matteo et al. (2005) is that feedback provided by the central black hole heats the galaxy’s gas and subsequently quenches star formation and further black hole accretion, transforming it into a ‘red and dead’  elliptical galaxy. In this framework this co-evolution leads the galaxy, i.e. its velocity dispersion (e.g., Ferrarese & Merritt, 2000; Gebhardt et al., 2000) or its bulge mass (e.g., Kormendy & Richstone, 1995; Häring & Rix, 2004), being related to the black hole mass. However, some works (e.g., Jahnke & Macciò, 2011) argue that correlations between the galaxy and its central black hole do not imply a physical connection of these two.

Numerical simulations have been proven very successful in investigating the formation and evolution of galaxies, leading to a number of projects that have been developed over the last years. E.g., Illustris (Vogelsberger et al., 2014), IllustrisTNG (Pillepich et al., 2018), Magneticum Pathfinder (Dolag et al., 2016), MassiveBlack-II (Khandai et al., 2015) and EAGLE (Schaye et al., 2015) are hydrodynamic simulations of cosmological volumes, and FIRE (Hopkins et al., 2014) and Auriga (Grand et al., 2016) consist of zoom-in simulations of individual galaxies.

The NIHAO project (Wang et al., 2015) is a suite of hydrodynamical cosmological zoom-in simulations of galaxies and is unique in combining (i) a high resolution of $\sim 10^{6}$ particles per halo, (ii) a large range of halo masses from dwarf to Milky Way masses ($\sim 5\times 10^{9}$ to $\sim 2\times 10^{12}\,\mathrm{M}_{\sun}$), and (iii) a large sample size of $\sim$ 100 galaxies. NIHAO successfully reproduces the stellar mass versus halo mass relation, one of the most fundamental constraints in galaxy formation and evolution, from redshift $z=0$ up to redshift $z=4$, and also matches the star formation rate versus stellar mass relation. However, up to this date the NIHAO suite did not include high mass, elliptical galaxies, because the stellar feedback implemented into NIHAO is insufficient to regulate star formation in elliptical galaxies (Dutton et al., 2015). It is now generally accepted (Croton et al., 2006) that black hole feedback is needed at higher masses to reproduce the sharp decline in the stellar mass function and to create red and dead galaxies.

The aim of this paper is to introduce algorithms for black hole formation, accretion and feedback into the NIHAO project, and to extend the NIHAO suite to include high mass, elliptical galaxies. We focus on the testing of the algorithms and the calibration of the free parameters on the observed stellar mass versus halo mass and black hole mass versus stellar mass relations, and also investigate how the scatter of these relations evolve with time and compare with observations. In section 2 we review the properties of the NIHAO project. We introduce the algorithms used for black hole formation, accretion and feedback in section 3 and the initial conditions in section 4. In section 5 we present our results and in section 6 we study the effect of the algorithms’ free parameters on the simulations. In section 7 we summarize our findings.

2 The NIHAO project

NIHAO uses an updated version (Keller et al., 2014) of the TreeSPH code Gasoline2 (Wadsley et al., 2004, 2017). The simulations use a flat LCDM cosmology with parameters from the Planck Collaboration et al. (2014).

The backbone of the NIHAO project consists of cosmological dark-matter-only simulations with box sizes of 60 and 20 Mpc/h from Dutton & Macciò (2014), and a new box of 15 Mpc/h with $400^{3}$ particles. These are evolved until $z=0$, then haloes from these boxes are selected and resimulated individually with a higher resolution and with gas particles.

All galaxies have the same relative resolution across the whole mass range, i.e. $\sim 10^{6}$ dark matter particles inside the virial radius at $z=0$. The particle masses and force softening lengths are chosen to resolve the mass profile at $\leq$ 1 per cent of the virial radius. The initial ratio of dark and gas particle mass equals the cosmological mass ratio of dark matter and baryons of $\Omega_{\mathrm{DM}}/\Omega_{\mathrm{b}}=5.48$. The softening length of gas and star particles is $(\Omega_{\mathrm{DM}}/\Omega_{\mathrm{b}})^{1/2}=2.34$ times smaller than the dark matter particle softening length. The free parameters of the stellar and supernova feedback model have been chosen to match the $M_{\star}$-$M_{{200}}$ relation for one Milky Way-like galaxy at $z=0$.

Cooling is provided via hydrogen, helium, and various metal-lines in a uniform ultraviolet ionizing background (Shen et al., 2010), this includes photoionization, UV background heating (Haardt & Madau, 2012), and Compton cooling.

Stars are formed from gas particles that pass a density and temperature threshold ($T<15000\,\mathrm{K}$, $n>10.3\,\mathrm{cm}^{-3}$) with a rate of $\dot{M_{\star}}=c_{\star}M_{\mathrm{gas}}t_{\mathrm{dyn}}^{-1}$, where $t_{\mathrm{dyn}}=(4\uppi G\rho)^{-1/2}$ is the gas particle’s dynamical time, $\rho$ its density, $M_{\mathrm{gas}}$ its mass and $c_{\star}=0.1$ the star formation efficiency.

Supernova feedback is modeled with the blastwave formalism of Stinson et al. (2006). Stars of mass $8<M_{\star}/M_{\sun}<40$ eject metals and energy to surrounding gas particles 4 Myr after their formation. For these gas particles cooling is delayed for $\sim 30\,\mathrm{Myr}$. Before they produce a supernova massive stars provide ‘early stellar feedback’  (Stinson et al., 2013), i.e. 13 per cent of the total stellar flux of $2\times 10^{50}\,\mathrm{erg}\,\mathrm{M}_{\sun}^{-1}$ is injected into the surrounding gas as thermal energy. No cooling delay is applied in this case. We refer to Wang et al. (2015) for more details of the NIHAO project.

NIHAO has been very successful in reproducing galaxy properties for halo masses of $M_{200}\leq 2\times 10^{12}\,\mathrm{M}_{\sun}$, e.g., the stellar mass versus halo mass relation (Wang et al., 2015), the galaxy velocity function (Macciò et al., 2016), the Tully-Fisher relation (Dutton et al., 2017) and the rotation curves of dwarf galaxies (Santos-Santos et al., 2018).

3 Computational methods for black hole physics

Black holes are modeled as sink particles (Bate et al., 1995) that only interact with their environment via gravitational forces and that can accrete matter from neighboring gas particles. For black hole accretion and feedback we choose the models introduced by Springel et al. (2005), because these are the most widely used and thus tested, and are known to be able to yield the correct relation between black hole mass and the stellar component (Di Matteo et al., 2005).

3.1 Black hole formation

We follow a common approach for modeling black hole formation: when a central halo111We do not seed subhaloes. exceeds a threshold mass $M_{\mathrm{h,t}}$ we convert the gas particle (or a part thereof) with the lowest gravitational potential into a black hole with seed mass $M_{\mathrm{BH,s}}$.222E.g., Sijacki et al. (2007) use $M_{\mathrm{BH,s}}=1\times 10^{5}\,\mathrm{M}_{\sun}$ and $M_{\mathrm{h,t}}=5\times 10^{10}$ M*☉, Di Matteo et al. (2008) and Schaye et al. (2015) use $M_{\mathrm{BH,s}}=1\times 10^{5}\,\mathrm{M}_{\sun}$ and $M_{\mathrm{h,t}}=1\times 10^{10}$ M☉, Schaye et al. (2010) use $M_{\mathrm{BH,s}}=9\times 10^{4}\,\mathrm{M}_{\sun}$ and $M_{\mathrm{h,t}}=4\times 10^{10}\,\mathrm{M}_{\sun}$ and Sijacki et al. (2015) use $M_{\mathrm{BH,s}}=1.4\times 10^{5}\,\mathrm{M}_{\sun}$ and $M_{\mathrm{h,t}}=7.1\times 10^{10}\,\mathrm{M}_{\sun}$. Haloes and their masses are found with the AMIGA Halo Finder (AHF, Gill et al., 2004; Knollmann & Knebe, 2009), see section 4 for more details. We follow Sijacki et al. (2007) and use the values $M_{\mathrm{BH,s}}=1\times 10^{5}$ M☉* and $M_{\mathrm{h,t}}=5\times 10^{10}\,\mathrm{M}_{\sun}$. An alternative model for black hole formation is used in the Romulus simulations (Tremmel et al., 2017), where black hole formation occurs when the gas meets specific thresholds for metallicity, density and temperature.

3.2 Black hole relocation

In galaxy-scale simulations two effects can lead the position of the black hole not to coincide with the halo centre: (i) The small ratio of black hole mass and gas/star/dark particle mass leads to stochastic motions of the black hole particle caused by the momentum of the accreted gas and gravitational interactions with nearby particles. (ii) In the case of two haloes merging dynamical friction would lead the black hole to sink to the centre of the newly forming halo. However, in galaxy-scale simulations dynamical friction is generally underestimated on small scales due to insufficient resolution.

To compensate for these effects a number of models have been developed in the last years: Debuhr et al. (2011, 2012) and Anglés-Alcázar et al. (2017) assign the black hole a ‘dynamical mass’  or ‘tracer mass’, which is several orders of magnitude larger than the actual black hole mass, to prevent the black hole from moving too far away from the halo centre. Johansson et al. (2009) relocate the black hole to the gas particle within the SPH smoothing length of the BH that has the lowest gravitational potential. Booth & Schaye (2009) also relocate the black hole to the neighboring gas particle with the lowest gravitational potential, but only if the relative velocity of black hole and gravitationally most bound gas particle is smaller than 25 per cent of the local speed of sound and if the black hole mass is smaller than ten times the initial gas particle mass. The IllustrisTNG simulations (Sijacki et al., 2015) force the black hole particle to be located at the potential minimum of its host halo. The EAGLE simulations (Schaye et al., 2015) force black holes with a mass smaller than 100 times the gas particle mass to migrate towards the minimum of the gravitational potential of its host halo. The Magneticum Pathfinder simulations (Hirschmann et al., 2014) and the the Romulus simulations (Tremmel et al., 2017; Tremmel et al., 2015) include a subgrid model to account for the dynamical friction force that is acting on the black hole.

However, relocating the black hole to the potential minimum of its host halo can lead it to move very large distances (the sum of the haloes virial radii, i.e., several $100\,\mathrm{kpc}$) in the case of two haloes merging. Furthermore, if these two haloes, whose black holes have just merged, disconnect again, the halo that just lost its black hole would get re-seeded. To avoid these problems every major time step333The simulation time of 13.8 Gyr is divided into 1024 major time steps of 13.4 Myr. we set the position and velocity of the black hole to the values of the dark matter particle within ten softening lengths that has the lowest gravitational potential. Using the star particle with the lowest gravitational potential might be problematic at high redshifts, when only few or no star particles are present. Placing the black hole at the gas particle with the lowest gravitational potential would artificially increase its kernel weight, thus its density and thus the black hole accretion rate. Furthermore, using a gas particle that has recently received feedback and is thus outflowing could lead to inappropriate black hole velocities and/or large jumps in the black hole position. (See Wurster & Thacker, 2013, for a discussion of a similar problem.).

3.3 Black hole merging

Two black holes are merging when their distance is smaller than the sum of their softening lengths. The resulting black hole inherits the position, velocity and acceleration of its predecessors’ barycentre, its mass is the sum of its predecessors’ masses, and its accretion rate is calculated in the next timestep.

3.4 Black hole accretion

We calculate the accretion rate of the black hole with the commonly used444E.g., in Di Matteo et al. (2005), Springel et al. (2005), Colberg & Di Matteo (2008), Di Matteo et al. (2008), Booth & Schaye (2009), Croft et al. (2009), Johansson et al. (2009), Choi et al. (2012). Bondi-Hoyle-Lyttleton parametrization (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952)

[TABLE]

where $M_{\mathrm{BH}}$ is the black hole mass and $\rho$, $c_{\mathrm{s}}$ and $v$ are the density, sound speed and velocity of the gas that surrounds the black hole. The parameter $\alpha$ was first introduced by Springel et al. (2005) to account for the limited resolution of these simulations, and is usually set to $\alpha=100$ (see, e.g., table 2 in Booth & Schaye, 2009). However, our resolution is higher than in earlier works that use the Bondi-Hoyle-Lyttleton parametrization, therefore we use an accretion parameter of $\alpha=70$ and justify this choice with a parameter study in section 6.

The black hole accretion rate is limited by the Eddington rate (Eddington, 1921)

[TABLE]

with the Salpeter time-scale $\tau_{\mathrm{S}}=4.5\times 10^{8}\,\mathrm{yr}$ (Salpeter, 1964) and the radiative efficiency $\epsilon_{\mathrm{r}}=0.1$ (Shakura & Sunyaev, 1973). The black hole accretion rate then is

[TABLE]

At each time step $\Delta t$ the black hole accretes the mass $\Delta t\,\dot{M}_{\mathrm{BH}}$ from the gas particle(s) that is (are) most gravitationally bound to the black hole. The momentum of the accreted mass is added to the black hole’s momentum. The black hole can accrete fractions of a gas particle, which might lead to gas particles with a very small mass. If a gas particle’s mass falls below 20 per cent of its initial mass it is deleted and its mass and momentum are distributed among the surrounding gas particles weighted with the SPH kernel.

To avoid too large accelerations and thus too small time steps of particles close to the black hole, as well as two body relaxation, we increase the black hole softening length as it grows. As the softening length of collisionless particles in our simulations is proportional to the square root of their mass, we multiply the black hole softening length with $(1+\Delta m/m_{\mathrm{BH}})^{1/2}$ when it increases by $\Delta m$ in mass.

Increasing the softening length of the black hole does not significantly affect its surroundings, as its gravitational force is usually much smaller than the gravitational force of the stars surrounding the centre. E.g., for the galaxy g7.92e12 that we use for our parameter study in section 6, the mass of the stars within $1\,\mathrm{kpc}$ of the black hole is at least 25 times larger than the mass of the black hole for all times.

We do not limit how far particles can be away in order to be accreted. However, usually in each timestep not more than one particle is accreted, and this one is the most bound one, and thus usually the closest to the black hole.

Alternative methods to model gas accretion onto black holes have been developed: Booth & Schaye (2009, 2010) assume that the accretion parameter $\alpha$ is not constant, but a function of the gas density, whereas in the EAGLE simulations (Schaye et al., 2015) and Rosas-Guevara et al. (2015) $\alpha$ is the ratio of Bondi time-scale and viscous time-scale. Tremmel et al. (2017) modify the Bondi-Hoyle-Lyttleton parametrization to take the angular momentum of the gas into account. Debuhr et al. (2011, 2012) calculate the black hole accretion rate based on the viscous evolution of an accretion disc that surrounds the black hole. Hopkins & Quataert (2011) and the FIRE simulations (Anglés-Alcázar et al., 2017) calculate the black hole accretion rate based on angular momentum and gravitational torques around the black hole.

3.5 Black hole feedback

Black hole accretion results in a luminosity of

[TABLE]

with speed of light $c$, and we assume that a fraction555Common values for the feedback efficiency are $\epsilon_{\mathrm{f}}=0.05$ (Di Matteo et al., 2005; Springel et al., 2005; Sijacki et al., 2007; Di Matteo et al., 2008; Johansson et al., 2009), $\epsilon_{\mathrm{f}}=0.15$ (Booth & Schaye, 2009, 2010; Schaye et al., 2015) and $\epsilon_{\mathrm{f}}=0.02$ (Tremmel et al., 2017). $\epsilon_{\mathrm{f}}=0.05$ of this luminosity is available as thermal energy for the gas that surrounds the black hole. Thus the gas receives an energy per time of

[TABLE]

that is distributed kernel weighted among the 50 nearest gas particles. To avoid too large sound speeds and thus too small time steps we limit the specific energy of a single gas particle to $(0.1c)^{2}$.

There are a few alternative models for black hole feedback. Sijacki et al. (2007) use an additional feedback mode operating at low black hole accretion rates that is modeled by injecting ‘bubbles’  into the host galaxy. Debuhr et al. (2011, 2012) and Choi et al. (2012) use ‘kinetic feedback’  that feeds momentum and mass to the gas surrounding the black hole.

4 Initial conditions, galaxy properties and parameters

Our galaxies and their properties are listed in Table 1 at the end of this paper. We take five galaxies from the original NIHAO suite that have a dark matter particle mass of $1.74\times 10^{6}\,\mathrm{M}_{\sun}$, a gas particle mass of $3.17\times 10^{5}\,\mathrm{M}_{\sun}$, a dark matter softening length of $931\,\mathrm{pc}$ and a gas particle softening length of $398\,\mathrm{pc}$. These originate from a cosmological dark-matter-only simulation with a box size of 60 Mpc and 6003 particles. From the same box we additionally take seven new galaxies with the same resolution that reach $z=0$ halo masses of $1-4\times 10^{12}\,\mathrm{M}_{\sun}$. We then take 38 new galaxies with a dark matter particle mass of $1.38\times 10^{7}\,\mathrm{M}_{\sun}$, a gas particle mass of $2.52\times 10^{6}\,\mathrm{M}_{\sun}$, a dark matter softening length of $1863\,\mathrm{pc}$ and a gas particle softening length of $782\,\mathrm{pc}$. These reach $z=0$ halo masses of $4\times 10^{12}$ to $4\times 10^{13}\,\mathrm{M}_{\sun}$, and originate from a cosmological dark-matter-only simulation with a box size of 90 Mpc and 4503 particles. A few galaxies were resimulated without black hole feedback, these have no entry for the black hole mass in Table 1.

The name of a galaxy refers to its halo mass at $z=0$ from the cosmological dark-matter-only simulations, i.e. its halo mass in the simulations presented in this paper might be slightly different. We use the AMIGA Halo Finder (AHF, Gill et al., 2004; Knollmann & Knebe, 2009) to identify haloes and their properties. AHF identifies overdensities on an adaptively refined grid as prospective halo centres. These haloes are then defined as a sphere with radius $R_{{200}}$ and mass $M_{{200}}$ such that their density equals 200 times the cosmic critical matter density. The stellar mass $M_{\star}$ of a galaxy is defined as the combined mass of all stellar particles within a radius of $r_{\mathrm{gal}}=0.2\,R_{{200}}$ (see Munshi et al., 2013, for a different approach in defining a galaxy’s stellar mass). Due to mergers some galaxies might contain more than one black hole, we define the central black hole as the one that is closest to the galaxy’s centre. All our haloes contain no intruder within 20 per cent of their virial radius and less than 20 intruders in total, where an intruder is a particle that originates from outside the zoom-in region. In total we have 50 galaxies with black hole feedback and 11 galaxies without black hole feedback in this paper.

5 Results

The magnitude of black hole feedback is determined by the black hole accretion rate, which in turn determines the black hole mass. As black hole feedback quenches star formation, one of its strongest effects is to reduce the stellar mass of a galaxy. Therefore we use the galaxy’s stellar mass and black hole mass to gauge our model, and thus use the stellar mass versus halo mass ($M_{\star}$-$M_{{200}}$) relation and the black hole mass versus stellar mass ($M_{\mathrm{BH}}$-$M_{\star}$) relation for calibration. We also analyze the star formation histories of three of our new galaxies to confirm black hole feedback having a quenching effect on them. In this section we use the reference parameters of our model, namely the accretion parameter $\alpha=70$, the feedback efficiency $\epsilon_{\mathrm{f}}=0.05$, the black hole seed mass $M_{\mathrm{BH,s}}=1\times 10^{5}\,\mathrm{M}_{\sun}$ and the halo threshold mass $M_{\mathrm{h,t}}=5\times 10^{10}\,\mathrm{M}_{\sun}$.

We compare our $M_{\star}$-$M_{{200}}$ relation to results from halo abundance matching of Moster et al. (2018), Behroozi et al. (2013) and Moster et al. (2013). These are based on the IMF of Chabrier (2003), but there is evidence for ‘heavier’  (bottom-heavy, i.e. more low mass stars) IMFs in massive galaxies (e.g., Conroy & van Dokkum, 2012; Dutton et al., 2013a; Dutton et al., 2013b). Therefore, following Dutton et al. (2013b), we add a correction

[TABLE]

for $m>9$ and $m=\log\left(M_{\star}/\mathrm{M}_{\sun}\right)$ to $\log M_{\star}$ of the $M_{\star}$-$M_{{200}}$ relations, which increases the stellar mass for high mass galaxies by a factor of $\sim$2. The Behroozi et al. (2013) and Moster et al. (2018) relations are a function of the halo’s virial mass, which we convert to $M_{200}$ by applying corrections suggested by Dutton & Macciò (2014).

Fig. 1 shows the $M_{\star}$-$M_{{200}}$ relation of our galaxies compared to the abundance matching results, our simulations with black hole feedback provide a good match. Only at very high and low halo masses our galaxies slightly deviate from the observed relations, but most of them are still within the 1-sigma scatter of the observations. Fig. 2 shows the $M_{\star}$-$M_{{200}}$ relation of our galaxies that were simulated without black holes, compared to their counterpart with black holes. Only the simulations with black holes provide a good match to the observed relations, demonstrating the importance of black hole feedback for the evolution of high mass galaxies. Fig. 3 shows the $M_{\mathrm{BH}}$-$M_{\star}$ relation of our galaxies for different redshifts compared with $z=0$ observations from Kormendy & Ho (2013) and Sani et al. (2011). At each redshift we fit a linear relation to our simulations and calculate the 1-sigma scatter. We apply a 5-sigma clipping to our data to remove the high-mass outlier that is visible in the $z=4$ subplot from the calculation of the scatter. Our results match the Kormendy & Ho (2013) relation well, only at the lower end the black hole masses of some galaxies are underestimated. We also show a relation between halo mass and black hole mass in Fig. 4.

Recent successes in reproducing the correct stellar masses and black hole masses in galaxy simulations include the EAGLE project (Schaye et al., 2015), which consists of cosmological simulations that reproduce the $M_{\star}$-$M_{200}$ relation in $\mathrm{log}\,M_{200}$ from 10.5 to 14.5 and the $M_{\mathrm{BH}}$-$M_{\star}$ relation in $\mathrm{log}\,M_{\star}$ from 8 to 12. The $M_{\star}$-$M_{200}$ relation of the Illustris project (Genel et al., 2014; Vogelsberger et al., 2014) is in good qualitative agreement with observations in $\mathrm{log}\,M_{200}$ from 10 to 14, also their $M_{\mathrm{BH}}$-$M_{\mathrm{bulge}}$ relation (in $\mathrm{log}\,M_{\mathrm{bulge}}$ from 8 to 12) is in good agreement with observations (Sijacki et al., 2015). The IllustrisTNG simulations (Weinberger et al., 2018) reproduce the $M_{\mathrm{BH}}$-$M_{\star}$ relation in $\mathrm{log}\,M_{\star}$ from 8 to 12. Also The ROMULUS cosmological simulations (Tremmel et al., 2017) match the $M_{\star}$-$M_{{\mathrm{vir}}}$ relation in $\mathrm{log}\,M_{\mathrm{vir}}$ from 10 to 13 and the $M_{\mathrm{BH}}$-$M_{\mathrm{\star}}$ relation in $\mathrm{log}\,M_{\mathrm{\star}}$ from 9 to 11.5. The Magneticum Pathfinder simulations (Hirschmann et al., 2014) match the $M_{\mathrm{BH}}$-$M_{\mathrm{stellar}}$ relation in $\mathrm{log}\,M_{\mathrm{stellar}}$ from 10.5 to 13.

Fig. 5 shows the scatter of the $M_{\star}$-$M_{{200}}$ relation versus time for our simulations and the Moster et al. (2018), Behroozi et al. (2013) and Moster et al. (2013) relations. The scatter of the Moster relations depends on the stellar mass, therefore we calculate the scatter for the range of stellar masses that are encompassed by our simulations in each redshift, its mean and standard deviation is shown in Fig. 5. We make sure that at each redshift there are at least 12 galaxies in the sample given our criteria of a minimum of $10^{4}$ particles per galaxy. The scatter of all relations is generally decreasing with time. Only at higher redshifts the scatter of our relation is rising, possibly caused by low number statistics. The scatter of our relation is always below the Moster et al. (2013) scatter, and below the Behroozi et al. (2013) scatter for z<1.5. The Moster et al. (2018) scatter is in good agreement with our simulations, with being only slightly below our scatter for z>1 and above for z<0.5.

Fig. 6 shows the scatter of the $M_{\star}$-$M_{{200}}$ relation versus halo mass for $z=0$ for the Behroozi et al. (2013), Moster et al. (2018) and Moster et al. (2013) relations, for the NIHAO galaxies with black holes presented in this paper (50 galaxies divided into 3 bins) and for 78 of the NIHAO galaxies (divided into 4 bins) without black holes presented in Wang et al. (2015) for four different redshifts. For $z=0$ our scatter is below the Moster et al. (2018) and Moster et al. (2013) scatter, except for halo masses around $10^{10}\,\mathrm{M}_{\sun}$, and below the Behroozi et al. (2013) scatter for halo masses above $5\times 10^{11}\,\mathrm{M}_{\sun}$. The scatter is generally increasing with decreasing halo mass. There is no clear trend as a function of redshift. Our results are in agreement with previous studies of the scatter by Wechsler & Tinker (2018) and Matthee et al. (2017). So far only the NIHAO simulations can show the scatter for halo masses covering five orders of magnitude, i.e., ranging from $5\times 10^{8}$ to $4\times 10^{13}\,\mathrm{M}_{\sun}$.

Fig. 7 shows the scatter of the $M_{\mathrm{BH}}$-$M_{\star}$ relation versus time for our simulations and the $z=0$ scatter of the Sani et al. (2011) and Kormendy & Ho (2013) relations. We only calculate the scatter for redshifts with at least 12 galaxies in our sample. At high redshifts the scatter is dominated by the black hole seed mass and possibly by low number statistics. The scatter then increases steeply, and is then declining for the remainder of the evolution. At redshift zero the scatter of the simulations is below the observed scatter. This is expected since we present here the intrinsic scatter in the simulations which does not account for observational biases and errors.

In Fig. 8 we show the star formation history, with and without black hole feedback, and the black hole accretion rate of three of our galaxies. Without feedback the galaxies show continuous star formation up to $z=0$, creating unrealistic and unobserved blue massive galaxies. With feedback the galaxies show an initial phase of high star formation, and are then quenched for the remainder of their evolution. The back hole accretion rate peaks shortly before the star formation rate starts declining, indicating that star formation is quenched by the black hole feedback.

All three galaxies in Fig. 8 show the black hole becoming active at about 2 to 4 gigayears. To determine the reason for this increase we look at Fig. 9, which shows the black hole accretion rate, the black hole mass and the gas density at the location of the black hole. The black hole accretion rate is given by the Bondi formula of eq. 1, and is determined mainly by the black hole mass and the gas density. In the first few time steps both quantities results in a low growth rate for all three galaxies, which would, if continuing until redshift zero, result in final black hole masses of less than $10^{7}\,\mathrm{M}_{\sun}$.

However, after the first few time steps (9 steps or 2.8 Gyr for g1.12e12, 18 steps or 4.3 Gyr for g7.92e12, 5 steps or 1.5 Gyr for g1.05e13) the central black hole experiences one or more merger events with other black holes. Thus the black hole mass roughly doubles, which roughly quadruples the black hole accretion rate. Now the black hole accretion rate is large enough to cause significant black hole growth, and the ever increasing black hole mass leads to an ever increasing accretion rate, resulting in a kind of ‘runaway growth’. The black hole gains most of its mass due to gas accretion, the black hole mergers serve merely as a trigger for higher accretion rates. For instance, the final central black hole mass of the galaxy g7.92e12 consists to 5 per cent of merged black holes and to 95 per cent of accreted gas.

When the black hole accretion rate, and thus the feedback, is high enough, it is capable of removing significant amounts of gas from the vicinity of the black hole, indicated by a drop in the gas density in Fig. 9. An increasing black hole mass and a decreasing gas density leads the black hole accretion rate to stabilize to a roughly constant value, at least for the next few gigayears.

Previous work (Anglés-Alcázar et al., 2017; Dubois et al., 2015) indicates that black hole growth is regulated by the gas inflow from the host galaxy, and that black hole growth is only efficient after the galaxy has reached a certain mass. However, in our work the gas density does not change significantly in the first stages of black hole growth, rather, black hole mergers are crucial in triggering significant black hole growth. This implies the existence of a critical black hole mass for the onset of black hole accretion of the order of $10^{5}\,\mathrm{M}_{\sun}$. Furthermore it outlines the importance of triggering efficient black hole growth in the first place, alongside with self-regulated black hole growth in the later stages of galaxy evolution. However, we point out that these effects are merely the result of the Bondi formula being proportional to the square of the black hole mass, and might change if other accretion schemes are used.

6 Parameter study

The reference parameters of our model are the accretion parameter $\alpha=70$, the feedback efficiency $\epsilon_{\mathrm{f}}=0.05$, the black hole seed mass $M_{\mathrm{BH,s}}=1\times 10^{5}\,\mathrm{M}_{\sun}$ and the halo threshold mass $M_{\mathrm{h,t}}=5\times 10^{10}\,\mathrm{M}_{\sun}$. In this section we explore how a variation of these model parameters affect the $M_{\star}$-$M_{{200}}$ and $M_{\mathrm{BH}}$-$M_{\star}$ relations, and show that our fiducial parameters are a reasonable choice. We restrict the parameter study to the galaxies g8.26e11 and g7.92e12, the first one is the same galaxy used in the first NIHAO paper Wang et al. (2015) to calibrate the stellar feedback, while the second one has been chosen since it is one order of magnitude larger in halo mass.

In this parameter study we vary only one parameter at the same time, thus it is possible that changing multiple parameters simultaneously gives an even better match to the observed scaling relations. Further improvement could be reached by varying parameters not associated to black hole feedback, e.g., parameters related to stellar feedback, or by changing the models for black hole formation, accretion and feedback as outlined in section 3.

6.1 Accretion parameter $\alpha$

For the galaxy g7.92e12 the accretion parameter $\alpha$ seems to have a threshold value of about $\alpha=60$. For $\alpha>60$ the variation in stellar mass (Fig. 10a) does not show a clear trend but is possibly of stochastic nature. The black hole mass does not change significantly (Fig. 10b), demonstrating the self-regulating nature of black hole feedback. However, for $\alpha<60$ the black hole feedback is insufficient to reduce the stellar mass to match the observed $M_{\star}$-$M_{{200}}$ relations (Fig. 10a) and the black hole mass rapidly drops below the observed $M_{\mathrm{BH}}$-$M_{\star}$ relations (Fig. 10b). For the galaxy g8.26e11 this threshold value seems to be smaller, possibly due to a higher resolution, as all simulations only show small variations in the $M_{\star}$-$M_{{200}}$ and $M_{\mathrm{BH}}$-$M_{\star}$ relations that are possibly of stochastic nature. Thus an accretion parameter of $\alpha=70$ seems to be a reasonable choice.

6.2 Feedback efficiency $e_{\mathrm{f}}$

A higher feedback efficiency quenches the accretion of gas by the black hole, and thus leads to a lower black hole mass (Fig. 10d). The total amount of feedback energy injected into the gas for the galaxy g7.92e12 is $E=e_{\mathrm{f}}e_{\mathrm{r}}c^{2}M_{\mathrm{BH}}=(8.9,9.6,9.3,9.1,7.6)\times 10^{53}\mathrm{J}$ for the values $e_{\mathrm{f}}=(0.01,0.03,0.05,0.10,0.15)$, i.e. even changing the feedback efficiency by a factor of 15 changes the total feedback energy by less than 20 per cent. Thus a higher feedback efficiency is compensated by a lower black hole mass, leading to approximately the same total feedback energy, again demonstrating the self-regulating nature of black hole feedback. The stellar mass is slightly decreasing with decreasing feedback efficiency (Fig. 10c), but this effect is small and could also be of stochastic nature. Feedback efficiencies of 0.10 and 0.15 far overpredict the $M_{\star}$-$M_{{200}}$ relation for the galaxy g7.92e12 (Fig. 10c), and feedback efficiencies of 0.01 and 0.03 far underpredict the $M_{\mathrm{BH}}$-$M_{\star}$ relation for the galaxy g8.26e11 (Fig. 10d), making a feedback efficiency of $e_{\mathrm{f}}=0.05$ an optimal choice.

6.3 Black hole seed mass $M_{\mathrm{BH,s}}$

Small seed masses of $5\times 10^{4}$ and $7\times 10^{4}\,\mathrm{M}_{\sun}$ provide insufficient feedback, leading the galaxy g7.92e12 to highly exceed the observed $M_{\star}$-$M_{{200}}$ relations (Fig. 10e). High seed masses of $3\times 10^{5}$ and $5\times 10^{5}\,\mathrm{M}_{\sun}$ bring the galaxy g7.92e12 to the edges of the observed $M_{\star}$-$M_{{200}}$ and $M_{\mathrm{BH}}$-$M_{\star}$ relations (Fig. 10e,f), leaving a seed mass of $M_{\mathrm{BH,s}}=1\times 10^{5}\mathrm{M}_{\sun}$ as optimal choice.

According to Fig. 10e,f the final black hole mass does not depend on the black hole seed mass (except for the lowest seed mass), but the final stellar mass does. This can be explained with Fig. 11, which shows the black hole mass and the stellar mass as a function of time for the galaxy g7.92e12 and for different black hole seed masses. The three lowest black hole seed masses show a black hole mass evolution as outlined in section 5: The growth rate is low at first, then black hole mergers trigger ‘runaway growth’  that is then stopped by feedback and ultimately leads to black hole masses in accordance with observations.

The two highest black hole seed masses show an entirely different evolution: Black hole mergers are not needed, the seed mass is already sufficiently high to trigger ‘runaway growth’.

Although the evolution of the black hole masses in both scenarios is entirely different at early times, they produce the same final black hole masses (except for the lowest seed mass), which can be attributed to the self-regulating nature of black hole feedback: Large black holes exert high feedback and allow for a low accretion rate, small black holes exert low feedback and allow for a high accretion rate.

Fig. 11 also shows the stellar mass as a function of time for different initial black hole seed masses. Until 2.5 gigayears the stellar masses evolve almost identical, then they start diverging: The higher the black hole seed mass the earlier the black hole accretion rate peaks, thus star formation is quenched earlier, which sets the stellar masses to ever lower values. The sudden jumps in the stellar mass at about 4-5 gigayears are caused by a merger event.

Different black hole seed masses lead to the same black hole mass, but to different stellar masses at redshift zero. However, calibration on scaling relations can be used to determine the optimal choice for the black hole seed mass.

6.4 Halo threshold mass $M_{\mathrm{h,t}}$

High halo threshold masses of $7\times 10^{10}$ or $7\times 10^{10}\,\mathrm{M}_{\sun}$ overpredict the $M_{\star}$-$M_{{200}}$ relation (Fig. 10g) for the galaxy g7.92e12, whereas decreasing the halo threshold mass to $1\times 10^{10}$ or $3\times 10^{10}\,\mathrm{M}_{\sun}$ moves the galaxy to the edge of the observed $M_{\mathrm{BH}}$-$M_{\star}$ relations (Fig. 10h). Thus a halo threshold mass of $M_{\mathrm{h,t}}=5\times 10^{10}\mathrm{M}_{\sun}$ is the most reasonable choice.

The results in Fig. 10g,h resemble the results in Fig. 10e,f: Increasing the halo threshold mass is equivalent to reducing the black hole seed mass, thus the same ratio of halo threshold mass to black hole seed mass produces the same results. This can be explained as follows: The halo threshold mass basically determines the time the black hole is seeded. If with our standard parameters a black hole with $10^{5}\,\mathrm{M}_{\sun}$ is seeded at time t, increasing the halo threshold mass means the black hole with $10^{5}\,\mathrm{M}_{\sun}$ is seeded at a time t+dt, and reducing the black hole seed mass means a black hole with $10^{4}\,\mathrm{M}_{\sun}$ is seeded at time t. The former is equivalent to the latter, as the black hole with $10^{4}\,\mathrm{M}_{\sun}$ seeded at time t will have grown to $10^{5}\,\mathrm{M}_{\sun}$ at the time t+dt.

7 Summary

We introduce and test algorithms for black hole formation, accretion and feedback to the NIHAO project. For black hole formation we place a black hole in the centre of a halo once it exceeds a threshold mass, for black hole accretion we use the Bondi-Hoyle-Lyttleton parametrization, and for black hole feedback we deposit thermal energy, which is proportional to the black hole accretion rate, into the gas that surrounds the black hole. This addition to the NIHAO project allows us to extend the NIHAO suite of galaxies to higher masses.

Our galaxies show good agreement with the observed $M_{\star}$-$M_{{200}}$ and $M_{\mathrm{BH}}$-$M_{\star}$ relations that we use to calibrate the free parameters of our model. We also investigate the scatter of these relations and their time evolution. The scatter of both relations is decreasing with time for $z<1$ (higher redshifts possibly suffer from low number statistics and seeding effects), and is also lower than the observed scatter, possibly because we measure the intrinsic scatter of these relations without any observational uncertainties.

In the high mass, elliptical galaxies the quenching of star formation occurs after an increase of the black hole accretion rate, confirming that star formation is quenched by the black hole feedback. A parameter study confirms that we have chosen the optimal parameters within the framework of our model. Our simulations provide a valuable tool to study the effect of black hole feedback on galaxy formation and evolution.

Acknowledgements

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). A part of this research was carried out on the High Performance Computing resources at New York University Abu Dhabi. We used the software package pynbody Pontzen et al. (2013) for our analyses. AO is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - MO 2979/1-1. We thank Tobias Buck for supplying the initial conditions for g7.92e12, g1.05e13 and g1.44e13.