High-redshift quasars and their host galaxies I: kinematical and dynamical properties and their tracers

TL;DR

This study uses high-resolution simulations to compare the properties of high-redshift quasar host galaxies with local relations, highlighting the importance of tracer choice in interpreting black hole-galaxy correlations.

Contribution

It provides a detailed simulation-based analysis of $z=7$ quasar hosts, comparing gas and stellar tracers to local black hole-galaxy relations, clarifying discrepancies in high-redshift observations.

Findings

High-redshift quasars follow local MBH-galaxy correlations when using stellar tracers.

Gas-based tracers suggest MBHs are over-massive at high redshift.

Tracer choice significantly affects interpretation of black hole growth in early galaxies.

Abstract

Observations of high-redshift quasars provide information on the massive black holes (MBHs) powering them and the galaxies hosting them. Current observations of $z \gtrsim 6$ hosts, at sub-mm wavelengths, trace the properties of cold gas, and these are used to compare with the correlations between MBHs and galaxies characterising the $z=0$ population. The relations at $z=0$, however, rely on stellar-based tracers of the galaxy properties. We perform a very-high resolution cosmological zoom-in simulation of a $z=7$ quasar including state-of-the-art non-equilibrium chemistry, MBH formation, growth and feedback, to assess the evolution of the galaxy host and the central MBH, and compare the results with recent ALMA observations of high-redshift quasars. We measure both the stellar-based quantities used to establish the $z=0$ correlations, as well as the gas-based quantities available in $z \gtrsim 6$ observations, adopting the same assumptions and techniques used in observational studies. The high-redshift studies argued that MBHs at high redshift deviate from the local MBH-galaxy correlations. In our analysis of the single galaxy we evolve, we find that the high-redshift population sits on the same correlations as the local one, when using the same tracers used at $z=0$. When using the gas-based tracers, however, MBHs appear to be over-massive. The discrepancy between local and high-redshift MBHs seems caused by the different tracers employed, and necessary assumptions, and not by an intrinsic difference. Better calibration of the tracers, higher resolution data and availability of facilities that can probe the stellar population will be crucial to assess precisely and accurately high-redshift quasar hosts.

Figures

Tables

Table 1. Table 1: Best-fit parameters from the velocity maps of [CII] line emission, resulting dynamical masses, and comparison with the intrinsic quantities. The first column identifies the tracer (and the resolution) used for the estimate and the second is the size of the emitting region. Columns 3-5 report the kinematical results; the third column is the FWHM of the Gaussian fit to the line-of-sight velocity (when the approximation is valid) or the distance between the two peaks, the fourth column is the mean of the Gaussian profile fit to the velocity dispersion, and the fifth column is the stellar velocity dispersion, measured on the stellar particles. Columns 6-9 report the mass measurements; the sixth column is the dynamical mass assuming a rotationally-supported system, and the seventh assumes a dispersion-dominated system. In the last two columns, we report the total gas and stellar masses enclosed within the emitting radius R x subscript 𝑅 𝑥 R_{x} .

Tracer/	$R_x$	FWHM	${\sigma }_{\mathrm{fit}}$	${\sigma }_{\ast }$	$M_{\mathrm{dyn},\mathrm{rot}}$	$M_{\mathrm{dyn},\mathrm{sph}}$
Resolution	(kpc)	(km s-1)	(km s-1)	(km s-1)	$\left({10}^{10}{M}_{\odot }\right)$	$\left({10}^{10}{M}_{\odot }\right)$	$\left({10}^{10}{M}_{\odot }\right)$	$\left({10}^{10}{M}_{\odot }\right)$
${\rho }_{\mathrm{gas},\mathrm{sim}}$	1.10	490	115	350	5.1	1.7	3.3	9.1
${\rho }_{\mathrm{gas},200\mathrm{pc}}$	1.10	460	140	350	4.5	1.5	3.3	9.1
${\rho }_{\mathrm{gas},2\mathrm{kpc}}$	1.15	350	175	350	2.7	0.9	3.4	9.1
$L_{\left[\mathrm{CII}\right],\mathrm{sim}}$	0.96	430	35	350	$3.5$	$1.1$	3.0	8.9
$L_{\left[\mathrm{CII}\right],200\mathrm{pc}}$	0.96	470	55	350	$4.1$	$1.3$	3.0	8.9
$L_{\left[\mathrm{CII}\right],2\mathrm{kpc}}$	1.04	390	90	350	$3.0$	$1.0$	3.2	9.0

Equations

$${\dot{\rho }}_{\mathrm{SF}}=\epsilon \frac{{\rho }_{\mathrm{g}}}{t_{\mathrm{ff}}},$$

$$\epsilon =\frac{{\epsilon }_{\star }}{2{\varphi }_{\mathrm{t}}}\exp \left(\frac{3}{8}{\sigma }_s^2\right)\left[1+\mathrm{erf}\left(\frac{{\sigma }_s^2-s_{\mathrm{crit}}}{\sqrt{2{\sigma }_s^2}}\right)\right],$$

$${\alpha }_{\mathrm{vir}}=\frac{5[\Vert \nabla \otimes \mathbf{v}{\Vert }^2+(c_{\mathrm{s}}/L{)}^2]}{\pi \mathrm{G}{\rho }_{\mathrm{gas}}},$$

$${\dot{M}}_{\mathrm{BHL}}=\frac{4\pi {\mathrm{G}}^2{\mathrm{M}}_{\bullet }^2{\rho }_{\mathrm{gas}}}{(v_{\mathrm{rel}}^2+c_{\mathrm{s}}^2{)}^{3/2}},$$

$${\dot{M}}_{\mathrm{accr}}=⟨\frac{4\pi {\mathrm{G}}^2{\mathrm{M}}_{\bullet }^2{\rho }_{\mathrm{gas}}}{(v_{\mathrm{rel}}^2+c_{\mathrm{s}}^2{)}^{3/2}}⟩.$$

$${\dot{M}}_{\mathrm{Edd}}=\frac{4\pi \mathrm{G}{\mathrm{m}}_{\mathrm{H}}}{{\eta }_{\mathrm{acc}}c{\sigma }_{\mathrm{T}}}{M}_{\bullet },$$

