#  Hyenas: X-ray Bubbles and Cavities in the Intra-Group Medium

Fred Jennings,1 Arif Babul,2,3,4 Romeel Dave,1,5,6 Weiguang Cui, 1,7,8 and
Douglas Rennehan 9
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh, EH9 3HJ, United Kingdom
2Leverhulme Visiting Prof., Institute for Astronomy, University of Edinburgh,
Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, United Kingdom
3Department of Physics & Astronomy, University of Victoria, BC V8X 4M6, Canada
4Infosys Visiting Chair Professor, Indian Institute of Science, Bangalore
560012, India
5University of the Western Cape, Bellville, Cape Town 7535, South Africa
6South African Astronomical Observatories, Observatory, Cape Town 7925, South
Africa
7Departamento de Física Teórica, Universidad Autónoma de Madrid, Módulo 15,
E-28049 Madrid, Spain
8Centro de Investigación Avanzada en Física Fundamental (CIAFF), Facultad de
Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain
9Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New
York, NY 10010, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the role of the Simba feedback model on the structure of the
Intra-Group Medium (IGrM) in the new Hyenas suite of cutting-edge cosmological
zoom-in simulations. Using 34 high-resolution zooms of halos spanning from
$10^{13}-10^{14}$ $M_{\odot}$ at $z=0.286$, we follow halos for 700 Myr, over
several major active galactic nuclei (AGN) jet feedback events. We use the
MOXHA package to generate mock Chandra X-ray observations, as well as
predictive mocks for the upcoming LEM mission, identifying many feedback-
generated features such as cavities, shock-fronts, and hot-spots, closely
mimicking real observations. Our sample comprises $105$ snapshots with
identified cavities, $50$ with single bubbles and $55$ with two, and spans
three orders of magnitude in observed cavity enthalpies, from
$10^{41}-10^{44}$ erg/s. Comparing semi-major axis length, midpoint radius,
and eccentricity to a matched sample of observations, we find good agreement
in cavity dimensions with real catalogues. We estimate cavity power from our
mock maps following observational procedures, showing that this is typically
more than enough to offset halo cooling, particularly in low-mass halos, where
we match the observed excess in energy relative to cooling. Bubble enthalpy as
measured with the usual midpoint pressure typically exceeds the energy
released by the most recent jet event, hinting that the mechanical work is
done predominantly at a lower pressure against the IGrM. We demonstrate for
the first time that X-ray cavities are observable in a modern large-scale
simulation suite and discuss the use of realistic cavity mock observations as
new halo-scale constraints on feedback models in cosmological simulations.

###### keywords:

galaxies:groups:general-galaxies:clusters:general–X-rays:galaxies:clusters.

††pubyear: 2015††pagerange: Hyenas: X-ray Bubbles and Cavities in the Intra-
Group Medium–Hyenas: X-ray Bubbles and Cavities in the Intra-Group Medium

##  1 Introduction

The physical process(es) driving the quenching of massive galaxies remains one
of the key mysteries in understanding the formation of today’s galaxy
population. Currently, models that are successful in creating a quenched
galaxy population all invoke a process referred to as “radio mode" or
“maintenance mode" feedback (Somerville & Davé, 2015) from active galactic
nuclei (AGN), which involves relativistic jets emanating from massive central
black holes that heat the surrounding circum-galactic medium and counteracts
the expected cooling of hot gas (Babul et al., 2002; Brüggen et al., 2005;
Roychowdhury et al., 2005; Sijacki & Springel, 2006; Nusser et al., 2006).
However, the exact mechanism by which jets transfer energy to the surrounding
gas to prevent cooling remains uncertain (Nusser et al., 2006; Fabian, 2012;
Babul et al., 2013; Cielo et al., 2018). Jets are seen to be highly
collimated, occasionally terminating in a lobe, yet somehow this must
counteract the cooling from a spherical halo of hot gas. The leading idea for
how this happens is that the jet lobes inflate cavities or “bubbles" of
overpressurised gas, generating pressure waves in the surrounding hot
atmosphere that then sphericalises the jet energy input over time. X-ray
observations indicate the presence of such bubbles in clusters and groups
(Bîrzan et al., 2004; Panagoulia et al., 2014; McDonald et al., 2018) lending
credence to this scenario, yet the link between jets, bubbles, halo gas
heating, and quenching remains tenuous.

Galaxy demographics indicate that widespread quenching must occur in the
regime of galaxy group (Rasmussen et al., 2012; Coenda et al., 2018; Salerno
et al., 2019). Isolated galaxies and those in poor groups such as the Local
Group tend to be star-forming, but moving to slightly larger halos finds a
growing fraction of quenched galaxies (Wetzel et al., 2012), and there exists
a similar trend also in the properties of their central galaxies (Weinmann et
al., 2006; Gozaliasl et al., 2016; Einasto et al., 2024). Satellites and
centrals also appear to be tightly correlated in terms of their star-forming
properties (Weinmann et al., 2006). Hence if black hole jets are the mechanism
that drives or maintains quenching, then it should be operating on this halo
mass scale. Indeed many galaxy groups are observed to contain a hot atmosphere
seen in the X-rays, i.e. an intragroup medium (IGrM) (Tully, 2015), which has
a complex, multiphase structure with cold, precipitated clouds and filaments
co-existing with a mass-dominant hot phase (Sharma et al., 2012; Prasad et
al., 2015, 2017; McCabe et al., 2021; McCabe & Borthakur, 2023; Saeedzadeh et
al., 2023).

The inferred cooling timescales in the centres of groups result in estimated
cooling rates of a few tens of $\mathrm{\ M}_{\odot}/$yr, in contrast to the
observed star formation rates and molecular gas masses in these systems which
puts the cooling rate at $\la 1-10\mathrm{\ M}_{\odot}/$yr (McDonald et al.,
2018). Hence despite the multi-phase nature of the IGrM, cooling appears to be
strongly suppressed throughout most of the hot gas.

To better understand how the jet energy suppresses cooling, many groups have
undertaken numerical simulations of this process. Early work focused on
constructing idealised halos of hot gas, and understanding how jets can
inflate bubbles (Vernaleo & Reynolds, 2006), excite gravitational modes as
cavities transform into overdensities (Omma et al., 2004; Reynolds et al.,
2005), and drive sound waves which heat the medium (Ruszkowski et al., 2004;
Fabian et al., 2005; Bambic & Reynolds, 2019). Subsequent work aimed at
simulating isolated halos but containing cosmologically motivated
substructure.

In these models, the plasma in the black hole accretion disk is heated and
uplifted via mechanical jets emanating from a central supermassive black hole
(SMBH) (Pope et al., 2010; Dubois et al., 2010; Gaspari et al., 2012, 2013;
Prasad et al., 2015, 2017, 2018). These jets then inflate large cavities, or
bubbles, of hot plasma in the IGrM, which shock the denser medium, heating the
gas and preventing cooling. Hence in these idealised situations, it appears
radio mode feedback is able to sphericalise jet energy input as long as the
jets are able to change direction on timescales that are short compared to the
dynamical and cooling timescales of the halo (Babul et al., 2013; Cielo et
al., 2018).

Separately, cosmological simulations of galaxy formation have attempted to
incorporate some version of jet feedback, with the primary motivation of
quenching massive galaxies. Sijacki et al. (2007) included the impact of jet-
driven bubbles via inflating bubbles by hand in the surrounding medium, which
showed some promise but was not able to fully quench galaxies. Dubois et al.
(2016); Kaviraj et al. (2017) directly modelled the jets accounting for black
hole spin within a cosmological setting in the HorizonAGN simulation (Dubois
et al., 2014), but likewise had trouble coupling the energy sufficiently to
produce a realistic quenched galaxy population. More recently, the effects of
AGN feedback have been included into cosmological simulations such as
IllustrisTNG (Weinberger et al., 2017) and Obsidian (Rennehan et al., 2023),
both of which have been able to produce a realistic quenched population. Both
models randomise the jet direction at every timestep in order to help
distribute energy in an isotropic sense. Locally isotropic active galactic
nucleus (AGN) feedback has been shown in the Romulus (Tremmel et al., 2017)
simulations to be able to still produce directional jets with varying
direction via interaction with an anisotropic gas profile in the cluster core
(Tremmel et al., 2019).

The Simba simulation successfully employs stably bipolar jets to produce a
quenched galaxy population in good agreement with observations (Davé et al.,
2019; Szpila et al., 2024). Simba stands apart in terms of its model as being
one that can achieve this whilst freeing the black hole jet axis to reorient
based on the angular momentum of gas within the accretion kernel, without
requiring enforced isotropy as several other simulations do. By comparing to
models run without jets, it is clear that galaxy quenching is primary driven
by the action of the jets, although the effect of an additional feedback mode
associated with X-ray photon pressure is not negligible (Cui et al., 2021).
Hence, Simba is at least one example of a full cosmological galaxy formation
simulation that is able to directly associate the quenching of massive
galaxies with the action of bipolar AGN jet feedback on the intracluster
medium (ICM) and ICM/IGrM, e.g. via heating and driving outflows. This begs
the question however on the exact nature of the energy transfer: Did Simba
achieve this by inflating bubbles in the surrounding gas and transferring
energy into the hot IGrM in a manner consistent with observations? This is the
question we seek to answer in this work.

Bubbles and cavities have now been observed in X-ray observations of many
groups and clusters , and allow one to make estimates of the jet power, the
time since the feedback event occurred, and the suppression of cooling within
the central regions of groups. There are however many aspects of AGN feedback
which we do not understand well. The energy output of individual feedback
events and their frequency, how exactly the jetted material propagates
through, and couples to, the hot IGrM, and to what extent cooling flows are
disrupted, are still not fully understood. One of the best ways to investigate
the nature of AGN feedback and its effects is through halo-scale, or, ideally,
cosmological simulations. Feedback is typically implemented in two modes; a
hot radiative, or "quasar", mode, operating at high Eddington fractions, and a
high-velocity "jet" mode (mechanical feedback) operating at low Eddington
ratios (Best & Heckman, 2012; Heckman & Best, 2014) (though see Rennehan et
al. 2023 for an example of an alternative implementation).

Simulations of AGN feedback in (generally idealised) clusters show that the
action of the jet on the IGrM or ICM gas is complex and non-linear. Feedback
has been shown to both shut-off star formation due to gas removal when active
(Sijacki & Springel, 2006; Sijacki et al., 2007; McCarthy et al., 2010; Le
Brun et al., 2014; Beckmann et al., 2017; Raouf et al., 2019; Cui et al.,
2021; Piotrowska et al., 2022), reducing stellar mass fractions in-line with
observation. However, in some circumstances, it can promote star formation
locally after the fact due to gas compression near shocks and through the
drawing out of low-entropy gas into filaments (Nayakshin & Zubovas, 2012;
Silk, 2013; Silk et al., 2014; Zubovas & Bourne, 2017; Huško et al., 2022).
Star formation is shown to be suppressed in and around radio lobes (Shabala et
al., 2011), and the promotion of star formation has been corroborated by
observations of sites of active star formation around X-ray cavities and jet-
induced shocks in multi-wavelength studies (Crockett et al., 2012; Vantyghem
et al., 2018).

Generally groups and clusters are simulated in idealised halos, either as
single systems or at best as a limited number of them. The computational cost
of running such simulations limits the sample size, and the idealised initial
conditions are in attempt to be able to draw general conclusions about
feedback, without findings being muddied by a complex halo thermodynamic and
velocity structure upon initialisation. In reality, groups are significantly
affected by their environment via gas inflow and accretion from the cosmic
web, tidal forces from neighbouring halos or subhalos, and mergers.
Furthermore, a representative sample of halos in terms of mass, assembly time,
and brightest cluster galaxy (BCG) morphology should be simulated for a true
statistical sample of feedback events. Finally, the feedback parameters in
these studies should be calibrated such that they produce the correct,
observed galaxy populations in cosmological simulations. The local properties
of jets can then be studied self-consistently and be predictive, as the
calibration is done on a global scale. This is what we aim to do.

An ideal platform for such a study is to use high-resolution zoom simulations
of a cosmological region encompassing a halo and its direct environment, with
a simulation model which is proven to reproduce realistic galaxy populations.
This gives the greatest balance between including large scale processes such
as halo accretion and subhalo mergers, as well as providing sufficiently high
resolution to be able to accurately simulate halo centers and resolve sub-
kiloparsec-scale gas physics. To this end we employ Hyenas, a cutting-edge
suite of zoom-in simulation boxes with halos selected based on virial mass and
formation time, simulated with the Simba model, which has a proven track-
record in reproducing observables in cosmological-scale boxes.

