© The Authors 2024

1. Introduction

Galaxy clusters are massive (M ∼ 10¹⁴−10¹⁵ M_⊙) structures resulting from the gravitational collapse of the highest over-densities 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⁷ − 10⁸ K) diffuse plasma sitting in hydrostatic equilibrium in the dark matter halo (Andreon 2010; Gonzalez et al. 2013; Laganá et al. 2013).

The intracluster medium (ICM) is expected to cool down radiatively through bremsstrahlung, recombination, and collisional processes (Boehringer & Hensler 1989). A third to a half of all clusters, usually designated as cool-core (CC) clusters, have a large central over-density 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 centre, with mass rates of up to 10³ M_⊙ yr⁻¹ (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. Observations of CC clusters shed light on the presence of large cavities apparently inflated by active galactic nucleus (AGN) jets (Fabian et al. 2006; Bîrzan et al. 2004). While it is generally accepted that AGNs 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 suggest that the gas inside the inner tens of kiloparsecs of clusters is multi-phase, with temperatures ranging from 10 to 10⁹ K (e.g. Allen 2000; Bonamente et al. 2001). Hα-emitting nebulae of cold, elongated filaments have been widely studied, and are thought to be ubiquitous in CC clusters (McDonald et al. 2012; Olivares et al. 2019). Multi-wavelength surveys of CC clusters show that the filaments have typical lengths ranging from a few kiloparsecs to tens of kiloparsecs, and thicknesses down to ∼50 pc (Fabian et al. 2008, 2016), and that they 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 thermal instability (TI) is a key process contributing to the formation of multi-phase gas (McCourt et al. 2012; Voit et al. 2017). Turbulence, induced either by AGN jets or by environmental effects, drives the formation of over-densities in the ICM. Due to the dependence of radiative cooling on density, these over-densities can grow non-linearly by drawing hot material out of the background ICM, resulting in structures of molecular and ionised gas with electron densities of ≳100 cm⁻³ (McDonald et al. 2012). These structures are then expected to decouple from the ICM and ‘rain’ towards the centre of the cluster, eventually fueling the supermassive black hole (SMBH) in a chaotic manner (Gaspari et al. 2013). It is also suspected that the AGN jets may be able to uplift cold material present in the inner kiloparsecs of the clusters up to altitudes of tens of kiloparsecs (Russell et al. 2017; Tremblay et al. 2018). This uplifted material could act as a seed for condensation and act as a precursor to the formation of massive, infalling filaments.

Various simulation set-ups have been employed to test this model. Hydrodynamic and magnetohydrodynamic (MHD) simulations of stratified or periodic boxes (McCourt et al. 2012; Valentini & Brighenti 2015; Mohapatra et al. 2022) and idealised CC clusters or massive elliptical galaxies (Li & Bryan 2014a; Beckmann et al. 2019; Qiu et al. 2019; Wang et al. 2021; Guo et al. 2023; Ehlert et al. 2023; Guo et al. 2024) have been able to reproduce some of the observed properties of cold filaments. In Li & Bryan (2014a), 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 hundred parsecs) 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 discs that are common in hydrodynamic set-ups (Li & Bryan 2014a; 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 parsecs) are still lacking, and consequences for AGN feedback remain largely unconstrained. In this paper, we present a set of high-resolution, 3D hydrodynamic 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 were performed on a static, nested mesh that was refined down to a maximum spatial resolution of 24.4 pc.

In Sect. 2, we present the numerical set-up 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 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 used the performance-portable ATHENAPK code based on the block-structured refinement framework PARTHENON (Grete et al. 2023). The fluxes were computed using Harten-Lax-van Leer-contact (HLLC) (Toro et al. 1994) and Harten-Lax-van Leer discontinuity (HLLD) (Miyoshi & Kusano 2005) solvers for our hydrodynamic and MHD runs, respectively. Overall, a second-order accurate method was used based on piecewise linear reconstruction, a predictor-corrector van-Leer-type integrator, and, in the MHD case, a generalised Lagrange multiplier-based divergence cleaning approach (Dedner et al. 2002). The basic system of equations is given by

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·[\rho \mathbf{u}]=0,\hfill \\ \hfill & \frac{\partial \rho \mathbf{u}}{\partial t}+\mathrm{\nabla }·[\rho \mathbf{u}\otimes \mathbf{u}+P_{\text{tot}}\mathbb{I}-\mathbf{B}\otimes \mathbf{B}]=0,\hfill \\ \hfill & \frac{\partial e_{\text{tot}}}{\partial t}+\mathrm{\nabla }·[\mathbf{u}(e_{\text{tot}}+P_{\text{tot}})-\mathbf{u}·\mathbf{B})\mathbf{B}]+\mathcal{L}(\rho ,T)=0,\hfill \\ \hfill & \frac{\partial \mathbf{B}}{\partial t}-\mathrm{\nabla }\times [\mathbf{u}\times \mathbf{B}]=0.\hfill \\ \hfill \end{array}\end{array}$$(1)

Here, ρ is the gas density, u is its velocity field, B is the magnetic field vector, and P_tot is the total pressure, including the thermal pressure, P = (1 − γ)e, and the magnetic field pressure, P_mag = B²/2. The gas is assumed to be non-relativistic, and thus we take γ = 5/3. e is the gas internal energy density and e_tot designates the total energy density, e_tot = e + 0.5ρu² + B²/2.

ℒ(ρ, 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 was set to T_floor = 10⁴ K. Finally, given the non-relativistic set-up of our simulations, we also limited the temperature to T <  5 × 10⁹ 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 adopted a set-up akin to that employed by Wang et al. (2021) and Grete et al. (in prep.). A root grid of 512³ cells covering a volume of (800 kpc)³ was statically refined in its central region. The width of cells at refinement level ℓ is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_{\ell }=\frac{L_{\mathrm{box}}}{512\times 2^{\ell }}.\end{array}$$(2)

In our highest-resolution runs, we used ℓ = 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 ∼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 was modelled using a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997):

$$\begin{array}{c}\hfill g_{\text{NFW}}(r)=\frac{\mathcal{G}}{r^2}\frac{M_{\text{NFW}}[ln(1+\frac{r}{R_{\text{NFW}}})-\frac{r}{r+R_{\text{NFW}}}]}{ln(1+c_{\text{NFW}})-\frac{c_{\text{NFW}}}{1+c_{\text{NFW}}}},\end{array}$$(3)

