How AGN feedback drives the size growth of the ﬁrst quasars

Quasars at z = 6 are powered by accretion onto supermassive black holes with masses M BH ∼ 10 9 M (cid:12) . Their rapid assembly requires eﬃcient gas inﬂow 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 ≈ 0 . 5 kpc . Despite their high binding energy, we ﬁnd that these compact bulges expand at z < 6 , with their half-mass radii reaching ≈ 5 kpc by z = 3 . The circular velocity drops by factors ≈ 2 from their initial values to 200 − 300 km s − 1 at z ≈ 3 and the stellar proﬁle undergoes a cusp-core transformation. By tracking individual stellar populations, we ﬁnd that the gradual expansion of the stellar component is mainly driven by ﬂuctuations in the gravitational potential induced by bursty AGN feedback. We also ﬁnd that galaxy size growth and the development of a cored stellar proﬁle does not occur if AGN feedback is ineﬀective. Our ﬁndings suggest that AGN-driven outﬂows 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 BH − σ (cid:63) seen at high velocity dispersions σ (cid:63) .


INTRODUCTION
Detections of bright quasars at z > 6 show that supermassive black holes with masses M BH 10 9 M 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 BH 10 10 M (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 activegalactic-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 re-E-mail: dvdvlugt@strw.leidenuniv.nl quired 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 DM 10 12 M (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 BH 10 9 M 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 σ * . In their high resolution cosmological simulations, Dubois et al. (2012) find galaxies with stellar bulges with σ * > 600 km s −1 measured within ≈ 500 pc embedded within haloes with virial mass M vir 10 12 M 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 6. Costa et al. (2015) show, using similar simulations, that 'compact bulges' with peak circular velocities 500 km s −1 form even in the presence of supernova and AGN feedback. Using ∼ 10 pc resolution simulations, Curtis & Sijacki (2016) find peak circular velocities as high as ≈ 900 km s −1 for the 'compact bulge' at z ≈ 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-σ 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 vir 3 × 10 12 M 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.

SIMULATIONS
We perform cosmological, 'zoom-in', hydrodynamic simulations targeting the five most massive haloes found in the dark matter-only Millennium simulation ) at z = 6.2 and follow their evolution down to z = 3. These haloes represent high-σ 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 ΛCDM cosmology with parameters h = 0.73, Ω m = 0.25, Ω Λ = 0.75, Ω b = 0.041, identical to the parent Millennium simulation, but a lower σ 8 value of 0.8, in better accordance with the Planck cosmology (Planck Collaboration et al. 2014). 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.

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 ) 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 gas = 1.8 × 10 6 h −1 M (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 ≈ 1.5% of the cells within the virial radius of Halo1 can be classified as 'small cells' (≤ 2 × ∆x, see Table 1). Individual high-resolution dark matter particles have a mass m dm = 6.8 × 10 6 h −1 M and the stellar particles have m = 1.3 × 10 5 h −1 M . We use a fixed comoving gravitational softening soft = 1h −1 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 timedependent UV background with reionization occurring at z ≈ 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 cm −3 and with temperatures T < 10 5 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 M = η M , where M is the star formation rate and η 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 η = 1 and kick-velocity of 1200 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 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 ∼ 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 where α 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, ρ gas and c 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 where m p is the proton mass, σ T the cross-section for Thomson scattering, c the speed of light and r the radiative efficiency. We set r to the standard value of 0.1 We seed haloes with a virial mass M 200 = 10 10 h −1 M with seed black hole particles with mass 10 5 h −1 M . 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 where f is the feedback efficiency, set to 0.05, which reproduces the normalisation of the M BH − σ 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).

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.

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.

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 M . Their central black holes have grown to masses 5×10 8 −10 9 M 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 ∼ 10 13 M , the central black holes masses at this redshift lie in the range 7 × 10 9 − 9 × 10 9 M (see second part of Table 2). The stellar mass within the stellar half mass radius, as predicted by SUBFIND (Springel et al. 2001), is ≈ 2.3×10 10 M at z = 6.2 and reaches ≈ 1.5 × 10 11 M at z = 3.3. For comparison, in the simulation without AGN feedback, the targeted galaxy's stellar mass grows from a comparable ≈ 2.2 × 10 10 M at z = 6.2 to a significantly higher value of ≈ 6.2 × 10 11 M at z = 3.3.
The typical cosmological environment of the targeted   Figure 1. Illustration of the typical cosmological environment of the targeted haloes with AGN feedback in the simulation. The first row shows mass-weighted density and temperature projected along a slab of thickness 2h −1 cMpc. The second row is a zoom-in of the first row projected along a slab of thickness 200h −1 ckpc. The left-hand panels show the environment at z = 6.2 and at the right-hand panels show the environment z = 3.3. The stellar maps of the central galaxy in the simulation with and without AGN feedback at the bottom of the figure are projected along a slab of thickness 100h −1 ckpc, sizes are in physical coordinates. The contours represent the surface density of bulge stars (see text for definition) and correspond to 5, 100, 700 M pc −2 surface density levels, where the thickest contour corresponds to the highest surface density. Table 2. Main properties of the sample of the 6 haloes followed by our simulations, evaluated at z = 6.2 and 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). haloes is illustrated in Fig. 1, where we show the massweighted density and temperature projected along a slab of thickness 2h −1 cMpc and 200h −1 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 ∼ kpc radius at z = 6.2, it has a size 5 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 × 10 6 M pc −2 which is ≈ 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 M 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 ≈ 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 ≈ 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 ≈ 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.

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.    Moster et al. (2013). There is good agreement between the simulated galaxies and the abundance matching constraints, though the simulated galaxies lie somewhat above the Moster et al. (2013) relation at z = 3.3. If AGN feedback is excluded (open symbols), then the discrepancy between simulations and abundance matching constraints widens significantly. Note that at z = 6.2 the symbols of Halo1 and Halo1-NoAGN overlap.

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.

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 ≈ 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 M BH − galaxy relationships. Merloni et al. (2010) propose that M BH − M stars evolves with redshift as ∝ (1 + z) 0.68 , while Decarli et al. (2010) suggest a different scaling ∝ (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 BH − M relation evolves in various other state-of-the-art simulations (Barber et al. 2016;Sijacki et al. 2015).
M BH and M stars estimates for various observed galaxies with over-massive BHs are also shown in Fig CID-947, NGC 1271, and NGC 1277, respectively 1 . 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.

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 M placed in the massive halo at around z = 15. By z = 10.1, the central black holes have grown to masses 3.1 × 10 5 − 5.6 × 10 6 M 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 M growing at the   Figure 3. The left-hand panel shows the black hole mass versus total stellar mass for all haloes compared to the relation for elliptical galaxies of Reines & Volonteri (2015). The blue-shaded region indicates the intrinsic scatter around the observed relation. We also re-scale this relation to higher redshift with the (1 + z) 0.68 scaling suggested by Merloni et al. (2010). The orange line indicates the re-scaled relation at z = 2.8. Various observed M BH − M stars outlier galaxies are indicated with black stars (see text for details and references). The black open symbols refer to other black holes selected from every simulation at every snapshot. They give the values of 7 other massive black holes in every simulated environment taken at the same redshifts as the colour-coded ones. The right-hand panel shows the evolution mass of the most massive black holes in each simulation (left axis). The black line shows the black hole mass expected from constant Eddington-limited accretion, starting from an initial mass of 10 7 M at z = 10. The dashed lines show the ratio between the total gravitational potential, measured within 1 kpc, and the total feedback energy released (right axis).
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 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 × 10 9 M .

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 κ rot parameter, defined as where the sum is performed over all stellar particles within a spherical radius of 30 kpc centred at the position of the most massive halo, m i is the stellar particle mass, i 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 κ co is a reasonable measure of visual morphology and established that κ co = 0.4 can be used to separate disc-like galaxies (κ co ≥ 0.4) from elliptical galaxies (κ co < 0.4). We thus use also κ co = 0.4 to separate the two populations. Fig. 4 shows the redshift evolution of κ 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 κ 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  . The solid black line shows the mean of all galaxies with AGN feedback, whereas the grey shaded region shows the standard deviation. The red solid line shows the evolution of the galaxy without AGN feedback. All quasar host galaxies are consistent with an elliptical morphology throughout the simulation. In the absence of AGN feedback, however, the massive galaxy develops a prominent stellar disc. From the standard deviation, we see that disc-like stellar structures form in some quasar host galaxies from time to time. These galaxies are, however, unable to maintain this structure. 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 κ 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 disclike 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 6 contain spheroidal, dispersiondominated 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).

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 where M 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 max as a proxy for bulge size, we find R max ≈ 0.2 − 0.5 kpc with a mean of ∼ 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 ∼ 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 ∼ 340 km s −1 for the halo without feedback, comparable to the mean of ∼ 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 max ≈ 2.5 − 4.0 kpc, an order of magnitude higher than the gravitational softening length and R 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 ∼ 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 feed-  Figure 6. The circular velocity profile of the stellar component for Halo1. The profile is shown at different redshifts with the coloured solid lines. The vertical grey area marks the location of the gravitational softening length for z = 5.2 to z = 3.8. As redshift decreases, the peak circular velocity drops and its location shifts to larger radii, indicating a gradual expansion of the quasar host galaxy.
back maintains its compact bulge. In the following section, we investigate whether AGN feedback is directly responsible for the expansion of the bulge.

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 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 eff as proxy for radius which, obviously, encloses 50 % of the stellar mass within the subhalo.

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 ≈ 0.2 kpc to ≈ 0.4 kpc at z ≈ 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  Figure 7. The left-hand panel shows the location of the maximum stellar circular velocity as a function of redshift for the six investigated galaxies. The solid lines show the evolution in our simulations with AGN feedback, while the dashed line shows the results for the simulation without black holes and AGN feedback. The grey band indicates the range of stellar gravitational softening lengths, which, since they are fixed in comoving coordinates, change with redshift in physical coordinate. At z ≈ 4, the massive galaxy in Halo4 undergoes a major merger, which results in a double-peaked circular velocity profile; we have here taken the radius associated with the second peak.
The right-hand panel shows the redshift evolution of the stellar half mass radius as a function of the maximum stellar circular velocity. We see that for z 6, all quasar host galaxies experience size growth, while their maximum circular velocity drops.
grows little, from ≈ 0.5 kpc to ≈ 0.6 kpc; it has roughly 10 % of the mean size of the other galaxies at z = 3.3 which is ≈ 5 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 kpc at z = 3.3. While galaxy growth is particularly rapid at z 4, it slows down at lower redshift, when the maximum circular velocity reaches values of≈ 250 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 2000 km s −1 , while the galaxy size remains close to ≈ 1 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  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 ≈ 0.02 spacing, which is comparable to the dynamical time-scale measured at ≈ 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. Fig. 8(a) 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 pc. At z 6, we see rapid, ≈ 1 Myr, and strong fluctuations in the enclosed mass. The pronounced peak seen at z ≈ 6 is caused by a merger which temporarily supplies the galaxy with fresh gas. The gas mass inside 1 kpc then starts to decrease below z = 5.7 more systemat- Clearly, the enclosed gas mass becomes highly variable after z ≈ 5.7 within ∼ kpc scales, approximately when the released AGN feedback energy becomes comparable to the galaxy binding energy. From this time onward, the central stellar and dark matter mass also starts to drop. The redshift evolution of the total gravitational potential for bulge stars (selected at z = 6.2) within different radii is shown in (c). In the absence of AGN feedback (dashed lines), the gravitational potential of 'bulge stars' increases, in contrast to all other simulations.