$${\tilde{L}}_{[\mathrm{CII}]}\left(\frac{{\mathrm{L}}_{\odot }}{{\mathrm{M}}_{\odot }}\right)=\{\begin{array}{ll}0.1\left(\frac{Z}{Z_{\odot }}\right)\left(\frac{n_{\mathrm{gas}}}{100\mathrm{c}{\mathrm{m}}^{-3}}\right)& n_{\mathrm{gas}}<10^3\mathrm{c}{\mathrm{m}}^{-3}\\ 1.0\left(\frac{Z}{Z_{\odot }}\right)& \mathrm{otherwise}\end{array}$$

$$M_{x,\mathrm{dyn}}=\{\begin{array}{ll}\frac{R_x}{\mathrm{G}}{\left(\frac{0.75\mathrm{FWHM}}{\sin i}\right)}^2& \mathrm{rotationally}-\mathrm{supported}\\ \frac{3R_x{\sigma }_{\mathrm{FWHM}}^2}{2\mathrm{G}}& \mathrm{dispersion}-\mathrm{dominated}\end{array}$$

$${\dot{M}}_{\mathrm{BHL}}=\frac{4\pi {\mathrm{G}}^2{\mathrm{M}}_{\bullet }^2⟨{\rho }_{\mathrm{gas}}⟩}{(⟨v_{\mathrm{rel}}{⟩}^2+⟨c_{\mathrm{s}}{⟩}^2{)}^{3/2}},$$

$${\dot{M}}_{\mathrm{accr}}=⟨\frac{4\pi {\mathrm{G}}^2{\mathrm{M}}_{\bullet }^2{\rho }_{\mathrm{gas}}}{(v_{\mathrm{rel}}^2+c_{\mathrm{s}}^2{)}^{3/2}}⟩.$$

$${\dot{M}}_{\mathrm{accr}}=-{\int }_{\Delta V}\nabla \cdot ({\rho }_{\mathrm{gas}}{\mathbf{v}}_{\mathrm{rel}}).$$

$$\Delta E=L_{\mathrm{BH}}\Delta t={\eta }_{\mathrm{fbk}}{\eta }_{\mathrm{rad}}{\dot{M}}_{\mathrm{BH}}{c}^2{\left(\frac{3N_{\mathrm{ngb}}{m}_{\mathrm{gas}}}{4\pi \rho c_s^3}\right)}^{1/3}.$$

$$\Delta T=\frac{\Delta E}{N_{\mathrm{ngb}}{m}_{\mathrm{gas}}}\frac{(\gamma -1)m_{\mathrm{p}}}{k_{\mathrm{B}}}=\zeta \frac{{\eta }_{\mathrm{c}}}{(N_{\mathrm{ngb}}{m}_{\mathrm{gas}}{)}^{2/3}},$$

$$r_{\Delta T}=\frac{\Delta T_{\mathrm{low}}}{\Delta T_{\mathrm{high}}}=\frac{{\eta }_{\mathrm{fbk},\mathrm{low}}}{{\eta }_{\mathrm{fbk},\mathrm{high}}}{\left(\frac{N_{\mathrm{ngb},\mathrm{high}}{m}_{\mathrm{gas},\mathrm{high}}}{N_{\mathrm{ngb},\mathrm{low}}{m}_{\mathrm{gas},\mathrm{low}}}\right)}^{2/3},$$

$$r_{\Delta T}=10{\left(\frac{32}{64}\frac{1.5\times 10^4}{2\times 10^6}\right)}^{2/3}=\mathrm{0.24.}$$

Full text

††pagerange: High-redshift quasars and their host galaxies I: kinematical and dynamical properties and their tracers–16††pubyear: 2018

1 Introduction

Massive black holes (MBHs) are ubiquitously observed in the centre of massive galaxies at all redshifts. By spanning a mass range that goes from about $10^{5}\,\rm{M_{\odot}}$ up to $10^{9-10}\,\rm{M_{\odot}}$, they also represent the most massive compact objects in the Universe. Observationally, they are usually detected during accretion events, when a fraction of the energy of the accreting gas is released as radiation, making them shine as active galactic nuclei (AGN), or produces powerful collimated jets. According to Soltan’s argument (Soltan, 1982), these MBHs must have gained most of their mass via gas accretion, and they should have formed as smaller objects (called ‘seeds’) in the early Universe.

Recent observations of quasars at high redshift have shown that MBHs with masses around 10${}^{9}\,\rm{M_{\odot}}$ were already present when the Universe was less than 1 Gyr old (Fan et al., 2006; Mortlock et al., 2011; Bañados et al., 2018), posing tight constraints on the mechanisms proposed to explain their formation and early growth. A crucial requirement to improve our understanding of how these MBHs formed and managed to grow so quickly is to properly resolve the environment around these objects and the properties of their galaxy hosts. Unfortunately, while in the nearby Universe the contribution of the galaxy and the central MBH to the total emission can be easily separated, at high redshift we are mainly limited to the most luminous sources, which, because of the intrinsic angular size of the galaxy and the limited resolution, are dominated by the central source emission, preventing us from properly resolving the galaxy host.

The advent of current facilities like ALMA has been crucial for our understanding of these objects. By probing far-infrared lines that are uncontaminated by the quasar light and unaffected by dust extinction, in particular the [CII] line at 158 $\mu$m, i.e., the main coolant of the cold (30 K $<T<$ 3000 K) interstellar medium (ISM) and the CO line emission, ALMA is giving important information about the kinematics and dynamics of molecular gas within the quasar hosts. In the near future, the James Webb Space Telescope (JWST) will extend these capabilities to the near-infrared band, providing us with a unique set of tools that will finally help us to shed light on the peculiar conditions in which these MBHs form and evolve.

Compared to the local population, MBHs at high redshift seem to be overmassive by about an order of magnitude at a given galaxy mass (Walter et al., 2004; Venemans et al., 2017; Decarli et al., 2018), and they appear to deviate from the low-redshift $M$-$\sigma$ relation (Shields et al., 2006; Coppin et al., 2008; Wang et al., 2010). The trend of increasing ratio between MBH and galaxy mass with redshift appears also in lower redshift samples (e.g., Merloni et al., 2010; Decarli et al., 2010). However, intrinsic limitations in observational techniques and selection criteria can affect the estimates of the different properties or their statistical significance in particular at high redshift, where only the most massive/luminous objects in relatively common galaxies can be observed (e.g. Lauer et al., 2007; Vestergaard et al., 2008; Volonteri & Stark, 2011).

Recently, many authors have investigated this peculiar class of objects with semi-analytic models (e.g. Valiante et al., 2011, 2014; Pezzulli et al., 2016) and numerical simulations (Costa et al., 2014; Richardson et al., 2016; Di Matteo et al., 2017; Smidt et al., 2018; Barai et al., 2018), employing different techniques and addressing different questions. For instance, Di Matteo et al. (2017) focus on the conditions for the efficient growth of these MBHs, while Costa et al. (2014) and Smidt et al. (2018) assess the role of the AGN feedback and X-ray AGN radiation in suppressing star formation (SF) in the MBH vicinity and in creating an HII region around the galaxy. Richardson et al. (2016) and Barai et al. (2018) investigate the impact of AGN feedback in a $z=5$ proto-cluster of galaxies and in a $z=6$ halo of $10^{12}\,\rm{M_{\odot}}$, respectively.

In this study, we investigate the evolution of a quasar host at high-redshift through very-high resolution numerical simulations of a massive galaxy and its central MBH ($M_{\rm halo}=1.5\times 10^{12}\,\rm{M_{\odot}}$ at $z=7$), employing a detailed sub-grid modelling that also includes non-equilibrium chemistry of the primordial species, to follow directly the kinematics and dynamics of molecular hydrogen in the ISM.

This is the first of a series of paper addressing properties of high-redshift quasar hosts and their MBHs. In this paper, we present and discuss the main evolution of the target galaxy and its central MBH, focussing on the stellar and gas tracers (total gas and [CII] emission). In paper II (Lupi et al. in prep.) we will focus on the molecular hydrogen, directly traced in the simulation, and in Paper III (Lupi et al. in prep.) we will discuss the evolution of the entire MBH population forming during the simulation, and possible outflows from the MBHs.

2 Simulation setup

We study the evolution of a high-redshift quasar host by means of a zoom-in cosmological simulation performed with the hydrodynamic code gizmo (Hopkins, 2015).

2.1 The hydrodynamic code GIZMO

gizmo (Hopkins, 2015) is a particle-based code that descends from Gadget3 and Gadget2 (Springel, 2005). It implements a new method to solve the hydrodynamic equations that exhibits at the same time the intrinsic adaptivity and almost perfect conservation of angular momentum typical of smoothed particle hydrodynamics codes, and the excellent shock-capturing properties of grid-based codes.

The code is based on a partition scheme to sample the simulated volume via a set of tracer ‘particles’ corresponding to unstructured cells. Compared to moving-mesh codes (e.g. arepo; Springel, 2010), in gizmo the effective volume of each cell is computed via a kernel smoothing procedure. Hydrodynamics is then treated like in mesh-based codes, with the Riemann problem solved with a Godunov-like method. Here, we employ the mesh-less finite-mass method, i.e. the mass preserving one available in the code, with a cubic-spline kernel, for which we set the desired number of neighbours to 32. Gravity is solved via a Barnes-Hut tree, as in Gadget3 and Gadget2. We keep the gravitational softening of dark matter and stars fixed at 40 and 10 pc h*-1*. For gas we employ the fully adaptive gravitational softening, where the softening length is maintained equal to the kernel size of the particle in an intrinsically adaptive way. This avoids the main issues arising when two different resolutions are used (Hopkins, 2015). The maximum spatial resolution for gas, determined by the minimum allowed softening, is 2.5 pc h*-1*, corresponding to an inter-particle spacing of $\sim 5$ pc. However, this value is reached only in the very high density clumps just before a SF event.

We now summarise the sub-grid prescriptions employed in the simulation.

• •

The chemical network and cooling/heating processes, modelled via krome (Grassi et al., 2014), are the same as in Lupi et al. (2018), and include non-equilibrium chemistry for 9 primordial species (H,H+,He,He+,He*++*,H-,H2, H${}_{2}^{+}$, and e-) with H2 formation via H- associative detachment and on dust (see Bovino et al., 2016, for details), a metagalactic ultraviolet (UV) background Haardt & Madau (2012), and look-up metal cooling tables tabulated by Shen et al. (2013) and obtained with Cloudy (Ferland et al., 2013). The metal cooling tables take into account the first 30 elements in the periodic table, and assume photo-ionisation equilibrium of an optically thin gas with the extragalactic UV background as a function of $z$ (Shen et al., 2010). Because of the limited resolution, we also include a clumping factor $C_{\rho}=\exp(\sigma_{s}^{2})=1+b^{2}\mathcal{M}^{2}$ in the H2 formation rate on dust, where $\mathcal{M}$ is the Mach number.

• •

SF is implemented using a stochastic prescription that converts gas particles into stellar particles, where the star formation rate (SFR) density is defined as

[TABLE]

with $\varepsilon$ the SF efficiency parameter, $t_{\rm ff}=\sqrt{3{\rm\pi}/(32G\rho_{\rm g})}$ the free-fall time, $\rho_{\rm g}$ the local gas density, and $G$ the gravitational constant. $\varepsilon$ is computed according to the theoretical studies of turbulent magnetised clouds by Padoan & Nordlund (2011), based on the assumption of a Log-Normal gas density distribution with average density $\rho_{\rm g}$ and the width of the underlaying gaussian distribution $\sigma_{s}=\ln(1+b^{2}\mathcal{M}^{2})$. The $b$ parameter, representing the ratio between solenoidal and compressive modes, is set to 0.4, that corresponds to a statistical mixture of the two modes (Federrath & Klessen, 2012). $\varepsilon$ can be expressed as (Federrath & Klessen, 2012)

[TABLE]

where $\varepsilon_{\star}=0.5$ is the normalisation calibrated against observations (Heiderman et al., 2010), $1/\phi_{\rm t}=0.49$ is a fudge factor that takes into account the free-fall time-scale uncertainty, and $s_{\rm crit}(\alpha_{\rm vir},\mathcal{M})$ is the critical logarithmic density for SF and depends on the virial parameter $\alpha_{\rm vir}$, that determines how bound the cloud is, and $\mathcal{M}$ (Federrath & Klessen, 2012). As in Lupi (2019), we define

[TABLE]

where $L\approx 0.5h$ is the particle grid-equivalent size, with $h$ the particle smoothing length, $\sigma_{\rm eff}\equiv L\sqrt{\|\nabla\otimes\mathbf{v}\|^{2}+(c_{\rm s}/L)^{2}}$ the total support against gravitational collapse111We employ the velocity gradient before the slope-limiting procedure is applied, as in Lupi (2019)., and $c_{\rm s}$ the sound speed. SF is allowed for gas matching two criteria222As discussed in Lupi et al. (2018), the density threshold does not play any role and is kept only for numerical reasons, to avoid wasting time computing the SF rate in low-density, unbound regions.: i) $\rho_{\rm g}>m\rm_{H}\,cm^{-3}$ and ii) $\mathcal{M}>2$.

• •

Stellar radiation is implemented as in model b of Lupi et al. (2018), by collecting all the stellar sources in the gravity tree.

• •

