© The Authors 2025

1. Introduction

As the culmination of the structure formation process, galaxy clusters are sensitive cosmological probes through their dependence on the high-mass end of the halo mass function (e.g., Allen et al. 2011; Pierre et al. 2011; Clerc & Finoguenov 2023). Surveys of increasing size and depth such as eROSITA (Ghirardini et al. 2024), Euclid (Euclid Collaboration: Adam et al. 2019; Euclid Collaboration: Bhargava et al. 2025), SPT (Bocquet et al. 2024), or ACT (Hilton et al. 2021) are now detecting many thousands of clusters, which can potentially lead to groundbreaking discoveries. At the present time, the main bottleneck for these surveys lies in the characterization of the detected clusters, and specifically on the relation between the survey observables and the mass of the underlying halo (see Pratt et al. 2019, for a review). In this respect, observables based on the properties of the intracluster medium (ICM) such as temperature, gas mass, or Compton-y parameter are of particular interest, as they are expected to be tightly correlated with the halo mass (e.g., Pratt et al. 2009; Ettori 2015; Lovisari et al. 2020; Sereno et al. 2020). The primary unknown is the impact of structure formation processes and baryonic physics on these observables (McCarthy et al. 2010; Le Brun et al. 2014; Truong et al. 2018; Pop et al. 2022; Braspenning et al. 2024), which currently prevent observables from being derived from the first principles in a straightforward manner.

In this respect, measurements of the gas velocity field in galaxy clusters hold a key role. Numerical simulations have indeed long predicted that a fraction of the energy in the ICM is not thermalized (Rasia et al. 2006; Lau et al. 2009; Nelson et al. 2014; Biffi et al. 2016; Ayromlou et al. 2024). Given the low density of the medium, the mean-free path in the ICM is expected to be on the order of a kiloparsec, although small-scale instabilities may reduce the mean-free path substantially (Schekochihin 2022). As such, the kinetic energy injected by successive merging events takes several gigayears to thermalize (see Simionescu et al. 2019, for a review). In addition, outflows from central active galactic nuclei (AGN) introduce perturbations in the neighboring medium by generating shocks, carving cavities, and driving turbulence, which affects the thermodynamic properties of the ICM and its gas distribution (Gaspari et al. 2011; Gaspari 2015; Bourne & Sijacki 2017; Wittor & Gaspari 2020; Bourne & Sijacki 2021). Numerical simulations predict that most of the remaining kinetic energy should be in the form of random gas motions (Lau et al. 2009; ZuHone et al. 2013; Vazza et al. 2017; Angelinelli et al. 2020; Pearce et al. 2020; Barnes et al. 2021; Bennett & Sijacki 2022; Groth et al. 2025). Perturbations injected on large scales generate a turbulent cascade and dissipate on small scales where thermalization becomes effective (e.g., Gaspari et al. 2014; Zhuravleva et al. 2014). The statistical properties of the velocity field can be characterized by a fluctuation power spectrum that peaks on large scales at the so-called injection scale and decreases on smaller scales as a power law with an index close to the classical Kolmogorov slope of 11/3 (Schuecker et al. 2004; Churazov et al. 2012; Gaspari & Churazov 2013; Vazza et al. 2017; Zhuravleva et al. 2018; Simonte et al. 2022).

Observational constraints on the velocity fluctuation power spectrum are currently scarce due to the lack of instruments with appropriate capabilities. Studies of line profiles in XMM-Newton/RGS observations of cool-core clusters only yielded mild constraints on the turbulence of the ICM (Sanders et al. 2011; Sanders & Fabian 2013; Pinto et al. 2015; Ogorzalek et al. 2017). More recently, Sanders et al. (2020) proposed a new technique to accurately calibrate the energy gain of the XMM-Newton/EPIC-pn CCD, which resulted in measurements of bulk velocities in a select number of nearby clusters (Gatuzz et al. 2022, 2023a, 2024), although the accuracy of the method is limited to ∼100 km/s. Alternatively, several works used the fluctuation field of thermodynamic quantities (density, pressure, and temperature) as an indirect tracer of the velocity field (Schuecker et al. 2004; Churazov et al. 2012; Gaspari et al. 2014; Zhuravleva et al. 2014; Hofmann et al. 2016; Khatri & Gaspari 2016). These studies imply that the turbulent pressure support in galaxy clusters is relatively mild (P_turb/P_tot ≲ 0.2, Zhuravleva et al. 2015; Eckert et al. 2017; Dupourqué et al. 2023, 2024; Romero et al. 2023, 2024; Heinrich et al. 2024; Marin-Gilabert et al. 2024; Lovisari et al. 2024). However, the relation between the thermodynamic fluctuations and the underlying velocity field is uncertain, most notably in merging clusters (Simonte et al. 2022), in large part due to contamination by other unrelated processes such as subhalos, gas sloshing, or ellipticity (Zhuravleva et al. 2023).

The field is expected to advance significantly thanks to the advent of X-ray micro-calorimeters. The first observation of a galaxy cluster with a micro-calorimeter was obtained in 2016 by the SXS instrument on board Hitomi, which targeted the Perseus cluster as its first-light observation (Hitomi Collaboration 2016, 2018). In the cluster core, the Fe XXV emission line complex was found to be broadened by 160 km/s, which, if interpreted as evidence for random isotropic motions, amounts to a small fraction (∼4%) of the thermal pressure (Hitomi Collaboration 2016). While Hitomi was unfortunately lost shortly thereafter, its successor XRISM (Tashiro et al. 2025) was launched in September 2023 and has been successfully operating since then. XRISM carries the Resolve instrument, an array of 6 × 6 individual pixels featuring a spectral resolution of 4.5 eV. XRISM has observed several clusters thus far (XRISM Collaboration 2025a,b,c), revealing moderate line widths (100−200 km/s) but substantial bulk motions as large as 700 km/s in Coma (XRISM Collaboration 2025a, hereafter X25). However, linking these measurements to the underlying velocity power spectrum is far from trivial, which makes the comparison with predictions from numerical simulations difficult.

In this paper, we propose a simulation-based inference (SBI) method to determine the power spectrum of velocity fluctuations in galaxy clusters based on X-ray micro-calorimeter data. We generate simulations of the velocity field for any given power spectrum shape, and fold in observational effects (emissivity weighting, projection, point spread function (PSF) smearing, measurement uncertainties) to generate mock bulk velocities and velocity dispersions for any observational setup. We then train a neural network to learn the relation between the parameters describing the 3D power spectrum and the mock datasets. The neural network can then be used in reverse to infer the model parameters given the observed data. We validate our method using extensive simulations and apply it to the measured bulk velocities and velocity dispersions in the Coma cluster. The paper is organized as follows. In Sect. 2.1 we describe the simulation pipeline and the optimization procedure. We validate our method in Sect. 2.3 and apply it to the XRISM/Resolve data of the Coma cluster in Sect. 3.1, comparing the results with the constraints obtained from the same data in X25. We discuss our results in Sect. 4.

