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

# The Role of Black Hole Feedback on Size and Structural Evolution in Massive Galaxies

Ena Choi<fnref n="1" />,<fnref n="2" /> ID, Rachel S. Somerville<fnref n="2" />,<fnref n="3" />, Jeremiah P. Ostriker<fnref n="1" />,<fnref n="4" /> ID, Thorsten Naab<fnref n="5" /> ID, and Michaela Hirschmann<fnref n="6" /> ID

Received 2018 February 14; revised 2018 September 4; accepted 2018 September 6; published 2018 October 16

## Abstract

We use cosmological hydrodynamical simulations to investigate the role of feedback from accreting black holes in the evolution of the size, compactness, stellar core density, and specific star formation of massive galaxies with stellar masses of $M_{*} > 10^{10.9} M_{\odot}$. We perform two sets of cosmological zoom-in simulations of 30 halos to $z = 0$: (1) without black holes and active galactic nucleus (AGN) feedback and (2) with AGN feedback arising from winds and X-ray radiation. We find that AGN feedback can alter the stellar density distribution, reduce the core density within the central 1 kpc by 0.3 dex from $z = 1$, and enhance the size growth of massive galaxies. We also find that galaxies simulated with AGN feedback evolve along tracks similar to those characterized by observations of specific star formation rate versus compactness. We confirm that AGN feedback plays an important role in transforming galaxies from blue compact galaxies into red extended galaxies in two ways: (1) it effectively quenches the star formation, transforming blue compact galaxies into compact quiescent galaxies, and (2) it also removes and prevents new accretion of cold gas, shutting down in situ star formation and causing subsequent mergers to be gas-poor or mixed. Gas-poor minor mergers then build up an extended stellar envelope. AGN feedback also puffs up the central region through fast AGN-driven winds as well as the slow expulsion of gas while the black hole is quiescent. Without AGN feedback, large amounts of gas accumulate in the central region, triggering star formation and leading to overly massive blue galaxies with dense stellar cores.

**Key words:** galaxies: evolution – galaxies: formation – galaxies: structure – methods: numerical – quasars: supermassive black holes

## 1. Introduction

In the local universe, galaxies show a bimodal color distribution: red quiescent galaxies (QGs) with old stellar populations, and young, blue, and star-forming galaxies (SFGs; e.g., Kauffmann et al. 2003; Bell et al. 2004). The color bimodality has been observed at early epochs up to $z \sim 2-3$ (Faber et al. 2007; Brammer et al. 2009). Quiescence seems to be strongly correlated with structural properties, such as galaxy size and concentration, at all epochs since $z \sim 3$. At a given stellar mass, QGs typically have smaller half-light radii, and denser cores, than SFGs (Shen et al. 2003; Williams et al. 2010; Wuyts et al. 2011; Bell et al. 2012; Fang et al. 2013; Omand et al. 2014; van der Wel et al. 2014; Brennan et al. 2017; Pandya et al. 2017). The fraction of QGs in the SDSS is found to be correlated with the mass and velocity dispersion of the bulge (Bluck et al. 2014, 2016).

The sizes of QGs are observed to be much smaller at fixed stellar mass at $z \sim 2$ compared to local elliptical galaxies (Ferguson et al. 2004; Trujillo et al. 2004, 2007; Daddi et al. 2005; Longhetti et al. 2007; Toft et al. 2007; Buitrago et al. 2008; Cimatti et al. 2008; Franx et al. 2008; van Dokkum et al. 2008; Szomoru et al. 2010; van der Wel et al. 2014). These compact QGs (cQGs), also known as “red nuggets” (Damjanov et al. 2009), have sizes of about $R_{\text{eff}} = 1$ kpc, about four times smaller than those of local QGs of the same mass. Moreover, the number density of compact galaxies is much lower in the local universe than in the past (Taylor et al. 2010; Cassata et al. 2011), leading to much speculation as to the fate of compact galaxies and how these

different populations of compact, extended, SF, and quiescent galaxies are related to one another at different cosmic epochs.

Barro et al. (2013) suggested that massive galaxies evolve through a characteristic track in the plane of the specific star formation rate (sSFR) versus compactness (defined as $\Sigma_{1.5} \equiv m_*/r_e^{1.5}$, where $m_*$ is the galaxy stellar mass and $r_e$ is the half-light radius; see their Figure 6): beginning at around redshift $z \sim 2.5-3$, compact SFGs (cSFGs) begin to quench, building up a population of cQGs. Below $z \sim 1$, galaxies increase in mass and (even more so) in size to form the extended QGs that are common today. The results of Barro et al. (2013) furthermore suggest that the quenching time is relatively short (300 Myr to 1 Gyr; see also Pandya et al. 2017). Barro et al. (2014) showed that the radii and stellar mass surface densities of cSFGs quantitatively matched those of cQGs, supporting the picture in which cSFGs rapidly quench into cQGs.

A similar picture was suggested by van Dokkum et al. (2015), who also noted the concurrent decline of the cSFG population and rise of the cQG population. They proposed a toy model in which SFGs evolve along parallel tracks in the size-mass plane with $\Delta r \propto \Delta m^{0.3}$, until they reach a stellar density or velocity dispersion threshold and quench their star formation. After quenching, galaxies evolve along a steeper track in the size-mass plane, with $\Delta r \propto \Delta m^2$ as observed (van Dokkum et al. 2010; Patel et al. 2013).

These observationally inferred formation scenarios pose several major questions to theoretical galaxy formation models:

<footnote n="1">¹Department of Astronomy, Columbia University, New York, NY 10027, USA; ena.choi@columbia.edu</footnote>

<footnote n="2">²Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854, USA</footnote>

<footnote n="3">³Simons Center for Computational Astrophysics, New York, NY 10010, USA</footnote>

<footnote n="4">⁴Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="5">⁵Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85741 Garching, Germany</footnote>

<footnote n="6">⁶Sorbonne Universités, UPMC-CNRS, UMR7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>
which physical processes drive the characteristic evolution of the structural properties of star-forming galaxies and QGs? What physical process is responsible for the quenching of star formation? Are these two sets of processes causally connected or the result of a common cause, or do they occur together only coincidentally? What causes the change in slope of the evolutionary path in the size-mass plane, or put another way, why do QGs grow so much more rapidly in radius relative to their mass growth than SFGs?

Several potential reasons for the early development of compact galaxies have been addressed in many papers, such as intense starbursts in galactic nuclei triggered by gas-rich major mergers (Barnes & Hernquist 1991; Hopkins et al. 2006; Wellons et al. 2015), strong gas inflows within the disk driven by violent disk instabilities (Gammie 2001; Dekel et al. 2009), and early dissipative assembly (Oser et al. 2010). Different authors have applied different labels to this phase of galaxy formation, but all pictures involve rapid inflow of gas into central star-forming regions (see reviews by Somerville & Davé 2015; Naab & Ostriker 2017). Some works have suggested a causal connection between these dissipative processes (“compaction”) and quenching via a combination of gas exhaustion, and stellar and active galactic nucleus (AGN) driven winds (Dekel & Burkert 2014; Zolotov et al. 2015). Observationally, strong galactic outflows are often associated with high star formation and gas densities (e.g., Diamond-Stanic et al. 2012).

The Λ Cold Dark Matter (ΛCDM) hierarchical structure formation model ubiquitously predicts that galaxies experience many mergers over the course of more than 11 Gyr from $z \sim 3$ to the present day (Conselice et al. 2009; Stewart et al. 2009; Lotz et al. 2011; López-Sanjuan et al. 2013; Rodriguez-Gomez et al. 2016). Most of these would be termed “minor mergers” with $\Delta m_*/m_* < 1/4$. However, gas-rich (“wet”) mergers can drive gas into the nucleus of galaxies, making the centers denser and more compact; gas-poor (“dry”) mergers produce remnants that are more extended than the progenitors (Hopkins et al. 2010; Covington et al. 2011). Several works have pointed out that a promising way to drive rapid size growth with relatively little growth in mass is for a galaxy to experience multiple dry, minor mergers (Bezanson et al. 2009; Naab et al. 2009; Oser et al. 2012; Hilz et al. 2013; Porter et al. 2014).