Supernova feedback is based on the mechanical feedback prescription in Lupi (2019), a variation of the publicly available implementation by Hopkins et al. (2018), that is able to reproduce the terminal momentum of state-of-the-art high-resolution simulations (Kim & Ostriker, 2015; Martizzi et al., 2015) independent of resolution. However, a consensus on the terminal momentum from clustered SNe has not been reached yet, with different models varying from almost no variation with respect to the single SN case (Kim et al., 2017; Gentry et al., 2018) up to one order of magnitude higher values (Keller et al., 2014; Gentry et al., 2017). Here, to make SN feedback more effective, we assume a terminal momentum twice that in Lupi (2019). Stellar particles represent an entire stellar population following a Chabrier initial mass function (IMF; Chabrier, 2003). Here we consider both type II and type Ia SNe. SN explosions are modelled as discrete events, and not as a continuous source, as in Lupi (2019). We assume that stars in the range $8-40\,\rm{M_{\odot}}$ explode as type II SNe, and release $10^{51}$ erg of energy and an IMF-averaged ejecta mass $M_{\rm ej}=15.1452-M_{\rm NS}=13.7452\,\rm{M_{\odot}}$, with $M_{\rm NS}=1.4\,\rm{M_{\odot}}$ the mass of the remnant neutron star, an oxygen mass $M_{\rm oxy}=1.2403\,\rm{M_{\odot}}$ and an iron mass $M_{\rm iron}=0.10422\,\rm{M_{\odot}}$. Assuming the solar ratio for the Oxygen and Iron groups, the total metal mass injected by every SN is $M_{\rm Z}=2.09M_{\rm oxy}+1.06M_{\rm iron}=2.7028\,\rm{M_{\odot}}$. The type Ia SN rates are instead based on the delay-time distribution by Maoz et al. (2012), and cover the stellar age interval 0.1-10 Gyr. Every type Ia SN is assumed to inject $10^{51}$ erg, $M_{\rm ej}=1.4\,\rm{M_{\odot}}$ and $M_{\rm Z}=0.9604\,\rm{M_{\odot}}$ ($M_{\rm oxy}=0.14\,\rm{M_{\odot}}$ and $M_{\rm iron}=0.63\,\rm{M_{\odot}}$). Every time a SN explodes, we distribute mass, metals, momentum, and thermal energy among all the gas neighbours within the star particle kernel (defined as the sphere enclosing 64 neighbours) and those whose kernel size encompasses the star particle, as in Hopkins et al. (2018). Every particle receives a fraction of the total SN mass, metals and energy, according to the solid angle covered by the particle around the star, which is further symmetrised via a tensor renormalization to ensure the explosion is spherical in the stellar reference frame. The proper amount of thermal vs kinetic energy to be given to each neighbour is then determined according to the gas properties, so that the feedback is mostly thermal if the cooling mass is resolved and mostly kinetic if not (see Lupi (2019) for details).

• •

For low-mass stars, we only distribute the mass lost in slow winds and the initial momentum to the neighbours, without accounting for possible shocks with the interstellar medium.

The total mass recycling of our stellar feedback prescription is about 42 per cent of the initial stellar mass over a Hubble time (see Kim et al., 2014, and references therein). We also include an additional time-step limiter for the stellar particles in the simulation, to ensure single SN events are resolved. In fact, the time-step is limited to $\sim 5\times 10^{4}$ yr, i.e. 1/100th of a $40\,\rm{M_{\odot}}$ star lifetime, for particles younger than 100 Myr (dominated by Type II SN events). This limit is instead increased to 1/10th of the stellar population age when the stellar particle is older than 100 Myr.

2.2 Black hole physics

Here, we describe the prescriptions we implemented in gizmo to model BH seeding, accretion, feedback, and the merger of BH binaries.

2.2.1 BH seeding

To follow BH evolution in a cosmological context, we need BHs to consistently form during the simulation, similarly to stars. The focus of this study is on the properties of the host galaxies, not on the properties of seed BHs, therefore we devise a scheme that ensures numerical stability and results in BH masses consistent with the available observations by $z=7$.

We seed BHs in galaxies with a stellar mass $M_{\rm star}>10^{8}\,\rm{M_{\odot}}$, unless another BH is already present. Unlike for SF, a strictly local process, our prescription for BH seeding requires global information about entire galaxies. In order to collect this information on-the-fly in the simulation, we employ the Friends-of-Friends (FoF) algorithm available in gizmo. We first identify DM groups assuming a linking length $d_{\rm FoF}=0.05$,333The linking length is defined in units of the average inter-particle distance of the high-resolution DM particles in the box $\Delta x_{\rm DM}=[M_{\rm DM}/(N_{\rm DM}\rho_{z,\rm DM})]^{1/3}$, where $M_{\rm DM}$ is the total mass in high-resolution DM particles, $N_{\rm DM}$ is their number, and $\rho_{z,\rm DM}$ is the DM cosmological density at redshift $z$. Although the commonly assumed linking length for DM halo identification is $d_{\rm FoF}=0.2$, we found that a smaller value is a better choice to avoid contamination from sub-haloes and to limit the selection to the stellar-populated region. and then attach gas, stellar and BH particles to the group of the closest DM particle. Finally, we flag only the groups that match our criteria for seed BH formation. In all the flagged groups, a new seed is finally spawn at the position of the star surrounded by the highest gas density (typically close to the galaxy centre), with an initial mass $M_{\bullet,\rm init}=10^{6}\,\rm{M_{\odot}}$. The gravitational softening of BH particles is set to 5 pc.

Compared to other studies, our initial BH mass here is ten to hundred times higher. This choice is motivated by the requirement of avoiding spurious scattering of BHs out of the galaxy centre. This is a common issue in numerical simulations, when the BH mass is comparable or smaller than that of the other particles of the simulation, and can be solved by means of ad-hoc sub-grid prescriptions like i) an instantaneous re-centering (Springel et al., 2005; Sijacki et al., 2007; Booth & Schaye, 2009; Schaye et al., 2015; Barai et al., 2018), ii) a correction for unresolved dynamical friction (Dubois et al., 2014b; Tremmel et al., 2015) or iii) an initially larger dynamical mass decoupled from the BH mass used for accretion (e.g. Anglés-Alcázar et al., 2017; Biernacki et al., 2017). Unfortunately, all these prescriptions can introduce additional undesired effects in the simulation, like superluminal motions when the BH is far away (i), dynamical acceleration (ii) when the real dynamical friction is properly solved (Beckmann et al. 2018, see Tremmel et al. 2015 and Pfister et al. 2019 for how to avoid this effect), or influence the galaxy dynamics (iii). In our case, instead, we took a simpler approach and we opted for an already large initial BH mass, $\approx 12$ times larger than the DM particle mass and $\approx 60$ times larger than the gas/star particle mass. Tremmel et al. (2015) also advocate that small gravitational softening is needed, along with a sufficiently high ratio of MBH to particle mass, to properly resolve dynamics; in our simulation the extremely high spatial resolution, hence the small gravitational softening, is enough to correctly account for the dynamics throughout the evolution, except at most for the very early times immediately after BH formation.

2.2.2 BH accretion

BH accretion is implemented via the commonly adopted Bondi-Hoyle-Lyttleton formula (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952, BHL hereafter), where the accretion rate for an homogeneous medium in relative motion around the BH $\dot{M}_{\rm BHL}$ is defined as

[TABLE]

where $M_{\bullet}$ is the BH mass, $\rho_{\rm gas}$ is the gas density around the BH and $v_{\rm rel}$ is the gas-BH relative velocity. Because of the limited resolution, many previous studies, increased the BH accretion rate adopting a ‘boost’ factor $\alpha$, either constant ($\alpha=100$; Di Matteo et al., 2005) or scaling with density ($\alpha\propto\rho_{\rm gas}^{2}$; Booth & Schaye, 2009). Here, thanks to the very-high resolution achieved, we can simply use the ‘real’ BHL accretion rate. Following Bellovary et al. (2010), Choi et al. (2012), and Tremmel et al. (2017), instead of the commonly employed kernel-weighted properties around the BH, we kernel-average the individual accretion rate for the particles within the BH kernel, obtaining:

[TABLE]

For the BH, we consider a cubic spline with 32 neighbours. However, in order to prevent the BH accreting low density material very far away from it, the BH smoothing length is not allowed to exceed a maximum physical radius $h_{\rm max}=150\rm\,pc\,h^{-1}$. When $h_{\rm max}$ is reached, only the particles within it are taken into account for the accretion.

In addition, we cap the accretion rate at the Eddington limit, i.e. the maximum accretion rate for which radiation luminosity does not overcome the gravitational pull, defined as

[TABLE]

where $\eta_{\rm acc}$ is the accretion radiative efficiency, that we set to 0.1 (Soltan, 1982), $c$ is the speed of light, and $\sigma_{\rm T}=6.65\times 10^{-25}\rm\,cm^{2}$ is the Thomson scattering cross section.

2.2.3 BH feedback

Every time a BH accretes material, a fraction $\eta_{\rm acc}=0.1$ of the energy is converted into radiation, with only the remaining mass–energy actively contributing to the BH mass growth. Moreover, only part of it actually couples with the surrounding gas, heating or expelling it. Unfortunately, we do not have a unique physically-motivated value for this feedback efficiency, and different studies, depending on resolution, on the number of neighbouring particles/cells considered, and on the numerical technique employed, made different choices, with values ranging from a few $10^{-3}$ (Ostriker et al., 2010) up to 0.05 (Di Matteo et al., 2005) or 0.15 (Dubois et al., 2012). Here, we assume a moderately weak feedback, where only a fraction $\eta_{\rm fbk}=0.005$ of the energy produced is dumped on to the gas as thermal energy in a kernel-weighted fashion, among the 32 neighbours defining the kernel size around the BH. In Appendix A, we validate our BH model on an isolate Milky Way-galaxy simulation, to highlight how the choice of the accretion prescription and the feedback efficiency affect the BH growth with respect to the galaxy SFR. In addition, in Appendix B, we also test the effect of different feedback efficiencies in the cosmological run by rerunning our simulation for $\sim 100$ Myr from $z=8$ to $z=7.2$, when the BH is already massive.

2.2.4 BH mergers

When two galaxies merge, the BHs hosted in their centres spiral towards the centre under the effect of dynamical friction on DM, gas and stars, until they bind in a BH binary. Then, three-body scattering or dynamical friction on gas kick in and harden the binary, down to the regime where gravitational wave emission leads to the coalescence of the two BHs. Cosmological simulations are unable to follow the entire binary formation and coalescence process, because of the huge dynamic range involved, and have to rely on sub-grid prescriptions to model BH mergers.

The most common algorithm to treat BH binary mergers checks whether the BH separation approaches or hits the resolution limit, i.e. they are in the same cell (or in close enough cells) in mesh-based codes or they share the kernel in particle-based codes, and merges them into a new single more massive BH (e.g. Dubois et al., 2012). Sometimes, a relative velocity/binding criterion is also employed, to prevent BHs to merge when they approach each other at very high speeds (e.g. Sijacki et al., 2009; Schaye et al., 2015). In our case, we employ both a distance and binding criteria, merging the BHs only when their relative speed is lower than the escape velocity from the binary system. When these criteria are matched, the mass of the least massive BH is added to that of the primary BH, which is moved to the centre of mass of the pre-existing binary, and linear momentum conservation is enforced.

