# Formation of quenched massive galaxies in FIRE cosmological zoom-in
simulations
with multi-channel AGN feedback

Lindsey Byrne Department of Physics and Astronomy and CIERA, Northwestern
University, Evanston, IL 60201, USA  Claude-André Faucher-Giguère Department
of Physics and Astronomy and CIERA, Northwestern University, Evanston, IL
60201, USA  Sarah Wellons Department of Astronomy, Van Vleck Observatory,
Wesleyan University, 96 Foss Hill Drive, Middletown, CT 06459, USA  Philip F.
Hopkins California Institute of Technology, TAPIR, Mailcode 350-17, Pasadena,
CA 91125, USA  Daniel Anglés-Alcázar Department of Physics, University of
Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA  Center
for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New
York, NY 10010, USA  Imran Sultan Department of Physics and Astronomy and
CIERA, Northwestern University, Evanston, IL 60201, USA  Nastasha Wijers
Department of Physics and Astronomy and CIERA, Northwestern University,
Evanston, IL 60201, USA  Jorge Moreno Department of Physics and Astronomy,
Pomona College, Claremont, CA 91711, USA  Center for Computational
Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
Sam Ponnada California Institute of Technology, TAPIR, Mailcode 350-17,
Pasadena, CA 91125, USA

###### Abstract

Feedback from supermassive black holes (SMBHs) is believed to be a critical
driver of the observed quenching of star formation and color bimodality of
galaxies above the Milky Way mass scale. In recent years, various forms of
SMBH feedback have been implemented as subgrid models in galaxy formation
simulations, but most implementations have involved simplified prescriptions
or coarse-grained models of the interstellar medium (ISM). We present the
first set of FIRE-3 cosmological zoom-in simulations with AGN feedback evolved
to $z\sim 0$, examining a set of galaxies with halos in the mass range
$10^{12}-10^{13}\,{\rm M_{\odot}}$. These high-resolution simulations combine
detailed stellar and ISM physics with multi-channel AGN feedback including
radiative feedback, mechanical outflows, and in some simulations, cosmic rays.
We find that massive (>L*) galaxies in these simulations can match local
scaling relations including the stellar mass-halo mass relation, the $M_{\rm
BH}$-$\sigma$ relation, the size-mass relation, and the Faber-Jackson
relation. Many of the massive galaxies in the simulations with AGN feedback
have quenched star formation and elliptical morphologies, in qualitative
agreement with observations. In contrast, simulations at the massive end
without AGN feedback produce galaxies that are too massive and form stars at
too high rates, are order-of-magnitude too compact, and have velocity
dispersions well above Faber-Jackson. Despite these successes, the AGN physics
models analyzed do not necessarily produce uniformly realistic galaxies across
the full mass range studied when the feedback parameters are held constant,
indicating that further refinements of the black hole modeling may be
warranted.

Supermassive black holes (1663), Active galactic nuclei(16), Galaxy quenching
(2040), AGN host galaxies (2017), Galaxy evolution (594), Hydrodynamical
simulations (767)

##  1 Introduction

Over the past two decades, observations of galaxies hosting luminous quasars
have found wide-angle kiloparsec-scale outflows of molecular gas (Feruglio et
al., 2010; Fischer et al., 2010; Rupke & Veilleux, 2013; Moe et al., 2009; Liu
et al., 2013; Fabian, 2012). These outflows provide crucial evidence for the
impact of active galactic nuclei (AGN) on their host galaxies: the outflow
rates can exceed the galaxy’s star formation rate, and the velocities of the
outflows are correlated to the luminosity of the AGN (Sturm et al., 2011;
Cicone et al., 2014). Other properties of supermassive black holes (BHs) have
also been demonstrated to correlate with properties of their host galaxies,
including scaling relations between black hole mass and stellar velocity
dispersion and between black hole mass and stellar bulge mass (Ferrarese &
Merritt, 2000; Magorrian et al., 1998; Gebhardt et al., 2000). These scaling
relations indicate that black holes and galaxies develop together. This co-
evolution can likely be explained through the feedback energy produced by
accretion onto the SMBH (Silk & Rees, 1998; Di Matteo et al., 2005; Hopkins et
al., 2007). In addition to driving winds, observational evidence suggests that
AGN can heat the interstellar/circumgalactic/intracluster medium (ISM/CGM/ICM)
with powerful, relativistic radio jets (Fabian, 2012).

Feedback from supermassive black holes is thought to be key to understanding
how galaxies quench, or stop actively forming stars (Di Matteo et al., 2005;
Springel et al., 2005; Sijacki et al., 2007). Substantial observational
evidence indicates that star formation becomes increasingly likely to cease in
galaxies above a stellar mass of approximately $10^{10.5}\,{\rm M_{\odot}}$,
corresponding to a halo mass $M_{\rm halo}\sim 10^{12}\,{\rm M_{\odot}}$,
creating a population of ‘red and dead’ massive galaxies (Peng et al., 2010;
Muzzin et al., 2013; Tomczak et al., 2014). Understanding the physical
processes which lead to this quenching is a key question in galaxy formation
theory. While several possible hypotheses have been proposed, feedback from
AGN has emerged as a leading candidate, with some studies showing that without
AGN feedback, simulations fail to reproduce the observed galaxy mass function
(Somerville et al., 2008; Harrison, 2017) and massive galaxy properties (Su et
al., 2019).

As the physics of SMBH accretion and feedback occurs on scales below the
resolution of such cosmological simulations, it is necessary to rely on
subgrid models. However, much uncertainty remains as to how best to represent
the physics involved due to open questions about the fundamental physics.
Various forms of AGN feedback—including various combinations of thermal energy
injection, winds, and jets—have been implemented as sub-grid models in galaxy
formation simulations (see reviews in Somerville & Davé, 2015; Naab &
Ostriker, 2017; Crain & Voort, 2023; Di Matteo et al., 2023). Large-volume
cosmological simulations (e.g Weinberger et al., 2017, 2018; Schaye et al.,
2015; Tremmel et al., 2017; Davé et al., 2019) have had some success in
producing quenched populations of massive galaxies. These large-volume
simulations utilize a coarse-grained treatment of the ISM with baryonic mass
resolutions in the range of $m_{\rm b}\sim 10^{5}-10^{7}\,{\rm M_{\odot}}$,
and BH models consisting of idealized phenomenological prescriptions which are
fine-tuned to produce galaxy properties that match galaxy-scale observations.
High-resolution simulations of isolated massive ellipticals with AGN feedback,
such as those studied by e.g Ciotti & Ostriker (2007); Ciotti et al. (2010);
Su et al. (2019); Ciotti et al. (2022); Yoon et al. (2019), have also produced
promising results with regards to reproducing the global properties of those
galaxies, though these were idealized galaxies not run with a full
cosmological context and evolution.

Recently, zoom-in simulations with a full cosmological context and detailed
physics models, including a detailed star-formation prescription and multi-
channel stellar feedback and AGN feedback prescriptions (Choi et al., 2017,
2018), have successfully reproduced observed scaling relations in massive
galaxies in simulations with a baryonic mass resolution $m_{\rm b}\sim
10^{6}\,{\rm M_{\odot}}$. The simulations we present in this paper build on
what has been done previously in several key ways, by combining high
resolution ($m_{\rm b}\sim 10^{3}-10^{5}\,{\rm M_{\odot}}$) with a full
cosmological context, a multiphase ISM, and a stellar feedback model that has
been extensively studied and shown to produce realistic galaxy properties for
a wide range of stellar masses and redshifts.

In this paper, we use cosmological zoom-in simulations of galaxy formation
from the Feedback In Realistic Environments (FIRE) project111See the FIRE
project website at: http://fire.northwestern.edu. to investigate the impact of
multi-channel AGN feedback in relatively massive galaxies with halos in the
mass range of $10^{12}-10^{13}\,{\rm M_{\odot}}$. These correspond to Milky
Way-mass ($\sim$ L*) galaxies and intermediate-mass (>L*) elliptical galaxies.
In contrast to most large-volume simulations, the high resolution of the FIRE
simulations allows us to resolve the multi-phase nature of the ISM and more
accurately capture interactions between AGN feedback and the ISM. All
simulations are evolved using the FIRE-3 code (Hopkins et al., 2023). FIRE-3
represents an update from FIRE-2 (Hopkins et al., 2018), including updates to
the treatment of stellar evolution physics, nucleosynthetic yields, and low-
temperature ISM cooling. Additionally, we use a detailed, multi-channel AGN
feedback model which combines multiple forms of radiative feedback (e.g.
radiation pressure, photoionization, and Compton heating) with kinetic
outflows and, in some simulations, cosmic rays. This is a more complete set of
feedback mechanisms than those used in current large-scale simulations. Our
new cosmological simulations expand previous, more idealized studies of AGN
feedback by our team in analytic models (Faucher-Giguère & Quataert, 2012;
Faucher-Giguère et al., 2012; Richings & Faucher-Giguère, 2018) and idealized
simulations (Richings & Faucher-Giguère, 2018; Richings et al., 2021) of AGN
winds, as well as in simulations including the FIRE ISM physics but neglecting
the cosmological environment (Hopkins et al., 2016; Torrey et al., 2020). This
work is complementary to Mercedes-Feliz et al. (2023) and Cochrane et al.
(2023), which analyzed the impact of AGN winds during a luminous quasar phase
in a massive FIRE-2 star-forming galaxy, and to Wellons et al. (2023), which
explored the impact of various models of SMBH accretion and feedback in a
large suite of FIRE-2 simulations and found several physically plausible
models. We apply AGN feedback models broadly consistent with those favored by
Wellons et al. (2023) to a suite of cosmological zoom-in simulations with
updated FIRE-3 physics. Moreover, about half of our new simulations are of
order-of-magnitude higher mass resolution than the simulations of comparable-
mass galaxies studied in Wellons et al. (2023).

