# Seeking Self-Regulating Simulations of Idealized Milky Way-Like Galaxies

Claire Kopenhafer Department of Physics & Astronomy, Michigan State University
Department of Computational Mathematics, Science, & Engineering, Michigan
State University  Brian W. O’Shea Department of Physics & Astronomy, Michigan
State University  Department of Computational Mathematics, Science, &
Engineering, Michigan State University  Facility for Rare Isotope Beams,
Michigan State University  G. Mark Voit Department of Physics & Astronomy,
Michigan State University

###### Abstract

Precipitation is potentially a mechanism through which the circumgalactic
medium (CGM) can regulate a galaxy’s star formation. Here we present idealized
simulations of isolated Milky Way-like galaxies intended to examine the
ability of galaxies to self-regulate their star formation, particularly via
precipitation. Our simulations are the first CGM-focused idealized models to
include stellar feedback due to the explicit formation of stars. We also
examine the impact of rotation in the CGM. Using six simulations, we explore
variations in the initial CGM $t_{\rm cool}/t_{\rm ff}$ ratio and rotation
profile. Those variations affect the amount of star formation and gas
accretion within the galactic disk. Our simulations are sensitive to their
initial conditions, requiring us to gradually increase the efficiency of
stellar feedback to avoid destroying the CGM before its gas can be accreted.
Despite this gradual increase, the resulting outflows still evacuate large,
hot cavities within the CGM and even beyond $r_{200}$. Some of the CGM gas
avoids interacting with the cavities and is able to feed the disk along its
midplane, but the cooling of feedback-heated gas far from the midplane is too
slow to supply the disk with additional gas. Our simulations illustrate the
importance of physical mechanisms in the outer CGM and IGM for star formation
regulation in Milky Way-scale halos.

††software: NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020),
Matplotlib (Hunter, 2007), Pandas (McKinney, 2010), yt (Turk et al., 2011),
and GNU parallel (Tange, 2020).

##  1 Introduction

Galaxies obey a number of scaling relations that suggest they have a mechanism
for regulating their star formation rates (SFRs). These scaling relations
include the star-forming main sequence (e.g., Popesso et al., 2022; Sherman et
al., 2021; Renzini & Peng, 2015) and the stellar mass–halo mass (SMHM)
relation (Somerville & Davé, 2015; Behroozi et al., 2019). Gas must be
continually supplied to a star-forming galaxy as its own interstellar medium
(ISM) only contains enough gas to fuel star formation for a few gigayears. The
circumgalactic medium (CGM) has gained much attention as a possible locale for
self-regulation as it mediates both inflowing and outflowing gas in addition
to functioning as a substantial baryonic reservoir in its own right (Tumlinson
et al., 2017).

One potential self-regulation mechanism is known as “precipitation” (Voit et
al., 2015). Under this framework, feedback (either stellar or from an AGN)
works to maintain the median $t_{\rm cool}/t_{\rm ff}$ ratio around 10.
Throughout this paper we’ll refer to the median value $\tau$ of $t_{\rm
cool}/t_{\rm ff}$ in the CGM as the “precipitation limit parameter.”
Importantly, $\tau$ is the median of a $t_{\rm cool}/t_{\rm ff}$ distribution
for the CGM gas that may have a broad dispersion (Voit, 2021). When $\tau\sim
10$, the $t_{\rm cool}/t_{\rm ff}$ distribution allows for the formation of a
multiphase medium. If $\tau$ declines, additional cold dense gas is able to
accrete onto the galaxy, driving a burst of star formation and feedback. This
feedback in turn heats and inflates the CGM, lowering its density and reducing
its ability to cool. Conversely, a CGM driven above $\tau\sim 10$ should feed
less cold, dense gas into the galaxy. This reduction of the cold gas supply
then stalls feedback until the CGM is again able cool and contract, thereby
lowering $\tau$.

Precipitation as a regulation mechanism was originally proposed for galaxy
cluster cores heated by active galactic nuclei (AGN). It was motivated by
observations suggesting that that the black hole fueling rate depends on the
development of a multiphase medium (Pizzolato & Soker, 2005; Cavagnolo et al.,
2008; Voit et al., 2008), which happens in simulations when a CGM in
approximate thermal balance becomes thermally unstable (Pizzolato & Soker,
2010; McCourt et al., 2012; Sharma et al., 2012). In simulated multiphase
media, cold clouds are able to precipitate out of the hot, ambient medium when
$t_{\rm cool}\,\lesssim\,10\,t_{\rm ff}$ (Gaspari et al., 2012, 2013, 2015; Li
& Bryan, 2014a, b; Prasad et al., 2015). If the medium is sufficiently
turbulent, these cold clouds are then able to accrete chaotically toward the
center, and in the systems originally considered, they trigger a strong AGN
feedback response. This feedback then heats and expands the CGM, raising
$t_{\rm cool}$ and diminishing accretion. See Donahue & Voit (2022) for an in-
depth review of this process.

Regulation via precipitation is observationally supported in both clusters
(Voit & Donahue, 2015) and elliptical galaxies (Frisbie et al., 2020). The
physical principles that underpin precipitation regulation are agnostic to the
source of feedback in a galaxy, and so could extend to smaller galaxies not
dominated by AGN feedback. This possibility is supported by recent theoretical
(Voit et al., 2019) and observational (Babyk et al., 2018) developments.

In this work, we explore the viability of precipitation as a regulatory
mechanism for star formation in Milky Way-like galaxies lacking AGN activity.
Self-regulation is broadly defined as a balance between gas accretion on the
one hand and stellar feedback on the other (Schaye et al., 2010; Bouché et
al., 2010; Davé et al., 2012; Zaragoza-Cardiel et al., 2019). This balance
might ultimately be what places galaxies on the stellar mass-halo mass
relation (Voit et al., 2015; Mitchell & Schaye, 2021).

Spiral galaxies like the Milky Way tend to live in complex environments where
they are continuously interacting with companion galaxies and cosmological
filaments. Such environments make it difficult to isolate the interactions of
the CGM and its host galaxy and to understand the long-term behavior of an
undisturbed Milky Way-like halo. We therefore use idealized simulations of an
isolated galaxy for our investigation, seeking to understand what a Milky Way-
like galaxy would look like without external interference. Idealized
simulations also allow us to maintain a constant dark matter halo mass, as
changing the dark matter mass should eventually change the stellar mass of the
galaxy according to the SMHM relation. With idealized simulations, we are also
able to achieve higher spatial resolution in the CGM for reasonable
computational cost. High CGM resolution broadens the distributions of gas
densities and temperatures (Corlies et al., 2020) and increases the overall
amount of neutral hydrogen (van de Voort et al., 2019). Both of these effects
should have an impact on precipitation (Voit, 2021).

Historically, idealized simulations of isolated Milky Way-like galaxies have
fallen into two types: simulations of a star-forming gaseous disk with no CGM,
and simulations of a CGM with no explicit star formation. The former include
the AGORA simulations of Kim et al. (2016), works using the AGORA initial
conditions such as Butsky & Quinn (2018) and Shin et al. (2021), and others
(Benincasa et al., 2016). Instead of a CGM, the gas disk of these simulations
is surrounded by a very low density medium with very long cooling times and
low total mass. In the latter category are the simulations from Fielding et
al. (2017) and Li & Tonnesen (2020). These CGM-focused works include stellar
feedback, but this feedback is tied to the gas flow rate through an inner
boundary rather than the explicit formation of stars in a gas disk. This work
is the first to combine both approaches, modeling star formation and the
resulting feedback in the context of both a gas disk and its surrounding non-
uniform CGM to create a fully self-consistent picture of how the galaxy and
its CGM interact. This is particularly important for modeling the inner CGM
and the disk-halo interface.

We detail our simulation setup in Section 2. Then, in Section 3, we examine
how the galactic disk is affected by our simulation setup and its variations.
In Section 4, we follow the movement of gas between the disk and the CGM by
tracing inflows, outflows, and the CGM’s gas supply. Section 5 looks at the
structure of the CGM in our simulations. We discuss various aspects of our
simulations in Section 6. Finally, we summarize and conclude our work in
Section 7.

##  2 Simulation Setup

We perform idealized simulations of isolated, Milky Way-like galaxies and
their circumgalactic media. Our simulations are performed with the Eulerian
astrophysical hydrodynamics code Enzo (Bryan et al., 2014; Brummel-Smith et
al., 2019). Enzo models stellar populations with particles and we include
feedback from Type II supernova, as will be discussed more in Section 2.2.

We employ a 9-species, non-equilibrium primordial chemistry network through
the Grackle cooling and chemistry library (Smith et al., 2017). This network
includes all ionization states of H, He, and $\rm H_{2}$, including $\rm
H^{-}$, as well as the density of free electrons. Grackle provides cooling
rates for H, He, and metals through precomputed Cloudy111A slightly modified
version of Cloudy 10 is used which saves outputs with more precision. tables
(Ferland et al., 2013). These tables include the extragalactic UV background
of Haardt & Madau (2012), which is fixed at $z=0$ for the purposes of our
simulations.

We adopt a static NFW dark matter potential (Navarro et al., 1996) with
$M_{\mathrm{200}}=10^{12}\ \mathrm{M_{\odot}}$ and a concentration of 10. We
define the edge of the dark matter halo using $r_{200}\approx 206$ kpc, which
is the radius enclosing a density that is 200 times the critical density.

On top of this we also impose a static stellar disk potential from Miyamoto &
Nagai (1975). For this potential we adopt $M_{\ast}=5.8\times 10^{10}\
\mathrm{M_{\odot}}$, which is the uncertainty-weighted average of observations
compiled in Table 6 of Côté et al. (2016). We also set the radial and vertical
scale heights of the stellar potential to match those of the gas for
simplicity (see Table 1). This potential is included because, though the
simulation is not initialized with any star particles, we want to model the
gravitational effects that the pre-existing stellar population of a $z\sim 0$
Milky Way-like galaxy would provide. In particular, this potential influences
the circular velocity of the disk gas. Star particles are allowed to form as
the simulation evolves, effectively allowing the stellar component of the disk
to grow from $M_{\ast}=5.8\times 10^{10}\ \mathrm{M_{\odot}}$.

Our galaxy and its background potentials are fixed to the center of a (1,638.4
kpc)3 box with periodic boundary conditions. While periodic boundary
conditions are not physically motivated, they are an easy way to ensure the
total mass of the box is conserved. Additionally, the width of the box is
about 8 times $r_{200}\approx 206$ kpc, ensuring that the periodic conditions
will not lead to any unexpected consequences.

The initial construction of the disk and CGM are covered in Sections 2.1.1 and
2.1.2, respectively. Section 2.2 covers the details of our star formation and
feedback models. In Section 2.3 we discuss our simulations’ approach to AMR.
Finally, Section 2.4 covers the multiple simulation variants we run.

###  2.1 Initial Conditions

The bulk values used to initialize our simulation can be found in Table 1.
Initial mass-averaged profiles are shown in Figure 1. The profiles are split
into disk and CGM components at small radii.

Our stellar mass is taken as the uncertainty-weighted average of the
observations in Table 6 of Côté et al. (2016). We then use the findings of
Peeples et al. (2014) to derive the ISM mass and metallicity based on this
stellar mass. The disk metallicity is $Z_{\rm ISM}=2.3\mathrm{\ Z_{\odot}}$,
assuming $M_{\mathrm{dust}}=M_{Z,\mathrm{dust}}$. Peeples et al. (2014)
suggests $M_{\rm ISM}=9.8\times 10^{9}\ \mathrm{M_{\odot}}$, but we lower this
to $7\times 10^{9}\ \mathrm{M_{\odot}}$ to help minimize the initial burst of
star formation experienced by the simulation; see Section 2.2.1 for more
discussion.

The CGM is initialized to a constant uniform metallicity of $0.3\
\mathrm{Z_{\odot}}$. This is the median value from Prochaska et al. (2017).
The metallicity of the CGM changes, however, as stellar populations inject
supernova feedback and drive outflows (see Section 2.2).

Quantity | Name | Value
---|---|---
$M_{200}$ | virial mass | $1.0\times 10^{12}\ \mathrm{M_{\odot}}$
$C_{\rm NFW}$ | dark matter concetration | 10
$r_{200}$ | virial radius | 206 kpc
$M_{\ast}$ | background stellar mass | $5.8\times 10^{10}\ \mathrm{M_{\odot}}$
$M^{\prime}_{ISM}$ | effective222The actual initial disk mass ends up being slightly higher than the value of the parameter used in Equation 1, which is $M_{\rm ISM}=5.0\times 10^{9}$. The actual disk mass that results is determined from a temperature cut, as the initial disk is isothermal and distinctly cooler than the surrounding CGM.  initial disk gas mass | $7.0\times 10^{9}\ \mathrm{M_{\odot}}$
$Z_{ISM}$ | initial disk metallicity | 2.3$\ \mathrm{Z_{\odot}}$
$R_{\rm s}$ | disk radial scale height | 3.5 kpc
$z_{\rm s}$ | disk vertical scale height | 0.325 kpc
$M_{CGM}$ | initial CGM mass (fiducial) | $2.3\times 10^{10}\ \mathrm{M_{\odot}}$
$Z_{CGM}$ | initial CGM metallicity | 0.3$\ \mathrm{Z_{\odot}}$
Table 1:  Important simulation parameters. See text for references and
derivations. Figure 1: In rows from top left, as a function of radius: initial
mass-weighted average density, temperature, entropy, pressure, and cooling
time, and total enclosed mass. The independent contributions of the disk
(dotted) and inner CGM (dashed) are shown for small radii, except for the
enclosed mass. For the enclosed mass, the entire CGM is shown as a dashed
line. Though the smallest cell size in the simulation is 100 pc, the x-axis
does not extend this far because of the low number of cells. The vertical
dashed line marks $r_{200}\approx 206$ kpc. The gray dot-dashed line shows
where $t_{\rm cool}/t_{\rm ff}=10$.

####  2.1.1 Gaseous Disk

The gas of the disk follows the softened profile of Tonnesen & Bryan (2009):

$\rho(R,z)=\frac{M_{\rm ISM}}{8\pi R_{\rm s}^{2}z_{\rm
s}}\operatorname{sech}\left(\frac{R}{R_{\rm
s}}\right)\operatorname{sech}\left(\frac{z}{z_{\rm s}}\right),$ (1)

where $R$ is the cylindrical radius, and $R_{\rm s}=3.5$ kpc and $z_{\rm
s}=0.325$ kpc are the scale heights of the gas. For $R>24\ \mathrm{kpc}$,
Equation 1 is multiplied by a smoothing factor of

$0.5\left[1+\cos\left(\pi\frac{R-24\ \mathrm{kpc}}{7.2\
\mathrm{kpc}}\right)\right]$ (2)

that tapers the radial edges of the disk. The disk gas is given a circular
velocity prescribed by the combination of our NFW and stellar potentials, with
$v_{\phi}=\sqrt{R\cdot g(r,R,z)}$.