To summarize, the suggested physical sequence (as described by Barro et al. 2013) is the following: (1) star-forming, extended galaxies are gas-rich and undergo dissipative processes to form cSFGs. (2) This compaction is followed by the rapid quenching of star formation to form cQGs. Finally, (3) the galaxy sizes are gradually expanded via multiple, minor dry mergers. Although this picture is plausible and appealing, it is still an open question whether fully self-consistent numerical cosmological hydrodynamical simulations can qualitatively or quantitatively reproduce the observed evolution.

Recent cosmological hydrodynamic simulations have demonstrated that AGN feedback can reduce galaxy stellar mass by preventing cooling flows (cf. Fabian 1994) and quenching subsequent star formation in group and cluster-sized halos (see reviews by Somerville & Davé 2015 and Naab & Ostriker 2017 for references). Several recent studies have examined predictions for the evolution of the size-mass relations for both star-forming galaxies and QGs in hydrodynamic simulations of large cosmological volumes that include various implementations of AGN feedback (Genel

et al. 2014, 2018; Crain et al. 2015; Wellons et al. 2015; Furlong et al. 2017).

While these studies found qualitative agreement with the observations, predictions for galaxy size and structure are known to be extremely sensitive to the details of how processes such as star formation, stellar feedback, and AGN feedback are implemented on subresolution scales (“subgrid” physics; see discussion and references in Somerville & Davé 2015). For example, the predicted size-mass relation for SFGs is quite different in the original Illustris simulation and the revised Illustris-TNG simulation (Genel et al. 2018), but the physical reasons for this are not fully understood (S. Genel 2018, private communication).

It has been shown that simulations that include AGN feedback produce considerably better agreement with the observed size-mass relation for massive QGs, while simulations that omit AGN feedback produce massive galaxies that are too compact (Dubois et al. 2013; Crain et al. 2015; Choi et al. 2017). However, once again, the physical processes at play are not fully understood. AGN feedback clearly reduces the gas fraction in early massive galaxies, leading to more dry merging, which, as discussed above, can drive rapid size growth. Furthermore, AGN feedback shuts off the late in situ star formation that tends to rebuild a compact stellar core.

In addition, a potentially important dynamical impact of AGN feedback on the dark matter and stellar distributions has been discussed recently. Martizzi et al. (2012, 2013) found that AGN feedback can produce a flat stellar and dark matter density core in a galaxy-cluster-scale zoom-in simulation, possibly via the coupling of several mechanisms such as orbital energy transfer by black holes through dynamical friction and “central revirialization” (Pontzen & Governato 2012) after AGN-driven gas outflows.

In this study, we make use of the high-resolution simulations of massive galaxies evolving in a full ΛCDM context presented by Choi et al. (2017) to study the evolution of the size, stellar mass, star formation rate, and compactness of galaxies during the critical transition phase $0 < z < 3$. These simulations utilize the “zoom-in” technique to attain mass and spatial resolution higher than those for large cosmological volumes. Furthermore, these simulations include a more detailed and physically motivated treatment of the feedback from radiatively efficient accretion onto supermassive black holes (see Section 2.5 in Choi et al. 2017) than previous cosmological simulations that have studied galaxy structural evolution. Our feedback model is motivated by observations of powerful winds seen in broad absorption line (BAL) AGNs (Gibson et al. 2009; Moe et al. 2009; Crenshaw & Kraemer 2012; Arav et al. 2013). We include mechanical feedback from the BAL winds, as well as radiative feedback via Compton and photoionization heating (Choi et al. 2012, 2015). Furthermore, we analyze a matched suite of simulations with and without AGN feedback, and delve in more detail into the physical processes responsible for driving the predicted structural evolution.

We describe our simulations in Section 2, present our results in Section 3, and conclude in Section 4.

## 2. Simulations

We performed two sets of high-resolution, cosmological zoom-in hydrodynamic simulations of 30 massive halos, with present-day total masses of $1.4 \times 10^{12} M_\odot \le M_{vir} \le 2.3 \times 10^{13} M_\odot$ and present-day stellar masses of $8.2 \times 10^{10} M_\odot \le M_* \le$
1.5 × 10¹² M☉ (for the central galaxies). The physics implemented in our simulation models includes star formation, supernova feedback, wind feedback from massive stars, asymptotic giant branch (AGB) stars, and metal cooling and diffusion. Our simulation code also incorporates a new treatment of mechanical and radiative AGN feedback, which is implemented in a self-consistent way, launching high-velocity mass outflows. The simulation suite used in this study is presented in Choi et al. (2017), and we refer the reader to that paper for further details. In the following, we briefly summarize our simulations.

## 2.1. The Hydrodynamic Simulation Code: SPHGal

We use SPHGal (Hu et al. 2014), a modified version of the parallel smoothed particle hydrodynamics (SPH) code GADGET-3 (Springel 2005). This code incorporates a density-independent pressure-entropy SPH formulation (Ritchie & Thomas 2001; Hopkins 2013; Saitoh & Makino 2013) to overcome the numerical fluid-mixing problems of classical SPH codes (e.g., Agertz et al. 2007). It also includes an improved artificial viscosity implementation (Cullen & Dehnen 2010), an artificial thermal conductivity (Read & Hayfield 2012), and a Wendland $C^4$ kernel with 200 neighboring particles (Dehnen & Aly 2012). In order to ensure proper treatment of shock propagation and feedback distribution, we employ a time-step limiter that makes neighboring particles have similar time steps (Saitoh & Makino 2009; Durier & Dalla Vecchia 2012).

## 2.2. Star Formation, Chemical Enrichment, and Stellar Feedback Model

Following Aumer et al. (2013), stars are stochastically formed within overdense regions when the gas density exceeds a density threshold for the Jeans gravitational instability of the enclosed mass. The star formation rate is computed as $d\rho_*/dt = \eta\rho_{gas}/t_{dyn}$, where $\rho_*$, $\rho_{gas}$, and $t_{dyn}$ are the stellar and gas densities, and local dynamical time of gas particle, respectively, and we set the star formation efficiency as $\eta = 0.025$.

The evolution of each star particle contributes to the chemical enrichment of the interstellar medium (ISM) in our simulation during various mass-loss events. We allow chemical enrichment via winds driven by Type I supernovae (SNe), Type II SNe, and asymptotic giant branch (AGB) stars with the chemical yields adopted from Iwamoto et al. (1999), Woosley & Weaver (1995), and Karakas (2010). We explicitly trace the mass in 11 chemical elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe), both for star and gas particles. We also allow the metal-enriched gas particles to mix their metals with neighboring metal-poor gas particles via turbulent diffusion, using the standard SPH neighbor searches following Aumer et al. (2013). The net cooling rate is calculated based on the individual chemical abundances, temperatures, and densities of gas particles following Wiersma et al. (2009), accounting for a redshift-dependent metalalactic UV/X-ray and cosmic microwave background with a modified Haardt & Madau (2001) spectrum.

The stellar feedback model is adopted from Núñez et al. (2017) and includes UV heating within Strömgren spheres around young massive stars, three-phase SN feedback by both Type I and Type II SNe, and winds from dying low-mass AGB stars. Each of the mass-loss events of the stellar particles

explicitly contributes mass, metals, momentum, and energy to the surrounding gas. First, the young star particles gradually heat the neighboring gas to $T = 10^4$ K within their H II region limit (Strömgren spheres; Strömgren 1939) before they explode as SNe.

We assume that a single SN event ejects mass in an outflow with a velocity $v_{SN} = 4500$ km s⁻¹, a typical velocity of outflowing SN ejecta (Janka 2012). Depending on the physical distance from the SN, each adjacent gas particle is affected by one of the three successive SN phases: (i) momentum-conserving free expansion phase, (ii) energy-conserving Sedov–Taylor phase where SN energy is transferred with 30% as kinetic and 70% as thermal, and (iii) the snowplow phase where radiative cooling becomes important. In this “snowplow” SN feedback model, each SN remnant launches standard Sedov–Taylor blast waves carrying energy as 30% kinetic and 70% thermal, and both dissipate with distance from the SN in its final pressure-driven snowplow phase. In addition, the old stellar particles still contribute their mass, energy, and metal output via slow winds during an AGB phase. The energy and metal-enriched mass output from the old star particles are transferred to the adjacent gas particles in a momentum-conserving manner.