We examine the impact of AGN feedback in the FIRE simulations by making a
detailed comparison between two possible models for AGN feedback, simulations
without BHs, and observations. We present a set of massive (>L*) galaxies
which broadly match observed scaling relations at low redshifts, which along
with some of the FIRE-2 test simulations studied by Wellons et al. (2023)
constitute the first such galaxies produced by the FIRE project.

In section 2 we describe the suite of zoom-in simulations analysed for this
study, along with the three SMBH physics models tested. In section 3 we
present the outcomes of the simulations with comparisons to observations of
galaxy masses and star-formation rates (section 3.2) and structural properties
(section 3.3). In section 4 we discuss our results, and in section 5 we
summarize our conclusions.

Throughout, we assume a standard flat $\Lambda$CDM cosmology consistent with
recent measurements (Planck Collaboration et al., 2020).

##  2 Methods

Parameter | NoBH | BH | BH+CR
---|---|---|---
vwind [km s-1] | – | 3000 | 3000
$\epsilon_{\rm wind}^{\rm BH}$ | – | $5.0\times 10^{-5}$ | $5.0\times 10^{-5}$
$\dot{p}_{\rm rad}\;[L/c]$ | – | $\sim 1$ | $\sim 1$
$\epsilon_{\rm CR}^{\rm BH}$ | – | – | $10^{-3}$
Table 1: BH feedback parameters and energetics for the physics models compared in this paper. vwind is the velocity of mechanical winds. $\epsilon^{\rm BH}_{\rm wind}=(v_{\rm wind}/c)^{2}/2$ is the efficiency of mechanical winds (eq. 5). $\dot{p}_{\rm rad}$ is the momentum flux deposited by radiative feedback. $\epsilon_{\rm CR}$ is the efficiency of cosmic ray feedback (eq. 6).  Halo | Model | Final $z$ | Halo Mass [$\,{\rm M_{\odot}}$] | Stellar Mass [$\,{\rm M_{\odot}}$] |  $m_{\rm b}$ [$\,{\rm M_{\odot}}$] |  $\epsilon_{\rm star}$ [$h^{-1}$ kpc]
---|---|---|---|---|---|---
| NoBH | 0.00 | $1.19\times 10^{12}$ | $1.45\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
m12f | BH | 0.00 | $1.15\times 10^{12}$ | $1.46\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
| BH+CR | 0.00 | $1.06\times 10^{12}$ | $8.59\times 10^{9}$ | $6\times 10^{4}$ | 0.013
| NoBH | 0.00 | $1.35\times 10^{12}$ | $3.45\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
m12q | BH | 0.20 | $1.15\times 10^{12}$ | $1.44\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
| BH+CR | 0.00 | $1.19\times 10^{12}$ | $7.05\times 10^{9}$ | $6\times 10^{4}$ | 0.013
m12b | BH | 0.00 | $1.10\times 10^{12}$ | $2.09\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
m12i | BH | 0.00 | $8.52\times 10^{11}$ | $1.36\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
| BH+CR | 0.00 | $7.78\times 10^{11}$ | $9.93\times 10^{9}$ | $6\times 10^{4}$ | 0.013
m12m | BH | 0.00 | $1.27\times 10^{12}$ | $4.47\times 10^{10}$ | $7\times 10^{3}$ | 0.0066
| BH+CR | 0.00 | $9.46\times 10^{11}$ | $1.64\times 10^{10}$ | $6\times 10^{4}$ | 0.013
m12r | NoBH | 0.00 | $8.52\times 10^{11}$ | $7.94\times 10^{9}$ | $7\times 10^{3}$ | 0.0068
| BH | 0.00 | $8.87\times 10^{11}$ | $8.77\times 10^{9}$ | $7\times 10^{3}$ | 0.0068
m12w | NoBH | 0.00 | $8.54\times 10^{11}$ | $9.24\times 10^{9}$ | $7\times 10^{3}$ | 0.0068
| BH | 0.00 | $8.90\times 10^{11}$ | $1.37\times 10^{10}$ | $7\times 10^{3}$ | 0.0068
| NoBH | 0.00 | $9.10\times 10^{12}$ | $4.56\times 10^{11}$ | $3\times 10^{5}$ | 0.024
m13h113 | BH | 0.08 | $8.07\times 10^{12}$ | $1.70\times 10^{11}$ | $3\times 10^{4}$ | 0.011
| BH+CR | 0.00 | $6.90\times 10^{12}$ | $8.89\times 10^{10}$ | $3\times 10^{5}$ | 0.024
| NoBH | 0.00 | $6.72\times 10^{12}$ | $3.30\times 10^{11}$ | $3\times 10^{5}$ | 0.024
m13h206 | BH | 0.20 | $5.84\times 10^{12}$ | $1.36\times 10^{11}$ | $3\times 10^{4}$ | 0.011
| BH+CR | 0.00 | $5.65\times 10^{12}$ | $5.22\times 10^{10}$ | $3\times 10^{5}$ | 0.024
m13h002 | BH+CR | 0.00 | $1.93\times 10^{13}$ | $9.61\times 10^{10}$ | $3\times 10^{5}$ | 0.024
m13h007 | NoBH | 0.00 | $2.32\times 10^{13}$ | $6.41\times 10^{11}$ | $3\times 10^{5}$ | 0.024
| BH+CR | 0.00 | $1.95\times 10^{13}$ | $1.48\times 10^{11}$ | $3\times 10^{5}$ | 0.024
m13h009 | BH+CR | 0.00 | $2.76\times 10^{13}$ | $7.98\times 10^{10}$ | $3\times 10^{5}$ | 0.024
m13h029 | NoBH | 0.14 | $1.35\times 10^{13}$ | $6.15\times 10^{11}$ | $3\times 10^{5}$ | 0.024
| BH+CR | 0.00 | $1.08\times 10^{13}$ | $1.20\times 10^{10}$ | $3\times 10^{5}$ | 0.024
m13h037 | BH+CR | 0.00 | $2.04\times 10^{13}$ | $1.28\times 10^{11}$ | $3\times 10^{5}$ | 0.024
m13h223 | NoBH | 0.06 | $4.09\times 10^{13}$ | $9.39\times 10^{11}$ | $3\times 10^{5}$ | 0.024
m13h236 | NoBH | 0.20 | $8.94\times 10^{12}$ | $3.15\times 10^{11}$ | $3\times 10^{5}$ | 0.024
| BH+CR | 0.00 | $1.48\times 10^{13}$ | $6.68\times 10^{10}$ | $3\times 10^{5}$ | 0.024
Table 2: Parameters of the FIRE-3 simulation suite analyzed in this work.
Milky Way-mass halos (‘m12s’) are listed in the upper half of the table, and
massive galaxies (‘m13s’) are listed in the lower half. Bolded simulations are
those for which the initial conditions were run with all three primary physics
models discussed in this paper. Columns are: physics model, final redshift,
halo mass and stellar mass at final redshift, the baryonic mass resolution
$m_{\rm b}$, and the stellar gravitational softening length $\epsilon_{\rm
star}$ in physical units. Figure 1: Edge-on mock Hubble images for two example
Milky Way-mass simulated galaxies (top two rows) and two example m13 massive
galaxies (bottom two rows) at $z\sim 0$. For each set of initial conditions,
we compare a model with no BHs or CRs (left) to a model with BHs but no CRs
(middle) and to a model with both BHs and CRs (right). Milky Way-mass galaxies
tend to be disk-like at $z=0$ with the NoBH and BH models, while the massive
galaxies produce a mix of disks and spheroidals for those two models. The
BH+CR model consistently produces elliptical galaxies. Figure 2: Stellar mass-
halo mass (SMHM) relation (left) and $M_{\rm BH}-\sigma$ relation (right) for
$z\sim 0$. We compare the NoBH (black, unfilled), BH (green), and BH+CR
(purple) models. Left panel: The SMHM relation from (Behroozi et al., 2019)
and its intrinsic scatter. Right panel: The observational $M_{\rm BH}-\sigma$
relation (Greene et al., 2020) and its intrinsic scatter.

###  2.1 FIRE Simulations

We analyze a set of FIRE-3 cosmological zoom-in simulations, all of which were
run with the hydrodynamics code GIZMO222See the GIZMO project website at:
http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html. in meshless finite
mass mode (Hopkins, 2015). Details of the FIRE-3 methods and physics are
explained in Hopkins et al. (2023). The simulations include all standard
FIRE-3 physics, including multiple forms of stellar feedback, including
feedback from supernovae of Type I and II, stellar winds from OB and AGB
stars, photoionization, and radiation pressure. Star formation occurs in
dense, self-gravitating gas. The high resolution of these simulations allows
us to resolve the multiphase ISM. All simulations include magnetic fields,
using the magneto-hydrodynamics methods described in Hopkins & Raives (2016)
and Hopkins (2016).

