© The Authors 2025

1. Introduction

Around 10% of the total mass of galaxy clusters is thought to be contained in a hot plasma of ionized gas, usually designated as the intracluster medium (ICM; e.g., Chiu et al. 2015). Cooling via Bremsstrahlung gives rise to the emission of radiations and makes clusters strong X-ray sources (Felten et al. 1966; Mitchell et al. 1976; Serlemitsos et al. 1977; Sarazin 1986). Mass estimates of clusters, as well as density or temperature profiles, can be derived from observations of this emission, assuming that the ICM is in hydrostatic equilibrium. However, both simulations and observations have shown that several processes are likely to cause deviations to the hydrostatic equilibrium (Rasia et al. 2006; Lau et al. 2013; Biffi et al. 2016; Barnes et al. 2021). In particular, it is expected that non-thermal sources of pressure, such as turbulence, might add to the thermal pressure support of the ICM against gravity (Fusco-Femiano & Lapi 2017; Angelinelli et al. 2020; Vazza et al. 2011).

Cluster mergers are thought to be the main drivers of turbulence on the scale of entire clusters. In a typical merger event, two clusters collide with relative velocities of ∼10³ km s⁻¹. As the ICM relaxes in the post-merger phase, shocks form and propagate through the ICM, carrying an energy of up to 10⁶⁵ erg and providing heat to the plasma on scales of up to tens of Mpc (e.g., Markevitch & Vikhlinin 2007). The subsequent dissipation of turbulence heats the ICM on longer timescales. Results from hydrodynamical simulations suggest that turbulent kinetic energy in perturbed clusters resulting from mergers can account for up to 20–30% of the ICM’s thermal energy (Dolag et al. 2005; Vazza et al. 2011). A common consequence of mergers in galaxy clusters is a temporary offset of the gas core with respects to its dark matter counterparts. In such cases, the ICM falls back toward the center of the cluster, and the resulting gas sloshing can drive turbulence on scales of up to 100 kpc (ZuHone et al. 2013). In the innermost regions of galaxy clusters, active galactic nuclei (AGNs) are suspected to be the main drivers of turbulence (Scannapieco & Brüggen 2008; Li et al. 2020). By accelerating jets of gas up to relativistic speeds, AGNs inflate cavities in the ICM, which rise buoyantly and propagate out to distances of a few to hundreds of kiloparsecs (Bîrzan et al. 2004; Dunn et al. 2005; Prunier et al. 2025). As they rise, these cavities displace the ICM, mix with the surrounding gas, and inject turbulent energy (Brüggen et al. 2009; Gatuzz et al. 2024). Constraining turbulence with observations is challenging and typically relies on indirect measurements. X-ray surface brightness (SB) fluctuations can be computed by subtracting smoothed models from the SB data (Churazov et al. 2012; Zhuravleva et al. 2014). The velocity power spectrum of the gas can then be inferred from the power spectrum of the SB fluctuations, but it is subject to a range of assumptions. X-ray calorimeters can also provide measurements of the typical turbulent velocities in cluster cores through line broadening measurements (Hitomi Collaboration 2016; XRISM Collaboration 2025). More recently, Li et al. (2020) have proposed using observations of ionized filaments found in the inner tens of kiloparsecs of cool-core clusters to probe the turbulent motion of the ICM. Such filamentary structures are common in cool-core clusters (Olivares et al. 2019) and are thought to originate from thermal instabilities in the ICM (McCourt et al. 2012), likely seeded by AGN outflows (Li & Bryan 2014; Fournier et al. 2024). The line-of-sight velocity of these filaments can be obtained from optical spectrometers with angular resolutions of less than ∼10⁻¹ arcsec². By measuring the velocity structure function (VSF) of the filaments, Li et al. (2020) found deviations from a Kolmogorov scaling characterized by a steeper power-law index than the expected ℓ^1/3, where ℓ designates the spatial scale. The origin of this steepening is debated in the literature. One hypothesis is that it is related to the injection of turbulence caused by the central AGN. As the bubbles rise in the ICM, they could drive turbulence and increase the amplitude of the VSF at scales on the order of their typical size (i.e., ∼10 kpc; Hillel & Soker 2020). Other possibilities include the contribution of gravitational acceleration (Wang et al. 2021) or the condensation of cold gas within supersonic outflows (Hu et al. 2022), which have both been found to result in VSFs steeper than the Kolmogorov scaling law. The interpretation of these observations could be affected by projection effects or limited angular resolution. Some of these effects have been investigated in turbulent box simulations (Mohapatra et al. 2022) and idealized cluster setups (ZuHone et al. 2016; Wang et al. 2021; Sotira et al. 2025).

In this paper, we use a high-resolution exascale magneto-hydrodynamical (MHD) simulation from the XMAGNET project¹ to characterize the VSFs of the multi-phase ICM. We explore the effects of projection, emission weighting, atmospheric seeing, and viewing angle on the VSF. In Sect. 2, we describe our numerical method and our simulation parameters. In Sect. 3, we quantify the turbulent motions of both the hot and cold phases using a first-order VSF. In Sect. 4, we investigate up to what extent projection effects might affect the measurement of the VSF before discussing our results in Sect. 5 and concluding in Sect. 6.

2. Method

2.1. Simulation suite and initial conditions

The XMAGNET simulation suite used in this study is introduced and described in Grete et al. (2025). Hereafter, we summarize the key points of the methodology presented in that paper.

We ran hydrodynamical and MHD simulations of an idealized cool-core cluster on the Frontier supercomputer². We used ATHENAPK³, a performance-portable version of the ATHENA++ astrophysical magnetohydrodynamics code (Stone et al. 2020). It is based on the block-structured refinement framework PARTHENON (Grete et al. 2023a) and includes magnetic fields, radiative cooling, and AGN and stellar feedback. The simulations used an overall second-order accurate, shock-capturing, finite-volume scheme consisting of second-order Runge Kutta (RK2) time integration, piecewise-linear reconstruction, and an HLLD Riemann solver (Miyoshi & Kusano 2005).

Gravity is treated as an acceleration term parametrized by a superposition of a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997), a Hernquist potential representing the brightest central galaxy’s mass (Hernquist 1990), and a central point mass modeling the effect of 3C 84, the supermassive black hole (SMBH) in the Perseus cluster’s BCG. The initial gas distribution is calculated numerically to obey hydrostatic equilibrium assuming the entropy profile of the Perseus cluster from the ACCEPT catalog (Cavagnolo et al. 2009). Optically thin cooling is calculated by the exact integration method introduced by Townsend (2009) and uses solar metallicity cooling tables by Schure et al. (2009) down to 10^4.2 K. A temperature floor of 10⁴ K was also enforced. Purely rotational velocity and magnetic perturbations are seeded at t = 0 to break the symmetry of the system, with respective dispersion of δv = 75 km s⁻¹ and δB = 1 μG. These perturbations are generated using an inverse parabolic spectrum with a peak scale of 100 kpc. AGN feedback is modeled via source terms for thermal, kinetic and magnetic energy. Cold material (defined as gas cells of temperature T ≤ 5 × 10⁴ K) is accreted at a rate Ṁ_acc related to the accretion power by ${\dot{E}}_{\mathrm{acc}}=\eta {\dot{M}}_{\mathrm{acc}}{c}^2$, where η = 0.001 is the accretion efficiency and c the speed of light. This accretion power is then split into the three feedback channels by their corresponding feedback fractions f_i, where i designates the kinetic, thermal and magnetic feedback channels. In the absence of star particles, stellar feedback is calculated analytically, by extracting gas denser than a density threshold $n_{\mathrm{thresh}.,\star }=50{\mathrm{cm}}^{-3}$ and redistributing it into thermal energy with an efficiency of 5 × 10⁻⁶. The structure of the grid is identical in all our runs and kept constant across time. The root mesh grid is made of 1024³ cells covering a total volume of (6.4 Mpc)³. The root grid is refined a further ℓ_max = 6 cubic and nested refinement levels. As a result, the maximally refined region has a volume of (250 kpc)³ and is covered by 2560³ cells of size Δ_x = 97.7 pc. The entire analysis presented in this paper is based on data extracted from this maximally refined region, i.e. for gas cells located at radii r ≤ 125 kpc. We use periodic boundary conditions. Considering the total width of the simulated volume, we can exclude the possibility that re-entering sound or shock waves impact our analysis, entirely performed in the inner (200 kpc)³ of the box.