The disk is set to an initial uniform temperature of $10^{5}$ K. In reality,
galaxy disks are usually at $10^{4}$ K or below, but gas cools very
efficiently around $10^{5}$ K. Our disk is initialized at this higher
temperature so that gas would be well above our temperature threshold for star
formation ($3\times 10^{3}$ K; see Section 2.2) but also rapidly cool. The
variation in disk density means that the temperature does not uniformly, which
helps spread out the initial burst of star formation across time.

The initial conditions for the disk and CGM are blended together based on
density; if $\rho_{\mathrm{disk}}(r)>\rho_{\mathrm{CGM}}(r)$ in a cell, that
cell is considered part of the disk. Otherwise, it is initialized to be part
of the CGM.

####  2.1.2 Circumgalactic Medium

The density and temperature structure of the CGM is set following Voit (2019),
who describes an entropy profile that precipitation-regulated galaxies in NFW
halos would follow under the assumption of hydrostatic equilibrium. This
assumption is reasonable for a $M_{\mathrm{200}}=10^{12}\ \mathrm{M_{\odot}}$
halo (Oppenheimer, 2018; Lochhaas et al., 2020). The profile therein assumes
the median $t_{\mathrm{c}}/t_{\mathrm{ff}}=10$, but we make $\tau\equiv
t_{\mathrm{c}}/t_{\mathrm{ff}}$ a free parameter in our model of the entropy
profile, so that

$\displaystyle K(r)$ $\displaystyle=(39\ \mathrm{keV\
cm^{2}})v_{200}^{2}\left(\frac{r}{r_{200}}\right)^{1.1}$ $\displaystyle+(2\mu
m_{p})^{1/3}\left[\frac{2n_{i}\tau}{3n}r\Lambda\left(2T_{\phi}(r),Z_{\rm
CGM}\right)\right]^{2/3},$ (3)

for entropy quantified in terms of $K\equiv kTn_{e}^{-2/3}$. This way we can
test halos that deviate from the apparent precipitation limit of $\tau\approx
10$. We assume $n_{e}=n_{i}$ such that $2n_{i}/n=1$. The initial metallicity
is $Z_{\rm CGM}=0.3$. The term $\Lambda$ is the cooling rate as a function of
of temperature and metallicity, which we calculate with the Grackle chemistry
and cooling library (Smith et al., 2017). The term $T_{\phi}(r)$ is the
gravitational temperature $kT_{\phi}(r)=\mu m_{p}v_{\rm c}^{2}(r)/2$.

To calculate the initial density and temperature profiles of the CGM, we
combine Equation 2.1.2 with $dP/dr$ from hydrostatic equilibrium. We adopt the
boundary condition $kT(r_{200})=0.25\mu m_{p}v^{2}_{c,\mathrm{max}}$ suggested
by Voit (2019). Beyond $r_{200}$, the temperature begins to plummet
dramatically. To mitigate this, we abandon a functional form for the entropy
and instead fix $d\log P/d\log r$ to its value just inside $r_{200}$. We then
set the temperature with a sigmoid function that smoothly blends $T(r_{200})$
with a floor of $4\times 10^{4}$ K. The resulting mass-weighted profiles for
the initial density, temperature, entropy, pressure, cooling time, and
enclosed mass are shown in Figure 1. For $\tau=10$, the actual cooling time
deviates slightly from the target value, falling between 7–10$t_{\rm ff}$
within $\sim 0.5r_{200}$ and rises to $\sim 20t_{\rm ff}$ at the virial
radius.

For simplicity, we use a modified NFW profile to calculate the circular
velocity $v_{\rm c}$ used in Equation 2.1.2. This modification approximates
the presence of a galactic disk at the center of the potential well: $v_{\rm
c}^{2}(r)=v_{\rm c,max}^{2}$ for $r\leq 2.163r_{s}$, where $r_{s}$ is the
scale radius of the NFW profile, and

$v_{\rm c}^{2}(r)=v_{\rm c,max}^{2}\cdot
4.625\left[\frac{\ln(1+r/r_{s})}{r/r_{s}}-\frac{1}{1+r/r_{s}}\right]$ (4)

otherwise. Using this modified NFW profile to represent the inner contents of
the halo is slightly incongruous with the construction of our gas disk and
static stellar potential; however, the modified NFW profile is spherically
symmetric, making it easier to calculate the CGM’s density and temperature
profiles. This modified NFW profile is only used to set up the initial CGM
conditions, while the unmodified NFW potential, the Miyamoto & Nagai (1975)
stellar potential, and the gaseous disk’s own gravity are what are actually
applied to the simulation during execution. Any incongruities introduced into
the initial hydrostatic equilibrium through our use of two slightly differing
potentials will quickly fade as the simulation evolves.

We give the CGM an initial azimuthal velocity, the strength of which is
determined by

$v_{\phi}(r,\theta)=v_{0}\sin^{2}\theta\left(\frac{r}{r_{0}}\right)^{\beta}.$
(5)

Here, $r$ is the spherical radius and $\theta$ is the polar angle. This
function replaces the disk’s Keplerian velocity profile wherever the CGM
density dominates, as described above. We do not consider Equation 5 when
placing the CGM in hydrostatic equilibrium, meaning the CGM may be slightly
oversupported against gravity in the cylindrically radial direction. We use
Hodges-Kluck et al. (2016) to choose $v_{0}=180$ km/s and $r_{0}=10$ kpc,
which is roughly the radius where their halo model’s specific angular momentum
matches that of the disk (see their Figure 5). It should be noted that Hodges-
Kluck et al. (2016) adopt a constant $v_{\phi}$ with their model, which they
assume is reasonable within $\sim 50$ kpc. We leave $\beta$ as a free
parameter; for our fiducial simulation, $\beta=-1/2$ (see Table 2).

###  2.2 Star Formation & Feedback

Our star formation and feedback algorithms are modified from the
implementation of Cen & Ostriker (1992) as described in Section 2.1 of Oh et
al. (2020). We require that gas in a grid cell have

  1. 1.

$\nabla\cdot v<0$

  2. 2.

either $t_{\rm cool}\leq t_{\rm dyn}$ or $T<3\times 10^{3}$ K

  3. 3.

$n\geq 10$ particles/cm-3

  4. 4.

$m_{\ast}\equiv f_{\ast}m_{\rm cell}\frac{\Delta t}{t_{\rm dyn}}>10^{4}\
\mathrm{M_{\odot}}$

where the local gas dynamical time is $t_{\rm dyn}=\sqrt{3\pi/(32G\rho)}$. If
all these criteria are met, a “star particle” is created representing a
population of individual stars with total mass $m_{\ast}$. An equivalent
amount of gas is also removed from the host cell, and the particle is given a
velocity such that momentum is conserved. The minimum stellar mass of $10^{4}\
\mathrm{M_{\odot}}$ is chosen as a balance between resolving the star
formation with more particles and the computational expense of tracking and
managing these particles. As in Smith et al. (2011), we ignore the Jeans
instability criterion as it is always met by the star forming gas in our
simulations. We set $f_{\ast}=0.2$ and impose a minimum dynamical time of one
million years.

Stellar feedback proceeds as described in Section 2.2 of Oh et al. (2020) and
is the same algorithm as used in Peeples et al. (2019). Though a star particle
is formed immediately, for the sake of stellar feedback, the accumulation of
stellar mass is assumed to be a drawn-out process that peaks in efficiency
after one dynamical time (Oh et al., 2020, Equation 3). The amount of mass,
momentum, and energy that are returned in timestep $\Delta t$ is then tied to
this extended star formation model through the mass of stars $\delta{M}_{\rm
SF}$ that would form in that timestep; e.g.,

$\Delta E=\epsilon_{\rm FB}\cdot c^{2}\delta{M}_{\rm SF}$ (6)

We adopt the same efficiencies as Oh et al. (2020) for the returned fraction
of total and metal masses. Energy, momentum, mass, and metals are then
deposited into the cube of 27 cells centered on the star particle’s host cell.

It should be noted that there are no star particles present in the simulation
initial conditions. This means that there are no pre-existing stellar
populations affecting the gas in the earliest moments of the simulation. This
artificiality has interesting consequences that are discussed in the following
subsection.

####  2.2.1 Ramping Stellar Feedback Efficiency

One of the consequences of idealized, isolated galaxy simulations like ours is
the behavior of the cold gas disk at early times. The dense gas of the
idealized disk is rotationally supported in the radial direction, but does not
have sufficient support against gravity in the vertical direction. It
therefore immediately starts to compress as the simulation begins to evolve,
allowing stars to form in a large volume of the disk as the density rises.
This leads to a large burst of star formation and feedback at very early
times, which disrupts the CGM before it can cool significantly. This means the
CGM has no time in which to influence the disk through accretion. Moreoever,
feedback raises cooling times and evacuates the gas above and below the disk
by successive supernova-driven shocks, preventing the CGM from interacting
with the disk on a reasonable timescale. The previous isolated galaxy
simulations listed in Section 1 did not have this problem because they either
lacked a significant CGM or did not explicitly model the star-forming disk.

Galaxies in the real universe, as well as in cosmological simulations, are
slowly built up over time. They do not “begin” as a $\sim 7\times 10^{9}\
\mathrm{M_{\odot}}$ disk of cold, dense gas. This assembly history means there
have been successive generations of both stars and stellar feedback. Real
galaxies will have effects such as stellar velocity dispersions and ISM
turbulence that support against vertical collapse. Additionally, the CGM’s
dynamical and thermodynamic structure will be impacted by the accretion of
matter over time and by historical outflows.

In spite of its artificiality, the value of an idealized galaxy system is that
it provides a controlled environment in which to explore the interactions of
the CGM and star formation. Ideally the modelled system can quickly “relax”
into a state where the oversimplification of the initial conditions has been
erased. We could minimize the collapse of the disk by adding, for instance, an
initial stellar distribution (so that feedback is active at the very start of
the simulation to limit the burst of star formation) or a 3D velocity
dispersion; however, both of these require multiple further parameter choices
to specify their distributions.

Instead, we choose to modulate the efficiency of stellar feedback,
$\epsilon_{\rm FB}$ (Equation 6). This dimensionless parameter is linearly
increased between $t=1$ Gyr and $t=2$ Gyr from $\epsilon_{\rm FB}=5\times
10^{-8}$ to $5\times 10^{-6}$. Our final $\epsilon_{\rm FB}$ is 0.5 dex lower
than both the cosmological simulations of (Oh et al., 2020) and FOGGIE
(Peeples et al., 2019). The efficiency ramp begins at 1 Gyr of simulation time
in order to give CGM gas with $t_{\rm cool}\lesssim 1$ Gyr time to interact
with the disk before being disrupted by feedback.

###  2.3 Resolution

Enzo uses block-structured Cartesian Adaptive Mesh Refinement (AMR; Berger &
Colella, 1989) to control the resolution of its grid cells. The base grid of
our simulation is $128^{3}$ cells (12.8 kpc per cell side) and the resolution
of a region is refined by a factor of two if the cell exceeds a mass threshold
of $2.67\times 10^{5}\ {\rm M_{\odot}}\times 2^{(-0.5l)}$, where $l$ is the
zero-based level index. The mass threshold for refinement therefore decreases
with level in a super-Lagrangian way, preventing excessive refinement on the
lowest levels which would increase computational cost. The highest levels are
concentrated in the galactic disk, which, while not the focus of this study,
is important to resolve for the purposes of star formation.

Given the low density of the CGM, a mass-based refinement criterion is not
enough for the CGM gas to become well-resolved. The importance of good CGM
resolution has been demonstrated by Hummels et al. (2019), Peeples et al.
(2019), Suresh et al. (2019), and van de Voort et al. (2019). Smaller clouds
are allowed to develop with higher CGM resolution, and turbulent structures
are resolved to smaller scales. We therefore also define six nested
rectangular regions of fixed minimum resolution (the first four of which are
cubic). Cells within these regions are allowed to refine further based on
mass. An example of this is shown in Figure 2. We allow up to seven levels of
refinement (a minimum spatial resolution of 100 pc).

Figure 2: Demonstration of our AMR scheme with edge-on (top) and zoomed face-
on (bottom) density projections at $t=3$ Gyr. Orange contours enclose regions
of uniform resolution at the center of the projection. Our refinement scheme
defines nested boxes that must be at a fixed minimum resolution, but cells are
allowed to refine beyond this minimum if they exceed a mass threshold. This
behavior can be seen with the 3,200, 1,600, and 200 pc contours. The 100 pc
contour is based exclusively on mass refinement and is the highest resolution
allowed.

###  2.4 Simulation Variants

Name | FB | $t_{\mathrm{c}}/t_{\mathrm{ff}}$ | $\beta$
---|---|---|---
Fiducial | Y | 10 | -1/2
CoolFlow | N | 10 | -1/2
LowRatio | Y | 5 | -1/2
HighRatio | Y | 20 | -1/2
LinRot | Y | 10 | -1
NoRot | Y | 10 | N/A
Table 2:  Our set of simulations and their varied parameters. Variations from
the fiducial model are highlighted in bold. The parameter $\beta$ refers to
the power law index of the initial azimuthal rotation profile.

The goal of our simulations is not only to test if we can create a
(precipitating) self-regulating system, but also to explore the robustness of
self-regulation. To that end, we explore five variations on our fiducial
simulation, as laid out in Table 2. The first variation, CoolFlow, is
identical to the fiducial run in its initial conditions but has stellar
feedback completely disabled (star formation is still allowed in order to
remove cold gas and preserve numerical stability). In this way, it functions
as a control to demonstrate the importance of feedback.

The next two variants, LowRatio and HighRatio, modify the precipitation limit
parameter $\tau$ in Equation 2.1.2. Their values of $\tau$ are chosen to be
the approximate lower and upper bounds experienced, on average, by a
precipitating system (Voit et al., 2017; Voit, 2018, 2021). We note that,
because the CGM’s radial density profile is defined by the initial $\tau$
through Equation 2.1.2, these variants start with CGM masses different by
about a factor of 2. Compared to the fiducial CGM mass given in Table 1, the
HighRatio variant starts with $1.3\times 10^{10}\ \mathrm{M_{\odot}}$ of CGM
gas and the LowRatio with $3.9\times 10^{10}\ \mathrm{M_{\odot}}$.

Though the CGM of galaxies likely rotates (Hodges-Kluck et al., 2016;
DeFelippis et al., 2020), the shape of its rotation profile is not well-
constrained. The LinRot variation modifies the initial rotation profile of the
CGM, changing the index $\beta$ in Equation 5, while NoRot is a control for
the impact of CGM rotation on the evolution of the system. This latter variant
has no initial rotation in the CGM, but still maintains Keplerian rotation in
the disk.

##  3 The Galactic Disk

