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

# Physics of Galactic Metals: Evolutionary Effects due to Production, Distribution, Feedback, and Interaction with Black Holes

Ena Choi¹, Jeremiah P. Ostriker²,³, Thorsten Naab⁴, Rachel S. Somerville¹,⁵,

Michaela Hirschmann⁶, Alejandro Núñez², Chia-Yu Hu⁴,⁵, and Ludwig Oser⁴

¹Department of Physics and Astronomy, Rutgers, The State University of New Jersey, New Brunswick, NJ 08854, USA

²Department of Astronomy, Columbia University, New York, NY 10027, USA

³Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

⁴Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85741 Garching, Germany

⁵Simons Center for Computational Astrophysics, New York, NY 10010, USA

⁶Sorbonne Université's, UPMC-CNRS, UMR7095, Institut d' Astrophysique de Paris, F-75014 Paris, France

Received 2016 October 28; revised 2017 May 11; accepted 2017 June 6; published 2017 July 19

## Abstract

We ask how the inclusion of various physical heating processes due to the metal content of gas affects the evolution of central massive galaxies and compute a suite of cosmological hydrodynamical simulations that follow these systems and their supermassive black holes. We use a smoothed particle hydrodynamics code with a pressure-entropy formulation and a more accurate treatment of the metal production, turbulent diffusion, and cooling rate based on individual element abundances. The feedback models include (1) active galactic nucleus (AGN) feedback via high-velocity broad absorption line winds and Compton/photoionization heating; (2) stellar feedback from multiple processes, including powerful winds from supernovae, stellar winds from young massive stars, and AGB stars, as well as radiative heating within Strömgren spheres; and (3) additional heating effects due to the presence of metals, including grain photoelectric heating and metallicity-dependent X-ray heating by nearby accreting black holes and from the cosmic X-ray background. With a suite of zoom-in simulations of 30 halos with $M_{vir} \sim 10^{12.0-13.4}$, we show that energy and momentum budgeted from all feedback effects generate realistic galaxy properties. We explore the detailed role of each feedback model with three additional sets of simulations with varying input physics. We show that the metal-induced heating reduces the fraction of accreted stellar material but overall has a relatively minor effect on the massive central galaxies. The inclusion of AGN feedback significantly improves the ability of our simulations to yield realistic gas and stellar properties of massive galaxies with a reasonable accreted star fraction from other galaxies.

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

## 1. Introduction

The physical treatment of galaxy formation and evolution has steadily improved over recent years as the cosmological setting within which this process occurs becomes more definite and confirmed by multiple observations. Thus, the gravitational accumulation of mass on galactic scales from well-defined cosmological perturbations is well understood. Recent progress has focused on the “feedback,” how the stars (e.g., Springel & Hernquist 2003; Oppenheimer & Davé 2006; Hopkins et al. 2014; Agertz & Kravtsov 2015) and black holes (e.g., Springel et al. 2005a; Sijacki et al. 2007; Booth & Schaye 2009; Dubois et al. 2012) formed from the infalling gas can initiate processes that in turn alter that infall, drive outflows, and change the rate of star formation from what would be guessed at from the simple picture of cosmological infall (see Somerville & Davé 2015, for a literature review). The two main drivers are supernovae (SNe) and active galactic nuclei (AGNs), but even normal star formation of massive and lower-mass stars can have significant effects (Agertz et al. 2013; Renaud et al. 2013; Stinson et al. 2013; Hopkins et al. 2014).

The first generation of physics treatments for massive galaxies allowing for the incorporation of feedback concentrated on the mass, energy, and momentum input by massive stars (e.g., Oppenheimer & Davé 2006; Dalla Vecchia & Schaye 2008) and black holes (e.g., Debuhr et al. 2010; Ostriker et al. 2010; Choi et al. 2012) found that the effects were dramatic. The ratio of baryons incorporated to those

available was reduced by a factor of over 10, while the size of moderate-mass galaxies and the duration of galaxy formation for these systems increased dramatically as kinetic feedback expelled gas. Some of this material returns at later times with increased angular momentum (Naab et al. 2014; Genel et al. 2015). For those galaxies and for more massive ones the production of a circumgalactic medium due to outflow was found to have dramatic physical and observational consequences (e.g., Ford et al. 2013; Faucher-Giguère et al. 2015; Oppenheimer et al. 2016).

Recently the feedback effects of the chemical evolution have been included, in addition to the mechanical effects. While most “metals” (everything heavier than H and He) and particularly α-process elements such as O and Mg are produced by young, massive stars, some elements, such as Fe, are known to be produced primarily (but not exclusively) by older, lower-mass stars exploding as Type I SNe, and some elements, e.g., N, are significantly contributed by normal evolution of solar-type stars. The line emission from these elements greatly increases cooling rates as compared to an (H, He) plasma, and recent papers such as Finlator & Davé (2008), Oppenheimer & Davé (2008), Wiersma et al. (2009a, 2009b), and Aumer et al. (2013) highlighted the dramatic consequences for the evolution of normal galaxies, with the emphasis on moderate-mass spiral systems. The net effect of the extra cooling processes is to again increase star formation rates somewhat, mitigating the mechanical effects of feedback. However, allowing only for the
<table><caption>**Table 1**
Synoptic Table of Galaxy Model Development</caption><thead><tr><th rowspan="2">Input Physics</th><th rowspan="2">Oser et al. (2010, 2012)</th><th rowspan="2">Choi et al. (2015)<br>ThSNnoMetal</th><th colspan="3">This Work</th></tr><tr><th>Fiducial</th><th>NoAGN</th><th>NoZHeating</th></tr></thead><tbody><tr><td>Improved SPH<sup>a</sup></td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Stellar feedback</td><td>Thermal<sup>b</sup></td><td>Thermal<sup>b</sup></td><td>Mechanical<sup>c</sup></td><td>Mechanical<sup>c</sup></td><td>Mechanical<sup>c</sup></td></tr><tr><td>AGN feedback</td><td>No</td><td>Mechanical<sup>d</sup></td><td>Mechanical<sup>d</sup></td><td>No</td><td>Mechanical<sup>d</sup></td></tr><tr><td>Metal production and cooling</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Metal heating<sup>e</sup></td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>No</td></tr></tbody></table>

**Notes.** A comparison of present and past treatment.

<sup>a</sup> Our improved SPH scheme includes density-independent pressure-entropy SPH formulation, improved artificial viscosity and thermal conductivity, and a time-step limiter for shock ambient particles (Hu et al. 2014).

<sup>b</sup> Thermal SN feedback model from Springel & Hernquist (2003).

<sup>c</sup> Our mechanical stellar feedback model includes various processes such as “snowplow” SN winds, stellar winds from young massive stars, and AGB stars (Núñez et al. 2017). See Section 2.3 for details.

<sup>d</sup> Our mechanical AGN feedback model includes high-velocity (10,000 km s<sup>-1</sup>) broad absorption line winds, as well as photoionization and Compton heating and associated radiation pressure from AGNs (Choi et al. 2015). See Section 2.5 for details.

<sup>e</sup> Metal-induced heating prescription includes the photoelectric heating and enhanced Compton heating by X-ray background. See Section 2.6 for details.

cooling effects of heavy elements is too one-sided and thus not fully realistic.

There are also substantial heating effects due to the presence of metals (Draine 2010), and this paper attempts to spell out the most important ones. As we shall see, this reverses some of the cooling effects of metals recently included as we proceed, by successive approximations, toward a physically based and hopefully accurate treatment of galaxy formation. With successive iterations the situation becomes more complex, so, for example, the metal component of interstellar gas is heated by X-rays both from the black holes within the galaxy and from the extragalactic background, and photons from massive stars can heat gas via the photoelectric effect on metals distributed in the interstellar medium (ISM). All of these processes tend to increase the heating rate, causing a reduction of star formation, i.e., causing “negative feedback.”

In this paper, we describe the implementation of the metal-induced heating mechanisms, as well as our more standard treatments of star formation and other physical processes. We present in some detail the black-hole-driven AGN feedback, as it is different from that implemented in most codes in that it allows for broad absorption line winds and X-ray and UV output of both energy and momentum, rather than explicitly postulating “radio” feedback (e.g., Vogelsberger et al. 2013) or turning off cooling (e.g., Tremmel et al. 2016). In addition, our SN feedback allows for three stages of SN remnant propagation (ejecta dominated, Sedov-Taylor, and snowplow) with a relatively realistic treatment (see Núñez et al. 2017, for details), without kinematically decoupling the outflow particles or disabling cooling for some fixed time interval as in many of the SN feedback implementations in the literature.

