Black Hole Growth and Feedback in Isolated ROMULUS25 Dwarf Galaxies

The Romulus suite of cosmological simulations was run using the Tree+SPH code ChaNGa (Menon et al. 2015), which inherits baryonic prescriptions from Gasoline and Gasoline2 (Wadsley et al. 2004, 2008; Stinson et al. 2006; Shen et al. 2010; Wadsley et al. 2017). In this work we analyze Romulus25, the 25 Mpc-per-side uniform box with periodic boundary conditions. We analyze Romulus25 because of its large, uniform sample of low-mass galaxies with BHs run to z = 0.

Romulus25 has comparable mass resolution to recent cosmological simulations such as IllustrisTNG (Springel et al. 2018) and Horizon-AGN (Volonteri et al. 2016), as well as force resolution comparable to the highest-resolution runs of EAGLE (Schaye et al. 2015). Romulus25 resolves gravity with a Plummer equivalent force softening of _g = 250 pc. Typically, the number of gas particles in a simulation is equal to the number of dark matter particles, while the relative masses are set according to the cosmic baryon fraction. Romulus25 instead contains 3.375× more dark matter particles than gas particles, with dark matter particles of mass 3.39 × 10⁵M_⊙ and gas particles of mass 2.12 × 10⁵M_⊙. This "oversampling" of dark matter provides Romulus25 with better-resolved BH dynamics (Tremmel et al. 2015). The Romulus suite of simulations was run with a Planck 2014 ΛCDM cosmology, with Ω_m = 0.3086, Ω_Λ = 0.6914, h = 0.6777, and σ₈ =0.8288 (Planck Collaboration et al. 2014).

Halos were identified using the Amiga Halo Finder (Knollmann & Knebe 2009) and analyzed using the simulation analysis code Pynbody (Pontzen et al. 2013). Key properties were organized into a Tangos database (Pontzen & Tremmel 2018).

We calculate halo mass, M₂₀₀, such that

$$\begin{eqnarray}&&{M}_{200}=4\pi {\rho }_{c}{{\rm{\Delta }}}_{h}{R}_{200}^{3},\end{eqnarray} \tag{ 1 }$$

where ρ_c is the critical density, Δ_h = 200 is the overdensity threshold, and R₂₀₀ is the halo radius.

Star formation in Romulus25 is regulated by the star formation efficiency, the efficiency of SN energy injection into the ISM, and the density/temperature threshold beyond which stars are allowed to form. SN feedback follows the blast wave prescription from Stinson et al. (2006).

Parameters governing star formation were constrained based on a series of 80 "zoom-in" (Governato et al. 2009) simulations of four galaxies with halo masses 10^10.5−10¹²M_⊙, run without BH physics, with the aim of reproducing a set of observed z = 0 scaling relations for galaxies. The parameter spaces were explored using the Kriging algorithm and graded by agreement with the stellar mass–halo mass relation (Moster et al. 2013); the relationship between stellar mass, angular momentum, and morphology (Obreschkow & Glazebrook 2014); and the H i gas–stellar mass relation (Cannon et al. 2011; Haynes et al. 2011). The Kriging algorithm efficiently traverses the parameter space and penalizes parameter selections that lead to deviations from the observed scaling relations. The parameters governing BH growth were afterward constrained in a similar fashion (see Section 2.4). Romulus25 adopts the following:

• 1.  
  SF efficiency, c_* = 0.15.
• 2.  
  Gas density threshold, n_* = 0.2m_p cm⁻³.
• 3.  
  SN coupling efficiency, _SN = 0.75.

Romulus25 incorporates prescriptions for metal and thermal diffusion from Shen et al. (2010) and Governato et al. (2015) and low-temperature radiative cooling as in Guedes et al. (2011). The SN feedback is based on a Kroupa (2001) initial mass function. Specifics of the physics and calibration process for Romulus25 are detailed in Tremmel et al. (2017).

Simulations commonly seed BHs by choosing a halo mass threshold above which galaxies are allowed to form a BH (Sijacki et al. 2014; Anglés-Alcázar et al. 2017). BHs are then formed from star-forming gas and placed in the halo center (Agarwal et al. 2014). Instead, Romulus25 seeds BHs based on conditions of the pre-collapse gas using characteristics of direct-collapse seeding (Haiman 2013). In this work we qualitatively define massive BHs as BHs with masses relevant to direct-collapse seeding models (M_BH ≳ 10⁵M_⊙). Romulus25 seeds massive BHs at high redshift ($z\gtrsim 5$) without assumptions of the BH occupation fraction. A gas particle is marked to form a BH if it has the following:

• 1.  
  Low metallicity, Z < 3 × 10⁻⁴.
• 2.  
  High density, 15× higher than the SF threshold.
• 3.  
  Temperature between 9500 and 10,000 K.

In other words, a gas particle will form a BH if the gas is set to collapse quickly and cool slowly, following predicted sites of direct-collapse seeding (Begelman et al. 2006; Johnson et al. 2012; Volonteri 2012; Haiman 2013; Reines & Comastri 2016). In particular, a low metallicity threshold prevents premature gas fragmentation and pushes seed formation to early times when gas has undergone little metal enrichment (Greene 2012). The values were chosen to restrict BH growth to the highest-density regions in the early universe, where BHs can undergo rapid accretion. Choosing lower metallicity or colder temperature thresholds was found to bias formation away from the densest regions (Tremmel et al. 2017).