In this study we analyze the data from the fiducial simulation of the XMAGNET suite. This run comprises magnetic fields as well as solar metallicity cooling. The simulation is run for a total duration of 4 Gyr. A detailed analysis of the cluster’s evolution in this simulation is presented in Grete et al. (2025). We note that the feedback power is rather constant for times t ≥ 500 Myr, with no clear duty cycle. The total cold gas mass varies between 10⁹ and 10¹¹ M_⊙. A visualization of the inner 80 kpc at t = 3.2 Gyr is presented in Fig. 1. It shows a projected map of the cold gas line-of-sight velocity, weighted by Hα emission, overlaid on top of a projected temperature map of the hot ICM. A zoom on the innermost region shows a comparison between the original projection, and a smoothed version produced to mimic the effects of atmospheric seeing.

Fig. 1. Hα emission – weighted projection of the line-of-sight velocity for the inner 80 kpc of our simulated box at t = 3.2 Gyr. A zoom-in on the innermost region showcases the effect of atmospheric seeing on the map. The background image is a projection of the temperature field. The figure can be compared to Fig. 2 from Li et al. (2020). ⊙ and ⊗ markers indicate the directions pointing toward (away) from the viewer, associated with negative (positive) line-of-sight velocity (respectively). More detail on the emissivity weighting are provided in Sect. 4.1. A movie showing the three dimensional structure of the cold phase is available online.

2.2. Velocity structure functions

A useful tool to quantify the scaling of turbulence is the pth-order VSF. It is a two-point correlation function defined as

$$\begin{array}{c}\hfill S_p(\ell )={⟨|\mathbf{u}(\mathbf{x}+\ell ·\mathbf{e})-\mathbf{u}(\mathbf{x})|}^p⟩,\end{array}$$(1)

where u is the velocity field evaluated at a pair of points of positions x and x + ℓ ⋅ e, and e is a randomly oriented unit vector. It is expected that the VSF scales as S_p(ℓ)∝ℓ^p/3 in the inertial range, assuming that the turbulence is isotropic, homogeneous and that the fluid is incompressible (Kolmogorov 1941). To facilitate direct comparison with observations, we focused on the first-order VSFs (i.e., S₁(ℓ)). It is worth noting that expressing the second-order VSF S₂(ℓ) in spectral space leads to the regular Kolmogorov scaling law:

$$\begin{array}{c}\hfill E(k)\propto {\epsilon }^{2/3}{k}^{-5/3},\end{array}$$(2)

where E(k) is the turbulent energy spectrum, ε is the turbulent dissipation rate, and k is the wavenumber. In Appendix C, we also present second-order VSFs for comparison.

We compute the VSF by randomly picking pairs of points separated by a distance ℓ, calculate their associated velocity difference, and then average them into bins of logarithmically increasing size. Since we compute the VSFs in a region containing a relatively large number of cells (i.e., 1024³), we employ a parallelized algorithm adapted from Federrath et al. (2010). Sampling pairs of points with separations on the order of a few cell sizes can be inefficient. When randomly selecting two cells in the grid, the probability of choosing a pair with separation ℓ is proportional to ℓ². Moreover, since the grid is discretized, the set of possible separations is not continuous. For separations smaller than 500 pc, we adjust the bin sizes to account for the discrete nature of the grid. Specifically, the separation between two cells is given by

$$\begin{array}{c}\hfill \ell =\sqrt{i^2+j^2+k^2}{\mathrm{\Delta }}_x,\end{array}$$(3)

where Δ_x is the cell size and i, j, k are integers representing the displacement between the two cells along the three Cartesian directions of the grid. When the grid is two-dimensional (for instance when computing projected VSFs), Eq. (3) reduces to $\ell =\sqrt{i^2+j^2}{\mathrm{\Delta }}_x$. To increase the efficiency of the point sampling, we randomly sampled N points on the grid, representing the first pixel of each pair. We then randomly sample a neighboring cell for each of these points following a uniform probability distribution of distance, up to 500 pc. The VSF above 500 pc is computed using a naive sampling of points (i.e., by uniformly sampling both pixels of each pair). The value of 500 pc has been chosen to maximize the efficiency of the code. We have verified that our results are equivalent to the VSF obtained when entirely computed from naive pair sampling. For projected VSF of the cold phase (see Sect. 4), the limited number of pixels allowed us to explicitly calculate the velocity difference of all possible pairs of pixels. Convergence tests of our VSF algorithms are presented in Appendix B.

3. Quantifying turbulence in the intracluster medium

3.1. Velocity dispersion

The time-averaged velocity probability density function (PDF) of the hot and cold gas is presented in Fig. 2. As the kinematics of the cold phase at the beginning of the active period of the simulation might carry some imprints from the initial conditions, we perform the averaging over all snapshots with time t ≥ 0.9 Gyr, i.e. 400 Myr after the first cold clumps appears. Darker to brighter colored curves represent annuli with increasing distances from the center, namely r ≤ 25 kpc, 25 ≤ r ≤ 50 kpc and 50 ≤ r ≤ 100 kpc. The dashed gray curves indicate the data for the gas located at radii r ≤ 100 kpc. The standard deviation σ(v_x) for each curve is derived by fitting a Gaussian function to each distribution. The colored areas represent the standard deviation resulting from time variation. The colored dashed lines indicate the median velocity. The standard deviation of the velocity of the hot gas decreases with increasing radius by roughly a factor of two between the inner and the outer regions of the maximally refined region of our simulated box, with a value of ∼60 km s⁻¹ in the innermost region. This value is relatively low compared to measurements of X-ray calorimeters (Hitomi Collaboration 2016) of order ∼ 150 km s⁻¹, but consistent with other idealized setups similar to ours (Wang et al. 2021; Ehlert et al. 2021). This could indicate that the contribution of other sources of turbulence, such as mergers, sloshing, or accretion cannot be ignored to reproduce the measured velocity dispersions. We emphasize that our upper temperature cutoff of 10⁸ K excludes the jet material while still including the (shocked) ICM. The velocity distribution of the cold phase is much wider, with dispersions on the order of ∼150 km s⁻¹ in the innermost region. The turbulence within the cold phase is likely supersonic, as the typical velocity fluctuations within the filaments are on the order of 10² km s⁻¹, roughly an order of magnitude above the sound speed. We have verified that including emission-weighting does not significantly impact the results for either the hot or the cold phase.

Fig. 2. Probability density function (PDF) of the velocities’ x-component for the hot phase (106 ≤ T ≤ 108 K, left panel) and cold phase (T ≤ 105 K, right panel) for various radius bins and the PDF of all gas located at radii r ≤ 100 kpc, averaged between 0.9 and 4 Gyr. Vertical dashed lines indicate the median values. Colored areas represent the standard deviation resulting from time variation. We have verified that the y-component leads to similar results.

3.2. Radial velocity

Similarly, we present the distribution of the radial velocity of the hot and cold gas in Fig. 3. The hot phase (left panel) distribution is characterized by a negative mean radial velocity, as large as ∼ − 50 km s⁻¹ in the innermost region and ∼ − 20 km s⁻¹ when including the whole maximally refined region. The velocity distribution of the cold phase is more complex and radially dependent. In the outermost region, it follows a broadly Gaussian distribution with a mean radial velocity of ∼ − 400  km s⁻¹. In the innermost region, the distribution is broader and more asymmetric. At the lower end of the distribution the velocity goes down to ∼ − 10³ km s⁻¹. A higher tail of positive radial velocity is also visible and is likely related to uplifting of warm and cold gas resulting from the coupling with the outflowing jets (see also Fournier et al. 2024). Overall, the motions of the hot and cold phases are mostly decoupled. The cold phase is dominated by a free-falling motion, while the hot phase only slightly deviates from hydrostatic equilibrium. The resulting turbulent to thermal pressure ratio for the hot phase is typically of ∼1%.

Fig. 3. Same as Fig. 2 but for radial velocity, vr.

3.3. Velocity structure functions

3.3.1. Overall intracluster medium velocity structure functions