Throughout the paper, we assume a flat ΛCDM cosmology with Ω_m = 0.3 and H₀ = 70 km/s/Mpc. At the mean redshift of the Coma cluster (z = 0.0233, Bilton & Pimbblet 2018), this corresponds to 1′ = 28.5 kpc.

2. Methodology

In this section we present the theoretical framework that we use to forward model the velocity measurements and our SBI implementation. We then validate the technique using a large set of simulations and study the ability of our method to recover the parameters of the underlying velocity power spectrum. We distribute our code in the form of an easy-to-use Python package¹.

Fig. 1. Example of fluctuation field generation. Left: Model 3D power spectrum described by an injection at kinj and a classical Kolmogorov cascade on smaller scales (Eq. 4). The vertical dashed line shows the position of kinj = (500 kpc)−1. Right: Slice of the real-space velocity fluctuation field along the z axis generated from this power spectrum (Eq. 3) for a velocity dispersion σv = 300 km/s. The color code shows the velocity in km/s. The box size is 1.7 Mpc on a side.

2.1. Simulation setup

We consider a three-dimensional velocity field v(x) = (v_x(x, y, z),v_y(x, y, z),v_z(x, y, z)), with x, y the Cartesian coordinates on the plane of the sky and z the coordinate along the line of sight. The three-dimensional power spectrum of the velocity field is defined as the square modulus of the velocity field in Fourier space,

$$\begin{array}{c}\hfill P_{3\mathrm{D}}(\mathbf{k})={\left|{\displaystyle \int \mathbf{v}}(\mathbf{x})e^{-i2\pi \mathbf{k}·\mathbf{x}}{\mathrm{d}}^{3}x\right|}^2,\end{array}$$(1)

with k = (k_x, k_y, k_z) the wave number. Assuming isotropy, the power spectrum can be reduced to the one-dimensional function P_3D(k), with k = |k| = 1/ℓ.

For a simulation box of size N³, we set linearly spaced wave numbers with the origin set at the center of the box, such that each axis (k_x, k_y, k_z) is linearly spaced and spans values in the range k ∈ ( − N/2, N/2). Given a parametric model for P_3D, we calculate the power P_3D(k) for every point on the grid and generate a random Gaussian realization of the Fourier transform $\widehat{\mathit{v}}$ of the velocity field as

$$\begin{array}{c}\hfill \widehat{\mathit{v}}(\mathbf{k})=P_{3\mathrm{D}}{(\mathbf{k})}^{1/2}{\mathcal{N}}_C,\end{array}$$(2)

with 𝒩_C a random complex number with normally distributed real and imaginary parts. The real-space velocity field can then be obtained as the real part of the inverse Fourier transform of $\widehat{\mathit{v}}$,

$$\begin{array}{c}\hfill \mathbf{v}(\mathbf{x})=\mathrm{Re}\left[{\displaystyle \int \widehat{\mathit{v}}(\mathbf{k})e^{i2\pi \mathbf{k}·\mathbf{x}}{\mathrm{d}}^{3}k}\right].\end{array}$$(3)

This procedure generates a 3D velocity field centered on zero with Gaussian fluctuations representing isotropic turbulent eddies. As an example, in Fig. 1 we show a realization of the velocity field from a power spectrum given by the classical expression for a turbulent cascade generated by injection on large scales (e.g., Zhuravleva et al. 2012; ZuHone et al. 2016),

$$\begin{array}{c}\hfill P_{3\mathrm{D}}(k)={\sigma }_v^2\frac{{(1+{(k/k_{\mathrm{inj}})}^2)}^{-\alpha /2}}{{\displaystyle \int 4\pi }{k}^2{(1+{(k/k_{\mathrm{inj}})}^2)}^{-\alpha /2}\mathrm{d}k}·\end{array}$$(4)

Here k_inj describes the wave number corresponding to the energy injection scale, α is the slope of the turbulent cascade, and the model is normalized such that the amplitude of the spectrum is equal to the variance σ_v² of the velocity field. The left-hand panel of Fig. 1 shows the power spectrum described by Eq. (4) for an injection scale ℓ_inj = (k_inj)⁻¹ = 500 kpc, a slope of 11/3 corresponding to the classical Kolmogorov slope, and an arbitrary normalization. In the right-hand panel we show a slice of the 3D velocity field along the z axis for a realization of a Gaussian random field characterized by this power spectrum, with a velocity dispersion σ_v = 300 km/s, a pixel size of 8 kpc, and a box size of 1.7 Mpc on a side. We can see that the procedure generates fluctuations on all scales with a maximum around ℓ_inj = 500 kpc, as expected from the model of choice.

Given a paramterization for P_3D(k) and a set of model parameters, we generate a real-space fluctuation field v(x) according to the procedure described above. Since measurements of Doppler shifts and velocity dispersions are obtained along the line of sight (z axis), we need to project the velocity field to predict the measured quantities. Let EM_3D(x) be the three-dimensional emissivity distribution, which is proportional to the squared gas density. The projected bulk velocity is given by the mean of the 3D velocity along the z direction, weighted by the local gas emissivity,

$$\begin{array}{c}\hfill v_b(x,y)=\frac{\int {\mathrm{EM}}_{3\mathrm{D}}(\mathbf{x})v_z(\mathbf{x})\mathrm{d}z}{\int {\mathrm{EM}}_{3\mathrm{D}}(\mathbf{x})\mathrm{d}z}·\end{array}$$(5)

Similarly, the line-of-sight averaged velocity dispersion σ(x, y) can be estimated as the quadratic mean of the velocity fluctuations along the z axis, weighted by the local emission measure,

$$\begin{array}{c}\hfill \sigma (x,y)={\left[\frac{\int {\mathrm{EM}}_{3\mathrm{D}}(\mathbf{x}){(v_z(\mathbf{x})-v_b(x,y))}^2\mathrm{d}z}{\int {\mathrm{EM}}_{3\mathrm{D}}(\mathbf{x})\mathrm{d}z}\right]}^{1/2}\end{array}$$(6)

Truong et al. (2024) demonstrated with mock spectra from simulated clusters from the TNG-Cluster suite that Eqs. (5) and (6) provide a good description of the quantities estimated from X-ray spectroscopy (see also Roncarelli et al. 2018).

To describe the emission measure distribution EM_3D(x), we assume that the system of interest is spherically symmetric with a gas distribution described by a spherical βmodel (Cavaliere & Fusco-Femiano 1976),

$$\begin{array}{c}\hfill {\mathrm{EM}}_{3\mathrm{D}}(r)={\mathrm{EM}}_0{(1+{\left(\frac{r}{r_c}\right)}^2)}^{-3\beta }.\end{array}$$(7)

This model provides a good representation of the surface brightness profile of the Coma cluster (e.g., Churazov et al. 2012). Fitting Eq. (7) to the XMM-Newton surface brightness profile of the system yields r_c = 266 kpc and β = 0.66. The choice of the βmodel is motivated by the fact that the relation between 2D and 3D distributions is analytic (e.g., Eckert et al. 2020). However, any model of choice can be substituted for cases in which the βmodel is not sufficiently accurate.