Fluctuations in the gravitational potential
ically, while the fluctuations in the enclosed gas mass start becoming noticeable at even larger radii of up to 3 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(b). The enclosed collisionless component mass also show a significant and gradual decline below z = 6 which occurs mainly within scales 3 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 Φ = − i φ i , where the sum is performed over the gravitational potential φ i of each bulge stellar particle, within different radii. As can be seen in Fig.  8(c), 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(c)), the gravitational potential of 'bulge stars' increases, in contrast to the simulation with AGN feedback.

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   Figure 9. The left-hand panel shows the contribution to the stellar density profile of in-, ex-situ and of the bulge 'relic' stars that existed at z = 6.7 already. Also shown is the total stellar density profile at z = 6.7 (dashed black line) and at z = 4.5 (solid black line). The right-hand panel shows the stellar radial distance travelled since z = 6.7 measured for star particles that reside within r < 1 kpc at z = 6.7 for the galaxy in Halo1 (black histogram) and the galaxy in the re-run simulation of Halo1 without AGN feedback (red histogram). The values give the percentage of star particles with D greater or smaller than 1 kpc. When AGN feedback operates, most stars that make-up the initially compact quasar host galaxy migrate outwards by several kpc. If AGN feedback is neglected, outwardly stellar migration is far less significant.
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 thresh = 3× 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 thresh . This radial threshold corresponds to ≈ 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 thresh that already existed at z = 6.7 within R 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 10 kpc, making up ≈ 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 ≈ 23% of the total stellar mass within R thresh , while the relic stars make up ≈ 18% of the final mass within R thresh . ≈ 4% of the stellar mass is not found using any these definitions, because they lie outside R 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 ≈ 60% to the total stellar mass within R 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 ≈ 4.5 (∆R = formation radius -final radius). We see that ≈ 39% of the in-situ stars migrate inwards and have a greater formation radius than their final radius. In contrast, ≈ 61% of the in-situ formed stars migrate outwards where ≈ 1% of the outward migrated stars have a significantly bigger final radius (∆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 atz = 6.7, and the green solid line, the profile at z = 4.5. About 86% of the relic stars move to larger radii and ≈ 5% of the outward-migrating stars have a significantly bigger final radius (∆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.

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 pc at z ≈ 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(Peirani et al. , 2019Costa 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 ≈ 500 pc stellar bulge. At this point, the velocity dispersions of the stellar bulges are high ≈ 500 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 ≈ 5 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 BH −σ relation seen for σ 270 km s −1 (e.g. Kormendy & Ho 2013). In our simulations, we find a drop in the peak circular velocity from 500 km s −1 to 300 km s −1 by z ≈ 3. Thus, while the black hole mass doesn't change significantly, the velocity dispersion drops considerably, resulting in an upturn in the M BH −σ relation.
The observed M BH − σ saturates at σ 270 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 BH −σ . One prediction of our scenario is that the characteristic stellar velocity dispersion at which the M BH − σ saturates should evolve with redshift. At z 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 BH − σ should persist even to velocity dispersions σ 400 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 σ 400 km s −1 should then start thinning out and an upturn in the M BH − σ 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(Debuhr et al. , 2012. Improving the assumed model is highly desirable to capture AGN feedback more realistically and test whether AGN-driven galaxy growth will still occur.

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 e f f ∝ M ). According to this scenario, the highredshift systems would be uniformly denser than their lowredshift 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 Myr timescales. The final mean size of ≈ 5 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(b)). 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 mi-  Figure 10. The left-hand panel shows the stellar half mass radius for the simulation with a constant physical softening length (Halo1-Soft), the simulation where AGN feedback is switched off at z = 5.2, (Halo1-weakAGN) and our fiducial simulation Halo1, as a function of redshift. Switching from constant comoving to constant physical softening lengths does not prevent the massive quasar host galaxy from growing in size. The right-hand panel shows the stellar radial distance travelled since z = 5.2 measured for star particles within the central region r < 1 kpc at z = 5.2 for Halo1 (black histograms) an d Halo1-weakAGN (red histograms). When AGN feedback is switched off, stellar outward migration ceases, showing that AGN feedback is directly responsible for the outward drift of stars in our simulations. After AGN feedback is switched-off, the tendency is for stars to migrate inwards as gas concentrates in the nucleus, deepening the potential well.  Figure 11. The left-and right-hand panel show the enclosed mass for different radii for the simulation where AGN feedback is switched off at z = 5.2 (dashed lines) and our fiducial simulation (solid lines). The left-hand panel shows the gas component and the right-hand panel shows the stellar component. Clearly, the enclosed gas mass stays highly variable after z ≈ 5.2 within ∼ kpc scales in the fiducial simulation. The simulation where AGN feedback is switched off stops showing this high variability. The same can be seen in the stellar component although the enclosed stellar mass shows smaller fluctuations. This shows that even the fluctuations within 300 pc (close to the resolution of the simulation) are physically meaningful and driven by AGN driven outflows. grate 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 ∼ 10 11 M . They subsequently grow in size by ≈ 1.5 dex down to z = 0 as well as in mass, in agreement with the r ∝ M 2 relation expected from minor dry mergers. Genel et al. (2018), however, find that quenched galaxies that are less massive (∼ 10 10.5 M at z = 0) follow a sizemass relation which is significantly steeper than r ∝ M 2 , which may be partly explained by AGN feedback-induced size growth.

Cusp/core transformation
The flattening in density profiles of both the stellar and dark matter component can be seen in Fig. 8(b) 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 ≈ 1 kpc (see Quillen  . Stellar and dark matter density profiles for Halo1 indicated with dashed and solid lines, respectively. The profiles are shown at z = 6.2 and z = 3.3 in red and green, respectively. The vertical grey region marks the range of the gravitational softening lengths found between z = 6.2 and z = 3.3. The stellar density profile develops a prominent ∼ kpc sized core by z = 3.3 if AGN feedback is efficient. 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 ∼ 10 10 M 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 BH = 10 9 M 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. (2017Peirani et al. ( , 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 × 10 11 M -1 × 10 13 M ) and the associated galaxies. They find that the density profiles significantly evolve with time. At z ≈ 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 ≤ 2).

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 km s −1 and scale radii of ≈ 300 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 ≈ 300 pc at z ≈ 6 to ≈ 5 kpc at z ≈ 3.
• Size growth is accompanied by a drop in the circular velocity, which falls from high values in the range 300-−500 km s −1 to 200−300 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 2 are direct tests to the ability of AGNdriven outflows to expel the central ISM and may ultimately be used to narrow down on the coupling efficiency of AGN feedback.