11institutetext: IRAP, Université de Toulouse, CNRS, UPS, CNES, Toulouse,
France
11email<EMAIL_ADDRESS>22institutetext: CNES, 18 Avenue Edouard
Belin, 31400 Toulouse, France

# Towards mapping turbulence in the intra-cluster medium

I. Sample variance in spatially-resolved X-ray line diagnostics

Nicolas Clerc 11 Edoardo Cucchetti 11 Etienne Pointecouteau 11 Philippe Peille
22

(Received April 2019; accepted July 2019.)

###### Abstract

Context.  X-ray observations of galaxy clusters provide insights on the nature
of gaseous turbulent motions, their physical scales and on the fundamental
processes they are related to. Spatially-resolved, high-resolution spectral
measurements of X-ray emission lines provide diagnostics on the nature of
turbulent motions in emitting atmospheres. Since they are acting on scales
comparable to the size of the objects, the uncertainty on these physical
parameters is limited by the number of observational measurements, through
sample variance.

Aims.  We propose a different and complementary approach for the computation
of sample variance to repeating numerical simulations (i.e. Monte-Carlo
sampling) by introducing new analytical developments for lines diagnosis.

Methods.  We consider the model of a ”turbulent gas cloud”, consisting in
isotropic and uniform turbulence described by a universal Kolmogorov power-
spectrum with random amplitudes and phases in an optically thin medium.
Following a simple prescription for the 4-term correlation of Fourier
coefficients, we derive generic expressions for the sample mean and variance
of line centroid shift, line broadening and projected velocity structure
function. We perform a numerical validation based on Monte-Carlo simulations
for two popular models of gas emissivity based on the $\beta$-model.

Results.  Generic expressions for the sample variance of line centroid shifts
and broadening in arbitrary apertures are derived and match the simulations
within their range of applicability. Generic expressions for the mean and
variance of the structure function are provided and verified against
simulations. An application to the _Athena_ /X-IFU (_Advanced Telescope for
High-ENergy Astrophysics_ /X-ray Integral Field Unit) and _XRISM_ /Resolve
(_X-ray Imaging and Spectroscopy Mission_) instruments forecasts the potential
of sensitive, spatially-resolved spectroscopy to probe the inertial range of
turbulent velocity cascades in a Coma-like galaxy cluster.

Conclusions.  The formulas provided are of generic relevance and can be
implemented in forecasts for upcoming or current X-ray instrumentation and
observing programs.

##  1 Introduction

Galaxy clusters form by accretion of matter along filaments of the cosmic web,
either continuously or episodically through major and minor merger events. The
baryonic gas flowing along filamentary structures and falling into their deep
gravitational wells acquires kinetic energy that is transformed into thermal
energy, magnetic field amplification and cosmic ray acceleration in the intra-
cluster medium, by a succession of shocks, large-scale motions and dissipation
by turbulent processes (Ryu et al., 2008; Zhuravleva et al., 2014; Gaspari et
al., 2014; Miniati & Beresnyak, 2015; Gaspari et al., 2018; Vazza et al.,
2018). Observational signatures of these phenomena are rare and difficult to
obtain. The most promising and efficient diagnostics are issued from
spectroscopic observations of the hot intra-cluster gas which permeates the
entire volume of massive halos and emits copious amounts of X-ray light.

Focusing mainly on X-ray emission lines extracted along a single line-of-
sight, Inogamov & Sunyaev (2003) demonstrated that departures from Gaussian
line shapes carry important indications on the nature of large-scale
turbulence in the intra-cluster medium. The authors extended formalism to two-
dimensional diagnostics by introducing the correlation function of the
projected velocity field and calculating its scaling relative to fundamental
parameters such as the turbulent injection and dissipation scales. Applying
these findings to simple, but realistic configurations of the intra-cluster
medium, Zhuravleva et al. (2012) calculated exact expressions for emission
line diagnostics such as centroid shift, broadening and two-dimensional
correlation function. They evaluated the associated sampling uncertainty (also
called ’sample variance’ or ’sampling variance’) by multiple Monte-Carlo
realisations of the velocity field and showed that it can dominate the overall
error budget in presence of large-scale turbulence. ZuHone et al. (2016) could
evaluate the contribution of sample variance and statistical errors for the
well-defined observational case of the Coma cluster, thereby demonstrating the
impact of the observational strategy on this source of uncertainty. Using
numerically simulated clusters instead, Roncarelli et al. (2018) performed
end-to-end simulations to derive expected values of the indicators of
turbulence issued from emission line measurements, postponing calculation of
sample variance to a later stage by means of multiple realisations.

In this work we propose a formal approach to the problem of sample variance by
considering the ideal case of an arbitrary, optically-thin gas distribution in
which uniform and isotropic turbulent motions take place. This study is
motivated by the intent to obtain reliable and fast estimates of this specific
class of uncertainties and to identify key parameters impacting them. We will
consider three popular diagnostics extracted from a continuum-free, isolated
spectral line in the X-ray wavebands: line centroid shift (hereafter $C$),
line broadening ($S$) and projected velocity structure function ($SF$). The
latter is defined as the squared difference of projected velocities averaged
among all points separated by a distance $s$ on sky. Instrumental
characteristics and signal-to-noise considerations related to, e.g. the
exposure time or the energy resolution, are deliberately excluded and
addressed in a separate work (Cucchetti et al., in press, hereafter paper II)
The results of the present work are therefore instrument-independent to some
extent.

Among these indicators, the structure function appears as a very promising
diagnostic of turbulence since it takes advantage of spatially-resolved
spectroscopic observations, as enabled by Integral Field Units. It is also the
least intuitive of all three. Effects such as heterogenous sampling, non-
stationarity, anisotropies, etc. reflect diversely in the modelling of $SF$.
Interestingly, the structure function as a mathematical tool has received
extensive interest in multiple fields of research involving spatial
statistics, notably geostatistics and Earth science, under the name
’variogram’ (e.g. Matheron, 1965, 1973; Cressie, 1985; Haslett, 1997;
Armstrong, 1998; Corstanje et al., 2008). In the field of astronomy and
astrophysics where its use is comparatively less widespread, it is involved in
various works under both terms ’structure function’ and ’variogram’, to
analyse data either in one dimension (e.g. Roelens et al., 2017, for stellar
variability), two (e.g. Cayón, 2010, for Cosmic Microwave Background) or three
and more dimensions (e.g. Martínez et al., 2010, for galaxy clustering).

We first introduce the derivation of the average and variance of the centroid
shift and line broadening for measurements of an X-ray spectral line along a
single line-of-sight, together with a numerical validation (Section 2).
Section 3 generalises these results to the case of three-dimensional turbulent
fields and the extension of the line diagnostics to two dimensions, thereby
treating the case of the structure function. We perform a numerical validation
of these results in Section 4. We finally discuss our results in Section 5 and
highlight two specific cases matching the future X-ray instruments _XRISM_
/Resolve (Ishisaki et al., 2018) and _Athena_ /X-IFU (Barret et al., 2016,
2018). We report most of the details on calculations and their discussions in
the appendices, to which the reader can refer for more details.

The convention in our notations is as follows: the line-of-sight direction is
denoted by $x$ and the plane-of-sky coordinate is $\@vec{\theta}=(y,z)$. Units
of these coordinates are physical (kpc) since in practice the angular distance
at the redshift of the object is known. Three-dimensional vectors are
underlined to differentiate them from two-dimensional vectors. The velocity
$v$ (units km s-1) is the component of the gas velocity projected along the
line-of-sight. All following definitions and derivations (e.g. turbulent
velocity dispersion, power-spectrum, etc.) are relative to this line-of-sight
component. We denote with brackets $\langle.\rangle$ the sample average of the
estimators and random variables. We will decompose the velocity field in
Fourier coefficients with discrete indices (involving the discrete summation
sign $\sum_{k}$). The emissivity and geometrical shapes will be treated with
their continuous Fourier transforms (involving the continuous summation sign
$\int{\rm d}k$). This distinction will often be purely formal: this is the
choice made for clarifying the calculations. One- and three-dimensional
Fourier transforms are indicated with a tilde (e.g. $\widetilde{\rho}$), two-
dimensional transforms with a hat (e.g. $\widehat{\mathcal{W}}$).

##  2 Measured velocity dispersion along single line-of-sight

In this section we assume the velocity structure diagnostics are issued from
the measurement of an emission line profile (e.g. iron XXV at $\sim 6.5$ keV)
along a given line of sight. Measuring a line profile is a complex task
involving tools and methods developed under a certain set of observational
conditions (binning of the spectra, level of noise, background subtraction,
continuum subtraction, etc.) In order to illustrate our findings, we adopt a
simplified approach where the analysis applies to a continuum- and background-
subtracted spectrum with no source of noise nor uncertainty and no systematic
(corresponding to a virtually infinite exposure time with a perfectly
calibrated instrument) . Only one emission line is investigated, thus we
neglect blending with neighbouring lines. Importantly, we do not provide a
prescription for the measurement process itself (Gaussian or more complex fit,
non-parametric fit, full spectrum fitting, etc.) Instead we model such
measurement as a calculation of the zero-th, first and second moments of the
line energy distribution $I_{l}(E)$ integrated along a line-of-sight
$\@vec{\theta}_{0}$, such that:

$\displaystyle F(\@vec{\theta}_{0})$ $\displaystyle=\int
I_{l}(\@vec{\theta}_{0};E){\rm d}E$ $\displaystyle\delta E(\@vec{\theta}_{0})$
$\displaystyle=\left(F^{-1}\int EI_{l}(\@vec{\theta}_{0},E){\rm
d}E\right)-E_{0}$ $\displaystyle\Sigma^{2}(\@vec{\theta}_{0})$
$\displaystyle=F^{-1}\int\left(E-\delta
E-E_{0}\right)^{2}I_{l}(\@vec{\theta}_{0};E){\rm d}E$

where $\delta E$ and $\Sigma$ are the observed centroid shift (relative to
reference energy $E_{0}$) and width (or broadening) of the line. $F$ is a
normalization factor, namely the flux in the line. We introduce the gas
velocity field along the line-of-sight fixed by $\@vec{\theta}_{0}$ with
$v(x)\equiv v(x,\@vec{\theta}_{0})$. We define $C=c\delta E/E_{0}$ and
$S^{2}=c^{2}\Sigma^{2}/E_{0}^{2}$ the observed centroid shift and width in
velocity space.

###  2.1 Emission along the line of sight

At microscopic level (i.e. below the turbulent dissipation scale in the
medium), emission is assumed to follow a thermally broadened line profile. We
neglect any additional broadening such as natural (Lorentzian) and assign each
ion a rest-frame Gaussian emission profile. Assuming purely collisional origin
of the emission line, the amplitude of the line is assumed to scale with
emissivity $\epsilon(\@vec{r})\propto n_{Fe}n_{e}=n_{e}^{2}$. The coefficient
of proportionality may depend on the local property of the medium
(metallicity, temperature, etc.)

In the following we assume that the turbulent dissipation scale $L_{diss}$ is
large enough for each volume of size $(L_{diss})^{3}$ to contain a
significantly large number of line emitters: it is practically always
fulfilled, even for the tenuous intra-cluster medium with typical density
$10^{-3}$ cm-3 and kpc-scale injection scales. Accounting for the Doppler
shift in energy $v(x)/c$ and emissivity $\epsilon(x)$, we therefore model the
total emission at each point by:

$\frac{{\rm d}I_{l}(E)}{{\rm
d}x}=\frac{\epsilon(x)c}{E_{0}\sigma_{th}(x)\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left[\frac{E-E_{0}(1+v(x)/c)}{E_{0}\sigma_{th}(x)/c}\right]^{2}\right)$
(1)

Assuming an optically thin medium and reordering the summation over velocities
we can write:

$\displaystyle F(\@vec{\theta}_{0})$ $\displaystyle=\int\epsilon(x){\rm d}x$
$\displaystyle C(\@vec{\theta}_{0})$
$\displaystyle=F^{-1}\int\epsilon(x)v(x){\rm d}x$ $\displaystyle
S^{2}(\@vec{\theta}_{0})$
$\displaystyle=F^{-1}\int\epsilon(x)\sigma_{th}^{2}(x){\rm
d}x+\frac{F^{-2}}{2}\iint G(x,x^{\prime}){\rm d}x{\rm d}x^{\prime}$ (2)
$\displaystyle G(x,x^{\prime})$
$\displaystyle=\epsilon(x)\epsilon(x^{\prime})\left[v(x)-v(x^{\prime})\right]^{2}$

These integrals extend along the line-of-sight indexed by $x$, that we assume
to range in a large interval $[-L,L]$. In principle $v(x)$ encapsulates the
effect of hydrodynamical motions (turbulence) and bulk motions of the gas.
Without loss of generality, we assume bulk motions over a small area are known
and subtracted from the measurements; we therefore set its contribution to
zero.

###  2.2 Turbulent velocity field

We describe the turbulent velocity field $v$ by its Fourier series expansion,
with positive and negative values of $k$:

$v(x)=\sum_{k}V_{k}\exp{\left(\frac{i2\pi kx}{L}\right)}$

The coefficients $V_{k}$ are complex random variables, defined as
$V_{k}=V_{-k}^{*}=|V_{k}|e^{i\psi_{k}}$. Here $\psi_{k}$ is a random phase and
$|V_{k}|$ a random modulus, supposed independent from each other. In the
following, we define for convenience $\omega=2\pi/L$. We note that $\langle
v\rangle=0$, leading to $V_{0}=0$. Averaging over multiple random realisations
provides:

$\langle V_{j}V_{k}\rangle=\delta_{j,-k}P(k)$ (3)

where $P(k)=P(-k)=\langle|V_{k}|^{2}\rangle$ is the power-spectrum of the
turbulent velocity field and $\delta_{ij}=1$ if $i=j$, 0 if $i\neq j$.

Our hypothesis of uniform turbulence implies that the normalization of the
power spectrum matches the square of the turbulent velocity dispersion
$\sigma_{turb}$ at any given point $x$. It is defined through the calculation
of the second moment $\sigma_{turb}=\sqrt{\langle v^{2}\rangle}$, where
averaging occurs over random phases and moduli. Therefore (taking e.g. $x=0$)
$\sigma_{turb}^{2}=\sum_{k}P(k)$. This definition does not involve the profile
of the emissivity, in contrast to e.g. ZuHone et al. (2016).

One simple assumption for the distribution of moduli (e.g. Inogamov & Sunyaev,
2003) consists in non-random coefficients, leaving only phases as random.

Another popular and physically motivated assumption (although not
systematically required in the rest of this paper) is the Rayleigh
distribution (e.g. ZuHone et al., 2016):

$|V_{k}|\sim\mathcal{P}(\nu){\rm d}\nu=\frac{2\nu}{P(k)}e^{-\nu^{2}/P(k)}{\rm
d}\nu$

which reflects that the turbulent velocity is a Gaussian random field.
Introducing $R_{k}=P(k)^{2}-\mathrm{Var}(|V_{k}|^{2})$, we obtain under the
assumption of Rayleigh-distributed moduli: $R_{k}=0$ for all $k$.

###  2.3 Statistics of the centroid shift and line broadening

Calculations in App. A provide the following formulas for the centroid shift:

$\langle C\rangle=0$ (4)

and its variance:

$\mathrm{Var}(C)=\langle C^{2}\rangle=F^{-2}\sum_{k}P_{\epsilon}(k)P(k)$ (5)

This expression is identical to Eq. (A9) in Zhuravleva et al. (2012), since
$P_{\epsilon}(k)=|\widetilde{\epsilon}(k)|^{2}$ is the Fourier power spectrum
of the (unnormalised, one-dimensional) emissivity $\epsilon$.

As for the line broadening, we obtain in App. A:

$\langle S^{2}\rangle=\overline{\sigma_{th}^{2}}+\sigma_{turb}^{2}-\langle
C^{2}\rangle$ (6)

The horizontal bar denotes averaging of the thermal component along the line-
of-sight. Again this expression is similar to Eq. B4 in Zhuravleva et al.
(2012). Also interesting is the contribution of the last term, indicating that
averaging many (independent) measurements of broadening measurements generally
provides a biased estimate of the (thermal+turbulent) broadening. The bias is
zero only in cases where the turbulent power-spectrum and the emissivity power
spectrum act on distinct spatial scales. Finally, the variance can be written:

$\mathrm{Var}(S^{2})=2\sum_{j,k}P(k)P(j)\left|\frac{\widetilde{\epsilon}(j+k)}{F}-\frac{\widetilde{\epsilon}(j)\widetilde{\epsilon}(k)}{F^{2}}\right|^{2}\\\
-\sum_{k}R_{k}\times\left\\{\left|\frac{\widetilde{\epsilon}(2k)}{F}-\frac{\widetilde{\epsilon}(k)^{2}}{F^{2}}\right|^{2}+2\left[1-\frac{P_{\epsilon}(k)}{F^{2}}\right]^{2}\right\\}$
(7)

Conveniently, if the moduli $|V_{k}|$ are Rayleigh-distributed, the term under
the second $k$-sum vanishes.

###  2.4 Numerical validation

A verification of the equations previously derived is a relatively quick task
with modern computing resources. We considered a power spectrum in the form
$P(k)\propto 1/k^{2\alpha+1}$ in the inertial range $[k_{min},k_{max}]$ and
zero outside of it. The slope is $\alpha=1/3$. We simulated the line profile
resulting from the projection through a (isothermal, isometallic)
$\beta$-model density profile (Cavaliere & Fusco-Femiano, 1978) with
$\beta=2/3$, core-radius $r_{c}$ and distance $\theta$ to the cluster centre:

$\epsilon(x)\propto\left(1+\frac{x^{2}+\theta^{2}}{r_{c}^{2}}\right)^{-3\beta}$

where $L=2$, $r_{c}=0.2$ and $\theta=0.2$. The moduli $V_{k}$ may be constant
(non-random) or follow a Rayleigh distribution. The Fourier coefficients of
the emissivity are given by ZuHone et al. (2016) (see also App. B):
$\widetilde{\epsilon}(k)/F=\exp(-\omega|k|c)(1+\omega|k|c)$ with
$c^{2}=r_{c}^{2}+\theta^{2}$. An example of such realisation is shown in Fig.
1. Only phases are random in this illustration. This configuration and random
realisation are specifically chosen to highlight the possible discrepancy
between the value of the broadening $S^{2}$ and the simple estimate
$\sigma_{turb}^{2}+\overline{\sigma_{th}^{2}}$. Indeed, purely by chance, most
of the points where velocity is high are located in low-emissivity regions,
therefore their contribution to line broadening is weak.

