AN ENHANCED MASSIVE BLACK HOLE OCCUPATION PREDICTED IN CLUSTER DWARF GALAXIES

The occupation fraction of massive black holes (MBHs) in low mass galaxies o ff ers interesting insights into initial black hole seeding mechanisms and their mass assembly history, though disentangling these two e ff ects remains challenging. Using the R omulus cosmological simulations we examine the impact of environment on the occupation fraction of MBHs in low mass galaxies. Unlike most modern cosmological simulations, R omulus seeds MBHs based on local gas properties, selecting very dense, pristine, 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 ⋆ < 10 9 M ⊙ in cluster environments are ∼ 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 di ff erence 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.


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 ∼ 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 & Corresponding Author: mtremmel@ucc.ieVolonteri 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 ⋆ < 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. 2015Baldassare 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. 2008Matteo et al. , 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. 2021Koudmani et al. , 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(Tremmel et al. , 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(Bellovary et al. , 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(Tremmel et al. , 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(Sharma et al. , 2022b)), something which is seen in other simulations as well (Koudmani et al. 2021(Koudmani et al. , 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-perside 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.

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. (2017Tremmel et al. ( , 2019)).

Overview of the Romulus Simulations
The Romulus Simulations (Tremmel et al. 2017(Tremmel et al. , 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 (Romu-lusC).With a dark matter mass resolution of 3.4 × 10 5 M ⊙ , both Romulus simulations are able to resolve halos as small as ∼ 3×10 9 M ⊙ with more than 10,000 particles.With typical gas and star particle masses of 2.1 × 10 5 and 6 × 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 ∼ 10 7 M ⊙ in stellar mass with hundreds of resolution elements.Both simulations are run with the same cosmology (Ω 0 = 0.3086, Λ = 0.6914, h = 0.6777, σ 8 = 0.8288; Planck Collaboration et al. 2016) and physics.

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: where m star = 0.3m gas , c ⋆ = 0.15, and t 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 'ultradiffuse' 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.

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: where G is the gravitational constant; n * is the star formation threshold (0.2 cm −3 ; see previous section); v bulk is the bulk velocity of the gas; v θ is its rotational velocity; c s is its sound speed; ρ is its mass density; and β is set to 2. The radiative efficiency (ϵ 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 where ϵ 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(Tremmel et al. , 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(Sharma et al. , 2022a,b) ,b) 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).

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, ϵ g , during each black hole timestep.
Here v • is the velocity of the black hole relative to nearby star and dark matter particles, ρ(< v • ) is the local density of star and dark matter particles that are moving slower than the black hole relative the background particles, and ln Λ is the Coulomb logarithm and equal to ln(ϵ g /r 90 ).Here r 90 is the 90 • 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 × ϵ 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).

Black Hole Seeding
A critical aspect of the Romulus simulations to this work is the prescription for MBH seeding.Like most modern, largescale 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 × 10 −4 Z ⊙ .
• Temperature is between 9500 K and 10000 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 ×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.

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 (R vir ) and mass (M vir ) 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 R 200 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 × 10 9 M ⊙ , such that each halo included in our analysis has at least ∼ 10, 000 particles.An exception is made for our analysis of RomulusC where we include halos down to 3 × 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 (∼ 2R 200 ) 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 ⋆ < 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 R 200 of any more massive halo.Additionally, for low mass galaxies with M ⋆ < 10 10 M ⊙ , they must be further than 1.5 Mpc from any galaxy with stellar mass greater than 2.5 × 10 10 M ⊙ .

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 (Romu-lusC) and in the field (Romulus25).For Romulus25, only galaxies within halos of mass > 3 × 10 9 M ⊙ are included.For RomulusC, only galaxies in halos of mass > 3 × 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 Romu-lus25 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 ⋆ < 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.
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 Edd ) accretion: This assumes the same constant radiative efficiency (ϵ 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 ϵ 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 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 ⋆ < 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 × 10 38 erg 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.
emitted light couples to nearby gas as thermal energy, which would no longer be true when using equation 5.
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 Xray 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 offcenter and indeed often are (Tremmel et al. 2018a;Tremmel et al. 2018b;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.

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 (∼ 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 ∼ 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 ∼ 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.

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 R 200 ) 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 2R 200 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;Tremmel et al. 2018b;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.

EXPLORING THE ORIGIN OF THE ENVIRONMENTAL DEPENDENCE OF DWARF GALAXY MBH OCCUPATION
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 R 200 , 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 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 ∼ 3. The large error bars are due to small number statistics at the highest mass bins at earlier times.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. 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 ultracompact 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).

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).
Focusing only on MBHs hosted in cluster dwarf galaxies (M ⋆ < 10 10 M ⊙ ) at z = 0.05, we find that the median MBH host galaxy has only experienced a net loss of ∼ 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 ∼ 3. Looking at figure 2, a typical stellar mass loss of a factor of ∼ 3 − 5 is needed to bring the occupation fractions in line with the field.In Romu-lusC such extreme mass loss is too rare to fully explain the overabundance of MBHs in low mass galaxies.

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 concentration1 at the time they reach maximum virial mass (t max ).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 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 10 9 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 2R 200 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.
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.

An Overabundance of Early-forming Dwarfs in Clusters
As discussed in Sharma et al. (2020), dwarf galaxies in Ro-mulus25 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 ∼ 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 lateforming 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.
Once again, we face the problem discussed in the previous section whereby low mass (M ⋆ < 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 Fig. 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.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(Sharma et al. , 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.

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(Bellovary et al. , 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. 2017Regan et al. , 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.

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 re-Fig.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.sults, 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.

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 oc-cupation 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.

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).

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. 2021Fahrion 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. 2018Nguyen et al. , 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.

CONCLUSIONS
We use the Romulus simulations to predict an enhanced MBH occupation fraction in cluster dwarf galaxies (M ⋆ < 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 ∼ 3 in the cluster environment relative to the field.This 'freezing out' of the occupation fraction results in a factor of ∼ 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. 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. 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 ⋆ < 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 Ro-mulusC 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.
Fig. 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.

Fig. 2 .
Fig.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 inMiller et al. (2015) andGreene (2012) respectively.The black dashed line is the predicted occupation of MBHs above 10 5 M ⊙ fromAskar 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 ⋆ < 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 × 10 38 erg 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 fromPacucci 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.

Fig. 3 .
Fig. 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 ∼ 3. The large error bars are due to small number statistics at the highest mass bins at earlier times.

Fig. 4 .
Fig. 4.-The Evolution of the MBH Occupation Fraction Stalls in the Cluster Environment.Similar to Figure3here 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.

Fig. 5 .
Fig. 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 R 200 of cluster center while the red points are for galaxies in the cluster outskirts (0.75 < D/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 10 9 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 2R 200 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 fromCarlsten et al. (2022) andHoyer et al. (2021) as hatched regions.