where 𝒢 is the gravitational constant, M_NFW is the NFW profile virial mass, c_NFW its concentration parameter, and R_NFW is a scale radius function of the virial mass and the concentration parameter, defined by the following relations:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill R_{\mathrm{NFW}}& ={\left(\frac{M_{\mathrm{NFW}}}{4\pi {\rho }_{\mathrm{NFW}}[ln(1+c_{\mathrm{NFW}})-c_{\mathrm{NFW}}/(1+c_{\mathrm{NFW}})]}\right)}^{1/3},\hfill \\ \hfill {\rho }_{\mathrm{NFW}}& =\frac{200H_0^2}{8\pi \mathcal{G}}\frac{c_{\mathrm{NFW}}^3}{ln(1+c_{\mathrm{NFW}})-c_{\mathrm{NFW}}/(1+c_{\mathrm{NFW}})}.\hfill \end{array}\end{array}$$(4)

Secondly, a Hernquist potential was used to model NGC1275, the brightest cluster galaxy (BCG), with M_BCG and R_BCG its mass and radius (Hernquist 1990):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill g_{\text{BCG}}(r)& =\frac{\mathcal{G}{M}_{\text{BCG}}}{R_{\text{BCG}}^2}{(1+\frac{r}{R_{\text{BCG}}})}^{-2}.\hfill \end{array}\end{array}$$(5)

Finally, we added a central point mass modelling the effect of the SMBH of mass M_BH. Density, pressure, and temperature profiles were then calculated from this gravitational field so that the ICM obeys a hydrostatic equilibrium. We assume an entropy profile, K(r), defined by

$$\begin{array}{c}\hfill K\equiv \frac{k_{\mathrm{B}}T(r)}{n_{\mathrm{e}}^{2/3}(r)},\end{array}$$(6)

where k_B is the Boltzmann constant, r is the radial distance to the centre 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):

$$\begin{array}{c}\hfill K(r)=K_0+K_{100}{\left(\frac{r}{100\text{kpc}}\right)}^{{\alpha }_K}.\end{array}$$(7)

In order to close the system of equations defining the hydrostatic equilibrium, we fixed the electron density of the cluster to a defined value, n_e(r_ref), at a defined radius, r_ref = 10 kpc.

In order to break the symmetry of the system, we added velocity perturbations in each run. 40 wave vectors with length scales between 12.5 kpc and 50 kpc were sampled using a random seed and their relative contribution was computed from an inverse parabolic spectrum. We took a peak wavelength of 25 kpc and a velocity dispersion of 75 km s⁻¹. We used a similar spectrum to initialise a perturbed magnetic field, whose average local amplitude follows the density profile

$$\begin{array}{c}\hfill ⟨B(r)⟩\propto \rho {(r)}^{2/3}.\end{array}$$(8)

The strength of the initial magnetic fields was chosen such that its average intensity in the inner few kiloparsecs of the box is ∼10 μG (Taylor et al. 2006). The same random seed was used in each of our runs so that the perturbed fields were identical from one simulation run to another. The values adopted for each of these parameters are summarised in Table 1.

2.4. Active galactic nuclei feedback

We modelled 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 centre of the cluster. Accordingly, we followed the chaotic cold accretion model developed in Gaspari et al. (2013) and later implemented in cluster-scale simulations (Li & Bryan 2014a; Meece et al. 2017). We used an AGN triggering model based on the accretion of cold material. We defined the accretion region as a sphere of radius R_acc = 500 pc, and selected the gas with a temperature below a threshold value of $T_{\mathrm{cold},\mathrm{acc}}=5\times {10}^4\mathrm{K}$ as the accretion material. The accretion rate at each time step was then derived by summing mass over all cells containing gas beneath T_cold, acc and within the accretion region:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}={\displaystyle \sum _{r_i<R_{\mathrm{acc}}}}\frac{\rho ({\widehat{\mathbf{r}}}_i)\mathrm{d}V({\widehat{\mathbf{r}}}_i)}{t_{\mathrm{acc}}}\delta (T({\widehat{\mathbf{r}}}_i)<T_{\mathrm{cold},\mathrm{acc}}),\end{array}$$(9)

where $\rho ({\widehat{\mathbf{r}}}_i)$ is the gas density in cell i, and dV is that cell’s volume. Cold gas in these cells was depleted over the accretion timescale t_acc = 5 Myr. This accretion rate was then used to compute the AGN feedback power, given by

$$\begin{array}{c}\hfill {\dot{E}}_{\text{AGN}}=\eta {\dot{M}}_{\mathrm{acc}}{c}^2,\end{array}$$(10)

where η = 10⁻² 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

$$\begin{array}{c}\hfill {\dot{E}}_{X,\text{AGN}}=f_X{\dot{E}}_{\text{AGN}},\end{array}$$(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 modelled the radiative feedback from the AGN by depositing a portion of the AGN feedback power, $f_T·{\dot{E}}_{\mathrm{AGN}}$, through thermal energy within a sphere of radius $R_T=500\mathrm{pc}$. The corresponding power is related to the AGN power (Eq. (10)) through

$$\begin{array}{c}\hfill {\dot{e}}_T(r\le R_T)=\frac{3}{4}\frac{f_T{\dot{E}}_{\mathrm{AGN}}}{\pi R_T^3}.\end{array}$$(12)

Mass was also injected within the same region of radius R_T:

$$\begin{array}{c}\hfill {\dot{\rho }}_T(r\le R_T)=\frac{3}{4}\frac{f_T{\dot{M}}_{\mathrm{acc}}}{\pi R_T^3}.\end{array}$$(13)

2.4.3. Kinetic feedback

We modelled the jets by redistributing material in two cylinders of radius r_jet = 4 Δx_min, thickness h_jet = 2 Δx_min, and with an offset of L_jet = 4 Δx_min with respect to the centre of the box. We positioned these cylinders in a static orientation along a fixed axis. The velocity of the jet was determined by

$$\begin{array}{c}\hfill v_{\mathrm{jet}}=\sqrt{2\eta }c,\end{array}$$(14)

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

$$\begin{array}{c}\hfill {\rho }_{\mathrm{jet}}=\frac{(1-\eta )f_K{\dot{M}}_{\mathrm{acc}}}{2\pi r_{\mathrm{jet}}^2{h}_{\mathrm{jet}}}.\end{array}$$(15)

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

2.4.4. Magnetic feedback

The importance of magnetic fields in the acceleration and collimation of AGN jets has been emphasised 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 Faraday rotation measurements in AGN jets leads to strengths ranging from ${10}^2\mathrm{G}$ in launching regions down to typical ICM strengths of a few μG (Kino et al. 2015; Park et al. 2019; Gómez et al. 2022).

We included 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_jet ≤ |z|≤L_jet + h_jet), following:

$$\begin{array}{c}\hfill {\mathcal{B}}_{\theta }(r)=2{\mathcal{B}}_0\left(\frac{r}{L_B}\right)exp[-{\left(\frac{r}{L_B}\right)}^2],\end{array}$$(16)

in cylindrical coordinates. L_B = 2.5 Δx_min is the characteristic length scale of the toroidal field set so that the injected fields are mostly confined to the launching region and ℬ₀ was dynamically adjusted so that the resulting injected magnetic energy matches the magnetic fraction of the target AGN feedback.

2.5. Stellar feedback

2.5.1. Type Ia supernova feedback

Type Ia supernova (SNIa) feedback was included, following the method described in Prasad et al. (2020). The energy density and mass density injected into a cell at a distance, r, from the cluster centre is

$$\begin{array}{c}\hfill \begin{array}{cc}& {\dot{e}}_{\mathrm{SNIa}}={\mathrm{\Gamma }}_{\mathrm{SNIa}}{E}_{\mathrm{SNIa}}{\rho }_{\mathrm{BCG}}(r),\hfill \\ \hfill & {\dot{\rho }}_{\mathrm{SNIa}}=\alpha {\rho }_{\mathrm{BCG}}(r),\hfill \end{array}\end{array}$$(17)

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

2.5.2. Stellar formation

In the absence of separate star particles, we also instantaneously converted 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{\mathrm{cm}}^{-3}$, a temperature below T <  2 × 10⁴ K, and within a radius of R_acc <  r <  25 kpc was locally converted into heat with an efficiency of 5 × 10⁻⁶.

2.6. Simulation parameters

We ran two high-resolution simulations, whose parameters are summarised in Table 2. In addition to our high-resolution run with a minimum cell width of 24.4 pc, we performed three hydrodynamic 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 ∼400 Myr of our simulation. During that time, the gas properties evolve due to radiative cooling and the initial velocity turbulence. Over-densities build up in the centre 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 the order of ${10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}$. Phases of low AGN power down to $\sim {10}^{42}\mathrm{erg}{\mathrm{s}}^{-1}$ lasting for a few tens of millions of years are visible in the hydrodynamic 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 converges to typically $\sim {10}^9{\mathrm{M}}_{\odot }$ after 1 Gyr, consistent with observational constraints (Salomé et al. 2006; Fabian et al. 2008) and previous simulations (Beckmann et al. 2019; Li & Bryan 2014a; Ehlert et al. 2023). We find that the specific angular momentum of the cold gas is less than in our MHD set-up. A massive cold disc quickly forms in our hydrodynamic run, while we find no such structure in our MHD runs. Such behaviour is 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 hydrodynamic simulations, from the start of the first AGN phase at t ∼ 400 Myr to the end of our runs. Observational constraints for the cooling time profile is also presented by the grey 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 hydrodynamic and MHD set-ups regarding the overall properties of the ICM. The cooling profiles remain broadly stable over time, with values of $2-3\times {10}^8\mathrm{yr}$ in the inner 10 kpc, consistent with observations. The electron density increases in both our runs up to $\sim {10}^{-1}{\mathrm{cm}}^{-3}$. At late times, the MHD runs show an elevated cooling time.