2.3 Initial conditions

According to recent results, high-redshift quasars are hosted in haloes with typical masses in the range $3\times 10^{11}-\times 10^{13}$ (Di Matteo et al., 2017; Tenneti et al., 2018b), they do not live in the most overdense regions of the Universe (Mazzucchelli et al., 2017; Uchiyama et al., 2018), and they are powered by BHs with masses above a few $10^{8}\,\rm{M_{\odot}}$ up to $10^{10}\,\rm{M_{\odot}}$.

We follow the evolution of a typical quasar host halo, with a virial mass $M_{\rm vir}\approx 3\times 10^{12}\,\rm{M_{\odot}}$ at redshift $z=6$, adopting the Planck Collaboration et al. (2016) cosmological parameters, with $\Omega_{\rm m}=0.3089$, $\Omega_{\Lambda}=0.6911$, $\Omega_{\rm b}=0.0489$, $\sigma_{8}=0.8159$, $n_{\rm s}=0.9667$, and $H_{0}=67.74\,\rm km\,s^{-1}Mpc^{-1}$, and we assume negligible contribution from both radiation and curvature.

Because of the rarity of such massive haloes, i.e. $\sim$1 in a 100 Mpc $h^{-1}$ comoving box, where $h=H_{0}/(100\rm\,km\,s^{-1}Mpc^{-1}$), a very large box is necessary to guarantee the presence of at least one of these objects (see, e.g., Barai et al., 2018). However, this can be avoided by creating a constrained realisation in a smaller box (Bertschinger, 1987; Hoffman & Ribak, 1991). Here, we consider a cosmological volume of 75 Mpc $h^{-1}$ comoving, where we impose the collapse of a halo with $\sim 10^{12}\,\rm{M_{\odot}}$ at $z=6$. The initial conditions are generated at $z=100$ by music (Hahn & Abel, 2013), with a coarse resolution grid of 2563, that gives a DM mass resolution of $\approx 3\times 10^{9}\,\rm{M_{\odot}}$, enough to resolve the target halo with about 1000 particles.

However, since we also want to enforce the overdensity criterion, a single realisation is not enough. Due to the constraint imposed to the random realization, the target halo is always the most massive halo in the simulation, and forms close to the centre of the box We have run several coarse-level DM-only simulations of the same box, with different random seeds, and we have identified all the haloes resolved by amiga halo finder (Knollmann & Knebe, 2009, ahf hereafter) with at least 50 particles. Then, we have used the stellar-to-halo mass relation by Behroozi et al. (2013) on the selected haloes to derive the rest-frame UV magnitude of the corresponding galaxies. Since at $z\sim 5-6$, the UV emission is shifted to the local $i$-band, we followed Uchiyama et al. (2018) and filtered the selected haloes with $i<25$. Finally, we have computed the galaxy overdensity significance $\sigma=(N-\bar{N})/\sigma_{\rm N}$, with $N$ the number of galaxies in the aperture, and $\bar{N}$ and $\sigma_{\rm N}$ the average and the standard deviation of $N$ in all the considered region, in apertures of $0.75$  proper Mpc ($\sim 1.8$ arcmin) around the target haloes of all realisations, and compared them with the results by Uchiyama et al. (2018). The chosen realisation corresponds to the one with $\sigma\sim 2-3$, reported in Fig. 1, with the target halo shown as a black star.

2.3.1 Refinement strategy for the high-resolution region

To refine the zoom-in region up to the desired resolution, we follow the approach reported in Fiacconi et al. (2017):

From the results of the coarse-level DM-only simulation, we identify the halo using ahf. We flag all the DM particles enclosed in a sphere centred on the target halo of radius 2.5$r_{\rm vir}$, with $r_{\rm vir}$ the halo virial radius, and we trace them back in time to the initial conditions. 2. 2.

We add one refinement level around the Lagrangian box surrounding the previously flagged DM particles, and we rerun the DM-only simulation. 3. 3.

We repeat step (i), but this time by adding two refinement levels in a convex hull region surrounding the flagged DM particles, and re-run the DM-only simulation.

This step-by-step approach guarantees that the contamination in the high-resolution region is well below 1 per cent, and is completely absent within $r_{\rm vir}$. In the zoom-in region, the mass resolution achieved is $1.53\times 10^{4}\,\rm{M_{\odot}}$ for gas and stellar particles, and $8.19\times 10^{4}\,\rm{M_{\odot}}$ for DM, with a total number of particles of 67.5M in gas and 90.2M in DM.

3 Results

We now present the main results of our simulation down to $z=7$. For all the following analyses, we identify the halo using ahf, and only select the particles within a sphere of radius $0.2r_{\rm vir}$, to avoid contamination from satellite galaxies.

3.1 Galaxy host and MBH evolution across cosmic time

Here, we discuss the evolution of the target galaxy and of its central MBH during the first 800 Myr of the Universe. We defer the discussion about the other MBHs existing in the simulation to Paper III (in preparation).

3.1.1 Halo and galaxy growth

In Fig. 2 we show the halo mass $M_{\rm halo}$ (black solid), the stellar mass $M_{\rm star}$ computed using the actual mass of the stellar particles in the simulation (cyan dashed) and gas mass $M_{\rm gas}$ (green dotted) as a function of time. We also show as a comparison the stellar mass estimate based on the UV luminosity, as derived in Song et al. (2016), as a red dot-dashed line (see Trebitsch et al., 2018, for a discussion). Despite the short time available, the typically high densities and the stronger turbulent motions in the interstellar medium result in a very rapid build up of the quasar host galaxy mass, up to $M_{\rm star}\sim 10^{11}\,\rm{M_{\odot}}$ at $z=7$. The stellar mass at $z=7$ corresponds to a conversion efficiency of about $f_{\star}=M_{\rm star}/M_{\rm halo}\approx 0.06$, broadly consistent with the extrapolation of the stellar-to-halo mass relation from Moster et al. (2018).

In the first 300 Myr (z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{>}}11), the main progenitor of the quasar host is not easily discernible from other haloes in the same region that show a comparable mass ($M_{\rm halo}\sim 2-3\times 10^{10}\,\rm{M_{\odot}}$). Then, around $z\sim 11$, a significantly massive halo emerges from the simultaneous coalescence of three progenitors, reaching $M_{\rm halo}\approx 10^{11}\,\rm{M_{\odot}}$. After this major event, a phase of steady growth begins, driven by smooth accretion or minor mergers, i.e., with mass ratio between the two interacting haloes less than 1:10.

During the first 400 Myr, most of the baryons remain in the form of gas, because of the low efficiency of SF and the powerful effect of SN feedback that heats the gas up and sweeps it away. The burstiness of the star formation history in these stages is also reflected in the discrepancy between the actual mass and the UV-based diagnostic (which is a better tracer of continuous SF histories). This trend changes as soon as the galaxy reaches $M_{\rm star}\sim 10^{10}\,\rm{M_{\odot}}$ (Dubois et al., 2015; Habouzit et al., 2017), when SNe become less effective and gas accumulates in the halo, settling in a disc-like structure. At this stage, self-regulation between gas cooling, SF, and SN feedback is able to stabilise the gas consumption and replenishing and yield a similar evolution, so that the gas-star ratio remains roughly constant. Similarly, below $z=10$, the UV-based stellar mass approaches the actual stellar mass, with only stronger fluctuations due to the evolution of the average age of the stellar population.

The quick build-up of a large stellar mass observed for these galaxies is the result of the massive inflows on to the halo and the very high SFRs shown in Fig. 3. During the first 200 Myr, multiple dips can be observed, corresponding to the ejection of gas by SN feedback. At later times, when the galaxy has become massive enough, the SFR becomes more stable with values ranging from 100 up to 800 $\rm\rm{M_{\odot}}yr^{-1}$. This is in good agreement with observations of high redshift quasars by Decarli et al. (2018, D18 hereafter), with SFRs in the range 100-2000 $\rm\rm{M_{\odot}}yr^{-1}$.

At $z\sim 7$, the galaxy exhibits a SFR that is consistent with the main sequence according to the empirical model by Behroozi et al. (2013) (extrapolated up to the stellar mass we have), see the green shaded area and the dotted green line in Fig. 3. On the other hand, these SFRs are well above the best-fit of observed main sequence galaxies (Salmon et al., 2015), suggesting that these galaxies are more likely starburst rather than normal galaxies as shown by the shaded cyan area and the dotted cyan line in Fig. 3.

Given the large number of SNe in the galaxy, the metallicity of the gas (and of the newly formed stars) rapidly increases to about solar values, as shown in Fig. 4. The stellar metallicity (red solid line) evolves smoothly with time, saturating around $Z_{\rm star}\approx 2\rm\,Z_{\odot}$, with small oscillations resulting from mergers with smaller lower-metallicity galaxies. The gas metallicity (blue dashed line) exhibits much stronger oscillations, mainly associated with massive inflow events of lower-metallicity gas from the filaments. Such large metallicities are consistent with the idea that quasar hosts at high redshift are already dust-rich, and produce a strong far-infrared emission of about $\sim 10^{12-13}\rm\,L_{\odot}$ (Venemans et al., 2017; Decarli et al., 2018).

As already discussed, after the first 400 Myr, the galaxy mainly grows via smooth accretion from the filaments around the halo, which continuously provide fresh material, as visible in Fig. 5, where we show the gaseous (left-hand panels) and stellar (right-hand panels) distribution within 400 comoving kpc from the central galaxy at $z\approx 10,8.5,$ and $7$. The maps correspond to the column density integrated along the z-direction in the 400 comoving kpc box, with the colormap spanning column densities in the range of $10^{7}\,\rm{M_{\odot}}kpc^{-2}<\Sigma_{\rm gas}<10^{11}\,\rm{M_{\odot}}kpc^{-2}$ for gas and $5\times 10^{6}\,\rm{M_{\odot}}kpc^{-2}<\Sigma_{\rm star}<5\times 10^{10}\,\rm{M_{\odot}}kpc^{-2}$ for stars. The left-hand panels show a large-scale gaseous disc extending over a few kpc, whose orientation changes with time following the angular momentum evolution of the massive gas inflows. The stellar distribution also exhibits a disc-like structure, with a denser concentration in the central region, where the MBH, shown as a red dot, is located.

