11institutetext: Universität Hamburg, Hamburger Sternwarte, Gojenbergsweg 112,
21029 Hamburg, Germany
11email<EMAIL_ADDRESS>22institutetext: Department of
Computational Mathematics, Science, and Engineering, Michigan State
University, East Lansing, MI 48824, USA 33institutetext: Department of Physics
and Astronomy, Michigan State University, East Lansing, MI 48824, USA
44institutetext: Facility for Rare Isotope Beams, Michigan State University,
East Lansing, MI 48824, USA

# The formation and survival of cold gas in a magnetized cool-core galaxy
cluster

M. Fournier 11 P. Grete 11 M. Brüggen 11 F. W. Glines and B. W. O’Shea
2233223344

(Received )

###### Abstract

Context. Filaments of cold gas ($T\leq 10^{4}$ K) are found in the inner
regions of many cool-core clusters. These structures are thought to play a
major role in the regulation of feedback from active galactic nuclei (AGN).

Aims. We study the morphology of the filaments, their formation, and their
impact on the propagation of the outflowing AGN jets.

Methods. We present a set of GPU-accelerated 3D (magneto)hydrodynamical
simulations of an idealized Perseus-like cluster using the performance
portable code AthenaPK. We include radiative cooling, and a self-regulated AGN
feedback model that redistributes accreted material through kinetic, thermal
and magnetic feedback.

Results. We confirm that magnetic fields play an important role in both the
formation and evolution of the cold material. These suppress the formation of
massive cold disks and favor magnetically supported filaments over clumpy
structures. Achieving resolutions of $25-50$ pc, we find that filaments are
not monolithic as they contain numerous and complex magnetically–supported
substructures. We find that the mass distribution of these clumps follows a
power law for all investigated filaments, consistent with previous cloud-
crushing simulations of individual clumps. Studying the evolution of
individual filaments, we find that their formation pathways can be diverse. We
find examples of filaments forming through a combination of cold gas uplifting
and condensation, as well as systems of purely infalling clumps condensing out
of the intracluster medium. The density contrast between the cold gas and the
outflowing hot material leads to recurring deflections of the jets, favoring
inflation of bubbles.

Conclusions. Filaments in cool–core clusters are clumpy and contain numerous
substructures, resulting from a complex interplay between magnetic fields,
thermal instability, and jet-cloud interaction. Frequent deflections of the
AGN outflows suppress jet collimation and favor the formation of large X–ray
bubbles, and smaller off-axis cavities.

###### Key Words.:

galaxies: clusters: intracluster medium – galaxies: jets – galaxies: clusters:
general – methods: numerical – magnetohydrodynamics

##  1 Introduction

Galaxy clusters are massive ($M\sim 10^{14}-10^{15}\ \rm{M}_{\odot}$)
structures resulting from the gravitational collapse of the highest
overdensities in the early Universe and their subsequent accretion of matter.
Dark matter accounts for roughly 85% of the cluster halo mass with baryons
accounting for only about 15%. Over 90% of the galaxy cluster’s baryonic mass
is in the form of a hot ($10^{7}-10^{8}$ K) diffuse plasma sitting in
hydrostatic equilibrium in the dark matter halo. (Andreon 2010; Gonzalez et
al. 2013; Laganá et al. 2013).

The ICM is expected to cool down radiatively through Bremsstrahlung,
recombination, and collisional processes (Boehringer & Hensler 1989). A third
to half of all clusters, usually designated as ”cool–core” (CC) clusters, have
a large central overdensity with an associated cooling time as small as two
orders of magnitudes below the Hubble time (Peres et al. 1998). This suggests
that in the absence of heating sources, cooling gas should flow towards the
cluster’s center, with mass rates of up to $10^{3}\,\rm{M}_{\odot}\
\rm{yr}^{-1}$ (Fabian 1994; White et al. 1997). The absence of specific
spectral signatures of such cooling flows from X-ray surveys or comparably
large star formation rates in the clusters’ central galaxies leads to the
conclusion that such gas flows are absent in most clusters (Peterson et al.
2001) and that a heating source must be counteracting the radiative cooling.
Following observations of CC clusters shed light on the presence of large
cavities apparently inflated by AGN jets (Fabian et al. 2006; Bîrzan et al.
2004). While it is generally accepted that active galactic nuclei (AGN) are
the main heating source balancing radiative cooling in galaxy clusters, the
way this energy is transferred and distributed in the surrounding ICM remains
an open question (Böhringer & Werner 2010).

Observations of various clusters suggests that the gas inside the inner tens
of kiloparsec (kpc) of clusters is multiphase, with temperatures ranging from
$10$ to $10^{9}$ K (Allen 2000; Bonamente et al. 2001, e.g.).
H$\alpha$-emitting nebulae of elongated, cold filaments have been widely
studied, and are thought to be ubiquitous in cool-core (CC) clusters (McDonald
et al. 2012; Olivares et al. 2019). Multi-wavelength surveys of cool core
clusters show that the filaments have typical lengths ranging from a few kpc
to tens of kpc, thicknesses down to $\sim 50$ pc (Fabian et al. 2008, 2016),
and are typically found below or at the edge of the X-ray cavities (Salomé et
al. 2006; Tremblay et al. 2018; Olivares et al. 2019), suggesting a tight
connection between AGN feedback and the formation of the filaments. It is
commonly assumed that the thermal instability (TI) is a key process
contributing to the formation of multiphase gas (McCourt et al. 2012; Voit et
al. 2017). Turbulence, either induced by AGN jets or by environmental effects
drives the formation of overdensities in the ICM. Due to the dependence of
radiative cooling on density, these overdensities can grow non-linearly by
drawing hot material out of the background ICM, resulting in structures of
molecular and ionized gas with electron densities $\gtrsim 100\,\rm{cm}^{-3}$
(McDonald et al. 2012). These structures are then expected to decouple from
the ICM and ”rain” towards the center of the cluster, eventually fueling the
SMBH in a choatic manner (Gaspari et al. 2013). It is also suspected that the
AGN jets may be able to uplift cold material present in the inner kpc of the
clusters up to altitudes of tens of kpc (Russell et al. 2017; Tremblay et al.
2018). This uplifted material could act as a seed for condensation and act as
a precursor for the formation of massive, infalling filaments.

Various simulation setups have been employed to test this model.
Hydrodynamical and magnetohydrodynamical (MHD) simulations of stratified or
periodic boxes (McCourt et al. 2012; Mohapatra et al. 2022) and idealised CC
clusters (Li & Bryan 2014; Beckmann et al. 2019; Qiu et al. 2019; Wang et al.
2021; Ehlert et al. 2023) have been able to reproduce some of the observed
properties of cold filaments. In Li & Bryan (2014), the interaction of AGN
outflows with the ICM leads to the formation of cold clouds along the jet
axis. However, the resulting structures are very clumpy and lack the coherent
and elongated morphologies observed in surveys of CC clusters. Other
simulations with coarser spatial resolutions (typically a few hundreds of pc)
have been performed to include jet re-orientation (Beckmann et al. 2019),
radiative transfer (Qiu et al. 2019), magnetic fields (Wang et al. 2021;
Ehlert et al. 2023), cosmic rays (Wang et al. 2020; Beckmann et al. 2022) and
anisotropic thermal conduction (Beckmann et al. 2022). Magnetic fields are
thought to play a significant role in the formation of elongated threads of
cold gas, suppressing the formation of cold clumps as well as massive cold
disks that are common in hydrodynamical setup (Li & Bryan 2014; Qiu et al.
2019). However, MHD simulations of idealised CC clusters with spatial
resolutions resolving the inner structure of the filaments (i.e., down to
scales of a few tens of pc) are still lacking, and consequences for AGN
feedback remain largely unconstrained. In this paper we present a set of high-
resolution, three-dimensional hydrodynamical and MHD simulations of an
idealised Perseus-like galaxy cluster with a turbulent, radiatively cooling
ICM and self-regulated AGN feedback. Our GPU-accelerated simulations are
performed on a static, nested mesh that is refined down to a maximum spatial
resolution of 24.4 pc.

In Sect. 2 we present the numerical setup of our simulations and describe the
parameters used for each of our runs. In Sect. 3 we discuss the overall
evolution of the ICM properties over the entire duration of our simulations.
In Sect. 4 we study the morphology and main physical properties of the cold
gas filaments found in our high–resolution MHD run. In Sect. 5 we take a
closer look at the formation of several filaments and discuss the relative
contribution of cold gas uplifting and condensation. In Sect. 6, we assess the
impact of cold gas on the morphology of AGN outflows, before discussing our
results and assessing the importance of missing physical process in our
simulation in Sect. 7 and concluding in Sect. 8.

##  2 Simulations and method

###  2.1 Magnetohydrodynamics with radiative cooling

We use the performance-portable AthenaPK code based on the block-structured
refinement framework Parthenon (Grete et al. 2023). The fluxes are computed
using Harten-Lax-van Leer Contact (HLLC) (Toro et al. 1994) and Harten-Lax-van
Leer discontinuities (HLLD) (Miyoshi & Kusano 2005) solvers for our
hydrodynamical and MHD runs, respectively. Overall, a second-order accurate
method is used based on piecewise linear reconstruction, a predictor-corrector
Van Leer-type integrator, and, in the MHD case, a generalized Lagrange
multiplier based divergence cleaning approach (Dedner et al. 2002). The basic
system of equations is given by:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot[\rho\mathbf{u}]=0,$
(1) $\displaystyle\frac{\partial\rho\mathbf{u}}{\partial
t}+\nabla\cdot\left[\rho\mathbf{u}\otimes\mathbf{u}+P_{\text{tot}}\mathbb{I}-\mathbf{B}\otimes\mathbf{B}\right]=0,$
$\displaystyle\frac{\partial e_{\text{tot}}}{\partial
t}+\nabla\cdot\left[\mathbf{u}(e_{\text{tot}}+P_{\text{tot}})-\mathbf{u}\cdot\mathbf{B})\mathbf{B}\right]+\mathcal{L}(\rho,T)=0,$
$\displaystyle\frac{\partial\mathbf{B}}{\partial
t}-\nabla\times[\mathbf{u}\times\mathbf{B}]=0.$

Here, $\rho$ is the gas density, $\mathbf{u}$ its velocity field, $\mathbf{B}$
is the magnetic field vector, $P_{\text{tot}}$ is the total pressure including
the thermal pressure $P=(1-\gamma)e$ and the magnetic field pressure
$P_{\text{mag}}=B^{2}/2$. The gas is assumed to be non-relativistic and thus
we take $\gamma=5/3$. $e$ is the gas internal energy density and
$e_{\text{tot}}$ designates the total energy density $e_{\text{tot}}=e+0.5\rho
u^{2}+B^{2}/2$.