Figure 1: One realisation of a unidimensional turbulent velocity field (middle
panel) along the spatial axis $x$ (in arbitrary units) with parameters
$\alpha=1/3$, $k_{min}=1$ ($L_{inj}=2$), $k_{max}=20$ ($L_{diss}=0.1$),
$\sigma_{turb}=160$ km s-1 and $\sigma_{th}=100$ km s-1 (materialized by the
yellow shading). The emissivity profile (top panel) of the gas corresponds to
a $\beta$ density model with core-radius $r_{c}=0.2$ and at a distance
$\theta=0.2$ from the centre. The lower panel shows the resulting line profile
as a thick black line. The ”best-fit” Gaussian centred on $C$ (vertical line)
of width $S$ is shown as a dashed red line. The green thin curve shows the
Gaussian centred on the line centroid of width the geometrical mean of the
thermal and turbulent broadening.

Figure 2 shows the excellent agreement between the theoretical and simulated
quantities after 5000 random realisations of the velocity field. It is
important to notice the strong non-gaussianity of $S^{2}$ and $C$ in general.
The configuration chosen for this simulation clearly illustrates that $\langle
S^{2}\rangle<\overline{\sigma_{th}^{2}}+\sigma_{turb}^{2}$. This is a
consequence of the injection scale ($L_{inj}\sim L$) being much larger than
the typical cluster characteristic scale (here $r_{c}=L/10$), leading to non-
gaussian line shapes (Inogamov & Sunyaev, 2003).

|
---|---
Figure 2: The line centroid $C$ versus line width (squared) $S^{2}$ of an
emission line along a single line-of-sight and 5000 realisations of a
1-dimensional turbulent velocity field ($\sigma_{turb}=160$ km/s, $k_{min}=1$,
$k_{max}=20$, $\sigma_{th}=100$ km/s). Left panel assumes only random phases
while right panel also includes randomly distributed moduli (Rayleigh
distribution). Points show measurements, red cross is the measured mean and
standard deviations along both axes. The blue lines represent the results
obtained from analytical calculations (Eq. 4, 5, 6 and 7). The plain green
line shows the location of the geometric mean of the turbulent and thermal
dispersions: the presence of turbulent motions on scales comparable to that of
the cluster makes such estimate a biased one.

##  3 Two-dimensional characterisation of the velocity field

Observations and diagnostics of the intra-cluster medium rarely rely on single
line-of-sight measurements. Instead, due to instrumental resolution limits and
signal-to-noise considerations, line centroid shifts and broadening are
measured from spectra collected over well-defined 2-dimensional regions,
sometimes denoted as ”bins” or ”pixels”. A popular diagnostic tool in the
field of astrophysical turbulence (e.g. Lis et al., 1998; Esquivel et al.,
2007; Anorve-Zeferino, 2019) is the two-dimensional structure function,
loosely speaking a 2-d correlation function analysis of the line centroid
shift map. Obtaining analytical expressions of the sample variance of these
estimators requires the formalism above to be extended and to account for the
3-dimensional structure of the velocity field and the emissivity field. A
supplementary difficulty one has to face is the non-stationarity of the
projected velocity field: even though the 3-d velocity field is homogeneous
(stationary) in the medium, spatial variations of the emissivity in general
break this property.

###  3.1 The tridimensional velocity field

Similarly as in previous section, we define the centroid shift and line width
measured over a spectral line being a sum over all individual line-of-sights
selected within a region. We introduce the window function
$\mathcal{W}(\@vec{\theta})=\mathcal{W}(y,z)$ equals to 1 (one) for selected
line-of-sights and zero elsewhere. The measured spectral parameters of the
line write:

$\displaystyle F(\mathcal{W})$ $\displaystyle=\int I_{l}^{\mathcal{W}}(E){\rm
d}E$ $\displaystyle\delta E(\mathcal{W})$ $\displaystyle=\left(F^{-1}\int
E\times I_{l}^{\mathcal{W}}(E){\rm d}\@vec{\theta}{\rm d}E\right)-E_{0}$
$\displaystyle\Sigma^{2}(\mathcal{W})$ $\displaystyle=F^{-1}\int\left(E-\delta
E-E_{0}\right)^{2}I_{l}^{\mathcal{W}}(E){\rm d}E$

by defining:

$I_{l}^{\mathcal{W}}(E)=\int
I_{l}(\@vec{\theta},E)\mathcal{W}(\@vec{\theta}){\rm d}\@vec{\theta}$

All results from previous sections are obviously recovered with
$\mathcal{W}(\@vec{\theta})=\delta(\@vec{\theta}-\@vec{\theta}_{0})$.

Introducing $v(\@vec{\underline{x}})\equiv v(x,y,z)$ the line-of-sight
component of the velocity, and operating the following substitutions in Eq. 1:
$\epsilon(x)\rightarrow\epsilon(\@vec{\underline{x}})$, $v(x)\rightarrow
v(\@vec{\underline{x}})$, we can rewrite the observed ”aperture” flux,
centroid and velocity dispersion as:

$\displaystyle F_{\mathcal{W}}$ $\displaystyle=\int{\rm
d}\@vec{\theta}\,\mathcal{W}(\@vec{\theta})\int{\rm
d}x\,\epsilon(x,\@vec{\theta})=\int\mathcal{W}F{\rm d}\@vec{\theta}$
$\displaystyle C_{\mathcal{W}}$ $\displaystyle=F_{\mathcal{W}}^{-1}\int{\rm
d}\@vec{\theta}\,\mathcal{W}(\@vec{\theta})\int\epsilon(x,\@vec{\theta})v(x,\@vec{\theta}){\rm
d}x$ $\displaystyle=F_{\mathcal{W}}^{-1}\int\mathcal{W}FC{\rm d}\@vec{\theta}$
$\displaystyle S^{2}_{\mathcal{W}}$
$\displaystyle=F_{\mathcal{W}}^{-1}\int{\rm
d}\@vec{\underline{x}}\,\sigma_{th}^{2}(\@vec{\underline{x}})\mathcal{W}(\@vec{\theta})\epsilon(\@vec{\underline{x}})+\frac{F_{\mathcal{W}}^{-2}}{2}\int{\rm
d}\@vec{\underline{x}}{\rm
d}\@vec{\underline{x}^{\prime}}\,G(\@vec{\underline{x}},\@vec{\underline{x}^{\prime}})$
$\displaystyle=F_{\mathcal{W}}^{-1}\int\mathcal{W}F\left(S^{2}+C^{2}\right){\rm
d}\@vec{\theta}-C_{\mathcal{W}}^{2}$ $\displaystyle
G(\@vec{\underline{x}},\@vec{\underline{x}^{\prime}})$
$\displaystyle=\mathcal{W}(\@vec{\theta})\epsilon(\@vec{\underline{x}})\mathcal{W}(\@vec{\theta}^{\prime})\epsilon(\@vec{\underline{x}^{\prime}})\left[v(\@vec{\underline{x}})-v(\@vec{\underline{x}^{\prime}})\right]^{2}$

Integrations over $x$ range over an arbitrary large interval $[-L,L]$ and over
all possible values of the plane-of-sky position $\@vec{\theta}$. The velocity
can be written in terms of its Fourier decomposition with
$\@vec{\underline{k}}=(k_{x},k_{y},k_{z})=(k_{x},\@vec{\xi})$:

$v(\@vec{\underline{x}})=\sum_{\@vec{\underline{k}}}V_{\@vec{\underline{k}}}\exp\left(i\omega\@vec{\underline{k}}\cdot\@vec{\underline{x}}\right)$

and we note
$P_{3D}(\@vec{\underline{k}})=\langle|V_{\@vec{\underline{k}}}|^{2}\rangle=P_{3D}(k)$.
We have the following relations:

$\displaystyle V_{-\@vec{\underline{k}}}$
$\displaystyle=V_{\@vec{\underline{k}}}^{*}$ $\displaystyle\langle
V_{\@vec{\underline{k}}}V_{\@vec{\underline{k}^{\prime}}}\rangle$
$\displaystyle=\delta_{\@vec{\underline{k}};-\@vec{\underline{k}^{\prime}}}P_{3D}(k)$

Similarly as in the 1-dimensional case (Sect. 2), we introduce
$R_{\@vec{\underline{k}}}=P_{3D}(\@vec{\underline{k}})^{2}-Var\left(|V_{\@vec{\underline{k}}}|^{2}\right)$,
such that $R_{\@vec{\underline{k}}}=0$ for Rayleigh-distributed moduli. A
minimal assumption on the 4-term bracket $\langle
V_{\@vec{\underline{j}}}V_{\@vec{\underline{k}}}V_{\@vec{\underline{l}}}V_{\@vec{\underline{m}}}\rangle$
is necessary and our ansatz is explicitly provided in App. F.

###  3.2 Statistics of the aperture line centroid

We find that the average of the velocity shift measurements over several
realisations is 0:

$\langle C_{\mathcal{W}}\rangle=0$ (8)

The calculations are actually very similar to the one-dimensional case and we
refer to App. A for details. By using the Fourier decomposition of the
velocity field, the variance in centroid shifts measurements reads:

$\langle
C_{\mathcal{W}}^{2}\rangle=\frac{1}{F_{\mathcal{W}}^{2}}\sum_{\@vec{\underline{k}}}P_{3D}(k)|c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})|^{2}$
(9)

Here $c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})$ is the Fourier
coefficient of the product $\epsilon(\@vec{\underline{x}})\mathcal{W}(y,z)$.
This expression differs from Eq. (E7) in Zhuravleva et al. (2012) because we
do not assume $\epsilon$ being independent of the line-of-sight. If instead
$\epsilon(x,y,z)=\epsilon(x)$ in the domain of $\mathcal{W}\neq 0$, then
$c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})=\widetilde{\epsilon}(k_{x})\widehat{\mathcal{W}}(\@vec{\xi})$;
therefore we can rewrite our finding under the factorised form:

$\langle
C_{\mathcal{W}}^{2}\rangle=\frac{1}{F_{\mathcal{W}}^{2}}\sum_{\@vec{\underline{k}}}P_{3D}(k)P_{\mathcal{W}}(\@vec{\xi})P_{\epsilon}(k_{x})$

with $P_{\mathcal{W}}$ being the power-spectrum of the window function
$\mathcal{W}$. This expression is applicable considering for instance small,
pencil-beam, window functions or, equally interesting, narrow annular window
functions, if the emissivity shows a circular symmetry. We provide in App. B a
detailed calculation of the function $c_{\epsilon.\mathcal{W}}$ for the case
of the isothermal, isometallicity $\beta$-model gas density.

###  3.3 Statistics of the aperture line broadening

The calculation of the average of $S_{\mathcal{W}}^{2}$ over multiple
realisations of the turbulent field follows similar steps as in the
1-dimensional case presented before and we find:

$\langle
S_{\mathcal{W}}^{2}\rangle=\overline{\sigma_{th}^{2}}+\sigma_{turb}^{2}-F_{\mathcal{W}}^{-2}\sum_{\@vec{\underline{k}}}P_{3D}(k)|c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})|^{2}$

The bar indicates the average of the thermal broadening over the cluster
volume defined by $\mathcal{W}$. With these notations the relation found in
the 1-d case still holds:

$\langle S_{\mathcal{W}}^{2}\rangle+\langle
C_{\mathcal{W}}^{2}\rangle=\overline{\sigma_{th}^{2}}+\sigma_{turb}^{2}$ (10)

Finally the variance of the line broadening writes:

$\mathrm{Var}(S_{\mathcal{W}}^{2})=\\\
2\sum_{\@vec{\underline{k}},\@vec{\underline{k}^{\prime}}}P_{3D}(\@vec{\underline{k}})P_{3D}(\@vec{\underline{k}^{\prime}})\left|\frac{c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}}+\@vec{\underline{k}^{\prime}})}{F_{\mathcal{W}}}-\frac{c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}^{\prime}})}{F_{\mathcal{W}}^{2}}\right|^{2}\\\
-\sum_{\@vec{\underline{k}}}R_{\@vec{\underline{k}}}\times\left\\{\left|\frac{c_{\epsilon.\mathcal{W}}(2\@vec{\underline{k}})}{F_{\mathcal{W}}}-\frac{c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})^{2}}{F_{\mathcal{W}}^{2}}\right|^{2}\right.\\\
\left.+2\left[1-\frac{|c_{\epsilon.\mathcal{W}}(\@vec{\underline{k}})|^{2}}{F_{\mathcal{W}}^{2}}\right]^{2}\right\\}$
(11)

which reduces to the first term only in case of Rayleigh-distributed moduli.

###  3.4 Statistics of the structure function

We define the structure function as the integral:

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

This expression simply describes an average over all pairs of regions (called
’bins’ or ’pixels’) $(\mathcal{W},\mathcal{W}^{\prime})$ separated by a
distance111There is quite a latitude in choosing the definition of distance,
either considering geometrical centres of each region or flux-weighted
barycentres, etc. $d(\mathcal{W},\mathcal{W}^{\prime})=s$. Here $N_{p}(s)$ is
the number of such pairs of regions. For instance, considering single line-of-
sight measurements and the Euclidian distance between two points on sky, i.e.
$\mathcal{W}(\@vec{\theta})\equiv\delta(\@vec{\theta}-\@vec{\theta}_{0})$, we
recover the standard formulation (e.g. ZuHone et al., 2016):

$SF(s)=\frac{1}{N_{p}(s)}\int_{\@vec{\theta}_{0},|\@vec{r}|=s}\left|C(\@vec{\theta}_{0}+\@vec{r})-C(\@vec{\theta}_{0})\right|^{2}{\rm
d}N_{p}$

The integration runs over an arbitrary large, but bounded region of sky
$\mathcal{A}$ of total area $\mathcal{S}_{\mathcal{A}}$. In the following we
consider $\mathcal{A}(\@vec{\theta})$ as a function taking value 1 in the
analysis region and 0 outside. There, $N_{p}(s)$ needs to be interpreted as
the integral $\int{\rm d}N_{p}$ for all $(\@vec{\theta}_{0},|\@vec{r}|=s$).

Such defined, $SF(s)$ is a random variable that depends on the particular
realisation of the velocity field and we can therefore compute its mean and
variance across several realisations, hereafter called $\mathrm{sf}(s)$ and
$\sigma^{2}_{\mathrm{sf}}(s)$.

####  3.4.1 Expected value $\mathrm{sf}(s)$

Under the assumption that finite-size (i.e border) effects are negligible, we
find for the most general expression of the emissivity field (see App. D) :

$\displaystyle\mathrm{sf}(s)$ $\displaystyle=$ $\displaystyle
2K\sum_{\@vec{\underline{k}}}P_{3D}(k)\int{\rm
d}\@vec{\xi}^{\prime}P_{\rho}(k_{x},\@vec{\xi}^{\prime})\left[1-J_{0}\left(\left|\@vec{\xi}+\@vec{\xi}^{\prime}\right|\omega
s\right)\right]$
$\displaystyle=2\left(\frac{\omega}{2\pi}\right)^{2}\int\left[1-J_{0}\left(\omega\left|\@vec{\xi}\right|s\right)\right]P_{2D}(\@vec{\xi}){\rm
d}\@vec{\xi}$ (12)

with $\rho(x,y,z)=\epsilon(x,\@vec{\theta})/F(\@vec{\theta})$,
$K=\omega^{2}/(4\pi^{2}\mathcal{S}_{\mathcal{A}})$, $J_{0}$ being the Bessel
function of the first kind and order 0. The function $P_{2D}$ is the 2D power
spectrum of the centroid shift map, which expressions are properly defined and
derived in App. C. One must be careful that the power-spectrum $P_{\rho}$
involved in these expressions is that of the normalised emissivity field
$\rho$, i.e. the 3-d emissivity $\epsilon$ divided by the ”flux map”
$F(\@vec{\theta})$. It is strongly dependent on the choice of analysis domain
$\mathcal{A}$.

In the special case where the two-dimensional spectrum is isotropic this
expression takes the following form:

$\mathrm{sf}(s)\simeq
4\pi\left(\frac{\omega}{2\pi}\right)^{2}\int_{0}^{+\infty}P_{2D}(\xi)\left(1-J_{0}(\omega\xi
s)\right)\xi{\rm d}\xi$

The latter equation resembles ZuHone et al. (2016), their eq. 29. However this
result implicitly includes the shape of the domain of analysis through
$P_{2D}$, which effectively acts as a high-pass spatial filter. We recall here
the assumptions leading to this result: i) centroid shift measurements are
performed along individual line-of-sights, ii) isotropy of the two-dimensional
power-spectrum $P_{2D}$, iii) the averaging domain allows all possible
orientations of the pair vector $\@vec{r}$ and iv) the sum over modes
$\@vec{\xi}$ can be written as an integral.

We provide in App. D a generic formula to correct the above expressions for
border effects. This involves calculation of the number of pairs enclosed
within the analysis domain and those crossing its frontier, both dependent on
the separation length $s$ and the exact shape of the domain. These are easily
calculated for a circular domain of analysis of radius $R$ and we provide the
equations in the appendix.

We also provide in App. E a prescription to account for pixelization of the
centroid map with pixels of arbitrary size and shapes. Provided such pixels
are small with respect to the typical scales of the surface brightness
fluctuations, a correction is obtained by multiplying $P_{2D}$ by the two-
dimensional power-spectrum of the pixel shape $P_{\ell}$. As shown in
appendix, this prescription should not be used in combination with the
correction formula for border effects, especially if pixels are of sizeable
length compared to the analysis domain. We do not provide here a complete
analytical formulation accounting simultaneously for border effects and
pixelization; it may be more advantageous in such case to numerically estimate
the average structure function from its primary definition involving the
$C_{\mathcal{W}}$’s.