3.1.2 The central MBH

In Fig. 6 we show the mass (red solid) and the average accretion rate (blue dashed) in units of the Eddington luminosity for the MBH, for $\eta_{\rm acc}=0.1$, as a function of time. Immediately after its formation, the gas density around the MBH is high enough to trigger a short accretion phase. This initial short accretion phase is a consequence of the seeding criterion employed here, which creates a new MBH at the highest gas density peak. However, after a few Myr, SN events (together with the AGN feedback) completely evacuate the region around the MBH, halting its growth. After this initial burst, the MBH experiences a long phase where accretion is suppressed by SN feedback (the galaxy stellar mass is below $\sim 10^{9-10}\rm\,\rm{M_{\odot}}$, Dubois et al., 2014c), and only later it starts to grow at the Eddington limit until it reaches a few $10^{8}\,\rm{M_{\odot}}$. At this point, because of the larger energy injection powered by accretion, the gas around the MBH is heated up to more than $10^{6}$ K and the growth slows down, with an average accretion rate oscillating between 0.5 and 0.8 times the Eddington limit, and the MBH reaches $M_{\rm BH}=8.6\times 10^{8}\,\rm{M_{\odot}}$ by $z\approx 7$. The behaviour observed for the MBH growth is consistent with several previous studies (Dubois et al., 2014c; Habouzit et al., 2017; Prieto et al., 2017; Anglés-Alcázar et al., 2017), where the MBH was found to starve until the galaxy reached a critical mass between $10^{9}\,\rm{M_{\odot}}$ and $10^{10}\,\rm{M_{\odot}}$, or a halo mass/temperature about $10^{11}-10^{12}$ and $10^{5.6}$ K (Bower et al., 2017; McAlpine et al., 2018). This result implies that our initial seed mass (see Section 2.2.1) is not crucial for the MBH evolution, which is instead regulated by the galaxy’s ability to funnel gas towards the centre. 444To further confirm this statement we have rerun our simulation with a $10^{5}\,\rm{M_{\odot}}$ initial BH mass, and evolved it down to $z=9$. The results are reported in Appendix C. whose results Finally, the lack of MBH growth in low-mass galaxies enables to reconcile different observations of the high-z MBH and quasar populations, as initially suggested in Volonteri & Stark (2011).

In Fig. 6, we also report the MBH mass estimates for the two most distant quasars known to date, J1120+0641 at $z=7.085$ (Mortlock et al., 2011) and J1342+0928 at $z=7.541$ (Bañados et al., 2018). Our result is in reasonable agreement with observations of the most distant quasars to date. Compared to a previous similar study by Smidt et al. (2018), our MBH forms later and, despite a larger seed mass, its growth is suppressed for the first 100 Myr, unlike in Smidt et al. (2018), who use a less efficient thermal SN feedback.

If we measure the accretion rate on the MBH directly computed by the code using the BHL formula before applying the cap at the Eddington limit, in many cases we get super-Eddington accretion rates, with Eddington ratios up to about 100. Since our sub-grid modelling for BH feedback does not take into account the super-critical regime where powerful jets are produced (e.g. Sa̧dowski et al., 2016), we opted for capping the accretion rate at the Eddington limit. Nevertheless, if such accretion rates are reached in massive galaxies at high redshift, even for a short time, the initial MBH growth could be accelerated, increasing the final mass and possibly alleviating the tight constraints on the seed BH masses.

In Fig. 7, we compare our simulation (as blue dots connected by a solid line) with local data from Reines & Volonteri (2015). The red squares are broad line AGN with single-epoch mass estimates, whereas the green stars correspond to MBHs with measurements from stellar or gas dynamics, masers and reverberation mapping. Interestingly, after the initial phase where SN feedback hampers MBH growth and its growth lags the galaxy’s, the MBH in our simulation grows along with its host, always remaining within the scatter of the observations of low-redshift MBHs. This suggests that MBHs in high-redshift galaxies are not different from their local counterparts. Analyses of observations at high redshift argued that MBHs grow faster than their galaxy host, resulting in MBHs overmassive with respect to those hosted in galaxies with the same masses at $z=0$. In Section 3.2.2, we measure the dynamical masses from our simulation, in order to directly compare them with the observed data, and we discuss this apparent discrepancy.

The relation between MBH and galaxy mass is built through how MBHs and galaxies gain their mass. The relation between MBH accretion rate and SFR is key to interpret this relation. In observations, MBH accretion rate and SFR appear to be uncorrelated on a source-by-source basis, but the ratio of the mean MBH accretion rate is about $10^{-3}$ the mean SFR (Chen et al., 2013, and references therein). For high redshift quasars, the observed ratios are consistently much higher, $\sim 10^{-2}$ (Wang et al., 2011; Netzer et al., 2014; Willott et al., 2015), or even larger, although we stress that these are not average values, and that they are obtained on sources selected for being quasars, i.e., where the MBH luminosity outshines the galaxy. In Fig. 8, we compare these results with those in our simulation, shown as filled circles. The size of the circle corresponds to the MBH mass, whereas the colour is $z$. The dashed curves, from bottom to top, correspond to $\dot{M}_{\rm BH}/\dot{M}_{\rm star}=10^{-3},10^{-2},$ and $10^{-1}$. At very high redshift, when the MBH is not growing much, the accretion rate is low, similar to, or lower than, $z\sim 1-2$ galaxies (Mullaney et al., 2012). At later times, after the MBH has entered the strong accretion regime, the data points move upwards, close to high redshift quasar measurements, although the observations are at somewhat lower redshift ($z=6$ and $z=4.8$) (Willott et al., 2015; Wang et al., 2011; Netzer et al., 2014). If we assume a baseline of about $10^{-3}$ for the average ratio of BH accretion rate to SFR, the interpretation is therefore that initially the MBH lags behind its galaxy ($z=10-11$), then it catches up and even surpasses it ($z<9$). Eventually, when BH feedback becomes important (at $t=700$ Myr, $z=7.5$), BH accretion slows down.

Because of the large luminosity associated to the central MBH in the optical/UV band, that exceeds by far their stellar luminosity, quasar hosts at high redshift can currently only be detected via their molecular and dust component in the sub-millimetre band. In Fig. 9, we show the evolution of the BH bolometric luminosity relative to the rest-frame UV luminosity of the galaxy as a blue line. The far-UV flux is computed from the up-to-date stellar population synthesis models by Bruzual & Charlot (2003) also employed for the stellar radiation sub-grid model. The UV luminosity has not been corrected for dust extinction, which is expected to be important, given the high metallicity of the gas (see also Tenneti et al., 2018a). The dashed region corresponds to the early phase when the MBH was not growing, hence only the galaxy would have been observable. For $z\lesssim 10$, the MBH has accreted enough mass to become more luminous than its host, i.e. a quasar. At $z\sim 7$, its luminosity has become almost a hundred times larger than that of the galaxy, completely dominating the emission. The advent of facilities that can observe the rest-frame optical/near-IR for high-redshift quasars with high angular resolution that permit good subtraction of the point spread function, like JWST, will soon allow us to measure the stellar properties of the host.

3.2 Quasar hosts at $z\lesssim 7$

In this section we focus on the properties of the MBH and its host at $z\sim 7$, and specifically how observational tracers and diagnostics compare to quantities directly measured in the simulation. Our aim is to aid interpretation of observations and identify possible biases caused by selection effects or limited spatial resolution.

3.2.1 Gas and stellar distribution

Thanks to the very high resolution of the simulation, we can properly resolve the gas/star distribution in the galaxy. At $z\sim 7$, metal rich gas can quickly cool within the halo, and it settles in a well defined dense disc extending out to $\sim 1$ kpc. Continuous inflows and SN explosions within the system create a massive reservoir of mildly lower density gas extending out to several kpc with filaments and clumps (see top panels in Fig. 10, where we report the total gas surface density face-on and edge-on). The dense disc is relatively thin, and shows a well defined spiral structure, also observable in the stellar component (whose far-UV flux is reported in the bottom panels). A low density region can also be seen close to the central MBH (red dot), due to the effect of the BH feedback. However, because of the high gas densities around the MBH, BH feedback is not effective enough to affect the gas at larger distances. Around the central stellar disc, many stellar clusters and some small galaxies can be distinctly seen in the far-UV maps. Although not visible in the young stellar component, a central bulge is already present in the galaxy. Following Vogelsberger et al. (2014), we measure the B/T ratio from the circularity distribution of the stars, that results in $\sim 0.45$, suggesting that an already massive bulge has formed in the centre of the system.

In the middle panels of Fig. 10, we show the expected luminosity from [CII] emission, assuming the [CII] emissivity defined in Pallottini et al. (2017) as

[TABLE]

with $n_{\rm gas}=\rho_{\rm gas}/\rm m_{p}$. As expected, [CII] traces very well the star-forming part of the galaxy, corresponding to the dense disc in the central kpc, and is present at larger distances only in the dense filaments or clumps flowing on to the galaxy or ejected from it. By integrating the [CII] flux in the disc, we obtain $L_{\rm[CII]}\sim\times 10^{9.5}\rm\,L_{\odot}$, in good agreement with the observed luminosities by D18. We also evaluated $L_{\rm[CII]}$ at different redshifts, and we found a reasonably good agreement with the $L_{\rm[CII]}$-SFR relation employed by D18, but with significant scatter, up to half a decade.

3.2.2 Gas/stellar kinematics and dynamical masses

At $z\sim 7$ and above, we are still unable to properly resolve the galaxy structure in observations, and furthermore for quasars the stellar distribution is hidden by the quasar emission in UV-based observations. This limitation results in the impossibility of accurately measuring the galaxy stellar masses. Dynamical masses are estimated from the gas properties in sub-mm observation, from the emitting region size and the line width by assuming virial equilibrium and either a rotationally supported disc or a dispersion dominated system.