Fig. 1. Evolution of the AGN power (top panel) as a function of time for our highest-resolution hydrodynamic (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 set-up to an other, but shows weaker fluctuations in MHD set-up. While similar quantities of cold gas are formed in each run, its overall specific angular momentum is found to be larger in the hydrodynamic case due to the formation of a massive disc, consistent with previous work (Li & Bryan 2014a; Qiu et al. 2019; Ehlert et al. 2023).

Fig. 2. Evolution of the radial profiles for cooling time and electron density evaluated at various times in our highest-resolution MHD (top) and hydrodynamic (bottom) runs for the hot ICM gas (i.e. T ≥ 105 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

4.1. Cluster visualisation

To begin our exploration of the cold gas properties in our simulated cluster, we present a series of visualisations 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 the order of ∼100 μG are visibly interacting with the outflowing AGN jets, deflecting them and inflating several cavities. In the top and centre right of the figure, we present mock Hα images of these cold structures. These have been produced by integration along the line of sight of the Hα emission rate, j_Hα, from Dong & Draine (2011):

$$\begin{array}{c}\hfill j_{\mathrm{H}\alpha }=2.83\times {10}^{-26}{T}_4^{-0.942-0.031ln(T_4)}{n}_{\mathrm{e}}n({\mathrm{H}}^{+})\mathrm{erg}{\mathrm{cm}}^3{\mathrm{s}}^{-1}{\mathrm{sr}}^{-1},\end{array}$$(18)

Fig. 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α and X-ray images are also shown in the lower right of the image. The top Hα image is untouched, while the centre 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, the 0.3–1.2 keV (red), 1.2–keV (green), and 2–7 keV (blue) bands, and is comparable to Fig. 1 from Fabian et al. (2006). The dashed white line is a contour map from the temperature field, showing the position of the outflowing bubbles.

where $T_4=T/{10}^4\mathrm{K}$ and T, n_e, and n(H⁺) are the temperature, electron density, and ionised hydrogen density of each gas cell in the simulated box, respectively. The second Hα flux has been smoothed with a Gaussian filter 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 represent the 0.3–1.2 keV, 1.2–2 keV, and 2–7 keV bands, respectively.

4.2. Morphology

Previous research has shown that magnetic fields affect the overall morphology of the cold gas. While high-resolution hydrodynamic simulations can produce cold gas in quantities consistent with observations, the morphology of the structures is usually very clumpy (Li & Bryan 2014a) and lacks a coherent, filamentary morphology. Magnetic fields have been found in both cluster and cloud-crushing set-ups to inhibit the formation of clumps and massive discs, 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 hydrodynamic run and the right panel the MHD run. In agreement with previous work (e.g. Ehlert et al. 2023), we find that the magnetic fields have an important impact on the morphology of the cold gas. In the hydrodynamic case, the cold gas condenses into thousands of cloudlets typically consisting of just 1–10 cells. In our MHD set-up, the gas gathers along extended filamentary structures, preferentially infalling towards the centre of the cluster. Each of these structures is dynamic, reaches maximum elongations up to ∼15 kpc, and contains up to $\sim {10}^{8.5}{\mathrm{M}}_{\odot }$ of cold gas.

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

The morphologies observed in our hydrodynamic 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 discs 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.

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, which is indicated in the top left of each panel. The typical filament ranges between 5 and 15 kpc in length, and its morphology can be complex. In particular, all filaments exhibit numerous inner structures of cold clumps and sub-filaments of various thicknesses spanning between ∼800 pc and the resolution limit of the simulation. To study quantitatively the distribution of these sub-structures in and around these five filaments, we used a clump-finding algorithm integrated in the visualisation 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 set-up, we used a minimum density threshold of $n_{\mathrm{clump},\min }=1{\mathrm{cm}}^{-3}$ and a step factor of s = 4 corresponding to the multiplication factor applied to the minimum density threshold at each iteration of the algorithm. The impact of the clump-finding algorithm is discussed in Appendix B. For each system, we limited the search for dense structures to all gas cells contained in a cylinder of radius r_cyl = 2 kpc manually aligned with each filament’s main axis. We also excluded all clumps with fewer than ten cells and with a mean temperature above T_clump, max = 10⁵ K to exclude possible artefacts from the hot ICM. Several cell volume-weighted quantities were calculated for each filament by summing over all gas cells belonging to a given clump. In particular, we derived the magnetic field strength, the ratio of the thermal pressure to magnetic pressure, β (later on referred as plasma, β), and the radial velocity of each clump belonging to each of the five filaments.

Fig. 5. Projection of the gas density for five different filaments in our fiducial MHD run. Each projection was performed at the time when each filament reached 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.

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 in our analysis; that is, each one containing 10 or more cells. The corresponding PDFs are represented by solid blue lines. The number in parentheses corresponds to the total number of clumps with no constraint on the minimum number of cell for each clump. The associated PDFs are also represented by dashed blue lines for comparison.

Fig. 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 is indicated in the legend of the left panel. The first number indicates the number of clumps included in our statistical analysis, while the number in parentheses indicates the total clumps population with no constraint on the 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 dashed bright blue lines. The dashed orange line indicates a dN/dM power law distribution of slope α = −0.6 for comparison. The shaded grey area in the top right plot corresponds to the region of positive radial velocity (i.e. clumps flowing out from the cluster centre).

As is 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 of $M_{\mathrm{clump}}\le {10}^7{\mathrm{M}}_{\odot }$ contributes ∼25% of the overall mass of each filament, while 50% of each filament’s mass is carried by the few clumps with masses of $M_{\mathrm{clump}}\ge {10}^{7.5}{\mathrm{M}}_{\odot }$. The number of such massive clumps 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}^5{\mathrm{M}}_{\odot }$ and ${10}^7-{10}^8{\mathrm{M}}_{\odot }$. The corresponding slope of the clump number distribution, $\mathrm{d}log\mathrm{N}/\mathrm{d}log\mathrm{M}$, is typically of α = −0.6 (see dashed orange line for comparison). While this is consistent with several previous simulation-based works (Jennings et al. 2023; Ramesh & Nelson 2024), steeper spectra with slopes down to ∼ − 2 are also reported in several previous CC cluster simulations (Li & Bryan 2014a, 2014b) and turbulent box set-ups (Gronke et al. 2022; Das & Gronke 2024).

The left and central panels support a picture in which the filaments are clumpy structures, largely dominated in number by light clumps with masses as low as ${10}^4{\mathrm{M}}_{\odot }$, and largely dominated in mass by the contribution of only a few massive clumps with individual masses above ${10}^7-{10}^8{\mathrm{M}}_{\odot }$. We emphasise 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 ten cells, suggesting that the lower mass end of the clump statistics has not yet converged. This is expected and consistent with results from 2D and 3D cloud-crushing simulations that conclude 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. As is visible, the bulk motion of the filaments is preferably infalling, with average velocities of between −100 and −200 km s⁻¹. All five filaments contain clumps characterised by a positive radial velocity, up to +400 km s⁻¹. 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).

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, β. The median values of the magnetic field strength in the clumps sample of our five filaments are typically of Med(B_clumps)∼80 − 100 μG. Some of the clumps have a particularly high inner magnetic field strength, ranging up to ∼200 μG. As is 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 half the clumps.

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

We find a correlation between the magnetic field lines orientation within the filaments and each 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, β, in the mid-plane 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 behaviour is also found in the other filaments we have inspected (see also Fig. 9). 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 ∼130° to reach alignment with the filament.

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

Fig. 9. Sequence of images spanning the whole lifetime of filament system #2 (see Sect. 4.2). The panels show, from left to right: the gas column density, density-weighted radial velocity, density of jet scalar, and density-weighted thermal-to-magnetic pressure ratio, β, overlaid with the magnetic field lines. In each case, the plotting region is aligned with the filament’s main axis. The SMBH is located to the right of the plot and the jet axis has an angle of ∼10° with the main filament axis. Time has been re-scaled such that t = 0 coincides with the last snapshot before the first clumps starts to move upwards along the subsequent filament’s main axis. The times of the snapshots have been selected such that the main steps of the filament evolution are most apparent, and are thus voluntarily unevenly spaced.

5. Formation of the cold filaments

5.1. Upward-moving gas and uplifting

We find evidence of ubiquitous upward-moving gas throughout the entire course of our fiducial MHD simulation. In Fig. 10, we show the evolution of M_cold, up, the total mass of upward-moving cold gas with radial velocities of a certain threshold, defined as

$$\begin{array}{c}\hfill M_{\mathrm{cold},\mathrm{up}}(v_{\mathrm{rad}})={\displaystyle \sum _i}{m}_i\delta (v_{\mathrm{rad},\mathrm{i}}\ge v_{\mathrm{rad}})\delta (T_i\le {10}^5\mathrm{K}),\end{array}$$(19)

Fig. 10. 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 solid orange line indicates the evolution of the AGN feedback power.

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

After about 400 Myr, the total upward-moving gas mass converges to ${10}^{8.4}{\mathrm{M}}_{\odot }$, which is about 15% of the total cold gas 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 less than a million years, while the typical variations in the cold gas mass are smoother and typically of a few million years 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 upward-moving cold mass start to be noticeable. In particular, some peaks and drops in AGN power are correlated with peaks and drops in M_cold, up, especially for radial velocity thresholds ≥300 km s⁻¹.