We start by looking at the gas content of the galactic disk and its star
formation. In later sections, we will look in detail at the exchange of gas
between the disk and CGM, but the effects of CGM accretion are readily
apparent within the disk.

###  3.1 Appearance

Figure 3: Slices of density, temperature, and radial velocity through the disk midplane in each simulation variant at 3 Gyr. Cells are 100 pc thick in the center of the galaxy, 200 pc thick out to $\pm 25$ kpc, and 400 pc thick out to the edge of the field of view at $\pm 50$ kpc from the center. Solid white circles show the average radius within which stars form over 2–4 Gyr. Dotted white circles mark the extent of the initial gas disk ($\sim 28$ kpc). Note that the colormap for radial velocity covers $\pm 3\times 10^{3}$ km/s.  Variant | Mean (kpc) | Std Dev (kpc)
---|---|---
Cflow | 15.4 | 2.1
LowRatio | 11.7 | 4.2
LinRot | 10.4 | 3.7
NoRot | 10.4 | 5.1
Fid | 9.0 | 3.9
HighRatio | 8.1 | 3.7
Table 3:  Mean and standard deviation of the maximum star formation radius
over 2–4 Gyr of simulation time, in kpc. The maximum SF radius is defined as
the largest cylindrical radius at which a star formed in the last 50 Myr
(which is the time between simulation outputs).

Figure 3 compares the disks of our simulation variants at a simulation time of
3 Gyr. We show the density, temperature, and radial velocity within a slice
through the disk midplane. The solid white circles represent the average
radius within which stars form, which is determined as follows: for each
simulation snapshot ($\Delta t=50$ Myr), we determine the maximum cylindrical
radius of star particles formed since the last snapshot. The average of these
maxima over $t=2$–4 Gyr give the white circles in Figure 3. The numerical
values of these radii, as well as their standard deviations over time, are
reported in Table 3. The dotted circles show the edge of each variant’s
initial gas disk. Due to differences in initial CGM structure and the way the
disk and CGM are blended together (see Section 2.1.2), the LowRatio and
HighRatio simulations have slightly different initial radii than the other
variants, but all are around 28 kpc.

From these slices we see that the gas in the disk is around $10^{4}$ K on
average, reaching $10^{3}$ K and lower in the densest spiral arms and near the
center. Star formation is most prevalent in the very center of the disk. There
are hotspots of $\sim 10^{5}$ K gas within the star-forming center of each
variant’s disk that are the result of feedback. The exception is the CoolFlow
simulation, which has more prominent hotspots but no feedback. For this
variant, the gas heating is due to a “pile-up” of accreting gas. This pile-up
will be discussed further in Section 4.1.

The radial velocity slices show us that the hottest gas at $T\sim 10^{6}$ K is
not strictly outflowing. Clear examples of this can be seen in the right side
of the fiducial slices, in the bottom corner of the LowRatio images, and
throughout the NoRot slices. The fiducial simulation, in particular, shows
evidence of a radial velocity gradiant across a hot, low density cloud. This,
along with the highly variable gas temperature outside the galactic disk,
suggests that hot gas is mixing with cooler gas in the disk midplane. The
prevalence of radial inflow in hot gas also suggests that gas is cooling as it
moves towards the disk. Though highlighted in Figure 3 by LowRatio and NoRot,
cooling inflow is present for $t>2$ Gyr in all variations.

The size of the disk’s star forming region appears impacted by the CGM’s
$t_{\rm cool}/t_{\rm ff}$ properties, as seen from the solid white circles in
Figure 3 and the values in Table 3. For both of these, we use the average
maximum radius within which stars form as a measure of the star forming
region. The LowRatio variant has an SF region that is statistically larger
than the fiducial and the HighRatio variants (adopting $\alpha=0.05$ or a 95%
confidence interval), though we note that the fiducial’s region is not
signficantly larger than HighRatio’s. This suggests a small dependence on the
simulation’s initial $t_{\rm cool}/t_{\rm ff}$ ratio. Moreover, the fiducial,
LinRot, and NoRot simulations are statistically indistinguishable when tested
against each other333 A slight difference between the fiducial, LinRot, and
NoRot simulations appears when testing each of these simulations against other
variants. While the LowRatio simulation has a statistically larger star
forming region than the fiducial, it fails to test as significantly larger
than the LinRot and NoRot simulations. Additionally, while the fiducial
simulation does not test as significantly larger than the HighRatio
simulation, the LinRot and NoRot variants do. . These simulations all have the
same initial $t_{\rm cool}/t_{\rm ff}$ ratio. The CoolFlow simulation has a
star forming region that is statistically significantly larger than all the
other simulation variants, but this follows from its lack of feedback and
aggressive star formation rather than a difference in its CGM.

###  3.2 Mass Growth

Figure 4: Disk mass (top) and SFR (bottom) as a function of simulation time.
Disk mass (solid) is broken into gas (dotted) and stellar (dashed) components.
The stellar mass reflects the cumulative mass of stars formed and neglects
stellar mass loss. To improve clarity, the SFR is shown as the rolling mean
with a 60 Myr rectangular window. The region between the two dotted grey lines
is when the stellar feedback efficiency is ramped (Section 2.2.1). From left
to right, the Fiducial model (blue) is compared with CoolFlow, LowRatio and
HighRatio, and LinRot and NoRot. Compared to the CoolFlow variant, the mass
growth in each component remains fairly stable after the onset of feedback at
1 Gyr. The final stellar mass varies for each variant, along with the early-
time ($t<1$ Gyr) SFR and gas consumption. Figure 5: Mass of the stellar (top)
and gaseous (bottom) disk components relative to $t=2$ Gyr. The legend is
split between both panels. The stellar mass grows similarly among the non-
CoolFlow runs, but the gas mass varies widely.

Figure 4 shows the disk mass and star formation rate over the full simulation
time (4 Gyr) for each variant. The disk mass is split into gaseous and stellar
components. The disk gas is defined by a cylindrical region with $R\approx
5.7R_{\rm s}=20$ kpc and thickness $8z_{\rm s}=2.6$ kpc. This cutoff in
cylindrical radius is somewhat arbitrary, but corresponds roughly to a
transition point in the profiles of density and temperature that is visible
from 2–4 Gyr. The disk’s stellar mass is defined as the total mass of all
stars formed, and neglects stellar mass loss. Including stellar mass loss
lowers the stellar masses at $t>1$ Gyr by $\sim 1.5\times 10^{9}\
\mathrm{M_{\odot}}$. The accompanying SFR has been smoothed using a trailing
moving average over a $\sigma=60$ Myr rectangular window. Vertical dotted
lines denote the beginning and end of the feedback efficiency ramp laid out in
Section 2.2.1. The fiducial simulation is recreated in all three panels.

In the leftmost panel, the fiducial simulation is compared to the CoolFlow
run. For the fiducial model, the growth in stellar mass slows past 1 Gyr, when
the stellar feedback efficiency begins to ramp. In the CoolFlow model,
however, the stellar mass continues to grow significantly as the simulation
progresses, depleting the disk’s gas mass more appreciably than in the
fiducial variant. The SFR is also generally higher in the CoolFlow simulation.
This behavior is expected as the CoolFlow model has stellar feedback disabled.
This also indicates the flattening of the fiducial model’s stellar growth is
due to feedback, which is confirmation of an expected result.

In the middle panel, variations in the initial $t_{\rm cool}/t_{\rm ff}$ ratio
are explored by comparing the fiducial simulation to the LowRatio and
HighRatio runs. The rightmost panel explores the rotation variants LinRot and
NoRot. All of the variants with stellar feedback demonstrate remarkably
similar gas mass curves. Stellar growth and gas mass loss flatten out around 1
Gyr, when the feedback efficiency begins to increase. The variants are
primarily distinguished by the total amount of stellar mass they form, with
the largest difference observed between the $t_{\rm cool}/t_{\rm ff}$
variants.

For the first $\sim 0.7$ Gyr, the fiducial and LinRot simulations have the
same stellar mass curves. This is also true for the fiducial and CoolFlow
simulations before $\sim 0.5$ Gyr. Since these simulations all start with an
identical initial disk, deviations in stellar mass can be attributed to gas
accretion from the CGM. The NoRot variant has little apparent overlap in
stellar mass growth with the fiducial simulation, indicating an earlier
deviation in the amount of CGM accretion. We attribute NoRot’s early
differentiation to the lack of CGM angular momentum, which should make it
easier for gas to settle onto the disk and form stars. The LowRatio and
HighRatio simulations start with slightly smaller/larger initial disk radii,
respectively, due to how the disk and CGM are blended (Section 2.1.1). This
has a negligible effect on their disk mass, while their CGM masses deviate
from the fiducial by about a factor of two with LowRatio being the most
massive. Given that the LowRatio simulation goes on to form the most stars (of
the variants with feedback), this is further evidence that CGM accretion is a
major contributor to stellar growth.

The fiducial model and all other variants with $\tau=10$ start with about
$8\times 10^{8}\ \mathrm{M_{\odot}}$ of CGM gas with $t_{\rm cool}\leq 1$ Gyr,
which is the timescale on which stellar feedback is kept inefficient. The
LowRatio variant has $\sim 4\times 10^{9}\ \mathrm{M_{\odot}}$ of this short-
cooling time gas ($\sim 5$ times greater than the fiducial), and forms the
most stars before 1 Gyr. On the other hand, the HighRatio variant forms the
fewest stars by 1 Gyr, and its CGM only has $\sim 1\times 10^{8}\
\mathrm{M_{\odot}}$ of gas with $t_{\rm cool}\leq 1$ Gyr (a factor of 8 lower
than the fiducial). These differences support the interpretation that star
formation within the first gigayear is fueled not only by disk mass, but also
by the accretion of gas from the CGM. Additionally, there is evidence for
continued CGM inflow as the gas mass continues to increase past 1 Gyr. The
inflow of gas will be discussed more in Section 4.2.

Though each variant experiences a different amount of star formation,
particularly at early times, their gas masses remain both relatively constant
and similar across variants. Additionally, the smoothed star formation rates
are also quite consistent once stellar feedback reaches its full strength at
$t=2$ Gyr. We attribute the constancy of the disk mass among simulation
variants to the Toomre criterion (Toomre, 1964). A mass of $\sim 4\times
10^{9}\ \mathrm{M_{\odot}}$ appears to be the value at which the disk is
marginally stable against gravitational instability. As gas is added to an
isothermal disk, the disk gas is able to fragment and form stars until the
surface density falls back below the threshold for instability. This effect
can also explain why the LowRatio variant experiences more star formation from
1–2 Gyr while maintaining a fairy stable gas mass in the disk. The role of the
Toomre instability will be discussed more in Section 6.3.

In Figure 5 we focus on the late time ($t\geq 2$ Gyr) stellar and gas masses.
These quantities are shown relative to each simulation’s gas and stellar
masses at $t=2$ Gyr. Once again, stellar mass loss is ignored in order to show
the integrated star formation. Including stellar mass loss drops the final
stellar mass increase by $\sim 1\times 10^{8}\ \mathrm{M_{\odot}}$ for the
simulations with feedback and $\sim 4\times 10^{8}\ \mathrm{M_{\odot}}$ for
the CoolFlow variant.

The most dramatic stellar mass increase can again be seen in the CoolFlow
simulation, while the variants with feedback experience more gradual increases
in their stellar masses. The fiducial and LinRot simulations grow at a similar
rate, while the LowRatio and NoRot variants “taper off” in their stellar mass
increase after 2.5 Gyr. The HighRatio variant appears to be a “late bloomer,”
with stellar mass increases smaller than the other variants with feedback
until around 3.6 Gyr. All of the variants with feedback form between
3–$7\times 10^{8}\ \mathrm{M_{\odot}}$ more stars after $t=2$ Gyr.

These final increases in stellar mass are partially compensated by overall
decreases in gas mass, which, ignoring the NoRot and CoolFlow variants, are
between 1.5 and $5.5\times 10^{8}\ \mathrm{M_{\odot}}$ at $t=4$ Gyr. Each
simulation variant, however, took a very different path to reach this point.
Except for the LinRot and CoolFlow simulations, most variants see an increase
in their disk’s gas mass at some point before $\sim 2.3$ Gyr. The NoRot
simulation sees the largest increase in disk gas mass, and is the only
simulation to end with more gas at $t=4$ Gyr than it had at $t=2$ Gyr. Despite
this, it is not the variant with the largest increase in stellar mass. The
differences in disk gas mass between the variants highlight that the story of
accretion is complex. We will explore it more in Section 4.2.

##  4 Outflows and Inflows

Feedback processes in the disk are able to drive large-scale galactic
outflows. In our simulations, these “feedback processes” are limited to Type
II supernovae. Indeed, self-regulation requires that inflowing gas must
balance gas loss through outflows and star formation. Given that our idealized
models ignore galaxy mergers, gas filaments, and other cosmological processes,
their evolution is a story of inflowing and outflowing CGM gas. Therefore, in
this section we examine the exchange of gas between the CGM and the disk as
our simulations evolve.

###  4.1 System Dynamics

Figure 6: Average line-of-sight density, temperature, entropy, pressure, and
radial velocity of the fiducial run at four different times. Averages are
taken from a thin slab two radial scale heights $2R_{\rm s}=7$ kpc thick and
600 kpc square centered on the disk, which is shown edge-on. The white circle
shows $r_{200}\approx 206$ kpc. The entropy color bar separates gas with
$K>10^{2}$ keV cm2 in pink, which should not cool within a few Gyr. The
efficiency of stellar feedback is ramped between 1 and 2 Gyr (Section 2.2.1).

=1.75in

=1in

Figure 7: Same as Figure 6 but for each simulation variant at 3 Gyr. The
leftmost column is reproduced from Figure 6.

Figure 6 shows the evolution of the fiducial simulation in four thermodynamic
quantities—density, temperature, entropy, and pressure—and radial velocity.
Average values along the line of sight are shown at four times (in gigayear
intervals) using an edge-on thin slab. These slabs are two radial scale
heights ($2R_{\rm s}=7$ kpc) thick and 600 kpc square on their sides, centered
on the disk. By averaging through a thin subdomain, we can show a more
representative visualization of the domain (as compared to a slice) while not
oversmoothing turbulent gas structures. The white circle denotes $r_{200}$.
The ramp in stellar feedback efficiency, as described in Section 2.2.1, occurs
from 1–2 Gyr.

At 1 Gyr, the initial burst of star formation begins to subside as the
feedback efficiency starts to increase. The initial disk has collapsed as
described in Section 2.2.1, pushing out gas in the radial direction. Sound
waves were driven by the disk’s initial collapse, including a spherical wave
that has propagated beyond the virial radius at 1 Gyr. Its propagation has
resulted in a positive radial velocity throughout the domain. Gas behind this
wave has cooled from its initial temperature profile and started falling
inward. Net inward flow can be seen in blue in the bottom row of Figure 6 and
is approximately 50 kpc in radius at $t=1$ Gyr. Further evidence of sound
waves can be seen faintly in the pressure above and below the disk. These
intersecting wave patterns will eventually be wiped away by stellar feedback,
which we can see starting to happen at 2 Gyr.