To consistently compare simulation results with the current observations, and to infer whether the dynamical mass estimates are good tracers of the entire system masses, we then need to mimic the observational approach in our numerical simulations. We summarise here the procedure we follow:

we estimate the galaxy angular momentum and we rotate the system by 55 degrees perpendicular to it; 2. 2.

we produce 2-dimensional average maps of the line-of-sight velocity $v_{\rm los}$ and the velocity dispersion $\sigma_{v,\rm los}$, weighted by the gas density (we limit to $\rho_{\rm gas}>1.67\times 10^{-24}\rm\,g\,cm^{-3}$ only to avoid considering low-density hot gas) (a) and the [CII] line intensity (b) at a resolution of 50 pc, the target resolution for observations to be able to resolve the Bondi region around the MBH; 3. 3.

we remove the low-density/luminosity outer regions by reducing the opacity of the map according to the gas column density (a) or the [CII] line intensity (b), in logarithmic steps over three orders of magnitude; 4. 4.

we degrade the spatial resolution of 50 pc by convolving the maps (and the opacity mask) with a Gaussian kernel with full width half maximum (FWHM) 200 pc and 2 kpc, respectively; 5. 5.

we compute the 1-dimension line profile from each map and estimate the width of the line; if the profile can be fitted by a Gaussian function, we estimate the FWHM from the best-fit, otherwise we assume the FWHM corresponds to the distance between the two peaks of the profile; for the velocity dispersion profile, we fit using a Gaussian function and take the mean value for $\sigma_{\rm fit}$; 6. 6.

we compute the emitting region size $R_{x}$ from the opacity mask, converting the 2-dimension image into a 1-dimension radial profile and taking the semi-major axis of the fitted Gaussian profile; 7. 7.

we estimate the dynamical mass from each channel (gas or [CII] luminosity) according to the virial theorem, as D18

[TABLE]

where the $x$ subscript identifies the tracer employed (gas density or [CII] emission), $\rm G$ is the gravitational constant, $\sigma_{\rm FWHM}=\rm FWHM/2.35$ and $i$ is the inclination angle (55 degrees in the case of rotationally-supported systems). Since the gas in the galaxy is distributed mainly in a disc, the best estimate for the dynamical mass is that of the rotationally-supported case, but for sake of completeness, we also compute the expected dynamical mass assuming the FWHM as a measure of the dispersion of the system.

As an example, we show in Fig. 11 the velocity maps and profiles obtained at $z\approx 7$ from the [CII] line intensity at resolutions of 50 pc (left column), 200 pc (middle column), and 2 kpc (right column), the latter two representing standard resolutions for ALMA observations of high-z quasars. In the top row we show the opacity mask, in the second row the line-of-sight velocity maps, and in the third row the line-of-sight velocity dispersion. In the bottom two rows we show the line-of-sight velocity and velocity dispersion profiles extracted from the [CII] maps, linearly sampled in bins 30 $\rm km\,s^{-1}$ wide. As we degrade the resolution, the velocity pattern of the spiral arms gets smoothed, resulting in a decrease of the peak velocity of a factor of a few relative to the original value. Similarly, at lower resolution the measured gas velocity dispersion is also suppressed. However, the fits to the profiles show that the FWHM of the Gaussian fit to the line-of-sight velocity is not significantly affected by resolution, although the approximation of a Gaussian distribution improves at lower resolution, while a double-peak is more clearly observed when the distribution is properly resolved. The mean of the velocity dispersion profile instead slightly rises with decreasing resolution, due to the increasing importance of the large-velocity tail resulting from the smoothing procedure.

In table 1 we report the relevant quantities extracted from the velocity maps of the target galaxy at $z\approx 7$ and compare them to the quantities measured directly from the simulation.

Regarding gas kinematics, we find that rotational support is dominant. If we measure the actual rotation curve and velocity dispersion from the gas particle distribution, the velocity dispersion in the central kpc is about 35 per cent of the rotational velocity, and it is consistent with the values inferred from the velocity maps. We repeat the kinematical analysis detailed above on the stellar distribution: for the three resolutions, the FWHM is 430, 390, and 105 $\rm km~{}s^{-1}$, and $\sigma_{\rm fit,*}$ (measured as described in point (v) above, but for stars) is 350, 330, and 300 $\rm km~{}s^{-1}$, much larger than that of the gaseous disc, because of the presence of a central bulge and the outer stellar halo, that are dispersion-dominated. The “true” stellar velocity dispersion, measured on the stellar particles, $\sigma_{*}$, is 350 $\rm km~{}s^{-1}$, while the virial circular velocity of the halo $v_{\rm circ,vir}\approx 385\rm\,km~{}s^{-1}$. In summary, $\sigma_{\rm FWHM}$ from the gas maps underestimates the rotational velocity of the gas, the velocity dispersion of the stellar components as well as the circular velocity of the halo. The estimate of $\sigma_{\rm fit,*}$ from the stellar maps is instead consistent with the “true" stellar velocity dispersion, and it is also consistent with the circular velocity of the halo.

Regarding dynamical masses, our estimates of about $2~{}-~{}5\times~{}10^{10}\,\rm{M_{\odot}}$ from the velocity maps are in good agreement with the observations by D18, but they underestimate the actual mass in the system (within the observed radius) by about a factor of 2 to 5, for rotationally-supported systems, with an even larger discrepancy for dispersion-dominated systems. The total mass within the emitting radius is $\sim 1.2\times 10^{11}\,\rm{M_{\odot}}$, mostly in stars, with gas contributing to about 25 per cent.

3.2.3 BH-galaxy correlations

Benchmark correlations in local galaxies are based on bulge stellar masses and stellar velocity dispersions (Kormendy & Ho, 2013, and references therein). The observations of high redshift quasars, instead, provide dynamical masses and gas-based velocity dispersions, which, as detailed in the previous section, may differ from (bulge) stellar masses and stellar velocity dispersions. We therefore perform two comparisons: with observations at $z>6$, using the high-redshift observables and with observations at $z=0$, using the low-redshift observables. This allows us on the one hand to verify that the simulation agrees with the observed properties of high-redshift quasars and validate our results, and on the other hand to compare correctly to the $z=0$ correlations.

MBHs at these redshift are commonly considered overmassive compared to their local counterparts, based on the ratio between MBH and dynamical mass being much above the “canonical” value of $10^{-3}$. Although these mass ratios are well above the local values for low-mass MBHs, that settle around $\sim 10^{-4}-10^{-3}$, they are consistent with those around $\sim 10^{-2}$ for high-mass MBHs (see Fig. 18 in Kormendy & Ho, 2013 and Fig. 11 in Reines & Volonteri, 2015), as shown in Fig. 12, where we report the stellar (or dynamical) mass to MBH mass ratio as a function of MBH mass. This result is also consistent with Lyu et al. (2016), who suggested that the BH-galaxy correlations could already be in place at z\mathrel{\hbox to0.0pt{\lower 3.0pt\hbox{\sim}\hss}\raise 2.0pt\hbox{>}}5.

As in Fig. 7, we report the data by Reines & Volonteri (2015) as red squares and green stars. Since at $z\sim 7$ bulge masses cannot be properly estimated in observations because of the limited angular resolution, even if there was no bright quasar, we compare to relations in local galaxies that use the total stellar mass, as in section 3.1.2. High-mass BHs, above $\sim 10^{8}\,\rm{M_{\odot}}$, at $z=0$ are hosted in elliptical galaxies, and therefore the difference between total and bulge mass is negligible in that mass range. We also include the measures by D18 as cyan plus signs (rotationally-supported) and orange crosses (dispersion-dominated). Our dynamical mass measurements are shown as black (rotationally-supported) and magenta circles (dispersion-dominated); using the dynamical mass estimated through the gas tracers, the BH mass - dynamical mass ratio at $z\approx 7$ is $\sim 0.017$, in agreement with D18. The actual ratio, for the stellar mass in the galaxy, is shown with a blue diamond.

The data points by D18, as well as our simulation results, are perfectly consistent with those of nearby massive ellipticals, without any sign of clear deviation. The observations would be slightly above the local data, but still consistent with them, if the galaxies were dispersion-dominated. Furthermore, for the simulation, the estimate from the dynamical mass is above the “true” value, based on the stellar mass. Since in the simulation, analysed in the same way observations are analysed, the dynamical mass underestimates the stellar mass in the system, and therefore even more the total mass, the inference is that MBHs at high redshift are in reality not different from their local counterparts, and they only represent the upper end of the $M_{\rm BH}-M_{\rm star}$ correlation, where massive elliptical galaxies lie today. This result is also in agreement with the semi-analytic model results by Valiante et al. (2011, 2014), who also suggested that dynamical masses could underestimate the actual galaxy mass.

In Fig. 13, we show the $M_{\rm BH}$-$\sigma$ correlation. Observations cannot access the stellar velocity dispersion in high-redshift quasars. The velocity dispersion is estimated through the Gaussian fit to the line profile of the gas component, although it is a priori correct only for dispersion-dominated galaxies. CO-detected quasars at $1.5<z<5$ from Shields et al. (2006) and Coppin et al. (2008) are shown as blue stars, high-redshift $z\sim 6$ quasars from Wang et al. (2010) as filled and empty green squares, with the empty squares based on the $\sigma$ estimate method by Ho (2007), and the D18 quasars as cyan crosses. For all these samples, $\sigma\approx\rm FWHM/2.35$ (Nelson, 2000; Shields et al., 2006). The results of the simulation, also based on $\sigma_{\rm FWHM}\equiv\rm FWHM/2.35$ as velocity dispersion (see Table 1), i.e., using the analogue of the quantity used for the high-redshift quasars, is shown with empty red diamonds. Independent of the method used to estimate $\sigma$ through gas tracers, the simulation agrees with observations, showing a shift with respect to the benchmark correlation in local galaxies, shown with black dashed line, $\log(M_{\rm BH}/\rm\rm{M_{\odot}})=8.13+4.02\log(\sigma/200\rm\,km\,s^{-1})$ (Tremaine et al., 2002) and a gray area for the 0.3 decades nominal scatter.

