# An Enhanced Massive Black Hole Occupation Predicted in Cluster Dwarf
Galaxies

M. Tremmel1 A. Ricarte2,3 P. Natarajan2,4,5 J. Bellovary6,7 R. Sharma8 T. R.
Quinn9 1 School of Physics, University College Cork, Cork, Ireland  2 Black
Hole Initiative, 20 Garden Street, Cambridge, MA 02138, USA  3 Center for
Astrophysics – Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  4 Department of Astronomy, Yale University, 52 Hillhouse Avenue New
Haven, CT 06511, USA
5Department of Physics, Yale University, P.O. Box 208121, New Haven, CT 06511,
USA  6 Department of Physics, Queensborough Community College, City University
of New York, 222-05 56th Ave, Bayside, NY, 11364  7 Department of
Astrophysics, American Museum of Natural History, Central Park West at 79th
Street, New York, NY 10024, USA  8Department of Physics & Astronomy, Rutgers,
The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ
08854, USA  9 Astronomy Department, University of Washington, Box 351580,
Seattle, WA, 98195-1580  Corresponding Author<EMAIL_ADDRESS>
###### Abstract

The occupation fraction of massive black holes (MBHs) in low mass galaxies
offers interesting insights into initial black hole seeding mechanisms and
their mass assembly history, though disentangling these two effects remains
challenging. Using the Romulus cosmological simulations we examine the impact
of environment on the occupation fraction of MBHs in low mass galaxies. Unlike
most modern cosmological simulations, Romulus seeds MBHs based on local gas
properties, selecting dense ($n>3$ cm-3, 15 times the threshold for star
formation), pristine ($Z<3\times 10^{-4}Z_{\odot}$), and rapidly collapsing
regions in the early Universe as sites to host MBHs without assuming anything
about MBH occupation as a function of galaxy stellar mass, or halo mass, a
priori. The simulations predict that dwarf galaxies with M${}_{\star}<10^{9}$
M⊙ in cluster environments are $\sim 2$ times more likely to host a MBH
compared to those in the field. The predicted occupation fractions are
remarkably consistent with those of nuclear star clusters. Across cluster and
field environments, dwarf galaxies with earlier formation times are more
likely to host a MBH. Thus, while the MBH occupation function is similar
between cluster and field environments at high redshift ($z>3$), a difference
arises as late-forming dwarfs – which do not exist in the cluster environment
– begin to dominate in the field and pull the MBH occupation fraction down for
low mass galaxies. Additionally, prior to in-fall some cluster dwarfs are
similar to progenitors of massive, isolated galaxies, indicating that they
might have grown to higher masses had they not been impeded by the cluster
environment. While the population of MBHs in dwarf galaxies is already widely
understood to be important for understanding MBH formation, this work
demonstrates that environmental dependence is important to consider as future
observations search for low mass black holes in dwarf galaxies.

###### keywords:

galaxies:dwarf – black hole physics – galaxies:clusters

##  1 Introduction

The origin and evolution of massive black holes (MBHs) of mass $>10^{5}$ M⊙,
which are ubiquitously found in the centers of massive galaxies (Tremaine et
al., 2002; Kormendy & Richstone, 1995; Kormendy & Ho, 2013), remains an
important open question. There are different models for how the seeds of MBHs
form in the early Universe (Volonteri, 2010; Natarajan, 2014) which predict
different formation rates and initial masses. Seeds formed through the stellar
evolution of population III stars would have initial masses on the order of
$\sim 100$ M⊙ and would exist in virtually every galaxy, with a fraction
capable of attaining much larger masses through, e.g. super-Eddington
accretion events (Volonteri & Rees, 2005; Alexander & Natarajan, 2014;
Inayoshi et al., 2016; Sassano et al., 2023). Other formation models require
more specific conditions at high redshift, such as very high densities and a
lack of both metals and molecular gas, but can produce seed black holes of
mass $10^{4}-10^{6}$ M⊙ (Natarajan, 2011; Alexander & Natarajan, 2014). The
collapse of a dense star cluster can produce a very massive star through
runaway stellar mergers that then forms a massive black hole (Devecchi &
Volonteri, 2009; Davies et al., 2011). Alternatively, if fragmentation into a
star cluster is prevented, gas can collapse directly into a very massive
(quasi-) star type object, which then forms a massive black hole (Lodato &
Natarajan, 2006, 2007). The origins of MBHs are notoriously difficult to
constrain observationally, as the early seeding epochs are currently
inaccessible (though JWST is expected to change that soon); and typical
detection methods are best only at finding very massive and/or rapidly growing
MBHs, which effectively have their initial conditions mostly erased (Volonteri
et al., 2008; Volonteri & Gnedin, 2009). JWST and next-generation telescopes,
such as Euclid, have the potential to observe some of the earliest phases of
MBH growth (Sesana & Khan, 2015; Natarajan et al., 2017; Pacucci et al., 2019)
and could shine new light into their formation physics and environments.
Electromagnetic observations will be complemented by new gravitational wave
detectors like LISA (Laser Interferometer Space Antenna), with the ability to
detect merging, low-mass black holes out to $z>20$ (Colpi et al., 2019; Amaro-
Seoane et al., 2023). Constraining the MBH population at high redshift will
help to constrain the relative efficiency of different seeding mechanisms
(Sesana, 2013; Sesana & Khan, 2015; Ricarte & Natarajan, 2018) while providing
new insight into the mechanisms driving their growth and co-evolution with
their host (proto-)galaxies.

However, clues to MBH formation exist more locally as well. Dwarf galaxies
(with stellar mass M${}_{\star}<10^{10}$ M⊙) may host black holes that have
not grown significantly, neither through accretion nor mergers, and therefore
may maintain their connection to their initial conditions (Volonteri &
Natarajan, 2009; van Wassenhove et al., 2010; Greene, 2012; Ricarte &
Natarajan, 2018; Regan et al., 2023). There is a growing sample of MBHs
detected in dwarf galaxies using a variety of methods and wavelengths,
typically observing MBHs as low luminosity active galactic nuclei (AGN; e.g.
Reines et al., 2011; Reines & Deller, 2012; Reines et al., 2020; Baldassare et
al., 2015; Baldassare et al., 2016; Mezcua et al., 2018; Nguyen et al., 2019;
Woo et al., 2019; Birchall et al., 2020; Molina et al., 2021; Latimer et al.,
2021; Cann et al., 2021; Burke et al., 2022). Observations like these have
been used to estimate the true occupation fraction - the fraction of galaxies
hosting any MBH, regardless of whether it is detected as an AGN (Greene, 2012;
Miller et al., 2015; Nguyen et al., 2019; Burke et al., 2022; Askar et al.,
2022). While there is evidence that the occupation fraction does not
dramatically change down to stellar masses as low as $10^{9}$ M⊙ (Miller et
al., 2015; Baldassare et al., 2020; Burke et al., 2022), mapping detected AGN
fractions to the underlying population of MBHs in low mass galaxies remains
difficult and highly uncertain. JWST and next-generation telescopes like Vera
Rubin will also prove useful in expanding our view of MBHs in dwarf galaxies
(Baldassare et al., 2018; Cann et al., 2021; Burke et al., 2022).

Cosmological simulations, which self-consistently model the collapse and
merger history of dark matter halos with the baryonic evolution (gas
accretion, star formaiton, MBH growth) of galaxies, have proven to be an
invaluable tool to study the co-evolution of MBHs and galaxies (e.g. Di Matteo
et al., 2008, 2017; Okamoto et al., 2008; Sijacki et al., 2015; Rosas-Guevara
et al., 2016; Dubois et al., 2016; Nelson et al., 2019; Blank et al., 2019;
Habouzit et al., 2021; Koudmani et al., 2021, 2022; Ni et al., 2022). While
important work has been done to study the evolution of MBHs in dwarf galaxies
using large-scale simulations (e.g. Haidar et al., 2022), difficulties arise
due to a wide range of model assumptions which often have strong effects on
the results. Even the most modern large-scale simulations lack the resolution
to correctly model galaxies much lower than $10^{9}$ M⊙ in stellar mass. Many
smaller-scale simulations that can attain very high resolution, such as TNG50
(Nelson et al., 2019) or FABLE (Henden et al., 2018), utilize simplistic
prescriptions for seeding MBHs whereby all galaxies residing in dark matter
halos above a certain mass are seeded with a MBH at their centers. While this
may allow for predictions of the accretion rates and luminosities of MBHs in
low mass galaxies (i.e the active fraction), the underlying occupation
fraction of MBHs is an explicitly assumed prior. This type of prescription
will mean that for many low mass galaxies black holes are seeded at rather
late times instead of at high redshift like most theoretical models (although
recent works by Natarajan (2021) and Mayer et al. (2023) have noted that MBH
formation could continue throughout cosmic time and may not be limited to
metal poor regions). In addition, when considering satellite or backsplash
galaxies, simplistic schemes to force MBHs to the centers of galaxies can
generate unrealistic numerical effects and artificially impact the occupation
fraction in some environments (Borrow et al., 2022).

Recent high-resolution simulations that resolve dwarf galaxies while also
incorporating more predictive models for MBH formation, such as the Romulus
simulations (Tremmel et al., 2017, 2019), the Obelisk Simulations (Trebitsch
et al., 2021), the New Horizon simulations (Dubois et al., 2021; Beckmann et
al., 2023), and the MARVELous and DC Justice League simulations (Bellovary et
al., 2019, 2021; Munshi et al., 2021; Applebaum et al., 2021) are valuable
tools in predicting MBH evolution and occupation fraction within low mass
galaxies. In this Paper we utilize the Romulus suite of simulations to examine
the environmental dependence of MBH occupation in dwarf galaxies. Romulus
forms MBHs from dense, pristine gas in the early Universe ($z>5$) without any
a priori assumptions about which halos should or should not host a MBH. The
dynamics of MBHs is followed realistically down to sub-kpc separations, so
they are not always forced to the centers of galaxies (Tremmel et al., 2015).
The physics models governing MBH growth, star formation, and supernovae
feedback that have been incorporated into Romulus have been optimized to
broadly reproduce observed scaling relations across nearly five orders of
magnitude in halo mass (Tremmel et al., 2017, 2019; Ricarte et al., 2019;
Sharma et al., 2020).

MBHs in field dwarf galaxies have been extensively studied in Romulus, which
produces a realistic population of dwarf galaxy AGN consistent with
observations (Sharma et al., 2022a) that can also affect their evolution
(Sharma et al., 2020, 2022b), something which is seen in other simulations as
well (Koudmani et al., 2021, 2022). We expand on previous analysis in this
work by examining how the dwarf galaxy MBH occupation fraction evolves with
both time and environment. To do this we compare results from Romulus25, a 25
Mpc-per-side uniform volume simulation, with RomulusC, a zoom-in simulation of
a $10^{14}$ M⊙ galaxy cluster. Between these two simulations, we have several
hundred simulated dwarf galaxies with stellar masses in the range
$10^{7}-10^{10}$ M⊙ that are resolved with at least 200 baryonic resolution
elements (and $10^{4}$ for dark matter).

In Section 2 we provide a brief overview of the Romulus simulations and the
relevant physics models incorporated in them. In Section 3 we present our
results of the dependence of occupation fraction on environment and its
evolution with time. In Section 4 we explore the origins of the enhanced MBH
occupation fraction in $z=0$ cluster dwarf galaxies that we find. Section 5
discusses our results and Section 6 provides a summary of our conclusions.

##  2 The Romulus Simulations

In this section we briefly describe the relevant properties of the
simulations. For a more detailed discussion, including how the parameters were
chosen, we point the reader to Tremmel et al. (2017, 2019).