At $t=2$ Gyr we see outflows have disrupted the inflowing gas above and below
the disk. Some gas inside these outflows has started falling back towards the
galaxy. By 3 Gyr, these outflows have expanded past the virial radius. Gas
near the edge of older outflows still exhibits some inward radial motion, but
is disrupted by younger outflows. The interior of these recent outflows has
very low densities and pressures; lower than at 2 Gyr. At both $t=2$ and 3
Gyr, cold dense gas is feeding the disk along its midplane.

By 4 Gyr, we can see that the history of successive feedback-driven outflows
has created a complex shock structure within 100–200 kpc directly above and
below the disk. Visible at the very top of the projection is an inflowing
front with high temperature and entropy, but very low density ($<10^{-32}$
g/cm3). Indeed, there is very little gas above and below the disk. Meanwhile,
accreting material continues to flow inward along the midplane. There is gas
infall at large radii that is outside the disk midplane, but this is prevented
from reaching the disk by younger outflows.

In Figure 7, we compare the fiducial simulation at 3 Gyr to each of the
variants, still using the line-of-sight average through a thin slab. With no
feedback to disrupt its CGM, the CoolFlow variant has a rough sphere of
inflowing gas around its disk. The high pressure region at the disk’s center
steadily grows in radius from the continually infalling gas. This increase in
pressure also drives an increase in temperature.

All other variants have outflow “plumes” that vary in size and shape. The
NoRot variant has a taller, narrower outflow envelope than the fiducial and
LinRot models, indicating that rotation in the CGM may be important for
distributing outflowing material such as metals evenly through the CGM.
Visible in the LinRot model is another high-entropy, low-density front of gas
falling back towards the disk. These features develop in many of the
simulations, but are prevented from reaching the disk by futher feedback.

Clearly visible in the CoolFlow variant, but present in all simulations except
NoRot, is the remnant of the initial disk’s radial spread due to vertical
collapse (Section 2.2.1). This gas is very low entropy, and though not shown,
has metallicity associated with the initial disk. This initial disk spreads
out the least for the NoRot simulation (reaching about half as far as the
fiducial simulation at $t=2$ Gyr), indicating the initial expansion of the
disk is encouraged by the CGM’s own angular momentum. For NoRot, this smaller
spreading feature is easily disrupted by feedback, which is why this feature
is not as prominent in Figure 7 and why hot gas appears much closer to the
spiral disk in Figure 3.

###  4.2 Gas Availability & Accretion

Figure 8: Net rate of change in disk mass $\dot{M}_{\rm disk}$ (left axis;
thick lines) compared to the star formation rate (right axis; gray and pink).
The rate of disk mass change is calculated as $\dot{M}_{\rm
disk}=\dot{M}_{\ast}+\dot{M}_{\rm gas}$. The dashed horizontal line marks the
zero-point for this quantity. Note that all displayed quantities are sampled
at different cadences: $\Delta t=50$ Myr for the $\dot{M}_{\rm disk}$ and
$\Delta t=0.6$ Myr for the gray SFR. The 60 Myr moving-averaged SFR from
Figure 4 is shown in pink for comparison, and the horizontal lines indicate
average net $\dot{M}_{\rm disk}$ over 1 Gyr.  Figure 9: Cumulative mass
distribution of CGM gas as a function of cooling time at five simulation
times: 0, 1, 2, 3, and 4 Gyr. The CGM is defined by $r<r_{200}=206$ kpc, with
a cylinder 40 kpc in diameter and $8z_{\rm s}=2.6$ kpc thick excised to remove
the disk. The top row compares the Fiducial model (blue) with the CoolFlow
variant, the middle row with LowRatio (green) and HighRatio (orange), and the
bottom row with LinRot (purple) and NoRot (brown). The initial distribution
(leftmost column) is identical for all but the $t_{\rm cool}/t_{\rm ff}$
variants.

Figure 10: The same as the middle (left) and bottom (right) rows of Figure 9,
but the mass distributions have been normalized in order to highlight
differences in shape at each timestamp. The simulations with varying initial
$t_{\rm cool}/t_{\rm ff}$ maintain cooling time distributions with similar
shapes. On the other hand, the variants with different CGM rotation profiles
vary considerably from each other across time.

In our simulations, there are two major sources of gas for star formation:
cold gas that was present in the disk from the initial conditions, and CGM gas
that accretes onto the galaxy. Dying stars also return gas to the disk, but
their contribution over a $\lesssim 1$ Gyr timescale isnegligible in
comparison to these other two sources (unlike in ellipticals; see Voit &
Donahue, 2011). The story of gas accretion can also be thought of as the story
of gas availability: what gas can cool and fall inward toward the disk? How
much is gas in a rotating CGM able to shed angular momentum? In Figure 4, the
differences in stellar mass formed by 1 Gyr are a clear indicator that our
simulation variants experience different rates of gas accretion and star
formation. In Section 3.2 we noted that these differences depend primarily on
the mass of CGM gas with $t_{\rm cool}<1$ Gyr and on the presence of rotation.
Figure 5 suggests that the simulation variants have unique and variable
patterns of gas inflow over time. In this section, we delve more deeply into
the ability of the CGM to contribute gas to the disk.

Figure 8 shows the net rate of change in disk mass for $t>2$ Gyr, which is
when stellar feedback is at its full efficiency. The disk is restricted to a
cylinder with $R=20$ kpc and thickness $8z_{\rm s}=2.6$ kpc at the center of
the simulation domain. Net $\dot{M}_{\rm disk}$ is calculated as the sum of
the star formation rate $\dot{M}_{\ast}$ and the rate of gas change
$\dot{M}_{\rm gas}$. If $\dot{M}_{\rm disk}=0$, gas consumption by star
formation precisely accounts for all the gas lost from the disk. If
$\dot{M}_{\rm disk}<0$, star formation accounts for only part of the gas loss,
meaning that outflows have removed some of the disk gas. And if $\dot{M}_{\rm
disk}>0$, then consumption of gas by star formation is more than compensated
by CGM accretion and the gas shed by stars, although the latter process
contributes on the order of $0.1\ \mathrm{M_{\odot}}\ \mathrm{yr^{-1}}$ or
less over $t=2$–4 Gyr.

We compare $\dot{M}_{\rm disk}$ to the star formation rate on the right axis,
but the SFR has been binned at a higher rate than was used for the net
$\dot{M}_{\rm disk}$ calculation ($\Delta t=0.6$ Myr for the SFR instead of
$\Delta t=50$ Myr for $\dot{M}_{\rm disk}$). A pink line shows the smoothed
SFR that was used in Figure 4. Horizontal bars show the average $\dot{M}_{\rm
disk}$ over 1 Gyr.

This figure makes it clear that the CoolFlow variant’s persistent star
formation (as seen in Figure 4) is accompanied by a persistent inflow of gas.
This is the expected behavior for a simulation without feedback. Intriguingly,
the NoRot simulation also sees a net growth in its disk gas from 2–3 Gyr. This
corresponds to the large growth (and then steady decline) in gas mass seen in
Figure 5.

The HighRatio variant has the most stable disk mass of the variants, being the
simulation with $\dot{M}_{\rm disk}$ that is consistently closest to zero. The
other variants exhibit much larger fluctuations (even CoolFlow, though its net
$\dot{M}_{\rm disk}$ is almost always positive). Yet for the variants with
feedback, averaging the net $\dot{M}_{\rm disk}$ over Gyr timescales brings
the net rate of change in disk mass closer to zero.

Visually, there is no clear time correlation between net $\dot{M}_{\rm disk}$
and SFR. Sometimes $\dot{M}_{\rm disk}$ decreases following a burst of star
formation, such as the fiducial simulation at $t\approx 2.4$ and 3.1 Gyr,
LinRot at 2.6 Gyr, NoRot at 2.2 Gyr, and HighRatio at 3.4 Gyr. Other times, a
dip in $\dot{M}_{\rm disk}$ precedes a burst of star formation, as in NoRot at
$t\approx 3.1$ Gyr and LinRot at 3.8 Gyr. Bursts of high instantaneous SFR and
more prolonged periods of steady star formation both tend to be close in time
to large changes in $\dot{M}_{\rm disk}$, but the sign of this change is not
consistent, nor is the magnitude or whether the star formation precedes or
follows a $\dot{M}_{\rm disk}$ change. This lack of a clear correlation will
be discussed more in Section 6.4.

In Figure 9 we look at the cooling time distribution of CGM gas. We define the
CGM as a sphere with radius $r_{200}=206$ kpc, with a cylinder 40 kpc in
diameter and $8z_{\rm s}=2.6$ kpc thick excised to remove the disk. This is
consistent with the region used for the disk in Figures 4, 5, and 8. The
cumulative mass distributions are presented in intervals of 1 Gyr, starting
from the initial conditions. The stellar feedback efficiency is ramped from
1–2 Gyr as explained in Section 2.2.1. The second column therefore represents
the cooling time structure of the CGM after it has had a chance to evolve and
interact with the disk, but before the impact of feedback. Dashed lines show
the cooling time mass distribution when the disk is not removed (but we still
restrict to gas with $r<r_{200}$).

The LowRatio and HighRatio variants have higher and lower CGM masses than the
fiducial simulation, respectively, as noted in Section 2.4. Indeed, the
initial $t_{\rm cool}/t_{\rm ff}$ ratio is the only thing that affects the
initial cooling time distribution.

The overall mass of the CGM drops over time in all variants, irrespective of
feedback. The majority of the mass loss (all but $\sim 2\times 10^{5}\
\mathrm{M_{\odot}}$) is an artifact of the initial conditions: as explained in
Section 4.1, the collapse of the initial disk pushes back on the CGM,
inflating it and pushing gas beyond $r_{200}$. Indeed, by 4 Gyr, the CGM mass
within $r_{200}$ has dropped by over an order of magnitude, with the fiducial
CGM mass dropping to $\sim 1.8\times 10^{9}\ \mathrm{M_{\odot}}$.

At 1 Gyr, we can see that each variant has increased its amount of CGM gas
with $t_{\rm cool}\sim 0.4$ Gyr. Given the difference between the solid and
dashed lines at $t_{\rm cool}\lesssim 0.4$ Gyr and the flatness of the disk-
included profile in this region, most of the very low cooling time gas has
been accreted onto the disk. The NoRot simulation is the one variant with
notably more gas with $t_{\rm cool}\lesssim 0.4$ Gyr left in its CGM at 1 Gyr.
This is intriguing given its rapid growth in stellar mass seen in Figure 4.

From 2–3 Gyr, the CoolFlow simulation has less CGM gas at short cooling times
($t_{\rm cool}<1$ Gyr) than in the fiducial, but a roughly equivalent amount
of gas at longer cooling times. This is the imprint of the CoolFlow
simulation’s constant inflow of gas (Figure 8), which depletes gas with short
cooling times. The fiducial simulation experiences gas inflow as well of
course, but in a reduced capacity thanks to stellar feedback: at the 0 and 1
Gyr snapshots—before the feedback efficiency is increased—the fiducial and
CoolFlow variants have essentially identical cooling time-mass distributions.
The distributions (both with and without the disk) deviate at 2 Gyr after
feedback has become effective. This also means that any variations in cooling
time distribution before 2 Gyr are due to CGM differences. Interestingly, the
fiducial and CoolFlow CGM distributions come roughly back into agreement at 4
Gyr (though the fiducial simulation has a more massive disk).

Figure 10 highlights the differences in shape between the CGM cooling time-
mass distributions of our variants by normalizing them. On the left we show
the $t_{\rm cool}/t_{\rm ff}$ variants: the fiducial, LowRatio, and HighRatio
simulations. On the right are our rotation variants: the fiducial, LinRot, and
NoRot runs.

We first focus on the $t_{\rm cool}/t_{\rm ff}$ variants. At 1 Gyr, the
LowRatio variant has more of its mass in gas with $1<t_{\rm cool}<10$ Gyr than
the other two, and slightly less in gas with $t_{\rm cool}<1$ Gyr. The excess
of gas with $t_{\rm cool}\sim 1$ Gyr seems consistent with the large, positive
$\dot{M}_{\rm disk}$ exhibited by the LowRatio simulations from $t=2$–2.5 Gyr
in Figure 8. The fiducial simulation also has positive $\dot{M}_{\rm disk}$
during this time, though of less magnitude than LowRatio. Furthermore, the
fiducial simulation has an intermediate amount of CGM gas with $1<t_{\rm
cool}<10$ Gyr. These differences in shape are accompanied by differences in
absolute mass of gas with $1<t_{\rm cool}<10$ Gyr. It nevertheless seems an
intuitive result that the (relative and absolute) amount of gas with $t_{\rm
cool}$ of a few gigayears affects the amount of gas accreted by a disk a few
gigayears in the future.

By 3 Gyr, we see that the three simulations have evolved cooling time-mass
distributions with the same general shape. These similarities are for the most
part retained by 4 Gyr. This suggests that, despite initial differences in the
CGM’s structure, feedback has impacted their cooling-time mass distributions
in a similar manner. This is despite the morphological differences in their
feedback (Figure 7).

The similarity between the fiducial and CoolFlow simulations in Figure 9 tells
us that the LowRatio, HighRatio, and fiducial cooling time distributions have
the same shape because feedback does not affect the cooling time mass
distribution much at all. This seems strange given the dramatic outflows in
Figure 7, but most of that outflowing gas has $t_{\rm cool}$ longer than a
Hubble time and therefore should not be expected to return to the disk. The
gas that is visible in Figure 9, particularly that with $t_{\rm cool}\lesssim
10$ Gyr, is gas that has necessarily not been heated by outflows. This is
likely gas near the disk midplane which avoids being heated.

What does have an impact on the CGM’s cooling-time structure is rotation, as
seen in the right side of Figure 10. Not only is there a difference between
NoRot and the other variants with rotating CGMs, but there are also variations
due to the different rotation profiles. We adopt a very straightforward
prescription for the CGM’s rotation (Equation 5), but these differences
highlight the importance of understanding the angular momentum of the CGM.

##  5 The Circumgalactic Medium

The CGM has a complex, multiphase structure that can be difficult to
encapsulate. The digestability of information must be balanced against the
loss of detail. We can see in Figures 6 and 7 that our simulations do not have
a spherically symmetric CGM and are instead closer to having cylindrical
symmetry. Trying to encapsulate the structure of the CGM with, e.g., mass-
weighted spherical profiles therefore constitutes a loss of information.
Furthermore, averages are biased towards the highest values, even when mass-
weighted.