In Fig. 4, we present the three-dimensional VSF (i.e., obtained taking the three-dimensional distance between pairs of cells) of the hot (left panel) and cold (right panel) phases. To ease comparison with observations, we restrict our analysis to spheres of varying radius, centered on the AGN. While a more thorough analysis of individual cold filaments’ kinematics might provide deeper insights into their dynamics and evolution, we postpone such analysis to future work. The solid (dashed) lines indicate the VSFs obtained from the magnitude of the three-dimensional velocity vectors (x-component of the velocity vectors) of each pairs of cells, respectively. Shaded areas indicate time variation throughout the whole active part of the simulation (i.e., starting from the formation of the first cold structures at t ∼ 0.9 Gyr). The different colors indicate VSFs obtained including gas located from various radial distances from the center of the box. At low spatial separations (typically up to ∼1 kpc), the power-law index of the VSFs is likely dominated by effects of numerical viscosity, leading to artificially increasing power-law indices (see also Grete et al. 2023b, for a more extensive discussion on this effect). At larger separations, both the time-averaged three-dimensional and line-of-sight VSFs are steeper than the Kolmogorov scaling. There is no significant variations of this scaling when including gas at larger radii. The cold phase VSFs follow a power law with a constant index between ∼500 pc and ∼10 kpc. The time-averaged VSFs are typically steeper than the Kolmogorov scaling index of m = 1/3, and on the order of m ∼ 0.4. In both phases, only including the x-component of the velocity vector decreases the amplitude of the VSFs by a factor of ∼2 − 3, which is larger than the $\sqrt{3}$ ratio expected in the case of isotropic turbulence. Several previous studies have discussed the possibility of using Hα-emitting gas or surface brightness fluctuations of the X–ray luminosity as a proxy to trace the hot gas VSF. In Fig. 5, we present the relationship between the power-law index of the line-of-sight (m_hot, x) and unprojected (m_hot, 3D) hot phase VSFs (top panels) and between that of the line-of-sight, cold phase (m_cold, x) and unprojected hot phase VSFs (m_hot, 3D) (middle panels) for each snapshot. The lower panels present the relationship between the amplitude v₀ of the line-of-sight, cold phase VSFs and the unprojected hot phase VSF. The left, center and right panels indicate results for VSFs obtained from all gas cells contained within a limit radius of 12, 25 and 50 kpc, respectively. Each power-law index m and amplitude v₀ were measured by fitting the VSFs using

$$\begin{array}{c}\hfill {log}_{10}{S}_1(\ell )={log}_{10}(v_0)+m{log}_{10}(\ell ),\end{array}$$(4)

where ℓ is the separation scale and v₀ is defined as the VSF amplitude. The fit is performed between 2 and 4 kpc, avoiding the region influenced by numerical viscosity and staying below the window effect caused by the limited domain size. The gray dashed line indicates the identity relationship.

While the power-law index of the line-of-sight, hot gas VSF shows a good correlation with its three-dimensional counterpart, there is no visible correlation between the power-law index of the line-of-sight, cold gas VSFs and that of the three-dimensional hot gas VSFs. The power-law index of the hot phase VSF is relatively constant throughout the whole duration of the simulation, and is mostly centered around m ∼ 0.4, with no significant variation when including gas at larger radii (see middle and right columns). The power-law index of the cold phase shows more fluctuations, ranging from m ∼ 0.2 to m ∼ 0.5. The power-law index of the hot phase VSFs are preferentially steeper than their unprojected hot phase counterparts. Similarly, we find no clear correlation between the amplitude of the cold and hot phases VSFs. We emphasize that the cold phase VSFs can vary significantly from one snapshot to an other and thus that many internal details of individual VSFs are smoothed out by time-averaging. Cold phase VSFs at individual snapshots are presented later on (see Sect. 4.6).

3.3.2. Effect of rising cavities

In Li et al. (2020) and Ganguly et al. (2023), it has been suggested that features observed in the line-of-sight Hα emitting gas VSF could be related to the presence and size of the rising cavities related to the AGN activity. In the previous subsection, we have reported the absence of strong correlation between the parameters of the VSFs for the gas located in the inner tens to hundred of kpc of the simulated box. To evaluate the possibility that the AGN activity still might affect the VSF locally, we re-run our cluster setup simulation and save a larger number of snapshots, with an output saving interval of Δt_snap = 2 Myr, starting from t = 1.5 Gyr.

Fig. 4. Velocity structure functions of the hot ICM (left panel) and the cold phase (right panel) for all gas located less than 12, 25, and 50 kpc away from the center of the simulated box. The dashed colored lines indicate the line-of-sight VSFs taking the x-component of the velocity field. We emphasize that the latter are not equivalent to projected VSFs, as the distance between pairs of cells is still three-dimensional (projected VSFs are discussed in Sect. 4). Colored area is showing standard deviation over time. The solid gray lines in the background indicate the Kolmogorov scaling law (Kolmogorov 1941). As our setup closely resembles the fiducial MHD run from Wang et al. (2021), we present the time-averaged line-of-sight VSFs for the hot and cold phases from that work, for radii r ≲ 10 kpc (dashed gray line labeled W21). A comparison of our results with this study is discussed in Sect. 5.2. We emphasize that VSFs below separations of ℓ ∼ 800 pc are likely to be steepened by the effect of numerical viscosity.

Fig. 5. Top panels: Relationship between the power-law index of the line-of-sight hot gas VSF and that of its three-dimensional counterpart. Middle (bottom) panels: Relationship between the power-law index m (amplitude v0) of the line-of-sight, cold phase VSF, and the three-dimensional hot phase VSFm, respectively. The color of the points represents time, with values ranging from 0.9 to 4.0 Gyr. The columns show VSF values for all gas within radii of 12, 25, and 50 kpc (left to right). The dashed gray line indicates the identity relationship and should not be interpreted as an attempt to fit the data. The light blue shaded area represents the region of the parameter space where the measured power-law index of the VSF is steeper than the m = 1/3 power-law index predicted by the Kolmogorov theory.

In Fig. 6 we present the evolution of the three-dimensional hot phase VSF evaluated in a sphere of radius r_sphere = 10 kpc between t₀ = 1.4 Gyr and t_f = t₀ + 100 Myr. The sphere was positioned 12 kpc below the SMBH along the z-axis, thus allowing it to encompass the whole volume of a cavity inflated by AGN activity at t − t₀ ∼ 65 Myr. The upper panels of Fig. 6 show projections of the temperature field with a thickness of 1 kpc. The second row shows the corresponding VSFs of the hot (red curves) and cold (gray curves) phases, computed within that sphere. The thin solid gray lines indicate the Kolmogorov scaling law of index m = 1/3, while the dashed orange line indicates the fit of the hot phase VSFs. In the lower panel, we present the time evolution of the corresponding power-law index for the cold and hot phases. Before the cavity starts forming, the power-law index of the VSF is around m = 0.3. As the jets inflate the cavity, the power-law index of the VSF becomes steeper and increases by a factor of ∼2 to reach ∼0.5 − 0.6. At later times, the cavity rises through the ICM, loses its spherical shape, and exits the sphere of cells under consideration. The power-law index of the VSF decreases again and converges to m = 0.4. The VSFs of the cold phase show wider variations and no definitive conclusion can reasonably be drawn from the time evolution of its power-law index. A peak in the VSF at scales of ∼10 kpc is visible during the phase when the cavity structure is clearly visible.

Fig. 6. Top panels: temperature projection of the gas. The dashed circle indicates the sphere within which the VSF are computed, whose position is kept fixed throughout the run. The round marker at the top of the temperature projection indicates the position of the SMBH particle and the arrow shows the jet injection axis. The projection depth is 1 kpc. Middle panels: Velocity structure functions of the cold (gray curves) and hot (red curves) phases. The dashed orange lines indicate the fit of the power-law index of the hot phase VSF. Lower panel: Time evolution of the cold and hot phases power-law index. Here, t0 refers to the restarting time of the simulation.

4. Effects of projection and comparison to observations

4.1. Mock Hα observations

In Li et al. (2020) and Ganguly et al. (2023), the optical imaging Fourier transform spectrometer (SITELLE) and the Multi Unit Spectroscopic Explorer (MUSE) are used to create maps of the optical emission from the warm ionized phase of several BCGs. The resolved detail in such maps depends on angular resolution and atmospheric seeing, where the latter is standardly modeled by a Gaussian smoothing kernel. For the Perseus cluster, SITELLE’s angular resolution corresponds to δx_pxl ∼ 250 pc, and atmospheric seeing induces a FWHM smoothing of ∼0.42 kpc (Li et al. 2020). To evaluate the effects of spatial resolution and seeing, we produce mock Hα emission-weighted velocity maps. We filtered the gas to keep only the cold component of temperature T ≤ 10⁵ K. The emissivity j_Hα of the gas in each cell is given by

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

where T₄ = T/10⁴ K and T, n_e, and n(H⁺) are the temperature, electron density, and ionized hydrogen density of each gas cell, respectively (Dong & Draine 2011). The line-of-sight two-dimensional velocity maps were then obtained by projecting the x-component of the velocity field on the (y, z) plane. The VSF of the projected map is then obtained by sampling all available pairs of pixels. The separation between each pair of pixels is computed from their projected, two-dimensional distance. We also compute the VSFs obtained from smoothed projected maps to mimic the effect of atmospheric seeing. These maps are obtained by applying a Gaussian kernel to the original Hα map (see Fig. 1). In Sect 4.3.1, we compute VSFs out of projected velocity maps with degraded resolution. The method employed to generate these resolution-downgraded maps is described in Appendix E.