$\mathcal{L}(\rho,T)$ is a function modelling the optically thin cooling
processes taking place in the gas component. We include radiative cooling of
gas down to $10^{4.2}$ K using the exact integration scheme introduced in
Townsend (2009) and the cooling table from Schure et al. (2009), assuming half
solar metallicity (Fabian et al. 2011). A temperature floor is set to
$T_{\rm{floor}}=10^{4}$ K. Finally, given the non-relativistic setup of our
simulations we also limit the temperature to $T<5\times 10^{9}$ K, and gas and
Alfvén velocities to $0.05$ c in the central region of the simulation ($r<20$
kpc).

###  2.2 Grid structure

We adopt a setup akin to that employed by Wang et al. (2021) and Grete et al.
(2024, in prep.). A root grid of $512^{3}$ cells covering a volume of (800
kpc)3 is statically refined in its central region. The width of cells at
refinement level $\ell$ is given by:

$\Delta x_{\ell}=\frac{L_{\rm{box}}}{512\times 2^{\ell}}.$ (2)

In our highest resolution runs we use $\ell=6$ refinement levels. The three
finest refinement levels, with cell sizes of 24.4 pc, 48.8 pc, and 97.6 pc,
cover the inner 6.4 kpc, 12.6 kpc, and 25 kpc of the box, respectively. Such
grid geometry allows us to simulate $\sim 80\%$ of the total cold gas mass
formed with cells of 48.8 pc at most in our highest resolution runs (see
Appendix A).

###  2.3 Initial conditions

In all our simulations we impose a constant gravitational potential. The dark
matter halo of the cluster is modeled using a Navarro-Frenk-White (NFW)
profile (Navarro et al. 1997):

$\displaystyle g_{\text{NFW}}(r)$
$\displaystyle=\frac{\mathcal{G}}{r^{2}}\frac{M_{\text{NFW}}\
\left[\ln\left(1+\frac{r}{R_{\text{NFW}}}\right)-\frac{r}{r+R_{\text{NFW}}}\right]}{\ln(1+c_{\text{NFW}})-\frac{c_{\text{NFW}}}{1+c_{\text{NFW}}}},$
(3)

where $\mathcal{G}$ is the gravitational constant, $M_{\text{NFW}}$ is the NFW
profile virial mass, $c_{\text{NFW}}$ its concentration parameter and
$R_{\text{NFW}}$ is a scale radius function of the virial mass and the
concentration parameter, defined by the following relations:

$\displaystyle
R_{\rm{NFW}}=\left(\frac{M_{\rm{NFW}}}{4\pi\rho_{\rm{NFW}}\left[\ln{\left(1+c_{\rm{NFW}}\right)}-c_{\rm{NFW}}/\left(1+c_{\rm{NFW}}\right)\right]}\right)^{1/3},$
(4)
$\displaystyle\rho_{\rm{NFW}}=\frac{200\,H_{0}^{2}}{8\pi\mathcal{G}}\frac{c_{\rm{NFW}}^{3}}{\ln{\left(1+c_{\rm{NFW}}\right)}-c_{\rm{NFW}}/\left(1+c_{\rm{NFW}}\right)}.$

Secondly, a Hernquist potential is used to model NGC1275, the brightest
cluster galaxy (BCG), with $M_{\text{BCG}}$ and $R_{\text{BCG}}$ its mass and
radius (Hernquist 1990):

$\displaystyle g_{\text{BCG}}(r)$
$\displaystyle=\frac{\mathcal{G}M_{\text{BCG}}}{R_{\text{BCG}}^{2}}\left(1+\frac{r}{R_{\text{BCG}}}\right)^{-2}.$
(5)

Finally, we add a central point mass modeling the effect of the supermassive
black hole (SMBH) of mass $M_{\rm{BH}}$.

Density, pressure, and temperature profiles are then calculated from this
gravitational field so that the ICM obey hydrostatic equilibrium. We assume an
entropy profile, $K(r)$, defined by:

$K\equiv\frac{k_{B}T(r)}{n_{e}^{2/3}(r)},$ (6)

where $k_{B}$ is the Boltzmann constant, $r$ is the radial distance to the
center of the cluster, $T(r)$ is the temperature at a given distance, $r$, and
$n_{e}(r)$ is the electron density. The profile follows the empirical power-
law introduced in the ACCEPT database (Cavagnolo et al. 2009):

$K(r)=K_{0}+K_{100}\ \left(\frac{r}{100\text{ kpc}}\right)^{\alpha_{K}}.$ (7)

In order to close the system of equations defining the hydrostatic
equilibrium, we fix the electron density of the cluster to a defined value
$n_{e}(r_{\rm{ref}})$ at a defined radius $r_{\rm{ref}}=10$ kpc.

Parameter | Value | Unit | Reference
---|---|---|---
$H_{0}$ | $75.1$ | $\rm{km}\,\rm{Mpc}^{-1}\,\rm{s}^{-1}$ | a
$M_{\rm{NFW}}$ | $6.6\times 10^{14}$ | $\rm{M}_{\odot}$ | b
$c_{\rm{NFW}}$ | 5.0 | - |
$R_{\rm{BCG}}$ | $10$ | $\rm{kpc}$ | c
$M_{\rm{BCG}}$ | $2.4\times 10^{11}$ | $\rm{M}_{\odot}$ | d
$M_{\rm{BH}}$ | $4\times 10^{8}$ | $\rm{M}_{\odot}$ | e
$K_{0}$ | $10.0$ | $\text{ keV}\ \text{cm}^{2}$ | f
$K_{100}$ | $150.0$ | $\text{ keV}\ \text{cm}^{2}$ | f
$\alpha_{K}$ | $1.1$ | - | f
$n_{e}(r_{\rm{ref}})$ | $0.05$ | $\rm{cm}^{-3}$ | g
Table 1: Values adopted for the various initial condition parameters in our
simulation runs. References: a) Anand et al. (2022), b) Simionescu et al.
(2011), c) de Vaucouleurs et al. (1991), d) Mathews et al. (2006), e) Park &
Trippe (2017), f) Cavagnolo et al. (2009), e Churazov et al. (2003).

In order to break the symmetry of the system, we add velocity perturbations in
each run. 40 wave vectors with lengthscales between 12.5 kpc and 50 kpc are
sampled using a random seed and their relative contribution is computed from
an inverse parabolic spectrum. We take a peak wavelength of 25 kpc and a
velocity dispersion of 75 km s-1. We use a similar spectrum to initialize a
perturbed magnetic field, whose average local amplitude follows the density
profile:

$\langle B(r)\rangle\propto\rho(r)^{2/3}.$ (8)

The strength of the initial magnetic fields is chosen such that its average
intensity in the inner few kpc of the box is $\sim 10\,\mu\rm{G}$ (Taylor et
al. 2006). The same random seed is used in each of our runs so that the
perturbed fields are identical from one simulation run to another.

The values adopted for each of these parameters are summarized in Table 1.

###  2.4 Active galactic nuclei feedback

We model the effect of the central SMBH on its environment by including AGN
feedback through three different channels, namely thermal, kinetic and
magnetic feedback.

####  2.4.1 Cold gas accretion

Since the ICM is thermally unstable, turbulence is expected to trigger the
formation of cold clouds which decouple from the surrounding hot gas and fall
chaotically towards the center of the cluster. Accordingly, we follow the
chaotic cold accretion model developed in Gaspari et al. (2013) and later
implemented in cluster scale simulations (Li & Bryan 2014; Meece et al. 2017).
We use an AGN triggering model based on the accretion of cold material. We
define the accretion region as a sphere of radius $R_{\rm{acc}}=500$ pc, and
select the gas with a temperature below a threshold value
$T_{\rm{cold,acc}}=5\times 10^{4}\,\rm{K}$ as the accretion material. The
accretion rate at each time step is then derived by summing mass over all
cells containing gas beneath $T_{\rm{cold,acc}}$ and within the accretion
region:

$\dot{M}_{\rm{acc}}=\sum_{r_{i}<R_{\rm{acc}}}\frac{\rho(\mathbf{\hat{r}}_{i})\
dV(\mathbf{\hat{r}}_{i})}{t_{\rm{acc}}}\
\delta(T(\mathbf{\hat{r}}_{i})<T_{\rm{cold,acc}}),$ (9)

where $\rho(\mathbf{\hat{r}}_{i})$ is the gas density in cell $i$, $d\rm{V}$
is that cell’s volume. Cold gas in these cells is depleted over the accretion
timescale $t_{\rm{acc}}=5$ Myr. This accretion rate is then used to compute
the AGN feedback power, given by:

$\dot{E}_{\text{AGN}}=\eta\ {\dot{M}}_{\rm{acc}}c^{2},$ (10)

where $\eta=10^{-2}$ is the cold mass triggering efficiency (Gaspari et al.
2011). The energy injected in each feedback channel is then a fraction of
$\dot{E}_{\text{AGN}}$, given by:

$\dot{E}_{X,\text{AGN}}=f_{X}\ \dot{E}_{\text{AGN}},$ (11)

where $X$ designates the feedback mode (i.e., $f_{T}$, $f_{K}$ and $f_{M}$ for
thermal, kinetic and magnetic feedback respectively, see Table 2).

####  2.4.2 Thermal feedback

We model the radiative feedback from the AGN by depositing a portion of the
AGN feedback power, $f_{T}\cdot\dot{E}_{\rm{AGN}}$ through thermal energy
within a sphere of radius $R_{T}=500\,\rm{pc}$. The corresponding power is
related to the AGN power (Eq. 10) through:

$\dot{e}_{T}(r\leq R_{T})=\frac{3}{4}\ \frac{f_{T}\ \dot{E}_{\rm{AGN}}}{\pi
R_{T}^{3}}.$ (12)

Mass is also injected within the same region of radius $R_{T}$:

$\dot{\rho}_{T}(r\leq R_{T})=\frac{3}{4}\ \frac{f_{T}\ \dot{M}_{\rm{acc}}}{\pi
R_{T}^{3}}.$ (13)

####  2.4.3 Kinetic feedback

We model the jets by redistributing material in two cylinders of radius
$r_{\rm{jet}}=4\,\Delta x_{\rm{min}}$, thickness $h_{\rm{jet}}=2\,\Delta
x_{\rm{min}}$ and with an offset $L_{\rm{jet}}=4\,\Delta x_{\rm{min}}$ with
respect to the center of the box. We position these cylinders in a static
orientation along a fixed axis. The velocity of the jet is determined by

