# Escape of ionizing radiation from high-redshift dwarf galaxies: role of AGN feedback

Maxime Trebitsch, Marta Volonteri, Yohan Dubois, Piero Madau

<figure><img src="image_1.png" /></figure>

## To cite this version:

Maxime Trebitsch, Marta Volonteri, Yohan Dubois, Piero Madau. Escape of ionizing radiation from high-redshift dwarf galaxies: role of AGN feedback. Monthly Notices of the Royal Astronomical Society, 2018, 478 (4), pp.5607-5625. 10.1093/mnras/sty1406. hal-01872280

HAL Id: hal-01872280

https://hal.science/hal-01872280v1

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Escape of ionizing radiation from high-redshift dwarf galaxies: role of AGN feedback

Maxime Trebitsch,<fnref n="1" />* Marta Volonteri,<fnref n="1" /> Yohan Dubois<fnref n="1" /> and Piero Madau<fnref n="1" />,<fnref n="2" />

Accepted 2018 May 23. Received 2018 May 22; in original form 2017 December 8

## ABSTRACT

While low-mass, star-forming galaxies are often considered as the primary driver of reionization, their actual contribution to the cosmic ultraviolet background is still uncertain, mostly because the escape fraction of ionizing photons is only poorly constrained. Theoretical studies have shown that efficient supernova feedback is a necessary condition to create paths through which ionizing radiation can escape into the intergalactic medium. We investigate the possibility that accreting supermassive black holes in early dwarf galaxies may provide additional feedback and enhance the leakage of ionizing radiation. We use a series of high-resolution cosmological radiation hydrodynamics simulations where we isolate the different sources of feedback. We find that supernova feedback prevents the growth of the black hole, thus quenching its associated feedback. Even in cases where the black hole can grow, the structure of the interstellar medium is strongly dominated by supernova feedback. We conclude that, in the dwarf galaxy regime, supermassive black holes do not appear to play a significant role in enhancing the escape fraction and in contributing to the early ultraviolet background.

**Key words:** radiative transfer – galaxies: formation – galaxies: high-redshift – quasars: supermassive black holes – dark ages, reionization, first stars.

## 1 INTRODUCTION

Cosmic reionization is the epoch during which the universe transitioned from the so-called Dark Ages, where most of the hydrogen in the intergalactic medium (IGM) was neutral, to the present epoch, where the IGM is fully ionized. It started as the first luminous objects in the universe formed (around $z \sim 15$), and observations of high-redshift quasars suggest that it proceeded until $z \lesssim 6$ (e.g. Becker et al. 2001; Fan et al. 2001) with some large-scale variations at the very end (Becker et al. 2015). The latest results from the Planck mission suggest that the reionization mid-point (when half of the volume of the universe is ionized) occurs around $z \sim 8.5$, only 600 Myr after the Big Bang (Planck Collaboration XLVI 2016b). While the study of the imprints on the IGM of the early stages of the Epoch of Reionization (EoR) is one of the major endeavours of upcoming 21-cm experiments such as HERA and SKA, there is a strong need to understand the nature of the sources that reionized the universe, especially as the James Webb Space Telescope (JWST) will be able to observe and characterize galaxies deep into the reionization era.

Two kinds of sources have predominantly been invoked as major contributors to the ionizing budget of the EoR since the first evidence

for an ionized IGM was found by Gunn & Peterson (1965): hot, massive stars in star-forming galaxies, and active galactic nuclei (AGNs) powered by supermassive black holes (SMBHs), with the exact contribution of each still highly uncertain. Even though they are extremely bright, the standard picture emerging from high-redshift observations is that quasars are too rare and cannot dominate the ionizing ultraviolet (UV) background during the EoR (e.g. Willott et al. 2010; Fontanot, Cristiani & Vanzella 2012; Grissom, Ballantyne & Wise 2014). In this paradigm, the main driver of reionization must be stellar light from star-forming galaxies (see e.g. Robertson et al. 2010; Haardt & Madau 2012; Kuhlen & Faucher-Giguère 2012). One of the main features of this family of models is that they all require a large number of faint galaxies, with a UV luminosity function (LF) that must be extrapolated down to magnitudes $M_{1500} \lesssim -13$ in order to match the observed integrated Thomson optical depth derived from CMB anisotropies (e.g. Robertson et al. 2013). Such extrapolation has been validated by extremely deep surveys that have shown no sign of a turnover in the steep faint end slope of the UV LF (Atek et al. 2015; Ishigaki et al. 2018; Livermore, Finkelstein & Lotz 2017), consistent with recent simulations of the EoR (e.g. O’Shea et al. 2015; Gnedin 2016; Ocvirk et al. 2016). This galaxy-dominated picture has recently been challenged by recent measurements of the quasar UV LF by Giallongo et al. (2015), who find many more faint AGN at high redshift than previously thought (but see Weigel et al. 2015; Vito et al. 2016; Parsa, Dunlop

<footnote n="1">¹Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="2">²Department of Astronomy and Astrophysics, University of California, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote>* E-mail: maxime.trebitsch@iap.fr</footnote>
McLure 2018, for a different point of view). These recent obser-
vations have revitalized some models where AGNs are contributing
a fraction (Volonteri & Gnedin 2009) or most of the ionizing pho-
tons (Madau et al. 2004; Madau & Haardt 2015). Other models and
seminumerical works appear to exclude an AGN-dominated sce-
nario (Mao & Kim 2016; Hassan et al. 2018; Onoue et al. 2017;
Qin et al. 2017; Ricci et al. 2017; Yoshiura et al. 2017).

Besides the intrinsic luminosity or abundance of either kind of
sources, the escape fraction *f*<sub>esc</sub> of ionizing photons (defined as the
fraction of photons that manage to escape the host halo and ionize
the IGM) is a key parameter in understanding reionization. Despite
more than 20 yr of effort (Leitherer et al. 1995; Deharveng et al.
2001), direct observations of the Lyman continuum (LyC) below
912 Å are still very difficult, and there are only a few galaxies for
which *f*<sub>esc</sub> has been directly measured (see e.g. Bergvall et al. 2013
for a review). At high redshift, the increasing opacity of the IGM
makes this even harder, and foreground contamination becomes a
major issue (Siana et al. 2015). Moreover, observations yield a large
scatter on the value of *f*<sub>esc</sub> between galaxies of different types and at
different redshift, with no clear pattern: Izotov et al. (2016) recently
reported LyC leakage in all of the four galaxies they observed at z ~
0.3, with *f*<sub>esc</sub> ≃ 8–13 per cent, while at intermediate redshift, Bridge
et al. (2010) present stringent upper limit of 1 per cent on a stacked
sample. Higher redshift observations of Vanzella et al. (2016) report
a very high *f*<sub>esc</sub> ~ 50 per cent for a compact star-forming galaxy at
z ~ 3, while Grazian et al. (2017) find only *f*<sub>esc</sub> < 1.7 per cent for
a stacked sample of bright galaxies at similar redshift, and Vasei
et al. (2016) find a much lower value of *f*<sub>esc</sub> < 8 per cent at z ~ 2.38.
This large variation from galaxy to galaxy might be explained by
variations in the luminosity (Japelj et al. 2017) or in the physical
properties of the galaxy (Izotov et al. 2016). Numerical simula-
tions suggest a similar large scatter in the value of *f*<sub>esc</sub> at fixed halo
mass (Kimm & Cen 2014; Paardekooper, Khochfar & Dalla Vec-
chia 2015), and this large variability is correlated with the intrinsic
burstiness of star formation in low-mass galaxies and the result-
ing supernova (SN) feedback (Wise et al. 2014; Ma et al. 2015;
Trebitsch et al. 2017). AGN-dominated reionization models usu-
ally assume that *f*<sub>esc</sub><sup>AGN</sup> ~ 100 per cent, meaning that all the emitted
ionizing photons contribute to the ionizing background (e.g. Gial-
longo et al. 2015), consistent with observations of bright quasars
(Cristiani et al. 2016). However, for fainter AGN, the observations
of Micheva, Iwata & Inoue (2017) suggest *f*<sub>esc</sub><sup>AGN</sup> well below unity.

At faint luminosities, it becomes increasingly difficult to separate
the AGN contribution from the starlight of its host, and conversely,
ionizing radiation coming from a source identified as an AGN might
actually be partially of stellar origin (see e.g. Volonteri et al. 2017).
For instance, the LyC emission coming from two AGNs out of the
four LyC candidates targeted by Micheva et al. (2017) could be of
stellar origin. Furthermore, even when they are undetected, AGNs
can still have an impact on the interstellar medium (ISM) of high-
redshift galaxies, providing an additional source of feedback in ad-
dition to SNs. This AGN feedback could in turn have a strong effect
on the escape of ionizing radiation from galaxies. Using Chandra
data, Kaaret et al. (2017) find some evidence for the presence of a
compact X-ray source in the LyC leaking galaxy Tol 1247-232 that
might be powered by a low-luminosity AGN. They further suggest
that the AGN could be responsible for outflows helping the escape
of ionizing radiation.

While signs of nuclear activity have been detected for three
decades in nearby dwarf galaxies (Kunth, Sargent & Bothun 1987;
Filippenko & Sargent 1989), systematic searches for black holes
(BHs) in the nuclei of dwarf galaxies are still fairly recent. Reines,

Greene & Geha (2013) (see also Moran et al. 2014) present a large spectroscopic sample of 136 dwarf galaxies 10$^{8.5}$ M$_{\odot} \lesssim M_{\star} \lesssim 10^{9.5}$ M$_{\odot}$ with BH signature, among which 35 have their spectra dominated by the AGN. This suggests that BHs, with masses sometimes as low as $5 \times 10^4$M$_{\odot}$ (Baldassare et al. 2015), can be present even in the lowest mass galaxies, such as NGC 4395, a galaxy with stellar mass below $10^8$M$_{\odot}$ (Filippenko & Sargent 1989). Searches for BHs in low-mass galaxies are, however, extremely difficult due to the contamination by stellar light at all wavelengths (see Reines & Comastri 2016, for a review of AGN searches in low-mass galaxies). At $z \lesssim 0.3$, Chen et al. (2017) find a sample of 10 more low-mass galaxies harbouring AGN activity. Mezcua et al. (2016) recently uncovered a population of BHs in dwarfs at $z \lesssim 1.5$ in the COSMOS survey. As dwarf galaxies have a simpler evolution than massive galaxies, their massive BH may be relics of seeds formed at very high redshift (e.g. Volonteri, Lodato & Natarajan 2008; Volonteri 2010), and the fraction of dwarf galaxies hosting massive BH could provide valuable clues to help understanding the formation of the progenitors of massive BH.

