Environmental Dependence of Self-regulating Black Hole Feedback in Massive Galaxies

The simulation domain is a (4 Mpc)³ box with a 64³ root grid and up to 10 levels of refinement. The central (128 kpc)³ region enforces static regions of grid refinement ranging from level 6 to level 9, with the refinement level increasing toward the center. In the central (2 kpc)³, the mesh is fixed to be at the highest level of refinement. This design ensures that the CGM is highly resolved at all times with a minimum cell size Δl ≈ 61 pc and a maximum cell size Δl < 1 kpc.

The gravitational potentials in our simulations have three components, and they do not change with time. We use an NFW form (Navarro et al. 1997) for the dark matter halo, with mass density ${\rho }_{\mathrm{DM}}\propto {r}^{-1}{(1+r/{r}_{{\rm{s}}})}^{-2}$, where r_s is the NFW scale radius and c₂₀₀ = r₂₀₀/r_s is a halo concentration parameter. For the central galaxy, we use a Hernquist profile (Hernquist 1990), with a stellar mass density ${\rho }_{* }\propto {r}^{-1}{(1+r/{r}_{{\rm{H}}})}^{-3}$, where r_H is the Hernquist scale radius. The potential of the central supermassive black hole of mass M_BH follows a Paczynski–Witta form (Paczyńsky & Wiita 1980).

Table 1 provides the parameter values for each model. Two of those models are based on particular galaxies. The multiphase galaxy (MPG) model is intended to resemble NGC 5044, which has a central stellar velocity dispersion⁵ ${\sigma }_{v}\approx 225\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and represents the population of massive elliptical galaxies with extended multiphase gas. The single-phase galaxy (SPG) model is intended to resemble NGC 4472, which has a central stellar velocity dispersion ${\sigma }_{v}\approx 282\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and represents the population of massive elliptical galaxies without much multiphase gas beyond the central kiloparsec. The brightest cluster galaxy (BCG) model is meant to represent the central region of a typical massive galaxy cluster. The smaller elliptical galaxy (SEG) model is designed to test how the AGN feedback mechanism used in the larger halos operates when the halo mass and central stellar velocity dispersion are reduced.

Table 1.  Galactic Environmental Parameters

Galaxy Model	M200	c200	r200	M*	${r}_{{\rm{H}}}$	σv	MBH	Analog Galaxy
	(M⊙)		(kpc)	(M⊙)	(kpc)	(km s−1)	(108M⊙)
BGC	8 × 1014	4.9	1920	3 × 1012	10	230	50	...
MPG	4.4 × 1013	9.5	730	1.2 × 1011	1.2	230	4.6	NGC 5044
SPG	4.0 × 1013	7.5	700	2 × 1011	1.6	280	26	NGC 4472
SEG	2 × 1012	5	275	1 × 1011	1.5	150	0.7	⋯

Download table as:  ASCIITypeset image

A simulation cell forms a star particle if its gas satisfies several criteria based on the Cen & Ostriker (1992) prescription:

• 1.  
  The baryon density must exceed a threshold density ($\sim 1\,{\mathrm{cm}}^{-3}$).