$v_{\rm{jet}}=\sqrt{2\eta}\ c,$ (14)

with $c=3\times 10^{8}\,\rm{km}\,\rm{s}^{-1}$ the speed of light. The density
of the deposited material is:

$\rho_{\rm{jet}}=\frac{(1-\eta)f_{K}\dot{M}_{\rm{acc}}}{2\pi
r_{\rm{jet}}^{2}h_{\rm{jet}}}.$ (15)

In order to track the jet material as passive scalar is evolved whose
concentration is set to unity within the jet launching regions.

Figure 1: Evolution of the AGN power (top panel) as a function of time for our
highest resolution hydrodynamical (dark blue) and MHD (light blue) runs. The
bottom panel presents the evolution of the total cold gas mass formed (solid
line), as well as its associated specific angular momentum (dashed line). The
average AGN power is broadly constant from one setup to an other, but shows
weaker fluctuations in MHD setup. While similar quantities of cold gas are
formed in each run, its overall specific angular momentum is found to be
larger in the hydrodynamical case due to the formation of massive disk,
consistently with previous work (Li & Bryan 2014; Qiu et al. 2019; Ehlert et
al. 2023).

####  2.4.4 Magnetic feedback

The importance of magnetic fields in the acceleration and collimation of AGN
jets has been emphasized by numerous theoretical and observational studies.
General relativistic MHD simulations have shown that these are dynamically
important in both the launching and propagation of the jets (Liska et al.
2020; Lalakos et al. 2024), while Farady rotation measurements in AGN jets
leads to strength ranging from $10^{2}\,\rm{G}$ in lauching regions down to
typical ICM strengths of a few $\mu$G (Kino et al. 2015; Park et al. 2019;
Gómez et al. 2022).

We include magnetic feedback in our MHD runs by depositing purely toroidal
magnetic field loops in all cells at the same height as the kinetic feedback
cylinders (i.e. within $L_{\rm{jet}}\leq|z|\leq L_{\rm{jet}}+h_{\rm{jet}}$),
following:

$\mathcal{B}_{\theta}(r)=2\mathcal{B}_{0}\ \left(\frac{r}{L_{B}}\right)\
\exp\left[-\left(\frac{r}{L_{B}}\right)^{2}\right],$ (16)

in cylindrical coordinates. $L_{B}=2.5\,\Delta x_{\rm{min}}$ is the
characteristic length scale of the toroidal field set so that the injected
fields are mostly confined to the launching region and $\mathcal{B}_{0}$
dynamically adjusted so that the resulting injected magnetic energy matches
the magnetic fraction of the target AGN feedback.

###  2.5 Stellar feedback

SNIa-type feedback is included to account for star-formation processes
following Prasad et al. (2020). The energy density and mass density injected
into a cell at a distance, $r$, from the cluster center is:

$\displaystyle\dot{e}_{\rm{SNIa}}=\Gamma_{\rm{SNIa}}\ E_{\rm{SNIa}}\
\rho_{\rm{BCG}}(r),$ (17) $\displaystyle\dot{\rho}_{\rm{SNIa}}=\alpha\
\rho_{\rm{BCG}}(r),$

where $\Gamma_{\rm{SNIa}}=3\times 10^{-14}\,\rm{yr}^{-1}\ \rm{M}_{\odot}^{-1}$
is the SNIa rate per unit mass, $E_{\rm{SNIa}}=10^{51}\rm{erg}$ is the amount
of energy released per supernova event, and $\alpha=10^{-19}\ \rm{s}^{-1}$ is
the mass injection rate.

In the absence of separate star particles, we also instantaneously convert gas
that would end up in stars into thermal energy. More specifically, a fraction
of the gas above a number density threshold of $n>50\,\rm{cm}^{-3}$, a
temperature below $T<2\times 10^{4}$ K, and within a radius of
$R_{\mathrm{acc}}<r<25$ kpc is locally converted into heat with an efficiency
of $5\times 10^{-6}$.

###  2.6 Simulation parameters

Table 2: Overview of the parameters used in our runs. Columns indicate (from left to right): simulation name; width of cells in the maximally refined region of the box; feedback fractions for kinetic, thermal and magnetic energy; average initial magnetic field amplitude at the center of the box. Name |  | Max. res.
---
[pc]
$f_{K}$ | $f_{T}$ | $f_{M}$ |  | $\langle B(r=0)\rangle$
---
[$\mu$G]
hydro | 24.4 | 0.75 | 0.25 | - | -
mhd | 24.4 | 0.74 | 0.25 | 0.01 | 10

We run two high-resolution simulations whose parameters are summarized in
Table 2. In addition to our high-resolution run with minimum cell width of
24.4 pc, we perform three hydrodynamical and three MHD simulations with fewer
refinement levels in order to study the impact of numerical resolution on our
results. The dependence of our results on spatial resolution is discussed in
Appendix A.

##  3 Cluster evolution

Because our AGN triggering mode is only sensitive to cold gas, the SMBH is
inactive during the first $\sim$400 Myr of our simulation. During that time
the gas properties evolve due to radiative cooling and the initial velocity
turbulence. Overdensities build up in the center of the box and end up forming
cold clumps after roughly 400 Myr. In Fig. 1 we present the evolution of the
AGN feedback power (top panel), total cold gas mass (bottom panel, solid
line), and its associated specific angular momentum (dashed line) as a
function of time for our high-resolution MHD and HD simulations. The average
AGN power is roughly constant across runs and is of order $10^{44}\
\rm{erg}\,\rm{s}^{-1}$. Phases of low AGN power down to $\sim 10^{42}\
\rm{erg}\,\rm{s}^{-1}$ lasting for a few tens of Myr are visible in the
hydrodynamical run, with a typical period of 100 – 200 Myr. The overall
feedback power remains much steadier in the MHD run, with weaker fluctuations
that last for shorter periods. The total cold gas mass formed throughout our
simulations converge to typically $\sim 10^{9}\ \rm{M}_{\odot}$ after $\rm{1}$
Gyr, consistent with observational constraints (Salomé et al. 2006; Fabian et
al. 2008) and previous simulations (Beckmann et al. 2019; Li & Bryan 2014;
Ehlert et al. 2023). We find that the specific angular momentum of the cold
gas is less than in our MHD setup. A massive cold disk quickly forms in our
hydrodynamical run, while we find no such structure in our MHD runs. Such
behavior are in agreement with previous simulation work (Wang et al. 2020;
Ehlert et al. 2023). In Fig. 2 we show the cooling time and electron density
radial profile of the hot ICM at various times in our MHD and hydrodynamical
simulations, from the start of the first AGN phase at $t\sim 400$ Myr to the
end of our runs. Observational constraints for the cooling time profile is
also presented by the gray lines and points (Sanders et al. 2004; Dunn &
Fabian 2006; Zhuravleva et al. 2014; Fabian et al. 2017). We find no major
differences between our hydrodynamical and MHD setups regarding the overall
properties of the ICM. The cooling profiles remain broadly stable over time,
with values of $2-3\,\times\,10^{8}\,\rm{yr}$ in the inner 10 kpc,
consistently with observations. The electron density increases in both our
runs up to $\sim 10^{-1}\,\rm{cm}^{-3}$. At late times, the MHD runs show an
elevated cooling time.

Figure 2: Evolution of the radial profiles for cooling time and electron
density evaluated at various times in our highest resolution MHD (top) and
hydrodynamical (bottom) runs for the hot ICM gas (i.e., $T\geq 10^{5}$ K).
Observational constraints for the Perseus cluster are also indicated for
cooling time for comparison (Sanders et al. 2004; Dunn & Fabian 2006;
Zhuravleva et al. 2014; Fabian et al. 2017).

##  4 Properties of the cold gas

Figure 3: Composite image showing the density-weighted projection of the
magnetic field strength (blue) and the temperature (orange) for the inner 25
kpc of our high-resolution MHD run at $t=1025$ Myr. Mock H$\alpha$ and X–ray
images are also shown in the lower right of the image. The top H$\alpha$ image
is untouched, while the center right one has been smoothed with a Gaussian
filter to match the angular resolution of the Muse instrument assuming our
cluster is at the same distance as the Perseus cluster, and is comparable to
Fig. 2 in Olivares et al. (2019). The X–ray image is a superposition of three
filters, namely 0.3–1.2 keV (red), 1.2–2 keV (green) and 2–7 keV (blue) bands,
and is comparable to Fig. 1 from Fabian et al. (2006). The white dashed line
is a contour map from the temperature field, showing the position of the
outflowing bubbles.

###  4.1 Cluster visualization

To begin our exploration of the cold gas properties in our simulated cluster,
we present a series of visualizations from our high-resolution MHD simulation
in Fig. 3. The main image shows the magnetic field strength (blue), overlaid
with a projection of the temperature field (orange). Elongated structures of
gas with magnetic field strengths of order $\sim 100\,\mu$G are visibly
interacting with the outflowing AGN jets, deflecting them and inflating
several cavities. In the top and center right of the figure, we present mock
H$\alpha$ image of these cold structures. The second H$\alpha$ flux has been
smoothed with a gaussian filte to match the angular resolution of the MUSE
instrument assuming that our simulated cluster is at the same distance from
the observer as the Perseus cluster, and is comparable to Fig. 2 from Olivares
et al. (2019). Our simulations show clumpy internal sub-structure inside the
filaments that is unresolved by MUSE. The image in the lower right is a
composite X–ray mock image imitating Fig. 1 from Fabian et al. (2006), and
traces the outflowing bubbles carved in the ICM by the AGN jets. The red,
green and blue components represents the 0.3–1.2 keV, 1.2–2 keV and 2–7 keV
bands, respectively.

###  4.2 Morphology

Previous research have shown that magnetic fields affect the overall
morphology of the cold gas. While high-resolution hydrodynamical simulations
can produce cold gas in quantities consistent with observations, the
morphology of the structures is usually very clumpy (Li & Bryan 2014) and
lacks a coherent, filamentary morphology. Magnetic fields have been found in
both cluster and cloud–crushing setups to inhibit the formation of clumps and
massive disks, and to enhance the formation of elongated structures (Wang et
al. 2020; Ehlert et al. 2023; Das & Gronke 2024).

In Fig. 4, we present a gas density projection of the inner 20 kpc in our two
high-resolution simulations. The left panel shows the hydrodynamical run and
the right panel the MHD run. In agreement with previous work (Ehlert et al.
2023, e.g.) we find that the magnetic fields have an important impact on the
morphology of the cold gas. In the hydrodynamical case the cold gas condenses
into thousands of cloudlets typically consisting of just 1–10 cells. In our
MHD setup, the gas gathers along extended filamentary structures
preferentially infalling towards the center of the cluster. Each of these
structures are dynamic, reach maximum elongations up to $\sim 15$ kpc, and
contain up to $\sim 10^{8.5}\,\rm{M}_{\odot}$ of cold gas.

