How AGN feedback drives the size growth of the first quasars

TL;DR

This study uses cosmological simulations to show that AGN feedback causes the expansion of quasar host galaxies' stellar bulges, influencing their size growth and internal structure from redshift 6 to 3.

Contribution

It demonstrates that AGN feedback-driven outflows are key to the size growth and core formation in massive galaxies, a novel insight into galaxy evolution.

Findings

Compact stellar bulges expand significantly over time.

AGN feedback fluctuations induce potential changes driving expansion.

Size growth and core formation are suppressed without AGN feedback.

Abstract

Quasars at $z \,=\, 6$ are powered by accretion onto supermassive black holes with masses $M_{\rm BH} \sim 10^9 \rm \, M_{\odot}$. Their rapid assembly requires efficient gas inflow into the galactic nucleus, sustaining black hole accretion at a rate close to the Eddington limit, but also high central star formation rates. Using a set of cosmological 'zoom-in' hydrodynamic simulations performed with the moving mesh code Arepo, we show that $z \,=\, 6$ quasar host galaxies develop extremely tightly bound stellar bulges with peak circular velocities $300$ - $500$ km s$^{-1}$ and half-mass radii $\approx 0.5 \, \rm kpc$. Despite their high binding energy, we find that these compact bulges expand at $z \, < \, 6$, with their half-mass radii reaching $ \approx 5$ kpc by $z \, = \, 3$. The circular velocity drops by factors $\approx 2$ from their initial values to $200$ - $300$ km s$^{-1}$ at $z \, \approx \, 3$ and the stellar profile undergoes a cusp-core transformation. By tracking individual stellar populations, we find that the gradual expansion of the stellar component is mainly driven by fluctuations in the gravitational potential induced by bursty AGN feedback. We also find that galaxy size growth and the development of a cored stellar profile does not occur if AGN feedback is ineffective. Our findings suggest that AGN-driven outflows may have profound implications for the internal structure of massive galaxies, possibly accounting for their size growth, the formation of cored ellipticals as well as for the saturation of the $M_{\rm BH}$ - $\sigma_{\star}$ seen at high velocity dispersions $\sigma_{\star}$.

Tables

Table 1. Table 1: Main numerical parameters in our simulations. The first, second and third column give the mass of high-resolution dark matter particles, the target mass of gas cells and the mass of stellar particles, respectively. The gravitational softening length is given in comoving units in the fourth column. In the redshift range investigated here, this corresponds to 150 - 400 physical pc. The minimum cell size is given in the fifth column.

$m_{\mathrm{dm}}$	$m_{\mathrm{gas}}$	$m_{\star }$	${ϵ}_{\mathrm{soft}}$	$\Delta x$
($h^{-1}{M}_{\odot }h$)	($h^{-1}{M}_{\odot }$)	($h^{-1}{M}_{\odot }$)	($h^{-1}$ckpc)	(pc)
$6.8\times {10}^6$	$1.8\times {10}^6$	$1.3\times {10}^5$	1	22

Table 2. Table 2: Main properties of the sample of the 6 haloes followed by our simulations, evaluated at z = 6.2 𝑧 6.2 z=6.2 and z = 3.3 𝑧 3.3 z=3.3 in, respectively, the top and bottom table. We list the mass of the most massive black hole present in each volume (second column), the total mass within the virial radius (third column), virial radius (fourth column), the virial temperature (fifth column) of its parent halo, the total star formation rate within the virial radius (sixth column) and the stellar half mass radius as predicted by SUBFIND (seventh column).

Simulation at $z=6.2$	$M_{\mathrm{BH}}$	$M_{200}$	$R_{200}$	$T_{200}$	$SF{R}_{200}$	$V_{200}$	$R_{\mathrm{eff}}$
	($\times {10}^8{M}_{\odot }$)	($\times {10}^{12}{M}_{\odot }$)	(kpc)	($\times {10}^6$ K)	($M_{\odot }{\mathrm{yr}}^{-1}$)	($\mathrm{km}$ $s^{-1}$)	(kpc)
Halo1	8.81	3.97	70.0	4.01	352.82	494.14	1.16
Halo1-NoAGN	-	3.97	70.0	4.24	396.59	494.02	1.05
Halo2	9.12	1.81	59.85	4.03	215.0	422.61	0.80
Halo3	17.8	3.73	68.54	3.25	237.56	483.92	0.59
Halo4	7.43	3.18	64.96	3.24	520.49	458.68	0.91
Halo5	9.66	1.95	55.21	3.02	70.81	389.78	0.76

Equations

$${\dot{M}}_{\mathrm{BH}}=\min \left[\frac{4\pi \alpha G^2{M}_{\mathrm{BH}}^2{\rho }_{\mathrm{gas}}}{(c_{\mathrm{s}}^2+v^2{)}^{3/2}},{\dot{M}}_{\mathrm{edd}}\right],$$

$${\dot{M}}_{\mathrm{edd}}=\frac{4\pi G{M}_{\mathrm{BH}}{m}_{\mathrm{p}}}{{ϵ}_{\mathrm{r}}{\sigma }_{\mathrm{T}}c},$$

$${\dot{E}}_{\mathrm{feed}}={ϵ}_f{L}_{\mathrm{bol}}={ϵ}_{\mathrm{f}}{ϵ}_{\mathrm{r}}(1-{ϵ}_{\mathrm{r}}{)}^{-1}{\dot{M}}_{\mathrm{BH}}{c}^2,$$

$${\kappa }_{\mathrm{rot}}=\frac{K_{\mathrm{rot}}}{K}=\frac{1}{K}\sum _i^{r<30\mathrm{kpc}}\frac{1}{2}{m}_i[L_{z,i}/m_i{R}_i{]}^2,$$

$$V_{\mathrm{circ}}=\sqrt{\frac{G{M}_{\mathrm{stars}}(<R)}{R}},$$

Full text

††pubyear: 2019††pagerange: How AGN feedback drives the size growth of the first quasars–How AGN feedback drives the size growth of the first quasars

1 Introduction

Detections of bright quasars at $z\,>\,6$ show that supermassive black holes with masses $M_{\rm BH}\,\gtrsim\,10^{9}\,\rm M_{\odot}$ assemble in less than a Gyr after the Big Bang (e.g. Fan et al., 2004; Willott et al., 2010; Mortlock et al., 2011; Venemans et al., 2013; Gallerani et al., 2017). The most distant quasar known has $z\,=\,7.5$ (Bañados et al., 2018), while the most massive black hole at $z>6$ has $M_{\rm BH}\gtrsim 10^{10}\,\rm M_{\odot}$ (Wu et al., 2015). The extreme conditions which are likely required to account for such rapid black hole growth provide interesting testbeds for models of black hole accretion and its associated active-galactic-nucleus (AGN) feedback.

If the main growth channel is accretion of interstellar gas, the black hole accretion rate must be close to the Eddington limit for a significant fraction of the Hubble time at $z>6$ even if black hole seeds are massive (e.g Sijacki et al., 2009). For lower mass seed black holes, prolonged episodes of super-Eddington accretion rates are likely required to achieve rapid enough black hole growth. Such constraints require supermassive black holes to reside in galaxies where gas in-fall towards their nuclei is sustained, a condition which is most easily satisfied if their host dark matter haloes are massive with $M_{\rm DM}\,\gtrsim\,10^{12}\,\rm M_{\odot}$ (Efstathiou & Rees, 1988; Volonteri & Rees, 2006; Costa et al., 2014).

Cosmological hydrodynamic simulations following black hole growth as well as the impact of energy deposition by the quasar (AGN feedback) have succeeded in forming massive black holes with $M_{\rm BH}\,\gtrsim\,10^{9}\,\rm M_{\odot}$ by $z\,=\,6$ (Sijacki et al., 2009; Di Matteo et al., 2012; Costa et al., 2014). The efficient in-fall of gas into the galactic nucleus necessary to sustain black hole growth, however, also likely triggers powerful central starbursts and the formation of compact stellar bulges characterised by extreme stellar velocity dispersions $\sigma_{\rm*}$. In their high resolution cosmological simulations, Dubois et al. (2012) find galaxies with stellar bulges with $\sigma_{\rm*}>600\,\rm km\,s^{-1}$ measured within $\approx 500\,\rm pc$ embedded within haloes with virial mass $M_{\rm vir}\,\gtrsim\,10^{12}\,\rm M_{\odot}$ at $z>6$. The compact sizes are attributed to their special location at the intersection between multiple large-scale cold gas filaments, a configuration which leads to efficient angular momentum cancellation.