In the presence of a non-negligible instrumental PSF, as is the case for XRISM/Resolve, the measured bulk velocities and velocity dispersions need to be averaged out also along the plane-of-the-sky coordinates (x, y), again weighted by the local emission measure. Equation (5) then becomes

$$\begin{array}{c}\hfill v_b(x\prime ,y\prime )=\frac{\int \mathrm{PSF}(x-x\prime ,y-y\prime ){\mathrm{EM}}_{3\mathrm{D}}(\mathbf{x})v_z(\mathbf{x}){\mathrm{d}}^{3}x}{\int \mathrm{PSF}(x-x\prime ,y-y\prime ){\mathrm{EM}}_{3\mathrm{D}}(\mathbf{x}){\mathrm{d}}^{3}x}·\end{array}$$(8)

Eq. (6) can be modified in a similar way to predict the velocity dispersion in the presence of a relatively broad instrumental PSF. In the case of XRISM/Resolve, we model the PSF as a two-dimensional Gaussian with σ = 34″, which gives a reasonable description of the telescope’s PSF (Tashiro et al. 2025). However, any PSF shape can be considered and folded into Eq. (8).

Finally, the observational configuration can be introduced, i.e. the regions for which bulk velocity and velocity dispersion measurements were obtained can be defined as a function of the plane-of-the-sky coordinates (x, y). The measurements of bulk velocities and velocity dispersions within each region of total area Ω are then obtained by integrating Eq. (8) over the defined area, i.e.

$$\begin{array}{c}\hfill v_{1\mathrm{D}}=\frac{1}{\mathrm{\Omega }}{\displaystyle \int v_b}(x,y)\mathrm{d}\mathrm{\Omega }.\end{array}$$(9)

The final values of v_1D and σ_1D are drawn from Gaussian distributions according to the measurement uncertainties in each region to generate a mock observational dataset which mimics the real procedure as closely as possible.

2.2. Optimization

As our test configuration, we consider the two XRISM/Resolve pointings of the Coma cluster presented in X25. The dataset consists in a deep (400 ks) pointing located close to the X-ray emission peak and a second slightly shallower pointing (160 ks) located 170 kpc south of the cluster center. Each pointing was split into four approximately independent sub-arrays of 1.5 arcmin (42 kpc) on a side, such that the considered dataset includes 8 individual bulk velocity measurements with an uncertainty of 30 km/s for the central pointing and 60 km/s for the south pointing. The bulk velocities are set with respect to the cluster rest frame defined as the mean redshift of ∼600 member galaxies. We consider the velocity dispersions from the four sub-arrays for the deep central pointing, while for the shallower south pointing we only include a single velocity dispersion measurement for the entire array. As a result, our simulation pipeline generates a data vector that contains a set of 8 bulk velocities and 5 velocity dispersions.

The procedure described in Sect. 2.1 is inherently random, as an infinity of real-space fluctuation fields can be described by the same power spectrum. In particular, for two realizations generated from the same parameters, the local bulk velocities vary substantially, as the distribution of local bulk velocities is randomly distributed around zero. This effect is known as cosmic variance (Clerc et al. 2019). For this reason, the data cannot be compared to the model using an analytic likelihood. To alleviate this issue, we adopt simulation-based inference as our fitting technique, as was done by Dupourqué et al. (2023, 2024) in the context of surface brightness fluctuations. Specifically, we use the sbi Python package (Tejero-Cantero et al. 2020) and the sequential neural posterior estimation (SNPE) method to fit the data. For a given parameterization of P_3D(k) (e.g., Eq. 4), we set a uniform prior range for the model parameters and sample parameter sets from the prior distribution. We then generate mock datasets using the simulation pipeline outlined in Sect. 2.1 and train a neural network to describe the relation between the model parameters and the multidimensional posterior density distribution. We use the automatic posterior transformation method (Greenberg et al. 2019) and the masked autoregressive flow (MAF) normalizing flow method (Papamakarios et al. 2017) for density estimation. We refer to these papers for the details on the architecture of the network and the training method. We infer the posterior distribution of the model parameters by sampling the trained density distribution given the observed data vector. For more details on the optimization technique, we refer the reader to Regamey et al. (2025).

Alternatively, following Molin et al. (2025) we consider a summary statistic for the bulk velocity field in the form of the velocity structure function (VSF, e.g., ZuHone et al. 2016; Roncarelli et al. 2018; Gatuzz et al. 2023b). The VSF is defined as the mean squared difference of all velocity measurements separated by a distance r,

$$\begin{array}{c}\hfill \mathrm{VSF}(r)=\langle {(v_b(x,y)-v_b(x+r_x,y+r_y))}^2\rangle ,\end{array}$$(10)

with (r_x, r_y) any 2D vector of modulus r. The VSF can be analytically related to the power spectrum (Zhuravleva et al. 2012; ZuHone et al. 2016; Clerc et al. 2019), such that a likelihood function can be constructed and optimized, as was done in X25 for the Coma cluster. The main drawback is that the VSF only depends on relative bulk velocities rather than absolute ones, thus it is insensitive to bulk motions occurring on scales larger than the covered observational footprint. In the following we consider both the case where we directly forward fit the bulk velocities and the case where we fit a summary statistic in the form of the VSF. In the latter case, for each simulated mock dataset, we calculate the VSF using Eq. (10). The fitted data vector is then made of the combination of the VSF points and the velocity dispersion measurements.

2.3. Recovery tests

We ran recovery tests with random points to test the accuracy of our method and its ability to recover the true underlying power spectrum. We considered the two-pointing configuration described in Sect. 2.2 and an emission measure profile that is well matched to the emissivity distribution in the Coma cluster (Eq. 7, see Appendix A). We considered boxes with N = 200 pixels on a side covering a range of 1.7 Mpc, which is much larger than the projected distance between the two pointings (170 kpc) and the effective depth along the line of sight (∼600 kpc). This configuration also provides a spatial resolution of 8 kpc, which substantially oversamples the instrumental PSF (σ = 34″ ≈ 16 kpc). We assumed that the power spectrum can be described by the three-parameter model shown in Eq. (4) and Fig. 1. The model does not include dissipation on small scales, as the resolution of XRISM/Resolve (1.5′∼43 kpc) largely exceeds the Coulomb mean-free path in the medium (∼10 kpc). The model parameters are the 3D velocity dispersion σ_v, the injection wave number k_inj, and the slope of the turbulent cascade α. Instead of the velocity dispersion, we define the model as a function of the Mach number ℳ_3D = σ_v/c_s, with c_s = (γk_BT/μm_p)^1/2 the sound speed in the medium. For a temperature of 8.37 keV in the core of Coma (see X25), the sound speed amounts to 1508 km/s. We set uniform priors on the model parameters spanning the range ℳ_3D ∼ 𝒰(0, 2), log k_inj ∼ 𝒰(−4, −2), α ∼ 𝒰(2, 8). These prior ranges are expected to be broad enough to encompass any reasonable configuration.

