© The Authors 2024

1. Introduction

Turbulence is one of the non-gravitational processes ruling the dynamics and thermodynamics of the ionized plasma found in groups and clusters of galaxies, contributing to its virialization and to the subtle balance between cooling and heating processes (see e.g., Simionescu et al. 2019; Kunz et al. 2022). Various sources can inject turbulence in the intragroup and intracluster medium (ICM hereafter).

One of the sources is the activity of the super-massive black holes (SMBHs) hosted at the center of galaxies within the cluster. These SMBHs eject large amount of energy via powerful jets. This process, known as AGN feedback, balances the gas cooling at the center of groups and clusters of galaxies, quenches star formation in galaxies, and displaces gas. This generates turbulent motions in massive halos out to scales of hundreds of kpc (see e.g., McNamara & Nulsen 2012; Eckert et al. 2021; Hlavacek-Larrondo et al. 2022).

A second source of turbulence is the accretion and merger of sub-halos, namely, massive galaxies and small groups, falling within the main halo potential well. Depending on their initial dynamics and the accretion and merger configuration, they tend to fall in a spiralling motion around the center of the main halo in sloshing motions that displace gas and induce turbulence (see e.g., Markevitch & Vikhlinin 2007; ZuHone & Su 2022).

A third notable source of turbulence is the growth of groups and clusters via the continuous accretion of dark matter, gas, and galaxies fed by the filamentary structures of the cosmic web to which they are connected. It also happens through merger events with other massive halos.

All these processes whether it is feedback, accretion, sloshing, generate shocks at all scales, where the released gravitational energy is transferred to the accreted gas mainly as kinetic energy, heating it to the temperatures of 10⁷ − 10⁸ K. However, a fraction is channeled as non-thermal energy in turbulent motions in the gas which cascade in eddies down to viscous scales, thereby contributing to the gas heating (see e.g., De Gasperin et al. 2017; Gaspari & Sdowski 2017; Voit 2018; Biffi et al. 2022). Numerical simulations of structure formation predict that the fraction of non-thermal pressure induced by turbulent motions in the total budget of the intragroup and intracluster gas is 10–30% (Lau et al. 2009; Vazza et al. 2009; Nelson et al. 2014; Shi et al. 2015; Biffi et al. 2016). This makes it a non negligible contribution to the overall gas pressure budget, which has to be accounted for, for instance, in estimation of the halo mass via the hydrostatic equilibrium (see e.g., Pratt et al. 2019).

As turbulent motions displace gas, they produce a velocity distribution of the gas particles within the halo referential, inducing a well-known Doppler-shift and broadening of spectral emission (Inogamov & Sunyaev 2003). Although emission lines from many chemical elements in the ICM have been observed for decades at X-ray wavelengths (see e.g., Mernier et al. 2018; Mernier & Biffi 2022), their proper characterization requires high spectral resolution. It is needed to disentangle the line broadening and centroid shift due to turbulent and gas motions from other contributors (e.g., thermal broadening, instrument response). To date, this has been achieved only with measurements from gratings instruments, such as RGS on board XMM-Newton (Pinto et al. 2020). However, the lack of spatial resolution and sensitivity of gratings mixes spatial scales and restricts measurement to the very core of the brightest clusters of galaxies. The only actual spectral high resolution measurements of lines centroid and broadening came from the SXS instrument on board the Hitomi satellite, in the core of the Perseus clusters (Hitomi Collaboration 2016, 2018).

Despite these instrumental limitations, indirect measurements of turbulence within the ICM have been obtained through various methods, such as the statistics of surface brightness fluctuations. These are assumed to trace the density fluctuations induced by turbulent gas motions (Churazov et al. 2012; Zhuravleva et al. 2014, 2018; Simionescu et al. 2019; Dupourqué et al. 2022). Another method relies on the inter-calibration of the Cu-Kα line between data from the instrumental background and the observations (Sanders et al. 2020; Gatuzz et al. 2022). As initially stressed by Churazov et al. (2012) the interpretation of density fluctuations as due to turbulent motions is one of several, as other processes could produce similar fluctuations (e.g., perturbations of the gravitational potential or entropy fluctuations due to infalling cold gas).

As a consequence, direct measurements (i.e., an observable directly linked to the speed of the medium) are needed to characterize the properties of the turbulence in groups and clusters of galaxies. The closest direct measurements available are high-resolution spectra at X-ray wavelengths, which provide the only way to properly measure the associated velocities. The next generation of X-ray observatories will embark integral field units with the necessary high spectral resolution of the order of a few eV (to be compared to the ∼150 eV resolution of the currently operating X-ray spectro-imager, that is EPIC on board XMM-Newton or ACIS on board Chandra). The Resolve instrument (Ishisaki et al. 2022) on board the XRISM mission (Ueda & Tashiro 2022) will soon provide measurements of the ICM emission line shifts and broadenings thanks to its ∼5 eV resolution over its 6 × 6 pixels matrix (over 3 × 3 arcmin²), albeit with a limited spatial resolution of ∼1 arcmin FWHM. In the mid-2030s, the European Space Agency flagship mission Athena (Barret 2022), with its X-ray Integral Field Unit instrument (X-IFU; Barret et al. 2023), will bring the spectral resolution to 2–3 eV and a spatial resolution of 5–10 arcsec. Associated with a large collecting area, providing sufficient sensitivity and a decent field of view (∼5 arcmin diameter), it will allow for a realistic and detailed line-of-sight-integrated velocity mapping of the ICM in groups and clusters of galaxies.

To fold such measurements into actual constraints on the properties of the turbulent motion and its associated power spectrum (e.g., velocity dispersion, the injection and dissipation scales) a statistical treatment of the relevant diagnostics, such as the line broadening or centroid shift, is required. Structure functions are related to the two point correlation functions of the projected velocity field underlining the 3D power spectrum of turbulent motions (Inogamov & Sunyaev 2003; Zhuravleva et al. 2012). They can be used to characterize of the turbulence in the ICM. ZuHone et al. (2016) used this diagnostic tool to demonstrate the capacity of Hitomi to map turbulence at the center of the Coma cluster. Following these authors, Roncarelli et al. (2018) applied this formalism to Athena/X-IFU simple mock observations giving a hint of how Athena/X-IFU will allow for breakthrough measurements of the velocity structure induced by bulk motions and turbulence.

Moreover, the treatment of the stochastic nature of the turbulent process raises the need to account for the limited number of observed realisation (i.e., one per cluster). This implies an associated sample (or cosmic) variance, that is the variance due to a finite number of available observations and samples, which is a key component to the error budget, although this is difficult to estimate. The aforementioned works partially accounted for this error component, making use of a Monte Carlo approach that involves the estimation of errors from a large number of time consuming mock observations.

In our first two papers (Clerc et al. 2019; Cucchetti et al. 2019), we have developed an analytical approach to provide relatively fast estimates of the structure function and its associated comprehensive budget of errors, including the statistical, systematics and sample variance for an ideal case of uniform and isotropic turbulence.

With this third instalment in the series, we aim to propagate this formalism down to the reconstruction of the turbulent power spectrum to quantitatively assess the ability to constrain its properties from X-ray measurement. In the present work, we focus on the measurement of centroid shifts of X-ray emission line from the ICM. We apply our method to the future Athena/X-IFU instrument in order to quantify the feasibility of this key scientific objective to the mission as a function of the instrument performances.

The paper is organized as follow: in Sect. 2, we introduce the basic formalism on the line shape and structure function. We present our method to mock X-ray observations as well as their post processing analysis in Sect. 3. In Sect. 4, we provide the details of the structure function model and complete error budget computation. We discuss the results of our end-to-end simulations and analysis down to quantitative constraints on the turbulent power spectrum in Sects. 5 and 6. We present our conclusions in Sect. 7.

The work presented here has been done with the current X-IFU instrumental configuration. However, a reformulation of the Athena mission has been led in an effort to meet budgetary requirements. A new instrumental configuration is yet to be confirmed and future work will need to be done to assess the impact of this new configuration on the work presented here. We assume a ΛCDM cosmology, with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7. With this setup, at a redshift of z = 0.1, 1 kilo-parsec (kpc) corresponds to an angular extend of 0.54 arcseconds (arcsec).

2. Structure function and its uncertainty

The two most common direct diagnostics available to study turbulence from X-ray spectroscopic observations, via gas motions along the line of sight (LoS), are the line centroid shift from its rest frame value and the line broadening (Cucchetti et al. 2019). In this work we focus on the centroid shift measurements over a given observed area. We assume that the statistics of this observable and, thus, of the underlying velocity, can be characterized using the second-order structure function, which is closely related to the power spectrum. In other words, we assume the turbulent field to be closely approximated by a Gaussian random field.

2.1. Second-order structure function

To characterize the gas turbulence as a function of different scales, we use the second order structure function SF(r). It is equivalent to a spatial two points correlation function and appears when observing the distribution of the velocity field over all the points in space. We assume that the velocity field, v, is isotropic and for any given component of the velocity field, we define a structure function such that:

$$\begin{array}{c}\hfill SF(\mathit{r})=\overline{{(v(\mathit{x}+\mathit{r})-v(\mathit{x}))}^2},\end{array}$$(1)

where the over-line denotes averaging over the different pairs of sky points x and x + r, separated by a distance, r = |r|.

In practice, only the LoS component of the velocity field is accessible from spectroscopy measurements, modulated by emissivity effects. We can, however, transpose the structure function definition to its 2D projected equivalent by using the centroid shift as a velocity proxy.

As developed in Clerc et al. (2019), in the case where line diagnostics are obtained from spectra measured over 2D bins (simply due to the pixel size or because of spatial binning over larger bins), the structure function can be written as :

$$\begin{array}{c}\hfill \mathrm{SF}(s)=\frac{1}{N_{\mathrm{p}}(s)}{\displaystyle \sum _s{|C_{{\mathcal{W}}^{\prime }}-C_{\mathcal{W}}|}^2,}\end{array}$$(2)