Extremely compact galaxies appear to be the general outcome of galaxy formation in massive haloes at $z\,\gtrsim\,6$. Costa et al. (2015) show, using similar simulations, that ’compact bulges’ with peak circular velocities $\gtrsim 500\,\rm km\,s^{-1}$ form even in the presence of supernova and AGN feedback. Using $\sim 10\,\rm pc$ resolution simulations, Curtis & Sijacki (2016) find peak circular velocities as high as $\approx 900\,\rm km\,s^{-1}$ for the ’compact bulge’ at $z\,\approx\,5$, while Di Matteo et al. (2017) show that extremely compact galaxies assemble already by $z\,=\,8$.

Based on the extreme binding energy of these systems, we are led to question whether they may be able to survive down to lower redshift. Using the Illustris simulations, Wellons et al. (2016) identify various evolutionary paths for high-z massive, compact galaxies. These include (i) the survival of the core within a subsequently accreted stellar envelope, (ii) its destruction in a major merger, (iii) its consumption by a more massive galaxy and (iv) its undisturbed survival down to $z\,=\,0$. In addition, Wellons et al. (2016) find that the compact galaxies grow both in mass and size, which they attribute to accretion of ex-situ stars, a process identified in earlier work (Naab et al., 2009; Oser et al., 2010; Johansson et al., 2012). According to these studies, most descendants of high-z massive, compact galaxies are no longer compact at $z\,=\,0$. Whether and how this picture may change for the extreme systems which form in high-$\sigma$ peaks at $z>6$, however, remains unexplored.

Previous studies have proposed quasar feedback itself as a mechanism for galactic growth (e.g. Fan et al., 2008). In this scenario, successive episodes of rapid ejection of gas from a massive galaxy lead to gravitational potential fluctuations and the expansion of the collisionless component (see also Navarro et al. 1996; Ragone-Figueroa & Granato 2011; Teyssier et al. 2013; Pontzen & Governato 2014; Ogiya & Mori 2014; Penoyre et al. 2017). Others have proposed ‘positive feedback’ as a possible channel for massive galaxy growth (Ishibashi et al., 2013).

In this paper, we investigate the fate of the compact galaxies hosting bright $z>6$ quasars down to $z\,=\,3.3$, using a suite of zoom-in simulations targeting massive $M_{\rm vir}\gtrsim 3\times 10^{12}\,\rm M_{\odot}$ haloes at $z\,=\,6$. We start by describing our suite of simulations in Section 2 and then present our main findings in Section 3. The implications of our results are discussed in Section 4 and our main conclusions are summarised in Section 5.

2 Simulations

We perform cosmological, ‘zoom-in’, hydrodynamic simulations targeting the five most massive haloes found in the dark matter-only Millennium simulation (Springel et al., 2005) at $z=6.2$ and follow their evolution down to $z\,=\,3$. These haloes represent high-$\sigma$ fluctuations of the cosmic density field (see Table 2 for a list of their main properties) and are massive enough to ensure rapid black hole growth (Costa et al., 2014). We assume a $\Lambda$CDM cosmology with parameters $h=0.73$, $\Omega_{m}=0.25$, $\Omega_{\Lambda}=0.75$, $\Omega_{b}=0.041$, identical to the parent Millennium simulation, but a lower $\sigma_{8}$ value of $0.8$, in better accordance with the Planck cosmology (Planck Collaboration et al. 2014).

2.1 Numerical Setup

The numerical setup of our simulations is identical to that of Costa et al. (2015) and is here only briefly summarised.

Our cosmological simulations are performed with the moving-mesh code Arepo (Springel, 2010) and consist of ‘zoom-ins’ of the most massive haloes found in the Millennium volume (Springel et al., 2005) at $z\,=\,6$.

Besides the collisionless dynamics of dark matter and stellar populations, our simulations follow the hydrodynamic evolution of gas on an unstructured Voronoi mesh using a second-order finite volume scheme. Explicit refinement and de-refinement is employed in order to maintain the masses of gas cells within a factor of two of a target mass, here chosen to be $m_{\rm gas}\,=\,1.8\times 10^{6}\,h^{-1}\rm M_{\odot}$ (see Table 1 for a list of the main numerical parameters of our simulations). The geometry and size of the fluid patches is adapted to higher density regions thereby concentrating resolution elements in the densest regions. For example, in Halo1 cells can be as small as 22 pc at $z\,=\,6.2$ and $\approx$ 1.5% of the cells within the virial radius of Halo1 can be classified as ‘small cells’ ($\leq 2\times\Delta x$, see Table 1). Individual high-resolution dark matter particles have a mass $m_{\rm dm}\,=\,6.8\times 10^{6}h^{-1}\,\rm M_{\odot}$ and the stellar particles have $m_{\star}\,=\,1.3\times 10^{5}h^{-1}\,\rm M_{\odot}$. We use a fixed comoving gravitational softening $\epsilon_{\rm soft}\,=\,1h^{-1}\,\rm ckpc$, which means their physical gravitational softening increases with time. We also perform a simulation, in which we maintain a constant physical gravitational softening below $z\,=\,6$. For the gas component, the gravitational softening length scales with the gas cell size, but is never allowed to be smaller than the collisionless gravitational softening.

We follow radiative cooling for a hydrogen and helium gas in collisional ionization equilibrium as in Katz et al. (1996), subjected to a spatially constant, but time-dependent UV background with reionization occurring at $z\approx 6$ (Faucher-Giguère et al. 2009). Metal-line cooling is treated as in Vogelsberger et al. (2013) with the minimum temperature down to which gas is allowed to cool radiatively set to $10^{3}$ K.

We follow Springel & Hernquist (2003) to model star formation and associated feedback. Stellar particles are spawned stochastically out of gas particles with number density exceeding a threshold $n_{0}\,=\,0.13\,\rm cm^{-3}$ and with temperatures $T<10^{5}\,\rm K$. Star-forming gas follows an effective equation of state resulting from the balance of radiative cooling and cold cloud evaporation within the unresolved interstellar medium.

Supernova-driven winds are modelled following Springel & Hernquist (2003); Vogelsberger et al. (2013), i.e. a gas is launched from star forming regions at a rate of $\dot{M}=\eta\dot{M}_{\star}$, where $\dot{M}_{\star}$ is the star formation rate and $\eta$ is the mass-loading factor. Gas cells ejected from a star-forming region are converted into wind particles and launched isotropically. We assume a mass loading factor of $\eta=1$ and kick-velocity of $1200\,\rm km\,s^{-1}$ (as in Costa et al. 2015). Wind particles interact with other matter components only gravitationally until they travel to a region where the local density drops below $0.05n_{\rm 0}$ or a maximum travel time is reached. When either of these criteria are met, the wind particle’s mass, metalicity, momentum and energy are deposited into its host gas cell.

Black holes are modelled as sink particles. Black hole accretion occurs at scales much smaller than cosmological simulations typically can resolve. Resolving the radius at which accretion takes place requires pc resolution which is $\sim\,2$ orders of magnitude lower than the resolution of our simulations. This means we cannot resolve transfer of angular momentum, gas collapse or even the structure of infalling material. The unresolved physics of accretion is therefore also modelled in a subgrid fashion following Springel et al. (2005) and Di Matteo et al. (2005). The accretion is assumed to be of a Bondi-Hoyle-Lyttleton type. Assuming this type of prescription has the advantage that it is a relatively simple parametrisation of black hole accretion. In addition, this type of accretion is commonly assumed in almost every state-of-the-art cosmological simulation (e.g., in EAGLE (Schaye et al., 2015), HORIZON-AGN (Dubois et al., 2014),ILLUSTRIS-TNG (Weinberger et al., 2017), MASSIVE-BLACK (Di Matteo et al., 2012), BLUE-TIDES (Feng et al., 2016)) with the exception of the SIMBA simulations (Davé et al., 2019). This allows us to connect our results with other studies.

The accretion rate is capped at the Eddington rate for each black hole. The accretion rate is given by

[TABLE]

where $\alpha$ is a dimensionless parameter, set to a value of 100, that is introduced to recover a volume coverage of the Bondi rates for the cold and hot interstellar medium phases, $\rho_{\rm gas}$ and $c_{\rm s}$ are the density and sound speed of the gas surrounding the black hole, respectively, $v$ is the speed of the black hole particle relative to the local gas speed and $G$ is the gravitational constant. The sound speed and density are estimated locally at the position of the black hole particle by taking an SPH average over its 64 nearest gas cell neighbours.

The Eddington accretion rate used to find the accretion rate is defined as

[TABLE]

where $m_{\rm p}$ is the proton mass, $\sigma_{\rm T}$ the cross-section for Thomson scattering, $c$ the speed of light and $\epsilon_{\rm r}$ the radiative efficiency. We set $\epsilon_{\rm r}$ to the standard value of 0.1

We seed haloes with a virial mass $M_{200}=10^{10}h^{-1}\rm M_{\odot}$ with seed black hole particles with mass $10^{5}h^{-1}\rm M_{\odot}$. Only one such particle is allowed to form per halo. Besides accretion, every black hole is allowed to grow by merging with other black holes that fall within its smoothing length, evaluated by averaging over the properties of the 64 nearest gas cell neighbours, and have a relative velocity lower or comparable to the local sound speed.