In this paper, we use high-resolution radiation hydrodynamics
simulations of a galaxy in a cosmological context to investigate
the growth of BHs in dwarf galaxies during the EoR. We take an
agnostic view with respect to the seeding mechanism, and we focus
only on the coevolution of high-redshift dwarfs and their putative
central BH. Our main goal is to understand the interplay between the
BH growth and the various feedback processes involving massive
stars and AGN. We then discuss the importance of the different
sources of feedback on the escape of ionizing radiation from the
galaxy.

The paper is structured as follows. We present in Section 2 the
numerical methods employed in this work, and we then discuss
the connection between galaxy and BH growth in Section 3. We
investigate how this affects the production and escape of ionizing
radiation in Section 4, and we discuss our results in Section 5 before
concluding in Section 6.

## 2 DESCRIPTION OF THE SIMULATIONS

The simulations presented here are run with *RAMSES-RT* (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), a public multigroup radiative transfer (RT) extension of the adaptive mesh refinement (AMR) code *RAMSES<fnref n="1" />* (Teyssier 2002). The collisionless dark matter (DM) and stellar particles are evolved with a particle-mesh solver with cloud-in-cell interpolation. We follow the gas evolution by solving the Euler equations with a second-order Godunov scheme using the HLLC Riemann solver (Toro, Spruce & Speares 1994), with a MinMod total variation diminishing scheme to reconstruct the intercell fluxes. We impose a Courant factor of 0.8 for all simulations. We use a standard quasi-Lagrangian refinement strategy, in which a cell is refined if $\rho_{\text{DM}}\Delta x^3 + \frac{\Omega_{\text{DM}}}{\Omega_b} \rho_{\text{gas}}\Delta x^3 + \frac{\Omega_{\text{DM}}}{\Omega_b} \rho_* \Delta x^3 > 8 m_{\text{DM}}^{\text{HR}}$, where $\rho_{\text{DM}}$, $\rho_{\text{gas}}$, and $\rho_*$ are, respectively, the DM, gas, and stellar densities in the cell, $\Delta x$ is the cell size, and $m_{\text{DM}}^{\text{HR}}$ is the mass of the highest resolution DM particle; this would allow refinement as soon as there are at least eight high-resolution DM particles in a cell in a DM-only run.

The RT module propagates radiation emitted by stars in three
energy intervals between 13.6–24.59 eV, 24.59–54.4 eV, and over
54.4 eV, respectively, describing the H, He, and HeII photon field.
The equation of RT is solved on the AMR grid using a first-order

<footnote n="1">¹https://bitbucket.org/rteyssie/ramses/</footnote>
Godunov method with the M1 closure (Levermore 1984; Dubroca & Feugeas 1999) for the Eddington tensor. The radiation is coupled to the hydrodynamical evolution of the gas through the non-equilibrium thermochemistry for hydrogen and helium. Since we are mainly interested in the absorption of radiation within the galaxy (rather than the expansion of cosmological H II regions in the IGM), we can use the reduced speed of light approximation, $\bar{c} = 0.01c$. Because RAMSES-RT can handle an arbitrarily large number of sources, we emit radiation from all star particles at each timestep using the GALAXEV model of Bruzual & Charlot (2003), assuming a Chabrier (2003) initial mass function (IMF). Finally, we use the on-the-spot approximation: any ionizing photon emitted by recombination is assumed to be absorbed locally, and we thus ignore direct recombinations to the ground level and the associated emission of ionizing radiation from the gas. We also ignore the radiation pressure from the UV field to the gas, which has been shown to have only a subdominant effect (Rosdahl et al. 2015).

## 2.1 Initial conditions

To select the halo of interest, we first run for 1 Gyr a DM-only simulation of a region of $10 h^{-1}$ comoving Mpc on a side, sampled with 512³ particles ($m_{DM} \simeq 6.4 \times 10^5 h^{-1} M_{\odot} \simeq 9.4 \times 10^5 M_{\odot}$) initialized at $z = 150$. We use the MUSIC<fnref n="2" /> code (Hahn & Abel 2011) to generate the initial conditions, and assume 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 XIII 2016a)

The target halo is selected to have a total mass of $5 \times 10^9 M_{\odot}$ in the final output at $z \sim 5.7$ with HALOMAKER (Tweed et al. 2009) using the ADAPTAHOP algorithm (Aubert, Pichon & Colombi 2004). We generate high-resolution initial conditions for that halo by refining the mass distribution such that the DM particle mass is $m_{DM}^{HR} \sim 1.5 \times 10^4 M_{\odot}$ in the zoomed-in region, equivalent to a 2048³ particles simulation. This corresponds to an AMR coarse grid level of $\ell = 11$ in the zoomed-in region, and we allow for 10 more levels of refinement, resulting in a most refined cell size of $\Delta x = 10 h^{-1} Mpc/2^{21} = 7$ pc. The gas is assumed to be initially neutral, and with an initial gas phase metallicity $Z = 10^{-3} Z_{\odot} = 2 \times 10^{-5}$ everywhere.

## 2.2 Galaxy formation physics

The subgrid recipes used in this work to model star formation, the subsequent feedback from massive stars, and the various gas cooling and heating processes are very similar to those of Trebitsch et al. (2017). In this section, we will summarize them briefly and highlight the few differences in their implementation.

### 2.2.1 Star formation

We model star formation using the model presented in Trebitsch et al. (2017) and Kimm et al. (2017) (see also Devriendt et al. in preparation). The approach is similar to that of Rasera & Teyssier (2006): in star-forming cells, $M_{sf} = \epsilon \rho_{gas} \Delta x^3 \Delta t / t_{ff}$ of gas is converted into star particles during one timestep $\Delta t$, where $G$ is the gravitational constant, $\epsilon$ is the local star formation efficiency, and $t_{ff} = \sqrt{3\pi/32G\rho_{gas}}$ is the gas free fall time. The actual number $N$

of star particles formed in one timestep $\Delta t$ is drawn from a Poisson distribution $P(N) = (\lambda^N/N!)\exp(-\lambda)$ of parameter $\lambda = M_{sf}/m_{\star}$, where $m_{\star} \simeq 2000 M_{\odot}$ is the minimum mass of a star particle in our simulations.

This model describes star-forming regions as clouds characterized by their average gas density $\rho_0$, sound speed $c_s$, and local velocity dispersion $\sigma_{gas}$ (which accounts for the turbulence in the cloud). The local star formation efficiency $\epsilon$ is computed following the ‘multi-ff PN’ model of Federrath & Klessen (2012), based on Padoan & Nordlund (2011).

$$ \epsilon \propto e^{\frac{3}{8}\sigma_s^2} \left[ 1 + \operatorname{erf}\left(\frac{\sigma_s^2 - s_{\operatorname{crit}}}{\sqrt{2\sigma_s^2}}\right) \right], \quad (1) $$

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

The primary difference with Trebitsch et al. (2017) is in the way we identify star-forming regions: instead of selecting cells that are Jeans-unstable (even with the additional turbulent support), we use a more relaxed criterion by simply enforcing the gas to be super-sonic. The motivation for this modified approach is twofold. First, the models of Federrath & Klessen (2012) break down when the gas becomes trans- or subsonic, so we require a Mach number $M \ge 2$ in star-forming regions. This is well motivated by observations of molecular clouds, in which the Mach number is usually far greater than unity. Note that we only allow star formation in cells where $\rho > \rho_{th} = 5 \text{ cm}^{-3}$ in order to avoid forming stars outside of the high-resolution region. Second, we drop the turbulent Jeans criterion to select star-forming regions. Indeed, enforcing the ‘turbulent Jeans length’ to be smaller than four cells directly translates into a constraint on the virial parameter of the cloud to be $\alpha_{vir} \lesssim 5-6$, while observations show a large scatter with values sometimes above this threshold (see for instance Kauffmann, Pillai & Goldsmith 2013). In our model, the dependence of the star formation efficiency on the virial parameter is already encoded in equation (1), such that a cloud with high virial parameter can form stars, but at a very low efficiency, consistent with the cloud-scale, high-resolution simulations of e.g. Padoan et al. (2017). We however checked that this change in the star formation model does not affect the global star formation history of the galaxy. Indeed, the fact that $\epsilon$ is non-uniform in the galaxy results by itself in a clumpy star formation process.

### 2.2.2 Feedback from massive stars

We implement two channels of stellar feedback: radiative feedback from stars resulting from photoionization heating and type II SNs explosions. Photoionization heating is the main channel of radiative feedback on galactic scale (Rosdahl et al. 2015), and is naturally included in the simulations as a consequence of the ionizing RT. For SN feedback, we use a mechanical feedback model similar to that first described by Kimm & Cen (2014) and Kimm et al. (2015), and later refined by Rosdahl et al. (2017). We release momentum in the neighbouring cells (and energy in the SN host cell) in a single event 5 Myr after the birth of the star particle. This scheme was developed to capture correctly the momentum transfer at all stages of the Sedov blast wave, and the amount of momentum deposited in the simulation depends on the local properties (density and metallicity) of the gas. We use a modified version of this scheme, implemented by Kimm et al. (2017), which takes into account the argument of Geen et al. (2015) that the photoionization pre-processing of