Our simulation suite consists of sixteen sets of initial conditions: seven
Milky Way-mass galaxies with final halo mass $M_{\rm halo}\sim 10^{12}\,{\rm
M_{\odot}}$ (which we hereafter refer to as ‘m12’ galaxies), and nine galaxies
reaching $M_{\rm halo}\sim 10^{13}\,{\rm M_{\odot}}$ at $z=0$ (which we
hereafter refer to as ‘m13’ galaxies). Note that while we will refer to the
m13 galaxies as ‘massive’ in comparison to the Milky Way-mass galaxies, they
are what would generally be considered intermediate-mass relative to >L*
ellipticals. As shown in Table 2, each set of initial conditions is evolved
with between one and three physics models. The three primary models we analyze
for this paper are as follows: one without any supermassive black holes at
all, one which models BHs with both radiative and mechanical feedback, and one
with both BHs and cosmic ray physics. These models, which we hereafter refer
to as the “NoBH”, “BH”, and “BH+CR” models, are described in more detail below
and defined by the parameters in Table 1. Four sets of initial conditions, two
m12s and two m13s, were evolved with all three of these models; these galaxies
are bolded in Table 2. We tested a range of mass resolutions, from a baryonic
mass of $m_{\rm b}=7\times 10^{3}\>\,{\rm M_{\odot}}$ for our highest-
resolution simulations to a baryonic mass of $m_{\rm b}=3\times 10^{5}\>\,{\rm
M_{\odot}}$ for our lowest-resolution simulations. The final redshift is $z=0$
for most of these simulations, with a few ending at redshifts of $z\sim 0.1$
or $z\sim 0.2$; these minor differences in final redshift do not affect our
conclusions, and for simplicity, we refer to this range as $z\sim 0$
throughout our analysis.

For each simulation snapshot, we define the halo center using the “shrinking
sphere” method of Power et al. (2003) to find the iterative center of mass of
the star particles in the zoom-in region. Throughout this work, we consider
the virial radius $R_{\rm vir}$ to be $R_{200c}$, the radius which encloses a
sphere whose average density is 200 times the critical density. We consider
the halo mass $M_{\rm halo}$ to be $M_{200c}$, the total mass within that
region.

###  2.2 Black Hole Physics

In simulations which include BHs, the BH seeding, dynamics, accretion, and
feedback follow the methods described in Hopkins et al. (2023). Black holes
with seed mass $M_{\rm seed}=100\,{\rm M_{\odot}}$ are probabilistically
generated from gas cells and form preferentially in areas of high surface
density and low metallicity. On average, one BH will seed per $M_{0}$ stellar
mass, with $M_{0}=(1-3)\times 10^{6}\,{\rm M_{\odot}}$. They move under the
influence of gravity, and following the method in Wellons et al. (2023), are
artificially accelerated towards the local particle with the highest binding
energy, which encourages them to move towards the center of the galaxy. If two
black hole particles become gravitationally bound to one another within the
same interaction kernel, and have overlapping force softenings, they merge
within a single time-step.

Each black hole, with mass $M_{\rm BH}$ is surrounded by an accretion disk of
mass $M_{\rm acc}$. Gas from the accretion kernel, which includes the 256
nearest gas resolution elements up to a maximum radius of 5 kpc, accretes onto
the disk at a rate

$\dot{M}_{\rm acc}=\eta_{\rm acc}(1-f_{w,*})\,M_{\rm gas}\,\Omega$ (1)

where the efficiency parameter

$\eta_{\rm acc}\equiv\frac{0.01[(M_{\rm BH}+M_{\rm acc})/M_{\rm
disk}]^{1/6}}{1+3(M_{\rm gas}/M_{\rm disk})(M_{\rm disk}/10^{9}\,{\rm
M_{\odot}})^{1/3}}$ (2)

is calibrated from high-resolution simulations to represent the effects of
gravitational torques, which lead to angular momentum loss for gas in the
galaxy’s nucleus and drive gas towards the BH (Hopkins & Quataert, 2011;
Anglés-Alcázar et al., 2017b, a; Anglés-Alcázar et al., 2021). The gas mass
$M_{\rm gas}$, disk mass $M_{\rm disk}$, and dynamical frequency
$\Omega=\sqrt{GM_{\rm enc,tot}}/R^{3}$ are evaluated within the radius of the
accretion kernel. The total gas and stellar disk mass, $M_{\rm disk}$, is
estimated following Hopkins & Quataert (2011), decomposing the stars in the
accretion kernel into isotropic and thin disk components based on their
angular momentum, and assuming cold gas forms a disk. The parameter

$f_{w,*}\equiv\left[1+\frac{\bar{a}_{\rm
grav}}{\langle\dot{p}_{*}/m_{*}\rangle}\right]^{-1}$ (3)

is a correction for gas lost by ejection from stellar feedback (Hopkins et
al., 2022b), where $\bar{a}_{\rm grav}$ is the gravitational acceleration
inwards and $\langle\dot{p}_{*}/m_{*}\rangle\sim 10^{-7}$ cm s-2 is the the
momentum flux per unit mass from stellar feedback.

Accreting BHs inject energy into their surroundings in the forms of both
radiative and mechanical feedback. Radiative feedback includes radiation
pressure, photo-ionization, and Compton heating. Radiation transport is
modelled using the same numerical method as is used for stellar radiative
feedback in the FIRE simulations (Hopkins et al., 2020). The accretion disk
emits radiation with a bolometric luminosity of

$\dot{E}^{\rm BH}_{\rm rad}\equiv L_{\rm bol}=\epsilon^{\rm BH}_{\rm
r}\,\dot{M}_{\rm BH}c^{2}$ (4)

where $\epsilon^{\rm BH}_{\rm r}=0.1$, and the total photon momentum flux is
$\dot{p}_{\rm rad}=L_{\rm abs}/c$ where $L_{\rm abs}$ is the absorbed photon
luminosity in a given gas element. The total momentum deposited by radiation
in the ISM in these simulations is typically $\sim L_{\rm bol}/c$ because
single scattering dominates.

Figure 3: Star formation histories and BH accretion histories for the same
four example simulated galaxies shown in Figure 1: two Milky Way-mass galaxies
(upper) and two massive galaxies (lower). We compare three different physics
models: NoBH (black), BH (green), and BH+CR (purple). The inclusion of BH
feedback reduces the star formation rate at later times compared to the SFRs
in simulations without BH feedback. This effect is particularly strong for the
BH+CR model in simulations of massive (m13) galaxies, in which the SFR is
consistently lower than for the other models starting early in the galaxies’
lifetimes. In addition to having lower SFRs than the BH model, the BH+CR model
also produces lower BH accretion rates.

Mechanical winds are modeled by the spawning of new resolution elements, or
wind cells, in the vicinity of the black hole (see also Richings & Faucher-
Giguère, 2018; Torrey et al., 2020; Su et al., 2021; Cochrane et al., 2023;
Mercedes-Feliz et al., 2023). These cells, which have masses of $100\,{\rm
M_{\odot}}$ (for m12 halos) or $400\,{\rm M_{\odot}}$ (for m13 halos), are
launched with a mass-loading $\dot{M}_{\rm wind}={\dot{M}}_{\rm BH}$ and a
velocity of $v_{\rm wind}=$ 3000 km s-1, generating kinetic energy

$\dot{E}^{\rm BH}_{\rm wind}=\dot{M}_{\rm BH}\,v_{\rm
wind}^{2}/2=\epsilon^{\rm BH}_{\rm wind}\,\dot{M}_{\rm BH}c^{2}$ (5)

where $\epsilon^{\rm BH}_{\rm wind}=(v_{\rm wind}/c)^{2}/2\approx 5\times
10^{-5}$. The outflows are initially jet-like: collimated along the angular
momentum axis of the BHs at the point of injection, though after meeting the
ISM and shock heating, they tend to spread along paths of least resistance,
with wide solid-angle coverage.

We also test the impact of cosmic ray feedback. The NoBH and BH models do not
include any cosmic rays. In the “BH+CR” model, cosmic rays are injected from
both stellar and black hole sources. The CRs are spectrally resolved for
protons and electrons, and 10% of the initial kinetic energy from all stellar
mechanical feedback processes (e.g. SNe and fast stellar winds) is injected in
the form of CRs. Additionally, whenever mechanical outflows are injected from
BHs, CRs are also injected with an efficiency of

$\dot{E}^{\rm BH}_{\rm CR}=\epsilon^{\rm BH}_{\rm CR}\,\dot{M}_{\rm BH}c^{2}$
(6)

where $\epsilon^{\rm BH}_{\rm CR}=10^{-3}$. For the CR scattering rate, we
adopt the modified external driving model described in §5.3.2 of Hopkins et
al. (2022a), which is calibrated to reproduce the Voyager and AMS-02
observations from MeV to TeV energies in Milky Way-like simulations. In this
model, the CR scattering rate depends on local plasma properties. Under local
ISM conditions, this model produces similar scattering rates to common
phenomenological models.