The purpose of this paper is to study the respective role of various feedback models implemented in our simulation on the physical properties of central massive galaxies. In Section 3, we present the 30 zoom-in simulations of halos with $M_{vir} \sim 10^{12-13.4}$ with four different feedback models outlined in Section 2, pointing out the changes induced by including the various physical processes. Finally, Section 4 summarizes our conclusions.

## 2. The Simulations

The galaxy formation model and input physics adopted in the simulations have been improved from the previous papers. In this section we describe the new ingredients of the current model, and we refer to the previous papers for the unchanged part. A comparison of present and past treatment is summarized in Table 1.

### 2.1. Hydrodynamics Code

We use a modified version of the parallel smoothed particle hydrodynamics (SPH) code GADGET-3 (Springel 2005). In order to avoid the numerical artifacts of the classical SPH code (e.g., Agertz et al. 2007), we use SPHGal (Hu et al. 2014), a modified version of GADGET-3 that includes a density-independent pressure-entropy SPH formulation (Ritchie & Thomas 2001; Hopkins 2013; Saitoh & Makino 2013). To further improve over standard SPH, we adopt the Wendland $C^4$ kernel with 200 neighboring particles (Dehnen & Aly 2012), the improved artificial viscosity implementation (Cullen & Dehnen 2010), and the artificial thermal conductivity following Read & Hayfield (2012). Lastly, we employ a time-step limiter to ensure that neighboring particles have a similar time step, such that ambient particles do not remain inactive when a shock travels through (Saitoh & Makino 2009; Durier & Dalla Vecchia 2012). We refer the readers to Hu et al. (2014) for the performance of the new SPH schemes in the test problems.

### 2.2. Star Formation and Chemical Enrichment

We adopted the star formation and chemical evolution model described in Aumer et al. (2013), which allows chemical enrichment by winds driven by Type I 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) respectively. We explicitly trace the mass in 11 different species, H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe, for both star and gas particles. Then the net cooling rate is calculated based on individual element abundances, temperature, and density of gas (see Aumer et al. 2013, for details). We adopted the cooling rate from Wiersma et al. (2009a) for optically thin gas in ionization equilibrium. We
also include a redshift-dependent UV/X-ray and cosmic microwave background with a modified Haardt & Madau (2001) spectrum.

We include a model for turbulent diffusion of gas-phase metals. Following Aumer et al. (2013), we allow the metal-enriched gas particles to mix their metals with neighboring gas particles using the standard SPH neighbor searches.

We stochastically form stars if the gas density exceeds a density threshold, which we defined as $n_{th} \equiv n_0(T_{gas}/T_0)^3(M_0/M_{gas})^2$, where critical threshold density and temperature are $n_0 = 2.0 \text{ cm}^{-3}$ and $T_0 = 12,000$ K, respectively, and $M_0$ is the gas particle mass in fiducial resolution. We require that the gas density should be higher than the value for the Jeans gravitational instability of a mass $M_{gas}$ at temperature $T_{gas}$. The star formation rate is calculated 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 for gas particles, respectively. The star formation efficiency $\eta$ is set as $\eta = 0.025$. In our model, star formation is feedback regulated and star-forming regions grow in mass until sufficient new stars have formed to stop further collapse.

## 2.3. Early Stellar Feedback and Supernova Feedback

We use the early stellar and SN feedback model adopted in Núñez et al. (2017; see Section 2 in their paper for further details). We include the winds from young stars, UV heating from young stars, Strömgren spheres, three-phase SN remnant input from both Type I and Type II SN feedback, and outflow and metals from dying low-mass AGB stars. Each of the processes is included as explicit physical processes communicating mass, metals, momentum, and energy from stellar particles to surrounding gaseous particles.

In the early stellar feedback model, we include the effect of stellar winds, as well as heating by the ionizing radiation from young massive stars before they explode as SNe. We distribute the momentum of winds from massive stars to the closest gas particles with the same amount of ejected mass and momentum as those of Type II SN explosions evenly spread in time before the moment of SN explosion $t_{SN} = 3$ Myr. We also add the effects of ionizing radiation from massive stars. The cold gas with $T < 10^4$ K within a Strömgren (1939) radius from each star particle is gradually heated to $T = 10^4$ K and is not allowed to cool below this temperature until it is no longer in an H II region, within the Strömgren radius (see also Hopkins et al. 2012).

In the SN feedback model, a single SN event is assumed to eject mass in an outflow with a velocity $v_{out,SN} = 4500 \text{ km s}^{-1}$, a typical velocity of outflowing materials in SNe. We distribute SN energy and momentum to the surrounding ISM from the SN event. We assume the SN outflows to be one of three characteristic SN remnant phases that determine how much SN energy is transferred to gas particles as kinetic and thermal energy. Depending on the distance from the SN events, each neighboring gas particle is affected by one of the three successive phases: (i) the momentum-conserving free expansion phase; (ii) the 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 significant. In this model, the SN remnant initiates standard Sedov–Taylor blast waves carrying energy as 30% kinetic and 70% thermal, and both amounts dissipate with distance from the SN, as described by the pressure-driven snowplow phase of SN remnants. See Appendix A in Núñez

et al. (2017) for a detailed description of the implementation of the SN feedback model. This model has provisions for treating the ISM as multiphase with most of the volume hot but most of the mass in warm or cool phases. The volume fraction in the hot phase is estimated by the volume-weighted average temperature following the detailed high-resolution simulations of Li et al. (2015).

Lastly, we include the feedback from the stars of low and intermediate initial mass via slow winds during an AGB phase. We transfer energy and momentum from old star particles to the neighboring gas particles in a momentum-conserving way. The outflowing wind velocity of AGB stars is assumed to be $v_{out,AGB} = 10 \text{ km s}^{-1}$, corresponding to typical outflowing velocities of AGB stars (e.g., Nyman et al. 1992).

The metal-enriched gas from dying stars is returned continuously to the ISM via winds from young stars, SNe, and AGB stars. Over 30% of total mass in stars is ultimately ejected via winds in our initial mass function assumption (Kroupa 2001) and is able to induce the late star formation and quasi-stellar object (QSO) activity by feeding the central supermassive black holes (see also Ciotti et al. 2010).

## 2.4. Black Hole Formation and Growth

In the simulations, the black holes are treated as collisionless sink particles and are seeded in newly forming dark matter halos. The dark matter halos are identified on the fly during a simulation by a friends-of-friends algorithm. The new black holes are seeded with mass of $10^5 \text{ h}^{-1} M_{\odot}$ such that any halo above $1 \times 10^{11} \text{ h}^{-1} M_{\odot}$ contains one black hole at its center if it does not already have a black hole. The dark matter halo threshold mass and black hole seed mass are set to roughly follow the Magorrian et al. (1998) relation and the theoretical calculations of Stone et al. (2016). The chosen seed mass makes a negligible contribution to the final black hole mass.

The black hole mass can grow via two channels: mergers with other black holes, and accretion of gas. We allow the mergers between two black hole particles only when they fall within their local SPH smoothing lengths and their relative velocities are smaller than the local sound speed. As we cannot directly resolve the accretion disk of the black holes on subparsec scales in the cosmological galaxy group-scale simulation, we estimate the rate of the gas infall onto the black hole with a Bondi–Hoyle–Lyttleton parameterization (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952). Following Springel et al. (2005b), the gas accretion rate onto the central region around the black hole is estimated as

$$ \dot{M}_{\text{inf}} = \frac{4\pi G^2 M_{\text{BH}}^2 \rho}{(c_s^2 + v^2)^{3/2}}, \qquad (1) $$

where $\rho$, $c_s$, and $v$ denote the density, the sound speed, and the velocity of the gas relative to the black hole, respectively. Several works adopting Bondi accretion in cosmological simulations often boost the accretion rate by a factor of $\alpha \sim 100$ (see Booth & Schaye 2009, for a literature review); here we do not employ an additional “boost” factor $\alpha$ with regard to the accretion rate.

We also include the soft Bondi criterion introduced in Choi et al. (2012) to avoid the unphysical accretion of unbound gas from outside the Bondi radius of the black hole. This criterion statistically limits the accretion to the gas within the Bondi radius. It also accounts for the size of the gas particle as the
physical properties of each gas particle are smoothed within the kernel size in the smoothed particle hydrodynamics simulations. The full accretion is only allowed when the total volume of a gas particle is included within the Bondi radius. If a gas particle volume is partially included within the Bondi radius, its probability of being absorbed by the black hole is reduced. Finally, we include the free-fall timescale in the accretion rate in order to account for the time that it takes a particle to be accreted to a black hole (Choi et al. 2012).

## 2.5. Feedback from Black Holes