We assume a Kroupa (2001) initial mass function (IMF), and we note that over 30% of the total mass in stars will be continuously ejected into the ISM via winds from SNe and AGB stars within ~13 Gyr of stellar evolution. This metal-enriched mass output from stars not only fuels late star formation but also feeds the central supermassive black hole, inducing AGN activity (see also Ciotti et al. 2010).

## 2.3. Black Hole Formation, Accretion, and AGN Feedback Model

In the simulations, new black hole seeds of mass $10^5 h^{-1} M_{\odot}$ are placed in the center of galaxies when they initially reach a dark matter halo mass of $10^{11} h^{-1} M_{\odot}$. The black hole seed mass and dark matter halo threshold mass are chosen to approximately follow the Magorrian et al. (1998) relation and the theoretical calculations of black hole seed formation (e.g., Volonteri 2010; Stone et al. 2017). We note that black hole seed mass makes a negligible contribution to the final black hole mass of the central galaxies in our simulations.

The black hole can grow by merging with other black holes and by gas accretion. When two black hole particles get closer than their local SPH smoothing lengths and their relative velocities are less than the local sound speed, they merge together. The black hole particle also grows in mass via the accretion of surrounding gas with a Bondi–Hoyle–Lyttleton accretion rate (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), $\dot{M}_{inf} = (4\pi G^2 M_{BH}^2 \rho) / (c_s^2 + v^2)^{3/2}$, where $\rho$, $c_s$, and $v$ are the density, the sound speed, and the velocity of the gas relative to the black hole, respectively. In order to prevent the unphysical accretion of unbound gas outside the Bondi radius, we incorporate the soft Bondi criterion first introduced in Choi et al. (2012). We statistically limit the accretion of the gas by the volume fraction of the gas particle lying within the Bondi radius; for example, the full accretion of the particle is allowed only when the total volume of a smoothed gas particle is contained within the Bondi radius.

The AGN feedback model we use is adopted from Choi et al. (2012, 2014) and consists of two main components: (1) mechanical feedback via winds, which carry energy, mass, and
momentum into the neighboring gas, and (2) radiative feedback via Compton and photoionization heating from the X-ray radiation from the accreting black hole. We also incorporate the radiation pressure associated with the X-ray heating, and the Eddington force. The emergent AGN spectrum and metal-line heating rate are taken from Sazonov et al. (2004). The winds are driven by radiation pressure on gas and dust on scales below those that we can simulate explicitly; we therefore treat them using a subgrid model as outlined below.

For the mechanical AGN feedback, the gas mass inflowing to the central region contributes to the accretion onto the black hole and AGN-driven winds. AGN winds carry mass with the kinetic energy rate given as $E_{\text{outf,AGN}} = 0.5\dot{M}_{\text{outf}}v_{\text{outf,AGN}}^2$, where the AGN outflowing wind velocity is assumed to be $v_{\text{outf,AGN}} = 10,000 \text{ km s}^{-1}$, motivated by observations of BAL winds (e.g., Arav et al. 2013). This kinetic energy rate is proportional to the gas mass accreted onto the black hole as $\dot{E}_{\text{outf,AGN}} \equiv \epsilon_w \dot{M}_{\text{acc}} c^2$, where $\epsilon_w$ indicates the AGN feedback efficiency and is set to be $\epsilon_w = 0.005$, similar to the values used in other AGN feedback model implementations (e.g., Di Matteo et al. 2005).

Based on our selection of the feedback efficiency $\epsilon_w$ and the wind velocity $v_{\text{outf,AGN}}$, a fraction of the gas particles entering the accretion region is stochastically selected as wind particles. The selected wind particles are kicked in a direction parallel or antiparallel to the direction of angular momentum of each gas particle; therefore, the wind tends to be oriented perpendicular to the disk plane (Proga & Kallman 2004) when the black holes are surrounded by a rotating gas disk. The ejected wind particle shares its energy and momentum with two adjacent gas particles and produces a shock-heated momentum-driven flow with a ratio of kinetic to thermal energy similar to that in the standard Sedov-Taylor blast wave.

We also incorporate the heating via hard X-ray radiation from the accreting black hole following Sazonov et al. (2005). At the position of each gas particle, we calculate the net luminosity flux from all black holes in the simulated zoom-in area. The calculated flux is then converted to the net volume heating rate $\dot{E}$ via Compton and photoionization heating using the Sazonov et al. (2005) formulae. The radiation pressure from the X-ray flux is also included, as every gas particle absorbing energy $\Delta E$ from X-ray radiation is given an additional momentum $\Delta p = \Delta E/c$ directed away from the black hole.

Finally, we include the Eddington force on electrons in the neighboring gas directed radially away from black holes instead of artificially limiting the gas accretion rate onto the black hole to not exceed the Eddington rate. In our simulations, we allow for super-Eddington gas accretion to occur occasionally. When this happens, the corresponding feedback effects naturally reduces the inflow and increases the outflow.

Choi et al. (2015) showed that the mechanical AGN feedback via BAL winds produces stronger effects on the galaxy compared to the thermal feedback treatment (e.g., Springel et al. 2005). The feedback energy injected via mechanical winds rather than that in thermal form more effectively drives the gas out of the galaxy via galactic outflows (Brennan et al. 2018) and quenches the star formation for a longer timescale.

## 2.4. Initial Conditions

We use the cosmological “zoom-in” initial conditions described in Oser et al. (2010). A subvolume for the zoom-in

<table><caption>**Table 1**
Summary of the Simulation Resolution</caption><thead><tr><th>Mass Ratio</th><th>Softening Length <i>&epsilon;</i> (pc h<sup>-1</sup>)</th><th>Particle Mass (<i>M</i><sub>&odot;</sub>)</th></tr></thead><tbody><tr><td>Dark matter</td><td>890</td><td>3.5 &times; 10<sup>7</sup></td></tr><tr><td>Baryon</td><td>400</td><td>5.8 &times; 10<sup>6</sup></td></tr></tbody></table>

initial condition is extracted from a larger volume, dark-matter-only simulation using a flat cosmology with parameters obtained from WMAP3 (Spergel et al. 2007, $h = 0.72$, $\Omega_b = 0.044$, $\Omega_{\text{dm}} = 0.216$, $\Omega_\Lambda = 0.74$, $\sigma_8 = 0.77$, and $n_s = 0.95$). All dark matter particles close to the halos of interest are traced from redshift zero and then replaced with higher resolution gas and dark matter particles.

We resimulate the new high-resolution initial conditions from redshift $z = 43$ to $z = 0$. The simulation sets used in this study have mass resolution for the baryonic particles (star and gas) of $m_{*,\text{gas}} = 5.8 \times 10^6 M_{\odot}$ and the dark matter particles of $m_{\text{dm}} = 3.5 \times 10^7 M_{\odot}$. We adopt the co-moving gravitational softening lengths $\epsilon_{\text{gas,star}} = 400 \text{ pc h}^{-1}$ for the baryonic particles and $\epsilon_{\text{halo}} = 890 \text{ pc h}^{-1}$ for the dark matter. The adopted resolution is summarized in Table 1.

The simulated halo masses range from $1.4 \times 10^{12} M_{\odot} \lesssim M_{\text{vir}} \lesssim 2.3 \times 10^{13} M_{\odot}$ at $z = 0$, and the stellar masses of central galaxies are $8.2 \times 10^{10} M_{\odot} \lesssim M_{*} \lesssim 1.5 \times 10^{12} M_{\odot}$ at present day. These galaxies are well resolved with $\approx 2.5 \times 10^4-4.8 \times 10^5$ stellar particles within the virial radius ($R_{200}$, the radius where the spherical overdensity drops below 200 times the critical density of the universe at a given redshift). The halo virial mass $M_{200}$ is the mass contained within a sphere with radius $R_{200}$.

In order to study the effects of AGN feedback on the star formation, size, and compactness of the galaxies, we run the full set of simulations with two different sets of physics:

(1) **NoAGN:** No black holes and no AGN feedback. This model is comparable to the results of Oser et al. (2010, 2012), but note that we use an alternative density-independent formulation of SPH designed to treat contact discontinuities more accurately and include an artificial viscosity and an energy diffusion implementation which results in more cooling. We also include metal enrichment and metal line cooling, which also enhance the cooling process. Moreover, we adopt a different recipe for stellar feedback (see Choi et al. 2017, for a further exploration of the effects of these changes).