where s is the distance separating any pair of bins, 𝒲 and 𝒲′, while N_p is the number of pairs separated by the distance, s, and C_𝒲 is the observed centroid shift in velocity space measured over a given emission line or spectrum, obtained from the flux-weighted sum of all the single lines of sight contributing to a given region (pixels or bins).

The measurement of the structure function provides a view of the underlying turbulent velocity power spectrum, and can be used to estimate some characteristic parameters of the turbulence such as injection and dissipation scales, the slope of the turbulent cascade, or its normalization factor.

A key point however, is that the realization of a turbulent flow is naturally stochastic. Thus, a given turbulence regime, characterized by its power spectrum, can produce different velocity fields and thus different measured structure functions if the observed field of view is limited. This principle is called sample (or cosmic) variance. The structure function SF is a random variable, which depends on the particular realization of the velocity field. We are interested in its mean value and variance across multiple realizations, which are denoted as⟨SF⟩ and ${\sigma }_{\text{sf}}^2$, respectively.

In addition to sample variance, the limited exposure time of observation and uncertainties in the fit or in the calibration lead to statistical and systematic errors in X-ray observations of the centroid shifts. We use σ_stat to denote the standard deviation of the statistical error for C_𝒲. This uncertainty adds scatter to the velocity measurements on top of the sample variance. ZuHone et al. (2016) showed that any measurement of the structure function is systematically biased by the statistical error, provided that statistical errors are not correlated with the value of C_𝒲, such that:

$$\begin{array}{c}\hfill {⟨\mathrm{SF}⟩}^{\prime }=⟨\mathrm{SF}⟩+2{\sigma }_{\mathrm{stat}}^2.\end{array}$$(3)

One way to understand this result is that the structure function provides a measurement of the power of fluctuations at different scales, in the same way the power spectrum does. In that sense, a uniform measurement error can be seen as another source of fluctuations and adds bias to the SF as such. This formula holds, as long as all bins have independent and identically distributed measurement uncertainties.

2.2. Estimation of the full error budget

Cucchetti et al. (2019) showed that the variance of the structure function is also systematically biased by terms related to the statistical error on the centroid shifts measurements. They obtained:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{sf},\mathrm{tot}}^2={\sigma }_{\mathrm{sf}}^2+{\sigma }_{\mathrm{vf}}^2+{\sigma }_{\mathrm{sff}}^2+{\sigma }_{\mathrm{statf}}^2.\end{array}$$(4)

Let us detail the terms in the equation above. ${\sigma }_{\text{sf}}^2$ the intrinsic variance of the structure function. This is the residual variance in the absence of statistical error, namely, if σ_stat = 0.

The velocity field fluctuation term, ${\sigma }_{\text{vf}}^2=4{\sigma }_{\text{D}}^2{\sigma }_{\text{stat}}^2$, is related to the velocity fluctuations in the field of view. The larger the pixel-to-pixel variations in the velocities, the larger this term. This term is small for dissipation scales larger than the pixel or binned region size. Also, ${\sigma }_{\text{D}}^2$ is the variance of the first order structure function (as defined in Eq. (B.1)).

The structure function fluctuation term, written such that ${\sigma }_{\mathrm{sff}}^2=\frac{4}{N_{\mathrm{p}}(s)}\mathrm{SF}(s){\sigma }_{\mathrm{stat}}^2$, is related to the uncertainty with which the structure function is computed when the turbulent velocities are affected by statistical errors. This term is negligible if the number of pairs, N_p(s), used to estimate the structure function is large.

The statistical fluctuation term, which is written such that ${\sigma }_{\mathrm{statf}}^2=\frac{4}{N_{\mathrm{p}}(s)}(N_{\mathrm{nei}}(s)+1){\sigma }_{\mathrm{stat}}^4$, is purely a contribution of the statistics to the overall variance of the structure function. This term is also negligible if the number of pairs used to estimate the structure function is large, but otherwise it is most important at small spatial separations. Then, N_nei(s) is the number of neighboring regions at a distance s of any given point.

3. Mock X-ray observations

One of the key science objectives of Athena/X-IFU is to understand the role of turbulence in the gas dynamical assembly and relaxation in massive halos, namely, groups and clusters of galaxies. In order to assess this capability, we perform end-to-end simulations using the SIXTE software (Dauser et al. 2019). The following sections describe our procedure, as a continuation of the work of Cucchetti et al. (2019), with the difference that the present simulations are based on a physically motivated input cluster model to Athena/X-IFU mock observations.

3.1. Galaxy cluster model

3.1.1. General characteristics

Nearby clusters are targets of choice in order to grasp the whole range of spatial scales spanned by turbulent cascades. We modeled a nearby galaxy cluster with a cosmological redshift z = 0.1 and a scale radius of R₅₀₀ = 1300 kpc, corresponding to an angular radius of 11.7 arcmin (where R₅₀₀ is the radius encompassing a mean density of the gas 500 times that of the critical density of the Universe at the cluster redshift).

We assumed a non-cool core cluster model. Systems in this class are more likely to have undergone major events such as cluster mergers or collision and are more likely to show a higher level of bulk motions and turbulence. Our cluster gas model departs from the popular, albeit simplistic, β-model. We remain under the assumption of a spherically symmetric body. We only include instrumental background as an additional component in our simulations, neglecting astrophysical back- and foreground emission as well as any second-order instrumental systematic.

3.1.2. Profile

The cluster density and temperature profiles are taken from the joint analysis of eight non-cool core clusters in the XMM-Newton Cluster Outskirts Project (X-COP, Ghirardini et al. 2019). This temperature profile follows the parametrisation given by:

$$\begin{array}{c}\hfill \frac{T(x)}{T_{500}}=T_0\frac{T_{\min }{T}_0+{(x{r}_{\mathrm{cool}})}^{a_{\mathrm{cool}}}}{1+{(x{r}_{\mathrm{cool}})}^{a_{\mathrm{cool}}}}\frac{1}{{(1+{(x{r}_t)}^2)}^{c2}},\end{array}$$(5)

where x = r/R₅₀₀, and {T₀, T_min, r_coola_cool, r_t, c/2} are six parameters, and T₅₀₀ follows a parametrization on the redshift, see Ghirardini et al. (2019). Their values are provided in Table 1.

The density profile follows:

$$\begin{array}{c}\hfill n_{\mathrm{e}}^2(x)=n_0^2\frac{{(x{r}_c)}^{-\alpha }}{{(1+{(x{r}_c)}^2)}^{3\beta -\alpha /2}}\frac{1}{{(1+{(x{r}_s)}^{\gamma })}^{ϵ\gamma }},\end{array}$$(6)

where x = r/R₅₀₀ and {γ, n₀, r_c, r_s, α, β, ϵ} are six parameters. Their values are provided in Table 1. We also set an arbitrary upper limit on the electron density value at 0.05 cm⁻³, in order to avoid divergence as we go closer to the center of the cluster.

For the metallicity profile, we followed the shape derived by Mernier et al. (2017), which was obtained from the study of a sample of 44 nearby cool-core clusters, groups, and ellipticals observed with XMM-Newton:

$$\begin{array}{c}\hfill \mathrm{Fe}(r)=0.21{(r+0.021)}^{-0.48}-6.54\times exp(-\frac{{(r+0.0816)}^2}{0.0027}),\end{array}$$(7)

where r is the distance from the center of the cluster in kiloparsecs. The solar abundances are fixed to those of Anders & Grevesse (1989).

We assume an emission spectrum for the intra-cluster gas described by the APEC model (Smith et al. 2001) implemented under XSPEC (Arnaud 1996). It assumes an ionized, collisional, diffuse, and optically thin plasma. As we assume a fully ionized gas with ∼10% Helium in atom number (and weakly contributing heavier elements), electron density, n_e and hydrogen density, n_H, are related by n_e ∼ 1.2 n_H. We imposed the emissivity of the gas to vanish at radii beyond 5R₅₀₀.

3.1.3. Turbulent velocity

We assumed that the turbulent field is homogeneous and isotropic. Following ZuHone et al. (2016), we modeled the 3D velocity power spectrum with a standard functional, with k_inj and k_diss as the cut-off frequencies of the Kolmogorov cascade at low and high wave numbers, $\underline{\mathit{k}}$ (with k = 1/r the conversion between scales and wave numbers), α as the slope of the spectrum, and σ as the 3D velocity dispersion (units km s⁻¹). With $k=|\underline{\mathit{k}}|$, our model can be expressed as:

$$\begin{array}{c}\hfill P_{3\mathrm{D}}(\underline{\mathit{k}})={\sigma }^2\frac{{\displaystyle k^{-\alpha }{e}^{-{(k_{\mathrm{inj}}/k)}^2}{e}^{-{(k/k_{\mathrm{diss}})}^2}}}{{\displaystyle \int 4\pi k^2{k}^{-\alpha }{e}^{-{(k_{\mathrm{inj}}/k)}^2}{e}^{-{(k/k_{\mathrm{diss}})}^2}\mathrm{d}k}}.\end{array}$$(8)

In our simulations, the injection scale is set to 300 kpc. While it has not been determined by observations, we expect similar values if turbulence is driven by mergers or accretion from the large-scale structure (as found in Dupourqué et al. 2022). The dissipation scale is set to 10 kpc. The exact range of it can depend on the details of the ICM dynamics and properties and can vary between 1 and 100 kpc (Cavaliere et al. 2011). It was chosen to be 10 kpc arbitrarily with the perspective of testing X-IFU capabilities in this regime. The slope is α = −11/3, corresponding to the usual index of the Kolmogorov cascade in the energy spectrum. The normalization constant, σ, is set such that the expected velocity broadening of the turbulent field at the center of the cluster equals (ℳc_s)². Here, ℳ is the Mach number, set to 0.3 following ZuHone et al. (2016), and c_s is the speed of sound, set to 833 km s⁻¹. This gives σ ≃ 250 km s⁻¹.