Figure 11: Entropy, pressure, $t_{\rm cool}/t_{\rm ff}$, and gas mass radial
profiles of the fiducial CGM as a function of polar angle $\theta_{\rm cen}$
at late times ($t=3$–4 Gyr). The disk midplane corresponds to
$\theta=90^{\circ}$. Solid lines are the time average of the 50th percentiles
(medians), and shaded regions are of the 16th and 80th percentiles; see text
for more details on our analysis. The vertical dashed line at 20 kpc
represents the approximate extent of the disk adopted throughout this text,
while the horizontal dashed lines show the $t_{\rm cool}/t_{\rm ff}$ range of
5–20 predicted by precipitation theory.

Therefore, in an attempt to encapsulate the radial structure of the CGM, we
subdivide the gas within the virial radius based on its polar angle $\theta$.
This creates a set of cones (at the poles) and circular wedges. We choose to
have seven regions subtending $\Delta\theta=30^{\circ}$; for example, the
region with $\theta\in(0^{\circ},15^{\circ})$ covers the top pole of the
simulation domain, and $\theta\in(75^{\circ},105^{\circ})$ encompasses the
disk midplane. To keep our notation concise, we will refer to these regions
based on their central polar angle, $\theta_{\rm cen}$. For the previous
examples, these would be $\theta_{\rm cen}=0^{\circ}$ and $\theta_{\rm
cen}=90^{\circ}$.

Within each region, we radially bin quantities of interest from $r=2$ to
$r=206\ \mathrm{kpc}\,\approx r_{200}$ using 51 bins. This gives a bin width
of $\Delta r\approx 4$ kpc. For each radial bin, we find the 16th, 50th, and
84th percentiles. This analysis is done for each simulation output (whose
cadence is $\Delta t=50$ Myr). Figures 11 and 12 show the time-average of
these percentiles from $t=3$ to $t=4$ Gyr. We restrict the average to the last
gigayear as the feedback ramp is at its full strength and the profiles are
visually the most stable over time.

Figure 11 shows the radial profiles of entropy, pressure, $t_{\rm cool}/t_{\rm
ff}$, and gas mass for the fiducial simulation. A vertical dashed line at
$r=20$ kpc marks the rough edge of the disk used throughout our analysis. The
horizontal dashed lines in the $t_{\rm cool}/t_{\rm ff}$ panel indicate the
5–20 range predicted by precipitation (Voit et al., 2017; Voit, 2018, 2021).
For $\theta_{\rm cen}=90^{\circ}$ (the disk midplane), the 16th percentile of
entropy is cut off at $r\sim 30$ kpc, reaching values of $\sim 10^{-4}$ keV
cm-2 at the smallest radii. This is the influence of the galactic disk, which,
unlike in much of our analysis, we do not excise. The disk is also evident in
the 84th percentile of the mass profile.

It is clearly visible that all four quantities vary with polar angle: entropy
and $t_{\rm cool}/t_{\rm ff}$ decrease towards the disk midplane at
$90^{\circ}$ while cell mass and pressure rise. The simulation is roughly
symmetric about the disk plane. The higher entropy, lower pressure, and lower
mass near the poles are consistent with the bipolar outflows seen in Figure 7,
while the opposite trends at the disk midplane are consistent with the inflow
seen along the midplane. All simulation variants demonstrate these trends.
Near the poles (far from the disk midplane), the time-averaged $t_{\rm
cool}/t_{\rm ff}$ ratio tends toward a smoothly declining profile. At these
high angles, the time-averaged cooling time is approximately constant at
$t_{\rm c}\sim 2\times 10^{3}$ Gyr. The shape of $t_{\rm cool}/t_{\rm ff}$ is
set primarily by the freefall time and, by extension, the NFW halo profile. On
the other hand, only the $t_{\rm cool}/t_{\rm ff}$ profile centered on
$90^{\circ}$ (the plane of the disk) is remotely consistent with the $t_{\rm
cool}/t_{\rm ff}$ range of 5–20 predicted by precipitation. This is also the
polar angle with the highest cell mass, and is the angle at which most of the
cold gas inflow is located.

Figure 12: Entropy, $t_{\rm cool}/t_{\rm ff}$, and radial velocity profiles
for the cones at $\theta_{\rm cen}=90^{\circ}$ (i.e., covering the disk
midplane), averaged over $t=3$–4 Gyr, for the simulation variants with stellar
feedback. Solid and shaded lines are the same as in Figure 11. The vertical
dashed line at 20 kpc marks the approximate edge of the disk, and the
horizontal dashed lines mark $t_{\rm cool}/t_{\rm ff}=5$–20 and $v_{r}=0$
km/s. The initial entropy profiles are shown as dotted gray lines, and dotted
colored lines highlight the time-averaged 16th percentiles in radial velocity.

The mass dominance of the $\theta_{\rm cen}=90^{\circ}$ profiles, as well as
their association with the inflowing gas along the disk midplane in Figure 7,
motivates isolating these profiles for comparison across the simulation
variants. This is done in Figure 12, where the time-averaged profiles of
entropy, $t_{\rm cool}/t_{\rm ff}$, and radial velocity for the
$\theta\in(75^{\circ},105^{\circ})$ wedge are shown for all of the simulation
variants with stellar feedback. The time averaging is again restricted to the
last 1 Gyr of simulation evolution. A vertical line at 20 kpc again marks the
nominal edge of the disk, and horizontal dashed lines mark $t_{\rm
cool}/t_{\rm ff}=5$–20 and $v_{r}=0$ km/s. The initial entropy profiles are
shown as dotted gray curves. Colored dotted lines are used to highlight the
16th percentile of the radial velocity, which is predominantly negative. The
entropy profile within the disk is once again cut off to better demonstrate
the CGM.

For all simulations, the spread in entropy and $t_{\rm cool}/t_{\rm ff}$
between the time-averaged 16th and 84th percentiles crosses several orders of
magnitude. Though this spread is still less than the overall dependence on
polar angle (Figure 11), we can see from the radial velocity that our wedges
encompass both inflowing and outflowing gas. Figures 6 and 7 show that the
inflow region is very thin. We chose a large $\Delta\theta$ in order to
capture potential warping of this region (most evident in the fiducial
simulation; see the 3rd and 4th columns of Figure 6) and so inevitably we
capture a mix of gas phases.

Focusing on the $t_{\rm cool}/t_{\rm ff}$ ratio, we see that the median of gas
near the midplane is predominantly higher than predicted from precipitation
theory, though the 16th percentile does extend into the 5–20 range for all but
the NoRot simulation. Precipitation requires a local $t_{\rm cool}/t_{\rm
ff}\sim 1$ for gas to actually condense, which may not be captured by the 16th
percentile. Yet this Figure is strong evidence that our simulated galaxies are
not being regulated by precipitation, even if we search for precipitation only
in the disk midplane.

##  6 Discussion

We simulated a suite of idealized galaxies that are similar to the Milky Way
in order to explore the conditions under which self-regulating feedback might
arise. Our hypothesis was that these galaxies would naturally regulate their
star formation rates according to the predictions of precipitation theory if
initialized with a CGM having $t_{\rm cool}/t_{\rm ff}\sim 10$ (Voit, 2019).
This work therefore complements that of Prasad et al. (2020), who studied the
precipitation-regulation of larger AGN-dominated systems with SNIa feedback
(i.e., massive central ellipitcal galaxies in groups and clusters), in that it
explores the ability of less massive galaxies to regulate themselves through
stellar feedback alone.

In light of the results from Prasad et al. (2020), the question motivating our
work is: can we also create a self-regulating galaxy-CGM system but with
feedback coupled to star formation instead of AGN and older stellar
populations? We broadly take self-regulation to mean that feedback tunes the
net inflow of CGM gas to match the disk’s time-averaged star formation rate.
Since our only feedback mechanism is Type II SNe, this also implies that the
star formation rate of self-regulated systems would be tuned to the inflow of
cold gas.

Ultimately, we do not consider the galaxies we have simulated to be self-
regulating. Instead, the CGM in our simulations experiences large-scale
disruption due to outflows. Only gas along the disk midplane is of low enough
entropy to be able to cool and accrete onto the disk. Even this limited
accretion mode does not appear to be precipitation-regulated based on the most
general measure, the median ratio of cooling and freefall times. Our star
formation rates drop to very low values of order $0.1\ \mathrm{M_{\odot}}$/yr
after the onset of our feedback efficiency ramp. This limited gas accretion
and star formation keep the average disk growth, $\dot{M}_{\rm disk}$, near
zero. Rather than creating a system that maintains a moderate SFR, we have
created a system in which feedback essentially shuts off star formation.

We start our discussion with Section 6.1, where we examine the impacts of
variations in the simulation parameters. Then, in Section 6.2, we compare the
structural features of our simulations to other works, both observational and
theoretical. Section 6.3 highlights the role the Toomre instability plays in
the evolution of our galaxies. Our simulations do not have self-regulated star
formation, but their failure to self-regulate is illuminating. In Section 6.4
we interrogate the definition of self-regulation and highlight the ambiguity
as to the expected timescales involved. Finally, in Section 6.5 we discuss the
physical effects missing in our simulations that could better couple gas
accretion and stellar feedback.

###  6.1 Impact of CGM Variations

The overall behavior of our simulations is unaffected by variations in the
CGM: once stellar feedback is at its full strength the CGM becomes disrupted
by wind-driven bubbles and the SFR drops by an order of magnitude. The
exception is, of course, the CoolFlow variant which completely lacks stellar
feedback. Yet despite observing no impact on the overall behavior, we do see
important differences manifest in our simulation variations. The variations in
initial conditions exert their biggest influence at $t<1$ Gyr, when feedback
is made artificially weak. In the following subsections we will discuss the
differences resulting from changes to the initial $t_{\rm cool}/t_{\rm ff}$
ratio (Section 6.1.1) and from alterations to the CGM’s initial rotation
profile (Section 6.1.2).

####  6.1.1 Variation in Initial $t_{\rm cool}/t_{\rm ff}$

The LowRatio and HighRatio variants are distinguished by modifying the
condensation criterion $\tau$ in Equation 1. This parameter is essentially the
initial $t_{\rm cool}/t_{\rm ff}$ ratio of the CGM, though as discussed in
Section 2.1.2 the actual $t_{\rm cool}/t_{\rm ff}$ ratio deviates slightly.
Including the fiducial simulation, we sample $\tau\in[5,10,20]$.

Because the CGM’s density structure is set by $\tau$, these variants have
different starting CGM masses. The LowRatio variant has the most mass at
$3.9\times 10^{10}\ \mathrm{M_{\odot}}$, followed by the fiducial with
$2.3\times 10^{10}\ \mathrm{M_{\odot}}$ and the HighRatio with $1.3\times
10^{10}\ \mathrm{M_{\odot}}$.444 The CGM is defined as a sphere with
$r=r_{200}\approx 206$ kpc with a cylinder excised for the disk. This disk has
height $z=4z_{\rm s}=1.3$ kpc from the midplane and radius $R=[27.5,28.5,29]$
kpc for the LowRatio, fiducial, and HighRatio variants respectively.
Therefore, a higher $\tau$ results in an overall less massive CGM.

Moreover, a higher $\tau$ lowers the relative amount of gas with initial
$t_{\rm cool}<1$ Gyr. This is seen in Figure 9. For a real multiphase CGM,
this can be understood via the framework of Voit (2021): the median $t_{\rm
cool}/t_{\rm ff}$ of the CGM defines the center of a distribution in $t_{\rm
cool}/t_{\rm ff}$. Moving the median to higher values means that less of the
distribution covers low cooling times, and therefore less gas is able to
efficiently cool, condense, and reach the galaxy. Our simulations, however, do
not start with a multiphase CGM; rather, the temperature and density are
smooth, spherically symmetric functions of radius. In fact, the initial
temperature profile is very similar between the fiducial, LowRatio, and
HighRatio simulations. The differing masses of gas with $t_{\rm cool}<1$ Gyr
are then a result of both variations in the initial density profiles and the
different initial $t_{\rm cool}/t_{\rm ff}$.

The initial stellar mass growth of the HighRatio and LowRatio simulations is
easily described by the different masses of gas with $t_{\rm cool}<1$ Gyr.
This is seen in Figure 4, where at $t<1$ Gyr the LowRatio simulation has the
highest disk gas mass, greatest growth in stellar mass, and highest SFR.
Conversely, the HighRatio variant has the least growth in stellar mass and
lowest SFR, although its gas mass is not as distinct from the fiducial
simulation as the LowRatio simulation is.

The difference in low cooling time gas also seems to have a slight impact on
the physical size of the disks, as seen in Figure 3 and Table 3. The LowRatio
simulation has the most gas with $t_{\rm cool}<1$ Gyr and correspondingly has
the largest radius within which star formation occurs. The converse is not
quite true for the HighRatio simulation, whose average radius of star
formation is not statistically significantly different from the fiducial
simulation’s. The average star formation radii in Table 3 are determined over
$t=2$–4 Gyr, after the feedback ramp has ended and the bulk of star formation
has occurred. This therefore suggests that the initial differences in CGM gas
accretion between the LowRatio, HighRatio, and fiducial simulations have
lingering effects on the structure of the galactic disk. This effect would
appear to be on the same scale as natural variation in the maximum radius of
star formation. These differences are likely not a result of the difference in
initial disk radius, as the initial LowRatio disk is about 1 kpc smaller than
the fiducial’s gas disk.

Feedback does not generally seem to have much impact on the $t_{\rm cool}$
distribution function for CGM gas. This is seen between the $t_{\rm
cool}/t_{\rm ff}$ variants in Figure 10. Though the $t_{\rm cool}/t_{\rm ff}$
variants have different star formation histories and stellar masses, feedback
neither amplifies nor diminishes any pre-existing differences in the $t_{\rm
cool}$ distribution function. The similarity between the fiducial and CoolFlow
distributions at $t=4$ Gyr further indicates that feedback has minimal impact
on the CGM’s overall availability of low cooling time gas.

Even though the cooling time distribution is largely unaffected by feedback,
the $t_{\rm cool}/t_{\rm ff}$ variants do not evolve identically. In Figure 7
we can see that outflows in the LowRatio simulation have traveled less far
than in either the fiducial or HighRatio simulations. For $t<1$ Gyr, the
LowRatio simulation has the highest CGM pressure of the $t_{\rm cool}/t_{\rm
ff}$ variants, and HighRatio the lowest. This generally tends to remain true
for the outer CGM ($r\approx 150$–200 kpc) even as outflows drop the overall
CGM pressure.

The $\dot{M}_{\rm disk}$ in Figure 8 also shows interesting differences. The
HighRatio variant experiences the smallest $\dot{M}_{\rm disk}$ fluctuations
of the simulations considered here, suggesting that relatively little mass is
involved in the cycle of accretion and star formation. The fluctuations of
$\dot{M}_{\rm disk}$ in the LowRatio simulation are on par with and
occasionally larger than those experienced by the fiducial simulation.
Generally, the fluctuations in $\dot{M}_{\rm disk}$ get smaller as the
simulations go on, as seen by the average from $t=3$–4 Gyr.