Figure 4: Gas density projection of the inner 20 kpc of our simulated cluster
at a $t=920$ Myr for our hydrodynamical (left panel) and MHD runs (right
panel). The gas mostly gathers into tiny cloudlets of a few tens of pc, while
magnetic fields promote the formation of elongated filaments.

The morphologies observed in our hydrodynamical run are hard to reconcile with
observations. Optical and radio surveys have shown that cold gas in CC
clusters is concentrated in filaments and that massive cold disks are rare
(Olivares et al. 2019; Tremblay et al. 2018; Russell et al. 2017). In the
following, and unless stated otherwise, we thus focus on the results from our
high–resolution MHD run.

Figure 5: Projection of the gas density for five different filaments in our
fiducial MHD run. Each projection is performed at the time when each filament
reaches its maximum elongation. The projected position of the SMBH and the
projected orientation of the jets is indicated by the round marker and the
dotted line, respectively. Figure 6: Main statistical properties of the clump
sample obtained from the five filaments presented in Fig. 5. The left panel
shows the cumulative mass fraction of the clumps as a function of their mass.
The PDF of the clumps for the five filaments as a function of their mass and
radial velocity are presented in the central and right panel, respectively.
The number of clumps for each filament are indicated in the legend of the left
panel. The first number indicates the number of clumps accounted for our
statistical analysis while the number in parentheses indicates the total
clumps population with no constraint on minimum number of cells per clump. The
associated clump mass distribution with no constraint on the minimum number of
cell is also indicated by the bright blue dashed lines. The orange dashed line
indicates a $d\rm{N}/d\rm{M}$ power law distribution of slope $\alpha=-0.6$
for comparison. The vertical dashed lines on the right panel indicates the
mass weighted average of each filament’s distribution.

Among the cold structures observed over the course of our fiducial MHD
simulation, we identify five filaments that lend themselves to a morphological
study because their relative isolation helps identify their boundaries. In
Fig. 5, we present the gas density projections of these five distinct
filaments at the time they reach their maximum elongation. These are referred
to by their index indicated on the top left of each panel. The typical
filaments length ranges between 5 and 15 kpc, and their morphology can be
complex. In particular, all filaments exhibit numerous inner structures of
cold clumps and sub-filaments of various thicknesses spanning between $\sim
800$ pc and the resolution limit of the simulation. To study quantitatively
the distribution of these substructures in and around these five filaments, we
use a clump-finding algorithm integrated in the visualization and data
analysis package Yt (Turk et al. 2011, 2024). The algorithm operates by first
defining a minimum density threshold. It then employs a contouring method to
recursively construct a hierarchical structure of clumps and sub-clumps. In
our setup, we use a minimum density threshold
$n_{\rm{clump,min}}=1\,\rm{cm}^{-3}$ and a step factor $s=4$ corresponding to
the multiplication factor applied to the minimum density threshold at each
iteration of the algorithm. For each system, we limit the search for dense
structures to all gas cells contained in a cylinder of radius $r_{\rm{cyl}}=2$
kpc manually aligned with each filament’s main axis. We also exclude all
clumps with fewer than 10 cells and with mean temperature above
$T_{\rm{clump,max}}=10^{5}$ K to exclude possible artifact from the hot ICM.
Several cell volume-weighted quantities are calculated for each filament by
summing over all gas cells belonging to a given clump. In particular, we
derive the magnetic field strength, the ratio of the thermal pressure to
magnetic pressure $\beta$ (later on referred as plasma $\beta$) and the radial
velocity of each clump belonging to each of the five filaments.

In Fig. 6 we present the cumulative mass fraction of the clumps belonging to
each filament, as well as the probability distribution function (PDF) of the
clumps as a function of their gas mass (central panel) and radial velocity
(right panel). $N$ in the left panel indicates the number of clumps identified
in each filament included for our analysis, i.e., each containing 10 or more
cells. The corresponding PDF are represented by the blue solid lines. The
number in parenthesis corresponds to the total number of clumps with no
constraint on the minimum number of cell for each clump. The associated PDF
are also presented by the dashed blue lines for comparison.

As visible in the left panel, a large fraction of the mass of each filament is
carried by massive clumps. In particular, the cumulative mass of all clumps
with individual masses $M_{\rm{clump}}\leq\leavevmode\nobreak\ 10^{7.5}\
\rm{M}_{\odot}$ contributes to 25 % at most of the overall mass of each
filament, while 50 % of each filament’s mass is carried solely by the few
clumps with masses $M_{\rm{clump}}\geq 10^{8}\ \rm{M}_{\odot}$. The number of
such massive clump is low and around 1–5 for each filament system.

The overall shape of the mass distributions (central panel) is remarkably
constant across filaments, and seems to follow a power law between
$10^{4.7}\,\rm{M}_{\odot}$ and $10^{8}-10^{9}\,\rm{M}_{\odot}$. The
corresponding slope of the clump number distribution
$\rm{d}\log\rm{N}/\rm{d}\log\rm{M}$ is typically of $\alpha=-0.6$ (see orange
dashed line for comparison), consistently with cloud–crushing simulations
(Jennings et al. 2023).

The left and central panel support a picture in which the filaments are clumpy
structures, largely dominated in number by light clumps with masses as low as
$10^{3}-10^{4}\,\rm{M}_{\odot}$, and largely dominated in mass by the
contribution of only a few massive clumps with individual masses above
$10^{8}\,\rm{M}_{\odot}$. We emphasize that most of the lightest clumps remain
spatially unresolved, and that increasing spatial resolutions might impact the
shape of the mass distribution, especially at its lower mass end. Indeed, the
position of the peak of the distribution is shifted to lower masses when
including the excluded clumps with fewer then 10 cells, suggesting that the
lower mass end of the clump statistics in not yet converged. This is expected
and consistent with results from 2D and 3D cloud-crushing simulations which
concluded that the fragmentation of cooling gas produces cloudlets with
typical radii reaching down to sub-parsec scales, well below our resolution
limit of 25 pc (McCourt et al. 2018; Jennings et al. 2023).

In the right panel of Fig. 6, we show the distribution of clumps as a function
of radial velocity. The vertical dashed lines indicates the mass weighted
average value for each filament’s distribution. As visible, the bulk motion of
the filaments is preferably infalling with average velocites between –100 and
–200 km s-1. All five filaments contain clumps characterized by a positive
radial velocity, up to +400 km s-1. These clumps are likely uplifted due to
their interaction with the AGN jets, and might play a role in the mass build-
up of the filaments (see Sect. 5).

Figure 7: Distribution of the clumps as a function of their volume-weighted
magnetic field strength (left panel) and their thermal-to-magnetic pressure
ratio $\beta$ (right panel), for the five filaments presented in Fig. 5.

###  4.3 Magnetic support

In Fig. 7 we present the distribution of the clumps as a function of their
volume-weighted magnetic field strength (left panel) and plasma $\beta$. The
median values of the magnetic field strength in the clumps sample of our five
filaments are typically of $\rm{Med}(B_{\rm{clumps}})\sim 80-100\ \mu$G. Some
of the clumps have a particularly high inner magnetic field strength, ranging
up to $\sim 200\ \mu$G. As shown in the right panel of Fig. 7, the magnetic
pressure is the main source of pressure supporting the filaments against the
thermal pressure of the ICM for more than a half the clumps.

We find a correlation between the magnetic field lines orientation within the
filaments and the filament’s main axis. As an example, we show the magnetic
field topology of the filament system #3 in Fig. 8 (see also associated
movies). The left panel shows the gas column density, highlighting the overall
structure of the clump system. The right panel shows a slice of the plasma
$\beta$, in the midplane of the system, overlaid with the magnetic field
lines. The dashed grey line indicates the contour of equipartition between
magnetic and thermal pressure. The magnetic field lines within the regions
dominated by magnetic pressure are aligned with the main axis of the filament.
This behavior is also found in the other filaments we have inspected (see also
Fig. 11). The magnetic field lines outside of the filament are typically
uncorrelated with the filament’s main axis. In certain regions, the lines are
found to flow towards the filament with an orthogonal orientation, and are
strongly bent by up to $\sim$130∘ to reach alignment with the filament.

Figure 8: Projection of the gas density (left) and slice of the plasma
$\beta$, overlaid with the magnetic field lines (right) for the filament #3.
The projection depth of the density plot is 4 kpc. The grey dashed line around
the filament core represents the contour of equi-partition between magnetic
and thermal pressure. (Movies available)

Figure 9: Total mass of upward-moving cold gas with radial velocities above 0,
150, 300, 400 and 500 km s-1 (from dark to bright blue). The orange solid line
indicates the evolution of the AGN feedback power.

##  5 Formation of the cold filaments

###  5.1 Gas uplifting and condensation

Figure 10: Evolution of the average radial velocity of three filament systems
exhibiting phases of uplift, as a function of their total cold gas mass
budget. The grey shaded area indicates positive radial velocities. At early
times, the bulk of the filament is uplifted due to its interaction with the
AGN outflows. The radial velocity of the system typically reverses $\sim$
20–40 Myr after the beginning of the uplift phase. The filament then continue
building up mass while falling back towards the center of the clusters for
$\sim$ 20–100 Myr, before being shredded. The time is rescaled such that $t=0$
coincides with the formation of each filament system. Figure 11: Sequence of
images spanning the whole lifetime of the filament system #2 (see Sec. 4.2).
The panels show, from left to right: gas column density, density-weighted
radial velocity, density of jet scalar and density-weighted thermal-to-
magnetic pressure ratio $\beta$, overlaid with the magnetic field lines. In
each case, the plotting region is aligned with the filament’s main axis. The
SMBH is located right of the plot and the jet axis has an angle of $\sim
10^{\circ}$ with the main filament axis. Time has been rescaled such that
$t=0$ coincides with the snapshot corresponding do the start of the clumps
uplifting phase. The times of the snapshots have been selected such that the
main steps of the filament evolution is most apparent, and are thus
voluntarily unevenly spaced.

We find evidence of ubiquitous cold gas uplifting during the entire course of
our fiducial MHD simulation. In Fig. 9, we show the evolution of
$M_{\rm{cold,up}}$, the total mass of upward-moving cold gas with radial
velocities a certain threshold, defined as:

$M_{\rm{cold,up}}(v_{\rm{rad}})=\sum_{i}m_{i}\,\delta(v_{\rm{rad,i}}\geq
v_{\rm{rad}})\,\delta(T_{i}\leq 10^{5}\,\rm{K}),$ (18)