We generated a realization of the line of-sight turbulent velocity field to match such a power spectrum by drawing a random field in the Fourier domain, V_x(k) = ve^iϕ. To do so, we used a Rayleigh distribution for the amplitude, v, and a uniform distribution for the phase, ϕ. The inverse Fourier transform of V_x(k) gives the velocity field (Cucchetti et al. 2019). The turbulent motions of the gas are then included in our cluster model by converting the LoS component of the velocity field into a redshift correction applied to the emitted spectra.

3.2. Post-processing

3.2.1. Mosaics

To fully map our toy model cluster at z = 0.1 out to R₅₀₀, a coverage of 19 separate pointings arranged around the center of the cluster is needed. We set the initial exposure time of each pointing to 500 ks. Each pointing mock observation is simulated separately. This relatively high exposure time is not meant to be performed on 19 different exposures during the mission. From this main simulation, we are able to pick the number of pointings and their individual exposure time by downsampling the observations to the desired length and by assembling the final event list suited to a given observing strategy. From the list of events from the SIXTE simulation, we reconstructed raw brightness maps in counts (later referred to as count maps).

3.2.2. Spatial binning

In order to reconstruct the line-of-sight-projected velocity field, we binned the 2D counts distribution of our mosaic mock observations into bins of multiple pixels. The binning was performed over multiple pointings of the same exposure times. We aim to define regions of equal signal-to-noise ratios (S/Ns) to extract individual centroid shift measurements with similar levels of statistical quality.

In practice, this binning procedure operates on the assembly of count maps with similar exposure time, dominated by the ICM Bremsstrahlung emission, and the bins (denoted as 𝒲) are created following a Voronoi tessellation (Cappellari 2009) with a constant S/N of 200.

3.2.3. Spectral fitting

For each of the Voronoi bins, we extracted a spectrum generated from the list of events. Each spectra is fitted under XSPEC, using the APEC model. More precisely, the spectra were fitted to an absorbed single-temperature thermal plasma model, including a velocity broadening component of the lines to account for turbulent velocities and accounting for the absorption along the LoS (i.e., phabs*bapec). The rebinning of the spectra was done before the fitting to increase the computation speed, following the optimal binning proposed by Kaastra & Bleeker (2016).

The column density, n_abs, of the phabs model is kept fixed to the value of 0.03 × 10²² cm⁻². The free parameters are the temperature, the abundance, the normalization, the redshift, and the line velocity broadening. All the parameters are fit simultaneously, over the full energy band (0.2–12 keV), and relying on C-statistics for the goodness of fit.

In order to account for the vignetting effects in our processing, a specific ARF is created for each bin, following the method by Cucchetti et al. (2018b). This is done by multiplying the response of each pixel within the bin by the vignetting function implemented in SIXTE, over all energies, and averaging, in each energy channel, the response of each pixel weighted by their respective number of counts.

3.2.4. Output and input parameter maps

After fitting, we used the best fit values and errors inside each bin to create output maps for each parameter and its associated error. Each spatial bin is associated to a value of the centroid shift, C, thus creating a spatial map of this parameter. The centroid shift, C, (tracing the turbulent motions), however, needs to be converted from the best fit value of the spectral redshift, z_fit, to velocity domain with :

$$\begin{array}{cc}& z_v=\frac{z_{\mathrm{fit}}-z}{1+z},\hfill \end{array}$$(9)

$$\begin{array}{cc}& C=\frac{{(1+z_v)}^2-1}{{(1+z_v)}^2+1}{c}_{\mathrm{light}},\hfill \end{array}$$(10)

where z_v is the redshift equivalent to the integrated LoS velocity induced by the injected turbulent motion in our simulations, z the cosmological redshift of the cluster, and c_light the speed of light. Equation (10) can be approximated by C = z_vc_light as the turbulent velocities are non relativistic. Our procedure relies on the assumption that the cosmological redshift of the cluster is perfectly known.

A map of the true centroid shift is created and serves as a reference to assess the goodness of fit of our spectral analysis. We obtain this map by projecting the velocity shift along the LoS, weighted by the emissivity of the ICM (see Sect. 5.4 in Cucchetti et al. 2018b).

3.2.5. Measured structure function

The second-order structure function was computed from the output map of measured centroid shifts, C, following Eq. (2). The associated statistical errors arising from the propagation of the best-fit parameter uncertainties were evaluated with the difference between the input and output bulk motion maps. This statistical error is then saved for subsequent application in the computation of the structure function variance.

4. Model and fit of the structure function

In order to model the measured second order SF, we applied and extended the formalism developed by Clerc et al. (2019) and Cucchetti et al. (2019) to our toy model. This formalism provides a semi-analytical way to estimate the full error budget over the structure function without having to reproduce the full simulations process a large number of times to estimate the SF mean and scatter (including the sample variance) at each step of a Monte Carlo fitting procedure.

4.1. Estimation of the SF

From the formulation of the structure function in Eq. (2), we demonstrate in Appendix A that for any emissivity distribution of the sources, projected over a field of view grouped into bins of any shape and at any position, the expected structure function can be given by:

$$\begin{array}{cc}\hfill ⟨\mathrm{SF}(s)⟩& ={\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(\underline{\mathit{k}})}\hfill \\ \hfill & \times \left(\frac{1}{N_{\mathrm{p}}(s)}{\displaystyle \sum _{d(\mathcal{W},{\mathcal{W}}^{\prime })=s}}{|\frac{c_{ϵ·\mathcal{W}}(\underline{\mathit{k}})}{F_{\mathcal{W}}}-\frac{c_{ϵ·{\mathcal{W}}^{\prime }}(\underline{\mathit{k}})}{F_{{\mathcal{W}}^{\prime }}}|}^2\right).\hfill \end{array}$$(11)

In the equation above, the second sum runs over all pairs, (𝒲, 𝒲′),  separated by a distance, s, which is fixed by the geometry of the bins. c_{ϵ ⋅ 𝒲} is the Fourier coefficient of the product, ϵ(x, θ) 𝒲(θ), the emissivity at a location, x, along a LoS labeled as θ, multiplied by the window function of the bin, which equals one inside the bin and zero outside. Then, θ is the 2D vector on the sky plane, designating a unique LoS. F_𝒲 is the total flux in the bin denoted as 𝒲.

In order to speed up the computation, we calculated the term in parenthesis in Eq. (11) only once for the given geometry and binning configuration. It was then saved on the disk and then multiplied by P_3D each time a new turbulent model was generated, thereby providing the ⟨SF⟩.

For a general emissivity model, it is necessary to compute the 3D Fourier transform of ϵ(x, θ) 𝒲(θ) for each bin 𝒲. In order to accelerate this step, we rely on the following approximation (see details in Appendix A): we assume that both the flux and the emissivity are independent from the LoS chosen within the bin. We can then choose any of those LoSs and define it as θ_eff, representative of the bin 𝒲 (but changing from bin to bin) such that:

$$\begin{array}{c}\hfill \forall (\mathit{\theta }\in \mathcal{W}):ϵ(x,\mathit{\theta })\approx ϵ(x,{\mathit{\theta }}_{\mathrm{eff}})\equiv {\epsilon }_{{\mathit{\theta }}_{\mathrm{eff}}}(x).\end{array}$$(12)

Hence, the simplified and faster computation of c_{ϵ ⋅ 𝒲} through:

$$\begin{array}{c}\hfill c_{ϵ·\mathcal{W}}\approx {\mathrm{FT}}_{1\mathrm{D}}({\epsilon }_{{\mathit{\theta }}_{\mathrm{eff}}}(x)){\mathrm{FT}}_{2\mathrm{D}}(\mathcal{W}(\mathit{\theta })).\end{array}$$(13)

This methodology allows for a fast computation of the structure function for any emissivity model and any observation and geometry. We implemented it to compute the modeled structure function in our MCMC procedure to model the structure function.

4.2. Estimation of the errors

The mean structure function, the sample variance, σ_sf, the velocity field fluctuations, σ_D, can be estimated numerically for any emissivity model, with a few simplifying assumptions on the detector and observation geometry. The other terms needed to assess the full error budget are themselves intrinsically related to the observational setup. These terms are N_p, N_nei, and σ_stat. We note that it is simple to derive N_p and N_nei, as they are related to the binning and geometry of the detector and can be obtained analytically with high accuracy.

We assumed a prior knowledge of σ_stat, which we obtained through a comparison of the output maps with the input maps. For a real observation, this quantity would be hard to estimate and would probably require a complete end-to-end simulation to get a proper estimate. Figure 1 shows the distribution of the measurement errors on the centroid shift, for two observation configurations: the “uniform exposure” and “mix exposure”, which are explicited later. The distributions are close to Gaussian, which justifies the correction explicited in Eq. (3).

Fig. 1. Histograms of the centroid shift measurement error calculated for the mosaic of 19 pointings for both mixed and uniform exposures. The plain lines show the Gaussian profile plotted with a standard deviation corresponding to that of the sample of the errors, that is, 34.1 km s−1 for the “mix exposure” and 34.5 km s−1 for “uniform exposure”.

4.2.1. Sample variance σ_sf

Computing the analytic sample variance involves integrals over the 3D Fourier space that make the full computation numerically untractable, at least in the context of multiple evaluations as in an MCMC procedure. However, Clerc et al. (2019) show a simplifying assumption that allows for a fast computation of the sample variance term. It requires the assumption of constant emissivity over the sky, as in Schuecker et al. (2004), for instance, with a so-called Xbeta model for the emissivity (defined in Appendix C). This allows for a decorrelation of the emissivity field and the Fourier transform of the emissivity-weighted centroid shift along the LoS. Under this approximation, and assuming a regular grid of fixed size bins, the variance of the structure function is expressed as:

$$\begin{array}{c}\hfill \mathrm{Var}[\mathrm{SF}(s)]=16\pi {\displaystyle {\int }_0^{\infty }\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{P}_l^2(\xi )P_{2\mathrm{D}}^{\infty 2}(\xi ){(1-J_0(2\pi \xi s))}^2\xi \mathrm{d}\xi }\end{array}$$(14)