We use mechanical and radiative AGN feedback models described in Choi et al. (2012, 2014). Our AGN feedback model consists of two main components: (1) mechanical feedback as in the broad absorption line winds, which carry energy, mass, and momentum into the surrounding gas; and (2) radiative feedback via the Compton and photoionization heating from the X-ray radiation from the accreting black hole, the radiation pressure associated with the heating, and the Eddington force. The emergent AGN spectrum and metal-line heating are taken from Sazonov et al. (2004). In this section we summarize the main aspects of the AGN feedback model; for more details we refer the reader to Choi et al. (2012).

### 2.5.1. Mechanical AGN Feedback

In the presence of significant AGN winds, only a small fraction of the gas mass inflowing to the central region ultimately accretes to the black hole. AGN winds carry a mass with the outflowing rate given as

$$ \dot{M}_{\text{outf}} = \dot{M}_{\text{inf}} - \dot{M}_{\text{acc}}, \qquad (2) $$

where $\dot{M}_{\text{outf}}$, $\dot{M}_{\text{inf}}$, and $\dot{M}_{\text{acc}}$ denote the outflowing and inflowing mass rate and the mass accretion rate to the black hole, respectively. We assume a constant velocity for AGN wind $v_{\text{outf,AGN}} = 10,000 \text{ km s}^{-1}$, corresponding to a typical broad absorption line wind velocity (e.g., Crenshaw et al. 2003). The momentum flux carried by the wind is $\dot{p} = \dot{M}_{\text{outf}} v_{\text{outf,AGN}}$, and the kinetic energy rate carried by the outflow will be given as

$$ \dot{E}_w \equiv \epsilon_w \dot{M}_{\text{acc}} c^2, \qquad (3a) $$

$$ = \frac{1}{2} \dot{M}_{\text{outf}} v_{\text{outf,AGN}}^2, \qquad (3b) $$

where $\epsilon_w$ denotes the feedback efficiency and we assume $\epsilon_w = 0.005$. The dimensionless quantity $\psi$ is defined as the ratio of the outflow rate to the accreted rate and is given as

$$ \psi \equiv 2 \epsilon_w c^2 / v_{\text{outf,AGN}}^2 = \dot{M}_{\text{outf}} / \dot{M}_{\text{acc}}. \qquad (4) $$

Thus, the black hole accretion rate depends on the mass inflow rate and the dimensionless parameter $\psi$, which is determined by mass and energy conservation as

$$ \dot{M}_{\text{acc}} = \dot{M}_{\text{inf}} \frac{1}{1 + \psi}. \qquad (5) $$

We have $\psi \sim 9$ with our choice of the feedback efficiency $\epsilon_w = 0.005$, and the wind velocity $v_{\text{outf,AGN}} = 10,000 \text{ km s}^{-1}$; thus, the bulk of mass entering the central region $\dot{M}_{\text{inf}}$ is expelled, 10% ($f_{\text{acc}} = \dot{M}_{\text{acc}} / \dot{M}_{\text{inf}} = 1/(1+\psi)$) of the inflowing mass is accreted onto the black hole, and 90%

($f_{\text{outf}} = \dot{M}_{\text{outf}} / \dot{M}_{\text{inf}} = \psi / (1 + \psi)$) is expelled in an outflowing wind.

Among all gas particles entering the central region, the wind particles are stochastically selected with $f_{\text{outf}} = \psi / (1 + \psi)$. We deposit the wind mass and momentum by giving kicks to the selected wind particles. The direction of the wind is set to be parallel or antiparallel to the direction of angular momentum of each gas particle, so that it results in a wind 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 momentum and energy with two nearby gas particles and reproduces the shock-heated momentum-driven flows. The residual energy increases the temperature of the impacted gas particles; therefore, the total energy and momentum are conserved. This prescription gives a ratio of kinetic to thermal energy in the outflowing particles similar to that in the standard Sedov–Taylor blast wave.

### 2.5.2. Radiative AGN Feedback

We also include the heating via X-ray radiation from the accreting black hole. Strong X-ray radiation can be coupled to the surrounding gas according to an approximation described in Sazonov et al. (2005). The net luminosity flux from all black holes in the simulated zoom-in area is calculated at the position of each gas particle, and the flux is converted to the net volume heating rate $\dot{E}$ by adopting the formulae that include Compton and photoionization heating from Sazonov et al. (2005). Finally, the radiation pressure from the X-ray flux absorbed is included. Every fluid element absorbing energy $\Delta E$ from quasar radiation is given an extra momentum $\Delta p = \Delta E/c$ directed away from the black hole.

In many galaxy simulations with black holes, the gas accretion rate onto the black hole is limited to the Eddington rate. In our simulations, we do not limit the accretion rate onto the black holes; instead, we include the Eddington force acting on electrons in the neighboring gas through the hydrodynamic equations, directed radially away from black holes. In this way we allow that super-Eddington gas accretion occasionally occurs, so that the corresponding feedback effect naturally reduces the inflow and increases the outflow.

### 2.5.3. Main Aspects of the AGN Feedback Model

The two AGN feedback mechanisms we have included originate from different parts of a quasar’s SED, i.e., mechanical feedback from UV radiation, and radiative feedback from X-ray radiation. Many previous high-resolution hydrodynamical simulations (e.g., Proga et al. 2000) showed that the flux from strong UV radiation in the AGN spectra dominates for the momentum-driven winds. The region where the winds are generated is very close to the accreting black holes; thus, it is impossible to resolve this scale in cosmological simulations. Therefore, we resort to the subgrid modeling of “mechanical feedback” via the broad absorption line winds that have been accelerated by metal-line trapping (Proga & Kallman 2004). Meanwhile, the heating and associated pressure from AGN radiation are dominated by the moderately hard X-ray region (~50 keV), which is nearly independent of obscuration (Sazonov et al. 2004). We therefore include this feedback effect following a standard atomic physics treatment.

In Choi et al. (2015), we found that AGN feedback has a very strong effect on the star formation and X-ray properties of
galaxies and what matters most is the kind of feedback one includes. The mechanical AGN feedback via broad absorption line winds has more dramatic effects than the traditional thermal feedback treatment in many papers (e.g., Di Matteo et al. 2005; Springel et al. 2005b; Johansson et al. 2009a; Steinborn et al. 2015). The two treatments—thermal versus mechanical—put in the same total energy for a given accretion rate and given efficiency, but putting some fraction of the energy into mechanical rather than thermal increases the effectiveness in driving gas out of the galaxy, as discussed in many recent papers (e.g., Gaspari et al. 2012; Simpson et al. 2015; Barai et al. 2016; Weinberger et al. 2017).

Choi et al. (2014) showed that the mechanical and radiative AGN feedback can drive strong nuclear outflows as observed in many luminous quasars (e.g., Greene et al. 2011; Zakamska et al. 2016), which can further progress on a galaxy-wide scale removing a significant fraction of gas from the host galaxies. Moreover, the relic AGN-driven winds from the height of recent activity found in the simulated galaxies provide a promising explanation for the moderate-velocity outflows (500–1000 km s⁻¹) observed in post-starburst galaxies at intermediate redshift (e.g., Tremonti et al. 2007) and in quiescent nearby galaxies (e.g., Cheung et al. 2016). The outflowing wind properties of our galaxy samples will be further discussed in a forthcoming paper (R. Brennan et al. 2016, in preparation).

## 2.6. Metallicity-dependent Heating Effect

### 2.6.1. Photoelectric Heating

A dust grain can absorb an energetic far-UV photon, excite it to higher energy, and reemit as a “photoelectron” (Spitzer 1948). This photoelectric emission from dust grains dominates the heating of the diffuse ISM and thus can have a significant effect on dynamical evolution of the ISM, as well as on the subsequent star formation (Draine 1978; Weingartner & Draine 2001). The precise modeling of photoelectric heating may require the estimation of the radial distance dependency from OB stars (e.g., Wolfire et al. 2003; Tasker 2011), self-shielding and shielding of the radiation field (Walch et al. 2015; Hu et al. 2016), and a detailed model of dust physical properties such as grain sizes, compositions, and charge state (Bakes & Tielens 1994). However, it would be beyond the scope of this paper to fully account for all these factors that may affect the photoelectric heating rate. We therefore take the following simple prescription of the photoelectric heating rate that is only dependent on the global star formation rate and the metallicity of gas that is related to the abundance of the smallest dust particles, which may primarily consist of polycyclic aromatic hydrocarbons.

From the numerical coefficient of the photoelectric heating rate for the Galaxy from Draine (2010), we added a temperature cutoff $T_{cut} = 10^4$ K so that only cool or cold gas is heated. Then, since the heating rate is proportional to the OB stars, i.e., the star formation rate, we added a factor for (SFR/SFR$_{MW}$). Then the factor (SFR/SFR$_{MW}$) is estimated by

$$ (\text{SFR}/\text{SFR}_{\text{MW}}) = \frac{Z H(t)}{Z_{\odot} H_0}, \qquad (6) $$