Nevertheless, an exact solution for $\mathrm{sf}(s)$ is obtained in case of a
stationary 2-dimensional velocity field – e.g. if
$\epsilon(x,y,z)=\epsilon(x)$ – by replacing $P_{2D}$ by $P_{2D}^{\infty}$ in
Eq. 3.4.1, that is the power spectrum computed in the limit of an infinitely
extended analysis domain (see App. C for details.) Such formulation then
matches exactly that proposed by ZuHone et al. (2016). The above prescription
for pixel binning then also becomes exact and raises no issue due to a finite
region of analysis. These properties are used in App. D and E to validate our
correction formulas and to stress their limitations.

####  3.4.2 Variance $\sigma^{2}_{\mathrm{sf}}(s)$

A full calculation of the covariance

$\Sigma_{ij}=\mathrm{Cov}\left(SF(s_{i}),SF(s_{j})\right)$

between structure functions measured at different scales is provided in App. F
under the assumption of negligible finite-size effects. The complete formula
is given in Eq. 35 and involves integrals of the Fourier transform
$\widetilde{\rho}$ of the normalised emissivity field. Because of the relative
position of the analysis region and the emissivity distribution, it is in
general not possible to factor their respective contributions in the
expression of the variance. However, we can study a simpler, practical case
where the emissivity is independent on the line-of-sight direction within the
given analysis field-of-view. This is for instance the case for an observation
pointing at the outskirts of a nearby galaxy cluster or towards the core of a
”flat” galaxy cluster (e.g. Coma). This particular case writes
$\rho(x,\@vec{\theta})=\epsilon(x)/F$. This leads to decoupling the
calculation of ”geometrical” terms (i.e. the shape and location of the
instrumental field-of-view) and ”fluctuation” terms (the coupling between the
cluster emissivity and the turbulent velocity spectrum). In App. 35 we obtain
the following simple formula, under the supplementary hypothesis of a very
large analysis region, i.e. for
$\mathcal{S}_{\mathcal{A}}^{1/2}\gg(s,L_{inj},...$):

$\Sigma_{ij}\simeq
16\pi\left(\frac{\omega}{2\pi}\right)^{2}\int\left[\frac{1}{\mathcal{S}_{\mathcal{A}}}P_{2D}^{\infty}(\xi)^{2}-\left(\frac{\omega}{2\pi}\right)^{2}Q_{2D}^{\infty}(\xi)^{2}\right]\\\
\times\left(1-J_{0}(\omega\xi s_{i})\right)\left(1-J_{0}(\omega\xi
s_{j})\right)\xi{\rm d}\xi$ (13)

where $Q_{2D}^{\infty}=0$ for Rayleigh-distributed moduli. The diagonal term
of this quantity is then: $\Sigma_{ii}=\sigma^{2}_{\mathrm{sf}}(s)$. Similarly
as for the calculation of the average $\mathrm{sf}(s)$ in case of pixelized
data, one has to multiply $P^{\infty}_{2D}$ by the power-spectrum of the
elementary pixel shape. Eq. 13 in the Rayleigh regime is the expression we
will validate in the next section, keeping in mind the series of assumption
made to obtain this simple formulation.

##  4 Numerical validation

We performed a set of numerical experiments to validate the equations derived
previously. This requires generation of multiple velocity boxes in three
dimensions. Given the high computational demand, only a selected set of cases
are treated.

###  4.1 Dataset of velocity cubes

We created a series of velocity boxes with characteristics indicated in Table
1. The smooth and continuous velocity power-spectrum takes a form similar to
that of ZuHone et al. (2016), namely:

$P_{3D}(\@vec{\underline{k}})=C_{n}e^{-(k/k_{diss})^{2}}k^{\alpha}e^{-(k_{inj}/k)^{2}}$
(14)

with $C_{n}$ a normalization constant having units such that $\int
P_{3D}(\@vec{\underline{k}}){\rm d}\@vec{\underline{k}}=\sigma_{turb}^{2}$ and
$\alpha=-11/3$ typical of a Kolmogorov turbulence spectrum (Kolmogorov, 1941).
Scales $k_{diss}=1/L_{diss}$ and $k_{inj}=1/L_{inj}$ represent the dissipation
and injection frequencies, respectively. Moduli of the Fourier coefficients
are drawn from a Rayleigh distribution.

3-d box size | Inject. scale | Dissip. scale | Slope ($\alpha$) | $C_{n}$ | $\sigma_{turb}$ | N. realisations
---|---|---|---|---|---|---
(pixel)3 | Mpc$\times$(kpc)2 | (kpc) | (kpc) |  | (*) | (km/s) |
$1936\times 242^{2}$ | $4.2\times 520^{2}$ | 100 | 10 | $-11/3$ | 807.9 | 448.3 | 100
$1936\times 242^{2}$ | $4.2\times 520^{2}$ | 200 | 10 | $-11/3$ | 428.8 | 443.7 | 100
$1936\times 242^{2}$ | $4.2\times 520^{2}$ | 300 | 10 | $-11/3$ | 307.1 | 442.3 | 100
Table 1: Numerical realisations of a 3-dimensional velocity cube used for
validating the analytic calculations of line centroid shift, line broadening
and structure function sample variances. The size of the box in the line-of
sight direction is 8 times larger than the transverse (plane-of-sky) box size.
(*: units km2 s-2 kpcα+3.)

The computationally demanding Fast Fourier transforms (FFT) were distributed
across 10 processors using 2DECOMP&FFT222http://www.2decomp.org (Li & Laizet,
2010). The histograms of the (3-d) velocity standard deviation in each of the
three configurations is shown in Fig. 3, this is an indicator of the goodness
of the simulated field. Clearly, as the injection scale increases the box
becomes too small for the periodic boundary condition to apply during the FFT.
The third panel indicates an additional ”noise” of order 5-10 km/s in the
$L_{inj}=300$ kpc run, which we attribute to aliasing effects. This extra
numerical scatter needs to be reminded while comparing analytic results to
simulations.

Figure 3: The 3-d velocity dispersion $\Delta V_{3D}$ in each of the 100
numerical realisations for the three configurations in Table 1 is shown as
blue histograms. Vertical dashed line indicates the exact value of
$\sigma_{turb}$ from integration of the input turbulent power-spectrum (Eq.
14). The increased numerical dispersions in those values as $L_{inj}$
increases, as a consequence of the finite simulation box size.

The emissivity of the galaxy cluster gas is taken as the square of a
(isothermal, isometallic) gas density considered either as a spherical
$\beta$-model (hereafter _beta_) or as a $\beta$-model along the line-of-sight
and constant over the plane of the sky (hereafter _Xbeta_). The $\beta$
parameter is held at a value $2/3$ while the core-radius takes value in
$\\{4,21,54,107,215,429\\}$ kpc. The normalization of the emissivity plays no
objective role in this study, since no signal-to-noise consideration is made.
Figure 4 shows one example of the line centroid and line width maps for a 200
kpc injection scale and the various _beta_ emissivity models, free of any
uncertainty other than numerical noise. In all numerical simulations there is
no thermal broadening ($\sigma_{th}=0$ km s-1 hereafter). In the following we
present the comparisons with the $L_{inj}=100$ kpc simulation only. Validation
of the other two runs is extensively presented in App. H for completeness.

Figure 4: Projection of a single realisation of a 3-d velocity field
(injection scale at 200 kpc) with several emissivity models (top row). All but
the last column correspond to spherical $\beta$-models with core-radii 4, 21,
54, 107, 215 and 429 kpc (from left to right). The rightmost column
corresponds to a constant emissivity in the entire simulation box. The size of
each panel is 520 kpc on a side. Middle row shows the centroid shift ($C$) and
bottom row shows the line width ($\sqrt{S^{2}}$). Particularly noticeable is
the decrease in contrast (or power) as the core-radius increase and the small
line broadening seen through a small cluster core (e.g. bottom-left figure).

###  4.2 Centroid and line broadening

We first carry out the validation of Equations 8, 9, 10 and 11, which provide
analytical representations of the sample average and variance of the line
centroid shift and line broadening (more specifically, the square of the line
width) measured in arbitrary apertures. We limit this validation exercise to
circular apertures centred on a galaxy cluster and allow their sizes to vary.

These analytical expressions involve calculation of the 3-d function
$c_{\epsilon.\mathcal{W}}$: we provide in App. B the analytical formulas for
both considered emissivity models and for circular apertures. Calculation of
this function for spherical $\beta$-models demands slightly more computing
time than for the _Xbeta_ model.

Equation 11 requires integration over 6 scalar variables. Taking advantage of
the isotropy of the velocity power spectrum and the 2-d rotational invariance
of this specific configuration, this can be reduced to five integration
variables only (e.g. $k_{x},k_{x}^{\prime},\xi,\xi^{\prime}$ and one angle
$\phi$). This integral is evaluated by Monte-Carlo sampling distributed over
40 computing cores by means of the MCQUAD library333Available in package
SciKit-Monaco, https://pypi.org/project/scikit-monaco/. The number of
samplings is $2.10^{6}$ and $2.10^{5}$ for the _Xbeta_ and _beta_ emissivity
models respectively and we monitor and store the statistical uncertainties out
of the numerical sampler.

Figures 5 and 6 show the comparison between analytical calculations (plain
lines) and numerical simulations (dots with error bars). They correspond to
the _Xbeta_ and _beta_ emissivity models respectively, using the same 100
velocity boxes with $L_{inj}=100$ kpc. Each dot corresponds to a calculation
using the 100 velocity realisations and a given core-radius size and a given
aperture size. Error bars are derived from bootstrap resampling. Because we
always used the same 100 simulations, the deviations to the expected trend
appear correlated: this is for instance striking in the left-most panel
showing $\langle C\rangle$. This behaviour is likely to disappear with a
higher number of realisations.

Figure 5: Numerical validation of equations 8, 9, 10 and 11, i.e. the expected
value and sample variance of the line centroid shift $C_{\mathcal{W}}$ (first
and second panel) and the expected value and sample variance of the line width
$\sqrt{S^{2}_{\mathcal{W}}}$ (third and fourth panel). Plain lines show the
analytical calculations, data points are measured on 11 numerical realisations
of a turbulent field (errors estimated via bootstrap) with $L_{inj}=100$ kpc.
The calculations are performed assuming measurements in circular apertures
$\mathcal{W}$ of various radii (x-axis). The emissivity model is _Xbeta_ with
core-radii indicated in legend. The uncertainty on the analytical results for
$\sigma(S^{2})$ (materialized by the line widths in the last panel) is due to
limitations of the numerical integrator used to evaluate Eq. 11. Figure 6:
Similar figure as Fig. 5 for a spherical $\beta$-model emissivity (_beta_).
The numerical uncertainties are slightly larger (in the 4th panel, compared to
Fig. 5) due to a lower accuracy in the numerical integration of Eq. 11.

In any case these figures demonstrate a very good agreement between analytical
calculations and simulations. The only exception is the case of very large
core-radii ($r_{c}=429$ kpc), for both emissivity models. This is a
consequence of the simulation box being too small in the x-direction (4240 kpc
along the line-of-sight). This causes a non-negligible sharp cut-off in the
simulated $\beta$-emissivity profile, not accounted for by the analytical
equations.

The sample variance of the centroid shift and the line width exhibit large
variations with respect to the aperture radius. Emission line diagnostics in
growing apertures for a selection of ’look-alike’ galaxy clusters has
interesting potential to reveal the properties of the underlying turbulent
power-spectrum. This is illustrated in Fig. 7 where we vary the injection
scale from 100 to 300 kpc for a given $\beta$-model ($r_{c}=107$ kpc). It is
out of the scope of this paper to provide forecasts on the constraining power
of this method, which must also include measurement uncertainties and
limitations related to the availability of samples.

Figure 7: The ”sample variance” associated to line measurements – centroid
shift (left) and broadening (right) – in apertures of growing sizes for
clusters presenting identical emissivity models (a spherical $\beta$-model
with core-radius 107 kpc.) These curves are predicted analytically by Eq. 9
and 11. They have been normalized to the value of the 3-d turbulent velocity
dispersion $\sigma_{turb}$. The different shapes of the curves as the aperture
radius is growing can be used as a diagnostic to discriminate between various
injection scales.

###  4.3 Structure function

We restrict the numerical validation to that of Eq. 3.4.1 and 13 for
computational reasons. First, we consider an emissivity model of type _Xbeta_
and a large enough analysis region compared to the typical separation and
pixelization of the line centroid map. The same $3\times 100$ simulated boxes
are projected and pixelized in square regions of size $\ell\times\ell$ where
$\ell=4,9,17,34,69$ kpc. Since the two-dimensional velocity field is
stationary, it is fine to use $P_{\ell}P_{2D}^{\infty}$. This replacement is
justified because in this configuration the surface brightness is constant
over the analysis domain. In what follows the analysis domain is a circular
aperture of diameter 520 kpc. Figure 18 illustrates how pixelization acts on a
simulated centroid shift map: it indeed is very close to a convolution or
’smoothing’. Similar maps are created for all pixel sizes and core-radii for
all simulated boxes. The geometrical centres of the pixels are used to compute
the structure function, as the arithmetic mean of the squared centroid
gradients at pre-defined separations $s$ (within a range $\delta s$). The
average of the structure functions and the standard deviations at each $s$
provide the numerical indicators to be compared to Eq. 3.4.1 and 13.

Analytical calculation of $P_{2D}$ is performed according to Eq. 26, by
2-dimensional convolution444Making use of the FFT convolution implemented in
the signal.convolve function of Numpy/Scipy. of the 3-d velocity power-
spectrum $P_{3D}$ and the function $P_{\rho}$ at each frequency $k_{x}$ and
eventually summing over those frequencies. The whole procedure is distributed
over 40 processors working in parallel. Analytical expressions for $P_{\rho}$
are given in App. G (particularly Eq. 37) for the emissivity models relevant
to our validation procedure. The calculation of $P_{2D}^{\infty}$ is much more
straightforward, see Eq. 27.

A comparison between the analytical and numerical results is illustrated in
Fig. 8 for a given turbulent power spectrum (injection scale at 100 kpc) and
various values for the core-radius and pixel size. The results from the 100
realisations are displayed as thin grey lines and their distribution at each
$s$ is likely not Gaussian.

The analytical and numerical values for the sample variance are in very good
agreement for all 9 configurations, which is a very encouraging result given
the various assumptions involved in both cases.

Figure 8: Comparison of numerical and analytical structure functions and their
sample variance for various cluster sizes (i.e. various core-radii ($r_{c}$)
of the $\beta$-model) and various pixel sizes ($\ell$). The data points and
thick error bars show the sample mean and standard deviation of the 100
realisations (individually represented as thin grey lines) for each of the
considered configurations. The coloured curves and shaded areas represent the
analytical calculations following Eqs. 3.4.1 and 13. The emissivity model is
_Xbeta_ and the region of analysis is a circle of diameter 520 kpc (as
diplayed in Fig. 18). The turbulent power-spectrum is that of Table 1 with
injection scale 100 kpc.

A thorough assessment of the agreement between the analytical calculations and
the numerical validation is summarised in Figures 9 and 10. They show the
relative difference (expressed in percent) between the analytical and the
numerical computations for the expected value of the structure function and
its variance respectively. The injection scale is $L_{inj}=100$ kpc, as in
Fig. 8. Three separations are illustrated: $s=20$ kpc (close to the
dissipation scale), $s=60$ kpc (within the inertial range) and $s=300$ kpc
(past the inertial range).

Regarding the sample mean, the analytical model (Eq. 3.4.1) performs well
within 20 % of the numerical experiment. Keeping in mind the limited number of
realisations (100 samples), this result appears satisfactory. A degradation of
the prediction accuracy arises as the binning size increases. This is
attributed to numerical approximation in computing $P_{\ell}$ and higher
levels of sampling noise in the simulation (larger pixels imply fewer
$s$-pairs). The slight decrease in accuracy at larger core-radii was already
pointed out in Sect. 4.2, as a result of the simulation box size.

Figure 9: Representation of the absolute relative difference between the
numerical and analytical estimates of the sample mean of the structure
function, $\langle SF\rangle$, at three distinct separations $s$. Each
coloured square corresponds to one experiment based on the same 100
realisations of the velocity field with $L_{inj}=100$ kpc and various binning
sizes (y-axis) and $\beta$-models core-radii (x-axis). Small red crosses
indicate locations where the binning size is larger than $s$.

As for the sample variance (Fig. 10), the relative differences between
analytical and numerical results show somewhat higher values, as is expected
for second order statistics. In general our formula tends to overpredict by a
few tens of percent the observed variance of the structure function at small
separations ($s=20$ kpc) as a result of numerical uncertainties both in the
simulations and the evaluation of integrals. At large separations ($s=300$
kpc) and for the largest pixel size, the analytic formula underpredicts the
variance by up to 80%. The analysis region $\mathcal{A}$ indeed cannot be
considered as infinitely large any longer, making simplification of Eq. 35
into Eq. 13 less accurate.

Figure 10: Similar as Fig. 9, but for the sample variance of the structure
function, $\mathrm{Var}(SF)$. Positive values indicate higher predicted
variance compared to that measured in the numerical validation procedure.
Although some of these numbers are high at face value, it is important to
recall the assumptions leading to the chosen analytical formula and the noise
inherent in our set of numerical simulations (see text).

We finally relax the assumption of a constant emissivity and we show in Fig.
11 a comparison of the structure functions obtained for a spherical
$\beta$-model density (_beta_ emissivity model). The analytical formula Eq.
3.4.1 recovers the mean structure function, despite the spatial non-
stationarity of the projected velocity field. We corrected for border effects
using Eq. 30, the region analysis being a circle of diameter 520 kpc. As
highlighted in App. E, we are not able to use the simple prescription for
significantly large pixel binnings. Moreover we did not carry the full
evaluation of the sample variance using Eq. 35 for this figure. Instead, we
computed the variance according to Eq. 13 assuming an effective core radius
$c=\sqrt{r_{c}^{2}+\theta_{\mathrm{eff}}}$. Such approximation of the complex
emissivity field by an effective emissivity extracted at a radius
$\theta_{\mathrm{eff}}=80$ kpc from the cluster centre makes the calculation
more tractable. It shows a good agreement with results obtained from the
Monte-Carlo simulations.