with P_l the power spectrum of the bins of characteristic size, l (square of the Fourier transform of the window function), $P_{\text{2D}}^{\infty }$ as the 2D power spectrum of the centroid shift, extrapolated for an infinite analysis region, 𝒮_𝒜 as the area of the observed region, J₀ as the 0-th order Bessel function. One expression for $P_{\text{2D}}^{\infty }$ is :

$$\begin{array}{c}\hfill P_{2\mathrm{D}}^{\infty }=2{\displaystyle {\int }_0^{\infty }{P}_{3\mathrm{D}}(\sqrt{k_x^2+{\xi }^2}){(1+2\pi k_x{R}_c)}^2{e}^{-4\pi k_x{R}_c}\mathrm{d}{k}_x,}\end{array}$$(15)

with P_3D as described in Eq. (8) and R_c the core radius of the Xbeta model. The integrals are computed using the double exponential quadrature, allowing for a fast computation of each integral with a few hundreds of integrand evaluations (Takahasi & Mori 1974).

Because the observed velocity field is not collected over a regular grid, but in spatial bins of varying sizes and shapes, we made another assumption. We assumed that to each separation, s, we can associate an equivalent bin geometry and Xbeta core radius to be the input of Eq. (14). This equivalent geometry is obtained by identifying the bins involved in the computation of the structure function at the separation defined by s and assuming that they can be approximated by a grid of identical pixels of fixed size. Their size is then computed with the mean of the sizes of the pixel used in the separation, s. The Xbeta core radius is found by fitting the structure function provided for the Xbeta model in Clerc et al. (2019), to the structure function computed for the full toy model, as described in Sect. 4.1 at the true parameters. In that sense, we adapted the Xbeta model to provide the best possible estimation of structure function (and, hence, the variance) around the true value of the turbulent parameters. Finally, the area, 𝒮_𝒜, associated to each separation, s, is taken as the area covered by all bins entering the calculation of SF at this separation. This approach assumes that there is no covariance between the different scales of the structure function. Implementing this effect in our current framework would require us to estimate the covariance over the structure function points, each being estimated on a different grid of separations. This issue is beyond the scope of this study.

To gain insights into the actual dependence of sample variance errors on the various model parameters, we provide in Appendix D an approximate derivation of the fractional uncertainty σ_sf/⟨SF⟩. This shows in particular that relative uncertainties are inversely proportional to the largest scale spanned by the observation and roughly proportional to the injection scale in the regime of small spatial bin sizes. The binning size is also relevant, as larger bins will wash out the signal and tend to increase the relative uncertainties. Exact derivations are intractable analytically and would require a numerical integration.

4.2.2. Velocity field fluctuation term, σ_D

The computation of this term calls for simplifying assumptions to be adopted as well. The detailed computation is explicited in Appendix B. We can summarize it in the following way: σ_D is computed once for a given observation geometry (i.e., binning) for the true turbulent parameters. It is then saved and later approximated as a term proportional to the structure function when called in the computation of the errors in the MCMC. This is a crude approximation reflecting the similarity between the expressions for ⟨SF⟩ and ${\sigma }_{\text{D}}^2$, both depending linearly on the 3D power spectrum. This can be explicitly expressed in the following way:

$$\begin{array}{c}\hfill {{\sigma }_{\mathrm{D}}^2(s)|}_X={{\sigma }_{\mathrm{D}}^2(s)|}_{X_0}\frac{{⟨\mathrm{SF}(s)⟩|}_X}{{⟨\mathrm{SF}(s)⟩|}_{X_0}},\end{array}$$(16)

where we call X the set of turbulent parameters and X₀ the set of the true turbulent parameters.

5. Results

5.1. Outputs from cluster simulation

For the simulations of 19 pointings, two different mosaics have been produced: one where all the pointings have the same exposure time of 125 ks (for a total exposure of 2.3 Ms) and one where the exposure time was increased toward the outskirts of the cluster by having the central pointing at 125 ks, the first ring at 250 ks, and the second ring at 500 ks (meaning a total exposure of 7.6 Ms, see Fig. 2). The reasoning for the mix exposure version is an attempt to mitigate the decrease in count rate and, hence, an increase in bin sizes, as the emissivity decreases toward the outskirts. This allows us to obtain output maps with more homogeneous bin sizes (see Fig. 3 for the centroid shift).

Fig. 2. Illustration of the ‘mix exposure’ pointing strategy with exposure times in ks at the center of each pointing.

Fig. 3. Output maps for mosaics with a mix exposure (top) versus constant exposure at 125 ks per pointing (bottom), for the centroid shifts. A black circle indicates the R500 radius. All the spatial bins were constructed such as to provide identical signal to noise levels in the resulting spectra.

Such exposure times are obviously prohibitive and not realistic, so they would not be used for actual observations with the X-IFU. They have been used in order to provide a first configuration with low statistical errors to validate this simulation framework. Further tests, looking at optimizing the observing strategy for future X-IFU observations and, thus, with more realistic total exposure times, are planned in the near future. Furthermore, the hypothesis of isotropy of turbulence would allow for the use of only a quarter or a slice of our 19 pointings, offering a significant time reduction.

From the output and input maps of the centroid shifts, we generated the histogram of the absolute errors from individual regions (see Fig. 1). The standard deviation correspond to the σ_stat from Eq. (3), used for the computation of the expected structure function and error budget. For both exposure configurations, we find a statistical deviation of ∼30 km s⁻¹.

5.2. Simulated SF versus modeled

From the output maps of the measured centroid shifts, we extracted the corresponding structure functions. From these, we subtracted the appropriate bias correction dependent on σ_stat.

Figure 4 shows a comparison of structure functions from the mock observations (dashdotted lines) and from modeling (plain lines with error bars), for the uniform exposure configuration (in red) and mixed exposure (in blue). We stress that computation of the uncertainties are entirely part of the model. This explains why we prefer representing error bars on top of the model SF, rather than the measured SF. The overall shape and normalization of the modeled structure functions match with those from the mock observations. The central part of the SF, between 70 and 300 kpc separations, shows that the model underestimates this specific realization of the structure function. It is partly due to this specific realization of the turbulent velocity field, as from its stochastic nature another realization would show a different deviation; and partly to the hypotheses underlying the computation of the theoretical structure function. In the uniform exposure, where the bins are larger, the approximation made in Eqs. (12)–(13) is more likely to fail and to result in a stronger mismatch with respect to the data.

Fig. 4. Structure functions and residuals for our simulation of 19 X-IFU pointings. The blue curves show the results for the “uniform exposure” configuration and the red curves refer to those of the “mix exposure” configuration. The upper plot shows the structure functions, observed in the mock observations in dash-dotted lines, and the modeled in plain lines. The total error is over-plotted on the modeled structure function with error bars, as its computation is part of the modeling, as described in Sect. 4.2. The lower panel shows the value of χ = (data-model) and the errors.

We illustrate the different contributions of the error terms to the total error in Fig. 5. The sample variance dominates at all scales. The other terms have comparable contributions at lower scales. At higher scales, the statistical fluctuations become comparable with the sample variance. All of the error terms but the cosmic variance include a contribution of the measurement error σ_stat. This means that higher measurement errors would lead to these other terms having a larger contribution compared to the sample variance. We recall that, in principle, ${\sigma }_{\mathrm{stat}}\propto \frac{1}{\mathrm{S}/\mathrm{N}}\propto \frac{1}{\sqrt{t_{\exp }}}$. This is true if we keep the bin size constant. In our specific case however, the high spectral resolution of X-IFU requires a rather high (i.e., approx. 200) S/N to populate the numerous X-IFU spectral channel in order to characterize emissions lines. This forces us to keep the S/N constant and to adapt the binning accordingly. Then, even with lower exposure times, the σ_stat term should remain constant. As a result, the relative contributions of the error should remain relatively unchanged. The resulting structure function would span a different range of distances because of the change in bin size to accommodate for the lower exposure time. The overall error level should increase, as all the terms are dependent on the bin size, but their relative contributions should remain on the same order.

Fig. 5. Error terms contributing to the total error computation as a function of separation. The top panel represents the error models for the “mix exposure” computed with the input turbulent parameters. The lower panel represents these same errors computed for the “uniform exposure”. Four terms are represented, going from darker to lighter shades: the sample variance, the velocity field fluctuation, the structure function fluctuation, and the statistical fluctuation.

In Fig. 6 we show how the modeled structure function varies with the turbulent parameters. The range given to the parameters is roughly representative to the expected range in which they should be observed in the ICM (e.g., ZuHone et al. 2016; Simionescu et al. 2019; Vazza et al. 2011). The variations of the parameters have a relatively similar impact on the shape of the structure function. This lead us to expect strong degeneracies in the posterior distributions of the parameters. The dissipation scales ranging between 1 and 10 kpc have little to no effect on the structure function. This is due to the finite size of the bins, which, at this redshift, are 8.72 kpc in width. From this we can expect the dissipation scale to be hardly constrained in this observational configuration.

Fig. 6. Structure function models for our simulation of 19 X-IFU pointings. The red curves show results for the “uniform exposure” configuration.

5.3. Retrieving the turbulent parameter distribution

To obtain the constraints on the turbulent parameters of the ICM, we used the emcee python package to perform a Monte Carlo sampling of the posterior distribution of the parameters. We ran 32 walkers for 10 000 steps each. We checked proper convergence of the chains with the $\widehat{R}$ statistic (Vehtari et al. 2021), which is $1.01<\widehat{R}<1.06$ for all parameters. This is not the optimal value but was deemed acceptable with respect to the computation time needed for a greater number of samples. We used a Gaussian likelihood, with errors depending on the turbulent parameters, resulting in the following expression:

$$\begin{array}{cc}\hfill log(p({\mathrm{SF}}_{\mathrm{obs}}|X))& =-0.5{\displaystyle \sum _s[log(2\pi {\sigma }_{\mathrm{sf},\mathrm{tot}}^2)}\hfill \\ \hfill & +\frac{{(⟨\mathrm{SF}(s,X)⟩-{\mathrm{SF}}_{\mathrm{obs}}(s))}^2}{{\sigma }_{\mathrm{sf},\mathrm{tot}}^2}],\hfill \end{array}$$