We sampled 10 000 parameter sets from the prior distributions and generated mock data vectors for each parameter set using the simulation pipeline defined in Sect. 2.1. We then used the SNPE method in sbi to estimate the posterior density distribution. The training was achieved by maximizing a loss function using a maximum-likelihood algorithm, which converged after 148 epochs. We generated another set of 1000 simulations and applied the model to the mock data vectors to infer the best fitting parameters. In each case, we adopted the median of the posterior distribution as our point estimate, as it is a robust estimator in the presence of skewed, non-Gaussian posteriors. We compared the resulting parameter values to the input values that were used to generate the simulations. In Fig. 2 we show the resulting coverage tests for the model defined in Eq. (4) and the two-pointing configuration. The gray points show the median of the recovered posterior distributions for each individual simulation, whereas the smooth color scale shows the point density estimated using a Gaussian kernel density estimator. We can see that the method accurately recovers the 3D Mach number, with an uncertainty ranging from ∼0.1 at low Mach numbers to ∼0.25 for ℳ_3D > 1.5. The recovered points are unbiased throughout most of the considered range, with the exception of the very high Mach number regime where the posterior distribution hits the edge of the prior. Conversely, for the slope α the recovered values are close to the median of the prior (α = 5), which shows that the parameter is poorly constrained. This is not surprising, as the two-pointing configuration only probes a limited range of independent scales (from ∼2 to ∼7 arcmin) and the available local measurements are affected by cosmic variance due to the limited azimuthal coverage. Additional pointings covering a wider range of scales are needed to set constraints on the power spectrum slope. The recovery test for the injection scale k_inj shows that the two-pointing configuration is sensitive to this parameter, as the distribution of the recovered values lies close to the one-to-one relation, although the measurement is not as accurate as for the Mach number. The injection scale can be estimated with a modest precision of 0.2−0.3 dex.

Fig. 2. Coverage tests for the two-pointing configuration and the model defined in Eq. (4). All the panels show the median of the inferred posterior distribution as a function of the true input value. The panels show the recovery tests for the Mach number ℳ3D = σv/cs (left), the slope of the turbulent cascade α (middle), and the logarithm of the injection wave number kinj, in units of kpc−1 (right). The gray dots show the results of individual simulations, whereas the color code indicates the point density inferred using a Gaussian kernel density estimator. The dotted black lines and the red dashed lines indicate the 1:1 relation and the running median of the data points, respectively.

We performed similar tests for the case where the fitted data points include the VSF and the velocity dispersions instead of the individual bulk velocities (see Molin et al. 2025, and Sect. 2.2). Fitting for the VSF yields noticeably poorer constraints. The posterior Mach number distributions are highly skewed, with a peak that is biased toward low values and a long tail toward high Mach numbers. This is due to the VSF measurement being sensitive only to relative motions (see Eq. 10) whereas the addition of the absolute bulk velocities with respect to the cluster rest frame yields much tighter constraints on the Mach number and the injection scale. Nonetheless, the VSF remains a very useful summary statistic to compare the best fitting model to the measurements, as it is less affected by the random nature of the process than the individual measurements.

3. Results

3.1. XRISM/Resolve observations of the Coma cluster

We applied the technique devised in Sect. 2 to the measurements of bulk velocities and velocity dispersions obtained by XRISM/Resolve for the Coma cluster (see X25). The position of the two XRISM pointings is highlighted in Fig. 3 and superimposed onto an XMM-Newton/EPIC map of the cluster in the [0.7−1.2] keV band. Each of the two pointings was subdivided into four approximately independent quadrants, which in the end provides us with a set of eight bulk velocity measurements. The retrieved bulk velocities with respect to the cluster rest frame z = 0.0233, set as the mean of the redshift distribution of ∼600 individual galaxy redshifts (Bilton & Pimbblet 2018), are shown in Fig. 3. Coma is known to contain two dominant galaxies, NGC 4874 and NGC 4889, that are separated by ∼200 kpc on the plane of the sky and exhibit a velocity difference of ∼700 km/s. This observation has been used as evidence that the cluster is currently in a post-merger phase. Using XMM-Newton/EPIC-pn, Sanders et al. (2020) found that the bulk velocities of the cluster show a bimodal behavior, with the central, southern and western regions aligning with the redshift of NGC 4889, whereas the remaining regions tend to align with NGC 4874. The Resolve measurements confirmed this picture with a much higher precision, showing that the gas in the core of the cluster exhibits considerable bulk velocities with respect to the cluster rest frame. The centroid of the lines was found to be blueshifted by 450 km/s and 730 km/s in the central and southern pointings, respectively. The gas velocities align with the redshift of NGC 4889, which identifies NGC 4889 as the dominant galaxy of the primary merging component. The velocity dispersion of the gas was found to be 210 km/s for the central pointing and 200 km/s for the southern region, which is substantially lower than the measured bulk velocities and possibly implies a steep velocity power spectrum. For more details on the analysis of the XRISM/Resolve data, we refer the reader to X25.

Fig. 3. Location of the XRISM/Resolve measurements superimposed on a XMM-Newton/EPIC surface brightness map of the Coma cluster. The white squares show the footprint of the two Resolve pointings split into four individual quadrants, with the number showing the bulk velocity of each region with respect to the cluster rest frame. The blue arrows show the positions of the two dominant galaxies NGC 4874 and NGC 4889 with their respective peculiar velocities. Figure reproduced with modifications from X25.

3.2. Power spectrum inference

We applied the model trained for the two-pointing configuration tuned to the Coma cluster observations and defined in Sect. 2.2 to the data vector made of the bulk velocity and velocity dispersion measurements reported in X25. In Fig. 4 we show a corner plot extracted from the inferred probability distributions for the three-parameter model (Eq. 4).

Fig. 4. Posterior distributions for the three-parameter model (Eq. 4) for the Coma cluster inferred from the XRISM/Resolve bulk velocity and velocity dispersion data.

The model retrieves a fairly high Mach number ${\mathcal{M}}_{3\mathrm{D},\mathrm{tot}}=0.{72}_{-0.22}^{+0.28}$, which is necessary to reproduce the large bulk velocities measured in the south pointing. If these large bulk velocities are interpreted as isotropic fluctuations, the measured 1D velocities in the southern field correspond to 3D velocities of ∼1250 km/s, i.e ∼80% of the sound speed in the medium. The model prefers a large injection scale ${\ell }_{\mathrm{inj}}=2.2_{-1.0}^{+2.0}$ Mpc. A large injection scale is required to match at the same time the large bulk velocities and the relatively low velocity dispersions. Reducing the injection scales to lower values (< 500 kpc) increases the velocity dispersion at fixed Mach number, as the power is shifted to smaller scales. The retrieved injection scale is very large compared to the projected distance between the two pointings (170 kpc). The measurement is rendered possible by our prior knowledge of the mean cluster redshift from the galaxies and the dependence of the velocity dispersion on the line-of-sight depth of the system.

The slope of the power spectrum is poorly constrained, as expected from the coverage test presented in Fig. 2. As already noted in X25, the model shows a preference for a steep slope, although the classical Kolmogorov slope is well within the 1σ contours. We note that the slope and the injection wave number are positively correlated. Indeed, models featuring a smaller injection scale but a steep slope also place most of the power on large scales, which is required to generate large bulk velocities without contributing to the velocity dispersion.