Figure 11: Comparison of the simulated and calculated structure function in a
similar way as Fig. 8, except the emissivity model is of type _beta_
(spherical $\beta$-model gas density). This induces non-stationarity of the
projected velocity field, noticeable by the drop at large $s$. The coloured
curves represent the analytical calculation of the mean structure function
following Eqs. 3.4.1 and 30. For simplicity, the variance remains calculated
according to Eq. 13, i.e. assuming _Xbeta_ emissivity with an effective core-
radius $c=\sqrt{r_{c}^{2}+\theta_{\mathrm{eff}}^{2}}$ with
$\theta_{\mathrm{eff}}=80$ kpc.

##  5 Discussion

This work provides an extension of earlier studies, among others Inogamov &
Sunyaev (2003), Churazov et al. (2012), Zhuravleva et al. (2012) and ZuHone et
al. (2016). We extended the formalism presented in these papers by: i)
addressing the case of an arbitrary emissivity field; ii) computing second
order statistics beyond expected values (i.e. the sample variance) and iii)
identifying limiting cases in which these studies coincide.

###  5.1 Validity of hypotheses and range of applicability

Our study remain formal and rely on strong hypotheses such as uniform, ergodic
and isotropic turbulent velocity fields throughout the intra-cluster medium,
whose physics is encapsulated in a universal Kolmogorov power-spectrum,
meaning that turbulence follows an identical physical description from cluster
to cluster. The latter assumption is often implicitly made and mirrors an
intent to concentrate all the unknown physics of turbulence in a single
mathematical description. It is clear that this hypothesis may fail if widely
distinct mechanisms produce turbulent motions: for instance large-scale matter
accretion and central AGN feedback.

More specifically, our assumption of isotropic turbulent motions may break
down in the stratified intra-cluster medium where buoyancy-restoring forces
tend to suppress motions along the radial direction. Numerical simulations
indicate a change in the morphology of turbulent fields in regions showing
strong density gradients (e.g. Shi & Zhang, 2019), even in cluster cores
(Valdarnini, 2019). According to these findings, one can therefore expect our
model to become less representative as larger and larger cluster radii enter
the emission line analysis, or in presence of strong cool-core clusters. A
study of anisotropic turbulence is out of the scope of this paper: for
instance, one would undertake similar derivation steps as shown in appendices,
dropping the assumption $P_{3D}(\@vec{\underline{k}})=P_{3D}(k)$.

The existence of several (two) drivers of turbulence acting at different
scales may change the shape of the velocity power-spectrum and more generally
the statistical relations between Fourier coefficients of the velocity field.
ZuHone et al. (2016) proposed to rewrite the resulting $P_{3D}$ as a sum of
two Kolmogorov-like power-spectra with different injection scales, based on
the simulations and results of Yoo & Cho (2014). Such a prescription enters
the framework presented in this paper, because our results are independent on
the exact shape of $P_{3D}$. However, it remains to be checked whether the
decomposition of $\langle
V_{\@vec{\underline{j}}}V_{\@vec{\underline{k}}}V_{\@vec{\underline{l}}}V_{\@vec{\underline{m}}}\rangle$
proposed in App. F holds under such conditions, which most likely can be
addressed through numerical simulations.

Our hypothesis also assumes full decoupling of the gas emissivity and the
local behaviour of turbulent motions. This simplifying assumption may largely
fail if gas motions are induced by merging of an external galaxy group which
shows high emissivity in its vicinity. An interesting perspective of the
present calculations would be the coupling between $P(k)$, the turbulent power
spectrum and $\epsilon$, the emissivity; however it is likely that
calculations would become more complex and the gain over Monte-Carlo
simulations would become less obvious. Moreover, density fluctuations, hence
emissivity fluctuations, are thought to be directly linked to the turbulent
power spectrum based on theoretical grounds (Churazov et al., 2012). At first
order though, the broad-scale emissivity of the galaxy cluster gas is the
dominant component modulating the Doppler shift in the integrated line
profiles and our approach remains a reasonable one in this regard.

Finally, our work deliberately neglects measurement uncertainties and
instrumental noises. We address this assumption in a subsequent study (paper
II, Cucchetti et al., in press) by propagating the impact of measurement
uncertainties on the line diagnostics, in particular the structure function.
An interesting conclusion of this study is that sample variance effects
dominate on large scales the error budget for observations based on next-
generation X-ray instruments such as _Athena_ /X-IFU, while statistics
dominate at small scales.

###  5.2 An application: forecasting line shift and width profiles

The formulas derived in Sect. 3 provide the sample mean and variance of both
the line centroid shift and width in arbitrary apertures. As such, they can be
used to predict measurements in concentric annuli centred on a galaxy cluster,
i.e. a radial profile. Formally, an annular aperture mask is defined as the
difference between two concentric circular apertures. Thanks to the linear
behaviour of the Fourier transform, the coefficient $c_{\epsilon.\mathcal{W}}$
is also the difference between the two corresponding coefficients, both easily
computed following App. B. Interestingly, the emissivity in each annulus can
be considered as constant, i.e. $\epsilon(x,\@vec{\theta})=\epsilon(x)$, if
the gas density shows spherical symmetry and the annuli are thin enough. As
already noted, this property drastically reduces the computing time needed to
integrate the equations. An example of the profiles of centroid shift
variance, line width average and line width variance are displayed on Fig. 12
for a turbulent power spectrum with injection scale 100 kpc and
$\sigma_{turb}=448$ km s-1. No thermal broadening is included in this
exercise. As expected, the centroid shift (whose average value is zero) shows
larger variance in the central bins than the outskirts and the larger the core
radius, the smaller the effect. The average broadening shows the reverse
behaviour with smaller widths in the central parts and reaching a plateau
(corresponding to $\sigma_{turb}$) in the outskirts. The line width variance
shows diverse behaviours but here again, the general trend is a decrease
towards the outskirts.

Figure 12: Model predictions for radial profiles of line properties, i.e
measurements in spectra collected in circularly concentric annuli of equal
width (21 kpc). Shown are the sample variance of the centroid shift (left
panel), the sample average of the line broadening (middle panel) and the
sample variance of the line broadening. The injection scale is $L_{inj}=100$
kpc, the emissivity model is a spherical $\beta=2/3$ model with core radii
indicated in legend. Shaded rectangles indicated bin widths (horizontally) and
numerical uncertainties (vertically).

###  5.3 Forecasting the structure function from upcoming instrumentation

One particularly interesting perspective consists in inverting the formulas
presented here to evaluate the power of future astronomical X-ray micro-
calorimeters in constraining the nature of turbulent motions in galaxy
clusters. By properly selecting the samples (typically, the number of objects
and their core-radii and distances) and the observing strategy (mapping,
exposure times, etc.) one is able to focus the constraints on, e.g. the slope
of the power-spectrum or the injection scale. This assumes that turbulence has
identical characteristics throughout the sample considered, which hopefully is
a reasonable guess. We postpone the complete exercise to later investigation.
Rather we compute the expected structure functions for a simplified Coma-like
galaxy cluster, following a setup similar to ZuHone et al. (2016) and the
associated uncertainties due to sample variance only (statistical errors are
disregarded). We consider two instruments: i)  _XRISM_ /Resolve with a
resolution element of $1.5\arcmin$ and a field-of-view of $3.4\arcmin$
equivalent diameter and ii)  _Athena_ /X-IFU with a resolution element of
$5\arcsec$ and a field-of-view $5\arcmin$ equivalent diameter (Barret et al.,
2018). We consider two observing strategies: either one single pointing
towards the cluster centre, or the mapping of a $\sim 15\arcmin\times
15\arcmin$ area with multiple pointings. We also consider the case of a
$15\arcsec$ pixelization rebinned images for _Athena_ /X-IFU, such that the
signal-to-noise ratio of each spectrum is increased (e.g. Roncarelli et al.,
2018, also Cucchetti et al., in press). At the redshift of Coma, $1\arcmin$ on
sky corresponds roughly to 27 kpc physical separation. We consider a turbulent
power spectrum with $\sigma_{turb}=438$ km/s, $\alpha=-11/3$, injection scale
at 200 kpc and dissipation scale at 20 kpc. Given the proximity of Coma and
its apparent size, using the _Xbeta_ emissivity model is amply justified, as
already noted by Churazov et al. (2012); ZuHone et al. (2016). Figure 13 shows
the result outcome of our model. For identical sky coverages, X-IFU provides
smaller relative variance in comparison to Resolve, thanks to its better
angular resolution. Even in one single pointing X-IFU can provide a
measurement of the structure function up to $\sim 100$ kpc separation scales.
The associated variance is larger though, due to a smaller number of pairs
entering the structure function.

This example provides the basis in view of optimising an observational
strategy for a given instrumental setup. Our formalism involves Fourier
transforms of window functions (denoted $\mathcal{A}$ and $\mathcal{W}$) and
therefore accounts for arbitrary instrumental shapes and pointing strategies,
by taking advantage of standard properties of the Fourier transform. For
instance, a window function made of multiple non-overlapping pointings can be
considered as a sum of identical, translated window functions; linearity then
makes the computation of its Fourier transform straightforward.

Figure 13: Model predictions for structure functions and their associated
sample variances under two instrumental setups: _XRISM_ /Resolve (assuming
$1.5\arcmin$ resolution elements) and _Athena_ /X-IFU (assuming $5\arcsec$ and
$15\arcsec$ resolution elements for high and low signal to noise ratios
respectively). The left panel shows predictions for a $\sim 15\arcmin\times
15\arcmin$ contiguous mapping of the Coma cluster while the right panel shows
the result for a single X-IFU pointing. A single Resolve pointing ($3\arcmin$
on a side) would be too small for a useful derivation of the structure
function. See text for details on the input turbulent power spectrum and gas
density model.

##  6 Conclusions

In this paper we have derived analytical expressions for the sample mean and
variance of three indicators of turbulence in X-ray emitting, optically-thin,
plasmas under the hypothesis of homogeneous and isotropic Kolmogorov
turbulence. These are the line centroid shift $C$, the line broadening $S$ and
the structure function $SF$.

  1. 1.

We obtained exact expressions for the mean and variance of $C$ and $S$
obtained from single line-of-sight measurements through arbitrary gas
emissivity: Eqs. 4, 5, 6 and 7. We numerically validated the results with
Monte-Carlo simulations of turbulent velocities with Gaussian or constant
amplitudes.

  2. 2.

We generalised these expressions for measurements in apertures of arbitrary
shapes and sizes and for arbitrary 3-dimensional emissivity fields: Eqs. 8, 9,
10 and 11. We provided in App. B useful formulas for the common $\beta$-model
and for circular apertures. We numerically validated the formulas using Monte-
Carlo simulations of 3-d velocity fields in a range of emissivity and power-
spectrum configurations.

  3. 3.

We derived an expression for the mean structure function under the assumption
of negligible border effects (Eq. 3.4.1). Notably, this formula does not
assume constant (’flat’) emissivity in the plane-of-sky direction. It involves
a specific definition for the two-dimensional power-spectrum of the projected
velocity field, introduced in App. C. We provided in App. G useful formulas
for the common $\beta$-model and for circular domains of analysis.

  4. 4.

In App. D we provided a correction formula for border effects (Eq. 30) valid
for non-binned maps of the projected velocity and for domains of arbitrary
shapes. We explicitly computed the case of a circular field-of-view.

  5. 5.

We provided in App. E a simple prescription to account for binning (or
pixelisation) on the mean structure function. It is valid as long as pixels
are smaller than the typical scale of flux variations and much smaller than
the domain of analysis.

  6. 6.

We derived a fairly generic expression for the sample variance of $SF$ under
assumption of negligible border effects and for arbitrary emissivity fields
(Eq. 35). This equation takes a tractable form in case of flat emissivity and
very large domain of analysis (Eq. 13).

  7. 7.

We numerically validated our results for the sample mean and variance of $SF$
in the case of ’flat’ emissivity fields ($\beta$-models with a range of core
radii) and various binnings.

  8. 8.

We numerically validated our results for the sample mean of $SF$ in the case
of non-flat emissivity fields (spherical $\beta$-models with a range of core
radii) and negligible binning.

  9. 9.

We discussed our results and presented forecasts for observations of the core
of the Coma cluster with the integral field units X-ray calorimeters planned
to embark onboard _XRISM_ (Resolve) and _Athena_ (X-IFU).

###### Acknowledgements.

The authors thank the referee for useful comments that helped to improve the
quality and perspectives of this work. The authors thank D. Barret for
fruitful discussions that led to improve this manuscript. NC thanks A. Marin-
Laflèche for help in the selection of an FFT algorithm suited to this
research, A. Proag for useful discussions on geometrical calculations and L.
Le Briquer for commenting a preliminary version of this manuscript. This
research made use of Astropy,555http://www.astropy.org a community-developed
core Python package for Astronomy (Astropy Collaboration et al., 2013; Price-
Whelan et al., 2018). This research made use of Matplotlib (Hunter, 2007).

## References

  * Anorve-Zeferino (2019) Anorve-Zeferino, G. A. 2019, MNRAS, 483, 704
  * Armstrong (1998) Armstrong, M. 1998, Basic linear geostatistics (Springer Science & Business Media)
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Barret et al. (2016) Barret, D., Lam Trong, T., den Herder, J.-W., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, 99052F
  * Barret et al. (2018) Barret, D., Lam Trong, T., den Herder, J.-W., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, 106991G
  * Cavaliere & Fusco-Femiano (1978) Cavaliere, A. & Fusco-Femiano, R. 1978, A&A, 70, 677
  * Cayón (2010) Cayón, L. 2010, MNRAS, 405, 1084
  * Churazov et al. (2012) Churazov, E., Vikhlinin, A., Zhuravleva, I., et al. 2012, MNRAS, 421, 1123
  * Corstanje et al. (2008) Corstanje, R., Grunwald, S., & Lark, R. 2008, Geoderma, 143, 123
  * Cressie (1985) Cressie, N. 1985, Journal of the International Association for Mathematical Geology, 17, 563
  * Esquivel et al. (2007) Esquivel, A., Lazarian, A., Horibe, S., et al. 2007, MNRAS, 381, 1733
  * Gaspari et al. (2014) Gaspari, M., Churazov, E., Nagai, D., Lau, E. T., & Zhuravleva, I. 2014, A&A, 569, A67
  * Gaspari et al. (2018) Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167
  * Haslett (1997) Haslett, J. 1997, Journal of the Royal Statistical Society. Series D (The Statistician), 46, 475
  * 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. (2018) Ishisaki, Y., Ezoe, Y., Yamada, S., et al. 2018, Journal of Low Temperature Physics, 193, 991
  * Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Li & Laizet (2010) Li, N. & Laizet, S. 2010, in Proceedings of Cray User Group 2010 conference, Edinburgh
  * Lis et al. (1998) Lis, D. C., Keene, J., Li, Y., Phillips, T. G., & Pety, J. 1998, The Astrophysical Journal, 504, 889
  * Martínez et al. (2010) Martínez, V., Arnalte-Mur, P., & Stoyan, D. 2010, A&A, 513, A22
  * Matheron (1965) Matheron, G. 1965, Les Variables régionalisées et leur estimation: une application de la théorie des fonctions aléatoires aux sciences de la nature (Masson et Cie)
  * Matheron (1973) Matheron, G. 1973, Advances in Applied Probability, 5, 439?
  * Miniati & Beresnyak (2015) Miniati, F. & Beresnyak, A. 2015, Nature, 523, 59
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123
  * Roelens et al. (2017) Roelens, M., Eyer, L., Mowlavi, N., et al. 2017, MNRAS, 472, 3230
  * Roncarelli et al. (2018) Roncarelli, M., Gaspari, M., Ettori, S., et al. 2018, A&A, 618, A39
  * Ryu et al. (2008) Ryu, D., Kang, H., Cho, J., & Das, S. 2008, Science, 320, 909
  * Shi & Zhang (2019) Shi, X. & Zhang, C. 2019, MNRAS, 487, 1072
  * Spiegel (2003) Spiegel, M. 2003, Manual de fórmulas y tablas matemáticas: 2400 fórmulas y 60 tablas, Serie de compendios Schaum (McGraw-Hill/Interamericana)
  * Valdarnini (2019) Valdarnini, R. 2019, ApJ, 874, 42
  * Vazza et al. (2018) Vazza, F., Angelinelli, M., Jones, T. W., et al. 2018, MNRAS, 481, L120
  * Yoo & Cho (2014) Yoo, H. & Cho, J. 2014, ApJ, 780, 99
  * Zhuravleva et al. (2012) Zhuravleva, I., Churazov, E., Kravtsov, A., & Sunyaev, R. 2012, MNRAS, 422, 2712
  * Zhuravleva et al. (2014) Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2014, Nature, 515, 85
  * ZuHone et al. (2016) ZuHone, J. A., Markevitch, M., & Zhuravleva, I. 2016, ApJ, 817, 110

##  Appendix A Calculations in the 1-dimensional case

This Appendix details the calculation leading to results presented in Sect. 2,
namely the expected values for the centroid shift $C$ and line broadening
$S^{2}$ and their variances, in the case of measurements along a single line-
of-sight (Eq. 4, 5, 6 and 7). These calculations are generalised in Sect. 3
for 2-dimensional diagnostics of the velocity field.

###  A.1 Statistics of the centroid

The finding $\langle C\rangle=0$ is a direct consequence of random
uncorrelated phases.

We calculate $\mathrm{Var}(C)=\langle C^{2}\rangle$ by noting that:

$\langle C^{2}\rangle=F^{-2}\iint\epsilon(x)\epsilon(x^{\prime})\langle
v(x)v(x^{\prime})\rangle{\rm d}x{\rm d}x^{\prime}$ (15)

From eq. 3, the term within brackets reads:

$\langle
v(x)v(x^{\prime})\rangle=\sum_{k}P(k)\exp{\left(ik\omega(x^{\prime}-x)\right)}$

It is then easily shown that:

$\langle
C^{2}\rangle=F^{-2}\sum_{k}P(k)\left|\int\epsilon(x)\exp\left(ik\omega
x\right){\rm d}x\right|^{2}$

The term within modulus is $\widetilde{\epsilon}(k)$, namely the $k^{th}$
Fourier coefficient of $\epsilon$. This identification leads to the expression
shown in Eq. 5.

###  A.2 Statistics of the dispersion