Once these gas conditions are met, a BH is seeded at a mass of M_BH = 10⁶M_⊙. Seeding accretes mass from nearby gas particles to conserve total mass and simulate rapid, unresolved growth thought to exceed 0.1M_⊙ yr⁻¹ (Hosokawa et al. 2013; Schleicher et al. 2013). This seed mass is high relative to some other simulations, such as IllustrisTNG (Nelson et al. 2019). Dynamical BH estimates (e.g., Nguyen et al. 2019) also indicate that BHs in dwarf galaxies can fall to masses well below 10⁶M_⊙. However, seeding BHs at a higher mass ensures that both BH dynamics and gas accretion are well resolved throughout their lifetimes (Tremmel et al. 2015). Further, the formation mechanism ensures that BHs are only seeded in regions with dense, collapsing gas that is unlikely to form stars, and hence are more likely to grow BHs rapidly. For more detailed caveats related to the BH seeding mechanism, see Section 3.5.

Dynamical friction between BHs and their hosts causes BHs to sink toward the galaxy center (Kazantzidis et al. 2005; Pfister et al. 2017). Dynamical friction interactions occur at both large scales (Colpi et al. 1999) and scales below the resolution of most cosmological simulations. Common practice in cosmological simulations is to reposition the BH along local potential gradients as it begins to migrate, forcefully recentering it. However, this method suppresses BH motion around the galaxy and artificially inflates BH growth rates (Tremmel et al. 2017). Romulus25 instead incorporates a dynamical friction subgrid recipe shown to reproduce realistic BH sinking timescales. This recipe allows BHs to naturally "wander" within galaxies (Tremmel et al. 2018; Bellovary et al. 2019; Reines et al. 2020). Tremmel et al. (2015) find that the spatial resolution in Romulus25, combined with its oversampling of dark matter, also avoids unrealistic numerical heating of BHs found in lower-resolution simulations.

Assuming an isotropic velocity distribution of particles within a gravitational softening length, ${\epsilon }_{g}$, from the BH, we can write the acceleration due to dynamical friction (Chandrasekhar 1943):

$$\begin{eqnarray}&&{a}_{\mathrm{DF}}=-4\pi {G}^{2}{M}_{\mathrm{BH}}\rho (\lt {v}_{\mathrm{BH}})\mathrm{ln}{\rm{\Lambda }}\displaystyle \frac{{v}_{\mathrm{BH}}}{{v}_{\mathrm{BH}}^{3}},\end{eqnarray} \tag{ 2 }$$

where M_BH is the mass of the BH, ρ(<v_BH) is the density of nearby particles moving slower than the BH, $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithm, and v_BH is the velocity of the BH. Velocities are calculated relative to the local center of mass within the smoothing kernel. The Coulomb logarithm is approximated by $\mathrm{ln}{\rm{\Lambda }}\sim \mathrm{ln}\tfrac{{b}_{\max }}{{b}_{\min }}$, where b_max and b_min are, respectively, the maximum and minimum impact parameters of the surrounding particles. The maximum impact parameter is restricted to ${b}_{\max }={\epsilon }_{g}$ to avoid double counting of resolved dynamical frictional forces. The minimum impact parameter is restricted to the 90° deflection radius with a lower limit of the Schwarzschild radius of the BH. Acceleration is calculated from the nearest 64 collisionless particles. Mergers occur when BHs fall within 2_g of one another and have low enough relative velocity to be considered gravitationally bound.

BH accretion is handled through a Bondi–Hoyle prescription modified to incorporate angular momentum support on resolved scales. The accretion rate is driven by mass flux across a resolved accretion radius, defined as the radius where gravitational potential balances with the minimally resolved thermal energy of the surrounding gas. Accretion rates are averaged over the smallest simulation time step at a given time, typically 10⁴−10⁵ yr. We can write the accretion rate depending on whether the dominant gas motion is rotational, v_θ, or bulk flow, v_bulk:

$$\begin{eqnarray}\dot{M}=\alpha \times \left\{\begin{array}{lc}\displaystyle \frac{\pi {G}^{2}{M}_{\mathrm{BH}}^{2}\rho }{{\left({v}_{\mathrm{bulk}}^{2}+{c}_{s}^{2}\right)}^{3/2}} & \mathrm{if}\ {v}_{\mathrm{bulk}}\gt {v}_{\theta }\\ \displaystyle \frac{\pi {G}^{2}{M}_{\mathrm{BH}}^{2}\rho {c}_{s}}{{\left({v}_{\theta }^{2}+{c}_{s}^{2}\right)}^{2}} & \mathrm{if}\ {v}_{\mathrm{bulk}}\lt {v}_{\theta },\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray*}\alpha =\left\{\begin{array}{lc}{\left(\displaystyle \frac{{n}_{\mathrm{gas}}}{{n}_{* }}\right)}^{\beta } & \mathrm{if}\ {n}_{\mathrm{gas}}\geqslant {n}_{* }\\ 1 & \mathrm{if}\ {n}_{\mathrm{gas}}\leqslant {n}_{* },\end{array}\right.\end{eqnarray*}$$

is the density-dependent boost factor that corrects for underestimates of accretion rate due to resolution limitations (Booth & Schaye 2009), β = 2 is the corresponding boost coefficient, n_gas is the number density of the surrounding gas, n_* is the star formation density threshold, ρ is the mass density of the surrounding gas, c_s is the local sound speed, v_θ is the rotational velocity of the surrounding gas at the smallest resolved scales, and v_bulk is the bulk velocity of the surrounding gas. This calculation is performed over the 32 nearest particles. BH accretion in Romulus25 is Eddington limited.

BH feedback is handled through a subgrid recipe similar to the blast wave SN feedback model. Thermal energy from accretion is isotropically transferred to the nearest 32 gas particles, weighted by the SPH kernel. To ensure realistic dissipation of feedback energy, gas particles that receive energy are stopped from cooling for roughly the gas dynamical time step over which the accretion is calculated (Tremmel et al. 2017). The energy coupled to the surrounding gas is

$$\begin{eqnarray}&&E={\epsilon }_{r}{\epsilon }_{f}\dot{M}{c}^{2}{dt},\end{eqnarray} \tag{ 4 }$$

where _r = 0.1 is the assumed radiative efficiency of the BH and _f = 0.02 is the coupling efficiency of the thermal energy to the surrounding gas (see below). This form of BH feedback has been shown to efficiently quench galaxies with halo masses above a few × 10¹¹M_⊙ (Pontzen et al. 2017).

The accretion efficiency, β, and coupling efficiency, _f, were constrained through 48 zoom-in simulations of the same four galaxies used for constraining star formation parameters, run with BH physics. The star formation parameters were left unchanged from the initial parameter search, and instead the BH parameters were allowed to change. The parameter space was explored using the Kriging algorithm and graded by agreement with an empirical z = 0 BH mass–stellar mass relation (Schramm & Silverman 2013). The results of this parameter search were also used in the high-resolution cosmological hydrodynamic, galaxy cluster simulation, RomulusC (Tremmel et al. 2019).

Figure 1 shows the z = 0.05 M_BH−M_star relation for all isolated Romulus25 galaxies with central BHs. Galaxies below M_star < 10¹⁰M_⊙ exhibit a high degree of scatter in M_BH, which is a phenomenon that has been predicted and observed before. Using semianalytic models, Volonteri & Natarajan (2009) find that the evolution of BHs toward the local M_BH−σ relation is dependent on both the BH growth history and seed BH mass. They find that low-mass BHs exhibit a higher amount of scatter on the M_BH−σ relation than higher-mass BHs. Using cosmological hydrodynamical simulations, Barai & de Gouveia Dal Pino (2019) similarly find high scatter in the M_BH−M_star relation at low stellar masses. Reines & Volonteri (2015) find that the M_BH−M_star relation observed in high-mass galaxies breaks down in low-mass galaxies. They find that low-mass star-forming galaxies may instead follow a relation with lower BH masses than expected. In high-mass galaxies, Shankar et al. (2016) find that dynamical estimates of the local M_BH−M_star relation are biased by the requirement that the BH sphere of influence be fully resolved. Shankar et al. find that corrections to this bias place dynamical M_BH−M_star relations in closer agreement with AGN-derived relations and eliminate much of the perceived scatter in M_BH at high M_star.

Zoom In Zoom Out Reset image size

Galaxies around M_star ∼ 10⁸M_⊙ tend to clump up at the BH seed mass and twice the seed mass. The simulation BH seeding mechanism introduces a floor that likely inflates BH masses in the lightest dwarf galaxies relative to observations. We do not remove hosts of seed mass BHs from the sample since it is unclear what additional biases would be introduced relative to a higher-resolution seeding mechanism. For some analyses that follow we separate the sample into two mass bins since resolution may impact results at the lowest masses.

We compare with observed relations from Schramm & Silverman (2013; SS13) and Greene et al. (2019). Romulus25 agrees with the Greene et al. relation for early-type galaxies at stellar masses M_star > 10¹⁰M_⊙, as well as with the SS13 relation. Above M_star ≳ 10¹⁰M_⊙, galaxies follow the SS13 relation in part because the BH physics is tuned to match with the low-mass end of the SS13 relation. For similar reasons, Romulus25 does not agree with the Greene et al. relation for late-type galaxies.

The divergence between the various relations may be due to a mixture of a few effects: morphological differences in the host galaxy, differences in accretion efficiencies between BHs, and uncertainties in observational BH mass estimators. Reines & Volonteri (2015) find significant differences between their M_BH−M_star relations for bulge-dominated galaxies and AGNs. Their bulge-dominated galaxy sample has dynamical BH mass estimates, while their AGN sample contains a mixture of pseudobulges and classical bulges with broad Hα virial + reverberation-mapped BH masses. They find that the scatter in the combined M_BH−M_star relation is closely related to the morphology of the host galaxy, where bulge-dominated galaxies sit on a relation with similar slope but higher normalization than spiral galaxies. SS13 do not make morphological distinctions in their sample and base their relation on X-ray-selected AGNs with broad Mg II virial BH masses. Similarly, Davis et al. (2018) find that late-type galaxies follow a steeper M_BH−M_star relation than early-type galaxies down to M_star ∼ 10^10.5M_⊙, which Sahu et al. (2019) find is fundamentally a difference in bulge morphology of the host galaxies.

Trump et al. (2011) find that broad emission line regions are only observed in rapidly accreting AGNs with luminosities greater than 10⁻² L_Edd, where L_Edd is the Eddington luminosity. By analyzing Romulus25, Ricarte et al. (2019) find that high Eddington ratios are associated with BHs that exhibit systematically lower masses than expected for their host stellar mass.

Indirect measurements of BH masses within AGNs are typically thought to yield large uncertainties, dependent on their luminosity and redshift (e.g., Shen & Kelly 2012; Kelly & Shen 2013). Virial BH mass estimates of AGNs require assumptions of the geometry and orientation of the broad-line region (Denney et al. 2010; Barth et al. 2011) and may be biased by nongravitational contributions to broad-line widths (Krolik 2001). However, Reines & Volonteri (2015) argue that the difference between M_BH−M_star relations for AGNs and elliptical galaxies is far larger than can be explained by uncertainties in virial BH masses. They further find that their virial BH masses would need an unreasonably high virial factor in order to match the elliptical galaxy relation.

We quantify scatter in the M_BH−M_star relation by the residual away from the median M_BH in bins of M_star. Within each bin we classify BHs with masses in the bottom 25% as "undermassive," while BHs within the top 25% are classified as "overmassive." BHs falling between the two quartiles are classified as "median." We summarize our sample sizes in Table 1. We show this classification scheme at all stellar masses, though the classification is most meaningful in dwarf galaxies where the central BHs show a high degree of scatter in M_BH. A similar classification scheme is built by Li et al. (2020) for M_star > 10¹⁰M_⊙ Illustris galaxies and M_star > 10⁹M_⊙ TNG100 galaxies, though they instead define overmassive and undermassive BHs on the M_BH−σ relation.

It should be emphasized that current observations find that BHs in dwarf galaxies can have masses lower than the resolution limit of Romulus25. Baldassare et al. (2015) find a BH in the dwarf galaxy RGG 118 with a virial mass of M_BH = 5 × 10⁴M_⊙. Nguyen et al. (2017) calculate a dynamical BH mass of M_BH = 1.5 × 10⁵M_⊙ in the nearby dwarf galaxy NGC 404. Nguyen et al. (2019) improve dynamical BH masses for three nearby low-mass galaxies, where all three show BH masses below 1 million solar masses and one shows a dynamical mass of M_BH = 6.8 × 10³M_⊙. Graham & Soria (2019) and Graham et al. (2019) use BH scaling relations to predict BH masses as low as M_BH ∼ 10⁴M_⊙ in low-mass galaxies in the Virgo Cluster.

Following Ricarte et al. (2019), we can partially compensate for the effects of the BH seeding mechanism and define the total mass a given BH has grown via accretion, M_acc. This definition completely excludes the contributions of BH seeding and only counts accretion onto every progenitor within the BH merger tree. We are able to trace M_acc well below the BH seed mass because of the high resolution of gas accretion in Romulus25. Figure 2 shows the M_acc−M_star relation for all isolated Romulus25 galaxies with central BHs, compared with observed M_BH−M_star relations. We find that the M_acc−M_star relation continues linearly below the BH seed mass and shows better agreement with the observed relations. Galaxies below M_star < 10¹⁰M_⊙ still show a high degree of scatter in M_acc, while higher mass galaxies do not. Overmassive BHs in galaxies below M_star < 10¹⁰M_⊙ tend to fall above the observed relations. Notably, BHs that are considered overmassive or undermassive in M_BH are typically overmassive or undermassive in M_acc as well. Although it is unclear how much the BH seed mass affects growth, Romulus25 is capable of producing many BHs with growth consistent with current observational constraints by z = 0.05.

Zoom In Zoom Out Reset image size

Understanding the source of scatter in the M_BH−M_star relation requires first understanding how the BHs evolved to the present day. BHs grow through BH–BH mergers, as well as through accretion of gas particles. We trace the growth history of each BH and calculate the growth through mergers onto the main progenitor, M_mergers.

Figure 3 shows the fraction of total BH mass grown via mergers by z = 0.05 versus the host stellar mass. By z = 0.05, only 8% of undermassive BHs have ever undergone a BH merger, and instead the vast majority grow solely through accretion. On the other hand, 95% of overmassive BHs have undergone at least one BH merger and grow through a mixture of BH mergers and accretion at higher rates than undermassive BHs. The primary mode of growth for overmassive BHs depends on host stellar mass, where those found in hosts below M_star ≲ 10⁹M_⊙ grow primarily through BH mergers and those found in hosts above ${M}_{\mathrm{star}}\gtrsim {10}^{9}{M}_{\odot }$ grow primarily through accretion onto the main BH progenitor.

Zoom In Zoom Out Reset image size

We next turn to how BHs evolve relative to the stellar mass of their hosts. Figure 4 shows the evolution of undermassive, overmassive, and median BHs and their hosts onto both the z = 0.05 M_BH−M_star relation and M_acc−M_star relation. The black line indicates the SS13 locally observed relation, where the dashed portion indicates a linear extrapolation. Regardless of how we frame BH growth, we find fundamental differences in growth histories between undermassive and overmassive BH hosts. Undermassive BHs tend to evolve onto the z = 0.05 relation by building up M_BH only after the host has built up its stars. On the other hand, overmassive BH hosts tend to either grow their BHs before their stellar mass or grow both in tandem. The differences in M_BH−M_star evolution strongly suggest that overmassive and undermassive BH hosts may build up their stars and dark matter in fundamentally different ways.

Zoom In Zoom Out Reset image size

With the growth histories of the BHs in hand, we now study the environments in which they formed and reside. We trace the structural evolution of stars and dark matter in the BH hosts across cosmic time. In the following analysis we consider the 197 isolated dwarf galaxies in Romulus25 that do not host any BHs alongside the 205 that do host BHs, in order to better understand the role of BHs in structural evolution.

Munshi et al. (2013) identify a systematic overestimate in dark-matter-only (DMO) simulation halo masses when compared to simulations that include baryon physics and outflows. When comparing with results from DMO simulations, we adjust halo masses such that M_200,sim = 0.8M_200,DMO for halos with masses M₂₀₀ = 10⁸−10¹²M_⊙.

3.3.1. Stellar Mass–Halo Mass Relation

Figure 5 shows the z = 0.05 stellar mass–halo mass (SMHM) relation for all isolated Romulus25 galaxies, with masses adjusted with corrections from Munshi et al. (2013). Points are colored by whether the hosted BH is overmassive, undermassive, or median, or if the galaxy does not host a BH. A similar figure of the SMHM relation in Romulus25 for all central halos can be found in Tremmel et al. (2017). We mark the M_star = 10¹⁰M_⊙ dwarf galaxy boundary and limit the axes to focus on low-mass galaxies. We include abundance matching estimates from Moster et al. (2013) and Kravtsov et al. (2018) for reference. Above M₂₀₀ > 10¹¹M_⊙, overmassive BHs tend to be found in halos with lower stellar masses than expected for their halo mass. Undermassive BHs instead tend to be found in halos with higher stellar masses than expected for their halo mass. Undermassive and median BH hosts in particular tend to sit along or above abundance matching estimates of the SMHM relation. Galaxies without BHs follow a similar relation to undermassive BH hosts, and above M₂₀₀ > 10¹¹M_⊙ they sit at higher M_star than overmassive BH hosts at a given halo mass.

Zoom In Zoom Out Reset image size

To better quantify the connection between scatter in M_acc and scatter in M_star, we define the residual quantities ${\rm{\Delta }}\mathrm{log}$ M_acc and ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{star}}$. We fit a smoothing spline to the median SMHM relation for all isolated Romulus25 galaxies and then calculate ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{star}}$ as the residual from the median $\mathrm{log}{M}_{\mathrm{star}}$ for a given halo mass. We similarly fit a spline to the median M_BH−M_star relation and use the previous fit to find the expected M_acc for a given halo mass. We then calculate ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{acc}}$ as the residual from the median. Figure 6 shows ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{acc}}$ versus ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{star}}$, split between galaxies with stellar mass 10⁸M_⊙ < M_star < 10⁹M_⊙ and 10⁹M_⊙ <M_star < 10¹⁰M_⊙. We split our sample in this way to better isolate resolution effects at the lowest masses. We distinguish between hosts of overmassive, undermassive, and median BHs. Positive/negative ${\rm{\Delta }}\mathrm{log}{M}_{\mathrm{acc}}$ roughly correspond with overmassive/undermassive BHs.