It may be that the LowRatio simulation, being able to accrete more cool gas
initially, was set onto a cycle of large amounts of gas accretion followed by
productive periods of star formation, rather like an oscillator with a large
initial perturbation. Also like most oscillators found in nature, the LowRatio
simulation experiences damping in its cycle of accretion and star formation.
Though it is the simulation with the highest stellar mass, Figure 5 shows that
it had only moderate gains in stellar mass after 2 Gyr. This figure also shows
that the LowRatio simulation continued to accrete gas from $t=2$–3 Gyr,
despite the stronger feedback. This may be the LowRatio simulation “refueling”
in order to maintain a more moderate SFR at the end of the simulation. On the
other hand, the fiducial simulation gains the most stellar mass over $t=2$–4
Gyr in Figure 5, while the HighRatio simulation is with LowRatio near the
bottom of the pack, although it is something of a “late bloomer.” It may be
that the fiducial simulation has the most ideal conditions of these variants
for sustained star formation growth. Whether or not we consider it “self
regulating” is a question we defer to Section 6.4.

To summarize, the LowRatio, HighRatio, and fiducial simulations all exhibit a
number of differences that are a consequence of how the condensation criterion
$\tau$ determines both the initial $t_{\rm cool}/t_{\rm ff}$ ratio and density
structure of the CGM. These density differences result in differences in total
mass. The LowRatio simulation not only starts with more gas, but a larger
fraction of it is able to cool efficiently. The converse is true for the
HighRatio simulation. These differences are persistent even in the face of
feedback. The structure of the disk and its star formation therefore retain
their early differences throughout the simulation runtime.

####  6.1.2 Variation in CGM Rotation

We now consider variations that arise from the CGM’s rotation. Apparent from
Figure 4 is that both the presence of rotation, as well as its variation with
radius, has an impact on the ability of the CGM to supply gas to the disk.

We have indicated throughout this work that the NoRot simulation should be
able to accrete gas more efficiently because its CGM gas does not have to shed
angular momentum to reach the disk. This is evident in the stellar growth of
Figure 4 at $t<1$ Gyr. It may also explain the large growth in gas mass seen
over $t\sim 2$–3 Gyr in Figures 5 and 8. In Section 4.1, we noted a remnant of
the initial disk conditions that is present in the CGM: due to insufficient
vertical support against gravity, the initial disk collapses and spreads
outwards into the CGM. The NoRot variant’s initial disk has the smallest
amount of radial spread due to the lack of rotation in its CGM.

The mere inclusion of rotation in the CGM has a large effect on gas accretion
and simulation evolution, but the radial profile of that rotation also has an
impact. The LinRot simulation forms more stars than the fiducial model within
the first 1 Gyr, before stellar feedback has much effect. The LinRot
simulation also has lower angular momentum at all CGM radii, making it easier
for gas to accrete than in the fiducial simulation. Rotation is also able to
change the shape of the cooling time-mass distribution in Figure 10, affecting
the availability of gas that can be accreted in addition to the ease of
accretion.

Su et al. (2020) also performed isolated galaxy simulations that included a
rotating CGM. These simulations targeted cool-core clusters with halo masses
between $10^{12}$ and $10^{14}\ \mathrm{M_{\odot}}$. Rotation followed a
$\beta$-profile with $\beta=1/2$, set to be twice the net dark matter spin.
Rotation therefore comprised 10-15% of the CGM’s support against gravity, with
the rest being supplied by thermal energy. This is in contrast to our
simulations, where rotation was not considered in the calculation of
hydrostatic balance. The simulations of Su et al. (2020) were included in the
Fielding et al. (2020) meta analysis, where it is noted that these simulations
had an enhanced cold phase at $r<0.2r_{200}$ compared to the other isolated
galaxy simulations in the analysis. This enhancement is attributed to CGM
rotation—which was not included in Fielding et al. (2017) and Li & Tonnesen
(2020)—in agreement with our results.

DeFelippis et al. (2020) studied the CGM angular momentum of Milky Way-mass
galaxies in IllustrisTNG. They split these galaxies into samples based on the
stellar specific angular momentum, resulting in high- and low-momentum
populations. Both populations have inflowing cold CGM gas near the disk plane.
In the sample with high specific angular momentum, this inflow is well aligned
with the disk; however, the radial inflow disappears for $r<0.2r_{200}$. Much
like with our simulations, it is suggested that this is due to the presence of
strong rotation. We note that our rotation speeds (derived from Hodges-Kluck
et al., 2016) are higher than the rotational speeds measured by DeFelippis et
al. (2020).

###  6.2 Comparison of Structural Features to Other Works

It is worth comparing some of the structures seen in our simulations with
those seen or predicted in other works. These works include observations,
analytic predictions, and cosmological simulations. We’ll begin with the gas
disk, expand outward to the accretion flows seen along the disk midplane, and
finally consider the structure of our feedback-driven outflows.

####  6.2.1 Extended Gas Disk

Starting with the disk, we can see in Figure 3 that there is cold, dense,
rotating gas that extends beyond the average radius of star formation in all
of our variants. This is consistent with observations of spiral galaxies in
the THINGS survey (Leroy et al., 2008). These observations show that spiral
galaxies have high star formation efficiency in their H2-dominated cores. This
efficiency declines with radius in the extended neutral hydrogen disk.

The structure of our galactic disk is also consistent with the findings of
Lopez et al. (2020) and Tejos et al. (2021). These papers present two
instances of extended Mg II disks ($r\sim 20$–40 kpc) that are co-rotating
with the interstellar medium. Our extended cool disks are of order $r\sim
20$-30 kpc. We note that, while an extended cool disk is an emergent feature
of our simulations, we purposely aligned the angular momentum vectors of the
CGM and the disk. It is unclear what impact misaligned angular momentum
vectors might have on our extended gas disk.

####  6.2.2 Accretion Along the Midplane

Figures 6 and 7 make it clear that, while inflowing gas can be found
throughout the CGM, only inflowing gas near the disk midplane is able to reach
the disk. Inflowing gas not along the midplane is disrupted by outflows. The
exception to this is of course the CoolFlow variant, which lacks stellar
feedback.

Inflow within the disk plane is not unprecedented in cosmological simulations.
Trapp et al. (2022) find this inflow mode in the FIRE-2 simulations and note
that it is the dominant source of accretion there just as it is in our
simulations. Gas accretion in the EAGLE simulations is also anisotropic,
favoring low heights relative to the disk and inflow speeds of 20–60 km/s (Ho
et al., 2019). This inflow is predominantly cold gas found near the disk or in
low-angular momentum streams, and estimates of the mass rate are enough to
meet or exceed the SFR.

Within IllustrisTNG, DeFelippis et al. (2020) also see this inflow structure
in $10^{11.75}$–$10^{12.25}\mathrm{M_{\odot}}$ halos. Truong et al. (2021) see
a global anisotropy in the CGM, with density enhanced parallel to the disk
plane, and temperature and metallicity enhanced along the orthogonal minor
axis. Their anisotropies are more subtle than what we see in our own
simulations, being of only 0.1–0.3 dex. These anisotropies peak in Milky Way-
like galaxies with $M_{\ast}\sim 6\times 10^{10}\ \mathrm{M_{\odot}}$, which
is the mass at which super massive black hole feedback turns on in TNG.
Generally, the anisotropies are larger in galaxies with SMBH feedback, except
for metallicity: the metallicity anisotropy is more pronounced in star-forming
and disky galaxies.

The inflowing gas spans a temperature range from $\sim 5\times
10^{3}$–$5\times 10^{5}$ K, with the coldest gas living in denser filaments.
Gas condenses and cools while being part of a rotating inflow. This is
reminiscent of the accretion mode observed in Hafen et al. (2022) and
described in detail in Stern et al. (2021). Gas inflow is hot until it reaches
$\sim 20$ kpc scales, at which point angular momentum slows its inward motion.
Radiative cooling then exceeds heating due to compression, and the gas
temperature drops to $\sim 10^{4}$ K or below. Importantly, Hafen et al.
(2022) note that the angular momentum of inflowing gas aligns itself with the
galaxy before cooling, but because the angular momentum of our CGM and disk
are constructed to be aligned, we cannot make any comparison on this point.

####  6.2.3 Outflow Structure

Our simulations are dominated by outflows to an unrealistic degree. Very low
density cavities extend out to the virial radius by $t=4$ Gyr. Outflows have a
very wide opening angle of essentially $180^{\circ}$, covering the face of the
disk. The CGM essentially becomes obliterated above and below the disk.

That said, the structure of our outflows matches well to analytic models from
Lochhaas et al. (2018). These models describe how relatively slow-moving gas
bubbles can be inflated by fast galactic winds. In their model, winds drive a
forward shock through the CGM. Their interaction leads to the development of a
reverse shock, and between these two fronts lies a contact discontinuity. This
discontinuity separates material driven out by the winds from the swept-up CGM
gas. Lochhaas et al. (2018) refer to the region between the reverse shock and
the contact discontinuity as the “shocked wind.” It is this shocked wind that
is responsible for the growth of the gas bubble. Immediately behind this
shocked wind (inside the reverse shock) lies cool, unshocked gas, followed by
hot unshocked gas immediately next to the galaxy.

We can clearly see the shocked wind in the density projections of Figure 6,
often with winds inside winds. The models of Lochhaas et al. (2018) use a
continuous wind, while our simulations have multiple discrete winds driven by
episodic Type II supernova feedback. At 2 Gyr, a band of shocked winds is
visible both above and below the disk. Subsequent winds are able to travel
faster due to the evacuated region that follows behind the first (Lochhaas et
al., 2018). The thickness of the outermost shocked wind grows as its reverse
shock propagates. It travels with speeds on the order of $\sim 100$ km/s,
which is an order of magnitude lower than the feedback-driven winds,
corroborating the major result of Lochhaas et al. (2018). The density of the
shocked wind is also relatively constant across its width as assumed by the
analytic models, though that density drops as the shocked wind expands.

Instead of a layer of cold unshocked gas, the gas behind the shocked wind is
very hot ($\gtrsim 5\times 10^{6}$ K) and chaotic thanks to successive
feedback events. Additionally, these continued outflows prevent the outermost
shocked wind from stalling and falling back onto the galaxy, as seen in the
additional tests run for Section 5.3 of Lochhaas et al. (2018). They also
appear to disrupt any dense, cooling material that may fall inward towards the
galaxy. Lochhaas et al. (2018) note that such gas could be seeded by a
Rayleigh-Taylor or Vishniac (Vishniac, 1983; Vishniac & Ryu, 1989) thin-shell
instability. While we see plenty of turbulence, we see no evidence of
infalling material reaching the galactic disk outside of the disk midplane.

###  6.3 The Role of the Toomre Criterion

Figure 13: Toomre parameter $Q$ as a function of radius for our simulation
initial conditions. We calculate $Q$ analytically based on the static
background gravitational potentials and the initial disk density profile
(Section 2.1.1). The disk is initially isothermal with $T=10^{5}$ K (orange),
but other profiles are shown assuming the disk cools uniformly to $10^{4}$ K
(green) and $10^{3}$ K (blue). The latter is just below the minimum
temperature for star particles to form ($3\times 10^{3}$ K). The vertical
dashed line marks the edge of the fiducial disk, based on how it blends with
the CGM.

One of the most striking features of Figure 4 is that the disk gas mass is
almost identical between our variants with feedback. This is true at all
times, both during the $t<1$ Gyr star formation burst and late into the
simulations, despite vast differences in stellar mass. This feature is
consistent with the Toomre criterion $Q$ for disk fragmentation

$Q=\frac{c_{s}\kappa}{\pi G\Sigma_{\rm g}}$ (7)

where $c_{s}$ is the sound speed, $\kappa$ is the epicyclic frequency of the
rotating disk, and $\Sigma_{g}$ is the gas surface density of the disk’s face
(Toomre, 1964). If $Q>1$ at a given disk radius, the gas is stable against
fragmentation. If $Q$ locally drops below 1, then at that radius, the disk is
able to fragment into clumps that lead to star formation.

Important for our discussion is the dependence of $Q$ on the gas surface
density, $\Sigma_{\rm g}$. As gas accretes, the surface density rises and thus
$Q$ drops. Star formation lowers $\Sigma_{\rm g}$, bringing $Q$ back towards 1
and the disk back towards stability. Gas accretion and fragmentation therefore
conspire to keep $Q\sim 1$, which results in a fairly stable gas surface
density and, by extension, overall mass. All disks are constructed to have the
same initial $\Sigma_{g}$ profile through Equations 1 and 2 (though the
HighRatio and LowRatio disks extend to larger or smaller radii, respectively).
They also all evolve within the same static gravitational potential, which
sets $\kappa$, and remain roughly isothermal with $T\sim 4\times 10^{3}$ K
(although $T\sim 10^{2}$ K in the densest clumps), setting $c_{s}$. The
similarity in disk temperature and sound speed between all of the variants
indicates that the Toomre criterion is involved in keeping all variants near
the same disk mass. The exception to this is the CoolFlow variant, whose disk
mass notably deviates from the other variants despite having only a slightly
warmer disk overall ($T\sim 5\times 10^{3}$ K). The lower disk gas mass in the
CoolFlow variant may occur because there is no supernova feedback to generate
turbulence, which would raise the effective value of $c_{s}$ in Equation 7.

Figure 13 shows the $Q$ parameter calculated as a function of radius for the
simulation initial conditions (Section 2.1.1). We perform this calculation
analytically (i.e., we do not estimate $Q$ from the simulation itself) using
the background stellar and dark matter potentials and the initial disk density
profile (Equations 1 and 2) to calculate the epicyclic frequency $\kappa$ and
surface density $\Sigma_{\rm g}$, respectively. Since the initial disk is
isothermal at $10^{5}$ K, the sound speed $c_{s}$ is constant. We also include
$Q$ calculated for lower temperatures, assuming the disk cools uniformly:
$10^{4}$ and $10^{3}$ K. The latter is close to the minimum temperature
required for stars to form ($3\times 10^{3}$ K) and is therefore an
approximate representation of $Q(r)$ when stars begin to form. We only have
$Q\sim 1$ within the inner $\sim 10$ kpc, which is commensurate with the
average star formation radii in Table 3 (with exception of the CoolFlow
variant).

###  6.4 The Timescale of Self-Regulation

Broadly speaking, self-regulation is a balance between gas inflow and
feedback. This definition is used in, e.g., Mitchell & Schaye (2021) and
Prasad et al. (2020). Self-regulation can be understood through analogy to a
thermostat, which works to keep the temperature of a room at a specified “set
point.” Self-regulation in a galaxy should work to keep the star formation
rate near a particular set point. What this set point is, exactly, is not set
externally (as with a thermostat), but would be set by the conditions of the
galaxy and its environment. Thermostats turn off and on a furnace, which in
our analogy, is the star formation. The product of a furnace is heat, and the
product of our galactic “furnace” is stellar feedback.