with X as the set of turbulent parameters, ⟨SF(s, X)⟩ as the modeled SF described in Sect. 4.1 and ${\sigma }_{\text{sf},\text{tot}}^2$ as the modeled variance of the SF described in Sect. 4.2.

We ran the MCMC sampler for different configurations of priors and free parameters, for both the mix exposure and uniform exposure observations. First, we investigated the possibility of constraining all the parameters at once, including the dissipation. We used uniform priors for the following parameters, such that π(k_inj)∼𝒰(0, k_diss),π(α)∼𝒰(−10, 0),π(σ²)∼𝒰(0, 500²). Because k_diss can span the entire range [0;+∞[, we set the prior to the exponential distribution, such that π(k_diss) = λe^−λk_diss. This follows from the exponential distribution being the maximum entropy distribution on the [0;+∞[ semi-open interval (Park & Bera 2009). Moreover, setting a uniform prior would be equivalent to assuming a specific range of dissipation scale and we found that in doing so the prior bounds were systematically hit. The exponential distribution has also the advantage of giving equal probability to all decades, that is, the dissipation scale is equally likely to be found in the range [1; 10] kpc and in [10; 100] kpc.

The resulting posterior distributions are shown in Fig. 7. For clarity, we show the results in terms of lengths (L_inj, diss), rather than wavenumbers (k_inj, diss). The posterior distributions are plotted with the 1 and 2-σ confidence regions, and the marginalized distributions show the 1-σ confidence interval as a shaded region. The blue distribution represents the result obtained in the uniform exposure configuration and the red stands for the mixed exposure configuration. The input values are recovered within 1-σ in the distributions. However, they are not centered on the input values, displaying a slight bias in the recovery of the parameters, which we argue could be due to the sample variance affecting this particular realization. As expected from Fig. 6, the parameters are strongly degenerated. This is due to the nature of our choice of observable SF, as this range of covered scales does not allow us to fully break the degeneracy on the shape of the structure function between the parameters. The normalization parameter in particular is strongly underestimated. This can again be explained by the strong degeneracy between the parameters and the very large range allowed for the slope of the spectrum.

Fig. 7. Posterior distribution of the turbulent parameters obtained through the MCMC sampling. The blue contour show the distribution for the mixed exposure observation and the red for the uniform exposure. The inner contour delineates the 1σ level of the posterior probability, the second contour is at 2σ. Dotted lines show the input values of the simulations. Uniform priors were used for all parameters, except Ldiss, for which the prior is plotted in black.

One could argue that a uniform prior on the slope of the power spectrum, especially a wide one, is a naive expectation, and could call for a more informative prior, on the basis of simulations or previous observations. For instance, Zhuravleva et al. (2012) used a dimensional analysis to argue that −5 < α < −3. Figure 8 shows the result of the same sampling procedure, however, with π(α)∼𝒩(−11/3, 0.5). This is a translation of the assumption that it is unlikely to find a turbulent power spectrum with a slope too far from that of a Kolmogorov cascade. In such a scenario, the prior of the slope is very informative or, conversely, the data are very uninformative with respect to this prior. This is again a translation of the difficulty of constraining this parameter. This results mainly in a slightly tighter distribution for the injection scale, as shown by the marginalized uncertainties shown in Table 2. Adding a Gaussian prior to the slope of the power spectrum is only one option and one could equivalently do the same for the other parameters on the same assumptions of prior belief and/or knowledge.

Fig. 8. Similar to Fig. 7, except that a Gaussian prior is set on parameter α, instead of an uniform prior.

As the influence of the dissipation scale on the SF shape is very limited on the range of scales we investigate here, it is reasonable to set the dissipation scale to the simulation input value. The result of this is shown in Fig. 9. The Gaussian prior was kept on the slope parameter to restrict its range to the most reasonable values. The marginalized uncertainties of all our MCMC procedures are shown in Table 2. The resulting uncertainties are not significantly better than in a configuration where the dissipation scale is left free with an exponential prior. This is due to the prior being strongly constraining on the dissipation scale, and its aforementioned weak dependency on the SF.

Fig. 9. Posterior distribution obtained with setting the dissipation scale to the input value and letting the injection scale, norm and slope free. A Gaussian prior was set on the slope with π(α)∼𝒩(−11/3, 0.5), whereas the other parameters had uniform priors. The blue contour show the posterior distribution for the Mix exposure observation and the red for the uniform exposure. The prior is plotted in black. Dotted lines show the input values of the simulations.

6. Discussion

We have accounted for the total error budget in our analysis of mock observations of a galaxy cluster with the X-IFU instrument, including the sample variance. To do so, we have formulated simplifying assumptions in the computation of the errors on the structure function. As such, our total error budget computation is an approximation. It is likely that the full computation of the error would lead to varying results. However, it is unlikely that it would significantly alter the precision on the recovery of the turbulent parameters.

We purposely chose not to include any astrophysical background in our mock observations. The addition of astrophysical background should, in principle, only increase the uncertainty in the spectral fitting and hence on the recovered velocities. As such, it should translate into an increase of σ_stat, thus marginally affecting the nature of our results. The main impact should be on the width of the posterior distributions of the turbulent parameters which would be larger, as a result of the larger error. In addition, with our current exposure and signal-to-noise binning setup, the background is expected to be subdominant over all observed regions and should not influence the line velocity measurements significantly (Cucchetti et al. 2018a). The astrophysical background will be included in a forthcoming study investigating the optimal realistic observing strategy.

Our correction of the structure function and computation of the error assumes a perfect knowledge of the statistical error, σ_stat (see Sect. 3). Quantifying this term could in principle be done with bootstrapping from observations. The possible error on the estimation would have to be properly propagated to the final distribution of the parameters. Equivalently, we could make use of complete simulations to assess σ_stat given an observational setup.

Our toy model cluster is obviously not plainly representative of real clusters. In a complete analysis, modeling the structure function would have to take into account the departure to the spherical symmetry of all physical quantities governing the emissivity of the ICM. Same goes for the underlying turbulent field, that was assumed to be isotropic and homogeneous. This is not likely to be the case in real clusters, given the complex dynamics and substructures altering the turbulent velocity field, as well as the gravitational field breaking isotropy.

In addition to the previous comments, we should emphasize the other possible contributions to the velocity field within galaxy clusters. Indeed, bulk motion within the cluster, or gravitational redshift gradients (Molin et al. 2023) will be measured through the centroid shift as well and will bias the structure function and the subsequent turbulent parameters. Their effect will need careful handling and comparison with cosmological simulations.

The numerical implementation of the analytic model requires a number of simplifying assumptions to adapt the formalism to the non-ideal situation of a realistic cluster profile with irregular regions. This pushes the edge of using analytical treatment in the prediction of such observables. In a forthcoming study, we shall investigate methods such as simulation based inference (SBI, Tejero-Cantero et al. 2020), which requires a great number of simulations, to derive a full error budget without any approximations. It would also naturally solve the issue of the inclusion of the structure function’s covariance. This would allow us to include other observables, such as the line broadening, which was not included in the present study, as we deliberately focused on the velocity structure function. Such added observables would help break the degeneracy between the turbulent parameters.

7. Conclusions

With this study, we have further investigated the ability to characterize the turbulent motions in galaxy clusters with integral field spectroscopy in the X-ray domain. We developed a framework incorporating a comprehensive error budget and we applied it to mock observations (with unrealistic exposure time) of a case study cluster with the X-IFU instrument on board the Athena mission.

• Using end-to-end simulations (on the one hand) and modeling of structure functions for a physically sound cluster model (on the other), we performed an MCMC sampling with a semi-analytical model to the structure function in order to retrieve the parameters from a simple Kolmogorov-like turbulent model (injection scale, slope, and normalization factor i.e., velocity dispersion).

• The error model presented here relies on a careful account of all errors terms (i.e., the sample variance and additional statistical terms as defined in Cucchetti et al. 2019), and on the further development of the formalism to compute the structure function and errors from Clerc et al. (2019). Such a comprehensive account is mandatory in any analysis aiming at a proper physical characterisation of the ICM turbulence. As a drawback, it affects the final constraints on the turbulence parameters. Overcoming the remaining assumptions made for the error computations is the goal to reach in order to plainly assess the full potential of the future XIFU observations to characterize the physical properties of the ICM turbulence. This will constitute the basis for future work.

• Under the observational conditions of this paper, the constraints on turbulent parameters are strongly degenerate. Breaking such degeneracies needs either prior knowledge of some of the parameters or different strategies to access a wider range of spatial scales with sufficiently high resolution and S/N values. Solving this problem likely requires optimizing the strategy by exploring a relatively wide configuration space.

• In the most favorable observational configuration and assuming prior knowledge of the dissipation scale and power spectrum slope, we were able to constrain the injection scale to +16%/−13%, velocity dispersion to +10%/−7%, and slope to +10%/−12% of the input values.

In the upcoming work in this series of papers, we will look into the optimal observing strategy. We will define the best combination of X-IFU pointings, in terms of geometrical configuration and exposure time to extract the most stringent constraints on the turbulent parameters. We will also address the possibility of stacking multiple observed clusters, taking into consideration the stochastic nature of the turbulent process. Co-adding independent observations of clusters together should reduce the errors on each parameter by a factor of $\sqrt{N}$, with N the number of clusters (see e.g., Dupourqué et al. 2022).

Finally, we stress that this study was conducted with the Athena/X-IFU mission configuration before the mission reformulation by the European Space Agency¹. The newly issued specifications of the Athena mission affect top-level requirements, such as the effective area, spectral resolution, or X-instrument field-of-view. These changes will likely quantitatively impact the results presented here. Our future work will focus on and integrate these new instrumental specifications.

Acknowledgments

SB, AM, EP and NC acknowledge the support of CNRS/INSU and CNES. The following python packages have been used throughout this work : astropy (Astropy Collaboration 2013); astropy:2018; astropy:2022, chainconsumer (Hinton 2016), emcee (Foreman-Mackey et al. 2013), matplotlib (Hunter 2007) and cmasher (van der Velden 2020).

References

Appendix A: Modeling ⟨SF⟩ in arbitrary spatial binning configurations

The model for the expectation value ⟨SF⟩ of the structure function SF requires knowledge of the underlying turbulent velocity field and the emissivity field. Both these quantities are defined on a 3D spatial grid (x, y, z) = (x, θ). An additional important input is the geometrical design describing the position and shape of 2D bins within which X-ray spectra are extracted and from which a line centroid shift is measured. The formalism of Clerc et al. (2019) does not provide a closed formula for ⟨SF⟩ in the very general case of arbitrary bin shapes and arbitrary emissivity field. This section intends to extend their approach with this purpose in mind; we derive ⟨SF⟩ in a form that relies on pre-computed Fourier transform of the emissivity field.

We start from the definition of the structure function SF(s) computed from line centroid shifts C_𝒲, each associated to a spatial bin 𝒲:

$$\begin{array}{c}\hfill SF(s)=\frac{1}{N_{\mathrm{p}}(s)}{\displaystyle \sum _{d(\mathcal{W},{\mathcal{W}}^{\prime })=s}{|C_{{\mathcal{W}}^{\prime }}-C_{\mathcal{W}}|}^2.}\end{array}$$

In this formula, the sum runs over all pairs of bins separated by a distance s, the total number of them reads N_p(s). Given a 2D geometry, the list of pairs and of pairwise distances is known and it is a fixed input in the problem.

Developing this expression and averaging over many realizations of the turbulent velocity field, we are left with computing $⟨C_{\mathcal{W}}{C}_{\mathcal{W}\prime }⟩$ for any two bins 𝒲 and 𝒲′ (including the case of 𝒲 = 𝒲′). We write the line centroid shift in bin 𝒲 as a flux-weighted average of centroid shifts projected along the LoS (i.e., along coordinate x):

$$\begin{array}{cc}\hfill C_{\mathcal{W}}& =F_{\mathcal{W}}^{-1}{\displaystyle \int }\mathrm{d}\mathit{\theta }\mathcal{W}(\mathit{\theta }){\displaystyle \int }\mathrm{d}xϵ(x,\mathit{\theta })v(x,\mathit{\theta })\hfill \\ \hfill & =F_{\mathcal{W}}^{-1}\int \mathrm{d}\mathit{\theta }\mathcal{W}(\mathit{\theta })F(\mathit{\theta })C(\mathit{\theta }),\hfill \end{array}$$

where v(x, θ) stands for the LoS velocity, ϵ is the (3D) emissivity, F(θ) is the line flux along a single LoS labeled by the sky coordinate, θ, and F_𝒲 is the line flux summed within bin, 𝒲. Decomposing the velocity field into its Fourier components, $V_{\underline{\mathit{k}}}$, we obtain an expression for the velocity shift along a single LoS direction:

$$\begin{array}{c}\hfill C(\mathit{\theta })={\displaystyle \sum _{\underline{\mathit{k}}}{V}_{\underline{\mathit{k}}}{e}^{i\omega \mathit{\xi }·\mathit{\theta }}\int e^{i\omega k_{x}x}\rho (x,\mathit{\theta })\mathrm{d}x,}\end{array}$$(A.1)

where $\underline{\mathit{k}}=(k_x,\mathit{\xi })$ is the Fourier coordinate, ω = 2π/L, L is a large characteristic length and ρ(x, θ)≡ϵ(x, θ)/F(θ).

Using the definition of the velocity power-spectrum P_3D and writing complex conjugation with symbol ^*, we find:

$$\begin{array}{c}\hfill ⟨C_{\mathcal{W}}{C}_{{\mathcal{W}}^{\prime }}⟩=F_{\mathcal{W}}^{-1}{F}_{{\mathcal{W}}^{\prime }}^{-1}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3D}(\underline{\mathit{k}})c_{ϵ·\mathcal{W}}(\underline{\mathit{k}})c_{ϵ·{\mathcal{W}}^{\prime }}^{\ast }(\underline{\mathit{k}}).}\end{array}$$

The complex quantity $c_{ϵ·\mathcal{W}}(\underline{\mathit{k}})$ depends solely on the bin geometry, size, position and on the emissivity. Its value can be pre-calculated as the (3D) Fourier transform of the quantity (ϵ ⋅ 𝒲), namely:

$$\begin{array}{c}\hfill c_{ϵ·\mathcal{W}}(\underline{\mathit{k}})={\displaystyle \int }{e}^{i\omega \underline{\mathit{k}}·{\underline{\mathit{x}}}^{\mathbf{\prime }}}\mathcal{W}(\mathit{\theta })F(\mathit{\theta })\rho ({\underline{\mathit{x}}}^{\mathbf{\prime }})\mathrm{d}{\underline{\mathit{x}}}^{\mathbf{\prime }},\end{array}$$

where the function 𝒲(θ) takes value 1 within bin 𝒲 and zero outside of it. We may then reassemble the terms and use the fact that SF is a real quantity to obtain the formula shown in Eq. 11.

While varying the characteristics of the turbulent field, the term between parentheses in Eq. 11 remains constant. Therefore it can be computed only once, which significantly accelerates the exploration of the turbulent parameters through MCMC sampling. The main burden then resides in pre-computing c_{ϵ ⋅ 𝒲} for each bin 𝒲. In general², this calculation requires numerical evaluations of the (3D) Fourier transforms of the field (ϵ ⋅ 𝒲). These can be simplified and accelerated if the emissivity is constant within a bin 𝒲, regardless of the LoS θ, i.e. if ϵ(x, θ) = ε(x) and F(θ) = F. This is a relevant approximation if the bins are small in comparison to typical emissivity gradients. There, we obtain:

$$\begin{array}{cc}\hfill c_{ϵ·\mathcal{W}}(k_x,\mathit{\xi })& \simeq \stackrel{\sim }{\epsilon }(k_x)\times \widehat{\mathcal{W}}(\mathit{\xi }),\hfill \\ \hfill F_{\mathcal{W}}& \simeq F\times S_{\mathcal{W}.}\hfill \end{array}$$

In these expressions, S_𝒲 is the surface area of a bin, tilde and hat stand for 1D and 2D Fourier transforms respectively. In this paper, we make use of this approximation in order to make calculations tractable.

Appendix B: Computation of ${\sigma }_D^2=\text{Var}(D)$

Cucchetti et al. (2019) introduced the quantity D(s) whose variance intervenes in the calculation of uncertainties of the structure function. They defined:

$$\begin{array}{c}\hfill D(s)=\frac{1}{N^{\mathit{in}}(s)}{\displaystyle \sum _{\mathit{in}}[C(\mathit{\theta })-C({\mathit{\theta }}^{\prime })].}\end{array}$$(B.1)

Similarly as in Appendix A, C(θ) represents the line centroid shift along a LoS labeled by the sky coordinate θ. We arbitrarily set C(θ) = 0 for locations outside of the region of analysis. The sum runs over all pairs of coordinates (θ, θ′) separated on sky by a distance, s, both being located inside the region of analysis (hence, the label "in"). In the following, we assume the region of analysis is circular of radius, R. The number of such pairs is Nⁱⁿ and its formal derivation is provided in Appendix D of Clerc et al. (2019). The above definition for D is however ambiguous and it hints to a necessary revision of the Cucchetti et al. (2019) formula, which is out of the scope of the present study. We disambiguate the expression by requiring that the pairs under the sum sign are ordered such that θ always lies southwards of θ′. The following calculation introduces an extra set of pairs separated by a distance s, labeled "ext," which cross the border of the region of analysis (see Fig. B.1). We first compute the sum over all the pairs (i.e., both the "in" and "ext" pairs) and we then subtract an explicit derivation of the contribution from the "ext" pairs.

Fig. B.1. Geometrical quantities used in counting ‘ext’ pairs of length s from a point θ located within the circular region of analysis, 𝒜, of radius, R, at a distance θ = |θ| from the center. Whether pairs are pointing northward or southward, their number depends on the angle φ+ or φ− and we have φ = φ+ + φ−. Symmetrical relationships arise for the number of ‘inner’ pairs, with angles ϕ+ and ϕ−.

Let us introduce the new quantity D′ through:

$$\begin{array}{c}\hfill D^{\prime }(s)=\frac{1}{N^{\mathit{in}}(s)}{\displaystyle \sum _{in+ext}[C(\mathit{\theta })-C({\mathit{\theta }}^{\prime })],}\end{array}$$

in which the sum runs over all pairs (always assuming the conventional orientation toward the north). Rearranging the summation terms and using that C(θ) = 0 outside of the region of analysis, we find D′(s) = 0.

We are thus left with calculating the sum over all pairs of kind ‘ext’. Those pairs having their southern end θ within the analysis region contribute a term C(θ); while those having their southern end outside of the region of analysis contribute a term −C(θ′). To each direction, θ, within the field of analysis, 𝒜, we associate the number of "ext" pairs: φ_s(θ). This number can be decomposed into pairs that point northwards (${\phi }_s^{+}$) and pairs that point southwards (${\phi }_s^{-}$). We write:

$$\begin{array}{c}\hfill {\displaystyle \sum _{\mathit{ext}}[C(\mathit{\theta })-C({\mathit{\theta }}^{\prime })]=\sum _{\mathit{\theta }\in \mathcal{A}}({\phi }_s^{+}(\mathit{\theta })-{\phi }_s^{-}(\mathit{\theta }))C(\mathit{\theta })}\end{array}$$

Introducing the number of "in" pairs ϕ_s(θ) = 2π − φ_s(θ), it is clear that the following relations hold (e.g., see Fig. B.1):

$$\begin{array}{ccc}\hfill {\phi }_s={\phi }_s^{+}+{\phi }_s^{-}& ;\hfill & \hfill {\varphi }_s={\varphi }_s^{+}+{\varphi }_s^{-};\\ \hfill {\varphi }_s^{+}=\pi -{\phi }_s^{+}& ;\hfill & \hfill {\varphi }_s^{-}=\pi -{\phi }_s^{-}.\end{array}$$

An explicit expression for ϕ_s is given in Appendix D of Clerc et al. (2019). Hence, we only need to calculate ${\varphi }_s^{+}$ from which we will deduce the other quantities (${\phi }_s^{-}$, ${\phi }_s^{+}$, and ${\varphi }_s^{-}$).

Let us describe the position of the point θ by two coordinates: θ = |θ| its distance to the center of the circular shape of 𝒜 ; and ψ is the angle θ makes with the horizontal line. Geometrical considerations lead to distinguish five mutually exclusive cases:

$$\begin{array}{c}\hfill {\varphi }_s^{+}(\theta ,\psi )=\{\begin{array}{cc}\pi & \mathrm{if}{\tau }_s(\theta ,\psi )>\pi ,\hfill \\ {\varphi }_s(\theta )-\pi & \mathrm{if}{\tau }_s(\theta ,\psi )<{\varphi }_s(\theta )-\pi ,\hfill \\ 0& \mathrm{if}{\tau }_s(\theta ,\psi )<0,\hfill \\ {\varphi }_s(\theta )& \mathrm{if}{\tau }_s(\theta ,\psi )>{\varphi }_s(\theta ),\hfill \\ {\tau }_s(\theta ,\psi )& \mathrm{otherwise},\hfill \end{array}\end{array}$$

where τ_s(θ, ψ) = ϕ_s(θ)/2 − ψ and τ_s ∈ [ − π/2, 3π/2]. Figure B.2 shows one example (s = 0.4R) representation of the function $f_s(\mathit{\theta })={\phi }_s^{+}(\mathit{\theta })-{\phi }_s^{-}(\mathit{\theta })$. This purely geometrical function associates to each point within the field of analysis the balance between the number of "ext" pairs pointing northwards and those pointing southwards.

Fig. B.2. Representation of the functional fs(θ). It takes value in [ − π, π] for each point θ = (y, z) located within the circular region of analysis. This region of radius, R, is bounded by the plain black line.

Using that D′ = 0, we write:

$$\begin{array}{c}\hfill D(s)=-\frac{1}{N^{\mathit{in}}(s)}\sum _{\theta \in \mathcal{A}}{f}_s(\mathit{\theta })C(\mathit{\theta }).\end{array}$$

Since ⟨C⟩ = 0, we have that ⟨D⟩ = 0 and Var(D) = ⟨D²⟩. Developing the calculation involves the quantity ⟨C(θ₁)C(θ₂)⟩ and Eq. A.1 may be used to obtain:

$$\begin{array}{c}\hfill Var(D)=\frac{1}{{(N^{\mathit{in}})}^2}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3D}(k){\left|{\int }_{\mathit{\theta }\in \mathcal{A}}\mathrm{d}\mathit{\theta }{f}_s(\mathit{\theta })e^{i\omega \mathit{\xi }·\mathit{\theta }}\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)\right|}^2.}\end{array}$$