The variations of $S^{2}$ are due to the second term of equation 2, therefore
we will focus now on studying the statistics of the double integral $\iint G$.

####  A.2.1 Average of $\iint G$

First we write:

$A\equiv\langle\iint G\rangle=\iint\langle
G\rangle=\iint\epsilon(x)\epsilon(x^{\prime})\langle\left[v(x)-v(x^{\prime})\right]^{2}\rangle$

And:

$v(x)-v(x^{\prime})=\sum_{k}V_{k}\left(e^{ik\omega x}-e^{ik\omega
x^{\prime}}\right)$

Therefore, using Eq. 3:

$\left\langle\left[v(x)-v(x^{\prime})\right]^{2}\right\rangle=\sum_{k}P(k)\left|e^{ik\omega
x}-e^{ik\omega
x^{\prime}}\right|^{2}=2\sum_{k}P(k)\left[1-\cos(k\omega(x^{\prime}-x))\right]$

This leads to:

$A=2\sum_{k}P(k)\iint{\rm d}x{\rm
d}x^{\prime}\epsilon(x)\epsilon(x^{\prime})\times\left[1-\cos(k\omega
x^{\prime})\cos(k\omega x)-\sin(k\omega x^{\prime})\sin(k\omega x)\right]$

This expression again can be rewritten using the power-spectrum of the
emissivity:

$A=2\sum_{k}P(k)\left[F^{2}-P_{\epsilon}(k)\right]$ (16)

The average of the measured line width then reads:

$\langle S^{2}\rangle=\frac{1}{F}\int\epsilon(x)\sigma_{th}^{2}(x){\rm
d}x+\sum_{k}P(k)\left[1-\frac{P_{\epsilon}(k)}{F^{2}}\right]$

Which we write under the simple form:

$\langle
S^{2}\rangle=\overline{\sigma_{th}^{2}}+\sigma_{turb}^{2}-F^{-2}\sum_{k}P(k)P_{\epsilon}(k)$

where an horizontal bar denotes averaging of the thermal component along the
line-of-sight.

####  A.2.2 Variance of $\iint G$

We define $B$ such that $\mathrm{Var}(S^{2})=F^{-4}(B^{2}-A^{2})/4$. We note
that:

$\left\langle\left(\iint
G\right)^{2}\right\rangle=B^{2}=\int\epsilon(x)\epsilon(y)\epsilon(z)\epsilon(t)\left\langle\left[v(x)-v(y)\right]^{2}\left[v(z)-v(t)\right]^{2}\right\rangle{\rm
d}x{\rm d}y{\rm d}z{\rm d}t$ (17)

The term within brackets reads:

$\langle[.]^{2}\times[.]^{2}\rangle=\sum_{j,k,l,m}\langle
V_{j}V_{k}V_{l}V_{m}\rangle\times\left(e^{ik\omega x}-e^{ik\omega
y}\right)\left(e^{ij\omega x}-e^{ij\omega y}\right)\times\left(e^{il\omega
z}-e^{il\omega t}\right)\left(e^{im\omega z}-e^{im\omega t}\right)$ (18)

Since phases are two-by-two independent, we assume the simplest possible
expression for the 4-term product of Fourier coefficients $V_{k}$, namely:

$\langle
V_{j}V_{k}V_{l}V_{m}\rangle=\left\\{\begin{array}[]{cl}P(k)P(l)&\mathrm{if\
}(k=-j);(l=-m);(k\neq\pm l)\,\\{A\\}\\\ P(k)P(j)&\mathrm{if\
}(k=-l);(j=-m);(k\neq\pm j)\,\\{B\\}\\\ P(k)P(j)&\mathrm{if\
}(k=-m);(j=-l);(k\neq\pm j)\,\\{C\\}\\\ \langle|V_{k}|^{4}\rangle&\mathrm{if\
}(k=-j=l=-m)\,\\{D\\}\\\ \langle|V_{k}|^{4}\rangle&\mathrm{if\
}(k=-j=-l=m)\,\\{E\\}\\\ \langle|V_{k}|^{4}\rangle&\mathrm{if\
}(k=j=-l=-m)\,\\{F\\}\\\ 0&\mathrm{else}\\\ \end{array}\right.$

All cases are mutually exclusive. Remarking that conditions $B$ and $C$ lead
to identical expressions under transformation $l\leftrightarrow m$, and
similarly for $D$ and $E$, we can rewrite the sum, hence the triple integral,
with a sum of 4 terms:

$B^{2}=b_{A}+2b_{B}+2b_{D}+b_{F}$

##### First term:

If $(k=-j)$ and $(l=-m)$ the bracket writes:

$\langle.\rangle_{A}=\sum_{k\neq\pm l}P(k)P(l)\left|e^{ik\omega x}-e^{ik\omega
y}\right|^{2}\left|e^{il\omega z}-e^{il\omega t}\right|^{2}$

The integration over $x,y,z$ and $t$ provides after some algebra:

$b_{A}=4\sum_{j\neq\pm
k}P(k)P(j)\left[F^{4}-F^{2}P_{\epsilon}(k)-F^{2}P_{\epsilon}(j)+P_{\epsilon}(k)P_{\epsilon}(j)\right]=4\sum_{j\neq\pm
k}P(k)P(j)\left[F^{2}-P_{\epsilon}(k)\right]\left[F^{2}-P_{\epsilon}(j)\right]$
(19)

where we have used the same trigonometric decomposition as for deriving Eq.
16.

##### Second term:

The symetries in the expression lead to:

$b_{B}=\sum_{k\neq\pm
j}P(k)P(j)\left|\iint\epsilon(x)\epsilon(y)f_{k}(x,y)f_{j}(x,y){\rm d}x{\rm
d}y\right|^{2}$

introducing the complex function:

$f_{k}(x,y)=f_{-k}^{*}(x,y)=e^{ik\omega x}-e^{ik\omega y}$

Developing the product $f_{k}f_{j}$ we can rewrite the term under the modulus
as:

$\left|\iint...\right|^{2}=4\left|F\widetilde{\epsilon}(j+k)-\widetilde{\epsilon}(j)\widetilde{\epsilon}(k)\right|^{2}$

##### Third term:

The term within brackets is rewritten as:

$\langle.\rangle_{D}=\sum_{k}\langle|V_{k}|^{4}\rangle\left|e^{ik\omega
x}-e^{ik\omega y}\right|^{2}\left|e^{ik\omega z}-e^{ik\omega t}\right|^{2}$

which, using similar calculations as for $b_{A}$, leads to:

$b_{D}=4\sum_{k}\langle|V_{k}|^{4}\rangle\left[F^{2}-P_{\epsilon}(k)\right]^{2}$

##### Fourth term:

Similarly as for the computation of $b_{F}$ above, we find:

$b_{F}=4\sum_{k}\langle|V_{k}|^{4}\rangle\left|F\widetilde{\epsilon}(2k)-\widetilde{\epsilon}(k)^{2}\right|^{2}$

Therefore, combining previous expressions we obtain:

$\mathrm{Var}(S^{2})=2\sum_{j\neq\pm
k}P(k)P(j)\left|\frac{\widetilde{\epsilon}(j+k)}{F}-\frac{\widetilde{\epsilon}(j)\widetilde{\epsilon}(k)}{F^{2}}\right|^{2}+\sum_{k}\langle|V_{k}|^{4}\rangle\left|\frac{\widetilde{\epsilon}(2k)}{F}-\frac{\widetilde{\epsilon}(k)^{2}}{F^{2}}\right|^{2}\\\
+2\sum_{k}\left(\langle|V_{k}|^{4}\rangle-P(k)^{2}\right)\left[1-\frac{P_{\epsilon}(k)}{F^{2}}\right]^{2}$

which is rearranged so to provide Eq. 7.

##  Appendix B Fourier transform of the emissivity field $\epsilon$
($\beta$-model)

We provide here calculations of $c_{\epsilon.\mathcal{W}}$, the 3-d power-
spectrum of the emissivity $\epsilon(x,y,z)$ in a case of a $\beta$-model,
seen through a sky aperture $\mathcal{W}(y,z)$. This is particularly useful
for deriving the line centroid and broadening statistics. We assume that
$\epsilon\propto n_{e}^{2}$ where $n_{e}$ is the gas density and effectively
follows a $\beta$-model profile, as in a isothermal, isometallic intra-cluster
medium.

###  B.1 Spherical model

For a spherical $\beta$-model density with core-radius $r_{c}$ centred on
$\@vec{\theta}=0$, the emissivity is expressed as:

$\epsilon(x,\@vec{\theta})=\epsilon(0)\left(1+\frac{x^{2}+\theta^{2}}{r_{c}^{2}}\right)^{-3\beta}$

The flux integrated along the line-of-sight writes:

$\displaystyle F(\@vec{\theta})$
$\displaystyle=\epsilon(0)r_{c}u_{\beta}\left(1+\frac{\theta^{2}}{r_{c}^{2}}\right)^{1/2-3\beta}$
(20) $\displaystyle\mathrm{with:\ }u_{\beta}$
$\displaystyle=2\int_{0}^{\pi/2}\cos^{6\beta-2}(t){\rm
d}t=\sqrt{\pi}\frac{\Gamma(3\beta-1/2)}{\Gamma(3\beta)}.$

The value of $c_{\epsilon.\mathcal{W}}$ is defined as follows:

$c_{\epsilon.\mathcal{W}}(k_{x},\@vec{\xi})=\iint\mathcal{W}(\@vec{\theta})\epsilon(x,\@vec{\theta})e^{-i\omega(k_{x}x+\@vec{\xi}\cdot\@vec{\theta})}{\rm
d}x{\rm d}\@vec{\theta}=\epsilon(0)\int{\rm
d}\@vec{\theta}\mathcal{W}(\@vec{\theta})e^{-i\omega\@vec{\xi}\cdot\@vec{\theta}}\int{\rm
d}xe^{-i\omega
k_{x}x}\left(1+\frac{x^{2}+\theta^{2}}{r_{c}^{2}}\right)^{-3\beta}$

This Fourier transform is calculated first along the x-axis:

$\displaystyle\int_{-\infty}^{+\infty}\left(1+\frac{x^{2}+\theta^{2}}{r_{c}^{2}}\right)^{-3\beta}e^{-i\omega
k_{x}x}{\rm d}x$
$\displaystyle=2r_{c}^{6\beta}(\omega|k_{x}|)^{6\beta-1}\int_{0}^{+\infty}\frac{\cos(t){\rm
d}t}{\left[(\omega k_{x})^{2}(\theta^{2}+r_{c}^{2})+t^{2}\right]^{3\beta}}$
$\displaystyle=\frac{2^{3/2-3\beta}\sqrt{\pi}}{\Gamma(3\beta)}r_{c}^{6\beta}\left(\frac{\sqrt{\theta^{2}+r_{c}^{2}}}{\omega|k_{x}|}\right)^{1/2-3\beta}K_{3\beta-1/2}\left(\omega|k_{x}|\sqrt{\theta^{2}+r_{c}^{2}}\right)$
(21)

where $K_{n}$ is the modified Bessel function of the second kind666The case
$k_{x}=0$ is recovered using $\lim_{x\rightarrow 0}x^{n}K_{n}(x)=\pi
2^{n-1}/\left[\sin(n\pi)\Gamma(1-n)\right]$ (Spiegel 2003).. In the special
case of $\beta=2/3$ this formula is equivalent777Because
$\mathcal{F}_{3/2}(x)=x^{3/2}K_{3/2}(x)=\sqrt{\pi/2}(1+x)\exp(-x)$ to equation
23 in ZuHone et al. (2016), namely:
$\widetilde{\epsilon}/F=(1+\omega|k_{x}|c)\exp(-\omega|k_{x}|c)$, with
$c^{2}=\theta^{2}+r_{c}^{2}$.

Introducing $\mathcal{F}_{n}(x)=x^{n}K_{n}(x)$, the integration over the
plane-of-sky coordinates $\@vec{\theta}$ provides:

$c_{\epsilon.\mathcal{W}}(k_{x},\@vec{\xi})=\epsilon(0)r_{c}\frac{2^{3/2-3\beta}\sqrt{\pi}}{\Gamma(3\beta)}\int\left(1+\frac{\theta^{2}}{r_{c}^{2}}\right)^{1/2-3\beta}\mathcal{F}_{3\beta-1/2}\left(\omega|k_{x}|\sqrt{\theta^{2}+r_{c}^{2}}\right)\mathcal{W}(\@vec{\theta})e^{-i\omega\@vec{\xi}\cdot\@vec{\theta}}{\rm
d}\@vec{\theta}$ (22)

The unknown normalization factor $\epsilon(0)$ is unimportant in this paper,
since the Fourier transform always appears divided by the aperture flux
$F_{\mathcal{W}}$ defined by:

$F_{\mathcal{W}}=\int F(\@vec{\theta})\mathcal{W}(\@vec{\theta}){\rm
d}\@vec{\theta}=\epsilon(0)r_{c}u_{\beta}\int\left(1+\frac{\theta^{2}}{r_{c}^{2}}\right)^{1/2-3\beta}\mathcal{W}(\@vec{\theta}){\rm
d}\@vec{\theta}$

An usual practical case is for a circular aperture $\mathcal{W}$ of radius
$R_{ap}$ centred on $\@vec{\theta}=0$. For this particular case:

$\frac{c_{\epsilon.\mathcal{W}}}{F_{\mathcal{W}}}(k_{x},\@vec{\xi})=\frac{2^{5/2-3\beta}(3\beta-3/2)}{\Gamma(3\beta-1/2)}\left(1-\left[1+\frac{R_{ap}^{2}}{r_{c}^{2}}\right]^{3/2-3\beta}\right)^{-1}\times\mathcal{I}_{(R_{ap}/r_{c});(3\beta-1/2)}\left(\omega|k_{x}|r_{c},\omega\xi
r_{c}\right)$

which uses the special integral defined below and represented in Fig. 14 for
$n=3/2$ ($\beta=2/3$):

$\mathcal{I}_{p;n}(u,v)=\int_{0}^{p}\frac{tJ_{0}(vt)}{(1+t^{2})^{n}}\mathcal{F}_{n}\left(u\sqrt{1+t^{2}}\right){\rm
d}t$

Figure 14: Numerical calculations of $\mathcal{I}_{p;n}(u,v)$ for various
values of $p$ and $n=3/2$. Logarithmically spaced contours (identical in all
panels) indicate the value of the function. This function is involved in the
calculation of the Fourier transform of a spherical $\beta$-model
($n=3\beta-1/2$) observed through a concentric circular aperture of radius $p$
times the core radius.

###  B.2 Plane-constant model

The integral 22 can be simplified if the core-radius $r_{c}$ is much larger
than the typical size of the window function $\mathcal{W}$. In such case, it
is equivalent to consider an emissivity that is independent of the line-of-
sight direction $\@vec{\theta}$, i.e. $\epsilon(x,y,z)=\epsilon(x)$. An
effective impact parameter $\theta_{\mathrm{eff}}$ is introduced so that:

$\epsilon(x)=\epsilon(0)\left(1+\frac{x^{2}+\theta_{\mathrm{eff}}^{2}}{r_{c}^{2}}\right)^{-3\beta}$

The calculations above then become:

$\displaystyle c_{\epsilon.\mathcal{W}}(k_{x},\@vec{\xi})$
$\displaystyle=F\left(\theta_{\mathrm{eff}}\right)\frac{2^{3/2-3\beta}}{\Gamma(3\beta-1/2)}\mathcal{F}_{3\beta-1/2}\left(\omega|k_{x}|c\right)\widehat{\mathcal{W}}(\@vec{\xi})$
(23) $\displaystyle F_{\mathcal{W}}$
$\displaystyle=\mathcal{S}_{\mathcal{W}}F\left(\theta_{\mathrm{eff}}\right)$

where we introduced $c^{2}=\theta_{\mathrm{eff}}^{2}+r_{c}^{2}$ and
$\mathcal{S}_{\mathcal{W}}=\int\mathcal{W}$ is the area of the aperture on sky
and $\widehat{\mathcal{W}}$ its (2-d) Fourier transform.

An usual practical case is for a circular aperture $\mathcal{W}$ of radius
$R_{ap}$ and an emissivity $\epsilon(x)$ in form of a $\beta=2/3$-model
independent of the line of sight. For this particular case:

$\frac{c_{\epsilon.\mathcal{W}}(k_{x},\@vec{\xi})}{F_{\mathcal{W}}}=2e^{-\omega
c|k_{x}|}\left(1+\omega c|k_{x}|\right)\frac{J_{1}(\omega\xi
R_{ap})}{\omega\xi R_{ap}}$

with $J_{1}$ the Bessel function of the first kind and order 1.

##  Appendix C Two-dimensional power-spectrum for generic emissivity field

For a pencil-beam aperture
$\mathcal{W}(\@vec{\theta}^{\prime})=\delta(\@vec{\theta}-\@vec{\theta}^{\prime})$,
the expression for the centroid shift writes:

$C(\@vec{\theta})=\int\rho(x,\@vec{\theta})v(x,\@vec{\theta}){\rm d}x$

Since $\chi^{\@vec{\theta}}(x)=\rho(x,\@vec{\theta})$ we obtain:

$C(\@vec{\theta})=\sum_{\@vec{\underline{k}}}V_{\@vec{\underline{k}}}e^{i\omega\@vec{\xi}\cdot\@vec{\theta}}\widetilde{\chi^{\@vec{\theta}}}(k_{x})$
(24)

where tilde indicates one-dimensional Fourier transform along direction $x$.
Indeed,

$\widetilde{\chi^{\@vec{\theta}}}(k_{x})=\int{\rm d}xe^{i\omega
k_{x}x}\rho(x,\@vec{\theta})=\left(\frac{\omega}{2\pi}\right)^{2}\int{\rm
d}\@vec{\xi}^{\prime}e^{-i\omega\@vec{\xi}^{\prime}\cdot\@vec{\theta}}\widetilde{\rho}(k_{x},\@vec{\xi}^{\prime})$
(25)

with $\widetilde{\rho}$ the (3D) Fourier transform of $\rho$. The last
equality derives from the definition of the inverse 3-dimensional Fourier
transform.

We remark that $\rho(x,y,z)=\epsilon/F$ is defined in a domain of space
$\mathcal{A}$ (area $\mathcal{S}_{\mathcal{A}}$) where the flux $F(y,z)$ of
the source is non-zero (which in practice is a bounded region). If the source
is infinitely extended (as in the formal case of a $\beta$-model) the boundary
is imposed by the domain of analysis $\mathcal{A}$, e.g. the instrument field
of view. We therefore consider that $\rho$ is defined over the entire
3-dimensional space by filling regions outside of the bounded domain with
zeros, which ensures the existence of $\widetilde{\rho}$. The 2D Fourier
transform $\widehat{C}(\@vec{\xi})$ of $C(\@vec{\theta})$ is used to define:

$P_{2D}(\@vec{\xi})=\frac{1}{\mathcal{S}_{\mathcal{A}}}\left\langle\left|\widehat{C}(\@vec{\xi})\right|^{2}\right\rangle$

The weighting by the total area ensures that the total ’energy’ does not
diverge as $\mathcal{A}$ becomes large. We provide later the expression for
the limiting case of an infinitely extended analysis domain.

Equation 25 shows that at any given $k_{x}$, $\widetilde{\chi^{...}}(k_{x})$
is the 2-dimensional inverse Fourier transform of
$\widetilde{\rho}(k_{x},...)$ and then:

$\int_{\@vec{\theta}}\widetilde{\chi^{\@vec{\theta}}}(k_{x})e^{i\omega\@vec{\xi}\cdot\@vec{\theta}}{\rm
d}\@vec{\theta}=\widetilde{\rho}(k_{x},\@vec{\xi})$

Therefore:

$\widehat{C}(\@vec{\xi})=\sum_{\@vec{\underline{k}}=(k_{x},\@vec{\alpha})}V_{\@vec{\underline{k}}}\int_{\@vec{\theta}}e^{i\omega(\@vec{\alpha}+\@vec{\xi})\cdot\@vec{\theta}}\widetilde{\chi^{\@vec{\theta}}}(k_{x}){\rm
d}\@vec{\theta}$

which leads to:

$\left|\widehat{C}(\@vec{\xi})\right|^{2}=\sum_{\@vec{\underline{k}}_{1},\@vec{\underline{k}}_{2}}V_{\@vec{\underline{k}}_{1}}V^{*}_{\@vec{\underline{k}}_{2}}\int_{\@vec{\theta}_{1},\@vec{\theta}_{2}}e^{i\omega\@vec{\xi}\cdot\left(\@vec{\theta}_{1}-\@vec{\theta}_{2}\right)}e^{i\omega\left(\@vec{\alpha}_{1}\cdot\@vec{\theta}_{1}-\@vec{\alpha}_{2}\cdot\@vec{\theta}_{2}\right)}\widetilde{\chi^{\@vec{\theta}_{1}}}(k_{x_{1}})\widetilde{\chi^{\@vec{\theta}_{2}}}^{*}(k_{x_{2}}){\rm
d}\@vec{\theta}_{1}{\rm d}\@vec{\theta}_{2}$

Averaging over all possible realisations provides:

$\displaystyle P_{2D}(\@vec{\xi})$
$\displaystyle=\frac{1}{\mathcal{S}_{\mathcal{A}}}\sum_{\@vec{\underline{k}}}P_{3D}(k)\int_{\@vec{\theta}_{1},\@vec{\theta}_{2}}e^{i\omega\left(\@vec{\xi}+\@vec{\alpha}\right)\cdot\left(\@vec{\theta}_{1}-\@vec{\theta}_{2}\right)}\widetilde{\chi^{\@vec{\theta}_{1}}}(k_{x})\widetilde{\chi^{\@vec{\theta}_{2}}}^{*}(k_{x}){\rm
d}\@vec{\theta}_{1}{\rm d}\@vec{\theta}_{2}$
$\displaystyle=\frac{1}{\mathcal{S}_{\mathcal{A}}}\sum_{\@vec{\underline{k}}}P_{3D}(k)\left|\int_{\@vec{\theta}}\widetilde{\chi^{\@vec{\theta}}}(k_{x})e^{i\omega(\@vec{\xi}+\@vec{\alpha})\cdot\@vec{\theta}}{\rm
d}\@vec{\theta}\right|^{2}$
$\displaystyle=\frac{1}{\mathcal{S}_{\mathcal{A}}}\sum_{\@vec{\underline{k}}=(k_{x},\@vec{\alpha})}P_{3D}(k)P_{\rho}\left(k_{x},\@vec{\alpha}+\@vec{\xi}\right)$

Which is equivalent to:

$P_{2D}(\@vec{\xi})=\frac{1}{\mathcal{S}_{\mathcal{A}}}\sum_{k_{x},\@vec{\xi}^{\prime}}P_{3D}\left(\sqrt{k_{x}^{2}+|\@vec{\xi}^{\prime}|^{2}}\right)P_{\rho}\left(k_{x},\@vec{\xi}-\@vec{\xi}^{\prime}\right)$
(26)

We note that $P_{2D}$ is in general non isotropic, as the emissivity and the
shape of the analysis domain are arbitrary. However, in the particular case
where the normalized line-of-sight emissivity is independent of the line-of-
sight, i.e. following our previous notations
$\chi^{\@vec{\theta}}(x)\equiv\epsilon(x)/F$, we find that:

$P_{\rho}(k_{x},\@vec{\xi})=\frac{1}{F^{2}}P_{\epsilon}(k_{x})P_{\mathcal{A}}(\@vec{\xi})$

with:

$\int P_{\mathcal{A}}(\@vec{\xi}){\rm
d}\@vec{\xi}=\left(\frac{2\pi}{\omega}\right)^{2}\mathcal{S}_{\mathcal{A}}$

This leads to the following expression:

$P_{2D}(\@vec{\xi})=\frac{1}{\mathcal{S}_{\mathcal{A}}}\sum_{\@vec{\xi}^{\prime}}P_{\mathcal{A}}(\@vec{\xi}-\@vec{\xi}^{\prime})\sum_{k_{x}}\frac{P_{\epsilon}(k_{x})}{F^{2}}P_{3D}\left(\sqrt{k_{x}^{2}+|\@vec{\xi}^{\prime}|^{2}}\right)=\frac{1}{\mathcal{S}_{\mathcal{A}}}\left(\frac{\omega}{2\pi}\right)^{2}\left(P_{\mathcal{A}}\otimes
P_{2D}^{\infty}\right)(\@vec{\xi})$

with $\otimes$ representing the discrete convolution product. The power-
spectrum $P_{2D}^{\infty}$ is defined such as it matches $P_{2D}$ for an
extremely large domain of analysis. Indeed, $P_{\mathcal{A}}$ then becomes a
very peaked function around $\@vec{\xi}=0$ and we obtain (see also Zhuravleva
et al. 2012):

$P_{2D}^{\infty}(\@vec{\xi})\equiv\lim_{\mathcal{A}\rightarrow\infty}P_{2D}(\@vec{\xi})\simeq\left(\frac{2\pi}{\omega}\right)^{2}\sum_{k_{x}}P_{3D}\left(\sqrt{k_{x}^{2}+\xi^{2}}\right)\frac{P_{\epsilon}(k_{x})}{F^{2}}$
(27)

##  Appendix D Structure function for generic emissivity field

###  D.1 Formal derivation neglecting border effects

For convenience, we introduce the $\mathcal{W}$-normalized emissivity:
$\chi^{\mathcal{W}}(\@vec{\underline{x}})=F_{\mathcal{W}}^{-1}\epsilon(\@vec{\underline{x}})$.
By extension, we define
$\chi^{\@vec{\theta}}(x)=\epsilon(x,\@vec{\theta})/F(\@vec{\theta})=\rho(x,\@vec{\theta})$.
We recall that $\rho=0$ outside of the domain of analysis by construction,
this is equivalent to imposing the centroid shift to vanish outside of this
region. Using the decomposition of the velocity field in Fourier series and
the definition of the velocity power spectrum, one obtains:

$\left\langle\left|C_{\mathcal{W}}-C_{\mathcal{W}^{\prime}}\right|^{2}\right\rangle=\sum_{\@vec{\underline{k}}}P_{3D}(|\@vec{\underline{k}}|)\left|\int{\rm
d}\@vec{\underline{x}}\,e^{i\omega\@vec{\underline{k}}\cdot\@vec{\underline{x}}}\left(\mathcal{W}(\@vec{\theta})\chi^{\mathcal{W}}(\@vec{\underline{x}})-\mathcal{W}^{\prime}(\@vec{\theta})\chi^{\mathcal{W}^{\prime}}(\@vec{\underline{x}})\right)\right|^{2}$

which for the most common ”pencil-beam” window function,
$\mathcal{\@vec{\theta}}=\delta(\@vec{\theta}_{0}-\@vec{\theta})$, reduces to:

$\left\langle\left|C(\@vec{\theta}_{0}+\@vec{r})-C(\@vec{\theta}_{0})\right|^{2}\right\rangle=\sum_{\@vec{\underline{k}}}P_{3D}(|\@vec{\underline{k}}|)\left|C_{\@vec{\theta}_{0},\@vec{r}}(\@vec{\underline{k}})\right|^{2}$
(28)

with:

$C_{\@vec{\theta}_{0},\@vec{r}}(\@vec{\underline{k}})=\int{\rm
d}xe^{ik_{x}\omega
x}\left[\chi^{\@vec{\theta}_{0}+\@vec{r}}(x)e^{i\omega\@vec{\xi}\cdot\@vec{r}}-\chi^{\@vec{\theta}_{0}}(x)\right]$

The expected value for the structure function therefore writes:

$\mathrm{sf}(s)=\frac{1}{N_{p}(s)}\sum_{\@vec{\underline{k}}}P_{3D}(k)I_{s}(\@vec{\underline{k}})$

with:

$I_{s}(\@vec{\underline{k}})=\int_{\@vec{\theta},|\@vec{r}|=s}|C_{\@vec{\theta},\@vec{r}}(\@vec{\underline{k}})|^{2}$

Following notations in previous appendix, we can rewrite
$C_{\@vec{\theta},\@vec{r}}(\@vec{\underline{k}})$ into:

$C_{\@vec{\theta},\@vec{r}}(\@vec{\underline{k}})=e^{i\omega\@vec{\xi}\cdot\@vec{r}}\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})-\widetilde{\chi^{\@vec{\theta}}}(k_{x})$