as the average value of this factor is just proportional to the metals produced locally and inversely proportional to the time over which they were produced. The photoelectric heating rate

is thus estimated by

$$ \Gamma = 1.4 \times 10^{-26} e^{-T/10^4} \frac{Z H(t)}{Z_{\odot} H_0} \text{ erg s}^{-1} \text{ per hydrogen atom.} \qquad (7) $$

By construction it matches the local photoelectric heating rate given by Draine (2010) for the solar neighborhood.

### 2.6.2. Cosmic X-Ray Background Heating

Energetic X-ray radiation from accreting black holes with a long mean free path can make a significant contribution to the heating of the early universe, as it can easily escape the galaxies and impact the intergalactic medium at long range. We self-consistently include the heating effect from X-ray radiation of individual black holes within our simulation volume. However, as we use zoom-in calculation with the limited volume, it is necessary to account for the effect of X-ray background, which may be able to pre-ionize large volumes of the intergalactic medium unreachable by UV sources. Thus, we include the effects of the low-temperature gas heating by the cosmic X-ray background from accreting black holes in background galaxies as follows.

First, we compute the cosmic X-ray background flux based on the information of the cosmic black hole mass accretion rates derived from the bolometric AGN luminosity functions. The comoving number density of luminous quasars reaches a maximum at $z_{max} \sim 2$ (e.g., Schmidt et al. 1995). We take the global black hole growth rates $\frac{d\rho_{BH}}{dt}$ derived from the bolometric AGN luminosity functions of Shankar et al. (2009) as

$$ \frac{d\rho_{\mathrm{BH}}(z)}{dt} = \frac{1 - \epsilon}{\epsilon c^2} \int_{0}^{\infty} \Phi_{L}(L) L d \log L, \qquad (8) $$

for example, $d\rho_{BH}/dt(z=0) = 3 \times 10^{-6} M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$. Then, we take the cosmic evolution of black hole mass growth rates from Madau & Dickinson (2014) and calculate the emissivity of AGNs in a unit comoving volume of the universe at redshift z as

$$ \varepsilon(z) = \epsilon \frac{d\rho_{\mathrm{BH}}(z)}{dt} c^2, \qquad (9) $$

where the radiative efficiency $\epsilon = 0.1$. For example, $\varepsilon(z=0) = 1.7 \times 10^{40} \text{ erg s}^{-1} \text{Mpc}^{-3}$. The bolometric flux from the cumulative background AGN light at redshift z is then

$$ F_{\mathrm{bk}}(z) = (1+z)^2 \frac{c}{4\pi H_0} \int_{z}^{\infty} \frac{\varepsilon(z') dz'}{(1+z') [\Omega_M (1+z')^3 + \Omega_{\Lambda}]^{1/2}}. \qquad (10) $$

The first $(1+z)^2$ term represents how physical area scales with redshift. We perform the integration of Equation (10) numerically and calculate the metagalactic background flux level as a function of redshift. Then the net flux at the position of each gas particle is calculated as a sum of the global flux from the background AGN at a given redshift and the flux from all black holes in the simulated zoom-in area based on their accretion rate. All gas particles in the simulation are heated accordingly following the heating rate described in Sazonov et al. (2005).
### 2.6.3. Metallicity-dependent X-Ray Heating

For cold gas exposed to the X-rays from an AGN the gas absorbs a certain number of X-ray photons that is proportional to the flux (see Sazonov et al. 2004, for details). Most of this absorption is due to iron and similar heavy elements, and the ejected electrons cause further ionizations and deposit much of the energy as heat. Thus, the heating rate per H atom will be proportional to the X-ray flux, as well as the metallicity. Thus, in order to allow for metal-line absorption, we include an additional metallicity-dependent factor ($1 + 10Z/Z_\odot$) in the X-ray heating rate by AGN radiation over and above Compton heating from Sazonov et al. (2005) as

$$S_1 = 4.1 \times 10^{-35}(1.9 \times 10^7 - T)\xi(1 + 10Z/Z_\odot), \quad (11)$$

where $\xi$ denotes the ionization parameter. Then the volume heating rate $\dot{E}$ by X-ray heating in cgs units is estimated as $\dot{E} = n^2 S_1$, where $n$ is the proton number density. For example, the radiative total X-ray heating for gas at solar metallicity is larger by a factor of ~11 than that of metal-free gas.

### 2.7. Cosmological “Zoom-in” Initial Conditions

The cosmological “zoom-in” initial conditions used in this study are described in detail in Oser et al. (2010). To achieve sufficiently high enough resolution to robustly model the evolution of galaxies, a subvolume 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_{dm} = 0.216$, $\Omega_\Lambda = 0.74$, $\sigma_8 = 0.77$, and $n_s = 0.95$). At any given snapshot we trace back all particles close to the halos of interest from redshift zero. We replace those particles with higher-resolution gas and dark matter particles. Then new, high-resolution initial conditions are resimulated from redshift $z = 43$ to $z = 0$.

The simulations have been performed at two resolutions: (1) The reference resolution has a mass resolution for the star and gas particles of $m_{*,gas} = 4.2 \times 10^6 h^{-1} M_\odot$, and the dark matter particles have $m_{dm} = 2.5 \times 10^7 h^{-1} M_\odot$. We use the comoving gravitational softening lengths $\epsilon_{gas,star} = 400$ pc $h^{-1}$ for the gas and star particles and $\epsilon_{halo} = 890$ pc $h^{-1}$ for the dark matter, which are scaled with the square root of the mass ratio following Dehnen (2001). (2) The high-resolution simulations have been performed with eight times better mass resolution than the reference resolution, with $m_{*,gas} = 5.3 \times 10^5 h^{-1} M_\odot$ and $m_{dm} = 3.1 \times 10^6 h^{-1} M_\odot$, and twice better spatial resolution with $\epsilon_{gas,star} = 200$ pc $h^{-1}$ and $\epsilon_{halo} = 450$ pc $h^{-1}$. In the Appendix we discuss the resolution convergence of the galaxy physical quantities.

The simulated halo masses are in the range $1.4 \times 10^{12} M_\odot \lesssim M_{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_{vir} \equiv r_{200}$, the radius where the spherical overdensity drops below 200 times the critical density of the universe at a given redshift).

In order to study the effects of each of the physical modules we included, we investigate the galaxy properties by running all 30 zoom-in simulations in the reference resolution with four different models:

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Fraction of baryons converted to stars at $z = 0$ relative to the universal baryon fraction as a function of halo mass for central galaxies in different models: fiducial model (red circles), model without AGN feedback effect (NoAGN; black triangles), model without metal heating effect (NoZHeating; green stars), and model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue squares). The predictions of the fraction of stars inferred from observations via the abundance-matching estimates are shown in black and blue solid lines from Moster et al. (2013) and Kravtsov et al. (2014), respectively, with the $1\sigma$ scatter region shown by dotted lines. Compared to the fiducial model, the NoAGN model predicts significantly higher conversion efficiencies, especially at high halo masses, while the effect of metal-induced heating mechanisms on the final stellar mass seems negligible.</figcaption></figure>

(1) Fiducial: the reference model that includes all physical modules listed above, i.e., mechanical and radiative AGN feedback, stellar feedback with snowplow SN feedback, metal cooling and enrichment, and the metal heating effect from photoelectric heating and cosmic X-ray background. Feedback-related numerical parameters were calibrated against the black hole mass-stellar velocity dispersion ($M_{BH} - \sigma_*$) relation and the baryonic conversion rate at $z = 0$ (see Figures 1 and 2).

(2) NoAGN: without black hole and AGN feedback. This model isolates the effect of the AGN feedback.

(3) NoZHeating: same as the fiducial model, but without the new ingredients of metal heating effects listed in Section 2.6, including metallicity-dependent Compton heating, photoelectric heating, and cosmic X-ray background heating.

(4) ThSNnoMetal: from a previous paper (Choi et al. 2015), this model uses thermal SN feedback (Springel & Hernquist 2003) instead of ejective SN feedback described in Section 2.3. This model does not include metal enrichment, metal-induced heating/cooling, and early stellar feedback. Note that we presented 20 halos in Choi et al. (2015), but we have performed 10 more zoom-in simulations for a fair comparison.

A comparative summary of the input physics for each model is given in Table 1. Throughout the paper, we focus on the central galaxies and discuss the respective role of various feedbacks on the formation of massive galaxies.

## 3. Variations of Galaxy Properties with Physical Models

In the following we will present simulated physical properties of the main four models and compare them to observations of various relations: stellar mass-halo mass relation, the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Black hole mass–stellar velocity dispersion ($M_{BH} - \sigma_{*}$) relation at $z=0$ of the fiducial model (red circles), the model without a metal heating effect (NoZHeating; green stars), and the model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue squares). The observed $M-\sigma_{*}$ relations of the elliptical galaxies are shown with dotted lines: the blue dotted line is from Kormendy & Ho (2013), and the black dotted line is from McConnell & Ma (2013).</figcaption></figure>