The function $\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)$ may be considered either as the (1D) Fourier transform of the function χ^θ(x) = ρ(x, θ); or, if seen as a function of θ, it is the (2D) inverse Fourier transform of the function $u_{k_x}(\mathit{\xi })=\stackrel{\sim }{\rho }(k_x,\mathit{\xi })$ (up to a constant factor). Hence, we obtain the following numerically tractable expression:

$$\begin{array}{c}\hfill Var[D(s)]=\frac{1}{N^{\mathit{in}}{(s)}^2}{\left(\frac{\omega }{2\pi }\right)}^4{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3D}(k){\left|(u_{k_x}\otimes \widehat{f_s})(\mathit{\xi })\right|}^2,}\end{array}$$(B.2)

involving a series of 2D convolutions indicated with the sign ⊗. We remind that a hat indicates 2D Fourier transform.

In practical applications, the computation of D(s) is performed in spatial bins of finite size. The derivation above assumes infinitely small bins. Similarly as in Clerc et al. (2019, in particula, Appendix E), we introduce an approximate scheme for dealing with spatial bins of finite size, ℓ, all supposedly of identical shape and regularly spaced. This approximation replaces C(θ) above by a filtered version:

$$\begin{array}{c}\hfill C(\mathit{\theta })\to {\displaystyle \int }\mathrm{d}\mathit{\mu }{\mathcal{F}}_{\ell }(\mathit{\mu })C(\mathit{\theta }-\mathit{\mu }),\end{array}$$