We model AGN feedback at all mass accretion regimes through an isotropic coupling of a small fraction of the AGN bolometric luminosity to the surrounding gas cells as in Di Matteo et al. (2005); Sijacki et al. (2009) or as in the ‘quasar mode’ feedback used in ILLUSTRIS (Vogelsberger et al., 2013) and ILLUSTRIS-TNG (Weinberger et al., 2017) simulations. Thermal energy is injected continuously into the 64 nearest gas cells at a rate given by

[TABLE]

where $\epsilon_{f}$ is the feedback efficiency, set to 0.05, which reproduces the normalisation of the $M_{\rm BH}-\sigma_{\star}$ relation, the correlation between the black hole mass and the stellar velocity dispersion, in fully cosmological simulations (Di Matteo et al., 2005; Sijacki et al., 2007; Di Matteo et al., 2008).

2.2 Simulation sample

We perform a total of 8 simulations, in which we explore the evolution of the targeted galaxies after their $z\,=\,6$ quasar phase in different environments. Our fiducial simulations are named Halo1 to Halo5 and these haloes are simulated till $z\,=\,3.3$. We also re-run simulations for Halo1 with a modified setup: (i) switching off black hole accretion and AGN feedback at $z\,=\,5.2$, (ii) excluding black hole growth and AGN feedback altogether and (iii) switching to fixed physical (as opposed to comoving) softening lengths at $z\,=\,6$. These simulations are named, respectively, Halo1-weakAGN, Halo1-NoAGN and Halo1-Soft. Halo1-NoAGN is also simulated till $z\,=\,3.3$ whereas Halo1-weakAGN and Halo1-Soft, simulations with finer time-steps, are only run till $z\,=\,4.5$.

3 Results

Here we present the main results from our simulations and analysis. We start with a qualitative description of the evolution of the quasar host galaxies in our simulations.

3.1 Overview

The properties of the galactic haloes targeted by our simulations at $z\,=\,6$ are summarised in Table 2. The total halo masses are all on the order of a few $10^{12}\,\rm M_{\odot}$. Their central black holes have grown to masses $5\times 10^{8}\--10^{9}\rm M_{\odot}$ as has been seen in previous simulations (Sijacki et al. 2009; Di Matteo et al. 2012; Costa et al. 2014).

By $z\,=\,3.3$ the total halo masses all grow to $\sim 10^{13}\,\rm M_{\odot}$, the central black holes masses at this redshift lie in the range $7\times 10^{9}\--9\times 10^{9}\rm M_{\odot}$ (see second part of Table 2). The stellar mass within the stellar half mass radius, as predicted by SUBFIND (Springel et al. 2001), is $\approx 2.3\times 10^{10}M_{\odot}$ at $z\,=\,6.2$ and reaches $\approx 1.5\times 10^{11}M_{\odot}$ at $z\,=\,3.3$. For comparison, in the simulation without AGN feedback, the targeted galaxy’s stellar mass grows from a comparable $\approx 2.2\times 10^{10}M_{\odot}$ at $z\,=\,6.2$ to a significantly higher value of $\approx 6.2\times 10^{11}M_{\odot}$ at $z\,=\,3.3$.

The typical cosmological environment of the targeted haloes is illustrated in Fig. 1, where we show the mass-weighted density and temperature projected along a slab of thickness $2h^{-1}\,\rm cMpc$ and $200h^{-1}\,\rm ckpc$, in the top and the second row respectively, at $z\,=\,6.2$ (left-hand panels) and at $z\,=\,3.3$ (right-hand panels). The massive galaxy resides at the intersection of a number of filaments which feed it with cold gas, as seen in previous simulations (Sijacki et al. 2009; Di Matteo et al. 2012; Dubois et al. 2012; Costa et al. 2014). Remarkably, the orientation of the filaments changes little between $z\,=\,6.2$ and $z\,=\,3.3$, though it is clear that a diffuse component becomes prominent at $z\,=\,3.3$, when the density field is overall smoother and much of the gas settles into a stable, hot atmosphere. This hot component forms through accretion shocks as well as through repeated cycles of AGN heating. The cold filaments that connect the cosmic web to the central galaxy at $z\,=\,6.2$ are more clearly disrupted at $z\,=\,3.3$. We find that such filaments become disrupted between $z\,=\,5.3$ and $z\,=\,3.9$ for all haloes including the halo without AGN feedback.

Due to its compact size, the quasar host galaxy is barely noticeable in the density projection at $z\,=\,6.2$. In contrast, at $z\,=\,3.3$, it clearly takes the shape of a larger disc, seen edge-on in the inset. The galaxy be seen as a cold, compact, disc surrounded by a large mass of cold, dense gas.

Just through visual inspection, it is clear that the targeted galaxy experiences substantial size growth, as illustrated by the stellar surface density maps in the bottom panels of Fig. 1. In the top row, we show the stellar distribution in our default setup (with AGN feedback) at $z\,=\,6.2$ on the left-hand side, and at $z\,=\,3.3$ on the right-hand side. While the galaxy appears as a flattened spheroid of $\sim\rm kpc$ radius at $z\,=\,6.2$, it has a size $\gtrsim 5\,\rm kpc$ at $z\,=\,3.3$.

The bottom row shows the results for the simulation without AGN feedback. In contrast to our default simulation, the galaxy appears to experience less significant size growth and its surface density remains high and centrally concentrated. The maximum stellar surface density for the galaxy without AGN feedback is $2.3\times 10^{6}\,\rm M_{\odot}pc^{-2}$ which is $\approx 66$ times higher than the maximum found for the galaxy with AGN feedback.

The contours in the projected stellar maps show the surface density of the stars that make up the ‘compact bulge’ at $z\,=\,6.2$. Each stellar particle has a unique ID that allows it to be traced back- and forward in time in the simulation. We select the stars inside the stellar half mass radius (1.2 kpc) at $z=6.2$ as proxy for ‘compact bulge’ stars, extract their IDs and identify their location at $z=3.3$. The contours correspond to $5,100,700\,\rm M_{\odot}\,pc^{-2}$ surface density levels, where the thickest contour corresponds to the highest surface density. At $z\,=\,6.2$, the ‘compact bulge’ stars are, by construction, concentrated within the stellar half mass radius of $\approx 1$ kpc and there is little difference between simulations with and without AGN feedback. While, at $z=3.3$ the bulge stars in the galaxy without AGN feedback are less concentrated than at $z\,=\,6.2$, they remain within a radius of $\approx 3$ kpc, a factor three larger than at $z\,=\,6.2$. More strikingly, in the galaxy with AGN feedback, the ‘compact bulge’ stars are significantly more spatially extended. The radius of the outer contour is $\approx 8$ kpc.

Thus, despite the remarkable compactness of quasar host galaxies at $z\,=\,6.2$, some process leads to the expansion of the ‘compact bulge’ and the migration of its stellar populations to the outskirts of the galaxy. We have shown that this process is much more significant in the presence of AGN feedback, which, as explained in Section 2.1, operates solely via thermal, ‘quasar mode’ feedback at all regimes. We thus now explore how quasar-driven outflows shape the evolution of the quasar host galaxies.

3.2 The impact of AGN feedback on the properties of the quasar host galaxy

Here, we investigate how strongly AGN feedback influences the evolution of quasar host galaxies and their host haloes.

3.2.1 Stellar mass vs. halo mass

We link the stellar mass in the halo to the halo mass and compare these with existing abundance matching constraints. In Fig. 2 we plot the stellar mass-halo mass ratio against halo mass. The colour-coded symbols show the evolution of the haloes with redshift and the coloured lines indicate the redshift dependent relation found by Moster et al. (2013) for $z\,=\,6.2$, $z\,=\,4.5$ and $z\,=\,3.3$. For the stellar mass and total halo mass, we use the values as predicted by SUBFIND.

We see that the simulations match the observational constraints well at $z\,=\,6.2$ and $z\,=\,4.5$. This agreement, however, breaks down at $z\,=\,3.3$ (blue coloured symbols) and it is clear that our massive galaxies form stars too efficiently, when compared to the abundance matching constraints. However, the stellar mass fractions of the simulated haloes are in much better agreement with the observed relation than in the simulation without AGN feedback. In the latter, the halo forms an order of magnitude more stars than expected based on the abundance matching constraints, while we see a disparity of only a factor two when we include AGN feedback.

Thus, up until the later stages of our simulations, galaxy stellar masses appear reasonable and, if anything, AGN feedback appears to be too weak. This is not surprising as continuous thermal energy injection is known to suffer from numerical radiative losses at the resolution typically achieved by even ‘zoom-in’ cosmological simulations (Booth & Schaye, 2009; Weinberger et al., 2018). The main point here, however, is that if AGN feedback is indeed responsible for the expansion of the galaxies’ ‘compact bulges’, this is unlikely to be because AGN feedback is unrealistically strong.