##  3 Results

###  3.1 Morphologies

The inclusion of BH feedback affects the morphology of the simulated galaxies.
Using the open-source visualization software FIREstudio (Gurvich et al.,
2022), we generated edge-on mock Hubble renderings of the stellar
distributions of two Milky Way-mass simulated galaxy, m12q and m12f, and two
massive (>L*) galaxies, m13h113 and m13h206, for each of the three models
(Figure 1). Milky Way-mass galaxies produce disk-like structures at $z\sim 0$
in simulations without BHs, and in simulations with BHs but no cosmic rays,
while the massive galaxies produce a mix of disks and spheroidals for those
two models. The addition of cosmic rays in simulations with BHs produces
galaxies which appear spherical regardless of mass. The simulated galaxies
with cosmic rays appear smaller than their counterparts with the same initial
conditions, most likely as a result of the reduced stellar mass caused by that
model as demonstrated in section 3.2.

###  3.2 Masses and Star-Formation Rates

Simulations of massive galaxies with AGN feedback can reproduce basic scaling
relations at $z\sim 0$. Figure 2 shows how the simulations compare to the
stellar mass-halo mass and $M_{\rm BH}-\sigma$ relations at $z\sim 0$. Points
are colored by physical model: unfilled black points represent simulations
without BHs (the NoBH model), green points represent simulations with BHs but
without cosmic rays (the BH model), and purple points represent simulations
with both BHs and cosmic ray physics (the BH+CR model).

The left panel shows the stellar mass-halo mass (SMHM) relation. Stellar
masses for the simulated galaxies are calculated as all stars within
$0.1R_{\rm vir}$ of the galaxies’ centers, and are compared to the relation
presented in Behroozi et al. (2019). At $z\sim 0$, massive galaxies with no
AGN feedback have become increasingly overmassive, indicating under-suppressed
star formation. By contrast, galaxies with AGN feedback lie much closer to the
relation. The two massive galaxies simulated with the BH model lie slightly
(around $1\sigma$) above the relation. Massive galaxies with the BH+CR model
lie either along the relation or slightly below it; six of the eight galaxies
in this subsample are within $1\sigma$, and the remaining two are just below
$1\sigma$ beneath the relation. The inclusion of cosmic rays as another AGN
feedback channel therefore appears to have a significant effect on regulating
the stellar mass. However, a few of the Milky Way-mass galaxies with the BH+CR
model are significantly undermassive relative to observations, indicating
possibly over-suppressed star formation, in contrast to the other two models
which more closely reproduce the observed relation at this mass scale. As this
model successfully reproduces the SMHM relation for the massive galaxies, this
result may suggest that fixed AGN feedback efficiencies like the ones in our
model are unable to reproduce the differential quenching at different mass
scales implied by observations.

The relationship between stellar mass and stellar velocity dispersion ($M_{\rm
BH}-\sigma$) is another key scaling relation, and both BH models (with and
without CRs) reproduce it well. The right panel shows the $M_{\rm BH}-\sigma$
relation for the two models which included black holes, compared to the
observational relation from Greene et al. (2020). The stellar velocity
dispersion $\sigma$ is calculated as the standard deviation of the stellar
velocities along the stellar angular momentum axis, for stars within the
projected half-mass radius. Almost all of the simulations remain within the
intrinsic scatter of the observed relation at $z\sim 0$; only one simulation,
an m12 galaxy with the BH model, falls very slightly below it.

In these simulations, BH feedback appears to regulate star formation and self-
regulate BH growth. In Figure 3, we examine the star formation histories and
black hole accretion rates of the four simulated galaxies for which we can
compare all three primary models. The star-formation histories are constructed
using an ‘archaeological’ method (e.g. as in Flores Velázquez et al., 2021;
Gurvich et al., 2023) in which we take all stars within 0.1 $R_{\rm vir}$ at
the final snapshot and estimate their formation masses by accounting for
stellar mass-loss according to rates from STARBURST99 (Leitherer et al.,
1999). We then reconstruct the time-dependent SFH from the distribution of
star particle formation times (SFTs) and star formation masses. Without BH
physics included, the star formation rates remain steady and relatively
constant over the course of the galaxy’s lifetime. The inclusion of AGN
feedback leads to reductions in the star formation rates after a few Gyr of
the galaxies’ lifetimes: stronger feedback (i.e. the BH+CR model) produces
lower star-formation rates at late times than the BH model, and both produce
lower rates than the NoBH model. The BH accretion rates appear to follow a
similar pattern: stronger feedback in the form of our BH+CR model leads to
reduced BH accretion rates, relative to the model without cosmic rays. At
times, the BH accretion rates appear to track the SF rates: remaining steady
when the SF is steady, and decreasing at approximately the same time as major
SF quenching episodes. While not every feature of the star formation histories
appears in the BH accretion histories, this general correlation is indicative
that both stars and BHs in the galactic nuclei are drawing on a common gas
supply—and BH feedback reducing that gas supply can regulate both star-
formation and the growth of the BH itself.

Figure 4: Specific star formation rate vs. stellar mass of the galaxy at
$z\sim 0$. We compare the NoBH (black, unfilled), BH (green), and BH+CR
(purple) models. Triangles indicate upper limits. Observed specific star
formation rates from Salim et al. (2007) are shown for the star-forming
sequence (dotted red line) along with the intrinsic observed scatter (red
shaded region). We indicate the stellar mass at which 50% of galaxies are
expected to be quenched with a black dashed line (Behroozi et al., 2019); the
grey shaded region represents the inter-quartile range of the quenched
fraction from 25%-75%. The NoBH model produces galaxies which are consistently
star-forming, while the BH+CR model produces galaxies which are consistely
quenched, and the BH model produces a mixture of quenched and star-forming
galaxies.  Figure 5: Projected stellar half-mass radius versus stellar mass
for galaxies at redshift $z=0$, compared to the observations of projected
half-light radii found by Baldry et al. (2012). We compare the NoBH (black,
unfilled), BH (green), and BH+CR (purple) models. Massive galaxies with the
NoBH model are consistently overly compact, while all but one of the
simulations with the BH+CR model are consistent with observed early-type
galaxies. Massive galaxies with the BH model, without cosmic rays, are also
overly compact, though not to the same degree as the NoBH galaxies.

Figure 4 shows the specific star formation rate for each simulation at $z\sim
0$, averaged over 10 Myr, along with observed specific star formation rates
for the star-forming sequence from Salim et al. (2007). We indicate the
stellar mass at which 50% of galaxies are expected to be quenched at $z=0$
(Behroozi et al., 2019), along with the inter-quartile range of the quenched
fraction from 25%-75%. Star formation rates for massive m13 galaxies with no
AGN feedback, or for the BH model without cosmic rays, remain at or above the
star-forming sequence. By contrast, all galaxies with the BH+CR model are well
below the star-forming sequences. Many of these massive galaxies with AGN
feedback are quenched entirely, while a few retain small amounts of star
formation, consistent with observations. For Milky Way-mass m12 galaxies,
however, the stronger feedback model (BH+CR) appears to quench all galaxies,
while the BH and NoBH models produce a more realistic mixture of star-forming
and quenched galaxies.

###  3.3 Structural Properties

Simulations without AGN have struggled to produce realistic galaxy sizes at
the massive end, as without AGN feedback, massive galaxies can become overly
compact by an order of magnitude (e.g. Choi et al., 2017, 2018; Wellons et
al., 2020; Parsotan et al., 2021). We consider the impact of AGN feedback on
galaxy sizes in our simulations. Figure 5 shows the relationship between the
stellar mass and the stellar half-mass radius, which is calculated as the
radius at which the cumulative stellar mass profile reaches half the total
stellar mass, taking the average after projecting along five randomly-selected
angles. Galaxy sizes are compared to the observational $R_{e}-M_{*}$ relations
found by Baldry et al. (2012) and Shen et al. (2003), assuming a constant
mass-to-light ratio, which is a reasonable approximation for old stellar
populations. Simulated massive galaxies without AGN feedback (the NoBH model)
are too compact by an order of magnitude; given their high stellar masses,
this indicates that these galaxies are dramatically more compact than observed
galaxies. By contrast, simulations with the BH+CR model are almost all
consistent with observations: all four m12s, and seven out of eight m13s, are
within the scatter of the observed relation. The BH+CR model performs better
by this metric than the BH model, as several m12 galaxies and both m13
galaxies with the latter model are also compact relative to observations.

Another important scaling relation is the relationship between stellar
velocity dispersion $\sigma$ and stellar mass—the Faber-Jackson relation
(Faber & Jackson, 1976). Figure 6 shows this relation for all simulations at
$z\sim 0$. While all but a few of the simulated Milky Way-mass m12 galaxies,
both with and without AGN feedback, fall within the intrinsic scatter of the
observed relation found by Gallazzi et al. (2006), the inclusion of AGN
feedback appears to have a profound impact on massive galaxies. All but one of
the simulated massive galaxies without AGN feedback have overly high velocity
dispersions, 0.3-0.5 dex above the intrinsic scatter of the observations. By
contrast, simulated massive galaxies with the BH+CR model lie on or near the
relation, and massive galaxies with the BH model lay slightly (within < 0.2
dex) above it. The inclusion of BH feedback therefore appears to be key to
successfully reproducing this relation for massive galaxies: the lower stellar
masses and larger sizes of the galaxies with BH feedback reduce the velocity
dispersions to more realistic levels, since $\sigma\varpropto\sqrt{M/R}$.