However, for local galaxies, shown as black dots, $\sigma$ is the stellar velocity dispersion, $\sigma_{*}$. The correct comparison to assess whether MBHs at high redshift are on the correlation is therefore with $\sigma_{\rm fit,*}$ derived from the stellar velocity dispersion maps (filled red diamonds) and $\sigma_{*}$ directly derived from the stellar particles (filled magenta diamond). As noted in Section 3.2.2, the velocity dispersion obtained from the gas maps underestimates the stellar velocity dispersion. When we use $\sigma_{*}$ as $\sigma$, shown as filled diamonds in Fig. 13, the position occupied by the simulation’s MBH is in very good agreement with local observations, which also adopt the stellar velocity dispersion as a tracer.

In summary, the MBH-galaxy relations appear to be consistent with $z=0$ values when the analogs of $z=0$ quantities are measured, and MBHs appear to be overmassive only because of the different tracers used. This result is in agreement with simulations having lower resolution but larger statistics (Huang et al., 2018), which also find relations fully consistent with local relations. To this, one has to add the bias caused by selection of quasars with respect to quiescent galaxies (Lauer et al., 2007). Since the galaxy mass function is steep, at fixed quasar luminosity it is more probable to select a MBH more massive than the average in a galaxy of modest mass than a MBH with average mass for its host in a much rarer high-mass galaxy. For examples of the consequences of this bias, see Volonteri & Stark (2011) and Volonteri & Reines (2016).

4 Discussion and conclusions

We have investigated the evolution of a high-redshift quasar host with an extremely high resolution cosmological zoom-in simulation performed with the code gizmo (Hopkins, 2015). In order to ensure that the target halo would be consistent with up-to-date constraints from both simulations (Di Matteo et al., 2017; Tenneti et al., 2018b, a) and observations (Mazzucchelli et al., 2017; Uchiyama et al., 2018), we carefully selected the initial conditions to reproduce the desired halo properties and environmental conditions. In our simulation, we employed a state-of-the-art physically-motivated sub-grid modelling, including non-equilibrium chemistry for primordial elements, turbulence-based SF and stellar feedback (mechanical feedback from type II/Ia SNe and radiation from young stars), as in Lupi (2019). Thanks to the extremely high resolution achieved (down to 5 pc for the gas), we were able to follow in great detail the evolution of the quasar host, and its central MBH, and accurately compare the results at $z\approx 7$ with recent observations of high-redshift quasars by Decarli et al. (2018).

Two main caveats are present in the study. First, the lack of additional processes like HII regions, radiation pressure, stellar winds could result in a slightly less effective suppression of SF, as already shown by Rosdahl et al. (2017); Lupi (2019). Our choice of a fiducial boost factor to the radial momentum by SNe was motivated by this limitation, and also by the recent results by Gentry et al. (2017, 2018), and Semenov et al. (2018), but it does not necessarily represent the “true” solution. Second, the MBH feedback model and the feedback efficiency employed could result in a slightly faster growth of the MBH and a subsequent less effective feedback on to the galaxy, that again could lead to a moderately higher stellar mass. Nonetheless, we do not expect these two limitations to significantly affect our conclusions: while a stronger SN/AGN feedback could reduce the stellar mass in the system, the resulting suppressed growth of the MBH due to a stronger AGN feedback could plausibly counterbalance this decrease, resulting in small differences compared to our results. Finally, our conclusions are based on one simulation of a particular halo, and additional simulations are needed before generalizing this result.

Our simulation showed that quasar hosts very quickly build up their stellar mass, exceeding $M_{\rm star}\sim 10^{10}\,\rm{M_{\odot}}$ already at $z=8$, in agreement with previous results (Smidt et al., 2018; Barai et al., 2018). To reach these conditions, very high SFRs are required, in very good agreement with observations, and these SFRs also result in a great amount of metals produced (about twice the solar value), and a strong dust emission. Despite gas being very turbulent at these high redshift, such massive galaxies are able to settle in a disc-like structure with well defined spiral arms, and the system we have studied is more likely rotationally-supported rather than dispersion-dominated. To assess the presence of a bulge, we computed the disc to total ratio, D/T, from the circularity of the stellar particles (Vogelsberger et al., 2014; Tenneti et al., 2016; Di Matteo et al., 2017). We placed the galaxy face-on, determined the specific angular momentum $j_{z}$ of the star particles along the z-direction and computed the circularity $c=j_{z}/j_{\rm circ}$, with $j_{\rm circ}$ the specific circular angular momentum. The disc component is then assumed to be composed of all the stellar particles with $c>0.7$ (Di Matteo et al., 2017). This procedure produces D/T$\sim 0.55$, the galaxy, therefore presents a distinct bulge. Di Matteo et al. (2017) suggest a predominance of spheroids among the hosts of the most massive black holes, with two of the hosts having D/T$\sim 0.2$ and two having $\sim 0.55$; our MBH, hosted in a massive spiral galaxy with an important central bulge, has in fact a similar D/T to BH3 and BH4 in Di Matteo et al. (2017), but our disc structure appears better resolved, likely because of the higher resolution in our simulation, leading to a better treatment of the disc vertical structure.

The large amount of gas available is also crucial for the growth of the central MBH, that accretes at the Eddington limit for most of the time, until its feedback becomes powerful enough to overcome the accretion, as already found by Di Matteo et al. (2017). However, a crucial parameter for the growth of the MBH is the galaxy mass, as already pointed out by Dubois et al. (e.g. 2014a), because of the impact of SN feedback. Indeed, in our simulation, MBH growth is strongly suppressed until the galaxy stellar mass reaches about $M_{\rm star}\sim 10^{10}\,\rm{M_{\odot}}$, in good agreement with expectations. For instance, our result agrees with that of Barai et al. (2018), that employ an effective kinetic SN feedback, but is very different from Smidt et al. (2018), where the weaker thermal SN feedback probably allows the MBH to efficiently grow already at very early times (see, e.g., Rosdahl et al. 2017 and Habouzit et al. 2017 for discussions about the relative impact of different SN feedback models on the galaxy and the central MBHs).

Although not explicitly considered here, we also found that the central MBH could have experienced super-Eddington accretion phases during its life (see, also, Pezzulli et al., 2016, where they found similar results with semi-analytic models). If these phases really occur in the Universe, we could expect a more massive MBH, that would most likely be overmassive. However, this faster growth would also strongly depend on how the feedback from a super-Eddington accreting MBH would couple with the surrounding gas, as recently shown by Regan et al. (2018).

We have examined the relative growth of the stellar component, via SFR, and of the MBH, via its accretion rate. If we consider a “canonical” value of $\sim 10^{-3}$ for keeping the relative growth at the pace needed for a symbiotic evolution maintaining the typical ratio between MBH and bulge mass, we find that at first the galaxy grows faster (ratio $<10^{-3}$), while later the MBH grows faster (ratio $>10^{-3}$). Towards the end of the simulation, once MBH feedback becomes important, the system is positioned close to the $<10^{-3}$ value. This suggests a three-phase scenario. At first, when the galaxy is small, SN feedback regulates SF in the galaxy and suppresses MBH accretion. At intermediate galaxy mass SN feedback becomes inefficient at suppressing MBH growth, while AGN feedback cannot stop a runaway phase of MBH accretion. Finally, when the MBH has become sufficiently massive, it is its own feedback that limits growth.

We have compared to observed MBH-galaxy correlations, local and at high redshift. We showed that gas-derived dynamical mass estimates typically underestimate the actual mass of the system, possibly explaining why these objects seem to be “overmassive” with respect to the host mass at high redshift, although the intrinsic relation is in good agreement with local measurements. Similarly, when considering the relation with velocity dispersion, gas-based estimates underestimate the stellar velocity dispersion, leading to the MBHs being overmassive at a given velocity dispersion with respect to the local population. In our simulation, when the “true” quantities are used, the discrepancy in the MBH-galaxy correlations completely disappears, simply leaving a population of “normal” MBHs in massive (quickly evolving) galaxy.

In summary, our main results are the following:

• •

The host galaxy quickly builds up a stellar mass reaching $10^{11}\,\rm{M_{\odot}}$ at $z=7$ (Fig. 2). The galaxy itself is a very large galaxy at these redshifts.

• •

Initially, the BH growth lags behind the galaxy growth, in relative terms, while later the BH grows faster than the galaxy, also in relative terms (Fig. 8).

• •

The BH-stellar mass evolution can be mapped on the sample of local MBHs, showing a similar behaviour as a function of galaxy mass (Fig. 7)

• •

Using the gas-based tracers for galaxy dynamical mass and velocity dispersion currently available for high-redshift quasars, MBHs appear overmassive with respect to the $z=0$ correlations (Fig. 12 and Fig. 13).

• •

Using the stellar-based tracers for galaxy stellar mass and velocity dispersion used to derive the $z=0$ correlation the tension disappears, and no difference is found between the high-redshift and the $z=0$ MBHs (Fig. 12 and Fig. 13).

Because current telescopes do not have high angular resolution at infrared wavelenghts (restframe optical), the actual stellar masses and velocity dispersion cannot be properly measured. JWST will improve our capabilities, but angular resolution will remain limited to a few hundred parsec. Gas-based estimates of the velocity dispersion are a viable alternative, but they can easily lead to an underestimation of the actual velocity dispersion or dynamical mass: our results can be used to optimize how to obtain better calibrated correction factors.

Acknowledgements

We thank the anonymous referee for constructive comments that helped us to improve the quality of the paper. We acknowledge support from the European Research Council (Project No. 267117, ‘DARK’, AL, JS; Project no. 614199, ‘BLACK’, AL, MV). SB thanks for funding through Fondecyt Iniciacion 11170268, and PCI Redes Internacionales project number REDI170093. SB also thanks for funding via Conicyt PIA ACT172033 and BASAL Centro de Astrofısica y Tecnologıas Afines (CATA) PFB-06/2007. This work was granted access to the High Performance Computing resources of CINES under the allocations A0020406955, and A0040406955 by GENCI, and it has made use of the Horizon Cluster, hosted by Institut d’Astrophysique de Paris, for the analysis of the simulation results. The maps reported in this work have been created using pynbody (Pontzen et al., 2013).

Appendix A Validation of the black hole accretion and feedback model