Real thermostats do not keep a room at the set point all the time. When the
temperature drops, the furnace is turned on. When the room gets sufficiently
warm (either reaching the set point exactly or maybe slightly above to
compensate for cooling), the furnace is shut off. We expect the same thing of
galaxies: gas accretion should lead to star formation. The feedback from this
star formation will influence the galaxy’s ability to accrete gas, modulating
future star formation.

The connection between gas accretion and star formation is not temporally or
spatially immediate, however. As we discussed in Section 6.2.2, gas is often
seen to accrete along the plane of the disk, but star formation is often
concentrated at the center of the disk where densities are highest. Recent
estimates put the inward radial mass flow rate at less than 1
$\mathrm{M_{\odot}}/\mathrm{yr}$ (Di Teodoro & Peek, 2021). During a feedback
burst (when the furnace is on), it is therefore probably not the case that
inflow is perfectly balanced by star formation. Rather, the two balance each
other on average over some timescale.

With Figure 8, we attempted to identify by eye correlations between the SFR
and the net rate of change in disk mass, $\dot{M}_{\rm disk}$, but ultimately
found no consistent relation. While there are certainly rigorous statistical
tools for addressing time correlations, this analysis is sufficient to raise
questions about what precisely one should be looking for in order to decide if
a galaxy is self-regulating its star formation.

The ambiguity over timescale is exacerbated when examining self-regulation
with simulations. Because we are limited in resolution and computational
power, simulations of Milky Way-mass galaxies must use particles to represent
entire stellar populations rather than individual stars. This modeling
limitation introduces shot noise into the SFR. It is therefore also worth
considering the timescales over which we are concerned with changes in the
SFR.

The timescale uncertainties are in large part due to the temporal and spatial
separation between gas accretion and feedback. Within our simulations, CGM gas
that accretes onto the galaxy is more likely to contribute to a star particle
if it reaches the denser disk center. The formation of a star particle
obfuscates the physical and temporal scales involved in the actual formation
of molecular clouds and protostars, but these are assumed to be below the
temporal resolution of the simulation. Connections between inflow and star
formation in simulations seem clearer when made on smaller physical scales,
such as gas inflow along dust lanes at the very nucleus of a barred spiral
galaxy (Moon et al., 2022).

Our efforts to directly explore the self-regulation of galaxy star formation
has highlighted the need to refine our questions. Specifically, if we as a
community are to define self-regulation as a balance between gas accretion and
star formation, we should consider on what timescale we expect this balance to
be achieved. This way, we may better identify which physical processes
contribute to the self-regulation of galaxies.

###  6.5 Missing Model Components

One of the key takeaways from Prasad et al. (2020) and Voit et al. (2020) is
that feedback from Type Ia supernovae assists an AGN in self-regulating its
own feedback. Similarly, it seems that our simulations are missing one or more
important features that work in concert with Type II supernova feedback to
regulate a Milky Way-like galaxy. While there may be evidence for weak
coupling between gas inflow and star formation, as discussed in Section 6.4,
we believe that our simulations are likely not capturing physical features or
effects that would make this coupling more obvious and/or tighter.

The biggest evidence for this is the large feedback-driven bubbles that are
present in all of our simulation variants with stellar feedback (i.e.,
excluding the CoolFlow variant). Though the ramp in feedback efficiency,
described in Section 2.2.1, delays the onset of these bubbles, we don’t expect
further tweaks to the ramp (such as modifications to its start and end points
in both time and efficiency values) to significantly affect the presence of
these large bubbles. This is because the dramatic cavities seen in at late
times in Figures 6 and 7 are due to successive bubbles being blown into the
same volume. Earlier feedback winds clear away material, making it easier for
later winds to travel faster and farther, and to further remove material from
the bubble’s cavity.

These bubbles travel well beyond the virial radius by the end of the
simulations at $t=4$ Gyr, leaving cavities of high entropy gas that will not
cool within a Hubble time. The density of the ambient medium continues to
decrease beyond the virial radius (Figure 1), providing nothing for the
feedback-driven bubbles to meaningfully collide with and halt against.
Additionally, the cooling time at large radii is high enough that material
cannot “backfill” the bubble cavities. The behavior of these outflows suggest
that we are missing either material that the winds would stall against, a
mechanism which would fill the cavities in with low entropy gas, or both.

The idealized nature of our simulations may be working against us in this
regard. Fielding et al. (2020) compared idealized and cosmological simulations
of Milky Way-like galaxies and found that the outer CGM structure of the
latter ($\gtrsim 0.5r_{200}$) was highly impacted by cosmological effects such
as nonspherical gas accretion and the presence of satellites. No attempt to
model cosmological accretion was included in our simulations, but these could
address the large wind-driven bubbles in our simulations. This is especially
true considering that the FOGGIE simulations (Peeples et al., 2019) use the
same feedback algorithm and a slightly higher feedback efficiency but do not
see the same long-term disruption due to outflows. Appropriately modeling
inflow in an idealized simulation is a challenge, however. Typically,
cosmological accretion in idealized simulations is treated as being
spherically symmetric, as in Fielding et al. (2017) (albeit with added density
fluctuations), but as stated above, Fielding et al. (2020) emphasized that
this spherical treatment is not sufficient. That work highlighted the strength
of idealized simulations for studying the inner CGM ($\lesssim 0.5r_{200}$),
but our work indicates that for self-regulation the outer CGM and perhaps the
nearby IGM exert an important influence.

Our initial CGM is also smooth and spherically symmetric. A more realistic CGM
would be properly multiphase, with fluctuations in density, temperature, and
velocity (both radial and tangentially). We initially expected to see such
fluctuations develop over the course of the simulation as a result of
outflows, but including these from the start would likely have a profound
effect. A nonuniform CGM will have a distribution of cooling times at a given
radius, better mixing, and may be better able to disrupt early outflows. In
the simulations presented in this work, the earliest outflows are able to
expand through the CGM nearly uniformly.

Initial density perturbations are the easiest way to disrupt this spherical
symmetry, but require a choice of power spectrum. A more natural way of
introducing perturbations at large radii may occur from adjusting our CGM
beyond $r_{200}$. Our current treatment for this “outer” gas was motivated by
practical considerations rather than observations. Notably, our entropy
initial profile increases monotonically with radius. More realistic halos are
likely to have low-entropy gas beyond $r_{200}$: gas that has either not yet
passed through a cosmological accretion shock, or has passed through the shock
with enough density for its current entropy to be less than the mean at
$r_{200}$. The relaxation of the simulation’s initial state could push this
gas to become Rayleigh-Taylor unstable, and either over time or with the
assistance of instabilities seeded in the initial conditions, could grow dense
clumps of gas that may affect the outflows and overall precipitation within
our simulations.

Cosmological inflow would further assist in seeding and/or maintaining CGM
fluctuations, naturally creating them as infalling material combines with the
galaxy. This is true both for filamentary inflow as well as the presence of
and mergers with small satellites and dwarf galaxies.

The failure of our simulations to “close the feedback loop” is not entirely
novel. Prasad et al. (2020) simulated an elliptical galaxy of similar mass
($2\times 10^{12}~{}\mathrm{M_{\odot}}$), dominated by AGN feedback, and also
saw a highly disrupted CGM emerge. Generally, lower mass galaxies have a
weaker gravitational potential and lower CGM pressure. Outflows can more
easily escape the halo. This reinforces our earlier observation that some
physical mechanism at the edge of the halo—cosmological inflows or simply more
generic density perturbations, to name two candidates—appears necessary for
lower mass halos to maintain their CGM structure.

Though we suspect the lack of cosmological effects is the biggest missing
piece from our simulations, there are other physical processes to consider.
Our current simulations use Type II supernovae as their only feedback source.
Though our simulations do not currently suffer from a lack of outflows,
efforts to counteract their current behavior may highlight the need to include
other forms of feedback, such as Type Ia supernovae and AGN. We adopt a fairly
straightforward prescription for rotation in the CGM and assume its angular
momentum is aligned with that of the disk, but the angular momentum of the CGM
is likely quite complicated (Cadiou et al., 2021). This could facilitate
better mixing between outflows and the ambient medium, disrupting the
structures that develop in our simulations.

Finally, our simulations omit two important plasma components: magnetic fields
and cosmic rays. These two influences would significantly alter the behavior
of outflows. Magnetic fields can slow outflows and raise the density of the
inner CGM, though they also hinder metal mixing (van de Voort et al., 2021).
Cosmic rays provide a form of non-thermal pressure support, allowing cold gas
to occupy more CGM volume (Ji et al., 2020). Cosmic rays also lead to larger,
lower density cold clouds, and keep cold gas in the CGM for longer (Butsky et
al., 2020). All of these changes have implications for the galaxy’s accretion
of CGM gas.

Disruption is not a universal feature of CGM-focused idealized simulations.
Fielding et al. (2017) and Li & Tonnesen (2020) do not include a cold gas disk
nor explicitly model star formation. Instead, they tie outflows to the amount
of inflow through an inner boundary. In both cases, outflows do not cause a
large-scale disruption of the CGM. Li & Tonnesen (2020) see a clear net
balance between gas inflow and outflow in the latter half of their simulation,
and both works observe clear cold gas condensation near the galaxy ($r\lesssim
100$ kpc). Yet tying outflows directly to inflows, rather than depending on
the intermediary process of star formation, may under-predict the strength of
stellar feedback because it neglects star formation due to gas already present
in the disk. It also inputs feedback energy into a region that is spatially
removed from star formation, where gas is more likely to be lower density.

##  7 Conclusions

We have run a suite of isolated, idealized Milky Way-like galaxy simulations
in order to examine the ability of galaxies to self-regulate their star
formation. They were specifically designed to explore the precipitation theory
of self-regulation (Voit et al., 2015). The circumgalactic medium (CGM) in our
galaxies was initialized in hydrostatic equilibrium with entropy profiles set
by expectations from precipitation (Voit, 2019). The CGM was also given an
initial azimuthal rotation scaled off the estimates of Hodges-Kluck et al.
(2016). We explored variations in the entropy profile through the
precipitation limit parameter $\tau$ in Equation 2.1.2 as well as variations
in the rotation profile of the CGM (Equation 5).

Ours are the first idealized simulations to include both the CGM and explicit
star formation in a rotating disk. Previous idealized CGM simulations
(Fielding et al., 2017; Li & Tonnesen, 2020) tied outflow directly to inflow
through an inner boundary. Explicitly modeling star formation brings
additional challenges, as the cold gas disk typically used as an initial
condition in isolated galaxy models leads to large initial burst of star
formation. Outflows from the resulting feedback quickly disrupt the CGM. This
has historically not been a challenge for isolated galaxy simulations (e.g.,
Kim et al., 2016; Benincasa et al., 2016) because they have had an essentially
non-existent CGM. To prevent the initial star formation burst from disrupting
our CGM before its gas can accrete onto the disk, we implement a ramp in the
stellar feedback’s efficiency parameter (Section 2.2.1). This ramp minimizes
the impact of feedback for the first 1 Gyr. As a result, the impacts of our
CGM variations become apparent.

Our work highlights that including explicit star formation in idealized
simulations is crucial to understanding the galaxy-CGM connection. Alternative
solutions, such as tying outflow rates directly to gas accretion rates,
obfuscate important steps connecting accretion and feedback. We ultimately
fail to produce isolated galaxies that are able to self-regulate their star
formation, but this failure is illuminating in several ways. The primary
results of our simulations are as follows:

  * •

Idealized galaxy simulations are highly sensitive to their initial conditions.
Our simulation setup is not unusual among isolated disk galaxy simulations,
but complications arise when including the CGM. Chief among these is a large
initial burst of star formation that can disrupt the CGM.

  * •

Even after mitigating the initial star formation burst, our simulations
contain outflows that are very disruptive to the CGM. Disruptive feedback
seems to be a common feature of isolated galaxy simulations at this halo mass
(Prasad et al., 2020). This indicates that current idealized simulations are
missing important features that would constrain, disrupt, or backfill these
outflows, such as cosmological inflow.

  * •

Our galaxies continue to accrete gas along the midplane of the disk, despite
the disruption of the CGM by outflows. Though this accretion channel contains
gas with the lowest $t_{\rm cool}/t_{\rm ff}$ ratio, it is still higher than
expected for a precipitation-regulated system. This accretion is able to
maintain low SFRs of $\sim 0.1\ \mathrm{M_{\odot}}$/yr.

  * •

Rotation in the CGM impacts the ability of gas to accrete onto the disk. This
accretion also varies with the rotation profile. Rotation is therefore an
important component of idealized CGM studies, and better understanding angular
momentum in the CGM is an important prerequisite to informative modeling.

Understanding the balance between accretion rate and star formation first
requires an understanding of the time scales over which we expect these
processes to balance. While our work is not able to answer this question, we
find it important for the community to consider as studies of the CGM’s impact
on star formation continue.

Future work will incorporate some of the physical features we have identified
as potentially important for mitigating the dominance of outflows in our
idealized simulations. We also wish to investigate the importance of magnetic
fields and cosmic rays, particularly since the latter can dramatically affect
the CGM’s cold gas fraction (Ji et al., 2020; Butsky & Quinn, 2018).

CK was supported in part by a Department of Energy Computational Science
Graduate Fellowship under Award Number DE-FG02-97ER25308. CK and BWO
acknowledge support from NSF grant #1908109, and BWO acknowledges further
support from NASA ATP grants NNX15AP39G and 80NSSC18K1105. GMV and BWO
acknowledge support from NSF grant #2106575. This work used the Extreme
Science and Engineering Discovery Environment (XSEDE) under allocation TG-
AST090040, as well as the resources of the Michigan State University High
Performance Computing Center (operated by the Institute for Cyber-Enabled
Research).