Zoom In Zoom Out Reset image size

We find that overmassive BHs tend to be found in halos with fewer stars than expected from the median, while undermassive BHs are found in halos with more stars than expected. For galaxies with 10⁹M_⊙ < M_star < 10¹⁰M_⊙, overmassive BHs tend to be found at lower ${\rm{\Delta }}\mathrm{log}$ M_star and higher ${\rm{\Delta }}\mathrm{log}$ M_acc than their undermassive counterparts. Galaxies with 10⁸M_⊙ < M_star < 10⁹M_⊙ show little difference in ${\rm{\Delta }}\mathrm{log}$ M_star between overmassive and undermassive BHs. This result implies that BH accretion, and hence feedback, may suppress star formation in isolated galaxies between 10⁹M_⊙ < M_star < 10¹⁰M_⊙.

3.3.2. Structural Evolution

We now turn to the impact of BHs on the structural evolution of dark matter and stars. We find that the undermassive/overmassive nature of a BH is directly tied to the structural evolution of its host. Figure 7 shows the evolution of M₂₀₀ of the main halo progenitor, scaled by M₂₀₀ at z = 0.05, distinguishing between undermassive, overmassive, and median BH hosts. We also include a comparison between galaxies hosting overmassive BHs and galaxies without BHs. We distinguish between galaxies with z = 0.05 stellar mass 10⁸M_⊙ < M_star < 10⁹M_⊙ and 10⁹M_⊙ < M_star <10¹⁰M_⊙. Overmassive BH hosts tend to build up their halos earlier than both hosts of undermassive BHs and galaxies without BHs. The delay in halo growth is most pronounced in dwarf galaxies above M_star > 10⁹M_⊙ at z = 0.05.