4.2. Mock X–ray observations

We also evaluated the impact of projection on the hot phase. To this end, we used the same method as for the optical emission. We separated the hot phase by only projecting gas whose temperature is above 10⁶ K and below 10⁸ K (to exclude the jets) and applied either no weighting or a square density weighting. A visualization of such emission-weighted projection is presented in Fig. 7. Although the cold phase is filtered out from these projections, the continuity of density between the hot ICM and the denser, cold filaments results in over densities in the hot phase enveloping the cold clumps. The square density weighting of the projected map enhance the contribution of these over densities, which are visible in the projected velocity maps. Their impact on the projected hot VSFs is discussed in Sect. 4.4.

Fig. 7. Emission-weighted (w) line-of-sight velocity projection of the hot phase within the inner 50 kpc of our fiducial MHD run at t = 3.1 Gyr. The marker and dashed line at the center indicate the position of the SMBH particle and the axis of the jets, respectively.

4.3. Effects of projection and smoothing

4.3.1. Cold phase

In Fig. 8, we show the effects of various projection methods on the resulting VSFs. All curves corresponds to the time-average of all VSFs between 0.9 and 4 Gyr for our fiducial MHD run and for all cold gas located at r ≤ 25 kpc. The colored areas represent the standard deviation resulting from time variation across all snapshots. The darkest curve indicates the VSF of the x-component of the velocity field, computed by taking the three-dimensional distance between each pair of cells. All the other curves are computed based on projected maps, implying that separations between pairs of pixels are computed taking their two-dimensional distances. The third curve is computed taking projected maps of the velocity field weighted by Hα emission. The last and brightest curve includes smoothing resulting from atmospheric seeing, assuming a smoothing kernel with a standard deviation consistent with the values expected for the Perseus cluster (i.e., a FWHM of 0.42 kpc, see Li et al. 2020). The blue vertical line indicates the full width at half maximum (FWHM) of the Gaussian kernel used to smooth the projected velocity maps. The effect of smoothing on the VSFs power-law index taking other filters’ kernel sizes is discussed in Appendix D.

Fig. 8. Line-of-sight VSFs for the cold component and for the projected radius r ≤ 25 kpc averaged between 0.9 and 4 Gyr. The darkest curve (labeled ℓ3D) was obtained by taking the three-dimensional distance between pairs of points, similarly to Wang et al. (2021). For all the other curves, the VSFs were computed out of projected maps, and the corresponding separations scales are two dimensional (ℓ2D). The term w(Hα) indicates weighting of the map by Hα emission, and 𝒮 is the additional smoothing effect resulting from atmospheric seeing (see Li et al. 2020). The shaded areas indicate the time variation of each VSF.

Taking into account the three-dimensional distribution of the cold gas, the time-averaged VSF is slightly steeper than the Kolmogorov scaling law, with a mean power-law index of m ∼ 0.4. When calculating the VSF extracted from projected maps, the time-averaged VSF flattens and its power-law index reduces to m ∼ 0.2. Including weighting by Hα emission has little impact on the power-law index, although leading to a minor additional flattening with respect to the unweighted projected VSF. This is consistent with results of previous studies (see, e.g., Mohapatra et al. 2022), and is likely related to the low volume filling fraction of the cold phase which results in limited cancellations along the line of sight. Finally, when taking into account the effect of atmospheric seeing on our maps by applying a Gaussian smoothing kernel, the power-law index of the VSF is dramatically increased. Below scales on the order of the smoothing kernel size, the power-law index is close to m = 1. Above this typical scale, the effect of smoothing on the power-law index becomes less significant, although it still affect the overall amplitude VSF. Since we use a smoothing kernel assuming the distance to the Perseus cluster, the expected effect of smoothing on the line-of-sight velocity maps of galaxy clusters located further away could be more important.

To evaluate the effect of smoothing on individual VSFs for scales on the order of the smoothing kernel size, we present in Fig. 9 the relationship between the power-law index $m_{\mathrm{cold}}^{\mathrm{us}}$ of the unsmoothed (indicated here and in the figure with “us”) projected VSFs (horizontal axis) and Δm, which represents the increase in the power-law index resulting from smoothing. All power-law indices are obtained by fitting the VSFs for separations of ℓ = 0.4 ± 0.1 kpc. As a comparison, lines of equal smoothed VSFs power-law index $m_{\mathrm{cold}}^{\mathrm{s}}$ are represented with the solid gray lines. Intrinsically shallower VSFs are more affected by smoothing than the steeper ones, as Δm varies broadly linearly with $m_{\mathrm{cold}}^{\mathrm{us}}$, with a slope close to –1. Consequently, the power-law index of the smoothed VSFs remains relatively constant across snapshots, with $m_{\mathrm{cold}}^{\mathrm{s}}=0.86\pm 0.06$.

Fig. 9. Relationship between the power-law index $m_{\mathrm{cold}}^{\mathrm{us}}$ of the unsmoothed (“us”) projected VSFs (horizontal axis) and Δm, which represents the increase in the power-law index resulting from smoothing. Lines of equal power-law index for smoothed VSFs ($m_{\mathrm{cold}}^{\mathrm{s}}$) are also indicated for comparison (gray lines).

Instruments used to obtain line-of-sight velocity maps of galaxy clusters have limited angular resolution, translating into limited spatial resolutions depending on the distance of the observed galaxy. To evaluate possible consequences on the VSFs, we present in Fig. 10 the influence of varying spatial resolution on the velocity maps of the cold phase. While the darkest curve is computed using projected maps with pixels of width Δ_x, 0 = 97.7 pc (corresponding to the maximum resolution of our simulation setup), all brighter curves correspond to VSF obtained from maps with degraded resolution, resulting in pixels of width Δ_x = k Δ_x, 0, where k is an integer. Downgrading the resolution leads to a steepening of the VSFs. The typical VSF amplitude decreases for separations scales lower than the domain size (i.e., for ℓ ∼ 25 kpc).

Fig. 10. Velocity structure functions obtained from projected maps taking various spatial resolutions for the cold gas located at a projected radius of r ≤ 25 kpc. The projection method for all curves is equivalent to the third curve in Fig. 8 (i.e., labeled as ℓyz + w(Hα)). Consequently, the VSFs presented here are shallower than their three-dimensional counterparts (labeled as ℓ3D in Fig. 8) due to projection and emission weighting. The darkest curve corresponds to the maximum spatial resolution Δx, 0 of our simulation. All brighter curves are obtained from coarse-grained maps with pixel size Δx. The colored areas represent the standard deviation resulting from time variation across all snapshots.

4.4. Hot phase

In Fig. 11, we study the effects of projection on the VSF of the hot phase in a similar way as we did in Fig. 8, but for the cold phase. w indicates that the VSFs were obtained from emission – weighted velocity maps. The amplitudes of the projected VSFs at any scale are typically lower by a factor of ∼2 compared to the three-dimensional VSFs. Both projected VSFs (weighted and unweighted) are steeper than the three-dimensional VSF for separations above ℓ ∼ 1 kpc, and have power-law index of m ∼ 1/2. Including emission weighting makes the VSF shallower than the unweighted ones for lower separations. This is likely resulting from the presence of the over dense hot clumps enveloping the cold filaments (see Fig. 7). We also note that the steepening of the weighted and projected VSF depends on the choice for the lower temperature cutoff used for the projection. Lower temperature cutoff results in stronger over densities in the hot phase, projected velocity maps with higher contrast, and consequently flatter VSFs. We emphasize that varying the temperature cutoff between 10⁶ and 10⁷ K has no impact on the unprojected or unweighted VSFs. These results are broadly consistent with that reported from turbulent box simulations (see the projected and unprojected VSFs of the radiative MHD run f0.10magHR in Fig. 7 from Mohapatra et al. 2022, for a direct comparison).

Fig. 11. Same as Fig. 8 but for the hot phase. The brightest curves indicate the VSF obtained from weighted and projected velocity maps, taking three different minimum temperature cutoffs. The projected VSFs are steeper than the unprojected one for separations above ℓ ∼ 1 kpc.

4.5. Effects of the observer orientation

In Fig. 12, we present a set of 2000 VSFs computed from Hα emission-weighted maps of the cold gas line-of-sight velocity in the inner 50 kpc of our simulated cluster at t = 3.2 Gyr. Each map is obtained by projecting the line-of-sight velocity of each cold gas cell along a randomly selected orientation identified by its azimuthal and polar angles (ϕ,  θ). As discussed earlier, the VSF do not follow a power law with a fixed index within the whole interval of separations considered here, as numerical viscosity results in artificially steeper VSF for spatial separations of a few cells. Here, we thus fit the VSF between 2 and 4 kpc, above the region affected by numerical viscosity and below the window effect resulting from the limited domain size. Varying these boundaries does not result in major changes as the power law typically keeps a constant index between 2 and 10 kpc. The values for v₀ and m are then presented following an Aitoff projection. Since we do not include effects of opacity during the projection of the velocity map, and the left and right half of the Aitoff projection are symmetrical for both v₀ and m. The injection axis of the jets corresponds to θ = ±90°–that is, the vertical axis at ϕ = 0°.