Our paper is organised as follows; In §2, we detail the Hyenas simulation set-
up as well as the group population we work with. In §3 we outline the
selection of cavities, and our analysis of them. We also present measurements
of bolometric cooling luminosities and accretion rates, linking them both to
cavity properties. §4 contains a few brief case-studies of groups selected
from our sample, and provides comparison between Chandra and LEM observations.
We conclude this work with a discussion of our findings in §5

##  2 The Simulations

###  2.1 Hyenas

The simulations in this paper are drawn from the Hyenas Level-1 suite of
simulation boxes (Cui et al., 2024), with the dark matter particle masses of
$8.1\times 10^{6}$ $M_{\odot}/h$ and gas particle masses of $1.5\times 10^{6}$
$M_{\odot}/h$(almost one order of magnitude resolution better than Simba). The
initial conditions are selected from a large dark matter box (L =
$200\,\mathrm{cMpc}\,h^{-1}$) spanning a range in mass and halo formation
time, with the aim of choosing a representative sample of halos. The galaxy
formation model, including the feedback prescription, is essentially identical
to Simba (Davé et al., 2019). We now outline the most important features of
the Hyenas (Simba) model in the context of this project, with a main focus on
the active galactic nuclei (AGN) feedback implementation. For complete details
of all aspects of the model, we direct the reader to the original code paper
(Davé et al., 2019).

The Simba model is unique within the set of current cosmological simulations
in its treatment of accretion onto black holes. Accretion is split by the
local gas properties into a Bondi-Hoyle-Lyttleton (BHL) mode and a torque-
limited accretion mode (Hopkins & Quataert, 2010, 2011) depending on the gas
temperature.

The hot gas ($T>10^{5}$K) is accreted onto black holes using a BHL accretion
model under the assumption of spherical symmetry. The mass accretion rate is
computed Davé et al. (2019);

$\dot{M}_{\text{Bondi}}=\varepsilon_{m}\frac{4\pi
G^{2}M^{2}_{\text{BH}}\rho}{(v^{2}+c_{s}^{2})^{3/2}},$ (1)

where $\rho$ is the density of the hot gas within the accretion kernel,
$c_{s}$ is an averaged sound speed of this gas, and $v$ is an average relative
velocity compared to the black hole itself. $\varepsilon_{m}$ is a tunable
efficiency factor, set to $0.1$ in Simba.

For cold gas ($T<10^{5}$K, $n_{\mathrm{H}}>0.13\,\mathrm{cm}^{-3}$), the
torque-limited accretion model (Anglés-Alcázar et al., 2017), yields an
accretion rate

$\dot{M}_{\text{Torque}}\approx\varepsilon_{T}f_{d}^{5/2}\times\left(\frac{M_{\text{BH}}}{10^{8}M_{\odot}}\right)^{1/6}\left(\frac{M_{\text{enc}}(R_{0})}{10^{9}M_{\odot}}\right)\left(\frac{R_{0}}{100\text{pc}}\right)^{-3/2}\times$
$\left(1+\frac{f_{0}}{f_{\text{gas}}}\right)^{-1}M_{\odot}\text{ yr}^{-1},$
(2)

where $f_{d}$ is the disk mass fraction (including stellar and gas
components), $M_{\text{enc}}(R_{0})$ is the total sum of these two components,
$f_{\text{gas}}$ is the gas mass fraction in the disk, and $R_{0}$ is the
distance from the black hole which encloses the nearest $256$ distinct gas
elements. $f_{0}$ is a factor that describes the amplitude of epicyclic
oscillations of a test particle around a perturbed circular orbit in the
central accretion disk - see (Hopkins & Quataert, 2011) for a full definition.
A kinematic decomposition is used during the simulation to separate the disk
and spheroidal components, each comprised of gas and star particles. The
interstellar medium (ISM) in the Simba galaxy model is artificially
pressurized in order to resolve the Jeans mass, so for accretion we treat all
ISM gas ($n_{\mathrm{H}}>0.13\,\mathrm{cm}^{-3}$) as "cold". We assume
$\varepsilon_{T}=0.1$ which is tuned to match the amplitude of the black hole
mass – stellar mass (or velocity dispersion) relation at $z\approx 0$.

The accretion rate is modulated by a radiative efficiency factor $\eta=0.1$,
and is given in its final form as
$\dot{M}_{\mathrm{BH}}=(1-\eta)\times(\dot{M}_{\text{Torque}}+\dot{M}_{\text{Bondi}})$.
We impose upper limits on the rates of both components of the total accretion
rate based on the Eddington accretion rate
$\dot{M}_{\text{Edd}}(M_{\text{BH}})$ for each black hole.

The feedback itself is also implemented in two modes; kinetic and X-ray. The
kinetic mode is implemented through a sub-grid model which is split into two
sub-modes, determined by the Eddington fraction
$f_{\text{Edd}}\equiv\dot{M}_{\text{BH}}/\dot{M}_{\text{Edd}}$ of the
accreting black hole. For high Eddington ratios ($\gtrapprox 0.2$) the AGN
drives winds with velocities on the order of $10^{3}$
$\mathrm{km}\,\mathrm{s}^{-1}$, which consist of molecular and warmly ionized
gas. We call this the radiative mode. At lower Eddington fractions a jet mode
is instead active, with hot gas being driven at velocities of order $10^{4}$
kms-1. Jet-mode feedback is effective at evacuating gas from the central
regions through entrainment (Borrow et al., 2020). Both modes of kinetic
feedback are implemented in a bipolar manner, with the jet axis calculated
according to the angular momentum of the inner gas disk.

We parameterise the radiative wind velocity as a function of the black hole
mass

$v_{\text{w,EL}}=500\left(1+\frac{\log_{10}M_{\text{BH}}-6}{3}\right)\text{ km
s}^{-1},$ (3)

with the gas being ejected at the ISM temperature. If the Eddington fraction
$f_{\mathrm{Edd}}<0.2$, the feedback instead occurs via the jet mode and the
velocity is boosted to

$v_{\text{w,jet}}=v_{\text{w,EL}}+7000\log_{10}(0.2/f_{\text{Edd}})\text{ km
s}^{-1}.$ (4)

We impose an upper limit of $7000$ kms-1 at $f_{\text{Edd}}<0.02$.
Additionally, jets will not activate unless the black hole mass reaches
$\gtrsim 10^{7.5}M_{\odot}$, in order to stop small, but slowly-accreting
black holes from generating jets. Gas ejected in jets is automatically raised
to the virial temperature of the halo in Simba.

The final AGN feedback mode present in Simba is in the form of X-ray feedback,
which models the radiation pressure from X-rays off the accretion disk. This
feedback mode is only active when full-velocity jet feedback is occurring, and
so is again determined by the mass accretion rate. Additionally, X-ray heating
is only active if $f_{gas}\equiv M_{gas}/M_{*}<0.2$. Simba heats gas within
the kernel of the BH, with the energy input scaled by the inverse square of
the distance from the black hole. X-rays heat the non-ISM gas directly,
whereas for ISM gas, which would be expected to just cool quickly, only half
the energy is added as heat, with the other half being used to give the gas a
kick radially outwards. The gas heating rates due to the X-rays is given in
Choi et al. (2012) as

$\dot{E}=n^{2}(S_{1}+S_{2}),$ (5)

where $n$ is the local proton number density, $S_{1}$ is the power
contribution from Compton heating, and $S_{2}$ is the contribution coming from
the sum of photoionization heating. The radiation-pressure change in momentum
of the kicked gas is straightforwardly calculated by
$\dot{p}=\dot{E}^{\prime}/c$ where $E^{\prime}$ is half the total energy
available.

The free parameters in the feedback model are broadly tuned to match
observations of the galaxy stellar mass function evolution and the stellar
mass–black hole mass relation. The X-ray feedback mode is included primarily
to obtain enough fully-quenched galaxies at $z=0$. Simba has been tested
against a large number of other galaxy, halo, and intergalactic medium
properties, with good alignment with available observations (e.g. Davé et al.,
2019) at a level comparable or better than other similar cosmological
simulations. For example, the X-ray scaling relations and baryonic content in
groups and clusters (Jennings & Davé, 2023), the gas content of the
circumgalactic medium (CGM) including observed $\mathrm{HI}$ equivalent width
and $\mathrm{OVI}$ absorbers (including the observed dichotomy around quenched
versus star-forming galaxies) (Appleby et al., 2021; Bradley et al., 2022),
the fraction of baryons locked into the intergalactic medium (IGM) at low
redshift (Sorini et al., 2022), black hole scaling relations, the satellite
galaxy stellar mass function and color-magnitude relationship, and stellar
properties of the BCG (Cui et al., 2021), richness of environments of radio
galaxies (Thomas et al., 2021), general agreement in the fraction of rapidly
quenched quiescent galaxies at $z=1$ (Zheng et al., 2022), and good agreement
with the number density of quenched galaxies between $z=3-4$ (Szpila et al.,
2024).

###  2.2 Resimulation technique

In order to capture the short timescales of the jet and bubble-inflation
processes, we select a Hyenas snapshot produced at $z=0.286$ and restart the
simulation from that point. This allows us to increase the snapshot cadence to
every $10$ Myr. We run the simulation for an additional $700$ Myr in total. We
check that the gas, stellar, and black hole properties are consistent after
restarting, yet still allow a grace period of a few tens of Myr before the
earliest analysis is done. This is to mitigate any potential transient effects
which may be present due to restarting from snapshot initial conditions within
Gizmo.

In Fig. 1 we plot the distribution of initial halo masses at the first
snapshot at redshift $z=0.286$. The range of halo masses that we simulate (and
find bubbles) is $10^{13}M_{\odot}$ to $10^{14}M_{\odot}$.

Figure 1: The initial mass distribution of halos for which we find bubbles as
measured at the first snapshot at $z=0.286$.

###  2.3 X-ray Maps

To make comparisons to existing catalogues of X-ray cavities, we produce mock
Chandra observations using the MOXHA package (Jennings & Davé, 2023). We use
the Chandra ACIS-I cycle-0 specification111Response files as known in SOXS:
acisi_aimpt_cy0.rmf, acisi_aimpt_cy0.arf, which does not suffer from the
degredation of effective area at low energies due to contamination build-up
over operating time compared to later cycles (Plucinsky et al., 2018). Our
motivation for using Chandra is that there is a wealth of observed bubbles in
the literature, and the instrument has excellent angular resolution
($0.492^{\prime\prime}$ per pixel) which is needed to spatially resolve
cavities. We artificially remove the Chandra chip gaps before observation, so
as to emulate a mosaic-ed observation. This means that our full field of view
is $16.6^{\prime}\times 16.6^{\prime}$, with approximately $4$ times the
single-chip area.

In order to offer insights into X-ray cavity detection using the next
generation of X-ray observatories, we produce mock observations using the Line
Emission Mapper (LEM) probe, with the 2eV spectral resolution
specification222Response files as known in SOXS:: lem_2ev_110422.rmf,
lem_300522.arf. We present some comparison X-ray maps in §5 and also perform
spectral fits for all halos for both Chandra and LEM, comparing the accuracy
of their recovered source luminosities when fitting to a jet-disturbed halo.

We create the X-ray maps using PyXSIM to generate photons from an APEC model
(Smith et al., 2001) and project them along the line of sight. We include the
effects of neutral hydrogen absorption from the Milky Way, and use SOXS to
observe these photons via convolution with the Chandra and LEM instrumental
responses. Furthermore, we select an observation redshift for each halo that
results in the field of view on the Chandra mosaic being $2\times 0.6R_{500}$
from edge-to-edge, where $R_{500}$ is measured at the initial redshift
$z=0.286$. This is to standardise the observations over the whole order-of-
magnitude range of halo masses, and prevent biasing against observations of
small X-ray bubbles in low-mass groups due to reduced resolution which would
occur if we instead held redshift constant. Our field of view choice results
in a range of observational redshifts: $0.0377-0.0172$. We use an exposure
time of $75$ks, which is consistent with real observations of clusters with
identifiable X-ray cavities.

###  2.4 Thermodynamic Profiles

We require radial thermodynamic profiles so that we can sample the
temperature, pressure, and (electron) density at given radius within a halo.
To this end, we generate emissivity-weighted electron density and temperature
profiles from the ISM-filtered gas particle data and fit the parameterised
models below.

For the temperature we follow Vikhlinin et al. (2006):

$T_{3D}(r)=T_{0}t_{\text{cool}}(r)t(r)$ (6)

where outside the central cooling region the profile is dominated by

$t(r)=\frac{(r/r_{\mathrm{t}})^{-a}}{\Big{[}1+(r/r_{\mathrm{t}})^{b}\Big{]}^{c/b}},$
(7)

and the central regions are described by

$t_{\text{cool}}(r)=\frac{(x+T_{\text{min}}/T_{0})}{x+1}\ \ \ \ \ \ \ \ ,\ \ \
\ \ \ \\\ x=\Bigg{(}\frac{r}{r_{\text{cool}}}\Bigg{)}^{a_{\text{cool}}}.$ (8)