## References

  * Babyk et al. (2018) Babyk, I. V., McNamara, B. R., Nulsen, P. E. J., et al. 2018, The Astrophysical Journal, 862, 39, doi: 10.3847/1538-4357/aacce5
  * 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
  * Benincasa et al. (2016) Benincasa, S. M., Wadsley, J., Couchman, H. M. P., & Keller, B. W. 2016, Monthly Notices of the Royal Astronomical Society, 462, 3053, doi: 10.1093/mnras/stw1741
  * Berger & Colella (1989) Berger, M. J., & Colella, P. 1989, Journal of Computational Physics, 82, 64, doi: 10.1016/0021-9991(89)90035-1
  * Bouché et al. (2010) Bouché, N., Dekel, A., Genzel, R., et al. 2010, The Astrophysical Journal, 718, 1001, doi: 10.1088/0004-637X/718/2/1001
  * Brummel-Smith et al. (2019) Brummel-Smith, C., Bryan, G., Butsky, I., et al. 2019, Journal of Open Source Software, 4, 1636, doi: 10.21105/joss.01636
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, The Astrophysical Journal Supplement Series, 211, 19, doi: 10.1088/0067-0049/211/2/19
  * Butsky et al. (2020) Butsky, I. S., Fielding, D. B., Hayward, C. C., et al. 2020, The Astrophysical Journal, 903, 77, doi: 10.3847/1538-4357/abbad2
  * Butsky & Quinn (2018) Butsky, I. S., & Quinn, T. R. 2018, The Astrophysical Journal, 868, 108, doi: 10.3847/1538-4357/aaeac2
  * Cadiou et al. (2021) Cadiou, C., Dubois, Y., & Pichon, C. 2021, doi: 10.48550/arXiv.2110.05384
  * Cavagnolo et al. (2008) Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2008, The Astrophysical Journal, 683, L107, doi: 10.1086/591665
  * Cen & Ostriker (1992) Cen, R., & Ostriker, J. P. 1992, The Astrophysical Journal Letters, 399, L113, doi: 10.1086/186620
  * Corlies et al. (2020) Corlies, L., Peeples, M. S., Tumlinson, J., et al. 2020, The Astrophysical Journal, 896, 125, doi: 10.3847/1538-4357/ab9310
  * Côté et al. (2016) Côté, B., Ritter, C., O’Shea, B. W., et al. 2016, The Astrophysical Journal, 824, 82, doi: 10.3847/0004-637X/824/2/82
  * Davé et al. (2012) Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, Monthly Notices of the Royal Astronomical Society, 421, 98, doi: 10.1111/j.1365-2966.2011.20148.x
  * DeFelippis et al. (2020) DeFelippis, D., Genel, S., Bryan, G. L., et al. 2020, The Astrophysical Journal, 895, 17, doi: 10.3847/1538-4357/ab8a4a
  * Di Teodoro & Peek (2021) Di Teodoro, E. M., & Peek, J. E. G. 2021, The Astrophysical Journal, 923, 220, doi: 10.3847/1538-4357/ac2cbd
  * Donahue & Voit (2022) Donahue, M., & Voit, G. M. 2022, Physics Reports, 973, 1, doi: 10.1016/j.physrep.2022.04.005
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Revista Mexicana de Astronomia y Astrofisica, 49, 137.  http://adsabs.harvard.edu/abs/2013RMxAA..49..137F
  * Fielding et al. (2017) Fielding, D., Quataert, E., McCourt, M., & Thompson, T. A. 2017, Monthly Notices of the Royal Astronomical Society, 466, 3810, doi: 10.1093/mnras/stw3326
  * Fielding et al. (2020) Fielding, D. B., Tonnesen, S., DeFelippis, D., et al. 2020, The Astrophysical Journal, 903, 32, doi: 10.3847/1538-4357/abbc6d
  * Frisbie et al. (2020) Frisbie, R. L. S., Donahue, M., Voit, G. M., et al. 2020, The Astrophysical Journal, 899, 159, doi: 10.3847/1538-4357/aba8a8
  * Gaspari et al. (2015) Gaspari, M., Brighenti, F., & Temi, P. 2015, Astronomy & Astrophysics, 579, A62, doi: 10.1051/0004-6361/201526151
  * Gaspari et al. (2013) Gaspari, M., Ruszkowski, M., & Oh, S. P. 2013, Monthly Notices of the Royal Astronomical Society, 432, 3401, doi: 10.1093/mnras/stt692
  * Gaspari et al. (2012) Gaspari, M., Ruszkowski, M., & Sharma, P. 2012, The Astrophysical Journal, 746, 94, doi: 10.1088/0004-637X/746/1/94
  * Haardt & Madau (2012) Haardt, F., & Madau, P. 2012, The Astrophysical Journal, 746, 125, doi: 10.1088/0004-637X/746/2/125
  * Hafen et al. (2022) Hafen, Z., Stern, J., Bullock, J., et al. 2022, arXiv e-prints.  https://arxiv.org/abs/2201.07235
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Ho et al. (2019) Ho, S. H., Martin, C. L., & Turner, M. L. 2019, The Astrophysical Journal, 875, 54, doi: 10.3847/1538-4357/ab0ec2
  * Hodges-Kluck et al. (2016) Hodges-Kluck, E. J., Miller, M. J., & Bregman, J. N. 2016, The Astrophysical Journal, 822, 21, doi: 10.3847/0004-637X/822/1/21
  * Hummels et al. (2019) Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, The Astrophysical Journal, 882, 156, doi: 10.3847/1538-4357/ab378f
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ji et al. (2020) Ji, S., Chan, T. K., Hummels, C. B., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 4221, doi: 10.1093/mnras/staa1849
  * Kim et al. (2016) Kim, J.-h., Agertz, O., Teyssier, R., et al. 2016, The Astrophysical Journal, 833, 202, doi: 10.3847/1538-4357/833/2/202
  * Leroy et al. (2008) Leroy, A. K., Walter, F., Brinks, E., et al. 2008, The Astronomical Journal, 136, 2782, doi: 10.1088/0004-6256/136/6/2782
  * Li & Tonnesen (2020) Li, M., & Tonnesen, S. 2020, The Astrophysical Journal, 898, 148, doi: 10.3847/1538-4357/ab9f9f
  * Li & Bryan (2014a) Li, Y., & Bryan, G. L. 2014a, The Astrophysical Journal, 789, 153, doi: 10.1088/0004-637X/789/2/153
  * Li & Bryan (2014b) —. 2014b, The Astrophysical Journal, 789, 54, doi: 10.1088/0004-637X/789/1/54
  * Lochhaas et al. (2020) Lochhaas, C., Bryan, G. L., Li, Y., Li, M., & Fielding, D. 2020, Monthly Notices of the Royal Astronomical Society, 493, 1461, doi: 10.1093/mnras/staa358
  * Lochhaas et al. (2018) Lochhaas, C., Thompson, T. A., Quataert, E., & Weinberg, D. H. 2018, Monthly Notices of the Royal Astronomical Society, 481, 1873, doi: 10.1093/mnras/sty2421
  * Lopez et al. (2020) Lopez, S., Tejos, N., Barrientos, L. F., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 4442, doi: 10.1093/mnras/stz3183
  * McCourt et al. (2012) McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, Monthly Notices of the Royal Astronomical Society, 419, 3319, doi: 10.1111/j.1365-2966.2011.19972.x
  * McKinney (2010) McKinney, W. 2010, in Python in Science Conference, Austin, Texas, 56–61, doi: 10.25080/Majora-92bf1922-00a
  * Mitchell & Schaye (2021) Mitchell, P. D., & Schaye, J. 2021, arXiv e-prints.  https://arxiv.org/abs/2103.10966
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, Publications of the Astronomical Society of Japan, 27, 533.  http://adsabs.harvard.edu/abs/1975PASJ...27..533M
  * Moon et al. (2022) Moon, S., Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2022, The Astrophysical Journal, 925, 99, doi: 10.3847/1538-4357/ac3a7b
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, The Astrophysical Journal, 462, 563, doi: 10.1086/177173
  * Oh et al. (2020) Oh, B. K., Smith, B. D., Peacock, J. A., & Khochfar, S. 2020, Monthly Notices of the Royal Astronomical Society, 497, 5203, doi: 10.1093/mnras/staa2318
  * Oppenheimer (2018) Oppenheimer, B. D. 2018, Monthly Notices of the Royal Astronomical Society, 480, 2963, doi: 10.1093/mnras/sty1918
  * Peeples et al. (2014) Peeples, M. S., Werk, J. K., Tumlinson, J., et al. 2014, The Astrophysical Journal, 786, 54, doi: 10.1088/0004-637X/786/1/54
  * Peeples et al. (2019) Peeples, M. S., Corlies, L., Tumlinson, J., et al. 2019, The Astrophysical Journal, 873, 129, doi: 10.3847/1538-4357/ab0654
  * Pizzolato & Soker (2005) Pizzolato, F., & Soker, N. 2005, The Astrophysical Journal, 632, 821, doi: 10.1086/444344
  * Pizzolato & Soker (2010) —. 2010, Monthly Notices of the Royal Astronomical Society, 408, 961, doi: 10.1111/j.1365-2966.2010.17156.x
  * Popesso et al. (2022) Popesso, P., Concas, A., Cresci, G., et al. 2022, arXiv e-prints.  https://arxiv.org/abs/2203.10487
  * Prasad et al. (2015) Prasad, D., Sharma, P., & Babul, A. 2015, The Astrophysical Journal, 811, 108, doi: 10.1088/0004-637X/811/2/108
  * Prasad et al. (2020) Prasad, D., Voit, G. M., O’Shea, B. W., & Glines, F. 2020, The Astrophysical Journal, 905, 50, doi: 10.3847/1538-4357/abc33c
  * Prochaska et al. (2017) Prochaska, J. X., Werk, J. K., Worseck, G., et al. 2017, The Astrophysical Journal, 837, 169, doi: 10.3847/1538-4357/aa6007
  * Renzini & Peng (2015) Renzini, A., & Peng, Y.-j. 2015, The Astrophysical Journal Letters, 801, L29, doi: 10.1088/2041-8205/801/2/L29
  * Schaye et al. (2010) Schaye, J., Vecchia, C. D., Booth, C. M., et al. 2010, Monthly Notices of the Royal Astronomical Society, 402, 1536, doi: 10.1111/j.1365-2966.2009.16029.x
  * Sharma et al. (2012) Sharma, P., McCourt, M., Quataert, E., & Parrish, I. J. 2012, Monthly Notices of the Royal Astronomical Society, 420, 3174, doi: 10.1111/j.1365-2966.2011.20246.x
  * Sherman et al. (2021) Sherman, S., Jogee, S., Florez, J., et al. 2021, Monthly Notices of the Royal Astronomical Society, 505, 947, doi: 10.1093/mnras/stab1350
  * Shin et al. (2021) Shin, E.-J., Kim, J.-H., & Oh, B. K. 2021, The Astrophysical Journal, 917, 12, doi: 10.3847/1538-4357/abffd0
  * Smith et al. (2011) Smith, B. D., Hallman, E. J., Shull, J. M., & O’Shea, B. W. 2011, The Astrophysical Journal, 731, 6, doi: 10.1088/0004-637X/731/1/6
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, Monthly Notices of the Royal Astronomical Society, 466, 2217, doi: 10.1093/mnras/stw3291
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, Annual Review of Astronomy and Astrophysics, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Stern et al. (2021) Stern, J., Faucher-Giguère, C.-A., Fielding, D., et al. 2021, The Astrophysical Journal, 911, 88, doi: 10.3847/1538-4357/abd776
  * Su et al. (2020) Su, K.-Y., Hopkins, P. F., Hayward, C. C., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 1190, doi: 10.1093/mnras/stz3011
  * Suresh et al. (2019) Suresh, J., Nelson, D., Genel, S., Rubin, K., & Hernquist, L. 2019, Monthly Notices of the Royal Astronomical Society, 483, 4040, doi: 10.1093/mnras/sty3402
  * Tange (2020) Tange, O. 2020, GNU Parallel 20200822 (’Beirut’), Zenodo, doi: 10.5281/zenodo.3996295
  * Tejos et al. (2021) Tejos, N., López, S., Ledoux, C., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 663, doi: 10.1093/mnras/stab2147
  * Tonnesen & Bryan (2009) Tonnesen, S., & Bryan, G. L. 2009, The Astrophysical Journal, 694, 789, doi: 10.1088/0004-637X/694/2/789
  * Toomre (1964) Toomre, A. 1964, The Astrophysical Journal, 139, 1217, doi: 10.1086/147861
  * Trapp et al. (2022) Trapp, C. W., Kereš, D., Chan, T. K., et al. 2022, Monthly Notices of the Royal Astronomical Society, 509, 4149, doi: 10.1093/mnras/stab3251
  * Truong et al. (2021) Truong, N., Pillepich, A., Nelson, D., Werner, N., & Hernquist, L. 2021, Monthly Notices of the Royal Astronomical Society, 508, 1563, doi: 10.1093/mnras/stab2638
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, Annual Review of Astronomy and Astrophysics, 55, 389, doi: 10.1146/annurev-astro-091916-055240
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, The Astrophysical Journal Supplement Series, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * van de Voort et al. (2021) van de Voort, F., Bieri, R., Pakmor, R., et al. 2021, Monthly Notices of the Royal Astronomical Society, 501, 4888, doi: 10.1093/mnras/staa3938
  * van de Voort et al. (2019) van de Voort, F., Springel, V., Mandelker, N., van den Bosch, F. C., & Pakmor, R. 2019, Monthly Notices of the Royal Astronomical Society, 482, L85, doi: 10.1093/mnrasl/sly190
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vishniac (1983) Vishniac, E. T. 1983, The Astrophysical Journal, 274, 152, doi: 10.1086/161433
  * Vishniac & Ryu (1989) Vishniac, E. T., & Ryu, D. 1989, The Astrophysical Journal, 337, 917, doi: 10.1086/167161
  * Voit (2018) Voit, G. M. 2018, The Astrophysical Journal, 868, 102, doi: 10.3847/1538-4357/aae8e2
  * Voit (2019) —. 2019, The Astrophysical Journal, 880, 139, doi: 10.3847/1538-4357/ab2bfd
  * Voit (2021) —. 2021, The Astrophysical Journal, 908, L16, doi: 10.3847/2041-8213/abe11f
  * Voit et al. (2015) Voit, G. M., Bryan, G. L., O’Shea, B. W., & Donahue, M. 2015, The Astrophysical Journal Letters, 808, L30, doi: 10.1088/2041-8205/808/1/L30
  * Voit et al. (2008) Voit, G. M., Cavagnolo, K. W., Donahue, M., et al. 2008, The Astrophysical Journal, 681, L5, doi: 10.1086/590344
  * Voit & Donahue (2011) Voit, G. M., & Donahue, M. 2011, The Astrophysical Journal, 738, L24, doi: 10.1088/2041-8205/738/2/L24
  * Voit & Donahue (2015) —. 2015, The Astrophysical Journal, 799, L1, doi: 10.1088/2041-8205/799/1/L1
  * Voit et al. (2019) Voit, G. M., Donahue, M., Zahedy, F., et al. 2019, The Astrophysical Journal Letters, 879, L1, doi: 10.3847/2041-8213/ab2766
  * Voit et al. (2017) Voit, G. M., Meece, G., Li, Y., et al. 2017, The Astrophysical Journal, 845, 80, doi: 10.3847/1538-4357/aa7d04
  * Voit et al. (2020) Voit, G. M., Bryan, G. L., Prasad, D., et al. 2020, The Astrophysical Journal, 899, 70, doi: 10.3847/1538-4357/aba42e
  * Zaragoza-Cardiel et al. (2019) Zaragoza-Cardiel, J., Fritz, J., Aretxaga, I., et al. 2019, Monthly Notices of the Royal Astronomical Society, 487, L61, doi: 10.1093/mnrasl/slz093