(2) **WithAGN:** all of the same physical models as in the NoAGN model, with the addition of black holes and mechanical and radiative AGN feedback.

In this study, we only consider the central galaxies within the simulated halos. The stellar mass and the star formation rate of the galaxy are measured within 10% of the virial radius $r_{10} = 0.1 \times R_{200}$. Then, we determine the effective radius of the galaxy $r_{\text{eff}}$ by determining the mean values of the half-mass radii of stars within $r_{10}$ projected along 20 randomly chosen directions of the main stellar body.

## 3. Results

In Choi et al. (2017), we showed that AGN feedback, as implemented in our simulations, can change the fate of massive galaxies, effectively quenching late star formation, turning the majority of massive galaxies into QGs, and reducing their final
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Size evolution of simulated galaxies with AGN feedback (green circles) and without AGN feedback (purple triangles) from $z=3$ to $z=0$. The observed size–mass relation of present-day quiescent galaxies from Nipoti et al. (2009) is shown by the red solid line (with the $1\sigma$ scatter in red dotted lines). The size–mass relation for high-redshift galaxies ($z=2.25$) from CANDELS (van der Wel et al. 2014) is shown by the red dotted line. (Left) Projected stellar half-mass radii vs. stellar mass for simulated galaxies shown at $z=3, 2, 1$, and 0. Quiescent galaxies with low specific star formation rates (sSFR $\le 0.3 t_{\text{H}}$) are shown by the solid symbols, while star-forming galaxies (sSFR $\ge 0.3 t_{\text{H}}$) are shown by the open symbols. Galaxies simulated with AGNs have their star formation quenched starting at $z \sim 2$, and then start to evolve on a steeper track on the size–mass relation. (Right) Average tracks of simulated galaxies on the size–mass plane from $z=3$ to $z=0$ are shown for two models, with AGN feedback (green arrow) and without AGN feedback (purple arrow). The points show the size–mass relation of the simulated galaxies at $z=3$ with open symbols and at $z=0$ with filled symbols, respectively.</figcaption></figure>

stellar mass by almost a factor of 4. AGN feedback effectively reduces the in situ star formation in the central galaxy where the most massive black holes are located. Our simulated galaxies without AGNs show much higher rates of in situ star formation in the innermost regions of the galaxies, as the gas reservoir is constantly refilled by “recycled gas” from SN and AGB winds from the old stellar population. In this section, we show that AGN feedback also affects various structural properties of galaxies, including galaxy sizes, stellar core density, and compactness.

### 3.1. Evolution in the Size–Mass Plane

In left panels of Figure 1, we show the evolution of the galaxies simulated with and without AGN feedback in the projected half-mass–radius (“size”) versus stellar mass plane from redshifts $z=3$ to $z=0$. We show the sizes and stellar masses of the 30 central galaxies at $z=0$ and their most massive progenitors at $z=1, 2$, and 3. In order to exclude outliers caused by galaxy interactions (e.g., ongoing mergers), we use the median value of the sizes that we obtained for three snapshots around the target redshift. We divide our galaxies into star-forming galaxies (open symbols) and QGs (solid symbols) based on the commonly adopted dividing line at sSFR = $0.3/t_{\text{H}}$, where $t_{\text{H}}$ is the age of the universe at each redshift (e.g., Franx et al. 2008). We adopt this method as this divider based on sSFR is roughly equivalent to the commonly adopted method that separates QGs and SFGs in color–color space up to $z \sim 2.5$ as described in Williams et al. (2009, 2010) and Muzzin et al. (2013). Note that we only consider the central galaxies in the simulated halos. The central galaxies in group environments may experience accelerated size growth

compared to galaxies in the field via minor mergers, which are more common in high-density environments (see Cappellari 2013; Chan et al. 2018).

At $z=3$, all progenitor galaxies simulated with and without AGN feedback are very compact ($r_{\text{eff}} < 1$ kpc), and the two models do not show a significant difference in size while NoAGN galaxies are slightly more massive than WithAGN galaxies. At $z=2$, both WithAGN and NoAGN galaxies are still compact, compared to the size–mass relation for high-redshift galaxies ($z=2.25$) from CANDELS (van der Wel et al. 2014). However, by $z=1$, the quenched galaxies among those simulated with AGN feedback start to diverge and form a size–mass relation clearly separated from that of NoAGN feedback galaxies. By $z=0$, all galaxies are quenched in the simulations with AGN feedback, and they closely follow the observed size–mass relation of the local QGs from the SLACS sample (Nipoti et al. 2009).

On the other hand, NoAGN galaxies show much less size growth since $z \sim 2$, evolving with a much shallower slope in the size–mass plane. The present-day sizes of galaxies simulated without AGN feedback are ~5 times smaller compared to the observed relation at a given stellar mass, consistent with Crain et al. (2015). In the absence of AGN feedback, we have continuous star formation in the central regions of galaxies, which results in a concentrated stellar mass profile as discussed in Martizzi et al. (2012). In addition, the “outer” size growth via the addition of stars from minor dry mergers is less significant in the NoAGN simulations, as in situ star formation always dominates over the accreted star component (see Figure 9 in Choi et al. 2017, and also Lackner et al. 2012; Dubois et al. 2013, 2016, and Hirschmann et al. 2013).
Since $z \sim 2$, the galaxies without AGN feedback mainly accumulate their stellar mass via in situ star formation (stars formed from gas within the main progenitor), while the galaxies with AGN feedback mainly grow their stellar mass by accreting stellar mass via minor mergers, after their in situ star formation is quenched. This difference in the evolutionary paths of galaxies in the two models in the size-mass plane is illustrated in the right panel of Figure 1, where we show the average tracks of simulated galaxies for the two models. The two models with and without AGN feedback follow very different evolutionary tracks. Until $z \sim 2$, galaxies in both models grow in mass and gradually increase their density, showing little growth in size. Then, after the in situ star formation is quenched (as shown in left panels), the galaxies simulated with AGN feedback show a much steeper evolution in the size-mass plane.

We note that the galaxies simulated with AGN feedback follow a very similar track to that suggested based on observations by van Dokkum et al. (2015), as discussed in the introduction. SFGs in our models evolve along the shallower track in the size-mass plane, $\Delta r \propto \Delta M^{0.3}$, until their star formation is quenched by the AGN-driven winds. After quenching, galaxies grow along a steeper track in the size-mass plane, with $\Delta r \propto \Delta M^2$. On the other hand, galaxies simulated without AGN feedback remain gas-rich and continue to evolve along the much shallower track with $\Delta r \propto \Delta M$.

Recently, Genel et al. (2018) also found such a turnover in the size-mass plane in the evolution of massive galaxies in the Illustris-TNG simulation. They found that galaxies have this steep size evolution as a function of the added stellar mass after they quench their star formation—mainly induced by their black hole driven kinetic and thermal feedback (Weinberger et al. 2017) as in this study—and the size growth of the most massive galaxies occurs mostly during their quiescent phase.

To quantify the size evolution of our simulated galaxies as a function of redshift, we fit a power law in $(1+z)^\alpha$ as is frequently done in observational studies (e.g., Newman et al. 2010; Damjanov et al. 2011). We show the size evolution of galaxies with $M_* > 6.3 \times 10^{10} M_\odot$ from $z=2$ to the present day in Figure 2. The median sizes of galaxies in each redshift bin are shown as squares, and we show the power-law fit to the median sizes with black lines. Overall, galaxy sizes are smaller in the NoAGN model at all redshifts, and the difference becomes even larger by $z=0$ as WithAGN galaxies show more rapid size growth per added stellar mass during the quenched phase as shown in Figure 1. Between redshift 2 and 0, our simulated galaxies grow on average by a factor of 5. Our power-law fit to galaxy sizes simulated with AGNs is $\alpha = -1.49$, which is in good agreement with observations (van der Wel et al. 2014).

In Figure 3, we show the projected stellar half-mass radii at $z=0$ as a function of the fraction of stellar mass formed in other galaxies and accreted onto the main progenitor, to the stellar mass formed in situ, for galaxies simulated with and without AGN feedback. There is a clear correlation between the sizes and the relative amount of accreted and in situ formed stars for WithAGN as well as NoAGN feedback models. This indicates that the accretion of stellar mass drives the growth of the effective radius of massive galaxies for both models. For galaxies whose mass is dominated by the accreted stellar component $M_{\text{acc}} > M_{\text{in situ}}$, the two models show similar sizes at a fixed accreted star fraction. However, for galaxies