Note that in our groups sample, the central regions are generally only at-best
partially captured on orders of $\sim$ a few kpc scales and below.

For the electron density $n_{e}$, we use the modified $\beta$ model from
Vikhlinin et al. (2006):

$n_{\mathrm{e}}^{2}=n_{\mathrm{e,0}}^{2}\frac{(r/r_{\mathrm{c}})^{-1}}{(1+r^{2}/r_{\mathrm{c}}^{2})^{3\beta}}\frac{1}{\left(1+(r/r_{\mathrm{s}})^{\gamma}\right)^{\frac{\varepsilon}{\gamma}}}.$
(9)

We fit both the temperature and density profiles individually to the
corresponding emissivity-weighted radial profiles out to $100$ kpc. For this
work, we use the emissivity between $0.2-3.0$ keV (source-frame) for the
profiles. We reject fitted profiles if the Coefficient of Determination value
(for model values $f_{i}$ and data $y_{i}$ this is
$R=1-\frac{\sum(y_{i}-f_{i})^{2}}{\sum(y_{i}-\bar{y})^{2}}$, which we
evaluated in logarithmic space) is $<0.5$ for either or both fits. In this
case, the halo is discounted from further quantitative analysis. We
furthermore obtain a radial pressure profile by multiplying our best-fit
temperature and electron density profiles. In doing this we assume a mean
molecular weight per electron of $\mu_{e}=1.155$ in order to mimic the usual
assumptions made when these profiles are derived from observations.

##  3 Cavity Analysis

###  3.1 Selecting Cavities

The first step in our study of X-ray cavities is finding them. They are mainly
characterised by depressions in surface brightness on scales of a kiloparsec
to a few tens-of-kpc close to the cluster or group centre. Several
observational techniques have become standard in the field (see for example
Hlavacek-Larrondo et al., 2015), and we use them in this work to match the
observations as closely as possible.

We search for cavities in lightly-smoothed exposure-corrected X-ray flux maps.
We focus on the region within $0.3R_{500}$ (measured from the minimum
potential position, which coincides with the central black hole position),
where the intensity contrast is likely to be greatest. Furthermore, we
generate a residual X-ray map, extracting the mean pixel flux of a lightly-
smoothed X-ray flux map in radial bins about the cluster center. We fit a King
model (King, 1966) to this mean radial flux profile

$I(r)=I_{0}\left(1+(1/r_{0})^{2}\right)^{-\beta}+C_{0},$ (10)

where $I_{0},r_{0},C_{0}$ are the normalization, core radius, and a constant
variable, respectively — all are free parameters in the fit. We then subtract
off the King model counts from the raw X-ray counts map, and divide by the
same King model, before lightly smoothing.

Next, we generate unsharped maps. These are designed to highlight features on
a particular scale, by subtracting a strongly-smoothed image from a lightly-
smoothed counterpart. We use Gaussian smoothing on several different pairs of
scales. The main scales we use to search for cavities are $(60,2)$kpc and
$(80,5)$kpc.

Finally, we run the machine-learning cavity detection pipeline cadet
333https://github.com/tomasplsek/CADET (Plšek et al., 2023) on the exposure-
corrected flux maps. cadet comprises a convolutional neural network that has
been trained on Chandra images to make pixel-wise predictions for the
identification of cavities. We apply cadet to the flux image subjected to
various different scales of cropping in order to test the robustness of a
given detection and to highlight cavities on different scales. We find that
cadet does well in identifying cavities in the vast majority of cases, and in
several cases can recover a very good match to the true region of low density
from comparison with the true projected maps even when the corresponding
surface brightness depression is not clear in the smoothed and/or filtered
X-ray maps.

In all we therefore have four different types of spatial map used for cavity
detection, that would be available to an observer; X-ray flux ($0.5-7$ keV),
King model-subtracted residual, unsharped, and cadet pixel-wise predictions.

We selected cavities based predominantly on mock X-ray quantities, but
occasionally we use the true thermodynamic maps (and their time-evolution) to
inform the direction and rough size and shape of bubbles, in order to reduce
false positive detections and poor shape fits to the surface brightness
depressions. Therefore, our search is not a truly blind test. However, we
point out that due to computational and time limitations, we are not afforded
access to other signatures of bubbles such as (a) tesselated spectrally-fitted
spatial maps of temperature and density over the group, or (b) a multi-
wavelength analysis looking for radio emission in the directions of cavities.

For each cavity where we observe double-cavities, we select only bubbles
launched at approximately the same time. Older, relic bubbles are ignored when
they appear in snapshots where younger bubbles/s are observed for the purposes
of calculating the recent AGN power; these relic bubbles may still be analysed
earlier in their lifetimes corresponding to a preceding snapshot and period of
jet activity, if they are clear at that earlier time.

In Fig. 2 and Fig. 3 we plot our sample of observed cavities in our catalogue.
Shown are the Unsharp-masked and King-model subtracted maps, respectively,
over which we overlay cadet confidence contours and our fitted cavity
ellipses. Qualitatively, we see that the cavities produced in Hyenas are
similar to those seen in real systems (e.g. Hlavacek-Larrondo et al., 2015),
and are largely well-identified by cadet. Later on in this section we shall
demonstrate that they are also quantitatively similar in terms of size and
location.

(a)  Figure 2: Unsharp-masked maps are show for a selection of our galaxy
groups. They are generated by subtracting a strongly-smoothed X-ray map from a
lightly-smoothed one; we use $\sigma=60,2$ here respectively. The orange and
red contours show the cadet predictions at a significance level of 90% when
run on two different cutout scales ($1024$ and $1280$ Chandra pixels
respectively), and we show the fitted ellipses that we identify in blue. The
scale bars denote $50$ kpc. Note that the colormap normalisation varies
between each plot and is designed to maximise contrast in the radius range
occupied by the bubble/s. Some apparent cavities are not designated so as they
are ruled out via consultation of maps and movies of the jets in the true
thermodynamic maps, or are relic cavities not associated with the most recent
bouts of jet activity, but possibly included in the analysis of a previous
snapshot.  (a)  Figure 3: Same as for Fig. 2 but for the King-model subtracted
X-ray images.

###  3.2 Cavity Uncertainties

Figure 4: An example (Zoom252_halo0 (snapshot 188)) of the posterior
distribution of the derived enthalpy obtained from our uncertainty estimation
method, showing it’s dependence on the uncertainty in the bubble measurements
and inclination.

The uncertainties on observed X-ray cavity sizes and, therefore, enthalpies,
are generally large. Unfortunately, it is difficult to identify their exact
dimensions and, simultaneously, it is difficult to measure the thermodynamic
radial profiles of the halo (Hlavacek-Larrondo et al., 2015). Therefore, we
must quantify these uncertainties accurately.

First, we assign a typical uncertainty of $20\%$ to the minor and major axes
of the fitted ellipse, and also to the radial distance of the bubble center
from the halo center. These are in line with observational studies (see e.g.
Hlavacek-Larrondo et al. 2015). We assume that the errors on these quantities
are Gaussian distributed.

Second, we consider uncertainties arising from bubble inclination. In general,
a nonzero, unknown inclination will move the true bubble center to a larger
radius than the projected bubble midpoint radius. We assign a uniform
distribution to this angle $\theta$, under the assumption that for us to
clearly identify a bubble the inclination can be no more than $30^{\circ}$
inclined from the plane perpendicular to the line of sight (i.e. the plane of
observation). We take $10,000$ enthalpy samples
${H_{\text{tot}}(\boldsymbol{\Theta}_{i})}$ with parameter values drawn from
these aforementioned distributions;

$H_{\text{tot}}(\boldsymbol{\Theta}_{i})=\frac{5}{2}*p\Big{(}r_{\text{proj,i}}/\sin(\theta_{i})\Big{)}*\frac{4}{3}\pi
a_{i}^{2}b_{i},$ (11)

where $r_{\text{proj}}$ is the projected bubble midpoint radial distance. We
assume cylindrical symmetry of the ellipse about the radial unit vector, so
$a$ is the tangential ellipse axis and $b$ the radial axis of the ellipse,
with the origin at the halo center. Therefore, we have

$\displaystyle
r_{\text{proj,i}}\sim\mathcal{N}\big{(}r_{\text{proj}},(0.2r_{\text{proj}})^{2}\big{)},$
$\displaystyle a_{i}\sim\mathcal{N}\big{(}a,(0.2a)^{2}\big{)},$ $\displaystyle
b_{i}\sim\mathcal{N}\big{(}b,(0.2b)^{2}\big{)},\mathrm{and}$
$\displaystyle\theta_{i}\sim\mathcal{U}\big{(}\frac{\pi}{3},\frac{\pi}{2}\big{)}.$

We assume that the change in the bubble volume due to inclination $\theta_{i}$
is small compared to the measurement uncertainties, and we also assume that we
know the pressure profile $P(r)$ perfectly. In reality, there is a non-
negligible uncertainty attached to the 3D de-projected radial profile.

Using this method, we can account for most of the uncertainties in the
measurements, and the usual assumption of zero inclination. Most importantly,
we account for the uncertainty in the change in the radius at which the
pressure is sampled. We show an example corner plot of the enthalpy posterior
of one of our cavities in Fig. 4. We find the resulting posterior
distributions for the logarithm of the enthalpy for our cavities are
approximately normally-distributed, and so we use the $16^{th}$ and $84^{th}$
percentiles as the $1\sigma$ error for each individual cavity. If we have two
bubbles present, the error on the final enthalpy is the sum in quadrature of
the individual enthalpy uncertainties.

###  3.3 Cavity Dimensions

In this Section, we investigate the physical characteristics of the X-ray
cavities in our Hyenas sample. In order to examine whether the simulated gas
properties are accurate in terms of where the bubbles are injected and how
large they grow once seeded, we compare to the observational results of
Rafferty et al. (2006) who present a sample of $62$ cavities in galaxy
clusters, and Shin et al. (2016), who present a study of 148 detected cavities
from a sample of 69 clusters, groups, and elliptical galaxies.

Figs 5 and 6 show that our cavity sample matches the Rafferty et al. (2006)
and Shin et al. (2016) data in both the linear trend and the scatter. Hyenas
cavities tend to be found at the larger end of the linear trend that is traced
by the observational data, but note this observational data includes also
individual elliptical galaxies. The resolution of the simulation could impact
these results since it would prevent clear, small cavities from being observed
very close to the centres of small groups on sub-$10$kpc scales. An
alternative explanation is that Simba’s strong feedback scheme is very
effective at quickly driving winds to large radii and inflating large
cavities. It is possible that Hyenas bubbles are slightly too large or
injected slightly too far from the centre of the halo, since our points have a
tendency to lie slightly below the Rafferty et al. (2006) and Shin et al.
(2016) data.

In Fig. 7, we show the eccentricity of Hyenas cavities from our mock
observation technique, comparing to data from (Rafferty et al., 2006) and
(Shin et al., 2016), matched to the midpoint distance range of this work. We
find promising agreement in the position of the peak of the eccentricity
distribution and in the general scatter, though we tend to have a long-tail
towards low-eccentricity compared to observations, as we show in the
normalised histogram presented alongside. Overall, we demonstrate that Hyenas
cavities are found in the right locations, with the correct sizes and shapes
when compared to real bubbles.

Simba’s feedback model, including the re-coupling timescale of jet outflows,
as well as the energy and momentum transfer rates, is tuned to reproduce both
the galaxy stellar mass function (GSMF) and the quenched fraction at low $z$.
Simba was not calibrated to match the general hot gas properties of halos,
let-alone the smaller scales of bubble injection. Therefore, it is encouraging
that a realistic X-ray cavity population emerges from a bipolar feedback
subscription that is calibrated on galaxy population properties. We suggest
that a bipolar feedback model such as Simba’s is able, and possibly required,
to match observable signatures of AGN feedback on the IGrM and ICM on the
scales of tens of kiloparsecs, whilst also producing realistic galaxy
populations.

Figure 5: The bubble surface area as seen in the plane of observation compared
to the observational data of Rafferty et al. (2006) and Shin et al. (2016).
Our bubbles match the relation and spread between area and midpoint distance
well. For clarity, no error bars are shown we ascribe an error of $20\%$ to
the midpoint distance and$\sqrt{0.2^{2}+0.2^{2}}\approx 28\%$ to the bubble
volume. Highlighted in blue is the range of midpoint distances we measure.
Figure 6: Same as Fig. 5 but for the semi-major axis length. Again we match
the observed trend of Rafferty et al. (2006) and Shin et al. (2016) well.
Figure 7: We plot the (projected) eccentricity of the Hyenas bubbles (black)
against those of the published data of Rafferty et al. (2006) (red) and Shin
et al. (2016) (blue), from which we restrict to those corresponding to our
range of midpoint-distances (blue band in Figs 5, 6). Contours correspond to
$50$th, $70$th, and $90$th percentiles. Also plotted is the normalised
histogram. We match observations reasonable well, with an average eccentricity
at $\sim 0.7$, only slightly lower than the observed trends. We tend to have a
slight overabundance of cavities at low eccentricities of $0.2-0.4$.