where ℱ_ℓ is a top-hat convolution filter describing the shape of bins. This approximation works under the assumption of sufficiently small bins with respect to the size of the analysis bin and to typical variations of the emissivity. The corrected equation for Var(D) is obtained by replacing $\widehat{f_s}$ by the product $\widehat{f_s}{\widehat{{\mathcal{F}}_{\ell }}}^{\ast }$ in Eq. B.2.

Appendix C: Plane-constant emissivity model

An useful simplified emissivity model consists of a β-model along the LoS direction x, identical for any selected LoS θ. We call this model Xbeta, it is given by:

$$\begin{array}{c}\hfill ϵ(x,\mathit{\theta })=\epsilon (x)=\epsilon (0){(1+\frac{x^2}{r_c^2+{\theta }_0^2})}^{-3\beta }.\end{array}$$(C.1)

This expression involves a core-radius, r_c, and a fixed impact parameter, θ₀. The characteristic size of the model along the LoS direction is $R_c={(r_c^2+{\theta }_0^2)}^{1/2}$. This emissivity model does not depend on the sky coordinate θ, hence, it produces uniform surface brightness. Analytical expressions for its Fourier transform make the computation of the 2D power spectrum of the centroid shift easily tractable (Churazov et al. 2012; ZuHone et al. 2016; Clerc et al. 2019).

Appendix D: Approximate scaling of sample variance

This section details the calculation leading to approximate scalings of the sample variance in the structure function. To this end, we restrict the following discussion to emissivity field described by Eq. C.1, regularly binned centroid shift maps and a very large region of analysis. Under these assumptions, we have (see Eq. 12 and 13 in Clerc et al. 2019):

$$\begin{array}{c}\hfill \mathrm{sf}(s)=4\pi {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle {\int }_0^{\infty }{P}_{\ell }{P}_{2D}^{\infty }(1-J_0(\omega \xi s))\xi \mathrm{d}\xi ,}\end{array}$$(D.1)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{sf}}^2(s)=\frac{16\pi }{{\mathcal{S}}_{\mathcal{A}}}{\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle {\int }_0^{\infty }{P}_{\ell }^2{P_{2D}^{\infty }}^2{(1-J_0(\omega \xi s))}^2\xi \mathrm{d}\xi ,}\end{array}$$(D.2)

where sf = ⟨SF⟩, ${\sigma }_{\text{sf}}^2$ = ⟨(SF − sf)²⟩. In these expressions, ω = 2π/L with L an arbitrary length setting the length unit system ; P_ℓ(ξ) is the power spectrum of the 2D bins, each of characteristic size, ℓ ; $P_{2D}^{\infty }(\xi )$ is the power spectrum of the centroid shift map, extrapolated to an unbounded region of analysis; 𝒮_𝒜 denotes the area of the region of analysis; J₀ is the Bessel function of the first kind of order 0. Furthermore, we specify the emissivity model following Eq. C.1, involving a characteristic size, R_c.

In what follows we will be interested in the relative sample variance, as given by σ_sf/sf. We formulate a number of additional assumptions in order to derive useful approximations:

• The 3D velocity power spectrum is a power-law on the interval bounded by the injection and dissipation scales and vanishes outside of this interval. We note this is a simplification of the actual power-spectrum used in our study (Eq. 8).

• All spatial bins are of identical size ℓ and the Fourier transform of their shape is well approximated by a top-hat function over the interval: −1/ℓ < ωξ < 1/ℓ.

• Bins are small with respect to the projected injection length and larger than the dissipation length of turbulent motions: k_inj < 1/ℓ < k_diss.

The expression for $P_{2D}^{\infty }$ provided in Eq. C.4 of Clerc et al. (2019) makes evident the role played by the emissivity profile as a low-pass filter, with characteristic cut-off frequency 1/R_c:

$$\begin{array}{c}\hfill P_{2D}^{\infty }(\xi )\propto {\displaystyle \int }{P}_{3D}\left(\sqrt{k_x^2+{\xi }^2}\right)P_{\epsilon }(k_x)\mathrm{d}{k}_x,\end{array}$$(D.3)

where P_ε represents the 1D power spectrum of the emissivity profile. Two distinct cases appear, according to R_c being larger (case A) or smaller (case B) than the turbulent injection scale.

D.1. Study of case A (1/R_c < k_inj)

In case A, the term under the integral in Eq. D.3 contributes non-zero terms only for values of k_x very close to zero. We then have $P_{2D}^{\infty }(\xi )\propto P_{3D}(\xi )$ and we approximate the result by writing:

$$\begin{array}{c}\hfill P_{2D}^{\infty }(\xi )\simeq \{\begin{array}{cc}\hfill A{\left(\frac{\omega \xi }{k_{\mathit{inj}}}\right)}^{\alpha }& \mathrm{if}\omega \xi \in [k_{\mathit{inj}};k_{\mathit{diss}}],\hfill \\ \hfill 0& \mathrm{otherwise}.\hfill \end{array}\end{array}$$

Here, A is a constant that cancels out in forming the ratio σ_sf/sf. The integrals involved in expressions D.1 and D.2 extend from 0 to +∞. The expressions for P_ℓ and P_2D show it is only necessary to compute the integrals over the domain [k_inj, 1/ℓ]. Two regimes are then worth investigating.

At small separations s ≪ ℓ, the second-order Taylor expansion of the Bessel function J₀ leads to:

$$\begin{array}{c}\hfill \frac{{\sigma }_{\mathrm{sf}}}{\mathrm{sf}}\simeq 2\sqrt{\pi }\frac{L_{\mathit{inj}}}{{\mathcal{S}}_{\mathcal{A}}^{1/2}}\frac{\mathcal{G}{(r,2\alpha +6)}^{1/2}}{\mathcal{G}(r,\alpha +4)}\mathrm{if}s\ll \ell .\end{array}$$(D.4)

Here, we define r = ℓk_inj ∈ [0, 1] and the specific function:

$$\begin{array}{c}\hfill \mathcal{G}(r,n)=\frac{1}{n}(r^{-n}-1).\end{array}$$

At large separations, s ≫ L_inj, assuming that the Bessel function vanishes toward large values of its argument:

$$\begin{array}{c}\hfill \frac{{\sigma }_{\mathrm{sf}}}{\mathrm{sf}}\simeq 2\sqrt{\pi }\frac{L_{\mathit{inj}}}{{\mathcal{S}}_{\mathcal{A}}^{1/2}}\frac{\mathcal{G}{(r,2\alpha +2)}^{1/2}}{\mathcal{G}(r,\alpha +2)}\mathrm{if}s\gg L_{\mathit{inj}}.\end{array}$$(D.5)

D.2. Study of case B (k_inj < 1/R_c < k_diss)

In case B, the integral in Eq. D.3 decomposes into a sum over three intervals: [0, k_inj], [k_inj, 1/c] and [1/R_c, +∞]. A close examination of the integrals for various values of ξ leads to the following approximation:

$$\begin{array}{c}\hfill P_{2D}^{\infty }(\xi )\simeq \{\begin{array}{cc}\hfill B& \mathrm{if}\omega \xi <k_{\mathit{inj}},\hfill \\ \hfill B{\left(\frac{\omega \xi }{k_{\mathit{inj}}}\right)}^{\alpha }& \mathrm{if}\omega \xi \in [k_{\mathit{inj}};k_{\mathit{diss}}],\hfill \\ \hfill 0& \mathrm{otherwise}.\hfill \end{array}\end{array}$$

Here, B is a constant of no relevance in the discussion. The expressions for P_ℓ and $P_{2D}^{\infty }$ show that integrals D.1 and D.2 must be calculated over domain [0, 1/ℓ]. Similarly as in case A, we identify two limiting regimes, allowing for a Taylor expansion of the J₀ Bessel function under the integrals. It follows that:

$$\begin{array}{c}\hfill \frac{{\sigma }_{\mathrm{sf}}}{\mathrm{sf}}\simeq 2\sqrt{\pi }\frac{L_{\mathit{inj}}}{{\mathcal{S}}_{\mathcal{A}}^{1/2}}\frac{{[\frac{1}{6}+\mathcal{G}(r,2\alpha +6)]}^{1/2}}{\frac{1}{4}+\mathcal{G}(r,\alpha +4)}\mathrm{if}s\ll \ell ,\end{array}$$(D.6)

$$\begin{array}{c}\hfill \frac{{\sigma }_{\mathrm{sf}}}{\mathrm{sf}}\simeq 2\sqrt{\pi }\frac{L_{\mathit{inj}}}{{\mathcal{S}}_{\mathcal{A}}^{1/2}}\frac{{[\frac{1}{2}+\mathcal{G}(r,2\alpha +2)]}^{1/2}}{\frac{1}{2}+\mathcal{G}(r,\alpha +2)}\mathrm{if}s\gg L_{\mathit{inj}}.\end{array}$$(D.7)

D.3. Interpretation

Under the assumptions stated in this Appendix, the relative uncertainties due to sample variance scale as provided in Eq. D.4, D.5, D.6, and D.7, each corresponding to a limiting case obtained by comparing together the values of R_c (typical size of the galaxy cluster), of L_inj, L_diss (typical scales of the turbulent field) and of ℓ (typical size of the spatial binning scheme).

In all four limiting cases, it clearly appears that relative uncertainties scale with the ratio $L_{\mathit{inj}}/\sqrt{{\mathcal{S}}_{\mathcal{A}}}$ at fixed r = ℓ/L_inj. Larger regions of analysis that sample multiple (projected) injection lengths are less prone to sample variance, as one might expect. Figure D.1 is obtained by fixing the bin size ℓ and letting only the injection scale vary. The relative uncertainty scales linearly, with L_inj in the regime of large separations, if L_inj is sufficiently larger than the bin size. The dependence with the bin size ℓ and the slope of turbulence, α, through the functional 𝒢 is less intuitive. Generally, such larger bin sizes are known to induce larger relative variance; whereas steeper power spectra also induce larger relative sample variance. Figure D.2 shows the dependency of relative uncertainties upon the bin size ℓ, all other quantities held fixed. In the regime of large separations, the relative uncertainties does not depend much on the bin size, provided bins are small enough.

Fig. D.1. Approximate relative uncertainty (sample variance only) in the two cases, A and B, and for the limiting regimes of small and large separations, s. In this figure, we fix the bin size, ℓ, the area, 𝒮𝒜, and the slope, α, letting only the injection scale vary (x-axis, in units of the bin size, ℓ).

Fig. D.2. Approximate relative uncertainty (sample variance only) in the two cases, A and B, and for the limiting regimes of small and large separations, s. In this figure we fix the injection scale, Linj, the area, 𝒮𝒜, and the slope, α, letting only the bin size, ℓ,  vary (x-axis, in units of Linj).

D.4. Exact calculations

The integrals in Eqs. D.1 and D.2 (and in Sect 4.2) for the sample variance computation can be computed exactly, namely, without approximations. This is the approach used in the core of this work. In Fig. D.3, we show the relative error obtained for different values of the turbulent parameters as a function of pair separation, s. We employ the same method as described in Sect. 4.2 to associate each value of s with a neighboring, best-matching Xbeta emissivity model, with its associated geometry (i.e., bin size and total area). As expected, the resulting relative error is independent of the norm of the power spectrum (bottom panel). Similarly as for the structure function, dissipation scales smaller than the bin size do not affect the relative error and even above this size, they have a limited influence. The injection scale has a significant impact on the relative error at large separations, with larger injection scales resulting in larger errors.

Fig. D.3. Exact relative error computation of the structure functions for the mix exposure observation, as a function of separation. Only the sample variance is included in the error budget.

All Tables

All Figures

	Fig. 1. Histograms of the centroid shift measurement error calculated for the mosaic of 19 pointings for both mixed and uniform exposures. The plain lines show the Gaussian profile plotted with a standard deviation corresponding to that of the sample of the errors, that is, 34.1 km s−1 for the “mix exposure” and 34.5 km s−1 for “uniform exposure”.
In the text

	Fig. 2. Illustration of the ‘mix exposure’ pointing strategy with exposure times in ks at the center of each pointing.
In the text

	Fig. 3. Output maps for mosaics with a mix exposure (top) versus constant exposure at 125 ks per pointing (bottom), for the centroid shifts. A black circle indicates the R500 radius. All the spatial bins were constructed such as to provide identical signal to noise levels in the resulting spectra.
In the text

Fig. 4. Structure functions and residuals for our simulation of 19 X-IFU pointings. The blue curves show the results for the “uniform exposure” configuration and the red curves refer to those of the “mix exposure” configuration. The upper plot shows the structure functions, observed in the mock observations in dash-dotted lines, and the modeled in plain lines. The total error is over-plotted on the modeled structure function with error bars, as its computation is part of the modeling, as described in Sect. 4.2. The lower panel shows the value of χ = (data-model) and the errors.

In the text

Fig. 5. Error terms contributing to the total error computation as a function of separation. The top panel represents the error models for the “mix exposure” computed with the input turbulent parameters. The lower panel represents these same errors computed for the “uniform exposure”. Four terms are represented, going from darker to lighter shades: the sample variance, the velocity field fluctuation, the structure function fluctuation, and the statistical fluctuation.

In the text

	Fig. 6. Structure function models for our simulation of 19 X-IFU pointings. The red curves show results for the “uniform exposure” configuration.
In the text

Fig. 7. Posterior distribution of the turbulent parameters obtained through the MCMC sampling. The blue contour show the distribution for the mixed exposure observation and the red for the uniform exposure. The inner contour delineates the 1σ level of the posterior probability, the second contour is at 2σ. Dotted lines show the input values of the simulations. Uniform priors were used for all parameters, except Ldiss, for which the prior is plotted in black.

In the text

	Fig. 8. Similar to Fig. 7, except that a Gaussian prior is set on parameter α, instead of an uniform prior.
In the text

Fig. 9. Posterior distribution obtained with setting the dissipation scale to the input value and letting the injection scale, norm and slope free. A Gaussian prior was set on the slope with π(α)∼𝒩(−11/3, 0.5), whereas the other parameters had uniform priors. The blue contour show the posterior distribution for the Mix exposure observation and the red for the uniform exposure. The prior is plotted in black. Dotted lines show the input values of the simulations.

In the text

	Fig. B.1. Geometrical quantities used in counting ‘ext’ pairs of length s from a point θ located within the circular region of analysis, 𝒜, of radius, R, at a distance θ = |θ| from the center. Whether pairs are pointing northward or southward, their number depends on the angle φ+ or φ− and we have φ = φ+ + φ−. Symmetrical relationships arise for the number of ‘inner’ pairs, with angles ϕ+ and ϕ−.
In the text

	Fig. B.2. Representation of the functional fs(θ). It takes value in [ − π, π] for each point θ = (y, z) located within the circular region of analysis. This region of radius, R, is bounded by the plain black line.
In the text

	Fig. D.1. Approximate relative uncertainty (sample variance only) in the two cases, A and B, and for the limiting regimes of small and large separations, s. In this figure, we fix the bin size, ℓ, the area, 𝒮𝒜, and the slope, α, letting only the injection scale vary (x-axis, in units of the bin size, ℓ).
In the text

	Fig. D.2. Approximate relative uncertainty (sample variance only) in the two cases, A and B, and for the limiting regimes of small and large separations, s. In this figure we fix the injection scale, Linj, the area, 𝒮𝒜, and the slope, α, letting only the bin size, ℓ,  vary (x-axis, in units of Linj).
In the text

	Fig. D.3. Exact relative error computation of the structure functions for the mix exposure observation, as a function of separation. Only the sample variance is included in the error budget.
In the text