dominated by in situ formed stars, $M_{\text{in situ}} > M_{\text{acc}}$, galaxies simulated without AGNs show smaller sizes at a fixed fraction. This is partly because of the fact that completely dry mergers are relatively rare due to the high gas fractions in satellite galaxies without AGN feedback, while merging satellite galaxies also have black holes and corresponding AGN feedback in the WithAGN feedback model. Moreover, galaxies with AGNs are also more efficiently puffed up via gas outflows (Fan et al. 2008). The slow expulsion of gas while the AGN is quiescent can induce adiabatic expansion of the stellar component (e.g., Hills 1980), and the rapid mass loss driven by fast AGN winds can further elevate this puffing-up process (Biermann & Shapiro 1979). We discuss and illustrate this puffing-up process of the stellar system in the next subsection.

## 3.2. Central Stellar Density Evolution

Cheung et al. (2012) proposed the quantity $\Sigma_1$, the projected stellar mass surface density within the central 1 kpc, as a robust measure of galaxies’ core density, and showed that at $z \sim 0.5-0.8$, $\Sigma_1$ is more strongly correlated with quiescence than stellar mass or the effective surface density $\Sigma_{\text{eff}}$. Fang et al. (2013) showed that $\Sigma_1$ shows a tight scaling relation with stellar mass for QGs at $z \sim 0.1$, and that a stellar-mass-dependent $\Sigma_1$ threshold is a better predictor of quiescence than stellar mass alone. Barro et al. (2017) studied the structural properties of massive galaxies in the CANDELS/GOODS-S field and found that this tight scaling relation in $m_*$ versus $\Sigma_1$ holds for QGs out to $z \sim 2.5$, but the zero point of this relation has declined by $\sim 0.3$ dex since $z=2$. Van Dokkum et al. (2014) also found that the number density of galaxies with dense cores has decreased from $z=2.5$ to today. Both papers suggested that stellar mass loss and subsequent adiabatic expansion can explain this observed decrease in core density.

In order to study the evolution of the stellar core density of our simulated galaxies, we measure the projected central stellar mass surface density within $r<1$ kpc, $\Sigma_1 \equiv M_*(<1 \text{ kpc})/\pi(1 \text{ kpc})^2$, and show the evolution of all galaxies in the stellar mass versus $\Sigma_1$ plane from $z=3$ to $z=0$ in Figure 4. The core densities $\Sigma_1$ indicated here are the mean values of the surface stellar densities projected along the three principal axes. We note that we only focus on the differential effect as this quantity is affected by gravitational softening.

In the left panels of Figure 4, we show $\Sigma_1$ as a function of stellar mass for 30 central galaxies at $z=0$ and their most massive progenitors at $z=1, 2$, and 3 as in Figure 1. We also separately show the star-forming galaxies and QGs (open and filled symbols, respectively) to study the effect of star formation quenching on the evolution of the core stellar density. Until $z \sim 2$, galaxies in both models showed a tight correlation between $\Sigma_1$ versus stellar mass and did not show significant differences in the relationship. By $z=1$, the quenched galaxies among those simulated with AGN feedback started to show decreased core densities and diverged from the previously established relation, while galaxies without AGN feedback continued to increase their core densities as well as stellar mass.

All galaxies are quenched in the WithAGN model by the present day, and they form a clearly separated relation from that of the NoAGN galaxies. While WithAGN galaxies are in better agreement with the observed relation from Barro et al. (2017) and Fang et al. (2013), NoAGN galaxies show $\sim 0.5$ dex higher
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Redshift evolution of the projected stellar half-mass radii of galaxies with $M_* > 6.3 \times 10^{10} M_\odot$. The galaxies simulated with AGNs are shown by the solid green (quiescent) and open green (star-forming) circles, and the galaxies simulated without AGN feedback are shown by the solid purple (quiescent) and open purple (star-forming) triangles. The green and purple squares show the mean sizes at a given redshift for the WithAGN and NoAGN models, respectively. NoAGN model galaxies are offset by 0.1 in redshift for clarity. The black lines show the result of a power-law fit for the NoAGN (dashed) and WithAGN (dotted) simulations.</figcaption></figure>

core densities at a given stellar mass than observed. This high core density in NoAGN galaxies is attributed to the continuous in situ star formation of the galaxies, which is due to the lack of star formation quenching mechanisms in the central region. NoAGN galaxies have extended star formation until the present day, especially in the innermost region of the galaxies, as the gas reservoir in the galactic center is constantly refilled by “recycled gas” from SN and AGB winds in the old stellar population. Although WithAGN galaxies show much lower core densities than NoAGN galaxies, they are still 0.2–0.3 dex denser than the observations (Fang et al. 2013) at low redshift. This is presumably due to overcooling at high redshift as galaxies in both models are already too dense in the cores at $z=2$ and 3.

In the right panel of Figure 4, we show the average tracks of simulated galaxies in the Σ₁-stellar mass plane for the WithAGN and NoAGN models. Until $z \sim 2$, galaxies in both models grow their stellar mass as well as core densities following a tight correlation. NoAGN galaxies keep following this evolutionary track until $z=0$, increasing their core densities, though the slope becomes a bit shallower by the present day. The core densities of galaxies simulated with AGN feedback, however, start to decline when the in situ star formation is suppressed (as shown in the left panels). This transition seems to happen when galaxies reach $M_* \sim 10^{11} M_\odot$, where a majority of galaxies start having their star formation quenched by AGN feedback. Overall, the galaxies with AGN feedback show almost a 0.3 dex decrease in core density from redshift 1 to 0, consistent with recent observations by Barro et al. (2017).

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Projected stellar half-mass radii at $z=0$ as a function of the fraction of accreted stellar mass to the stellar mass formed in situ for the WithAGN model (green circles) and NoAGN model (purple triangles). Quiescent galaxies with low specific star formation rates (sSFR $\le 0.3\,\text{f}_\text{H}$) are shown by the solid symbols, and star-forming galaxies (sSFR $\ge 0.3\,\text{f}_\text{H}$) are shown by the open symbols. For both WithAGN and NoAGN models, we find a correlation between galaxy size and the fraction of accreted stars, i.e., galaxies with a higher fraction of accreted stars have larger sizes. This suggests that the accretion of stellar mass through mergers is a main driver of the size growth of massive galaxies.</figcaption></figure>

We find that, in our simulations, the decrease in core density is a result of a combination of several physical processes. One contributor is stellar mass loss, i.e., old stars gradually lose their mass by stellar evolution via SN and AGB winds. Without additional star formation, the central stellar density naturally decreases due to this effect. In our simulations, star particles gradually return mass to the adjacent gas particles via SN- or AGB-driven winds. In order to check how much mass stellar particles have lost due to stellar evolution, we trace the evolution of all central star particles within $r < 1$ kpc at $z=2$ over time until $z=0$. On average, these central star particles have lost $\sim 10\%$ of their mass since $z=2$ via stellar evolution. We note that with our assumed IMF, over 30% of the total mass in newly formed stars will eventually be returned to the gas phase via winds over $\sim 13$ Gyr of stellar evolution, but the bulk of the mass loss occurs when the stellar populations are very young, via SN explosions (see also Leitner & Kravtsov 2011). Therefore, the relatively old stellar population already in place in the central region of the galaxies only shows moderate stellar mass loss since $z \sim 2$, and this cannot fully explain the $\sim 50$ decrease in central stellar densities in WithAGN galaxies.

Second, the expansion of the collisionless particles in the central region after gas ejection by outflows can “puff up” the central region due to a decrease in the gravitational potential, leading to a more diffuse stellar core. The effect of outflows on the central stellar density depends on the amount of mass ejected as well as on the timescale of ejection. We have two channels for this “puffing-up” process, depending on the timescale of the mass loss: an adiabatic expansion after slow
<figure><img src="image_5.png" /></figure>

