key: cord-0426912-pxcju0bt
authors: cCatmabacak, Onur; Feldmann, Robert; Angl'es-Alc'azar, Daniel; Faucher-Giguere, Claude-Andr'e; Hopkins, Philip F.; Science, Duvsan Kerevs Institute for Computational; Zurich, University of; Zurich,; Switzerland,; Physics, Department of; Connecticut, University of; Storrs,; CT,; USA,; Astrophysics, Center for Computational; Institute, Flatiron; York, New; NY,; Astronomy,; CIERA,; University, Northwestern; Evanston,; IL,; TAPIR,; Technology, California Institute of; Pasadena,; CA,; Astrophysics, Center for; Sciences, Space; Diego, University of California at San; Jolla, La
title: Black hole -- galaxy scaling relations in FIRE: the importance of black hole location and mergers
date: 2020-07-23
journal: nan
DOI: 10.1093/mnras/stac040
sha: 5145586d0b8d90dc15156c6f6c87fb9788ac0c72
doc_id: 426912
cord_uid: pxcju0bt

The concurrent growth of supermassive black holes (SMBHs) and their host galaxies remains to be fully explored, especially at high redshift. While often understood as a consequence of self-regulation via AGN feedback, it can also be explained by alternative SMBH accretion models. Here, we expand on previous work by studying the growth of SMBHs with the help of a large suite of cosmological zoom-in simulations (MassiveFIRE) that are part of the Feedback in Realistic Environments (FIRE) project. The growth of SMBHs is modelled in post-processing with different black hole accretion models, placements, and merger treatments, and validated by comparing to on-the-fly calculations. Scaling relations predicted by the gravitational torque driven accretion (GTDA) model agree with observations at low redshift without the need for AGN feedback, in contrast to models in which the accretion rate depends strongly on SMBH mass. At high redshift, we find deviations from the local scaling relations in line with previous theoretical results. In particular, SMBHs are under-massive, presumably due to stellar feedback, but start to grow efficiently once their host galaxies reach $M_* sim 10^{10} M_{odot}$. We analyse and explain these findings in the context of a simple analytic model. Finally, we show that the predicted scaling relations depend sensitively on the SMBH location and the efficiency of SMBH merging, particularly in low-mass systems. These findings highlight the relevance of understanding the evolution of SMBH-galaxy scaling relations to predict the rate of gravitational wave signals from SMBH mergers across cosmic history.

cally well constrained, their observational status at higher redshift is less clear, with different authors suggesting both redshift-dependent (Treu et al. 2004; Walter et al. 2004; Merloni et al. 2010; Targett et al. 2012; Netzer & Trakhtenbrot 2014; Bongiorno et al. 2014) and redshift-independent relations (Shields et al. 2003; Jahnke et al. 2009; Cisternas et al. 2011; Ding et al. 2017) .

Even though SMBHs and galaxies follow relatively tight scaling relations in the local Universe, it is currently unknown whether such tight relations hold in the early Universe (Huang et al. 2018; Trakhtenbrot et al. 2017; Delvecchio et al. 2019; Shirakata et al. 2016; Izumi et al. 2018) . In particular, SMBHs at high redshift may be over-or under-massive compared to their host galaxies or could grow in lock-step with each other (Volonteri 2012) . The redshift evolution and the scatter of various SMBH-galaxy scaling relations may provide critical insights into the physics of black hole and galaxy growth. Which physical processes might be responsible for reproducing the local scaling relations? Is it possible to produce the local scaling relations without self-regulating black hole feedback? If so, how do the SMBH-galaxy scaling relations evolve at high redshift? These are the questions we would like to address in this paper.

The standard approach to model the growth of SMBHs is via the spherical accretion approximation (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952) . In its basic form, Bondi-like models assume radial accretion of non-self-gravitating gas onto a point-like source to estimate the accretion rate from large scales to black holes at the centre of galaxies. However, there are caveats to these prescriptions. Springel et al. (2005) and Booth & Schaye (2009) introduced an ad hoc boost factor of the Bondi model to avoid underestimating the accretion rate. Furthermore, the Bondi ansatz does not account for the angular momentum of the inflowing gas (Hopkins & Quataert 2010 Angles-Alcazar et al. 2020) . However, modifications of the Bondi model that include rotation have been proposed by, e.g. Hobbs et al. (2012) ; Tremmel et al. (2017) .

On their own, Bondi-like models tend to overpredict the SMBH mass since they scale with 2 BH . Black hole feedback is thus critical as it avoids overly massive SMBHs relative to the local scaling relations by regulating both black hole growth and star formation Di Matteo et al. 2005; Sijacki et al. 2007 ). This idea has been widely used to investigate the evolution of galaxies and SMBHs in cosmological simulations such as I (Vogelsberger et al. 2014; Genel et al. 2014 ), H -AGN (Dubois et al. 2014; Volonteri et al. 2016; Kaviraj et al. 2017 ), E (Schaye et al. 2015) , M B (Khandai et al. 2015) , B T (Feng et al. 2016 ), R (Tremmel et al. 2017) , and I -TNG (Springel et al. 2018 ).

On the other hand, alternative models for the gas accretion onto SMBHs have been proposed. Hopkins & Quataert (2010) performed nested simulations of star-forming galaxies to understand how gas can accrete from galactic scales (∼ 10 − 100 kpc) to smaller scales (< 1 pc). Non-axisymmetric features that result in gravitational torques caused by galaxy mergers, spiral instabilities and eccentric disc modes efficiently remove angular momentum of the gas and drive it further down to the sub-pc scales (Hopkins & Quataert 2011) . This model has been successfully used in galaxy simulations to reproduce the local scaling relations without the need for selfregulatory black hole feedback (Anglés-Alcázar et al. 2013 Davé et al. 2019 ; Thomas et al. 2019) .

The present work studies the SMBH growth in a fully cosmological context with the help of high-resolution, zoom-in, hydrodynamical simulations. The simulations used in this paper (M FIRE) resolve scales down to tens of pc in a cosmological environment. High resolution is essential to properly trace the flow of gas into the centres of galaxies. Furthermore, the relatively large number of well-resolved galaxies in our sample (43 different galaxies at = 6 and 132 at = 2) compared to full cosmological simulations allows us to address the questions listed above with a statistically significant set of simulated galaxies over a wide range of redshifts (2 ≤ ≤ 12) and halo masses (10 < log( h / ) < 13.5). The outline of the paper is as follows; section §2 introduces the simulation properties. Section §3 lays out the details of our postprocessing analysis. We present our main results in the following section §4. Specifically, section §4.4 presents a toy model to explain the physical origin of the * − BH scaling relation. We discuss the caveats of our post-processing analysis in §5 and give our summary in §6.