Fig. 12. Aitoff projection of the values of the VSF amplitude v0 and the power-law index m, measured for 2000 randomly selected orientations at fixed time t = 3.2 Gyr within r = 25 kpc. Each VSF was computed from projected Hα emission weighted maps. The power-law fitting was performed between 2 and 4 kpc. As we do not include effects of opacity when generating the projection of the velocity maps used to obtain each VSF, the left and right half of the plot are symmetrical for both v0 and m. The coordinates of the jet injection axis is at θ = ±90°.

While the typical VSF amplitudes does not vary strongly along the azimuthal angle ϕ, it varies of nearly a factor of two between viewing angles parallel to the midplane of the simulated box (θ = 0°) and those parallel to the jet axis (θ = ±90°). This is likely related to the orientation of the filaments, which tend to align with the AGN jet axis. Since they are preferentially infalling structures, their line-of-sight velocity amplitude is larger along the z axis than the perpendicular directions with respect to the jets. As roughly half of the gas is falling along the +z direction and the other half along the −z direction, this bipolarity results in larger velocity differences at separations on the order ℓ ≤ 10 kpc. The power-law index of the VSF also varies by more than a factor of 2 depending on the observer’s orientation. Its value is typically found between m = 0.15 and m = 0.35, and is on average shallower than the Kolmogorov value of 1/3. We attribute this flattening to projection (see Sect. 4.3). The polar and azimuthal dependence of the power-law index is likely influenced by the three dimensional orientation of the filaments, and is likely to vary from one snapshot to an other.

4.6. Comparison with optical observations

In order to compare our results to those in Li et al. (2020), we show in Fig. 13 time-averaged VSFs of the cold gas obtained from Hα emission-weighted maps. The maps have been coarse-grained to a spatial resolution of 291 pc (the actual resolution of the SITELLE instrument is approximately 255 pc, but the coarse-graining method used here only permits the resolution to be reduced by integer multiples of the original map’s spatial resolution, namely 97 pc). The labels “inner”, “outer”, and “all” (see legend) designate the VSFs obtained from all gas cells located at projected radii less than 12 kpc, at projected radii more than 12 kpc, and for all gas cells, respectively. Gray curves correspond to the VSF inferred from optical observations and are taken from Li et al. (2020). The left and right panels show results for VSFs calculated from unsmoothed and smoothed, projected velocity maps, respectively. The unsmoothed, coarse-grained map is typically shallower than the Kolmogorov scaling law m = 1/3. The amplitude at low spatial separations is also greater than observations by a factor of two. Smoothing results in time-averaged VSFs that have power-law index and amplitude more consistent with observations. We note that numerical viscosity can impact the dynamics of gas for separations lower than ℓ∼ 1 kpc (see Sect. 5).

Fig. 13. Time-averaged VSFs of the cold phase obtained from mock Hα emission-weighted maps projected along the x-axis (red lines) compared to observational data (gray lines, from Li et al. 2020). The images have been coarse-grained to a spatial resolution of 290 pc to mimic the spatial resolution of instruments assuming Perseus cluster distance. The right panel was obtained from smoothed maps to mimic the effects of atmospheric seeing. The labels “inner”, “outer”, and “all” (see legend) respectively designate the VSFs obtained from all gas cells located at a projected radii of less than 12 kpc, at a projected radius of more than 12 kpc, and for all gas cells irrespective of their radius.

Time-averaging VSFs over nearly 4 Gyr (i.e., approximately 170 snapshots) naturally smoothes out many details visible in individual VSFs calculated at fixed times. Although the steepening effect due to smoothing is visible in all VSFs, the relative curvature of the time-averaged smoothed VSFs presented in the right panel of Fig. 13 does not accurately represent the typical shape of individual VSFs. To provide context for the time-averaged VSFs, we present in Fig. 14 individual projected and smoothed VSFs calculated at three different times, namely: 1.28, 3.28, and 3.48 Gyr. For comparison, we also include the 11 VSFs of nearby BCGs from Ganguly et al. (2023). Our individual VSFs were selected to highlight typical trends observed throughout our simulation, namely: (i) VSFs characterized by a low amplitude and a relative flatness for separation scales above a few kiloparsecs (t = 1.28 Gyr), (ii) VSFs exhibiting a distinct peak structure at separation scales beyond 10 kpc (t = 3.28 Gyr), and (iii) VSFs closely matching the observed VSFs both in terms of steepening and amplitude (t = 3.48 Gyr). Although a visual inspection of the hot and cold phases suggests that peaks at large separations correlates with the presence and size of X-ray cavities, we were unable to rigorously establish this correlation. This is because the number of pixel pairs at separations approaching the domain size is low, leading to increased statistical uncertainty and significant scatter in the VSF values.

Fig. 14. Projected cold gas VSFs including smoothing from atmospheric seeing at three different times in the simulation: 1.28 Gyr, 3.28 Gyr, and 3.48 Gyr (orange, red, and purple lines). For comparison, we also include the 11 VSFs of nearby BCGs from Ganguly et al. (2023) (thick blue lines and red dashed line).

5. Discussion

5.1. On the steepening of the velocity structure functions

Although our results suggest that smoothing from atmospheric seeing may artificially steepen VSFs, they do not rule out the possibility of an intrinsic steepening in either the cold or hot phase VSFs. As shown in Fig. 9, the impact of smoothing depends on the intrinsic power-law index of the unsmoothed VSF. Since VSFs with intrinsically steeper slopes experience weaker steepening from atmospheric seeing, it is possible that VSFs measured from observations of nearby BCGs are intrinsically steeper than the Kolmogorov scaling law of ℓ^1/3, artificially steepened by smoothing, or a combination of both.

Also, several arbitrary parameter choices in our model could significantly impact gas dynamics and, consequently, the resulting VSFs. As discussed in Sect. 2, our radiative cooling algorithm assumes a constant and homogeneous metallicity of the gas. This likely leads to an overestimation of the cooling rate at radii beyond ∼25–50 kpc, which in turn may result in unrealistically extended cold gas structures and eventually excessive radial velocities. Additionally, our feedback parameters–such as jet density–could lead to a different jet-ICM coupling than seen in other simulations. For instance, Beckmann et al. (2019) demonstrated that precessing low-density jets can strongly interact with the cold phase, producing effects not observed in our setup, such as the coherent outward motion of entire cold gas structures. However, it seems unlikely that the specific choices of our setup alone can fully explain the effects reported in, for example, Figs. 12 and 13, which we attribute to observational limitations. The effect of smoothing on projected maps has been investigated in recent observational studies, which conclude that atmospheric seeing causes a global steepening effect on the VSFs (see notably Chen et al. 2022 and Li et al. 2025). The absence of magnetic fields in a purely hydrodynamical setup might also lead to steeper VSFs in the cold phase due to the preferentially disky morphologies reported in many previous studies (e.g., Li & Bryan 2014; Beckmann et al. 2019; Qiu et al. 2019; Fournier et al. 2024). A qualitative comparison between projected cold gas VSFs from the purely hydrodynamical and the fiducial MHD run of the XMAGNET suite is presented in Appendix A. We note that recent observational work focusing on jellyfish galaxies has shown that removing rotational components and large-scale motion can lead to a significant reduction in the observed steepness of the VSFs (Chen et al. 2022; Ignesti et al. 2024). We also emphasize that our results are unlikely to be affected by our choice for the initial velocity perturbations. We have performed several low resolution test simulation and found out that the kinematics of the cold phase is impacted by the initial conditions only during the first ∼100 Myr of the active phase (i.e. starting 500 Myr after the start of the simulation). We thus intentionally included only data from snapshots with times t ≥ 900 Myr. It is, however, likely that lower initial magnetic field strengths will result in quantitatively different behavior for the cold phase, closer to purely hydrodynamical setup (see results for the weakB run in Grete et al. 2025, and associated discussions).

