11institutetext: IRAP, Université de Toulouse, CNRS, CNES, UT3-PS, Av. du
Colonel Roche 9, 31400, Toulouse, France 22institutetext: Centre National
d’Etudes Spatiales, Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401
Toulouse Cedex 9, France

# Toward mapping turbulence in the intracluster medium

III. Constraints on the turbulent power spectrum with Athena/X-IFU

S. Beaumont 1122 A. Molin Corresponding author : A. Molin -
<EMAIL_ADDRESS>N. Clerc 11 E. Pointecouteau 11 M. Vanel  E.
Cucchetti 22 P. Peille 22 F. Pajot 11

(Received …)

###### Abstract

Context. Future X-ray observatories with high spectral resolution and imaging
capabilities will enable measurements and mappings of emission line shifts in
the intracluster medium (ICM). Such direct measurements can serve as unique
probes of turbulent motions in the ICM. Determining the level and scales of
turbulence will improve our understanding of the galaxy cluster dynamical
evolution and assembly, together with a more precise evaluation of the non
thermal support pressure budget. This will allow for more accurate constraints
to be placed on the masses of galaxy clusters, among other potential benfits.

Aims. In this view, we implemented the methods presented in the previous
instalments of our work to characterising the turbulence in the intra-cluster
medium in a feasibility study with the X-ray Integral Field Unit (X-IFU) on
board the future European X-ray observatory, Athena.

Methods. From idealized mock observations of a toy model cluster, we
reconstructed the second-order structure function built with the observed
velocity field to constrain the turbulence. We carefully accounted for the
various sources of errors to derive the most realistic and comprehensive error
budget within the limits of our approach. With prior assumptions on the
dissipation scale and power spectrum slope, we constrained the parameters of
the turbulent power spectrum model through the use of Markov chain Monte Carlo
(MCMC) sampling.

Results. With a very long exposure time, a favourable configuration, and a
prior assumption of the dissipation scale, we were able to retrieve the
injection scale, velocity dispersion, and power spectrum slope, with 1$\sigma$
uncertainties for better than $\sim$15% of the input values. We demonstrated
the efficiency of our carefully set framework to constrain the turbulence in
the ICM from high-resolution X-ray spectroscopic observations, paving the way
for more in-depth investigation of the optimal required observing strategy
within a more restrictive observational setup with the future Athena/X-IFU
instrument.

###### Key Words.:

Galaxy clusters – Turbulence – Athena/X-IFU