###  2.1 Overview of the Romulus Simulations

The Romulus Simulations (Tremmel et al., 2017, 2019) are a suite of
cosmological hydrodynamic simulations that includes a 25 Mpc-per-side uniform
volume simulation (Romulus25) and a zoom-in simulation of a $10^{14}$ M⊙
galaxy cluster (RomulusC). With a dark matter mass resolution of $3.4\times
10^{5}$ M⊙, both Romulus simulations are able to resolve halos as small as
$\sim 3\times 10^{9}$ M⊙ with more than 10,000 particles. With typical gas and
star particle masses of $2.1\times 10^{5}$ and $6\times 10^{4}$ M⊙
respectively, and a spline softening of 350 pc (equivalent to 250 pc plummer
softening), the two simulations are also able to resolve the the baryonic
structure of dwarf galaxies as small as $\sim 10^{7}$ M⊙ in stellar mass with
hundreds of resolution elements. Both simulations are run with the same
cosmology ($\Omega_{0}=0.3086,\Lambda=0.6914,h=0.6777,\sigma_{8}=0.8288$;
Planck Collaboration et al., 2016) and physics.

The Romulus simulations were run using the N-body+Smoothed particle
hydrodynamics code, ChaNGa (Menon et al., 2015), which incorporates standard
physics modules previously used in GASOLINE (Wadsley et al., 2004; Wadsley et
al., 2008; Wadsley et al., 2017). This includes a cosmic UV background (Haardt
& Madau, 2012) with self-shielding (Pontzen et al., 2008), star formation with
‘blastwave’ supernovae feedback (Stinson et al., 2006), low temperature metal
cooling (Guedes et al., 2011), and thermal and metal diffusion (Shen et al.,
2010; Governato et al., 2015). It also includes recent improvements, such as
an SPH force calculation that uses the geometric mean density (Ritchie &
Thomas, 2001; Menon et al., 2015; Governato et al., 2015), an updated
turbulent diffusion implementation (Wadsley et al., 2017), and a time-
dependent artificial viscosity and time-step adjustment system (Saitoh &
Makino, 2009; Wadsley et al., 2017).

###  2.2 Star Formation, Supernovae Feedback and Gas Cooling

Star formation occurs within dense (n$>0.2$ cm-3), cold (T$<10^{4}$ K) gas.
Each gas particle that meets these criteria is allowed to form star particles
on a characteristic timescale of $10^{6}$ years with the following
probability:

$p=\frac{m_{gas}}{m_{star}}\left(1-e^{-c_{\star}(10^{6}/t_{\mathrm{dyn}})}\right),$
(1)

where $m_{\mathrm{star}}=0.3m_{\mathrm{gas}}$, $c_{\star}=0.15$, and
$t_{\mathrm{dyn}}$ is the dynamical time of the gas particle. Energy from
supernovae couples thermally to nearby gas with an efficiency of 75%.
Supernova feedback uses the ‘blastwave’ implementation (Stinson et al., 2006),
where gas cooling is shutoff for a period of time to avoid numerical
overcooling. This implementation of feedback and star formation produces dwarf
galaxies that lie on the observed stellar mass-halo mass relations (Tremmel et
al., 2017) and includes a realistic population of ‘ultra-diffuse’ galaxies
(Tremmel et al., 2020; Wright et al., 2021; Van Nest et al., 2022).

An important limitation of Romulus is the lack of high temperature metal-line
cooling. This can affect the accretion history of gas onto massive galaxies
(van de Voort et al., 2011). This choice, discussed in more detail in Tremmel
et al. (2019), is motivated by a variety of previous simulations showing that
the inclusion of metal-line cooling in simulations without adding molecular
hydrogen physics and more detailed star formation prescriptions results in
unrealistic dwarf galaxies (Christensen et al., 2014). Despite being among the
highest resolution simulations of its class, even Romulus is unable to resolve
the multiphase ISM so it is not possible to include these more detailed
physical processes while maintaining realistic low mass galaxies. It has been
shown that RomulusC, our most massive halo using these input physics,
maintains a realistic intracluster medium (Tremmel et al., 2019). Because this
paper is focused on dwarf galaxies and on the formation of MBHs at high
redshift from very metal-poor gas (see below), this choice does not impact our
results.

###  2.3 Black Hole Accretion and Feedback

Massive black holes are allowed to grow by accreting nearby gas via a modified
Bondi-Hoyle formalism (Bondi & Hoyle, 1944) that accounts for angular momentum
support and includes a density-dependent boost factor (Booth & Schaye, 2009)
that is meant to account for unresolved, multiphase gas:

$\dot{M}_{\bullet}=\left(\frac{n}{n_{*}}\right)^{\beta}\begin{cases}\frac{(GM_{\bullet})^{2}\rho}{(v^{2}_{\mathrm{bulk}}+c^{2}_{s})^{3/2}}&\text{if
$v_{\mathrm{bulk}}>v_{\theta}$}\\\ \frac{(GM_{\bullet})^{2}\rho
c_{s}}{(v^{2}_{\theta}+c^{2}_{s})^{2}}&\text{if
$v_{\mathrm{bulk}}<v_{\theta}$},\\\ \end{cases}$ (2)

where $G$ is the gravitational constant; $n_{*}$ is the star formation
threshold (0.2 cm-3; see previous section); $v_{\mathrm{bulk}}$ is the bulk
velocity of the gas; $v_{\theta}$ is its rotational velocity; $c_{s}$ is its
sound speed; $\rho$ is its mass density; and $\beta$ is set to 2. The
radiative efficiency ($\epsilon_{r}$) of accreting MBHs is assumed to be 0.1
and all accretion is capped at the Eddington rate.

Growing MBHs produce feedback at a rate

$\displaystyle\dot{E}_{BH}=\epsilon_{r}\epsilon_{f}\dot{M}_{\bullet}c^{2},$
(3)

where $\epsilon_{f}$ is set to 0.02 and $c$ is the speed of light. This energy
is distributed instantaneously as thermal energy to the surrounding 32 nearest
gas particles. Feedback from MBHs has been shown to successfully regulate the
star formation in simulated massive galaxies in Romulus (Tremmel et al., 2017,
2019; Chadayammuri et al., 2021). However, too many massive galaxies remain
star forming and disk-dominated in Romulus at $z=0$ (Jung et al., 2022),
indicating that stronger modes of feedback are still needed for higher mass
galaxies. On the other hand, at the lower mass end, feedback from MBHs has
also been shown to influence the evolution of dwarf galaxies in Romulus
(Sharma et al., 2020, 2022a, 2022b) and has been shown to over-quench star
formation at low masses, indicating that feedback may be too efficient at low
masses. Because this paper focuses on the formation of black holes and less on
the detailed evolution of their host galaxies, this should not affect our
results. The majority of low mass galaxies in RomulusC are quenched by the
cluster environment and not by internal processes (Tremmel et al., 2019).

###  2.4 Black Hole Dynamics

Unlike most cosmological simulations that include MBHs, Romulus does not force
black holes to the centers of galaxies. Rather, they are allowed to move
realistically within their host galaxies. This is done through the
implementation of a sub-grid routine to account for unresolved dynamical
friction (Tremmel et al., 2015). This estimates and applies the necessary
force each MBH would feel due to dynamical friction by integrating the
Chandrasekhar formula (Chandrasekhar, 1943) out to the gravitational softening
length, $\epsilon_{g}$, during each black hole timestep.

$\boldsymbol{a}_{\mathrm{DF}}=-4G^{2}M_{\bullet}\rho(<v_{\bullet})\ln\Lambda\frac{\boldsymbol{v}_{\bullet}}{v_{\bullet}^{3}},$
(4)

Here $v_{\bullet}$ is the velocity of the black hole relative to nearby star
and dark matter particles, $\rho(<v_{\bullet})$ is the local density of star
and dark matter particles that are moving slower than the black hole relative
the background particles, and $\ln\Lambda$ is the Coulomb logarithm and equal
to $\ln(\epsilon_{g}/r_{90})$. Here $r_{90}$ is the $90^{\circ}$ deflection
radius based on the black hole’s mass and local relative velocity. As shown in
Tremmel et al. (2015), this method is able to produce realistically decaying
orbits. The result of this model is that MBHs take non-negligible time to
reach the halo center after galaxy mergers (Tremmel et al., 2018a), and
sometimes their orbits fail to decay and they end up as ‘wandering’, off-
center MBHs (Tremmel et al., 2018b; Bellovary et al., 2021; Ricarte et al.,
2021a, b). Two MBHs are allowed to merge when they are within 700 pc
($2\times\epsilon_{g}$) and mutually bound to one another.

The technique used here to model dynamical friction is similar to that
employed in a growing number of other simulations, including Magneticum
(Hirschmann et al., 2014), Obelisk (Trebitsch et al., 2021; Pfister et al.,
2019) and Astrid (Chen et al., 2022; Ni et al., 2022). The high resolution of
Romulus allows for the dynamical evolution of MBHs to be accurately tracked
down to sub-kpc scales with dynamical friction applied only locally, sampling
densities of stars and dark matter near each MBH (distances $<350$ pc). Other
simulations like Horizon-AGN (Dubois et al., 2016) have models that account
only for dynamical friction from gas (Ostriker, 1999) while Romulus only
accounts for dynamical friction from stars and dark matter (i.e. collision-
less particles). While it is possible that gas may dominate the density of a
galaxy at times, it is difficult to fully account for the effects of torques
due to unresolved structure and turbulence within the ISM which can be at
least as important (Roškar et al., 2015; Bortolas et al., 2022; Lescaudron et
al., 2022).

Figure 1: Distribution of Black Hole Formation Times. The cumulative
distribution of formation times for MBHs in Romulus25 and RomulusC. Despite
sampling different environments, both simulations seed MBHs at very similar
times.

###  2.5 Black Hole Seeding

A critical aspect of the Romulus simulations to this work is the prescription
for MBH seeding. Like most modern, large-scale cosmological simulations, the
actual formation sites for MBH seeds are unresolved in Romulus. To approximate
the regions most likely to create an early, massive black hole, we form MBH
seeds in gas that would otherwise be forming stars (see above) but with the
following additional criteria:

  * •

Metallicity less than $3\times 10^{-4}\ Z_{\odot}$.

  * •

Particle density is greater than $3\ cm^{-3}=15n_{\star}$.

  * •

Temperature is between $9500\ \mathrm{K}$ and $10000\ \mathrm{K}$.

This prescription will naturally pick out regions in the very early Universe
that have yet to be polluted with metals from star formation, yet are still
able to collapse to high densities (15 times greater than the threshold for
star formation) on timescales shorter than that for star formation (i.e. the
gas is able to reach densities well beyond the star formation threshold before
forming stars) and without cooling much below $10^{4}$ K (densities are
growing faster than the gas can effectively cool). While these criteria are
most directly reminiscent of the direct collapse formation scenario (gas
collapsing into a single massive object; Lodato & Natarajan, 2006, 2007; Regan
et al., 2017; Wise et al., 2019; Regan et al., 2020b), or the collapse of a
dense star cluster (runaway stellar collisions produce a single massive
object; Devecchi & Volonteri, 2009; Davies et al., 2011), they are also
consistent with regions where a low mass black hole seed (e.g. from a
population III star) may be able to grow rapidly in a very short amount of
time thanks to the presence of of dense, rapidly collapsing gas (e.g.
Volonteri & Rees, 2005; Alexander & Natarajan, 2014; Inayoshi et al., 2016;
Sassano et al., 2023).