relation between black hole mass and velocity dispersion, galaxies sizes, and stellar mass.

### 3.1. The Relation between Stellar Mass and Halo Mass

Figure 1 compares the simulated galaxies to the abundance-matching results of Moster et al. (2013) and Kravtsov et al. (2014). The mass of star particles within 10% of the virial radius $r_{10}$ is defined as the stellar mass of simulated central galaxies. Note that we use the total baryonic and dark matter mass within the $r_{vir}$ for $M_{vir}$, while the abundance-matching models use the dark-matter-only simulations. As the baryonic physics can alter the dark matter distribution and reduce the halo masses, $M_{vir}$ can be overestimated by 10% in the abundance-matching models, especially for the small halo mass range $M_{vir} < 10^{10} M_{\odot}$ (e.g., Sawala et al. 2013; Martizzi et al. 2014). For our galaxy mass range, massive galaxies with $M_{vir} > 10^{12} M_{\odot}$, the stellar masses measured in the abundance-matching model become increasingly sensitive to the aperture used to measure the stellar mass. Recently, Kravtsov et al. (2014) demonstrated that improved photometric techniques used to measure stellar mass lead to a significant effect on the stellar mass–halo mass relation predicted by abundance-matching models. Their prediction is shown by the blue solid line. We see that the baryonic conversion efficiency in our models is very close to that inferred from abundance matching when AGN feedback is included. Models without AGN feedback overpredict the baryon abundances by a factor of 2–3 at $z=0$.

### 3.2. The Relation between Black Hole Mass and Stellar Velocity Dispersion

Figure 2 shows the relation between the mass of central supermassive black holes and the velocity dispersion of galaxies and its comparison with observational data from

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top: averaged star formation rate over time for the 30 central galaxies for different models: fiducial model (red), model without AGN feedback effect (NoAGN; black), model without metal heating effect (NoZHeating; green), and model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue). The solid lines show the average value, and shaded regions illustrate the $1\sigma$ scatter. For clarity of display, the $1\sigma$ region of the NoZHeating model is not shown. Bottom: same as in the top panel, but for the median specific star formation rates. The dotted black line indicates the specific star formation rates equal to $0.3/t_{H}$, a commonly used criterion separating quiescent and star-forming galaxies (e.g., Franx et al. 2008). The NoAGN feedback model (black) stays above this criterion, constantly star-forming throughout the evolution.</figcaption></figure>

Kormendy & Ho (2013) and the black dotted line from McConnell & Ma (2013). We first determine the half-mass radii of stars within 10% of the virial radius $r_{10}$ projected along the 20 randomly chosen directions. The mean value is taken as a representative effective radius $r_{e}$ of each galaxy. Then, the line-of-sight velocity dispersion $\sigma_{*}$ is calculated by considering all the stars within the half of the effective radius (0.5 $r_{e}$) along the three principal axes. Our new simulations show higher black hole masses and higher velocity dispersions compared to our previous work, shown by the blue squares (Choi et al. 2015). This is because we have more gas available for both star formation and black hole accretion via enhanced gas metal-line cooling and recycled gas ejected from the death of stars in the old stellar population.

### 3.3. Star Formation Rates

In Figure 3, we address “quenching” of star formation and show the dramatic effects that the various physical processes
have on the mean star formation rate within 10% of the virial radius $r_{10}$ (top panel) and the specific star formation rate (bottom panel) of 30 central galaxies in our four models as a function of redshift. The star formation rate of all models peaks at $z \sim 3$ and drops rapidly afterward, except for the NoAGN model. With mechanical and radiative AGN feedback included, the star formation is very rapidly quenched by the effective removal of gas via AGN-driven large-scale winds (see Pandya et al. 2016, for detailed discussion on the quenching timescale of the fiducial galaxies).

Compared to the NoZHeating model, metal-dependent heating in the fiducial run has the effect of reducing the star formation rate around the redshift of quenching at $2 < z < 4$. This is because both photoelectric heating and cosmic X-ray background heating from accreting black holes mainly affect diffuse and small stellar systems, effectively suppressing their growth.

With metal-induced heating, the formation of low-mass stellar systems is efficiently delayed, but the heated gas will later cool down and come back to the central system; therefore, the star formation rate later compensates at $z < 1$.

In the bottom panel, we show the median specific star formation rates of four models. We also show the commonly used criteria separating quiescent and star-forming galaxies (Franx et al. 2008), with the specific star formation rates equal to $0.3/t_H$ shown by the dotted black line. The NoAGN feedback model (black) stays above this criterion constantly, star-forming throughout the simulation. By this criterion the fiducial models cross over to the red sequence and become quiescent between redshifts $z = 1.5$ and $z = 0.3$. Again compared to our previous models (ThSNnoMetal, Choi et al. 2015), allowing for the additional metal cooling from recycled gas increases star formation rates (red curves) but leaves them at $z = 0$, still more than two orders of magnitude below the NoAGN models.

### 3.4. Galaxy Sizes and Velocity Dispersions

In Figure 4, we show the projected half-mass radius $r_e$, i.e., the radius that encloses 50% of the stellar mass in projection, as a function of the galaxy mass at $z = 0$. As described above, we calculate the half-mass radii of stars within $r_{10}$ projected along the 20 randomly chosen directions of the main stellar body, and the mean value is taken as a representative effective radius $r_e$ of each galaxy. The sizes of the fiducial galaxies at high stellar masses simulated with AGN feedback are in good agreement with observations of early-type galaxies from Nipoti et al. (2009), but we tend to produce typically too small sizes for lower-mass galaxies.

The sizes of galaxies in the NoAGN feedback models are even smaller; their effective radii are $\sim 5$ times smaller compared to the observed ones at a given stellar mass. The continuous star formation without AGN feedback in the central regions of galaxies leads to a concentrated stellar mass profile (Martizzi et al. 2012). Moreover, the galaxies are less puffed up by minor mergers as the absence of AGN feedback decreases the fraction of accreted material (Dubois et al. 2013, 2016). In the absence of AGN feedback, in situ star formation dominates over accreted star formation at all times (Lackner et al. 2012; Hirschmann et al. 2013); therefore, relatively fewer stars are added to the outskirts of galaxies.

The inclusion of physical heating processes due to the metal content of gas slightly decreases the galaxy sizes. The metal-dependent heating mainly affects the smaller systems,

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Projected stellar half-mass radii of the simulated galaxies vs. stellar masses at redshifts $z = 0$ for four models: fiducial model (red), model without AGN feedback effect (NoAGN; black), model without metal heating effect (NoZHeating; green), and model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue). The red solid line indicates the observed size-mass relation of the SLACS sample of local early-type galaxies from Nipoti et al. (2009), with the $1\sigma$ scatter given by the dotted lines. The blue solid line is from Shen et al. (2003).</figcaption></figure>

decreasing the total mass of accreted stars. Therefore, the accreted star fraction decreases when we include metal-dependent heating, and this leads to smaller galaxy sizes.

Similarly in the ThSNnoMetal model, even more small stellar systems are accreted to the central galaxy, due to insufficient stellar feedback to suppress star formation in the small building blocks, leading to a significant size growth of galaxies at late times. In addition, the lack of metal cooling in this model prevents gas from concentrating in the central region of galaxy, thus contributing to the size increase.

From Oser et al. (2012), our galaxy evolution code has been improved with the addition of the density-independent SPH formulation, the metal enrichment, and the subsequent metal-line cooling. All of these changes enhance the gas cooling and star formation especially for the central galaxies. In addition, the SN feedback model is changed from thermal to mechanical, which drives more powerful winds. However, this change also challenges the late-time star formation quenching in central galaxies, as the material removed from lower-mass progenitors can be accreted to more massive ones at later times (e.g., Oppenheimer & Davé 2008; Hirschmann et al. 2013). The previous version of our physical model was able to reproduce the quiescent galaxies via gravitational heating (e.g., Johansson et al. 2009b) and sufficient size growth afterward even in the absence of AGN feedback, but this seems to be due to some missing physics. The gravitational heating helps but by itself is insufficient to quench extended star formation in massive galaxies. Overall, appropriately predicting the final effective radius and its evolution provides a strong test of the accuracy of the physical modeling. A further study of the physical mechanisms determining the size evolution of galaxies
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Projected velocity dispersion measured within $0.5r_e$ as a function of stellar mass at $z=0$ for central galaxies in different models: fiducial model (red), model without AGN feedback effect (NoAGN; black), model without metal heating effect (NoZHeating; green), and model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue). The observational relations for local early-type galaxies from Nipoti et al. (2009) and Taranu et al. (2013) are shown by the black and blue solid lines, respectively, with the dotted lines indicating the $1\sigma$ scatter.</figcaption></figure>