3.2.2 Black hole mass vs. stellar mass

Another relation that tests our implemented feedback model is the black hole mass - stellar mass relation. Fig. 3 shows the black hole mass as a function of total stellar mass in our various simulated haloes. The colour-coded symbols show the redshift evolution of the most massive black holes located in the nuclei of the most massive galaxies in each simulation. The black symbols give the values of 7 other massive black holes in every simulated environment and these are taken at the same redshift as the colour-coded ones.

While the slope is in good agreement with the local relation, the normalisation appears to be higher by a factor of $\approx 3$. It has been claimed that observed high redshift black holes are, typically, over-massive at fixed galaxy mass (e.g., Volonteri & Stark 2011) and that there can be significant redshift evolution in various $\rm M_{\rm BH}-$ galaxy relationships. Merloni et al. (2010) propose that $M_{\rm BH}-\rm M_{\rm stars}$ evolves with redshift as $\propto(1+z)^{0.68}$, while Decarli et al. (2010) suggest a different scaling $\propto(1+z)^{0.28}$. If we rescale the relation found by Reines & Volonteri (2015) using the scaling found by Merloni et al. (2010) (shown in orange in Fig. 3), we obtain better agreement between the simulations and the observations. We note that the normalisation of the $M_{\rm BH}-\rm M_{\star}$ relation evolves in various other state-of-the-art simulations (Barber et al., 2016; Sijacki et al., 2015).

$M_{\rm BH}$ and $M_{\rm stars}$ estimates for various observed galaxies with over-massive BHs are also shown in Fig. 3 for reference. Values were taken from Seth et al. (2014), Saglia et al. (2016), McConnell & Ma (2013), van Loon & Sansom (2015), Trakhtenbrot et al. (2015), Walsh et al. (2015) and Walsh et al. (2016) for M60-UCD1, NGC 4486B, NGC 4342, S536, CID-947, NGC 1271, and NGC 1277, respectively111Note that for NGC 4486B, NGC 4342 and NGC 1271, we plot stellar bulge mass since total stellar mass is not available.. We see that the simulated black holes compare well with the observed over-massive black holes and suggest that the simulated environments could thus be the birth place of such nearby observed over-massive black holes.

3.2.3 Black hole growth

Here, we investigate how the central black holes evolve before and after their bright quasar phase at $z\,=\,6$, In the right-hand panel of Fig. 3, we show the central black hole mass as a function of redshift, where the colour-coded symbols indicate the evolution for the different haloes.

Each black hole starts from an initial seed of $10^{5}h^{-1}\rm M_{\odot}$ placed in the massive halo at around $z\,=\,15$. By $z\,=\,10.1$, the central black holes have grown to masses $3.1\times 10^{5}\--5.6\times 10^{6}$ $\rm M_{\odot}$ through a combination of mergers and accretion (see Sijacki et al., 2009). In the redshift range $z\,=\,10-6$, we see exponential growth as the accretion rate is close to the Eddington limit almost all the time. This can be seen from the comparison with the black line which shows the expected mass for a black hole with a mass of $10^{7}\rm M_{\odot}$ growing at the Eddington rate starting at $z\,=\,10$. It is clear that the slopes of the black hole masses in all our simulations are similar to that expected for Eddington-limited accretion at $z\gtrsim 7$.

Black hole growth slows down significantly below $z\,=\,6$. This indicates that AGN feedback has become efficient. This can also be seen from the dashed lines in the right-hand panel of Fig. 3 which show ratio between the evolution of the total gravitational potential of gas within 1 kpc divided and the cumulative amount of feedback energy released with redshift. We can see that below $z\,=\,6$ the energy released by the AGN becomes comparable to the total gravitational potential in the inner kpc of the galaxy. AGN feedback has become strong enough to remove gas from the inner kpc. The feedback prevents the gas from accreting onto the black hole, thus slowing down the black hole growth.

AGN feedback is a circular process; after gas is removed from the central region, new gas will cool and move into the galactic nucleus and eventually accrete onto the black hole, starting a new feedback cycle (Costa et al. 2014). These short bursts of accretion allow the black hole to continue growing, though at modest rates. By $z\,=\,3.3$, the central black holes in our simulations reach masses in the range $6.9\--8.6\,\times 10^{9}\rm M_{\odot}$.

3.2.4 Morphology

In order to quantify the morphology of the compact stellar component that forms at the centre of the quasar host galaxies, we evaluate the fraction of kinetic energy invested in motion which is co-rotating with respect to the halo (Correa et al. 2017). We evaluate the $\kappa_{\rm rot}$ parameter, defined as

[TABLE]

where the sum is performed over all stellar particles within a spherical radius of $30\,\rm kpc$ centred at the position of the most massive halo, $m_{i}$ is the stellar particle mass, $K=\sum_{i}^{r<30\,\rm kpc}\frac{1}{2}m_{i}v_{i}^{2}$ the total kinetic energy, $L_{z,i}$ the particle angular momentum along the direction of the total angular momentum of the stellar component inside the radius of 30 kpc and $R_{i}$ the distance to the centre of the galaxy. Following Correa et al. (2017) we only take corotating stellar particles ($L_{z,i}>0$) into account which is found to be a better measure of ordered rotation. Correa et al. (2017) find that $\kappa_{\rm co}$ is a reasonable measure of visual morphology and established that $\kappa_{\rm co}=0.4$ can be used to separate disc-like galaxies ($\kappa_{\rm co}\geq 0.4$) from elliptical galaxies ($\kappa_{\rm co}<0.4$). We thus use also $\kappa_{\rm co}=0.4$ to separate the two populations.

Fig. 4 shows the redshift evolution of $\kappa_{\rm co}$ where the black line shows the mean evolution for all quasar host galaxies and the solid red line shows the evolution for the galaxy without AGN feedback. We confirm that all investigated galaxies start as dispersion-dominated ellipticals (as in Dubois et al., 2012). At $z\,=\,3.3$, the galaxies with AGN feedback have a mean value of $\kappa_{\rm co}=0.26$ and are thus in agreement with an elliptical morphology.

In addition, we find that the evolution of the mean value for all galaxies indicates that the quasar host galaxies never develop pronounced discs. From the standard deviation shown as the grey shaded region in Fig. 4, we see that disc-like stellar structures form in some quasar host galaxies. Specifically, the galaxy in Halo3 becomes disc-like at $z\,=\,7.2$ and $z\,=\,6.2$. Form $z\,=\,6.2$, this galaxy stays elliptical. The galaxy in Halo1 also becomes disc-like at $z\,=\,4.8$ but remains elliptical from $z\,=\,4.5$. Quasar host galaxies are thus able to form disc-like stellar structures but do not maintain the disc-like structure. This in contrast to the galaxy without AGN feedback, which becomes disc-like at $z\,=\,5.2$ and is able to maintain this structure reaching $\kappa_{\rm co}=0.59$ at $z\,=\,3.3$. The disc-like structure can also be seen from the stellar surface density maps, shown in the bottom panels of Fig. 1. The galaxy without AGN feedback gains this disc-like configuration at $z\,=\,5.3$ and remains disc-like from then onwards. A similar behaviour is seen in the HORIZON-AGN simulations in the case in which AGN feedback is excluded; late-time accretion of gas results in the formation of a disc while, if AGN feedback is active, the formation of the disc is prevented (Dubois et al. 2016).

AGN feedback thus appears to ensure that the galaxies hosting quasars at $z\gtrsim 6$ contain spheroidal, dispersion-dominated stellar components. We emphasize that the gas component, which evolves on Myr timescales, does have a disc-like morphology, but we find that this changes orientation rapidly as the central galaxy is bombarded by satellite galaxies as well as by a smoother inflow component. The resulting angular momentum cancellation ultimately promotes the formation of a compact, dispersion-dominated bulge (Dubois et al., 2012).

3.2.5 Circular velocity

In Fig. 5 we show stellar circular velocities profiles for all our simulations at $z\,=\,6.2$ and $z\,=\,3.3$. The circular velocity is defined as

[TABLE]

where $M_{\rm stars}(<R)$ is the stellar mass enclosed within a radius $R$.

At $z\,=\,6.2$, all quasar host galaxies are, as we have seen in Fig. 1, very compact. If we consider the radial position of the maximum (stellar) circular velocity $R_{\rm max}$ as a proxy for bulge size, we find $R_{\rm max}\approx 0.2\--0.5\,\rm kpc$ with a mean of $\sim 0.4$ kpc which is roughly twice the size of the gravitational softening length. The compact bulges have peak circular velocities between 300 km s*-1* and 500 km s*-1*. Stars make up between 40% and 57% of the total mass within the inner kpc where dark matter and gas make up between 10% and 33% of the total mass. The black hole only makes up $\sim 3$% of the total mass.