we then obtain:

$|C_{\@vec{\theta},\@vec{r}}(\@vec{\underline{k}})|^{2}=|\widetilde{\chi^{\@vec{\theta}}}(k_{x})|^{2}+|\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})|^{2}-2\times\mathrm{Re}\left[e^{i\omega\@vec{\xi}\cdot\@vec{r}}\widetilde{\chi^{\@vec{\theta}}}^{*}(k_{x})\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})\right]$

In a first approximation, let us perform summation over all pairs, including
those fully comprised within the domain of analysis $\mathcal{A}$ (”inner”
pairs on Fig. 15) and those with only one end in $\mathcal{A}$ (”Ext” pairs).
By construction $\widetilde{\chi^{\@vec{\theta}}}=0$ for $\@vec{\theta}$
outside of $\mathcal{A}$. This approximation is equivalent to neglecting
border effects and correction terms are discussed in the following subsection.

Using the relation between $\widetilde{\chi}$ and $\widetilde{\rho}$
identified previously, we obtain:

$\int_{p}|\widetilde{\chi^{\@vec{\theta}}}(k_{x})|^{2}+|\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})|^{2}=2\pi\int_{\@vec{\theta}}|\widetilde{\chi^{\@vec{\theta}}}(k_{x})|^{2}=2\pi\left(\frac{\omega}{2\pi}\right)^{2}\int
P_{\rho}(k_{x},\@vec{\xi}){\rm d}\@vec{\xi}$

Using Eq. 25 and some algebra leads to:

$\widetilde{\chi^{\@vec{\theta}}}^{*}(k_{x})\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})=\left(\frac{\omega}{2\pi}\right)^{2}\int{\rm
d}\@vec{\xi}^{\prime}{\rm
d}\@vec{\xi}^{\prime\prime}\widetilde{\rho}^{*}(k_{x},\@vec{\xi}^{\prime})\widetilde{\rho}(k_{x},\@vec{\xi}^{\prime\prime})e^{i\omega(\@vec{\xi}^{\prime}-\@vec{\xi}^{\prime\prime})\cdot\@vec{\theta}}e^{-i\omega\@vec{\xi}^{\prime\prime}\cdot\@vec{r}}$

Summing the term under the $\mathrm{Re}$ function over all pairs (without
double-counting) we write:

$\displaystyle\int_{p}\left[e^{i\omega\@vec{\xi}\cdot\@vec{r}}\widetilde{\chi^{\@vec{\theta}}}^{*}(k_{x})\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})\right]=$
$\displaystyle\frac{1}{2}\int_{\@vec{\theta}}\int_{|\@vec{r}|=s}\left[e^{i\omega\@vec{\xi}\cdot\@vec{r}}\widetilde{\chi^{\@vec{\theta}}}^{*}(k_{x})\widetilde{\chi^{\@vec{\theta}+\@vec{r}}}(k_{x})\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{\omega}{2\pi}\right)^{2}\int{\rm
d}\@vec{\xi}^{\prime}P_{\rho}(k_{x},\@vec{\xi}^{\prime})\int_{0}^{2\pi}e^{i\omega|\@vec{\xi}+\@vec{\xi}^{\prime}|s\cos\phi}{\rm
d}\phi$ $\displaystyle=$
$\displaystyle\pi\left(\frac{\omega}{2\pi}\right)^{2}\int{\rm
d}\@vec{\xi}^{\prime}P_{\rho}(k_{x},\@vec{\xi}^{\prime})J_{0}(\omega|\@vec{\xi}+\@vec{\xi}^{\prime}|s)$

We finally obtain:

$I_{s}(\@vec{\underline{k}})=2\pi\left(\frac{\omega}{2\pi}\right)^{2}\int{\rm
d}\@vec{\xi}^{\prime}P_{\rho}(k_{x},\@vec{\xi}^{\prime})\left[1-J_{0}\left(\left|\@vec{\xi}+\@vec{\xi}^{\prime}\right|\omega
s\right)\right]$ (29)

Dividing by the total number of pairs
$N_{p}^{tot}(s)\simeq\frac{1}{2}N_{p}(\@vec{\theta})\times
N_{p}(r=s)=\pi\mathcal{S}_{\mathcal{A}}$ provides the general expression for
the expected value of the structure function (see Eq. 3.4.1). These
calculations assume that integration over all pairs $(\@vec{\theta},\@vec{r})$
is continuous. Appendix E describes the effect of pixelized and filtered data.

###  D.2 Finite-size effects (circular domain of analysis)

Previous calculation neglects border effects in the integration over pairs of
points. We have set $\rho=0$ outside of the domain of analysis, which implies
$C=0$. Consequently a number of extra pairs are erroneously included in this
derivation, translating into extra terms
$\langle|C(\@vec{\theta}+\@vec{r})-C(\@vec{\theta})|^{2}\rangle=\langle|C(\@vec{\theta})|^{2}\rangle$
in the numerator $I_{s}$ and the number of pairs entering the denominator
$N_{p}(s)$ needs to be corrected (see Fig. 15).

Figure 15: Sketch illustrating the counting of pairs within a circular domain
of analysis of radius $R$ represented by the large black circle. Within this
domain, the centroid shift $C(\@vec{\theta})$ takes values determined by the
stochastic turbulent field, while we set $C=0$ outside. Counting inner pairs
(materialized with green and red sticks) separated by a distance $s$ is
performed by computing the range of accessible angles
$\phi_{s}(\@vec{\theta})$ for a given position in the domain, then dividing by
two. External pairs have only one end within the domain of analysis and the
range of accessible angles is $2\pi-\phi_{s}$ at a given position.

The integrals shown previously run over all pairs separated by $s$ with at
least one extremity within the field of view. There are $N_{p}^{tot}(s)$ such
pairs and $N_{p}^{ext}(s)$ pairs with only one end within the field of view.
Naturally we denote $N_{p}^{in}=N_{p}^{tot}-N_{p}^{ext}$ the pairs fully
comprised within the domain analysis. A given point $\@vec{\theta}$ in the
analysis domain belongs to $\phi_{s}(\@vec{\theta})\in[0,2\pi]$ pairs in the
field of view. Let us first compute the exact number of pairs, assuming an
infinitely fine tessellation:

$N_{p}^{tot}(s)=\pi\mathcal{S}_{\mathcal{A}}+\frac{1}{2}N_{p}^{ext}(s)$

We write:

$\left(N_{p}^{tot}(s)-\frac{1}{2}N_{p}^{ext}(s)\right)\mathrm{sf}(s)=\int_{p}\langle.\rangle=\int_{\mathrm{in}}\langle.\rangle^{in}+\int_{\mathrm{ext}}\langle.\rangle^{\mathrm{ext}}=N_{p}^{in}\mathrm{sf}^{\rm
corr}(s)+\int_{\mathrm{ext}}\langle.\rangle^{ext}$

denoting by $\mathrm{sf}^{\rm
corr}(s)=1/N_{p}^{in}\int_{\mathrm{in}}\langle.\rangle$ the value of the
structure function corrected from finite-size effects. We obtain:

$\int_{\mathrm{ext}}\langle.\rangle^{ext}=\int_{\@vec{\theta}\in\mathcal{A}}(2\pi-\phi_{s}(\@vec{\theta}))\langle|C(\@vec{\theta})|^{2}\rangle{\rm
d}\@vec{\theta}$

As demonstrated in Sect. 2:

$\langle|C(\@vec{\theta})|^{2}\rangle=\sum_{\@vec{\underline{k}}}P_{3D}(\@vec{\underline{k}})\left|\widetilde{\chi^{\@vec{\theta}}}(k_{x})\right|^{2}$

Reassembling terms, we obtain the corrected mean structure function
$(\mathrm{sf}^{\mathrm{corr}})$:

$\mathrm{sf}^{\mathrm{corr}}(s)=\left(\frac{N_{p}^{ext}(s)}{2N_{p}^{in}(s)}+1\right)\mathrm{sf}(s)-\frac{1}{N_{p}^{in}(s)}\sum_{\@vec{\underline{k}}}P_{3D}(\@vec{\underline{k}})\int_{\@vec{\theta}\in\mathcal{A}}(2\pi-\phi_{s}(\@vec{\theta}))\left|\widetilde{\chi^{\@vec{\theta}}}(k_{x})\right|^{2}{\rm
d}\@vec{\theta}$ (30)