In Fig. 11, 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, filaments #1, #2, #4, and #5) exhibit an initial phase where clumps move upwards, with average radial velocities ranging from 20–100 km s⁻¹. Since not all clumps are moving upwards and because the region within which the physical quantities were 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 move with radial velocities far beyond these values, up to 400 km s⁻¹. For 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 phase during which the bulk of the filament is characterised by a positive radial velocity lasts for 20–40 Myr in these four filaments, after which the system starts falling back towards the centre of the cluster. Filament #3 exhibits a different behaviour as no clear evidence of upward-moving cold gas is found prior to the formation of the filament through visually inspection. It is likely that this system forms solely from condensation of the ICM.

Fig. 11. Evolution of the average radial velocity of our five filament systems, as a function of their total cold gas mass budget. The shaded grey area indicates positive radial velocities. For three of our filament systems (namely, #1, #2, #4, and #5), the bulk of the filament is moving upwards at early times due to its interaction with the AGN outflows. The radial velocity of the system typically reverses ∼20–40 Myr after the beginning of the positive radial velocity phase. The filaments then continue building up mass while falling back towards the centre of the clusters for ∼20–100 Myr, before being shredded. The time is re-scaled based on visual inspection such that t = 0 coincides with the last snapshot before the formation sequence of each filament starts. For filaments #1, #2, #4, and #5, this corresponds to when the first clumps are visibly starting to move upwards in the direction of what will later on become the filament. For filament #3, there is no visible upward bulk motion as the filament solely forms out of condensation; thus, t = 0 is defined here as the last snapshot before infalling clumps visibly gather to form the filament.

5.2. Evolution of the filaments structure

A sequence of images showing the formation and evolution of filament system #2 is shown in Fig. 9. The panel shows, from left to right: column gas density, density-weighted radial velocity, and a mid-plane slice of the plasma, β, 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 ∼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 centre of the cluster, and other ones exhibiting a positive radial velocity above 200 km s⁻¹, 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 ∼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 centre of the SMBH. Interestingly, the clouds located at the tip of the expanding structure have velocities above 300 km s⁻¹ at t = 40 Myr. Considering the gravitational field imposed in our simulation (see Sect. 2.3), a simple integration of the equations of motion shows that the clumps should come to rest at t ∼ 40 Myr if their motion was ballistic. Hence, the presence of some clumps with radial velocities above 300 km s⁻¹ after several tens of millions of years suggests that the interaction between the AGN outflows and the cold clouds lasts for several tens of millions of years, sustaining momentum transfer from the jets to the cold clouds.

At t ∼ 100 Myr, a clear filamentary structure starts to appear, which extends up to distances of ∼15 kpc from the centre of the cluster. The bulk of this thread is magnetically supported, as β drops by nearly two orders of magnitude between the ICM and the core of the filament, where its reaches values of β ∼ 10⁻¹ − 10⁻². At later times, the filament falls back towards the centre of the cluster with radial velocities down to – (200–300) km s⁻¹. 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.

The last column of Fig. 9 suggests that the formation of the filament leads to a reshaping of the magnetic field’s topology, which seems to transition between a relatively disorganised state to a topology where field lines in and around the filament are aligned with the main filament axis. To study this phenomenon more quantitatively, we present Fig. 12, which depicts the cell volume-weighted distribution of the angle between the magnetic field, B, and the unit vector, ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$, for filament #2, where ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$ indicates the main axis of the filament, oriented outwards from the centre of the cluster. The dashed curve indicates the distribution averaged on the snapshots corresponding to the phase where the average radial velocity of the cold gas is positive (see Fig. 11). The solid curves indicate the distribution in the later phase corresponding to the infall of the filament. The purple (orange) curves indicate the distribution for the cold (hot) phase, respectively. All of the considered cells are restrained here to the same region as for the analysis performed in Fig. 6: a cylinder aligned with ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$ of radius r_cyl = 2 kpc. In the first phase, there is no clear trend visible in the distribution of the magnetic field orientation. At later times (solid curves), the distribution moves towards $cos(\widehat{\theta })=+1$ due to a preferential orientation of the magnetic field in both the hot and cold phases towards the main filament axis.

Fig. 12. Cell volume-weighted distribution of the angle between the magnetic field, B, and the unit vector, ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$, for filament #2, with ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$ indicating the main axis of the filament, pointing outwards from the cluster’s centre. The dashed line indicates the distribution averaged over the snapshots corresponding to the phase where the average radial velocity of the filament is positive, while the solid curved is averaged on the subsequent infalling part (see Fig. 11). The thin curves in the background indicate the successive individual curves for the cold phase, with darker colours associated with later times.

In Fig. 13, we show the evolution of several physical quantities characterising this same filament system, #2, presented in Fig. 9 (see also left panel of Fig. 11). 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 centre of mass to the centre of the cluster (dashed line). The other panels show 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 referred here with respect to the begin of the formation of the filament; that is, when the first clumps start to move upwards (see Fig. 9). 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 centre 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. We find no clear evolution in the magnetic pressure inside the clumps. The median value of β remains remarkably constant, around ∼10⁻¹, and the proportion of the number of clumps with a dominant magnetic pressure (i.e. with β <  1) is > 90% at all times.

Fig. 13. Evolution of the properties of the filament system (#2) (see also Fig. 9) at five equally spaced times. Top left panel: Evolution of the total cold mass budget of the system (solid line) and position of the centre-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 shaded grey 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, β, 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–ten–of cells for a clump to be included in the sample.

As is 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\mathrm{km}{\mathrm{s}}^{-1}$. Even when the overall bulk motion motion of the filaments falls back towards the centre of the cluster (i.e. as of t ∼ 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 was discussed above, the formation of filaments in our MHD run results from the combination of both upward-moving cold clumps and in situ condensation of the ICM. Since our current set-up does not include tracer particles, we cannot unambiguously disentangle the contribution of these two phenomena and we postpone such analysis to a subsequent study.

6. Impact of cold gas on the active galactic nucleus outflows

Cavities in galaxy clusters have been widely studied using X-ray observations. Typically, pairs of cavities with radii of ∼10 kpc are found up to a few tens of kiloparsecs away from the central galaxy. These are sometimes filled with radio emission associated with AGN radio jets that 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 works, it has been discussed that high-resolution simulations (typically with maximally refined cells smaller than ∼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 kiloparsecs. 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 was discussed in Li & Bryan (2014a) and Vernaleo & Reynolds (2006)). In previous works, this problem has been addressed by including jet re-orientation, either by enforcing jet precession around a fixed axis (Li & Bryan 2014a; 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 ∼6 orders of magnitude, the light jets are easily deflected by the filaments.

In Fig. 14, we show a volume-rendering of the inner (20 kpc)³ of our cluster. The yellow and white components represent the isodensity surfaces for $3\times {10}^{-3}{\mathrm{cm}}^{-3}$ and $2\times {10}^{-2}{\mathrm{cm}}^{-3}$ and trace the outflowing AGN jet material. The blue component traces the cold gas of density $20{\mathrm{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).

Fig. 14. 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}{\mathrm{cm}}^{-3}$ and $2\times {10}^{-2}{\mathrm{cm}}^{-3}$, respectively, and trace the outflowing AGN jet material. The blue component traces the cold gas of density $20{\mathrm{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 centre of each panel represents the position of the SMBH particle (movies available online).

In Fig. 15, we show a projection of the density-weighted temperature (top panels) and passive jet scalar (bottom panels). The image sequence in Fig. 14 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 kiloparsecs of the box, allowing us to unambiguously visualise 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 is shown in the temperature maps, several cold filaments converge towards the upper jet (visible in blue in the first column). The top left filament starts to interact with the jet in the second column and leads to a splitting of the jet. 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 columns. By t = 720 Myr, all the remaining clumps in the way of the jet have been destroyed.

Fig. 15. 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 online).

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 behaviours are presented in Fig. 16, showing a density-weighted projection of the temperature field in a narrow, 500 pc-deep region centred on the mid-plane of the simulation. In the left panel, the jets exhibit a strongly collimated morphology at distances of up to ∼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 their possible impact on the ICM heating is beyond the scope of this paper.

Fig. 16. Projection of the temperature field in the mid-plane 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 can 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 kiloparsecs to extreme values of hundreds of kiloparsecs, 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 produced mock images with the YT visualisation package. In Fig. 17, 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 emphasises 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 were selected to showcase the diversity of X-ray structures found in our simulation.

Fig. 17. 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 GGM, revealing shocks, ripples, and cavities.

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. It should be noted that the formation of this small cavity has been presented in Sect. 6.1 and Fig. 15. 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.

7. Discussion

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

7.1. Upward-moving gas and uplifting

Various channels for the formation of the filaments have been discussed in the literature. Li & Bryan (2014b) and Voit et al. (2017) discuss a picture in which low-entropy hot gas is uplifted through interaction with the AGN jets. This uplifted hot gas then condenses and builds up mass, while decelerating and falling back towards the centre of the accretion zone. In Prasad et al. (2015), the jets are found to effectively uplift clumps of cold gas, which might rise up to altitudes of a few tens of kiloparsecs. Since our simulations do not include tracer particles, we cannot unambiguously quantify the relative importance of these two formation channels and our interpretation can only rely on simple considerations. Throughout the whole duration of our fiducial MHD run, the inner 3–4 kpc of the simulated cluster constantly contains a myriad of clumps and short filaments (with a typical length of a few hundred parsecs) of cold (T <  10⁵ K) gas. The first few kiloparsecs of the path of the jets thus consist of a strongly multi-phase medium, and interactions between the jets and cold clumps of gas are frequent, as is discussed in Sect. 6. Looking at the kinematics of cold gas in the direction of propagation of the jets, we found ubiquitous cold gas structures with positive radial velocities of up $\sim {10}^3\mathrm{km}{\mathrm{s}}^{-1}$ at distances as small as a few hundred parsecs away from the kinetic jet injection region. This is an indication that the presence of fast outflows of cold gas is likely the result of a physical uplifting by the jets. While there remains a possibility that the origin of these fast upward-moving clumps could be explained by the condensation of low-entropy hot gas, this condensation then must happen on distances as small as a few hundred parsecs, roughly an order of magnitude below results from previous work (Li & Bryan 2014b).

7.2. Length of the filaments

We note that our filaments have relatively short lengths compared to the extents of optical nebulae in some CC clusters. In the Perseus cluster, Hα-emitting filaments are observed at distances of up to 50 kpc from the centre of the BCG (Salomé et al. 2006). The typical extents of the cold features observed in our simulation range between ∼5 and ∼15 kpc, closer to the Hα filaments observed in the Centaurus cluster (Olivares et al. 2019). One possible explanation is 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 modelling the relativistic nature of the jets and resolving gas with increased spatial resolution near the core of the cluster will result in a more effective coupling between the jets and the ambient gas, leading to higher launching velocities and higher altitudes being reached from the centre of the cluster, potentially triggering condensation at larger radii.

7.3. Magnetic field strength in the filaments

It has been a matter of discussion whether the kiloparsec-scale filaments of cold gas observed in the ICM of clusters and in the circumgalactic medium of some elliptical galaxies are magnetically supported. In Fabian et al. (2008), a minimum magnetic field strength of 25 μG is found to be required to stabilise 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 kiloparsecs 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 ∼100 μ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 μG, we also find individual cells or small clumps with values up to 200–250 μ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 favour values of up to ∼50 − 100 μG. It is possible, though, that external sources of pressure such as cosmic rays may reduce the overall magnetic field strength in the threads.

7.4. Jet propagation and off-axis cavities

As is discussed in Sect. 6, our high-resolution MHD run produces many cavities that are misaligned with the jets injection axis. Although a similar behaviour has been observed and/or discussed in previous works (Li & Bryan 2014a; Qiu et al. 2019; Wang et al. 2020, 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 ∼10 kpc. A possible explanation for this behaviour 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 behaviour closer to that observed in Ehlert et al. (2023), with larger cavities fully deflected in nearly any direction. We emphasise that we assume a sub-relativistic propagation of our jets, a hypothesis that might not hold for typical length scales of ∼1 kpc. Previous studies have shown that relativistic jets with powers above $\sim {10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}$ are less affected by the presence of a clumpy interstellar medium in disc 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 modelling colder gas to further investigations. Cavities such as those presented in the first two columns of Fig. 17 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 artefact of our set-up, 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.5. Additional physics

Anisotropic thermal conduction can change the dynamics and morphology of the multi-phase gas. Thermal conduction is mediated by electrons from the hot surrounding medium that penetrate into the cold filaments. This has been studied analytically (e.g. Cowie & McKee 1977) 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 observations 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 AGNs 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 hydrodynamic set-up and find large and fast variations in the jet orientations, up to ${\dot{\theta }}_{\mathrm{jet}}\sim 10°{\mathrm{Myr}}^{-1}$. As a consequence, the ICM is propelled more isotropically by the jets and the feedback power exhibits strong variations with peaks up to ${10}^{49}\mathrm{erg}{\mathrm{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 are achievable in reality is a matter of debate (Nixon & King 2013).

It has been shown that threads of molecular gas spatially correlate with the location of the Hα-emitting nebula observed near the BCG of the Perseus cluster, and it is believed that gas in those regions could reach temperatures as low as 10–100 K (Salomé et al. 2006). In our current set-up, we impose a temperature floor of 10⁴ K, which is likely to limit the total amount of cold gas formed and its maximum density (McDonald et al. 2010). Below such a temperature threshold, it is suspected that photoionisation from young stars could play a role in the powering of the filaments, and should thus be included in the cooling calculation, either with explicit 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α filaments observed in CC clusters on kiloparsec scales, such as the Centaurus or Perseus clusters. In the absence of magnetic fields, simulations predominantly produce clumps instead of filaments with a lower specific angular momentum.

• The jet lifts up gas and promotes in situ cooling. This leads to cold gas formation at the edges of the hot, outflowing bubbles or jets. As much as a few tens of percent of the mass of some filaments is carried by upward-moving cold clumps with velocities of up to $\sim {10}^3\mathrm{km}{\mathrm{s}}^{-1}$, while the remainder condenses out of the hot phase during the fall back onto the AGN.

• Magnetic pressure largely dominates thermal pressure in the inner regions of the filaments. The values for plasma-β go as low as 10⁻¹ − 10⁻², and the magnetic field strengths in the filament core are typically around 100 μ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 parsecs. Using a clump-finding algorithm, we isolated several distinct filaments and studied their internal compositions 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 relatively shallow power law with a slope of −0.6, in agreement with more highly idealised cloud-crushing simulations.

• We have identified the cyclic life of some filaments. We find that small clumps are initially moving upwards, likely through jet interaction. The filament structure becomes visible as the clumps reach their highest altitude. The filamentary structures are preferably found to be infalling towards the centre of the cluster.

• The formation of filaments lead to a re-organisation of the magnetic field topology in the hot and cold phases. While the topology of the magnetic field has no preferential orientation prior to the filaments’ formation, field lines gradually align with the main filament axis during the formation sequence and point outwards from the cluster’s centre.

• 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, favouring 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 kiloparsecs and its morphology resembles that of 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 kiloparsecs. 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 hydrodynamic runs, cavities are oriented along the jet axis.

Data availability

Movies associated to Figs. 8, 14, and 15 are available at https://www.aanda.org

Acknowledgments

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 visualisation 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. (http://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 https://github.com/parthenon-hpc-lab/athenapk, which makes use of the https://github.com/kokkos/kokkos (Carter Edwards et al. 2014) library and the https://github.com/parthenon-hpc-lab/parthenon adaptive mesh refinement framework (Grete et al. 2023). All data analysis was performed with https://yt-project.org/ (Turk et al. 2011, 2024), https://matplotlib.org/ (Hunter 2007), https://numpy.org/ (Harris et al. 2020), https://seaborn.pydata.org/ (Waskom 2021) and https://github.com/JBorrow/swiftascmaps. We thank their authors for making these software and packages publicly available.

References

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 summarised in Table A.1. 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 centre of the box, and thus vary from one run to an other (see Sect. 2.4).

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. A.1, 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 hydrodynamic runs.

Fig. A.1. 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{\mathrm{M}}_{\odot }$ around t ∼ 1 Gyr.

In Fig. A.2, we show the evolution of r₈₀, 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₈₀ converging to ∼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.

Fig. A.2. Evolution of r80, 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_jet = 4 Δx_min and h_jet = 2 Δx_min. As this might result in jet injection velocities above the velocity ceiling value v_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 emphasise that the jets remain relatively light in every run, with densities in the injection region always of the order of $\sim {10}^{-26}\mathrm{g}{\mathrm{cm}}^{-3}$ ($\sim {10}^{-2}{\mathrm{cm}}^{-3}$) at most.

Fig. A.3. Gas density projection in the inner 30 kpc of our simulated cluster for each of our simulations.

In Fig. A.3, 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 sub-structures 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 hydrodynamic 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 hydrodynamic runs.

Appendix B: The impact of the clump-finding algorithm

As stated in Sect. 4.2, the clump mass distribution around filaments in our MHD simulations is found to be characterised by a power law ranging from ∼10⁵ to $\sim {10}^8{\mathrm{M}}_{\odot }$, with a slope α ∼ −0.7. While this result is consistent with results from some previous simulation based work (Jennings et al. 2023; Ramesh & Nelson 2024), other work also find steeper distribution, with slopes between α ∼ −1.5 and α ∼ −2 (Li & Bryan 2014b; Gronke et al. 2022; Das & Gronke 2024). In particular, the clump-finding algorithm used in Gronke et al. (2022) and Das & Gronke (2024) is different than ours, as it uses Scipy’s ndimage (Virtanen et al. 2020) library to extract temperature contours, while our analysis is performed with yt clump-finding algorithm which performs a hierarchical search for structures recursively between two density extrema, similar to the analysis performed in Li & Bryan (2014b). To assess the sensitivity of our results on the choice of the clump-finding algorithm, we present in Fig. B.1 the clump distribution for the inner 10 kpc of our simulated box for our fiducial MHD run. The blue and orange curves are obtained with the yt clump finder. For the blue one, the standard hierarchical search of disconnected structures is performed such as described in Sect. 4.2. For the orange curve, the selection is performed with yt based on a temperature criteria, such that clumps are defined as disconnected structures with a T = 10⁵ K contour. This second method is theoretically identical to the third one which is comparable to the method used in (Gronke et al. 2022) and Das & Gronke (2024). Here, the gas of temperature T ≤ 10⁵ K is isolated and Scipy’s ndimage (Virtanen et al. 2020) is used to find the clumps. As visible, the second and last method are in a nearly perfect agreement, while the number of clumps is a factor of ∼2 larger using the method employed for the analysis in this paper. This is due to the fact that clumps are not defined the same way between these two methods. Using a density criteria and a hierarchical arrangement of clumps allow to distinguish structures which are not distinguished by the scipy method. The slope of the distribution is consistent across methods.

Fig. B.1. Comparison of the clump distribution in the inner 10 kpc of our fiducial MHD simulation at t = 960 Myr with three different method. The blue and orange curves shows the distribution of the clumps with yt (Turk et al. 2011) based on density and temperature criteria. The red curve indicates the clump distribution using Scipy’s ndimage (Virtanen et al. 2020), and is comparable to the method used in (Gronke et al. 2022) and Das & Gronke (2024). The inclination of a power law of slope α ∼ −0.6 is indicated for comparison.

All Tables

All Figures

Fig. 1. Evolution of the AGN power (top panel) as a function of time for our highest-resolution hydrodynamic (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 set-up to an other, but shows weaker fluctuations in MHD set-up. While similar quantities of cold gas are formed in each run, its overall specific angular momentum is found to be larger in the hydrodynamic case due to the formation of a massive disc, consistent with previous work (Li & Bryan 2014a; Qiu et al. 2019; Ehlert et al. 2023).

In the text

	Fig. 2. Evolution of the radial profiles for cooling time and electron density evaluated at various times in our highest-resolution MHD (top) and hydrodynamic (bottom) runs for the hot ICM gas (i.e. T ≥ 105 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).
In the text

Fig. 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α and X-ray images are also shown in the lower right of the image. The top Hα image is untouched, while the centre 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, the 0.3–1.2 keV (red), 1.2–keV (green), and 2–7 keV (blue) bands, and is comparable to Fig. 1 from Fabian et al. (2006). The dashed white line is a contour map from the temperature field, showing the position of the outflowing bubbles.

In the text

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

	Fig. 5. Projection of the gas density for five different filaments in our fiducial MHD run. Each projection was performed at the time when each filament reached 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.
In the text

Fig. 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 is indicated in the legend of the left panel. The first number indicates the number of clumps included in our statistical analysis, while the number in parentheses indicates the total clumps population with no constraint on the 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 dashed bright blue lines. The dashed orange line indicates a dN/dM power law distribution of slope α = −0.6 for comparison. The shaded grey area in the top right plot corresponds to the region of positive radial velocity (i.e. clumps flowing out from the cluster centre).

In the text

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

	Fig. 8. Projection of the gas density (left) and slice of the plasma, β, overlaid with the magnetic field lines (right) for the filament #3. The projection depth of the density plot is 4 kpc. The dashed grey line around the filament core represents the contour of equi-partition between magnetic and thermal pressure (movies available online).
In the text

Fig. 9. Sequence of images spanning the whole lifetime of filament system #2 (see Sect. 4.2). The panels show, from left to right: the gas column density, density-weighted radial velocity, density of jet scalar, and density-weighted thermal-to-magnetic pressure ratio, β, overlaid with the magnetic field lines. In each case, the plotting region is aligned with the filament’s main axis. The SMBH is located to the right of the plot and the jet axis has an angle of ∼10° with the main filament axis. Time has been re-scaled such that t = 0 coincides with the last snapshot before the first clumps starts to move upwards along the subsequent filament’s main axis. The times of the snapshots have been selected such that the main steps of the filament evolution are most apparent, and are thus voluntarily unevenly spaced.

In the text

	Fig. 10. 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 solid orange line indicates the evolution of the AGN feedback power.
In the text

Fig. 11. Evolution of the average radial velocity of our five filament systems, as a function of their total cold gas mass budget. The shaded grey area indicates positive radial velocities. For three of our filament systems (namely, #1, #2, #4, and #5), the bulk of the filament is moving upwards at early times due to its interaction with the AGN outflows. The radial velocity of the system typically reverses ∼20–40 Myr after the beginning of the positive radial velocity phase. The filaments then continue building up mass while falling back towards the centre of the clusters for ∼20–100 Myr, before being shredded. The time is re-scaled based on visual inspection such that t = 0 coincides with the last snapshot before the formation sequence of each filament starts. For filaments #1, #2, #4, and #5, this corresponds to when the first clumps are visibly starting to move upwards in the direction of what will later on become the filament. For filament #3, there is no visible upward bulk motion as the filament solely forms out of condensation; thus, t = 0 is defined here as the last snapshot before infalling clumps visibly gather to form the filament.

In the text

Fig. 12. Cell volume-weighted distribution of the angle between the magnetic field, B, and the unit vector, ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$, for filament #2, with ${\widehat{\mathbf{r}}}_{\mathrm{fil}}$ indicating the main axis of the filament, pointing outwards from the cluster’s centre. The dashed line indicates the distribution averaged over the snapshots corresponding to the phase where the average radial velocity of the filament is positive, while the solid curved is averaged on the subsequent infalling part (see Fig. 11). The thin curves in the background indicate the successive individual curves for the cold phase, with darker colours associated with later times.

In the text

Fig. 13. Evolution of the properties of the filament system (#2) (see also Fig. 9) at five equally spaced times. Top left panel: Evolution of the total cold mass budget of the system (solid line) and position of the centre-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 shaded grey 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, β, 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–ten–of cells for a clump to be included in the sample.

In the text

Fig. 14. 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}{\mathrm{cm}}^{-3}$ and $2\times {10}^{-2}{\mathrm{cm}}^{-3}$, respectively, and trace the outflowing AGN jet material. The blue component traces the cold gas of density $20{\mathrm{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 centre of each panel represents the position of the SMBH particle (movies available online).

In the text

Fig. 15. 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 online).

In the text

Fig. 16. Projection of the temperature field in the mid-plane 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 can jet bending. The circular marker represents the position of the SMBH particle, and the jet injection axis is the vertical axis of the image.

In the text

	Fig. 17. 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 GGM, revealing shocks, ripples, and cavities.
In the text

	Fig. A.1. Total amount of cold gas contained in the inner 40 kpc of the simulated box.
In the text

	Fig. A.2. Evolution of r80, the radius containing 80% of the cold gas in our various runs.
In the text

	Fig. A.3. Gas density projection in the inner 30 kpc of our simulated cluster for each of our simulations.
In the text

Fig. B.1. Comparison of the clump distribution in the inner 10 kpc of our fiducial MHD simulation at t = 960 Myr with three different method. The blue and orange curves shows the distribution of the clumps with yt (Turk et al. 2011) based on density and temperature criteria. The red curve indicates the clump distribution using Scipy’s ndimage (Virtanen et al. 2020), and is comparable to the method used in (Gronke et al. 2022) and Das & Gronke (2024). The inclination of a power law of slope α ∼ −0.6 is indicated for comparison.

In the text