• 2.  
  The flow must be converging (${\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{b}\lt 0$).
• 3.  
  The gas must be cold (T < 1.1 × 10⁴ K).
• 4.  
  The gas mass of the cell (m_b) must exceed 10³ M_⊙.
• 5.  
  The cooling time of the gas must be less than the dynamical time for that cell's gas, ${t}_{\mathrm{dyn}}=\sqrt{3\pi /(32G{\rho }_{\mathrm{gas}})}$.

A star particle is then formed with mass ${m}_{* }={f}_{* ,\mathrm{eff}}({\rm{\Delta }}t/{t}_{\mathrm{dyn}}){m}_{{\rm{b}}}$, where the star formation efficiency parameter is set to f_*,eff = 0.1 in this simulation suite.

The central galaxy's stars heat the gas through two separate channels:

• 1.  
  New stars forming during the course of the simulation produce Type II SNe (SNe II) that impart both thermal energy and momentum.
• 2.  
  Old stars with the density distribution of the Hernquist potential (Section 3.2) add heat through SN Ia explosions and thermalization of stellar kinetic energy.

Feedback from stars formed during the simulation follows the prescription from Bryan et al. (2014). After a star particle of mass m_* forms, a fraction f_m,* of its mass is added back to the cell, along with thermal energy E_SN = f_SNm_*c². Our simulations adopt the parameter values ${f}_{\mathrm{SN}}=1\times {10}^{-5}$ and f_m,* = 0.25. The returned gas formally has a metallicity f_Z,* = 0.02, which can be used as a passive tracer but is not included in our radiative cooling calculations. This process starts immediately after the formation of the star particle and decays exponentially, with a time constant of 1 Myr.

SN Ia heating is modeled with steady, spherically symmetric injection of thermal energy into the simulation domain at a rate proportional to the stellar mass density. The total energy ejected from SN Ia explosions assumes 10⁵¹ erg per SN Ia at a specific rate of $3\times {10}^{-14}\,\mathrm{SNIa}\,{\mathrm{yr}}^{-1}\,{M}_{\odot }^{-1}$, following Voit et al. (2015c). At this rate, an old stellar population of mass ∼10¹¹ M_⊙ adds $\sim {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ of thermal energy. The old stellar population also injects kinetic energy as it sheds gas mass in the form of stellar winds and SN Ia explosions. Our simulations assume that this kinetic energy immediately thermalizes. We assume a specific gas ejection rate ${\alpha }_{* }={10}^{-19}\,{{\rm{s}}}^{-1}$, such that the net ejected matter per unit time per unit volume is ${\alpha }_{* }{\rho }_{* }$. To simplify the calculation of thermalized kinetic energy, we assume an isotropic 1D stellar velocity dispersion of 300 km s⁻¹ at all radii in all of our runs, following Wang et al. (2019). The difference between this uniform value of ${\sigma }_{v}$ and the actual one is inconsequential, because energy input from SN Ia heating is several times (≳5) greater than the kinetic energy injection in all cases.

AGN feedback is introduced into the simulation using a feedback zone attached to the AGN particle (Meece et al. 2017), which is always located at the geometric center of the halo. We drive AGN feedback in the form of a bipolar outflow by putting source terms for mass, momentum, and energy into the fluid equations as follows:

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{u}})={S}_{\rho }\\ \displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}\otimes {\boldsymbol{v}})=-{\rm{\nabla }}P-\rho {\rm{\nabla }}{\rm{\Phi }}+{S}_{{\boldsymbol{p}}}\\ \displaystyle \frac{\partial (\rho e)}{\partial t}+{\rm{\nabla }}\cdot [(\rho e+P){\boldsymbol{v}}]=\rho {\boldsymbol{v}}\cdot {\boldsymbol{g}}-{n}_{e}{n}_{i}{\rm{\Lambda }}(T,Z)+{S}_{e},\end{array}\end{eqnarray*}$$

where ${S}_{\rho }$, ${S}_{{\boldsymbol{p}}}$, and S_e are the density, momentum, and energy source terms, respectively. The specific energy e includes kinetic energy, and the corresponding equation of state is $P=(\gamma -1)\rho (e-{v}^{2}/2)$.

3.5.1. Accretion and AGN Feedback Efficiency

The accretion rate ${\dot{M}}_{\mathrm{acc}}$ onto the central supermassive black hole is calculated by assuming that all the cold gas ($T\lt {10}^{5}$ K) within $r\lt 0.5\,\mathrm{kpc}$ accretes onto the central black hole on a 1 Myr timescale. This mass accretion rate fuels AGN feedback at an energy output rate given by

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{AGN}}={\epsilon }_{\mathrm{AGN}}{\dot{M}}_{\mathrm{acc}}{c}^{2},\end{eqnarray} \tag{ 1 }$$

where c is the speed of light and the feedback efficiency parameter ${\epsilon }_{\mathrm{AGN}}$ is taken to be 10⁻⁴ for all our runs with AGN feedback. A cold gas mass equal to ${\dot{M}}_{\mathrm{acc}}{\rm{\Delta }}t$ is removed from the spherical accretion zone ($r\lt 0.5$ kpc) by subtraction of gas mass from each cell in that zone with a temperature below 10⁵ K, and the amount of gas mass subtracted from the cell is proportional to its total gas mass.

3.5.2. Feedback Energy Deposition

The AGN output energy is partitioned into kinetic and thermal parts and introduced using the source terms described above. The source regions are cylinders of radius 0.5 kpc extending along the jet axis from $r=0.5\,\mathrm{kpc}$ to r = 1 kpc in each direction. Each jet therefore subtends 1 rad at r = 1 kpc.

In each cell of volume ${\rm{\Delta }}{V}_{\mathrm{cell}}$ within the source region of volume ${V}_{\mathrm{jet}}$, the density source term is ${S}_{\rho }={\rm{\Delta }}{m}_{\mathrm{cell}}/{\rm{\Delta }}{V}_{\mathrm{cell}}$, where ${\rm{\Delta }}{m}_{\mathrm{cell}}=({\rm{\Delta }}{V}_{\mathrm{cell}}/{V}_{\mathrm{jet}}){\dot{M}}_{\mathrm{acc}}{\rm{\Delta }}t$ and ${\rm{\Delta }}t$ is the time step. Within the source region, the corresponding energy source term in each cell is ${S}_{e}=({\rm{\Delta }}{e}_{\mathrm{cell}}+{\rm{\Delta }}{\mathrm{KE}}_{\mathrm{cell}})/{V}_{\mathrm{cell}}$, where