Moreover, it has been shown that numerical viscosity affects the motion of the simulated fluid (Grete et al. 2023b) and that velocity differences computed between pairs of cells located less than 8–10 cells apart (i.e., up to a physical scale of ∼1 kpc in our setup) are thus likely to be underestimated. This is likely to impact the unprojected VSFs presented in Sect. 3 since the separations between pairs of cells are computed in 3D, and we attribute the significant steepening of the hot VSFs at separations lower than 1 kpc to this effect (see Fig. 4). Evaluating how much this might affect the projected VSFs presented in Sect. 4 is not straightforward. Each pixel of the underlying projected velocity maps corresponds to a single velocity value obtained from the averaging of up to hundreds of cells along the line-of-sight, potentially located tens of kpc away from each other. The two-dimensional distance between any given pair of pixels from a projected map is thus not indicative of the underlying three-dimensional distribution of gas cells. This may result in a dilution of the numerical viscosity, and we could expect projected VSFs to be less sensitive to this bias. This effect although remains hard to quantify and we postpone such analysis to future work.

5.2. Comparison to other studies

In Mohapatra et al. (2022), periodic box simulations with turbulence driving are used to study the effects of radiative cooling, magnetic fields, as well as projection and emission weighting on first- and second-order VSF of the hot and cold phases of the ICM. One key result of their work is that the three-dimensional hot phase VSF is always steeper than the Kolmogorov prediction. Our analysis of the first-order VSF supports such conclusions, as we find a scaling close to ℓ^1/2 near the driving scale, and even steeper at smaller scales (i.e., below ∼1 kpc). However, we find that the cold phase moves differently owing to the radial stratification in the gravitational field of the cluster. Consequently, it is less coupled to the hot phase, and dominated by a preferentially infalling motion toward the center. A direct consequence is that we find the cold phase VSF to be always shallower than the hot phase VSF, and slightly steeper than the Kolmogorov scaling of ℓ^1/3.

A simulation setup closer to ours is presented in Wang et al. (2021). In this work, an idealized cool-core cluster with AGN feedback, radiative cooling, and magnetic field is followed. Static grid refinement is used and the VSFs are calculated for gas located at distances less than 10 kpc away from the center of the box. The two main differences with respect to our setup are: (i) the grid structure: our maximally refined region has a width of 250 kpc versus ∼20 kpc in Wang et al. (2021), allowing us to quantify turbulence at much larger radii without refinement level transition effects or varied scales of numerical viscosity across levels, and (ii) the AGN feedback mixture: we used a three feedback channels, namely thermal, kinetic, and magnetic, while the AGN feedback is purely kinetic in Wang et al. (2021). In this work, the authors conclude that the motion of the cold filaments is mostly driven by gravitational acceleration, possibly slowed down by magnetic fields in the MHD case. The authors also find a correlation between the AGN activity and the power-law index of the hot phase VSF. Quiescent periods give rise to shallower VSFs as compared to active periods. Our results are broadly consistent with the conclusions of this study. The amplitudes and scalings of our line-of-sight VSFs are close to that of their fiducial MHD run (see Fig. 4 in Wang et al. 2021 and the dashed lines in Fig. 4 of our study for a direct comparison). As reported in Grete et al. (2025), the AGN activity in our setup does not exhibit significant duty cycles as it is mostly switched on throughout the entire duration of the simulation, with a power of P_AGN ∼ 10⁴⁵ erg s⁻¹. Thus, we do not observe the same alternation between steeper and flatter VSFs correlated with active and quiescent phases, as shown in Fig. 7 of Wang et al. (2021). However, we find evidence that a time-dependent steepening of the VSF is present when the selection of pixel pairs is biased toward regions around rising cavities (see Fig. 6).

5.3. Missing physics

Our work presents a simplified and idealized modeling of a cool-core cluster. We have ignored several phenomena and components, which could impact our results. These include the following:

• The contribution of orbiting cluster galaxies and sloshing to the local turbulence in the core of the cluster: in our setup, we only model the BCG, although the motion of the other galaxies of the cluster is known to contribute to the overall level of turbulence inside the cluster (Fielding et al. 2020; Ayromlou et al. 2024).

• The contribution of cosmic rays, which may affect the self-regulation of cool-core clusters (Jacob & Pfrommer 2017; Beckmann et al. 2022) and possibly the formation and morphology of cold structures (Huang et al. 2022).

• The impact of a radially dependent metallicity profile: we consider an ICM with a constant metallicity of 1.0 Z_⊙. While this reproduces well metallicities measured in the inner tens of kpc of cool-core clusters (Sanders & Fabian 2007; Mernier et al. 2017; McDonald et al. 2019), we might overestimate cooling at larger distances (typically at radius r ≥ 25 − 50 kpc), resulting in cold structures that are more elongated than measured in real cool-core clusters.

• Uncertainties on the AGN feedback parameters: many jet parameters remain largely unconstrained, such as the jet material density and temperature. Changing our jet parameters might impact the way AGN outflows propagate (e.g., Weinberger et al. 2023), and consequently the kinematics of the cold phase and the development of turbulence (e.g., Beckmann et al. 2019).

While these simplifications allow us to isolate and study specific physical mechanisms, they also introduce limitations to the general applicability of our results. Future work incorporating these additional processes will be essential to build a more complete and realistic model of cool-core cluster dynamics.

6. Conclusions

We have performed MHD simulations of AGN feedback in a cool-core cluster with an unprecedented combination of high resolution (Δ_x = 97 pc) and large grid size (2560³ cells in the inner region). This simulation was run on Frontier, the first exascale supercomputer available to academic researchers (Atchley et al. 2023), and performed with the AthenaPK code. Our setup includes gravity, radiative cooling, and a model of cold gas-triggered AGN feedback that converts accreted gas into bipolar magnetized kinetic jets and thermal energy. The dimensions and resolution of the innermost maximally refined grid allowed us to study turbulence while minimizing numerical effects from transitions between refinement levels. Our main conclusions are summarized as follows:

• Based on quantifying the typical line-of-sight and radial velocities of the hot and cold phases, we find that the two phases have significantly different kinematic properties. The hot phase is close to hydrostatic equilibrium, with the velocity dispersions lower than what is measured in clusters with X-ray calorimeters by more than a factor of two. The cold phase follows a ballistic motion and is characterized by infalling and outflowing motions with velocities of up to ±10³ km s⁻¹.

• In unprojected VSFs, we find no clear correlation between the VSF properties of the hot and cold phases. The power-law index of the hot phase VSF varies only weakly over time and remains around m ∼ 0.4. In contrast, the power-law index of the cold-phase VSF fluctuates significantly over time, ranging from m ∼ 0.2 to m ∼ 0.6, with a time-averaged value of m ∼ 0.4.

• We find that the VSF of the hot phase around a rising cavity exhibits a transient steepening, likely related to the injection of turbulence. The power-law index increases by a factor of two, reaching up to m = 0.6. This steepening lasts for a period of ∼50 Myr, after which the VSF reverts to m = 0.4.

• Based on our analysis using mock optical images, we conclude that the cold phase VSF is significantly impacted by projection effects. Projection causes substantial flattening of the VSF, and weighting by Hα emission further enhances this effect.

• Gaussian smoothing of the projected maps modeling the effect of atmospheric seeing significantly affects the projected VSFs at any scale. For scales below and on the order of the smoothing kernel size, the VSFs significantly steepen. Intrinsically shallower VSFs are more impacted by smoothing than steeper VSFs. For all separation scales, the amplitude of the VSFs is decreased by a factor of about two.

• Changing the viewing angle leads to a variation of the amplitude and power-law indices of the projected VSFs by more than a factor of two.

• Projection effects reduce the overall amplitude of the hot-phase VSF compared to its unprojected counterpart and cause a slight steepening.

We conclude that projection effects, atmospheric seeing, and the dependence on the viewing angle must be taken into account when interpreting the VSF of the multiphase ICM. Possible extensions of this work include examining the influence of AGN feedback parameters on the VSFs of both the hot and cold phases, exploring potential relationships between VSF features, viewing angle, and the orientation of cold filaments. Incorporating additional sources of turbulence, such as sloshing, would be also be useful to evaluate the relative contribution of the AGN-induced turbulence.

Data availability

Movie associated to Fig. 1 is available at https://www.aanda.org

Acknowledgments

MF thanks Ricarda Beckmann, Tirso Marin-Gilabert and Yuan Li for insightful discussions, as well as Jacob Shen for providing help with visualisation packages. The authors thank Muzi Li and Brian McNamara for sharing a pre-publication draft of their manuscript on velocity structure functions in multiphase galaxy cluster cores. MB and MF acknowledge 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”). This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). BWO acknowledges support from NSF grants #1908109 and #2106575, NASA ATP grants NNX15AP39G and 80NSSC18K1105, and NASA TCAN grant 80NSSC21K1053. 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 W. O’Shea). DP is supported by the Royal Society through UKRI grant RF-ERE-210263 (PI: Freeke van de Voort). The authors also gratefully acknowledge the Gauss Centre for Supercomputing e.V. (https://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). All simulations were performed using the public MHD code ATHENAPK, which makes use of the KOKKOS (Trott et al. 2022) library and the PARTHENON adaptive mesh refinement framework (Grete et al. 2023a). All data analysis was performed with YT (Turk et al. 2011, 2024), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SEABORN (Waskom 2021) and CMASHER (van der Velden 2020). We thank their authors for making these software and packages publicly available.