In addition, we note that, at $z\,=\,6.2$, the circular velocity profile of the simulation without AGN feedback is comparable to the circular velocity profiles found in the other five simulations. We find a peak circular velocity of $\sim 340$ km s*-1* for the halo without feedback, comparable to the mean of $\sim 350$ km s*-1* for haloes with AGN feedback. This, along with Figs 1 and 2, indicates that AGN feedback has not yet significantly impacted the structure of the galaxy by $z\,=\,6.2$, despite the presence of vigorous quasar-driven outflows (Costa et al., 2015; Curtis & Sijacki, 2016).

At $z\,=\,3.3$, all quasar host galaxies have, as we have seen for one of our simulations in Fig. 1, become less compact. The radial position of the maximum circular velocity grows to $R_{\rm max}\approx 2.5\--4.0\,\rm kpc$, an order of magnitude higher than the gravitational softening length and $R_{\rm max}$ at $z\,=\,6.2$. Thus, while a compact bulge is initially found in every of our quasar host galaxies, it expands significantly in all our simulations. The evolution of the stellar circular velocity for a single one of our systems can be seen in Fig. 6. The peak of the circular velocity profile shifts to larger radii and decreases as the redshift decreases.

Importantly, the above statements only hold in simulations with AGN feedback. In the simulation without AGN feedback, the compact bulge seen at $z\,=\,6.2$ becomes more compact and more massive by $z\,=\,3.3$; the maximum circular velocity increases from $300$ km s*-1* to $2000$ km s*-1*. The peak value of the circular velocity stays at roughly the same radius of $\sim 0.3$ kpc. From Fig. 1 we do, however, observe some stellar migration which is likely due to numerical effects. This effect is further explored in Section 4.2.

In conclusion, we thus confirm that our implemented AGN feedback model is reasonable and does a reasonable job in preventing the overproduction of stars in massive galaxies (but is still likely too weak, especially at later times). The simulated black holes also compare well with the observed over-massive black holes. When examining the morphology of the quasar host galaxies, we find that these stay elliptical throughout the simulation. The galaxy without AGN feedback, on the other hand, gets a disc-like configuration. We also find that, while all quasar host galaxies host highly compact stellar bulges at $z\,=\,6$, they all experience significant expansion with time. The galaxy without AGN feedback maintains its compact bulge. In the following section, we investigate whether AGN feedback is directly responsible for the expansion of the bulge.

3.3 Growing galaxies with AGN-driven outflows

Size is not a well-defined parameter for massive galaxies with extended stellar mass distributions, which makes it hard to compare our size measurements to other simulations and observations. The effective radius used in observations depends, for example, on resolution, filter and the adopted model for the light profile. We have therefore adopted the position of the peak of the stellar component circular velocity ($R_{\rm max}$) as a simple proxy for radius. This radius encloses approximately $15\--30$ % of the stellar mass within the virial radius. We also use the stellar half-mass radius as defined by SUBFIND $R_{\rm eff}$ as proxy for radius which, obviously, encloses $50$ % of the stellar mass within the subhalo.

3.3.1 Size growth across time

The redshift evolution of the size of the targeted galaxies, here defined as the radius of the stellar circular velocity peak is shown on the left-hand panel of Fig. 7 with solid lines for simulations including black holes and AGN feedback and a dashed line for the simulation without black holes and AGN feedback. The grey band in the figure indicates the range over which the gravitational softening lengths change over time in our simulations.

The targeted galaxies start with roughly the same size, as we have seen in the previous section, ranging from $\approx 0.2$ kpc to $\approx 0.4$ kpc at $z\,\approx\,7$. In simulations which include AGN feedback, the massive galaxies start growing in size after $z\,=\,7$. The sizes of the galaxies in Halo2 and Halo4 start to grow rapidly after $z\,=\,7$ whereas the other galaxies initially grow only slowly. The galaxy in Halo3 starts to grow faster after $z\,=\,5.7$ and the galaxy in Halo5 starts to grow rapidly in size after $z\,=\,4.8$. The galaxy in Halo1-NoAGN grows little, from $\approx 0.5$ kpc to $\approx 0.6$ kpc; it has roughly 10 % of the mean size of the other galaxies at $z\,=\,3.3$ which is $\approx 5\,\rm kpc$.

In the right-hand panel of Fig. 7, we plot the stellar half mass radius as a function of the maximum value of the circular velocity for all simulations with AGN feedback. We see that, as the quasar host galaxies grow in size, their maximum stellar circular velocity drops, indicating that the galaxies gradually become less tightly bound. The peak circular velocity and the stellar half mass radius have similar values for all haloes at $z\,=\,6.2$ and all galaxies end up having roughly the same half mass radius of $6\--9$ $\rm kpc$ at $z\,=\,3.3$. While galaxy growth is particularly rapid at $z\gtrsim 4$, it slows down at lower redshift, when the maximum circular velocity reaches values of$\approx 250\,\rm km\,s^{-1}$. Halo5 is simulated down to $z\,=\,2.6$ and continues to show size growth at a rate similar to that seen in all targeted haloes at $z<4$. One striking feature in the right-hand panel of Fig. 7 is the sudden increase of the peak circular velocity in Halo1 between $z\,=\,6.2$ and $z\,=\,5.7$. Stellar images and circular velocity profiles show that a merger is happening at this point which increases the circular velocity. Afterwards, the galaxy shows the same evolution as the other quasar hosting galaxies, the galaxy grows in size and its maximum stellar circular velocity drops gradually. We also note that for Halo1-NoAGN, the maximum circular velocity increases rapidly to extreme values $\lesssim 2000\,\rm km\,s^{-1}$, while the galaxy size remains close to $\approx 1\,\rm kpc$ (see Fig. 5).

While it is clear that the ‘puffing-up’ of the stellar component appears to be connected to the presence of an AGN, as found in various previous studies (Peirani et al., 2017; Costa et al., 2018; Choi et al., 2018), we have not yet shown that AGN feedback plays a direct role on the expansion of the compact bulges. For instance, it could act simply by reducing the gas fraction of the progenitor galaxies and thereby increasing the number of dry mergers undergone by the progenitors. As mentioned in Section 1, this is the prevailing explanation for the size growth of high redshift massive compact ellipticals. In the following sections, we show that, in our simulations, AGN feedback plays quite a direct role in galaxy size growth.

In order to understand how AGN feedback may provoke galaxy growth, we have performed a new simulation identical to Halo1, but using a larger number of simulation outputs with a $z\approx 0.02$ spacing, which is comparable to the dynamical time-scale measured at $\approx 1$ kpc. The simulation is run from $z\,=\,6.7$ to $z\,=\,4.5$, the time range where we saw the most intense size growth. In addition, we also perform another simulation like Halo1, but for which we ‘switch-off’ black hole accretion and AGN feedback at $z\,=\,5.2$. We analyse the results of these simulations in the next section.

3.3.2 Fluctuations in the gravitational potential

Fig. 8 shows the enclosed gas mass in Halo1 at different radial scales, as indicated by different colours. AGN activity becomes strong enough to regulate the central gas reservoir already at $z\,=\,6.5$, when we see a drop in the enclosed mass within a radius of $500\,\rm pc$. At $z\lesssim 6$, we see rapid, $\approx$ $1\,\rm Myr$, and strong fluctuations in the enclosed mass. The pronounced peak seen at $z\,\approx\,6$ is caused by a merger which temporarily supplies the galaxy with fresh gas. The gas mass inside $1\,\rm kpc$ then starts to decrease below $z\,=\,5.7$ more systematically, while the fluctuations in the enclosed gas mass start becoming noticeable at even larger radii of up to $3\,\rm kpc$.

The enclosed dark matter and stellar masses also undergo small fluctuations in response to the gas motions as can be seen from Fig. 8. The enclosed collisionless component mass also show a significant and gradual decline below $z\,=\,6$ which occurs mainly within scales $\lesssim 3\,\rm kpc$ and is not seen at larger scales. This suggests the gradual expansion of stellar and dark matter within those scales, in good match to what was seen in the circular velocity profiles discussed in the previous section.

We also measure the actual gravitational potential of ‘compact bulge’ stars, which we here select as those stellar particles located within a radius of 3 kpc at $z\,=\,6.2$. We track the stellar particles using their IDs and calculate their total gravitational potential $\Phi\,=\,-\sum_{i}\phi_{i}$, where the sum is performed over the gravitational potential $\phi_{\rm i}$ of each bulge stellar particle, within different radii. As can be seen in Fig. 8, we find the expected fluctuations starting from $z\,=\,5.5$ in the gravitational potential of the bulge stars within radii of 1 kpc, and from $z\,=\,5.2$ onwards, we also see the fluctuations for radii < 3 kpc. As the collisionless component migrates outwards, the gravitational potential also drops with time.