$$\begin{eqnarray}&&{\rm{\Delta }}{e}_{\mathrm{cell}}=(1-{f}_{\mathrm{kin}}){\dot{E}}_{\mathrm{AGN}}{\rm{\Delta }}t\cdot \displaystyle \frac{{\rm{\Delta }}{V}_{\mathrm{cell}}}{{V}_{\mathrm{jet}}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{KE}}_{\mathrm{cell}}={f}_{\mathrm{kin}}{\dot{E}}_{\mathrm{AGN}}{\rm{\Delta }}t\cdot \displaystyle \frac{{\rm{\Delta }}{V}_{\mathrm{cell}}}{{V}_{\mathrm{jet}}}.\end{eqnarray} \tag{ 3 }$$

For all our runs the ratio of kinetic to total AGN output energy is fixed at ${f}_{\mathrm{kin}}=0.9$.

The momentum source term in each cell is ${S}_{{\boldsymbol{p}}}={\rm{\Delta }}{{\boldsymbol{p}}}_{\mathrm{cell}}/{\rm{\Delta }}{V}_{\mathrm{cell}}$, where

$$\begin{eqnarray}&&{\rm{\Delta }}{{\boldsymbol{p}}}_{\mathrm{cell}}={\rm{\Delta }}{m}_{\mathrm{cell}}\sqrt{\displaystyle \frac{2{\rm{\Delta }}\mathrm{KE}}{{\rm{\Delta }}{m}_{\mathrm{cell}}}}\hat{{\boldsymbol{n}}}\end{eqnarray} \tag{ 4 }$$

and $\hat{{\boldsymbol{n}}}$ is a unit vector pointing away from the origin along the jet injection axis.

To see how quickly star formation begins and what its rate would be without AGN feedback, we ran simulations of the three smaller halos with ${\epsilon }_{\mathrm{AGN}}=0$. We did not perform a similar simulation of the BCG, because stellar feedback is obviously insufficient to limit star formation in such a massive halo. Figure 2 shows the resulting star formation rates as functions of time. Unsurprisingly, the MPG model begins to form stars almost immediately, because radiative cooling exceeds SN Ia heating at all radii. More than ${10}^{9}\,{M}_{\odot }$ of cold gas ($T\lt {10}^{5}$ K) accumulates by $t\approx 250\,\mathrm{Myr}$, and star formation then proceeds at a steady rate $\sim 25\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. According to Voit (2011), the steady cooling flow rate associated with an entropy profile $K(r)\propto r$ is

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{8\pi }{3}\mu {m}_{p}{({kT})}^{2}{\rm{\Lambda }}(T){\left(\displaystyle \frac{K}{r}\right)}^{-3}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{rcl} & \approx & 24\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\ {\left(\displaystyle \frac{{kT}}{1\mathrm{keV}}\right)}^{2}\ \left[\displaystyle \frac{{\rm{\Lambda }}(T)}{{10}^{-23}\,\mathrm{erg}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}}\right]\\ & & \times {\left(\displaystyle \frac{K/r}{1\mathrm{keV}{\mathrm{cm}}^{2}{\mathrm{kpc}}^{-1}}\right)}^{-3}\end{array}\end{eqnarray} \tag{ 6 }$$

in an isothermal potential. The asymptotic star formation rate observed in the MPG simulation without AGN feedback is therefore consistent with its initial entropy profile, which is $K/r\approx 1\,\mathrm{keV}\,{\mathrm{cm}}^{2}\,{\mathrm{kpc}}^{-1}$ at ∼10 kpc.

Zoom In Zoom Out Reset image size

More surprisingly, star formation also begins promptly in the SEG simulation without AGN feedback, even though SN Ia heating initially exceeds radiative cooling out to ∼5 kpc from the center. Figure 2 shows that star formation rises to $\sim 8\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ within the first 200 Myr. The amount of cold gas that accumulates during that same time period is comparable to the steady-state amount in the MPG simulation and is similar to the amount of hot gas that starts the simulation with a cooling time ${t}_{\mathrm{cool}}\lesssim 200\,\mathrm{Myr}$. The asymptotic star formation rate is again consistent with Equation (6), but with ${kT}\sim 0.3\,\mathrm{keV}$.

Cooling and star formation begin within ∼100 Myr in the SEG model despite the central SN Ia heating, because the weight of its CGM prevents the gas ejected by old stars from leaving the galaxy. The initial gas-mass density at ∼1 kpc is $\rho \approx 5\times {10}^{-26}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. At similar radii, the stellar mass density is ${\rho }_{* }\approx {10}^{-22}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, meaning that stellar ejecta can double the gas-mass density there on a timescale $\rho {({\alpha }_{* }{\rho }_{* })}^{-1}\sim 150\,\mathrm{Myr}$ if the gas cannot be pushed outward. Some of the gas is pushed outward, but not enough to prevent a buildup of gas within the central 5 kpc. Meanwhile, radiative cooling of gas just beyond ∼5 kpc produces an entropy inversion that makes the atmosphere convectively unstable and promotes thermal instability. Condensing gas clouds then sink to the center and initiate star formation. The resulting SN II explosions briefly suppress additional star formation but result in uplift of low-entropy ambient gas that precipitates at ∼10 kpc, forming new cold clouds. Those clouds then rain down toward the galaxy's center, and during the next 50 Myr the ambient gas settles into a steady cooling flow.