References

Appendix A: Comparison with hydrodynamical setup

For brevity, the results section of this paper does not discuss VSFs calculated from the purely hydrodynamical run of the XMAGNET suite. However, it is worth noting that because the structure of the cold phase differs significantly depending on the presence or absence of magnetic fields (see Grete et al. 2025), the resulting VSFs are also markedly different. In Fig. A.1, we present the time-averaged projected VSFs of the cold gas for our hydrodynamical (red curves) and MHD (gray curve) runs. In the hydrodynamical case, the cold gas resides in a massive disk with a radius of r_disk ∼ 5 − 10 kpc, which persists until t ∼ 2.7 Gyr. After this time, the disk disappears as thousands of cold clumps condenses along the jet axis, increasing the total mass of cold gas in the simulation by nearly two orders of magnitude, up to ∼10¹¹ M_⊙ (see Fig. 3 in Grete et al. 2025). The later phase is characterized by ballistic motion closer to what observed in the fiducial MHD run. We distinguish these two phases using two time-averaged VSFs, with the solid line corresponding to the disky phase and the dashed line to the non-disky phase.

Fig. A.1. Comparison between the time-averaged projected VSFs of our hydrodynamical (red lines) and MHD (gray line) runs. In the hydrodynamical run, the gas is nearly entirely contained in a disk of radius rdisk ∼ 5 − 10 kpc for times t ≤ 2.7 Gyr. We thus distinguish this specific phase (solid line) from the remaining part of the run (dashed line) where no disk is visible (i.e., for t ≥ 2.7 Gyr).

Because the disk is rotating along an axis broadly aligned with the z-axis of the simulation’s coordinate system, orthogonal to the projection axis, the resulting projected velocity maps exhibit a distinct pattern where half of the disk has positive line-of-sight velocities and the other half has negative line-of-sight velocities. This results in a characterized peak in the VSFs, visible in the solid line at separation scales broadly consistent with the disk’s size (i.e., r_disk ∼ 5 − 10 kpc). The time-averaged VSF at later times (i.e., when no disk is visible) does not exhibit any peak. Its scaling is broadly consistent with that of the MHD run (gray curve), although its mean amplitude is larger by a factor of ∼2.

Appendix B: Velocity structure function algorithm convergence

Figure B.1 presents the three-dimensional VSFs of the hot and cold phases at t = 3.28 Gyr, considering all gas particles within a radius of 12 kpc. For the hot phase, we compare the VSF obtained using a naive sampling method—where all cell pairs are randomly sampled—with our optimized approach, which prioritizes alternative sampling for separations below 500 pc (as detailed in Sect. 2.2). As visible, the two methods are in good agreement. We find that VSF convergence is achieved with 10⁷ cell pairs for the hot phase and 10⁸ for the cold phase. However, the volume filling fraction of the cold phase varies significantly throughout the simulation, particularly within the inner 12 kpc, making the number of pairs required for convergence time-dependent. This number also changes with the selected domain size. After testing convergence across multiple snapshots, we adopt 10⁷ (10⁸) pairs for the hot (cold) phase when analyzing an r < 12 kpc region. For larger domains (r < 25 or r < 50 kpc), we use 10⁸ (10⁹) pairs for the hot (cold) phase, respectively.

Fig. B.1. Convergence of the VSF algorithm taking various numbers of pairs of cells for both the hot (left panel) or cold (right panel) phases and for the gas located at radii r ≤ 12 kpc at t = 3.28 Gyr. The labels “os” and “ns” designate the optimized and naive sampling methods.

Appendix C: Comparison with second-order velocity structure functions

In Fig. C.1, we present a comparison between the first-order (blue curve) and second-order (purple) VSFs for the hot (left panel) and cold (right panel) phases, calculated for the gas cells within r ≤ 12 kpc. All curves have been rescaled such that a pth-order VSF following Kolmogorov scaling, ie. with a power-law index m = p/3, appears as a horizontal line in the plot. The Kolmogorov theory of turbulence states that S₂(ℓ) and the longitudinal VSF S₂^∥(ℓ) should scale similarly. For comparison, we also present the time-averaged second-order longitudinal VSF (red curve), defined as

Fig. C.1. Comparison of the time-averaged first-order VSF (blue curve), second-order VSF (purple curve), and second-order longitudinal VSF (red curve) for the hot phase (left panel) and the cold phase (right panel). The separations scales have been normalized by the domain size L, and the VSFs have been normalized by their time-averaged root mean square velocity vrms. The colored areas represent the standard deviation resulting from time variation across all snapshots.

$$\begin{array}{c}\hfill S_2^{\Vert }(\ell )=⟨\{\mathbf{v}({\text{r}}_j+\ell ·\mathbf{e})-\mathbf{v}({\text{r}}_j)\}·\mathbf{e}⟩,\end{array}$$(C.1)

where e is the unit vector pointing from cell j to cell i. The scaling of the first- and second-order VSFs are in good agreement.

Appendix D: Projection depth and smoothing

In Fig. D.1, we present the power-law index of VSFs measured at separation scales of 600 pc (left panel) and 3 kpc (right panel), as a function of the projection depth and smoothing length. Each pixel in the figure represents the power-law index of a VSF computed from a projected velocity map. The projection is smoothed using a 2D Gaussian filter with a standard deviation σ_filter. The velocity maps include all cold gas cells within a projection depth Δx_proj, defined by the condition −Δx_proj/2 ≤ x ≤ Δx_proj/2, where x is the coordinate along the line-of-sight (with the SMBH particle being located at x_SMBH = 0.0). The horizontal gray dashed line indicates the minimum depth of projection beyond which at least 95 % of all the gas cells contained in the simulated box are included in the projection.

Fig. D.1. Various power-law indices measured by fitting the VSFs over ℓ = 0.6 ± 0.1 kpc (left panel) and ℓ = 3 ± 1 kpc (right panel) as a function of the Gaussian smoothing kernel size (expressed in pixels in the lower axis and in physical length for the upper axis, assuming the distance to the Perseus cluster for the later) on the horizontal axis and the depth of projection Δxproj along the vertical axis. The vertical orange dashed line indicates the standard deviation σPerseus of the smoothing effect for the Perseus cluster. The inner figure within the left panel is a 1D visualization of the uppermost row.

As the projection depth increases more cold structures are included in the projection, modifying the power-law index of the VSFs. This effect is particularly visible in the right panel—that is, for the power-law index measured at 3 kpc. The power-law index measured at 600 pc, the effect of the smoothing is nearly independent of the projection depth, and more significant than at 3 kpc. A smoothing length of 2–3 pixels results in a doubling of the VSFs’ power-law index. The inner figure in the left panel is a 1D visualization of the top row (including all cold gas cells of the box), showing the evolution of the VSFs’ power law as a function of the smoothing kernel size. For smoothing lengths larger than 4 pixels, the index converges to m ∼ 1.

Appendix E: Coarse-grained velocity maps

In Sect. 4.3.1, we compute the VSF extracted from projected velocity maps with degraded resolution. To degrade the resolution of the velocity maps, we employ a downsampling technique where the maps are coarsened by a factor N. Specifically, we group adjacent pixels into N × N blocks, compute the mean velocity for each block, and assign the result to a single pixel in the downsampled map. During this process, we ensured that blocks with too many empty pixels (i.e., pixels corresponding to the projection of a column containing no cold gas cell) were excluded from the mean computation by applying a threshold s = 0.5 on the proportion of non-empty pixels per block. A comparison of the resulting maps is presented in Fig. E.1. This threshold method allows one to ignore small clumps of one or a few cells, which are visible in the high-resolution maps but expected to be spatially unresolved in the degraded map.

Fig. E.1. Comparison of three projected line-of-sight velocity maps of the same snapshot and with the same orientation. The panels indicate, from left to right, maps with pixel widths of 0.1 kpc (the maximum resolution of our simulation), 0.6 kpc, and 1.1 kpc.

All Figures

Fig. 1. Hα emission – weighted projection of the line-of-sight velocity for the inner 80 kpc of our simulated box at t = 3.2 Gyr. A zoom-in on the innermost region showcases the effect of atmospheric seeing on the map. The background image is a projection of the temperature field. The figure can be compared to Fig. 2 from Li et al. (2020). ⊙ and ⊗ markers indicate the directions pointing toward (away) from the viewer, associated with negative (positive) line-of-sight velocity (respectively). More detail on the emissivity weighting are provided in Sect. 4.1. A movie showing the three dimensional structure of the cold phase is available online.