**Figure 4.** Projected stellar core density, $\Sigma_I = M_*(<1 \text{ kpc})/\pi(1 \text{ kpc})^2$ vs. stellar mass of simulated galaxies with AGN feedback (green circles) and without AGN feedback (purple triangles) from $z=3$ to $z=0$. The observed $\Sigma_I - M_*$ relation of quiescent galaxies is shown for $1.4 < z < 2.2$ (red dotted line), for $0.5 < z < 1.0$ (red solid line) from Barro et al. (2017), and for $z \sim 0.1$ (red dashed line) from Fang et al. (2013), respectively. (Left) Projected stellar core density within 1 kpc vs. stellar mass for simulated galaxies shown at $z=3, 2, 1$, and 0. Quiescent galaxies with low specific star formation rates (sSFR $\le 0.3 f_{\text{H}}$) are shown by the solid symbols, while star-forming galaxies (sSFR $\ge 0.3 f_{\text{H}}$) are shown by the open symbols. Quenched galaxies in the AGN feedback model decrease their core density with time since $z \sim 2$, while galaxies simulated without AGNs keep increasing their core density. (Right) Average tracks of simulated galaxies on the $\Sigma_I - M_*$ plane from $z=3$ to $z=0$ shown for two models, with AGN feedback in green and without AGN feedback in purple. The points show the core density–mass relation of the simulated galaxies at $z=3$ with open symbols and at $z=0$ with filled symbols respectively.

gas mass loss and an impulsive expansion consequent to the rapid gas mass loss. The former includes the slow stellar mass loss associated with the death of old stars, i.e., AGB winds. In addition, the slow expulsion of gas during the quiescent phases of AGNs can also adiabatically expand the central region. The latter rapid mass loss involves SN-driven as well as AGN-driven winds (Fan et al. 2008), with a shorter ejection timescale than the dynamical timescale.

If we define the fractional changes in mass and radius of a self-gravitating system as $\delta_m \equiv (m_1 - m_0)/m_0$ and $\delta_r \equiv (r_1 - r_0)/r_0$, where $m_0$ and $m_1$ are the initial and final masses after mass loss, and $r_0$ and $r_1$ are the initial and final radii, we have adiabatic expansion with a rate of $\delta_r = -\delta_m/(\delta_m + 1)$ after mass loss on a timescale longer than a dynamical timescale. But when the gas mass is ejected with a timescale shorter than the dynamical timescale of the system, the expansion proceeds at a higher rate as $\delta_r = -\delta_m/(2\delta_m + 1)$. For example, for a 25% increase in radius, $\delta_r = 0.25$, which leads to a ~50% decrease in density, we require a mass loss of $\delta_m = 1/3$ for the adiabatic case, but only $\delta_m = 1/6$ for the impulsive mass-loss case. Outflowing gas in our galaxies simulated with AGN-driven wind feedback show a much higher characteristic velocity (500–1000 km s⁻¹; Brennan et al. 2018).

Finally, core scouring by binary black hole systems can heat and expel collisionless matter from the central region: the gravitational slingshot effect (Milosavljević & Merritt 2001). The black hole scouring effect is able to produce the observed light profiles of galaxy cores (Faber et al. 1997), and its impact is found to be imprinted on recently observed massive elliptical galaxies (Thomas et al. 2016), which show a tight correlation between the sizes of the core in the

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

**Figure 5.** Stellar radial migration distance since $z=2$ measured for all star particles within the central region $R < 1 \text{ kpc}$ for 30 galaxies simulated with an AGN.

observed light profile and the radii of the black hole's sphere of influence. Our simulation incorporates black hole mergers and naturally includes the heating effect of the corresponding black hole binary orbital decays. However, the resolution in
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of the galaxies simulated with AGN feedback (green circles) and without AGN feedback (purple triangles) on the sSFR-$\Sigma_{\text{eff}}$ ($=M_*/\pi r_{\text{eff}}^2$) plane from $z=3$ to $z=0.5$. (Left) Specific star formation rate vs. projected stellar density measured within the effective radius for simulated galaxies are shown at $z=3, 2, 1$, and $0.5$. As in Figures 1 and 4, we show quiescent galaxies with low specific star formation rates (sSFR $\le 0.3\,h_{\text{0}}^{-1}$) in solid symbols and star-forming galaxies (sSFR $\ge 0.3\,h_{\text{0}}^{-1}$) in open symbols. We indicate the galaxies with specific star formation rates of log sSFR $<-4$ with downward arrows showing their effective stellar densities. Galaxies simulated with AGN feedback start to quench their star formation at $z \sim 2$, and gradually move downwards on the sSFR-$\Sigma_{\text{eff}}$ plane, and move toward the left afterward. (Right) Average tracks of simulated galaxies on the sSFR-$\Sigma_{\text{eff}}$ plane from $z=3$ to $z=0.5$ shown for two models, with AGN feedback in green and without AGN feedback in purple. The points show the location of the simulated galaxies on the sSFR-$\Sigma_{\text{eff}}$ plane at $z=3$ with open symbols and at $z=0$ with filled symbols.</figcaption></figure>

our simulation is not sufficient to capture this effect of black hole scouring in detail, as the black hole binary scouring is considered to be mainly important on length scales of $\lesssim 100$ pc (see Rantala et al. 2017, 2018).

In order to check the contribution of the “puffing-up” process to the stellar core decrease, we show the histogram of the stellar radial migration distance from $z=2$ to $z=0$ for all star particles which were within the central region $r<1$ kpc at $z=2$ in 30 WithAGN galaxies in Figure 5. Although there are a few star particles that travelled inward (radial travel distance $<0$), most of stars have been migrating radially outward since $z=2$. Over 70% of star particles have been migrating radially outward more than 1 kpc since $z=2$, that is, the majority of star particles that constituted the core at $z=2$ are no longer within the central region. This is primarily due to the expansion after gas ejection as discussed above. Furlong et al. (2017) also recently found such a migration of stars in their compact galaxies selected from the EAGLE simulation. They showed that this star migration from the central region of galaxies can result in the size growth of compact galaxies, but its contribution is modest, when renewed star formation and mergers dominate the size growth.

### 3.3. Evolution of Galaxies in the sSFR-$\Sigma_{\text{eff}}$ Plane

In order to study the structural evolution of galaxies and its connection to the quenching of star formation, in Figure 6 we show the evolution of simulated central galaxies in the plane of the sSFR versus the effective stellar density $\Sigma_{\text{eff}}(=M/\pi r_{\text{eff}}^2)$, following Barro et al. (2013). In the left panels of Figure 6, we show the sSFR and $\Sigma_{\text{eff}}$ of the 30 central galaxies simulated

with and without AGN feedback at $z=0.5$, and their most massive progenitors at $z=1, 2$, and 3. In the right panel of Figure 6, we show the average tracks of simulated galaxies in the sSFR-$\Sigma_{\text{eff}}$ plane for the WithAGN and NoAGN models, showing the schematic evolutionary tracks with the arrows. At $z=3$, all progenitor galaxies simulated with and without AGN feedback are distributed horizontally in the sSFR-$\Sigma_{\text{eff}}$ plane, with almost constant sSFR over a broad range in $\Sigma_{\text{eff}}$. By $z=2$, galaxies in both models have moved toward the right, becoming more compact. Note that we see this “compaction” trend even before galaxies start to quench, presumably due to dissipative processes such as gas-rich mergers and disk instabilities (Zolotov et al. 2015; Tacchella et al. 2016). By $z=0.5$, the galaxies simulated with AGN feedback show a rapid decline in sSFR, moving downwards in the diagram rapidly, as the quenching timescale is very short (see Pandya et al. 2017). NoAGN galaxies also show a decline in their sSFR, but the changes are much smaller and slower, and they do not decline below $\log \text{sSFR} = -1.5$. After star formation is quenched, galaxies simulated with AGNs decrease their stellar densities, moving toward the left on the sSFR-$\Sigma_{\text{eff}}$ plane, while galaxies simulated without AGN feedback show a relatively small decrease in their stellar density.