We use 34 high-resolution, cosmological zoom-in simulations from the M FIRE suite (Feldmann et al. 2016 Feldmann 2017; Anglés-Alcázar et al. 2017c ) that is part of the Feedback in Realistic Environments (FIRE 1 ) project (Hopkins et al. 2014 (Hopkins et al. , 2018 . Simulations were run with the gravity-hydrodynamics solver GIZMO 2 in Pressure-Energy Smoothed Particle Hydrodynamics (P-SPH, FIRE-1) and Meshless Finite Mass (MFM, FIRE-2) mode. Table 1 provides an overview of the simulations used in the present work.

The selection of the zoom-in regions for runs from series A, B, and C is described in detail in Feldmann et al. (2017) . In brief, isolated halos are selected from a low-resolution DM-only run of an = 100 Mpc h −1 comoving cosmological volume. The halos are selected based on their = 2 masses (we consider 3 narrow mass bins corresponding to 2.5 − 3.6 × 10 12 , 0.9 − 1.1 × 10 13 , and 2.5 − 3.6 × 10 13 ) and local environmental densities (based on the enclosed mass within a 1.8 pMpc radius). In total, 18 haloes are selected with a range of masses (10, 5, and 3 haloes from the low, intermediate, and high mass bins) and environmental densities. Initial conditions for the zoom-in runs were generated using the multi-scale initial condition tool MUSIC (Hahn & Abel 2011 ) using a convex hull for all particles within 3 × vir at = 2.

Additional zoom-in simulations (series D and E) are created in a similar fashion from low-resolution DM-only simulations of = 400 Mpc h −1 and = 762 Mpc h −1 comoving cosmological volumes. The 3 most massive halos at = 6 are selected from each volume. In addition, 5 zoom-in regions are created based on the = 400 Mpc h −1 volume by selecting halos with halo ( = 6) > 10 12.5 from a range of local environmental densities. Another zoom-in region is selected from the same volume based on having a halo mass at = 2 of approximately 10 14 . FIRE-1 simulations use a quintic spline kernel with 60-62 neighbors for gravitational softening (Morris 1996; Dehnen & Aly 2012) , see Hopkins et al. (2014) . The gravitational softening lengths of dark matter and star particles are fixed at 143 and 21 pc (physical) respectively, while the softening length of gas particles is adaptive and reaches a minimum value of 9 pc (physical) in the dense interstellar medium. FIRE-2 simulations use a cubic spline kernel with 32 neighbors (Morris 1996 FIRE-2 100 1 13.10 5 Table 1 . List of simulations used in this work. Column 1 refers to the name of the simulation, see Feldmann et al. (2017) . Column 2 lists whether simulations were run with FIRE-1 or FIRE-2 physics. Column 3 provides the box sizes from which the zoom-in simulations were selected (in comoving units). (2018). The gravitational softening lengths of dark matter and star particles are 57 and 7 pc. The minimum softening length of gas particles is 0.7 pc. In all simulations, the gas softening lengths are chosen sufficiently small to capture gas densities well above the star formation threshold. All simulations have a mass resolution of 1.7 × 10 5 for dark matter particles and 3.3 × 10 4 for gas and star particles.

Star formation takes place only in self-gravitating, dense, molecular gas with a density above 5 and 1000 atoms per cm 3 for FIRE-1 and FIRE-2 simulations, respectively. The simulations include various stellar feedback channels such as energy, momentum, and mass injection from stellar winds and supernovae, local and long-range momentum flux from radiative pressure, a uniform UV background using the model from (Faucher-Giguère et al. 2009 ) and photo-ionization and photo-electric heating (Hopkins et al. 2014 (Hopkins et al. , 2018 .

The growth of black holes is modelled fully in post-processing. Our FIRE-1 simulations do not directly account for black hole physics, while the FIRE-2 runs include live black hole sink particles but do not model AGN feedback. However, we adopt the same post-processing approach for FIRE-1 and FIRE-2 simulations in the present study. A comparison between the prediction of our post-processing model and the on-the-fly calculation is shown in appendix A.

We refer the reader to Feldmann et al. (2016 Feldmann et al. ( , 2017 ; Anglés-Alcázar et al. (2017c); Hopkins et al. (2014 Hopkins et al. ( , 2018 for more detailed information about the simulations and the properties of the simulated galaxies.

In this section, we will introduce our post-processing approach and describe the various SMBH accretion models studied in this work. By considering different models of black hole accretion, we can analyze how the resultant SMBH masses change in the absence of self-regulating AGN feedback. As we will show in subsequent sections, the choice of the accretion model affects both the mass evolution of SMBHs as well as the resulting * − BH scaling relation. We follow the exact same post-processing approach in both our FIRE-1 and FIRE-2 simulations.

We use the publicly available Amiga Halo Finder 3 (AHF, Knollmann & Knebe 2009 ) to identify dark matter haloes and to find their centres in the M FIRE simulations. The identified virialized structures contain at least 100 particles ( halo ∼ 10 7 ℎ −1 ). SMBHs are placed either at the centre of mass (COM) or the maximum density centre (MAX) of each halo as provided by AHF. The former is defined as the centre of mass of the gas, star, and dark matter particles on the finest level of refinement of the host halo. Therefore, the COM often represents a typical environment in the central region of the host halo. The MAX is calculated as the position of the highest-density cell within the halo by AHF. These positions typically correspond to dense star clusters or gas clouds. Clearly, the SMBH placement can have a substantial impact on the early SMBH growth, in addition to the role played by stellar feedback (Anglés-Alcázar et al. 2017c) .

The difference between the two centering approaches is shown in physical units in Figure A3 for simulation A1 with FIRE-2 physics. The center positions differ by 0.5−1.0 kpc at relatively high redshift ( 5) and the difference becomes much smaller (< 100 pc) at later times. This behavior is potentially linked to the transition in galactic structure from an irregular morphology at high redshift to well-settled disc galaxies at later times (Sparre et al. 2017; Stern et al. 2020) .

Our sample consists of the most massive halo in each simulation and of all the haloes above halo = 10 10 at the final simulation snapshot in the 6 simulations that run to 1 (FIRE-1 A1 and A4, and FIRE-2 A1, A2, A4, and A8). Our selection is subject to the following contamination condition. The mass fraction of high-resolution dark matter particles is required to be larger than 98% so that haloes are not significantly polluted with low-resolution dark matter particles. Progenitors of selected halos are traced back Table 1 at = 2 and = 6. Our sample includes the central galaxy of the primary halo from 34 different zoom-in simulations and all galaxies in halos above halo = 10 10 from 6 zoom-in simulations evolved to 1. The total number of distinct galaxies followed across cosmic time varies with redshift (132 at = 2, 43 at = 6). Our sample includes a broad range of low to moderately massive galaxies and halos.

in time with the AHF M T tool. We only consider progenitors with a stellar mass of at least 10 7 within 10% of their virial radii (Moore et al. 1998; Devriendt et al. 2010) . Progenitors of lower stellar mass do not host a SMBH in our model and are thus neglected in this study. Figure 1 shows the distribution of stellar masses in our sample at = 2 and = 6. The stellar mass range of M FIRE galaxies is 7.0 < log * / < 11.5 at ∼ 2 and 7.0 < log * / < 10.6 at ∼ 6 while the halo mass range is 8.3 < log * / < 12.1 at ∼ 2 and 8.4 < log * / < 11.6 at ∼ 6. In total, we have 4681 galaxy data sets corresponding to ∼ 100 distinct galaxies traced across cosmic time ( < 12). We refer the reader to Feldmann et al. (2016) ; Price et al. (2017); Feldmann et al. (2017) ; Cochrane et al. (2019) ; Wellons et al. (2020) ; Parsotan et al. (2021) for the general properties of galaxies in M FIRE simulations. The radii of central galaxies, gal , are defined as 10% of the virial radii of their parent halos, e.g., Price et al. (2017) . The total stellar mass of a given central galaxy is subsequently calculated as the mass of stellar particles within gal . We separate the stellar mass of the galaxy into a bulge and a disc component following Anglés-Alcázar et al. (2014) . Specifically, the bulge mass in an enclosed region with a given radius is calculated as twice the mass of all star particles that are counter-rotating ( < 0) according to the unit stellar angular momentum vector of the galaxy. The stellar disc mass within is then defined as the difference between the stellar mass within and the bulge mass within . The total disc mass is defined as the sum of the stellar disc mass and the gas mass within . We calculate the velocity dispersion as the velocity dispersions of all star particles that belong to the stellar bulge. The comparison of velocity dispersion within different radial apertures ( 0 ) is shown in the appendix A.

The evolution and growth of SMBHs are treated in postprocessing. A (virtual) SMBH seed is placed at the centre of each progenitor halo. This choice establishes early SMBHs on the local * − BH relation when 10 4 SMBH seed masses are adopted (our fiducial option). Lower or higher seed masses result in highredshift SMBHs that start either below or above the local * − BH scaling relation. Also, the velocity of each SMBH is defined as the average velocity of the 100 youngest star particles in its vicinity (within 0 , see §3.2.1). The post-processing analysis is performed in the rest-frame of the SMBH of the given galaxy.

We define the stellar growth rate (SGR = dM * /dt) as the past change in stellar mass over a time interval avg in a three-dimensional sphere of fixed physical size divided by avg . In practice, we measure the difference in stellar mass in the considered physical volume between the current snapshot at time and an earlier snapshot near time − avg . We thus do not calculate SGRs when avg is much shorter than the time interval between snapshots (10 − 25 Myr). In contrast, the star formation rate (SFR) is defined as the stellar mass belonging to recently (within time avg ) formed stars currently present in a three-dimensional sphere divided by avg . It can be thus be calculated for both long and short avg . Both SFRs and SGRs are subject to stellar mass loss.

We need to differentiate between mergers of SMBHs and mergers of the dark matter halos in which they reside as not every halo merger results in the merger of their central galaxies and not every galaxy merger results in a prompt merger of their SMBHs. We thus analyse three different SMBH merger scenarios for the remainder of this paper.

Our first model (" ") assumes that SMBHs merge as soon as their parent dark matter haloes merge, regardless of the mass ratio of the two haloes. This scenario results in the maximum possible number of SMBH mergers. Hence, this extreme scenario explores the most optimistic contribution of BH merging to SMBH growth.

Our second model (" ") poses more stringent constraints on SMBH mergers and is based on the distance sep of the two parent halos when they are identified as separate (sub-)halos for the last time ( merge ). Specifically, the SMBHs at the centres of these halos are assumed to have merged by the next snapshot if (i) their sep is smaller than 10% vir of the more massive of the two halos or (ii) their dynamical friction timescale df is smaller than the Hubble time at merge . We adopt the following analytic estimate of the dynamical friction time (Binney & Tremaine 1987) :

Here, h1 and h2 are the virial masses of the more massive and the less massive halo, respectively. Furthermore, Λ is Coulomb logarithm where ln(Λ) = ln(1 + h1 / h2 ). vir is the virial radius and vir is the circular velocity of the more massive of the merging haloes. The less massive halo is often affected by tidal stripping at the late stages of the halo merger. To mitigate the effect of stripping -which may extend the inferred merging timescale-the maximum mass of the progenitors of the less massive halo is used instead of its current halo mass. A halo typically reaches its maximum mass 60 − 100 Myrs before the halo merger.

Our third option (" ") is to neglect SMBH mergers altogether so that SMBHs only grow via mass accretion throughout cosmic time.

Hopkins & Quataert (2011) introduced a model of gas accretion from kpc to sub-pc scales driven by gravitational torques. The functional form of the black hole accretion rate (BHAR) in the gravitational torque-driven accretion (GTDA) model is

In the GTDA model, the accretion rate is calculated based on the properties of gas particles inside of a sphere with a radial aperture ( 0 ). All the terms in Equation 2 have a radial aperture ( 0 ) dependency except the SMBH mass, m , and T . Here, T is a function of nuclear star formation law (see equations 39 and 65 in Hopkins & Quataert (2011) ), however for practical reasons we set T to 5 (see Figure 10 in Hopkins & Quataert (2011) ). The mass retention rate m is a fudge factor between 0 and 1 which reduces the analytically-derived accretion rate of the GTDA model. Physically, m captures the reduction in accretion rate due to unresolved winds. The product of m and T is an overall normalisation that covers the effects on gas dynamics (i.e. stellar and BH feedback) at unresolved scales (Anglés-Alcázar et al. 2017a).

The disc fraction d is the ratio of the disc mass (stellar disc mass + gas mass) to total baryonic mass (stellar mass + gas mass). The accretion rate scales super linearly with the disc fraction, ∝ 5/2 d . BH is the mass of the SMBH and tot is the total mass within the radial aperture ( 0 ). The total mass is the sum of dark matter and baryonic matter inside 0 . However, especially at the central regions, baryonic matter dominates over dark matter. Although the accretion rate scales only linearly with the total mass inside 0 , the total mass is the determining factor for the accretion rate onto SMBH. Amongst the remaining parameters, gas is the ratio of the gas mass to the total baryonic mass inside 0 and 0 ≈ 0.31 2 d ( tot /10 9 ) −1/3 . We update black hole masses in our post-processing analysis iteratively for each accretion model, as briefly explained in §3.2, BH,i+1 = BH,i + Δ × BH,i . We repeat the same analysis for the densest (MAX) and average density (COM) centres, different SMBH merger treatments (with or without SMBH mergers) and radial apertures ( 0 ) changing from 1 kpc to 100 pc. The BHAR is limited to the Eddington rate for the Bondi-like models (see below) and to ten times the Eddington limit for all other accretion models.

Spherical accretion onto a point object has a solution known as Bondi accretion, (Bondi & Hoyle 1944; Bondi 1952) and the corresponding accretion rate is

In the equation above, B is the boost factor introduced by Springel et al. (2005) , is the gravitational constant, and is the volume density of gas particles within 0 . The bulk velocity of gas and the sound speed are denoted by bulk and s , respectively.

The Bondi model is valid for the case of hot virialised gas with negligible angular momentum and radiative cooling. It does not account for SMBH growth via accretion of high density gas which tends to cool efficiently. Hobbs et al. (2012) proposed a modification to the Bondi model to account for the contribution of the halo to the gas dynamics. They replaced the relative velocity with the velocity dispersion for the external potential, ∼ √︁ enc ( )/ , and the SMBH mass with the enclosed mass of the external potential.

Observational and theoretical evidence for a roughly constant ratio between BHAR and SFR rate has led to the idea that the growth of SMBHs and galaxies are coupled, especially at the nuclear scales (Hopkins & Quataert 2010; Volonteri et al. 2015; Dai et al. 2018a; Yang et al. 2017) . A simple ansatz is to model the SMBH accretion rate as a linear function of the SFR inside R 0 :

Finally, we also use a simple accretion model for comparison where the BHAR scales with the free-falling gas inside R 0 (Anglés-Alcázar et al. 2017c):

Here, is a scaling factor that controls the percentage of free-falling gas accreted onto SMBHs. The typical values of changes from 0 to 100%. The free-fall timescale of the gas is dyn = 3 0 /(2 tot (< 0 )).

In this section, we analyse how the choice of our model parameters affects the predicted SMBH growth. Our post-processing analysis includes several key parameters that influence the growth of SMBHs such as the SMBH seed mass ( seed ), the black hole merger treatment and, most importantly, the accretion model (Equation 2 -6) and its free parameters, e.g. the mass retention rate ( m ) and the radial aperture ( 0 ) for the GTDA model. We will use the fiducial values for our model parameters in Table 2 together with the MAX centring method and the model to see the precise impact of our model parameters. Throughout this paper, these settings (the fiducial parameters from Table 2 , the MAX centring method and the model) will be our primary choice unless stated otherwise.

We tested various values for the radial aperture ( 0 ) in the absence of self-regulating AGN feedback. As the top panel of Figure 2 shows, the growth history of SMBHs is only mildly affected if 0 is varied between 100 pc and 1 kpc. This finding is perhaps surprising, given that the GTDA model has a strong dependence on the radial aperture ( 0 ), ∝ −3/2 0 . However, various other terms in Equation 2 also depend on the radial aperture ( 0 ), largely cancelling the overall dependence on 0 in agreement with similar tests in Anglés-Alcázar et al. (2015) .

As shown in Hopkins & Quataert (2010) , a smaller 0 results in a more precise prediction (lower scatter) of the instantaneous accretion rate on small scales. Hence, to mimic the accretion from galactic scales to sub-pc scales as accurately as we can, we adopt 0 = 100 pc as the fiducial value.

The mass retention rate is a normalisation of the overall gas accretion rate. As the bottom panel of Figure 2 highlights, varying the mass retention rate creates a noticeable shift in the normalisation of the SMBH mass, at any redshift. Hence, we can adjust the normalisation of the predicted M * − M BH scaling relation by choosing an appropriate value of the mass retention rate (Anglés-Alcázar et al. 2013). We adopt a mass retention rate of 10% as our fiducial value. The fiducial values in Table 2 are chosen so that SMBH masses in our post-processing analysis are in approximate agreement with the local * − BH scaling relation. Figure 3 shows the insensitivity of the GTDA model predictions to the black hole seed mass choice (Anglés-Alcázar et al. 2013 . Specifically, SMBH seeds with masses 10 − 10 5 result in a similar SMBH mass by ∼ 4.5 if accretion is limited to less than ten times Eddington and by ∼ 3 in the case of Table 2 . SMBHs are placed in the densest centres (MAX), and model is used to model SMBH growth. A larger 0 slightly increases the SMBH mass, but overall the choice of 0 has little impact on the growth history of the SMBH. In contrast, the SMBH scales approximately linearly with m .

Eddington-limited accretion. This convergence is a consequence of the SMBH accretion rate being only a weak function of black hole mass in the GTDA model. This figure suggests that observations of SMBH masses in the progenitors of massive galaxies at > 3 − 5 may provide useful constraints on the masses of the first black hole seeds. Figure 4 compares the SMBH growth via gas accretion with the contribution from black hole seed masses for the A-series simulations run with FIRE-2 physics (Anglés-Alcázar et al. 2017c). For low-mass SMBHs ( BH < 10 6 ) at 3 < < 6, the seed mass contribution can be significant if relatively heavy seeds (here 10 4 ) are chosen. However, most of the SMBH mass is acquired via gas accretion either in-situ or by merging at lower redshifts. We refer the reader to Figure A2 for the same analysis with different seed masses. Figure 5 compares the results for SMBH growth via GTDA . Effect of the SMBH seed mass choice and a limit on the growth rate on the mass evolution of SMBHs for our fiducial settings with the model and SMBHs located in the densest centres (MAX). The top (bottom) panel shows the prediction if the growth rates are limited to ten times the Eddington rate (to the Eddington rate). The masses of SMBHs converge to ∼ 10 7 after 1.5-2 Gyr of cosmic time (by z ∼ 3 − 4) independent of the initial seed mass. The influence of the seed mass on the growth of SMBHs is thus limited to high redshift in the progenitors of massive galaxies, and it is smaller if super-Eddington accretion rates are possible.

with the predictions of various other accretion models employed in the literature (see section 3.2). These models can be divided into two types. The first group includes models in which the accretion rate depends strongly on the black hole mass, such as in Bondi accretion. In this case, black holes may grow extremely fast in the absence of AGN feedback. As a consequence, the predicted black holes are overly massive compared to the local scaling relations.

In the second group of models, accretion rates scale weakly with SMBH mass. These models use the properties of the host galaxy such as stellar mass, disc fraction, free-fall timescale of gas and SFR. Such models result in a much more steady black hole growth in massive galaxies from an initial seed mass of 10 4 at > 7 to SMBHs of 10 8 or more by = 1. The most significant difference between these two groups of accretion models is the necessity of self-regulating AGN feedback to reproduce the (10 4 ). Overall, gas accretion is the driving force for black hole growth throughout much of cosmic history. In the absence of self-regulating AGN feedback, the Bondi accretion model and its modification by Hobbs et al. (2012) result in overly massive black holes at < 4. In contrast, SMBHs grow more steadily in GTDA, SFR-based, and dynamical free-fall models, and, for appropriate choices of the overall normalisation, they lead to similar mass growth histories.

local scaling relations. The Bondi-like accretion models have a high dependence on the SMBH mass, thus require strong AGN feedback to regulate SMBH growth to provide reasonable SMBH masses (Anglés-Alcázar et al. 2015) . Other models can produce matching results without the need for expelling material from the centre of the host galaxy. Figure 6 compares the gas fraction within the central 100 pc region and the whole galaxy for the simulation A1 with FIRE-2 physics. At high redshift, gas fraction demonstrates a bursty behaviour until stellar feedback becomes inefficient to remove gas from the central region of galaxy (Muratov et al. 2015; Anglés-Alcázar et al. 2017c) . Gas fraction in the centre region follows the increase in the galactic gas fraction. This suggests that gas inflow at high redshift can reach to the central regions more easily compared to the galaxies in the local Universe. The change in the galaxy size could play an important role for the bursty behaviour of gas fraction (Torrey et al. 2017; Faucher-Giguère 2018) . Furthermore, peaks in gas fraction decrease rapidly at early times. The replenished gas reservoir triggers expeditious SF and feedback from the newly formed stars evacuates gas from the star forming regions. We refer the reader to Figure A4 for a comparison of the gas fraction for different centring methods. Figure 7 shows the relation between the stellar mass of the host galaxies and SMBH masses for our fiducial set of parameters as well as for different assumptions regarding the SMBH location and merger treatment. In all cases, SMBHs reach similar masses ( BH 10 7.5 ) in massive galaxies ( * 10 11 ). Different combinations of SMBH location and merger treatment result in different tracks ( see Volonteri 2012 for details) in the * − BH plane, notably at high redshift. For each case, efficient SMBH growth starts around a similar stellar mass threshold ∼ 10 10 . The growth trend of SMBHs shown in Figure 7 can be divided into three phases. During phase I ( * < 10 8.5

), the contribution of black hole mergers to the total SMBH mass is negligible. Instead, early SMBH growth is driven by accretion from the densest central region (upper panels in Figure 7 ). In this scenario, the growth of SMBHs and their host galaxies follow the local scaling relation for AGN galaxies in Reines & Volonteri (2015) . In contrast, the early growth of SMBHs is shifted if they reside in average density central regions, lower panels in Figure 7 .

During phase II (8.5 < log( * / ) < 10), SMBH mergers play an important role. In particular, SMBH mergers contribute more to the growth of the SMBH mass than gas accretion. In contrast, SMBH growth stalls during this phase if SMBH mergers are not considered. In our sample, the contribution from mergers peaks at 2 < < 6 (see Figure 4 ). The impact of mergers is more pronounced when the seed mass is heavier. We refer the reader to Appendix A for the comparison of different seed masses.

The SMBH merger treatment and the seed mass choice appear not to have a major impact on the SMBH mass at low redshift in massive galaxies. In each case, we find ∼ 10 8 SMBHs in ∼ 10 11 galaxies, see Figure 7 . While this result disagrees with other theoretical studies Shirakata et al. (2016) ; Park et al. (2016) , we suspect the origin of this difference is in the modelling of the SMBH accretion and stellar feedback. For any accretion model that has a black hole mass dependency in the form of BHAR ∝ p BH for > 0 (e.g. = 1/6 for the GTDA model), the specific BHAR (sBHAR ≡ BHAR/M BH ) scales with p−1 BH . In the case of Bondilike models, p is equal to 2 and sBHAR scales linearly with the black hole mass. This means that heavier black holes grow faster than smaller black holes in Bondi-like models for the same environmental conditions. SMBH masses are asymptotically insensitive to the seed mass if BH accretion is sublinear ( < 1) (Anglés-Alcázar et al. 2015 . In contrast, BH growth histories are highly sensitive to the initial seed mass in the case of super-linear accretion ( > 1), especially in the absence of AGN feedback and/or at early times when BH ∼ seed . We note that the presence of feedback may (Taylor & Kobayashi 2014) or may not (Dubois et al. 2015) affect the sensitivity of SMBH growth on the seed mass.

Phase III of the SMBH growth starts with accelerated SMBH growth when the stellar mass of the host galaxy reaches ∼ 10 10 , roughly coinciding with the time when the escape velocity of the central region becomes comparable to the velocity of galactic winds (Anglés-Alcázar et al. 2017b). The transition around the threshold mass is an empirical result from the simulations (Dubois et al. 2015; Habouzit et al. 2017; McAlpine et al. 2018) . Interestingly, the threshold stellar mass for efficient SMBH growth is similar to the divider between early and late type galaxies seen in both observations (Kauffmann et al. 2003 ) and simulations (Bower et al. 2017a; Taylor et al. 2017) . We refer the reader to Byrne et al. (in prep; see also the discussion in Stern et al. 2020 ) for a more detailed analysis of the physical drivers of delayed vs. efficient SMBH growth. Finally, once galaxies reach a stellar mass of * 10 11 ; this accelerated growth comes to an end. After this time, galaxies and SMBHs grow again at a similar rate.

The relation between the stellar bulge and SMBH mass shown in Figure 8 is qualitatively similar to the relation between galaxy stellar mass and SMBHS mass ( Figure 7) . Typically, the b − BH relation predicted by our sample at high redshift falls below the local b − BH scaling relation Häring & Rix (2004) . Our postprocessing analysis thus predicts that black hole masses are lower (or bulge masses and galaxy masses are larger) than expected from the local scaling relations.

We refer the reader to Figure A7 , 0 = 100 pc, and m = 10%), for different SMBH locations and merger treatments using the complete M FIRE sample. Colour reflects the number of galaxies in each pixel for all simulation snapshots. The panels show SMBHs from their time of seeding (z ∼ 6 − 12) until the final redshift of each simulation, see Table 1 . In each panel, the solid red line shows the * − BH scaling relation for disc galaxies in Reines & Volonteri (2015) . The first row of panels use maximum-density centre and the second row of panels use centre-of-mass as the halo centre. Columns indicate * − BH scaling relation for different merger treatments. Early growth of SMBHs is suppressed when SMBHs are placed at the centre of mass of the halo (bottom panels), which is more strongly affected by stellar feedback compared to the maximum-density centre. The maximum-density location (top panels) result in SMBH growth at high redshift, in line with the local scaling relation. SMBH mergers make a considerable contribution to the total SMBH mass in intermediate-mass galaxies (8.5 < log( * / ) < 10). Efficient SMBH growth starts when the stellar mass reaches ∼ 10 10 . SMBHs grow at a similar rate as their hosts in massive galaxies ( * > 10 11 ). In all six cases, SMBHs end up with similar masses. Panel A is the most optimistic scenario for SMBH growth and results from Panel E and F are in line with the findings of Anglés-Alcázar et al. (2017c), in which early SMBH growth is suppressed due to strong stellar feedback. We refer the reader to Figure A5 to see the effect of different seed mass choice on the * − BH scaling relation.

understood from Figure A3 which shows that at late times ( 4) the centre of the galaxy becomes well defined. Both choices lead to virtually identical SMBH placements.

Our standard approach to determine the velocity dispersion assumes non-rotating bulges (see §3 for details). We also show the effect of replacing the velocity dispersion of the bulge with the velocity dispersion of all star particles in 1 kpc, within the half stellar mass radius, and within the galactic radius in Figure A1 , finding little difference.

We convert our SFRs and BHARs into IR and bolometric luminosities via the following conversions: IR = SFR × 1.49 × 10 10 (Kennicutt 1998) and bol = /(1− ) × 2 , see Figure 9 . Here is the radiative efficiency of the accreting black hole which is generally taken as 10%. In general, we find good agreement between our model predictions and observations at < 2 and large differences at high redshift. Specifically, Figure 9 compares SFR and BHAR of our sample with the ones in galaxies hosting luminous AGN for ≤ 4 (left-hand Comparison of IR and bolometric luminosities derived from star formation and mass accretion rates with the observational data from literature Priddey et al. 2003; Wang et al. 2011; Xu et al. 2015; Fan et al. 2016; Gruppioni et al. 2016; Netzer et al. 2016; Duras et al. 2017; Bischetti et al. 2018; Díaz-Santos et al. 2018; Izumi et al. 2018 ). Data from observations are shown with black edge colours, coloured dots without edge colours are from post-processing analysis. Our predictions are consistent with the IR luminosity of low-redshift sources. However, our model is unable to reproduce the bolometric luminosities of the most luminous observed AGNs, falling short by several orders of magnitude. This suggests that most of the black holes at low luminosities in the early Universe are below the detection limit of current observational surveys.

panel) and ≥ 4 (right panel). The BHAR-SFR ratio scatters around a mean value of 1/10 3 , which is consistent with what Mullaney et al.

(2012) predict for main-sequence AGN at = 1. High SFRs and BHARs seem to be characteristic features of observed AGN beyond = 2, while M FIRE simulations produce relatively moderate SFRs and BHARs. Assuming our sample is a good representation of less luminous AGN and galaxies at high redshift, we predict a large number of yet unobserved low-luminosity AGN at high redshift.

Observation and theory show that the average BHAR and the average SFR correlate well (Mullaney et al. 2012; Calhau et al. 2017; Dai et al. 2018b; Volonteri et al. 2015) . While the cause of this correlation is not yet fully understood, several explanations have been proposed. The first explanation refers to a common cause. In particular, gas reservoirs in the galaxy are both the primary source for BH feeding and SF (Anglés-Alcázar et al. 2015 ). An alternative explanation is the idea that AGN activity may drive SF. For instance, the high-velocity outflows from AGN can sweep the gas away and pierce a cavity along its way but also trigger SF by induced pressure of the edges (Cresci et al. 2015) . Hence, AGN feedback may enhance SF and be responsible for its suppression at the same time Best & Heckman 2012; Ivison et al. 2012; Norris et al. 2012) . Figure 10 shows how the ratios of BHAR/SFR and BHAR/SGR scale with the total stellar mass of host galaxies in the A, B, and Cseries of M FIRE. We also investigate how the relation between BHAR and SFR (or SGR) depends on spatial scale (< gal and < 100 pc) as well as on the averaging timescale avg of the SFR and SGR (ranging from 5 Myr to 100 Myr). The description of how SFRs and SGRs are calculated is provided in section 3.1. We caution that our results are inferred by combining galaxies over a range of redshifts with lower masses being probed at higher redshift. These scaling relations may thus look different if measured for a sample of galaxies of different mass at a fixed redshift.

When employing galaxy-wide SFRs or SGRs (bottom row of Figure 10 ), we find that the slope of the * -BHAR/SFR relation predicted by our study (∼ 0.80 ± 0.22) is in good agreement with the slope (∼ 0.73[+0.22, −0.29]) inferred based on a compilation of observations of star-forming galaxies by Delvecchio et al. (2019) . However, we predict a normalization of the relation that is approximately an order of magnitude smaller than observed. Interestingly, the highest BHAR/SFR ratios that we obtain at a given stellar mass match observed values well suggesting that observational selection biases against low-luminosity AGN could be an explanation for the difference in normalization. If we aim to match the normalization by, e.g., boosting the mass retention rate from 5% to 25%, our predicted SMBHs become too massive at low redshift and thus inconsistent with the local * − BH scaling relation.

The analogous relations measured within the central 100 pc (top row of Figure 10 ) are significantly different. They have significantly larger normalization and a shallower slope. Furthermore, we find that the scatter of the BHAR/SFR and BHAR/SGR relations remains very substantial even when SFRs and SGRs are measured within the central regions of galaxies. In fact, we find that using central SFRs and SGRs result in a larger scatter than using galaxy-wide SFRs and SGRs. This latter result is perhaps somewhat unexpected. The BHAR should better correlate with the nuclear SFR rather than the total SFR, except during galaxy mergers, since the timescales of the nuclear star formation and the accretion onto SMBHs are close to the dynamical timescale of matter in the nuclear region (∼ 100 pc) Hopkins & Quataert (2010) . A similar conclusion was reached by Yang et al. (2019) . However, during mergers the total SFR of the host galaxy correlates well with BHAR because global dynamics becomes more important than the local processes in terms of angular momentum loss . Figure 10 shows that the difference in scatter between the central and the wholegalaxy BHAR/SFR relations is most noticeable at high redshift, when mergers are expected to be much more frequent (Netzer et al. 2016; Silva et al. 2021) . Hence, high merger rates at high provide a possible explanation of our result (Duncan et al. 2019 ).

The precise slope, normalization, and scatter of the relation between * and BHAR/SFR (or BHAR/SGR) in Figure 10 varies depending on the averaging timescale (the different columns of the figure show averaging times ranging from 5 to 100 Myr). First, the slope decreases when the averaging timescales are increased. Secondly, the scatter increases especially when comparing the * -BHAR/SFR relation on 100 pc scales for short and long averaging times. Overall, however, the choice of the averaging timescale plays a much smaller role than the choice of the spatial scale over which SFRs (or SGRs) are measured.

The figure also reveals that the BHAR/SGR relation at = 1−4 has a steeper slope (∼ 1.40 ± 0.25) than the BHAR/SFR relation. This finding implies that the SFR increases more strongly with increasing stellar mass than the SGR, while the opposite behavior would have been expected based on the increased galaxy merging activity in massive galaxies (Ferreras et al. 2014 (Ferreras et al. , 2016 Zahid et al. 2019 ). However, the SGR differs from the SFR not only by the additional merger contribution, but also by the decrease in stellar mass of stars already present at time − avg , i.e., those formed before − avg . This reduction in stellar mass can arise in multiple ways. First, supernovae and stellar winds return mass from the stellar component to the gas component. The contribution of this stellar mass loss to the SGR (when averaged over avg ) becomes less severe with increasing avg , in agreement with the results shown in Figure 10 , as a single stellar population loses much of its mass early on (Chabrier 2003) . Secondly, stellar mass may also be lost when stellar particles migrate outside the fixed physical radius used to calculate the SGRs. We thus conclude that the steeper scaling of the * -BHAR/SGR relation reflects the higher importance of stellar mass loss compared with galaxy mergers in our sample. Figure 11 presents how the * − BH scaling relation changes with redshift. To arrive at quantitative estimates of the slope of the scaling relation, we fit log stellar masses of all galaxies above 10 8 in our sample as well as their SMBHs at redshifts = 0, 1, 2, 4, 6 with a linear function. We also estimate the slopes for galaxies with stellar masses between 10 8 and 10 10.5 to mitigate biases due to the increasing stellar mass with redshift for the galaxies in our sample.

We fit the slope and normalization of the * − BH relation shown in Figure 11 via linear regression, see Table 3 . In the and models, the slope and normalization evolve only mildly with redshift, especially at ≤ 4. In contrast, the model shows a significant evolution of the * − BH relation with shallower slopes and lower normalizations at higher redshift. We perform the linear regression both on our complete galaxy sample and on galaxies with stellar masses between 10 8 and 10 10.5 , finding very similar results in either case. Clearly, the frequency of SMBH mergers play a critical role in shaping the * − BH relation at higher redshift. The left panel shows the * − BH scaling relation in the model. Overall, the predictions of this model are in line with the local * − BH scaling relation from Reines & Volonteri (2015) . The dotted dashed lines are the best fit lines for M FIRE data selected at certain redshifts. The slope of these fit lines increases from = 6 to = 0. The middle panel of Figure 11 show the predictions when we consider a more realistic model of SMBH merging (the model). Here, the * − BH scaling relation has a noticeably shallower slope, especially at > 2, compared to the ansatz. When SMBH mergers are not considered ( ), SMBH masses at = 2, 4, 6 are well below the local scaling relation but they still catch up to the local relation at < 1. In summary, SMBH mergers have a significant effect on the * − BH scaling relation, especially at high redshift. Hence, by accurately measuring the redshift evolution of the galaxy-SMBH scaling relation, it may be possible to constrain the rates of SMBH merging. Figure 12 contains a literature compilation of the correlation between the properties of host galaxies such as stellar (Merloni et SFRs of 0.5 < < 3 galaxies. It should be compared to the model predictions (dashed lines) in the lower panels. The BHAR normalized to SFR (or SGR) increases with the stellar mass of galaxies in M FIRE in qualitative agreement with observations. However, we predict a large number of galaxies with BHARs that are an order of magnitude lower than expected from the fit by Delvecchio et al. (2019) . Overall, slope, normalization, and amount of scatter can vary widely depending on how SFRs and SGRs are calculated.

Manne-Nicholas 2018), bulge (Savorgnan et al. 2016; Sahu et al. 2019) , and dynamical mass (Maiolino et al. 2005; Riechers et al. 2009; Venemans et al. 2012; Wang et al. 2013; Venemans et al. 2013; Kimball et al. 2015; Bañados et al. 2015; Willott et al. 2015; Bischetti et al. 2016; Wang et al. 2016; Venemans et al. 2016; Trakhtenbrot et al. 2017; Tsai et al. 2018; Decarli et al. 2018; Eilers et al. 2018; Izumi et al. 2018; Feruglio et al. 2018 ) and SMBH mass for 0 7.5. Here, bulge mass represents a lower limit for the stellar mass of the host galaxy while dynamical mass reflects an upper limit for the stellar mass. SMBH masses in this sample are determined via spectral lines for the AGN host galaxies and dynamics for the elliptical galaxies in the local Universe. Combining the various observational data we find that log BH ∼ (1.20 ± 0.06) log( gal /10 11 ) + (8.34 ± 0.04) where gal is the stellar mass for most observations, and bulge or dynamical mass for the remaining observations.

We also study how the offset of SMBH mass from this average relation, Δ log( BH ) obs , evolves with redshift. We fit the observational data in Figure 12 via linear regression finding Δ log( BH ) obs = (1.49 ± 0.10) × log(1 + ) − (0.28 ± 0.04) pointing to a super linear correlation. Measuring the offset with respect to the observed − BH or dyn − BH results in slopes which vary between 1.38 ± 0.11 to 1.57 ± 0.10. Merloni et al. (2010) re-port a positive slope of ∼ 0.68 ± 0.12 for their sample between 1 < < 2.2, as well as Ding et al. (2020) reporting 1.03 ± 0.25 for 0 < < 1.7.

The slope we obtain with the observational sample is ≈ 2.5 times larger than the slope reported by Merloni et al. (2010) . Also, the offset of SMBH mass appears to steadily increase with redshift. The difference between our re-analysis of observational data and the results by Merloni et al. (2010) are thus likely attributable to selection effect as many new sources at high redshift were discovered in the last decade.

The result of a similar analysis for the M FIRE sample is shown in Figure 13 for a variety of seed masses, centring methods and stellar mass bins. A detailed version of Figure 13 is shown in Figure A6 . We only include SMBHs with host galaxies that have stellar masses larger than 10 9 , to stay consistent with the observational data from the literature. The left panel of Figure 13 contains the best fit lines for different cases. The offset of BH from the best fit line (Δ log( BH ) pp ) is computed for postprocessing data following the same method in Figure 12 . However, our simulations show that the offset generally decreases with redshift on the right panel of Figure 13 , i.e., SMBHs at high redshifts tend to be undermassive compared to the local * − BH scaling relation, see Figure A6 . Consequently, our simulations predict Triangles and squares show data belonging to simulations that were run with FIRE-1 and FIRE-2 physics, respectively. Different colours represent different redshifts. The solid line stands for the fit function from Reines & Volonteri (2015) , while the dashed line is its extrapolation. Dotted dashed lines show linear fits to the log * −log BH at = 0−6. The fit only uses galaxies with * ≥ 10 8 to minimize the impact of the SMBH seed mass choice on the slope. The slope of the * − BH scaling relation is close to linear at = 0 (0.93±0.18, 0.96±0.18, 0.97±0.18) while the slope of Reines & Volonteri (2015) for disk galaxies is 1.05±0.11, but decreases with increasing redshift (e.g., 0.05 − 0.6 at = 6 depending on the SMBH merger model). A model with a higher number of SMBH mergers typically results in a steeper slope. The importance of mergers is particularly evident at high redshift, while the slope at the late times is not affected much, in agreement with our previous finding ( Figure 4 ) that SMBH mass at low redshift is primarily set by gas accretion and not SMBH mergers. a large number of low-luminosity AGN in high redshift galaxies. The result holds for all different post-processing models. Furthermore, the most negative slopes are obtained in low mass galaxies ( * < 10 10 ) when SMBHs are placed at the densest regions (MAX) and in massive galaxies ( * > 10 10 ) when a COM placing is used. The difference between our simulation results and observational data may be explained by selection effects that bias the samples of observed high redshift AGN. Efforts in the search for the low-luminosity AGNs in the early Universe (i.e. Subaru High-z Exploration of Low-luminosity Quasars Project, Matsuoka et al. 2016 ) are therefore paramount to better constrain the redshift evolution of Δ log( BH ) .

The slope of the * − BH scaling relation contains essential information about the growth trends of SMBHs and their host galaxies. In particular, it is set by the interplay of the stellar mass, the black hole mass, the black hole accretion rate, and the stellar growth rate, see Equation 7. In this section we will discuss the slope of the trajectory of black holes evolving in the * − BH plane. We note that this slope is not strictly identical to the slope of the * − BH relation of a population of galaxies with different stellar masses at a fixed redshift. However, for the latter we find only a small amount of redshift dependence at < 4 implying that the slope of trajectories in the * − BH plane will be similar to the The first and second columns shows the mergers model used in the post-processing analysis and the redshift. The third and the fourth columns show slope and normalization of the linear fit when we include only galaxies with 8 ≥ log * / < 10.5 from our sample. The lower table lists the slope and normalization that we obtain when we include all galaxies with * ≥ 10 8 . SMBH mergers have a significant impact on the slope and normalization of the * − BH scaling relation, especially at ≥ 2. In the and models the slope and normalization evolve much more gradually with redshift (with only mild changes since = 4) than in the model.

slope of a population of galaxies and SMBHs at fixed redshift for < 4.

According to Equation 7, the slope is larger (smaller) than unity when the sBHAR exceeds (is less than) the sSGR. The slope goes to zero when the black hole ceases to grow, and it reaches infinity for a growing SMBH in a non-growing galaxy. When the specific growth rates of the SMBH equals that of its host galaxy, the slope is unity.

Panels A and B of Figure 14 show how the trajectory and the slope of the most massive galaxy in the simulation A1 evolves in the model and with MAX centring. Here, we calculate the slope based on the (smoothed) trajectory of the galaxy and Red, green and blue data points indicate, respectively, bulge mass, stellar mass , and dynamical mass for both panels. Bulge masses put a lower limit to the stellar mass of the galaxies while dynamical mass is the upper limit for the galaxy stellar mass. Different lines correspond to the best fit lines taken from literature, HR04 ( Ding et al. (2020) found: 2 = 0.68 ± 0.12 for 1.1 < < 2.2 and 2 = 1.03 ± 0.25 for 0 < < 1.7, respectively. its SMBH in the * − BH plane. More precisely, the smoothed trajectory is calculated with the help of sliding bins in log stellar mass of width 0.1 and with subsequent shifts of 0.01. The slope is then calculated from the smoothed trajectory using a 0.1 dex in stellar mass.

Subsequently, we analyze the various terms in Equation 7 to see how they affect the slope. According to the panel C of Figure 14 , the sBHAR shows significant variations but not a strong evolutionary trend, except for a moderate decrease with increasing stellar mass when * 10 10. 3 . Overall, the sBHAR mostly lies between 10 −1 and 1 Gyr −1 . In contrast, the sSGR in panel D decreases steadily with time (and stellar mass). It starts at 10 1.5 Gyr −1 and reaches 10 −1 Gyr −1 when the galaxy becomes massive. Hence, at early times (when the galaxy mass is low), the sBHAR is often much lower than the sSGR resulting in a sub-unity slope. By the time the galaxy reaches a stellar mass of * ∼ 10 10. 3 , the sSGR has decreased sufficiently such that the sBHAR is now larger than the sSGR and the slope becomes very large. Subsequently, the sBHAR decreases to a similar level as the sSGR and the slope reaches unity.

Another perspective can be gained by comparing the sSGR and sBHAR to the inverse Hubble time −1 Hubble in panel E. When the galaxy has a relatively low stellar mass, the sSGR exceeds the inverse Hubble time indicating a quickly growing galaxy. In contrast, the sBHAR typically falls below −1 Hubble during this time, indicating slow SMBH growth. However, as soon as the galaxy reaches * ∼ 10 10.3 , the sBHAR approaches and then exceeds −1 Hubble . Hence, when the galaxy becomes massive, the SMBH grows as fast as or even faster than its host galaxy on a inverse timescale similar to −1 Hubble .

Equation 7 allows us to further understand the slope of the * − BH relation via the BHAR, SGR, and the * / BH ratio. In particular, the panel F of Figure 14 shows that the * / BH ratio does not change enough to affect the trend of the slope over much of the history of this galaxy. A change in * / BH is thus clearly not driving the slope of the * − BH relation. Instead, the slope is set by the BHAR to SGR ratio. Whenever this ratio exceeds the * / BH ratio, the slope becomes large, while for small values of the BHAR to SGR ratio (i.e., when BHAR < 10 −4 SGR), the slope is below unity.

Finally, Figure 14 also offers insights into which of the two terms, BHAR and SGR, plays a more important role in setting the slope of the * − BH relation. The slope is low at early times because the BHAR is there much lower than the SGR. A higher slope would require either faster SMBH growth or slower galaxy growth. Subsequently, when the galaxy grows its stellar mass from * ∼ 10 10.1 to * ∼ 10 10.4 , the BHAR increases and the SGR decreases. This combination results in a BHAR to SGR ratio that finally exceeds ∼ 10 −4 and thus a high value of the slope. Subsequently, the SGR slightly increases again while the BHAR remains nearly constant resulting in BHAR ≈ SGR and thus an approximately linear slope. Furthermore, a slope change at a later time is well correlated with a change in the BHAR while the SGR is approximately constant. We therefore conclude that both the BHAR and the SGR contribute in a significant manner to the slope evolution of the * − BH relation. Best fit lines of the * − BH relation in the model for various seed masses (10 2 in green, 10 3 in red, and 10 4 in blue) and centring methods (solid lines for MAX, dashed lines for COM). We calculate Δ log(M BH ) pp (see text) as the offset from these best fit lines. Right Panel: Change in slope ( 2 ) for different stellar mass bins in the " " model. Here, 2 represents the slope of the best fit lines in log(1 + z) vs Δ log(M BH ) pp , see Figure A6 . Different colors refer to different seed masses (see legend). Dots and triangles show results for MAX and COM centring methods, respectively. The error bars indicate half the size of the stellar mass bins, while error bars show the fit error of the slope, 2 . The slope is generally negative and tends to increase with increasing stellar mass, except in the COM centring model with seed ≥ 10 3 .

In this section, we introduce two simple models, a one-and a two-zone model, for the growth of SMBHs based on the stellar growth history of the host galaxy to offer additional insight into the origin of the * − BH scaling relation in the context of the GTDA model. In both models, the stellar mass of the galaxy is given by the integral of their star formation rates. However, the models differ in how they calculate the black hole accretion rate. In the one-zone model, the central mass which determines the accretion rate is set to a fixed fraction (5%) of the total stellar mass of the galaxy. In contrast, in the two-zone model the central to total stellar mass ratio is allowed to vary. We start by making the following assumptions:

• First, we assume that the total mass within 0 is dominated by the stellar mass, tot (< 0 ) ∼ * (< 0 ). This assumption holds in our simulations, see Figure A4 , as the stellar-to-total ratio is close to unity for 6. At higher redshift, this basic assumption may break down. We refer the reader to §5 for the caveats when modelling SMBH growth via the GTDA model in the early Universe.

• Secondly, we split Equation 2 into the term tot (< 0 ) ∼ * (< 0 ) and collect all other dependencies into a time-dependent function . In this approximation, the BHAR becomes proportional to the total stellar mass inside 0 . This ansatz allows us to directly tie the slope to the stellar growth history of the host galaxy.

• Thirdly, we ignore the merger contribution to SMBH growth.

With the assumptions listed above, we can link the growth of SMBHs to the growth of their host galaxies. The BHAR scales pseudo-linearly with the stellar mass of a galaxy,

where ( ) encapsulates the rather complex dependencies of Equation 2. Finally, we assume that can be modeled by a parametric function of the form

to simplify our analytic calculations. The values x 1 = 0.052 M yr −1 , x 2 = 0.002 M yr −1 , t * = 0.35 Gyr, and = 0.23 Gyr result in ( ) that is (at ≥ 1) in good agreement with our full post-processing analysis, see Figure 15 .

The SGR is simply expressed as the change in the stellar mass within (Δ * (< )/Δ ) between two adjacent snapshots in the post-processing analysis. Here, refers to either the radius of the central region 0 or to the size of the galaxy gal . The stellar mass ( * (< )) is obtained by integrating a function fitted (see Equation 10) to the central SGR data from the post-processing analysis.

SGR(< R) = a 1 tanh t − t * + a 2 (10) The stellar mass in fixed radius is connected to the stellar growth rate within the same radius via the integral:

Finally, the black hole mass is given as

with seed = 10 4 placed at = 20 ( 0 ∼ 180 Myr) as the fiducial case.

The one-zone model makes the simplifying assumption that * (< 0 ) ∝ * (< gal ). As a result, the BHAR is directly linked to the stellar mass of galaxies. We adopt a pre-factor of 5% to approximately match the observed normalization for the * − BH relation and to bring the predictions of the one-zone model in better agreement with the results of the two-zone model discussed in the next section. First, we demonstrate how different analytical SGR histories affect the evolution of the slope of the * − relation. For instance, let us assume a constant SGR evolution in the form of SGR = c 1 , where the unit of 1 is mass over time. Then, the host galaxy stellar mass becomes M * = c 1 t, and

By inserting Equation 13 into Equation 12, we arrive at the following estimate for the SMBH mass M BH = M seed + 0.05

x 1 c 1 10 9 M 2 Ω 1 (t) + 0.05 x 2 c 1 10 9 M t 2 2

A constant SGR thus yields a slope of

In particular, the slope is 0 when the black hole seed mass is much larger than the accretion contribution ( seed ∫ BHAR dt ). At early times the slope is thus rather shallow unless the seed mass is sufficiently small, see Figure 7 . Instead, when the accretion contribution matches surpasses the seed mass, the slope approaches ∼ 2 for a constant SGR. The slope takes various values for 0 < ∞, and sometimes the same value more than once. Therefore eq. 15 explains the slow-to-fast transition of SMBH growth seen, e.g., in Figure 7 . At early times the SMBH mass is dominated by its seed mass. Hence, the slope is shallow. At late times, the seed mass is small compared to the accreted mass and the slope is thus steep.

Another simple scenario is a galaxy in which the SGR increases linearly with time, SGR = 2 where the unit of 2 is mass over time squared, resulting in slope = c 2 t 2 2

An increasing SGR lowers the slope compared with the case of a constant SGR. However, the SMBH can still grow quickly compared to its host galaxy whenever the accretion contribution exceeds the black hole seed mass.

Finally, a galaxy may experience an epoch in which the SGR decreases with time, e.g., SGR( ) = SGR( 1 ) − 3 ( − 1 ) for ≥ 0, where the unit of 3 is inverse time. In this case,

.

(17) Figure 15 shows the predictions for the SMBH growth when applying the one-zone model to the aggregated stellar growth histories of the A, B, and C-series of M FIRE, see Table 1 . In agreement with the simple analytical examples discussed above, the one-zone model predicts a shallow-to-steep transition of the slope. The transition takes place between = 4 and = 2, when the galaxy transitions from being low mass to becoming a massive galaxy (near * < 10 10 ). Interestingly, at this time the SMBH mass already exceeds the seed mass by an order of magnitude, i.e., we are already in the seed ∫ 0 regime discussed in, e.g., eq:constant SGR one-zone model. We also confirmed that reducing the seed mass to, e.g., seed < 10 2 does not change the location of the slope transition, again indicating that this increase in slope is not a seed mass effect. Instead, Figure 15 shows that during the time of the transition ( = 2 − 4) the sBHAR declines much slower than the sSGR, resulting in an increasing slope. Hence, the increase in the slope of the * − BH relation in ∼ * < 10 10 galaxies is a natural consequence of the specific evolutionary history of the SGRs of such galaxies. Ultimately, their SGRs are set by the complex interplay of gas inflows, feedback, and galaxy mergers. Figure 7 shows that the * − BH relation flattens at late times in very massive galaxies. While this finding is only based on a small number of simulated galaxies at relatively low redshift, and thus tentative, we would like to explore the physical origin of this finding with the help of the one-zone model. To this end, we study 3 possible scenarios. In our first scenario, we assume that the SGR decreases linearly with increasing log stellar mass by half a dex between = 1 and = 0 while is kept constant at its = 1 value, ∼ 10 −3 −1 . In the second scenario, the SGR stays constant until = 0 but the ( ) parameter decreases by half a dex. The third scenario is identical to the first, except that the SGR increases by half a dex between = 1 and = 0.

For each of the three scenarios, we first fit the average behaviour of the galaxy-wide SGR and of ( ) until = 1 and then extrapolate to = 0. Next, we integrate the SGR to obtain the stellar mass and the SMBH mass as function of time as described in §4.4.1, see in particular Equation 8.

In the first scenario, the SGR decreases at low which results in a * that stays nearly constant after = 1. However, the average BHAR continues to grow since the BHAR is assumed to be propor- z Figure 15 . Results from the one-zone model. The first row show the average SGR measured from the A, B, and C-series of M FIRE at ≥ 1 (dark blue) and for the two simulations that continue below = 1 (light blue). The solid line shows 3 scenarios of how the > 1 SGR can be extended to low z. (From left-to-right) the SGR decreases by half a dex between = 1 and = 0, remains constant at the value at = 1, or increases by half a dex. The second row shows the average stellar mass (red) while the third row shows the average ( ) parameter, see Equation 2, which normalizes the pseudo-linear dependence of the BHAR on * . The ( ) parameter is assumed to remain constant at the = 1 value in the left and right columns but it decreases by half a dex between = 1 and = 0 in the middle column. The fourth row show the sSGR, the sBHAR while the fifth panel shows the slope of the trajectory in * − BH space derived from Equation 7. The trajectory in * − BH space is shown in the bottom row as a colored line indicating the redshift. These panels also show the * − BH scaling relation from the post-processing analysis as an histogram. The shaded region shows the result from post-processing analysis using fiducial settings, see Figure 7 . The one-zone model describes the overall evolution of simulated galaxies in the * − BH space well, despite its high degree of simplification compared to the full GTDA model. tional to the stellar mass. Consequently, the slope of the trajectory in * − BH space gradually steepens between = 1 and = 0 in line with Equation 17. This trend disagrees however with the results shown in Figure 7 leading us to exclude this scenario.

In second scenario the SGR is kept constant between = 1 and = 0, while ( ) decreases. In this case, the average * grows more rapidly, while the decreasing slows the growth of the SMBH resulting in a nearly constant BHAR at late times. ( ) is strongly dependent on the disc fraction ( ( ) ∝

) in the GTDA model and we expect a decrease in disc fraction at late times as massive galaxies transition from disc to early type morphology.

The final scenario assumes an increasing SGR for < 1. Subsequently, the stellar mass grows substantially (by about one dex) between = 1 and = 0. Even though the larger stellar mass also boosts the BHAR, the fast galaxy growth results in a slope that is only mildly super-linear and approaches unity at late times. Physically, an increasing SGR at low redshift may arise from the late assembly of massive galaxies via merging expected in a hierarchical Universe.

In contrast to the one-zone model, the two-zone model uses the stellar mass within 0 , and not within the galaxy, to calculate the BHAR. Furthermore, the link between BHAR and stellar mass (as opposed to SFR) within the central region as given by equation 2 has also observational support. Yang et al. (2019) find that the BHAR correlates better with * rather than SFR in the central region in non-bulge dominated galaxies which applies to the majority of the galaxies in our sample, see Figure A8 .

In the one-zone model, the slope of the * − BH relation depends mainly on the total SGR of a given galaxy. Here, however, the slope is affected both by the galaxy-wide SGR (which drives * ) and by the central SGR (which drives the BHAR). We can thus investigate how the slope changes for various options of a constant or increasing SGR, either galaxy-wide or in the central region. The slope is provided in Table 4 in the limits of seed exceeding, is equal to, or smaller than ∫ 0 BHAR . Similar to §4.4.1, we can apply also the two-zone model to the SGR histories of M FIRE galaxies. Given the lack of many galaxies at low , we again explore three possibilities for the SGR evolution at < 1. We will show that while these scenarios differ, they predict similar results for the slope of * − BH scaling relation in the local Universe. The three scenarios are as follows:

• The central SGR decreases by half a dex since = 1 while ( ) and the total SGR remain constant.

• ( ) decreases with time at < 1 while the central and total SGR remain constant.

• The central SGR and ( ) are constant while the total SGR increases by half a dex since = 1.

A decrease in the central SGR, as speculated in the first scenario, could originate, e.g., in inside-out quenching due to the AGN feedback (Tacchella et al. 2015; Ellison et al. 2018; Abdurro'uf 2018; Tacchella et al. 2018) . Strong stellar feedback (Cox et al. 2006) , merger quenching (Gabor et al. 2010) , or gravitational heating due to clumpy accretion (Birnboim et al. 2007; Dekel & Birnboim 2008; Dekel et al. 2009 ) could also reduce the central SGR by lowering the central star formation activity. The decrease in ( ), proposed in the second scenario, is expected from a disc-to-early type morphological transformation, while the increase in the total SGR (but not central SGR) could arise from galaxy merging.

As Figure 16 shows, all three scenarios result in a * − BH relation that is consistent with the aggregated M FIRE results even at < 1. In particular, it shows the shallow-to-steep transition of the slope when galaxies approach * ∼ 10 10

. In addition, all three scenarios show a reduction in slope in massive galaxies at late times. Scenario 2 and 3 leads to final slopes of order unity, while the first scenario lowers the slope to about 1.5.

We expect that all three scenarios are partly at play in the real Universe. The * − BH relation at late times may thus be especially susceptible to the differential growth of galaxies in their central region and on galaxy-wide scales. In addition, similar to the one-zone model, galaxy merging or a change in galaxy morphology can strongly affect the * − BH slope. Figure 17 compares the predictions of one-zone and twozone models. A main difference is the somewhat more pronounced change in the slope of the * − BH relation at ∼ 1 − 4 Gyr in the two-zone model compared to the one-zone model. However, overall the predictions are rather similar. The right panel of Figure 17 shows the BHAR-SGR ratio for both the one-zone and the two-zone model. Again, both models make rather similar predictions for the * − BH relation and for the stellar mass dependence of the ratio between the BHAR and the SGR. Overall, the two-zone model results in a slightly more accurate representation of the full post-processing analysis, see Figure 15 and Figure 16 , specifically in more pronounced changes from the shallow to the steep slope regime as well as from the steep to the approximately linear slope regime.

As we saw in Figure 10 , the SGR and the BHAR trace each other reasonably well both in observations and simulations albeit with significant scatter. In the context of the GTDA model, this link is facilitated via the amount of mass (stars and gas) in the central region of a galaxy which is directly driving the BHAR. In the previous sections, we discussed the consequences of simplified toy models based on this general idea. Specifically, in the one-zone model (two-zone model) the galaxy-wide stellar mass (the central stellar mass) is assumed to be proportional to the baryonic mass in the central region. Hence, the BHAR and thus SMBH growth can be calculated once the SGR history is known.

We have shown that there are in general (at least) 3 distinct z Figure 16 . Results of the two-zone model, similar to the one-zone model. We consider three different scenarios (columns left or right) for the late time evolution of SGR and . The first column considers a scenario in which central SGR decreases while the total SGR and ( ) remains constant. In the second column, central and total SGRs are constant, and ( ) decreases by half a dex after z = 1. The final column studies the impact of growing total SGR while central SGR and constant. The first row shows the average points of central (red) and total (blue) SGR of the A, B, and C-series of M FIRE where log( halo ) = 12.5 at = 2 and best fit lines to the average central and total SGR for 1. The light blue and light red colours represent the averaged data from A1 and A4 simulations with FIRE-1 physics that we do not include in the best fit. The solid red line in the second row shows the central stellar mass (integration of the red line in the first row) denoted as M * ,in . The solid black line in the third row shows BHAR/M * ,in in the MAX centring case. The solid blue line is the integration of the total SGR, blue line in the first row. In the fourth row, we show the specific growth rates of the galaxies (sSGR = SGR/M * ) and SMBHs (sBHAR = BHAR/M BH ). The fifth row consists of the slope of the scaling relation (sBHAR/sSGR) as a solid black line. The black dashed and dotted dashed lines are the slope for different seed masses, 10 2 and 10 5 . Further reduction in central SGR and ( ) give flatter slopes, and a boost in the total stellar growth also has a similar effect on the slope of the * − BH scaling relation. The final panel at the bottom shows the prediction of our two-zone model for the * − BH scaling relation colour-coded by redshift. The shaded region shows the result from post-processing analysis using fiducial settings, see Figure 7 . Compared to the one-zone model, the two-zone model predicts a more noticeable shallow-to-steep transition of the slope of the * − relation for log( * / ) = 10 galaxies. Reines & Volonteri (2015) . The red and blue solid lines are the predictions of our one-zone and two-zone models, where the BHAR is linked to the stellar growth history. The slopes predicted by the one-and two-zone models flatten at high masses log( * / ) ∼ 11 − 11.5 as in the transition from LTGs to ETGs in Sahu et al. (2019) . Right Panel: Predictions for how the BHAR-SGR ratio changes with the total stellar mass. This panel is the same as Figure 10 , but compares the BHAR-SGR ratio from the toy models with the observational result of Aird et al. (2019) , instead of being based on the post-processing data from M FIRE. The solid red line indicates the BHAR-SGR ratio from the one-zone model and the solid blue line represents the ratio between the BHAR and the total SGR in two-zone model. epochs as galaxies move through * − BH space. First, galaxies grow quickly while the SMBH mass does not (shallow slope). This epoch lasts until the total stellar mass reaches about 10 10 . This threshold is generally reached during cosmic noon ( ∼ 2 − 3) for our M FIRE sample. Secondly, the slope increases steadily marking efficient SMBH growth (steep slope). Finally, our analysis predicts that under certain assumptions the slope decreases again at late times (approximately linear slope).

During early times, galaxies have relatively low masses, form stars at high rates (Riechers et al. 2013; Finkelstein et al. 2013; Casey et al. 2014; Zavala et al. 2018; Bowler et al. 2018; Berta et al. 2020 ) and thus increase their stellar masses quickly. Specifically, the shallower potential wells in the centers of galaxies at earlier times may boost SN-driven mass ejections (Dubois et al. 2015) which plays a crucial role in suppressing the accretion onto SMBH. Also, stellar feedback may drive buoyant outflows of high-entropy gas from the central regions (Bower et al. 2017b ) which ceases to be effective in massive halos with hot gas coronas. In addition, a virialization of the circum-galactic medium down to the central galaxy ('inner CGM virialization', Stern et al. 2020) , which can stabilize disks against feedback driven outflows, has typically not yet taken place. Consequently, galaxies grow much quicker than their SMBHs at those early times resulting in a shallow slope for the trajectory in * − BH space.

As a galaxy becomes moderately massive, ∼ 10 10 , the escape speed of the galaxy exceeds the characteristic speed of SNdriven winds thus strongly reducing the effectiveness of galactic outflows (Anglés-Alcázar et al. 2017c) . In fact, a stellar compactness of 10 10 kpc −1 in the central region would be enough to keep the SN-driven winds within the host galaxy centre (Dubois et al. 2015) . During this time, SMBHs will be able to grow quickly resulting in a steep slope because (i) galaxies may contain significant reservoirs of gas in their centers and (ii) galaxies of this mass have often a major disc component allowing gravitational torques to operate efficiently (Querejeta et al. 2016; Anglés-Alcázar et al. 2017a; Blumenthal & Barnes 2018; Thomas et al. 2019 ).

Finally, we find that when galaxies become very massive, ∼ 10 11 , the masses of galaxies and of their SMBHs often grow at similar rates. The one-and two-zone models presented in the previous sections provide some insights into the origin of this behavior. Specifically, we pointed to two possibilities consistent with our analysis. First, the close to linear slope could originate in a reduction in the ( ) term in the GTDA model. ( ) depends super-linearly on the disc fraction, ( ) ∝ 5/2 disc . Hence, a transition from disc to elliptical morphology in the central regions of galaxies reduces the SMBH accretion rate. Such a transition is expected given the change in overall Hubble type when galaxies grow in mass (D'Onofrio et al. 2015; Tacchella et al. 2019; Cooke et al. 2019) . Supporting this scenario is also the observational finding that the slope of the * − BH relation depends on galaxy morphology Davis et al. (2019) . Secondly, the slope could flatten not because the SMBH accretion rate decreases, but because the SGR increases, e.g., due to a larger number of galaxy mergers in massive galaxies (Marchesini et al. 2014; Bellstedt et al. 2016; Buchan & Shankar 2016; Vulcani et al. 2016; Nipoti et al. 2018 ).

The results of this study are subject to a few potential caveats. First, there is the potential question whether the GTDA model is applicable to galaxies over a large range of masses and redshift. The GTDA model was developed for disc galaxies and estimates the accretion from circum-nuclear (∼ 100 pc) to sub-pc scales using the properties of the circum-nuclear region, including its stellar mass, baryonic disc fraction, and gas mass. However, especially at high redshift, the host galaxies may not always be disc galaxies (Cowie et al. 1995; van den Bergh et al. 1996; Tacconi et al. 2010; Genzel et al. 2011; Guo et al. 2012; Zanella et al. 2015 ) even though massive disc galaxies have certainly been found even at = 4 and beyond (Hodge et al. 2012; Neeleman et al. 2020) . Fortunately, the specific details of the GTDA model appear to be less important given that a very similar growth history for SMBHS can be obtained by using the much simpler dynamical accretion model (Equation 6) with = 10 −4 . Hence, our predictionns are likely robust as long as the accretion model results in SMBHs accreting only a small fraction (∼ 0.1%) of the available gas per free-fall time (Anglés-Alcázar et al. 2017b,c) . Models that may contribute to a faster, more efficient gas accretion onto SMBHs, such as chaotic accretion of hot gas (Thomas et al. 2019; Davé et al. 2019) or merger-triggered accretion (Capelo & Dotti 2017; Ricarte & Natarajan 2018) will be left to future work.

A second potential concern is that none of the SMBHs in our post-processing analysis are as massive as the most luminous AGN observed at 6 (Mortlock et al. 2011; Wu et al. 2015; Mazzucchelli et al. 2017; Bañados et al. 2018; Yang et al. 2020; Wang et al. 2021) . While some of our haloes are sufficiently massive enough ( 10 12 ) to potentially host a very luminous AGN, simple number density and clustering arguments show that our simulation volume is likely too small to contain even a single luminous AGN. In particular, the number density of luminous AGNs at ∼ 6 − 7 is ∼ 1 cGpc −3 for 1450 < −26 (Wang et al. 2019) . While some of our simulations are run in boxes with sizes of ∼ 0.8 cGpc, we only simulate a small number of (massive) halos selected from those boxes via the zoom-in approach. Hence, the chance of selecting the halo of even a single luminous AGN is very small. The most luminous AGNs may have an atypical formation path which leads to a larger seed mass, e.g. direct collapse black holes (Bromm & Loeb 2003; Volonteri 2010) . These high-mass SMBH seeds can grow quickly when the accretion rate is strongly dependent on the SMBH mass. Furthermore, the duty cycle of luminous AGN is close to the unity at 6 (Shankar et al. 2010b (Shankar et al. ,a, 2019 implying that analyzing snapshots at different times does not substantially increase the odds of reproducing a luminous AGN in our simulations.

Thirdly, a significant simplification of our model is the treatment of SMBH mergers. To robustly asses the impact of SMBH merging on our result we have considered two extreme scenarios in addition to our fiducial ( ) case. In the first of these extreme scenarios, no SMBH mergers take place. In the second extreme case, SMBH of galaxies merge as soon as their parent halos become sub-halos of each other. While the predictions for the low Universe are shown to be rather robust to the specifics of the SMBH merger model, the SMBH-galaxy scaling relation at early times are sensitive to details of SMBH merging, thus highlighting the importance of properly accounting for SMBH mergers especially in the young Universe (Ma et al. 2021 ).

Finally, most of the M FIRE simulations used in this paper do not include black hole physics on-the-fly, especially AGN feedback. This is by design and allows us to study SMBH growth and scaling relations in the absence of AGN feedback. Thus, this work provides a basis for the future comparison with the FIRE simulations including black hole physics.

Importantly, as we showed in Figure 7 , accretion models that are weakly dependent on the SMBH mass can reproduce the * − BH local scaling relation without AGN feedback (Anglés-Alcázar et al. 2013) . We speculate that the negative effect of AGN feedback should result in overall slower growth of SMBHs. If true, the current analysis provides an upper limit on how fast SMBHs can grow in the context of the GTDA model. We leave the study of BH growth in simulations with AGN feedback for the future.

We have carried out a post-processing analysis of 34 high resolution cosmological zoom-in simulations from the M FIRE suite to study the growth of SMBHs and their host galaxies across cosmic history in the absence of AGN feedback. In particular, we have analyzed the effect of SMBH placements, SMBH merger treatments, and specifics of the accretion models on the * − BH scaling relation, focusing in particular on the gravitational torque driven accretion model (GTDA) by Hopkins & Quataert (2011) and on the evolution at 2. Our main findings are as follows.

• The masses of galaxies and their central SMBHs co-evolve, even in the absence of AGN feedback, in the GTDA model approximately in line with the local * − BH scaling relation, see Figure 7 .

• While overall in line with the local * − BH scaling relation, we find clear evidence of a significant deviation from a simple power-law relationship (a "shallow-to-steep" transition of the slope) in low to moderately-massive galaxies. The strength of this deviation depends on the specific modeling assumptions. In particular, it is more pronounced if early SMBH growth is stunted by placing them on more typical (i.e., not the most gas rich) regions near the centers of galaxies and if SMBHs are not allowed to merge. SMBHs mergers and efficient early growth of SMBHs significantly reduces this deviation from the local * − BH scaling relation.

• Model assumptions, especially, about SMBH placements and mergers leave a clear imprint on the * − BH scaling relation at high redshift in the absence of AGN feedback (see Figure 11 ). Hence, we expect a link between the SMBH merger rate and their mass ratios with any deviations of the * − BH scaling relation at high from those of local galaxies.

• Different SMBH placement and merger models have no apparent efffect on the final SMBH mass at low redshift in the context of the GTDA model. The masses of SMBHs at late times are also largely independent of the BH seed mass.

• Aside from the GTDA model, we also study alternative accretion models. Figure 5 shows that large-scale accretion models can be divided into two major groups depending on whether the accretion rate scales superlinear (e.g., Bondi-like models) or sub-linear (e.g., the GTDA model) with the SMBH mass. The first class of models results typically in over-massive SMBHs strongly indicating the need for AGN feedback. The second class of models, however, is able to reproduce the local scaling relations without the inclusion of AGN feedback.

• Currently, none of the SMBHs in our post-processing analysis are as luminous as the billion solar mass SMBHs in the early Universe likely due to the limited volume probed by our simulations. However, we have considerable overlap in the IR luminosity.

• This study predicts a large number of low luminosity AGN at high redshift which may be potentially observable with JWST. Mergers between these SMBHs may be detectable by gravitational wave experiments.

• The offset of the SMBH mass from the local * − BH scaling relation, Δ log( BH ) obs , increases towards higher redshift. While our finding differs from observational data taken at face value, such a comparison does not account for observational selection biases. We thus predict that the discovery of dimmer AGNs at high redshift could decrease the slope of − Δ log( BH ) obs relation.

• We develop two variants of an analytical model in §4.4 that link the growth of SMBHs to the stellar growth history of the host galaxies within the frame of the GTDA model. These models capture the * − BH trajectory predicted by the full post-processing analysis remarkably well thus allowing us to understand the shape and normalization of the * − BH relation in terms of the stellar growth history of galaxies.

A high merging efficiency results in a close to linear slope of the * − BH scaling relation for all stellar masses, see Figure 13 . In contrast, the * − BH shows a clear non-linear scaling if BH mergers are rare. Consequently, the slope of the * − BH relation and the merger rate of SMBHs appear intricately linked. This link may be explored observationally by constraining the mass distribution of SMBHs residing in moderately massive galaxies and by measuring BH merger rates via gravitational wave signals.

We perform a bulge-disc decomposition following (Anglés-Alcázar et al. 2014), which is explained in §4.1 in detail. In Figure A1 , we compare what we found based on the recipe of Anglés-Alcázar et al. (2014) to the velocity dispersion of all star particles within different radii; 1 kpc, the half stellar mass radius, and galactic radius gal . The black hole masses are calculated using our fiducial settings and " " model. The different radii we used to estimate the velocity dispersion of star particles gives indifferent results from the bulge-disc decomposition method we used.

In Figure 4 , we show the contribution of both accretion and seed mass to the total SMBH mass for a seed of 10 4 . Figure A2 searches for the effect of different seed masses on the total SMBH mass for FIRE-2 A series simulations. As the seed mass increases, the offset between total SMBH mass and in-situ accretion contribution increases, especially at high redshifts. The effect of seed mass choice is more visible at 2 for seed masses equal to or above 10 4 , which points out that the effect of small seeds on SMBH mass is not observable. Figure A3 shows the difference between the central coordinates of two different centring methods available in AHF. The difference between MAX and COM centring methods are generally above half of a kpc in physical units for 4.5 for a selected simulation. The coordinates of maximum density centre and centre-of-mass of the host halo are roughly the same for the rest of the analysis. This finding suggests that the central galaxy settles at the host halo centre for the sample simulation around ∼ 4.5

The gas-to-total ( gas ) and stellar-to-total ( star ) mass ratios for different halo centring methods within 0 are shown in Figure A4 . Here, the total mass is the summation of the gas, stellar and dark matter mass of the galaxy. The gas and star are the ratios of gas mass and stellar mass to the total mass, respectively. The gas and stellar fractions within 100 pc track each other for different centring methods until ∼ 4.5. The difference between MAX and COM within 1 kpc track each more consistently compared to 100 pc case. For GTDA model, it is conceivable to assume that the total stellar mass within 0 is the dictating term in Equation 2 to determine the BHAR since the ratio of stellar and total mass within 0 is close to unity for 6 − 7. Moreover, a high gas fraction at high redshift is the one of the caveats discussed in §5 for the modelling of SMBH growth using GTDA model in the early Universe. Figure A5 produces the * − BH scaling relation for different seed masses, halo centring methods, and merger treatments. The effect of full merger treatment (" ") is distinguishable again for seed masses that are equal to or above 10 4 . The increase in seed mass boosts the importance of mergers on the * − BH scaling relation. The smaller seed masses that were born in gas-rich environments follow local scaling relation since the growth of the black hole is dominated by in-situ accretion. Figure A6 shows the different slopes following the same analysis method in Figure 12 for M FIRE data in post-processing analysis using fiducial settings. Each panel in Figure A6 represents the blue and the red data points in the right panel of Figure 13 . The larger seed mass and COM centring choices flattens the slope. /h. Figure A8 shows the bulge mass to total stellar mass ratio of the sample regarding the Figure 7 . In general, the galaxies in our sample are marginally disk dominated, even at high redshift. This finding makes the GTDA model suitable for modelling the SMBH growth in M FIRE galaxies. This paper has been typeset from a T E X/L A T E X file prepared by the author. . Gas-to-total ( gas ) and stellar-to-total ( star ) ratios for different halo centring methods within 100 pc and 1 kpc central regions for FIRE-2 A1 simulation. The solid black line shows the gas fraction for the densest central region, and the red dashed line shows the same quantities for the average density central region. The blue and cyan lines stand for stellar-to-total ratio for MAX and COM centring methods, respectively. Figure 7 . The red solid line shows the best fit line for the spiral galaxies in Reines & Volonteri (2015) . The effect of gas-rich early region on the early evolution of SMBH growth decreases as seed mass increases, especially for seeds heavier than 10 4 . . Redshift evolution of the offset of the * − BH scaling relation from our best fit line for different seed mass choices, centring methods, and the model in different stellar mass bins. We follow the same method as in Figure 12 and define the offset as the difference between the best fit line to the whole sample and the post-processing data. There is a strong negative correlation between the offset from our best fit line at high stellar masses. Seed mass and the slope are inversely proportional. The maximum density centring results in a stepper negative slope for log( * / ) < 10.0. The opposite is true for stellar masses greater than 10 10 . Finally, the slope tends to flatten above * = 10 10.5 . . Bulge to total stellar mass of the galaxies in our sample in Figure 7 with " " model, SMBH seed mass of 10 4 and MAX centres. Colorbar shows the number of galaxies and black dots with error bars represent the mean value of each bin. Most of the galaxies in the sample reach = 2 and some reach = 1 while only galaxies from two simulations reach = 0. Galaxies in our sample are marginally disk dominated (with a bulge to total ratio of 0.3-0.5).

Active Galactic Nuclei 12: A Multi-Messenger Perspective (AGN12)

Frontiers in Astronomy and Space Sciences

OÇ thanks Pedro R. Capelo, Alexander P. Hobbs and Mehmet Hakan Erkut for the discussion and their valuable feedback on the manuscript. OÇ also thanks his wife MÖÇ and two cats LÇ and KÇ for their tireless efforts to create a motivating working environment during COVID-19 lockdown. RF acknowledges financial support from the Swiss National Science Foundation (grant no 157591 and 194814). DAA acknowledges support by NSF grant AST-2009687 and by the Flatiron Institute, which is supported by the Simons Foundation. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grant 17-ATP17-0067; and by a Cottrell Scholar Award and Scialog Award #26968 from the Research Corporation for Science Advancement.Simulations were run with resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research centre, proposal SMD-14-5492. Additional computing support was provided by HEC allocations SMD-14-5189, SMD-15-5950, by NSF XSEDE allocations AST120025, AST150045, AST160048, by allocations s697, s698 at the Swiss National Supercomputing Centre (CSCS), and by S3IT resources at the University of Zurich. Numerical calculations were run on the Quest computing cluster at Northwestern University; XSEDE allocation TG-AST140023; and NASA HEC allocation SMD-16-7561 and SMD-17-1204.

Please contact the corresponding author if you have a sharing request for the data underlying this article. However, this small difference does not have any impact on the final SMBH mass.