The correction term depends both on the number extra pairs and on their
separation $s$ relative to the size of velocity fluctuations.

We have:

$N_{p}^{in}(s)=\frac{1}{2}\int_{\@vec{\theta}\in\mathcal{A}}\phi_{s}(\@vec{\theta}){\rm
d}\@vec{\theta}$
$N_{p}^{ext}(s)=\int_{\@vec{\theta}\in\mathcal{A}}(2\pi-\phi_{s}(\@vec{\theta})){\rm
d}\@vec{\theta}$

Estimating $\phi$ is easy under the assumption of a circular analysis region
of radius $R$. We find the following expressions, graphically represented in
Fig. 16, left:

$\phi_{s}(\@vec{\theta})=\left\\{\begin{array}[]{cl}2\pi&\mathrm{if\
}\theta<R-s\mathrm{\ and\ }s<R\\\ 0&\mathrm{if\ }\theta<s-R\mathrm{\ and\
}s>R\\\ 2\mathrm{Arccos}\left(\frac{\theta^{2}+s^{2}-R^{2}}{2\theta
s}\right)&\mathrm{if\ }(R-s<\theta<R\mathrm{\ and\ }s<R)\mathrm{\ or\
}(s-R<\theta<R\mathrm{\ and\ }s>R)\\\ \end{array}\right.$

We therefore rewrite $N_{p}^{in}(s)=s^{2}F(R/s)$ and
$N_{p}^{out}(s)=s^{2}G(R/s)$ with:

$F(x)=\left\\{\begin{array}[]{cl}\pi^{2}(x-1)^{2}+2\pi\int_{x-1}^{x}\gamma(u;x)u{\rm
d}u&\mathrm{if\ }x>1\\\ 2\pi\int_{1-x}^{x}\gamma(u;x)u{\rm d}u&\mathrm{if\
}0.5<x<1\\\ \end{array}\right.$

and

$G(x)=\left\\{\begin{array}[]{cl}4\pi\int_{x-1}^{x}\left[\pi-\gamma(u;x)\right]u{\rm
d}u&\mathrm{if\ }x>1\\\
2\pi^{2}(1-x)^{2}+4\pi\int_{1-x}^{x}\left[\pi-\gamma(u;x)\right]u{\rm
d}u&\mathrm{if\ }0.5<x<1\\\ \end{array}\right.$

having introduced
$\gamma(u;x)=\mathrm{Arccos}\left(\frac{u^{2}-x^{2}+1}{2u}\right)$.

The expressions for the number of pairs as a function of the separation
distance are represented on Fig. 16. As expected, the number of extra pairs is
negligible for small pair separations. It equals the number of regular
(”inner”) pairs for $s\simeq 0.5R$ and becomes dominant past this value.

|
---|---
Figure 16: _Left:_ two-dimensional representation of the function
$\phi_{s}(\@vec{\theta})$ for various values of the separation $s$ and the
position $\theta$ in a circular field of view of radius $R$. _Right:_ scaling
of $N_{p}(s)$, the number of pairs of points separated by a distance $s$
within a circular field-of-view of radius $R$. ”Inner” concerns those pairs
integrally contained within the circular domain, while ”Ext” concerns those
pairs with only one end within the domain. ”All” is the sum of the two
numbers. The curves actually show $\delta(sN_{p}(s))/\delta s$, that is the
differential number of pairs per interval of $s$ expressed in units of the
radius $R$. Here $\ell\ll R$ is the side length of an elementary pixel, so
that the total number of pixels in the circle is $\pi R^{2}/\ell^{2}$.

Finally, we note that for the flat emissivity field
($\rho(x,\@vec{\theta})=\epsilon(x)/F$) we have
$|\widetilde{\chi^{\@vec{\theta}}}|^{2}=P_{\epsilon}/F^{2}$ and then:

$\mathrm{sf}^{\mathrm{corr}}(s)=\left(\frac{N_{p}^{ext}(s)}{2N_{p}^{in}(s)}+1\right)\mathrm{sf}(s)-\frac{N_{p}^{ext}(s)}{N_{p}^{in}(s)}\times\left(\frac{\omega}{2\pi}\right)^{2}\sum_{\@vec{\xi}}P_{2D}^{\infty}(\@vec{\xi})$
(31)

Since in this case the projected velocity field is stationary and isotropic,
it may be easier and more exact to compute the mean structure function using
$P_{2D}^{\infty}$ in Eq. 3.4.1, instead of involving $P_{2D}$ and applying
this correction formula. This property is used to check the validity of the
correction formula in Eq. 31. Figure 17 shows the result of our calculation
for a flat emissivity field with core radius $r_{c}=400$ kpc and a turbulent
power spectrum with injection scale $L_{inj}=10L_{diss}=200$ kpc. The domain
of analysis is circular of radius $R=70,120$ or 500 kpc. The result for an
unbounded domain is also shown. For small analysis domains ($R=70$ kpc) the
uncorrected formula induces discrepancies at small separations, due to the
high-pass behaviour of the mask $\mathcal{A}$. As the field-of-view increases
($R=500$ kpc), border effect become negligible and all structure functions
match the exact one. Finally, as $N_{p}^{in}$ approaches zero for separation
length of size $s=2R$, the correction formula becomes numerically unstable at
large separation lengths.

Figure 17: Impact of finite-size corrections on the mean structure functions
computed according to our model. The domain of analysis is a circle of radius
$R$. Both panels show the same data, in logarithmic or linear scales. The
emissivity model is of type _Xbeta_ with a core radius $r_{c}=400$ kpc. The
turbulent power spectrum has injection and dissipation scales $L_{inj}$ and
$L_{diss}$ respectively. The thick dashed line is barely visible and shows the
exact result obtained assuming an infinitely extended analysis domain. Points
at large separations $s$ are subject to slight numerical instabilities.

##  Appendix E Structure function from pixelized and/or filtered data

Previous derivations assume that the centroid shift can be measured along
every line of sight. In general, real datasets are convolved by an
instrumental point-spread function and a pixel design is effectively grouping
line-of-sights within a single spectral line measurement. Both processes are
formally close to each other, since pixelization along a regular grid can be
reformulated as a top-hat filtering followed by the selection of points at the
centre of each pixel (Fig. 18).

Figure 18: Effect of a larger pixelization when computing the structure
function. Pixel size ranges from $\ell=4,17,34,69$ kpc (from left to right).
All four panels represent the same projected velocity field (injection scale
at 100 kpc) in a galaxy cluster represented by a _Xbeta_ model of core-radius
21 kpc. In general, larger pixels reduce the power in the 2-dimensional
velocity fluctuations and roughly act like a smoothing convolution filter on
the high-resolution centroid map.

We define a new map $D(\@vec{\theta})$ as:

$D=\frac{\mathcal{F}_{\ell}*\left(FC\right)}{\mathcal{F}_{\ell}*F}$ (32)

where $*$ represents the usual convolution product, $F(\@vec{\theta})$ and
$C(\@vec{\theta})$ are respectively the flux and centroid maps as defined in
Sect. 3.1. The filter $\mathcal{F}_{\ell}$ may represent the instrumental
point-spread function, or the pixel window function (previously noted
$\mathcal{W}$) or a combination of the two, we assume its characteristic scale
is $\ell$ (e.g. instrument FWHM, pixel size, etc.) and it is normalized to 1
by integrating over all values of $\@vec{\theta}$. It is clear that
$D(\@vec{\theta})$ is the value of the centroid shift measured after the
filtering process, resulting from a flux-weighted average of individual
centroid shifts.

This formula reduces to $D(\@vec{\theta})\simeq C(\@vec{\theta})$ for
components of $C$ varying on scales much larger than $\ell$ (equivalently, for
very sharp filters). For components of $C$ oscillating on tiny scales (much
smaller than the filter size) we find $D\simeq 0$: as expected the
pixelization or filtering process suppresses information on small scales.

A useful derivation can be carried out in case of a smooth flux map, varying
on scales much larger than the filter size $\ell$. Then $F$ can be considered
constant in the convolution products and one obtains:
$D\simeq\mathcal{F}_{\ell}*C$. All previous calculations now must incorporate
this convolution product. For instance the following replacement takes place:

$C_{\@vec{\theta}_{0},\@vec{r}}\rightarrow\int{\rm
d}\@vec{\mu}\mathcal{F}_{\ell}(\@vec{\mu})e^{i\omega\@vec{\mu}\cdot\@vec{\xi}}C_{\@vec{\theta}_{0}-\@vec{\mu},\@vec{r}}$

The calculation steps are similar to previous case, thanks to permutations of
the integrals over $\@vec{\mu}$ and other integrals. It leads to the expected
value of the structure function at each $s$:

$\displaystyle\mathrm{sf}(s)$
$\displaystyle=2\sum_{\@vec{\underline{k}}}P_{3D}(k)\int{\rm
d}\@vec{\xi}^{\prime}P_{\rho}(k_{x},\@vec{\xi}^{\prime})P_{\ell}(\@vec{\xi}+\@vec{\xi}^{\prime})\left[1-J_{0}\left(\left|\@vec{\xi}+\@vec{\xi}^{\prime}\right|\omega
s\right)\right]$
$\displaystyle=2\int\left[1-J_{0}\left(\omega\left|\@vec{\xi}\right|s\right)\right]P_{D}(\@vec{\xi}){\rm
d}\@vec{\xi}$

where $P_{D}=P_{\ell}P_{2D}$ is the power-spectrum of the map
$D(\@vec{\theta})$ and $P_{\ell}$ is the power-spectrum of the filter.
Therefore, in the case of small filter sizes (relative to the flux variation
scale) it is legitimate to replace in Eq. 3.4.1 the power-spectrum of the
centroid map, $P_{2D}$, by the power-spectrum of the filtered centroid map,
$P_{D}$. In case of larger pixels, this is generally no longer valid. This is
critical in presence of a finite domain of analysis of size comparable to the
pixel size, since then border effects must be treated more carefully. Figure
19 shows the result of applying the simple prescription $P_{2D}\rightarrow
P_{\ell}P_{2D}$ to Eq. 31, with a similar parametric setup as in Fig. 17.
Since in this case the velocity field is stationary, we also have an exact
computation of the structure function obtained by neglecting the finite-size
domain, i.e. by using $P_{\ell}P_{2D}^{\infty}$ in Eq. 3.4.1. It is then
obvious that, strictly speaking, the correction formula in Eq. 31 is valid
only for unbinned data.

Figure 19: Analytical implementation of binning with pixels of size $\ell$ in
the computation of the average structure function. The turbulent velocity
field and the cluster emissivity of type _Xbeta_ are both identical to Fig.
17. Plain lines show exact results using the stationary, unbounded
($R=\infty$) velocity field, effectively replacing $P_{2D}$ by
$P_{\ell}P_{2D}^{\infty}$ in Eq. 3.4.1. Dots are obtained by combining the
correction formula Eq. 31 for a circular domain of radius $R=250$ kpc, with
the prescription $P_{2D}\rightarrow P_{\ell}P_{2D}$. Since the latter is only
valid for slowly-varying flux maps, it fails at reproducing the true structure
function if pixels have sizeable length with respect to $R$.

Finally, in addition to this ’smoothing’ effect, pixelization induces a
discretization effect, or ’aliasing’. We do not develop a calculation for this
effect. For a given 2-dimensional frequency $\omega\@vec{\xi}$, aliasing
arises for very specific combinations of separations $s$ and pixel sizes,
matching integer multiples of the associated spatial scales. It therefore
strongly depends on the exact definition of the pixel grid. Since the power-
spectrum is continuous in the inertial range, aliasing effects are smoothly
distributed across the range of separations $s$ between the injection and
dissipation scales. It is expected to mostly affect the structure function
$SF(s)$ at separations close to the dissipation and the injection scales,
where significant discontinuities show up in the power-spectrum.

##  Appendix F Derivation of the variance of the structure function

###  F.1 General expression

Following similar notations, for a given pair indexed by $i$, we write
$n_{i}=N_{p}(s_{i})$ and we introduce:

$U_{i}=\left|C\left(\@vec{\theta}_{i}+\@vec{r}_{i}\right)-C\left(\@vec{\theta}_{i}\right)\right|^{2}=\sum_{\@vec{\underline{k}},\@vec{\underline{k}^{\prime}}}V_{\@vec{\underline{k}}}V_{\@vec{\underline{k}^{\prime}}}e^{i\omega(\@vec{\xi}+\@vec{\xi}^{\prime})\cdot\@vec{\theta}_{i}}C_{\@vec{\theta}_{i},\@vec{r}_{i}}(\@vec{\underline{k}})C_{\@vec{\theta}_{i},\@vec{r}_{i}}(\@vec{\underline{k}^{\prime}})$

We note in particular that:

$\langle
U_{i}\rangle=\sum_{\@vec{\underline{k}}}P_{3D}(k)\left|C_{\@vec{\theta}_{i},\@vec{r}_{i}}(\@vec{\underline{k}})\right|^{2}$
(33)

Even in case of an infinitely dense grid of pixels, there is a source of
uncertainty arising from the stochastic nature of the velocity field itself.
To see this more clearly, we compute:

$\left\langle
SF(s_{i})SF(s_{j})\right\rangle=\left\langle\frac{1}{n_{i}n_{j}}\int_{\@vec{\theta}_{i},\@vec{r}_{i},\@vec{\theta}_{j},\@vec{r}_{j}}U_{i}U_{j}\right\rangle=\frac{1}{n_{i}n_{j}}\int_{\@vec{\theta}_{i},\@vec{r}_{i},\@vec{\theta}_{j},\@vec{r}_{j}}\langle
U_{i}U_{j}\rangle$

Neglecting border effects, $n_{i}=n_{j}\propto 2\pi\mathcal{S}_{\mathcal{A}}$.

$\langle
U_{i}U_{j}\rangle=\sum_{\@vec{\underline{k}},\@vec{\underline{k}^{\prime}},\@vec{\underline{l}},\@vec{\underline{l}^{\prime}}}\langle
V_{\@vec{\underline{k}}}V_{\@vec{\underline{k}^{\prime}}}V_{\@vec{\underline{l}}}V_{\@vec{\underline{l}^{\prime}}}\rangle
e^{i\omega(\@vec{\xi}+\@vec{\xi}^{\prime})\cdot\@vec{\theta}_{i}}e^{i\omega(\@vec{\chi}+\@vec{\chi}^{\prime})\cdot\@vec{\theta}_{j}}C_{\@vec{\theta}_{i},\@vec{r}_{i}}(k_{x},\@vec{\xi})C_{\@vec{\theta}_{i},\@vec{r}_{i}}(k_{x}^{\prime},\@vec{\xi}^{\prime})C_{\@vec{\theta}_{j},\@vec{r}_{j}}(l_{x},\@vec{\chi})C_{\@vec{\theta}_{j},\@vec{r}_{j}}(l_{x}^{\prime},\@vec{\chi}^{\prime})$

We stress that $|\@vec{r}_{i}|=s_{i}$ and that $s_{i}$ and $s_{j}$ are not
necessarily equal. Integration runs over all pairs indexed by $i$ and $j$
within the field. Similarly to the 1-dimensional case we assume that:

$\langle
V_{\@vec{\underline{j}}}V_{\@vec{\underline{k}}}V_{\@vec{\underline{l}}}V_{\@vec{\underline{m}}}\rangle=\left\\{\begin{array}[]{cl}P_{3D}(k)P_{3D}(l)&\mathrm{if\
}(\@vec{\underline{k}}=-\@vec{\underline{j}});(\@vec{\underline{l}}=-\@vec{\underline{m}});(\@vec{\underline{k}}\neq\pm\@vec{\underline{l}})\,\\{A\\}\\\
P_{3D}(k)P_{3D}(j)&\mathrm{if\
}(\@vec{\underline{k}}=-\@vec{\underline{l}});(\@vec{\underline{j}}=-\@vec{\underline{m}});(\@vec{\underline{k}}\neq\pm\@vec{\underline{j}})\,\\{B\\}\\\
P_{3D}(k)P_{3D}(j)&\mathrm{if\
}(\@vec{\underline{k}}=-\@vec{\underline{m}});(\@vec{\underline{j}}=-\@vec{\underline{l}});(\@vec{\underline{k}}\neq\pm\@vec{\underline{j}})\,\\{C\\}\\\
\langle|V_{\@vec{\underline{k}}}|^{4}\rangle&\mathrm{if\
}(\@vec{\underline{k}}=-\@vec{\underline{j}}=\@vec{\underline{l}}=-\@vec{\underline{m}})\\{D\\}\\\
\langle|V_{\@vec{\underline{k}}}|^{4}\rangle&\mathrm{if\
}(\@vec{\underline{k}}=-\@vec{\underline{j}}=-\@vec{\underline{l}}=\@vec{\underline{m}})\\{E\\}\\\
\langle|V_{\@vec{\underline{k}}}|^{4}\rangle&\mathrm{if\
}(\@vec{\underline{k}}=\@vec{\underline{j}}=-\@vec{\underline{l}}=-\@vec{\underline{m}})\\{F\\}\\\
0&\mathrm{else}\\\ \end{array}\right.$

The decomposition of the product in brackets therefore involves six terms
$b_{A},b_{B},b_{C},b_{D},b_{E},b_{F}$ with $b_{B}=b_{C}$ and $b_{D}=b_{E}$.

##### Computation of $b_{A}$:

It corresponds to a case
$\@vec{\underline{k}}=-\@vec{\underline{k}^{\prime}}$,
$\@vec{\underline{l}}=-\@vec{\underline{l}^{\prime}}$ and
$\@vec{\underline{k}}\neq\pm\@vec{\underline{l}}$. It writes:

$\displaystyle b_{A}$
$\displaystyle=(n_{i}n_{j})^{-1}\sum_{\@vec{\underline{k}}\neq\pm\@vec{\underline{l}}}P_{3D}(k)P_{3D}(l)\int\left|C_{\@vec{\theta}_{i},\@vec{r}_{i}}(\@vec{\underline{k}})\right|^{2}\left|C_{\@vec{\theta}_{j},\@vec{r}_{j}}(\@vec{\underline{l}})\right|^{2}$
$\displaystyle=(n_{i}n_{j})^{-1}\left[\sum_{\@vec{\underline{k}}}P_{3D}(k)I_{s_{i}}(\@vec{\underline{k}})\right]\left[\sum_{\@vec{\underline{l}}}P_{3D}(l)I_{s_{j}}(\@vec{\underline{l}})\right]-2(n_{i}n_{j})^{-1}\sum_{\@vec{\underline{k}}}P_{3D}(k)^{2}I_{s_{i}}(\@vec{\underline{k}})I_{s_{j}}(\@vec{\underline{k}})$

with $I_{s}$ already defined in App. D. Therefore, the first term is simply
$\mathrm{sf}(s_{i})\mathrm{sf}(s_{j})$.

##### Computation of $b_{B}$:

It corresponds to a case $\@vec{\underline{k}}=-\@vec{\underline{l}}$,
$\@vec{\underline{k}^{\prime}}=-\@vec{\underline{l}^{\prime}}$ and
$\@vec{\underline{k}}\neq\pm\@vec{\underline{k}^{\prime}}$. It writes:

$\displaystyle b_{B}$
$\displaystyle=(n_{i}n_{j})^{-1}\sum_{\@vec{\underline{k}}\neq\pm\@vec{\underline{k}^{\prime}}}P_{3D}(k)P_{3D}(k^{\prime})\int
e^{i\omega(\@vec{\xi}+\@vec{\xi}^{\prime})\cdot(\@vec{\theta}_{i}-\@vec{\theta}_{j})}C_{\@vec{\theta}_{i},\@vec{r}_{i}}(k_{x},\@vec{\xi})C_{\@vec{\theta}_{i},\@vec{r}_{i}}(k_{x}^{\prime},\@vec{\xi}^{\prime})C_{\@vec{\theta}_{j},\@vec{r}_{j}}^{*}(k_{x},\@vec{\xi})C_{\@vec{\theta}_{j},\@vec{r}_{j}}^{*}(k_{x}^{\prime},\@vec{\xi}^{\prime})$
$\displaystyle=(n_{i}n_{j})^{-1}\left(\sum_{\@vec{\underline{k}},\@vec{\underline{l}}}P_{3D}(k)P_{3D}(l)J_{s_{i}}(\@vec{\underline{k}},\@vec{\underline{l}})J_{s_{j}}^{*}(\@vec{\underline{k}},\@vec{\underline{l}})-\sum_{\@vec{\underline{k}}}P_{3D}(k)^{2}J_{s_{i}}(\@vec{\underline{k}},\@vec{\underline{k}})J_{s_{j}}^{*}(\@vec{\underline{k}},\@vec{\underline{k}})-\sum_{\@vec{\underline{k}}}P_{3D}(k)^{2}I_{s_{i}}(\@vec{\underline{k}})I_{s_{j}}(\@vec{\underline{k}})\right)$

This involves the function:

$J_{s}(\@vec{\underline{k}},\@vec{\underline{l}})=\int_{\@vec{\theta},|\@vec{r}|=s}e^{i\omega(\@vec{\xi}+\@vec{\chi})\cdot\@vec{\theta}}C_{\@vec{\theta},\@vec{r}}(\@vec{\underline{k}})C_{\@vec{\theta},\@vec{r}}(\@vec{\underline{l}}){\rm
d}\@vec{\theta}{\rm d}\@vec{r}$

By developing the expression of $C_{\@vec{\theta},\@vec{r}}$ and using
Parseval’s theorem, we find:

$J_{s}(\@vec{\underline{k}},\@vec{\underline{l}})=4\pi\left(\frac{\omega}{2\pi}\right)^{2}\int\widetilde{\rho}(k_{x},\@vec{\xi}+\@vec{\kappa})\widetilde{\rho}(l_{x},\@vec{\chi}-\@vec{\kappa})\left(1-J_{0}(\omega\kappa
s)\right){\rm d}\@vec{\kappa}$ (34)

##### Computation of $b_{D}$:

It corresponds to
$\@vec{\underline{k}}=-\@vec{\underline{k}^{\prime}}=\@vec{\underline{l}}=-\@vec{\underline{l}^{\prime}}$.
It writes:

$b_{D}=(n_{i}n_{j})^{-1}\sum_{\@vec{\underline{k}}}\left\langle\left|V_{\@vec{\underline{k}}}\right|^{4}\right\rangle
I_{s_{i}}(\@vec{\underline{k}})I_{s_{j}}(\@vec{\underline{k}})$

##### Computation of $b_{F}$:

It corresponds to
$\@vec{\underline{k}}=\@vec{\underline{k}^{\prime}}=-\@vec{\underline{l}}=-\@vec{\underline{l}^{\prime}}$.
It writes:

$b_{F}=(n_{i}n_{j})^{-1}\sum_{\@vec{\underline{k}}}\left\langle\left|V_{\@vec{\underline{k}}}\right|^{4}\right\rangle
J_{s_{i}}(\@vec{\underline{k}},\@vec{\underline{k}})J_{s_{j}}^{*}(\@vec{\underline{k}},\@vec{\underline{k}})$

Reassembling all terms together we obtain the covariance term defined by:

$\Sigma_{ij}=\left\langle SF(s_{i})SF(s_{j})\right\rangle-\left\langle
SF(s_{i})\right\rangle\left\langle SF(s_{j})\right\rangle$

and which writes:

$\Sigma_{ij}=\frac{1}{(2\pi\mathcal{S}_{\mathcal{A}})^{2}}\left[2\sum_{\@vec{\underline{k}},\@vec{\underline{l}}}P_{3D}(k)P_{3D}(l)J_{s_{i}}(\@vec{\underline{k}},\@vec{\underline{l}})J_{s_{j}}^{*}(\@vec{\underline{k}},\@vec{\underline{l}})-\sum_{\@vec{\underline{k}}}R_{\@vec{\underline{k}}}\times\left\\{2I_{s_{i}}(\@vec{\underline{k}})I_{s_{j}}(\@vec{\underline{k}})+J_{s_{i}}(\@vec{\underline{k}},\@vec{\underline{k}})J_{s_{j}}^{*}(\@vec{\underline{k}},\@vec{\underline{k}})\right\\}\right]$
(35)

The expressions for $I$ and $J$ are given by equations 29 and 34 respectively.
Notably, the second term within brackets vanishes for Rayleigh-distributed
coefficients. In principle, finite-size corrections must also apply, similarly
as for the expected value of the structure function. We do not provide such
corrections here and keep in mind that our variance estimate neglects border
effects.

###  F.2 Case of an emissivity independent of the line-of-sight

The expression above can be further simplified if the emissivity is
independent of the line-of-sight direction. Introducing
$\widehat{\mathcal{A}}$ the Fourier transform of the analysis region and
$P_{\mathcal{A}}=|\widehat{\mathcal{A}}|^{2}$ its power-spectrum, we obtain
$\widetilde{\rho}(k_{x},\@vec{\xi})=\widetilde{\epsilon}(k_{x})\widehat{\mathcal{A}}(\@vec{\xi})/F$.
We define the following functions, whose principal interest resides in the
fact that they only depend on the definition of the analysis region and can be
precomputed numerically for any given instrumental field-of-view:

$\displaystyle\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\chi};s)$
$\displaystyle=K\int\widehat{\mathcal{A}}(\@vec{\xi}+\@vec{\kappa})\widehat{\mathcal{A}}(\@vec{\chi}-\@vec{\kappa})\left(1-J_{0}(\omega\kappa
s)\right){\rm d}\@vec{\kappa}$
$\displaystyle\mathcal{T}_{\mathcal{A}}(\@vec{\xi};s)$
$\displaystyle=\mathcal{U}_{\mathcal{A}}(\@vec{\xi},-\@vec{\xi};s)=K\int
P_{\mathcal{A}}(\@vec{\kappa})\left(1-J_{0}(\omega|\@vec{\xi}+\@vec{\kappa}|s)\right){\rm
d}\@vec{\kappa}$

The normalization constant is $K=(\omega/2\pi)^{2}/\mathcal{S}_{\mathcal{A}}$
throughout this paper.

Using these functions the variance then writes:

$\Sigma_{ij}=8\left(\frac{\omega}{2\pi}\right)^{4}\sum_{\@vec{\xi},\@vec{\chi}}P_{2D}^{\infty}(\xi)P_{2D}^{\infty}(\chi)\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\chi};s_{i})\mathcal{U}_{\mathcal{A}}^{*}(\@vec{\xi},\@vec{\chi};s_{j})\\\
-4\left(\frac{\omega}{2\pi}\right)^{4}\sum_{\@vec{\xi}}Q_{2D}^{\infty}(\xi)^{2}\left\\{2\mathcal{T}_{\mathcal{A}}(\@vec{\xi};s_{i})\mathcal{T}_{\mathcal{A}}(\@vec{\xi};s_{j})+\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\xi};s_{i})\mathcal{U}_{\mathcal{A}}^{*}(\@vec{\xi},\@vec{\xi};s_{j})\right\\}$
(36)