will be discussed in a forthcoming paper (E. Choi et al. 2016, in preparation).

Figure 5 shows the Faber–Jackson relation (Faber & Jackson 1976), the relation between the stellar velocity dispersion $\sigma_*$ and stellar mass for simulated galaxies at $z=0$. The line-of-sight velocity dispersions $\sigma_*$ are measured for all stars within $0.5 \times r_e$ along the three principal axes. The observed relations with $1\sigma$ ranges are from Nipoti et al. (2009) based on the SLACS sample of local early-type galaxies at $z=0$ and from Taranu et al. (2013) based on the Sloan Digital Sky Survey (SDSS) catalog. All simulations except the NoAGN feedback model are in good agreement with the data, but the fiducial model tends to have overall higher velocity dispersions. Compared to the ThSNnoMetal model, the fiducial model shows higher velocity dispersion owing to the increased central star formation driven by enhanced gas metal-line cooling and recycled gas ejected from the old stellar population. The NoAGN model significantly overpredicts the velocity dispersion of normal massive elliptical galaxies, due to the enhanced star formation in the central region of galaxies.

### 3.5. X-Ray Luminosity and Gas Mass Fraction

Recent X-ray observations have shown that hot gas halos are ubiquitous from galaxy clusters, groups, and galaxies with mass down to $M_* \sim 10^{11} M_\odot$ (e.g., Anderson et al. 2015). We explore the impact of various feedback models on hot gas in the simulated galaxy halos. In Figure 6, we show the X-ray luminosity of the hot gas versus the stellar velocity dispersion $\sigma_*$ and total mass within $5r_e$ of the simulated central galaxies at $z=0$ for four models. Following Cen et al. (1995), we calculate the X-ray luminosity for the simulated galaxies in Chandra bands (0.3–8 keV), including bremsstrahlung radiation and metal-line emission from all relevant species measured. For the ThSNnoMetal model, we assume solar

abundance. For $L_X - \sigma_*$, we show the observed relation from Boroson et al. (2011) and data from Goulding et al. (2016), who archived Chandra X-ray observation of local elliptical galaxies and combined the data from the MASSIVE galaxy survey (Ma et al. 2014) and the ATLAS3D survey. For $L_X - M_{\text{Total}}(<5r_{\text{eff}})$, we compare our simulated galaxy properties with those of Kim & Fabbiano (2013) and Forbes et al. (2017), who recently reexamined the dynamical masses of the galaxies with the globular cluster kinematics.

In the fiducial model, the mechanical energy output of the black hole at the peak epoch of galaxy formation drives galactic outflows and expels large amounts of hot gas from the galactic potential. These large-scale winds result in overall lower gas density and thus significantly reduce the gas X-ray luminosity as $L_X \propto n^2$. The fiducial model with mechanical AGN feedback reproduces X-ray luminosity within the observed range, as well as the observed scatter with the two orders of magnitude spread in $L_X$. Conversely, the NoAGN model shows much higher X-ray luminosity up to $\sim 1–2$ orders of magnitude higher compared to the observations elliptical galaxies. In our simulations, X-ray luminosity of hot halos of massive galaxies is most sensitive to AGN feedback, consistent with many previous works, including McCarthy et al. (2010) and Le Brun et al. (2014, 2015, 2016). In particular, Le Brun et al. (2014) showed that their AGN feedback model with a heating temperature of $\Delta T_{\text{heat}} = 10^8$ K produced the observed X-ray luminosity—mass relation in agreement with observations, and they also noted that X-ray luminosities of simulated galaxies are sensitive to the choice of the heating temperature parameter in their AGN feedback model.

In Figure 7 we show the gas mass fraction of all simulated halos and compare them to observations. We measure the fraction of gas mass to the total mass within $r_{500}$, the radius where the spherical overdensity drops below 500 times the critical density of the universe at $z=0$. The observed $f_{\text{gas,500}} - M_{500}$ relation derived from 43 nearby galaxy groups with $M_{500} = 10^{13-14} M_\odot$ from Sun et al. (2009) is shown by the black solid line, with $1\sigma$ scatter shown by dotted lines. We also overplot the observed gas mass fraction of galaxy groups from the same paper with black crosses. In the fiducial model, the AGN-driven strong wind results in much lower gas mass fraction compared to the NoAGN model, especially for the low-mass galaxies, because a lot of hot gas is swept up and removed from the halo by AGN-driven galactic outflows. The inclusion of metal-induced heating also lowers the gas fraction, but its effect is relatively small.

### 3.6. Mass–Metallicity Relation

The relations between stellar metallicity and mass of simulated galaxies are shown in Figure 8 and compared with observations from SDSS from Panter et al. (2008) and Gallazzi et al. (2005). The NoAGN feedback shows systemically higher metallicity by $\sim 0.05$ dex compared to the fiducial model simulated with AGN feedback. The higher metallicity in the NoAGN feedback model is driven by the excessive late-time star formation out of metal-enriched gas in the absence of an effective star formation quenching mechanism. A more detailed investigation of the impact of AGN feedback on individual metal abundances and metallicity gradients of stellar populations will be given in M. Hirschmann et al. (2017, in preparation). We also refer the reader to Hirschmann et al. (2017) for the further details on the cosmic evolution of gas
<figure><img src="image_7.png" /><figcaption>**Figure 6.** $L_X-\sigma_*$ relation (left) and $L_X-M_{\text{total}}(<5r_c)$ relation (right) at $z=0$ for four models: fiducial model (red), model without AGN feedback effect (NoAGN; black), model without metal heating effect (NoZHeating; green), and model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue). Observed relations and data for the normal early-type galaxies are from Boroson et al. (2011) and Goulding et al. (2016) (left) and from Kim & Fabbiano (2013) and Forbes et al. (2017) (right).</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Enclosed gas mass fraction within $r_{500}$, $f_{\text{gas,500}} - M_{500}$ relation for four models: fiducial model (red), model without AGN feedback effect (NoAGN; black), model without metal heating effect (NoZHeating; green), and model with thermal SN feedback and without metal enrichment (ThSNnoMetal; blue). The NoZHeating model shows a lower gas fraction compared to the fiducial model. The black crosses and black solid lines are the observed data and fitted $f_{\text{gas,500}} - M_{500}$ relation of the nearby galaxy groups with $M_{500} = 10^{13-14} M_{\odot}$ from Sun et al. (2009).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Mass–metallicity relation—stellar metallicity vs. galaxy stellar mass at $z=0$ for three models: fiducial model (red circles), model without AGN feedback effect (NoAGN; black triangles), and model without metal heating effect (NoZHeating; green stars). ThSNnoMetal model is not shown, as this model assumed a primordial metal abundance and did not include metal enrichment. The blue dotted line and black solid line indicate the local mass–metallicity relation for the SDSS galaxy population from Panter et al. (2008) and Gallazzi et al. (2005), respectively. The fiducial model is in good agreement with observations.</figcaption></figure>

metallicity and optical-emission line ratios in our simulated galaxies.

## 3.7. The Two Phases of Galaxy Formation

Recent observations have established that early-type galaxies form and become red and dead early but continue to grow in mass and size without much late star formation (e.g., Daddi et al. 2005; Trujillo et al. 2006; Buitrago et al. 2008; Szomoru et al. 2012). These developments challenge classical formation

models of early-type galaxies via monolithic collapse (Eggen et al. 1962) and equal-mass major mergers (e.g., Toomre & Toomre 1972), but favor a two-phase formation scenario (Naab et al. 2007; Oser et al. 2010). In this scenario, the cores of elliptical galaxies form early at $2 < z < 6$ by dissipational processes and cold gas flows (e.g., Dekel et al. 2009) and by merging of smaller structures of stars and gas. Then, the outer parts of elliptical galaxies grow by accretion of old stars and galaxy mergers of all mass ratios (Bezanson et al. 2009;
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Stellar origin diagram for all stars within $r_{10}$ of 30 galaxies at $z=0$ in each model: (a) fiducial, (b) NoAGN, (c) NoZHeating, and (d) ThSNnoMetal. Upper panels show the histogram of formation radii of stars, and the right panels show the star formation histories. The vertical red dotted line shows the 10% of virial radius, $r_{10}$.</figcaption></figure>

van Dokkum et al. 2010; Hirschmann et al. 2012; Sales et al. 2012). But this subsequent buildup of the stellar envelope is dominated by nondissipational processes via dry “minor” mergers, as major mergers are rare at late times. This two-phase formation scenario can explain the galaxy size growth and slight decline in velocity dispersion (Oser et al. 2012). The outer parts of ellipticals are found to be old and metal-poor (Greene et al. 2015), as would be expected if they had accreted at late times from low-mass dwarf companions (Lackner et al. 2012; Hirschmann et al. 2015).