Given the weak constraints on the power spectrum slope, we also fitted the data with a simpler model whereby the slope α is fixed to the Kolmogorov value of 11/3, which is consistent with the slope of the density fluctuation power spectrum within the uncertainties (Zhuravleva et al. 2019). While the Mach number is left unchanged (${\mathcal{M}}_{3\mathrm{D},\mathrm{tot}}=0.{73}_{-0.20}^{+0.25}$), the injection scale is pulled to even larger values (${\ell }_{\mathrm{inj}}=6.6_{-2.2}^{+2.5}$ Mpc) and the posterior hits the upper boundary of the prior (10 Mpc). This is expected from the correlation between α and k_inj shown in Fig. 4, where we can see that fixing α = 11/3 enforces lower values of k_inj close to the lower edge of the prior.

Finally, to determine the constraints that can be obtained from XRISM data only, we ran the power spectrum reconstruction by fitting the VSF instead of the individual bulk velocities as discussed in Sect. 2.2. In this case, we retrieve statistically consistent results, although the posterior distribution of the Mach number (${\mathcal{M}}_{3\mathrm{D}}=0.{36}_{-0.12}^{+0.50}$) is found to be highly skewed, with a peak at a substantially lower value and an extended tail toward high values essentially corresponding to a lower limit. This result can be explained by the lack of constraints on the absolute bulk velocities, as the VSF only depends on relative bulk velocities. The results on the injection scale and slope are consistent with the values obtained by fitting the individual bulk velocity values and independently confirm the need for a large injection scale.

Fig. 5. Goodness-of-fit tests applied to the Coma cluster XRISM/Resolve data. The left and middle panels show the comparison between the Resolve data (red data points) and the best fitting model for the VSF (left) and the velocity dispersion measurements in the four quadrants of the central pointings and the full southern pointing (middle). The blue curves represent the median of 1000 simulations generated from the posterior distribution, whereas the blue shaded areas show the 16th to 84th percentiles of the generated mock datasets. The right-hand panel shows the cumulative distribution of test statistic value (Eq. 11) for 1000 simulations generated from the posterior, with the test statistic value obtained for the data indicated as the dotted vertical line. The dashed black line indicates the fraction of simulations with a better test statistic value than the one obtained for the data.

3.3. Goodness of fit

To check whether the model provides an appropriate description of the data, we extracted 1000 random parameter sets from the posterior distribution of our fiducial three-parameter model (see Fig. 4) and generated mock data vectors for each parameter set. We then studied the distribution of the mock datasets and attempted to compare the simulated data with the observations. However, as discussed in Sect. 2.2 the individual bulk velocities fluctuate considerably from one simulation to another, as the simulated random velocity fields can generate both positive and negative fluctuations, which makes it difficult to directly compare the distribution of simulated bulk velocities with the measured values. To alleviate this issue, we computed the VSF (Eq. 10) from each simulation and compared the distribution of mock VSF values to the observed VSF.

In the left-hand panel of Fig. 5 we show the VSF obtained by combining the eight individual bulk velocity measurements shown in Fig. 3 and compare the result with the median and dispersion of the mock VSF generated from the posterior parameter sets. We see that the median of the mock datasets lies below the observed values. However, the generated VSF values span a broad range of values which encompass the measured VSF values within 1σ. As noted in X25, this is due to cosmic variance, as the existing measurements only cover a small fraction of the cluster’s volume and the velocity field substantially fluctuates from one region to another. In the middle panel of Fig. 5 we show a similar comparison for the velocity dispersions averaged over the five regions considered (four quadrants of the central pointing and the full southern pointing). Here we can see that the data points lie very close to the median of the 1000 models. The model predicts that the velocity dispersions could vary in the range 100−400 km/s for similar power spectrum parameters. Therefore, if the power spectrum parameters retrieved here accurately represent the true power spectrum, the measured velocity dispersions correspond to the typical expected values, whereas the bulk velocities fluctuate more than the typical value, although such fluctuations are not unlikely given the large cosmic variance.

To test whether the observed deviation of the VSF from the median value matches the expectations of the model, following von Wietersheim-Kramsta et al. (2025) we attempted to build a test statistic that quantifies the expected level of deviations between the model and the data. Since the uncertainties in the fitted velocities are approximately Gaussian, a natural choice for the test statistic is the χ² distribution. However, in the case of the VSF we found that the distribution of generated values follows a log-normal distribution rather than a Gaussian distribution, such that we build a χ² distribution in log space instead of real space. For each of the data points in Fig. 5 (left and middle panels), we calculated the mean value of the mock datasets, ⟨σ_1D⟩ and ⟨lnVSF⟩. We then defined our test statistic for each realization as

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i\in {\sigma }_{1\mathrm{D}}}}\frac{{({\sigma }_{1\mathrm{D},i}-⟨{\sigma }_{1\mathrm{D},i}⟩)}^2}{\mathrm{\Delta }{\sigma }_{1\mathrm{D},i}^2}+{\displaystyle \sum _{i\in \mathrm{VSF}}}\frac{{(ln{\mathrm{VSF}}_i-⟨ln{\mathrm{VSF}}_i⟩)}^2}{\mathrm{\Delta }ln{\mathrm{VSF}}_i^2},\end{array}$$(11)

where the summation is done over the velocity dispersion (σ_1D) and VSF points, and the denominators represent the statistical uncertainties in the velocity dispersion and VSF, respectively.

We computed the value of our test statistic for the 1000 individual realizations of the posterior distribution and built the cumulative distribution of the χ² values. The resulting distribution is shown in the right-hand panel of Fig. 5. We performed the same calculation for the observed dataset and compared the retrieved test statistic value with the distribution of the 1000 random datasets. We found that the measured χ² value is lower than the mock χ² value for 58% of the datasets, such that the observed deviations are fully consistent with the expectations. Thus, we conclude that the model provides an adequate description of the data at hand. Obviously, this conclusion only holds for the currently available dataset, and additional constraints providing a wider spatial coverage or a higher spatial resolution may still require a more complex model.

4. Discussion

4.1. Power spectrum reconstruction method

The recovery tests shown in Sect. 2.3 demonstrate that the method proposed here allows to retrieve interesting constraints on the velocity fluctuation power spectrum in galaxy clusters, even in the presence of a coarse observational sampling. Our results show that in the case of the Coma cluster, two XRISM/Resolve pointings are sufficient to get a reliable estimate of the turbulent Mach number and interesting constraints on the injection scale, although a wider coverage is required to constrain the slope of the turbulent cascade. Compared to surface brightness fluctuation techniques (e.g., Zhuravleva et al. 2014; Dupourqué et al. 2023), our study directly traces the underlying velocity field without requiring assumptions on the link between surface brightness and velocity fluctuations. The use of simulations allows us to include a number of important observational effects that are difficult to implement analytically, such as projection effects, the emissivity distribution of the source, and the PSF of the instrument. Our code also allows us to implement the exact spatial footprint of the available observational datasets.