where we introduced $P_{2D}^{\infty}$ the 2D power-spectrum of the centroid
map over an infinitely extended domain (App. C) and

$Q^{\infty}_{2D}(\xi)^{2}=\left(\frac{2\pi}{\omega}\right)^{4}\sum_{k_{x}}\frac{P_{\epsilon}(k_{x})^{2}}{F^{4}}R_{k_{x},\@vec{\xi}}$

If the moduli are Rayleigh-distributed it is evident that $Q_{2D}=0$ and the
expression for the variance (Eq. 36) depends only on the 2D power-spectrum of
the velocity. As already noted, the terms encapsulating the field-of-view
geometry ($\mathcal{U}_{\mathcal{A}},\mathcal{T}_{\mathcal{A}}$) are factored
out from the emissivity and turbulence part ($P_{2D},Q_{2D}$).

Further simplification can be made when the analysis region is extremely wide
compared to the separations $s$ and to the largest fluctuation scale of the
velocity map (still assuming a constant emissivity on sky). The limit
$\mathcal{A}\rightarrow\infty$ then applies and $\widehat{\mathcal{A}}$
becomes a strongly peaked function around 0. The above functions rewrite:

$\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\chi};s)\simeq
K\left(1-J_{0}(\omega\xi
s)\right)\int\widehat{\mathcal{A}}(\@vec{\xi}+\@vec{\kappa})\widehat{\mathcal{A}}(\@vec{\chi}-\@vec{\kappa}){\rm
d}\@vec{\kappa}=K\left(1-J_{0}(\omega\xi
s)\right)\left\\{(\widehat{\mathcal{A}}*\widehat{\mathcal{A}})(\@vec{\xi}+\@vec{\chi})\right\\}$

The sign $*$ indicates the convolution product. Since
$\mathcal{A}^{2}=\mathcal{A}$ we obtain:

$\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\chi};s)\simeq
K\left(1-J_{0}(\omega\xi
s)\right)\left(\frac{2\pi}{\omega}\right)^{2}\widehat{\mathcal{A}}(\@vec{\xi}+\@vec{\chi})$

This automatically shows that
$\mathcal{T}_{\mathcal{A}}(\@vec{\xi};s)=1-J_{0}(\omega\xi s)$ and
$\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\xi};s)=0$. Moreover:

$\displaystyle\sum_{\@vec{\chi}}P_{2D}^{\infty}(\chi)\mathcal{U}_{\mathcal{A}}(\@vec{\xi},\@vec{\chi};s_{i})\mathcal{U}^{*}_{\mathcal{A}}(\@vec{\xi},\@vec{\chi};s_{j})$
$\displaystyle\simeq\left(\frac{2\pi}{\omega}\right)^{4}K^{2}\left(1-J_{0}(\omega\xi
s_{i})\right)\left(1-J_{0}(\omega\xi
s_{j})\right)\sum_{\@vec{\chi}}P_{\mathcal{A}}(\@vec{\chi}+\@vec{\xi})P_{2D}^{\infty}(\chi)$
$\displaystyle=\left(\frac{2\pi}{\omega}\right)^{2}\mathcal{S}_{\mathcal{A}}^{-1}P_{2D}^{\infty}(\xi)\left(1-J_{0}(\omega\xi
s_{i})\right)\left(1-J_{0}(\omega\xi s_{j})\right)$

Grouping terms together in Eq. 36 leads to Equation 13.

##  Appendix G Fourier transform of the normalized emissivity field $\rho$
(spherical $\beta$-model)

We provide useful calculations for the 3-d power-spectrum of the normalised
emissivity $\rho$ in a case of a $\beta$-model. We also discuss the limiting
case of very extended sources (equivalently, very small field-of-views).

###  G.1 General case

Calculation of the 2-dimensional power-spectrum involves the calculation of
$P_{\rho}=|\widetilde{\rho}|^{2}$ with
$\rho(x,\@vec{\theta})=\epsilon(x,\@vec{\theta})/F(\@vec{\theta})$. As noted
above (App. C), it is compulsory to fill $\rho$ with zeros outside of its
domain of definition or outside of the analysis domain. We consider here an
arbitrarily large circular analysis domain $\mathcal{A}$ (of radius $R$,
centred on the source) to perform the following calculations.

Using the expression for $F(\@vec{\theta})$ derived in previous section (Eq.
20) we obtain:

$\rho(x,\@vec{\theta})=\frac{1}{r_{c}u_{\beta}}\left(1+\frac{\theta^{2}}{r_{c}^{2}}\right)^{3\beta-1/2}{\left(1+\frac{x^{2}+\theta^{2}}{r_{c}^{2}}\right)^{-3\beta}}$

The Fourier transform $\widetilde{\rho}(k_{x},\@vec{\xi})$ is calculated in
two steps. The integration over the $x$ axis is performed first and its result
is already displayed in Eq. 21. Integration over the second axis runs for all
$\@vec{\theta}\in\mathcal{A}$ and we find:

$\widetilde{\rho}(k_{x},\@vec{\xi})=\frac{2^{5/2-3\beta}\pi
r_{c}^{2}}{\Gamma(3\beta-1/2)}\times\mathcal{H}_{(R/r_{c});(3\beta-1/2)}\left(\omega|k_{x}|r_{c},\omega\xi
r_{c}\right)$

which uses the special integral defined below and represented in Fig. 20 for
$n=3/2$ ($\beta=2/3$):

$\mathcal{H}_{p;n}(u,v)=\int_{0}^{p}tJ_{0}(vt)\mathcal{F}_{n}\left(u\sqrt{1+t^{2}}\right){\rm
d}t$

recalling that $\mathcal{F}_{n}(x)=x^{n}K_{n}(x)$,

Figure 20: Numerical calculations of $\mathcal{H}_{p;n}(u,v)$ for various
values of $p$ and $n=3/2$. Logarithmically spaced contours (identical in all
panels, dashed lines for negative values) indicate the value of the function.
This function is involved in the calculation of the Fourier transform of the
normalized emissivity of a spherical $\beta$-model ($n=3\beta-1/2$) analysed
over a concentric circular aperture of radius $p$ times the core radius. The
strong oscillatory behaviour is particularly noticeable in the last panel.

###  G.2 Limit for small analysis domains ($R\ll r_{c}$)

If $R/r_{c}$ is small, those terms in $\sqrt{1+t^{2}}\simeq 1$ are
approximatively constant under the integral and:

$P_{\rho}(k_{x},\@vec{\xi})\simeq\frac{2^{3-6\beta}}{\Gamma(3\beta-1/2)^{2}}\left[\mathcal{F}_{3\beta-1/2}\left(\omega|k_{x}|r_{c}\right)\right]^{2}\times\left[2\pi
R^{2}\frac{J_{1}(\omega\xi R)}{\omega\xi R}\right]^{2}$

The rightmost factor involving the order 1 Bessel function $J_{1}$ is the
power spectrum of a circular pupil of radius $R$. Its integral over
$\@vec{\xi}$ equals $\pi R^{2}(2\pi/\omega)^{2}$ and rapidly falls to zero for
$\omega|\@vec{\xi}|\gtrsim 3.83R^{-1}$. This result can easily be generalised
to an analysis domain of arbitrary shape, introducing its power-spectrum
$P_{\mathcal{A}}$:

$P_{\rho}(k_{x},\@vec{\xi})\simeq\frac{P_{\epsilon}(k_{x};\theta=0)}{F(0)^{2}}\times
P_{\mathcal{A}}(\@vec{\xi})$ (37)

and we naturally recover the limiting case discussed several times in this
paper where the emissivity $\epsilon(x,\@vec{\theta})=\epsilon(x)$ does not
depend on the line-of-sight $\@vec{\theta}$ over the analysis domain. Finally,
$R(\ll r_{c})$ can still be very large and therefore
$P_{\mathcal{A}}(\@vec{\xi})\rightarrow\left(\frac{2\pi}{\omega}\right)^{2}\mathcal{S}_{\mathcal{A}}\delta(\@vec{\xi})$:
the component of $P_{\rho}$ in the $\@vec{\xi}$ plane can then be seen as a
sharp low-pass filter ; this corresponds to the case discussed in previous
works (e.g. ZuHone et al. 2016).

##  Appendix H Numerical validation results, continued

We show here the equivalent of Figures 5, 6, 9 and 10 for the two other sets
of 100 simulations (Table 1). All parameters remain the same apart from the
injection scale, respetively $L_{inj}=200$ kpc (Figs. 21-24) and $L_{inj}=300$
kpc (Figs. 25-28). These comparisons still demonstrate a good match between
the simulations and analytic results.

Figure 21: Figure similar to Fig. 5 for the simulation with injection scale
200 kpc. Figure 22: Figure similar to Fig. 6 for the simulation with injection
scale 200 kpc. Figure 23: Figure similar to Fig. 9 for the simulation with
injection scale 200 kpc. Figure 24: Figure similar to Fig. 10 for the
simulation with injection scale 200 kpc. Figure 25: Figure similar to Fig. 5
for the simulation with injection scale 300 kpc. Figure 26: Figure similar to
Fig. 6 for the simulation with injection scale 300 kpc. Figure 27: Figure
similar to Fig. 9 for the simulation with injection scale 300 kpc. Figure 28:
Figure similar to Fig. 10 for the simulation with injection scale 300 kpc.