We investigate the effect of various physics implemented in our simulation on the fundamental formation and assembly of galaxies. Figure 9 visualizes when and at which radius a star

was born for star particles ending up within 10% of the virial radius of a present-day galaxy. We have stacked all 30 simulated galaxies for each physics model: (1) fiducial (mechanical AGN and mechanical SN feedback model), (2) NoAGN (without AGN feedback), (3) NoZHeating (without metal-induced heating effect), and (4) ThSNnoMetal (with thermal SN feedback and without metal enrichment). For stars that reside in final galaxies, some are made in situ, within the $r_{10}$, while some are made ex situ, outside of the $r_{10}$, and later accreted. The vertical red dotted line indicates the 10% of virial radius, $r_{10}$, that clearly separates these two phases, as shown in the histogram of the formation radii in the upper panels of Figure 9. We found two peaks at relatively similar locations in the histogram of all four physics
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Stellar origin diagram for (a) low-mass ($8.5 \times 10^{10} < M_{*}/M_{\odot} < 1.25 \times 10^{11}$), (b) intermediate-mass ($1.25 \times 10^{11} < M_{*}/M_{\odot} < 2.1 \times 10^{11}$), and (c) high-mass ($2.1 \times 10^{11} < M_{*}/M_{\odot} < 5 \times 10^{11}$) galaxies from the fiducial model. The ratio of in situ and ex situ stars is found to be sensitive to the galaxy mass.</figcaption></figure>

models: for the in situ formed stars at $\log(r/r_{vir}) \sim -2.5$ and for the ex situ formed stars at $\log(r/r_{vir}) \sim 0.7$, respectively.

However, the ratio of in situ to ex situ stars is found to strongly vary with the feedback model we include. Our fiducial model galaxies on average have 40% of stars formed in situ (see panel (a)), but this fraction increases to 64% when we exclude the AGN feedback powered by the supermassive black holes. In addition, in the NoAGN feedback model the formation radius peak of in situ star formation is located further inside of galaxies ($\log(r/r_{vir}) \sim -3.0$) compared to the fiducial model. The stellar feedback solely is not strong enough to generate galactic outflows in high-mass galaxies and cannot quench in situ star formation. The central accretion of ambient gas from stellar mass loss continuously triggers star formation (Ciotti & Ostriker 1997), and this late-time star formation contributes significantly to the in situ fraction of stars.

In the NoZHeating model, when we exclude the metal-induced heating, the fraction of ex situ formed stars increases from 60% to 65%. These extra heating effects (photoelectric and cosmic X-ray background heating) mainly affect the small systems, which later accrete onto the central galaxies, preventing small blobs from condensing. This tends to decrease the accreted star fraction.

The ratio of in situ and ex situ stars also strongly depends on *how* the feedback prescription is implemented in the model. In the ThSNnoMetal model, where we use thermal SN feedback, which does not produce strong SN winds compared to the mechanical SN and early stellar feedback, the fraction of ex situ formed stars increases to 82%.

In Figure 10 we divide our simulated galaxies into three mass bins, high ($2.1 \times 10^{11} < M_{*}/M_{\odot} < 5 \times 10^{11}$), intermediate ($1.25 \times 10^{11} < M_{*}/M_{\odot} < 2.1 \times 10^{11}$), and low ($8.5 \times 10^{10} < M_{*}/M_{\odot} < 1.25 \times 10^{11}$), so that we have 10 galaxies in each mass bin. The ratio of in situ and ex situ stars is found to be sensitive to the galaxy mass, as found in Oser et al. (2010), Qu et al. (2017), and Rodriguez-Gomez et al. (2016), i.e., more stars form in situ for low-mass galaxies and more form ex situ for higher-mass galaxies. Also, there is a continuous trend of in situ star formation quenching redshift with stellar mass. While low-mass galaxies show high in situ star formation rate until redshift $z = 1$, high-mass galaxies quench early before $z = 1$. Some of our high-mass galaxies are rejuvenated shortly, due to the recycled gas ejected from the

old stellar population within the innermost region of the galaxies ($\log r_{10}/r_{vir} \sim -3$).

## 4. Summary and Discussion

In this work, we present a more comprehensive treatment of the physical processes consequent to the allowance of chemical evolution highlighting the extra heating processes due to the interaction of metals and radiation. We then used updated simulations of 30 massive galaxies with halo mass of $M_{vir} \sim 10^{12-13} M_{\odot}$ at $z=0$ to make comparisons between the simulation and observations. The feedback models account for (1) AGN feedback via broad absorption line winds and X-ray radiation heating; (2) stellar feedback via UV heating and stellar winds from young and old stars and SNe I and II, as well as associated mass and metal generation; and (3) an additional metal heating effect via photoelectric heating and cosmic X-ray background heating from accreting black holes in background galaxies. The overall energy and momentum budget from all feedback effects has been shown to generate realistic galaxy properties (fiducial model). We also explore the detailed role of separate feedback variables by running three additional sets of simulations without AGN feedback (NoAGN model), without metal heating effects (NoZHeating model), and without the suite of mechanical stellar feedback but with thermal SN feedback (ThSNnoMetal model), to better understand the separate effects of various physical processes. The input physics variation of the simulation suite is summarized in Table 1.

We find that the AGN feedback plays the most dominant role in reproducing the basic physical properties of observed massive early-type galaxies. The mechanical and radiative AGN feedback generates strong, large-scale outflows during early stages of galaxy evolution and efficiently quenches in situ star formation in massive galaxies. In the absence of AGN feedback, the conversion efficiency of baryons into stars in the central galaxies is increased by a factor of 3 (Figure 1), and all simulated galaxies never quench and keep forming stars throughout to $z=0$ (Figure 3). The formation of central galaxies is dominated by in situ star formation without AGN feedback (Figure 9(b)), and this leads to very compact stellar cores with a factor of five smaller effective radii than in the observed mass-size relation and high velocity dispersions
(Figures 4, 5). AGN feedback can efficiently drive out and remove gas from the galaxy halos and prevent large gas concentrations in halos, resulting in much lower X-ray luminosities and gas mass fractions (Figures 6, 7).

We show that the metal-induced heating mechanisms primarily affect diffuse and smaller galaxies at high redshift. They prevent baryons from condensing and accumulating in dense small blobs mainly at high redshift and therefore delay star formation. Star formation, however, later compensates, as the heated gas cools and comes back to the central galaxies. Therefore, metal-induced heating mechanisms have a negligible effect on the final stellar mass, as well as the final gas properties of central galaxies. However, they certainly reduce the mass of small stellar systems later accreted to the central galaxy. The fraction of accreted stars is reduced from 65% to 60% when we add metal-induced heating (Figure 9(c)), and this results in slightly smaller stellar galaxy sizes, as they tend to have fewer dissipationless low-redshift mergers. Although the metal-induced heating overall produces weak effects on this stellar mass scale ($M_* \sim M^{11-12} M_\odot$), they can have a dramatic impact on the evolution of dwarf galaxies with much smaller masses. For example, Forbes et al. (2016) recently showed that the photoelectric heating plays a more important role than SN feedback in regulating star formation in dwarf galaxies with a stellar mass of $10^7 M_\odot$.

Although AGN feedback seems to exert the dominant impact on massive galaxy formation, it is also crucial to properly include the stellar feedback to regulate the star formation in small building blocks and to suppress in situ star formation at high redshift (Hirschmann et al. 2013). We show that inclusion of thermal SN feedback instead of ejective SN feedback yields final massive galaxies dominated by hierarchically assembled small galaxies and therefore a much higher accreted fraction (Figure 9(d)).

Overall, we achieve fairly good agreement between our fiducial model and observational constraints in the local universe, although some small discrepancies still remain. To some extent, this may be linked to further physical mechanisms that have not been included so far, such as cosmic rays, magnetohydrodynamics, and effects of relativistic jets from AGNs. In particular, the simulations do not include cosmic-ray-driven winds (e.g., Booth et al. 2013; Hanasz et al. 2013; Salem & Bryan 2014; Wiener et al. 2016) and the effect of runaway OB stars that migrate into low-density regions (Li et al. 2015), and their absence may contribute to the excess of stellar masses in massive galaxies. In future work, we will extend our study to include several of these effects, all of which can act directly on the gas and reduce the star formation rate.