##  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 spiraling 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^{7}-10^{8}$ 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 & Sa̧dowski, 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 (Aharonian et al., 2016;
Aharonia et al., 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$\alpha$ 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
$\sim$150eV 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 $\sim$5 eV resolution over its 6$\times$6 pixels matrix (over
3$\times$3 arcmin2), albeit with a limited spatial resolution of $\sim$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., 2022), 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 ($\sim$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 $\Lambda$CDM
cosmology, with $H_{0}=70$ km s-1 Mpc-1, $\Omega_{m}=0.3$ and
$\Omega_{\Lambda}=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:

$SF(\mathbf{r})=\overline{(\emph{v}(\mathbf{x}+\mathbf{r})-\emph{v}(\mathbf{x}))^{2}},$
(1)

where the over-line denotes averaging over the different pairs of sky points x
and x+r, separated by a distance, $r=|\mathbf{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 :

$SF(s)=\frac{1}{N_{p}(s)}\sum_{s}|C_{\mathcal{W}^{\prime}}-C_{\mathcal{W}}|^{2},$
(2)

where $s$ is the distance separating any pair of bins, $\mathcal{W}$ and
$\mathcal{W}^{\prime}$, while $N_{p}$ is the number of pairs separated by the
distance, $s$, and $C_{\mathcal{W}}$ 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$\langle SF\rangle$ and $\sigma_{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
$\sigma_{stat}$ to denote the standard deviation of the statistical error for
$C_{\mathcal{W}}$. 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_{\mathcal{W}}$, such that:

$\langle SF\rangle^{\prime}=\langle SF\rangle+2\sigma_{stat}^{2}.$ (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:

$\sigma_{sf,{\rm
tot}}^{2}=\sigma_{sf}^{2}+\sigma_{vf}^{2}+\sigma_{sff}^{2}+\sigma_{statf}^{2},$
(4)

Let us detail the terms in the equation above. $\sigma_{sf}^{2}$ the intrinsic
variance of the structure function. This is the residual variance in the
absence of statistical error, namely, if $\sigma_{stat}$=0.
The velocity field fluctuation term,
$\sigma_{vf}^{2}=4\,\sigma_{D}^{2}\,\sigma_{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_{D}^{2}$ is the variance of the first order structure function (as
defined in Eq. 18)
The structure function fluctuation term, written such that
$\sigma_{sff}^{2}=\frac{4}{N_{p}(s)}\,SF(s)\,\sigma_{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_{statf}^{2}=\frac{4}{N_{p}(s)}\,(N_{nei}(s)+1)\,\sigma_{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_{500}=1300$ kpc, corresponding to an angular radius of $11.7$ arcmin (where
$R_{500}$ 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, $\beta$-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:

$\dfrac{T(x)}{T_{500}}=T_{0}\dfrac{\nicefrac{{T_{min}}}{{T_{0}}}+(\nicefrac{{x}}{{r_{cool}}})^{a_{cool}}}{1+(\nicefrac{{x}}{{r_{cool}}})^{a_{cool}}}\dfrac{1}{\left(1+(\nicefrac{{x}}{{r_{t}}})^{2}\right)^{\nicefrac{{c}}{{2}}}},$
(5)

where $x=r/R_{500}$, and {$T_{0},T_{min},r_{cool}a_{cool},r_{t},c/2$} are six
parameters, and $T_{500}$ follows a parametrization on the redshift, see
Ghirardini et al. (2019). Their values are provided in Table 1.

Parameter | Value
---|---
$T_{0}$ | 1.09
$T_{min}/T_{0}$ | $0.66$
$\log(r_{cool})$ | -4.4
$a_{cool}$ | 1.33
$r_{t}$ | 0.45
$c/2$ | 0.30
$\gamma$ | 3
$\log(n_{0})$ [$\text{cm}^{-3}$] | -4.9
$\log(r_{c})$ [kpc] | -2.7
$\log(r_{s})$ [kpc] | -0.51
$\alpha$ | 0.70
$\beta$ | 0.39
$\epsilon$ | 2.60
Table 1: Values of the parameters used for the temperature and emission
models, obtained from the best-fit of these profiles to observed clusters
(Ghirardini et al., 2019).

The density profile follows:

$n_{e}^{2}(x)=n_{0}^{2}\dfrac{(\nicefrac{{x}}{{r_{c}}})^{-\alpha}}{(1+(\nicefrac{{x}}{{r_{c}}})^{2})^{3\beta-\alpha/2}}\dfrac{1}{(1+(\nicefrac{{x}}{{r_{s}}})^{\gamma})^{\nicefrac{{\epsilon}}{{\gamma}}}},$
(6)

where x = r/$R_{500}$ and {$\gamma,n_{0},r_{c},r_{s},\alpha,\beta,\epsilon$}
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-3, 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:

$Fe(r)=0.21(r+0.021)^{-0.48}-6.54\times\exp\left(-\frac{(r+0.0816)^{2}}{0.0027}\right),$
(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 $\sim$10% Helium in atom number (and weakly
contributing heavier elements), electron density, $n_{e}$ and hydrogen
density, $n_{H}$, are related by $n_{e}\sim 1.2\>n_{H}$. We imposed the
emissivity of the gas to vanish at radii beyond $5R_{500}$.

####  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_{\text{inj}}$ and $k_{\text{diss}}$ as the cut-
off frequencies of the Kolmogorov cascade at low and high wave numbers,
$\boldsymbol{\underline{k}}$ (with $k=1/r$ the conversion between scales and
wave numbers), $\alpha$ as the slope of the spectrum, and $\sigma$ as the 3D
velocity dispersion (units km s-1). With $k=|\boldsymbol{\underline{k}}|$, our
model can be expressed as:

$P_{3D}(\boldsymbol{\underline{k}})=\sigma^{2}\frac{\displaystyle
k^{-\alpha}e^{-(k_{\text{inj}}/k)^{2}}e^{-(k/k_{\text{diss}})^{2}}}{\displaystyle\int
4\pi
k^{2}k^{-\alpha}e^{-(k_{\text{inj}}/k)^{2}}e^{-(k/k_{\text{diss}})^{2}}dk}.$
(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 to 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 $\alpha=-11/3$, corresponding to the usual index of the
Kolmogorov cascade in the energy spectrum. The normalization constant,
$\sigma,$ is set such that the expected velocity broadening of the turbulent
field at the center of the cluster equals $(\mathcal{M}c_{s})^{2}$. Here,
$\mathcal{M}$ 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 $\sigma\simeq
250$ km s-1.

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}(\mathbf{k})=ve^{i\phi}$. To do so, we used a Rayleigh distribution for
the amplitude, $v$, and a uniform distribution for the phase, $\phi$. The
inverse Fourier transform of $V_{x}(\mathbf{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_{500}$, 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 $\mathcal{W}$) 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\cdot 10^{22}\text{cm}^{-2}$. 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. (2018).
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 :

$\displaystyle z_{v}$ $\displaystyle=\frac{z_{fit}-z}{1+z},$ (9)
$\displaystyle C$
$\displaystyle=\frac{(1+z_{v})^{2}-1}{(1+z_{v})^{2}+1}\>c_{light},$ (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_{v}c_{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., 2018).

####  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:

$\langle
SF(s)\rangle=\sum_{\boldsymbol{\underline{k}}}P_{3D}(\boldsymbol{\underline{k}})\\\
\times\left(\frac{1}{N_{p}(s)}\sum_{d(\mathcal{W},\mathcal{W}^{\prime})=s}\left|\frac{c_{\epsilon\cdot\mathcal{W}}(\boldsymbol{\underline{k}})}{F_{\mathcal{W}}}-\frac{c_{\epsilon\cdot\mathcal{W}^{\prime}}(\boldsymbol{\underline{k}})}{F_{\mathcal{W}^{\prime}}}\right|^{2}\right),$
(11)

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

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 $\langle SF\rangle$.

For a general emissivity model, it is necessary to compute the 3D Fourier
transform of
$\epsilon(x,\boldsymbol{\theta})\,\mathcal{W}(\boldsymbol{\theta})$ for each
bin $\mathcal{W}$. 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 $\boldsymbol{\theta}_{\rm eff}$,
representative of the bin $\mathcal{W}$ (but changing from bin to bin) such
that:

$\forall\>(\boldsymbol{\theta}\,\in\,\mathcal{W})\>:\>\epsilon(x,\boldsymbol{\theta})\approx\epsilon(x,\boldsymbol{\theta}_{\rm
eff})\equiv\varepsilon_{\boldsymbol{\theta}_{\rm eff}}(x).$ (12)

Hence, the simplified and faster computation of $c_{\epsilon\cdot\mathcal{W}}$
through:

$\displaystyle
c_{\epsilon\cdot\mathcal{W}}\approx\text{FT}_{1D}(\varepsilon_{\boldsymbol{\theta}_{\rm
eff}}(x))\;\text{FT}_{2D}(\mathcal{W}(\boldsymbol{\theta})).$ (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, $\sigma_{sf}$, the velocity
field fluctuations, $\sigma_{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 $\sigma_{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 $\sigma_{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.

Figure 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 $\sigma_{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) shows 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:

$Var[SF(s)]=16\pi\int_{0}^{\infty}\frac{1}{\mathcal{S}_{\mathcal{A}}}P_{l}^{2}(\xi)P_{2D}^{\infty
2}(\xi)(1-J_{0}(2\pi\xi s))^{2}\xi{\rm d}\xi$ (14)

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

$P_{2D}^{\infty}=2\int_{0}^{\infty}P_{3D}(\sqrt{k_{x}^{2}+\xi^{2}})(1+2\pi
k_{x}R_{c})^{2}e^{-4\pi k_{x}R_{c}}{\rm d}k_{x},$ (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, $\mathcal{S}_{\mathcal{A}}$,
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 $\sigma_{sf}/\langle SF\rangle$. 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, $\sigma_{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: $\sigma_{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
$\langle SF\rangle$ and $\sigma_{D}^{2}$, both depending linearly on the 3D
power spectrum. This can be explicitly expressed in the following way:

$\left.\sigma_{D}^{2}(s)\right|_{X}=\left.\sigma_{D}^{2}(s)\right|_{X_{0}}\frac{\left.\langle
SF(s)\rangle\right|_{X}}{\left.\langle SF(s)\rangle\right|_{X_{0}}},$ (16)

where we call $X$ the set of turbulent parameters and $X_{0}$ 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).
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.

Figure 2: Illustration of the ’mix exposure’ pointing strategy with exposure
times in ks at the center of each pointing. Figure 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 $R_{500}$
radius. All the spatial bins were constructed such as to provide identical
signal to noise levels in the resulting spectra.

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 $\sigma_{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 $\sim$30 km s-1.

###  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 $\sigma_{stat}$.

Figure 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 $\chi=$(data-model) and the errors.

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 Equations 12–13 is more likely to fail and
to result in a stronger mismatch with respect to the data.
We illustrate the different contributions of the error terms to the total
error in Figure 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 $\sigma_{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_{stat}\propto\frac{1}{\mathrm{SNR}}\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 SNR constant and to
adapt the binning accordingly. Then, even with lower exposure times, the
$\sigma_{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.

Figure 5: Plot of the 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 Figure 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.

Figure 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 10000 steps each. We
checked proper convergence of the chains with the $\hat{R}$ statistic (Vehtari
et al., 2019), which is 1.01¡$\hat{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:

$\displaystyle\log\left(p(SF_{obs}|X)\right)$
$\displaystyle=-0.5\sum_{s}\biggl{[}\log(2\pi\sigma_{sf,{\rm tot}}^{2})$
$\displaystyle+\frac{(\langle SF(s,X)\rangle-SF_{obs}(s))^{2}}{\sigma_{sf,{\rm
tot}}^{2}}\biggr{]},$

with $X$ as the set of turbulent parameters, $\langle SF(s,X)\rangle$ as the
modeled SF described in Sect. 4.1 and $\sigma_{sf,{\rm 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
$\pi(k_{\mathrm{inj}})\sim\mathcal{U}(0,k_{\mathrm{diss}}),\pi(\alpha)\sim\mathcal{U}(-10,0),\pi(\sigma^{2})\sim\mathcal{U}(0,500^{2})$.
Because $k_{\mathrm{diss}}$ can span the entire range [0;+$\infty$[, we set
the prior to the exponential distribution, such that
$\pi(k_{\mathrm{diss}})=\lambda e^{-\lambda k_{\mathrm{diss}}}$. This follows
from the exponential distribution being the maximum entropy distribution on
the [0;+$\infty$[ 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 Figure 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-$\sigma$ confidence regions, and the marginalized distributions show the
1-$\sigma$ 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-$\sigma$ 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.

Figure 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$\sigma$ level of the posterior probability, the second
contour is at 2$\sigma$. Dotted lines show the input values of the
simulations. Uniform priors were used for all parameters, except
$L_{\mathrm{diss}}$, 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<\alpha<-3$. Figure 8 shows the result of the same sampling procedure,
however, with $\pi(\alpha)\sim\mathcal{N}(-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.

Table 2: Summary table of the marginalized parameter uncertainties for each of the result figures presented in this section. The uncertainties given are $\pm 1\sigma$. The first column refers to the model and its corresponding figure. The second and third column show the priors used for each configuration. Model | $\pi(1/L_{diss})$ | $\pi(\alpha)$ | $L_{diss}[\mathrm{kpc}]$ | $L_{inj}[\mathrm{kpc}]$ | $\alpha$ | $\sigma[\mathrm{km}/\mathrm{s}]$
---|---|---|---|---|---|---
True Value | - | - | 10 | 300 | $-11/3$ | 253
Uniform exp. (Fig. 7) | Exp($\frac{1}{10}$) | – | $6.2^{+12.2}_{-4.7}$ | $244^{+73}_{-74}$ | $-3.98^{+0.86}_{-1.64}$ | $209^{+30}_{-17}$
Mixed exp. (Fig. 7) | Exp($\frac{1}{10}$) | – | $6.4^{+12.4}_{-5.0}$ | $215^{+37}_{-75}$ | $-4.7^{+1.2}_{-2.1}$ | $200^{+21.5}_{-9.2}$
Uniform exp. (Fig. 8) | Exp($\frac{1}{10}$) | $\mathcal{N}(-\frac{11}{3},0.5)$ | $6.6^{+12.0}_{-5.0}$ | $296^{+48}_{-38}$ | $-3.84^{+0.41}_{-0.40}$ | $233^{+26}_{-18}$
Mixed exp. (Fig. 8) | Exp($\frac{1}{10}$) | $\mathcal{N}(-\frac{11}{3},0.5)$ | $6.2^{+13.6}_{-4.9}$ | $275^{+39}_{-31}$ | $-3.91^{+0.37}_{-0.42}$ | $231^{+21}_{-17}$
Uniform exp. (Fig. 9) | – | $\mathcal{N}(-\frac{11}{3},0.5)$ | – | $301^{+48}_{-40}$ | $-3.78^{+0.37}_{-0.45}$ | $241^{+24}_{-17}$
Mixed exp. (Fig 9) | – | $\mathcal{N}(-\frac{11}{3},0.5)$ | – | $277^{+38}_{-31}$ | $-3.89^{+0.35}_{-0.43}$ | $238^{+19}_{-16}$
Figure 8: Similar to Fig. 7, except that a Gaussian prior is set on parameter
$\alpha$, 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.

Figure 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
$\pi(\alpha)\sim\mathcal{N}(-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
$\sigma_{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.,
2018). 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, $\sigma_{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 $\sigma_{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
Agency111ESA just rendered public the new specification of the Athena mission
in november 2023.. 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.

###### Acknowledgements.

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 et al., 2013, 2018, 2022), chainconsumer (Hinton, 2016), emcee
(Foreman-Mackey et al., 2013), matplotlib (Hunter, 2007) and cmasher (van der
Velden, 2020).

## References

  * Aharonia et al. (2018) Aharonia, F., Aharonian, F., Akamatsu, H., et al. 2018, Publications of the Astronomical Society of Japan, 70, 1
  * Aharonian et al. (2016) Aharonian, F., Akamatsu, H., Akimoto, F., et al. 2016, Nature, 535, 117
  * Anders & Grevesse (1989) Anders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, apj, 935, 167
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Barret (2022) Barret, D. 2022, in 44th COSPAR Scientific Assembly. Held 16-24 July, Vol. 44, 2316
  * Barret et al. (2022) Barret, D., Albouys, V., den Herder, J.-W., et al. 2022, arXiv e-prints, arXiv:2208.14562
  * Biffi et al. (2016) Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112
  * Biffi et al. (2022) Biffi, V., ZuHone, J. A., Mroczkowski, T., Bulbul, E., & Forman, W. 2022, Astronomy and Astrophysics, 663, 1
  * Cappellari (2009) Cappellari, M. 2009, arXiv e-prints, arXiv:0912.1303
  * Cavaliere et al. (2011) Cavaliere, A., Lapi, A., & Fusco-Femiano, R. 2011, ApJ, 742, 19
  * Churazov et al. (2012) Churazov, E., Vikhlinin, A., Zhuravleva, I., et al. 2012, Monthly Notices of the Royal Astronomical Society, 421, 1123
  * Clerc et al. (2019) Clerc, N., Cucchetti, E., Pointecouteau, E., & Peille, P. 2019, Astronomy & Astrophysics, 629 [arXiv:1904.06248]
  * Cucchetti et al. (2019) Cucchetti, E., Clerc, N., Pointecouteau, E., Peille, P., & Pajot, F. 2019, Astronomy & Astrophysics, 629, 1
  * Cucchetti et al. (2018) Cucchetti, E., Pointecouteau, E., Barret, D., et al. 2018, in Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, Vol. 10699, SPIE, 1141–1151
  * Cucchetti et al. (2018) Cucchetti, E., Pointecouteau, E., Peille, P., et al. 2018, A&A, 620, A173
  * Dauser et al. (2019) Dauser, T., Falkner, S., Lorenz, M., et al. 2019, Astronomy & Astrophysics, 630, 1
  * De Gasperin et al. (2017) De Gasperin, F., Intema, H. T., Shimwell, T. W., et al. 2017, Science Advances, 3, 1
  * Dupourqué et al. (2022) Dupourqué, S., Pointecouteau, E., Clerc, N., & Eckert, D. 2022, in European Physical Journal Web of Conferences, Vol. 257, mm Universe @ NIKA2 - Observing the mm Universe with the NIKA2 Camera, 00015
  * Dupourqué et al. (2022) Dupourqué, S., Pointecouteau, E., Clerc, N., & Eckert, D. 2022, EPJ Web of Conferences, 257, 00015
  * Eckert et al. (2021) Eckert, D., Gaspari, M., Gastaldello, F., Le Brun, A. M., & O’Sullivan, E. 2021, Universe, 7, 1
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gaspari & Sa̧dowski (2017) Gaspari, M. & Sa̧dowski, A. 2017, The Astrophysical Journal, 837, 149
  * Gatuzz et al. (2022) Gatuzz, E., Sanders, J. S., Canning, R., et al. 2022, MNRAS, 513, 1932
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, Astronomy and Astrophysics, 621 [arXiv:1805.00042]
  * Hinton (2016) Hinton, S. R. 2016, The Journal of Open Source Software, 1, 00045
  * Hlavacek-Larrondo et al. (2022) Hlavacek-Larrondo, J., Li, Y., & Churazov, E. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, 5
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Inogamov & Sunyaev (2003) Inogamov, N. A. & Sunyaev, R. A. 2003, Astronomy Letters, 29, 791
  * Ishisaki et al. (2022) Ishisaki, Y., Kelley, R. L., Awaki, H., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12181, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. J.-W. A. den Herder, S. Nikzad, & K. Nakazawa, 121811S
  * Kaastra & Bleeker (2016) Kaastra, J. S. & Bleeker, J. A. M. 2016, A&A, 587, A151
  * Kunz et al. (2022) Kunz, M. W., Jones, T. W., & Zhuravleva, I. 2022, 1
  * Lau et al. (2009) Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129
  * Markevitch & Vikhlinin (2007) Markevitch, M. & Vikhlinin, A. 2007, Physics Reports, 443, 1
  * McNamara & Nulsen (2012) McNamara, B. R. & Nulsen, P. E. J. 2012, New Journal of Physics, 14, 055023
  * Mernier & Biffi (2022) Mernier, F. & Biffi, V. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, 12
  * Mernier et al. (2018) Mernier, F., Biffi, V., Yamaguchi, H., et al. 2018, Space Sci. Rev., 214, 129
  * Mernier et al. (2017) Mernier, F., De Plaa, J., Kaastra, J. S., et al. 2017, Astronomy & Astrophysics, 603, 1
  * Molin et al. (2023) Molin, A., Clerc, N., Pointecouteau, E., Pajot, F., & Cucchetti, E. 2023, A&A, 679, A24
  * Nelson et al. (2014) Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25
  * Park & Bera (2009) Park, S. Y. & Bera, A. K. 2009, Journal of Econometrics, 150, 219
  * Pinto et al. (2020) Pinto, C., Fabian, A. C., Sanders, J. S., & de Plaa, J. 2020, Astronomische Nachrichten, 341, 217
  * Pratt et al. (2019) Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Science Reviews, 215 [arXiv:1902.10837]
  * Roncarelli et al. (2018) Roncarelli, M., Gaspari, M., Ettori, S., et al. 2018, Astronomy & Astrophysics, 618, 1
  * Sanders et al. (2020) Sanders, J. S., Dennerl, K., Russell, H. R., et al. 2020, Astronomy and Astrophysics, 633 [arXiv:1911.13108]
  * Schuecker et al. (2004) Schuecker, P., Finoguenov, A., Miniati, F., Böhringer, H., & Briel, U. G. 2004, A&A, 426, 387
  * Shi et al. (2015) Shi, X., Komatsu, E., Nelson, K., & Nagai, D. 2015, MNRAS, 448, 1020
  * Simionescu et al. (2019) Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Science Reviews, 215 [arXiv:1902.00024]
  * Smith et al. (2001) Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91
  * Takahasi & Mori (1974) Takahasi, H. & Mori, M. 1974, Publications of the Research Institute for Mathematical Sciences, 9, 721
  * Tejero-Cantero et al. (2020) Tejero-Cantero, A., Boelts, J., Deistler, M., et al. 2020, Journal of Open Source Software, 5, 2505
  * Ueda & Tashiro (2022) Ueda, Y. & Tashiro, M. 2022, in 44th COSPAR Scientific Assembly. Held 16-24 July, Vol. 44, 2312
  * van der Velden (2020) van der Velden, E. 2020, The Journal of Open Source Software, 5, 2004
  * Vazza et al. (2011) Vazza, F., Brunetti, G., Gheller, C., Brunino, R., & Brüggen, M. 2011, Astronomy & Astrophysics, 529, 1
  * Vazza et al. (2009) Vazza, F., Brunetti, G., Kritsuk, A., et al. 2009, Astronomy & Astrophysics, 504, 33
  * Vehtari et al. (2019) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. 2019, arXiv e-prints, arXiv:1903.08008
  * Voit (2018) Voit, G. M. 2018, The Astrophysical Journal, 868, 102
  * Zhuravleva et al. (2018) Zhuravleva, I., Allen, S. W., Mantz, A., & Werner, N. 2018, ApJ, 865, 53
  * Zhuravleva et al. (2012) Zhuravleva, I., Churazov, E., Kravtsov, A., & Sunyaev, R. 2012, Monthly Notices of the Royal Astronomical Society, 422, 2712
  * Zhuravleva et al. (2014) Zhuravleva, I., Churazov, E. M., Schekochihin, A. A., et al. 2014, The Astrophysical Journal Letters, 788, 1
  * ZuHone & Su (2022) ZuHone, J. & Su, Y. 2022, Handbook of X-ray and Gamma-ray Astrophysics, 1
  * ZuHone et al. (2016) ZuHone, J. A., Markevitch, M., & Zhuravleva, I. 2016, The Astrophysical Journal, 817, 110

##  Appendix A Modeling $\langle SF\rangle$ in arbitrary spatial binning
configurations

The model for the expectation value $\langle SF\rangle$ 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,\boldsymbol{\theta})$. 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
$\langle SF\rangle$ 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 $\langle SF\rangle$ 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_{\mathcal{W}}$, each associated to a spatial bin
$\mathcal{W}$:

$SF(s)=\frac{1}{N_{p}(s)}\sum_{d(\mathcal{W},\mathcal{W}^{\prime})=s}\left|C_{\mathcal{W}^{\prime}}-C_{\mathcal{W}}\right|^{2}.$

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 $\langle
C_{\mathcal{W}}C_{\mathcal{W}^{\prime}}\rangle$ for any two bins $\mathcal{W}$
and $\mathcal{W}^{\prime}$ (including the case of
$\mathcal{W}=\mathcal{W}^{\prime}$). We write the line centroid shift in bin
$\mathcal{W}$ as a flux-weighted average of centroid shifts projected along
the LoS (i.e., along coordinate $x$):

$\displaystyle C_{\mathcal{W}}$ $\displaystyle=F_{\mathcal{W}}^{-1}\int{\rm
d}\boldsymbol{\theta}\mathcal{W}(\boldsymbol{\theta})\int{\rm
d}x\epsilon(x,\boldsymbol{\theta})v(x,\boldsymbol{\theta})$
$\displaystyle=F_{\mathcal{W}}^{-1}\int{\rm
d}\boldsymbol{\theta}\mathcal{W}(\boldsymbol{\theta})F(\boldsymbol{\theta})C(\boldsymbol{\theta}),$

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

$C(\boldsymbol{\theta})=\sum_{\boldsymbol{\underline{k}}}V_{\boldsymbol{\underline{k}}}e^{i\omega\boldsymbol{\xi}\cdot\boldsymbol{\theta}}\int
e^{i\omega k_{x}x}\rho(x,\boldsymbol{\theta}){\rm d}x,$ (17)

where $\boldsymbol{\underline{k}}=(k_{x},\boldsymbol{\xi})$ is the Fourier
coordinate, $\omega=2\pi/L$, $L$ is a large characteristic length and
$\rho(x,\boldsymbol{\theta})\equiv\epsilon(x,\boldsymbol{\theta})/F(\boldsymbol{\theta})$.

Using the definition of the velocity power-spectrum $P_{3D}$ and writing
complex conjugation with symbol ∗, we find:

$\langle
C_{\mathcal{W}}C_{\mathcal{W}^{\prime}}\rangle=F_{\mathcal{W}}^{-1}F_{\mathcal{W}^{\prime}}^{-1}\sum_{\boldsymbol{\underline{k}}}P_{3D}(\boldsymbol{\underline{k}})c_{\epsilon\cdot\mathcal{W}}(\boldsymbol{\underline{k}})c_{\epsilon\cdot\mathcal{W}^{\prime}}^{*}(\boldsymbol{\underline{k}}).$

The complex quantity
$c_{\epsilon\cdot\mathcal{W}}(\boldsymbol{\underline{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
$(\epsilon\cdot\mathcal{W})$, namely:

$c_{\epsilon\cdot\mathcal{W}}(\boldsymbol{\underline{k}})=\int
e^{i\omega\boldsymbol{\underline{k}}\cdot\boldsymbol{\underline{x}}}\mathcal{W}(\boldsymbol{\theta})F(\boldsymbol{\theta})\rho(\boldsymbol{\underline{x}}){\rm
d}\boldsymbol{\underline{x}},$

where the function $\mathcal{W}(\boldsymbol{\theta})$ takes value 1 within bin
$\mathcal{W}$ 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_{\epsilon\cdot\mathcal{W}}$ for each bin $\mathcal{W}$. In
general222The case of $\beta$-model emissivity fields benefits from a number
analytical expressions that may accelerate the process, see App. B in Clerc et
al. (2019)., this calculation requires numerical evaluations of the (3D)
Fourier transforms of the field $(\epsilon\cdot\mathcal{W})$. These can be
simplified and accelerated if the emissivity is constant within a bin
$\mathcal{W}$, regardless of the LoS $\boldsymbol{\theta}$, i.e. if
$\epsilon(x,\boldsymbol{\theta})=\varepsilon(x)$ and
$F(\boldsymbol{\theta})=F$. This is a relevant approximation if the bins are
small in comparison to typical emissivity gradients. There, we obtain:

$\displaystyle c_{\epsilon\cdot\mathcal{W}}(k_{x},\boldsymbol{\xi})$
$\displaystyle\simeq\widetilde{\varepsilon}(k_{x})\times\widehat{\mathcal{W}}(\boldsymbol{\xi}),$
$\displaystyle F_{\mathcal{W}}$ $\displaystyle\simeq F\times S_{\mathcal{W}.}$

In these expressions, $S_{\mathcal{W}}$ 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}={\rm 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:

$D(s)=\frac{1}{N^{in}(s)}\sum_{in}\left[C(\boldsymbol{\theta})-C(\boldsymbol{\theta}^{\prime})\right].$
(18)

Similarly as in Appendix A, $C(\boldsymbol{\theta})$ represents the line
centroid shift along a LoS labeled by the sky coordinate
$\boldsymbol{\theta}$. We arbitrarily set $C(\boldsymbol{\theta})=0$ for
locations outside of the region of analysis. The sum runs over all pairs of
coordinates $(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})$ 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^{in}$ 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 $\boldsymbol{\theta}$ always lies
southwards of $\boldsymbol{\theta}^{\prime}$. 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. 10). 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.

Figure 10: Geometrical quantities used in counting ‘ext’ pairs of length $s$
from a point $\boldsymbol{\theta}$ located within the circular region of
analysis, $\mathcal{A,}$ of radius, $R$, at a distance
$\theta=|\boldsymbol{\theta}|$ from the center. Whether pairs are pointing
northward or southward, their number depends on the angle $\varphi^{+}$ or
$\varphi^{-}$ and we have $\varphi=\varphi^{+}+\varphi^{-}$. Symmetrical
relationships arise for the number of ‘inner’ pairs, with angles $\phi^{+}$
and $\phi^{-}$.

Let us introduce the new quantity $D^{\prime}$ through:

$D^{\prime}(s)=\frac{1}{N^{in}(s)}\sum_{in+ext}\left[C(\boldsymbol{\theta})-C(\boldsymbol{\theta}^{\prime})\right],$

in which the sum runs over all pairs (always assuming the conventional
orientation toward the north). Rearranging the summation terms and using that
$C(\boldsymbol{\theta})=0$ outside of the region of analysis, we find
$D^{\prime}(s)=0$.

We are thus left with calculating the sum over all pairs of kind ‘ext’. Those
pairs having their southern end $\boldsymbol{\theta}$ within the analysis
region contribute a term $C(\boldsymbol{\theta})$; while those having their
southern end outside of the region of analysis contribute a term
$-C(\boldsymbol{\theta}^{\prime})$. To each direction, $\boldsymbol{\theta},$
within the field of analysis, $\mathcal{A,}$ we associate the number of ”ext”
pairs: $\varphi_{s}(\boldsymbol{\theta})$. This number can be decomposed into
pairs that point northwards ($\varphi^{+}_{s}$) and pairs that point
southwards ($\varphi^{-}_{s}$). We write:

$\sum_{ext}\left[C(\boldsymbol{\theta})-C(\boldsymbol{\theta}^{\prime})\right]=\sum_{\boldsymbol{\theta}\in\mathcal{A}}\left(\varphi^{+}_{s}(\boldsymbol{\theta})-\varphi^{-}_{s}(\boldsymbol{\theta})\right)C(\boldsymbol{\theta})$

Introducing the number of ”in” pairs
$\phi_{s}(\boldsymbol{\theta})=2\pi-\varphi_{s}(\boldsymbol{\theta})$, it is
clear that the following relations hold (e.g., see Fig. 10):

$\displaystyle\varphi_{s}=\varphi^{+}_{s}+\varphi^{-}_{s}$ ;
$\displaystyle\phi_{s}=\phi_{s}^{+}+\phi_{s}^{-};$
$\displaystyle\phi_{s}^{+}=\pi-\varphi_{s}^{+}$ ;
$\displaystyle\phi_{s}^{-}=\pi-\varphi_{s}^{-}.$

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

Let us describe the position of the point $\boldsymbol{\theta}$ by two
coordinates: $\theta=|\boldsymbol{\theta}|$ its distance to the center of the
circular shape of $\mathcal{A}$ ; and $\psi$ is the angle
$\boldsymbol{\theta}$ makes with the horizontal line. Geometrical
considerations lead to distinguish five mutually exclusive cases:

$\phi_{s}^{+}(\theta,\psi)=\left\\{\begin{array}[]{cl}\pi&{\rm if\
}\tau_{s}(\theta,\psi)>\pi,\\\ \phi_{s}(\theta)-\pi&{\rm if\
}\tau_{s}(\theta,\psi)<\phi_{s}(\theta)-\pi,\\\ 0&{\rm if\
}\tau_{s}(\theta,\psi)<0,\\\ \phi_{s}(\theta)&{\rm if\
}\tau_{s}(\theta,\psi)>\phi_{s}(\theta),\\\
\tau_{s}(\theta,\psi)&\mathrm{otherwise,}\\\ \end{array}\right.$

where $\tau_{s}(\theta,\psi)=\phi_{s}(\theta)/2-\psi$ and
$\tau_{s}\in[-\pi/2,3\pi/2]$. Figure 11 shows one example ($s=0.4R$)
representation of the function
$f_{s}(\boldsymbol{\theta})=\varphi^{+}_{s}(\boldsymbol{\theta})-\varphi^{-}_{s}(\boldsymbol{\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.

Figure 11: Representation of the functional $f_{s}(\boldsymbol{\theta})$. It
takes value in $[-\pi,\pi]$ for each point $\boldsymbol{\theta}=(y,z)$ located
within the circular region of analysis. This region of radius, $R,$ is bounded
by the plain black line.

Using that $D^{\prime}=0$, we write:

$D(s)=-\frac{1}{N^{in}(s)}\sum_{\theta\in\mathcal{A}}f_{s}(\boldsymbol{\theta})C(\boldsymbol{\theta}).$

Since $\langle C\rangle=0$, we have that $\langle D\rangle=0$ and
$Var(D)=\langle D^{2}\rangle$. Developing the calculation involves the
quantity $\left\langle
C(\boldsymbol{\theta}_{1})C(\boldsymbol{\theta}_{2})\right\rangle$ and Eq. 17
may be used to obtain:

$Var(D)=\frac{1}{(N^{in})^{2}}\sum_{\boldsymbol{\underline{k}}}P_{3D}(k)\left|\int_{\boldsymbol{\theta}\in\mathcal{A}}{\rm
d}\boldsymbol{\theta}f_{s}(\boldsymbol{\theta})e^{i\omega\boldsymbol{\xi}\cdot\boldsymbol{\theta}}\widetilde{\chi^{\boldsymbol{\theta}}}(k_{x})\right|^{2}.$

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

$Var[D(s)]=\frac{1}{N^{in}(s)^{2}}\left(\frac{\omega}{2\pi}\right)^{4}\sum_{\boldsymbol{\underline{k}}}P_{3D}(k)\left|\left(u_{k_{x}}\otimes\widehat{f_{s}}\right)(\boldsymbol{\xi})\right|^{2},$
(19)

involving a series of 2D convolutions indicated with the sign $\otimes$. 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, $\ell$, all
supposedly of identical shape and regularly spaced. This approximation
replaces $C(\boldsymbol{\theta})$ above by a filtered version:

$C(\boldsymbol{\theta})\rightarrow\int{\rm
d}\boldsymbol{\mu}\mathcal{F}_{\ell}(\boldsymbol{\mu})C(\boldsymbol{\theta}-\boldsymbol{\mu}),$

where $\mathcal{F}_{\ell}$ 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 ${\rm Var}(D)$ is
obtained by replacing $\widehat{f_{s}}$ by the product
$\widehat{f_{s}}\widehat{\mathcal{F}_{\ell}}^{*}$ in Eq. 19.

##  Appendix C Plane-constant emissivity model

An useful simplified emissivity model consists of a $\beta$-model along the
LoS direction $x$, identical for any selected LoS $\boldsymbol{\theta}$. We
call this model Xbeta, it is given by:

$\epsilon(x,\boldsymbol{\theta})=\varepsilon(x)=\varepsilon(0)\left(1+\frac{x^{2}}{r_{c}^{2}+\theta_{0}^{2}}\right)^{-3\beta}.$
(20)

This expression involves a core-radius, $r_{c}$, and a fixed impact parameter,
$\theta_{0}$. 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 $\boldsymbol{\theta}$, 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. 20, 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):

${\rm
sf}(s)=4\pi\left(\frac{\omega}{2\pi}\right)^{2}\int_{0}^{\infty}P_{\ell}P_{2D}^{\infty}\left(1-J_{0}(\omega\xi
s)\right)\xi{\rm d}\xi,$ (21)

${\sigma_{\rm
sf}^{2}}(s)=\frac{16\pi}{\mathcal{S}_{\mathcal{A}}}\left(\frac{\omega}{2\pi}\right)^{2}\int_{0}^{\infty}P_{\ell}^{2}{P_{2D}^{\infty}}^{2}\left(1-J_{0}(\omega\xi
s)\right)^{2}\xi{\rm d}\xi,$ (22)

where ${\rm sf}=\langle SF\rangle$, $\sigma_{\rm sf}^{2}=\langle(SF-{\rm
sf})^{2}\rangle$. In these expressions, $\omega=2\pi/L$ with $L$ an arbitrary
length setting the length unit system ; $P_{\ell}(\xi)$ is the power spectrum
of the 2D bins, each of characteristic size, $\ell$ ; $P_{2D}^{\infty}(\xi)$
is the power spectrum of the centroid shift map, extrapolated to an unbounded
region of analysis; $\mathcal{S}_{\mathcal{A}}$ denotes the area of the region
of analysis; $J_{0}$ is the Bessel function of the first kind of order 0.
Furthermore, we specify the emissivity model following Eq. 20, involving a
characteristic size, $R_{c}$.

In what follows we will be interested in the relative sample variance, as
given by $\sigma_{\rm sf}/{\rm 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 $\ell$ and the Fourier transform of
their shape is well approximated by a top-hat function over the interval:
$-1/\ell<\omega\xi<1/\ell$.

  * •

Bins are small with respect to the projected injection length and larger than
the dissipation length of turbulent motions: $k_{\rm inj}<1/\ell<k_{\rm
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}$:

$P_{2D}^{\infty}(\xi)\propto\int
P_{3D}\left(\sqrt{k_{x}^{2}+\xi^{2}}\right)P_{\varepsilon}(k_{x}){\rm
d}k_{x},$ (23)

where $P_{\varepsilon}$ 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. 23 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:

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

Here, $A$ is a constant that cancels out in forming the ratio $\sigma_{\rm
sf}/{\rm sf}$. The integrals involved in expressions 21 and 22 extend from $0$
to $+\infty$. The expressions for $P_{\ell}$ and $P_{2D}$ show it is only
necessary to compute the integrals over the domain $[k_{inj},1/\ell]$. Two
regimes are then worth investigating.

At small separations $s\ll\ell$, the second-order Taylor expansion of the
Bessel function $J_{0}$ leads to:

$\frac{\sigma_{\rm sf}}{\rm{sf}}\simeq
2\sqrt{\pi}\frac{L_{inj}}{\mathcal{S}_{\mathcal{A}}^{1/2}}\frac{\mathcal{G}\left(r,2\alpha+6\right)^{1/2}}{\mathcal{G}\left(r,\alpha+4\right)}{\
\ \ \ \ \ {\rm if\ }s\ll\ell.}$ (24)

Here, we define $r=\ell k_{inj}\in[0,1]$ and the specific function:

$\mathcal{G}(r,n)=\frac{1}{n}\left(r^{-n}-1\right).$

At large separations, $s\gg L_{inj}$, assuming that the Bessel function
vanishes toward large values of its argument:

$\frac{\sigma_{\rm sf}}{\rm{sf}}\simeq
2\sqrt{\pi}\frac{L_{inj}}{\mathcal{S}_{\mathcal{A}}^{1/2}}\frac{\mathcal{G}\left(r,2\alpha+2\right)^{1/2}}{\mathcal{G}\left(r,\alpha+2\right)}{\
\ \ \ \ \ {\rm if\ }s\gg L_{inj}}.$ (25)

###  D.2 Study of case B ($k_{inj}<1/R_{c}<k_{diss}$)

In case B, the integral in Eq. 23 decomposes into a sum over three intervals:
$[0,k_{inj}]$, $[k_{inj},1/c]$ and $[1/R_{c},+\infty]$. A close examination of
the integrals for various values of $\xi$ leads to the following
approximation:

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

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

$\frac{\sigma_{\rm sf}}{\rm{sf}}\simeq
2\sqrt{\pi}\frac{L_{inj}}{\mathcal{S}_{\mathcal{A}}^{1/2}}\frac{\left[\frac{1}{6}+\mathcal{G}\left(r,2\alpha+6\right)\right]^{1/2}}{\frac{1}{4}+\mathcal{G}\left(r,\alpha+4\right)}{\
\ \ \ \ \ {\rm if\ }s\ll\ell,}$ (26)

$\frac{\sigma_{\rm sf}}{\rm{sf}}\simeq
2\sqrt{\pi}\frac{L_{inj}}{\mathcal{S}_{\mathcal{A}}^{1/2}}\frac{\left[\frac{1}{2}+\mathcal{G}\left(r,2\alpha+2\right)\right]^{1/2}}{\frac{1}{2}+\mathcal{G}\left(r,\alpha+2\right)}{\
\ \ \ \ \ {\rm if\ }s\gg L_{inj}}.$ (27)

###  D.3 Interpretation

Under the assumptions stated in this Appendix, the relative uncertainties due
to sample variance scale as provided in Eq. 24, 25, 26, and 27, 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 $\ell$ (typical size of the
spatial binning scheme).

In all four limiting cases, it clearly appears that relative uncertainties
scale with the ratio $L_{inj}/\sqrt{\mathcal{S}_{\mathcal{A}}}$ at fixed
$r=\ell/L_{inj}$. Larger regions of analysis that sample multiple (projected)
injection lengths are less prone to sample variance, as one might expect.
Figure 12 is obtained by fixing the bin size $\ell$ 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 $\ell$ and the slope of
turbulence, $\alpha,$ through the functional $\mathcal{G}$ 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 13 shows the dependency of relative uncertainties upon the bin size
$\ell$, 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.

Figure 12: 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, $\ell$, the area,
$\mathcal{S}_{\mathcal{A}}$, and the slope, $\alpha$, letting only the
injection scale vary (x-axis, in units of the bin size, $\ell$). Figure 13:
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, $L_{inj}$, the area,
$\mathcal{S}_{\mathcal{A}}$, and the slope, $\alpha$, letting only the bin
size, $\ell,$ vary (x-axis, in units of $L_{inj}$).

###  D.4 Exact calculations

The integrals in Eqs. 21 and 22 (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. 14, 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.

Figure 14: 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.