We note that galaxies simulated with AGN feedback evolve along an evolutionary track strikingly similar to the one suggested by Barro et al. (2013, see their Figure 6). They characterized the structural evolution of galaxies in the $\Sigma_{1.5}(=M/r_{\text{eff}}^{1.5})$ plane and showed that galaxies follow three sequences: (1) star-forming and diffuse galaxies transformed into cSFGs via multiple dissipative processes, (2) cSFGs turned into cQGs by rapid quenching of star formation, and finally, (3)
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Center panel: evolutionary tracks of two example galaxies on the plane of sSFR–Σ<sub>eff</sub>. The redshifts z = 3, 2, and 1 are marked along the tracks by open squares. The fixed constant threshold log sSFR = -0.5 adopted in Barro et al. (2013) is shown by the black horizontal dotted line. Four left panels: projected stellar density maps at redshifts z = 3, 2, 1, and z = 0 of the example galaxy with M<sub>*,z=0</sub> = 3.2 × 10<sup>11</sup> M<sub>☉</sub> at z = 0 (shown in blue in the center panel). The snapshots are 100 kpc on a side, and a redder color indicates a higher density. Four right panels: same as left panels but for the example galaxy with M<sub>*,z=0</sub> = 9.3 × 10<sup>10</sup> M<sub>☉</sub> at z = 0 (shown in red in the center panel).</figcaption></figure>

in cQGs, the galaxy sizes are extended via gradual growth and expansion due to multiple dry mergers. The galaxies simulated with AGN feedback show a similar “clockwise” evolutionary track and similar transitions from extended SFGs, to cSFGs, to cQGs, and finally to extended QGs. By contrast, the galaxies in the NoAGN simulations show a much less dramatic evolution.

Figure 7 shows the evolution of two galaxies: we show the evolutionary tracks in sSFR versus Σ<sub>eff</sub> in the center panel and their projected stellar density snapshots at redshifts z = 3, 2, 1, and z = 0. Again, the galaxies follow tracks similar to those shown in Barro et al. (2013). The galaxy sizes become smaller at a roughly constant sSFR, and the onset of quenching again occurs at a nearly constant Σ<sub>eff</sub>.

**4. Summary and Discussion**

We have explored the role of AGN feedback in the evolution of galaxy sizes, compactness, and core densities using numerical cosmological hydrodynamical simulations with an observationally based novel subgrid treatment of winds driven by radiatively efficient accretion onto a supermassive black hole. The results can be summarized as follows.

1. We show that galaxies simulated with AGN feedback follow much steeper evolution with Δr ∝ ΔM² in the size–mass plane compared to galaxies without AGN feedback (Figure 1), and show rapid size evolution with redshift, r<sub>eff</sub> ∝ (1 + z)<sup>-1.6</sup> (Figure 2).

2. AGN feedback enhances the “inside out” evolution of massive galaxies via efficient star formation quenching mechanisms (e.g., Hopkins et al. 2010; van de Sande et al. 2013) by effectively quenching in situ star formation, therefore increasing the fraction of accreted stars in a galaxy (Figure 3).

3. We also show that galaxies simulated with AGN feedback decrease their core densities after the quenching of star formation, while galaxies without AGN feedback continue to increase their core densities as well as stellar mass (Figure 4). The decrease in core density is caused by stellar mass loss as well as the “puffing-up” process in response to gas mass loss (Figure 5).

We also study the evolution of galaxies in the star formation versus compactness (size) plane (Figure 6). We show the role of AGN feedback on the formation of cQGs at z ~ 2 and
present-day extended QGs through the evolutionary scenario suggested by Barro et al. (2013), which can be summarized as follows.

1. The formation of cQGs first follows an evolutionary sequence from extended SFGs (upper left) to cSFGs (upper right) due to a dissipative process, such as highly dissipative mergers between gas-rich progenitors, which are more common at high redshift (e.g., Covington et al. 2011).

2. Then, the AGN feedback quenches the star formation of cSFGs and lowers their sSFR at a roughly constant Σ<sub>eff</sub>, populating the cQG region (lower right) rapidly.

3. At later times, after in situ star formation has been quenched, the simulated galaxies increase their stellar mass primarily by the accretion of smaller stellar systems, which leads to a strong size growth (Naab et al. 2006, 2009; Hilz et al. 2012, 2013; Johansson et al. 2012; Oser et al. 2012; Oogi & Habe 2013; Bédorf & Portegies Zwart 2013, but see Newman et al. 2012; Nipoti et al. 2012). Minor mergers predominantly deposit materials in the outskirts of the galaxies and finally form extended QGs populating the lower-left region of the sSFR−Σ<sub>eff</sub> plane.

In summary, mechanical AGN feedback plays an important role in making massive galaxies red and dead, and it also plays an important role in making them extended. Including mechanical AGN feedback enhances the size growth in different ways. First of all, it has an indirect, supporting, but major impact, as with effective in situ star formation quenching, it increases the fraction of accreted stars (Hirschmann et al. 2012; Dubois et al. 2013), so the impact of dry mergers is much enhanced. Moreover, since completely dry mergers are less frequent due to the high gas fractions in satellite galaxies without AGN feedback, including AGN feedback can increase the frequency of dry mergers. Second, it has a direct but probably minor impact, as AGN-driven gas outflows can produce fluctuations of the gravitational potential and puff up the central region. It will also enhance the adiabatic expansion in response to the slow expulsion of gas from the central region during the quiescent mode of AGN activity.

The authors would like to thank the anonymous referee for helpful comments. We also thank Guillermo Barro, Renyue Cen, Sandra Faber, Shy Genel, and Sarah Wellons for helpful conversations and suggestions. Numerical simulations were run on the computer clusters of the Princeton Institute of Computational Science and engineering. This research was supported by NASA through grant number HST Cycle 23 AR-14287 from the Space Telescope Science Institute.

ORCID iDs

Ena Choi https://orcid.org/0000-0002-8131-6378

Jeremiah P. Ostriker https://orcid.org/0000-0002-6405-9904

Thorsten Naab https://orcid.org/0000-0002-7314-2558

Michaela Hirschmann https://orcid.org/0000-0002-3301-3321

References

Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380, 963

Arav, N., Borguet, B., Chamberlain, C., Edmonds, D., & Danforth, C. 2013, MNRAS, 436, 3286

Aumer, M., White, S. D. M., Naab, T., & Scannapieco, C. 2013, MNRAS, 434, 3142

Barnes, J. E., & Hernquist, L. E. 1991, ApJL, 370, L65

Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2013, ApJ, 765, 104

Barro, G., Faber, S. M., Pérez-González, P. G., et al. 2014, ApJ, 791, 52

Bédorf, J., & Portegies Zwart, S. 2013, MNRAS, 431, 767

Bell, E. F., van der Wel, A., Papovich, C., et al. 2012, ApJ, 753, 167

Bell, E. F., Wolf, C., Meisenheimer, K., et al. 2004, ApJ, 608, 752

Bezanson, R., van Dokkum, P. G., Tal, T., et al. 2009, ApJ, 697, 1290

Biermann, P., & Shapiro, S. L. 1979, ApJL, 230, L33

Bluck, A. F. L., Mendel, J. T., Ellison, S. L., et al. 2014, MNRAS, 441, 599

Bluck, A. F. L., Mendel, J. T., Ellison, S. L., et al. 2016, MNRAS, 462, 2559

Bondi, H. 1952, MNRAS, 112, 195

Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273

Brammer, G. B., Whitaker, K. E., van Dokkum, P. G., et al. 2009, ApJL, 706, L173

Brennan, R., Choi, E., Somerville, R. S., et al. 2018, ApJ, 860, 14

Brennan, R., Pandya, V., Somerville, R. S., et al. 2017, MNRAS, 465, 619

Buitrago, F., Trujillo, I., Conselice, C. J., et al. 2008, ApJL, 687, L61

Cappellari, M. 2013, ApJL, 778, L2

Cassata, P., Giavalisco, M., Guo, Y., et al. 2011, ApJ, 743, 96

Chan, J. C. C., Beifiiori, A., Saglia, R. P., et al. 2018, ApJ, 856, 8

Cheung, E., Faber, S. M., Koo, D. C., et al. 2012, ApJ, 760, 131

Choi, E., Naab, T., Ostriker, J. P., Johansson, P. H., & Moster, B. P. 2014, MNRAS, 442, 440

Choi, E., Ostriker, J. P., Naab, T., et al. 2017, ApJ, 844, 31

Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125

Choi, E., Ostriker, J. P., Naab, T., Oser, L., & Moster, B. P. 2015, MNRAS, 449, 4105

Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21

Ciotti, L., Ostriker, J. P., & Proga, D. 2010, ApJ, 717, 708

Conselice, C. J., Yang, C., & Bluck, A. F. L. 2009, MNRAS, 394, 1956