with the sum running over every cell of the box, designated by their index
$i$. The various curves indicates the evolution of $M_{\rm{cold,up}}$ for the
following radial velocity thresholds: 0, 150, 300, 400 and 500 km s-1 (from
dark to bright blue lines). For comparison, the orange curve indicates the
evolution of the AGN feedback power $\dot{E}_{\rm{AGN}}$.

After about 400 Myrs, the total uplifted gas mass converges to $10^{8.4}\
\rm{M}_{\odot}$ which is about 15% of the total cold gass mass (e.g., see Fig.
1). There is no obvious correlation between this quantity and the AGN power as
the latter flickers on timescales down to sub–Myr while the typical variations
of the cold gas mass are smoother and typically of a few Myr at least. As the
minimum radial velocity threshold is increased (i.e., for brighter blue
curves), correlations between peaks of AGN power and peaks of uplifted cold
mass starts to be noticeable. In particular, some peaks and drops of AGN power
are correlated with peaks and drops of $M_{\rm{cold,up}}$, especially for
radial velocity threshold $\geq 300$ km s-1.

In Fig. 10, we show the evolution of the total cold mass as a function of
mass-weighted average radial velocity for the five filament systems presented
in Fig. 5. Four of these filaments (namely the filaments #1, #2, #4 and #5)
exhibit an initial phase of uplifting, with average radial velocities ranging
from 20 – 100 km s-1. Since not all clumps are uplifted by jet-cloud
interactions and because the region within which the physical quantities are
computed also contains infalling clumps, the average radial velocity is a
lower estimate of the overall upward motion of the cold material as individual
clumps can uplifted with radial velocities much beyond these values, up to 400
km s-1. For the filament systems $\\#2$ and $\\#4$, the total amount of gas
present in the system during the uplifting phase is 30 % and 50 % of the
maximum mass reached by the filament in its lifetime, respectively. The
uplifting phase lasts for a time of 20 – 40 Myr in these four filaments after
which the system starts falling back towards the center of the cluster. The
filament #3 exhibits a different behavior as no clear traces of uplifting is
found when visually inspecting the filament over its formation time. It is
likely that this system is solely condensing out of the ICM.

###  5.2 Evolution of the filaments structure

Figure 12: Evolution of the properties of the filament system (#2) (see also
Fig. 11) at five equally spaced times. Top left panel: evolution of the total
cold mass budget of the system (solid line) and position of the center-of-mass
of the system (dashed line) as a function of time. Top right panel:
distribution of the radial velocity for various times, along with their
associated median values (dashed line). The grey shaded area represents the
positive radial velocity region. Lower left panel: mass distribution of the
clumps at various times. Lower right panel: distribution of the plasma $\beta$
at various times and associated median values (dashed line). The number $N$ in
the top right panel indicates the number of clumps found at each time, with
and without (in parenthesis) the constraint on the minimum number 10 of cells
for a clump to be included in the sample.

A sequence of images showing the formation and evolution of the filament
system #2 is shown in Fig. 11. The panel shows, from left to right: column gas
density, density-weighted radial velocity, and a midplane slice of the plasma
$\beta$, overlaid by magnetic field lines. At the initial state (first row),
the system consists of a distribution of small, magnetically supported
cloudlets and short filaments of length $\sim$ 1–2 kpc, and thicknesses close
to the resolution limit. The radial velocity of this distribution of clouds
shows strong internal variations depending on the location, with some clouds
falling towards the center of the cluster, and other ones exhibiting positive
radial velocity above 200 km s-1, likely resulting from their interaction with
the jet. After 40 Myr (second row), a more extended distribution of cloudlets
starts gathering along the main axis of the plot, on a length scale of $\sim$
10 kpc. The radial velocity profile varies smoothly along the distribution
from positive radial velocities at the tip of the structure, to steady and
infalling motions in regions closer to the center of the SMBH. Interestingly,
the clouds located at the tip of the expanding structure have velocities above
300 km s-1 at $t=40$ Myr. Considering the gravitational field imposed in our
simulation (see Sect. 2.3), a simple integration of the equations of motion
lead to the conclusion that the clumps should come to rest at $t\sim 40$ Myr
if their motion was ballistic. Hence, the presence of some clumps with radial
velocities above 300 km s-1 after several tens of Myr suggests that the
interaction between the AGN outflows and the cold clouds lasts for several
tens of Myr, sustaining momentum transfer from the jets to the cold clouds.

At $t\sim$ 100 Myr a clear filamentary structure starts to appear which
extends up to distances of $\sim 15$ kpc from the center of the cluster. The
bulk of this thread is magnetically supported, as $\beta$ drops by nearly two
orders of magnitude between the ICM and the core of the filament where its
reaches values of $\beta\sim 10^{-1}-10^{-2}$. At later times, the filament
falls back toward the center of the cluster with radial velocities down to -
(200 – 300) km s-1. At $t=170$ Myr, the filament structure of the system
vanishes. The resulting distribution of clouds looks similar to the initial
state, with dozens of clumps and small-scale filaments.

While the magnetic field topology at $t=0$ is uncorrelated with the subsequent
orientation of the filament main axis, the field lines in the inner region of
the structure are visibly aligned with the thread as of $t=100$ Myr, likely
due to flux freezing and magnetic field amplification during the growth of the
cloud. Interestingly, the magnetic field lines near the base of the filament
are strongly bent at its surface, as they are orthogonal to its main axis
outside of the structure. This specific configuration lasts for nearly 100 Myr
as the filament falls back towards the center of the cluster.

In Fig. 12, we show the evolution of several physical quantities
characterizing this same filament system #2 presented in Fig. 11 (see also
left panel of Fig. 10). The top left panel indicates the evolution of the cold
gas mass contained in the filament (solid line), as well as the distance of
its center-of-mass to the center of the cluster (dashed line). The other
panels shows the PDF of the clumps belonging to this filament system for
various physical quantities, namely: the radial velocity of the clumps (top
right panel), their mass distribution (lower left) and their plasma beta
(lower) right. Time is here referred with respect to the begin of the
formation of the filament, i.e. when the first clumps starts being uplifted
through interaction with the jet (see Fig. 11). The system of clumps reaches
its maximum altitude at $t=100$ Myr, while the cold gas mass peaks 30 Myr
later, suggesting that a large fraction of the gas (roughly 60 %) forms
through condensation when the filament falls back towards the center of the
cluster. Interestingly, both the mass and plasma beta PDF exhibit a constant
shape over the course of the simulation, even though the overall morphology of
the system varies widely over the life cycle of the filament. An excess of
clumps with mass $\geq 10^{7}\,\rm{M}_{\odot}$ is visible at $t\sim 125$ Myr,
while the system falls back towards the center of the cluster. We find no
clear evolution in the magnetic pressure inside the clumps. The median value
of $\beta$ remains remarkably constant, around $\sim 10^{-1}$, and the
proportion of the number of clumps with a dominant magnetic pressure (i.e.
with $\beta<1$) is $>$90 % at all times.

As shown in the top right panel, the system does not move as a rigid body
since the radial velocity PDF exhibits a large scatter, ranging from $-400$ to
$+300\,\rm{km}\,\rm{s}^{-1}$. Even when the overall bulk motion motion of the
filaments falls back towards the center of the cluster (i.e. as of $t\sim 100$
Myr), up to 30 % of all the clumps still have a positive radial velocity,
likely due to constant interactions between clumps located close to the
feedback injection region and the AGN jets.

As discussed above, the formation of filaments in our MHD run results from the
combination of both gas uplifting and condensation. Since our current setup
does not include tracer particles, we cannot unambiguously disentangle the
contribution of these two phenomenon and we postpone such analysis to
subsequent study.

Figure 13: Volume-rendering of the inner (20 kpc)3 of our fiducial MHD run at
three different times. The yellow and white components represent the
isodensity surfaces for $3\times 10^{-3}\ \rm{cm}^{-3}$ and $2\times 10^{-2}\
\rm{cm}^{-3}$ and traces the outflowing AGN jet material. The blue component
traces the cold gas of density $20\ \rm{cm}^{-3}$. The figure shows that, even
though the AGN jets are accelerated along the vertical axis of this image,
interaction with the dense infalling filaments leads to the formation of
cavities that move in directions strongly inclined relative to the main jet
axis. The round black and white marker at the center of each panel represents
the position of the SMBH particle. (Movies available)

##  6 Impact of cold gas on the AGN outflows

Cavities in galaxy clusters have been widely studied using X-ray observations.
Typically, pairs of cavities with radii of $\sim 10$ kpc are found up to a few
tens of kpc away from the central galaxy. These are sometimes filled with
radio emission associated with AGN radio jets which inflate the cavities
(Boehringer et al. 1993; Bîrzan et al. 2004). Understanding the processes
shaping these cavities is critical to assessing how AGN feedback can balance
the cooling of the ICM. In previous work it has been discussed that high-
resolution simulations (typically with maximally refined cells smaller than
$\sim$ 100 pc) tend to form morphologies closer to FR I galaxies, with
strongly collimated jets that propagate out to distances of up to hundreds of
kpc. A direct consequence is a less isotropic heating of the ICM and possibly
a failure of the AGN feedback to prevent the catastrophic cooling of the ICM
(as discussed in Li & Bryan (2014) and Vernaleo & Reynolds (2006)). In
previous work, this problem is addressed by including jet re-orientation,
either by enforcing jet precession around a fixed axis (Li & Bryan 2014; Meece
et al. 2017) or by computing the jet axis following a spin-driven approach
(Beckmann et al. 2019, 2022). In Ehlert et al. (2023), the cold material is
found to deflect the gas, inflating large bubbles that are eventually no
longer aligned with the injection axis. In this section, we explore the
effects of the small-scale jet-filament interaction on the overall morphology
of the outflowing jets and discuss the morphology of the outflows.

###  6.1 Deflection of the jets by cold gas

We find numerous occurrences of interactions between the infalling cold gas
and the outflowing AGN jets. Since the density contrast between the injected
hot material and the clouds is up to $\sim$ 6 orders of magnitude the light
jets are easily deflected by the filaments.

In Fig. 13 we show a volume-rendering of the inner (20 kpc)3 of our cluster.
The yellow and white components represent the isodensity surfaces for $3\times
10^{-3}\,\rm{cm}^{-3}$ and $2\times 10^{-2}\,\rm{cm}^{-3}$ and trace the
outflowing AGN jet material. The blue component traces the cold gas of density
$20\,\rm{cm}^{-3}$. Interactions of the cold gas with the hot material
accelerated by the AGN jet lead to the formation and propagation of cavities
in all directions relative to the vertical jet axis. In some cases the
distribution of the infalling cold gas splits both jets simultaneously,
leading to the formation of a second pair of cavities (right panel).