Zoom In Zoom Out Reset image size

Figure 8 similarly shows the median cumulative star formation history for each class of BH hosts. Overmassive BH hosts build up their stellar mass much more rapidly than both hosts of undermassive BHs and galaxies without BHs. As with halo mass, the differences in stellar mass evolution are most pronounced in dwarf galaxies above M_star > 10⁹M_⊙ at z = 0.05.

Zoom In Zoom Out Reset image size

Further analysis shows that the formation of overmassive/undermassive BHs can likely be attributed to differences in halo assembly times. We fit each halo density profile with a flexible five-parameter model (Dekel et al. 2017):

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{c}}{{x}^{\alpha }{\left(1+{x}^{1/\beta }\right)}^{\beta (\gamma -\alpha )}},\end{eqnarray} \tag{ 5 }$$

where α and γ are, respectively, the inner and outer profile slopes, β is the intermediate shape parameter, and the radius is scaled by the scale radius, R_s, which is related to the halo concentration, c,

$$\begin{eqnarray}&&x=\displaystyle \frac{r}{{R}_{s}},\quad c=\displaystyle \frac{{R}_{200}}{{R}_{s}}.\end{eqnarray} \tag{ 6 }$$

Figure 9 shows the evolution of halo concentration over time. Overmassive BH hosts have higher halo concentrations than their undermassive counterparts (see Macciò et al. 2008; Zhao et al. 2009, and references therein), as well as halos without BHs. A two-sample, two-sided Kolmogorov–Smirnov (K-S) test on overmassive versus undermassive BH host concentrations yields a K-S statistic D_n,m = 0.29 and allows us to reject the null hypothesis at the 0.02 level that the halos are drawn from the same distribution of concentrations. Similarly, a K-S test on overmassive BH hosts versus halos without BHs yields a K-S statistic D_n,m = 0.45 and allows us to reject the null hypothesis at the 8 × 10⁻⁸ level.