<footnote n="2">²https://bitbucket.org/ohahn/music/</footnote>
the ISM prior to the SN can augment the final radial momentum from an SN. While this should be taken into account by the RT in our simulation, Trebitsch et al. (2017) argued that a significant fraction of the radiation is actually emitted in regions where the Strömgren radius is not resolved, leading to a likely underestimation of this momentum increase. We thus follow the subgrid model of Kimm et al. (2017), which adds this missing momentum when the Strömgren radius is locally not resolved. We again checked that this change does not affect our results at an appreciable level.

### 2.2.3 Gas cooling and heating

RAMSES/GRAMSES-RT features non-equilibrium primordial cooling tracking the abundances of five species: H, H+, He, He+, He++. We include an extra cooling term for metals, using the standard tabulated rates included in and computed with CLOUDY<fnref n="3" /> (last described in Ferland et al. 2017) above $10^4$ K. We also account for energy losses via metal line cooling below $10^4$ K following the prescription of Rosen & Bregman (1995), and approximate the effect of the metallicity by scaling linearly the metal cooling enhancement. We do not take into account the impact of the local ionizing flux on metal cooling. Instead, metal cooling is computed assuming photoionization equilibrium with a redshift dependent Haardt & Madau (1996) UV background. We stress that we do not use this UV background for the hydrogen and helium non-equilibrium photochemistry, but rather take into account the local photoheating rate self-consistently as we follow the radiation transport on the fly, which would be the physically reasonable assumption for an isolated galaxy before the end of the EoR. We adopt an initial homogeneous metallicity floor of $Z = 10^{-3}Z_\odot$ in the whole box to compensate for the lack of molecular cooling and to allow the gas to cool down below $10^4$ K before the first stars have formed.

## 2.3 BH model and AGN feedback

We use the fiducial implementation of Dubois et al. (2012) for the BH seeding and growth and the associated feedback from the AGN. In the following, we recall the main features of the model.

BHs are represented by sink particles of initial mass $M_{\bullet,0} = 10^4 M_\odot$ that are created in cells that meet the following criteria: both the gas and stellar densities must exceed a density threshold $\rho_{\text{sink}}$, the cell must be Jeans-unstable, and there must be enough gas in the cell to form the sink. In order to avoid the creation of multiple BH in the same galaxies, we also enforce that there should be no other sink particle closer than an exclusion radius $r_{\text{excl}}$. In our simulation, we use a density threshold $\rho_{\text{sink}} = 400 \text{ cm}^{-3}$ (which corresponds to $10^4 M_\odot$ in a cell of 10 pc), and an exclusion radius $r_{\text{excl}} = 50 \text{ kpc}$ (comoving). Once the sink particle has been created, we spawn 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 averaging the properties of the gas around the BH.

BHs are allowed to grow both via gas accretion and via mergers. The accretion rate is computed using the Bondi–Hoyle–Lyttleton prescription (Bondi 1952):

$$ \dot{M}_{\text{BHL}} = 4\pi G^2 M_{\bullet}^2 \frac{\bar{\rho}}{(\bar{c}_s^2 + \bar{u}^2)^{3/2}} \quad (2) $$

where $M_{\bullet}$ is the BH mass, $\bar{\rho}$ and $\bar{c}_s$ are, respectively, the average gas density and sound speed, and $\bar{u}$ 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(-r^2/r_{\text{sink}}^2)$, where $r_{\text{sink}}$ is defined using the Bondi radius $r_{\text{BHL}} = \frac{GM_{\bullet}}{\bar{c}_s^2 + \bar{u}^2}$:

$$ r_{\text{sink}} = \begin{cases} \Delta x/4 & \text{if } r_{\text{BHL}} < \Delta x/4, \\ r_{\text{BHL}} & \text{if } \Delta x/4 \le r_{\text{BHL}} < 2\Delta x, \\ 2\Delta x & \text{if } 2\Delta x \le r_{\text{BHL}}. \end{cases} \quad (3) $$

The accreted gas is removed from the cells surrounding the BH using the same weighting. Note that we do not boost accretion as is usually done in kpc-resolution cosmological simulations (see e.g. Booth & Schaye 2009) because with our set of simulations we are able to resolve dense gas clumps. Even though we do not resolve the clumpy accretion onto the BH on the smallest scales, our ~10 pc resolution allows us to start to account for the multiphase structure of the ISM, and we are close to the regime where the non-boosted Bondi formalism captures correctly the accretion rate (Negri & Volonteri 2017). The BH accretion rate is capped at the Eddington rate, $\dot{M}_{\text{Edd}} = 4\pi GM_{\bullet}m_p/(\epsilon_r\sigma_Tc)$, where $m_p$ is the mass of a proton, $\sigma_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, taken to be the canonical value of $\epsilon_r = 0.1$. For numerical stability, we prevent the accretion to remove more than 25 per cent of the gas content of a given cell. Finally, following galaxy mergers, BH binaries are allowed to coalesce when they are closer than $4\Delta x$ from each other and when their relative velocity is smaller than the escape velocity of the binary.

In order to avoid spurious spatial oscillations of the BH, we include a drag force to model the unresolved dynamical friction from the gas on the BH as introduced in Dubois et al. (2013). The frictional force is proportional to $F_{\text{DF}} = \alpha f_{\text{gas}}4\pi\rho(GM_{\bullet}/\bar{c}_s^2)$, where $\alpha$ is an artificial boost introduced to account for unresolved small scale structures, with $\alpha = (\rho/\rho_{\text{th}})^2$ if $\rho > \rho_{\text{th}}$ and 1 otherwise, and $f_{\text{gas}}$ is a fudge factor varying between 0 and 2 which depends on the BH Mach number $M_{\bullet} = \bar{u}/\bar{c}_s$ (Ostriker 1999; Chapon, Mayer & Teyssier 2013).

We include the same dual mode of AGN feedback as in Dubois et al. (2012), with “radio mode” feedback when the Eddington ratio $\lambda_{\text{Edd}} = \dot{M}_{\text{BHL}}/\dot{M}_{\text{Edd}} < 0.01$ and “quasar mode” when $\lambda_{\text{Edd}} \ge 0.01$. For the Aquasar mode’, we release thermal energy at the coarse timestep in a sphere of radius $\Delta x$ around the BH at a rate $\dot{E}_{\text{AGN}} = \epsilon_f \epsilon_r 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. In the radio mode, the energy is deposited as a bipolar outflow with a jet velocity $u_J = 10^4 \text{ km s}^{-1}$, modelled as a cylinder of radius $\Delta x$ and height $2\Delta x$ weighted by a kernel function

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

with $r_{\text{cyl}}$ the cylindrical radius. We denote as $\Psi$ the integral of $\psi$ over the cylinder. Mass is removed from the central cell and deposited in the cells enclosed by the jet at a rate $\dot{M}_J$ with

$$ \dot{M}_J(r_{\text{cyl}}) = \frac{\psi(r_{\text{cyl}})}{\Psi} \eta_J \dot{M}_{\bullet} \quad (5) $$

where $\eta_J = 100$ is the mass-loading factor of the jet accounting for the interaction between the jet and the ISM at unresolved scales. Momentum is deposited with opposite outflowing directions on each

<footnote n="3">³http://www.nublado.org/</footnote>
<table><caption>**Table 1.** Summary of all the simulations presented in this work and the variations on the physical processes used.</caption><thead><tr><th>Simulation</th><th>SN</th><th>AGN</th><th>Accretion</th></tr></thead><tbody><tr><td>AGN_SN</td><td>√</td><td>√</td><td>Bondi</td></tr><tr><td>SN_only</td><td>√</td><td>–</td><td>–</td></tr><tr><td>AGN_only</td><td>–</td><td>√</td><td>Bondi</td></tr><tr><td>BH_SN</td><td>√</td><td>–</td><td>Bondi</td></tr><tr><td>AGNforced_SN</td><td>√</td><td>√</td><td>forced</td></tr></tbody></table>

side of the cylinder with a norm $q(r_{cyl})$ with

$$q_J(r_{cyl}) = \dot{M}_J(r_{cyl})u_J, \qquad (6)$$

and we inject the corresponding kinetic energy into the gas:

$$\dot{E}_J(r_{cyl}) = \frac{q_J(r_{cyl})^2}{2\dot{M}_J(r_{cyl})} = \frac{\psi(r_{cyl})}{\Psi} \eta_J \dot{E}_{AGN}. \qquad (7)$$

The jet is aligned with the total angular momentum of the accreted gas. As a result, in a turbulent environment such as the ISM of high-z dwarf galaxies, where no clear disc can be seen (see Section 3.1), the outflow is not required to be perpendicular to the angular momentum of the whole galaxy, but only to that of the clump in which the BH is hosted. The radio mode efficiency is assumed to be larger than in the Åquasar mode', with $\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.4 Summary of the simulations

In order to assess the exact role of BH growth and the associated feedback on the early evolution of low-mass galaxies, we rerun the simulation while varying the origin and amount of feedback. We present on Table 1 the list of the simulations used in this work.

For our fiducial model, which we will note AGN_SN, we use all the subgrid models described in the previous sections: mechanical feedback for SNs, dual-mode AGN feedback, and Bondi accretion. For comparison with Trebitsch et al. (2017), our SN_only simulation only includes SN feedback, and no BH physics at all. We then alternatively removed SN feedback (AGN_only) and AGN feedback (BH_SN) to isolate their respective effects. Finally, we perform a simulation where we allow the BH to accrete as much gas as is available (we prevent the BH to remove more than 25 per cent of the gas available) while keeping the Eddington limit as a maximal accretion rate rather than using the Bondi prescription, which we will denote by AGNforced_SN. Since we do not boost the accretion explicitly, this provides an extreme case to assess the robustness of our results, were we to underestimate the real accretion rate of our BH using the standard Bondi approach.

## 3 GALAXY-BLACK HOLE COEVOLUTION

In this section, we discuss how the main galaxy in our simulations assembles its stellar mass, and how this correlates with the growth of the central BH. We then discuss the role of the various feedback processes in controlling the accretion onto the BH.

### 3.1 Stellar mass assembly

We illustrate in Fig. 1 (left-hand column) the gas density around the main galaxy in our AGN_SN simulation projected along the x and