Covington, M. D., Primack, J. R., Porter, L. A., et al. 2011, MNRAS, 415, 3135

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Crenshaw, D. M., & Kraemer, S. B. 2012, ApJ, 753, 75

Cullen, L., & Dehnen, W. 2010, MNRAS, 408, 669

Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680

Damjanov, I., Abraham, R. G., Glazebrook, K., et al. 2011, ApJL, 739, L44

Damjanov, I., McCarthy, P. J., Abraham, R. G., et al. 2009, ApJ, 695, 101

Dehnen, W., & Aly, H. 2012, MNRAS, 425, 1068

Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870

Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

Diamond-Stanic, A. M., Moustakas, J., Tremonti, C. A., et al. 2012, ApJL, 755, L26

Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS, 433, 3297

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Durier, F., & Dalla Vecchia, C. 2012, MNRAS, 419, 465

Faber, S. M., Tremaine, S., Ajhar, E. A., et al. 1997, AJ, 114, 1771

Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265

Fabian, A. C. 1994, ARA&A, 32, 277

Fan, L., Lapi, A., De Zotti, G., & Danese, L. 2008, ApJL, 689, L101

Fang, J. J., Faber, S. M., Koo, D. C., & Dekel, A. 2013, ApJ, 776, 63

Ferguson, H. C., Dickinson, M., Giavalisco, M., et al. 2004, ApJL, 600, L107

Franx, M., van Dokkum, P. G., Schreiber, N. M. F., et al. 2008, ApJ, 688, 770

Furlong, M., Bower, R. G., Crain, R. A., et al. 2017, MNRAS, 465, 722

Gammie, C. F. 2001, ApJ, 553, 174

Genel, S., Nelson, D., Pillepich, A., et al. 2018, MNRAS, 474, 3976

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175

Gibson, R. R., Jiang, L., Brandt, W. N., et al. 2009, ApJ, 692, 758

Haardt, F., & Madau, P. 2001, in Proc. XXIst Morion Astrophysics Meeting, XXVth Reconcres de Moriond, ed. D. M. Neumann & J. T. T. Van (Paris: ARISF), 64

Hills, J. G. 1980, ApJ, 235, 986

Hilz, M., Naab, T., Ostriker, J. P., et al. 2012, MNRAS, 425, 3119

Hilz, M., Naab, T., & Ostriker, J. P. 2013, MNRAS, 429, 2924

Hirschmann, M., Naab, T., Davé, R., et al. 2013, MNRAS, 436, 2929

Hirschmann, M., Naab, T., Somerville, R. S., Burkert, A., & Oser, L. 2012, MNRAS, 419, 3200

Hopkins, P. F. 2013, MNRAS, 428, 2840

Hopkins, P. F., Bundy, K., Hernquist, L., Wuyts, S., & Cox, T. J. 2010, MNRAS, 401, 1099

Hopkins, P. F., Hernquist, L., Cox, T. J., Robertson, B., & Springel, V. 2006, ApJS, 163, 50

Hoyle, F., & Lyttleton, R. A. 1939, PCPS, 34, 405
Hu, C.-Y., Naab, T., Walch, S., Moster, B. P., & Oser, L. 2014, MNRAS, 443, 1173

Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439

Janka, H.-T. 2012, ARNPS, 62, 407

Johansson, P. H., Naab, T., & Ostriker, J. P. 2012, ApJ, 754, 115

Karakas, A. I. 2010, MNRAS, 403, 1413

Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54

Kroupa, P. 2001, MNRAS, 322, 231

Lackner, C. N., Cen, R., Ostriker, J. P., & Joung, M. R. 2012, MNRAS, 425, 641

Leitner, S. N., & Kravtsov, A. V. 2011, ApJ, 734, 48

Longhetti, M., Saracco, P., Severgnini, P., et al. 2007, MNRAS, 374, 614

López-Sanjuan, C., Le Fèvre, O., Tasca, L. A. M., et al. 2013, A&A, 553, A78

Lotz, J. M., Jonsson, P., Cox, T. J., et al. 2011, ApJ, 742, 103

Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

Martizzi, D., Teyssier, R., & Moore, B. 2013, MNRAS, 432, 1947

Martizzi, D., Teyssier, R., Moore, B., & Wentz, T. 2012, MNRAS, 422, 3081

Milosavljević, M., & Merritt, D. 2001, ApJ, 563, 34

Moe, M., Arav, N., Bautista, M. A., & Korista, K. T. 2009, ApJ, 706, 525

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJL, 699, L178

Naab, T., Khochfar, S., & Burkert, A. 2006, ApJL, 636, L81

Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59

Newman, A. B., Ellis, R. S., Bundy, K., & Treu, T. 2012, ApJ, 746, 162

Newman, A. B., Ellis, R. S., Treu, T., & Bundy, K. 2010, ApJL, 717, L103

Nipoti, C., Treu, T., Auger, M. W., & Bolton, A. S. 2009, ApJL, 706, L86

Nipoti, C., Treu, T., Leauthaud, A., et al. 2012, MNRAS, 422, 1714

Núñez, A., Ostriker, J. P., Naab, T., et al. 2017, ApJ, 836, 204

Omand, C. M. B., Balogh, M. L., & Poggianti, B. M. 2014, MNRAS, 440, 843

Oogi, T., & Habe, A. 2013, MNRAS, 428, 641

Oser, L., Naab, T., Ostriker, J. P., & Johansson, P. H. 2012, ApJ, 744, 63

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

Pandya, V., Brennan, R., Somerville, R. S., et al. 2017, MNRAS, 472, 2054

Patel, S. G., van Dokkum, P. G., Franx, M., et al. 2013, ApJ, 766, 15

Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464

Porter, L. A., Somerville, R. S., Primack, J. R., & Johansson, P. H. 2014, MNRAS, 444, 942

Proga, D., & Kallman, T. R. 2004, ApJ, 616, 688

Rantala, A., Johansson, P. H., Naab, T., Thomas, J., & Frigo, M. 2018, ApJ, 864, 113

Rantala, A., Pihajoki, P., Johansson, P. H., et al. 2017, ApJ, 840, 53

Read, J. I., & Hayfield, T. 2012, MNRAS, 422, 3037

Ritchie, B. W., & Thomas, P. A. 2001, MNRAS, 323, 743

Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371

Saitoh, T. R., & Makino, J. 2009, ApJL, 697, L99

Saitoh, T. R., & Makino, J. 2013, ApJ, 768, 44

Sazonov, S. Y., Ostriker, J. P., Ciotti, L., & Sunyaev, R. A. 2005, MNRAS, 358, 168

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

Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Spergel, D. N., Bean, R., Doré, O., et al. 2007, ApJS, 170, 377

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776

Stewart, K. R., Bullock, J. S., Barton, E. J., & Wechsler, R. H. 2009, ApJ, 702, 1005

Stone, N. C., Küpper, A. H. W., & Ostriker, J. P. 2017, MNRAS, 467, 4180
Strömgren, B. 1939, ApJ, 89, 526

Szomoru, D., Franx, M., van Dokkum, P. G., et al. 2010, ApJL, 714, L244

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242

Taylor, E. N., Franx, M., Glazebrook, K., et al. 2010, ApJ, 720, 723

Thomas, J., Ma, C.-P., McConnell, N. J., et al. 2016, Natur, 532, 340

Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285

Trujillo, I., Conselice, C. J., Bundy, K., et al. 2007, MNRAS, 382, 109

Trujillo, I., Rudnick, G., Rix, H.-W., et al. 2004, ApJ, 604, 521

van de Sande, J., Kriek, M., Franx, M., et al. 2013, ApJ, 771, 85

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

van Dokkum, P. G., Bezanson, R., van der Wel, A., et al. 2014, ApJ, 791, 45

van Dokkum, P. G., Franx, M., Kriek, M., et al. 2008, ApJL, 677, L5

van Dokkum, P. G., Nelson, E. J., Franx, M., et al. 2015, ApJ, 813, 23

van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018

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

Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291

Wellons, S., Torrey, P., Ma, C.-P., et al. 2015, MNRAS, 449, 361

Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99

Williams, R. J., Quadri, R. F., Franx, M., et al. 2010, ApJ, 713, 738

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Wuyts, S., Förster Schreiber, N. M., Lutz, D., et al. 2011, ApJ, 738, 106

Zołotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327