Importantly, these adopted criteria result in an epoch of MBH seeding that is
roughly 90% complete by $z=5$ and makes no assumptions as to which halos
should (or should not) host a MBH (Tremmel et al., 2017). This is in stark
contrast to the methods more typically utilized in cosmological simulations,
where MBHs are seeded based on a halo mass threshold (typically a few $\times
10^{10}$ M⊙), which preferentially prevent MBHs from populating low mass
galaxies early on, and will result in a prolonged period of MBH seeding. By
predicting the existence of MBHs based only on local gas properties, Romulus
can therefore predict the occupation of MBHs in galaxies in a way that even
simulations of similar resolution often cannot. Some cosmological simulations
have started to employ a similar MBH formation prescription. For example, the
New Horizon simulations have shown the importance of both resolving low mass
galaxies and allowing them to be seeded with MBHs in studying MBH mergers and
potential gravitational wave sources (Volonteri et al., 2020). However, it is
important to note that this method would not capture formation processes that
occur at later times in more metal-enriched gas (e.g. Regan et al., 2020a;
Natarajan, 2021; Mayer et al., 2023).

In Figure 1 we compare the distribution of black hole formation times in the
field ( Romulus25) with that of the galaxy cluster (RomulusC), confirming that
they are virtually indistinguishable in our simulations. This means that the
initial seeding of MBHs at high redshift, which only depends on local gas
properties, is not impacted directly by the environment. This will be an
important fact to keep in mind as we delve deeper into the differences
inferred in the MBH population between the two environments.

###  2.6 Halo Selection and Galaxy Properties

Halo finding is performed with the Amiga Halo Finder (Knollmann & Knebe,
2009), which uses a spherical top-hat collapse technique to define the virial
radius (Rvir) and mass (Mvir) of each halo and sub-halo. It also assigns all
of the baryonic content belonging to each dark matter halo. In our analysis we
use R200 which is the radius enclosing a mean density 200 times the critical
density of the Universe at a given redshift. We only include halos with
M${}_{vir}>3\times 10^{9}$ M⊙, such that each halo included in our analysis
has at least $\sim 10,000$ particles. An exception is made for our analysis of
RomulusC where we include halos down to $3\times 10^{8}$ M⊙ to account for the
fact that many halos have experienced tidal stripping of the outer regions of
their dark matter halos in the dense cluster environment. As discussed in
Tremmel et al. (2020), while the total halo mass is affected by the cluster
environment, dwarf galaxies often keep most of their stellar mass intact. This
choice to lower the halo mass threshold is to avoid discounting galaxies that
would otherwise have been previously considered well-resolved prior to falling
into the cluster. We have confirmed that our results are not sensitive to this
choice.

The position of each galaxy/halo is calculated using the shrinking spheres
approach (Power et al., 2003), which reliably extracts the centers of the
central galaxies within each halo. As RomulusC is a zoom-in simulation, it is
surrounded by a region of low resolution dark matter. Galaxies on the
outskirts may be contaminated by high mass (low resolution) dark matter
particles. We avoid including such galaxies in our analysis by requiring each
galaxy to have less than $5\%$ of its dark matter particles be contaminated by
low resolution elements. This cut removes only 3% of galaxies within 2 Mpc
($\sim 2$R200) of the cluster center, the vast majority of which are beyond
the virial radius.

When quoting the stellar mass of each galaxy, we estimate what the observed
stellar mass would be from typical techniques. Munshi et al. (2013) found that
typical observational techniques result in a systematic underestimate of
galaxy stellar mass. Based on those results, we apply a correction factor of
0.6 to the ‘raw’ stellar masses (the summed mass of all star particles
associated with a given halo), a conservative estimate given the results of
Munshi et al. (2013). Leja et al. (2019) also find a similar discrepancy when
comparing advanced spectroscopic techniques to estimating stellar masses with
common photometric techniques. This choice affects only the stellar masses we
present, but since all galaxies are treated equally in this analysis this has
no effect on the substance of our results comparing different simulation
regions.

Throughout this work we also classify some galaxies as ‘isolated’. We base our
definition of isolated on the results of Geha et al. (2012) that show the
onset of environmental effects for dwarf galaxies (M${}_{\star}<10^{10}$ M⊙)
occur when they are approximately 1.5 Mpc away from a massive galaxy. The same
environmental effects are seen on the quenched fraction of galaxies in Romulus
(Sharma et al., 2022b). Isolated galaxies in Romulus25 are therefore defined
so that they are not within R200 of any more massive halo. Additionally, for
low mass galaxies with M${}_{\star}<10^{10}$ M⊙, they must be further than 1.5
Mpc from any galaxy with stellar mass greater than $2.5\times 10^{10}$ M⊙.

##  3 The Predicted occupation fraction of MBHs in Romulus

The top panel of Figure 2 shows the fraction of galaxies hosting at least one
MBH in both the galaxy cluster (RomulusC) and in the field (Romulus25). For
Romulus25, only galaxies within halos of mass $>3\times 10^{9}$ M⊙ are
included. For RomulusC, only galaxies in halos of mass $>3\times 10^{8}$ M⊙
and within $2R_{200}$ of cluster center are included (see previous section for
more justification of this choice, though it does not affect our overall
results). The results are compared with multiple observationally derived
estimates for the underlying occupation fraction (Greene, 2012; Miller et al.,
2015; Askar et al., 2022) and the observed occupation fraction of nuclear star
clusters (NSCs) in dwarf galaxies within the Virgo and Fornax clusters
(Sánchez-Janssen et al., 2019; Muñoz et al., 2015), as well as the Local Group
(Hoyer et al., 2021). Note that observational estimates for the underlying MBH
occupation fraction are highly uncertain and rely on inference from active
galactic nuclei. In the context of this work, these are useful benchmarks to
compare our results while more detailed comparisons are not very useful. The
error bars in occupation fraction in this and all following figures represent
95% binomial confidence intervals (Cameron, 2011).

The occupation fractions in the field derived from Romulus25 are consistent
with estimates derived from observations and match quite well with observed
NSC occupation fractions in the local volume (Hoyer et al., 2021). We find a
significantly elevated occupation fraction among cluster dwarf galaxies
compared to field dwarfs with M${}_{\star}<10^{9}$ M⊙. Interestingly, the
occupation fraction in these low mass cluster galaxies matches remarkably well
with the occupation fraction of nuclear star clusters in Fornax and Virgo,
which also show enhanced occupation compared to lower density environments
(Sánchez-Janssen et al., 2019; Muñoz et al., 2015).

The bottom panel of Figure 2 shows the occupation fraction of central (within
1 kpc of the halo center) black holes that would be observed as active
galactic nuclei at $z=0.05$ as a function of stellar mass. We limit the
spatial offset to mimic the fact that observers often search for X-ray sources
associated with galactic centers specifically. We find that cluster dwarf
galaxies are less likely to host X-ray luminous AGN compared to the field. As
discussed in Tremmel et al. (2019) this is because the majority of cluster
dwarf galaxies have had their gas supply removed due to ram pressure
stripping, resulting in very low MBH growth. We note that it is possible that
ram pressure can also cause gas to compress to the center and, in some cases,
momentarily activate black hole growth (Poggianti et al., 2017). We see this
effect in RomulusC (Ricarte et al., 2020) but this process is transient and
only common among the more massive in-falling galaxy population. Low mass
galaxies quickly lose their gas and both star formation and black hole
accretion are shut off.

Figure 2: The excess of mostly quiescent massive black holes in cluster dwarf
galaxies. Top: The fraction of galaxies hosting at least one MBH as a function
of stellar mass for galaxies in Romulus25 (blue) and RomulusC (orange). Error
bars represent 95% binomial confidence intervals (Cameron, 2011). The grey
region and open diamonds are both estimates derived from observations in
Miller et al. (2015) and Greene (2012) respectively. The black dashed line is
the predicted occupation of MBHs above $10^{5}$ M⊙ from Askar et al. (2022)
assuming a 30 Myr mass doubling time for MBHs. The three hatched regions are
occupation fractions of nuclear star clusters (NSCs) observed in dwarf
galaxies in Virgo (green; Sánchez-Janssen et al., 2019), Fornax (magenta;
Muñoz et al., 2015), and the local volume (cyan; Hoyer et al., 2021). For the
local volume data, the hatched region represents 95% binomial confidence
intervals, which were calculated using the total number of observed galaxies
in each bin. While the field population in Romulus25 has occupation fractions
roughly consistent with observationally derived estimates, RomulusC has an
enhanced occupation fraction at low masses (M${}_{\star}<10^{9}$M⊙). The
occupation fraction in the RomulusC simulation is remarkably consistent with
that in the observed NSCs within galaxy clusters, while the occupation
fraction in field galaxies (Romulus25) is very consistent with local volume
NSCs. Bottom: The fraction of galaxies hosting a central ($D<1$ kpc), luminous
(L${}_{x}>2\times 10^{38}\ \mathrm{erg}\;\mathrm{s}^{-1}$) MBH as a function
of stellar mass. Also shown are two observational data sets examining AGN in
cluster and early-type galaxies with a similarly low luminosity threshold
(Gallo et al., 2010; Miller et al., 2015). Also plotted in the dashed black
line is the active fraction estimated from Pacucci et al. (2021). The cluster
dwarf population (orange), which are preferentially quenched galaxies, matches
well with the observations and are noticeably lower than the occupation of
luminous MBHs in field dwarfs (blue). Despite having a higher underlying MBH
occupation fraction, cluster dwarfs are less likely to host low luminosity AGN
compared to the field by $z=0$.

To estimate the X-ray luminosity of growing MBHs in the simulation, we first
estimated the bolometric luminosity. To do this we followed previous work
(Churazov et al., 2005; Habouzit et al., 2022; Sharma et al., 2022a) and
implemented a two-mode model that accounts for radiatively inefficient
accretion flows during low Eddington ratio ($f_{\mathrm{Edd}}$) accretion:

$L_{\rm bol}=\begin{cases}\epsilon_{r}\dot{M}_{\rm BH}c^{2},\quad&f_{\rm
Edd}\geq 0.1\\\ 10f_{\rm Edd}\epsilon_{r}\dot{M}_{\rm BH}c^{2},\quad&f_{\rm
Edd}<0.1.\end{cases}$ (5)

This assumes the same constant radiative efficiency ($\epsilon_{r}=0.1$) used
in the simulation to calculate Eddington ratio, applying a scale factor to
estimate an ‘effective’ radiative efficiency for low, radiatively inefficient
accretion rates ($f_{Edd}<0.1$). This is purely a post-processing calculation
to estimate the total escaping luminosity from each black hole. In the
simulation a constant $\epsilon_{r}=0.1$ was used to calculate both Eddington
ratio and feedback (see equation 3). While this estimate of bolometric
luminosity doesn’t directly contradict any of the physics implemented in the
simulation, it does change how feedback relates to (observed) luminosity, as
the original feedback model rationale was that some constant fraction of
emitted light couples to nearby gas as thermal energy, which would no longer
be true when using equation 5.

Figure 3: Comparing MBH Occupation Fraction Across Cosmic Time. The occupation
fraction of MBHs as a function of galaxy stellar mass at different redshifts.
Blue points show the results for galaxies in the field (Romulus25) and the
orange points show galaxies within 2 Mpc of the (proto-)cluster center at each
redshift. At $z=5$, when the vast majority ($>95\%$) of MBHs have formed in
both simulations, the occupation fractions look very similar. The differences
in the two distributions becomes more dramatic with cosmic time, particularly
after $z\sim 3$. The large error bars are due to small number statistics at
the highest mass bins at earlier times.