Figure 6: Stellar velocity dispersion versus stellar mass for galaxies at
redshift $z=0$, compared to the observed relation for early-type galaxies
found by Gallazzi et al. (2006). We compare the NoBH (black, unfilled), BH
(green), and BH+CR (purple) models. At the massive end, NoBH simulations have
overly high velocity dispersions, while massive galaxies with the BH+CR model
lie on or near the observed relation, and massive galaxies with the BH model
lie slightly above it. For Milky Way-mass galaxies, all galaxies produced with
the NoBH model and all but one with the BH model lie near the relation, while
two of the four galaxies with the BH+CR model lie slightly below it. Figure 7:
Stellar surface density profiles for massive galaxies with the NoBH (black,
unfilled), BH (green), and BH+CR (purple) models. Stellar surface densities
are calculated in radial annuli and shown as a function of $r^{1/4}$. Observed
profiles of local massive ellipticals (grey) are reproduced from Hopkins et
al. (2010), with observations from Lauer et al. (2007) separated into ‘cusp’
and ‘core’ systems. The shaded regions represent the $1\sigma$ scatter for the
observed sample.The NoBH galaxies are comparable in total stellar mass to the
observed ‘core’ galaxies, while the BH+CR galaxies are comparable in mass to
the ‘cusp’ galaxies, and the BH galaxies lie somewhere in between. Simulations
without BH feedback are overdense by up to two orders of magnitude in the
inner regions, while introducing BH feedback reduces stellar densities to more
realistic levels. Without cosmic rays, simulated galaxies with BH feedback are
overly dense in the central regions, but converge to observations by 1 kpc;
with cosmic rays, the central stellar surface densities are further reduced.
Figure 8: Mass-weighted normalized distributions of the circularity
$\epsilon=j_{z}/j_{c}(E)$ of star particles in each simulation. Each line
represents an individual simulation. Histograms are categorized by mass
scale—MW-mass m12s (top) and massive m13s (bottom)—and by model—NoBH (left),
BH (center), and BH+CR (right). Distributions with peaks around $\epsilon=0$
indicate elliptical morphologies, while peaks at $\epsilon=1$ indicate disk-
like morphologies.

In Figure 7, we examine the stellar surface density profiles for the massive
galaxies in our simulated sample. Projected stellar surface densities are
calculated in radial annuli and shown as a function of $r^{1/4}$ (de
Vaucouleurs et al., 1948). Observed profiles of local massive ellipticals
(grey) are reproduced from Hopkins et al. (2010), with observations from Lauer
et al. (2007) separated into ‘cusp’ and ‘core’ systems. Simulations without BH
feedback are overdense by up to two orders of magnitude in the inner regions.
Introducing BH feedback reduces stellar densities to more realistic levels.
Without cosmic rays, simulated galaxies with BH feedback are still overdense
in the central regions, but converge to observation at by 1 kpc. With cosmic
rays, the stellar surface densities are further reduced, more closely matching
observations in the inner regions of the galaxies. Stellar densities may be
artificially flat within a few times the stellar softening lengths listed in
Table 2.

We examine the effects of AGN feedback on the kinematic characteristics of
each galaxy in Figure 8, which shows the the mass-weighted normalized
distributions of the orbital circularity $\epsilon$ of every star particle
within $0.1R_{\rm vir}$ in each galaxy. Circularity $\epsilon=j_{z}/j_{c}(E)$
is defined as the ratio of each particle’s angular momentum in the $\hat{z}$
direction to that of a circular orbit with the same energy (Abadi et al.,
2002). Following Scannapieco et al. (2008), we define
$j_{c}(E)=rv_{c}=r\sqrt{GM_{\rm enc}/r}$ for a particle at radius $r$, a
reasonable approximation for galaxies with spherically symmetric potentials. A
spheroidal galaxy will have a circularity distribution centered around
$\epsilon=0$, while for a rotation-supported disk, the distribution will peak
at $\epsilon=1$. Galaxies with both a disk and a bulge will have multiple
peaks. We find that simulations without BHs, or with BHs but no CRs, tend to
produce a mixture of morphologies, including many galaxies with a significant
disk component. Simulations with both BHs and CRs, however, consistently
produce spheroidal galaxies with distributions peaking at $\epsilon=0$, as is
also visually suggested in Figure 1.

##  4 Discussion

###  4.1 Scaling Relations and the Quenching of Star Formation

We demonstrate that cosmological simulations with detailed stellar and ISM
physics combined with multi-channel AGN feedback can successfully produce
massive galaxies with little or no star formation, broadly in agreement to
observations which show galaxy populations dominated by quenched galaxies at
stellar masses of $M_{*}\gtrsim 10^{10.5}\,{\rm M_{\odot}}$, and in contrast
to simulations without AGN feedback which do not quench (Figures 3 and 4).
This provides further evidence that BH feedback is necessary to produce a
realistic bimodal population of galaxies. The reduction in star formation for
galaxies at the massive end when AGN feedback is included is in broad
agreement with many other simulations (e.g. Di Matteo et al., 2005; Dubois et
al., 2012; Schaye et al., 2015; Weinberger et al., 2017, 2018; Wellons et al.,
2023; Mercedes-Feliz et al., 2023; Weinberger et al., 2023).