Zoom In Zoom Out Reset image size

The buildup of stars prior to BH growth, the lack of BH mergers in undermassive BH hosts, and the lower halo concentrations all suggest that undermassive BHs initially formed in environments with a lower abundance of cold gas than is necessary to seed multiple BHs. In contrast, overmassive BHs were likely initially seeded in environments with an abundance of cold gas, where BHs had a higher likelihood to merge and accrete efficiently.

There is strong evidence that the central regions of massive galaxies are most affected by the presence of a central BH (Cheung et al. 2012; Barro et al. 2017; Choi et al. 2018). We trace the evolution of stars within the central regions of the host galaxy and search for a connection to the central BH.

We define the stellar mass surface density within the stellar half-light radius, r_e:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{e}=\displaystyle \frac{{M}_{\mathrm{star}}(\lt {r}_{e})}{\pi {r}_{e}^{2}},\end{eqnarray} \tag{ 7 }$$

where we calculate r_e by fitting a Sérsic profile to the projected face-on V-band surface brightness profile, with a surface brightness cutoff of 32 mag arcsec⁻² (Abraham & van Dokkum 2014).

Cheung et al. (2012) find that stellar density within the central 1 kpc, Σ₁, robustly follows quenching in local galaxies above M_star ≳ 10⁸M_⊙. Chen et al. (2019) build a schematic model that finds that galaxies begin to rapidly quench once they evolve over a boundary in Σ₁−M_star space. Many Romulus25 dwarf galaxies fall below r_e < 1 kpc at early times; hence, we use r_e to consistently define a central region across time. Both Franx et al. (2008) and Barro et al. (2013) find a strong relationship between Σ_e and both the star formation rate (SFR) and total stellar mass out to z ∼ 3.5. Hence, tracing the evolution of Σ_e while distinguishing between hosts of overmassive and undermassive BHs can give insight into the effects of BH growth on central star formation.

Figure 10 shows the evolution of Σ_e across time for galaxies with undermassive, overmassive, and median BHs, as well as for galaxies without BHs. We distinguish between galaxies with z = 0.05 stellar mass 10⁸M_⊙ < M_star < 10⁹M_⊙ and 10⁹M_⊙ < M_star < 10¹⁰M_⊙. Similar to the buildup of total stellar mass, overmassive BH hosts build up their central stellar density more rapidly than both undermassive BH hosts and galaxies without BHs. However, above M_star > 10⁹M_⊙, overmassive BH hosts show suppression of Σ_e growth starting at redshift z ∼ 2. By z ∼ 0.5, many hosts of undermassive BHs and galaxies without BHs reach 0.5 dex higher Σ_e than hosts of overmassive BHs. Hosts of undermassive BHs and galaxies without BHs may instead flatten in Σ_e at late times, around z ∼ 0.1. We confirm that these results would remain qualitatively the same were we to use surface densities calculated within the central 1 kpc rather than r_e.

Zoom In Zoom Out Reset image size

In short, overmassive BHs were first seeded in early-forming halos, building up their dark matter and stars rapidly in tandem with growth of the BH. Despite having higher halo concentrations, overmassive BH hosts have similar or lower central stellar density than their undermassive counterparts by late times, indicating a measure of star formation suppression within the central regions. Undermassive BH hosts instead follow nearly identical evolutionary tracks to galaxies without BHs, growing dark matter and stars later and exhibiting high central stellar densities at late times. Dickey et al. (2019) find a similar relationship in isolated galaxies with stellar mass 10⁹M_⊙ < M_star < 10^9.5M_⊙, where strong signatures of AGNs correlate with an older stellar population in the host galaxy. Li et al. (2020) similarly find that overmassive BH hosts form earlier and have lower present-day SFRs.

Choi et al. (2018) find similar evolution of Σ_e in zoom-in simulations of galaxies with M_star > 10^10.9M_⊙ run with and without AGN feedback. They find that galaxies with AGN feedback build up Σ_e until z ∼ 2, after which Σ_e turns over and begins to decrease. The stellar cores become diffuse primarily through AGN-driven stellar mass loss and gas mass loss "puffing up" the central region. They find that the turnover in Σ_e is concurrent with quenching of star formation. Galaxies run without AGN feedback indefinitely increase their central stellar densities and do not experience the same level of quenching. Both Guo et al. (2013) and Barro et al. (2017) similarly observe that stellar cores diffuse over time in CANDELS GOODS-S galaxies with stellar masses 10⁹M_⊙ < M_star < 10^11.5M_⊙.

Stagnation in central stellar density can come about in a few ways. Stellar evolution can eject mass from stars and return it to gas in the ISM (van Dokkum et al. 2014). Mass loss through stellar outflows directly competes with new star formation promoted by the increased gas mass (Kennicutt et al. 1994). As a result, Choi et al. (2018) find this effect to contribute little to central stellar density suppression.

Compaction events may occur through strongly dissipational, gas-rich mergers driving stars and gas toward the galaxy center. Conversely, gas-poor mergers can reduce the core stellar density by driving rapid size growth with little growth in mass (Hopkins et al. 2008; Nipoti et al. 2009; Covington et al. 2011; Oser et al. 2012; Hilz et al. 2013; Porter et al. 2014). Overmassive BH hosts in Romulus25 are often found to have up to 4 orders of magnitude lower gas fractions relative to galaxies without BHs (see Section 3.4.2) and hence may be subject to stellar diffusion through gas-poor mergers, though we find no evidence of major mergers driving rapid changes in stellar or gas profiles of overmassive or median BH hosts.