After calculating the bolometric luminosity, we apply the bolometric
correction from Shen et al. (2020) to estimate the 0.5-10 keV X-ray luminosity
of each MBH. The rationale and application of this approach is discussed
further in Sharma et al. (2022a) where it is shown that calculating luminosity
using equation 5 brings the simulated dwarf AGN fraction in Romulus much
closer to observed values. Note that here we are not making any additional
cuts based on the estimated X-ray luminosity of stars and gas in the galaxy,
as is done in Sharma et al. (2022a). Including this would likely bring down
the active fraction in the field and have less of an effect for the cluster
galaxies, bringing the two potentially closer together. Regardless, our
results would remain the same in that the enhanced occupation fraction of MBHs
in the cluster is not expressed in the active fraction, which only includes
MBHs that are likely to be detectable with observations.

MBHs in the Romulus simulations are allowed to be off-center and indeed often
are (Tremmel et al., 2018a, b; Ricarte et al., 2021a, b). Such off-center
black holes may even be preferentially more common in dwarf galaxies, which
host lower mass black holes that experience more inefficient dynamical
friction (Bellovary et al., 2021). Observations searching for MBHs in low mass
galaxies often focus on their centers and this would artificially lower the
observed active fraction. However, we confirm for both sides of Figure 2 that
the choice to include/exclude off-center ($D>1$ kpc) MBHs makes very little
difference in our predicted occupation fractions. Placing more strict criteria
for hosting a MBH will decrease the overall fraction of galaxies but the
effect is minor and of similar magnitude across environments. For luminous
MBHs, it is difficult for non-central MBHs to accrete enough material to
become luminous so the inclusion of off-center MBHs also has little effect
here. Importantly, the decision to include/exclude off-center MBHs does not
affect our main prediction: an enhanced occupation fraction of MBHs in low
mass cluster galaxies.

###  3.1 Evolution with Redshift

In Figure 3 we show the MBH occupation fraction as a function of stellar mass
at six snapshots at different redshifts. We only examine this out to $z=5$
because before this time many MBHs are still actively being seeded while, by
$z=5$, the vast majority ($\sim 90\%$) of MBHs have formed (see Figure 1).
From Figure 3 we can see that the occupation fraction in the different
environments begins to look very similar beyond $z\sim 3$. At $z=5$ the
occupation of MBHs as a function of stellar mass is nearly identical between
field and cluster environments. Therefore, both the formation times (Figure 1)
and host halos at high redshift (Figure 3; lower right panel) are similar
between RomulusC and Romulus25, indicating that the presence and location of
dense, metal-free gas is very similar between environments at early epochs.

In Figure 4 we compare the occupation fraction within each environment at
different redshifts. In the field (right), the occupation fraction decreases
steadily with redshift, particularly at lower masses. For the cluster (left)
the occupation fraction ceases to decrease significantly past $z\sim 3$. It is
this difference in evolution that results in the difference seen at $z=0.05$.
The occupation fraction in the cluster at $z=0.05$ is more similar to the
field at $z=3$, when the evolution stopped for cluster galaxies. In contrast,
the decline in MBH occupation fraction continues for field dwarfs throughout
cosmic time. As we will discuss further in Section 4, this decline seen in the
field is driven by late-forming dwarf galaxies which do not exist in the
cluster environment.

###  3.2 Dependence on Cluster-centric Distance and Halo Mass

We can explore in more detail the extent to which the occupation fraction is
dependent on environment. In the left-hand plot in Figure 5 we find evidence
that the occupation fraction evolves with cluster-centric distance. As one
looks more toward the cluster outskirts ($0.75<$D$<2$ R200) the MBH occupation
fraction in galaxies does decline, although it remains systematically higher
than that in isolated galaxies at stellar masses below $10^{9}$ M⊙. This
implies that while we might expect a gradual evolution from isolated to
cluster galaxies, the influence of the environment persists even in the
outskirts of clusters.

We can use Romulus25 to see if there are enhancements in the occupation
fraction in less dense environments, such as low mass groups. The right-hand
plot in Figure 5 plots the occupation fraction of isolated galaxies in blue,
cluster galaxies in orange, and galaxies within 2R200 of a more massive halo
with virial mass between $10^{12}$ \- $10^{13.3}$ M⊙ in red (i.e. the most
massive halos in the 25 Mpc volume). We find no evidence that dwarf galaxies
near these larger galaxies have any systematic enhancements to their MBH
occupation. We compare again to the observed nuclear star cluster occupation
in dwarfs from Hoyer et al. (2021), as well as Carlsten et al. (2022), and
find that the MBH occupation in dwarf galaxies associated with lower mass
halos match well with observed NSCs in the local volume.

We confirm that these results are not sensitive to the specific host halo mass
range used, though we are limited by the small volume of Romulus25. These
results imply that the enhancement in occupation fraction we see in the
Romulus simulations is isolated to very massive halos (massive groups and
above; M${}_{vir}>10^{13.3}$ M⊙). Lower mass halos assemble their mass
earlier, meaning that low mass galaxies in-fall at higher redshift when the
halo has a shorter characteristic dynamical time. These halos are more likely
to become disrupted before $z=0$ and contribute to the population of wandering
black holes in massive halos (Tremmel et al., 2018a, b; Ricarte et al., 2021a,
b). It may also be that these lower density environments have a MBH occupation
fraction enhancement only at galaxy masses that are currently unresolved in
Romulus.

##  4 Exploring the Origin of the Environmental Dependence of Dwarf Galaxy MBH
Occupation

Figure 4: The Evolution of the MBH Occupation Fraction Stalls in the Cluster
Environment. Similar to Figure 3 here we show the occupation fraction for
galaxies as a function of stellar mass at six different redshifts for cluster
galaxies (left) and field galaxies (right). While the occupation fraction in
field dwarf galaxies decreases steadily through time, the evolution is slower
in the cluster environment and plateaus at $z=2-3$.

The evolution of galaxies in dense environments like galaxy clusters is
different compared with galaxies in the field. For cluster member galaxies,
eventually, their ability to accrete new mass or form new stars will be shut
off by the cluster environment (Mistani et al., 2016). However, even before
galaxies become bound to the cluster, or in-fall to within R200, they are
growing within an over-dense region. This means that the galaxies that exist
within a cluster environment at a given stellar (or halo) mass will have very
different evolutionary histories compared with similar mass galaxies in the
field (Gao et al., 2005; Gao & White, 2007; Croton et al., 2007; Boylan-
Kolchin et al., 2009). Cluster dwarf galaxies are likely to have earlier
formation times because they must assemble their material before the cluster
environment shuts down their growth. It may also be true that cluster dwarfs
would have grown to much larger masses if their mass assembly was allowed to
continue unabated by their environment (Mistani et al., 2016). Finally,
cluster galaxies may also have their mass removed through tidal stripping as
they interact with the cluster potential. In some cases, this stripping can
unbind the majority of their stars, leaving behind only the dense, central
stellar component (such as a NSC), which is then detected as an ultra-compact
dwarf (Drinkwater et al., 2003; Bekki et al., 2003; Seth et al., 2014; Voggel
et al., 2016). In the following sections we evaluate the ability of each of
these evolutionary scenarios in explaining the unique MBH occupation fraction
in the simulated cluster (RomulusC) environment compared to the field
(Romulus25).

###  4.1 Cluster dwarfs as stripped remnants of more massive galaxies

The enhanced dwarf galaxy MBH occupation fraction in RomulusC could be
explained if low mass host galaxies in the cluster were once much more
massive. If a more massive galaxy were to become tidally stripped, it could
lose significant stellar mass while retaining its central MBH. If this is true
for a large number of low mass galaxies in RomulusC, then it would make sense
that their MBH occupation fractions would resemble that of more massive
galaxies (i.e. their original mass prior to tidal stripping). However, as
discussed in Tremmel et al. (2020), while significant dark matter mass is lost
due to interacting with the cluster potential, the majority of galaxies have
not lost much stellar mass. While much of the dark matter stripping occurs on
larger scales, stellar mass is confined within the galaxies themselves. Tidal
stripping of this much more compact component requires closer, more intense
interactions with the center cluster potential that are likely to result in
the complete disruption of the galaxy, as far as the simulation and halo
finder are concerned. It is important to note that some of this disruption
could be artificial due to the limited resolution of the simulation (van den
Bosch & Ogiya, 2018). It is also important to note that Romulus cannot resolve
extremely compact structures within galaxies that would be able to best
survive significant tidal interactions, such as NSCs.

The analysis from Tremmel et al. (2020) was done by tracing halos back in time
to compare their $z=0$ stellar mass with their maximum stellar mass. This
could induce a bias where the halos that pass closest to cluster center are
more likely to have time-steps where they are missed by the halo finder, a
well known issue (e.g. Knebe et al., 2011; Onions et al., 2012; Joshi et al.,
2016). This would cause us to be unable to fully trace their evolution through
time and these potentially heavily stripped galaxies would be ignored.
Focusing on MBH hosts, we can instead trace the MBHs themselves through time
without relying on halo finding and compare the final host stellar mass with
the maximum host stellar mass, excluding any intervening steps where they are
temporarily taken to be hosted by the main cluster halo (i.e. times where the
halo finder fails to extract their host sub-halo).

Figure 5: The Occupation Fraction As A Function of Environment. Left: The
occupation fraction of MBHs in galaxies in RomulusC at different cluster-
centric distances. The orange points represent galaxies within 0.75 R200 of
cluster center while the red points are for galaxies in the cluster outskirts
($0.75<\mathrm{D}/\mathrm{R}_{200}<2$). The blue points are for isolated
galaxies in Romulus25. Although a subtle effect, there is evidence that
galaxies in cluster outskirts do have less MBHs, though the low mass galaxies
below 109 M⊙ are still enhanced compared to isolated galaxies. Right: The
occupation fraction of galaxies in Romulus25 in different environments. The
red points show the occupation fraction of galaxies that are within 2R200 of a
halo of mass between $10^{12}$ and $10^{13.3}$ M⊙. There is little difference
between the occupation fraction at these intermediate masses and that for
isolated dwarf galaxies (blue points). We also show observational results for
local volume dwarfs from Carlsten et al. (2022) and Hoyer et al. (2021) as
hatched regions.

Focusing only on MBHs hosted in cluster dwarf galaxies
(M${}_{\star}<10^{10}$M⊙) at $z=0.05$, we find that the median MBH host galaxy
has only experienced a net loss of $\sim 20\%$ of its stellar mass as it
interacts with the cluster environment. Only one fifth of the MBH hosts have
seen their stellar mass decrease by more than a factor of $\sim 3$. Looking at
figure 2, a typical stellar mass loss of a factor of $\sim 3-5$ is needed to
bring the occupation fractions in line with the field. In RomulusC such
extreme mass loss is too rare to fully explain the overabundance of MBHs in
low mass galaxies.

###  4.2 Cluster Dwarfs as Failed Massive Galaxies

Galaxies in cluster environments will eventually stop accreting new material,
as both dark matter and gas will flow onto the primary halo instead. The lack
of replenishing gas supply combined with ram pressure removing the ISM will
eventually slow down or completely quench new star formation in the galaxy.
However, it is possible that, were these galaxies allowed to continue to grow
unimpeded, they would be more massive at $z=0$. If the progenitor galaxies to
our cluster dwarfs are more similar to progenitors of massive field galaxies
than they are to those of field dwarfs, this could explain the difference in
MBH occupation fraction. In other words, it may be that cluster dwarfs
actually represent progenitors to more massive field galaxies that were
instead ‘frozen’ in their growth by their environment. In this scenario it is
not required the galaxies lose stellar mass, just that they fail to reach the
same masses as their field counterparts.