In Fig. 14 we show a projection of the density-weighted temperature (top
panels) and passive jet scalar (bottom panels). The image sequence in Fig. 13
shows the formation of the cavity visible in the left panel. The projection
depth of 250 pc is set to the typical jet diameter in the inner few kpc of the
box, allowing us to unambiguously visualize the interaction between the cold
gas and the jets. The velocity stream lines are overlaid on the temperature
maps and the round marker indicates the position of the black hole. As shown
in the temperature maps, several cold filaments converge towards the upper jet
(visible in blue on the first first column). The top left filament starts to
interact with the jet in the second column and leads to a splitting of the
jet. A part of the hot material is strongly deflected and starts carving a
cavity orthogonal to the main jet axis. The filament is shattered from its
interaction with the jet, and breaks into smaller clumps, visible in the last
two column. By $t=720$ Myr all the remaining clumps in the way of the jet have
been destroyed.

Figure 14: Density-weighted projections of the temperature field (top panels)
and passive jet scalar (bottom panels), overlaid with the velocity streamlines
of the gas. The depth of projection is 250 pc. The circular marker indicates
the position of the SMBH particle and the injection axis of the jet is the
vertical axis of the image. The interaction of infalling filaments with the
AGN outflows leads to a deflection of the upper jet and the formation of a
cavity propagating in a direction orthogonal to the main jet axis. (Movies
available)

###  6.2 Morphology of the outflows

Cold gas is not always found to lie along the jet axis where it could perturb
the jets. Typical examples for each of these two behaviors are presented in
Fig. 15, showing a density-weighted projection of the temperature field in a
narrow, 500 pc–deep region centered on the midplane of the simulation. In the
left panel the jets exhibit a strongly collimated morphology at distances of
up to $\sim 45-50$ kpc. The right panel shows the same region of the
simulation nearly 100 Myr later. Both jets are bent due to their interaction
with the clumps and several hot spots that formed during previous AGN bursts
are visible. These are propagating in directions different from the jet axis.
Quantifying the relative duration of these two regimes and possible impact on
the ICM heating is beyond the scope of this paper.

Figure 15: Projection of the temperature field in the midplane of the
simulation at two different times. The left panel shows a typical example of
jets propagating in an unperturbed way, staying collimated up to a distance of
45 kpc. In the right panel, several hot spots resulting from previous AGN
bursts can be identified in various directions, as well as jet bending. The
circular marker represents the position of the SMBH particle, and the jet
injection axis is the vertical axis of the image.

Cavities seen in X–ray observations of galaxy clusters show diverse
morphologies and spatial scales. The typical diameter and cluster-centric
distance of cavities can range from a few kpc to extreme values of hundreds of
kpc, such as in MS 0735.6+7421 (Hlavacek-Larrondo et al. 2022). In some cases
more than one pair of bubbles is observed. These are probably associated with
successive AGN bursts, such as in the Perseus cluster (Fabian et al. 2006),
NGC 5813 (Hlavacek-Larrondo et al. 2022) or RBS 797 (Ubertosi et al. 2023). In
order to compare the appearance of our simulation in X-rays to observations we
produce mock images with the yt visualization package. In Fig. 16 we show a
series of X–ray mock images at four different times in our high-resolution MHD
simulation. The top row shows the X–ray luminosity in the 2–7 keV range for
the inner 50 kpc of the box, while the bottom row shows the corresponding
Gaussian Gradient Magnitude (GGM) in the inner 100 kpc of the box. The GGM
image emphasizes faint luminosity variations in the X–ray map and reveals
shocks, cavities and ripples (see Fig. 2 from Sanders et al. (2016) for
comparison). The specific times are selected to showcase the diversity of
X–ray structures found in our simulation.

The first column shows an example of a pair of cavities. Their typical height
and width is roughly 10 kpc and their aspect ratio is close to unity. Such a
morphology is qualitatively similar to the pair of cavities found in the
innermost region of the Perseus cluster (Fabian et al. 2006). There is an
additional, smaller cavity resulting from the deflection of the jet by
infalling filaments close to the AGN, with a typical diameter of 5 kpc. Its
position and relative size can be compared to the second pair of cavities
observed in RBS 797 (see Fig. 7 from Ubertosi et al. (2023)), although it is
unclear whether these have a similar physical origin. Note that the formation
of this small cavity has been presented in Sect. 6.1 and Fig. 14. The shock
wave enveloping the bubbles is also visible in the GGM figure. The third and
fourth columns show typical occurrences of jet propagation without jet
deflection. In that case, the jets propagate along a straight path and carve
narrow channels of under-dense gas into the ICM.

Figure 16: X–ray mock images of our cluster at the native resolution, i.e. not
matched to a particular X-ray instrument. The top row shows the X–ray
luminosity in the 2–7 keV range. The dashed lines indicate the position of the
cavities. The lower row is comparable to Fig. 2 from Sanders et al. (2016) and
shows the corresponding Gaussian Gradient Magnitude (GGM), revealing shocks,
ripples and cavities.

##  7 Discussion

We have achieved spatial resolutions close to a cloud–crushing setup in a
simulation of an idealized cluster with self-regulated AGN feedback. We have
shown that the inclusion of magnetic fields with high–resolution has a
significant impact on the morphology and dynamics of cold gas in a jet-
regulated ICM. Our study however remain imperfect and several question remain
open.

###  7.1 Length of the filaments

We note that our filaments have relatively short lengths compared to the
extents of optical nebulae of some CC clusters. In the Perseus cluster,
H$\alpha$-emitting filaments are observed at distances of up to 50 kpc from
the center of the BCG (Salomé et al. 2006). The typical extents of the cold
features observed in our simulation range between $\sim$5 and $\sim$15 kpc,
closer to the H$\alpha$ filaments observed in the Centaurus cluster (Olivares
et al. 2019). One hypothesis explaining this difference (considering our
initial conditions are meant to be similar to the Perseus cluster) is the
relatively short duration of our simulation compared to other work (Li & Bryan
2014; Qiu et al. 2019). Our simulation lasts for 1.4 Gyr, including the
initial 400 Myr required for gas to cool before the SMBH becomes active. It is
likely that running our simulation for longer will result in the formation of
more extended structures since cooling time increases with radii. For
comparison, the cooling time measured at the end of our simulation at a
distance of 50 kpc from the BCG is about 2 Gyr. It is also possible that the
use of static grid refinement leads to an under-estimate of cooling at larger
radii, where the cells are coarser. However, we note that other simulations
using adaptive mesh refinement also find similar filament lengths (Beckmann et
al. 2019). It is also likely that modeling the relativistic nature of the jets
and resolving gas with increased spatial resolution near the core of the
cluster might result in a more effective coupling between the jets and the
cold clumps, leading to higher clumps launching velocities and larger reached
altitudes from the center of the cluster, thus triggering condensation at
larger radii.

###  7.2 Magnetic field strength in the filaments

It has been a matter of discussion whether the kpc–scale filaments of cold gas
observed in the ICM of clusters and in the circumgalactic medium (CGM) of some
elliptical galaxies are magnetically supported. In Fabian et al. (2008), a
minimum magnetic field strength of 25 $\mu$G is found to be required to
stabilize the filaments observed in NGC 1275 against tidal shear. Similar
field strengths are also obtained from Faraday rotation measurements of the
magnetic field strength in the inner kpc of the Perseus cluster (Taylor et al.
2006). Measurements of the magnetic field strength of the other systems, such
as the cold gas filaments in the Virgo cluster, lead to inferred values of up
to $\sim 100\ \mu$G (Werner et al. 2013). In our MHD run, we find that the
magnetic field strength of the clumps and filaments is particularly high.
While the typical strength values for our filaments are around 80–100 $\mu$G,
we also find individual cells or small clumps with values up to 200–250
$\mu$G. While similar maximum values have been found in MHD cloud-crushing
simulations (Jennings et al. 2023) as well as in self-regulated cluster
simulations (Ehlert et al. 2023), it remains unclear whether such values could
be reached in real CC clusters since observations favor values up to $\sim
50-100\,\mu$G. It is although possible that external sources of pressure, such
as cosmic rays may reduce the overall magnetic field strength in the threads.

###  7.3 Jet propagation and off–axis cavities

As discussed in Sect. 6, our high-resolution MHD run produces many cavities
misaligned with the jets injection axis. Although similar behavior has been
observed and/or discussed in previous work (Wang et al. 2021; Ehlert et al.
2023), we find that these deflected cavities are relatively short-lived and
expand on much shorter distances than the outflows propagating in directions
aligned with the main jet axis, typically up to distances of $\sim 10$ kpc. A
possible explanation for this behavior is our spatial resolution, which allows
us to resolve the fragmentation of the cold filaments into smaller clumps. The
deflection of the jet is then only partial, with a significant fraction of the
material still propagating along the injection axis. We find our low
resolution runs to exhibit behavior closer to that observed in Ehlert et al.
(2023), with larger cavities fully deflected in nearly any direction. We
emphasize that we assume a sub-relativistic propagation of our jets,
hypothesis that might not hold for typical length scales of $\sim 1$ kpc.
Previous studies have shown that relativistic jets with powers above $\sim
10^{44}\ \rm{erg}\,\rm{s}^{-1}$ are less affected by the presence of a clumpy
interstellar medium in disk galaxies (Tanner & Weaver 2022), which might
reduce the deflection by the cold gas even further. However, denser, molecular
gas that would form if we allowed gas to cool below our temperature floor at
$10^{4.2}$ K, could lead to more deflection. We postpone the study of
relativistic effects and the impact of modeling colder gas to further
investigations. Cavities such as those presented in the first two columns of
Fig. 16 are common in our simulations, and are consistent with the typical
sizes of observed X–ray cavities in galaxy clusters (Panagoulia et al. 2014).
Whether the collimated jet morphology presented in the last two columns has a
correspondence in observations is less clear. This is perhaps an artifact of
our setup, however, since ambient bulk motions in real galaxy clusters may
prevent such elongated cavities. It could also indicate additional magnetic
support, a wider jet opening and/or jet re-orientation.

###  7.4 Additional physics

Anisotropic thermal conduction can change the dynamics and morphology of the
multiphase gas. Thermal conduction is mediated by electrons from the hot
surrounding medium that penetrate into the cold filaments. This was studied
analytically, (Cowie & McKee 1977, e.g.) and with simulations (Jennings et al.
2023; Brüggen et al. 2023). In the case of anisotropic conduction, the impact
of conduction has been shown to be limited and depends on the field
orientation and degree of tangling.