There are other physical processes that have been commonly found to reduce central stellar density but are unresolved in Romulus25. Binary BH systems may be capable of clearing out galaxy cores (Milosavljević & Merritt 2001; Kormendy et al. 2009), but BH scouring occurs on scales below the spatial resolution limit of Romulus25  (Rantala et al. 2017, 2018). Mass loss and heating driven by outflows and SNe may reduce the gravitational potential and in turn allow for outward stellar migration (Fan et al. 2008; Choi et al. 2018), but Romulus25 does not resolve these effects on the central potential.

Finally, heating and mass loss driven by outflows and SNe can in turn drive gas outward, restrict gas cooling, and hence suppress star formation (Somerville & Davé 2015). Below we show that the differences in structural evolution between Romulus25 galaxies with overmassive BHs and other galaxies are likely due to this BH feedback.

3.4.1. Energy Injection by BHs

We find that median and overmassive BHs are capable of injecting more energy than SNe into the surrounding gas of their hosts. We calculate total energy injection via BHs and SNe by first integrating their energy outputs across cosmic time. A set fraction of the energy output is injected into the surrounding ISM (see Section 2.1). Figure 11 shows the ratio of energy injected by BHs to SNe versus the host stellar mass. Many overmassive and median BHs are capable of injecting substantially more energy than SNe. On the other hand, undermassive BH hosts are all dominated by energy injection by SNe. At all stellar masses, overmassive BHs hosts are injected with two to three more combined BH + SN energy than undermassive BH hosts.

Zoom In Zoom Out Reset image size

A higher E_BH/E_SN suggests that outflows from BHs may more efficiently heat and drive gas than SNe outflows. It is important to note that while BHs produce and inject copious amounts of energy, it is ultimately the feedback prescription that determines the effect on the host. Feedback models that inject kinetic energy typically drive winds at higher velocities than those that inject purely thermal energy (Choi et al. 2018).

3.4.2. Cold Gas Depletion

Next, we turn to the relationship between BHs and the amount of cold gas in the host galaxy. Figure 12 shows the H i gas mass versus stellar mass relation for isolated dwarf galaxies. We compare with observations from Bradford et al. (2018) and Catinella et al. (2018). Bradford et al. combine a new set of 21 cm observations with the H i-selected ALFALFA survey (Haynes et al. 2011) to analyze gas depletion in local galaxies (0.002 < z < 0.055) with stellar mass 10⁷M_⊙ < M_star < 10^9.5M_⊙. Catinella et al. measure H i content of local (0.01 < z < 0.05) stellar-mass-selected xGASS galaxies, with stellar masses 10⁹M_⊙ < M_star < 10^11.5M_⊙.

Zoom In Zoom Out Reset image size

Hosts of undermassive BHs and galaxies without BHs follow the Catinella et al. relation at the high-mass end and are consistent with nondepleted galaxies from Bradford et al. across all stellar masses. Hosts of undermassive BHs and galaxies without BHs show little indication of significant cold gas depletion. On the other hand, many overmassive and median BH hosts exhibit lower M_{H i} than expected for their stellar mass, indicating a high degree of cold gas depletion. Although ALFALFA is most sensitive to H i above ${M}_{{\rm{H}}{\rm{I}}}\gtrsim {10}^{7}{M}_{\odot }$, it is worth noting that the extreme levels of cold gas depletion seen in overmassive BH hosts are often orders of magnitude below what is seen in either observational comparison sample. Bradford et al. similarly find that isolated galaxies with stellar mass 10^9.2M_⊙ < M_star < 10^9.5M_⊙ with strong signatures of AGNs tend to show a higher degree of gas depletion than similar galaxies with weaker AGN signatures, though they find that this effect does not extend to more massive galaxies. Bradford et al. do not rule out the ejection and heating of cold gas by unusually bursty and compact star formation activity.

3.4.3. Star Formation Quenching

Tremmel et al. (2019) fit the Romulus25 star formation main sequence by calculating median values of the SFR within 0.1 dex bins of stellar mass between 10⁸M_⊙ < M_star < 10¹⁰M_⊙, for galaxies considered relatively isolated by the same criteria we use in this work. They find a best fit of $\mathrm{log}(\mathrm{SFR})\,=1.206\times \mathrm{log}({M}_{\mathrm{star}})-11.7$ at z = 0. We define galaxies whose SFR falls a factor of 10 below the main sequence to have quenched star formation. We calculate SFRs averaged over the previous 250 Myr.

Figure 13 shows the z = 0.05 star formation main sequence for Romulus25 dwarf galaxies. We distinguish between hosts of overmassive, undermassive, and median BHs, as well as isolated galaxies without BHs. Quenched galaxies are marked with filled points. We mark the Tremmel et al. (2019) relation with a solid line.

Zoom In Zoom Out Reset image size

Quenching tends to occur in isolated dwarf galaxies that host overmassive or median BHs. We find 12 quenched dwarf galaxies, 7 of which host a BH. Above ${M}_{\mathrm{star}}\gtrsim {10}^{8.6}{M}_{\odot }$, quenching only occurs in dwarf galaxies that host BHs. Quenched galaxies that host a BH tend to host overmassive BHs, regardless of host stellar mass.

To help identify the source of quenching as internal or external, we track the interaction history of each galaxy of interest and estimate the tidal effects of close encounters. Following Karachentsev & Makarov (1999), we calculate tidal indices, Θ, such that

$$\begin{eqnarray}&&{\rm{\Theta }}=\max \{\mathrm{log}({M}_{k}/{D}_{k}^{3})\}-11.75,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&k=1,2,...,N,\end{eqnarray} \tag{ 9 }$$