In order to explore this scenario, we trace our cluster galaxies back in time,
finding the redshift, stellar mass, halo mass, and concentration111We define
concentration here as vmax/v200, where vmax is the maximum circular velocity
and v200 is the circular velocity at R200. at the time they reach maximum
virial mass (tmax). We exclude all galaxies that fail to trace backward to a
time prior to falling into the cluster. These progenitors are then matched to
galaxies at that same redshift in Romulus25 which are major progenitors to
isolated $z=0.05$ galaxies, requiring the stellar and halo masses be within
0.2 dex and the difference in concentration be less than 0.2. For each cluster
galaxy, we require they match with at least 4 field galaxies with this
criteria. We then recalculate the occupation fractions using for each cluster
galaxy the median $z=0.05$ stellar mass among the matched isolated galaxies.
The idea is that this should approximate the stellar mass they would have
attained were they allowed to continue to grow.

Figure 6 shows the results of this analysis, plotting the occupation fraction
using both the original $z=0$ stellar masses (open orange points) and the
stellar masses calculated by matching the progenitors to isolated galaxies
(solid orange points). We are only able to do this analysis with galaxies with
original stellar masses above $10^{8}$ M⊙. At lower masses too many galaxies
become excluded because we fail to trace them back before in-fall into the
cluster due to the halo finder failing to identify them at some point. We
confirm that these results are not sensitive to our specific matching
criteria, specifically which combination of halo mass, stellar mass, and
concentration were used, as well as the number of matches required for the
analysis.

If the majority of low mass MBH host galaxies are more like progenitors to
massive, isolated galaxies, the ‘corrected’ occupation fraction technique
would see many MBHs hosts move to larger (corrected) stellar masses,
decreasing the occupation fraction at low mass and bringing the results more
in-line with the field. While we see in Figure 6 that this matching technique
does result in lower occupation fractions at low masses, the results are still
systematically higher than the field. This implies that while some MBH host
galaxies may be considered ‘failed’ massive galaxies (i.e. they could have
grown larger had their mass accretion not been shut down by the cluster
environment) this scenario fails to fully explain the discrepancy in
occupation fraction between environments. Of course, because some galaxies had
to be excluded due to bad tracking, the error bars are larger and the
difference between the corrected and original occupation fractions for
individual mass bins are often only marginally significant.

###  4.3 An Overabundance of Early-forming Dwarfs in Clusters

As discussed in Sharma et al. (2020), dwarf galaxies in Romulus25 hosting
black holes, especially those that are more massive, tend to have earlier
formation times for both their stars and their overall halo mass. While
feedback from black holes could influence the assembly history of stars,
potentially quenching star formation even in low mass galaxies (Sharma et al.,
2022b; Koudmani et al., 2021), the fact that this trend is also seen in dark
matter halo assembly indicates that it is likely something more fundamental.
In cluster environments, the assembly of galaxies and the dark matter halos in
which they reside is stopped by the cluster environment, such that a galaxy of
a given mass in the cluster must have assembled that mass prior to in-fall.
This is an expected result of assembly bias in the formation of dark matter
halos (e.g. Gao et al., 2005; Gao & White, 2007; Croton et al., 2007; Boylan-
Kolchin et al., 2009) and is seen in other cosmological hydrodynamic
simulations (e.g. Mistani et al., 2016; Chaves-Montero et al., 2016).

Figure 7 plots the occupation fraction as a function of formation time for the
same two stellar mass bins of dwarf galaxies. The top panels show the halo
formation time (time to accumulate 50% of the maximum halo mass) and the
bottom panels show the stellar formation time (time to form 50% of the maximum
stellar mass). The orange and blue bands show the average occupation fraction
for successfully tracked galaxies in each mass bin. Note that the average
values are calculated only for the galaxies that have been successfully traced
back through time and included in the calculation of the individual data
points shown here.

In both the field and cluster environments, dwarf galaxies with earlier
formation times are more likely to host MBHs. This is true when examining
either stars or halo mass. The MBH occupation fraction for cluster dwarfs is
similar to field dwarfs when controlling for formation time, though this
connection is better illustrated by halo mass when considering higher mass
dwarfs. In the field, dwarfs of a given mass are allowed to form throughout
cosmic time, but those that accumulate their mass later are less likely to
host MBHs. This lack of late-forming dwarf galaxies in the cluster is what
causes the evolution of the occupation fraction to stop evolving after $z\sim
3$ (see Figures 3 and 4). In the field, ‘new’ dwarf galaxies grow at later
times, filing those lower mass bins with galaxies that lack MBHs. While this
is a function of the specific seeding criteria we use, these results show that
regions of very dense, rapidly collapsing, pristine gas are more likely to
exist in the progenitors to early-forming dwarf galaxies compared to late-
forming dwarfs. In late-forming dwarfs, such early phases of collapse occur
too slowly, allowing for the formation of stars and metal enrichment before
the required high densities (far beyond the threshold for star formation in
the simulation) are reached (if they ever are). In this scenario, the role of
environment is more to stop the formation of late-forming dwarf galaxies,
rather than influence on the formation sites of MBHs. As can be seen in
Figures 1, 3, and 4, the time and host halos of MBH seeding is very similar
between the two environments.

Figure 6: The Occupation Fraction for Galaxies Matched with Isolated Galaxies.
Here we test the extent to which the difference in occupation fraction can be
explained by cluster dwarf galaxies being ‘failed’ larger galaxies, i.e.
galaxies that, were they not in a cluster environment, would have grown to be
more massive. In blue we plot the occupation fraction for isolated galaxies in
Romulus25 at $z=0$. The solid orange points show the occupation fraction as a
function of galaxy mass after each cluster galaxy, where its mass has been
corrected to estimate what it might have attained had the galaxy not fallen
into a cluster. We do this by matching each cluster galaxy with a $z=0$
isolated galaxy based on its stellar and halo masses at $t_{50,halo}$ (the
time when the host dark matter halo reaches 50% of its maximum mass; see text
for details). The open, light orange points are these same galaxies but with
their original final masses. Note that we do not include the lowest mass
galaxies in this analysis because too many of them fail to be traced back in
time successfully (see text for details). While this matching process does
alleviate some of the differences, the ‘corrected’ occupation fraction for
cluster dwarfs remains systematically higher than isolated dwarf galaxies.
This implies that stunted mass growth is only a part of the explanation for
the different occupation fractions.

Once again, we face the problem discussed in the previous section whereby low
mass (M${}_{\star}<10^{8}$ M⊙) galaxies with MBHs are more likely to be
excluded because they fail to be tracked back in time successfully. These
galaxies form earlier and fall into the cluster sooner so they are more likely
to have passed closer to the cluster center and also be missed by the halo
finder. However, there is still a significant difference in the mean
occupation fraction in each mass bin relative to the field (comparing the
orange and blue bands).

Controlling for formation time results in a better match between the two
populations of galaxies when looking at halo mass, rather than stellar mass.
This makes sense, as many additional factors may play a role in the star
formation history of a galaxy, including the presence of feedback from a MBH
(Sharma et al., 2020, 2022b). The more massive dwarfs in clusters still appear
biased high relative to isolated galaxies with similar formation times. This
may indicate that a combination of effects are needed to fully explain this
enhanced MBH occupation population, i.e. some dwarfs could have assembled into
more massive galaxies were they isolated (see previous section) combined with
a lack of late-forming dwarf galaxies in clusters.

Figure 7: Early Forming Dwarfs are More Likely to Host MBHs. The occupation
fraction of MBHs as a function of $t_{50,halo}$ (top) and $t_{50,stars}$
(bottom) for dwarf galaxies in two stellar mass bins. The bands represent the
total occupation fraction in each mass bin (only for galaxies that can be
traced adequately far back in time; see text for details). The relationship
between formation time and MBH occupation fraction is similar for cluster
dwarfs (orange) and isolated dwarfs (blue), though galaxies in the higher mass
bin are still biased slightly high when controlling for either formation
timescale. Dwarf galaxies that form earlier are more likely to host a MBH. The
important difference between the two simulations, therefore, is that dwarf
galaxies within cluster environments form earlier than those that are more
isolated.

##  5 Discussion

Observational constraints on the underlying MBH occupation fraction in low
mass galaxies are uncertain because it is difficult to detect the low mass,
low luminosity black holes. Indeed while the Romulus simulation has been shown
to reproduce observed samples of dwarf galaxy AGN fractions and luminosities
(with specific assumptions made to convert between black hole accretion rate
and X-ray luminosity) the simulation predicts a large population of MBHs that
would go undetected by even the most sensitive modern X-ray surveys (Sharma et
al., 2022a). Still, as observations improve, evidence increasingly points to a
significant number of MBHs in low mass galaxies (Nguyen et al., 2018;
Baldassare et al., 2020; Burke et al., 2022). Much work is still needed from
the observational side, but upcoming time domain surveys from the Vera Rubin
Telescope may offer hope for dramatically increasing the completeness of the
observed MBH popuation through their intrinsic variability (Baldassare et al.,
2018; Burke et al., 2022) as well as tidal disruption events (Bricman &
Gomboc, 2020). JWST may also be a powerful tool for detecting AGN with low
X-ray luminosities (Cann et al., 2021). As observations continue to get better
at detecting MBHs in low mass galaxies, predictions like the ones made in this
work will be crucial in understanding and contextualising them with respect to
MBH formation models.

The challenge for simulations is, as always, resolution which comes at the
cost of the size and statistical sample of the data. Large-scale simulations
like Romulus reach a middle ground by being large enough to have many galaxies
while also capable of resolving dwarf galaxies. Still, smaller volumes means a
lack of environment diversity with only a handful of groups and a single low-
mass galaxy cluster. While newer, large-volume simulations are becoming better
in terms of both resolution and the black hole physics they implement (e.g.
Dubois et al., 2021; Trebitsch et al., 2021; Ni et al., 2022) it remains a
challenging balance. Even at the resolution of Romulus and these other state-
of-the-art simulations, the ISM remains largely unresolved, requiring
relatively simple prescriptions for black hole formation (see below for
further discussion). Zoom-in simulations are another viable path forward,
allowing for more detailed formation prescriptions (e.g. Dunn et al., 2018)
and more detailed analysis of both black hole dynamics and the internal
structure of dwarf galaxies (Bellovary et al., 2019, 2021). Very high
resolution simulations targeting the high redshift Universe have also proven
useful tools for examining the physics of MBH formation (Wise et al., 2019;
Regan et al., 2017, 2020b; Regan, 2023).

The seeding algorithm for MBHs implemented in Romulus is more predictive than
many previous large-scale cosmological simulations because it seeds MBHs based
on local gas properties (density, temperature, metalicity) without making any
a priori assumptions on which halos/galaxies should host a MBH. As discussed
below, our model is still simplistic because of limited resolution, but the
simulations still have significant predictive power. In particular, our
results demonstrate that the gas properties of galaxies at $z>5$ is connected
to their formation history and, therefore, so may be their likelihood of
hosting a MBH. Higher resolution simulations will be able to further test this
prediction, as will future observations. The environmental dependence of MBH
occupation that we predict here should be considered a potential way to
differentiate between MBH formation mechanisms. For example, an observed
environmental dependence of MBH occupation would support the theory that the
primary formation channel of MBHs occurs in the early ($z>5$) Universe from
metal poor gas. However, a lack of observed environmental dependence, based on
our results here, would indicate that other formation channels which have MBHs
grow at later times and from more metal polluted gas (e.g. Regan et al.,
2020a; Natarajan, 2021; Mayer et al., 2023) likely dominate.

###  5.1 The Effect of MBH Formation Model