z lines of sight of the simulation, with the virial radius shown as a large circle. The right-side column presents mock ugr images of the galaxy at $z \sim 5.7$, zooming on the inner region within 50 per cent of the virial radius. The galaxy resides in a large-scale filament and is well isolated from its nearest neighbours. The gas distribution is very clumpy and disturbed, and the galaxy has clearly not settled in a fully formed, well-defined disc, but rather shows a disturbed morphology. BHs are shown as the black dots, and we can see that while one of the BHs resides in the centre of the galaxy, the others have not yet fallen at the bottom of the potential well and are still attached to smaller clumps.

In order to compare our simulation to high-redshift observations, we derive the physical properties of the galaxy (stellar mass and SFR) using an observational approach. We compare these estimates with the quantities directly measured in the simulation in Appendix A. We estimate the stellar mass from the UV continuum at 1500 Å based the empirical power-law relation obtained by Song et al. (2016). We assume that each star particle in the simulation has the spectral energy distribution (SED) of a single stellar population (SSP) with the same age, mass, and metallicity as the particle. Using Bruzual & Charlot (2003) models, we compute the UV luminosity at $\lambda = 1500$ Å of each star particle, and measure the UV magnitude of the galaxy by summing over all the luminosities. This is likely to yield an overestimate of what would actually be measured in observations, since we assume no flux losses due to low surface brightness. However, given the very small size of the galaxy ($\lesssim 1$ kpc or $\lesssim 0.2$ arcsec at $z \sim 6$), it would be at most barely resolved even with the high-resolution imaging of Song et al. (2016). We also assume that there is no dust attenuation, which is a fairly reasonable approximation at first order in the high-redshift, low-mass regime. We assume that the stellar mass can be related to the UV magnitude through $\log M_* = a(z) + b \times M_{1500}$, and we linearly interpolate the best-fitting normalization of Song et al. (2016, Table 1) between $z = 5$ and $z = 8$ with fixed slope of $b = -0.5$. We note that Song et al. (2016) assume a Salpeter (1955) IMF, while we use a Chabrier (2003) IMF to estimate the luminosity of our stellar particles. At fixed age, mass and metallicity, an SSP populated with a Salpeter (1955) IMF will contain more low-mass stars, so its UV luminosity will be lower. Conversely, an SSP needs to be more massive to yield the same UV luminosity. We thus reduce the $\log M_* - M_{1500}$ scaling of Song et al. (2016) by 0.2 dex in $M_*$ to take this into account.

We present in Fig. 2 the resulting stellar mass-to-halo mass (SMHM) relation for the AGN_SN simulation. The solid, dashed, and dotted lines respectively indicate 100 per cent, 10 per cent, or 1 per cent of the universal baryon fraction. Using this UV estimate for the stellar mass, it appears that the galaxy converts around 1 per cent of its baryons into stars, with some increase at higher halo mass. While this is above the predictions from abundance matching techniques (e.g. Behroozi, Wechsler & Conroy 2013, shown here as dash-dotted lines), we note that such models are explicitly not designed for the mass regime we are studying in this work (see the discussion in Moster, Naab & White 2018). The large scatter at fixed halo mass mostly comes from the estimate of the stellar mass. Indeed, at low mass, galaxies assemble their mass by burst. Since the UV luminosity is dominated by massive, short-lived stars, a galaxy right after a star-forming episode will appear brighter than a few tens of Myr later. At that stage, even though the real stellar mass does not decrease, our UV-based estimation will yield a lower stellar mass, and should therefore be considered as a lower limit of the true stellar mass. We still prefer to use this quantity rather than the stellar mass measured in the simulation, since it suffers from the same uncertainties as high-redshift observations (e.g. degeneracy
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Views of the main galaxy in our AGN-SN simulation along the *x* and *z* lines of sight. Left: mass-weighted projection of the gas density, with the BHs shown as black dots and the virial radius of the halo as a white circle. Right: Mock ugr image, zoomed on the galaxy at *z* ~ 5.7.</figcaption></figure>

with the star formation history). By comparison, we also show in Fig. 2 the observations of Miller et al. (2014), who present a sample of 41 galaxies at much lower redshift ($0.2 \lesssim z \lesssim 1$) and for which they compute halo masses from dynamical modelling and stellar masses from photometry. They find a SMHM systematically above the models of Behroozi et al. (2013), with roughly 10 per cent of the universal baryon fraction converted into stars, very similar to our results in the mass range where observations and simulations overlap.

Still with the focus of comparing our simulation to observable quantities, we estimate star formation rates (SFR) using the far-UV calibration from Kennicutt & Evans (2012). While they assume a Kroupa & Weidner (2003), models with a Chabrier (2003) IMF yield identical results, so we include no correction for the IMF conversion.

We show in Fig. 3 where our galaxy falls on the star-forming main sequence (MS), with the UV-derived SFR as a function of the UV-derived stellar mass colour-coded with time since the beginning

of the simulation. The various dashed lines indicate observational fits of the MS at *z* ~ 6 from Salmon et al. (2015, in red), Santini et al. (2017, in pink), and Speagle et al. (2014, in grey), and the shaded areas represent the standard deviation on the MS, shown when available. Both for Salmon et al. (2015) and Speagle et al. (2014), the standard deviation is of order 0.2 dex. The grey dotted line shows the MS extrapolated to *z* ~ 8 from Speagle et al. (2014), and is almost indistinguishable from the MS at *z* ~ 6. We overplot in orange an extrapolation of the MS of Behroozi et al. (2013) as fitted by Salmon et al. (2015). We highlight that the apparent tight correlation between the stellar mass and the SFR comes from the fact that both quantities are derived from the UV luminosity, and we are effectively plotting against each other two functions of the same quantity. If we were to use the quantities directly measured in the simulation, we would find a much wider scatter (see Fig. A3). Nevertheless, as stated previously, we prefer to use the UV-derived properties in order to get a fairer comparison to the observations shown on the figure. Overall, the galaxy in our
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Stellar-to-halo mass relation for our simulated galaxy, using the UV estimate for the stellar mass. The solid, dashed, and dotted line corresponds, respectively, to 100 per cent, 10 per cent and 1 per cent of the baryons converted into stars.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Star formation MS based on UV properties for the AGN\_SN simulation, colour-coded by redshift. The lines correspond to SF MSs extrapolated from the models of Behroozi et al. (2013, orange) and from the observations of Salmon et al. (2015, red), Santini et al. (2017, pink), and Speagle et al. (2014, grey), with their standard deviation shown as a shaded area when available.</figcaption></figure>

simulation falls within the various lines, with a somewhat shallower slope than most observations. We note however that the uncertainty on the exact shape of the MS at high redshift exceeds by far the standard deviation of ~0.2 dex. As a result, extrapolations to low masses such as those we present in Fig. 3 are only trends around which some scatter is naturally expected due to intrinsic burstiness of star formation in low-mass galaxies. While we should not consider quantitatively the position and evolution of our galaxy on the MS, the broad consistency between our simulation and observations suggests that our simulated galaxy assembles its mass at a reasonable pace.

The other simulations, including SN feedback (SN_only, BH\_SN and AGN\_forced\_SN), present a very similar stellar mass assembly and star formation history, while the simulation without SN feedback (AGN_only) forms much more stars, as expected.

## 3.2 Black hole growth

We present in Fig. 4 the mass growth of the most massive BH of our simulations as a function of the stellar mass of the host galaxy; the colour coding indicates the redshift. Note that contrary to the previous section, we use here the stellar mass directly measured in the simulation. The various symbols represent various simulations: stars for the AGN\_SN case, circles for AGN_only, crosses for BH\_SN, and triangles for the forced accretion simulation (AGN-forced\_SN). We only show one in every five snapshots to make the figure easier to read, and join the symbols with thin lines to account for the complete growth history. As guide for the eye, we indicate constant BH to stellar mass ratio of $10^{-3}$ ($10^{-4}$) as a dashed (dotted) line, and the empirical relation derived by Reines & Volonteri (2015) for galaxies in the local universe as a dash-dotted line.

It appears very clearly that when SNe and AGN are included, the BH does not grow at all. In 1 Gyr, the most massive BH of the AGN\_SN simulation, which lies at the centre of the galaxy in Fig. 1, grows by less than 1000 $M_{\odot}$. By comparison, when SN feedback is removed, the most massive BH in AGN_only grows much faster and reaches a higher mass, roughly following $M_{\bullet} \simeq 10^{-3} M_{*}$. Interestingly, in the BH\_SN simulation (with only SN feedback), the BH reaches a slightly higher mass than in the fiducial simulation. Indeed, over the course of the whole simulation, the BH grows by ~3500 $M_{\odot}$. For the simulation where the back hole is forced to accrete (AGN\_forced\_SN, diamonds in Fig. 4), it grows more regularly but still much less than in the AGN_only case.

This seems to indicate that the main mechanism regulating the growth of BHs in small galaxies is SN feedback. Indeed, even when the BH is forced to accrete independently of the local gas properties, its growth is strongly limited by the lack of available gas. We develop this argument further in Fig. 5, which presents the local gas density around the main BH in each simulation. The accretion rate is shown as a grey background, with accretion rates below $10^{-2} M_{\text{Edd}}$ in white and a logarithmic scaling up to $\dot{M}_{\text{Edd}}$ in black. The two upper plots (AGN\_SN and BH\_SN, respectively, on panels 5a and 5b) display very similar features, with a low gas density on average, oscillating around $0.1-1 \text{ cm}^{-3}$. The rare peaks of high densities correspond to episodes of stronger accretion, but these only last for less than 5–10 Myr (the typical duration of a star formation episode). We interpret this as follows: at some point, the BH falls into a denser gas clump, which will start to form stars. For a few Myr, the BH will accrete efficiently, until the cloud is dispersed by SN feedback.