Another physical process that has received attention in recent years is
heating by cosmic rays, which can play an important role in clusters. Radio
observation show that cosmic rays pervade at least the inner regions of galaxy
clusters. For instance, Ruszkowski et al. (2017) found that cosmic-ray
streaming was essential for AGN to regulate the cooling of the ICM (see also
Ruszkowski & Pfrommer (2023) for a recent review on the theory of cosmic ray
feedback in galaxies and clusters of galaxies). While the modelling of cosmic
ray feedback is complicated by the largely unknown parameters that govern
cosmic-ray transport, there are clearly avenues for future research in this
area.

The overall morphology of the cold gas could also be impacted by a re-
orientation of the jet injection axis. Beckmann et al. (2019) follow a spin-
driven approach in a purely hydrodynamical setup and find large and fast
variations of the jet orientations, up to $\dot{\theta}_{\rm{jet}}\sim
10^{\circ}\,\rm{Myr}^{-1}$. As a consequence, the ICM is propelled more
isotropically by the jets and the feedback power exhibit strong variations
with peaks up to $10^{49}\,\rm{erg}\,\rm{s}^{-1}$, orders of magnitude above
the maximum AGN power found in our runs. Whether such extreme re-orientations
of the black hole spin axis is achievable in reality is a matter of debate
(Nixon & King 2013).

It has been shown that threads of molecular gas spatially correlates with the
location of the H$\alpha$ emitting nebula observed near the BCG of the Perseus
Cluster, and it is believed that gas in those region could reach temperature
down 10–100 K (Salomé et al. 2006). In our current setup, we impose a
temperature floor of $10^{4}$ K which is likely to limit the total amount of
cold gas formed and its maximum density (McDonald et al. 2010). Below such
temperature threshold, it is suspected that photoionization from young stars
could play a role in the powering of the filaments, and should thus be
included in the cooling calculation, either with explicite radiative transfer
or additional heating source terms in the radiative cooling algorithm.

##  8 Conclusions

The main conclusions of our work are as follows:

  * •

We confirm that magnetic fields have an important impact on the dynamics of
the cold gas. The structures obtained in our MHD simulation resemble the
morphology of the H$\alpha$ filaments observed in cool-core clusters on kpc
scales, such as the Centaurus or Perseus clusters. In the absence of magnetic
fields, simulations produce predominantly clumps instead of filaments with
lower specific angular momentum.

  * •

The jet lifts up cold gas and promotes in-situ cooling. This leads to cold gas
formation at the edges of the hot, outflowing bubbles or jets. In some cases,
up to a few tens of percents of the mass of cold gas is formed during the
uplift, while the rest condenses during the fallback onto the AGN.

  * •

Magnetic pressure largely dominates thermal pressure in the inner regions of
the filaments. The values for plasma-$\beta$ go as low as $10^{-1}-10^{-2}$,
and the magnetic field strengths in the filament core are typically around 100
$\mu$G.

  * •

Cold filaments in our high-resolution MHD run are not monolithic. Their
internal structure is complex and contains hundreds of sub-structures such as
clumps with thicknesses down to a few tens of pc. Using a clump finding
algorithm, we isolate several distinct filaments and study their internal
composition individually. While filaments vary in length and mass, the overall
distribution of their sub-structures is constant across filaments. In
particular, the mass distribution of clumps follows a power-law of slope –
0.6, in agreement with more highly idealized cloud-crushing simulations.

  * •

We identify the cyclic life of some filaments. We find that small clumps are
initially uplifted through jet interaction. The filament structure becomes
visible as the clumps reach their highest altitude. The filamentary structure
are preferably found to be infalling toward the center of the cluster.

  * •

Magnetic field lines inside filaments are aligned with the main filament axis,
while the topology of the magnetic fields in the hot ICM around the filaments
is highly disordered and not necessarily aligned with the filaments. The
properties of ICM turbulence in response to the AGN is the subject of a
separate, forthcoming publication.

  * •

Jet-cloud interactions contribute to the formation of large cavities.
Infalling threads and clumps of cold gas lying along the main jet axis
efficiently deflect and scatter the hot gas, favoring the inflation of
bubbles. While jet-cloud interactions are frequent, we also find transient AGN
bursts with no deflection. In that case the jet remains collimated out to
distances of tens of kpc and its morphology resembles an FR I galaxy.

  * •

Off-axis cavities are ubiquitous in our MHD runs. The large density contrast
between the cold gas and the jet material leads to the formation of cavities
propagating in any direction, including orthogonal to the jet axis. These
cavities are relatively small, with a typical diameter of a few kpc. This is
likely due to the fact that most of the filaments are too narrow to deflect
the jet entirely, and thus only a fraction of the injected material
participates in this process. In purely hydrodynamical runs cavities are
oriented along the jet axis.

###### Acknowledgements.

MF thanks Yu Qiu, Ricarda Beckmann, Marine Prunier, Deovrat Prasad and Ben
Wibking for insightful discussions, as well as Jacob Shen and Evan P. O’Connor
for providing help with visualization packages.
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 101030214.
MB acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum
Universe” – 390833306 and project number 443220636 (DFG research unit FOR
5195: ”Relativistic Jets in Active Galaxies”).
BWO acknowledges support from NSF grants #1908109 and #2106575, NASA ATP
grants NNX15AP39G and 80NSSC18K1105, and NASA TCAN grant 80NSSC21K1053.
The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V.
(www.gauss-centre.eu) for funding this project by providing computing time
through the John von Neumann Institute for Computing (NIC) on the GCS
Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).
This research used resources of the Oak Ridge Leadership Computing Facility at
the Oak Ridge National Laboratory, which is supported by the Office of Science
of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. These
resources were provided by as part of the DOE INCITE Leadership Computing
Program under allocation AST-146 (PI: Brian O’Shea).
All simulations were performed using the public MHD code AthenaPK, which makes
use of the Kokkos (Carter Edwards et al. 2014) library and the Parthenon
adaptive mesh refinement framework (Grete et al. 2023). All data analysis was
performed with yt (Turk et al. 2011, 2024), Matplotlib (Hunter 2007), Numpy
(Harris et al. 2020), Seaborn (Waskom 2021) and swiftascmaps. We thank their
authors for making these software and packages publicly available.