###  3.4 Cavity Power Estimation

In this Section, we investigate the simulated, X-ray determined cavity power
compared to the true jet power. Additionally, we compare the fraction of jet
power being contributed by the BHL and torque-limited accretion modes. For the
X-ray measured cavity power, we select a snapshot at which a detected bubble
appears clearly in the intensity map and/or in the cadet prediction. We assume
that the cavity is in approximate pressure equilibrium with the IGrM, such
that $P_{cav}\approx P_{\text{ICM}}\equiv p$, and therefore the total enthalpy
is (Bîrzan et al., 2004):

$H_{\text{tot}}=U_{\text{int}}+PV=\frac{1}{\gamma-1}PV+PV=\frac{5}{2}pV$ (12)

which gives the (minimum) energy that is needed to account for both the
internal energy and the energy required for an adiabatic bubble to expand
against the ICM (at a constant pressure), and we assume (as in Simba) that
$\gamma=5/3$. Note, this differs slightly from the usual observer assumption
that the adiabatic index $\gamma=4/3$ since the cavity gas in real systems is
thought to be relativistic (Graham et al., 2008; Fabian, 2012).

To determine the ambient pressure at the appropriate radius for each halo and
cavity event, we take the best-fit radial profiles of the electron density
$n_{e}$ and the plasma temperature $k_{B}T$ using the PyXSIM-derived,
emissivity-weighted filtered gas, as described in §2. We assume that the total
pressure follows $P=(\mu_{e}/\mu)*n_{e}k_{B}T$, computed at the bubble radius.
$\mu_{e}=1.155$ here is the mean molecular weight per electron. The volume of
the cavity is measured by associating an ellipse to each bubble based on the
size, location, and orientation of the features seen in the flux map and the
filtered images, as well as the cadet predictions. We assume symmetry about
whichever axis (semi-minor or semi-major) aligns closest to the radial vector
pointing out from the halo center. We calculate the enthalpy by computing the
pressure at the radius of the bubble centre; unless otherwise stated, this is
what we present as being used for the cavity enthalpy/power. For comparison,
we also calculate the enthalpy from the pressure at the radius of the
outermost point on the bubble perimeter. This gives a lower limit on any
enthalpy that could be calculated from a method that instead considers a
radially-varying pressure profile within the bubble as opposed to a
constant/midpoint value.

We calculate the AGN energy injected into the halo via kinetic jets directly
from the simulation data, with bubble gas being isolated via a cut on the hot
gas phase space and its change between the snapshots identified as containing
the jets contributing to the energy of the bubble. We identify cavity members
in a first sweep as any particle that between any two snapshots ($10$Myr
spacing) has satisfied all of the following: the particle has increased its
internal energy by at least $40\%$, reduced its Lagrangian density by at least
$40\%$, has a post-kick velocity and velocity increase of at least $200$ km/s,
and lays within $20$ kpc with respect to the halo center before being heated.
A second sweep is performed to also identify particles which may be initially
hot (so that their internal energy may not have increased by the above
percentage threshold) but which are also entrained in outflows. This sweep
finds outflowing particles which have had their internal energy raised by the
minimum energy increment found in the particles identified from the first
sweep. We visually inspect projections of thermodynamic properties that these
criteria result in bubble/jet particles being clearly identified.

Furthermore, we calculate the AGN energy via a second method by simply summing
the wind kinetic energies (which are saved in Simba at each timestep) over the
time range that the jets responsible for inflating the cavities are active. We
label these two schemes as "Gas Energies" and "Wind Energies" in Fig. 8.

(a)

(b)

Figure 8: On the left we plot the observed bubble enthalpies as calculated
using the midpoint pressure against the AGN energy, calculated using the two
methods described in the text. The scatter shows data for the observed
enthalpy versus the Wind Energies. On the right we plot the same but for using
the enthalpy derived with the pressure at the bubble’s outermost tip. Contours
correspond to $20$th, $50$th, and $80$th percentiles. For the midpoint energy
the bubbles are over-powered, but a 1:1 relation is recovered when using the
lower pressure.

The left plot in Fig. 8 shows that the cavity enthalpies (from a midpoint
pressure) generally exceed the apparent energy input by the AGN over the
corresponding jetting time-frame 444The duration of time that jets directly
responsible for the inflation of the bubble/s under cosnideration are actively
recoupling and driving outflows. The energy estimated from the change in gas
energies around the brightest group galaxy (BGG) is slightly lower than the
value obtained from the winds themselves. We posit that this is caused a
combination of cooling losses, and from gentle heating by weak shocks and
sound waves, which may not be captured through our tagging scheme. We leave
this for future work to quantify.

It is interesting that both of our AGN energy estimates appear too low in
power compared to the enthalpy contained in cavities. The fact that the cavity
dimensions are consistent with observations (if anything, we produce slightly
under-sized bubbles) implies that the cavity volumes are correct. There are
two possible explanations. First, there could be a significant decrement in
pressure along "channels" in the IGrM carved by older feedback episodes, that
can enhance the apparent energy contained within a freshly inflated bubble by
reducing the amount of $PdV$ work it must do in reality compared to inflating
within pristine halo gas. Second, the enthalpy prediction for the bubble is
too high. To test the latter option, we compute the enthalpy using the lowest
pressure that the bubble can "see" — the pressure at the outermost tip of the
bubbleplot. We show this result in the right plot in Fig. 8. This gives the
lowermost limit of the enthalpy within the bubble since, at these radii, the
pressure is monotonically decreasing with radius from the halo centre. There
is a clear linear trend between the bubble enthalpies and the AGN output, with
very little offset (but a large scatter). This indicates that the more
important physical quantity is the enthalpy calculated from the outermost
pressure, and not from the midpoint. We argue that the bubble inflation
process is dominated by work done along the radial axis into the IGrM at the
bubble tip, where pressure is lowest and where expansion requires the least
work. So long as the resistance provided by the "added mass" that the bubble
must uplift is comparatively small, this radial direction offers the "path of
least resistance" for the bubble to expand along, and it will do $PdV$ work
using the lowest pressure the bubble can "see". This preferential direction of
expansion also explains why both observed and simulated cavities are so
eccentric (Fig. 7).

To compute the cavity ages, we use the characteristic timescales based on the
instantaneous position and thermodynamic state of the bubble when observed,
and via a Direct estimate making use of the many snapshots that we actually
have access to from the simulation. To transform the observed cavity energy
into a power, we calculate the buoyancy timescale associated with each cavity
(Bîrzan et al., 2004)

$t_{\mathrm{buoy}}=r_{\text{mid}}/v_{t}=r_{\text{mid}}\sqrt{\frac{SC}{2gV}},$
(13)