The most important caveat to these results is that they will naturally rely on
our choice of MBH formation criteria. The main concern is whether our choices
directly influence our results, which is not the case here. The criteria in
Romulus are common sense requirements given any of the leading MBH formation
models and the requirement that each MBH should be able to grow to large
masses in a relatively short amount of time. In practice, our criteria will
pick out gas that is collapsing to very high densities on a timescale shorter
than the typical star formation timescale (assumed to be $10^{6}$ yr) and
faster than it can effectively cool. The additional criteria that this gas
must be (nearly) pristine means that such locations must form prior to or
simultaneously with the very first stars forming in the (proto-)galaxy.
Despite the simplicity of the model, the connection between the high-redshift
properties of (proto-)galaxies and their future assembly history and
environment remains a prediction of the model, rather than a direct
consequence of our choice in criteria. Still, we do not attempt to test
different criteria and it is very possible that this would influence our
results. For example, softening the metallicity or density requirements would
make MBHs much more common likely wash out any environmental dependence.

It should be noted that this model is only capable of capturing MBH formation
channels that take place in pristine (or near-pristine) gas. This is primarily
what results in an early formation epoch, as most gas becomes polluted as
stars form in the simulation. However, it may be possible to grow MBHs at
later times in metal enriched gas, either growing a low mass seed within star
clusters (Natarajan, 2021) or in massive merger events (Mayer et al., 2023).
The formation model implemented in Romulus would not capture such channels.

The fact that our theoretical results produce active fractions consistent with
observations (see Figure 2 and Sharma et al., 2022a) indicates that our model
parameters are reasonable, at least. However, Sharma et al. (2022b) find that
AGN feedback in dwarf galaxies is the primary cause for over-quenching low
mass galaxies. This could indicate that our occupation fractions are too high
in the field, though this is just as likely an issue with overly efficient MBH
accretion and/or feedback. In any case, a more strict formation criteria would
decrease occupation fractions of MBHs and could potentially increase the
divide between environments even further.

###  5.2 Halo Finder Limitations

Our analysis has been limited by our ability to extract halos in consecutive
timesteps. The difficulty of extracting substructure close to the centers of
dense structures like clusters is a well-known issue with halo finding
routines. While this may result in an artificial lack of low mass galaxies
deep within the cluster, this should not effect our overall results on the
occupation fraction of cluster dwarf galaxies. In fact, galaxies closer to the
center of the cluster are likely to have fallen in earlier and are therefore
more likely to have hosted a MBH, so including more central dwarfs could
increase our cluster occupation fractions further.

More important is the effect on our ability to trace halos backward in time.
This requires that a given halo is detected in all timesteps while it is in
the cluster, which may not happen if it passes close to the center at any
point. Given the wide mass range we examine and the fact that we are able to
successfully trace back the majority of even the smallest galaxies, these
missed galaxies should not affect our conclusions. In fact, we should
preferentially miss the earliest forming dwarf galaxies that fall into the
cluster first, which would only strengthen the effect that we see already.

###  5.3 Effect of Resolution

An important numerical affect caused by limited resolution is the artificial
disruption of dark matter halos (and galaxies). As discussed in detail in van
den Bosch & Ogiya (2018), limited particle count and gravity resolution
results in substructure becoming artificially disrupted. In the simulation,
most galaxies that do experience significant tidal stripping of their stars
are very soon completely disrupted while, in reality, it is possible they
should survive. This would mean that we underestimate the portion of MBHs that
exist in significantly stripped galaxies. This would likely further increase
the difference in occupation fraction we already see in the cluster, as we
would have an additional population of dwarf galaxies hosting MBHs that we
currently do not resolve. It would also mean that tidal stripping is more
important than we currently predict to the overall MBH population in cluster
dwarfs (see Voggel et al., 2019, for more discussion on this population of
MBHs based on observations).

Similarly, limited resolution means that Romulus will not resolve dense
stellar structures, such as nuclear star clusters, at the centers of galaxies.
These structures are more resilient to tidal effects, so even if the
disruption of subhalos were all real it is possible that some of these
structures would survive around the MBHs as ultra-compact dwarf galaxies
(Drinkwater et al., 2003; Bekki et al., 2003; Seth et al., 2014). Similar to
artificial disruption, the effect of this would be that there is a more
significant population of tidally stripped remnants that host MBHs. It would
also further increase the predicted environmental dependence of the occupation
fraction by, once again, adding a new population of dwarf galaxies hosting
MBHs to our sample. Further, this would create a population of galaxies with
significantly overmassive black holes. As discussed in Ricarte et al. (2019),
the accretion histories of MBHs in RomulusC dwarf galaxies is very similar to
that of field galaxies, and so galaxies in both environments exist on the same
stellar mass-black hole mass relation. This might not remain the case if the
number of artificially disrupted galaxies/nuclear star clusters is accounted
for, but we leave this question to future work more focused on explaining
observations of MBHs in ultra-compact dwarf galaxies (e.g. Seth et al., 2014;
Afanasiev et al., 2018; Voggel et al., 2019).

###  5.4 Connection to Nuclear Star Clusters

The MBH occupation fractions we predict with our relatively simple seed
formation model matches remarkably well with observations of NSCs in cluster
dwarf galaxies (Muñoz et al., 2015; Sánchez-Janssen et al., 2019), as well as
local group dwarfs (Hoyer et al., 2021; Carlsten et al., 2022). There is
reason to think that the formation of NSCs and MBHs could be connected. It may
be that MBHs are seeded as a result of the evolution of a dense nuclear star
cluster (Devecchi & Volonteri, 2009; Davies et al., 2011; Kroupa et al.,
2020). More broadly, the environment likely to form/grow a MBH (very dense,
pristine gas) is also the site of very dense, early star formation that seeds
an initial NSC (note that NSCs often include stars with variety of ages and
metalicities, indicated extended star formation histories Seth et al., 2006;
Carson et al., 2015; Kacharov et al., 2018). Some work suggests that NSCs and
MBHs form from entirely separate mechanisms (e.g. Scott & Graham, 2013) or
that NSCs may grow mostly from mergers of other star clusters, rather than in-
situ formation (Antonini et al., 2015; Fahrion et al., 2020). In reality, it
is likely that a combination of mechanisms are occurring (e.g. Fahrion et al.,
2021; Fahrion et al., 2022a).

Many NSCs co-exist with MBHs and their masses scale with the mass and
properties of their host galaxies in similar ways (Wehner & Harris, 2006;
Ferrarese et al., 2006; Seth et al., 2008; Georgiev et al., 2016; Nguyen et
al., 2018, 2019; Fahrion et al., 2022b). The potential well of a NSC can help
low mass MBHs grow (Natarajan, 2021; Askar et al., 2022). Conversely, feedback
from MBHs, as well as their dynamical interactions with nearby stars, can
hinder the growth of NSCs and even disrupt them completely (Antonini, 2013;
Antonini et al., 2015; Sánchez-Janssen et al., 2019). The presence of MBHs
within dense, nuclear regions are thought to explain the presence of
supermassive black holes in ultra-compact dwarf galaxies typically found in
groups and clusters (Seth et al., 2014; Afanasiev et al., 2018; Voggel et al.,
2019).

The formation criteria for MBHs in Romulus is relatively agnostic regarding
the exact physics as we are far from resolving the formation process itself.
Rather, it is meant to encapsulate the type of environment where one may
expect a MBH to form and grow rapidly in the early Universe, i.e. where cold,
low metalicity gas is collapsing on timescales much shorter than the star
formation timecale. Such conditions are required for all formation channels of
MBHs. Even starting as Population III remnants, the seeds would need a lot of
dense gas nearby to quickly grow large enough to become $10^{5}$ M⊙ in a short
amount of time (see Volonteri, 2012, for further discussion on this). While we
are far from being able to resolve NSCs in the simulation, the connection to
observed NSC occupation makes sense regardless of the details of NSC formation
and suggests that NSCs might serve as incubators for the formation and growth
of MBH seeds over cosmic time as noted by Natarajan (2021) . An in-situ
formation channel would require similar properties to MBH formation (very
dense gas in the early Universe) but a formation channel dominated by globular
cluster mergers (e.g. Antonini & Merritt, 2012; Antonini et al., 2015; Fahrion
et al., 2020) would also fit with a connection to MBHs. Romulus predicts that
the environment required for MBH formation is more likely in galaxies (and
dark matter halos) that form earlier. Such galaxies will be likely to form
more globular clusters that will then have more time to sink and merge and
form a NSC. Observations have shown an overabundance of globular clusters in
dwarf ellipticals residing in galaxy cluster environments (e.g. Miller et al.,
1998; Miller & Lotz, 2007; Jordán et al., 2007; Peng et al., 2008; Sánchez-
Janssen & Aguerri, 2012) and the cause of this can be attributed to their
earlier formation times (Mistani et al., 2016; Carleton et al., 2021).

The fact that our results match so well with observed NSC populations across
environments (Muñoz et al., 2015; Sánchez-Janssen et al., 2019; Hoyer et al.,
2021; Carlsten et al., 2022) while producing realistic black hole occupation
fractions in the field (though observational constraints are murky at best)
supports the notion that MBH formation could be connected with NSC formation.
While we cannot directly resolve the formation of NSCs, our results indicate
that their presence, like that of MBHs, may be connected not only to the
properties of gas at high redshift, but to the overall formation history of
galaxies and therefore, indirectly, their environment.

##  6 Conclusions

We use the Romulus simulations to predict an enhanced MBH occupation fraction
in cluster dwarf galaxies ($M_{\star}<10^{9}$M⊙). The Romulus simulations are
unique in their ability to both resolve low mass galaxies and implement a
model for MBH seeding that relies only on local gas properties (density,
temperature, metallicity) rather than requiring any a priori assumptions about
which halos should or should not host a MBH. Despite forming black holes at
similar times and within similar galaxy masses, we find that the cosmic
evolution of the MBH occupation fraction in galaxies is halted at $z\sim 3$ in
the cluster environment relative to the field. This ‘freezing out’ of the
occupation fraction results in a factor of $\sim 2$ enhancement to the
fraction of dwarf galaxies that host MBHs in clusters relative to the field at
$z=0.05$.

We investigate the cause of the enhancement in more detail and find that it
can likely be explained by a combination of two mechanisms:

  1. 1.

Early formation times of dwarf galaxies in cluster environments makes them
more likely to host a MBH. When controlling for formation time, cluster and
field galaxies have similar MBH occupation fractions, but late-forming dwarf
galaxies that do not host MBHs dillute the field population and pull down the
MBH occupation fraction, while these systems do not exist in clusters.

  2. 2.

Some cluster dwarf galaxies may be ‘failed’ massive galaxies. Were they
allowed to grow in the field, unimpeded by the cluster environment, they would
likely have attained much higher masses.

We do not find evidence that many MBH host galaxies experience tidal stripping
of their stars. However, future work will examine whether some MBHs in the
simulation may have unresolved, compact stellar structures around them.

The enhanced MBH occupation fraction in the cluster simulation appears to fall
with increased cluster-centric distance, but it remains in place out to
$2R_{200}$ for galaxies with $M_{\star}<10^{9}$M⊙. While we do not attempt to
model the detailed physics of MBH formation in these simulations, the
connection between galaxy assembly history, environment, and the properties of
gas at high redshift are important predictions. The presence of quickly
collapsing, high density regions of pristine gas are likely formation sites of
MBH seeds. While such environments are common in the progenitors of massive
galaxies, we predict that only early forming dwarf galaxies are typically able
to produce such environments. Such dwarf galaxies make up a much higher
fraction of the overall population in dense environments like galaxy clusters.
Our findings have important consequences for the origin and evolution of host
galaxy-MBH co-evolution.