In the text

Fig. 2. Probability density function (PDF) of the velocities’ x-component for the hot phase (106 ≤ T ≤ 108 K, left panel) and cold phase (T ≤ 105 K, right panel) for various radius bins and the PDF of all gas located at radii r ≤ 100 kpc, averaged between 0.9 and 4 Gyr. Vertical dashed lines indicate the median values. Colored areas represent the standard deviation resulting from time variation. We have verified that the y-component leads to similar results.

In the text

	Fig. 3. Same as Fig. 2 but for radial velocity, vr.
In the text

Fig. 4. Velocity structure functions of the hot ICM (left panel) and the cold phase (right panel) for all gas located less than 12, 25, and 50 kpc away from the center of the simulated box. The dashed colored lines indicate the line-of-sight VSFs taking the x-component of the velocity field. We emphasize that the latter are not equivalent to projected VSFs, as the distance between pairs of cells is still three-dimensional (projected VSFs are discussed in Sect. 4). Colored area is showing standard deviation over time. The solid gray lines in the background indicate the Kolmogorov scaling law (Kolmogorov 1941). As our setup closely resembles the fiducial MHD run from Wang et al. (2021), we present the time-averaged line-of-sight VSFs for the hot and cold phases from that work, for radii r ≲ 10 kpc (dashed gray line labeled W21). A comparison of our results with this study is discussed in Sect. 5.2. We emphasize that VSFs below separations of ℓ ∼ 800 pc are likely to be steepened by the effect of numerical viscosity.

In the text

Fig. 5. Top panels: Relationship between the power-law index of the line-of-sight hot gas VSF and that of its three-dimensional counterpart. Middle (bottom) panels: Relationship between the power-law index m (amplitude v0) of the line-of-sight, cold phase VSF, and the three-dimensional hot phase VSFm, respectively. The color of the points represents time, with values ranging from 0.9 to 4.0 Gyr. The columns show VSF values for all gas within radii of 12, 25, and 50 kpc (left to right). The dashed gray line indicates the identity relationship and should not be interpreted as an attempt to fit the data. The light blue shaded area represents the region of the parameter space where the measured power-law index of the VSF is steeper than the m = 1/3 power-law index predicted by the Kolmogorov theory.

In the text

Fig. 6. Top panels: temperature projection of the gas. The dashed circle indicates the sphere within which the VSF are computed, whose position is kept fixed throughout the run. The round marker at the top of the temperature projection indicates the position of the SMBH particle and the arrow shows the jet injection axis. The projection depth is 1 kpc. Middle panels: Velocity structure functions of the cold (gray curves) and hot (red curves) phases. The dashed orange lines indicate the fit of the power-law index of the hot phase VSF. Lower panel: Time evolution of the cold and hot phases power-law index. Here, t0 refers to the restarting time of the simulation.

In the text

	Fig. 7. Emission-weighted (w) line-of-sight velocity projection of the hot phase within the inner 50 kpc of our fiducial MHD run at t = 3.1 Gyr. The marker and dashed line at the center indicate the position of the SMBH particle and the axis of the jets, respectively.
In the text

Fig. 8. Line-of-sight VSFs for the cold component and for the projected radius r ≤ 25 kpc averaged between 0.9 and 4 Gyr. The darkest curve (labeled ℓ3D) was obtained by taking the three-dimensional distance between pairs of points, similarly to Wang et al. (2021). For all the other curves, the VSFs were computed out of projected maps, and the corresponding separations scales are two dimensional (ℓ2D). The term w(Hα) indicates weighting of the map by Hα emission, and 𝒮 is the additional smoothing effect resulting from atmospheric seeing (see Li et al. 2020). The shaded areas indicate the time variation of each VSF.

In the text

	Fig. 9. Relationship between the power-law index $m_{\mathrm{cold}}^{\mathrm{us}}$ of the unsmoothed (“us”) projected VSFs (horizontal axis) and Δm, which represents the increase in the power-law index resulting from smoothing. Lines of equal power-law index for smoothed VSFs ($m_{\mathrm{cold}}^{\mathrm{s}}$) are also indicated for comparison (gray lines).
In the text

Fig. 10. Velocity structure functions obtained from projected maps taking various spatial resolutions for the cold gas located at a projected radius of r ≤ 25 kpc. The projection method for all curves is equivalent to the third curve in Fig. 8 (i.e., labeled as ℓyz + w(Hα)). Consequently, the VSFs presented here are shallower than their three-dimensional counterparts (labeled as ℓ3D in Fig. 8) due to projection and emission weighting. The darkest curve corresponds to the maximum spatial resolution Δx, 0 of our simulation. All brighter curves are obtained from coarse-grained maps with pixel size Δx. The colored areas represent the standard deviation resulting from time variation across all snapshots.

In the text

	Fig. 11. Same as Fig. 8 but for the hot phase. The brightest curves indicate the VSF obtained from weighted and projected velocity maps, taking three different minimum temperature cutoffs. The projected VSFs are steeper than the unprojected one for separations above ℓ ∼ 1 kpc.
In the text

Fig. 12. Aitoff projection of the values of the VSF amplitude v0 and the power-law index m, measured for 2000 randomly selected orientations at fixed time t = 3.2 Gyr within r = 25 kpc. Each VSF was computed from projected Hα emission weighted maps. The power-law fitting was performed between 2 and 4 kpc. As we do not include effects of opacity when generating the projection of the velocity maps used to obtain each VSF, the left and right half of the plot are symmetrical for both v0 and m. The coordinates of the jet injection axis is at θ = ±90°.

In the text

Fig. 13. Time-averaged VSFs of the cold phase obtained from mock Hα emission-weighted maps projected along the x-axis (red lines) compared to observational data (gray lines, from Li et al. 2020). The images have been coarse-grained to a spatial resolution of 290 pc to mimic the spatial resolution of instruments assuming Perseus cluster distance. The right panel was obtained from smoothed maps to mimic the effects of atmospheric seeing. The labels “inner”, “outer”, and “all” (see legend) respectively designate the VSFs obtained from all gas cells located at a projected radii of less than 12 kpc, at a projected radius of more than 12 kpc, and for all gas cells irrespective of their radius.

In the text

	Fig. 14. Projected cold gas VSFs including smoothing from atmospheric seeing at three different times in the simulation: 1.28 Gyr, 3.28 Gyr, and 3.48 Gyr (orange, red, and purple lines). For comparison, we also include the 11 VSFs of nearby BCGs from Ganguly et al. (2023) (thick blue lines and red dashed line).
In the text

	Fig. A.1. Comparison between the time-averaged projected VSFs of our hydrodynamical (red lines) and MHD (gray line) runs. In the hydrodynamical run, the gas is nearly entirely contained in a disk of radius rdisk ∼ 5 − 10 kpc for times t ≤ 2.7 Gyr. We thus distinguish this specific phase (solid line) from the remaining part of the run (dashed line) where no disk is visible (i.e., for t ≥ 2.7 Gyr).
In the text

	Fig. B.1. Convergence of the VSF algorithm taking various numbers of pairs of cells for both the hot (left panel) or cold (right panel) phases and for the gas located at radii r ≤ 12 kpc at t = 3.28 Gyr. The labels “os” and “ns” designate the optimized and naive sampling methods.
In the text

Fig. C.1. Comparison of the time-averaged first-order VSF (blue curve), second-order VSF (purple curve), and second-order longitudinal VSF (red curve) for the hot phase (left panel) and the cold phase (right panel). The separations scales have been normalized by the domain size L, and the VSFs have been normalized by their time-averaged root mean square velocity vrms. The colored areas represent the standard deviation resulting from time variation across all snapshots.

In the text

Fig. D.1. Various power-law indices measured by fitting the VSFs over ℓ = 0.6 ± 0.1 kpc (left panel) and ℓ = 3 ± 1 kpc (right panel) as a function of the Gaussian smoothing kernel size (expressed in pixels in the lower axis and in physical length for the upper axis, assuming the distance to the Perseus cluster for the later) on the horizontal axis and the depth of projection Δxproj along the vertical axis. The vertical orange dashed line indicates the standard deviation σPerseus of the smoothing effect for the Perseus cluster. The inner figure within the left panel is a 1D visualization of the uppermost row.

In the text

	Fig. E.1. Comparison of three projected line-of-sight velocity maps of the same snapshot and with the same orientation. The panels indicate, from left to right, maps with pixel widths of 0.1 kpc (the maximum resolution of our simulation), 0.6 kpc, and 1.1 kpc.
In the text