## References

  * Allen (2000) Allen, S. W. 2000, MNRAS, 315, 269
  * Anand et al. (2022) Anand, G. S., Tully, R. B., Rizzi, L., Riess, A. G., & Yuan, W. 2022, ApJ, 932, 15
  * Andreon (2010) Andreon, S. 2010, MNRAS, 407, 263
  * Beckmann et al. (2019) Beckmann, R. S., Dubois, Y., Guillard, P., et al. 2019, A&A, 631, A60
  * Beckmann et al. (2022) Beckmann, R. S., Dubois, Y., Pellissier, A., et al. 2022, A&A, 665, A129
  * 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 & Hensler (1989) Boehringer, H. & Hensler, G. 1989, A&A, 215, 147
  * Boehringer et al. (1993) Boehringer, H., Voges, W., Fabian, A. C., Edge, A. C., & Neumann, D. M. 1993, MNRAS, 264, L25
  * Böhringer & Werner (2010) Böhringer, H. & Werner, N. 2010, A&A Rev., 18, 127
  * Bonamente et al. (2001) Bonamente, M., Lieu, R., & Mittaz, J. P. D. 2001, ApJ, 546, 805
  * Brüggen et al. (2023) Brüggen, M., Scannapieco, E., & Grete, P. 2023, ApJ, 951, 113
  * Carter Edwards et al. (2014) Carter Edwards, H., Trott, C. R., & Sunderland, D. 2014, Journal of Parallel and Distributed Computing, 74
  * Cavagnolo et al. (2009) Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, ApJS, 182, 12
  * Churazov et al. (2003) Churazov, E., Forman, W., Jones, C., & Bohringer, H. 2003, The Astrophysical Journal, 590, 225–237
  * Cowie & McKee (1977) Cowie, L. L. & McKee, C. F. 1977, ApJ, 211, 135
  * Das & Gronke (2024) Das, H. K. & Gronke, M. 2024, MNRAS, 527, 991
  * de Vaucouleurs et al. (1991) de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Herold G., J., et al. 1991, Third Reference Catalogue of Bright Galaxies
  * Dedner et al. (2002) Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645
  * Dunn & Fabian (2006) Dunn, R. J. H. & Fabian, A. C. 2006, MNRAS, 373, 959
  * Ehlert et al. (2023) Ehlert, K., Weinberger, R., Pfrommer, C., Pakmor, R., & Springel, V. 2023, MNRAS, 518, 4622
  * Fabian (1994) Fabian, A. C. 1994, ARA&A, 32, 277
  * Fabian et al. (2008) Fabian, A. C., Johnstone, R. M., Sanders, J. S., et al. 2008, Nature, 454, 968–970
  * Fabian et al. (2011) Fabian, A. C., Sanders, J. S., Allen, S. W., et al. 2011, Monthly Notices of the Royal Astronomical Society, 418, 2154–2164
  * Fabian et al. (2006) Fabian, A. C., Sanders, J. S., Taylor, G. B., et al. 2006, MNRAS, 366, 417
  * Fabian et al. (2016) Fabian, A. C., Walker, S. A., Russell, H. R., et al. 2016, MNRAS, 461, 922
  * Fabian et al. (2017) Fabian, A. C., Walker, S. A., Russell, H. R., et al. 2017, MNRAS, 464, L1
  * Gaspari et al. (2011) Gaspari, M., Melioli, C., Brighenti, F., & D’Ercole, A. 2011, MNRAS, 411, 349
  * Gaspari et al. (2013) Gaspari, M., Ruszkowski, M., & Oh, S. P. 2013, Monthly Notices of the Royal Astronomical Society, 432, 3401–3422
  * Gaspari et al. (2013) Gaspari, M., Ruszkowski, M., & Oh, S. P. 2013, MNRAS, 432, 3401
  * Gómez et al. (2022) Gómez, J. L., Traianou, E., Krichbaum, T. P., et al. 2022, ApJ, 924, 122
  * Gonzalez et al. (2013) Gonzalez, A. H., Sivanandam, S., Zabludoff, A. I., & Zaritsky, D. 2013, ApJ, 778, 14
  * Grete et al. (2023) Grete, P., Dolence, J. C., Miller, J. M., et al. 2023, The International Journal of High Performance Computing Applications, 37, 465
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359
  * Hlavacek-Larrondo et al. (2022) Hlavacek-Larrondo, J., Li, Y., & Churazov, E. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, 5
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Jennings et al. (2023) Jennings, F., Beckmann, R. S., Sijacki, D., & Dubois, Y. 2023, MNRAS, 518, 5215
  * Kino et al. (2015) Kino, M., Takahara, F., Hada, K., et al. 2015, ApJ, 803, 30
  * Laganá et al. (2013) Laganá, T. F., Martinet, N., Durret, F., et al. 2013, A&A, 555, A66
  * Lalakos et al. (2024) Lalakos, A., Tchekhovskoy, A., Bromberg, O., et al. 2024, ApJ, 964, 79
  * Li & Bryan (2014) Li, Y. & Bryan, G. L. 2014, ApJ, 789, 54
  * Liska et al. (2020) Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656
  * Mathews et al. (2006) Mathews, W. G., Faltenbacher, A., & Brighenti, F. 2006, ApJ, 638, 659
  * McCourt et al. (2018) McCourt, M., Oh, S. P., O’Leary, R., & Madigan, A.-M. 2018, MNRAS, 473, 5407
  * McCourt et al. (2012) McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, MNRAS, 419, 3319
  * McDonald et al. (2012) McDonald, M., Veilleux, S., & Rupke, D. S. N. 2012, The Astrophysical Journal, 746, 153
  * McDonald et al. (2010) McDonald, M., Veilleux, S., Rupke, D. S. N., & Mushotzky, R. 2010, ApJ, 721, 1262
  * Meece et al. (2017) Meece, G. R., Voit, G. M., & O’Shea, B. W. 2017, ApJ, 841, 133
  * Miyoshi & Kusano (2005) Miyoshi, T. & Kusano, K. 2005, Journal of Computational Physics, 208, 315
  * Mohapatra et al. (2022) Mohapatra, R., Jetti, M., Sharma, P., & Federrath, C. 2022, MNRAS, 510, 3778
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Nixon & King (2013) Nixon, C. & King, A. 2013, ApJ, 765, L7
  * Olivares et al. (2019) Olivares, V., Salome, P., Combes, F., et al. 2019, A&A, 631, A22
  * Panagoulia et al. (2014) Panagoulia, E. K., Fabian, A. C., Sanders, J. S., & Hlavacek-Larrondo, J. 2014, MNRAS, 444, 1236
  * Park et al. (2019) Park, J., Hada, K., Kino, M., et al. 2019, ApJ, 871, 257
  * Park & Trippe (2017) Park, J. & Trippe, S. 2017, ApJ, 834, 157
  * Peres et al. (1998) Peres, C. B., Fabian, A. C., Edge, A. C., et al. 1998, MNRAS, 298, 416
  * Peterson et al. (2001) Peterson, J. R., Paerels, F. B. S., Kaastra, J. S., et al. 2001, Astronomy & Astrophysics, 365, L104–L109
  * Prasad et al. (2020) Prasad, D., Voit, G. M., O’Shea, B. W., & Glines, F. 2020, ApJ, 905, 50
  * Qiu et al. (2019) Qiu, Y., Bogdanović, T., Li, Y., Park, K., & Wise, J. H. 2019, ApJ, 877, 47
  * Russell et al. (2017) Russell, H. R., McDonald, M., McNamara, B. R., et al. 2017, ApJ, 836, 130
  * Ruszkowski & Pfrommer (2023) Ruszkowski, M. & Pfrommer, C. 2023, A&A Rev., 31, 4
  * Ruszkowski et al. (2017) Ruszkowski, M., Yang, H. Y. K., & Reynolds, C. S. 2017, ApJ, 844, 13
  * Salomé et al. (2006) Salomé, P., Combes, F., Edge, A. C., et al. 2006, A&A, 454, 437
  * Sanders et al. (2004) Sanders, J. S., Fabian, A. C., Allen, S. W., & Schmidt, R. W. 2004, MNRAS, 349, 952
  * Sanders et al. (2016) Sanders, J. S., Fabian, A. C., Russell, H. R., Walker, S. A., & Blundell, K. M. 2016, Monthly Notices of the Royal Astronomical Society, 460, 1898–1911
  * Schure et al. (2009) Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751
  * Simionescu et al. (2011) Simionescu, A., Allen, S. W., Mantz, A., et al. 2011, Science, 331, 1576
  * Tanner & Weaver (2022) Tanner, R. & Weaver, K. A. 2022, AJ, 163, 134
  * Taylor et al. (2006) Taylor, G. B., Gugliucci, N. E., Fabian, A. C., et al. 2006, MNRAS, 368, 1500
  * Toro et al. (1994) Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25
  * Townsend (2009) Townsend, R. H. D. 2009, ApJS, 181, 391
  * Tremblay et al. (2018) Tremblay, G. R., Combes, F., Oonk, J. B. R., et al. 2018, ApJ, 865, 13
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, The Astrophysical Journal Supplement Series, 192, 9
  * Turk et al. (2024) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2024, Introducing yt 4.0: Analysis and Visualization of Volumetric Data, https://yt-project.github.io/yt-4.0-paper/
  * Ubertosi et al. (2023) Ubertosi, F., Gitti, M., Brighenti, F., et al. 2023, ApJ, 944, 216
  * Vernaleo & Reynolds (2006) Vernaleo, J. C. & Reynolds, C. S. 2006, ApJ, 645, 83
  * Voit et al. (2017) Voit, G. M., Meece, G., Li, Y., et al. 2017, ApJ, 845, 80
  * Wang et al. (2021) Wang, C., Ruszkowski, M., Pfrommer, C., Oh, S. P., & Yang, H. Y. K. 2021, MNRAS, 504, 898
  * Wang et al. (2020) Wang, C., Ruszkowski, M., & Yang, H. Y. K. 2020, MNRAS, 493, 4065
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021
  * Werner et al. (2013) Werner, N., Oonk, J. B. R., Canning, R. E. A., et al. 2013, ApJ, 767, 153
  * White et al. (1997) White, D. A., Jones, C., & Forman, W. 1997, MNRAS, 292, 419
  * Zhuravleva et al. (2014) Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2014, Nature, 515, 85

##  Appendix A The impact of spatial resolution

To study the effect of spatial resolution on our results, we run additional
simulations with fewer refinement levels whose parameters are summarized in
Table 3. The extent of each remaining refinement level, as well as the initial
conditions are unchanged from one run to an other. The size of the injection
region for the AGN feedback is adapted to the resolution of the cells at the
center of the box, and thus vary from one run to an other (see Sect. 2.4).

Table 3: Extended version of Table 2 including lower resolution test simulations. Name |  | Max. res.
---
[pc]
$f_{K}$ | $f_{T}$ | $f_{M}$ |  | $\langle B(r=0)\rangle$
---
[$\mu$G]
hydro_xlow | 200 | 0.75 | 0.25 | - | -
hydro_low | 100 | 0.75 | 0.25 | - | -
hydro | 25 | 0.75 | 0.25 | - | -
mhd_xlow | 200 | 0.74 | 0.25 | 0.01 | 10
mhd_low | 100 | 0.74 | 0.25 | 0.01 | 10
mhd_mid | 50 | 0.74 | 0.25 | 0.01 | 10
mhd | 25 | 0.74 | 0.25 | 0.01 | 10

Most MHD simulations of AGN feedback have achieved spatial resolutions of no
more than a few hundred pc (Beckmann et al. 2022; Wang et al. 2021; Ehlert et
al. 2023). We find that increasing spatial resolution to less than 100 pc
affects the morphology of the cold gas, its extent, the subsequent
interactions with the AGN jets, and the morphology of the outflows. In Fig.
17, we present the total cold gas mass in the inner 40 kpc of our box as a
function of time for our 7 runs. Considering the analysis performed in this
paper mostly relies on our highest resolution MHD simulation, we mostly focus
on discussing resolution dependence of our MHD runs. For comparison, we also
show the results for purely hydrodynamical runs.

Figure 17: Total amount of cold gas contained in the inner 40 kpc of the
simulated box.

As visible, lower resolution runs produce on average more cold gas than higher
resolution ones, and can vary by up to roughly one order of magnitude (e.g.,
comparing 25 pc and 200 pc MHD runs at $t=1.0$ Gyr). We also note that our 50
pc and 25 pc MHD runs are producing comparable quantities of cold gas with
variations less than a factor of two throughout the whole duration of the test
run, and converges to $\sim 10^{9}\,\rm{M}_{\odot}$ around $t\sim 1$ Gyr.

In Fig. 18, we show the evolution of $r_{80}$, the radius containing 80% of
the cold gas to quantify its spatial extent. Our 50 pc and 25 pc resolution
MHD runs are broadly consistent, with $r_{80}$ converging to $\sim 5$ kpc, in
agreement with previous work (see notably Fig. 8 from Beckmann et al. (2022)).
Interestingly, our 100 pc resolution MHD run seems to lead to a more extended
distribution of cold gas, which we find to be also preferably located along
the jet injection axis.

Figure 18: Evolution of $r_{80}$, the radius containing 80% of the cold gas in
our various runs.

Identifying the exact physical or numerical processes leading to these spatial
resolution dependencies is beyond the scope of this study as numerous possible
explanations could be investigated. In particular, jet properties are not
identical from one resolution run. As the spatial resolution increases, the
size of the injection region is decreased to fulfill $r_{\rm{jet}}=4\,\Delta
x_{\rm{min}}$ and $h_{\rm{jet}}=2\,\Delta x_{\rm{min}}$. As this might result
in jet injection velocities above the velocity ceiling value
$v_{\rm{ceiling}}=0.05\,c$, we enforce energy conservation by increasing the
density of the injected jet material. This leads to jets in our 25 pc runs to
be heavier than in our 200 pc runs by a factor of up to 2–3. This is likely
affecting their propagation properties and the way turbulence is deposited in
the ICM. We however emphasize that the jets remain relatively light in every
run, with densities in the injection region always of the order of $\sim
10^{-26}\,\rm{g}\,\rm{cm^{-3}}$ ($\sim 10^{-2}\,\rm{cm^{-3}}$) at most.

Figure 19: Gas density projection in the inner 30 kpc of our simulated cluster
for each of our simulations.

In Fig. 19, we show projections of the gas density in the inner 30 kpc of our
simulated box for our various runs. Numerous filaments are visible in each of
our MHD runs, although their substructures are unresolved for resolutions
coarser than 100 pc. No clearly filamentary structures such as the ones
presented in Sect. 4.2 are found throughout the duration of our purely
hydrodynamical simulations apart from a few short-lived and relatively small
filaments (i.e., with typical lengths shorter than 5 kpc, thicknesses close to
resolution limit) visible in our low resolutions hydrodynamical runs.