Black hole obscuration and duty cycles mediated by AGN feedback in high redshift galaxies

TL;DR

This study uses high-resolution simulations to explore how AGN feedback and gas obscuration influence black hole growth and visibility in high-redshift galaxies, revealing that feedback limits accretion and dust obscuration affects observed luminosity functions.

Contribution

It provides new insights into the role of gas and dust in obscuring AGN and regulating black hole growth during early galaxy evolution through detailed cosmological simulations.

Findings

AGN feedback controls black hole growth and limits accretion episodes.

Dust-enriched gas obscures AGN, reducing UV/optical visibility.

Obscuration from gas outside the nucleus is significant.

Abstract

Dense gas in the centre of galaxies feeds massive black holes, but can also become a source of obscuration and limit our ability to find faint Active Galactic Nuclei (AGN). We use a high resolution cosmological radiative hydrodynamics simulation to connect the properties of the gas in the central region (a few tens of parsecs) of a high redshift galaxy to the growth of a massive black hole during the first billion years of the Universe. We find that the feedback from the AGN efficiently controls the growth of the black hole and limits the duration of the high accretion episodes by emptying the gas reservoir. As the galaxy grows in mass, the production of metals results in the presence of dust-enriched gas in the galaxy centre that can obscure highly accreting black holes enough to strongly reduce their UV/optical visibility. We also find that the gas outside the very centre of the galaxy can contribute to the total column density and obscuration at a level at least comparable to the gas in the nuclear region. We suggest that this explains the different duty cycles required to explain the masses of high redshift quasars and the observed UV/optical luminosity functions: in our case, the AGN would be observed with an X-ray luminosity above $L_X = 10^{42}$ erg/s around 30% of the time, but with UV magnitude brighter than $M_{1450}$ = -23 only 4% of the time.

Figures

Equations

$$ϵ\propto {\text{e}}^{\frac{3}{8}{\sigma }_s^2}\left[1+\mathrm{erf}\left(\frac{{\sigma }_s^2-s_{\mathrm{crit}}}{\sqrt{2{\sigma }_s^2}}\right)\right],$$

$${\dot{M}}_{\mathrm{BHL}}=4\pi G^2{M_{\bullet }}^2\frac{\bar{\rho }}{{\left({\bar{c}}_s^2+{\bar{v}}_{\mathrm{rel}}^2\right)}^{3/2}}$$