We stress that our method naturally takes cosmic variance into account, since the large training sample includes fluctuations in the local measurements induced by the varying large-scale fluctuation fields. At this point, cosmic variance is the largest source of uncertainty in the reconstruction of the velocity power spectrum because of the limited spatial coverage. Additional XRISM/Resolve pointings are needed to constrain the velocity field over a wider spatial range and decrease the cosmic variance, which is the current limiting factor for the reconstruction of the turbulent Mach number and the injection scale (see Fig. 5). For instance, ZuHone et al. (2016) originally proposed a five-pointing strategy with one pointing in the center and four equally spaced offset pointings located 200 kpc from one another. This strategy greatly reduces cosmic variance with respect to the current configuration and probes a wider range of scales, from 40 to 400 kpc, which allows us to break the degeneracy between injection scale and slope. Three additional pointings would be sufficient to distinguish unambiguously between a small (few hundred kpc) and a large (≳1 Mpc) injection scale.

Our method naturally takes into account a potential offset between the local gas bulk velocities and the rest frame of the system set by the mean redshift of the member galaxies. This is crucial in the case of Coma, where large bulk velocities between the gas and the galaxies (up to ∼800 km/s) are observed. Conversely, approaches relying exclusively on the velocity structure function (e.g., ZuHone et al. 2016; Molin et al. 2025) require a wider spatial coverage of the system to provide unbiased estimates of power spectrum parameters, as the VSF only depends on relative bulk velocities rather than the absolute ones. Correspondingly, our recovery tests show that with a limited number of available measurements the VSF only provides lower limits to the true turbulent Mach number. In X25, we circumvented this problem by including an additional VSF point to represent the velocity difference between the gas and the galaxies, although the scale had to be fixed at an arbitrary value of 1 Mpc. Future instruments such as the X-ray Integral Field Unit (X-IFU, Peille et al. 2025) on board NewAthena (Cruise et al. 2025) will allow us to measure gas velocities at higher spatial resolution and over wider areas, which will provide constraints on the Mach number independently of any external measurement. These measurements will also probe a much wider range of spatial scales, which will allow us to determine the slope of the turbulent cascade and potentially detect the dissipation scale.

4.2. Model power spectrum

Applying our technique to the set of two XRISM/Resolve pointings of the Coma cluster presented in X25 (see Sect. 3.2), we found that the data can be well described by an isotropic turbulent cascade model with a 3D Mach number of ∼0.7 and a large injection scale (> 1 Mpc). Considering the expected level of cosmic variance in a two-pointing configuration, the observed data are consistent with the expected level of fluctuations around the mean value (Fig. 5), such that the model provides an adequate description of the data at hand. In Table 1 we provide the best fitting parameters for the fiducial 3-parameter model as well as for several other configurations (see below). Here we discuss the implications of these results for the case in which the underlying assumptions of a single clump with an isotropic velocity field are born out.

In Fig. 6 we show the posterior envelope of the model power spectra converted into fluctuation amplitudes (Churazov et al. 2012),

$$\begin{array}{c}\hfill A_{3\mathrm{D}}(k)=\frac{v_{3\mathrm{D}}(k)}{c_s}=\sqrt{4\pi k^3{P}_{3\mathrm{D}}(k)}\end{array}$$(12)

Fig. 6. Amplitude of velocity fluctuations as a function of wave number for the Coma cluster. The curves are calculated from the posterior distribution of model parameters for a model with free power-law slope (blue), a model with Kolmogorov slope α = 11/3 (light green), and a model with Kolmogorov slope and ℓinj = 1 Mpc (orange). The solid curves and the shaded areas show the median and the 1σ error envelope of the fitted models. The dashed red curve shows the X25 model fitted to the VSF only, whereas the dotted green curve represents the fit to the VSF and velocity dispersion data from X25. For comparison, the gray shaded area shows the range of scales probed by the bulk velocity data, whereas the gray dashed line indicates the Kolmogorov slope.

with P_3D defined as a function of the 3D Mach number rather than the velocity dispersion (see Fig. 4). As discussed in X25, the large bulk velocities observed in Coma exceed the values measured in other clusters observed by XRISM (Perseus, A2029, Centaurus), where the bulk gas velocities reach a maximum of 300 km/s with respect to the mean of the galaxy redshifts (XRISM Collaboration 2025b,d). Our modeling implies that Coma is much more dynamically active than the aforementioned clusters, which are all classified as morphologically relaxed systems with a central cool core (Campitiello et al. 2022).

In Fig. 6 we compare our best-fitting models with the results presented in X25, where the VSF and the velocity dispersions were calculated directly from the inverse Fourier transform of the model power spectrum. The injection scale was fixed to 1 Mpc and an additional point at an arbitrary scale of 1 Mpc was added to the VSF to represent the velocity difference between the gas and the galaxies. Our results are statistically consistent with the X25 results, which lie within the 1σ range. To allow for a direct comparison between this method and the results of X25, we ran an additional model reconstruction with ℓ_inj fixed to 1 Mpc and Kolmogorov slope, leaving only the normalization as a free parameter. The retrieved power spectrum is in excellent agreement with the inference of the same model from X25. However, the goodness-of-fit test (Sect. 3.3) shows that the single-parameter provides a substantially worse description of the data, which motivates keeping the injection scale or the power-law slope as a free parameter. We also considered the model power spectrum for a configuration where the injection scale is fixed to 1 Mpc and the slope is left free to vary. In this case, we retrieve a very steep slope ($\alpha =7.{45}_{-1.25}^{+0.46}$), which hits the boundary of the prior, in agreement with the discussion in X25.

4.3. Non-thermal pressure support

For cases in which the velocity field is fully isotropic, the turbulent Mach number can be linked to the non-thermal pressure fraction in the system as (Eckert et al. 2019)

$$\begin{array}{c}\hfill \frac{P_{\mathrm{NT}}}{P_{\mathrm{tot}}}=\frac{{\mathcal{M}}_{3\mathrm{D}}^2}{{\mathcal{M}}_{3\mathrm{D}}^2+\frac{3}{\gamma }},\end{array}$$(13)

with γ = 5/3 the adiabatic index. It is important to note that, given the large injection scale and/or steep turbulent cascade, most of the power in our best fitting model lies in the very large scales (ℓ ≳ 2 Mpc), which rival or even exceed the virial radius of the system. These large-scale motions are not directly probed by our model and, if present, they do not contribute to the kinetic pressure within the halo. We recall that the turbulent Mach number inferred in our model is defined as the integral of the power spectrum (Eq. 4) over all scales, which includes a significant contribution from scales larger than 2 Mpc that are not directly probed in this study and may not even contain much gas. Therefore, it is worthwhile estimating the Mach number only within the scales that are relevant to the pressure balance within R₅₀₀ and/or R₂₀₀. From the mass-temperature relation of Umetsu et al. (2020) and a temperature of 8.37 keV (X25), we retrieve R₅₀₀ = 1400 kpc and R₂₀₀ = 2153 kpc. These estimates agree within 5% with the values obtained from the PlanckY_SZ measurement (Planck Collaboration X 2013). The effective turbulent Mach number within R₅₀₀ can be estimated as the integral of the power spectrum over all wavenumbers greater than k₅₀₀ = 1/R₅₀₀, i.e.