In the absence of AGN feedback (dashed lines in Fig. 8), the gravitational potential of ‘bulge stars’ increases, in contrast to the simulation with AGN feedback.

3.3.3 Ex/in situ stars

The size growth we see in our simulations is, in part, driven by the expansion of the old stellar component which already formed before $z\,=\,6.2$ (see bottom panels in Fig. 1). The accretion of stars that formed outside the galaxy (‘ex-situ’ stellar populations) could, however, also contribute to the overall size growth, as it builds up the outer layers of the quasar host galaxy.

To quantify the contribution of different populations to the size growth, we follow the stellar particles with their IDs and classify them based on their formation radius. We distinguish between ‘ex-situ’, ‘in-situ’ and ‘relic’ stars. We define ex-situ stars as those which form outside $R_{\rm thresh}\,=\,3\times$ the stellar half mass radius, as defined by SUBFIND. This radius is measured at every snapshot and changes thus with redshift. The in-situ stars are defined as those which form within $R_{\rm thresh}$. This radial threshold corresponds to $\approx$ 5% of the virial radius which is somewhat smaller than the radius used in Oser et al. (2010) and Frigo et al. (2018), which is 10% of the virial radius. The choice of radial threshold, however, does not significantly affect our results as long as it lies in the range 5-15% of the virial radius. Stars at $z\,=\,4.5$ within $R_{\rm thresh}$ that already existed at $z\,=\,6.7$ within $R_{\rm thresh}$ are defined as relic stars. We select an initial redshift of $z\,=\,6.7$ here, because a highly tightly bound stellar bulge has already formed by that time.

The density profile evaluated at $z\,=\,4.5$ shown in Fig. 9 with the solid black line consists of three different components: in-situ stars, ex-situ stars and relic stars. We find that the in-situ stellar population dominates the stellar mass at scales $\lesssim 10\,\rm kpc$, making up $\approx 55$% of the final mass, and, in particular, within the innermost few kpc, where size growth takes place. In contrast, ex-situ stars make up $\approx 23$% of the total stellar mass within $R_{\rm thresh}$, while the relic stars make up $\approx 18$% of the final mass within $R_{\rm thresh}$. $\approx 4$% of the stellar mass is not found using any these definitions, because they lie outside $R_{\rm thresh}$ at $z\,=\,4.5$. The ex-situ stellar population is more concentrated than the in-situ stellar population. We find that early accreted ex-situ stars make up this concentrated central part. The contribution of ex-situ stars to the total final stellar mass is thus small compared to the contribution of in-situ stars, but comparable to the contribution of relic stars. Note that if we select relic stars to be those which are in place at $z\,=\,5.5$, instead of $z\,=\,6.7$, we find that they contribute $\approx 60$% to the total stellar mass within $R_{\rm thresh}$.

To see how different stellar populations contribute to the size growth, we check the difference between the formation radius and their radius at $z\approx 4.5$ ($\Delta R$ = formation radius - final radius). We see that $\approx 39$% of the in-situ stars migrate inwards and have a greater formation radius than their final radius. In contrast, $\approx 61$% of the in-situ formed stars migrate outwards where $\approx 1$% of the outward migrated stars have a significantly bigger final radius ($\Delta R$ < -5 kpc).

The relic stars also mostly move out to greater radii with time as can be seen in Fig. 9 from the difference between the black dashed line, the profile at$z\,=\,6.7$, and the green solid line, the profile at $z\,=\,4.5$. About $86$% of the relic stars move to larger radii and $\approx 5$% of the outward-migrating stars have a significantly bigger final radius ($\Delta R$ < -5 kpc).

In order to quantify how relic stars are affected by size growth, we use a similar approach as Choi et al. (2018). First, we select all stars within the stellar half mass radius at $z\,=\,6.7$. We then trace these stars down to $z\,=\,3.3$ and measure the radial distance $D$ they have traveled. The red and black histograms in the right-hand panel of Fig. 9 show the stellar radial migration distance from $z\,=\,6.7$ to $z\,=\,3.3$ for, respectively, the galaxy without and with AGN feedback in Halo1.

Most stars migrate radially outwards since $z\,=\,6.7$ for the galaxy with AGN feedback, 97.1% of the selected stars migrate more than 1 kpc outwards. The difference with the galaxy without AGN feedback is clear, only 17.4% of the stars migrate radially outward more than 1 kpc. 82.6% of star particles do not migrate further than 1 kpc since $z\,=\,6.7$. The majority of star particles which constituted the core at $z\,=\,6.7$ still are within the central region in the galaxy without AGN feedback.

In this section we investigated three different stellar populations and showed how they contributed to the size growth observed in Section 3.2.5. We find that the contribution of in-situ stars to the density profile at $z\,=\,4.5$ dominates. In-situ and relic stars are found to contribute most to the size growth as they migrate mostly outwards. We also show that relic stars migrate outwards in the galaxy in Halo1 which is not the case for the galaxy in the re-run simulation of Halo1 without AGN feedback. In addition, we observe that the density profile of the galaxy shown on the left-hand panel in Fig. 9 experiences a flattening because of the outward-migrating stars. The black dash-dotted line (density profile at $z\,=\,6.7$) shows an overall steep profile whereas the black solid line (density profile at $z\,=\,4.5$) shows the start of the formation of a core in the profile. For the galaxy without AGN feedback we observe that the dark matter and the stellar density profiles only steepen with decreasing redshift. This will be further discussed in Section 4.3.

4 Discussion

4.1 Growing massive galaxies through AGN feedback

We have unambiguously shown that AGN feedback is directly responsible for the gradual expansion of the compact stellar bulges that form in massive galaxies at $z>3$. Size growth driven by AGN feedback in our simulations begins roughly when the energy injected by the AGN becomes comparable to the binding energy of the compact bulges.

While they have a scale radius of a few $100\,\rm pc$ at $z\approx 6$, a time at which black hole growth is efficient and proceeds close to the Eddington rate, the compact bulges gain sizes an order of magnitude larger by $z\,=\,3$ as long as AGN feedback operates. The ‘puffing-up’ of the stellar component has been previously described in the literature (Fan et al., 2008; Ragone-Figueroa & Granato, 2011; Lapi et al., 2018) and has also been seen in cosmological simulations (Peirani et al., 2017, 2019; Costa et al., 2018; Choi et al., 2018) comparing the properties of galaxies with and without AGN feedback. We note that in those studies, fluctuations in the gravitational potential driven by outflows had been suggested as the main agent of galaxy growth, but this had not been explicitly shown. That AGN feedback drives sufficiently strong gravitational potential fluctuations has, however, been shown using idealised simulations (Martizzi et al., 2013). The above described mechanism will be further explored in Section 4.2.

Based on our results, we find that AGN feedback may result in significant size growth in galaxies that have recently hosted intense quasar activity. We suggest that massive galaxy growth occurs in different phases: (i) At high-z, these systems initially go through a phase in which gas is accreted onto the nucleus at high rates, powering black hole accretion but also the formation of a compact $\approx 500\,\rm pc$ stellar bulge. At this point, the velocity dispersions of the stellar bulges are high $\approx 500\,\rm km\,s^{-1}$, (ii) as enough AGN energy is coupled to the gas in the nucleus and the mass is expelled in brief AGN-driven outflows, the stellar bulges start to unbind and expand until they reach radii $\approx 5\,\rm kpc$, (iii) as the gas fractions drop, dry mergers start contributing significantly to size-growth. Note that another mechanism which can aid in the formation of stellar cores at the centre of massive galaxies is ‘scouring’ by binary black holes (Begelman et al., 1980), which we can, however, not follow with our simulations due to insufficient resolution. To some extent, the size growth of massive galaxies may be the result of multiple, independent physical processes operating at different stages in their evolution.

Another interesting implication of our findings relates to the possible upturn in the $M_{\rm BH}\--\sigma_{\star}$ relation seen for $\sigma_{\star}\gtrsim 270\,\rm km\,s^{-1}$ (e.g. Kormendy & Ho, 2013). In our simulations, we find a drop in the peak circular velocity from $\lesssim 500\,\rm km\,s^{-1}$ to $\lesssim 300\,\rm km\,s^{-1}$ by $z\approx 3$. Thus, while the black hole mass doesn’t change significantly, the velocity dispersion drops considerably, resulting in an upturn in the $M_{\rm BH}\--\sigma_{\star}$ relation.