$$r_{\mathrm{sink}}=\{\begin{array}{ll}\Delta x/4& \text{if}{r}_{\mathrm{BHL}}<\Delta x/4,\\ r_{\mathrm{BHL}}& \text{if}\Delta x/4\le r_{\mathrm{BHL}}<2\Delta x,\\ 2\Delta x& \text{if}2\Delta x\le r_{\mathrm{BHL}}.\end{array}$$

$${\dot{M}}_{\mathrm{Edd}}=\frac{4\pi G{M}_{\bullet }{m}_p}{{ϵ}_r{\sigma }_{\mathrm{T}}c},$$

$$\dot{M_{\bullet }}=\min \left({\dot{M}}_{\mathrm{BHL}},{\dot{M}}_{\mathrm{Edd}}\right),$$

$$\psi \left(r_{\mathrm{cyl}}\right)=\frac{1}{2\pi \Delta x^2}\exp \left(-\frac{r_{\mathrm{cyl}}^2}{\Delta x^2}\right),$$

$${\dot{M}}_{\mathrm{J}}(r_{\mathrm{cyl}})=\frac{\psi (r_{\mathrm{cyl}})}{\Psi }{\eta }_{\mathrm{J}}{\dot{M}}_{\bullet }$$

$$q_{\mathrm{J}}(r_{\mathrm{cyl}})={\dot{M}}_{\mathrm{J}}(r_{\mathrm{cyl}})u_{\mathrm{J}},$$

$${\dot{E}}_{\mathrm{J}}(r_{\mathrm{cyl}})=\frac{q_{\mathrm{J}}(r_{\mathrm{cyl}}{)}^2}{2{\dot{M}}_{\mathrm{J}}(r_{\mathrm{cyl}})}=\frac{\psi (r_{\mathrm{cyl}})}{\Psi }{\eta }_{\mathrm{J}}{\dot{E}}_{\mathrm{AGN}}.$$

$$\dot{M_{\bullet }}=\frac{4\pi G^2{m}_p{M_{\bullet }}^2{N}_{\mathrm{H}}}{R_{\mathrm{acc}}{\left({\bar{c}}_s^2+{\bar{v}}_{\mathrm{rel}}^2\right)}^{3/2}},$$

$$f_{\mathrm{Edd}}=7.58{ϵ}_{r,0.1}{M_{\bullet }}_{,6}⟨N_{\mathrm{H},22}⟩R_{\mathrm{acc},40}^{-1}(c_{\mathrm{s},10}^2+v_{\mathrm{rel},10}^2{)}^{-3/2},$$

$$\frac{E_b}{E_{\mathrm{inj}}}$$

$$\sim 1.9\times 10^{-4}\frac{⟨N_{\mathrm{H},22}⟩}{f_{\mathrm{Edd}}}\frac{M_{\mathrm{tot}}}{M_{\bullet }}\frac{R_{\mathrm{acc},40}}{{ϵ}_{r,0.1}{ϵ}_{f,0.15}\Delta t_{0.1}}$$

$$k(L_{\mathrm{bol}})=10.83{\left(\frac{L_{\mathrm{bol}}}{10^{10}{L}_{\odot }}\right)}^{0.28}+6.08{\left(\frac{L_{\mathrm{bol}}}{10^{10}{L}_{\odot }}\right)}^{-0.020}$$

$$\log \left(1450L_{1450}\right)=\frac{\log L_{\mathrm{bol}}-4.74}{0.91}.$$

$${\tau }_d(\lambda )={\int }_{\mathrm{ray}}{n}_d(\ell ){\sigma }_d(\lambda )d\ell ,$$

Full text

††pubyear: 2018††pagerange: Black hole obscuration and duty cycles mediated by AGN feedback in high redshift galaxies–Black hole obscuration and duty cycles mediated by AGN feedback in high redshift galaxies

1 Introduction

Active Galactic Nuclei (AGN) are sources emitting over a very large range of wavelengths and, as a family, they display numerous differences. Many of these differences are ascribed to the presence of material between the engine and the observer (Hickox & Alexander, 2018, and references therein). The absorbing material absorbs a part of the flux, dimming its brightness and modifying the transmitted spectral energy distribution. Broad emission lines, produced near the black hole (BH), can only be observed if no absorbing material is present, otherwise one observes narrow lines and infrared emission produced by reprocessing of the continuum emitted by the source.

The classic unification scheme (Antonucci, 1993; Urry & Padovani, 1995; Netzer, 2015, and references therein) explains the variety of observed AGN properties with absorbing material distributed non-uniformly, possibly in an axisymmetric or toroidal structure, “near” the BH, within the central few parsec. While this picture has had great success in explaining and interpreting most observations, gas in the host galaxy can further modulate the level of obscuration. In particular, high levels of host-wide obscuration are often advocated in conjunction with galaxy mergers (Sanders et al., 1988; Hopkins et al., 2006). In the case of powerful quasars the obscured fraction seems to increases with strong far-infrared luminosity, indicating high levels of star formation, while high-luminosity AGN do not seem, necessarily, to be characterized by non-nuclear obscuration (Merloni et al., 2014; Shimizu et al., 2018).

In local AGN, Ricci et al. (2017b) find that the fraction of obscured AGN is almost independent of luminosity, but has a strong dependence on Eddington ratio, with BHs accreting at high rates being less obscured. They advocate for the obscuring gas being close to the BH (within the sphere of influence) and for radiative feedback for the AGN being responsible for clearing out this gas when the BH reaches an effective dust-dependent Eddington rate.

The fraction of obscured AGN appears to increase with redshift (Treister & Urry, 2006; Vito et al., 2014), and this evolution is more evident at high luminosity, $L_{\rm X}\sim 10^{44}-10^{45}\,\mbox{erg.s}^{-1}$ for both optically and X-ray obscured AGN (Merloni et al., 2014; Vito et al., 2018). In Vito et al. (2018) they compared the obscured fraction as a function of luminosity of AGN at $z=3-6$ to the local one and found a strong evolution at high luminosities, with an increase of the fraction from $\sim 10-20\%$ to $\sim 70-80\%$.

High-redshift galaxies are often rich in gas (e.g. Tacconi et al., 2010, 2013; Santini et al., 2014), although not necessarily very dusty or metal rich in their earliest phases of evolution. This gas is not always available for feeding the BH in the nuclear region, because of local supernova feedback (Dubois et al., 2015; Bower et al., 2017; Anglés-Alcázar et al., 2017; Habouzit et al., 2017; Trebitsch et al., 2018) especially in low mass galaxies, but once the galaxy/halo become sufficiently massive and/or compact, the central regions can retain gas and feed the BH. This gas provides the reservoir to feed the accretion disc and then the BH, but it is also a source of obscuration.

In this paper we examine the role of gas in the near and around the Bondi radius in obscuring radiation from BHs accreting in high-$z$ galaxies. The region we study is intermediate between the “torus” (in its various nuances) and the “nuclear starburst”, beyond the dust sublimation radius and around the gravitational sphere of the BH (see Merloni et al., 2014; Hickox & Alexander, 2018, for a discussion of the various physical scales). Previous numerical investigations of obscuration around BHs focused on multi-scale re-simulations of idealized galaxies (Hopkins et al., 2012, 2016); the final re-simulations had a higher spatial and mass resolution than our study, but did not include the cosmological evolution of a young galaxy. We do not fully resolve the broad line region and the putative torus in our simulations, therefore their effects should be added as a correction. In other words, we estimate a lower limit to the total amount of obscuration in young galaxies that are actively star forming without being starbursts but live in very dynamic environments.

2 Description of the simulation

We perform a zoom-in simulation using the public, multigroup radiative transfer (RT) version of the adaptive mesh refinement (AMR) code Ramses 111https://bitbucket.org/rteyssie/ramses/ (Teyssier, 2002; Rosdahl et al., 2013; Rosdahl & Teyssier, 2015). The collisionless particles (stars and dark matter) are evolved using a particle-mesh method with a cloud-in-cell interpolation. For the gas, Ramses solves the Euler equations with the second-order MUSCL scheme (van Leer, 1979) using the HLLC Riemann solver from Toro et al. (1994) and a MinMod total variation diminishing scheme to reconstruct the intercell-fluxes. We impose a Courant factor of 0.8 to define the timestep.

The RT module propagates the radiation emitted by both stars and BHs (modelled as sink particles, see Sect. 2.3) in three frequency intervals describing the H i , He i and He ii photon field. The frequencies are defined by the ionization energies of these species, so the groups are 13.60–24.59 eV, 24.59–54.42 eV and we limit the final He ii -ionizing band to the interval 54.42–1000 eV. The radiation is evolved on the AMR grid using a first-order Godunov method to solve the first two moments of the RT equation and assuming the M1 closure (Levermore, 1984; Dubroca & Feugeas, 1999) for the Eddington tensor. To maintain the computational cost of the simulation manageable, we utilize the reduced speed of light approximation as implemented in Rosdahl et al. (2013), and we choose a reduced speed of light of $\tilde{c}=0.01c$, which is reasonable provided we do not focus on the propagation of ionization fronts in the IGM. The coupling to the hydrodynamical evolution of the gas is done through the non-equilibrium thermochemistry for hydrogen and helium, assuming the on-the-spot approximation (any photon emitted after a recombination is absorbed locally). We ignore the radiation pressure from the UV field to the gas, which plays only a subdominant role at galactic scales (Rosdahl et al., 2015). Radiation is emitted by star particles (representing fully sampled stellar population) as a function of their age and metallicity following the models of Bruzual & Charlot (2003). The AGN emit radiation when the accretion rate onto BH is high. The full description of the model is given in Sect. 2.3.5.

The AMR grid is refined using a quasi-Lagrangian criterion: a cell is selected for refinement if $\rho_{\rm DM}\Delta x^{3}+(\Omega_{\rm DM}/\Omega_{b})\rho_{\rm gas}\Delta x^{3}+(\Omega_{\rm DM}/\Omega_{b})\rho_{*}\Delta x^{3}>8\ m_{\rm DM}^{\rm HR}$, where $\rho_{\rm DM}$, $\rho_{\rm gas}$ and $\rho_{*}$ are respectively the DM, gas and stellar densities in the cell, $\Omega_{\rm DM}$ and $\Omega_{b}$ respectively the cosmic DM and baryon mass density, $\Delta x$ is the cell size, and $m_{\rm DM}^{\rm HR}$ is the mass of the highest resolution DM particle. In a DM-only run, this criterion would allow refinement as soon as there are at least 8 high-resolution DM particles in a cell.

2.1 Initial conditions

In this work, we zoom on a halo that ends up with a mass of $\sim 3\times 10^{11}{\rm M}_{\odot}$ at redshift $z\sim 5.7$ embedded in a cosmological volume of $40h^{-1}$ comoving Mpc on a side. The halo undergoes a $\sim$ 4:1 merger starting around $z\sim 7.1$ which translates into a galaxy merger around $z\sim 6.3$. The halo has been selected from a DM-only run of the full cosmological volume, initially run with $512^{3}$ particles ($m_{\rm DM}^{\rm LR}\simeq 5\times 10^{7}\,{\rm M}_{\odot}$). The initial conditions for the cosmological box as well as the zoom region have been created with Music222https://bitbucket.org/ohahn/music/ (Hahn & Abel, 2011) for a flat $\Lambda$CDM cosmology consistent with the Planck results (dark energy density $\Omega_{\Lambda}=0.692$, total matter density $\Omega_{m}=0.308$, Hubble parameter $h=0.6781$ and baryon matter density $\Omega_{b}=0.048$, Planck Collaboration et al., 2016). We select the target halo in the final output with HaloMaker (Tweed et al., 2009), which uses the AdaptaHOP algorithm (Aubert et al., 2004).

The initial conditions for the halo of interest are then regenerated at an effective resolution of $4096^{3}$ elements (level $\ell=12$), and the external region is degraded by one level so that the box is covered by a grid with $256^{3}$ elements (level $\ell=8$). This translates in a mass resolution of $m_{\rm DM}^{\rm LR}\simeq 4\times 10^{8}\,{\rm M}_{\odot}$ for the lowest resolution dark matter particles, and $m_{\rm DM}^{\rm HR}\simeq 10^{5}\,{\rm M}_{\odot}$ for the high-resolution particles. For the full radiative hydrodynamics run, we then allow for refinement down to a minimum cell size of $\Delta x=40h^{-1}\mbox{Mpc}/2^{23}\simeq 7\,\mbox{pc}$. The gas in the simulation is assumed to be neutral and homogeneously metal poor, with an initial gas phase metallicity $Z=5\times 10^{-3}Z_{\odot}=10^{-4}$.

2.2 Star formation and feedback

At the resolution of our simulation, we describe the stars forming in the galaxies as star particles with mass $m_{\star}\sim 1.8\times 10^{4}\,{\rm M}_{\odot}$.

2.2.1 Star formation

We model star formation with a Schmidt-like law (Schmidt, 1959), with an approach similar to that Rasera & Teyssier (2006). This thermo-turbulent model has been previously described in details in Kimm et al. (2017); Trebitsch et al. (2017), and we use the same version as implemented in Trebitsch et al. (2018).

During one timestep $\Delta t$, we convert a gas mass $M_{\rm sf}=\epsilon\rho_{\rm gas}\Delta x^{3}\Delta t/t_{\rm ff}$ in star-forming cells, where $G$ is the gravitational constant, $\epsilon$ is the local star formation efficiency and $t_{\rm ff}=\sqrt{3\pi/32G\rho_{\rm gas}}$ is the free fall time of the gas. The exact number $N$ of star particles formed in $\Delta t$ is drawn from a Poisson distribution $P(N)=(\lambda^{N}/N!)\exp(-\lambda)$ of parameter $\lambda=M_{\rm sf}/m_{\star}$.

Cells are considered star forming when the local density $\rho_{0}$ exceeds a threshold $\rho_{\rm threshold}=1\,\mbox{cm}^{-3}$, chosen as the typical ISM density, and when the local turbulent Mach number exceeds $\mathcal{M}\geq 2$. In each of these cells (assumed to correspond to a star forming cloud), the efficiency $\epsilon$ is computed following the ‘multi-ff PN’ model of Federrath & Klessen (2012), based on Padoan & Nordlund (2011).

[TABLE]

where $\sigma_{s}=\sigma_{s}(\sigma_{\rm gas},c_{\rm s})$ characterizes the turbulent density fluctuations, $s_{\rm crit}\equiv\mathrm{ln}\left(\frac{\rho_{\rm gas,crit}}{\rho_{0}}\right)$ is the critical density above which the gas will be accreted onto stars, and $\rho_{\rm gas,crit}\propto(\sigma_{\rm gas}^{2}+c_{\rm s}^{2})\frac{\sigma_{\rm gas}^{2}}{c_{\rm s}^{2}}$.

2.2.2 Feedback from massive stars

Our simulation models feedback from massive stars through two main channels: type II supernovae (SNe) and radiative feedback resulting from photoionization heating. Supernovae explosions are modelled using the subgrid implementation first detailed in Kimm & Cen (2014) and Kimm et al. (2015). We deposit mass and momentum in the cells neighbouring the SN host cell in a single event $t_{\rm SN}=5\,\mbox{Myr}$ after each star particle is formed. The amount of momentum released depends on the local density and metallicity of each neighbouring cell in order to capture correctly the momentum transfer at all stages of the Sedov blast wave. Photoionization heating is directly included through the coupling between the ionizing photon field and the gas. Following Kimm et al. (2017), when the Strömgren sphere of a star particle is not resolved, we increase the final radial momentum from SNe due to the pre-processing of the ISM by radiation before the SN explosion, as suggested by Geen et al. (2015).

2.3 BH model

Just as in Trebitsch et al. (2018), the BH seeding as well as the subsequent growth and associated feedback are modelled following the fiducial implementation of Dubois et al. (2012). We will now briefly recall the main features of this model, and describe two new features used in this work for the BH dynamics and for the radiation produced by the AGN.

2.3.1 BH seeding

We represent SMBH using sink particles with initial mass $M_{\bullet,0}=3\times 10^{4}\,{\rm M}_{\odot}$. These sink particles are created in cells where the following criteria are met: both the gas and stellar density exceeds a threshold that we choose to be $\rho_{\rm sink}=100\,\mbox{cm}^{-3}$; the gas must be Jeans-unstable, and there should be enough gas in the cell to form the sink particle. We also impose an exclusion radius $r_{\rm excl}=40\,\mbox{kpc}$ to avoid the formation of multiple SMBH in the same galaxy. Each sink particle is then dressed with a swarm of ‘cloud’ particles equally spaced by $\Delta x/2$ on a regular grid lattice within a sphere of radius $4\Delta x$. These cloud particles provide a convenient way of probing and averaging the properties of the gas around the BH.

2.3.2 Accretion

Once a sink particle is formed, it is allowed to grow via gas accretion and BH-BH mergers. Mergers are allowed to happen when two BH are closer than $4\Delta x$ from one another if their relative velocity is lower than the escape velocity of the binary system they would form. Gas accretion is implemented following the classical Bondi-Hoyle-Lyttleton formula (Bondi, 1952)

[TABLE]

where ${M_{\bullet}}$ is the BH mass, $\bar{\rho}$ and $\bar{c}_{s}$ are respectively the average gas density and sound speed, and $\bar{v}_{\rm rel}$ the relative velocity between the BH and the surrounding gas. The bar notation denotes an averaging over the cloud particles. The quantities are averaged using a Gaussian kernel $w\propto\exp\left(-r^{2}/r_{\rm sink}^{2}\right)$, where $r_{\rm sink}$ is defined using the Bondi radius $r_{\rm BHL}=\frac{G{M_{\bullet}}}{\bar{c_{s}}^{2}+\bar{v}_{\rm rel}^{2}}$:

[TABLE]

Just as in Trebitsch et al. (2018), we do not use any artificial boost for the gas accretion onto the BH. The accretion rate is capped at the Eddington rate:

[TABLE]

where $m_{p}$ is the proton mass, $\sigma_{\rm T}$ is the Thompson cross section, $c$ is the speed of light and $\epsilon_{r}$ is the radiative efficiency of the accretion flow onto the BH $\epsilon_{r}=0.1$. The final BH accretion rate is:

[TABLE]

and we further prevent the BH to accrete more than 25% of the gas content of the cell in one timestep for numerical stability.

2.3.3 Dynamics

As discussed e.g. in Pfister et al. (2017), dynamical friction plays a large role in setting the detailed dynamics of a BH moving in a galaxy. We however do not have the resolution in our simulation to resolve properly this phenomenon. We include a drag force to model the unresolved dynamical friction exerted by the gas lagging behind the BH as introduced in Dubois et al. (2013). This will mostly play a role at high redshift, when the galaxies are gas rich and clumpy, and will help to maintain the BH in dense gas clumps without resorting to regularly re-positioning the BH. The frictional force is proportional to $F_{\rm DF}=\alpha f_{\rm gas}4\pi\rho(G{M_{\bullet}}/\bar{c_{s}}^{2})$, where $\alpha$ is an artificial boost, with $\alpha=(\rho/\rho_{\rm DF,th})^{2}$ if $\rho>\rho_{\rm th}$ and 1 otherwise, and $f_{\rm gas}$ is a fudge factor varying between 0 and 2 which depends on the BH Mach number $\mathcal{M}_{\bullet}=\bar{v}_{\rm rel}/\bar{c}_{s}$ (Ostriker, 1999; Chapon et al., 2013). In this work, we take $\rho_{\rm DF,th}=50\,\mbox{cm}^{-3}$. Additionally, we include a second contribution to the dynamical friction, this time caused by the collisionless particles (stars and DM). The implementation, based on the analysis of Chandrasekhar (1943); Binney & Tremaine (1987) will be described in Pfister et al. (in prep), and is somewhat similar to that of Tremmel et al. (2015).

2.3.4 AGN feedback

AGN feedback resulting from the accretion onto a BH is modelled using the dual mode implementation of Dubois et al. (2012). At low Eddington ratio $\lambda_{\rm Edd}=\dot{M}_{\rm BHL}/\dot{M}_{\rm Edd}<0.01$, the AGN is in “radio mode”, and in “quasar mode” when $\lambda_{\rm Edd}\geq 0.01$. For the quasar mode, we release thermal energy in a sphere of radius $\Delta x$ centred on the BH at a rate $\dot{E}_{\rm AGN}=\epsilon_{f}\epsilon_{r}\dot{M}_{\bullet}c^{2}$, where $\epsilon_{f}=0.15$ is the fraction of the radiated energy that is transferred to the gas, and $\dot{{M_{\bullet}}}$ is the BH accretion rate. For the radio mode, we deposit momentum as a bipolar outflow with a jet velocity $u_{\rm J}=10^{4}\ \mbox{km\,s}^{-1}$, modelled as a cylinder of radius $\Delta x$ and height $2\Delta x$ weighted by a kernel function

[TABLE]

with $r_{\rm cyl}$ the cylindrical radius. Mass is removed from the central cell and deposited in the cells enclosed by the jet at a rate $\dot{M}_{\rm J}$ with

[TABLE]

where $\Psi$ is the integral of $\psi$ over the cylinder and $\eta_{\rm J}=100$ is the mass-loading factor of the jet accounting for the interaction between the jet and the ISM at unresolved scales. Gas is pushed away from the BH on each side of the cylinder with the norm of the momentum $q(r_{\rm cyl})$ given by

[TABLE]

and we inject the corresponding kinetic energy into the gas:

[TABLE]

The jet is aligned with the total angular momentum of the accreted gas, and the radio mode efficiency is assumed to be $\epsilon_{f}=1$. The feedback efficiencies in both the radio and quasar modes are empirically determined in order to reproduce the BH-to-bulge mass relations at $z=0$ (Dubois et al., 2012).

2.3.5 AGN radiation

Finally, we have developed a radiative feedback implementation for the AGN when it is in quasar mode, similar to the work of Bieri et al. (2017) and to the Appendix C of Trebitsch et al. (2018). We release radiation at each fine timestep, and the amount of radiation released in each frequency bin is given by the luminosity of the quasar in each band. Contrary to Bieri et al. (2017) and Trebitsch et al. (2018), we do not use the averaged Sazonov et al. (2004) spectrum for the BH but instead use a piecewise power-law corresponding to a Shakura & Sunyaev (1973) thin disc extended by a power-law at high energy. In each band, the AGN luminosity is therefore a function of ${M_{\bullet}}$ and $\dot{{M_{\bullet}}}$ (formally, there is also a spin dependence, but we do not follow the spin of the BH in this work). We further assume that $f_{\rm IR}=30\%$ of the bolometric luminosity of the disc is absorbed by dust and re-emitted as IR radiation, that we do not model here, and we keep the radiative efficiency of the accretion flow to $\epsilon_{r}=0.1$.

We ensured that the adopted AGN SED yielded an average spectrum similar to the AGN SED used in Volonteri et al. (2017). However, in the remainder of this work, we only use this SED to estimate the ionizing luminosity of the AGN.

2.4 Gas cooling and heating

Our simulation includes non-equilibrium cooling for hydrogen and helium by tracking the abundances of H, H+, He, He+, He*++*, complemented by metal cooling modelled using the standard rates tabulated in Ramses and computed with Cloudy333http://www.nublado.org/ (last described in Ferland et al., 2017) above $10^{4}$ K. Below $10^{4}$ K, we account for energy losses via metal line cooling following Rosen & Bregman (1995), where the effect of the metallicity is included by scaling linearly the metal cooling enhancement. We do not take into account the impact of the local ionizing flux on metal cooling, but instead assume photo-ionization equilibrium with a redshift dependent Haardt & Madau (1996) UV background for the metals. This UV background is not used for the hydrogen and helium non equilibrium photo-chemistry, which is done using the local photon field transported self-consistently by the RT solver.

3 Results

3.1 Global properties

We run AdaptaHOP on the star particles to identify galaxies in the simulation. The main galaxy is identified at the last timestep ($z\sim 5.7$, illustrated on Fig. 1), and we track back its progenitors using TreeMaker (Tweed et al., 2009). The central BH at that time is then matched to the its host galaxy at each timestep. We show in Fig. 2 the evolution of both the host galaxy (red solid line) and the BH (black dashed line), with the upper panel displaying the evolution of the mass of the galaxy and the BH, and the lower panel showing the star formation rate (SFR) and the BH accretion rate (BHAR). The vertical grey line indicates the last major merger for the host galaxy.

When the galaxy is settled and well identified ($z\lesssim 9$), the SFR slowly increases from $\dot{M}_{\star}\simeq 10\,{\rm M}_{\odot}\mbox{yr}^{-1}$ to $\dot{M}_{\star}\simeq 80\,{\rm M}_{\odot}\mbox{yr}^{-1}$ for a stellar mass growing from $10^{9}{\rm M}_{\odot}$ to almost $M_{\star}\simeq 2.5\times 10^{10}\,{\rm M}_{\odot}$. The specific SFR oscillates around $1-5\,\mbox{Gyr}^{-1}$, and strongly peaks at the time of the major merger ($z\sim 6.3$) as the SFR reaches $\dot{M}_{\star}\gtrsim 200\,{\rm M}_{\odot}\mbox{yr}^{-1}$. The dip in the BH accretion rate at $t\simeq 620-680\,\mbox{Myr}$ is due to the BH being ejected from the central region of the galaxy due to a minor merger. This also explains the stagnation of the BH mass around ${M_{\bullet}}\simeq 1.5\times 10^{5}\,{\rm M}_{\odot}$ during that time. When the BH falls back to the galaxy, it merges with another BH of similar mass to reach ${M_{\bullet}}\simeq 3\times 10^{5}\,{\rm M}_{\odot}$ and then steadily grows until the end of the simulation, where it reaches ${M_{\bullet}}\simeq 1.4\times 10^{7}\,{\rm M}_{\odot}$.

3.2 Self-regulation of the local column density

Let us now focus on the properties of the gas surrounding the BH in the central region of the galaxy in order to examine the accretion-feedback cycle at play when the BH is growing.

In the simulation, accretion is estimated through the Bondi formalism, therefore we can rewrite Eq. 2 as a relation between accretion rate, or Eddington ratio, and column density:

[TABLE]

where ${N_{\rm H}}$ is the local gas column density around the BH measured in the accretion kernel of radius $R_{\rm acc}=4\Delta x$. From this, we can write the dimensionless Eddington ratio ${f_{\rm Edd}}=\dot{{M_{\bullet}}}/\dot{M}_{\rm Edd}$ as

[TABLE]

where $\epsilon_{r,0.1}$ is the radiative efficiency normalized to 0.1, ${M_{\bullet}}_{,6}$ is the BH mass in units of $10^{6}\,{\rm M}_{\odot}$, $\langle N_{\rm H,22}\rangle$ is the average gas column density in units of $10^{22}\,\mbox{cm}^{-2}$ within the accretion radius, $R_{\rm acc,40}$ is the radius within which the accretion rate is estimated, in units of $40~{}{\rm pc}$, and $c_{\rm s,10}$ and $v_{\rm rel,10}$ are the sound speed and relative velocity between gas and BH in units of $10{\,\rm km\,s^{-1}}$. Note that in the simulation $f_{\rm Edd}$ is capped at unity.

From Eq. 11, it appears that a BH needs to be surrounded by dense gas (i.e., by a large average column density) in its vicinity to be able to accrete at high rate. If this were the only physics at play, one should expect a linear correlation between $f_{\rm Edd}$ and $\langle N_{\rm H}\rangle$ at fixed BH mass, temperature and dynamical state. Since temperature and density are not independent, with dense gas being colder than rarefied gas, and relative velocity between the gas and the BH are also a function of the dynamical state, the relation we measure in the simulation is close to $N_{\rm H,22}\propto{f_{\rm Edd}}^{0.50}$. In this analysis, we exclude BHs at ${f_{\rm Edd}}=1$ to avoid being biased by the cap. Including the BHs at ${f_{\rm Edd}}=1$ only marginally affects the trend, and yields $N_{\rm H,22}\propto{f_{\rm Edd}}^{0.56}$. In any case, a higher $f_{\rm Edd}$ corresponds to a higher $\langle{N_{\rm H}}\rangle$ and vice-versa.

We now move to the effect of feedback, through which an accreting BH injects energy in its surroundings. In our simulation, three types of feedback are included: mechanical energy injection at low $f_{\rm Edd}$ , thermal energy injection at high $f_{\rm Edd}$ , and radiation from the AGN at high $f_{\rm Edd}$ . Since we do not include Compton heating from X-rays in the simulation, we do not expect the radiation to play a huge role at the scales we resolve here. Indeed, the photoionization heating from the AGN will at most increase the local gas temperature to a few $10^{4}$ K, and would maintain the local sound speed above 10 km/s. In the simulation, the effective sound-speed $v_{\rm Bondi}=\sqrt{\bar{c}_{s}^{2}+\bar{v}_{\rm rel}^{2}}$ is most of the time above 10 km/s. We therefore limit ourselves to the case of direct energy injection, where we can estimate the relation with gas column density by taking the ratio of the gas binding energy to the injected energy:

[TABLE]

with $M_{\rm gas}=4\pi m_{p}\langle{N_{\rm H}}\rangle R_{\rm acc}^{2}/3$. Here $\Delta t\sim 0.1$ Myr is the coarse timestep in the simulation (hence $\Delta t_{0.1}$ is normalized to 0.1 Myr) and $\epsilon_{f,0.15}$ is the feedback efficiency normalized by 15%, the efficiency of the thermal energy injection. In this derivation, we have defined $M_{\rm tot}$ as the total enclosed mass within $R_{\rm acc}$, but in practice, we do not resolve the influence radius of the BH, and at the resolution of our simulation the gravitational potential is dominated by the stellar component in the vicinity of the BH. This quantity can be seen as a stability criterion for the central gas concentration: when $E_{b}/E_{\rm inj}<1$, the AGN is powerful enough to evacuate the gas from the central region and stopping the accretion onto the BH. At fixed Eddington ratio, Eq. 12 shows that the higher the column density the AGN has to lift, the less effective the feedback will be at unbinding the central gas concentration. We should also keep in mind that dense gas can cool rapidly and thus some of the injected energy is ineffective in unbinding the gas. Recall that for radiative cooling, dominant on the scales and densities of relevance here, $t_{\rm cool}\propto\rho^{-1}\propto{N_{\rm H}}^{-1}$. We do not include the correction here, but the injected energy should be rescaled by the fraction of gas that can cool in one timestep. We also note that Eq. 12 is formally only valid at ${f_{\rm Edd}}>0.01$, because the jet feedback model at lower accretion rate is, by definition, not isotropic.

In Fig. 3, we show the two-dimensional distribution444Note that in this figure as well as all the following (unless mentioned explicitly), we stack all the snapshots together. This effectively smooths any trend that would depend on redshift, but these trend would not be well captured by a single zoom simulation. of the column density $N_{\rm H}$ against the Eddington ratio $f_{\rm Edd}$ in a sphere of 40 pc around the BH, colour-coded by the ratio of the binding energy of the gas to the injected energy from AGN feedback (assuming that the feedback mode is always thermal for simplicity; the grey shaded region marks where this assumption formally breaks down). The yellow line shows the $N_{\rm H,22}\propto{f_{\rm Edd}}^{0.496}\sim{f_{\rm Edd}}^{0.5}$ fit discussed previously. In the red (blue) area, the feedback from the AGN is (not) strong enough to remove the gas from the centre. The general behaviour is perhaps unsurprising: at fixed column density, the BH is unable to unbind the gas when $f_{\rm Edd}$ is low, and conversely at fixed $f_{\rm Edd}$ , the higher the column density and the easier it is for the gas to resist the feedback. It is interesting to note that at fixed $f_{\rm Edd}$ , the average column density (indicated by the yellow line) is just slightly higher than that for which $E_{b}=E_{\rm inj}$ (the white region). This can be understood in terms of AGN feedback self-regulating the BH growth: if the BH is able to remove all of the gas that it feeds on ($E_{b}<E_{\rm inj}$), the accretion rate will drop, and as a result the injected energy from AGN feedback will drop as well. From Eq. 11, we see that AGN feedback reduces ${f_{\rm Edd}}$ faster than $N_{\rm H}$. Indeed, the AGN-driven wind will not only remove some gas, it will also heat whatever remains, therefore reducing further ${f_{\rm Edd}}$ For example, after a feedback event, if $N_{\rm H}$ decreases by a factor 10, $f_{\rm Edd}$ will decrease by more than a factor 10, because the sound speed and relative velocity will increase as well. If the BH falls in a state where $E_{b}>E_{\rm inj}$, then the surrounding gas will be able to cool and feed the BH again. This will be followed by a new episode of feedback. As a result of this feeding-feedback cycle, $E_{b}/E_{\rm inj}$ will hover around unity most of the time, naturally leading to an average behaviour of $N_{\rm H,22}\propto{f_{\rm Edd}}^{0.5}$. During the last phase of the simulation ($t\gtrsim 800$ Myr), when the BH self-regulate its growth, we have a series of episodes where $E_{b}\lesssim E_{\rm inj}$ followed by short episodes where the BH stops accreting, resulting in $E_{\rm inj}\rightarrow 0$ and therefore $E_{b}/E_{\rm inj}\rightarrow\infty$.

This overall paints a picture in which AGN feedback is effectively regulating the gas content in the vicinity of the central BH: when there is too much gas in the centre, the BH will accrete efficiently and release enough energy through feedback to reduce the density and slow down the BH growth. This translates into a modulation of the column density around the BH, that we discuss in the following.

3.3 Obscuration and local column density

We will now discuss the spatial distribution of obscuring material around the BH. To this end, we measure the column density of neutral gas around the BH integrating the local density along rays starting from the BH position and with a length 40 pc. We use a total of 192 rays, and the direction of each ray is drawn from a HealPix555https://healpix.sourceforge.io/ (Górski et al., 2005) tessellation of the sphere.

We show in Fig. 4 the distribution of the fraction of sightlines from the BH, above or below given column densities ($f_{>24}$ is the fraction above $10^{24}\mbox{cm}^{-2}$, $f_{>22}$ is the fraction above $10^{22}\mbox{cm}^{-2}$, $f_{<20}$ is the fraction below $10^{20}\,\mbox{cm}^{-2}$) when the central BH is accreting efficiently (in blue) or not (in red). At higher accretion rates, the BH is more often surrounded by a large fraction of dense gas. This is especially striking in the middle panel of Fig. 4: at high accretion rate, the BH is almost always surrounded by a column density larger than ${N_{\rm H}}\geq 10^{22}\,\mbox{cm}^{-2}$, and this is almost never the case for BH accreting inefficiently. This is in contrast with Ricci et al. (2017b), who analyse the Swift/BAT sample of local AGN and find an anti-correlation between column density and Eddington ratio. Ricci et al. (2017b) interpret their results as obscuration being caused by gas within a few tens of parsecs from the BH, with radiative feedback from dusty gas in this inner region responsible for clearing our gas around the most accretion BHs (see also Lawrence, 1991). They suggest that BHs will loop through a series of high and low accretion rates events, where the feedback from the AGN would clear out the gas obscuring the BH (reducing $N_{\rm H}$ ) while it is still accreting efficiently from a gas reservoir, therefore leading to a large fraction of unobscured AGN at high Eddington ratio. Once the gas reservoir is emptied, $f_{\rm Edd}$ decreases and the AGN will become obscured again as the gas falls back to the BH vicinity. While the feedback process they describe happens on scales smaller than our resolution, we confirm, numerically and analytically, that a BH moves to higher $N_{\rm H}$ as its Eddington ratio increases and that feedback modulates the column density of gas around the BH. In our simulation, however, accretion on the BH is determined instantaneously from the local gas density, without an additional accretion reservoir, therefore we cannot directly assess the time over which a BH can continue accreting after feedback has cleared its surroundings of gas.

Because our simulation includes the galactic and circum-galactic context, we can reproduce a similar analysis for the obscuring material not only in the inner region of the galaxy, but all the way to the circum-galactic medium. For this, we cast again 192 rays from the central BH of the galaxy at each snapshot of the simulation, and integrate the neutral column density over 10 kpc. Fig. 5 compares the column density measured in inner 40 pc around the BH (${N_{\rm H}}_{,\rm loc}$) and in a 10 kpc sphere (${N_{\rm H}}_{,\rm gal}$). The dotted lines marks the ${N_{\rm H}}_{,\rm loc}={N_{\rm H}}_{,\rm gal}$ relation. We can see that for a large part, the galaxy-scale column density (and so the obscuring material) is completely dominated by the small-scale $N_{\rm H}$ , especially in the ${N_{\rm H}}>10^{23}\,\mbox{cm}^{-2}$ regime. However, at lower central $N_{\rm H}$ (and notably below $10^{20}\,\mbox{cm}^{-2}$), there is a significant fraction of the sightlines where the obscuration comes from the gas beyond the centre of the galaxy. More quantitatively, around 14% of all the sightlines have $\xi=\frac{{N_{\rm H}}_{,\rm gal}}{{N_{\rm H}}_{,\rm loc}}>2$, and more than 1% have ${N_{\rm H}}_{,\rm gal}>8\,{N_{\rm H}}_{,\rm loc}$. If we further select only the sightlines with ${N_{\rm H}}_{,\rm loc}$ above (below) $10^{22}\,\mbox{cm}^{-2}$, the fraction of sightlines where $\xi>2$ becomes $\sim 8\%$ ($\sim 40\%$). Outside of the central region, the gas distribution is far from homogeneous (e.g. when the galactic disc forms, as show in Fig. 1). This explains the scatter in ${N_{\rm H}}_{,\rm loc}$ at fixed ${N_{\rm H}}_{,\rm gal}$: an AGN for which we would measure a (galaxy-scale integrated) $N_{\rm H}$ of order $10^{22}\,\mbox{cm}^{-2}$ could either correspond to a BH embedded in a dense enough clump (and would therefore accrete efficiently from a high ${N_{\rm H}}_{,\rm loc}$ reservoir) or to a BH that has cleared its local environment through AGN feedback (low ${N_{\rm H}}_{,\rm loc}$). In the first case, we would expect that the BH is surrounded by dense gas in all directions, while in the second case, the observed $N_{\rm H}$ should depend on the line of sight. We quantify this by defining the porosity of the medium around the BH, defined by the fraction of sightlines with ${N_{\rm H}}<{N_{\rm H}}_{,\rm avg}/5$, where ${N_{\rm H}}_{,\rm avg}$ is the average column density in the 10 kpc sphere around the BH. We illustrate the dependence of the porosity as a function of ${N_{\rm H}}_{,\rm avg}$ in Fig. 6, and the colour-coding indicate the redshift of each snapshot. The three insets correspond to the all-sky projections of the column density around the BH at three specific timesteps of the simulation, with $N_{\rm H}$ increasing from red to blue. The projection is oriented so that the angular momentum of the gas accreted onto the BH is aligned vertically. The figure present three distinct regions: at very high ${N_{\rm H}}_{,\rm avg}$, the porosity is almost always low, meaning that the BH is embedded in a fairly homogeneous clump, as we expected. The low end of the ${N_{\rm H}}_{,\rm avg}$ distribution corresponds to $z\geq 8$, when the galactic disc is not settled at all, and this even includes the event during which the BH spends some time out of the galaxy. The porosity is mostly high, meaning that the BH will be surrounded by a turbulent distribution of gas with a mixture of holes and clumps, as suggested by the inset. The third region is defined by the points at relatively high ${N_{\rm H}}_{,\rm avg}$, but also high porosity: this matches what we described earlier, where the AGN has cleared the local obscuring material, but in certain lines of sight, the galactic disc (the yellow horizontal structure in the inset) still contributes to the total $N_{\rm H}$ .

3.4 AGN luminosity

From the discussion in Sect. 3.2, we expect that the accretion rate onto the BH globally increases with the local column density around the BH. This also holds for the AGN bolometric luminosity, defined as ${L_{\rm bol}}=\epsilon_{r}\dot{M}_{\bullet}c^{2}$. Fig. 7 shows the bolometric luminosity in the simulation as a function of the local column density in a sphere of 40 pc (the accretion region). The background colours mark three regimes of obscuration: blue for an unobscured AGN (${N_{\rm H}}<10^{22}\,\mbox{cm}^{-2}$), yellow for an obscured AGN (${N_{\rm H}}<10^{24}\,\mbox{cm}^{-2}$), and red for Compton-thick AGN (${N_{\rm H}}>10^{24}\,\mbox{cm}^{-2}$). While we do not run any X-ray photon package to estimate directly the X-ray luminosity from the simulation, we use the Hopkins et al. (2007) bolometric correction to convert the bolometric luminosity measured in the simulation to the X-ray luminosity in the $[2-10]$ keV band, through ${L_{\rm bol}}=kL_{\rm X}$, with

[TABLE]

so that $L_{\rm X}={L_{\rm bol}}/k$.

With these definitions, we can quantify the AGN obscured fraction in our simulation, defined here as the fraction of AGN lines of sight with column densities above $10^{23}\,\mbox{cm}^{-2}$, following the work of Vito et al. (2018). Fig. 8 shows the obscured fraction as a function of the intrinsic X-ray luminosity using the column density measured within 10 kpc from the BH. On both panels, the grey shaded region indicates the incompleteness limit from Vito et al. (2018). Overall, we find that the obscured fraction increases with luminosity, in qualitative agreement with Vito et al. (2018), but in contradiction with the study of e.g. Burlon et al. (2011) who find the opposite. Vito et al. (2018) studies a sample of X-ray selected AGN at $z=3-6$ in a parameter space similar to that of the simulation: high-$z$ BHs down to low luminosities. They find, in agreement with our results, a high obscured fraction at high luminosity. Below $10^{43}\,\mbox{erg\,s}^{-1}$ their sample is incomplete, but after correcting for incompleteness, the obscured fraction appears to decrease at low-luminosity, as in our simulation. They note that at low redshift the obscured fraction anti-correlates with luminosity, as found by Ricci et al. (2017b) and Burlon et al. (2011). Vito et al. (2018) note also that the obscured fraction overall increases with redshift, and suggest that the differences with respect to the low-redshift Universe are caused by the larger gas fractions of galaxies at earlier cosmic epochs. Our results support this interpretation. For this simulation, we find that while the galaxy-wide gas contributes somewhat to the total column density of obscuring material (see Fig. 5), the obscuration is often dominated by the nuclear gas. In particular, the results of Fig. 8 are essentially unchanged if we consider only the obscuring material within 40 pc from the BH (shown as dotted lines). This is in part due to the fact that in this simulation, the AGN only reaches luminosities comparable to that of the sample of Vito et al. (2018) when the BH is in a regime where the two column densities (locally and on a galactic scale) are comparable or dominated by the nuclear gas. Recently, Circosta et al. (submitted), reported observations of obscured AGN at $z>2.5$ where the ISM of the galaxy is dense enough to contribute significantly to the obscuration of the central AGN. While our results cannot be extrapolated to higher luminosity, our findings are qualitatively consistent with their study: for most of their objects, they find that the ISM column density is at least comparable to the nuclear column density. We would also expect that the ISM contribution to the obscuration decreases at low redshift, as in the local Universe, only (U)LIRGs reach gas surface density comparable to high-$z$ normal galaxies (e.g. Daddi et al., 2010). This would explain the dichotomy between the results of Vito et al. (2018) and e.g. Burlon et al. (2011): the large gas content of high-$z$ massive galaxies would both significantly contribute to the AGN obscuration and easily replenish the nuclear reservoir once it has been emptied by AGN feedback. This is exactly what our simulation suggest, although at lower AGN luminosity.

We can use a similar approach to convert the AGN bolometric luminosity to a rest-frame UV magnitude at 1450 Å. Namely, we invert the recommended quasar bolometric correction of Runnoe et al. (2012):

[TABLE]

So far, we only considered the intrinsic UV or X-ray luminosity of the AGN: however, we have just show in Figs. 7 and 8 that the BH can be significantly obscured. For the X-rays, we can directly compute the obscuration via $\tau_{X}=\sigma_{\rm T}{N_{\rm H}}$. In the UV, the attenuation mainly comes from dust absorption. We estimate the corresponding wavelength-dependent optical depth, $\tau_{d}(\lambda)$, by integrating the dust column density along the rays used to compute $N_{\rm H}$ (see Sect. 3.3):

[TABLE]

with $\sigma_{d}$ the dust interaction cross section and $n_{d}$ the number density of dust grains. We follow the implementation of the Rascas code (Michel-Dansac et al., in prep) based on the approach of Laursen et al. (2009) to estimate these quantities: $\sigma_{d}$ is defined as the cross section per hydrogen atom, using the fits of Gnedin et al. (2008) for the wavelength dependence based on their Small Magellanic Cloud (SMC) model. The dust number density is then to be understood as a pseudo number density, given by $n_{d}=(n_{{\mathrm{H\,\textsc{i}}}}+f_{\rm ion}n_{{\mathrm{H\,\textsc{ii}}}})Z/Z_{0}$, where $f_{\rm ion}\sim 0.01$ accounts for the presence of dust in ionized gas, and $Z_{0}=0.005$ the mean metallicity of the SMC. We then evaluate the dust optical depth at 1450 Å as $\tau_{1450}=\tau_{d}(1450\,\mbox{\AA})$.

Fig. 9 shows the distribution of $\tau_{X}$ against $\tau_{1450}$, colour-coded by the stellar mass of the galaxy at each timestep. The dotted line indicate $\tau=1$, which corresponds to an attenuation of roughly a factor three ($\exp(-\tau)\sim 0.37$). As the colour goes from yellow to purple, the galaxy grows in mass and the gas around the BH will be more and more metal-rich: this explains why at fixed $\tau_{X}$, the dust optical depth $\tau_{1450}$ increases.

We can then compute an attenuated UV or X-ray luminosity. We illustrate this with Fig. 10, which shows the UV magnitude as a function of X-ray luminosity, colour-coded by the local $N_{\rm H}$ . The red dashed line shows the intrinsic emission (i.e. the bolometric luminosity directly converted in X or UV). The deviations from this line correspond to differential attenuation between the X and the UV due to the presence of dust, and are more significant at higher luminosity. This is due to the fact that on average, the AGN in our simulation is brighter when the galaxy is more massive (and therefore more metal rich), so that the same $N_{\rm H}$ will correspond to a higher dust obscuration. However, even at low $N_{\rm H}$ and low luminosity, the UV luminosity can be significantly reduced.

3.5 AGN duty-cycle

This brings us to distinguish between an “intrinsic duty-cycle”, i.e., the fraction of time the BH is actively growing, and an “observable duty-cycle”, i.e., the fraction of time that a growing BH can be observed as such. If the phases when the BH is most active are characterized by high obscuring column densities (see Figs. 4 and 7), then the measurable growth time will be lower than the intrinsic growth time. The upper panel Fig. 11 shows the fraction of time the BH is active, above a bolometric luminosity or Eddington ratio threshold, and compares the intrinsic duty-cycle to the one that would be inferred from observations. The lower panel shows the corresponding duty-cycle making use of the estimates of $L_{\rm X}$ and $M_{1450}$ discussed in the previous section. Ricci et al. (2017a) show that the optical AGN luminosity function is in good agreement with the X-ray luminosity function of unobscured AGN (${N_{\rm H}}<10^{21}\,\mbox{cm}^{-2}$), therefore we can use ${N_{\rm H}}$ ranges to estimate duty-cycles in the UV/optical and X-rays.

For this particular BH, it would be accreting above ${f_{\rm Edd}}>0.2$ for around 50% of the time, and $\sim 70$% of this time (i.e. $\sim 35\%$ of the total cosmic time) would be characterized by ${N_{\rm H}}\leq 10^{24}\,\mbox{cm}^{-2}$. In other words, as expected, X-ray surveys would miss less than $\sim$ 30% of the growth time. Conversely, fully unobscured (${N_{\rm H}}<10^{20}\,\mbox{cm}^{-2}$) growth above ${f_{\rm Edd}}>0.2$, that we can associate as discussed above to optical visibility, represents less than 2-3% of the total cosmic time. Requiring a growth above ${f_{\rm Edd}}>0.5$ pushes this fraction far below 1% of the cosmic time. This is maybe even more noticeable on the lower panel of Fig. 11: while this specific AGN would be observed with an X-ray luminosity above $3\times 10^{42}\,\rm erg\,s^{-1}$ for approximately a third of the cosmic time, it would be fully unobscured for only 10% of that time (around 3% of the total time). This translates into the UV: the AGN only spends around 2-3% of the time brighter than $M_{1450}\leq-19$, irrespective of the gas column density.

This explains why theoretical models which connect the halo mass function to the optical quasar luminosity function require very low duty-cycles, of order 0.01 (e.g., Haiman & Menou, 2000; Haiman & Hui, 2001; Shankar et al., 2013; Bolgar et al., 2018) while timescale arguments require duty-cycles of order of unity to explain the masses of $z\sim 6$ quasars (e.g., Haiman & Loeb, 2001). Our results are in agreement with, and confirm, their assumptions. Interestingly, the constraints on the quasar lifetimes derived by Chen & Gnedin (2018) are extremely sensitive to the kinematics and distribution of cold gas in quasar hosts.

4 Summary and conclusions

In this work, we have used a cosmological zoom-in simulation of a galaxy at $z\sim 6$, performed with the radiative hydrodynamics code Ramses, to investigate the link between AGN obscuration by gas on the scale of a few tens of parsecs and accretion onto BHs in high-$z$ galaxies.

Our main results are the following:

• •

AGN feedback efficiently regulates the growth of the central BH by controlling the gas reservoir available for accretion (Figs. 2 and 3).

• •

Black holes accreting efficiently tend to be buried under a large column of gas more often than the ones accreting slowly (Fig. 4).

• •

At high redshift, the gas in the galaxy and the circum-galactic medium contributes significantly to the total column density along the line of sight of an AGN (Figs. 5 and 6)

• •

Consistent with the findings of Vito et al. (2018), the obscured fraction (Fig. 8) seems to be very high at high redshift and shows no clear trend with luminosity.

• •

The UV emission from our high-$z$ AGN is more attenuated when the galaxy is more massive / metal rich (Fig. 9) and therefore at higher (observed) X-ray luminosity (Fig. 10).

• •

The fraction of time our BH is active and visible in the UV/optical is very low (of order 2%), while its intrinsic growth duty-cycle is rather of order 50%: the difference is due to the presence of highly obscuring gas in the galaxy (Fig. 11).

Despite the limited observability of an AGN in UV/optical based on its luminosity, the presence of an AGN in a galaxy can be revealed through emission line diagnostics (e.g. Baldwin et al., 1981; Juneau et al., 2011; Feltre et al., 2016), a promising tool as shown by recent simulations (Hirschmann et al., 2017). The combination of deep X-ray surveys, as we confirm here, and emission line diagnostics, especially in the upcoming JWST era, will contribute to derive a more complete census of BH growth in high redshift galaxies.

Acknowledgements

MT and MV acknowledge funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 614199, project ‘BLACK’). This work has made use of the Horizon Cluster hosted by Institut d’Astrophysique de Paris; we thank Stéphane Rouberol for running smoothly this cluster for us. This work was granted access to the HPC resources of CINES under the allocation A0040406955 made by GENCI. This work has made extensive use of the yt666https://yt-project.org/ analysis package (Turk et al., 2011), as well as the Matplotlib (Hunter, 2007), HealPy, Numpy/Scipy (Jones et al., 2001) and IPython (Perez & Granger, 2007) packages.