Star formation requires more time to reach a steady state in the SPG simulation without AGN feedback. A brief burst of star formation happens at $t\sim 50\,\mathrm{Myr}$, because radiative cooling initially exceeds SN Ia heating everywhere. SN II feedback from that initial burst then lowers the central density, which allows SN Ia heating to exceed radiative cooling out to ∼1 kpc. Star formation remains suppressed for the next ∼1.3 Gyr, while the central gas pressure gradually rises. The central pressure goes up during this period because SN Ia heating cannot push ejected stellar gas outward as fast as it accumulates and also because cooling of the overlying layers increases their weight. Eventually, the central gas density becomes great enough for radiative cooling to exceed SN Ia heating, and the resulting cooling flow boosts the star formation to a steady-state rate of ∼3.5 ${M}_{\odot }$ yr⁻¹ at $t\sim 2.0\,\mathrm{Gyr}$.

In all of our simulations with AGN feedback, star formation is highly suppressed relative to the respective no-AGN counterparts (see Figure 2). However, condensation of the ambient medium couples with AGN feedback differently, depending on both the depth of the central potential well and the atmospheric pressure at larger radii. In the SPG simulation, coupling between condensation and AGN feedback is remarkably tight and maintains a nearly steady feedback-regulated state. In contrast, the MPG and BCG simulations exhibit greater feedback bursts. This section examines how these three massive halos self-regulate, while Section 5.3 looks at what happens in the SEG simulation, which fails to self-regulate.

5.2.1. Radial Profiles

Figure 3 shows the median emissivity-weighted radial entropy profiles in the SPG, MPG, and BCG simulations with AGN feedback. A dashed magenta line indicates the initial entropy profile in each simulation. A solid red line traces the median profile for the entire period from 0.5 to 1.5 Gyr. Dark-gray shading shows the 20th–80th percentile range of the median entropy profile during that period, and light-gray shading shows the 1st–99th percentile range.

Zoom In Zoom Out Reset image size

In each case, the median entropy profile shifts from its initial state into a different self-regulated state. The BCG entropy profile settles from a flat-entropy core into a self-regulated state with $K\propto {r}^{2/3}$ in the central ∼30 kpc. The MPG entropy profile rises to a self-regulated state with a mean entropy at $\lt 10\,\mathrm{kpc}$ several times greater than the initial state. The SPG entropy profile also rises within the central ∼10 kpc by a factor of ∼2. However, the shading shows that tightly coupled feedback confines the SPG entropy profile to a much narrower range than in the BCG or MPG.

Comparing the data in Figure 3 with the simulation results reveals a significant discrepancy in the vicinity of 1 kpc, where specific entropy in the SPG and MPG simulations exceeds the observations by a factor of 2–3. In the middle panel showing the MPG simulation, some of that discrepancy might arise from a selection effect. The galaxies shown have particularly bright X-ray emission produced by atmospheres denser than average for their mass, meaning that they have lower-than-average specific entropy, but the data still remain within the fluctuation range of the simulation. However, the data in the SPG panel are well outside the fluctuation range of the SPG simulation. We therefore suspect that the entropy excess near 1 kpc results from a limitation of our simulation, which is the thickness of the AGN jet there. Section 6 discusses that limitation in more detail.

Figure 4 shows the median radial profiles of both ${t}_{\mathrm{cool}}$ (top panels) and the ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio (bottom panels) during the same time period for the same three simulations. The left panels show that the SPG simulation self-regulates with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\gt 20$ at $\gt 1\,\mathrm{kpc}$ and ${t}_{\mathrm{cool}}\gt 1\,\mathrm{Gyr}$ at 10 kpc. The middle panels show that the MPG self-regulates with most of its time spent in the $10\lt \min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\lt 20$ range, with ${t}_{\mathrm{cool}}\lt 1\,\mathrm{Gyr}$ at 10 kpc. The right panels show that the BCG self-regulates with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})$ fluctuating in and out of that range, also with ${t}_{\mathrm{cool}}\lt 1\,\mathrm{Gyr}$ at 10 kpc.

Zoom In Zoom Out Reset image size

5.2.2. Self-regulation and SN Ia Heating