The observed $M_{\rm BH}\--\sigma_{\star}$ saturates at $\sigma_{\star}\gtrsim 270\,\rm km\,s^{-1}$ (e.g. Kormendy & Ho, 2013), when the morphology of elliptical galaxies switches from predominantly coreless to mainly cored. This finding can be explained by our theoretical results; when black holes grow to sufficiently high masses, they release energies comparable to the binding energy of their host galaxies. As the stellar component expands in response to AGN-driven gas ejection, the massive galaxy develops a cored stellar profile and the stellar velocity dispersion in the central regions drops. AGN feedback is here responsible for both the formation of a core and the drop in the velocity dispersion necessary to result in an upturn in the $M_{\rm BH}\--\sigma_{\star}$. One prediction of our scenario is that the characteristic stellar velocity dispersion at which the $M_{\rm BH}\--\sigma_{\star}$ saturates should evolve with redshift. At $z\gtrsim 5$, the amount of AGN energy released is not sufficient and there has not been enough time to initiate strong galaxy size growth. Thus, at this point, we predict that the $M_{\rm BH}\--\sigma_{\star}$ should persist even to velocity dispersions $\sigma_{\star}\gtrsim 400\,\rm km\,s^{-1}$. Since those systems are likely to be the first in which supermassive black holes form and grow most rapidly, they should then be the first to be unbound by AGN feedback. The population of systems with $\sigma_{\star}\gtrsim 400\,\rm km\,s^{-1}$ should then start thinning out and an upturn in the $M_{\rm BH}\--\sigma_{\star}$ relation should become visible. Testing the importance of this scenario in simulations following large cosmological boxes will be important to develop more rigorous observational diagnostics.

We have shown that our adopted AGN feedback model is, if anything, too weak. It is thus unlikely that the size growth we find for the quasar host galaxies is driven by too strong AGN feedback. There are two possible ways, however, in which our AGN feedback model may amplify galaxy size growth. One is related to limited spatial resolution and the absence of a cold ISM phase from our simulations. As we increase resolution, we may expect the gas distribution to become clumpier, which may possibly permit AGN outflows to escape without coupling with such large volume (e.g. Bourne et al., 2015; Curtis & Sijacki, 2016). Determining whether this is the case will, however, require an accurate treatment of the cold ISM in massive, high redshift galaxies in future simulations. Another possible way in which galaxy size growth may become less important than seen in our simulations is if AGN feedback couples to halo gas directly rather than to gas in the galactic nucleus. For instance, if AGN feedback proceeds through jets which pierce through the ISM, but inflate hot bubbles in the diffuse halo, then gas in the galactic nucleus is likely to be less efficiently expelled. The heating of halo gas may also suppress cooling onto the central galaxy, which in turn may suppress black hole accretion. Thus, it will be important to establish whether the same AGN-driven galaxy growth is seen in simulations employing different AGN feedback models.

Another uncertainty in this study is the adopted Bondi-Hoyle-Lyttleton formalism to model the black hole accretion. Due to insufficient resolution, we are unable to resolve the accretion onto the black hole self-consistently and use therefore a subgrid model. The parameters of this model (gas density and sound speed) are estimated from scales that are resolved in the simulation. It has been shown that different choices on how these parameters are measured could lead to black hole masses that differ in two orders of magnitude (Curtis & Sijacki, 2015; Negri & Volonteri, 2017). In addition, the Bondi-Hoyle-Lyttleton formalism contains assumptions which are unlikely to be true for realistic black hole accretion. The model assumes, for example, spherical symmetry and does not consider the accretion of a turbulent, non-uniform gas flow (Negri & Volonteri, 2017). Alternative models have been developed where, for instance, the angular momentum of the acrreting gas is taken into account (Tremmel et al., 2017) or where the accretion rate is based on the viscous evolution of the accretion disk (Debuhr et al., 2011, 2012). Improving the assumed model is highly desirable to capture AGN feedback more realistically and test whether AGN-driven galaxy growth will still occur.

4.2 Other mechanisms of size growth

Other potential origins for massive galaxy size growth include physical effects such as major mergers, minor mergers, and adiabatic expansion due the sudden gas expulsion, but there are also potential numerical effects stemming from the non-constant physical size of the gravitational softening length.

In the case of dry major mergers, the effective radii and stellar mass both approximately double in the merger remnant ($R_{eff}\propto M_{\star}$). According to this scenario, the high-redshift systems would be uniformly denser than their low-redshift descendants (Hernquist et al. 1993; Boylan-Kolchin et al. 2006; Naab et al. 2009). This kind of growth is not observed for the galaxies in our simulation as can be seen from the right-hand panel in Fig. 7 showing the stellar half mass radius as a function of maximum circular velocity. We see a gradual size growth, that is not expected from a major merger scenario. In the minor merger scenario, the compact core remains intact and the overall size growth is then mainly driven by the acquisition of stellar material from less-bound galaxies. Dubois et al. (2016), for example, find that the accretion of ex-situ stars drive the galaxy growth at redshift $<4$. In our simulations, we, however, find a small fraction of ex-situ stars at least until $z\,=\,4.5$, when the galaxy has already grown in size significantly, and, in addition, show unambiguously that the compact bulge expands. Thus, we conclude that the increase in scale radius we find in our simulations is not driven by dry minor mergers.

Fan et al. (2008) and Fan et al. (2010) proposed a mechanism in which AGN feedback directly triggers size growth. In this scenario, which is akin to that which may operate via supernova explosions in dwarf galaxies (Navarro et al., 1996), AGN feedback removes large masses of gas from the central regions on the dynamical timescale or shorter. The centripetal force, holding orbiting collisionless particles in their orbit, drops substantially and the system relaxes into a new equilibrium configuration with a ‘puffed-up’ stellar and dark matter distributions. Repeating the process accentuates the effect, which allows a significant transformation to be accomplished by recycling a small amount of gas instead of expelling an unfeasibly large amount of gas in one episode. As we have seen in Section 3.3.2, AGN feedback in the central galaxy in Halo1 removes large amounts of cold gas from the central regions on $\rm Myr$ timescales. The final mean size of $\approx 5\,\rm kpc$ is in line with the expected final radius of $3-5$ kpc after ‘puffing up’ the stellar distribution by feedback processes as described by Lapi et al. (2018).

We explicitly investigate the influence of AGN feedback on the size growth by running the same simulation as described in Section 3.3, but now we switch AGN feedback off at $z\,=\,5.2$. In Fig. 10, we compare the stellar half mass radius as predicted by SUBFIND (see also Section 3.3.1) for this simulation with that for the simulation with AGN feedback. We see that the characteristic size of the targeted galaxies in these two simulations become very different at $z<5.2$, the galaxy where AGN feedback is switched off stops growing in size and gradually decreases in size. The galaxy with AGN feedback has 2.2 times the size of the galaxy with AGN feedback switched off at $z\,=\,4.6$.

Fig. 11 shows the enclosed gas mass (left-hand panel) and enclosed stellar mass (right-hand panel) at different radial scales, as indicated by the different colours, for our fiducial simulation (Halo1) and the simulation where AGN feedback is switched off at $z\,=\,5.2$ (Halo1 - weakAGN). We see that the fluctuations in the enclosed gas mass and enclosed stellar mass decrease and almost disappear after AGN feedback is switched off at $z\,=\,5.2$ for all radii. The fiducial simulation continues to show the fluctuations. In addition, we see that the enclosed stellar mass in the simulation where AGN feedback is switched off gradually increases whereas the enclosed stellar mass gradual decreases in the fiducal simulation (see also Fig. 8). This comparison shows that even the fluctuations at a radial scale of 300 pc, close to the resolution of the simulation, are physically meaningful and driven by AGN driven outflows. We also see that the gradual expansion of stellar matter stops when AGN feedback is switched off.

We also check how switching off AGN feedback affects the migration of relic stars. We take the same strategy as described in Section 3.3.3 and select all stars within the stellar half mass radius at $z\,=\,5.2$. We then trace these stars down to $z\,=\,4.6$ and measure the radial distance $D$ they have traveled. The histograms in the right-hand panel of Fig. 10 show the stellar radial migration distance for the galaxy with AGN feedback switched off at $z\,=\,5.2$ (red) and the galaxy with AGN feedback (black). As we have also seen in Fig. 9, most stars migrate radially outwards in the galaxy with AGN feedback. Since $z\,=\,5.2$, 70.3% of the selected stars migrate outwards. The difference with the galaxy with AGN feedback switched off is clear, most stars (65.2%) migrate inwards and if stars migrate outwards they do not migrate as far as is seen for the galaxy with AGN feedback. AGN feedback thus clearly promotes the outward migration of stars in the galaxy.

Using the ILLUSTRIS simulations, Wellons et al. (2016) discuss how numerical effects may lead to adiabatic expansion of massive, compact galaxies. If the softening length is fixed in comoving coordinates, the physical scale over which the gravitational potential is softened increases with time, and stellar particles may migrate outwards as a result. Wellons et al. (2016) find that such numerical relaxation can change the stellar half mass radius by about 15% from $z\,=\,2$ to $z\,=\,0$, much smaller than the size growth we find in our simulations.