The limited resolution in large-scale simulation does not allow a proper modelling of the accretion process on to MBHs. In many cases, even the influence radius of the MBH is not properly resolved. Therefore, to model gas accretion and the MBH growth, different sub-grid prescriptions have been proposed in the literature. These prescriptions exhibit different sensitivity to the mass/spatial resolution achieved and to the thermodynamic modelling of the gas.

When MBHs accrete mass, a fraction of the mass-energy is released in the form of radiation, with the accretion-powered luminosity $L_{\rm BH}=\eta_{\rm rad}\dot{M}_{\rm BH}c^{2}$, where $\eta_{\rm rad}$ is the radiative efficiency. The exact value of $\eta_{\rm rad}$ depends on the radius of the innermost stable circular orbit, which is tied to the MBH spin $a$, and ranges from $\sim 0.057$ (for a Schwarzchild MBH, i.e. $a=0$) up to $\sim 0.42$ for a maximally rotating Kerr MBH ($a=1.0$), but the commonly assumed value, that we also employ here, is 0.1 (Soltan, 1982). However, only a fraction of the produced radiation actually couples with the surrounding gas, resulting in a feedback effect from the MBH. In numerical simulations, this coupling efficiency $\eta_{\rm fbk}$ is sensitive to the resolution, the numerical technique employed and also depends on the way feedback is modelled on unresolved scales. The common approach is to employ ad-hoc values calibrated to reproduce the local galaxy-MBH correlations. The coupling efficiencies thus calibrated range from a few $10^{-3}$ up to 0.15 (e.g. Di Matteo et al., 2005; Booth & Schaye, 2009; Ostriker et al., 2010). Obviously, different choices result in different accretion histories and different feedback effects on the galaxy host.

The accretion prescription employed in this study has been chosen among several others after a detailed investigation in idealised conditions. In particular, we tested three different prescriptions on an idealised Milky-Way like galaxy with a MBH with $M_{\rm BH}=4.6\times 10^{6}\,\rm{M_{\odot}}$, evolved for 1 Gyr using the same sub-grid model employed in the cosmological run. The initial conditions are those of the AGORA isolated disc comparison (Kim et al., 2016).

In particular, we tested the following models:

Common BHL accretion, based on the kernel-weighted quantities around the MBH, where

[TABLE]

where $M_{\bullet}$ is the BH mass, $\langle\rho_{\rm gas}\rangle$ is the average gas density, $\langle v_{\rm rel}\rangle$ the average gas-BH relative velocity, and $\langle c_{\rm s}\rangle$ the average sound speed. 2. 2.

Weighted BHL accretion, i.e. our fiducial model, where the accretion rate is the kernel average of the local accretion rate from each gas neighbour (Choi et al., 2012), where

[TABLE] 3. 3.

Flux accretion, where the accretion rate comes from the continuity equation as the mass flux within the MBH accretion radius, i.e.

[TABLE]

For the BH feedback, we explore different values for $\eta_{\rm fbk}$, respectively 0.5, 2.0, and 5.0 per cent, to assess how sensitive the different models were to this choice.

All simulations produce galaxies that remain on the correlation between MBH and galaxy mass. We have therefore explored additional ways to test the model, by looking at the relation between MBH accretion rate and SFR. The results are shown in Fig. 14, where we report the ratio between MBH accretion rate and SFR for the different models, averaged over the last 500 Myr of each run (checking that different choices for the period do not affect our conclusions), and compare it with the observational constraints (grey shaded area) determined from Mullaney et al. (2012) (upper edge) and Chen et al. (2013) (lower edge). The red diamonds correspond to model (i), blue squares to model (ii), and green stars to model (iii). These results show that our fiducial model, applied on a local galaxy, produces a better agreement with local correlations of the $\dot{M}_{\rm BH}-\dot{M}_{\rm star}$ ratio.

Furthermore, as demonstrated by Negri & Volonteri (2017), the feedback efficiency is dependent on resolution, with higher resolutions requiring a lower efficiency. This is due to the fact that, for a fixed number of elements affected by BH feedback and energy injected, heating is more effective at higher resolution, because the mass to be heated and swept up is smaller.

This can be also confirmed via simple analytical arguments. Let’s assume that we have two simulations with different resolutions (a ’low’ resolution case and a ’high’ resolution case). We consider that both simulations have a MBH with mass $M_{\rm BH}$ surrounded by gas with the same properties, so that they measure the same $\dot{M}_{\rm BH}$. Assuming the BHL prescription, this ensures that the sound speed $c_{s}$ and the gas density $\rho$ are fixed for both cases.

In the two simulations, the accretion time-step can be defined by the Courant condition, i.e. $\Delta t=H/c_{S}$, with $H=\left(\frac{3N_{\rm ngb}m_{\rm gas}}{4\pi\rho}\right)^{1/3}$ the kernel size, with $N_{\rm ngb}$ the number of encompassed neighbours and $m_{\rm gas}$ the gas particle mass. The only parameters that change in the simulations are $N_{\rm ngb}$ and $m_{\rm gas}$.

In a single time-step, the accretion powered energy injected can be written as

[TABLE]

Now, assuming for simplicity that every neighbour in the kernel receives the same amount of energy (let’s keep in mind that the injection is usually done via a kernel-weighting procedure, which gives more weight to the closest particles), the average temperature increase is

[TABLE]

where $\zeta=\left(4/3\pi\rho\right)^{-1/3}[\eta_{\rm rad}\dot{M}_{\rm BH}c^{2}(\gamma-1)m_{\rm p}]/[c_{s}k_{\rm B}]$ is constant for both resolutions.

If we now compute the temperature variation ratio between low- and high-resolution simulations, we get

[TABLE]

where ‘low’ and ‘high’ identify the two resolutions.

If we require the same temperature increase observed in our isolated-galaxy, where $m_{\rm gas,low}=8.6\times 10^{4}\,\rm{M_{\odot}}$ and $0.005\lesssim\eta_{\rm fbk}\lesssim 0.02$, in the cosmological run ($m_{\rm gas,high}\sim 1.5\times 10^{4}\,\rm{M_{\odot}}$), then we should set $0.001\lesssim\eta_{\rm fbk,high}\lesssim 0.006$, consistent with our choice. This interval is also consistent with that proposed by Choi et al. (2012) for kinetic feedback models.

In general, if we compare our choice ($m_{\rm gas}\sim 1.5\times 10^{4}\rm M_{\odot}$, $N_{\rm ngb}\sim 32$, and $\eta_{\rm fbk}=0.005$) to other cosmological simulations (e.g. Schaye et al., 2015; Di Matteo et al., 2017), where $m_{\rm gas}\sim 2-3\times 10^{6}\rm\,M_{\odot}$, $N_{\rm ngb}\sim 64$ and $\eta_{\rm fbk}=0.05$ were employed, we get

[TABLE]

Since $r_{\Delta T}<1$, for fixed $M_{\rm BH}$, $c_{s}$, and $\rho$, our feedback typically heats up the gas more than that in low-resolution simulations. Obviously, the combination of the idealised simulations and the simple analytical arguments do not assure that the value we employed is the right one to reproduce the local correlations in a full cosmological simulation, but unfortunately a full box simulation evolved down to $z=0$ with the resolution we achieve is currently unfeasible.

Appendix B The effect of the BH coupling efficiency in the cosmological run

In our simulation, we employ a low efficiency value $\eta_{\rm fbk}=5\times 10^{-3}$, finding a good agreement with observations of high-redshift quasars. To assess how the star formation history of the galaxy and accretion history of the MBH would change if a more powerful feedback was used, we restarted our simulation at $z=8$ and evolved it for about 100 Myr with an efficiency four times higher, i.e. $\eta_{\rm fbk}=0.02$, as in Tremmel et al. (2015). This change in the injected energy results in a moderately larger low-density region around the MBH only, but its impact on the galaxy host as a whole is small, as shown in the top panel of Fig. 15, where we show the evolution of the galaxy host stellar mass in the two cases (the red dot-dashed line is the fiducial run and the blue dashed one the higher-efficiency case). The black solid line corresponds to the redshift range with only the fiducial case, with the cyan dotted vertical line identifying the initial redshift of the comparison. After 100 Myr, the difference in the stellar mass is about 10 per cent. In the bottom panel, we show instead the MBH mass in the two cases (with the same color/line-style scheme). Because of the stronger feedback, gas becomes hotter, and the growth is suppressed. However, during the first 45-50 Myr (about an $e$-folding time), when the MBH is less massive and the Eddington limit lower, only a small difference in accretion is observed, with the mass still increasing by about a factor of 2 and a difference of about 20 per cent. This is due to the balance between the energy injected by the MBH and the very high densities around it, that power strong radiative cooling, moderately inhibiting the MBH feedback effect. At later times, instead, the effect of the stronger feedback builds-up, keeping the gas hotter and quickly dropping the accretion rate to about 0.1 times the Eddington limit. Nevertheless, the difference in the MBH mass after 100 Myr is about a factor of two. According to these results, we can expect that, because of the very high densities in high-redshift quasar hosts, a higher $\eta_{\rm fbk}$ would only play a role at late times, between $z=6$ and $z=8$, when the MBH has already exceeded $M_{\rm BH}\sim 10^{8}\,\rm{M_{\odot}}$. The main growth phase, instead, would still be at the Eddington limit.

Appendix C The role of the initial BH mass

In order to test the impact of the initial BH mass ($M_{\rm BH,0}$) on the initial BH growth, we have performed an additional run with a seed mass of $10^{5}\,\rm{M_{\odot}}$ down to $z\sim 10.5$ and compared it with our fiducial case. In Fig. 16, we report the BH masses and accretion rates for the two runs. Immediately after formation, the lower mass seed accretes at a similar pace relative to the fiducial case. However, the situation changes when the BH in the fiducial case stops growing after a few tens of Myr (because of the combined effect of the SN feedback and the AGN heating), and the low-mass seed grows almost undisturbed. This confirms our conclusion that the initial accretion phases are mainly determined by the galaxy’s ability to funnel gas towards the BH, and not by the initial BH mass. At later times, when the galaxy mass has exceeded a critical mass, the growth would proceed at the Eddington rate for both BHs, until AGN feedback could take over. In the lower-mass case, the AGN-dominated regime would be reached below $z=8$, the redshift where it occurs in the fiducial run, because of the lower BH mass, and this would likely result in similar masses observed by $z=7$.