Self-regulation in these three simulations depends on how multiphase condensation couples AGN feedback with the state of the ambient medium. The top panels of Figure 5 show how injected jet power (${P}_{\mathrm{jet}}$) is related to X-ray luminosity (L_X) over the time period 0–1.5 Gyr. The X-ray luminosity is calculated in the "yt" package (Turk et al. 2011) using the "Cloudy" emission table for temperatures in the range of 0.5–7 keV with solar metallicity for all runs except the galaxy-cluster case, where a one-third solar metallicity emission table is used. After the initial AGN outburst of $\sim {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, the SPG settles into a nearly steady state with time-averaged AGN jet power that is similar to the X-ray luminosity from within the central 30 kpc. Frequent bursts fueled by fluctuations in the amount of cold gas within the central 0.5 kpc cause jet power to vary by a factor of ∼10, but the power output remains steady when smoothed over timescales $\gt 100\,\mathrm{Myr}$. However, the other two simulations experience much greater fluctuations in jet power.

Zoom In Zoom Out Reset image size

In the MPG, AGN feedback is bimodal. The simulation starts with an outburst of jet power $\sim {10}^{44}$ erg s⁻¹. Heat input from that outburst causes the galaxy's atmosphere to expand, lowering its density and significantly reducing radiative cooling of the central 30 kpc. The AGN then enters a low-power state with ${P}_{\mathrm{jet}}$ fluctuating on a ∼100 Myr timescale between $\lt {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and a few times ${10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Meanwhile, the atmosphere's X-ray luminosity climbs, because time-averaged AGN power is much less than L_X from within 30 kpc. Those radiative losses allow the weight of the CGM to compress the galactic atmosphere, gradually raising its density and pressure. The AGN remains in this low-power mode for $\approx 800\,\mathrm{Myr}$ but then reverts back to a high-power state, similar to the initial one, for another $\approx 200\,\mathrm{Myr}$.

The BCG simulation remains in a state similar to the high-power mode of the MPG simulation most of the time and does not have a low-power mode. It is either near ${P}_{\mathrm{jet}}\sim {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ or at ${P}_{\mathrm{jet}}\lt {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. The state of extremely low power is likely to be artificial, resulting from the fact that our feedback algorithm sets AGN power to zero if there is no cold gas within the central 0.5 kpc. A more realistic model would include AGN power resulting from Bondi-like accretion of hot ambient gas, but jet power in that mode would be far too low to significantly affect the surrounding atmosphere.

The bottom panels of Figure 5 show how the mode of AGN feedback is related to the ratio of SN Ia heating to radiative cooling in the central few kiloparsecs. The radiative cooling for each radial shell is calculated $\sim {n}_{e}{n}_{i}{\rm{\Lambda }}(T){dV}$, where dV is the volume of each radial shell and ${\rm{\Lambda }}(T)$ is the tabulated Sutherland–Dopita cooling function (Sutherland & Dopita 1993). Initially, radiative cooling in the SPG is slightly greater than SN Ia heating within ∼5 kpc of the center. Cooling of that gas fuels a ∼100 Myr burst of feedback that lowers the central gas density until SN Ia heating exceeds radiative cooling from ∼0.5 to ∼5 kpc. AGN feedback then enters the low-power mode, fueled only by cooling of gas within ∼0.5 kpc of the center. And that mode is sufficient to keep the SPG in a steady state, for at least ∼1.5 Gyr.

In the MPG, a larger initial burst of AGN power is needed to lower the atmosphere's density because its confining CGM pressure is greater. However, SN Ia heating becomes comparable to radiative cooling at ∼1 kpc by $t\approx 300\,\mathrm{Myr}$. The simulation then settles into the low-power feedback mode for nearly 1 Gyr. During that time, AGN power is less steady than in the SPG because near equality of SN Ia heating and radiative cooling at $r\lt 3\,\mathrm{kpc}$ allows larger condensation events to intermittently feed the AGN.

The BCG simulation, on the other hand, is in a cooling-dominated state everywhere during virtually the entire time period. AGN feedback cannot lower the atmosphere's central density enough for SN Ia heating to equal radiative cooling. Therefore, the mode of self-regulation connects AGN feedback to large condensation events, which occur well outside of the central kiloparsec.

5.2.3. Cold Gas and Star Formation

In these simulations, central accumulation of cold gas couples atmospheric conditions with AGN feedback, while accumulations of cold gas at larger radii facilitate star formation. Figure 6 shows how the masses of cold gas (${M}_{\mathrm{cold}}$) and new stars change with time. In the SPG, the accumulations of cold gas are always small (10⁴–${10}^{6}\,{M}_{\odot }$), and star formation is negligible. Note that the feedback algorithm described in Section 3.5.1 produces an AGN feedback power

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{AGN}}=6\times {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\ \left[\displaystyle \frac{{M}_{\mathrm{cold}}(\lt 0.5\,\mathrm{kpc})}{{10}^{6}\,{M}_{\odot }}\right],\end{eqnarray} \tag{ 7 }$$

given ${\epsilon }_{\mathrm{AGN}}={10}^{-4}$. The fluctuations in feedback power shown for this galaxy in Figure 5 are therefore consistent with the fluctuations in cold gas mass shown in Figure 6, as long as a large proportion of that gas ends up accreting onto the central black hole. Cold gas clouds forming through condensation are therefore consumed before they can form stars, linking the precipitation rate within 0.5 kpc directly to AGN feedback.

Zoom In Zoom Out Reset image size

In the MPG, the mass of cold gas is $\sim {10}^{7}\,{M}_{\odot }$ in the high-power mode and $\sim {10}^{5}$–${10}^{6}\,{M}_{\odot }$ in the low-power mode. Those amounts of cold gas are also largely consistent with the fluctuations in feedback power shown in Figure 5. Furthermore, the periods when cold gas extends beyond 1 kpc correlate with periods of greater star formation and AGN power. Feedback events that cause multiphase precipitation at larger radii therefore promote star formation in our simulations, because the cold gas clouds have time to form stars before sinking into the AGN accretion region in the central 0.5 kpc. Figure 2 shows that those star formation events briefly peak at a rate of $\sim 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, but the accumulated stellar mass in Figure 6 implies a time-averaged rate of $\sim 1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and a specific star formation rate of $\sim {10}^{-12}\,{\mathrm{yr}}^{-1}$.

The BCG experiences several condensation events that push ${M}_{\mathrm{cold}}$ up to $\sim {10}^{8}\,{M}_{\odot }$, producing surges of AGN power exceeding ${10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. Figure 5 shows that surges of this magnitude with a significant duty cycle are necessary to compensate for the radiative losses from within 30 kpc. Note that increasing the parameter ${\epsilon }_{\mathrm{AGN}}$ and accretion time in our feedback algorithm would result in greater accumulation of cold gas, because balance between AGN feedback and radiative cooling would require less cold gas to be reprocessed within the accretion zone. The amount of star formation in this simulation is therefore contingent on the choice of ${\epsilon }_{\mathrm{AGN}}$. For the choice ${\epsilon }_{\mathrm{AGN}}={10}^{-4}$, the time-averaged star formation rate is $\sim 0.4\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$.

Figure 7 shows the maximum radial extent of cold gas ($T\lt {10}^{5}$ K) as these different systems evolve. In the SPG simulation, cold gas remains concentrated within 1 kpc and usually within 0.5 kpc, except during the initial outburst. In contrast, cold gas in the MPG generally extends beyond 1 kpc (sometimes as far as ∼10 kpc), but cold gas in the BCG simulation tends to be less extended. Comparing cold gas radial extent to ${P}_{\mathrm{jet}}$ in Figure 5 shows that cold gas becomes most extended following periods of strong AGN feedback, indicating that uplift of central gas promotes condensation at greater altitudes (Revaz et al. 2008; McNamara et al. 2016; Voit et al. 2017).

Zoom In Zoom Out Reset image size

Our AGN feedback simulation in a lower-mass halo (the SEG with ${M}_{200}=2\times {10}^{12}{M}_{\odot }$) dramatically differs from the others. Figure 8 shows that the simulation produces a large AGN outburst and some star formation during the first 0.4 Gyr and then enters a state in which AGN feedback does not compensate for radiative losses. Instead, star formation and AGN power shut down, while ${L}_{{\rm{X}}}$ and ${M}_{\mathrm{cold}}$ steadily decline with time.

Zoom In Zoom Out Reset image size

As in the SEG simulation without AGN feedback (see Section 5.1), the action begins at ∼50 Myr when the central gas starts to condense. Those first cold clouds then trigger a self-exciting AGN feedback outburst. Uplift of low-entropy ambient gas simulates multiphase condensation that causes $\sim 5\times {10}^{8}\,{M}_{\odot }$ of cold gas to precipitate by t = 200 Myr. Much of that cold gas falls radially back into the accretion zone (<0.5 kpc), boosting the jet power by a few times—up to $\sim {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, more than an order of magnitude greater than the radiative losses from the ambient medium (${L}_{{\rm{X}}}$). This powerful feedback event blows out much of the hot gas atmosphere but does not destroy the cold gas clouds, which can continue to rain back down in the accretion zone. The result is a decaying AGN feedback mode, in which intermittent accretion events produce smaller feedback outbursts that gradually lower the X-ray luminosity.

Figure 9 shows the evolution of entropy (left panel), ${t}_{\mathrm{cool}}$ (middle panel), and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio (right panel) during the first 600 Myr of the SEG simulation with AGN feedback. During the first 300 Myr, strong AGN feedback pushes some of the low-entropy ambient gas from the central region to beyond ∼10 kpc and heats the rest. The radiatively cooling outflow becomes most unstable to condensation near 10 kpc, where ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ remains near unity for tens of megayears. That low-entropy gas then sinks inward and becomes the cold gas that sustains runaway feedback. A similar runaway does not happen in the more massive halos, because AGN feedback in those systems does not produce as much uplift and convection. Later in the simulation, we find ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\lt 10$ at ∼100 kpc, but that gas condenses slowly because the cooling time there is several gigayears.

Zoom In Zoom Out Reset image size

The SPG and MPG simulations presented here were designed, in part, to test the "black hole feedback valve" mechanism proposed by Voit et al. (2015c, 2020) and summarized in Section 2. Qualitatively, the SPG and MPG simulations with AGN feedback do indeed self-regulate as envisioned, with AGN feedback tuning itself so that local radiative cooling is similar to SN Ia heating out to distances several kiloparsecs from the galaxy's center. Figure 5 shows that the SPG simulation (${\sigma }_{v}\approx 280\,\mathrm{km}\,{{\rm{s}}}^{-1}$) begins with cooling exceeding heating everywhere and settles into a steady state with SN Ia heating exceeding radiative cooling from ∼1 to ∼5 kpc, as predicted by the feedback valve model for galaxies with ${\sigma }_{v}\gt 240\,\mathrm{km}\,{{\rm{s}}}^{-1}$. AGN feedback in this mode is fueled by cooling of gas within the central 0.5 kpc, while SN Ia heating sweeps much of the gas released by stars at larger radii out of the galaxy. This state can remain steady as long as AGN feedback prevents the confining CGM pressure from building up, and it succeeds for at least 1.5 Gyr because the time-averaged jet power roughly matches radiative losses from the inner 30 kpc. However, we have not yet tested whether a galaxy with the SPG potential but a higher-pressure CGM (like the initial state of the MPG) tunes itself to this same steady state.

The MPG simulation (${\sigma }_{v}\approx 230\,\mathrm{km}\,{{\rm{s}}}^{-1}$) also begins with cooling exceeding SN Ia heating everywhere. Kinetic AGN feedback with power $\sim {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ then lowers the atmospheric density and abates when radiative cooling becomes similar to SN Ia heating within ∼3 kpc. The galaxy remains in this low-power state for nearly 1 Gyr but cannot sustain it because radiative losses from the inner 30 kpc exceed the jet power. The CGM pressure there remains high and gradually increases, preventing SN heating from sweeping ejected stellar gas out of the galaxy. Gas density and radiative cooling within the galaxy therefore both rise until the gas density reaches a ceiling imposed by the condition $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ (see Figure 4), at which precipitation of cold clouds inevitably triggers a large increase in feedback power. This second burst of feedback again lowers the central gas density until SN Ia heating is comparable to local radiative cooling (see the yellow line in the bottom middle panel of Figure 5). Consequently, the configuration of the ambient medium fluctuates but is bracketed by the state in which SN Ia heating exceeds radiative cooling and the state with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$, in alignment with the black hole feedback valve model for galaxies with ${\sigma }_{v}\approx 230\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Quantitatively, however, the simulations do not match the black hole feedback valve model in detail. One of the model's key predictions is that the power-law slope of the entropy profile should exceed $K\propto {r}^{2/3}$ at $r\sim 1$–10 kpc in galaxies with ${\sigma }_{v}\gt 240\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The entropy profile at 1–10 kpc in our SPG simulation is much flatter than this prediction, which assumes that SN Ia heating exceeds AGN heating within the galaxy. While it is difficult to measure with precision how much of the heating at ∼1 kpc is resulting from thermalization of jet energy, the excess entropy at that radius in the SPG simulation, relative to both the data and the analytical steady-flow models of Voit et al. (2020), suggests that the discrepancy results from excessive thermalization of jet kinetic energy at ∼1 kpc in the simulation (see Section 6).

The predictions made in Voit et al. (2020) for galaxies like the SEG and BCG are less specific, but the results of our simulations of the SEG and BCG with AGN feedback generally conform to the model's expectations. According to the model, a galaxy that does not supply enough feedback power to lift its CGM and alleviate the confining pressure should remain in a precipitation-limited state that self-regulates through multiphase condensation. The BCG simulation is consistent with that expectation of the model. For the SEG, the model predicts that multiphase circulation should be inevitable, because feedback overturns the atmosphere's entropy gradient. And indeed, the feedback events observed in the simulation lift low-entropy ambient gas out of the center, catalyzing widespread condensation and production of cold gas, much of which falls back toward the center.

Several earlier numerical studies have explored the role of kinetic AGN feedback fueled by cold gas accretion in massive elliptical galaxies like the ones simulated in this paper. The efforts most closely related to our SPG and MPG simulations were published by Gaspari et al. (2011a, 2012a) and Wang et al. (2019). The ones most closely related to our BCG simulation were published by Gaspari et al. (2011b, 2012b), Li et al. (2015), Prasad et al. (2018), and Meece et al. (2017).

Gaspari et al. (2011a, 2011b, 2012a) performed the first suite of simulations to demonstrate that bipolar jets fueled by cold accretion can tune themselves to balance radiative cooling without overheating the central gas. Collectively, those three papers explored a range of halo mass similar to the range spanned by our SPG, MPG, and BCG models, but with substantially lower spatial resolution. Also, they adjusted their AGN feedback efficiency parameter, equivalent to our ${\epsilon }_{\mathrm{AGN}}$, to optimize agreement with observations, finding the best results for ${\epsilon }_{\mathrm{AGN}}\lesssim 3\times {10}^{-4}$ in lower-mass halos and ${\epsilon }_{\mathrm{AGN}}\gtrsim 5\times {10}^{-3}$ in cluster-scale halos. While mass and energy input from the old stellar population were included in these simulations, the role of the old stellar population in the overall feedback loop was not specifically analyzed.

The simulations of Wang et al. (2019), like ours, were motivated by the analysis of Voit et al. (2015c) and focused on distinguishing the roles of the central gravitational potential and the stellar mass and energy sources. Wang et al. (2019) performed two simulations similar to our SPG and MPG simulations with AGN feedback. The initial conditions in those simulations were not identical to ours but were inspired by the same two galaxies, with NGC 4472 representing single-phase elliptical galaxies and NGC 5044 representing MPGs. In alignment with our simulation results, Wang et al. (2019) found that AGN feedback in the galaxy similar to NGC 4472 maintained a relatively steady hot gas atmosphere with small amounts of centrally concentrated cold gas, while the same AGN feedback algorithm in the galaxy similar to NGC 5044 caused greater fluctuations in the hot gas atmosphere and produced larger quantities of extended cold gas. Their general findings therefore also support the black hole feedback valve model.

However, the details of our simulation results differ from those of Wang et al. (2019). First, the median AGN power in Wang et al. (2019) is much greater, with jet power often rising above ${10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the SPG and rarely dropping below ${10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the MPG. Their large AGN feedback efficiency (equivalent to ${\epsilon }_{\mathrm{AGN}}=5\times {10}^{-3}$) allows the same amount of cold gas accretion to produce much more power, but that cannot be the whole explanation for the power difference, because self-regulation over a 1.5 Gyr period requires time-integrated heat input within the central 10–30 kpc (with ${t}_{\mathrm{cool}}\lesssim 1.5$ Gyr) to balance radiative losses from the same region. Therefore, kinetic AGN power in Wang et al. (2019) must be thermalizing over a larger region, implying that it is propagating farther from the center. Second, specific entropy near ∼1 kpc in the SPG of Wang et al. (2019) remains below $10\,\mathrm{keV}\,{\mathrm{cm}}^{2}$ most of the time and is typically $\approx 5\,\mathrm{keV}\,{\mathrm{cm}}^{2}$, in better agreement with observations of SPGs than our SPG simulation. In Section 6, we hypothesized that the excess entropy at ∼1 kpc in our SPG simulation resulted from jets that were insufficiently narrow. And indeed, the jets implemented by Wang et al. (2019) are narrower, having a transverse momentum profile $\propto \exp (-{r}^{2}/2{r}_{\mathrm{jet}}^{2})$ with ${r}_{\mathrm{jet}}=183$ pc. With greater power and a smaller cross section, the jets in Wang et al. (2019) have a much greater momentum flux than ours and are capable of propagating to much greater distances, also accounting for why the simulations of Wang et al. (2019) require more AGN power to self-regulate.

Self-regulation of AGN feedback in our BCG simulation is broadly similar to what is observed in other simulations of its type (e.g., Gaspari et al. 2011b, 2012b; Li & Bryan 2014; Li et al. 2015; Prasad et al. 2015, 2018; Meece et al. 2017). Our BCG simulation's typical value of $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})$ is greater than most, with a median ratio of ∼25. In a future paper, we will show that the self-regulated K(r) profile of our BCG simulation, which has an inner slope $K\propto {r}^{2/3}$ (see Figure 3), is in excellent agreement with the observations of Hogan et al. (2017) and Babyk et al. (2018). However, unlike some of the other simulations, it produces less cold gas than is observed in cluster cores, and the cold gas it does produce rarely extends beyond 3 kpc.

The main reason for the lack of cold gas in our BCG simulation is the low feedback efficiency parameter and the small accretion time we have chosen. The small accretion time causes the cold gas forming within $r\lt 0.5\,\mathrm{kpc}$ to be quickly removed from the simulation domain. This results in a high accretion rate, producing powerful AGN jets before much cold gas accumulates in the cluster core, despite the low AGN feedback efficiency. For ${\epsilon }_{\mathrm{AGN}}={10}^{-4}$, the cold gas accretion rate required to sustain ${10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ of feedback power is $18\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Our algorithm converts all of that cold gas to hot gas and expels it from the central region in a jet. Over the course of 1 Gyr, more than ${10}^{10}\,{M}_{\odot }$ would otherwise have accumulated, and much of it would likely have formed stars at a rate of $\sim 10\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. In some of the other cluster-scale simulations (e.g., Gaspari et al. 2012b; Li & Bryan 2014; Prasad et al. 2015), much of the cold gas persists indefinitely in a torus orbiting outside of the accretion zone because of the stochastic angular momentum it gains during kinetic feedback bursts. Our simulation does not produce such a torus, and we will analyze what inhibits torus formation in a future paper.