We investigate the influence of numerical relaxation by running the same simulation as described in Section 3.3, but now keeping a constant physical softening length at $z<6$. In Fig. 10, we compare the stellar half mass radius as predicted by SUBFIND (see also Section 3.3.1) for this simulation with that for the simulation with constant comoving softening length. We see that the characteristic size of the targeted galaxies in these two simulations only becomes systematically different at $z\,<\,5.5$. The simulation with constant physical softening length results in a galaxy with a greater size as the size found in the simulation with comoving physical softening length. We find that black holes grow more by accretion in the simulation with fixed physical softening length and therefore the injected feedback energy is somewhat larger. This promotes slightly more efficient galaxy size growth. It is clear that there is still significant galaxy growth in the simulation with constant physical softening length. The change in physical softening length that occurs in our fiducial simulations thus contributes only weakly to the galaxy growth we find. In addition, we see in Fig. 7 that there is almost no size growth in the simulation without AGN feedback, for which numerical relaxation would be the main reason for size growth. The size growth in the galaxy without AGN feedback is about 13%. We therefore conclude that the softening length cannot explain the observed expansion.

AGN-driven expansion induced by AGN feedback thus seems the main explanation for the observed size growth in our galaxy sample, while other size growth mechanisms may contribute or even dominate at later stages in the evolution of the massive galaxies. Also, it is important to stress that we here analyse massive haloes at $z\,=\,6$, which form in extreme environments which favour rapid black hole growth. The size growth seen might be connected to this rare environment.

Another recent study by Genel et al. (2018) suggested a similar mechanism for size growth as we found. Genel et al. (2018) analysed scaling relations and evolution histories of galaxy sizes in Illustris-TNG. The most massive galaxies investigated in Genel et al. (2018) quench at a stellar mass of $\sim 10^{11}\rm M_{\odot}$. They subsequently grow in size by $\approx 1.5$ dex down to $z\,=\,0$ as well as in mass, in agreement with the $r\propto M^{2}$ relation expected from minor dry mergers.

Genel et al. (2018), however, find that quenched galaxies that are less massive ($\sim 10^{10.5}\rm M_{\odot}$ at $z\,=\,0$) follow a size-mass relation which is significantly steeper than $r\propto M^{2}$, which may be partly explained by AGN feedback-induced size growth.

4.3 Cusp/core transformation

The flattening in density profiles of both the stellar and dark matter component can be seen in Fig. 8 as a decrease in enclosed mass at small radii (< 1 kpc). Fig. 12 shows the actual density profiles of stars and dark matter for $z\,=\,6.2$ and $z\,=\,3.3$ and shows that these profiles flatten significantly over this redshift range. The flattened density profiles are more compatible with the observations of both the stellar and dark matter profiles. Massive elliptical galaxies exhibit very shallow slopes in the stellar surface brightness profiles within small radii of $\approx$ 1 kpc (see Quillen et al. 2000; Laine et al. 2003; Graham 2004; Trujillo et al. 2004; Lauer et al. 2005; Côté et al. 2007; Kormendy et al. 2009; Graham 2013).

Pontzen & Governato (2012) simulated two dwarf galaxies with a mass of $\sim 10^{10}\rm M_{\odot}$ and investigated their evolution under the influence of supernova feedback. The mass scale of these simulations is different to our scale mass, as our simulations investigate the high mass end of the galaxy mass function. They also found that the potential in the central kiloparsec changes on sub-dynamical time-scales over the redshift interval $4>z>2$. These changes were, in contrast to our simulations, induced by feedback from central bursts of star formation. They also show that the fluctuations irreversibly transfer energy into collisionless particles, thus generating the dark matter core, as had been shown much earlier by Navarro et al. (1996). Observations that suggest some dwarf galaxies contain AGN in their nuclei (e.g. Pardo et al., 2016; Penny et al., 2018). It will be interesting to explore whether AGN could shape the stellar and dark matter profiles also in dwarf galaxies in the future.

Martizzi et al. (2013) also find that AGN feedback can induce a cusp-core transformation, which they showed by using idealised, controlled simulations to follow the response of a massive, cluster-like dark matter halo to feedback from a central black hole of $M_{\rm BH}=10^{9}\rm M_{\odot}$ over 3.5 Gyr. They fitted a power law to the density profiles, in the region $0.4<R<8$ kpc, showing a similar flattening of the slope of the dark matter profile.

Peirani et al. (2017, 2019) used the HORIZON and HORIZON-AGN simulations to investigate how AGN feedback affects the evolution of the inner density profiles of massive dark matter haloes ($5\times 10^{11}\rm M_{\odot}$ - $1\times 10^{13}\rm M_{\odot}$) and the associated galaxies. They find that the density profiles significantly evolve with time. At $z\approx 1.5$, the mean halo density profiles flatten in the presence of AGN activity. At $z=0$, the halo (total) density profiles drift back to a more cusped shape, as the AGN feedback efficiency drops. It will be interesting to see if the dark matter profiles in high-z quasar host galaxies experience a similar evolution at lower redshift ($z\leq 2$).

5 Summary & Conclusions

We have performed a suite of cosmological, zoom-in hydrodynamical simulations using the moving mesh code AREPO in order to investigate how the stellar component of the galaxies hosting $z\,=\,6$ quasars evolve after their initial period of rapid black hole growth. Our main findings are:

• •

Stellar compact bulges with peak stellar circular velocities in the range $300\--500\,\rm km\,s^{-1}$ and scale radii of $\approx 300\,\rm pc$ form by $z\,=\,6$ even in the presence of strong AGN feedback. While strong outflows are generated (see e.g. Costa et al., 2015), their impact on the structure of the quasar host galaxy at $z\,=\ 6$ is marginal.

• •

The quasar host galaxies remain velocity dispersion dominated ellipticals in the presence of AGN feedback, while they become disc-like in its absence.

• •

Despite their high binding energy, all compact bulges in our quasar systems experience considerable size growth from $\approx 300\,\rm pc$ at $z\,\approx\,6$ to $\approx 5\,\rm kpc$ at $z\,\approx\,3$.

• •

Size growth is accompanied by a drop in the circular velocity, which falls from high values in the range $300\--500\,\rm km\,s^{-1}$ to $200\--300\,\rm km\,s^{-1}$. This expansion of the central stellar bulge is seen only if AGN feedback is included.

• •

The inner stellar and dark matter density profiles are flattened in the presence of AGN feedback. In the absence of AGN feedback, the compact bulge does not grow significantly in size and the inner dark matter and stellar density profiles instead become steeper.

• •

We show that bursty AGN feedback drives fluctuations in the gravitational potential in the innermost few kpc of the quasar host galaxies. We show that these fluctuations lead to the gradual expansion of the collisionless (stellar and dark matter) components.

• •

We explicitly verify that switching-off AGN feedback half-way through the simulation leads to an interruption in galaxy size growth and the fluctuations driven by AGN feedback, showing that AGN feedback is directly responsible for the generation of a stellar core.

We thus conclude that AGN-driven outflows may be a key ingredient in generating stellar cores at the centre of massive, elliptical galaxies. In fact, the generation of cored stellar profiles is a test to the ability of AGN-driven outflows to remove large amounts of the ISM on short time-scales. This mechanism is also a potential alternative (or may operate in addition) to ‘scouring’ by binary black holes (e.g. Rantala et al., 2018). AGN-driven outflows may also be responsible for kick-starting the size growth of massive galaxies in the high-redshift Universe, before size growth proceeds through dry mergers. AGN may thus not only suppress the ability of massive galaxies to form stars, but also shape the gravitational potential that facilitates black hole accretion and star formation to begin with.

A crucial next step will be to investigate the potential of AGN feedback to drive galaxy size growth in cosmological boxes with a large number of galaxies. In addition, higher resolution and additional physical processes (AGN and stellar radiation) as well as an exploration of the effect of different feedback mechanisms will prove essential. It is likely that different AGN feedback models will affect galaxy size growth differently, depending on where the energy is effectively injected (halo gas vs. ISM) and on how much gas AGN outflows are able to eject. This means that observational constraints of the inner structure of massive galaxies in particular at $z\gtrsim 2$ are direct tests to the ability of AGN-driven outflows to expel the central ISM and may ultimately be used to narrow down on the coupling efficiency of AGN feedback.

Acknowledgements

We thank the anonymous referee for a constructive and helpful report. We thank Benny Trakhtenbrot, Martin Haehnelt, Thorsten Naab, Joop Schaye, Debora Sijacki, Volker Springel and Simon White for helpful discussions and comments on the manuscript. This work was partially carried out on the Dutch national e-infrastructure with the support of SURF cooperative. We acknowledge that the results of this research have been achieved using the DECI resource Eagle, at PSNC in Poland with support from the PRACE aisbl.