We are grateful to the anonymous referee for very helpful comments on the manuscript. We also thank Volker Springel for making the GADGET-3 code available and Sergey Sazonov for helpful discussions. Support for this work is provided by NASA through grant no. HST Cycle 23 AR-14287 from the Space Telescope Science Institute. T.N. acknowledges support from the DFG cluster of excellence “Origin and Structure of the Universe.” M.H. acknowledges financial support from the European Research Council via an Advanced Grant under grant agreement no. 321323 NEOGAL. Numerical simulations were run on the computer facility of the Princeton Institute of Computational Science and engineering.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Cosmic evolution of virial radius, virial mass, effective radius, and stellar mass of the fiducial run of halo m0290 at two different resolutions (black: high resolution; red: fiducial resolution). We find no significant difference in the physical properties of central galaxies with different resolution.</figcaption></figure>

# Appendix
## Convergence with Resolution

We test the convergence of the physical properties of central galaxies with respect to resolution using two realizations of the halo m0290 with $M_{vir} \sim 4 \times 10^{12} M_\odot$. In Figure 11 we label the fiducial m0290 simulation “fiducial resolution” and run a higher-resolution simulation. The “high-resolution” run has a gas and dark matter particle mass eight times smaller and a gas softening length twice as small as the reference resolution. The high-resolution run is simulated with the number of SN events reduced by a factor of three but with increased SN feedback energy in order to account for clumping of star formation, but otherwise has identical parameters. The agreement between the high-resolution and low-resolution runs is very good.
References

Agertz, O., & Kravtsov, A. V. 2015, ApJ, 804, 18

Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25

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

Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806

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

Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822

Barai, P., Murante, G., Borgani, S., & Osterriker, J. P., & Ryu, D. 1995, ApJ, 451, 436

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

Bondi, H. 1952, MNRAS, 112, 195

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

Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJL, 777, L16

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

Boroson, B., Kim, D.-W., & Fabbiano, G. 2011, ApJ, 729, 12

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

Cen, R., Kang, H., Ostriker, J. P., & Ryu, D. 1995, ApJ, 451, 436

Cheung, E., Bundy, K., Cappellari, M., et al. 2016, Natur, 533, 504

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

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

Ciotti, L., & Ostriker, J. P. 1997, ApJL, 487, L105

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

Crenshaw, D. M., Kraemer, S. B., & George, I. M. 2003, ARA&A, 41, 117

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

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

Dalla Vecchia, C., & Schaye, J. 2008, MNRAS, 387, 1431

Debuhr, J., Quataert, E., Ma, C.-P., & Hopkins, P. 2010, MNRAS, 406, L55

Dehnen, W. 2001, MNRAS, 324, 273

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

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Natur, 457, 451

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

Draine, B. 2010, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Draine, B. T. 1978, ApJS, 36, 595

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

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

Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748

Faber, S. M., & Jackson, R. E. 1976, ApJ, 204, 668

Faucher-Giguère, C.-A., Hopkins, P. F., Kereš, D., et al. 2015, MNRAS, 449, 987

Finlator, K., & Davé, R. 2008, MNRAS, 385, 2181

Forbes, D. A., Alabi, A., Romanowsky, A. J., et al. 2017, MNRAS, 464, L26

Forbes, J. C., Krumholz, M. R., Goldbaum, N. J., & Dekel, A. 2016, Natur, 535, 523

Ford, A. B., Oppenheimer, B. D., Davé, R., et al. 2013, MNRAS, 432, 89

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

Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41

Gaspari, M., Brighenti, F., & Temi, P. 2012, MNRAS, 424, 190

Genel, S., Fall, S. M., Hernquist, L., et al. 2015, ApJL, 804, L40

Goulding, A. D., Greene, J. E., Ma, C.-P., et al. 2016, ApJ, 826, 167

Greene, J. E., Janish, R., Ma, C.-P., et al. 2015, ApJ, 807, 11

Greene, J. E., Zakamska, N. L., Ho, L. C., & Barth, A. J. 2011, ApJ, 732, 9

Haardt, F., & Madau, P. 2001, in Clusters of Galaxies and the High Redshift Universe Observed in X-rays, ed. D. M. Neumann & J. T. V. Tran (Saclay: CEA), 64

Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJL, 777, L38

Hirschmann, M., Charlot, S., Feltre, et al. 2017, MNRAS, submitted (arXiv:1706.00010)

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

Hirschmann, M., Naab, T., Ostriker, J. P., et al. 2015, MNRAS, 449, 528

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., Kereš, D., Óñorbe, J., et al. 2014, MNRAS, 445, 581

Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522

Hoyle, F., & Lyttleton, R. A. 1939, PCPS, 34, 405

Hu, C.-Y., Naab, T., Walch, S., Glover, S. C. O., & Clark, P. C. 2016, MNRAS, 458, 3528

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

Johansson, P. H., Naab, T., & Burkert, A. 2009a, ApJ, 690, 802

Johansson, P. H., Naab, T., & Ostriker, J. P. 2009b, ApJL, 697, L38

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

Kim, D.-W., & Fabbiano, G. 2013, ApJ, 776, 116

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Kravtsov, A., Vikhlinin, A., & Meshscheryakov, A. 2014, arXiv:1401.7329

Kroupa, P. 2001, MNRAS, 322, 231

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

Le Brun, A. M. C., McCarthy, I. G., & Melin, J.-B. 2015, MNRAS, 451, 3868

Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270

Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2016, MNRAS, 466, 4442

Li, M., Ostriker, J. P., Cen, R., Bryan, G. L., & Naab, T. 2015, ApJ, 814, 4

Ma, C.-P., Greene, J. E., McConnell, N., et al. 2014, ApJ, 795, 158

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

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

Martizzi, D., Mohammed, I., Teyssier, R., & Moore, B. 2014, MNRAS, 440, 2290

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

McCarthy, I. G., Schaye, J., Ponman, T. J., et al. 2010, MNRAS, 406, 822

McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Naab, T., Johansson, P. H., Ostriker, J. P., & Efstathiou, G. 2007, ApJ, 658, 710

Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357

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

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

Nyman, L.-A., Booth, R. S., Carlstrom, U., et al. 1992, A&AS, 93, 121

Oppenheimer, B. D., Crain, R. A., Schaye, J., et al. 2014, MNRAS, 460, 2157

Oppenheimer, B. D., & Davé, R. 2006, MNRAS, 373, 1265

Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577

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

Ostriker, J. P., Choi, E., Ciotti, L., Novak, G. S., & Proga, D. 2010, ApJ, 722, 642

Pandya, V., Brennan, R., Somerville, R. S., et al. 2016, arXiv:1611.03869

Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2008, MNRAS, 391, 1117

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

Proga, D., Stone, J. M., & Kallman, T. R. 2000, ApJ, 543, 686

Qu, Y., Helly, J. C., Bower, R. G., et al. 2017, MNRAS, 464, 1659

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

Renaud, F., Bournaud, F., Emsellem, E., et al. 2013, MNRAS, 436, 1836

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

Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312

Sales, L. V., Navarro, J. F., Theuns, T., et al. 2012, MNRAS, 423, 1544

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

Schmidt, M., Schneider, D. P., & Gunn, J. E. 1995, AJ, 110, 68

Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2009, ApJ, 690, 20

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

Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877

Simpson, C. M., Bryan, G. L., Hummels, C., & Ostriker, J. P. 2015, ApJ, 809, 69

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

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

Spitzer, L., Jr. 1948, ApJ, 107, 6

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Di Matteo, T., & Hernquist, L. 2005a, ApJL, 620, L79

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

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289
Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.-S. 2015, MNRAS, 448, 1504

Stinson, G. S., Brook, C., Macciò, A. V., et al. 2013, MNRAS, 428, 129

Stone, N. C., Kuepper, A. H. W., & Ostriker, J. P. 2016, arXiv:1606.01909

Strömgren, B. 1939, ApJ, 89, 526

Sun, M., Voit, G. M., Donahue, M., et al. 2009, ApJ, 693, 1142

Szomoru, D., Franx, M., & van Dokkum, P. G. 2012, ApJ, 749, 121

Taranu, D. S., Dubinski, J., & Yee, H. K. C. 2013, ApJ, 778, 61

Tasker, E. J. 2011, ApJ, 730, 11

Toomre, A., & Toomre, J. 1972, ApJ, 178, 623

Tremmel, M., Karcher, M., Governato, F., et al. 2016, arXiv:1607.02151

Tremonti, C. A., Moustakas, J., & Diamond-Stanic, A. M. 2007, ApJL, 663, L77

Trujillo, I., Förster Schreiber, N. M., Rudnick, G., et al. 2006, ApJ, 650, 18

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

Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031

Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238

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

Weingartner, J. C., & Draine, B. T. 2001, ApJS, 134, 263

Wiener, J., Pfrommer, C., & Oh, P. 2016, arXiv:1608.02585

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

Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009b, MNRAS, 399, 574

Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278

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

Zakamska, N. L., Hamann, F., Pâris, I., et al. 2016, MNRAS, 459, 3144