Our simulations of massive galaxies with AGN feedback can successfully match
local scaling relations including the stellar mass-halo mass relation (Figure
2, left), the $M_{\rm BH}$-$\sigma$ relation (Figure 2, right), the size-mass
relation (Figure 5, and the Faber-Jacksonn relation ($M_{*}-\sigma$, Figure
6). Previous simulations also have demonstrated that in the absence of AGN
feedback, simulations of massive galaxies can become too compact and overly
dense, falling below the observed size-mass relation unless specifically
calibrated to reproduce it (Wellons et al., 2020; Parsotan et al., 2021). Our
results (see Figures 5 and 7) are in agreement with other simulations which
have found that the inclusion of AGN feedback can increase the size of massive
galaxies and reduce their central stellar densities, producing more realistic
galaxy sizes (Dubois et al., 2016; Choi et al., 2018; Cochrane et al., 2023).

###  4.2 The Role of Cosmic Rays

We find that AGN feedback with cosmic rays is more efficient at quenching star
formation than AGN feedback with only radiative and kinetic feedback channels.
In the simulations presented in this paper, while both BH feedback models
produce massive galaxies which better match observed properties than the model
without BHs, we found that the BH+CR model was more successful than the BH
model at reproducing the size-mass and $M_{*}-\sigma$ relations (Figures 5 and
6). Galaxies simulated with the BH+CR model also display more drastic
reductions in their star formation rates than galaxies simulated with the BH
model. However, the BH+CR model may be too strong—it over-quenches Milky Way-
mass galaxies, and consistently produces very spherical ellipticals across the
mass range studied in this paper (Figures 1 and 4). Additionally, substantial
uncertainty remains as to how best to model CR transport. In future work, it
will therefore be important to explore the effect of alternate CR models on
our results.

These results are consistent with previous studies of FIRE simulations which
found that AGN feedback in the form of CRs could have a significant impact on
galaxy quenching. Wellons et al. (2023) found that the inclusion of CRs
significantly increases the potential of a BH to suppress star formation even
when the energy in CRs is relatively low. Su et al. (2020) found that CR jets
suppress cooling flows more effectively than thermal heating or momentum
injection, while Su et al. (2021) found that CR-dominated jets are able to
quench galaxies with order-of-magnitude lower energetics than jets dominated
by thermal or kinetic energy. It is worth noting, however, that the
simulations in those works used a CR model with constant diffusion, while we
use a model in which CR transport depends on local plasma properties (see
section 2.2 for more details). Sultan et al. (in prep) has found that, in the
set of simulations used in this work, CR pressure is sub-dominant in the
circumgalactic medium (CGM), in contrast to the constant diffusion model which
can produce CGMs dominated by CRs at the L* mass scale (e.g., Ji et al., 2020,
2021; Butsky et al., 2023).

As the BH+CR model includes cosmic rays from both stellar and BH feedback, it
is also important to determine whether these differences are indeed driven by
BH feedback, or if stellar cosmic rays could be the cause. To disentangle the
effects of stellar cosmic rays from BH cosmic ray feedback, we tested a fourth
physics model with cosmic rays from stellar feedback processes but without
black holes (hereafter referred to as the “NoBH+CR” model) on one m12 and one
m13 galaxy each. The m12 galaxy was integrated down to a redshift of $z=0$,
while we ended the simulation of the m13 galaxy with the NoBH+CRmodel at
$z=0.5$ after it became overly compact, due to computational expense. The m13
galaxy with the NoBH+CRmodel showed similar properties as the m13 galaxies
with the NoBH model: it failed to reproduce the observed SMHM, size-mass, and
$M_{*}-\sigma$ relations due to being overmassive, too compact, and having
overly high stellar velocity dispersions. It also was continuously star-
forming until the final snapshot at $z=0.5$, several gigayears after all
simulations with the BH+CR model had already quenched. For these reasons, the
success of the BH+CR model appears to be driven specifically by cosmic rays
from BH feedback, not cosmic rays from stellar sources.

It is worth noting that the total amounts of cosmic ray energy from stellar
feedback and BH feedback are comparable in the simulations when integrated
over the age of the universe. To see this, note that there is about one
supernova per 100 M⊙ of stars formed, so the integrated energy in CRs from
stars is $E_{\rm CR}^{\star}\approx 0.1\times 10^{51}~{}{\rm
erg}~{}(M_{\star}/{\rm 100~{}M_{\odot}})=10^{59}~{}M_{\star,11}~{}{\rm erg}$,
where $M_{\star,11}\equiv M_{\star}/10^{11}~{}{\rm M_{\odot}}$. Similarly, the
integrated energy in CRs from the BH is $E_{\rm CR}^{\rm BH}=\epsilon_{\rm
CR}^{\rm BH}M_{\rm BH}c^{2}\approx 1.8\times 10^{59}~{}{\rm erg}~{}M_{\rm
BH,8}\epsilon_{\rm CR,-3}^{\rm BH}$, where $M_{\rm BH,8}\equiv M_{\rm
BH}/10^{8}~{}{\rm M_{\odot}}$ and $\epsilon_{\rm CR,-3}^{\rm
BH}\equiv\epsilon_{\rm CR}^{\rm BH}/10^{-3}$. This suggests that total
feedback energy alone does not determine the impact on a galaxy’s star
formation. The lack of a CR-dominated halo found by Sultan et al. (in prep)
for these simulations indicates that our CR transport model produces high-
energy CR outbursts which can escape from the halo, rather than the smooth CR
injection into the halo seen in simulations with a constant diffusion model.
These results are therefore consistent with the findings of Wellons et al.
(2023), who found that the responsiveness of the AGN feedback model—that is,
how quickly and powerfully the BH will respond to changes in its environment,
for example by creating bursts of feedback when gas in the accretion kernel is
abundant—plays an important role in determining whether galaxies quench. This
suggests that CRs from AGN feedback may be more efficient at suppressing star
formation than CRs from stellar feedback because the AGN duty cycle produces
strong outbursts of CRs, while the CR injection from stellar sources is more
time-constant.

###  4.3 Evidence for a mass dependence of AGN feedback

In the simulations presented in this paper, the frequency with which galaxies
quenched is plausibly realistic for massive galaxies but perhaps too high for
Milky Way-mass galaxies (Figure 4). The reduced star formation in MW-mass
galaxies leads to final stellar masses which are in some cases significantly
below the observed SMHM relation. These results are broadly consistent with
other studies of both hydrodynamic and semi-analytic models which have found
that it is necessary for AGN feedback to become more efficient above a
critical mass scale of $\sim L*$ for simulations to reproduce the observed
bimodality (e.g. Croton et al., 2006; Nelson et al., 2018).

Rather than enforcing an explicit mass scale for AGN activity or star-
formation quenching, we chose to use fixed AGN feedback efficiencies. Our
simulations capture the multiphase ISM and CGM at high resolution, and self-
consistently resolve the formation of hot halos and the transition from bursty
to steady star formation, both of which appear to occur at approximately the
same galaxy mass as the critical mass scale for AGN feedback (Bower et al.,
2017; Anglés-Alcázar et al., 2017b; Habouzit et al., 2017; Stern et al., 2021;
Byrne et al., 2023). The lack of an explicitly-set mass scale in our models
allows us to investigate whether the mass scales set by these processes alone
could be responsible for the transition from inefficient to efficient AGN
feedback. Our findings suggest that this is not the case: resolving these
processes self-consistently does not appear to be sufficient to produce the
critical mass scale above which AGN feedback is effective at quenching star
formation, and below which it is not.

However, the subgrid treatment of BH accretion remains a source of
uncertainty, with only recent cosmological hyper-refinement simulations able
to explicitly model gas inflow rates down to scales $<0.1$ pc (Anglés-Alcázar
et al., 2021). In addition, our BH-feedback model is mass-invariant, in the
sense that energy proportional to accretion rate is being injected
independently of the Eddington ratio, the spin of the BH, the magnetic field
properties, etc. We have not yet determined to what extent these types of
small-scale accretion-flow physics may affect the large-scale properties of
AGN feedback in our simulations. In future work, we plan to update the
accretion model to incorporate insights from smaller-scale simulations about
these physical processes.

##  5 Summary and Conclusions

We present the first set of FIRE-3 simulations of massive galaxies evolved
with multi-channel AGN feedback. These high-resolution cosmological
simulations combine detailed stellar physics and a resolved multiphase ISM
with multi-channel AGN feedback in the forms of radiative feedback (photo-
electric, photo-ionization, Lyman-Werner, and Compton heating mechanisms, and
radiation pressure), mechanical outflows, and cosmic rays. The combination of
a highly-resolved ISM with detailed multi-channel stellar and AGN feedback
processes distinguishes this work from other large-volume simulations.

We simulate galaxies with halos in the mass range of $10^{12}-10^{13}\,{\rm
M_{\odot}}$ down to low redshift, and demonstrate that massive (>L*) FIRE
galaxies with realistic AGN feedback can produce reasonable matches to many
observations. We test two models of AGN feedback: the BH model, which includes
both radiative and kinetic feedback channels, and the BH+CR model, which
additionally includes feedback in the form of cosmic rays. These models are
both more successful at reproducing observations of massive galaxies than
simulations without AGN (the NoBH model), which can fail to match observed
properties such as galaxy sizes by an order of magnitude. Our main results are
as follows:

  1. 1.

The addition of BH feedback affects the structure of the galaxies. For the
same set of initial conditions, a galaxy that looks like a disk in the absence
of BH feedback can become spheroidal when feedback is included. Strong
feedback, such as our model in which the AGN feedback includes cosmic rays,
appears to consistently produce spheroidal galaxies, while the weaker feedback
model produced a mixture of disk and spheroidal galaxies. We demonstrate this
through both visual inspection of the stellar distributions, and by examining
the distributions of the orbital circularities of star particles.

  2. 2.

Massive galaxies with AGN feedback can reproduce several observed low-redshift
scaling relations, including the observed stellar mass-halo mass and $M_{\rm
BH}-\sigma$ relations. Without AGN feedback, simulated galaxies have stellar
masses which are consistently too high relative to their halo masses by a
factor of 3–5.

  3. 3.

The presence of AGN feedback can quench star formation in simulated massive
galaxies. We demonstrate that massive galaxies in simulations without BH
physics consistently remain actively star forming regardless of their mass:
their star-formation rates are roughly constant throughout their lifetimes,
and their specific star-formation rates remain at or above the star-forming
sequence as of $z\sim 0$. By contrast, simulations including BH feedback
produce quenched galaxies. Comparing the star formation histories of galaxies
with the same initial conditions evolved with each of our models, we find that
the inclusion of AGN feedback leads to reductions in the star formation rates
after a few Gyr of the galaxies’ lifetimes: stronger feedback (i.e. the BH+CR
model) produces lower star-formation rates at late times than the BH model,
and both produce lower rates than the NoBH model. By redshift $z\sim 0$, most
massive galaxies with AGN feedback are quenched.

  4. 4.

AGN feedback also produces galaxies with more realistic structural properties.
We demonstrate that the inclusion of AGN feedback can increase galaxy half-
mass radii to be consistent with the observed size-mass relation, reduce
stellar velocity dispersions to levels more consistent with observations of
the Faber-Jackson ($M_{*}-\sigma$) relation, and reduce stellar surface
densities in the central regions of the galaxies. By contrast, simulations
without AGN feedback produce galaxies which are inconsistent with observed
properties by up to an order of magnitude. The weaker BH model without cosmic
rays is somewhat less successful at reproducing some of these relations for
massive galaxies, though it still represents a significant improvement
relative to the model with no BHs at all.

With the addition of multi-channel AGN feedback to the FIRE-3 model, we are
able to produce realistic massive quenched galaxies. There are several
interesting directions for future work. First, the sample size of the set of
simulations studied in this work is relatively small. A future analysis
involving a larger sample would allow for more statistically rigorous
comparisons to observations. Secondly, the BH accretion and feedback models
used in this work are mass-invariant, and our results suggest that these
models do not accurately capture the differential quenching rates for galaxies
at different mass scales, even with a highly-resolved ISM and self-
consistently resolved hot halo formation and star formation burstiness. In
future work, we plan to explore additional physics that could potentially
introduce a mass dependence on the global effect of AGN feedback. For example,
the Eddington ratio, BH spin, or magnetic flux could affect the accretion disk
properties and therefore the intrinsic AGN feedback efficiencies, and we plan
to investigate these effects.

## Acknowledgements

LB was supported by the DOE Computer Science Graduate Fellowship through grant
DE-SC0020347. CAFG was supported by NSF through grants AST-2108230,
AST-2307327, and CAREER award AST-1652522; by NASA through grants
17-ATP17-0067 and 21-ATP21-0036; by STScI through grant HST-GO-16730.016-A;
and by CXO through grant TM2-23005X. Support for PFH and SP was provided by
NSF Research Grants 1911233, 20009234, 2108318, NSF CAREER grant 1455342, NASA
grants 80NSSC18K0562, HST-AR-15800. DAA acknowledges support by NSF grants
AST-2009687 and AST-2108944, CXO grant TM2-23006X, JWST grant GO-01712.009-A,
Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-
CSA-2023-028 by the Research Corporation for Science Advancement. JM is funded
by the Hirsch Foundation. N.A.W. was supported by a CIERA Postdoctoral
Fellowship. Numerical calculations were run on the Caltech computer cluster
Wheeler, the Northwestern computer cluster Quest, Frontera allocation FTA-
Hopkins/AST20016 supported by the NSF and TACC, XSEDE allocations ACI-1548562,
TGAST140023, and TG-AST140064, and NASA HEC allocations SMD-16-7561,
SMD-17-1204, and SMD-16-7592.

## Data Availability

The data supporting the plots within this article are available on reasonable
request to the corresponding author. A public version of the GIZMO code is
available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.
Additional data including simulation snapshots, initial conditions, and
derived data products are available at http://fire.northwestern.edu/data/.

## References

  * Abadi et al. (2002) Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2002, The Astrophysical Journal, Volume 597, Issue 1, pp. 21-34., 597, 21, doi: 10.1086/378316
  * Anglés-Alcázar et al. (2017a) Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F. 2017a, Monthly Notices of the Royal Astronomical Society, 464, 2840, doi: 10.1093/mnras/stw2565
  * Anglés-Alcázar et al. (2017b) Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017b, Monthly Notices of the Royal Astronomical Society: Letters, 472, L109, doi: 10.1093/mnrasl/slx161
  * Anglés-Alcázar et al. (2021) Anglés-Alcázar, D., Quataert, E., Hopkins, P. F., et al. 2021, The Astrophysical Journal, 917, 53, doi: 10.3847/1538-4357/AC09E8
  * Baldry et al. (2012) Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, Monthly Notices of the Royal Astronomical Society, 421, 621, doi: 10.1111/J.1365-2966.2012.20340.X/2/M_MNRAS0421-0621-MU37.GIF
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, Monthly Notices of the Royal Astronomical Society, 488, 3143, doi: 10.1093/MNRAS/STZ1182
  * Bower et al. (2017) Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, Monthly Notices of the Royal Astronomical Society, 465, 32, doi: 10.1093/mnras/stw2735
  * Butsky et al. (2023) Butsky, I. S., Nakum, S., Ponnada, S. B., et al. 2023, MNRAS, 521, 2477, doi: 10.1093/mnras/stad671
  * Byrne et al. (2023) Byrne, L., Faucher-Giguère, C. A., Stern, J., et al. 2023, Monthly Notices of the Royal Astronomical Society, 520, 722, doi: 10.1093/MNRAS/STAD171
  * Choi et al. (2017) Choi, E., Ostriker, J. P., Naab, T., et al. 2017, The Astrophysical Journal, 844, 31, doi: 10.3847/1538-4357/aa7849
  * Choi et al. (2018) Choi, E., Somerville, R. S., Ostriker, J. P., Naab, T., & Hirschmann, M. 2018, The Astrophysical Journal, 866, 91, doi: 10.3847/1538-4357/AAE076
  * Cicone et al. (2014) Cicone, C., Maiolino, R., Sturm, E., et al. 2014, Astronomy & Astrophysics, 562, A21, doi: 10.1051/0004-6361/201322464
  * Ciotti & Ostriker (2007) Ciotti, L., & Ostriker, J. P. 2007, The Astrophysical Journal, 665, 1038, doi: 10.1086/519833/XML
  * Ciotti et al. (2022) Ciotti, L., Ostriker, J. P., Gan, Z., et al. 2022, The Astrophysical Journal, 933, 154, doi: 10.3847/1538-4357/ac70c7
  * Ciotti et al. (2010) Ciotti, L., Ostriker, J. P., & Proga, D. 2010, The Astrophysical Journal, 717, 708, doi: 10.1088/0004-637X/717/2/708
  * Cochrane et al. (2023) Cochrane, R. K., Anglés-Alcázar, D., Mercedes-Feliz, J., et al. 2023, eprint arXiv:2303.12858, 000, arXiv:2303.12858, doi: 10.48550/ARXIV.2303.12858
  * Crain & Voort (2023) Crain, R. A., & Voort, F. v. d. 2023, https://doi.org/10.1146/annurev-astro-041923-043618, 61, doi: 10.1146/ANNUREV-ASTRO-041923-043618
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, Monthly Notices of the Royal Astronomical Society, 365, 11, doi: 10.1111/j.1365-2966.2005.09675.x
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, Monthly Notices of the Royal Astronomical Society, 486, 2827, doi: 10.1093/MNRAS/STZ937
  * de Vaucouleurs et al. (1948) de Vaucouleurs, G., Vaucouleurs, d., & Gerard. 1948, AnAp, 11, 247.  https://ui.adsabs.harvard.edu/abs/1948AnAp...11..247D/abstract
  * Di Matteo et al. (2023) Di Matteo, T., Angles-Alcazar, D., Shankar, F., et al. 2023, arXiv, arXiv:2304.11541, doi: 10.48550/ARXIV.2304.11541
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604, doi: 10.1038/nature03335
  * Dubois et al. (2012) Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, Monthly Notices of the Royal Astronomical Society, 420, 2662, doi: 10.1111/J.1365-2966.2011.20236.X/2/MNRAS0420-2662-F17.JPEG
  * Dubois et al. (2016) Dubois, Y., Peirani, S., Pichon, C., et al. 2016, Monthly Notices of the Royal Astronomical Society, 463, 3948, doi: 10.1093/MNRAS/STW2265
  * Faber & Jackson (1976) Faber, S. M., & Jackson, R. E. 1976, Astrophysical Journal, Vol. 204, p. 668-683 (1976), 204, 668, doi: 10.1086/154215
  * Fabian (2012) Fabian, A. C. 2012, Annu. Rev. Astron. Astrophys, 50, 455, doi: 10.1146/annurev-astro-081811-125521
  * Faucher-Giguère & Quataert (2012) Faucher-Giguère, C.-A., & Quataert, E. 2012, MNRAS, 425, 605, doi: 10.1111/j.1365-2966.2012.21512.x
  * Faucher-Giguère et al. (2012) Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2012, MNRAS, 420, 1347, doi: 10.1111/j.1365-2966.2011.20120.x
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, The Astrophysical Journal, 539, L9, doi: 10.1086/312838
  * Feruglio et al. (2010) Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, Astronomy and Astrophysics, 518, L155, doi: 10.1051/0004-6361/201015164
  * Fischer et al. (2010) Fischer, J., Sturm, E., González-Alfonso, E., et al. 2010, Astronomy and Astrophysics, 518, L41, doi: 10.1051/0004-6361/201014676
  * Flores Velázquez et al. (2021) Flores Velázquez, J. A., Gurvich, A. B., Faucher-Giguère, C. A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 501, 4812, doi: 10.1093/MNRAS/STAA3893
  * Gallazzi et al. (2006) Gallazzi, A., Charlot, S., Brinchmann, J., & White, S. D. 2006, Monthly Notices of the Royal Astronomical Society, 370, 1106, doi: 10.1111/J.1365-2966.2006.10548.X/2/MNRAS0370-1106-F19.JPEG
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, The Astrophysical Journal, 539, L13, doi: 10.1086/312840/FULLTEXT/
  * Greene et al. (2020) Greene, J. E., Strader, J., & Ho, L. C. 2020, https://doi.org/10.1146/annurev-astro-032620-021835, 58, 257, doi: 10.1146/ANNUREV-ASTRO-032620-021835
  * Gurvich et al. (2022) Gurvich, A. B., Gurvich, & B., A. 2022, ascl, ascl:2202.006.  https://ui.adsabs.harvard.edu/abs/2022ascl.soft02006G/abstract
  * Gurvich et al. (2023) Gurvich, A. B., Stern, J., Faucher-Giguère, C. A., et al. 2023, Monthly Notices of the Royal Astronomical Society, 519, 2598, doi: 10.1093/mnras/stac3712
  * Habouzit et al. (2017) Habouzit, M., Volonteri, M., & Dubois, Y. 2017, Monthly Notices of the Royal Astronomical Society, 468, 3935, doi: 10.1093/mnras/stx666
  * Harrison (2017) Harrison, C. M. 2017, Nature Astronomy, 1, 0165, doi: 10.1038/s41550-017-0165
  * Hopkins (2015) Hopkins, P. F. 2015, Monthly Notices of the Royal Astronomical Society, 450, 53, doi: 10.1093/mnras/stv195
  * Hopkins (2016) —. 2016, Monthly Notices of the Royal Astronomical Society, 462, 576, doi: 10.1093/MNRAS/STW1578
  * Hopkins et al. (2020) Hopkins, P. F., Grudić, M. Y., Wetzel, A., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 3702, doi: 10.1093/MNRAS/STZ3129
  * Hopkins et al. (2007) Hopkins, P. F., Hernquist, L., Cox, T. J., Robertson, B., & Krause, E. 2007, The Astrophysical Journal, 669, 45, doi: 10.1086/521590
  * Hopkins et al. (2010) Hopkins, P. F., Murray, N., Quataert, E., & Thompson, T. A. 2010, Monthly Notices of the Royal Astronomical Society: Letters, 401, L19, doi: 10.1111/j.1745-3933.2009.00777.x
  * Hopkins & Quataert (2011) Hopkins, P. F., & Quataert, E. 2011, Monthly Notices of the Royal Astronomical Society, 415, 1027, doi: 10.1111/j.1365-2966.2011.18542.x
  * Hopkins & Raives (2016) Hopkins, P. F., & Raives, M. J. 2016, Monthly Notices of the Royal Astronomical Society, 455, 51, doi: 10.1093/MNRAS/STV2180
  * Hopkins et al. (2022a) Hopkins, P. F., Squire, J., Butsky, I. S., & Ji, S. 2022a, Monthly Notices of the Royal Astronomical Society, 517, 5413, doi: 10.1093/mnras/stac2909
  * Hopkins et al. (2016) Hopkins, P. F., Torrey, P., Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2016, MNRAS, 458, 816, doi: 10.1093/mnras/stw289
  * Hopkins et al. (2022b) Hopkins, P. F., Wellons, S., Angles-Alcázar, D., Faucher-Giguère, C. A., & Grudić, M. Y. 2022b, Monthly Notices of the Royal Astronomical Society, 510, 630, doi: 10.1093/mnras/stab3458
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/sty1690
  * Hopkins et al. (2023) Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, Monthly Notices of the Royal Astronomical Society, 519, 3154, doi: 10.1093/mnras/stac3489
  * Ji et al. (2021) Ji, S., Kereš, D., Chan, T. K., et al. 2021, MNRAS, 505, 259, doi: 10.1093/mnras/stab1264
  * Ji et al. (2020) Ji, S., Chan, T. K., Hummels, C. B., et al. 2020, MNRAS, 496, 4221, doi: 10.1093/mnras/staa1849
  * Lauer et al. (2007) Lauer, T. R., Gebhardt, K., Faber, S. M., et al. 2007, The Astrophysical Journal, 664, 226, doi: 10.1086/519229
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, The Astrophysical Journal Supplement Series, 123, 3, doi: 10.1086/313233/FULLTEXT/
  * Liu et al. (2013) Liu, G., Zakamska, N. L., Greene, J. E., Nesvadba, N. P. H., & Liu, X. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2576, doi: 10.1093/mnras/stt1755
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, The Astronomical Journal, 115, 2285, doi: 10.1086/300353
  * Mercedes-Feliz et al. (2023) Mercedes-Feliz, J., Angles-Alcázar, D., Hayward, C. C., et al. 2023, Monthly Notices of the Royal Astronomical Society, 524, 3446, doi: 10.1093/mnras/stad2079
  * Moe et al. (2009) Moe, M., Arav, N., Bautista, M. A., & Korista, K. T. 2009, The Astrophysical Journal, 706, 525, doi: 10.1088/0004-637X/706/1/525
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, The Astrophysical Journal, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Naab & Ostriker (2017) Naab, T., & Ostriker, J. P. 2017, https://doi.org/10.1146/annurev-astro-081913-040019, 55, 59, doi: 10.1146/ANNUREV-ASTRO-081913-040019
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, Monthly Notices of the Royal Astronomical Society, 475, 624, doi: 10.1093/mnras/stx3040
  * Parsotan et al. (2021) Parsotan, T., Cochrane, R. K., Hayward, C. C., et al. 2021, Monthly Notices of the Royal Astronomical Society, 501, 1591, doi: 10.1093/MNRAS/STAA3765
  * Peng et al. (2010) Peng, Y. J., Lilly, S. J., Kovač, K., et al. 2010, The Astrophysical Journal, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, Astronomy and Astrophysics, 641, A6, doi: 10.1051/0004-6361/201833910
  * Power et al. (2003) Power, C., Navarro, J. F., Jenkins, A., et al. 2003, Monthly Notices of the Royal Astronomical Society, 338, 14, doi: 10.1046/J.1365-8711.2003.05925.X/2/M_338-1-14-EQ037.JPEG
  * Richings & Faucher-Giguère (2018) Richings, A. J., & Faucher-Giguère, C.-A. 2018, MNRAS, 478, 3100, doi: 10.1093/mnras/sty1285
  * Richings & Faucher-Giguère (2018) Richings, A. J., & Faucher-Giguère, C.-A. 2018, Monthly Notices of the Royal Astronomical Society, 474, 3673, doi: 10.1093/mnras/stx3014
  * Richings et al. (2021) Richings, A. J., Faucher-Giguère, C.-A., & Stern, J. 2021, MNRAS, 503, 1568, doi: 10.1093/mnras/stab556
  * Rupke & Veilleux (2013) Rupke, D. S. N., & Veilleux, S. 2013, The Astrophysical Journal, 775, L15, doi: 10.1088/2041-8205/775/1/L15
  * Salim et al. (2007) Salim, S., Rich, R. M., Charlot, S., et al. 2007, The Astrophysical Journal Supplement Series, 173, 267, doi: 10.1086/519218/FULLTEXT/
  * Scannapieco et al. (2008) Scannapieco, C., Tissera, P. B., White, S. D. M., & Springel, V. 2008, Mon. Not. R. Astron. Soc, 389, 1137, doi: 10.1111/j.1365-2966.2008.13678.x
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, Monthly Notices of the Royal Astronomical Society, 446, 521, doi: 10.1093/mnras/stu2058
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, Mon. Not. R. Astron. Soc, 343, 978.  https://academic.oup.com/mnras/article/343/3/978/1125459
  * Sijacki et al. (2007) Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, Monthly Notices of the Royal Astronomical Society, 380, 877, doi: 10.1111/J.1365-2966.2007.12153.X/2/MNRAS0380-0877-F17.JPEG
  * Silk & Rees (1998) Silk, J., & Rees, M. J. 1998, Astronomy and Astrophysics, v.331, p.L1-L4 (1998), 331, L1, doi: 10.48550/ARXIV.ASTRO-PH/9801013
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, Annu. Rev. Astron. Astrophys, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Somerville et al. (2008) Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, Monthly Notices of the Royal Astronomical Society, 391, 481, doi: 10.1111/j.1365-2966.2008.13805.x
  * Springel et al. (2005) Springel, V., Di Matteo, T., & Hernquist, L. 2005, The Astrophysical Journal, 620, L79, doi: 10.1086/428772
  * Stern et al. (2021) Stern, J., Faucher-Giguère, C.-A., Fielding, D., et al. 2021, The Astrophysical Journal, 911, 88
  * Sturm et al. (2011) Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, The Astrophysical Journal, 733, L16, doi: 10.1088/2041-8205/733/1/L16
  * Su et al. (2019) Su, K. Y., Hopkins, P. F., Hayward, C. C., et al. 2019, Monthly Notices of the Royal Astronomical Society, 487, 4393, doi: 10.1093/MNRAS/STZ1494
  * Su et al. (2020) —. 2020, Monthly Notices of the Royal Astronomical Society, 491, 1190, doi: 10.1093/MNRAS/STZ3011
  * Su et al. (2021) Su, K. Y., Hopkins, P. F., Bryan, G. L., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 175, doi: 10.1093/MNRAS/STAB2021
  * Tomczak et al. (2014) Tomczak, A. R., Quadri, R. F., Tran, K. V. H., et al. 2014, The Astrophysical Journal, 783, 85, doi: 10.1088/0004-637X/783/2/85
  * Torrey et al. (2020) Torrey, P., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 497, 5292, doi: 10.1093/MNRAS/STAA2222
  * Tremmel et al. (2017) Tremmel, M., Karcher, M., Governato, F., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 1121, doi: 10.1093/mnras/stx1160
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, Monthly Notices of the Royal Astronomical Society, 465, 3291, doi: 10.1093/mnras/stw2944
  * Weinberger et al. (2018) Weinberger, R., Springel, V., Pakmor, R., et al. 2018, Monthly Notices of the Royal Astronomical Society, 479, 4056, doi: 10.1093/mnras/sty1733
  * Weinberger et al. (2023) Weinberger, R., Su, K. Y., Ehlert, K., et al. 2023, Monthly Notices of the Royal Astronomical Society, 523, 1104, doi: 10.1093/MNRAS/STAD1396
  * Wellons et al. (2020) Wellons, S., Faucher-Giguère, C. A., Anglés-Alcázar, D., et al. 2020, Monthly Notices of the Royal Astronomical Society, 497, 4051, doi: 10.1093/MNRAS/STAA2229
  * Wellons et al. (2023) Wellons, S., Faucher-Giguère, C. A., Hopkins, P. F., et al. 2023, Monthly Notices of the Royal Astronomical Society, 520, 5394, doi: 10.1093/mnras/stad511
  * Yoon et al. (2019) Yoon, D., Yuan, F., Ostriker, J. P., Ciotti, L., & Zhu, B. 2019, The Astrophysical Journal, 885, 16, doi: 10.3847/1538-4357/ab45e8