By comparison, without SNe (Fig. 5c), the average density around the BH is much higher, typically between 30 and $200 \text{ cm}^{-3}$. In that case, the BH will accrete much more steadily, and as a result the AGN will be in Åquasar mode' for most of time. Because of this smooth energy injection, the gas surrounding the BH is continuously heated and pushed away from the very centre of the galaxy, and thus the BH does not grow at the Eddington limit all of the time. However, the BH sometimes falls into a massive clump where accretion becomes Eddington-limited, this happens for example from $t = 550$ Myr to $t = 620$ Myr in Fig. 5(c). The accretion will then only drop when the AGN feedback has injected enough energy to lower the density around the BH. Interestingly, the formation of this specific clump is not triggered by a galaxy merger
<figure><img src="image_7.png" /><figcaption>**Figure 4.** BH mass versus stellar mass colour-coded by redshift for all our simulations: stars for AGN_SN, circles for AGN_only, crosses for BH_SN and triangles for AGNforced_SN. For better legibility, we only show one point every three snapshots (so every 15 Myr) and use thin lines for the missing snapshots. SN feedback is by far the strongest mechanism to suppress BH growth in our low-mass galaxy.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Density in the vicinity of the most massive BH of each simulation. The background colour indicates the accretion rate, ranging from $10^{-2} M_{\text{Edd}}$ and below in white up to $M_{\text{Edd}}$ in dark grey. Note that the BHs form at slightly different times in the different simulation, as a result of a slightly different star formation history for the host galaxy.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Projected gas distribution for two snapshots of the AGN_only simulation, with the BH position highlighted by a black circle. The left-hand panel corresponds to a phase where the BH is embedded in a dense clump and accretes at the Eddington limit, while in the right-hand panel, the Eddington accretion phase just finished.</figcaption></figure>

(even though the central BH in that simulation does undergo a merger ~30 Myr earlier), but by the close encounter with a nearby satellite.

We illustrate this process in Fig. 6, which shows a projected map of the gas density for two snapshots of the AGN_only simulation, at $t \sim 600$ Myr (left) and $t \sim 640$ Myr (right). The black circle at the centre of the figure depicts the BH position. A quick visual inspection confirms that the drop in density and accretion rate seen in Fig. 5(c) is indeed due to feedback and not dynamics: the BH stays at the centre of the galaxy while a low-density hole is created nearby.

The situation is very different for the AGN_forced_SN simulation (Fig. 5d), where we allow for both AGN and SN feedback. Since the BH is forced to accrete as much gas as is available, the accretion rate is directly tied to the local gas density. At the same time, this artificially boosted accretion will result in an artificially strong AGN feedback, which will reinforce the effect of SNe to reduce the local gas density around the BH. This explains why the average density in Fig. 5(d) is approximately one order of magnitude lower than for the fiducial AGN_SN simulation.

We quantify this further in Fig. 7, which presents the probability distribution function of Eddington ratios $\lambda_{Edd} = \dot{M}/\dot{M}_{Edd}$ for all the simulations. As already expected from Fig. 5, the BH in the AGN_SN simulation (filled histogram in blue) accretes very little. All the other curves show episodes where the BH is accreting at the Eddington limit, but interestingly the AGN_forced_SN simulation (solid grey line) is the one for which these Eddington-limited episodes represent the shortest amount of time. For the BH_SN simulation (red dotted line), the BH accretes very little for most of the time, but experiences some high accretion episodes when it falls in a star-forming clump prior to its destruction via SN feedback. Finally, the BH in the simulation AGN_only spends a large amount of time accreting at the Eddington limit (mostly corresponding to the clump capture at $t \simeq 600$ Myr), but the $\lambda_{Edd}$ distribution has a broad secondary component around $\lambda_{Edd} \sim 10^{-2}$, corresponding to the phases where the BH is not embedded in a very dense clump.

# 4 ESCAPE OF IONIZING RADIATION

## 4.1 Escape of stellar radiation

We now discuss the effect of AGN feedback on the escape of H ionizing radiation from the galaxy where they are produced. We

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Probability distribution function of Eddington ratios for all the simulations: blue-filled histogram for AGN_SN, orange dashed line for AGN_only, red dotted line for BH_SN, and solid grey line for AGN_forced_SN.</figcaption></figure>

measure the fraction $f_{esc}$ of photons that manage to escape the halo following Trebitsch et al. (2017): we integrate the total ionizing flux escaping the halo and compare it to the ionizing photon production rate of the galaxy within the halo. We assume that all the stars are at the centre of the galaxy, and we correct the stellar ages for the light crossing time of the halo, $t_{cross} = R_{vir}/\tilde{c}$, where $R_{vir}$ is the virial radius and $\tilde{c}$ is the reduced speed of light. The photon production rate is then measured using SED models of Bruzual & Charlot (2003).

We show in Fig. 8 the escape fraction $f_{esc}$ as a function of time for all the simulations presented in this work. The BH formation times are annotated as arrows at the top of the figure. All the simulations with SN feedback (black dash-dot-dotted line for SN_only, blue solid line for AGN_SN, red dotted line for BH_SN, and grey solid line for AGN_forced_SN) behave very similarly, with some large spikes of high $f_{esc}$ episodes lasting for up to ~10–20 Myr. The two simulations without AGN feedback (SN_only in black and BH_SN in red) follow each other even more closely, but some differences seem to appear (e.g. a peak of $f_{esc}$ at 680 Myr for SN_only does not exist for BH_SN), mostly due to the intrinsic stochasticity of star formation. This tends to indicate that the mere presence of a BH accreting in the galaxy can slightly change its instantaneous star formation properties.

The fiducial AGN_SN simulation seems to undergo a very high escape fraction episode around the formation time of its BH, around 410 Myr, but it is actually due to a concurrent star formation event (see the jump in stellar mass around that time in Fig. 2 and Appendix B). These two events are of course related: since we impose a threshold on the local stellar density to spawn our sink particles, such a massive star formation event will necessarily result in the apparition of a new BH at that time, if not already present. The beginning of the high $f_{esc}$ episode corresponds to the end of the first (and strongest) BH accretion episode (see Fig. 5a), further suggesting that AGN feedback participates in clearing the ISM. However, the spike in $f_{esc}$ is also present (though reduced by a factor of ~3) for other simulations that include SN feedback. While this suggests that its origin is physical rather than a numerical artefact, it also means that AGN feedback is only helping the radiation to escape at that time. The simulation without SN feedback (AGN_only, the orange dashed line) exhibits a very different time evolution, which is expected based on the results of e.g. Trebitsch et al. (2017). Yet, Fig. 8
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Evolution of the escape fraction $f_{esc}$ with time for the different simulations presented in this work. Arrows mark BH formation times. When SN feedback is included, $f_{esc}$ has a very bursty behaviour, almost independently of the BH accretion and feedback. Without SN feedback however, the high $f_{esc}$ episodes are much less frequent and directly tied to AGN feedback.</figcaption></figure>

features two spikes of very high $f_{esc}$, which indeed correspond to episodes of high accretion rate, e.g. at $t \sim 300$ Myr and $t \sim 370$ Myr (see Fig. 5c). However, these episodes (at $z \sim 13-14$) correspond to a merger between the progenitor of the halo of interest at $z \sim 6$ and another halo that had already formed stars. The high $f_{esc}$ episodes result from this merger, which breaks our assumption that all stars lie at the centre of the galaxy, and strongly disturbs the ISM of both galaxies and therefore facilitate the escape of ionizing radiation. The long episode of very high accretion rate that lasts from $t \sim 550$ Myr to $t \sim 620$ Myr has apparently no effect on $f_{esc}$. It is only at the end of this episode that $f_{esc}$ rises again to around 10 per cent, when AGN feedback finally manages to pierce the dense clump in which the BH is embedded.

Because of the high variability of $f_{esc}$, it is difficult to quantitatively comment the differences between all the simulations including SN feedback, and isolate the effect of BH growth on the contribution of galaxies to reionization. A more practical quantity for that purpose is the total number of H I-ionizing photons produced by the galaxy and available to reionize the IGM. We compute the total number of photons produced as the integral of the photons production rate $\dot{N}_{em}$ between the time of formation of a galaxy $t_{form}$ and the current time $t$:

$$N_{em}(<t) = \int_{t_{form}}^{t} \dot{N}_{em}(t') dt'. \qquad (8)$$

Similarly, the number of photons that escape the halo to contribute to reionization is given by

$$N_{esc}(<t) = \int_{t_{form}}^{t} \dot{N}_{esc}(t') dt' = \int_{t_{form}}^{t} f_{esc}(t') \dot{N}_{em}(t' - t_{cross}) dt', \qquad (9)$$

where $t_{cross} = R_{vir}/\tilde{c}$ is the time it takes for radiation travelling at the reduced speed of light $\tilde{c}$ to reach the virial radius of the halo. The upper panel of Fig. 9 shows $N_{em}$ and $N_{esc}$ (respectively, in thick and thin lines) for the various simulations presented in this work, with the same legend as in Fig. 8. A valuable feature of Fig. 9 compared to

Fig. 8 is that it shows clearly that, provided SN feedback is present, the feedback from the BH has hardly any effect on either the production or the escape of ionizing radiation, for instance $N_{esc}$ seems to converge after $z \sim 9$ ($t \sim 550$ Myr), when the stellar mass of the galaxy reaches $M_{*} \sim 10^{7} M_{\odot}$. The simulation without SN feedback ends up producing more ionizing photons (a factor of 2 more from $z \sim 9$), consistent with a higher stellar mass, but a lower fraction manages to escape the halo. Around $t \sim 410$ Myr, the AGN-SN simulation exhibits a sudden increase in both $N_{em}$ and $N_{esc}$: this explains the high $f_{esc}$ episode already discussed in Fig. 8. The jump in $N_{em}$ is caused by a massive star formation event, which in turn causes a very strong SN feedback episode. In the AGNforced SN run, even though the BH is actively growing and provides feedback continuously, the stellar radiation does not seem to escape more easily than in runs where BH growth is stalled. This suggests that SN feedback remains the dominant mechanism regulating $f_{esc}$ in that galaxy.

In order to facilitate the comparison with reionization models, we show on the lower panel of Fig. 9 the averaged photon production rate (thick lines) and contribution rate (thin lines), defined as $\langle \dot{N}_{em} \rangle = \frac{N_{em}(<t)}{t-t_{form}}$ and $\langle \dot{N}_{esc} \rangle = \frac{N_{esc}(<t)}{t-t_{form}}$. On average, below $z \sim 9$ (when the halo is more massive than $M_{vir} \gtrsim 10^{9} M_{\odot}$ and has a stellar mass above $M_{*} \gtrsim 10^{7} M_{\odot}$), the galaxy contributes $\langle \dot{N}_{esc} \rangle \simeq 10^{51} \text{ s}^{-1}$ when SNs are included, very consistent with the findings of e.g. Kimm & Cen (2014) for haloes of similar mass.

While these two previous quantities are ultimately the most relevant ones to assess the contribution of an individual galaxy to the reionization of its local environment, the average escape fraction $\langle f_{esc} \rangle$ is a more common way to estimate the importance of a given class of galaxies to the reionization of the universe. We present in Fig. 10 the ratio of $N_{esc}$ and $N_{em}$, which gives an effective escape fraction averaged over the history of the galaxy for the different simulations, using the same colour-coding as in Fig. 8. While this has no statistical significance since it is still tied to the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top: Cumulative number of ionizing photons that are produced (thick lines) and that escape (thin lines) the halo as a function of time. Bottom: Averaged production rate (thick lines) and contribution rate (thin lines) of ionizing photons. This essentially reduces to the upper panel divided by the time since the formation of the galaxy. Both panels use the same legend as Fig. 8.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Escape fraction averaged over the lifetime of the galaxy. Legend is the same as Fig. 8. Without SNs, $f_{esc}$ plateaus around 1–2 per cent, except for the initial spike.</figcaption></figure>

past history of one single galaxy, it still gives a reasonable average once the galaxy becomes massive enough such that one single burst of star formation does not change significantly the total stellar mass.

## 4.2 Escape of radiation from the AGN

Even though we did not explicitly track the radiation from the AGN,<fnref n="4" /> it is nonetheless interesting to try to estimate its contribution to the total released ionizing radiation. We can compute this quantity in a very similar fashion $N_{esc}$ discussed in the previous section by estimating an AGN escape fraction $f_{esc}^{AGN}$ and the luminosity in the H I-ionizing band produced by the accretion onto the BH. This luminosity $L_{UV1}$ is directly proportional to the AGN bolometric luminosity

$$L_{UV1} = f_{UV1} L_{bol} = f_{UV1} \epsilon_r \dot{M}_{\bullet} c^2, \quad (10)$$

where the proportionality factor $f_{UV1}$ is given by the AGN spectrum integrated between 13.6 and 24.59 eV, the bounds of the H I-ionizing radiation used with RAMSES-RT. We use the composite AGN spectrum of Sazonov, Ostriker & Sunyaev (2004) for an unobscured AGN, and we find $f_{UV1} \simeq 7.9$ per cent (see also Bieri et al. 2017). The number of photons produced in that band is then simply given by

$$\dot{N}_{\text{em}}^{\text{AGN}} = \frac{L_{\text{UV1}}}{\langle \epsilon \rangle}, \quad (11)$$

where $\langle \epsilon \rangle \simeq 18$ eV is the luminosity-weighted average energy in the band. From this, we can finally compute the total number of escaped photons per unit time:

$$\dot{N}_{\text{esc}}^{\text{AGN}} = f_{\text{esc}}^{\text{AGN}} \dot{N}_{\text{em}}^{\text{AGN}}. \quad (12)$$

We then use a ray-tracing technique using the RASCAS code (Michel-Dansac et al., in preparation) to measure the column density seen by the central BH in order to estimate the escape fraction $f_{esc}^{AGN}$: we cast $N = 100000$ rays from the position of the BH and integrate the H I column density in each direction $N_{H,i}$ up to the virial radius, by assuming that the BH is always at the centre of the halo. We then estimate the average $f_{esc}^{AGN}$ as

$$f_{\text{esc}}^{\text{AGN}} = \frac{1}{N} \sum_{i=1}^{N} e^{-t_i} = \frac{1}{N} \sum_{i=1}^{N} e^{-\sigma_{\text{cross}} N_{\text{H},i}}, \quad (13)$$

where $\sigma_{\text{cross}}$ is the average cross-section within that energy band.

Since the BH is only growing in the AGN_only and AGN-forced_SN simulations, we only compute $f_{esc}^{AGN}$ for these two simulations. We present the results on the top panel of Fig. 11, where we compare the escape fraction for stellar radiation (noted $f_{esc}^{\star}$, in solid and dashed lines) to the escape fraction from the AGN ($f_{esc}^{AGN}$, shaded areas). We use the same colours as Fig. 8 for the different simulations. Interestingly, the two simulations present a very different behaviour.

For the simulation with SN feedback (AGNforced_SN), there is a relatively good correlation between the stellar and AGN escape fraction: most peaks of stellar $f_{esc}$ correspond to a peak of similar duration of $f_{esc}^{AGN}$. This can easily be understood in light of the results of Section 3.2. We have seen that the accretion rate is limited by the amount of available gas in the vicinity of the BH. Whenever the BH is in a clump and accretes efficiently, it will be in competition with star formation to access the local gas reservoir: after a few

<footnote n="4">⁴During the refereeing process, we ran two extra simulations where we followed the radiation from the AGN, corresponding to the AGN_SN and AGN_forced_SN simulations, both using the composite spectrum described here. The results of these simulations, discussed in Appendix C, are consistent with the findings of this section.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top: Instantaneous escape fraction for the radiation produced by the stars (lines) and from the AGN (shaded area) for the AGN_only and AGNforced_SN simulations. For the latter, there is a good correlation between $f_{esc}^*$ and $f_{esc}^{AGN}$, because the gas distribution in the ISM is dominated by SN feedback. Lower left: Ionizing emissivity of the AGN in the AGN_SN, AGN_only and AGNforced_SN simulations, which follows the accretion rate onto the BH. Lower right: Cumulative number of ionizing photons that are produced (thick lines) and that escape the halo (thin lines) as a function of time. All panels use the same legend as Fig. 8.</figcaption></figure>

Myr, when the cloud starts to be disrupted by SN explosions, the column density around the BH drops and $f_{esc}^{AGN}$ increases. The tight correlation between $f_{esc}^*$ and $f_{esc}^{AGN}$ confirms that the AGN has little impact on the galaxy. The ISM is already so stirred by SN feedback that even in the extreme case of forced accretion, AGN feedback cannot disturb it much more (see also Prieto et al. 2017). In this picture, we can explain the relative height of the $f_{esc}^{AGN}$ peaks compared to the peaks in $f_{esc}^*$: while the BH is either inside or outside of a clump, the stellar distribution in the galaxy is more complex. Trebitsch et al. (2017) showed that the stellar escape fraction is almost completely determined by the local properties around a stellar cluster, meaning that either radiation escapes, or it doesn't at all. If the average $f_{esc}^*$ over the galaxy is lower than $f_{esc}^{AGN}$, this indicates that some stars are still embedded in dense, neutral clouds. Further absorption of ionizing photons in the circum-galactic medium

accounts for the fact that the escape fraction ($f_{esc}^*$ or $f_{esc}^{AGN}$) do not reach 100 per cent. For the simulation without SN feedback (AGN_only), the situation is very different: while $f_{esc}^*$ is almost always close to zero, $f_{esc}^{AGN}$ is non-zero most of the times, especially after 600 Myr. This is again consistent with our finding: while stars are trapped in dense clumps with no feedback to create low-density channels, AGN feedback can heat and push the gas around the BH, for example at $t = 620$ Myr on the right-hand panel of Fig. 6. This creates easy escape routes for the ionizing radiation, and results in $f_{esc}^{AGN} \sim 60$ per cent. At the same time, some radiation produced by young stars in the clump will also find low-density channels, carved by AGN feedback, explaining the value of $f_{esc}^* \sim 10$ per cent at $t = 620$ Myr.

We show in the lower left-hand panel of Fig. 11 the instantaneous photon production rate, $\dot{N}_{em}^{AGN}$, resulting from the accretion onto the
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Ratio of the AGN feedback power to the ionizing luminosity emitted by stars and absorbed by the gas for the three simulations with AGN feedback.</figcaption></figure>

mode. The horizontal black line indicates equal energy injection from these two channels. Most of the time (after $t \gtrsim 400$ Myr), AGN feedback transfers much less energy to the gas compared to photoionization. This is especially true for the AGN\_SN run, where the accretion rate is very low.

## 5.3 Contribution of faint AGN to the reionization

In the analysis of Section 4.2, we assumed that the radiation emitted by the AGN followed an energy distribution based on the average non-obscured quasar spectrum of Sazonov et al. (2004). For an accretion disc following the Shakura & Sunyaev (1973) disc model, we expect the peak of the SED to shift towards higher frequency when the Eddington ratio $\lambda_{Edd}$ increases and when the BH mass decreases. When computing the ionizing luminosity of the AGN in equation (10) and more specifically the $f_{UV1}$, we should use a more physically motivated model for the spectrum, for instance the OXAF model of Thomas et al. (2016) based on Done et al. (2012). This is however at most a sub-dominant effect: even if we assume that the ionizing fraction of the bolometric luminosity is 10 times larger ($f_{UV1} \simeq 80$ per cent), the total production of ionizing photons by the AGN would still be two orders of magnitudes below the amount of ionizing photons produced by the stars. The inclusion of secondary ionizations from the hard AGN spectrum does not change this result either, as at most 35–37 per cent of their energy goes into secondary reionizations (with only a weak dependence on the photon energy, and the fraction dropping as the ionized hydrogen fraction increases; Shull & van Steenberg 1985; Venkatesan, Giroux & Shull 2001). This is because the dominant factor here is not so much the AGN spectrum, but the fact that in our most physical simulation (AGN\_SN), the accretion rate is low at all times, yielding a low bolometric luminosity.

The direct consequence of this is that BHs in the low-mass regime are unlikely to contribute much to the ionizing budget of the EoR. It is of course possible that we are underestimating the accretion onto the BH we are simulating here, for instance because our $\lesssim 10$ pc resolution can only resolve structures down to that scale, and we do not correct for any unresolved subgrid clumping of the gas. For instance, Negri & Volonteri (2017) discuss how the Bondi prescription for accretion can lead to both over- and under-estimation of the real accretion rate. However, even in our AGN-forced\_SN simulation (which corresponds to an extreme model

where the both the galaxy and the BH grow reasonably), the ionizing contribution from the AGN is negligible compared to the ionizing starlight. Going beyond this model and increasing even more the accretion rate would result in an overmassive BH compared to its host galaxy.

# 6 SUMMARY AND CONCLUSIONS

Motivated by the recent suggestions that faint AGN could play an important role in the reionization of the universe, we performed a series of cosmological zooms on a dwarf galaxy ($M_{vir} \simeq 5 \times 10^9 M_{\odot}$ at $z \simeq 6$) using a high resolution of $\sim 10$ pc, fully coupled radiative hydrodynamics simulation, and following the growth of massive BHs and the associated AGN feedback. In our reference run, we use a dual-mode AGN feedback describing both the radio and quasar modes, and we account for SN feedback. We compare this reference run to simulations where we alternately turned off SN and AGN feedback and we also tried an alternative prescription where the BH is forced to accrete as much gas as possible. The main goal of our study was to determine how and how much BH can affect the ISM of low-mass galaxies, and more specifically how this affects their contribution to the ionizing budget of the EoR. We summarize our main results below:

(i) SN feedback strongly hampers BH growth in low-mass galaxies. The resulting constant stirring of the ISM removes gas from around the BH, resulting in a very low accretion rate. Even when the BH is forced to accrete, the actual accretion is limited by the amount of gas locally available, and the BH grows at sub-Eddington rates.

(ii) BH growth has no effect on the escape of ionizing radiation. In the run without SN feedback, the integrated $f_{esc}$ is much lower than that in all other runs. There is no sign of amplification of the integrated $f_{esc}$ by AGN feedback, even in runs where the BH grows.

(iii) The instantaneous $f_{esc}$ is virtually unaffected by the presence of a BH or by its growth: AGN feedback does not provide an extra source of energy injection in the ISM, and therefore does not help ionizing radiation to escape before the destruction of the star-forming clouds by SN feedback.

(iv) When the BH grows in presence of SN feedback (AGN-forced\_SN simulation), the ionizing radiation produced by the accretion onto the BH escapes synchronously with the stellar ionizing radiation.

(v) Only a fraction $f_{esc}^{AGN}$ of the ionizing luminosity produced by the gas accretion onto the BH manages to escape in the IGM to contribute to reionization, with $f_{esc}^{AGN} < 50$ per cent and on average $\langle f_{esc}^{AGN} \rangle \simeq 10$ per cent.

# ACKNOWLEDGEMENTS

We are grateful to the anonymous referee for a constructive report that improved the manuscript. We thank Jérémy Blaizot and Léo Michel-Dansac for providing us with the RASCAS machinery. 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 Stephane Rouberol for running smoothly this cluster for us. This work was granted access to the HPC resources of CINES under the allocation A0020406955
made by GENCI. This work has made extensive use of the PYMSES<fnref n="5" /> analysis package (Guillet, Chapon & Labadens 2013), as well as the MATPLOTLIB (Hunter 2007), NUMPY/SCIPY (Jones et al. 2001) and IPYTHON (Perez & Granger 2007) packages.

## REFERENCES

Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109

Atek H. et al., 2015, ApJ, 814, 69

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Baldassare V. F., Reines A. E., Gallo E., Greene J. E., 2015, ApJ, 809, L14

Becker R. H. et al., 2001, AJ, 122, 2850

Becker G. D., Bolton J. S., Madau P., Pettini M., Ryan-Weber E. V., Vemmans B. P., 2015, MNRAS, 447, 3402

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

Bergvall N., Leitet E., Zackrisson E., Marquart T., 2013, A&A, 554, A38

Bieri R., Dubois Y., Rosdahl J., Wagner A., Silk J., Mamon G. A., 2017, MNRAS, 464, 1854

Bondi H., 1952, MNRAS, 112, 195

Booth C. M., Schaye J., 2009, MNRAS, 398, 53

Bower R. G., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2017, MNRAS, 465, 32

Bridge C. R. et al., 2010, ApJ, 720, 465

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Chabrier G., 2003, PASP, 115, 763

Chapon D., Mayer L., Teyssier R., 2013, MNRAS, 429, 3114

Chen C.-T. J. et al., 2017, ApJ, 837, 48

Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015, MNRAS, 450, 1514

Cristiani S., Serrano L. M., Fontanot F., Vanzella E., Monaco P., 2016, MNRAS, 462, 2478

Deharveng J.-M., Buat V., Le Brun V., Milliard B., Kunth D., Shull J. M., Gry C., 2001, A&A, 375, 805

Domínguez A., Siana B., Brooks A. M., Christensen C. R., Bruzual G., Stark D. P., Alavi A., 2015, MNRAS, 451, 839

Done C., Davis S. W., Jin C., Blaes O., Ward M., 2012, MNRAS, 420, 1848

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

Dubois Y., Pichon C., Devriendt J., Silk J., Haehnelt M., Kimm T., Slyz A., 2013, MNRAS, 428, 2885

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

Dubroca B., Feugeas J., 1999, Acad. Sci. Paris Comptes Rendus Ser. Sci. Math., 329, 915

Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267

Fan X. et al., 2001, AJ, 122, 2833

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Ferland G. J. et al., 2017, RMXAA, 53, 385

Filippenko A. V., Sargent W. L. W., 1989, ApJ, 342, L11

Fontanot F., Cristiani S., Vanzella E., 2012, MNRAS, 425, 1413

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248

Giallongo E. et al., 2015, A&A, 578, A83

Gnedin N. Y., 2016, ApJ, 825, L17

Grazian A. et al., 2017, A&A, 602, A18

Grissom R. L., Ballantyne D. R., Wise J. H., 2014, A&A, 561, A90

Guillet T., Chapon D., Labadens M., 2013, PyMSES: Python modules for RAMSES, Astrophysics Source Code Library, record ascl:1310.002

Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633

Haardt F., Madau P., 1996, ApJ, 461, 20

Haardt F., Madau P., 2012, ApJ, 746, 125

Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935

Hahn O., Abel T., 2011, MNRAS, 415, 2101

Hassan S., Davé R., Mitra S., Finlator K., Ciardi B., Santos M. G., 2018, MNRAS, 473, 227

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73

Izotov Y. I., Schaerer D., Thuan T. X., Worseck G., Guseva N. G., Orlitová I., Verhamme A., 2016, MNRAS, 461, 3683

Japelj J. et al., 2017, MNRAS, 468, 389

Jones E. et al., 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/

Kaaret P., Brorby M., Casella L., Prestwich A. H., 2017, MNRAS, 471, 4234

Kauffmann J., Pillai T., Goldsmith P. F., 2013, ApJ, 779, 185

Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531

Kimm T., Cen R., 2014, ApJ, 788, 121

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, 4826

Kroupa P., Weidner C., 2003, ApJ, 598, 1076

Kuhlen M., Faucher-Giguère C.-A., 2012, MNRAS, 423, 862

Kunth D., Sargent W. L. W., Bothun G. D., 1987, AJ, 93, 29

Leitherer C., Ferguson H. C., Heckman T. M., Lowenthal J. D., 1995, ApJ, 454, L19

Levermore C. D., 1984, J. Quant. Spec. Radiat. Transf., 31, 149

Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113

Ma X., Kasen D., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Kereš D., Murray N., 2015, MNRAS, 453, 960

Madau P., Haardt F., 2015, ApJ, 813, L8

Madau P., Rees M. J., Volonteri M., Haardt F., Oh S. P., 2004, ApJ, 604, 484

Mao J., Kim M., 2016, ApJ, 828, 96

McAlpine S., Bower R. G., Harrison C. M., Crain R. A., Schaller M., Schaye J., Theuns T., 2017, MNRAS, 468, 3395

Mezcua M., Civano F., Fabbiano G., Miyaji T., Marchesi S., 2016, ApJ, 817, 20

Micheva G., Iwata I., Inoue A. K., 2017, MNRAS, 465, 302

Miller S. H., Ellis R. S., Newman A. B., Benson A., 2014, ApJ, 782, 115

Moran E. C., Shahinyan K., Sugarman H. R., Vélez D. O., Eracleous M., 2014, AJ, 148, 136

Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822

Negri A., Volonteri M., 2017, MNRAS, 467, 3475

O'Shea B. W., Wise J. H., Xu H., Norman M. L., 2015, ApJ, 807, L12

Ocvirk P. et al., 2016, MNRAS, 463, 1462

Onoue M. et al., 2017, ApJ, 847, L15

Ostriker E. C., 1999, ApJ, 513, 252

Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Padoan P., Haugbølle T., Nordlund Å., Frimann S., 2017, ApJ, 840, 48

Parsa S., Dunlop J. S., McLure R. J., 2018, MNRAS, 474, 2904

Perez F., Granger B. E., 2007, Comput. Sci. Eng., 9

Planck Collaboration XIII, 2016a, A&A, 594, A13

Planck Collaboration XLVI, 2016b, A&A, 596, A107

Prieto J., Escala A., Volonteri M., Dubois Y., 2017, ApJ, 836, 216

Qin Y. et al., 2017, MNRAS, 472, 2009

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Reines A. E., Comastri A., 2016, PASA, 33, e054

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Reines A. E., Greene J. E., Geha M., 2013, ApJ, 775, 116

Ricci F., Marchesi S., Shankar F., La Franca F., Civano F., 2017, MNRAS, 465, 1915

Robertson B. E. et al., 2013, ApJ, 768, 71

Robertson B. E., Ellis R. S., Dunlop J. S., McLure R. J., Stark D. P., 2010, Nature, 468, 49

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Rosdahl J., Schaye J., Teyssier R., Agertz O., 2015, MNRAS, 451, 34

<footnote n="5">⁵https://bitbucket.org/dchapon/pymses/</footnote>
Rosdahl J., Schaye J., Dubois Y., Kimm T., Teyssier R., 2017, MNRAS, 466, 11

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

Salmon B. et al., 2015, ApJ, 799, 183

Salpeter E. E., 1955, ApJ, 121, 161

Santini P. et al., 2017, ApJ, 847, 76

Sazonov S. Y., Ostriker J. P., Sunyaev R. A., 2004, MNRAS, 347, 144

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shull J. M., van Steenberg M. E., 1985, ApJ, 298, 268

Siana B. et al., 2015, ApJ, 804, 17

Simpson C. M., Bryan G. L., Johnston K. V., Smith B. D., Mac Low M.-M., Sharma S., Tumlinson J., 2013, MNRAS, 432, 1989

Song M. et al., 2016, ApJ, 825, 5

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Teyssier R., 2002, A&A, 385, 337

Thomas A. D., Groves B. A., Sutherland R. S., Dopita M. A., Kewley L. J., Jin C., 2016, ApJ, 833, 266

Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

van Wassenhove S., Volonteri M., Walker M. G., Gair J. R., 2010, MNRAS, 408, 1139

Vanzella E. et al., 2016, ApJ, 825, 41

Vasei K. et al., 2016, ApJ, 831, 38

Venkatesan A., Giroux M. L., Shull J. M., 2001, ApJ, 563, 1

Vito F. et al., 2016, MNRAS, 463, 348

Volonteri M., 2010, A&A Rev., 18, 279

Volonteri M., Gnedin N. Y., 2009, ApJ, 703, 2113

Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079

Volonteri M., Reines A., Atek H., Stark D. P., Trebitsch M., 2017, ApJ, 849, 155

Weigel A. K., Schawinski K., Treister E., Urry C. M., Koss M., Trakhtenbrot B., 2015, MNRAS, 448, 3167

Willott C. J. et al., 2010, AJ, 139, 906

Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L., Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560

Yoshiura S., Hasegawa K., Ichiki K., Tashiro H., Shimabukuro H., Takahashi K., 2017, MNRAS, 471, 3713

<figure><img src="image_16.png" /><figcaption>**Figure A1.** Stellar-to-halo mass relation for our simulated galaxy, using the stellar mass directly measured in the simulation. The solid, dashed, and dotted lines correspond respectively to 100 per cent, 10 per cent, and 1 per cent of the baryons converted into stars. For the last 500 Myr of the simulation, the galaxy converts 10 per cent of its baryons in stars.</figcaption></figure>

that derived from the UV luminosity of the galaxy, we compare in Fig. A2 both stellar mass estimates. The open symbols denote snapshots at $z > 8$, where the log $M_\star - M_{1500}$ relation from Song et al. (2016) is extrapolated rather than interpolated. The observational estimate of $M_\star$ underestimates by sometimes one order of magnitude the actual stellar mass of the galaxy, and taking into account flux losses and dust would reduce even more the UV-based stellar mass. The large scatter in the UV-based $M_\star$ at fixed actual mass results from the burstiness of the star formation history that translate in a UV magnitude varying more rapidly than the stellar mass (see for instance Domínguez et al. 2015).

## APPENDIX A: INTRINSIC VERSUS OBSERVED GALAXY PROPERTIES

In Section 3.1, we discussed the assembly history of our simulated galaxy based on properties measured using its inferred UV magnitude. We now wish to discuss the impact of this approach, and how it compares to quantities directly measured in the simulation.

### A1 Stellar-to-halo mass relation

We show in Fig. A1 the stellar-to-halo mass relation, similar to Fig. 2, but this time using the stellar mass directly measured in the simulation. At very early times, the galaxy converts less than 10 per cent of its baryons into stars, and this fraction saturates after 500 Myr ($z \sim 9.5$). When the true stellar mass is used, our simulation lies well above the Behroozi et al. (2013) model, but gives results very similar to other simulations including similar subgrid models at equivalent spatial resolution, such as the *Renaissance* simulation (O’Shea et al. 2015) shown here in pink, with the fit to their stellar-to-halo mass relation shown as a solid line and an extrapolation to higher masses as a dashed line.

To further investigate the origin of the apparent discrepancy between the stellar mass directly measured in the simulation and

### A2 Star formation MS

In Fig. A3, we present the star formation MS using stellar masses and instantaneous SFR directly from the simulation. As expected, the scatter is much larger than with UV-derived properties. As star formation proceeds by bursts in low-mass galaxies, the SFR can range from very low to very high at fixed stellar mass, depending on whether the galaxy is caught during or after a star formation episode. This scatter is reduced when measuring the SFR from the UV luminosity as in Fig. 3, since the far UV at $\lambda \sim 1500$ Å traces timescales typically longer than the interval between two bursts of star formation.

We confirm this in Fig. A4, which shows again the MS for the AGN SN simulation, but using this time the direct estimate of the stellar mass and the SFR estimate based on Kennicutt & Evans (2012). The coloured lines correspond to the same observations as in Fig. 3. At fixed stellar mass, the scatter is strongly reduced compared to Fig. A3. By comparison, Fig. A5 shows the same MS, but with the stellar mass derived using the relation from Song et al. (2016) and the SFR directly measured in the simulation. This can be thought of as an SFR measured from emission lines such as $H\alpha$, which trace very recent star formation. There, the scatter in stellar mass is large and originates in the burstiness of the SFR. This
<figure><img src="image_17.png" /><figcaption>**Figure A2.** Stellar mass recovered from the UV magnitude using the Song et al. (2016) empirical relation versus the measured stellar mass in the simulation. The log $M_\star - M_{1500}$ is extrapolated above $z > 8$, as indicated by open symbols.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure A3.** Star formation MS, similar to Fig. 3, but with the stellar masses and SFRs directly measured in the simulation. The large scatter in the simulation is due to the intrinsic burstiness of star formation.</figcaption></figure>

explains why the simulation points match the range of observations relatively well.

This all suggests that great care should be taken when comparing simulations to high-redshift observations, probably even more than presented here, since we did not account for dust or nebular emission. Conversely, this points to the expectation that there should not generically be a perfect match between simulated quantities such as the stellar mass or the SFR directly measured in simulations and those derived from high-redshift observations, especially in the low-mass regime.

<figure><img src="image_19.png" /><figcaption>**Figure A4.** Star formation MS, similar to Fig. A3, but with only the stellar masses directly measured in the simulation.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure A5.** Star formation MS, similar to Fig. A3, but with only the SFR directly measured in the simulation.</figcaption></figure>

## APPENDIX B: STELLAR MASS ASSEMBLY

In this appendix, we present the star formation history of the main galaxy in all the runs discussed in this work. We compare the evolution of the instantaneous SFR in Fig. B1 and the evolution of the stellar masses in Fig. B2. The quantities presented here are those directly measured in the simulation, rather than what one would derive using the UV luminosity (see Section 3.1). Despite a large variability of the instantaneous SFR, it seems that the presence of an AGN has very little effect on the final stellar mass. Indeed, even though the galaxy in the run without AGN feedback (BH_SN) ends up with a slightly higher stellar mass than the reference run (AGN_SN), it also has roughly the same stellar mass as the run where the AGN is actually active (AGN_forced_SN). As expected for low-mass galaxies, SN feedback is the main regulator of star formation: the simulation without SN (AGN_only) forms many more stars. The sudden jump in stellar mass for the AGN_only simulation around $t = 400$ Myr is due to the mis-identification of the galaxy progenitor's host halo.
<figure><img src="image_21.png" /><figcaption>**Figure B1.** Evolution of the SFR for all the simulations presented in this work. The legend is the same as in Fig. 8.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure B2.** Stellar mass as a function of time for all the simulations presented in this work. The legend is the same as in Fig. 8.</figcaption></figure>

## APPENDIX C: RADIATION FROM THE AGN

In this appendix, we explore in more detail the impact of including the ionizing radiation produced by the AGN directly in the simulation. We performed two additional runs, AGNRT_SN and AGNRTforced_SN, which include the same subgrid models as the AGN_SN and AGNforced_SN runs (respectively), but further track radiation produced by the AGN itself. For consistency with Section 4.2, we emit radiation from the accreting SMBH using the Sazonov et al. (2004) spectrum. We stopped the AGNRT-forced_SN run after 700 Myr of evolution ($z \sim 7.5$).

We show in Fig. C1 the ionizing production in the two pairs of runs split in two components: stellar radiation on the upper panel, and AGN radiation on the lower panel. We can see that including AGN radiation on the fly has only a minor impact on the production of (stellar) ionizing photons. As expected, the difference is a bit more important for the forced accretion case (in grey, with a solid

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure C1.** Cumulative number of ionizing photons emitted by stars (top) and AGN (bottom) in the run as a function of time. Solid (dotted) lines indicate the runs without (with) radiation from the AGN.</figcaption></figure>

line for the AGNforced_SN run and dotted line for the AGNRT-forced_SN run. When focusing on the actual radiation produced by the AGN (lower panel), we see that there is a very good match between the AGNforced_SN run, where AGN radiation is only considered in post-processing, and the AGNRTforced_SN run, where it is fully tracked. There is a slightly larger difference between the AGN_SN and AGNRT_SN runs, but it is due to the fact that in these simulations, the SMBH grows very little. Indeed, a total mass accretion of $10^3 M_{\odot}$ corresponds to a production of less than $5 \times 10^5$ hydrogen-ionizing photons. Thus, the small difference between the AGN_SN and AGNRT_SN runs is directly due to a small difference in the SMBH accretion history.

Finally, we show in Fig. C2 the total number of ionizing photons escaping the halo as a function of time, for the same two pairs of runs. Here again, the differences between the runs with and without AGN radiation are very limited once the galaxy is settled enough, and are mostly driven by the stochasticity of the star formation history.
<figure><img src="image_25.png" /></figure>

**Figure C2.** Total cumulative number of ionizing photons escaping the halo as a function of time. Solid (dotted) lines indicate the runs without (with) radiation from the AGN.

This paper has been typeset from a TeX/LA TeX file prepared by the author.