where $v_{\mathrm{t}}$ is the terminal buoyant rising velocity (obtained
simply by equating the ram pressure force and the buoyancy force (Churazov et
al., 2001), $S$ is the surface area of the bubble, $C=0.75$ is the drag
coefficient (Churazov et al., 2001), $g$ is the gravitational acceleration at
the bubble midpoint, and $V$ is the bubble volume. We calculate $S$ by
assuming a circular cross-section with radius equal to the ellipse axis that
aligns best with the tangential unit vector around the group centre. We
calculate the gravitational acceleration $g$ by assuming hydrostatic
equilibrium and using our best-fit pressure profile. For the error
calculation, we calculate this quantity for each sample of our cavity
properties which are generated as described in Section 3.2. Note that we treat
the thermodynamic profile as exact in that we do not consider any uncertainty
around the best-fit fitting parameters (although the midpoint of the bubble at
which we sample the profile will be different for each sample) and so our
error estimates will be under-estimated compared to actual observational data.

Additionally, we require an estimate of the sonic timescale (Bîrzan et al.,
2004), which is the time-scale for the cavity to travel from the central AGN
to the bubble midpoint radius $r_{\mathrm{mid}}$. Therefore,

$t_{\mathrm{c_{s}}}\equiv
r_{\text{mid}}/c_{\text{s}}=r_{\text{mid}}\sqrt{\frac{\mu
m_{\mathrm{p}}}{\gamma k_{B}T}},$ (14)

where $c_{\mathrm{s}}$ is the speed of sound in the IGrM and we take
$\mu=0.62$. Compared to the free-fall time or buoyant-rising time, the sonic
time-scale does not depend on estimates of the local gravitational
acceleration, which can be very uncertain (Hlavacek-Larrondo et al., 2015).
Furthermore, it is more accurate for high-powered jets, where the bubbles are
driven outwards rather than buoyantly inflating (Omma et al., 2004). Because
Simba has a powerful jet model, this timescale may be the more appropriate
measure of the cavity lifetime.

For our Direct estimates, we find the time difference between the snapshot
where we observe the cavity, and the midpoint of the timescale we identify
associated jet activity. For example, if we identify jets contributing to a
given bubble over a timescale of $20-30$ Myr before we measure the bubble, we
associate a Direct cavity timescale of $25$ Myr. We also calculate a Maximum
Direct timescale, simply the time between the bubble observation and the
estimated start time of all associated bubble activity (this would be $30$ Myr
in the example above). These estimates are based on visual inspection of the
density projection movies, with cross-referencing to large spikes in the time-
binned wind kinetic energy.

Figure 9:  The Kernel-Density Estimate for the distribution of our estimated
cavity ages is shown, with the normalized distribution density on the y-axis
plotted against the cavity age. The best-fit Gaussian distribution is
calculated for the curves above a normalised density of $0.3$, and its mean
and standard deviation stated. We compare to the observational data of
Rafferty et al. (2006) that lies within our minimum and maximum midpoint
distances as plotted in Fig. 5.

We plot the kernel density estimates (KDEs) for the distributions of our age
estimates in Fig. 9 for comparison555Using the scipy statistics library..
Median and $1-\sigma$ errors were estimated from samples drawn from the KDE
distribution above a normalized density of $0.3$ and fit to a Gaussian. All
distributions are well-fit by a Gaussian above a normalised density of $0.3$.
The direct estimate is the shortest, with an average timescale of $30.3$ Myr
and a standard deviation of $15.0$ Myr. Next are the Maximum Direct and
buoyancy timescales, followed by the sound-speed/sonic timescale, with best-
fit mean values at $42.7,49.2,53.0$ Myr respectively. We see there is a
relatively large spread in the time taken for cavities to migrate from the
central few kpc, even for the most conservative estimate. The fact that the
sound speed timescale is significantly larger than the others may be a result
of the rather arbitrary definition of the sound-speed, where the bubble is
assumed to always travel at a Mach number of unity. In reality, as long as the
shock remains weak in order to reproduce the observational dearth of strong
cavity shocks and the gentle heating observed in systems, this cavity speed
can be larger and the sound-speed timescale can be brought more in line with
the buoyancy estimate. Indeed, observations have shown that lobes can possess
Mach numbers greater than unity, and of up to $\mathcal{M}\approx 2$
(Simionescu et al., 2009; Kraft et al., 2012; Sonkamble et al., 2015;
Vagshette et al., 2017). Modifying the formula for the sound-speed cavity age
by a universal conservative Mach number $\mathcal{M}=1.1-1.4$ would align
$t_{cs}$ with $t_{\mathrm{buoy}}$ and $t_{\textit{MaxDirect}}$ resepectively.

The buoyancy timescale and the maximum direct timescale distributions agree
reasonably well in the mean and scatter. The maximum direct timescale
corresponds to the bubble launched by the first jet during a period of
activity, and so it is this bubble that experience the maximal drag force,
which is set by the (undisturbed) hot atmosphere. Subsequent bubbles during
the same period of activity, blowing along the same axis of outflow, will
experience a lowered drag force and so a shorter cavity lifetime, because of
the fact that dense gas has already been evacuated from this channel.
Therefore, it makes sense that these two estimates for the cavity lifetime are
in good agreement and, furthermore, indicates that the buoyancy timescale is
likely close to being an unbiased estimator for the true cavity lifetime, at
least for the first bubble launched along an axis. Finally we note that the
Hyenas $t_{\text{buoy}}$ matches the buoyancy timescales as measured by
Rafferty et al. (2006) very well in terms of the average, although the scatter
on our values is smaller by around a third.

##  4 Case Studies of Individual Groups

We begin by considering a few case studies of selected individual galaxy
groups that display interesting behaviour. We particularly highlight features
in our X-ray maps that are reminiscent of real systems, to provide qualitative
support for the idea that the way in which Hyenas’s jets interact with the
IGrM is reasonably realistic.

In Fig. 10, we show one of our largest halos with an initial $M_{500}$ of
$10^{13.77}M_{\odot}$. This halo has some of the clearest, most well-defined
bubbles we see in our sample. A north and south cavity system is clearly
visible in both the Chandra and LEM mock images, which have been lightly
smoothed. The bipolar nature of the launching is evident, as seen in both the
radial velocity map and the thermodynamic projection (also lightly smoothed).
Large sound waves are emitted during these events, which appear as edges in
the pressure projection, and clearly reach at least several tenths of
$R_{500}$ ($100$s of kpc) before being significantly damped. The fact that we
see these features in our simulations is encouraging, because sound waves are
thought to be a significant mechanism for transporting energy (possibly up to
$25\%$ of the heating power; Bambic & Reynolds 2019) out of the group or
cluster core, though this is a topic of current debate and investigation
(Fabian et al., 2017; Bambic & Reynolds, 2019). Sound waves have been observed
in the cluster environment, notably in the Perseus cluster (Fabian et al.,
2003; Sanders & Fabian, 2007; Graham et al., 2008). It is left to future work
to quantify the dissipation of energy into the hot gas by compressive acoustic
waves in our suite of simulations.

Strong, bipolar cavities have been observed in several known clusters. For
example, in RBS 797 (Schindler et al., 2001; Doria et al., 2012), Abell 3017
(Parekh et al., 2017; Pandge et al., 2021), and also in the Perseus cluster
(Boehringer et al., 1993; Fabian et al., 2000, 2006). It is yet to be shown
whether such features are realistically reproducible in cosmological
simulations which do not employ a bipolar feedback model, although some work
has been done on this in the very recent past. For example, Truong et al.
(2023) examine Perseus-like clusters in the tng-cluster suite demonstrate that
bubbles can be formed via their isotropic feedback model in which the wind
direction at each event is random (Weinberger et al., 2017), although it is to
be seen whether these quantitatively match observed cavities.

In Figs. 11, 12, and 13 we show that even halos in the lowest mass bin show
clear x-shaped emission structure due to cavities displacing fluid as they
grow and rise. This suggests that groups in this respect are similar to
clusters, and we should expect to observe similar "wing"-like excess emission
around the centres of groups corresponding to the opening of cavities. Such
features, on scales of a few tens of kiloparsecs, are ubiquitous in
observations of cluster-scale halos, for example in Abell 133 (Fujita et al.,
2002, 2004; Randall et al., 2010), and in SDSS 1531 where they may coincide
with increased star formation (Omoruyi et al., 2024).

In Figs 13 and 14 we see examples of "hotspots" in temperature occurring at
the apexes of each jet lobe. This matches observations of, for example, Cygnus
A (Steenbrugge et al., 2008; Snios et al., 2020), and is due to the gas
shocking at the head of the jet (Smith et al., 1985; Smith & Donohoe, 2021).

Finally in Fig. 15, we show a time-series for the radial velocity evolution
spanning $380$ Myr, where the correspondence between outflows and inflated
bubbles which become visible in the X-ray maps. Hot "channels" are carved in
the IGrM over several tens of Myr before the jet axis re-orients, at which
point the jet starts to work against pristine halo gas, causing an X-ray
bright hotspot to appear at the site of the shock. Bubbles tend to live for
several tens of Myr and become dimmer as they migrate radially.

(a)

(b)

Figure 10: Zoom252_halo0 ($z_{obs}=0.032$): Strong bipolar jets are seen in
the temperature and density maps. Inset are mock Chandra (left) and LEM
(right) observations of the central regions, with X-ray cavities identified in
green dashed line.

(a)

(b)

Figure 11: Zoom825_halo0 ($z_{obs}=0.026$): A mid-range mass halo in our
sample, with an initial mass of $M_{500}=10^{13.49}M_{\odot}$. The cavity to
the south is clearest, and cooler, yet interestingly these two cavities were
inflated by jets launching at about the same time. The south cavity is larger
whereas the northern cavity is more compact and contains hotter gas. The
southern cavity is especially clear in both Chandra and LEM flux maps. This
group has a complex multi-phase structure, with filaments of cold gas
extending into the core, and cold blobs peppering the IGrM.

(a)

(b)

Figure 12: Zoom433_halo0 ($z_{obs}=0.0237$): A group with mass
$10^{13.39}M_{\odot}$, which shows the classic x-shaped region corresponding
to jet feedback as has been observed in real systems. This corresponds to
"wings" of emission which bound the cavities.

(a)

(b)

Figure 13: Zoom284_halo4 ($z_{obs}=0.0179$): One of the lowest-mass halos we
have in our sample at $M_{500}=10^{13.03}M_{\odot}$, this group still produces
clear bipolar X-ray cavities characterised by a dip in the density and a
boosted plasma temperature. Wing-like structures are still observable in the
mock X-ray maps. Figure 14: We show how the Simba model can re-produce hot-
spots and jet-associated X-ray excess that has been seen in observations, for
instance in Cygnus A (Duffy et al., 2018), as well as on smaller galactic
scales (Marshall et al., 2001; Hardcastle et al., 2016).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

(s)

(t)

Figure 15: We show the time-evolution of the radial velocity of all gas in
Zoom252_halo0 ($z_{obs}=0.032$) to demonstrate how a strongly-jetting system
evolves with time in Hyenas. In the earlier snapshots, the jet is North-South
and inflates two large bubbles which grow and rise with time. Sound waves in
the IGrM can also be clearly seen emanating for several $10$s of kpc from the
central regions, propagating perpendicular to the jet. After $~{}200$ Myr the
jet axis flips and the bubbles inflate along the East-West direction. The jets
themselves are sufficiently hot that the X-ray emission is maintained such
that in several cases we only observe cavities slightly adjacent to them, and
observe the jets themselves in the X-ray (for instance at a simulation time of
$10848$ Myr). This matches what is seen occasionally in real systems, for
instance in Cygnus A (Steenbrugge et al., 2008; Snios et al., 2020). After the
initial burst, the jet cools and cavities become apparent along the jet axes.

##  5 Bubbles Suppressing Cooling Flows

The relationship between the central cooling luminosity and the cavity power
is important because it indicates whether AGN feedback can sufficiently heat
the ICM and prevent large-scale cooling flows. If the cavity power is of order
the energy loss rate due to cooling, this is strong evidence towards jet
feedback being able to supply sufficient energy (though the exact method of
transferring this energy to the hot atmosphere is still poorly understood).

Observations have shown that in a large fraction of observed systems the
energy contained within hot cavities is sufficient to offset radiative losses
(Bîrzan et al., 2004; Hlavacek-Larrondo et al., 2015; Olivares et al., 2022).
Heating on sufficient or near-sufficient levels has also been reproduced in
idealised simulations of jet-mode feedback within cluster atmospheres,
suppressing large-scale cooling flows (Sijacki & Springel, 2006; Sijacki et
al., 2007; Puchwein et al., 2008; Cielo et al., 2018). It has yet to be
determined if AGN jets in cosmological simulations can effectively heat the
IGrM, given that their energies are tuned to reproducing larger-scale trends
in galaxy quenching and the galaxy stellar mass function, rather than to
reproduce hot gas properties in halos. In this Section, we whether this is
possible in the Simba model.

###  5.1 Calculating the Cooling Luminosity

(a)  Figure 16: The best fit to the $0.5-2$keV (observed-frame) spectrum
extracted from the central $7.7$ Gyr cooling radius, as observed by LEM with
$2$ eV spectral resolution for one of the larger halos in our sample
Zoom266_halo0 (snapshot 172). Posterior samples of the model parameter vector
are shown in red and the maximum-likelihood solution is shown in blue. (a)
Figure 17: The same as the above figure, but for the corresponding Chandra
observation.

First, we calculate a central cooling luminosity and plot this against cavity
power, to test the energy balance between heating and cooling. The central
cooling luminosity is defined as the luminosity within a radius
$r_{\text{cool}}$, where this radius is selected such that the cooling
timescale $t_{\text{cool}}$ within $r_{\text{cool}}$ is equal to $7.7$ Gyrs
(Macconi et al., 2022). We assume that the halo is approximately in
thermodynamic balance, since the luminosity of an actively cooling gas
component is typically fairly small compared to the total luminosity (
$\lesssim 20\%$ and frequently $<10\%$) (Rafferty et al., 2006), meaning that
$L_{\text{cool}}$ is the luminosity which must be matched by mechanical
heating. We calculate $r_{\text{cool}}$ from the emission-weighted radial
profiles, through the cooling timescale

$t_{\text{cool}}(r)=\frac{E_{\text{th}}(<r)}{L_{\text{X,bol}}(<r)},$ (15)

where $E_{\text{th}}(<r)$ is the sum of the internal energy (assuming ideal
gas with mean molecular weight $\mu=0.62$) as well as the work done to the gas
as it cools isobarically (Hlavacek-Larrondo et al., 2015). We calculate
$E_{\text{th}}(<r)$ with the ISM-free filtered field within radius $r$, and
$L_{X,bol}(<r)$ is the total bolometric ($0.1-100$ keV) luminosity within this
radius. We then select $r_{\text{cool}}$ such that
$t_{\text{cool}}(r_{\text{cool}})=7.7$ Gyrs.

We use two separate methods to calculate the central bolometric cooling
luminosity. First, we define $L_{\text{cool}}\equiv
L_{X,bol}(<r_{\text{cool}})$. Second, to better match the results and spread
of real halos, we calculate an observational bolometric luminosity with our
mock Chandra and LEM events, $L_{\text{cool}}^{obs}$. To do this we extract
the spectra within $r_{\text{cool}}$ using soxs, and fit with the Bayesian
X-ray parameter estimation code bxa (Buchner, 2016a), which connects the
spectral-fitting software (Py)XSPEC (Arnaud, 1996) to the nested-sampling code
Ultranest (Buchner, 2016b, 2019, 2021). We fit the spectra in the observed-
frame energy range $0.5-2.0$ keV with a bvapec * tbabs absorbed collisional
ionization equilibrium (CIE) emission model, using the Cash statistic as is
appropriate for spectra which may have low bin-counts. We choose
${kT,n_{e},O,Ne,Si,S,Fe}$ and the velocity as free parameters for the fit, and
freeze the other quantities to their known values for the Milky-Way foreground
column density and the redshift, and to $0.3Z_{\odot}$ with solar abundances
for the remaining elemental parameters. Using bxa, we generate a chain of
posterior parameter samples for the free parameters.

For each sample in the chain, we generate a bolometric luminosity for the
intrinsic, unabsorbed bvapec model for the halo, using the XSPEC commands
setEnergies and calcLumin. Therefore, we generate a posterior probability
distribution of bolometric luminosities, from which we calculate the median as
the maximum-likelihood value, and the $16^{th}$ and $84^{th}$ percentiles as
the $1-\sigma$ error.

It is important as a check to compare the summed intrinsic luminosity, as
calculated from the particle properties, with that obtained via a spectral fit
to the subsequently projected and observed photons. Since jet feedback has a
significant impact on the thermodynamic and emission properties of halos,
features like bubbles could compromise the accuracy of a finite (thermodynamic
and metal-) component fit. For example, the inflation of cavities can cause
asymmetrical heating within the lobes themselves and also compression of the
ambient atmosphere at the location of shock fronts, increasing the local
density and luminosity. In that case, there would be at least three major
components to the gas; hot bubble fluid, compressed gas in and around the
cavity shell, and the undisturbed halo atmosphere. In Figs 16 and 17 we show
example spectral fits for both LEM and Chandra respectively and in Fig. 18 we
plot the ground-truth bolometric luminosity as summed from the particles in
our datasets to the luminosities obtained via spectral fits using both of
these instruments. We see that even for our sample of jet-disrupted halos the
single-component fit recovers the true intrinsic source luminosity. LEM
luminosities tend to slightly overpredict, though this may be a coarse-grain
effect due to the cutting-out of a circular aperture with larger pixels
compared to Chandra. Alternatively, this may be a projection effect due to
added light from larger radii than $r_{\text{cool}}$, and in fact LEM is
capturing this better than Chandra. We proceed in our analysis mainly making
use of the intrinsic luminosity quantity with the knowledge that this value
will be essentially identical to that gained via a full observation and
spectral fit.

###  5.2 Does Cavity Power Offset Cooling?

Figure 18: We plot the intrinsic luminosity against the fit-obtained
luminosity from an bvapec*tbabs model applied to both Chandra and LEM within
$r_{\text{cool}}$. We plot only fits with $0.5<\chi^{2}/\text{d.o.f.}<5$. Even
for jet-disrupted halos the single-component fit recovers well the true
intrinsic source luminosity.

In Fig. 19, we show the observed cavity power against the luminosity within
the cooling radius, and compare to the observational results of Rafferty et
al. (2006) and Nulsen et al. (2009). To ensure a fair comparison with the
Rafferty et al. (2006) sample, we take their total cavity power for each halo,
and plot against their total luminosity within the cooling radius. For the
data of Cavagnolo et al. (2010), O’Sullivan et al. (2011), O’Sullivan et al.
(2017), O’Sullivan et al. (2018), and O’Sullivan et al. (2019), we also sum
individual cavity powers where appropriate, and plot against the cooling
luminosity (this data is compiled in Eckert et al. (2021)).

In all we have a reasonably good match to observational cavity powers. Both
the Hyenas and observational data of Nulsen et al. (2009), Cavagnolo et al.
(2010), O’Sullivan et al. (2011), O’Sullivan et al. (2017), O’Sullivan et al.
(2018), and O’Sullivan et al. (2019) appear too shallow below $L_{X,bol}\sim
10^{43}$ erg/s and, therefore, lie above the one-to-one heating/cooling
balance. Hence we find that halos in the groups regime are overheated. In this
scenario, the feedback from the core results in a net heating which over-
pressurises the halo, removing cooler gas and reducing the baryon fraction,
and quenching star formation (McNamara et al., 2008). Eckert et al. (2021)
showed that Simba matches the hot gas fraction in the group regime, and it is
encouraging that we find further evidence that the feedback model is matching
real halos (see also Robson & Davé, 2020).

From Fig. 19, we see that at the high-luminosity end the Hyenas cavities may
be under-luminous, which in this scenario could explain why the halo gas
fraction for high-mass groups and clusters in Simba is too high compared to
observational constraints; there is insufficient feedback in high-mass, high-
luminosity systems. However, we are limited by the small number of halos that
we have at these high luminosities and we have several examples at cavity
powers near $10^{44}$ erg/s, agreeing with the data and lying on or close to
the equality line where heating is more or less exactly balanced by cooling.
Given our Hyenas sample we cannot robustly tell whether, had we simulated more
massive halos, whether they would lie below the 1:1 line or whether they would
veer off and lie on the line as the observations do.

Overall, the Simba feedback model reproduces the $\sim 2-3$ orders of
magnitude dynamic range in cavity powers in galaxy groups, roughly from
$10^{41}-10^{44}$ erg/s, and matches the observed trend with cooling
luminosity. This indicates that the correct coupling of the mechanical and
X-ray luminosities is achieved through the Simba model, with heating from the
central AGN via jets providing a substantial fraction, if not the majority, of
the required energy to the group gas and preventing large-scale cooling flows.
In the low-mass, low-luminosity range, it is likely halos undergo a net
heating due to AGN activity. It is left to future work to conduct a similar
analysis in the clusters regime, but the highest luminosity groups we do have
in this sample have heating and cooling in balance, much like the
observational data.

Figure 19: We plot our Hyenas cavity powers (the sum of all cavities seen for
each snapshot) in black against the bolometric X-ray luminosity within the
cooling radius, and compare to the data of Rafferty et al. (2006), Nulsen et
al. (2009), Cavagnolo et al. (2010), O’Sullivan et al. (2011), O’Sullivan et
al. (2017), O’Sullivan et al. (2018), and O’Sullivan et al. (2019). Where
appropriate we have combined individual cavity powers and their errors for a
fair comparison to the Hyenas points. The Hyenas values are using the midpoint
pressure. We find very good agreement with real systems, with Hyenas able to
sufficiently heat halos at all luminosities, and apparently overheat at low
luminosities ($\lesssim 10^{43}$ erg/s) in agreement with observations.

###  5.3 Relationship Between Cavity Power and Accretion Rate

To understand how the cavities are powered and how this is related to the
dominant mode in which the gas is accreting onto the cental black hole, we
consider the relationship between the mass accretion rates (due to both BHL
and torque-limited accretion) and the cavity power inferred from our mock-
observation routine.

Fig. 20 shows the Bondi accretion power against the cavity power, which is a
correlation of interest since it is readily available to observers. There is a
large scatter in the Hyenas relationship, with no obvious correlation. This is
unlike the strong linear trend found in Allen et al. (2006), but we do see
overlap with the results for individual systems of Russell et al. (2013).

The overlap at the high-power end is encouraging, showing that Hyenas can
reproduce the observed range of accretion powers. However, we also predict
that there should be a large population of Bondi-dominated systems at low
accretion rates hosting powerful cavities. Such low Bondi powers are not seen
in observations, which could owe to observational sensitivity limits.

It may be that those low-accretion cavities are powerful enough to disrupt the
hot atmospheres and therefore the supply of hot gas to the central black hole.
However, it could also owe to the large stochasticity in Bondi accretion in
the Simba model (Thomas et al., 2019), as it is quite sensitive to the number
of hot gas fluid elements within the BH kernel which can vary strongly in hot
halos due to low resolution. Thomas et al. (2021) mitigated this effect by
smoothing the accretion rate over a fixed (short) timescale. In our case,
however, we have a timescale for the jet event, namely the same timescale used
for the Direct bubble lifetime as introduced in §3. For our purposes this
provides a more physically meaningful smoothing scale. We did however check
and find that our results are qualitatively the same if we instead used the
total accretion power averaged over a narrow timescale at the observation
snapshot (i.e. a value obtainable by observers).

In the left-hand plot of Fig. 21, we show the total accretion power (averaged
over the aforementioned jetting timescale) plotted against the observed cavity
power. We indicate the relative contributions of both Bondi and torque-limited
accretion rates on the plot.

Almost all cavities lie above the line at which the AGN can fuel them with
recent accretion activity. At the low power range, there are a few anomalies
in each plot that appear below the line $P_{\text{cav}}=P_{\text{accretion}}$.
In the right-hand plot of Fig. 21, we show the same data, except that for the
cavity enthalpy we use the enthalpy derived from the pressure at the bubble’s
outermost tip, instead of at the midpoint. We obtain an almost-perfect
partition of points into the upper-half plane, further corroborating our claim
made in §3.4 — that this is perhaps the more useful estimate of the bubble
enthalpy. Although artificially extending the timescale over which the AGN
winds energy is summed can bring the AGN energy in-line with the cavity
enthalpy in §3.4 (Fig. 8), this extension of the time-range considered does
not bring the accretion power in-line here, indicating that indeed it is the
cavity enthalpy that must be altered (through use of the bubble tip pressure)
and not the estimate of the timescale over which the jets are actively
inflating the cavity. With this estimate (cavity tip pressure), the accretion
rate is able to fully explain the rate of energy transfer to the bubble and
eventually to the IGrM. Where the bubble appears underpowered in these plots
indicates non-adiabatic losses, for instance from shocks, or from other
feedback modes being active as well as the jets.

Looking now at which accretion mode is most active when the bubbles are
produced, we see that the very highest-power jet-events ($P_{\dot{M}}\gtrsim
10^{44.5}$ ergs/s) are powered by Bondi accretion, but the remainder are
powered by a mix of Bondi and torque-limited accretion, with no apparent
relation to cavity power. There are many cavities inflated when the black hole
is powered predominantly by torque-limited accretion ($64\%$). This implies
that caution must be used in the groups regime when estimating accretion
rates; there may be a significant fraction of the accretion provided by the
disk, raising the total accretion rate above that of the Bondi rate and adding
to the powers of cavities. Cavities which appear significantly under-powered
based on their accretion rates are likely to be in a highly active quasar
state, with jets contributing less to the total power output at high Eddington
fractions.

Overall, we find that the simulated Hyenas groups display jet-driven bubbles
that are in good quantitative agreement with observations of bubbles in
similar-mass systems. This is somewhat remarkable given the crude nature of
the implementation of jets in the Simba model, necessitated by poor
resolution. The bubbles identified in mock X-ray maps using tools analogous to
those used by observers show similar cavity powers, and are generally
sufficient to counteract cooling in the halo. In fact, in lower-mass groups
they can significantly exceed cooling losses, which explains why Simba’s
groups are strongly evacuated of baryons (Appleby et al., 2021; Sorini et al.,
2022). The power injected by the jet exceeds that seen in the cavity, as
expected based on there being some radiative losses and/or energy deposition
outside the cavity, and this is only exposed in Hyenas if one measures the
pressure at the bubble’s leading tip rather than its midpoint.

Figure 20: We plot Bondi accretion powers against observed cavity power (using
$t_{buoy}$) and compare to the derived observational log-linear relation of
Allen et al. (2006) and the datapoints of Russell et al. (2013).

(a)

(b)

Figure 21: Left: We plot the power associated with the total accretion rate
using the Simba radiative efficiency of $\eta=0.1$ against the observed cavity
powers obtained with the midpoint pressure. The accretion power is calculated
as the mean over the jetting period, although the results are insensitive to
instead using the instantaneous accretion powers at the observation snapshot.
Most (but not all) cavity powers appear accounted for by the recent AGN
accretion. Right: As in the left-hand plot but now with the enthalpy
calculated from the pressure at the bubble’s outermost tip. Again we see a
more encouraging picture with this enthalpy value, with all cavities able to
be accounted for from the recent accretion rate

##  6 Conclusions

In this work, we present the first quantitative investigation of X-ray
cavities in a large-scale cosmological model suite of simulations, calibrated
solely to galaxy properties, that properly straddles simulation and
observation. By taking an observer-centric approach to analysing our Hyenas
suite of galaxy group zooms, we are able to make predictions and comparisons
realistically in the context of existing X-ray catalogues. Using mock X-ray
maps generated from the PyXSIM and SOXS software packages, we fully mimic the
analyses carried out by observers on real data. We use model-subtracted
residual maps, unsharp-masked images, and smoothed X-ray flux maps to identify
and measure cavity properties. In addition, we use the CADET machine-learning
cavity-identifying package and find that it is broadly very useful for
identifying potential cavities for follow-up. We create spectra within
$r_{\text{cool}}$ for each system using SOXS, and measure bolometric
luminosities through a Bayesian fitting approach. Our overarching conclusion
is that we demonstrate, for the first time within a fully cosmological
setting, how jets from black hole activity can inflate bubbles that deposit
sufficient energy in their environments to offset radiative cooling, thereby
suppressing cooling flows to yield quenched galaxies. This in a sense
"completes the loop" between galaxies, their black holes, and their circum-
galactic media in order to allow holistic studies of the process of galaxy
transformation.

Based on our analysis in this work, we draw the following conclusions:

  * •

Hyenas successfully produces many X-ray cavities. These coincide with
depressions in the density maps and with high outgoing velocities in the
radial velocity projection maps. Comparing to the sample of Rafferty et al.
(2006) and Shin et al. (2016), we find that the dimensions of the Hyenas
cavities in terms of semi-major axis, surface area, and eccentricity in
relation to the bubble midpoint distance match very closely both the trend and
scatter of observed cavities. Furthermore our cavity timescales also agree
well with observations, with an average buoyancy timescale of
$t_{\text{buoy}}=49.2\pm 19.5$ Myr.

  * •

We present King model-subtracted residual maps and unsharp-masked maps in §3
and also case-studies of a small number of selected halos in §4. We
demonstrate that the Simba model is in the very least qualitatively accurate
in reproducing the observed cavity populations of real systems. Through
comparison with the pressure maps and radial velocity maps obtained directly
from the particle data, we find that cavities are associated with pressure-
fronts, and large radial velocities relative to the ambient medium. We also
observe large sound waves emanating from the central AGN, similar to those
seen in the Perseus cluster. We leave it to future work to quantify these
various features, and most importantly the relative importance of each
mechanism in transferring energy from the AGN to the hot atmosphere.

  * •

Our bubble enthalpies (for a given snapshot) tend to appear over-powered
compared to the energy injected by recent AGN feedback. This is remedied
remarkably well if we instead simply calculate an enthalpy using the pressure
at the bubble’s outermost tip for the $PdV$ work estimate, hinting that this
may give a closer estimate of the "effective" pressure of the bubble over its
expansion. This prediction that the bubble should do work preferentially along
the ’axis of least resistance’ is supported by our observation of high bubble
eccentricities (of which $65\%$ have the semi-major axis aligned radially) in
agreement with observations. Our simulations show excellent agreement with the
observed data when we use the bubble midpoint pressure in terms of the
$P_{cav}=L_{X}$, indicating that this result (that the bubble tip pressure is
actually more accurate) should generalise to estimating the energy of the AGN
in real systems.

  * •

We compute the observed cavity powers and compare to the bolometric X-ray
luminosity within the $7.7$ Gyr cooling radius. We find good agreement at the
higher luminosity end with the observational data of Rafferty et al. (2006)
and excellent agreement at the lower-$L_{X}$ end with that of Nulsen et al.
(2009), Cavagnolo et al. (2010), O’Sullivan et al. (2011), O’Sullivan et al.
(2017), O’Sullivan et al. (2018), and O’Sullivan et al. (2019). This relation
lies on the line $P_{cav}=L_{X}$, except at lower luminosities
($L_{X,bol}\lesssim 10^{43}$ erg/s) where both the data and our results lie to
higher cavity powers. Therefore, we find that the Simba feedback model is
sufficiently powerful in terms of energy output to balance cooling in the IGrM
and prevent a large-scale cooling-flow, and potentially overheats the halo gas
at low halo masses. This may be the reason why Simba tends to have a lower hot
gas fraction in the groups regime compared to other simulations, which is in
line with observational data (Robson & Davé, 2020; Eckert et al., 2021).

  * •

Plotting the Bondi power against the cavity power we observe a very large
scatter around the one-to-one relationship. This is unlike the strong
correlation found in Allen et al. (2006) but similar to the results in Russell
et al. (2013). However, we predict a large population of powerful cavities in
systems with low Bondi accretion rates, hinting at a bias in observations, and
suggesting alternate modes of accretion onto the SMBH such as described in
Anglés-Alcázar et al. (2017); Prasad et al. (2017). In Hyenas we find that
Bondi accretion dominates in the very most massive systems, but outside of
those there are roughly comparable numbers of torque-dominated and Bondi-
dominated jet-driving black holes.

We demonstrate the feasibility of a large survey of X-ray halos in a
cosmological suite of simulations, and how they should be analysed in order to
properly compare with catalogues of real systems. We emphasise the importance
of this kind of analysis, and the need for feedback models to be tested on
small (group or cluster) scales and not just compared against summary galaxy
population properties. We envision the simulated cavity population can act as
a useful calibration barometer for large-scale simulations in the future,
especially as simulation resolution increases, and as the real cavity
population is better understood with upcoming X-ray missions. The main barrier
to this is the time required to manually identify and measure cavities.
However, this can at least partly be automated with machine-learning
pipelines, such as cadet.

## Software

This project made use of the following external software packages for
analysis; pyXSIM 4.3.0 (ZuHone & Hallman, 2016), SOXS 4.6.0 (ZuHone et al.,
2023), yT Project 4.2.2 (Turk et al., 2011), BXA 4.1.1 (Buchner, 2016a) and
UltraNest 3.6.2 (Buchner, 2021), Scipy 1.11.2 (Virtanen et al., 2020), Numpy
1.24.3 (Harris et al., 2020), mpi4py 3.1.4 (Dalcin et al., 2011), CAESAR 0.2b0
(Thompson, 2015), PyGadgetReader 2.6 (Thompson, 2014), AstroPy 5.3.3 (Astropy
Collaboration et al., 2013), Corner 2.2.2 (Foreman-Mackey, 2016), (Py)CADET
0.1.4 (Plšek et al., 2024), Matplotlib 3.7.3 (Hunter, 2007), smplotlib 0.0.9
(Li, 2023), cmasher 1.6.3 (van der Velden, 2021), pandas 2.1.0 (The pandas
development Team, 2023), h5py 3.9.0 (Collette et al., 2022), XSPEC 12.13.1
(Arnaud et al., 1999), and SAOImageDS9 8.4.1 (Smithsonian Astrophysical
Observatory, 2000). We are grateful to their respective authors for making
them available.

## Acknowledgements

FJJ would like to acknowledge the support of the Science and Technology
Facilities Council. He would also like to thank Ewan O’Sullivan and Nhut
Truong for useful conversations, Britton Smith for advice regarding yT and
John ZuHone for advice regarding the PyXSIM and SOXS software packages, as
well as for making these packages publicly available.

AB acknowledges support from the Natural Sciences and Engineering Research
Council of Canada (NSERC) through its Discovery Grant program. He also
acknowledges support from the Infosys Foundation via an endowed Infosys
Visiting Chair Professorship at the Indian Institute of Science, and from the
Leverhulme Trust via a Visiting Professorship at the Institute for Astronomy,
University of Edinburgh. RD is supported by the STFC RC grant ST/Y001117/1. WC
is supported by the STFC AGP Grant ST/V000594/1, the Atracción de Talento
Contract no. 2020-T1/TIC-19882 was granted by the Comunidad de Madrid in
Spain, and the science research grants were from the China Manned Space
Project. He also thanks the Ministerio de Ciencia e Innovación (Spain) for
financial support under Project grant PID2021-122603NB-C21 and HORIZON EUROPE
Marie Sklodowska-Curie Actions for supporting the LACEGAL-III project with
grant number 101086388. DR is supported by the Simons Foundation.

These original Hyenas simulations were performed using the DiRAC@Durham
facility managed by the Institute for Computational Cosmology on behalf of the
STFC DiRAC HPC Facility (www.dirac.ac.uk) with the DiRAC Project: ACSP252,
titled Simba Zoom Simulations of Galaxy Groups. The equipment was funded by
BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and
ST/S002502/1, Durham University, and STFC operations grant ST/R000832/1. DiRAC
is part of the National e-Infrastructure. Computations and analysis specific
to this work were made possible with support from SciNet and the Niagara
supercomputing cluster, operated by the Digital Research Alliance of Canada
(alliancecan.ca).

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## Data Availability

The relevant data will be shared upon reasonable request to the corresponding
author.

## References

  * Allen et al. (2006) Allen S. W., Dunn R. J. H., Fabian A. C., Taylor G. B., Reynolds C. S., 2006, MNRAS, 372, 21
  * Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840
  * Appleby et al. (2021) Appleby S., Davé R., Sorini D., Storey-Fisher K., Smith B., 2021, MNRAS, 507, 2383
  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Arnaud et al. (1999) Arnaud K., Dorman B., Gordon C., 1999, XSPEC: An X-ray spectral fitting package, Astrophysics Source Code Library, record ascl:9910.005
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Babul et al. (2002) Babul A., Balogh M. L., Lewis G. F., Poole G. B., 2002, MNRAS, 330, 329
  * Babul et al. (2013) Babul A., Sharma P., Reynolds C. S., 2013, ApJ, 768, 11
  * Bambic & Reynolds (2019) Bambic C. J., Reynolds C. S., 2019, ApJ, 886, 78
  * Beckmann et al. (2017) Beckmann R. S., et al., 2017, MNRAS, 472, 949
  * Best & Heckman (2012) Best P. N., Heckman T. M., 2012, MNRAS, 421, 1569
  * Bîrzan et al. (2004) Bîrzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen P. E. J., 2004, ApJ, 607, 800
  * Boehringer et al. (1993) Boehringer H., Voges W., Fabian A. C., Edge A. C., Neumann D. M., 1993, MNRAS, 264, L25
  * Borrow et al. (2020) Borrow J., Anglés-Alcázar D., Davé R., 2020, MNRAS, 491, 6102
  * Bradley et al. (2022) Bradley L., Davé R., Cui W., Smith B., Sorini D., 2022, arXiv e-prints, p. arXiv:2203.15055
  * Brüggen et al. (2005) Brüggen M., Ruszkowski M., Hallman E., 2005, ApJ, 630, 740
  * Buchner (2016a) Buchner J., 2016a, BXA: Bayesian X-ray Analysis, Astrophysics Source Code Library, record ascl:1610.011
  * Buchner (2016b) Buchner J., 2016b, Statistics and Computing, 26, 383
  * Buchner (2019) Buchner J., 2019, PASP, 131, 108005
  * Buchner (2021) Buchner J., 2021, The Journal of Open Source Software, 6, 3001
  * Cavagnolo et al. (2010) Cavagnolo K. W., McNamara B. R., Nulsen P. E. J., Carilli C. L., Jones C., Bîrzan L., 2010, ApJ, 720, 1066
  * Choi et al. (2012) Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125
  * Churazov et al. (2001) Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, ApJ, 554, 261
  * Cielo et al. (2018) Cielo S., Babul A., Antonuccio-Delogu V., Silk J., Volonteri M., 2018, A&A, 617, A58
  * Coenda et al. (2018) Coenda V., Martínez H. J., Muriel H., 2018, MNRAS, 473, 5617
  * Collette et al. (2022) Collette A., et al., 2022, h5py/h5py: 3.7.0, doi:10.5281/zenodo.6575970
  * Crockett et al. (2012) Crockett R. M., et al., 2012, MNRAS, 421, 1603
  * Cui et al. (2021) Cui W., Davé R., Peacock J. A., Anglés-Alcázar D., Yang X., 2021, Nature Astronomy, 5, 1069
  * Cui et al. (2024) Cui W., Jennings F., Dave R., Babul A., Gozaliasl G., 2024, arXiv e-prints, p. arXiv:2406.17012
  * Dalcin et al. (2011) Dalcin L. D., Paz R. R., Kler P. A., Cosimo A., 2011, Advances in Water Resources, 34, 1124
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Doria et al. (2012) Doria A., Gitti M., Ettori S., Brighenti F., Nulsen P. E. J., McNamara B. R., 2012, ApJ, 753, 47
  * Dubois et al. (2010) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409, 985
  * Dubois et al. (2014) Dubois Y., et al., 2014, MNRAS, 444, 1453
  * Dubois et al. (2016) Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948
  * Duffy et al. (2018) Duffy R. T., et al., 2018, MNRAS, 476, 4848
  * Eckert et al. (2021) Eckert D., Gaspari M., Gastaldello F., Le Brun A. M. C., O’Sullivan E., 2021, Universe, 7, 142
  * Einasto et al. (2024) Einasto M., Einasto J., Tenjes P., Korhonen S., Kipper R., Tempel E., Liivamägi L. J., Heinämäki P., 2024, A&A, 681, A91
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Fabian et al. (2000) Fabian A. C., et al., 2000, MNRAS, 318, L65
  * Fabian et al. (2003) Fabian A. C., Sanders J. S., Allen S. W., Crawford C. S., Iwasawa K., Johnstone R. M., Schmidt R. W., Taylor G. B., 2003, MNRAS, 344, L43
  * Fabian et al. (2005) Fabian A. C., Reynolds C. S., Taylor G. B., Dunn R. J. H., 2005, MNRAS, 363, 891
  * Fabian et al. (2006) Fabian A. C., Sanders J. S., Taylor G. B., Allen S. W., Crawford C. S., Johnstone R. M., Iwasawa K., 2006, MNRAS, 366, 417
  * Fabian et al. (2017) Fabian A. C., Walker S. A., Russell H. R., Pinto C., Sanders J. S., Reynolds C. S., 2017, MNRAS, 464, L1
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24
  * Fujita et al. (2002) Fujita Y., Sarazin C. L., Kempner J. C., Rudnick L., Slee O. B., Roy A. L., Andernach H., Ehle M., 2002, ApJ, 575, 764
  * Fujita et al. (2004) Fujita Y., Sarazin C. L., Reiprich T. H., Andernach H., Ehle M., Murgia M., Rudnick L., Slee O. B., 2004, ApJ, 616, 157
  * Gaspari et al. (2012) Gaspari M., Ruszkowski M., Sharma P., 2012, ApJ, 746, 94
  * Gaspari et al. (2013) Gaspari M., Brighenti F., Ruszkowski M., 2013, Astronomische Nachrichten, 334, 394
  * Gozaliasl et al. (2016) Gozaliasl G., Finoguenov A., Khosroshahi H. G., Mirkazemi M., Erfanianfar G., Tanaka M., 2016, MNRAS, 458, 2762
  * Graham et al. (2008) Graham J., Fabian A. C., Sanders J. S., 2008, MNRAS, 391, 1749
  * Hardcastle et al. (2016) Hardcastle M. J., et al., 2016, MNRAS, 455, 3526
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Heckman & Best (2014) Heckman T. M., Best P. N., 2014, ARA&A, 52, 589
  * Hlavacek-Larrondo et al. (2015) Hlavacek-Larrondo J., et al., 2015, ApJ, 805, 35
  * Hopkins & Quataert (2010) Hopkins P. F., Quataert E., 2010, MNRAS, 407, 1529
  * Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Huško et al. (2022) Huško F., Lacey C. G., Schaye J., Schaller M., Nobels F. S. J., 2022, MNRAS, 516, 3750
  * Jennings & Davé (2023) Jennings F., Davé R., 2023, arXiv e-prints, p. arXiv:2306.01397
  * Kaviraj et al. (2017) Kaviraj S., et al., 2017, MNRAS, 467, 4739
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kraft et al. (2012) Kraft R. P., et al., 2012, ApJ, 749, 19
  * Le Brun et al. (2014) Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, MNRAS, 441, 1270
  * Li (2023) Li J., 2023, AstroJacobLi/smplotlib: v0.0.9, doi:10.5281/zenodo.8126529
  * Macconi et al. (2022) Macconi D., Grandi P., Gitti M., Vignali C., Torresi E., Brighenti F., 2022, A&A, 660, A32
  * Marshall et al. (2001) Marshall H. L., et al., 2001, ApJ, 549, L167
  * McCabe & Borthakur (2023) McCabe T., Borthakur S., 2023, in American Astronomical Society Meeting Abstracts. p. 175.03
  * McCabe et al. (2021) McCabe T., Borthakur S., Heckman T., Tumlinson J., Bordoloi R., Dave R., 2021, ApJ, 923, 189
  * McCarthy et al. (2010) McCarthy I. G., et al., 2010, MNRAS, 406, 822
  * McDonald et al. (2018) McDonald M., Gaspari M., McNamara B. R., Tremblay G. R., 2018, ApJ, 858, 45
  * McNamara et al. (2008) McNamara B. R., Bîrzan L., Rafferty D. A., Nulsen P. E. J., Carilli C., Wise M. W., 2008, in Rector T. A., De Young D. S., eds, Astronomical Society of the Pacific Conference Series Vol. 386, Extragalactic Jets: Theory and Observation from Radio to Gamma Ray. p. 311 (arXiv:0708.0579), doi:10.48550/arXiv.0708.0579
  * Nayakshin & Zubovas (2012) Nayakshin S., Zubovas K., 2012, MNRAS, 427, 372
  * Nulsen et al. (2009) Nulsen P., Jones C., Forman W., Churazov E., McNamara B., David L., Murray S., 2009, in Heinz S., Wilcots E., eds, American Institute of Physics Conference Series Vol. 1201, The Monster’s Fiery Breath: Feedback in Galaxies, Groups, and Clusters. pp 198–201 (arXiv:0909.1809), doi:10.1063/1.3293033
  * Nusser et al. (2006) Nusser A., Silk J., Babul A., 2006, MNRAS, 373, 739
  * O’Sullivan et al. (2011) O’Sullivan E., Giacintucci S., David L. P., Gitti M., Vrtilek J. M., Raychaudhury S., Ponman T. J., 2011, ApJ, 735, 11
  * O’Sullivan et al. (2017) O’Sullivan E., et al., 2017, MNRAS, 472, 1482
  * O’Sullivan et al. (2018) O’Sullivan E., Kolokythas K., Kantharia N. G., Raychaudhury S., David L. P., Vrtilek J. M., 2018, MNRAS, 473, 5248
  * O’Sullivan et al. (2019) O’Sullivan E., Schellenberger G., Burke D. J., Sun M., Vrtilek J. M., David L. P., Sarazin C., 2019, MNRAS, 488, 2925
  * Olivares et al. (2022) Olivares V., Su Y., Nulsen P., Kraft R., Somboonpanyakul T., Andrade-Santos F., Jones C., Forman W., 2022, MNRAS, 516, L101
  * Omma et al. (2004) Omma H., Binney J., Bryan G., Slyz A., 2004, MNRAS, 348, 1105
  * Omoruyi et al. (2024) Omoruyi O., et al., 2024, ApJ, 963, 1
  * Panagoulia et al. (2014) Panagoulia E. K., Fabian A. C., Sanders J. S., Hlavacek-Larrondo J., 2014, MNRAS, 444, 1236
  * Pandge et al. (2021) Pandge M. B., Sebastian B., Seth R., Raychaudhury S., 2021, MNRAS, 504, 1644
  * Parekh et al. (2017) Parekh V., Durret F., Padmanabh P., Pandge M. B., 2017, MNRAS, 470, 3742
  * Piotrowska et al. (2022) Piotrowska J. M., Bluck A. F. L., Maiolino R., Peng Y., 2022, MNRAS, 512, 1052
  * Plucinsky et al. (2018) Plucinsky P. P., Bogdan A., Marshall H. L., Tice N. W., 2018, in den Herder J.-W. A., Nikzad S., Nakazawa K., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray. p. 106996B (arXiv:1809.02225), doi:10.1117/12.2312748
  * Plšek et al. (2024) Plšek T., Werner N., Topinka M., Simionescu A., 2024, MNRAS, 527, 3315
  * Plšek et al. (2023) Plšek T., Werner N., Topinka M., Simionescu A., 2023, CAvity DEtection Tool (CADET): Pipeline for automatic detection of X-ray cavities in hot galactic and cluster atmospheres (arXiv:2304.05457)
  * Pope et al. (2010) Pope E. C. D., Babul A., Pavlovski G., Bower R. G., Dotter A., 2010, MNRAS, 406, 2023
  * Prasad et al. (2015) Prasad D., Sharma P., Babul A., 2015, ApJ, 811, 108
  * Prasad et al. (2017) Prasad D., Sharma P., Babul A., 2017, MNRAS, 471, 1531
  * Prasad et al. (2018) Prasad D., Sharma P., Babul A., 2018, ApJ, 863, 62
  * Puchwein et al. (2008) Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53
  * Rafferty et al. (2006) Rafferty D. A., McNamara B. R., Nulsen P. E. J., Wise M. W., 2006, ApJ, 652, 216
  * Randall et al. (2010) Randall S. W., Clarke T. E., Nulsen P. E. J., Owers M. S., Sarazin C. L., Forman W. R., Murray S. S., 2010, ApJ, 722, 825
  * Raouf et al. (2019) Raouf M., Silk J., Shabala S. S., Mamon G. A., Croton D. J., Khosroshahi H. G., Beckmann R. S., 2019, MNRAS, 486, 1509
  * Rasmussen et al. (2012) Rasmussen J., Mulchaey J. S., Bai L., Ponman T. J., Raychaudhury S., Dariush A., 2012, ApJ, 757, 122
  * Rennehan et al. (2023) Rennehan D., Babul A., Moa B., Davé R., 2023, arXiv e-prints, p. arXiv:2309.15898
  * Reynolds et al. (2005) Reynolds C. S., Brenneman L. W., Stocke J. T., 2005, MNRAS, 357, 381
  * Robson & Davé (2020) Robson D., Davé R., 2020, MNRAS, 498, 3061
  * Roychowdhury et al. (2005) Roychowdhury S., Ruszkowski M., Nath B. B., 2005, ApJ, 634, 90
  * Russell et al. (2013) Russell H. R., McNamara B. R., Edge A. C., Hogan M. T., Main R. A., Vantyghem A. N., 2013, MNRAS, 432, 530
  * Ruszkowski et al. (2004) Ruszkowski M., Brüggen M., Begelman M. C., 2004, ApJ, 615, 675
  * Saeedzadeh et al. (2023) Saeedzadeh V., et al., 2023, MNRAS, 525, 5677
  * Salerno et al. (2019) Salerno J. M., Martínez H. J., Muriel H., 2019, MNRAS, 484, 2
  * Sanders & Fabian (2007) Sanders J. S., Fabian A. C., 2007, MNRAS, 381, 1381
  * Schindler et al. (2001) Schindler S., Castillo-Morales A., De Filippis E., Schwope A., Wambsganss J., 2001, A&A, 376, L27
  * Shabala et al. (2011) Shabala S. S., Kaviraj S., Silk J., 2011, MNRAS, 413, 2815
  * Sharma et al. (2012) Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420, 3174
  * Shin et al. (2016) Shin J., Woo J.-H., Mulchaey J. S., 2016, ApJS, 227, 31
  * Sijacki & Springel (2006) Sijacki D., Springel V., 2006, MNRAS, 366, 397
  * Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  * Silk (2013) Silk J., 2013, ApJ, 772, 112
  * Silk et al. (2014) Silk J., Di Cintio A., Dvorkin I., 2014, in Proceedings of the International School of Physics ’Enrico Fermi’ Course 186 ’New Horizons for Observational Cosmology’ Vol. 186. pp 137–187 (arXiv:1312.0107), doi:10.3254/978-1-61499-476-3-137
  * Simionescu et al. (2009) Simionescu A., Roediger E., Nulsen P. E. J., Brüggen M., Forman W. R., Böhringer H., Werner N., Finoguenov A., 2009, A&A, 495, 721
  * Smith & Donohoe (2021) Smith M. D., Donohoe J., 2021, MNRAS, 502, 423
  * Smith et al. (1985) Smith M. D., Norman M. L., Winkler K. H. A., Smarr L., 1985, MNRAS, 214, 67
  * Smith et al. (2001) Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91
  * Smithsonian Astrophysical Observatory (2000) Smithsonian Astrophysical Observatory 2000, SAOImage DS9: A utility for displaying astronomical images in the X11 window environment, Astrophysics Source Code Library, record ascl:0003.002
  * Snios et al. (2020) Snios B., Johnson A. C., Nulsen P. E. J., Kraft R. P., de Vries M., Perley R. A., Sebokolodi L., Wise M. W., 2020, ApJ, 891, 173
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Sonkamble et al. (2015) Sonkamble S. S., Vagshette N. D., Pawar P. K., Patil M. K., 2015, Ap&SS, 359, 21
  * Sorini et al. (2022) Sorini D., Davé R., Cui W., Appleby S., 2022, MNRAS, 516, 883
  * Steenbrugge et al. (2008) Steenbrugge K. C., Blundell K. M., Duffy P., 2008, MNRAS, 388, 1465
  * Szpila et al. (2024) Szpila J., Davé R., Rennehan D., Cui W., Hough R., 2024, arXiv e-prints, p. arXiv:2402.08729
  * The pandas development Team (2023) The pandas development Team 2023, pandas-dev/pandas: Pandas, doi:10.5281/zenodo.10107975
  * Thomas et al. (2019) Thomas N., Davé R., Anglés-Alcázar D., Jarvis M., 2019, MNRAS, 487, 5764
  * Thomas et al. (2021) Thomas N., Davé R., Jarvis M. J., Anglés-Alcázar D., 2021, MNRAS, 503, 3492
  * Thompson (2014) Thompson R., 2014, pyGadgetReader: GADGET snapshot reader for python, Astrophysics Source Code Library, record ascl:1411.001
  * Thompson (2015) Thompson R., 2015, SPHGR: Smoothed-Particle Hydrodynamics Galaxy Reduction, Astrophysics Source Code Library, record ascl:1502.012
  * Tremmel et al. (2017) Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T. R., Pontzen A., Anderson L., Bellovary J., 2017, MNRAS, 470, 1121
  * Tremmel et al. (2019) Tremmel M., et al., 2019, MNRAS, 483, 3336
  * Truong et al. (2023) Truong N., Pillepich A., Nelson D., Zhuravleva I., Lee W., Ayromlou M., Lehle K., 2023, arXiv e-prints, p. arXiv:2311.06334
  * Tully (2015) Tully R. B., 2015, AJ, 149, 54
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
  * Vagshette et al. (2017) Vagshette N. D., Naik S., Patil M. K., Sonkamble S. S., 2017, MNRAS, 466, 2054
  * Vantyghem et al. (2018) Vantyghem A. N., et al., 2018, ApJ, 863, 193
  * Vernaleo & Reynolds (2006) Vernaleo J. C., Reynolds C. S., 2006, ApJ, 645, 83
  * Vikhlinin et al. (2006) Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Van Speybroeck L., 2006, ApJ, 640, 691
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weinmann et al. (2006) Weinmann S. M., van den Bosch F. C., Yang X., Mo H. J., 2006, MNRAS, 366, 2
  * Wetzel et al. (2012) Wetzel A. R., Tinker J. L., Conroy C., 2012, MNRAS, 424, 232
  * Zheng et al. (2022) Zheng Y., Dave R., Wild V., Montero F. R., 2022, MNRAS, 513, 27
  * ZuHone & Hallman (2016) ZuHone J. A., Hallman E. J., 2016, pyXSIM: Synthetic X-ray observations generator, Astrophysics Source Code Library, record ascl:1608.002
  * ZuHone et al. (2023) ZuHone J. A., Vikhlinin A., Tremblay G. R., Randall S. W., Andrade-Santos F., Bourdin H., 2023, SOXS: Simulated Observations of X-ray Sources, Astrophysics Source Code Library, record ascl:2301.024
  * Zubovas & Bourne (2017) Zubovas K., Bourne M. A., 2017, MNRAS, 468, 4956
  * van der Velden (2021) van der Velden E., 2021, CMasher: Scientific colormaps for making accessible, informative plots, Astrophysics Source Code Library, record ascl:2102.008

  *[IGrM]: intragroup medium
  *[SMBH]: supermassive black hole
  *[AGN]: active galactic nucleus
  *[ICM]: intracluster medium
  *[BCG]: brightest cluster galaxy
  *[BHL]: Bondi-Hoyle-Lyttleton
  *[ISM]: interstellar medium
  *[CGM]: circumgalactic medium
  *[IGM]: intergalactic medium
  *[BGG]: brightest group galaxy
  *[CIE]: collisional ionization equilibrium