Finally, we find that the predicted MBH occupation fraction in the cluster is
remarkably consistent with the observed occupation fraction of nuclear star
clusters in Virgo and Fornax, while the field occupation fraction is similar
to that of dwarfs in the local volume. So, the dense, pristine gas at $z>5$
that the simulation attributes to MBH formation may also be connected to the
formation of nuclear star clusters, either directly (they form from similar
gas at similar times) or indirectly (e.g. connected to the earlier formation
times associated with the host galaxies).

These results show how MBH occupation may be not just a function of galaxy
mass but also environment, and that there are many connections between the
high redshift properties of galaxies and their overall assembly history. We
note that RomulusC is a relatively low mass cluster and more massive ones may
have even more enhancement in their occupation fractions. These results also
indicate that we must be cautious in how we utilize AGN observations to
constrain the underlying MBH occupation fraction, as the connection is heavily
dependent on environment and, more generally, the age of the galaxy (i.e. its
formation time). Further work and higher resolution simulations will be needed
to better understand how the environmental dependence of the underlying MBH
occupation fraction may be inferred from observations. More broadly, these
results demonstrate how simulations with more predictive models for MBH
physics are crucial to our understanding of MBH formation and evolution.

## Acknowledgements

This work used the pynbody (Pontzen et al., 2013) and Tangos (Pontzen &
Tremmel, 2018) software packages. MT was supported by an NSF Astronomy and
Astrophysics Postdoctoral Fellowship under award AST-2001810. PN gratefully
acknowledges support from the Black Hole Initiative at Harvard University,
which is funded by grants from the John Templeton Foundation and the Gordon
and Betty Moore Foundation. The authors thank Frank van den Bosch, Marta
Volonteri, Mark Kennedy, Paul Callanan, and Katja Fahrion for useful
discussions and comments.

## Data Availability

The Tangos database for both Romulus simulations are available upon email
request. Recreaction of all of the figures in this paper with processed data
(derived using Tangos databases) can be done using the jupyter notebook, code,
and data files that are available at
https://github.com/mtremmel/tremmel2023_mbh_occFrac