where M and D are the masses and 3D separations, respectively, of the kth closest halo. Negative values of Θ indicate no tidal disturbance of the main halo by nearby halos, while high positive values (we arbitrarily choose Θ > 5) indicate significant tidal disturbance.

Analyzing the encounter history of our quenched galaxies, there is one undermassive BH host that is particularly close to the quenched boundary, and further analysis reveals that it was stripped of H i gas and quenched immediately following tidal disturbance by a more massive halo. Similarly, the two quenched median BH hosts both show evidence of such encounters followed by H i depletion and quenching. Of the five quenched overmassive BHs, two exhibit high tidal indices from encounters with massive halos. Thus, quenching can occur in galaxies below M_star < 10^8.6M_⊙ even if they do not host a BH, though all but one show clear evidence of past encounters with a more massive galaxy and a subsequent drop in star formation.

As discussed in Tremmel et al. (2019), our definition of quenched is different from some observations, such as Wetzel et al. (2012), who adopt a flat threshold in specific SFR of 10⁻¹¹ yr⁻¹. Regardless, our results change little when we instead use a flat specific SFR threshold of 10⁻¹¹ yr⁻¹.

Although Romulus25 uses a purely thermal feedback model, simulations have found success in using feedback models that incorporate both thermal feedback and mechanical feedback that efficiently drives high-velocity winds. Choi et al. (2015) find that the inclusion of mechanical feedback more efficiently suppresses late-time star formation and produces AGN luminosities in line with observations. Choi et al. (2016) find that mixed thermal and mechanical AGN feedback models yield reasonable results for ex situ and in situ star formation and realistic gas and stellar structural properties. Weinberger et al. (2017) find that dual-mode AGN feedback for weakly/highly accreting BHs yields realistic star formation suppression, gas fractions, BH growth, and thermodynamic profiles.

Our findings are in line with results from Ricarte et al. (2019), who find that isolated Romulus25 galaxies with M_star > 10^9.5M_⊙ show signs of coevolution with their central BH. Ricarte et al. find that the BH accretion rate follows the SFR in star-forming galaxies, regardless of redshift, stellar mass, or large-scale environment. Further, they find that such BHs grow in tandem with their host galaxies, eventually being confined to a line of constant M_BH/M_star. They suggest that self-regulation of BH growth through feedback is a possible driver of coevolution seen in isolated Romulus25 galaxies.

Our analysis of the role of BHs in dwarf galaxy evolution has a few caveats. Although we find correlations between BH properties and properties of the host galaxy, the precise effect of BH activity on dwarf galaxy evolution is unclear. We have not directly traced the effects of BHs on the surrounding environment (e.g., tracing outflows, tracking heating, turning off BH physics, and rerunning the simulation), and hence we cannot say for certain how BHs may drive changes in their hosts within Romulus25. Further, observations from Mezcua et al. (2019) find that dwarf galaxies can host powerful jets whose mechanical feedback may strongly impact the host star formation history, an effect that is not accounted for in Romulus25. We find no evidence of starbursts in compact galaxies (Diamond-Stanic et al. 2012) as the dominant mechanism driving trends in star formation suppression, though the spatial resolution of Romulus25 likely restricts especially compact galaxies from forming. It appears that BH feedback plays a larger role than often thought within dwarf galaxies (Martín-Navarro & Mezcua 2018), and it is likely that both stellar feedback and BH feedback together drive suppression of central stellar density and overall star formation in Romulus25 dwarf galaxies. Future work will further analyze the role of AGNs in the evolution of dwarf galaxies, in particular how BH activity relates to suppression of star formation.

As seen in Equation (3), accretion onto BHs is sensitive to the BH mass. This is particularly important for two reasons: (1) some BHs in Romulus25 may unphysically merge immediately after seeding, and (2) the seed mass is likely too high in the lowest-mass galaxies. Some BHs effectively form at higher masses than the seed mass owing to seeding of multiple, clustered BHs and subsequent rapid merging. Bellovary et al. (2019) find a similar phenomenon in zoom-in simulations of dwarf galaxies. Bellovary et al. define "overmerging" to occur if either of the merging BHs were seeded <100 Myr prior to the merger event. They suggest that this time frame is long enough for BH feedback from existing BHs to suppress future BH formation. We find that approximately 35% of overmassive BHs and 15% of median BHs have experienced an overmerging event. Since BHs that undergo BH mergers have a subsequently higher accretion rate, such overmerging may unphysically contribute to the BH mass. However, we find that overmerging does not guarantee that a BH will become overmassive or grow to high M_BH.

Finally, current observations of BHs in dwarf galaxies indicate that BHs may have masses lower than the BH seed mass used in this work (Baldassare et al. 2015; Reines & Volonteri 2015; Nguyen et al. 2017, 2018, 2019; Schutte et al. 2019). Observed BH masses in M_star ∼ 10⁹ galaxies are typically a few × 10⁵M_⊙ but can reach as low as 6.8 ×10³M_⊙ (Nguyen et al. 2019). Mezcua et al. (2018) combine Chandra data of z ≲ 2.4 dwarf galaxies with the M_BH−M_star relation from Reines & Volonteri (2015) and calculate BH masses consistent with the undermassive but not the overmassive BH masses in this work. In particular, overmassive BHs in Romulus25 may be unrealistically massive and accrete too much, causing us to overestimate their ability to drive galaxy-scale changes through feedback. Moving to higher-resolution simulations in the future may give a more clear understanding of both how BHs grow within dwarf galaxies and how dwarf galaxies may have their star formation suppressed by the BH.