Roughly 50 years ago, Lynden-Bell (1969) suggested that the high mass-to-light ratio objects (active galactic nuclei, AGN) observed at the centre of galaxies are actually supermassive black holes (SMBHs). Today it is commonly accepted that the luminosity of AGN is powered by accretion onto SMBHs (Soltan 1982) and almost all massive galaxies host SMBHs of millions to billions of solar masses at their centres (Lynden-Bell & Rees 1971; Rees 1984; Kormendy & Richstone 1995) . The universal existence of SMBHs in the E-mail: onurc@physik.uzh.ch (UZH ICS) centres of galaxies has raised questions about their nature, formation mechanisms, and relation with the environment.
Observational studies in the past two decades have shown strong correlations between the SMBH mass and various properties of their host galaxies such as the bulge mass (Kormendy & Richstone 1995; Marconi & Hunt 2003; Häring & Rix 2004; McConnell & Ma 2013; Kormendy & Ho 2013; Graham & Scott 2015) , the stellar velocity dispersion of the bulge (Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Gültekin et al. 2009; Bennert et al. 2015) , and the stellar mass of the host galaxy (Reines & Volonteri 2015; Beifiori et al. 2012; Savorgnan et al. 2016) in the local Universe. A good understanding of the scaling relations 2Çatmabacak et al.
is essential, as they may reflect a co-evolution of SMBHs and galaxies. While the scaling relations at low redshift are empirically well constrainted, 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 undermassive 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 reproduce the local relations without selfregulating black hole feedback? If so, how do the SMBHgalaxy 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-selfgravitating 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 angular momentum of the inflowing gas (Hopkins & Quataert 2010 . 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 M 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 Illustris (Vogelsberger et al. 2014; Genel et al. 2014) , Horizon-AGN (Dubois et al. 2014; Volonteri et al. 2016; Kaviraj et al. 2017) , Eagle (Schaye et al. 2015) , MassiveBlack (Khandai et al. 2015) , BlueTides (Feng et al. 2016) , Romulus (Tremmel et al. 2017) , and Illustris-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 disk 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 sim-ulations to reproduce the local scaling relations without the need for self-regulatory black hole feedback (Anglés-Alcázar et al. 2013 , 2017a 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, hydrodynamical simulations. The simulations used in this paper (MassiveFIRE) 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 large number of galaxies in our sample minimizes statistical variance and allows us to address the questions listed above with a statistically meaningful set of simulated galaxies over a wide range of redshifts (12 <z< 2) and halo masses (10 < log(M h /M ) < 13.5).
The outline of the paper is as follows; section §2 introduces the simulation properties. Section §3 lays out the details of our post-processing analysis. We present our main results in the following section §4. Specifically, section §4.3 presents a toy model to explain the physical origin of the M * − M BH scaling relation.
We use 34 high-resolution, cosmological zoom-in simulations from the MassiveFIRE 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 . Initial conditions were generated using the multi-scale initial condition tool MU-SIC (Hahn & Abel 2011) . 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.
For FIRE-1 simulations (Feldmann et al. 2016 , the gravitational softening length of dark matter and star particles is 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. For FIRE-2 simulations (Anglés-Alcázar et al. 2017c), the gravitational softening lengths of dark matter and star particles are 57 and 7 pc. The minimum softening length of gas particles is minimum of 0.7 pc. All the simulations presented here have a mass resolution of 1.7×10 5 M for dark matter particles and 3.3 × 10 4 M for gas and star particles.
FIRE-1 simulations do not include black hole sink particles, i.e., the growth of black holes is fully modeled in postprocessing. FIRE-2 simulations include black hole particles with a gravitational softening of 7 pc (physical) but no black hole feedback. We adapt the same post-processing approach for the analysis of both FIRE-1 and FIRE-2 simulations.
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 photoionization and photo-electric heating (Hopkins et al. 2014 (Hopkins et al. , 2018 . We refer the reader to Feldmann et al. (2016 Feldmann et al. ( , 2017 ; Feldmann (2017) ; Anglés-Alcázar et al. (2017c); Hopkins et al. (2014 Hopkins et al. ( , 2018 for more detailed information about the simulations and their sub-grid models. Table 1 provides an overview of the simulations used in the present work. 
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 MassiveFIRE simulations. Identified virialized structures contain at least 100 particles (M halo ∼ 10 7 M h −1 ). In all simulations, black hole particles are added in post-processing. The positions of the SMBHs are defined using either the center of mass (COM) or the maximum density center (MAX) provided by AHF. While our FIRE-2 simulations include black hole particles, we follow the same post-processing methodology in both FIRE-1 and FIRE-2 simulations. The properties of the black hole particles in the FIRE-2 runs are only used for the purpose of validating the outlined post-processing approach.
AHF calculates COM as the centre of mass of gas, star, and dark matter particles on the finest refinement of the host halo. Therefore, COM represents the average density regions in the host halo's centre. MAX is calculated as the highestdensity cell in the end-leave grid. Maximum density central regions can be at gas-rich regions or dense star clusters. Consequently, the centre of mass of star, gas, and dark matter particles can be quite different from the maximum density centre of particles. The different centring choices thus serve as proxies for two different SMBH locations; the average density (COM) and the densest (MAX) central regions. The difference between the two central coordinates are shown in physical units in Figure A3 for an example simulation, A1 with FIRE-2 physics. The discrepancy between different methods is around 0.5 − 1.0 pkpc at relatively high redshift (z 5) and drops below 100 pc in physical units towards late times. The discrepancy of central coordinates may suggest a transition from a chaotic structure to well-settled disc galaxies (Sparre et al. 2017; Stern et al. 2020) .
In the early Universe, the centering methods we used can lead to gas poor central regions. Strong stellar feedback is the most likely physical mechanism responsible for such evacuated central regions. The dynamics of the SMBH determines its location in the galaxy. Therefore, it can have an important effect on early SMBH growth, in addition to the suppression by stellar feedback (Anglés-Alcázar et al. 2017c) which is seen in either case.
Once the halo centres are identified, we use AHF Merg-erTree to map the haloes between snapshots. Our sample consists of the most massive halo in each simulation 3 AMIGA Halo Finder: http://popia.ft.uam.es/AHF/ Download.html and all haloes above 10 10 M for 6 simulations that reach z ≤ 1. Selected samples of haloes and their progenitors satisfy the following selection criteria. The mass fraction of high-resolution dark matter particles is larger than 98% so that haloes are not significantly polluted with low-resolution dark matter particles. Progenitor haloes have at least 10 7 M of stellar mass within 10% of their virial radii. The radii of the central galaxies are defined as 10% of their virial radii. We define the velocity of each SMBH as the average velocity of the 100 youngest star particles in its vicinity (within R 0 , see §3.2.1). It is useful to note here that our post-processing analysis is performed in the rest-frame of the SMBH of the given galaxy.
Our approach of separating the stellar mass of the galaxy into a bulge and a disk component follows Anglés-Alcázar et al. (2014) . Specifically, the bulge mass in an enclosed region with a given radius r is calculated as twice the mass of all star particles that are counter-rotating (v φ < 0) according to the unit stellar angular momentum vector of the galaxy. The stellar disk mass within r is then defined as the difference between the stellar mass within r and the bulge mass within r. The total disk mass is defined as the sum of the stellar disk mass and the gas mass within r. 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 is shown in appendix A.
The evolution and growth of SMBHs is treated in postprocessing. A (virtual) SMBH seed is placed at the centre of each progenitor halo. This choice places early SMBHs on the local M * − M BH relation when 10 4 M SMBH seed masses are adopted (our fiducial option). Lower or higher seed masses result in high-redshift SMBHs that start either below or above the local M * − M BH scaling relation.
We need to differentiate between mergers of SMBHs and 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 merger of their SMBHs. We thus analyze three different SMBH merger scenarios for the remainder of this paper.
Our first model ("all mergers") 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, therefore the merger contribution is overly optimistic. As a result, this scenario serves as the most optimistic case for SMBH growth.
Our second model ("few mergers") poses more stringent constraints on SMBH mergers and is based on the distance d sep of the two parent halos when they are identified as separate (sub-)halos for the last time (z merge ). Specifically, the SMBHs at the centers of these halos are assumed to have merged if (i) their d sep is smaller than 10% R vir of the more massive of the two halos or (ii) their dynamical friction timescale is smaller than the Hubble time at z merge . We adopt the following analytic estimate of the dynamical friction time (Binney & Tremaine 1987) :
Here, M h1 and M h2 are the first and the second halo mass, calculated by AHF. Furthermore, Λ is Coulomb logarithm where ln(Λ) = ln(1 + M h1 /M h2 ). R vir is the virial radius and V vir is the circular velocity of the more massive of the merging haloes. It is important to be careful about how the mass of the second halo is calculated. Dark matter particles of the second halo begin to be stripped when the two haloes are sufficiently close. To mitigate the effect of stripping -which may extend the merging timescale-the maximum mass of the progenitors of the second halo is used instead of the current halo mass. The second halo typically reaches its maximum mass 60 − 100 Myrs before the halo merger.
Our third option ("no mergers") is to neglect SMBH mergers altogether, so that SMBH seeds only grow via mass accretion throughout cosmic time.
Hopkins & Quataert (2011) model gas accretion from kpc to sub-pc scales due to 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 R 0 . All the terms in Equation 2 have a radial aperture dependency except the SMBH mass and constant m . Therefore, a change in R 0 means a change in the other parameters, too. Here, α T is a function of nuclear star formation law (see equations 39 and 65 in Hopkins & Quataert (2011) ), however for practical reasons we treat α T is set to 5 (see Figure 10 in Hopkins & Quataert (2011) ). The mass retention rate m is an overall normalization that covers the effects on gas dynamics (i.e. stellar and BH feedback) at unresolved scales (Anglés-Alcázar et al. 2017a). The prefactor α T is a function of disk, bulge and gas surface densities and is equal to 1 − 10. We adopt α T = 5 to reflect an average value.
The disk fraction f d is the ratio of the disk mass (stellar disk mass + gas mass) to total baryonic mass (stellar mass + gas mass). The accretion rate scales super linearly with the disk fraction, M ∝ f 5/2 d . M BH is the mass of the SMBH and M tot is the total mass within the radial aperture R 0 . The total mass is the sum of dark matter and baryonic matter inside R 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 R 0 , the total mass is the determining factor for the accretion rate onto SMBH. Amongst the remaining parameters, f gas is the ratio of the gas mass to the total baryonic mass inside R 0 and f 0 ≈ 0.31 f 2 d (M tot /10 9 M ) −1/3 . We update black hole masses in our post-processing analysis iteratively for each accretion model, as briefly explained in §3.2, M BH,i+1 = M BH,i + ∆t × M 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 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) , G is the gravitational constant, and ρ is the volume density of gas particles within R 0 . The bulk velocity of gas and the sound speed are denoted by v bulk and c s , respectively.
Nonetheless, the Bondi model is only valid for the special case of hot virialised gas with negligible angular momentum and radiative cooling. 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, σ ∼ GM enc (r)/r, 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. 2015b; 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 test a model where the BHAR scales with the free-falling gas inside R 0 :
Here, γ is a scaling factor that controls the percentage of Table 2 . SMBHs are placed in the densest centres (MAX) and few mergers model is used to model SMBH growth. A larger R 0 slightly increases the SMBH mass but overall the choice of R 0 has little impact on the growth history of the SMBH. In contrast, the SMBH scales approximately linearly with m .
free-falling gas accreted onto SMBHs. The free-fall timescale of the gas is denoted by t dyn .
In this section, we analyze how the choice of our model parameters affects the predicted SMBH growth. Our postprocessing analysis includes several key parameters that influence the growth of SMBHs such as the SMBH seed mass (M seed ), the black hole merger treatment and, most importantly, the accretion model (Equation 2 -6) and its free Figure 2 . 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 few mergers 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 M 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.
parameters, e.g. the mass retention rate ( m ) and the radial aperture (R 0 ) for the GTDA model. We will use the fiducial values for our model parameters in Table 2 and vary one parameter at a time to see the clear impact of our model parameters.
We tested various values for the radial aperture. As top panel of Figure 1 shows, the growth history of SMBHs is only mildly affected if R 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, M ∝ R (2010), a smaller R 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 R 0 = 100 pc as the fiducial value.
The mass retention rate m is a normalization of the overall gas accretion rate. As the bottom panel of Figure 1 highlights, varying the mass retention rate creates a noticeable shift in the normalization of the SMBH mass, at any redshift. Hence, we can adjust the normalization 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. Figure 2 shows the effect of 'convergence' of black hole masses in the GTDA model (Anglés-Alcázar et al. 2013 , 2017a . Specifically, SMBH seeds with masses 10 − 10 5 M result in a similar SMBH mass by z ∼ 4.5 if accretion is limited to less than 10 times Eddington and by z ∼ 3 in the case of 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 z > 3 − 5 may provide useful constraints on the masses of the first black hole seeds. Figure 3 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 (M BH < 10 6 M ) at 3 < z < 6, the seed mass contribution can be significant if relatively heavy seeds (here 10 4 M ) 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 in Appendix for the same analysis with different seed masses. Figure 4 compares the results for SMBH growth via GTDA 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, disk fraction, freefall 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 M at z > 7 to SMBHs of 10 8 M or more by z = 1. The biggest difference between these two groups of accretion models is the necessity of self-regulating AGN feedback to reproduce the 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 5 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. 2017b) . Gas fraction in the central 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 centering methods.
4.2.1 M * , M b , and σ vs M BH Figure 6 shows the relation between 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 (M BH 10 7.5 M ) in massive galaxies (M * 10 11 M ). Different combinations of SMBH location and merger treatment result in different tracks ( see Volonteri 2012 for details) in the M * − M BH plane, notably at high redshift. For each case, efficient SMBH growth starts around a similar stellar mass threshold ∼ 10 10 M .
The growth trend of SMBHs shown in Figure 6 can be divided into three phases. During phase I (M * < 10 8.5 M ), 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 6 ). In this scenario, the average SMBH-galaxy coevolution follows the local scaling relation for AGN galaxies in Reines & Volonteri (2015) . In contrast, the early growth (lower panels in Figure 6 ) of SMBHs is shifted if they reside in average density central regions. Reines & Volonteri (2015) . Early growth of SMBHs is suppressed when SMBHs are placed at the center of mass of the halo, which is more strongly affected by stellar feedback compared to the maximum-density center (bottom panels). 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(M * /M ) < 10). Efficient SMBH growth starts when the stellar mass reaches ∼ 10 10 M . SMBHs grow at a similar rate as their hosts in massive galaxies (M * > 10 11 M ). In all four 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.
During phase II (8.5 < log(M * /M ) < 10), SMBH mergers play an important role. In particular, they increase the mass of the most massive central SMBH at fixed stellar mass of the host galaxy. In contrast, SMBH growth stalls during this phase if SMBH mergers are not considered. In our sample, the contribution from mergers peaks at 2 < z < 6 (see Figure 3 ). 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 M SMBHs in ∼ 10 11 M galaxies, see Figure 6 . While this result disagrees with other 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 feedback. For any accretion model that has a black hole mass dependency in the form of BHAR ∝ M p BH for p > 0 (e.g. p = 1/6 for the GTDA model), the specific BHAR (sBHAR ≡ BHAR/M BH ) scales with M p−1 BH . In the case of Bondi-like 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.
Furthermore, high values of p cause heavier SMBH seeds to grow faster compared to light SMBH seeds in the early Universe. Thus, BHAR with a super linear dependence on the SMBH mass (p > 1) may introduce a strong sensitivity of BH growth on the seed mass. On the other hand, SMBH masses are asymptotically insensitive to the seed mass if BH accretion is sublinear (p < 1) ical drivers of delayed vs. efficient SMBH growth. Finally, once galaxies reach a stellar mass of M * 10 11 M , this accelerated growth comes to an end. After this time, galaxies and SMBHs grow again at a similar rate. The relation between stellar bulge and SMBH mass shown in Figure 7 is qualitatively similar to the relation between galaxy stellar mass and SMBHS mass ( Figure 6 ). Typically, the M b − M BH relation predicted by our sample at high redshift falls below the local M b − M BH scaling relation Häring & Rix (2004) . Our post-processing 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 for a comparison of the M b -M BH scaling relation predicted by our postprocessing analysis with the corresponding on-the-fly and post-processing predictions by Anglés-Alcázar et al. (2017c). The early growth of SMBHs is suppressed if SMBHs are placed in typical regions in the center of galaxies instead of the densest regions (COM model vs MAX model). In either case, the resulting SMBH masses at low redshift are not affected by this SMBH placement choice. This result can be understood from Figure A3 which shows that at late times (z 4) the center of the galaxy becomes well defined and 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: L IR = SF R × 1.49 × 10 10 L (Kennicutt 1998) and L bol = /(1 − ) × Mc 2 , see Figure 8 . In general, we find good agreement between our model predictions and observations at z < 2 and large differences at high redshift. Assuming our sample is a good representation of less luminous AGN and galaxies at high redshift, we predict a large number of yet unobserved lowluminosity AGN at high redshift. However, we do not exclude the need for additional accretion mechanisms besides GTDA to model ultraluminous AGN at high redshift. Figure 8 compares SFR and BHAR of our sample with the ones in galaxies hosting luminous AGN for z ≤ 4 (lefthand panel) and z ≥ 4 (right panel). The BHAR-SFR ratio scatters around a mean value 1/10 3 which is consistent with what Mullaney et al. (2012) predict for main sequence AGN at z = 1. High SFRs and BHARs seem to be one of the characteristic features of the most luminous AGN beyond z = 2, while MassiveFIRE simulations produce relatively moderate SFRs and BHARs. Figure 9 shows how the BHAR/SFR ratio and the BHAR/ M * ,in ratio scale with the total stellar mass of the host galaxies. Here, M * ,in refers to the rate of change of stellar mass within R 0 and we show results for the A series of simulations. The average behaviour of the BHAR/SFR ratio at 4 ≥ z ≥ 1 as predicted by our simulations has a similar slope (slope= 0.79) as predicted by Delvecchio et al. (2019) (slope= 0.73) based on an empirical model. However, the normalizations differ by up to an order of magnitude. Our predicted BHAR/ M * ,in ratios agree well with the results presented in Delvecchio et al. (2019) .
Observation and theory show that the average BHAR (BH AR) and the average SFR (SF R) correlate well (Mullaney et al. 2012; Calhau et al. 2017; Dai et al. 2018b; Volonteri et al. 2015b) . While the cause of this correlation is not yet fully understood, several explanations have been proposed.
The first explanation refers to a common cause. Gas reservoirs in the galaxy are the primary source for BH feeding and SF (Anglés-Alcázar et al. 2015) . However, gas consumption may lead to the injection of energy and momentum to the ISM. On one hand, the negative feedback by AGN and supernovae regulates these gas reservoirs via the energy and the momentum injected to the surrounding (Granato et al. 2004; Di Matteo et al. 2005; Menci et al. 2006) . BHAR and nuclear SFR can be suppressed by a factor up to ∼ 30 as result of both AGN and stellar feedback .
Another explanation suggests that AGN activity can lead to enhanced SFR. 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) .
SN feedback can expel the surrounding gas and suppress BH feeding. However, the gas ejected by SN feedback can cool down at the outskirts of the galaxy and return to the galaxy, a process known as "Galactic fountain" (Shapiro & Field 1976; Bregman 1980; Melioli et al. 2008) . It is also possible to accomplish such a replenishment in the form of positive SN feedback through which the feedback steals part of the hot gas from the intergalactic medium (IGM) and cools it down to constitute a galactic fountain (Hobbs & Feldmann 2020) . These findings suggest that the refurbished cold gas budget may be crucial to sustain SF and black hole feeding, especially in the local Universe. Volonteri et al. (2015a,b) discuss the time variability of BHAR and SFR before and after the merger phase at different scales both theoretically and observationally, using a set of simulations from Capelo et al. (2015) . Their analyses suggest that BHAR better correlates with the nuclear SFR rather than the total SFR except during galaxy mergers, as also suggested by Hopkins & Quataert (2010) . The total SFR of the host galaxy correlates well with BHAR during galaxy mergers because global dynamics becomes more important than the local processes in terms of angular momentum loss (Volonteri et al. 2015b ). Timescales of the nuclear star formation and the accretion onto SMBHs show similar patterns that are close to the dynamical timescale of matter in the nuclear region (∼ 100 pc). Figure 10 presents how SMBH mergers alter the slope of the scaling relations. Mergers increase the mass of SMBHs thus increasing the normalization of the scaling relations. The slope of scaling relations also decreases with increasing redshift. Therefore, SMBH mergers can have a noticeable impact on the SMBH-galaxy scaling relations. Figure 11 contains a literature compilation of the corre- lation between the properties of host galaxies such as stellar, bulge, and dynamical mass and SMBH mass for 0 z 7.5 (see the caption of Figure 11 for the references). 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. We assume an evolution of the form ∆ log(M BH ) = δ 1 + δ 2 log(1 + z) for the z-∆ log(M BH ) relation. Here δ 1 is the normalization constant and δ 2 is the slope of the z − ∆ log(M BH ) relation which changes from 1.38 ± 0.11 to 1.57 ± 0.10 for different fit functions. For illustrive purposes, Figure 11 shows the SMBH mass offset calculated with respect to our best fit line for the M * − M BH scaling relation. The offset from our best fit line gives a slope of 1.50 ± 0.10, pointing to a super linear correlation. Merloni et al. (2010) report a positive slope of δ 2 ∼ 0.68 ± 0.12 for their sample between 1 < z < 2.2, as well as Ding et al. (2020) reporting δ 2 = 1.03 ± 0.25 for 0 < z < 1.7. If we make an estimation about the host galaxy stellar mass by assuming a factor of 2 increase in bulge masses and decrease in dynamical mass, this would lead to a slope of 1.75 ± 0.10.
The slope we obtain with the observational sample is ≈ 2.5 times bigger than the prediction of Merloni et al. (2010) . This difference would be caused by the amount of data points and the redshift range as well as the selection effects, as only the brightest sources at higher redshifts are detectable. There is a clear trend of increase in the offset of SMBH mass towards higher redshift. The upper limit of the SMBH mass offset seems to be in agreement at different times, although the lowest data points in Figure 11 increase with redshift. This result further suggests that selection effects can create biases in the high redshift AGN sample. 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 precious to have a better understanding of the redshift evolution of ∆ log(M BH ).
A similar analysis for our MassiveFIRE sample is shown in Figure 12 . We only include SMBHs with host galaxies that 14Çatmabacak et al.
have stellar mass larger than 10 9 M , to stay consistent with the observational data from the literature. The offset of M BH from the best fit line has a negative slope for all cases we studied in our post-processing analysis. Among these different SMBH locations and merger treatments choices, the strongest negative slope always belongs to the SMBHs located at the densest regions (MAX).
In this section, we introduce a toy model in order to offer additional insight into the physical origin of the M * − M BH scaling relation in the context of the GTDA model. We start by making the following assumptions: First, we assume that the total mass in R 0 is dominated by the stellar mass, M tot ∼ M * . Second, BHAR is assumed to be proportional to the normalized stellar mass in the central 100 pc region. In Equation 2, the total stellar mass is the most dominant term, therefore the remaining terms can be approximately treated as a normalization constant. With these simplifications, we can link the growth rates of SMBHs and their host galaxies via the following toy model.
In the equations above, M * is the change of stellar mass per unit time. M * includes contributions from star formation, stellar outflows, and galaxy mergers. Integration of the stellar growth rate at the galactic and nuclear scales (∼ 100 pc) allows us to calculate the total and central stellar masses of galaxies. The total mass within R 0 is the most dominant term in Equation 2 and we can use Equation 8 to derive BHAR as a function of the central stellar mass.
In Equation 10, we derived SMBH mass as a function of the stellar growth rate. The growth of SMBHs and their host galaxies are linked to the growth of central and total stellar mass in our toy model, as shown in Equation 7-10. Thus, it is possible to model the slope of the M * − M BH scaling relation as a function of central and total M * . The slope of the M * − M BH scaling relation can be written as the ratio of specific BHAR (sBHAR) and the specific M * (s M * ). The growth of host galaxies at different scales will determine the slope of the M * − M BH scaling relation.
Only two simulations were run to z = 0 (A1 and A4 with FIRE-1 physics) and they may not be representative of the average behaviour of the stellar growth rate at relatively late times. Therefore, we calibrate our toy model using the average stellar growth rates from MassiveFIRE simulations only until z = 1. At lower redshift we consider the following possibilities: 1) the central M * decreases ∼ 2 orders of magnitude while A is constant ; 2) "A" constantly decreases while the central M * stays constant; and 3) central M * and A are constant while total M * increases. The purpose of exploring various scenarios is to show that the three scenarios listed above can culminate in similar results for the slope of M * − M BH scaling relation at lower redshift, where galaxies start outgrowing their SMBHs.
The host galaxies grow faster than their SMBHs until z = 4 since the specific growth rate of galaxies is much bigger than the SMBH growth rate. However, in the first scenario of Figure 13 , the central stellar growth rate decays two orderof-magnitude after z = 1. Such a fast decay means that an inside-out quenching of the stellar growth rate is in process. Inside-out quenching may happen due to the strong AGN feedback as it was found earlier (Tacchella et al. 2015; Ellison et al. 2018; Abdurro'uf 2018; Tacchella et al. 2018 ). This scenario is ruled out for the MassiveFIRE sample since these simulations do not include AGN feedback.
For 4 > z > 1, the central stellar growth rate still increases while the total stellar growth stays constant. The accelerated stellar growth of the central region creates a denser galactic core and a deep gravitational potential which helps the host galaxy to build a steady disk. It is useful to note here that the stellar mass of the host galaxy reaches 10 10 M around z ∼ 3 − 4, which intersects with the time at which the slope of the scaling relation (≡sBHAR/s M * ) starts climbing over unity.
Until z = 1, the central stellar growth rate is the same for all columns in the figure. We consider three different scenarios to reproduce the decreasing slope of the M * −M BH scaling relation at late times. The left-hand column of Figure 13 addresses a scenario in which the central stellar growth rate decays while the total stellar growth rate stays constant beyond z = 1. The decaying central M * roughly keeps BHAR at a constant value for the late times together with the constant A. Thus, the specific SMBH growth rate slows down and the ratio of specific growth rates of the host galaxy and SMBH start to decay.
The central column portrays a similar decay trend of the slope of the M * − M BH scaling relation by keeping the central stellar growth rate constant while A decreases for z 1. BHAR converges to a steady accretion rate and the stellar mass at the galactic centre reaches 10 10 M , almost 20% of the total stellar mass of the host galaxy. A decrease in A or central stellar growth rates would make it possible for the accretion rate to stay constant for the late times of SMBH-galaxy co-evolution.
The disk fraction is the second most dominant term in Equation 2, A ∝ f 5/2 d . Here A decreases for z 1, which means that the disk fraction also decreases, which is consistent with a scenario in which galaxies make a transition from spiral to elliptical morphology. As Davis et al. (2019) discussed, the decrease in the slope of the M * − M BH scaling relation is consistent with the scenario in which AGN host galaxies evolve morphologically into elliptic galaxies.
The right-hand column keeps the central stellar growth rate and A constant, which results in an increasing accretion rate at late times. In the mean time, the total stellar growth rate increases, which points out the contribution from minor galaxy mergers. As a result, the host galaxies start to grow faster than their SMBHs, again. In all these three cases, it is possible to have a decaying slope for the M * − M BH scaling relation. It is likely that all these possibilities contribute to the evolution of the slope of the M * − M BH scaling relation at lower redshift.
A more realistic picture would suggest a combination of the three scenarios to explain the final stages of SMBHgalaxy co-evolution. In such an evolution scheme, the relative growth of galaxies compared to their SMBHs at late times can accelerate due to the inside-out quenching of SFR caused by strong AGN feedback, transition from spiral to elliptical morphologies and minor galaxy mergers.
Our post-processing analysis and toy model rely on a number of simplifications. Hence, they are subject to several caveats.
The GTDA model was developed for disk galaxies and estimates the accretion from circumnuclear (∼ 100 pc) to sub-pc scales using the properties of the circumnuclear region, including its stellar mass, baryonic disk fraction, gas mass, and so on. However, especially at high redshift, the host galaxies may not be very well settled structures. Hence, it is not obvious that the GTDA model is a good description of SMBH growth at high redshift. However, we find that a similar growth history of SMBHs can be obtained by using the dynamical accretion model (Equation 6) with γ = 10 −4 . This suggests that details of the accretion model are less important as long as only a small fraction (0.1%) of the available gas is accreted by the BH per free-fall time, see also Anglés-Alcázar et al. (2017c) .
None of the SMBHs in our post-processing analysis are as massive as the most massive, observed SMBHs at z 6 (Wu et al. 2015; Mortlock et al. 2011; Bañados et al. 2018) . While some of our haloes are sufficiently massive ( 10 12 M ) to potentially host a luminous AGN, simple number density and clustering arguments show that our simulation volume is likely too small to host even a single luminous AGN. In particular, the number density of luminous AGNs at z ∼ 6−7 is ∼ 1cGpc 3 for M 1450 < 26 (Wang et al. 2019 ) and the duty cycle is close to the unity for z 6 (Shankar et al. 2010b (Shankar et al. ,a, 2019 . Our largest box size is ∼ 1Gpc in extent. Hence, if we simulated the whole box, we would expect to find one luminous AGN on average. However, the actual number of luminous AGN in our sample will be much lower as our simulation suite includes only three halos from this box.
Gravitational torques are an important physical mechanism to remove angular momentum from gas, but they may not necessarily be the most dominant process at all scales to drift gas radially inward. There may be other mechanisms as well, such as chaotic accretion of hot gas (Thomas et al. 2019) or other possible effects of mergers as in mergertriggered accretion (Capelo & Dotti 2017; Ricarte & Natarajan 2018) . In the future, we aim to assess the importance of these accretion models by modelling the cosmic evolution of the accretion rate with two different SMBH growth modes, considering the accretion of hot, pressure-supported gas in addition to the accretion of cold, rotationally supported gas driven by gravitational torques Thomas et al. 2019) .
Most of the MassiveFIRE 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. Importantly, as we showed in Figure 6 , accretion models that are weakly dependent on the SMBH mass can reproduce the M * − M BH local scaling relation without AGN feedback (Anglés-Alcázar et al. 2013) . We leave similar studies of BH growth in simulations with AGN feedback for the future.
Our SMBH merger treatment considers an instantaneous black hole merger once the two haloes merge, regardless of the mass ratio of haloes. Therefore, our analysis overestimates the number of merging SMBHs. Such an overestimation of SMBH mergers boosts the effect of the mergers on the growth history of SMBHs and the SMBH-galaxy scaling relations at early times.
We study the effect of different SMBH locations, SMBH merger treatments, and accretion models on the SMBH growth, SMBH-galaxy scaling relations, and the redshift evolution of the M * − M BH scaling relation. Our findings are:
• We have analyzed the growth history of SMBHs and the co-evolution with their host galaxies at relatively high redshift (z 2) in a post-processing analysis using 34 high resolution cosmological zoom-in simulations from the Mas-siveFIRE suite. We find that different SMBH locations and SMBH mergers may change the slope of the M * −M BH scaling relation at high redshift in the absence of AGN feedback, although these aspects have no clear effect on the final SMBH mass in the local Universe.
• Large-scale accretion models can be divided into two major groups. The first group consists of Bondi-like models which have a super linear dependency on the SMBH mass. The second group of accretion models have almost no dependency on the SMBH mass. They model large-scale accretion using the large scale properties of the host galaxy. Contrary to the first group, the second group do not require self-regulating AGN feedback to reproduce the local scaling relations, at least for the AGN host galaxies. SMBHs lose their memory of seed mass choice in the second group where p < 1, while the seed mass choice is important for the first group to determine the final SMBH mass.
• It is possible to reproduce the different evolutionary phases of the M * − M BH scaling relation, considering the different SMBH locations and SMBH merger treatments. The SMBH merger treatment and the seed mass choice have inconsequential impact on the final SMBH mass for the massive galaxies that we focus on in this paper. Nonetheless, they are important to set the trend of the scaling relations in the early Universe and the redshift evolution of ∆ log(M BH ).
• Currently, none of the SMBHs in our post-processing analysis are as luminous as the billion solar mass SMBHs in the early Universe, although we have a considerable overlap in the IR luminosity. The most luminous high-redshift AGNs are quite rare objects (≈ 1 cGpc −3 ) and studying them requires not only large enough box sizes but also the highest density environments.
16Çatmabacak et al.
• We find a slope for the evolution of the SMBH mass offset with respect to the average trend of the M * − M BH local scaling relation between 1.38 − 1.57 using various fit functions. ∆ log(M BH ) increases towards higher redshift, although the maximum offset values are similar for both the early and the late times. We propose a method to anticipate what the SMBH seed mass and SMBH location should be in order to reproduce the slope of the z − ∆ log(M BH ) relation.
• We develop a toy model that links the growth of SMBHs to the central stellar growth rate of the host galaxies, motivated by the GTDA model. Our two-zone toy model may provide insights about the properties of SMBH seeds by reproducing the slope of the redshift evolution of ∆ log(M BH ).
Using the same toy model, it is also possible to come up with possible scenarios to speculate about the decay of the slope of the M * − M BH scaling relation for z 1 in our sample.
Assuming that SMBH seeds are born in average density central regions (COM), the merger of SMBHs would be important to put constraints on the SMBH seed mass. If the birth environments of SMBH seeds are the densest central regions, these objects can be even detected by LSST, WFIRST, and JWST to put new constraints on the availability of 10 8−9 M AGN for z 7.
The low luminosity AGN (LLAGN, M mag −25) predicted in our post-processing analysis (M mag ∼ −21 at z = 7) may tend to lie on or even below the local M * − M BH scaling relation, especially at the high dynamical masses (M dyn 4 × 10 10 M ) (Matsuoka et al. 2016) . The increasing number of LLAGN have a potential to change the slope of the M * − M BH scaling relation at high redshift (z 6).
Massive galaxies (∼ 10 11 M ) at z 1 can grow faster than their SMBHs via devouring their satellite galaxies which causes the slope of the M * − M BH scaling relation to decline. A combination of the physical mechanisms such as inside-out quenching, morphological transition from spiral to elliptical galaxies, and minor mergers may be necessary to explain the flattening of the slope of the M * − M BH scaling relation at low redshift.
The slope of the M * − M BH scaling relation changes with redshift because the evolution of M BH and the growth of M * are linked in a non-linear fashion, see §4.3. Furthermore, the efficiency of SMBH merging also changes the slope of the relation. In particular, a high merging efficiency results in a steep slope of the M * − M BH scaling relation at high z similar to the local scaling relation, see Figure 12 . Consequently, the observational finding of a steep M * − M BH relation at high z would be suggestive of a relatively high SMBH merger rate at those times, thus potentially boosting the generation of gravitational wave signals in the early Universe.
We plan to extend this work and compare our current findings with the simulations that include AGN feedback. The importance of AGN feedback on SMBH growth, SMBHgalaxy co-evolution, and scaling relations needs to be addressed in the future. 
Our bulge-disk decomposition method (Anglés-Alcázar et al. 2014) works well for non-rotating bulges. In Figure A1 , we compare velocity dispersion of all star particles within 1 kpc, half stellar mass radius and whole galaxy with the stars in the galactic stellar bulge. There is no clear difference between the velocity dispersion calculated using different approaches.
In Figure 3 , we show the contribution of both accretion and seed mass to the total SMBH mass for a seed of 10 4 M . Figure A2 searches for the effect of different seed masses on the total SMBH mass. As the seed mass increases, the offset between total SMBH mass and in-situ accretion contribution increase, especially at high redshifts. Effects of seed mass choice are observable beyond z = 2 for seed masses equal to or above 10 4 M , which points out that the effect of small seeds on SMBH mass is not observable. Figure A6 produces the M * − M BH scaling relation for different seed masses, birth environments, and merger treatments. The effect of full merger treatment is distinguishable again for seed masses that are equal to or above 10 4 M . Increase in seed mass choice boosts the importance of mergers on the trend of scaling relation. Smaller seed masses that were born in gas-rich environments follow local scaling relation since the growth of black hole is dominated by in-situ accretion. Figure A5 shows the accretion rates of the SMBHs in our sample according to our fiducial model introduced in Table 2 . The Eddington ratio can exceed 1 at early times (z > 3) and for low mass black holes (< 10 6 M ). More massive black holes or black holes at later times almost never . SMBH accretion rates normalized to the Eddington rate for our sample of MassiveFIRE galaxies at different redshifts. All model predictions are based on our fiducial GTDA model, see Table 2 . Points are colored according to black hole masses. The average Eddington ratio λ decreases with decreasing redshift by about one order of magnitude between z = 10 and z = 0. SMBHs of mass > 10 6 M almost never exceed the Eddington limit. 