$$\begin{array}{cc}\hfill {\mathcal{M}}_{3\mathrm{D},500}& ={\left[{\displaystyle {\int }_{k_{500}}^{\infty }4\pi k^2{P}_{3\mathrm{D}}(k)\mathrm{d}k}\right]}^{1/2}\hfill \end{array}$$(14)

$$\begin{array}{cc}& ={\mathcal{M}}_{3\mathrm{D},\mathrm{tot}}{\left[\frac{{\int }_{k_{500}}^{\infty }{k}^2{(1+{(k/k_{\mathrm{inj}})}^2)}^{-\alpha /2}\mathrm{d}k}{{\int }_0^{\infty }{k}^2{(1+{(k/k_{\mathrm{inj}})}^2)}^{-\alpha /2}\mathrm{d}k}\right]}^{1/2}.\hfill \end{array}$$(15)

For our fiducial model, the Mach number integrated within R₅₀₀ decreases to ${\mathcal{M}}_{3\mathrm{D},500}=0.{45}_{-0.13}^{+0.18}$. The non-thermal pressure contribution within R₅₀₀ can be estimated from the posterior chains on the power spectrum parameters by calculating the range of values of ℳ_3D, 500 and inferring the 16th and 84th percentiles of the non-thermal pressure ratio distribution (Eq. 13), which yields

$$\begin{array}{c}\hfill \frac{P_{\mathrm{NT}}}{P_{\mathrm{tot}}}(<R_{500})=0.{10}_{-0.04}^{+0.08}.\end{array}$$(16)

This estimate considers only the scales that effectively contribute to balancing gravity within the halo. In Table 1 we report the best-fit parameters for the three models considered here (fiducial, Kolmogorov slope, fixed injection scale) as well as the turbulent Mach numbers estimated within R₅₀₀ and R₂₀₀. Within R₅₀₀, all the models retrieve a Mach number in the range 0.4 − 0.6, which exceeds the value estimated directly from the broadening of the emission lines by a factor of ∼2 (ℳ_3D = 0.24 from small-scale motions in X25). Jointly modeling the velocity dispersions and the bulk velocities allows us to probe the full range of spatial scales, which leads to a larger kinetic energy fraction compared to the values retrieved from the line-of-sight velocity dispersion. Indeed, the estimated 1D velocity dispersion is weighted by the local gas emissivity (Eq. 6), which varies substantially along the line of sight. Even in the case of the Coma cluster, which features a flat emissivity distribution within its core, 90% of the emission originates from the innermost 600 kpc, such that large-scale modes do not contribute to the measured line broadening. We stress, however, that the concept of non-thermal pressure is only meaningful for cases in which the cluster is essentially virialized and the acceleration term in the Euler equation is negligible (e.g., Lau et al. 2009; Vazza et al. 2018; Angelinelli et al. 2020), which may not hold in the early stages of a major merger.

4.4. Limitations of the model

In the standard turbulent cascade scenario, the injection scale is expected to be comparable to the size of the main perturbing body. The very large injection scale inferred by our model is comparable to the size of the cluster itself, which makes it unrealistic in cosmological structure formation scenarios (e.g., Vazza et al. 2012, 2017). The large injection scale inferred for the Coma cluster may represent a velocity field that includes a major merger between two cluster-scale entities with characteristic sizes greater than a Mpc. This picture is consistent with the simulations of Vazza et al. (2012), who showed that velocity power spectra in clusters do not show a clear cut-off at low wave numbers (and can be formally described by a high injection scale) and are dominated by ongoing mergers on large scales and turbulence on smaller scales (ℓ ≲ 0.3 R₂₀₀). Coma is known to be a dynamically active system, as it features a pair of concentric shock fronts located ∼1 Mpc east and west of the cluster core (Planck Collaboration X 2013; Churazov et al. 2021) and a giant radio halo (Brown & Rudnick 2011). The large bulk velocities measured in the two Resolve pointings imply that there is a substantial line-of-sight component to the velocity field, although the presence of the concentric shock fronts requires plane-of-the-sky motions as well. Alternatively, if the merger occurs preferentially on the plane of the sky, the gas may be compressed along the merger axis and form perpendicular outflowing gas “tongues”, as observed for instance in Abell 754 (Botteon et al. 2024). Comparing the observed Coma velocity field and thermodynamic properties with a grid of merger simulations with different mass ratios and impact parameters is necessary to test this scenario (e.g., ZuHone et al. 2018).

Obviously, our method relies on the assumption that the velocity field can be accurately described by a Gaussian random field. This assumption is not necessarily satisfied for cases in which the merger is in its early phase, where the velocity difference between the two merging entities is the primary source of large-scale bulk motions. If so, we could be observing the system before the turbulent cascade is fully developed and the power spectrum may be temporarily much steeper than Kolmogorov (see the discussion in X25). In this case, a two-component model including the bulk motions on large scales and the previously generated turbulent cascade on smaller scales may be more physically motivated than the single-component isotropic model considered here. Additional pointings covering a wider area are required to test this scenario.

5. Conclusion

In this paper, we proposed a method to retrieve the shape of the velocity fluctuation power spectrum in galaxy clusters based on measurements of X-ray line shifts and line broadening. We then applied our model to the XRISM/Resolve observations of the Coma cluster to study the shape of the velocity power spectrum in this system, estimate the non-thermal pressure fraction, and constrain the merger scenario. Our results can be summarized as follows:

• We showed that the velocity fluctuation power spectrum can be estimated by generating simulations of a Gaussian random field and forward modeling the measured bulk velocities and velocity dispersions. This approach allows us to include several important observational effects (e.g., projection effects, emissivity weighting, PSF smearing) and generate realistic mock datasets.

• Given a parametric form for the velocity power spectrum, our method generates a large number of simulations from random parameter sets and trains a neural network to learn the mapping between the model power spectrum parameters and the mock data vectors. This is achieved through the SNPE algorithm implemented in the sbi package (Tejero-Cantero et al. 2020). We distribute the code in the form of an easy-to-use Python package which can be applied to any observational configuration and emissivity distribution.

• We considered a two-pointing XRISM/Resolve configuration tailored to the available observations of the Coma cluster (X25). We tested the ability of our code to recover the input power spectrum parameters for a simple paramterization of the power spectrum considering a turbulent velocity field characterized by an injection scale and a turbulent cascade toward smaller scales. We found that two pointings are sufficient to accurately determine the turbulent Mach number and provide interesting estimates of the injection scale. However, a finer sampling is needed to constrain the slope of the turbulent cascade.

• We compared two different configurations for the provided data vectors. In the first case, we fit directly for the individual bulk velocities, whereas in the other case we compute the VSF from each generated simulation and compare the results with the measured VSF. We showed that the former configuration provides tighter constraints on the turbulent Mach number, as the absolute bulk velocities include additional information on the motion of the gas with respect to the cluster’s rest frame.