## References

  * Afanasiev et al. (2018) Afanasiev A. V., et al., 2018, MNRAS, 477, 4856
  * Alexander & Natarajan (2014) Alexander T., Natarajan P., 2014, Science, 345, 1330
  * Amaro-Seoane et al. (2023) Amaro-Seoane P., et al., 2023, Living Reviews in Relativity, 26, 2
  * Antonini (2013) Antonini F., 2013, ApJ, 763, 62
  * Antonini & Merritt (2012) Antonini F., Merritt D., 2012, ApJ, 745, 83
  * Antonini et al. (2015) Antonini F., Barausse E., Silk J., 2015, ApJ, 812, 72
  * Applebaum et al. (2021) Applebaum E., Brooks A. M., Christensen C. R., Munshi F., Quinn T. R., Shen S., Tremmel M., 2021, ApJ, 906, 96
  * Askar et al. (2022) Askar A., Davies M. B., Church R. P., 2022, MNRAS, 511, 2631
  * Baldassare et al. (2015) Baldassare V. F., Reines A. E., Gallo E., Greene J. E., 2015, ApJ, 809, L14
  * Baldassare et al. (2016) Baldassare V. F., et al., 2016, ApJ, 829, 57
  * Baldassare et al. (2018) Baldassare V. F., Geha M., Greene J., 2018, preprint, p. arXiv:1808.09578 (arXiv:1808.09578)
  * Baldassare et al. (2020) Baldassare V. F., Geha M., Greene J., 2020, ApJ, 896, 10
  * Beckmann et al. (2023) Beckmann R. S., et al., 2023, MNRAS,
  * Bekki et al. (2003) Bekki K., Couch W. J., Drinkwater M. J., Shioya Y., 2003, MNRAS, 344, 399
  * Bellovary et al. (2019) Bellovary J., et al., 2019, BAAS, 51, 175
  * Bellovary et al. (2021) Bellovary J. M., et al., 2021, MNRAS, 505, 5129
  * Birchall et al. (2020) Birchall K. L., Watson M. G., Aird J., 2020, MNRAS, 492, 2268
  * Blank et al. (2019) Blank M., Macciò A. V., Dutton A. A., Obreja A., 2019, MNRAS, 487, 5476
  * Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273
  * Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
  * Borrow et al. (2022) Borrow J., Vogelsberger M., O’Neil S., McDonald M. A., Smith A., 2022, arXiv e-prints, p. arXiv:2205.10376
  * Bortolas et al. (2022) Bortolas E., Bonetti M., Dotti M., Lupi A., Capelo P. R., Mayer L., Sesana A., 2022, MNRAS, 512, 3365
  * Boylan-Kolchin et al. (2009) Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150
  * Bricman & Gomboc (2020) Bricman K., Gomboc A., 2020, ApJ, 890, 73
  * Burke et al. (2022) Burke C. J., et al., 2022, MNRAS, 516, 2736
  * Cameron (2011) Cameron E., 2011, Publications of the Astronomical Society of Australia, 28, 128
  * Cann et al. (2021) Cann J. M., et al., 2021, ApJ, 912, L2
  * Carleton et al. (2021) Carleton T., Guo Y., Munshi F., Tremmel M., Wright A., 2021, MNRAS, 502, 398
  * Carlsten et al. (2022) Carlsten S. G., Greene J. E., Beaton R. L., Greco J. P., 2022, ApJ, 927, 44
  * Carson et al. (2015) Carson D. J., Barth A. J., Seth A. C., den Brok M., Cappellari M., Greene J. E., Ho L. C., Neumayer N., 2015, AJ, 149, 170
  * Chadayammuri et al. (2021) Chadayammuri U., Tremmel M., Nagai D., Babul A., Quinn T., 2021, MNRAS, 504, 3922
  * Chandrasekhar (1943) Chandrasekhar S., 1943, ApJ, 97, 255
  * Chaves-Montero et al. (2016) Chaves-Montero J., Angulo R. E., Schaye J., Schaller M., Crain R. A., Furlong M., Theuns T., 2016, MNRAS, 460, 3100
  * Chen et al. (2022) Chen N., Ni Y., Tremmel M., Di Matteo T., Bird S., DeGraf C., Feng Y., 2022, MNRAS, 510, 531
  * Christensen et al. (2014) Christensen C. R., Governato F., Quinn T., Brooks A. M., Shen S., McCleary J., Fisher D. B., Wadsley J., 2014, MNRAS, 440, 2843
  * Churazov et al. (2005) Churazov E., Sazonov S., Sunyaev R., Forman W., Jones C., Böhringer H., 2005, MNRAS, 363, L91
  * Colpi et al. (2019) Colpi M., et al., 2019, arXiv e-prints, p. arXiv:1903.06867
  * Croton et al. (2007) Croton D. J., Gao L., White S. D. M., 2007, MNRAS, 374, 1303
  * Davies et al. (2011) Davies M. B., Miller M. C., Bellovary J. M., 2011, ApJ, 740, L42
  * Devecchi & Volonteri (2009) Devecchi B., Volonteri M., 2009, ApJ, 694, 302
  * Di Matteo et al. (2008) Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33
  * Di Matteo et al. (2017) Di Matteo T., Croft R. A. C., Feng Y., Waters D., Wilkins S., 2017, MNRAS, 467, 4243
  * Drinkwater et al. (2003) Drinkwater M. J., Gregg M. D., Hilker M., Bekki K., Couch W. J., Ferguson H. C., Jones J. B., Phillipps S., 2003, Nature, 423, 519
  * Dubois et al. (2016) Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS,
  * Dubois et al. (2021) Dubois Y., et al., 2021, A&A, 651, A109
  * Dunn et al. (2018) Dunn G., Bellovary J., Holley-Bockelmann K., Christensen C., Quinn T., 2018, ApJ, 861, 39
  * Fahrion et al. (2020) Fahrion K., et al., 2020, A&A, 634, A53
  * Fahrion et al. (2021) Fahrion K., et al., 2021, A&A, 650, A137
  * Fahrion et al. (2022a) Fahrion K., Leaman R., Lyubenova M., van de Ven G., 2022a, A&A, 658, A172
  * Fahrion et al. (2022b) Fahrion K., et al., 2022b, A&A, 667, A101
  * Ferrarese et al. (2006) Ferrarese L., et al., 2006, ApJ, 644, L21
  * Gallo et al. (2010) Gallo E., Treu T., Marshall P. J., Woo J.-H., Leipski C., Antonucci R., 2010, ApJ, 714, 25
  * Gao & White (2007) Gao L., White S. D. M., 2007, MNRAS, 377, L5
  * Gao et al. (2005) Gao L., Springel V., White S. D. M., 2005, MNRAS, 363, L66
  * Geha et al. (2012) Geha M., Blanton M. R., Yan R., Tinker J. L., 2012, ApJ, 757, 85
  * Georgiev et al. (2016) Georgiev I. Y., Böker T., Leigh N., Lützgendorf N., Neumayer N., 2016, MNRAS, 457, 2122
  * Governato et al. (2015) Governato F., et al., 2015, MNRAS, 448, 792
  * Greene (2012) Greene J. E., 2012, Nature Communications, 3, 1304
  * Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76
  * Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  * Habouzit et al. (2021) Habouzit M., et al., 2021, MNRAS,
  * Habouzit et al. (2022) Habouzit M., et al., 2022, MNRAS, 509, 3015
  * Haidar et al. (2022) Haidar H., et al., 2022, MNRAS, 514, 4912
  * Henden et al. (2018) Henden N. A., Puchwein E., Shen S., Sijacki D., 2018, MNRAS, 479, 5385
  * Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
  * Hoyer et al. (2021) Hoyer N., Neumayer N., Georgiev I. Y., Seth A. C., Greene J. E., 2021, MNRAS, 507, 3246
  * Inayoshi et al. (2016) Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS, 459, 3738
  * Jordán et al. (2007) Jordán A., et al., 2007, ApJS, 171, 101
  * Joshi et al. (2016) Joshi G. D., Parker L. C., Wadsley J., 2016, MNRAS, 462, 761
  * Jung et al. (2022) Jung S. L., et al., 2022, MNRAS, 515, 22
  * Kacharov et al. (2018) Kacharov N., Neumayer N., Seth A. C., Cappellari M., McDermid R., Walcher C. J., Böker T., 2018, MNRAS, 480, 1973
  * Knebe et al. (2011) Knebe A., et al., 2011, MNRAS, 415, 2293
  * Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kormendy & Richstone (1995) Kormendy J., Richstone D., 1995, ARA&A, 33, 581
  * Koudmani et al. (2021) Koudmani S., Henden N. A., Sijacki D., 2021, MNRAS, 503, 3568
  * Koudmani et al. (2022) Koudmani S., Sijacki D., Smith M. C., 2022, MNRAS, 516, 2112
  * Kroupa et al. (2020) Kroupa P., Subr L., Jerabkova T., Wang L., 2020, MNRAS, 498, 5652
  * Latimer et al. (2021) Latimer L. J., Reines A. E., Bogdan A., Kraft R., 2021, ApJ, 922, L40
  * Leja et al. (2019) Leja J., et al., 2019, ApJ, 877, 140
  * Lescaudron et al. (2022) Lescaudron S., Dubois Y., Beckmann R. S., Volonteri M., 2022, arXiv e-prints, p. arXiv:2209.13548
  * Lodato & Natarajan (2006) Lodato G., Natarajan P., 2006, MNRAS, 371, 1813
  * Lodato & Natarajan (2007) Lodato G., Natarajan P., 2007, MNRAS, 377, L64
  * Mayer et al. (2023) Mayer L., Capelo P. R., Zwick L., Di Matteo T., 2023, arXiv e-prints, p. arXiv:2304.02066
  * Menon et al. (2015) Menon H., Wesolowski L., Zheng G., Jetley P., Kale L., Quinn T., Governato F., 2015, Comp. Astrophysics and Cosmology, 2, 1
  * Mezcua et al. (2018) Mezcua M., Civano F., Marchesi S., Suh H., Fabbiano G., Volonteri M., 2018, MNRAS, 478, 2576
  * Miller & Lotz (2007) Miller B. W., Lotz J. M., 2007, ApJ, 670, 1074
  * Miller et al. (1998) Miller B. W., Lotz J. M., Ferguson H. C., Stiavelli M., Whitmore B. C., 1998, ApJ, 508, L133
  * Miller et al. (2015) Miller B. P., Gallo E., Greene J. E., Kelly B. C., Treu T., Woo J.-H., Baldassare V., 2015, ApJ, 799, 98
  * Mistani et al. (2016) Mistani P. A., et al., 2016, MNRAS, 455, 2323
  * Molina et al. (2021) Molina M., Reines A. E., Latimer L. J., Baldassare V., Salehirad S., 2021, ApJ, 922, 155
  * Muñoz et al. (2015) Muñoz R. P., et al., 2015, ApJ, 813, L15
  * Munshi et al. (2013) Munshi F., et al., 2013, ApJ, 766, 56
  * Munshi et al. (2021) Munshi F., Brooks A. M., Applebaum E., Christensen C. R., Quinn T., Sligh S., 2021, ApJ, 923, 35
  * Natarajan (2011) Natarajan P., 2011, Bulletin of the Astronomical Society of India, 39, 145
  * Natarajan (2014) Natarajan P., 2014, General Relativity and Gravitation, 46, 1702
  * Natarajan (2021) Natarajan P., 2021, MNRAS, 501, 1413
  * Natarajan et al. (2017) Natarajan P., Pacucci F., Ferrara A., Agarwal B., Ricarte A., Zackrisson E., Cappelluti N., 2017, ApJ, 838, 117
  * Nelson et al. (2019) Nelson D., et al., 2019, MNRAS, 490, 3234
  * Nguyen et al. (2018) Nguyen D. D., et al., 2018, ApJ, 858, 118
  * Nguyen et al. (2019) Nguyen D. D., et al., 2019, ApJ, 872, 104
  * Ni et al. (2022) Ni Y., et al., 2022, MNRAS, 513, 670
  * Okamoto et al. (2008) Okamoto T., Nemmen R. S., Bower R. G., 2008, MNRAS, 385, 161
  * Onions et al. (2012) Onions J., et al., 2012, MNRAS, 423, 1200
  * Ostriker (1999) Ostriker E. C., 1999, ApJ, 513, 252
  * Pacucci et al. (2019) Pacucci F., et al., 2019, BAAS, 51, 117
  * Pacucci et al. (2021) Pacucci F., Mezcua M., Regan J. A., 2021, ApJ, 920, 134
  * Peng et al. (2008) Peng E. W., et al., 2008, ApJ, 681, 197
  * Pfister et al. (2019) Pfister H., Volonteri M., Dubois Y., Dotti M., Colpi M., 2019, MNRAS, 486, 101
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Poggianti et al. (2017) Poggianti B. M., et al., 2017, Nature, 548, 304
  * Pontzen & Tremmel (2018) Pontzen A., Tremmel M., 2018, ApJS, 237, 23
  * Pontzen et al. (2008) Pontzen A., et al., 2008, MNRAS, 390, 1349
  * Pontzen et al. (2013) Pontzen A., Roškar R., Stinson G., Woods R., 2013, pynbody: N-Body/SPH analysis for python, Astrophysics Source Code Library, record ascl 1305.002 (ascl:1305.002)
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Regan (2023) Regan J., 2023, The Open Journal of Astrophysics, 6, 12
  * Regan et al. (2017) Regan J. A., Visbal E., Wise J. H., Haiman Z., Johansson P. H., Bryan G. L., 2017, Nature Astronomy, 1, 0075
  * Regan et al. (2020a) Regan J. A., Haiman Z., Wise J. H., O’Shea B. W., Norman M. L., 2020a, The Open Journal of Astrophysics, 3, E9
  * Regan et al. (2020b) Regan J. A., Wise J. H., Woods T. E., Downes T. P., O’Shea B. W., Norman M. L., 2020b, The Open Journal of Astrophysics, 3, 15
  * Regan et al. (2023) Regan J. A., Pacucci F., Bustamante-Rosell M. J., 2023, MNRAS, 518, 5997
  * Reines & Deller (2012) Reines A. E., Deller A. T., 2012, ApJ, 750, L24
  * Reines et al. (2011) Reines A. E., Sivakoff G. R., Johnson K. E., Brogan C. L., 2011, Nature, 470, 66
  * Reines et al. (2020) Reines A. E., Condon J. J., Darling J., Greene J. E., 2020, ApJ, 888, 36
  * Ricarte & Natarajan (2018) Ricarte A., Natarajan P., 2018, MNRAS, 481, 3278
  * Ricarte et al. (2019) Ricarte A., Tremmel M., Natarajan P., Quinn T., 2019, MNRAS, 489, 802
  * Ricarte et al. (2020) Ricarte A., Tremmel M., Natarajan P., Quinn T., 2020, ApJ, 895, L8
  * Ricarte et al. (2021a) Ricarte A., Tremmel M., Natarajan P., Zimmer C., Quinn T., 2021a, MNRAS, 503, 6098
  * Ricarte et al. (2021b) Ricarte A., Tremmel M., Natarajan P., Quinn T., 2021b, ApJ, 916, L18
  * Ritchie & Thomas (2001) Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743
  * Rosas-Guevara et al. (2016) Rosas-Guevara Y., Bower R. G., Schaye J., McAlpine S., Dalla Vecchia C., Frenk C. S., Schaller M., Theuns T., 2016, MNRAS, 462, 190
  * Roškar et al. (2015) Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J., 2015, MNRAS, 449, 494
  * Saitoh & Makino (2009) Saitoh T. R., Makino J., 2009, ApJ, 697, L99
  * Sánchez-Janssen & Aguerri (2012) Sánchez-Janssen R., Aguerri J. A. L., 2012, MNRAS, 424, 2614
  * Sánchez-Janssen et al. (2019) Sánchez-Janssen R., et al., 2019, ApJ, 878, 18
  * Sassano et al. (2023) Sassano F., Capelo P. R., Mayer L., Schneider R., Valiante R., 2023, MNRAS, 519, 1837
  * Scott & Graham (2013) Scott N., Graham A. W., 2013, ApJ, 763, 76
  * Sesana (2013) Sesana A., 2013, MNRAS, 433, L1
  * Sesana & Khan (2015) Sesana A., Khan F. M., 2015, MNRAS, 454, L66
  * Seth et al. (2006) Seth A. C., Dalcanton J. J., Hodge P. W., Debattista V. P., 2006, AJ, 132, 2539
  * Seth et al. (2008) Seth A., Agüeros M., Lee D., Basu-Zych A., 2008, ApJ, 678, 116
  * Seth et al. (2014) Seth A. C., et al., 2014, Nature, 513, 398
  * Sharma et al. (2020) Sharma R. S., Brooks A. M., Somerville R. S., Tremmel M., Bellovary J., Wright A. C., Quinn T. R., 2020, ApJ, 897, 103
  * Sharma et al. (2022a) Sharma R. S., Brooks A. M., Tremmel M., Bellovary J., Ricarte A., Quinn T. R., 2022a, arXiv e-prints, p. arXiv:2203.05580
  * Sharma et al. (2022b) Sharma R. S., Brooks A. M., Tremmel M., Bellovary J., Quinn T. R., 2022b, arXiv e-prints, p. arXiv:2211.05275
  * Shen et al. (2010) Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581
  * Shen et al. (2020) Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., Richards G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252
  * Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  * Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074
  * Trebitsch et al. (2021) Trebitsch M., et al., 2021, A&A, 653, A154
  * Tremaine et al. (2002) Tremaine S., et al., 2002, ApJ, 574, 740
  * Tremmel et al. (2015) Tremmel M., Governato F., Volonteri M., Quinn T. R., 2015, MNRAS, 451, 1868
  * Tremmel et al. (2017) Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T. R., Pontzen A., Anderson L., Bellovary J., 2017, MNRAS, 470, 1121
  * Tremmel et al. (2018a) Tremmel M., Governato F., Volonteri M., Quinn T. R., Pontzen A., 2018a, MNRAS, 475, 4967
  * Tremmel et al. (2018b) Tremmel M., Governato F., Volonteri M., Pontzen A., Quinn T. R., 2018b, The Astrophysical Journal Letters, 857, L22
  * Tremmel et al. (2019) Tremmel M., et al., 2019, MNRAS, 483, 3336
  * Tremmel et al. (2020) Tremmel M., Wright A. C., Brooks A. M., Munshi F., Nagai D., Quinn T. R., 2020, MNRAS, 497, 2786
  * Van Nest et al. (2022) Van Nest J. D., Munshi F., Wright A. C., Tremmel M., Brooks A. M., Nagai D., Quinn T., 2022, ApJ, 926, 92
  * Voggel et al. (2016) Voggel K., Hilker M., Richtler T., 2016, A&A, 586, A102
  * Voggel et al. (2019) Voggel K. T., Seth A. C., Baumgardt H., Mieske S., Pfeffer J., Rasskazov A., 2019, ApJ, 871, 159
  * Volonteri (2010) Volonteri M., 2010, A&A Rev., 18, 279
  * Volonteri (2012) Volonteri M., 2012, Science, 337, 544
  * Volonteri & Gnedin (2009) Volonteri M., Gnedin N. Y., 2009, ApJ, 703, 2113
  * Volonteri & Natarajan (2009) Volonteri M., Natarajan P., 2009, MNRAS, 400, 1911
  * Volonteri & Rees (2005) Volonteri M., Rees M. J., 2005, ApJ, 633, 624
  * Volonteri et al. (2008) Volonteri M., Lodato G., Natarajan P., 2008, MNRAS, 383, 1079
  * Volonteri et al. (2020) Volonteri M., et al., 2020, MNRAS, 498, 2219
  * Wadsley et al. (2004) Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy, 9, 137
  * Wadsley et al. (2008) Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427
  * Wadsley et al. (2017) Wadsley J. W., Keller B. W., Quinn T. R., 2017, MNRAS, 471, 2357
  * Wehner & Harris (2006) Wehner E. H., Harris W. E., 2006, ApJ, 644, L17
  * Wise et al. (2019) Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., Downes T. P., Xu H., 2019, Nature, 566, 85
  * Woo et al. (2019) Woo J.-H., Cho H., Gallo E., Hodges-Kluck E., Le H. A. N., Shin J., Son D., Horst J. C., 2019, Nature Astronomy, 3, 755
  * Wright et al. (2021) Wright A. C., Tremmel M., Brooks A. M., Munshi F., Nagai D., Sharma R. S., Quinn T. R., 2021, MNRAS, 502, 5370
  * van Wassenhove et al. (2010) van Wassenhove S., Volonteri M., Walker M. G., Gair J. R., 2010, MNRAS, 408, 1139
  * van de Voort et al. (2011) van de Voort F., Schaye J., Booth C. M., Haas M. R., Dalla Vecchia C., 2011, MNRAS, 414, 2458
  * van den Bosch & Ogiya (2018) van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066