• Applying our method to the Resolve measurements, we found that a power spectrum model with ${\mathcal{M}}_{3\mathrm{D},\mathrm{tot}}=0.{72}_{-0.22}^{+0.28}$ and ${\ell }_{\mathrm{inj}}=2.2_{-1.0}^{+2.0}$ Mpc provides a good description of the available data. Specifically, a relatively high Mach number together with a large injection scale are required to reproduce the large bulk velocities but low velocity dispersions observed in the two Resolve pointings. Such large injection scales are unexpected in the usual cosmological structure formation, which implies that the system is likely in a merging phase with a substantial line-of-sight component.

• Our power spectrum model is in full agreement with the one derived in X25 from the same data but using a different method, based on fitting the velocity structure function. In comparison, our method allows for a natural inclusion of the external information on galaxy velocities.

• The inferred turbulent Mach number within R₅₀₀ estimated with our model (${\mathcal{M}}_{3\mathrm{D},500}=0.{45}_{-0.13}^{+0.18}$) is noticeably larger than the value estimated from the line broadening alone ($\sqrt{3}{\mathcal{M}}_z=0.24-0.25$). The difference can be explained by the joint modeling of line shifts and line broadening, which given the large bulk velocities requires a large injection scale and/or a steep turbulent cascade. Additional pointings covering a wider area are required to reduce the cosmic variance and confirm our results.

Acknowledgments

DE and MR acknowledge support from the Swiss National Science Foundation (SNSF) under grant agreement 200021_212576. Support for JAZ was provided by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. IZ acknowledges partial support from the Alfred P. Sloan Foundation through the Sloan Research Fellowship and from NASA grant 80NSSC18K1684. NT acknowledges the support by NASA under award number 80GSFC24M0006. NO acknowledges JSPS KAKENHI Grant Number JP25K07368. DRW acknowledges support from NASA grant 80NSSC23K0740.

References

Appendix A: Surface brightness distribution

The results presented here make important assumptions on the true underlying emissivity distribution. That is, we assume that the emissivity distribution is spherically symmetric and can be described by a single βmodel (Eq. 7). This assumption is supported by the surface brightness distribution of the gas in Coma, which appears relatively symmetric out to ∼30 arcmin with a large central core. In Fig. A.1 we show the surface brightness profile of Coma extracted from ROSAT/PSPC data because of the instrument’s large FoV (Briel et al. 1992). The profile was extracted from a mosaic made of three pointed PSPC observations covering the cluster out to ∼1.5 degrees (∼2.5 Mpc) using the Python package pyproffit (Eckert et al. 2020). While the βmodel is too simplistic to reproduce all the features observed in the data, it provides a good description of the brightness distribution, with deviations at the level of ≤10% out to ∼50 arcmin (see the bottom panel of Fig. A.1). In the regions where the XRISM/Resolve observations are obtained, the model predicts that 90% of the emission originates from a region of ∼600 kpc length.

While the spherical βmodel provides an adequate description of the 1D profile, the isophotes are elongated along the E-W axis, such that an elipsoidal distribution likely provides a better description of the 3D cluster shape. To study the impact of ellipticity on the results presented in this work, we extracted the brightness profile in elliptical annuli, with the shape of the elliptical annuli matching the ellipticity ratio of ∼1.2 between the major and minor projected axes observed in the X-ray isophotes. We then assumed that the third (line-of-sight) dimension has a similar shape as the projected major axis. We then used the corresponding 3D distribution to generate a new set of cluster simulations and train another model. With this configuration, we obtain ${\mathcal{M}}_{3D}=0.{80}_{-0.24}^{+0.20}$, $\alpha =6.0_{-1.7}^{+1.4}$, and $log{k}_{\mathit{inj}}=-3.{34}_{-0.34}^{+0.19}$, which is consistent with the results obtained from a spherical emissivity distribution within 1σ.

Fig. A.1. ROSAT/PSPC surface brightness profiles of the Coma cluster (black data points) modeled with a spherical βmodel (Eq. 7, blue curve). The bottom panel shows the residuals quantified as the ratio of data to model.

All Tables

All Figures

Fig. 1. Example of fluctuation field generation. Left: Model 3D power spectrum described by an injection at kinj and a classical Kolmogorov cascade on smaller scales (Eq. 4). The vertical dashed line shows the position of kinj = (500 kpc)−1. Right: Slice of the real-space velocity fluctuation field along the z axis generated from this power spectrum (Eq. 3) for a velocity dispersion σv = 300 km/s. The color code shows the velocity in km/s. The box size is 1.7 Mpc on a side.

In the text

Fig. 2. Coverage tests for the two-pointing configuration and the model defined in Eq. (4). All the panels show the median of the inferred posterior distribution as a function of the true input value. The panels show the recovery tests for the Mach number ℳ3D = σv/cs (left), the slope of the turbulent cascade α (middle), and the logarithm of the injection wave number kinj, in units of kpc−1 (right). The gray dots show the results of individual simulations, whereas the color code indicates the point density inferred using a Gaussian kernel density estimator. The dotted black lines and the red dashed lines indicate the 1:1 relation and the running median of the data points, respectively.

In the text

Fig. 3. Location of the XRISM/Resolve measurements superimposed on a XMM-Newton/EPIC surface brightness map of the Coma cluster. The white squares show the footprint of the two Resolve pointings split into four individual quadrants, with the number showing the bulk velocity of each region with respect to the cluster rest frame. The blue arrows show the positions of the two dominant galaxies NGC 4874 and NGC 4889 with their respective peculiar velocities. Figure reproduced with modifications from X25.

In the text

	Fig. 4. Posterior distributions for the three-parameter model (Eq. 4) for the Coma cluster inferred from the XRISM/Resolve bulk velocity and velocity dispersion data.
In the text

Fig. 5. Goodness-of-fit tests applied to the Coma cluster XRISM/Resolve data. The left and middle panels show the comparison between the Resolve data (red data points) and the best fitting model for the VSF (left) and the velocity dispersion measurements in the four quadrants of the central pointings and the full southern pointing (middle). The blue curves represent the median of 1000 simulations generated from the posterior distribution, whereas the blue shaded areas show the 16th to 84th percentiles of the generated mock datasets. The right-hand panel shows the cumulative distribution of test statistic value (Eq. 11) for 1000 simulations generated from the posterior, with the test statistic value obtained for the data indicated as the dotted vertical line. The dashed black line indicates the fraction of simulations with a better test statistic value than the one obtained for the data.

In the text

Fig. 6. Amplitude of velocity fluctuations as a function of wave number for the Coma cluster. The curves are calculated from the posterior distribution of model parameters for a model with free power-law slope (blue), a model with Kolmogorov slope α = 11/3 (light green), and a model with Kolmogorov slope and ℓinj = 1 Mpc (orange). The solid curves and the shaded areas show the median and the 1σ error envelope of the fitted models. The dashed red curve shows the X25 model fitted to the VSF only, whereas the dotted green curve represents the fit to the VSF and velocity dispersion data from X25. For comparison, the gray shaded area shows the range of scales probed by the bulk velocity data, whereas the gray dashed line indicates the Kolmogorov slope.

In the text

	Fig. A.1. ROSAT/PSPC surface brightness profiles of the Coma cluster (black data points) modeled with a spherical βmodel (Eq. 7, blue curve). The bottom panel shows the residuals quantified as the ratio of data to model.
In the text