© N. Clerc et al. 2019

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, 2018; Miniati & Beresnyak 2015; 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 their formula 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 the presence of large-scale turbulence. ZuHone et al. (2016) were able to 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 the 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 that have an impact upon 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 amongst all points separated by a distance s on sky. Instrumental characteristics and signal-to-noise considerations, related for example to the exposure time or the energy resolution, are deliberately excluded and addressed in a separate work (Cucchetti et al. 2019, hereafter Paper II) The results of the present work are therefore instrument-independent to some extent.

Amongst 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 diagnostics. Effects such as heterogenous sampling, non-stationarity, and anisotropies 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 of the 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 (Sect. 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 Sect. 4. We finally discuss our results in Sect. 5 and highlight two specific cases matching the future X-ray instruments XRISM/Resolve (X-ray Imaging and Spectroscopy Mission, Ishisaki et al. 2018) and Athena/X-IFU (Advanced Telescope for High-ENergy Astrophysics/X-ray Integral Field Unit, 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 θ = (y, z). Units of these coordinates are physical (kiloparsec) 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 of km s⁻¹) 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 ⟨.⟩ 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 ∑_k). The emissivity and geometrical shapes will be treated with their continuous Fourier transforms (involving the continuous summation sign ∫dk). This distinction will often be purely formal: we made this choice to clarify the calculations. One- and three-dimensional Fourier transforms are indicated with a tilde (e.g. $\stackrel{\sim }{\rho }$), two-dimensional transforms with a hat (e.g. $\hat{\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 ∼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 or 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 a 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 θ₀, such that

$$\begin{array}{cc}& F({\mathit{\theta }}_0)={\displaystyle \int }{I}_l({\mathit{\theta }}_0;E)\mathrm{d}E\hfill \\ \hfill & \delta E({\mathit{\theta }}_0)={\displaystyle (F^{-1}\int E{I}_l({\mathit{\theta }}_0,E)\mathrm{d}E)}-E_0\hfill \\ \hfill & {\mathrm{\Sigma }}^2({\mathit{\theta }}_0)=F^{-1}{\displaystyle \int }{(E-\delta E-E_0)}^2{I}_l({\mathit{\theta }}_0;E)\mathrm{d}E,\hfill \end{array}$$

where δE and Σ are the observed centroid shift (relative to reference energy E₀) and width (or broadening) of the line, and F is a normalisation factor, namely the flux in the line. We introduce the gas velocity field along the line of sight fixed by θ₀ with v(x) ≡ v(x, θ₀). We define C = cδE/E₀ and $S^2=c^2{\mathrm{\Sigma }}^2/E_0^2$ as the observed centroid shift and width in velocity space.

2.1. Emission along the line of sight

At a 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 a purely collisional origin of the emission line, the amplitude of the line is assumed to scale with emissivity $ϵ(\mathit{r})\propto n_{\text{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)³ to contain a significantly large number of line emitters: it is practically always fulfilled, even for the tenuous intra-cluster medium with a typical density of 10⁻³ cm⁻³ and kiloparsec-scale injection scales. Accounting for the Doppler shift in energy v(x)/c and emissivity ϵ(x), we therefore model the total emission at each point by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{I}_l(E)}{\mathrm{d}x}=\frac{ϵ(x)c}{E_0{\sigma }_{\mathrm{th}}(x)\sqrt{2\pi }}exp(-\frac{1}{2}{\left[\frac{E-E_0(1+v(x)/c)}{E_0{\sigma }_{\mathrm{th}}(x)/c}\right]}^2)·\end{array}$$(1)

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

$$\begin{array}{cc}& F({\mathit{\theta }}_0)={\displaystyle \int }ϵ(x)\mathrm{d}x\hfill \\ \hfill & C({\mathit{\theta }}_0)=F^{-1}{\displaystyle \int }ϵ(x)v(x)\mathrm{d}x\hfill \\ \hfill & S^2({\mathit{\theta }}_0)=F^{-1}{\displaystyle \int }ϵ(x){\sigma }_{\mathrm{th}}^2(x)\mathrm{d}x+\frac{F^{-2}}{2}{\displaystyle \int \int }G(x,x^{\prime })\mathrm{d}x\mathrm{d}{x}^{\prime }\hfill \\ \hfill & G(x,x^{\prime })=ϵ(x)ϵ(x^{\prime }){[v(x)-v(x^{\prime })]}^2.\hfill \end{array}$$(2)

These integrals extend along the line of sight indexed by x, which 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:

$$\begin{array}{c}\hfill v(x)={\displaystyle \sum _k{V}_{k}exp\left(\frac{i2\pi kx}{L}\right)·}\end{array}$$

The coefficients V_k are complex random variables, defined as $V_k=V_{-k}^{*}=\left|V_k\right|e^{i{\psi }_k}$. Here ψ_k is a random phase and |V_k| a random modulus, supposedly independent from each other. In the following, we define for convenience ω = 2π/L. We note that ⟨v⟩=0, leading to V₀ = 0. Averaging over multiple random realisations provides

$$\begin{array}{c}\hfill ⟨V_j{V}_k⟩={\delta }_{j,-k}P(k),\end{array}$$(3)

where P(k)=P(−k)=⟨|V_k|²⟩ is the power spectrum of the turbulent velocity field and δ_ij = 1 if i = j, 0 if i ≠ j.

Our hypothesis of uniform turbulence implies that the normalisation of the power spectrum matches the square of the turbulent velocity dispersion σ_turb at any given point x. It is defined through the calculation of the second moment ${\sigma }_{\mathrm{turb}}=\sqrt{⟨v^2⟩}$, where averaging occurs over random phases and moduli. Therefore (taking e.g. x = 0), ${\sigma }_{\text{turb}}^2={\displaystyle {\sum }_{k}P}(k)$. This definition does not involve the profile of the emissivity, in contrast for example to 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),

$$\begin{array}{c}\hfill |V_k|\sim \mathcal{P}(\nu )\mathrm{d}\nu =\frac{2\nu }{P(k)}{e}^{-{\nu }^2/P(k)}\mathrm{d}\nu ,\end{array}$$

which reflects the fact that the turbulent velocity is a Gaussian random field. Introducing R_k = P(k)² − Var(|V_k|²), 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 Appendix A provide the following formulas for the centroid shift:

$$\begin{array}{c}\hfill ⟨C⟩=0\end{array}$$(4)

and its variance

$$\begin{array}{c}\hfill \mathrm{Var}(C)=⟨C^2⟩=F^{-2}{\displaystyle \sum _k{P}_{ϵ}(k)P(k).}\end{array}$$(5)

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

As for the line broadening, in Appendix A we obtain

$$\begin{array}{c}\hfill ⟨S^2⟩=\overline{{\sigma }_{\mathrm{th}}^2}+{\sigma }_{\mathrm{turb}}^2-⟨C^2⟩.\end{array}$$(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 plus 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 as

$$\begin{array}{cc}\hfill \mathrm{Var}(S^2)=& 2{\displaystyle \sum _{j,k}P(k)P(j){|\frac{\stackrel{\sim }{ϵ}(j+k)}{F}-\frac{\stackrel{\sim }{ϵ}(j)\stackrel{\sim }{ϵ}(k)}{F^2}|}^2}\hfill \\ \hfill & -{\displaystyle \sum _k{R}_k\times \{{|\frac{\stackrel{\sim }{ϵ}(2k)}{F}-\frac{\stackrel{\sim }{ϵ}{(k)}^2}{F^2}|}^2+2{[1-\frac{P_{ϵ}(k)}{F^2}]}^2\}·}\hfill \end{array}$$(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)∝1/k^2α + 1 in the inertial range [k_min, k_max] and zero outside of it. The slope is α = 1/3. We simulated the line profile resulting from the projection through a (isothermal, isometallic) β-model density profile (Cavaliere & Fusco-Femiano 1978) with β = 2/3, core-radius r_c, and distance θ to the cluster centre:

$$\begin{array}{c}\hfill ϵ(x)\propto {(1+\frac{x^2+{\theta }^2}{r_c^2})}^{-3\beta }\end{array}$$

where L = 2, r_c = 0.2, and θ = 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 Appendix B): $\stackrel{\sim }{ϵ}(k)/F=exp(-\omega |k|c)(1+\omega |k|c)$ with $c^2=r_c^2+{\theta }^2$. An example of such a 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² and the simple estimate ${\sigma }_{\mathrm{turb}}^2+\overline{{\sigma }_{\mathrm{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 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² and C in general. The configuration chosen for this simulation clearly illustrates that $⟨S^2⟩<\overline{{\sigma }_{\mathrm{th}}^2}+{\sigma }_{\mathrm{turb}}^2$. This is a consequence of the injection scale (L_inj ∼ L) being much larger than the typical cluster characteristic scale (here r_c = L/10), leading to non-Gaussian line shapes (Inogamov & Sunyaev 2003).

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 two-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 2D 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 three-dimensional structure of the velocity field and the emissivity field. A further difficulty one has to face is the non-stationarity of the projected velocity field: even though the 3D velocity field is homogeneous (stationary) in the medium, spatial variations of the emissivity in general break this property.

3.1. Tridimensional velocity field

Similarly to the previous section, we define the centroid shift and line width measured over a spectral line as a sum over all individual lines of sight selected within a region. We introduce the window function 𝒲( θ )=𝒲(y, z) as equal to 1 (one) for selected lines of sight, and as zero elsewhere. The measured spectral parameters of the line are written

$$\begin{array}{cc}& F(\mathcal{W})={\displaystyle \int }{I}_l^{\mathcal{W}}(E)\mathrm{d}E\hfill \\ \hfill & \delta E(\mathcal{W})=(F^{-1}{\displaystyle \int }E\times I_l^{\mathcal{W}}(E)\mathrm{d}\mathit{\theta }\mathrm{d}E)-E_0\hfill \\ \hfill & {\mathrm{\Sigma }}^2(\mathcal{W})=F^{-1}{\displaystyle \int }{(E-\delta E-E_0)}^2{I}_l^{\mathcal{W}}(E)\mathrm{d}E\hfill \end{array}$$

by defining

$$\begin{array}{c}\hfill I_l^{\mathcal{W}}(E)={\displaystyle \int }{I}_l(\mathit{\theta },E)\mathcal{W}(\mathit{\theta })\mathrm{d}\mathit{\theta }.\end{array}$$

All results from previous sections are recovered with 𝒲(θ) = δ(θ − θ₀).

Introducing $v(\underline{\mathit{x}})\equiv v(x,y,z)$ as the line-of-sight component of the velocity, and operating the substitutions in Eq. (1) $ϵ(x)\to ϵ(\underline{\mathit{x}})$, $v(x)\to v(\underline{\mathit{x}})$, we can rewrite the observed aperture flux, centroid, and velocity dispersion as

$$\begin{array}{cc}& F_{\mathcal{W}}={\displaystyle \int }\mathrm{d}\mathit{\theta }\mathcal{W}(\mathit{\theta }){\displaystyle \int }\mathrm{d}xϵ(x,\mathit{\theta })={\displaystyle \int }\mathcal{W}F\mathrm{d}\mathit{\theta }\hfill \\ \hfill & C_{\mathcal{W}}=F_{\mathcal{W}}^{-1}{\displaystyle \int }\mathrm{d}\mathit{\theta }\mathcal{W}(\mathit{\theta }){\displaystyle \int }ϵ(x,\mathit{\theta })v(x,\mathit{\theta })\mathrm{d}x\hfill \\ \hfill & =F_{\mathcal{W}}^{-1}{\displaystyle \int }\mathcal{W}FC\mathrm{d}\mathit{\theta }\hfill \\ \hfill & S_{\mathcal{W}}^2=F_{\mathcal{W}}^{-1}{\displaystyle \int }\mathrm{d}\underline{\mathit{x}}{\sigma }_{\mathrm{th}}^2(\underline{\mathit{x}})\mathcal{W}(\mathit{\theta })ϵ(\underline{\mathit{x}})+\frac{F_{\mathcal{W}}^{-2}}{2}{\displaystyle \int }\mathrm{d}\underline{\mathit{x}}\mathrm{d}{\underline{\mathit{x}}}^{\mathbf{\prime }}G(\underline{\mathit{x}},{\underline{\mathit{x}}}^{\mathbf{\prime }})\hfill \\ \hfill & =F_{\mathcal{W}}^{-1}{\displaystyle \int }\mathcal{W}F(S^2+C^2)\mathrm{d}\mathit{\theta }-C_{\mathcal{W}}^2\hfill \\ \hfill & G(\underline{\mathit{x}},{\underline{\mathit{x}}}^{\mathbf{\prime }})=\mathcal{W}(\mathit{\theta })ϵ(\underline{\mathit{x}})\mathcal{W}({\mathit{\theta }}^{\prime })ϵ({\underline{\mathit{x}}}^{\mathbf{\prime }}){[v(\underline{\mathit{x}})-v({\underline{\mathit{x}}}^{\mathbf{\prime }})]}^2.\hfill \end{array}$$

Integrations over x range over an arbitrarily large interval [ − L, L] and over all possible values of the plane-of-sky position θ. The velocity can be written in terms of its Fourier decomposition with $\underline{\mathit{k}}=(k_x,k_y,k_z)=(k_x,\mathit{\xi })$,

$$\begin{array}{c}\hfill v(\underline{\mathit{x}})={\displaystyle \sum _{\underline{\mathit{k}}}{V}_{\underline{\mathit{k}}}exp(i\omega \underline{\mathit{k}}·\underline{\mathit{x}}),}\end{array}$$

and we note $P_{3\mathrm{D}}(\underline{\mathit{k}})=⟨|V_{\underline{\mathit{k}}}{|}^2⟩=P_{3\mathrm{D}}(k)$. We have the relations

$$\begin{array}{cc}& V_{-\underline{\mathit{k}}}=V_{\underline{\mathit{k}}}^{\ast }\hfill \\ \hfill & \langle V_{\underline{\mathit{k}}}{V}_{{\underline{\mathit{k}}}^{\mathbf{\prime }}}\rangle ={\delta }_{\underline{\mathit{k}};-{\underline{\mathit{k}}}^{\mathbf{\prime }}}{P}_{3\mathrm{D}}(k).\hfill \end{array}$$

As in the one-dimensional case (Sect. 2), we introduce $R_{\underline{\mathit{k}}}=P_{3\mathrm{D}}{(\underline{\mathit{k}})}^2-\mathrm{Var}\left(|V_{\underline{\mathit{k}}}{|}^2\right)$, such that $R_{\underline{\mathit{k}}}=0$ for Rayleigh-distributed moduli. A minimal assumption on the four-term bracket $⟨V_{\underline{\mathit{j}}}{V}_{\underline{\mathit{k}}}{V}_{\underline{\mathit{l}}}{V}_{\underline{\mathit{m}}}⟩$ is necessary and our ansatz is explicitly provided in Appendix F.

3.2. Statistics of the aperture line centroid

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

$$\begin{array}{c}\hfill ⟨C_{\mathcal{W}}⟩=0.\end{array}$$(8)

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

$$\begin{array}{c}\hfill ⟨C_{\mathcal{W}}^2⟩=\frac{1}{F_{\mathcal{W}}^2}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){|c_{ϵ.\mathcal{W}}(\underline{\mathit{k}})|}^2.}\end{array}$$(9)

Here $c_{ϵ.\mathcal{W}}(\underline{\mathit{k}})$ is the Fourier coefficient of the product $ϵ(\underline{\mathit{x}})\mathcal{W}(y,z)$. This expression differs from Eq. (E7) in Zhuravleva et al. (2012) because we do not assume ϵ to be independent of the line of sight. If instead ϵ(x, y, z) = ϵ(x) in the domain of 𝒲 ≠ 0, then $c_{ϵ.\mathcal{W}}(\underline{\mathit{k}})=\stackrel{\sim }{ϵ}(k_x)\hat{\mathcal{W}}(\mathit{\xi })$; therefore we can rewrite our finding under the factorised form

$$\begin{array}{c}\hfill ⟨C_{\mathcal{W}}^2⟩=\frac{1}{F_{\mathcal{W}}^2}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k)P_{\mathcal{W}}(\mathit{\xi })P_{ϵ}(k_x),}\end{array}$$

with P_𝒲 being the power spectrum of the window function 𝒲. 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. In Appendix B we provide a detailed calculation of the function c_ϵ.𝒲 for the case of the isothermal, isometallicity β-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 to the one-dimensional case presented before and we find

$$\begin{array}{c}\hfill ⟨S_{\mathcal{W}}^2⟩=\overline{{\sigma }_{\mathrm{th}}^2}+{\sigma }_{\mathrm{turb}}^2-F_{\mathcal{W}}^{-2}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){|c_{ϵ.\mathcal{W}}(\underline{\mathit{k}})|}^2.}\end{array}$$

The bar indicates the average of the thermal broadening over the cluster volume defined by 𝒲. With these notations the relation found in the 1D case still holds:

$$\begin{array}{c}\hfill ⟨S_{\mathcal{W}}^2⟩+⟨C_{\mathcal{W}}^2⟩=\overline{{\sigma }_{\mathrm{th}}^2}+{\sigma }_{\mathrm{turb}}^2.\end{array}$$(10)

Finally the variance of the line broadening is written as

$$\begin{array}{cc}\hfill \mathrm{Var}(S_{\mathcal{W}}^2)=& 2{\displaystyle \sum _{\underline{\mathit{k}},{\underline{\mathit{k}}}^{\mathbf{\prime }}}{P}_{3\mathrm{D}}(\underline{\mathit{k}})P_{3\mathrm{D}}({\underline{\mathit{k}}}^{\mathbf{\prime }}){|\frac{c_{ϵ.\mathcal{W}}(\underline{\mathit{k}}+{\underline{\mathit{k}}}^{\mathbf{\prime }})}{F_{\mathcal{W}}}-\frac{c_{ϵ.\mathcal{W}}(\underline{\mathit{k}})c_{ϵ.\mathcal{W}}({\underline{\mathit{k}}}^{\mathbf{\prime }})}{F_{\mathcal{W}}^2}|}^2}\hfill \\ \hfill & -{\displaystyle \sum _{\underline{\mathit{k}}}{R}_{\underline{\mathit{k}}}\times \{{|\frac{c_{ϵ.\mathcal{W}}(2\underline{\mathit{k}})}{F_{\mathcal{W}}}-\frac{c_{ϵ.\mathcal{W}}{(\underline{\mathit{k}})}^2}{F_{\mathcal{W}}^2}|}^2}\hfill \\ \hfill & +2{[1-\frac{|c_{ϵ.\mathcal{W}}{(\underline{\mathit{k}})|}^2}{F_{\mathcal{W}}^2}]}^2\},\hfill \end{array}$$(11)

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

3.4. Statistics of the structure function

We define the structure function as the integral

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

This expression simply describes an average over all pairs of regions (called bins or pixels) (𝒲, 𝒲′) separated by a distance¹d(𝒲, 𝒲′) = 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, 𝒲(θ)≡δ(θ − θ₀), we recover the standard formulation (e.g. ZuHone et al. 2016)

$$\begin{array}{c}\hfill \mathrm{SF}(s)=\frac{1}{N_p(s)}{{\displaystyle \int }}_{{\mathit{\theta }}_0,|\mathit{r}|=s}{|C({\mathit{\theta }}_0+\mathit{r})-C({\mathit{\theta }}_0)|}^2\mathrm{d}{N}_p.\end{array}$$

The integration runs over an arbitrarily large, but bounded region of sky 𝒜 of total area 𝒮_𝒜. In the following we consider 𝒜(θ) as a function taking value 1 in the analysis region and 0 outside. There, N_p(s) needs to be interpreted as the integral ∫dN_p for all (θ₀, |r|=s).

Thus 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 sf(s) and ${\sigma }_{\text{sf}}^2(s)$.

3.4.1. Expected value of SF(s): sf(s)

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

$$\begin{array}{cc}\hfill \mathrm{sf}(s)& =2K{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }{P}_{\rho }(k_x,{\mathit{\xi }}^{\prime })[1-J_0\left(|\mathit{\xi }+{\mathit{\xi }}^{\prime }|\omega s\right)]}\hfill \\ \hfill & =2{\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }[1-J_0\left(\omega \left|\mathit{\xi }\right|s\right)]P_{2\mathrm{D}}(\mathit{\xi })\mathrm{d}\mathit{\xi },\hfill \end{array}$$(12)

with ρ(x, y, z) = ϵ(x, θ)/F(θ), K = ω²/(4π²𝒮_𝒜), J₀ 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 Appendix C. One must be careful that the power spectrum P_ρ involved in these expressions is that of the normalised emissivity field ρ, which is the 3D emissivity ϵ divided by the flux map F(θ). It is strongly dependent on the choice of analysis domain 𝒜.

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

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

The latter equation resembles Eq. (29) in ZuHone et al. (2016). 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 lines of sight, (ii) isotropy of the two-dimensional power spectrum P_2D, (iii) the averaging domain allows all possible orientations of the pair vector r, and (iv) the sum over modes ξ can be written as an integral.

We provide in Appendix D a generic formula to correct the above expressions for border effects. This involves the 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 Appendix E a prescription to account for pixelisation 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_ℓ. As shown in the appendix, this prescription should not be used in combination with the correction formula for border effects, especially if pixels are of a sizeable length compared to the analysis domain. We do not provide here a complete analytical formulation accounting simultaneously for border effects and pixelisation; it may be more advantageous in such a case to numerically estimate the average structure function from its primary definition involving the C_𝒲s.

Nevertheless, an exact solution for sf(s) is obtained in case of a stationary two-dimensional velocity field – for example if ϵ(x, y, z) = ϵ(x) – by replacing P_2D by $P_{2\mathrm{D}}^{\infty }$ in Eq. (12), which is the power spectrum computed in the limit of an infinitely extended analysis domain (see Appendix C for details.) Such a 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 Appendices D and E to validate our correction formulas and to stress their limitations.

3.4.2. Variance ${\sigma }_{\text{sf}}^2(s)$

A full calculation of the covariance

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}=\mathrm{Cov}(\mathrm{SF}(s_i),\mathrm{SF}(s_j))\end{array}$$

between structure functions measured at different scales is provided in Appendix F under the assumption of negligible finite-sized effects. The complete formula is given in Eq. (F.3) and involves integrals of the Fourier transform $\stackrel{\sim }{\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 into the expression of the variance. However, we can study a simpler, practical case where the emissivity is independent of the line-of-sight direction within the given analysis field of view. For instance, this is 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 is written ρ(x, θ) = ϵ(x)/F. This leads to a decoupling of 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 Appendix F.3 we obtain the following simple formula, under the supplementary hypothesis of a very large analysis region, that is for ${\mathcal{S}}_{\mathcal{A}}^{1/2}\gg (s,L_{\mathrm{inj}},\dots$):

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}& \simeq 16\pi {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }[\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{P}_{2\mathrm{D}}^{\infty }{(\xi )}^2-{\left(\frac{\omega }{2\pi }\right)}^2{Q}_{2\mathrm{D}}^{\infty }{(\xi )}^2]\hfill \\ \hfill & \times (1-J_0(\omega \xi s_i))(1-J_0(\omega \xi s_j))\xi \mathrm{d}\xi ,\hfill \end{array}$$(13)

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

4. Numerical validation

We performed a set of numerical experiments to validate the equations derived previously. This requires the 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

$$\begin{array}{c}\hfill P_{3\mathrm{D}}(\underline{\mathit{k}})=C_n{e}^{-{(k/k_{\mathrm{diss}})}^2}{k}^{\alpha }{e}^{-{(k_{\mathrm{inj}}/k)}^2},\end{array}$$(14)

with C_n a normalisation constant with units such that ${\displaystyle \int }{P}_{3\mathrm{D}}(\underline{\mathit{k}})\mathrm{d}\underline{\mathit{k}}={\sigma }_{\mathrm{turb}}^2$ and α = −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.

The computationally demanding fast Fourier transforms (FFT) were distributed across ten processors using 2DECOMP&FFT² (Li & Laizet 2010). The histograms of the (3D) velocity standard deviation in each of the three configurations are 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” on the order of 5−10 km s⁻¹ in the L_inj = 300 kpc run, which we attribute to aliasing effects. This extra numerical scatter needs to be remembered while comparing analytic results to simulations.

The emissivity of the galaxy cluster gas is taken as the square of a (isothermal, isometallic) gas density considered either as a spherical β-model (hereafter beta) or as a β-model along the line of sight and constant over the plane of the sky (hereafter Xbeta). The β parameter is held at a value of 2/3 while the core radius takes a value of {4, 21, 54, 107, 215, 429} kpc. The normalisation 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 (σ_th = 0 km s⁻¹ 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 Appendix H for completeness.

4.2. Centroid and line broadening

We first carry out the validation of Eqs. (8)–(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 the calculation of the 3D function c_ϵ.𝒲: in Appendix B we provide the analytical formulas for both considered emissivity models and for circular apertures. Calculation of this function for spherical β-models demands slightly more computing time than for the Xbeta model.

Equation (11) requires integration over six scalar variables. Taking advantage of the isotropy of the velocity power spectrum and the 2D rotational invariance of this specific configuration, this can be reduced to five integration variables only (k_x, $k_x^{\prime }$, ξ, ξ′, and one angle ϕ). This integral is evaluated by Monte-Carlo sampling distributed over 40 computing cores by means of the MCQUAD library³. The number of samplings is 2.10⁶ and 2.10⁵ for the Xbeta and beta emissivity models respectively, and we monitor and store the statistical uncertainties produced by 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 from the expected trend appear correlated: this is striking for instance in the left-most panel showing ⟨C⟩. This behaviour is likely to disappear with a higher number of realisations.

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 cutoff in the simulated β-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 β-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.

4.3. Structure function

We restrict the numerical validation to that of Eqs. (12) 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 pixelisation of the line centroid map. The same 3 × 100 simulated boxes are projected and pixelised in square regions of size ℓ × ℓ, where ℓ = 4, 9, 17, 34, 69 kpc. Since the two-dimensional velocity field is stationary, it is fine to use $P_{\ell }{P}_{2\mathrm{D}}^{\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 E.1 illustrates how pixelisation acts on a simulated centroid shift map: it 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 δs). The average of the structure functions and the standard deviations at each s provide the numerical indicators to be compared to Eqs. (12) and (13).

Analytical calculation of P_2D is performed according to Eq. (C.3), by 2D convolution⁴ of the 3D velocity power spectrum P_3D and the function P_ρ at each frequency k_x and eventual summing over those frequencies. The whole procedure is distributed over 40 processors working in parallel. Analytical expressions for P_ρ are given in Appendix G (particularly Eq. (G.1)) for the emissivity models relevant to our validation procedure. The calculation of $P_{2\mathrm{D}}^{\infty }$ is much more straightforward (see Eq. (C.4)).

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 nine configurations, which is a very encouraging result given the various assumptions involved in both cases.

A thorough assessment of the agreement between the analytical calculations and the numerical validation is summarised in Figs. 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. (12)) 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_ℓ 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.

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 𝒜 indeed cannot be considered as infinitely large any longer, making simplification of Eq. (F.3) into Eq. (13) less accurate.

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 β-model density (beta emissivity model). The analytical formula Eq. (12) recovers the mean structure function, despite the spatial non-stationarity of the projected velocity field. We corrected for border effects using Eq. (D.3), the region analysis being a circle of diameter 520 kpc. As highlighted in Appendix E, we are not able to use the simple prescription for significantly large pixel binnings. Moreover we did not carry out the full evaluation of the sample variance using Eq. (F.3) 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}}^2}$. Such approximation of the complex emissivity field by an effective emissivity extracted at a radius of θ_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.

5. Discussion

This work provides an extension of earlier studies, amongst others those of 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 remains formal and relies on strong hypotheses such as uniform, ergodic, and isotropic turbulent velocity fields throughout the intra-cluster medium, whose physics are 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 the presence of strong cool-core clusters. A study of anisotropic turbulence is out of the scope of this paper, but one could for instance undertake similar derivation steps as shown in the appendices, dropping the assumption $P_{3\mathrm{D}}(\underline{\mathit{k}})=P_{3\mathrm{D}}(k)$.

The existence of several 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 of the exact shape of P_3D. However, it remains to be checked whether the decomposition of $⟨V_{\underline{\mathit{j}}}{V}_{\underline{\mathit{k}}}{V}_{\underline{\mathit{l}}}{V}_{\underline{\mathit{m}}}⟩$ proposed in Appendix 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 the merging of an external galaxy group that shows high emissivity in its vicinity. An interesting perspective provided by the present calculations would be the coupling between P(k), the turbulent power spectrum, and ϵ, 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, and 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) 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. 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, that is, 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_ϵ.𝒲 is also the difference between the two corresponding coefficients, both easily computed following Appendix B. Interestingly, the emissivity in each annulus can be considered as constant, ϵ(x, θ) = ϵ(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 an injection scale of 100 kpc and σ_turb = 448 km s⁻¹. 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 σ_turb) in the outskirts. The line width variance shows diverse forms of behaviour but here again the general trend is a decrease towards the outskirts.

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 for example on 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 using 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′ and a field of view of 3.4′ equivalent diameter, and (ii) Athena/X-IFU with a resolution element of 5″ and a field of view of 5′ equivalent diameter (Barret et al. 2018). We consider two observing strategies: either one single pointing towards the cluster centre, or the mapping of a ∼15′×15′ area with multiple pointings. We also consider the case of a 15″ pixelisation rebinned image for Athena/X-IFU, such that the signal-to-noise ratio of each spectrum is increased (e.g. Roncarelli et al. 2018, also Paper II). At the redshift of Coma, 1′ on sky corresponds roughly to a 27 kpc physical separation. We consider a turbulent power spectrum with σ_turb = 438 km s⁻¹, α = −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) and ZuHone et al. (2016). Figure 13 shows the 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 ∼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 for further research, in view of optimising an observational strategy for a given instrumental setup. Our formalism involves Fourier transforms of window functions (denoted 𝒜 and 𝒲) 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.

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. We obtained exact expressions for the mean and variance of C and S obtained from single line-of-sight measurements through arbitrary gas emissivity. These expressions are Eqs. (4)–(7). We numerically validated the results with Monte-Carlo simulations of turbulent velocities with Gaussian or constant amplitudes.

2. We generalised these expressions for measurements in apertures of arbitrary shapes and sizes and for arbitrary three-dimensional emissivity fields, demonstrated in Eqs. (8)–(11). In Appendix B we provide useful formulas for the common β-model and for circular apertures. We numerically validated the formulas using Monte-Carlo simulations of 3D velocity fields in a range of emissivity and power-spectrum configurations.

3. We derived an expression for the mean structure function under the assumption of negligible border effects (Eq. (12)). 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 Appendix C. In Appendix G we provide useful formulas for the common β-model and for circular domains of analysis.

4. In Appendix D we provide a correction formula for border effects (Eq. (D.3)) 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. In Appendix E we provide 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. We derived a fairly generic expression for the sample variance of SF under the assumption of negligible border effects and for arbitrary emissivity fields (Eq. (F.3)). This equation takes a tractable form in case of flat emissivity and very large domain of analysis (Eq. (13)).

7. We numerically validated our results for the sample mean and variance of SF in the case of “flat” emissivity fields (β-models with a range of core radii) and various binnings.

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

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).

---

Acknowledgments

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 the improvement of 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 on a preliminary version of this manuscript. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013; Price-Whelan et al. 2018). This research made use of Matplotlib (Hunter 2007).

References

Appendix A: Calculations in the one-dimensional case

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

A.1. Statistics of the centroid

The finding ⟨C⟩=0 is a direct consequence of random uncorrelated phases. We calculate Var(C) = ⟨C²⟩ by noting that

$$\begin{array}{c}\hfill ⟨C^2⟩=F^{-2}{\displaystyle \int }{\displaystyle \int }ϵ(x)ϵ(x^{\prime })⟨v(x)v(x^{\prime })⟩\mathrm{d}x\mathrm{d}{x}^{\prime }.\end{array}$$(A.1)

From Eq. (3), the term within brackets reads

$$\begin{array}{c}\hfill ⟨v(x)v(x^{\prime })⟩={\displaystyle \sum _{k}P(k)exp(ik\omega (x^{\prime }-x)).}\end{array}$$

It is then easily shown that

$$\begin{array}{c}\hfill ⟨C^2⟩=F^{-2}{\displaystyle \sum _{k}P(k){|{\displaystyle \int }ϵ(x)exp\left(ik\omega x\right)\mathrm{d}x|}^2.}\end{array}$$

The term within modulus is $\stackrel{\sim }{ϵ}(k)$, namely the kth Fourier coefficient of ϵ. This identification leads to the expression shown in Eq. (5).

A.2. Statistics of the dispersion

The variations of S² are due to the second term of Eq. (2), therefore we will focus now on studying the statistics of the double integral ∬G.

A.2.1. Average of ∬G

First we write

$$\begin{array}{c}\hfill A\equiv \langle {\displaystyle \int }{\displaystyle \int }G\rangle ={\displaystyle \int }{\displaystyle \int }⟨G⟩={\displaystyle \int }{\displaystyle \int }ϵ(x)ϵ(x^{\prime })\langle {[v(x)-v(x^{\prime })]}^2\rangle \end{array}$$

and

$$\begin{array}{c}\hfill v(x)-v(x^{\prime })={\displaystyle \sum _k{V}_k(e^{ik\omega x}-e^{ik\omega x^{\prime }}).}\end{array}$$

Therefore, using Eq. (3)

$$\begin{array}{cc}\hfill \langle {[v(x)-v(x^{\prime })]}^2\rangle & ={\displaystyle \sum _{k}P(k){|e^{ik\omega x}-e^{ik\omega x^{\prime }}|}^2}\hfill \\ \hfill & =2{\displaystyle \sum _{k}P(k)[1-cos(k\omega (x^{\prime }-x))].}\hfill \end{array}$$

This leads to

$$\begin{array}{cc}\hfill A=& 2{\displaystyle \sum _{k}P(k){\displaystyle \int }{\displaystyle \int }\mathrm{d}x\mathrm{d}{x}^{\prime }ϵ(x)ϵ(x^{\prime })}\hfill \\ \hfill & \times [1-cos(k\omega x^{\prime })cos(k\omega x)-sin(k\omega x^{\prime })sin(k\omega x)].\hfill \end{array}$$

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

$$\begin{array}{c}\hfill A=2{\displaystyle \sum _{k}P(k)[F^2-P_{ϵ}(k)].}\end{array}$$(A.2)

The average of the measured line width then reads

$$\begin{array}{c}\hfill \langle S^2\rangle =\frac{1}{F}{\displaystyle \int }ϵ(x){\sigma }_{\mathrm{th}}^2(x)\mathrm{d}x+{\displaystyle \sum _{k}P(k)[1-\frac{P_{ϵ}(k)}{F^2}],}\end{array}$$

which we write under the simple form

$$\begin{array}{c}\hfill \langle S^2\rangle =\overline{{\sigma }_{\mathrm{th}}^2}+{\sigma }_{\mathrm{turb}}^2-F^{-2}{\displaystyle \sum _{k}P(k)P_{ϵ}(k),}\end{array}$$

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

A.2.2. Variance of ∬G

We define B such that Var(S²)=F⁻⁴(B² − A²)/4. We note that

$$\begin{array}{cc}\hfill \langle {\left({\displaystyle \int }{\displaystyle \int }G\right)}^2\rangle & =B^2\hfill \\ \hfill & ={\displaystyle \int }ϵ(x)ϵ(y)ϵ(z)ϵ(t)\hfill \\ \hfill & \langle {[v(x)-v(y)]}^2{[v(z)-v(t)]}^2\rangle \mathrm{d}x\mathrm{d}y\mathrm{d}z\mathrm{d}t.\hfill \end{array}$$(A.3)

The term within brackets reads

$$\begin{array}{cc}\hfill \langle {[.]}^2\times {[.]}^2\rangle =& {\displaystyle \sum _{j,k,l,m}⟨V_j{V}_k{V}_l{V}_m⟩\times (e^{ik\omega x}-e^{ik\omega y})(e^{ij\omega x}-e^{ij\omega y})}\hfill \\ \hfill & \times (e^{il\omega z}-e^{il\omega t})(e^{im\omega z}-e^{im\omega t}).\hfill \end{array}$$(A.4)

Since phases are two-by-two independent, we assume the simplest possible expression for the four-term product of Fourier coefficients V_k, namely:

$$\begin{array}{c}\hfill ⟨V_j{V}_k{V}_l{V}_m⟩=\{\begin{array}{cc}P(k)P(l)& \mathrm{i}f(k=-j);(l=-m);(k\ne \pm l)\{A\}\hfill \\ P(k)P(j)& \mathrm{i}f(k=-l);(j=-m);(k\ne \pm j)\{B\}\hfill \\ P(k)P(j)& \mathrm{i}f(k=-m);(j=-l);(k\ne \pm j)\{C\}\hfill \\ ⟨|V_k{|}^4⟩& \mathrm{i}f(k=-j=l=-m)\{D\}\hfill \\ ⟨|V_k{|}^4⟩& \mathrm{i}f(k=-j=-l=m)\{E\}\hfill \\ ⟨|V_k{|}^4⟩& \mathrm{i}f(k=j=-l=-m)\{F\}\hfill \\ 0& \mathrm{else}\hfill \end{array}\end{array}$$

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

$$\begin{array}{c}\hfill B^2=b_A+2b_B+2b_D+b_F.\end{array}$$

First term: If (k = −j) and (l = −m) the bracket is written

$$\begin{array}{c}\hfill {⟨.⟩}_A={\displaystyle \sum _{k\ne \pm l}P(k)P(l){|e^{ik\omega x}-e^{ik\omega y}|}^2{|e^{il\omega z}-e^{il\omega t}|}^2.}\end{array}$$

After some algebra, the integration over x, y, z, and t provides

$$\begin{array}{cc}\hfill b_A& =4{\displaystyle \sum _{j\ne \pm k}P(k)P(j)[F^4-F^2{P}_{ϵ}(k)-F^2{P}_{ϵ}(j)+P_{ϵ}(k)P_{ϵ}(j)]}\hfill \\ \hfill & =4{\displaystyle \sum _{j\ne \pm k}P(k)P(j)[F^2-P_{ϵ}(k)][F^2-P_{ϵ}(j)],}\hfill \end{array}$$(A.5)

where we have used the same trigonometric decomposition as in the derivation of Eq. (A.2).

Second term: The symmetries in the expression lead to

$$\begin{array}{c}\hfill b_B={\displaystyle \sum _{k\ne \pm j}P(k)P(j){\left|{\displaystyle \int }{\displaystyle \int }ϵ(x)ϵ(y)f_k(x,y)f_j(x,y)\mathrm{d}x\mathrm{d}y\right|}^2,}\end{array}$$

introducing the complex function

$$\begin{array}{c}\hfill f_k(x,y)=f_{-k}^{\ast }(x,y)=e^{ik\omega x}-e^{ik\omega y}.\end{array}$$

Developing the product f_kf_j we can rewrite the term under the modulus as

$$\begin{array}{c}\hfill {|{\displaystyle \int }{\displaystyle \int }\dots |}^2=4{|F\stackrel{\sim }{ϵ}(j+k)-\stackrel{\sim }{ϵ}(j)\stackrel{\sim }{ϵ}(k)|}^2.\end{array}$$

Third term: The term within brackets is rewritten as

$$\begin{array}{c}\hfill {⟨.⟩}_D={\displaystyle \sum _k\langle |V_k{|}^4\rangle {|e^{ik\omega x}-e^{ik\omega y}|}^2{|e^{ik\omega z}-e^{ik\omega t}|}^2,}\end{array}$$

which, using similar calculations as for b_A, leads to

$$\begin{array}{c}\hfill b_D=4{\displaystyle \sum _k\langle |V_k{|}^4\rangle {[F^2-P_{ϵ}(k)]}^2.}\end{array}$$

Fourth term: In a similar way to the computation of b_F above, we find

$$\begin{array}{c}\hfill b_F=4{\displaystyle \sum _k\langle |V_k{|}^4\rangle {|F\stackrel{\sim }{ϵ}(2k)-\stackrel{\sim }{ϵ}{(k)}^2|}^2.}\end{array}$$

Therefore, combining previous expressions we obtain

$$\begin{array}{cc}\hfill \mathrm{Var}(S^2)=& 2{\displaystyle \sum _{j\ne \pm k}P(k)P(j){|\frac{\stackrel{\sim }{ϵ}(j+k)}{F}-\frac{\stackrel{\sim }{ϵ}(j)\stackrel{\sim }{ϵ}(k)}{F^2}|}^2}\hfill \\ \hfill & +{\displaystyle \sum _k\langle |V_k{|}^4\rangle {|\frac{\stackrel{\sim }{ϵ}(2k)}{F}-\frac{\stackrel{\sim }{ϵ}{(k)}^2}{F^2}|}^2}\hfill \\ \hfill & +2{\displaystyle \sum _k(⟨|V_k{|}^4⟩-P{(k)}^2){[1-\frac{P_{ϵ}(k)}{F^2}]}^2,}\hfill \end{array}$$

which is rearranged to provide Eq. (7).

Appendix B: Fourier transform of the emissivity field ϵ ( β-model)

We provide here calculations of c_ϵ.𝒲, the 3D power spectrum of the emissivity ϵ(x, y, z) in the case of a β-model, seen through a sky aperture 𝒲(y, z). This is particularly useful for deriving the line centroid and broadening statistics. We assume that $ϵ\propto n_e^2$ where n_e is the gas density and effectively follows a β-model profile, as in an isothermal, isometallic intra-cluster medium.

B.1. Spherical model

For a spherical β-model density with core radius r_c centred on θ = 0, the emissivity is expressed as

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

The flux integrated along the line of sight is written

$$\begin{array}{cc}& F(\mathit{\theta })=ϵ(0)r_c{u}_{\beta }{(1+\frac{{\theta }^2}{r_c^2})}^{1/2-3\beta }\hfill \\ \hfill & \mathrm{w}ith:u_{\beta }=2{{\displaystyle \int }}_0^{\pi /2}{cos}^{6\beta -2}(t)\mathrm{d}t=\sqrt{\pi }\frac{\mathrm{\Gamma }(3\beta -1/2)}{\mathrm{\Gamma }(3\beta )}·\hfill \end{array}$$(B.1)

The value of c_ϵ.𝒲 is defined as

$$\begin{array}{cc}\hfill c_{ϵ.\mathcal{W}}(k_x,\mathit{\xi })& ={\displaystyle \int }{\displaystyle \int }\mathcal{W}(\mathit{\theta })ϵ(x,\mathit{\theta })e^{-i\omega (k_{x}x+\mathit{\xi }·\mathit{\theta })}\mathrm{d}x\mathrm{d}\mathit{\theta }\hfill \\ \hfill & =ϵ(0){\displaystyle \int }\mathrm{d}\mathit{\theta }\mathcal{W}(\mathit{\theta })e^{-i\omega \mathit{\xi }·\mathit{\theta }}{\displaystyle \int }\mathrm{d}x{e}^{-i\omega k_{x}x}{(1+\frac{x^2+{\theta }^2}{r_c^2})}^{-3\beta }·\hfill \end{array}$$

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

$$\begin{array}{cc}& {{\displaystyle \int }}_{-\infty }^{+\infty }{(1+\frac{x^2+{\theta }^2}{r_c^2})}^{-3\beta }{e}^{-i\omega k_{x}x}\mathrm{d}x\hfill \\ \hfill & =2r_c^{6\beta }(\omega |k_x{|)}^{6\beta -1}{{\displaystyle \int }}_0^{+\infty }\frac{cos(t)\mathrm{d}t}{{[{(\omega k_x)}^2({\theta }^2+r_c^2)+t^2]}^{3\beta }}\hfill \\ \hfill & =\frac{2^{3/2-3\beta }\sqrt{\pi }}{\mathrm{\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),\hfill \end{array}$$(B.2)

where K_n is the modified Bessel function of the second kind⁶. In the special case of β = 2/3, this formula is equivalent⁷ to Eq. (23) in ZuHone et al. (2016), namely $\stackrel{\sim }{ϵ}/F=(1+\omega |k_x|c)exp(-\omega |k_x|c)$, with $c^2={\theta }^2+r_c^2$.

Introducing ℱ_n(x)=xⁿK_n(x), the integration over the plane-of-sky coordinates θ provides

$$\begin{array}{cc}\hfill c_{ϵ.\mathcal{W}}(k_x,\mathit{\xi })& =ϵ(0)r_c\frac{2^{3/2-3\beta }\sqrt{\pi }}{\mathrm{\Gamma }(3\beta )}\hfill \\ \hfill & {\displaystyle \int }{(1+\frac{{\theta }^2}{r_c^2})}^{1/2-3\beta }{\mathcal{F}}_{3\beta -1/2}\left(\omega |k_x|\sqrt{{\theta }^2+r_c^2}\right)\mathcal{W}(\mathit{\theta })e^{-i\omega \mathit{\xi }·\mathit{\theta }}\mathrm{d}\mathit{\theta }.\hfill \end{array}$$(B.3)

The unknown normalisation factor ϵ(0) is unimportant in this paper, since the Fourier transform always appears divided by the aperture flux F_𝒲 defined by

$$\begin{array}{c}\hfill F_{\mathcal{W}}={\displaystyle \int }F(\mathit{\theta })\mathcal{W}(\mathit{\theta })\mathrm{d}\mathit{\theta }=ϵ(0)r_c{u}_{\beta }{\displaystyle \int }{(1+\frac{{\theta }^2}{r_c^2})}^{1/2-3\beta }\mathcal{W}(\mathit{\theta })\mathrm{d}\mathit{\theta }.\end{array}$$

An usual practical case is for a circular aperture 𝒲 of radius R_ap centred on θ = 0. For this particular case,

$$\begin{array}{cc}\hfill \frac{c_{ϵ.\mathcal{W}}}{F_{\mathcal{W}}}(k_x,\mathit{\xi })=& \frac{2^{5/2-3\beta }(3\beta -3/2)}{\mathrm{\Gamma }(3\beta -1/2)}{(1-{[1+\frac{R_{\mathrm{ap}}^2}{r_c^2}]}^{3/2-3\beta })}^{-1}\hfill \\ \hfill & \times {\mathcal{I}}_{(R_{\mathrm{ap}}/r_c);(3\beta -1/2)}(\omega |k_x|r_c,\omega \xi r_c),\hfill \end{array}$$

which uses the special integral defined below and represented in Fig. B.1 for n = 3/2 (β = 2/3):

$$\begin{array}{c}\hfill {\mathcal{I}}_{p;n}(u,v)={{\displaystyle \int }}_0^p\frac{t{J}_0(vt)}{{(1+t^2)}^n}{\mathcal{F}}_n\left(u\sqrt{1+t^2}\right)\mathrm{d}t.\end{array}$$

B.2. Plane-constant model

The integral B.3 can be simplified if the core radius r_c is much larger than the typical size of the window function 𝒲. In such a case, it is equivalent to consider an emissivity that is independent of the line-of-sight direction θ, that is, ϵ(x, y, z)=ϵ(x). An effective impact parameter θ_eff is introduced so that

$$\begin{array}{c}\hfill ϵ(x)=ϵ(0){(1+\frac{x^2+{\theta }_{\mathrm{eff}}^2}{r_c^2})}^{-3\beta }·\end{array}$$

The calculations above then become

$$\begin{array}{cc}& c_{ϵ.\mathcal{W}}(k_x,\mathit{\xi })=F\left({\theta }_{\mathrm{eff}}\right)\frac{2^{3/2-3\beta }}{\mathrm{\Gamma }(3\beta -1/2)}{\mathcal{F}}_{3\beta -1/2}\left(\omega |k_x|c\right)\hat{\mathcal{W}}(\mathit{\xi })\hfill \\ \hfill & F_{\mathcal{W}}={\mathcal{S}}_{\mathcal{W}}F\left({\theta }_{\mathrm{eff}}\right),\hfill \end{array}$$(B.4)

where we introduced $c^2={\theta }_{\text{eff}}^2+r_c^2$ and 𝒮_𝒲 = ∫𝒲 is the area of the aperture on sky and $\hat{\mathcal{W}}$ its (2D) Fourier transform.

An usual practical case is for a circular aperture 𝒲 of radius R_ap and an emissivity ϵ(x) in the form of a β = 2/3-model independent of the line of sight. For this particular case,

$$\begin{array}{c}\hfill \frac{c_{ϵ.\mathcal{W}}(k_x,\mathit{\xi })}{F_{\mathcal{W}}}=2e^{-\omega c|k_x|}\left(1+\omega c|k_x|\right)\frac{J_1(\omega \xi R_{\mathrm{ap}})}{\omega \xi R_{\mathrm{ap}}},\end{array}$$

with J₁ 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 𝒲(θ′) = δ(θ − θ′), the expression for the centroid shift is written

$$\begin{array}{c}\hfill C(\mathit{\theta })={\displaystyle \int }\rho (x,\mathit{\theta })v(x,\mathit{\theta })\mathrm{d}x.\end{array}$$

Since χ^θ(x)=ρ(x, θ), we obtain

$$\begin{array}{c}\hfill C(\mathit{\theta })={\displaystyle \sum _{\underline{\mathit{k}}}{V}_{\underline{\mathit{k}}}{e}^{i\omega \mathit{\xi }·\mathit{\theta }}\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x),}\end{array}$$(C.1)

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

$$\begin{array}{c}\hfill \stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)={\displaystyle \int }\mathrm{d}x{e}^{i\omega k_{x}x}\rho (x,\mathit{\theta })={\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }{e}^{-i\omega {\mathit{\xi }}^{\prime }·\mathit{\theta }}\stackrel{\sim }{\rho }(k_x,{\mathit{\xi }}^{\prime }),\end{array}$$(C.2)

with $\stackrel{\sim }{\rho }$ the (3D) Fourier transform of ρ. The last equality derives from the definition of the inverse 3D Fourier transform.

We note that ρ(x, y, z)=ϵ/F is defined in a domain of space 𝒜 (area 𝒮_𝒜) 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 β-model), the boundary is imposed by the domain of analysis 𝒜, for example the instrument field of view. We therefore consider that ρ is defined over the entire 3D space by filling regions outside of the bounded domain with zeros, which ensures the existence of $\stackrel{\sim }{\rho }$. The 2D Fourier transform $\hat{C}(\mathit{\xi })$ of C(θ) is used to define

$$\begin{array}{c}\hfill P_{2\mathrm{D}}(\mathit{\xi })=\frac{1}{{\mathcal{S}}_{\mathcal{A}}}\langle {\left|\hat{C}(\mathit{\xi })\right|}^2\rangle ·\end{array}$$

The weighting by the total area ensures that the total “energy” does not diverge as 𝒜 becomes large. We provide later the expression for the limiting case of an infinitely extended analysis domain.

Equation (C.2) shows that at any given k_x, $\stackrel{\sim }{{\chi }^{\dots }}(k_x)$ is the 2D inverse Fourier transform of $\stackrel{\sim }{\rho }(k_x,\dots )$ and then

$$\begin{array}{c}\hfill {{\displaystyle \int }}_{\mathit{\theta }}\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)e^{i\omega \mathit{\xi }·\mathit{\theta }}\mathrm{d}\mathit{\theta }=\stackrel{\sim }{\rho }(k_x,\mathit{\xi }).\end{array}$$

Therefore,

$$\begin{array}{c}\hfill \hat{C}(\mathit{\xi })={\displaystyle \sum _{\underline{\mathit{k}}=(k_x,\mathit{\alpha })}{V}_{\underline{\mathit{k}}}{{\displaystyle \int }}_{\mathit{\theta }}{e}^{i\omega (\mathit{\alpha }+\mathit{\xi })·\mathit{\theta }}\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)\mathrm{d}\mathit{\theta },}\end{array}$$

which leads to

$$\begin{array}{cc}\hfill {\left|\hat{C}(\mathit{\xi })\right|}^2& ={\displaystyle \sum _{{\underline{\mathit{k}}}_1,{\underline{\mathit{k}}}_2}{V}_{{\underline{\mathit{k}}}_1}{V}_{{\underline{\mathit{k}}}_2}^{\ast }{{\displaystyle \int }}_{{\mathit{\theta }}_1,{\mathit{\theta }}_2}{e}^{i\omega \mathit{\xi }·({\mathit{\theta }}_1-{\mathit{\theta }}_2)}{e}^{i\omega ({\mathit{\alpha }}_1·{\mathit{\theta }}_1-{\mathit{\alpha }}_2·{\mathit{\theta }}_2)}}\hfill \\ \hfill & \stackrel{\sim }{{\chi }^{{\mathit{\theta }}_1}}(k_{x_1}){\stackrel{\sim }{{\chi }^{{\mathit{\theta }}_2}}}^{\ast }(k_{x_2})\mathrm{d}{\mathit{\theta }}_1\mathrm{d}{\mathit{\theta }}_2.\hfill \end{array}$$

Averaging over all possible realisations provides

$$\begin{array}{cc}\hfill P_{2\mathrm{D}}(\mathit{\xi })& =\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){{\displaystyle \int }}_{{\mathit{\theta }}_1,{\mathit{\theta }}_2}{e}^{i\omega (\mathit{\xi }+\mathit{\alpha })·({\mathit{\theta }}_1-{\mathit{\theta }}_2)}\stackrel{\sim }{{\chi }^{{\mathit{\theta }}_1}}(k_x){\stackrel{\sim }{{\chi }^{{\mathit{\theta }}_2}}}^{\ast }(k_x)\mathrm{d}{\mathit{\theta }}_1\mathrm{d}{\mathit{\theta }}_2}\hfill \\ \hfill & =\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){\left|{{\displaystyle \int }}_{\mathit{\theta }}\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)e^{i\omega (\mathit{\xi }+\mathit{\alpha })·\mathit{\theta }}\mathrm{d}\mathit{\theta }\right|}^2}\hfill \\ \hfill & =\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{\displaystyle \sum _{\underline{\mathit{k}}=(k_x,\mathit{\alpha })}{P}_{3\mathrm{D}}(k)P_{\rho }(k_x,\mathit{\alpha }+\mathit{\xi }),}\hfill \end{array}$$

which is equivalent to

$$\begin{array}{c}\hfill P_{2\mathrm{D}}(\mathit{\xi })=\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{\displaystyle \sum _{k_x,{\mathit{\xi }}^{\prime }}{P}_{3\mathrm{D}}\left(\sqrt{k_x^2+{|{\mathit{\xi }}^{\prime }|}^2}\right)P_{\rho }(k_x,\mathit{\xi }-{\mathit{\xi }}^{\prime }).}\end{array}$$(C.3)

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 normalised line-of-sight emissivity is independent of the line of sight, and thus following our previous notations χ^θ(x)≡ϵ(x)/F, we find that

$$\begin{array}{c}\hfill P_{\rho }(k_x,\mathit{\xi })=\frac{1}{F^2}{P}_{ϵ}(k_x)P_{\mathcal{A}}(\mathit{\xi }),\end{array}$$

with

$$\begin{array}{c}\hfill {\displaystyle \int }{P}_{\mathcal{A}}(\mathit{\xi })\mathrm{d}\mathit{\xi }={\left(\frac{2\pi }{\omega }\right)}^2{\mathcal{S}}_{\mathcal{A}}.\end{array}$$

This leads to the expression

$$\begin{array}{cc}\hfill P_{2\mathrm{D}}(\mathit{\xi })& =\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{\displaystyle \sum _{{\mathit{\xi }}^{\prime }}{P}_{\mathcal{A}}(\mathit{\xi }-{\mathit{\xi }}^{\prime })\sum _{k_x}\frac{P_{ϵ}(k_x)}{F^2}{P}_{3\mathrm{D}}\left(\sqrt{k_x^2+{|{\mathit{\xi }}^{\prime }|}^2}\right)}\hfill \\ \hfill & =\frac{1}{{\mathcal{S}}_{\mathcal{A}}}{\left(\frac{\omega }{2\pi }\right)}^2(P_{\mathcal{A}}\otimes P_{2\mathrm{D}}^{\infty })(\mathit{\xi }),\hfill \end{array}$$

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

$$\begin{array}{c}\hfill P_{2\mathrm{D}}^{\infty }(\mathit{\xi })\equiv \underset{\mathcal{A}\to \infty }{\lim }{P}_{2\mathrm{D}}(\mathit{\xi })\simeq {\left(\frac{2\pi }{\omega }\right)}^2{\displaystyle \sum _{k_x}{P}_{3\mathrm{D}}\left(\sqrt{k_x^2+{\xi }^2}\right)\frac{P_{ϵ}(k_x)}{F^2}·}\end{array}$$(C.4)

Appendix D: Structure function for generic emissivity field

D.1. Formal derivation neglecting border effects

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

$$\begin{array}{cc}\hfill \langle {|C_{\mathcal{W}}-C_{{\mathcal{W}}^{\prime }}|}^2\rangle & ={\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(|\underline{\mathit{k}}|)}\hfill \\ \hfill & {\left|{\displaystyle \int }\mathrm{d}\underline{\mathit{x}}{e}^{i\omega \underline{\mathit{k}}·\underline{\mathit{x}}}(\mathcal{W}(\mathit{\theta }){\chi }^{\mathcal{W}}(\underline{\mathit{x}})-{\mathcal{W}}^{\prime }(\mathit{\theta }){\chi }^{{\mathcal{W}}^{\prime }}(\underline{\mathit{x}}))\right|}^2,\hfill \end{array}$$

which for the most common “pencil-beam” window function, 𝒲θ = δ(θ₀ − θ), reduces to

$$\begin{array}{c}\hfill \langle {|C({\mathit{\theta }}_0+\mathit{r})-C({\mathit{\theta }}_0)|}^2\rangle ={\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(|\underline{\mathit{k}}|){\left|C_{{\mathit{\theta }}_0,\mathit{r}}(\underline{\mathit{k}})\right|}^2,}\end{array}$$(D.1)

with

$$\begin{array}{c}\hfill C_{{\mathit{\theta }}_0,\mathit{r}}(\underline{\mathit{k}})={\displaystyle \int }\mathrm{d}x{e}^{i{k}_x\omega x}[{\chi }^{{\mathit{\theta }}_0+\mathit{r}}(x)e^{i\omega \mathit{\xi }·\mathit{r}}-{\chi }^{{\mathit{\theta }}_0}(x)].\end{array}$$

The expected value for the structure function therefore is written

$$\begin{array}{c}\hfill \mathrm{sf}(s)=\frac{1}{N_p(s)}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k)I_s(\underline{\mathit{k}}),}\end{array}$$

with

$$\begin{array}{c}\hfill I_s(\underline{\mathit{k}})={{\displaystyle \int }}_{\mathit{\theta },|\mathit{r}|=s}{|C_{\mathit{\theta },\mathit{r}}(\underline{\mathit{k}})|}^2.\end{array}$$

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

$$\begin{array}{c}\hfill C_{\mathit{\theta },\mathit{r}}(\underline{\mathit{k}})=e^{i\omega \mathit{\xi }·\mathit{r}}\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x)-\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x).\end{array}$$

We then obtain

$$\begin{array}{c}\hfill |C_{\mathit{\theta },\mathit{r}}{(\underline{\mathit{k}})|}^2=|\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x){|}^2+{|\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x)|}^2-2\times \mathrm{Re}\left[e^{i\omega \mathit{\xi }·\mathit{r}}{\stackrel{\sim }{{\chi }^{\mathit{\theta }}}}^{\ast }(k_x)\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x)\right].\end{array}$$

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

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

$$\begin{array}{c}\hfill {{\displaystyle \int }}_p|\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x){|}^2+|\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x){|}^2=2\pi {{\displaystyle \int }}_{\mathit{\theta }}{|\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)|}^2=2\pi {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }{P}_{\rho }(k_x,\mathit{\xi })\mathrm{d}\mathit{\xi }.\end{array}$$

Using Eq. (C.2) and some algebra leads to

$$\begin{array}{c}\hfill {\stackrel{\sim }{{\chi }^{\mathit{\theta }}}}^{\ast }(k_x)\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x)={\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }\mathrm{d}{\mathit{\xi }}^{″}{\stackrel{\sim }{\rho }}^{\ast }(k_x,{\mathit{\xi }}^{\prime })\stackrel{\sim }{\rho }(k_x,{\mathit{\xi }}^{″})e^{i\omega ({\mathit{\xi }}^{\prime }-{\mathit{\xi }}^{″})·\mathit{\theta }}{e}^{-i\omega {\mathit{\xi }}^{″}·\mathit{r}}.\end{array}$$

Summing the term under the Re function over all pairs (without double counting), we write

$$\begin{array}{cc}\hfill {{\displaystyle \int }}_p\left[e^{i\omega \mathit{\xi }·\mathit{r}}{\stackrel{\sim }{{\chi }^{\mathit{\theta }}}}^{\ast }(k_x)\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x)\right]& =\frac{1}{2}{{\displaystyle \int }}_{\mathit{\theta }}{{\displaystyle \int }}_{|\mathit{r}|=s}\left[e^{i\omega \mathit{\xi }·\mathit{r}}{\stackrel{\sim }{{\chi }^{\mathit{\theta }}}}^{\ast }(k_x)\stackrel{\sim }{{\chi }^{\mathit{\theta }+\mathit{r}}}(k_x)\right]\hfill \\ \hfill & =\frac{1}{2}{\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }{P}_{\rho }(k_x,{\mathit{\xi }}^{\prime }){{\displaystyle \int }}_0^{2\pi }{e}^{i\omega |\mathit{\xi }+{\mathit{\xi }}^{\prime }|scos\varphi }\mathrm{d}\varphi \hfill \\ \hfill & =\pi {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }{P}_{\rho }(k_x,{\mathit{\xi }}^{\prime })J_0(\omega |\mathit{\xi }+{\mathit{\xi }}^{\prime }|s).\hfill \end{array}$$

We finally obtain

$$\begin{array}{c}\hfill I_s(\underline{\mathit{k}})=2\pi {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }{P}_{\rho }(k_x,{\mathit{\xi }}^{\prime })[1-J_0\left(|\mathit{\xi }+{\mathit{\xi }}^{\prime }|\omega s\right)].\end{array}$$(D.2)

Dividing by the total number of pairs $N_p^{\mathrm{tot}}(s)\simeq \frac{1}{2}{N}_p(\mathit{\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. (12)). These calculations assume that integration over all pairs (θ, r) is continuous. Appendix E describes the effect of pixelised and filtered data.

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

Previous calculations neglect border effects in the integration over pairs of points. We have set ρ = 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 ⟨|C(θ + r)−C(θ)|²⟩=⟨|C(θ)|²⟩ in the numerator I_s and the number of pairs entering the denominator N_p(s) needs to be corrected (see Fig. D.1).

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^{\text{tot}}(s)$ such pairs and $N_p^{\text{ext}}(s)$ pairs with only one end within the field of view. Naturally we denote $N_p^{\text{in}}=N_p^{\text{tot}}-N_p^{\text{ext}}$ the pairs fully comprised within the domain analysis. A given point θ in the analysis domain belongs to ϕ_s(θ)∈[0, 2π] pairs in the field of view. Let us first compute the exact number of pairs, assuming an infinitely fine tessellation,

$$\begin{array}{c}\hfill N_p^{\mathrm{tot}}(s)=\pi {\mathcal{S}}_{\mathcal{A}}+\frac{1}{2}{N}_p^{\mathrm{ext}}(s).\end{array}$$

We write

$$\begin{array}{cc}\hfill (N_p^{\mathrm{tot}}(s)-\frac{1}{2}{N}_p^{\mathrm{ext}}(s))\mathrm{sf}(s)& ={{\displaystyle \int }}_p⟨.⟩={{\displaystyle \int }}_{\mathrm{in}}{⟨.⟩}^{\mathrm{in}}+{{\displaystyle \int }}_{\mathrm{ext}}{⟨.⟩}^{\mathrm{ext}}\hfill \\ \hfill & =N_p^{\mathrm{in}}{\mathrm{sf}}^{\mathrm{corr}}(s)+{{\displaystyle \int }}_{\mathrm{ext}}{⟨.⟩}^{\mathrm{ext}},\hfill \end{array}$$

denoting by ${\text{sf}}^{\text{corr}}(s)=1/N_p^{\text{in}}{{\displaystyle \int }}_{\text{in}}\langle .\rangle$ the value of the structure function corrected from finite-sized effects. We obtain

$$\begin{array}{c}\hfill {{\displaystyle \int }}_{\mathrm{ext}}{⟨.⟩}^{\mathrm{ext}}={{\displaystyle \int }}_{\mathit{\theta }\in \mathcal{A}}(2\pi -{\varphi }_s(\mathit{\theta })){⟨|C(\mathit{\theta })|}^2⟩\mathrm{d}\mathit{\theta }.\end{array}$$

As demonstrated in Sect. 2,

$$\begin{array}{c}\hfill {⟨|C(\mathit{\theta })|}^2⟩={\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(\underline{\mathit{k}}){\left|\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)\right|}^2.}\end{array}$$

Reassembling terms, we obtain the corrected mean structure function (sf^corr):

$$\begin{array}{cc}\hfill {\mathrm{sf}}^{\mathrm{corr}}(s)=& (\frac{N_p^{\mathrm{ext}}(s)}{2N_p^{\mathrm{in}}(s)}+1)\mathrm{sf}(s)\hfill \\ \hfill & -\frac{1}{N_p^{\mathrm{in}}(s)}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(\underline{\mathit{k}}){{\displaystyle \int }}_{\mathit{\theta }\in \mathcal{A}}(2\pi -{\varphi }_s(\mathit{\theta })){\left|\stackrel{\sim }{{\chi }^{\mathit{\theta }}}(k_x)\right|}^2\mathrm{d}\mathit{\theta }.}\hfill \end{array}$$(D.3)

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

We have

$$\begin{array}{cc}& N_p^{\mathrm{in}}(s)=\frac{1}{2}{{\displaystyle \int }}_{\mathit{\theta }\in \mathcal{A}}{\varphi }_s(\mathit{\theta })\mathrm{d}\mathit{\theta },\hfill \\ \hfill & N_p^{\mathrm{ext}}(s)={{\displaystyle \int }}_{\mathit{\theta }\in \mathcal{A}}(2\pi -{\varphi }_s(\mathit{\theta }))\mathrm{d}\mathit{\theta }.\hfill \end{array}$$

Estimating ϕ is easy under the assumption of a circular analysis region of radius R. We find the following expressions, graphically represented in Fig. D.2, left:

$$\begin{array}{c}\hfill {\varphi }_s(\mathit{\theta })=\{\begin{array}{cc}2\pi & \mathrm{i}f\theta <R-sands<R\hfill \\ 0& \mathrm{i}f\theta <s-Rands>R\hfill \\ 2\mathrm{Arccos}\left(\frac{{\theta }^2+s^2-R^2}{2\theta s}\right)& \mathrm{i}f(R-s<\theta <Rands<R)\hfill \\ & \mathrm{o}r(s-R<\theta <Rands>R).\hfill \end{array}\end{array}$$

We therefore rewrite $N_p^{\text{in}}(s)=s^{2}F(R/s)$ and $N_p^{\text{out}}(s)=s^{2}G(R/s)$ with

$$\begin{array}{c}\hfill F(x)=\{\begin{array}{cc}{\pi }^2{(x-1)}^2+2\pi {{\displaystyle \int }}_{x-1}^x\gamma (u;x)u\mathrm{d}u& \mathrm{i}fx>1\hfill \\ 2\pi {{\displaystyle \int }}_{1-x}^x\gamma (u;x)u\mathrm{d}u& \mathrm{i}f0.5<x<1\hfill \end{array}\end{array}$$

and

$$\begin{array}{c}\hfill G(x)=\{\begin{array}{cc}4\pi {{\displaystyle \int }}_{x-1}^x[\pi -\gamma (u;x)]u\mathrm{d}u& \mathrm{i}fx>1\hfill \\ 2{\pi }^2{(1-x)}^2+4\pi {{\displaystyle \int }}_{1-x}^x[\pi -\gamma (u;x)]u\mathrm{d}u& \mathrm{i}f0.5<x<1\hfill \end{array}\end{array}$$

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 in Fig. D.2. As expected, the number of extra pairs is negligible for small pair separations. It equals the number of regular (“inner”) pairs for s ≃ 0.5R and becomes dominant past this value.

Finally, we note that for the flat emissivity field (ρ(x, θ)=ϵ(x)/F) we have $|\stackrel{\sim }{{\chi }^{\mathit{\theta }}}{|}^2=P_{ϵ}/F^2$ and then

$$\begin{array}{c}\hfill {\mathrm{sf}}^{\mathrm{corr}}(s)=(\frac{N_p^{\mathrm{ext}}(s)}{2N_p^{\mathrm{in}}(s)}+1)\mathrm{sf}(s)-\frac{N_p^{\mathrm{ext}}(s)}{N_p^{\mathrm{in}}(s)}\times {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \sum _{\mathit{\xi }}{P}_{2\mathrm{D}}^{\infty }(\mathit{\xi }).}\end{array}$$(D.4)

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_{2\mathrm{D}}^{\infty }$ in Eq. (12), instead of involving P_2D and applying this correction formula. This property is used to check the validity of the correction formula in Eq. (D.4). Figure D.3 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 with a radius of 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 𝒜. As the field of view increases (R = 500 kpc), border effects become negligible and all structure functions match the exact one. Finally, as $N_p^{\text{in}}$ approaches zero for a separation length of size s = 2R, the correction formula becomes numerically unstable at large separation lengths.

Appendix E: Structure function from pixelised 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 effectively groups lines of sight within a single spectral line measurement. Both processes are formally close to each other, since pixelisation 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. E.1).

We define a new map D(θ) as

$$\begin{array}{c}\hfill D=\frac{{\mathcal{F}}_{\ell }\ast \left(FC\right)}{{\mathcal{F}}_{\ell }\ast F},\end{array}$$(E.1)

where represents the usual convolution product, and F(θ) and C(θ) are respectively the flux and centroid maps as defined in Sect. 3.1. The filter ℱ_ℓ may represent the instrumental point spread function, or the pixel window function (previously denoted 𝒲) or a combination of the two; we assume its characteristic scale is ℓ (e.g. instrument full width at half maximum, pixel size, etc.) and it is normalised to 1 by integrating over all values of θ. It is clear that D(θ) 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(θ)≃C(θ) for components of C varying on scales much larger than ℓ (equivalently, for very sharp filters). For components of C oscillating on tiny scales (much smaller than the filter size), we find D ≃ 0: as expected the pixelisation or filtering process suppresses information on small scales.

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

$$\begin{array}{c}\hfill C_{{\mathit{\theta }}_0,\mathit{r}}\to {\displaystyle \int }\mathrm{d}\mathit{\mu }{\mathcal{F}}_{\ell }(\mathit{\mu })e^{i\omega \mathit{\mu }·\mathit{\xi }}{C}_{{\mathit{\theta }}_0-\mathit{\mu },\mathit{r}}.\end{array}$$

The calculation steps are similar to the previous case, thanks to permutations of the integrals over μ and other integrals. This calculation leads to the expected value of the structure function at each s:

$$\begin{array}{cc}\hfill \mathrm{sf}(s)& =2{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){\displaystyle \int }\mathrm{d}{\mathit{\xi }}^{\prime }{P}_{\rho }(k_x,{\mathit{\xi }}^{\prime })P_{\ell }(\mathit{\xi }+{\mathit{\xi }}^{\prime })[1-J_0\left(|\mathit{\xi }+{\mathit{\xi }}^{\prime }|\omega s\right)]}\hfill \\ \hfill & =2{\displaystyle \int }[1-J_0\left(\omega \left|\mathit{\xi }\right|s\right)]P_D(\mathit{\xi })\mathrm{d}\mathit{\xi },\hfill \end{array}$$

where P_D = P_ℓP_2D is the power spectrum of the map D(θ) and P_ℓ 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. (12) the power spectrum of the centroid map, P_2D, with the power spectrum of the filtered centroid map, P_D. In the case of larger pixels, this is generally no longer valid. This is critical in the presence of a finite domain of analysis of a size comparable to the pixel size, since then border effects must be treated more carefully. Figure E.2 shows the result of applying the simple prescription P_2D → P_ℓP_2D to Eq. (D.4), with a similar parametric setup as in Fig. D.3. Since in this case the velocity field is stationary, we also have an exact computation of the structure function obtained by neglecting the finite-sized domain, that is, by using $P_{\ell }{P}_{2\mathrm{D}}^{\infty }$ in Eq. (12). It is then obvious that, strictly speaking, the correction formula in Eq. (D.4) is valid only for unbinned data.

Finally, in addition to this “smoothing” effect, pixelisation induces a discretization effect, or “aliasing”. We do not develop a calculation for this effect. For a given two-dimensional frequency ωξ, 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. This aliasing effect 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

$$\begin{array}{c}\hfill U_i={|C({\mathit{\theta }}_i+{\mathit{r}}_i)-C\left({\mathit{\theta }}_i\right)|}^2={\displaystyle \sum _{\underline{\mathit{k}},{\underline{\mathit{k}}}^{\mathbf{\prime }}}{V}_{\underline{\mathit{k}}}{V}_{{\underline{\mathit{k}}}^{\mathbf{\prime }}}{e}^{i\omega (\mathit{\xi }+{\mathit{\xi }}^{\prime })·{\mathit{\theta }}_i}{C}_{{\mathit{\theta }}_i,{\mathit{r}}_i}(\underline{\mathit{k}})C_{{\mathit{\theta }}_i,{\mathit{r}}_i}({\underline{\mathit{k}}}^{\mathbf{\prime }}).}\end{array}$$

We note in particular that

$$\begin{array}{c}\hfill ⟨U_i⟩={\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k){\left|C_{{\mathit{\theta }}_i,{\mathit{r}}_i}(\underline{\mathit{k}})\right|}^2.}\end{array}$$(F.1)

Even in the 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

$$\begin{array}{c}\hfill \langle \mathrm{SF}(s_i)\mathrm{SF}(s_j)\rangle =\langle \frac{1}{n_i{n}_j}{{\displaystyle \int }}_{{\mathit{\theta }}_i,{\mathit{r}}_i,{\mathit{\theta }}_j,{\mathit{r}}_j}{U}_i{U}_j\rangle =\frac{1}{n_i{n}_j}{{\displaystyle \int }}_{{\mathit{\theta }}_i,{\mathit{r}}_i,{\mathit{\theta }}_j,{\mathit{r}}_j}⟨U_i{U}_j⟩.\end{array}$$

Neglecting border effects, n_i = n_j ∝ 2π𝒮_𝒜. Therefore one writes

$$\begin{array}{cc}\hfill ⟨U_i{U}_j⟩& ={\displaystyle \sum _{\underline{\mathit{k}},{\underline{\mathit{k}}}^{\mathbf{\prime }},\underline{\mathit{l}},{\underline{\mathit{l}}}^{\mathbf{\prime }}}\langle V_{\underline{\mathit{k}}}{V}_{{\underline{\mathit{k}}}^{\mathbf{\prime }}}{V}_{\underline{\mathit{l}}}{V}_{{\underline{\mathit{l}}}^{\mathbf{\prime }}}\rangle e^{i\omega (\mathit{\xi }+{\mathit{\xi }}^{\prime })·{\mathit{\theta }}_i}{e}^{i\omega (\mathit{\chi }+{\mathit{\chi }}^{\prime })·{\mathit{\theta }}_j}}\hfill \\ \hfill & C_{{\mathit{\theta }}_i,{\mathit{r}}_i}(k_x,\mathit{\xi })C_{{\mathit{\theta }}_i,{\mathit{r}}_i}(k_x^{\prime },{\mathit{\xi }}^{\prime })C_{{\mathit{\theta }}_j,{\mathit{r}}_j}(l_x,\mathit{\chi })C_{{\mathit{\theta }}_j,{\mathit{r}}_j}(l_x^{\prime },{\mathit{\chi }}^{\prime }).\hfill \end{array}$$

We stress that |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. In a similar way to the one-dimensional case, we assume that

$$\begin{array}{c}\hfill \langle V_{\underline{\mathit{j}}}{V}_{\underline{\mathit{k}}}{V}_{\underline{\mathit{l}}}{V}_{\underline{\mathit{m}}}\rangle =\{\begin{array}{cc}{P}_{3\mathrm{D}}(k)P_{3\mathrm{D}}(l)& \mathrm{i}f(\underline{\mathit{k}}=-\underline{\mathit{j}});(\underline{\mathit{l}}=-\underline{\mathit{m}});(\underline{\mathit{k}}\ne \pm \underline{\mathit{l}})\{A\}\hfill \\ P_{3\mathrm{D}}(k)P_{3\mathrm{D}}(j)& \mathrm{i}f(\underline{\mathit{k}}=-\underline{\mathit{l}});(\underline{\mathit{j}}=-\underline{\mathit{m}});(\underline{\mathit{k}}\ne \pm \underline{\mathit{j}})\{B\}\hfill \\ P_{3\mathrm{D}}(k)P_{3\mathrm{D}}(j)& \mathrm{i}f(\underline{\mathit{k}}=-\underline{\mathit{m}});(\underline{\mathit{j}}=-\underline{\mathit{l}});(\underline{\mathit{k}}\ne \pm \underline{\mathit{j}})\{C\}\hfill \\ \langle |V_{\underline{\mathit{k}}}{|}^4\rangle & \mathrm{i}f(\underline{\mathit{k}}=-\underline{\mathit{j}}=\underline{\mathit{l}}=-\underline{\mathit{m}})\{D\}\hfill \\ \langle |V_{\underline{\mathit{k}}}{|}^4\rangle & \mathrm{i}f(\underline{\mathit{k}}=-\underline{\mathit{j}}=-\underline{\mathit{l}}=\underline{\mathit{m}})\{E\}\hfill \\ \langle |V_{\underline{\mathit{k}}}{|}^4\rangle & \mathrm{i}f(\underline{\mathit{k}}=\underline{\mathit{j}}=-\underline{\mathit{l}}=-\underline{\mathit{m}})\{F\}\hfill \\ 0& \mathrm{else}\hfill \end{array}\end{array}$$

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 the case in which $\underline{\mathit{k}}=-\underline{\mathit{k}}\mathbf{\prime }$, $\underline{\mathit{l}}=-\underline{\mathit{l}}\mathbf{\prime }$ and $\underline{\mathit{k}}\ne \pm \underline{\mathit{l}}$. It is written

$$\begin{array}{cc}\hfill b_A& ={(n_i{n}_j)}^{-1}{\displaystyle \sum _{\underline{\mathit{k}}\ne \pm \underline{\mathit{l}}}{P}_{3\mathrm{D}}(k)P_{3\mathrm{D}}(l){\displaystyle \int }{\left|C_{{\mathit{\theta }}_i,{\mathit{r}}_i}(\underline{\mathit{k}})\right|}^2{\left|C_{{\mathit{\theta }}_j,{\mathit{r}}_j}(\underline{\mathit{l}})\right|}^2}\hfill \\ \hfill & ={(n_i{n}_j)}^{-1}\left[{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}(k)I_{s_i}(\underline{\mathit{k}})}\right]\left[{\displaystyle \sum _{\underline{\mathit{l}}}{P}_{3\mathrm{D}}(l)I_{s_j}(\underline{\mathit{l}})}\right]\hfill \\ \hfill & -2{(n_i{n}_j)}^{-1}{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}{(k)}^2{I}_{s_i}(\underline{\mathit{k}})I_{s_j}(\underline{\mathit{k}}).}\hfill \end{array}$$

with I_s already defined in Appendix D. Therefore, the first term is simply sf(s_i)sf(s_j).

Computation of b_B: It corresponds to the case in which $\underline{\mathit{k}}=-\underline{\mathit{l}}$, $\underline{\mathit{k}}\mathbf{\prime }=-\underline{\mathit{l}}\mathbf{\prime }$ and $\underline{\mathit{k}}\ne \pm \underline{\mathit{k}}\mathbf{\prime }$. It is written

$$\begin{array}{cc}\hfill b_B& ={(n_i{n}_j)}^{-1}{\displaystyle \sum _{\underline{\mathit{k}}\ne \pm {\underline{\mathit{k}}}^{\mathbf{\prime }}}{P}_{3\mathrm{D}}(k)P_{3\mathrm{D}}(k^{\prime }){\displaystyle \int }{e}^{i\omega (\mathit{\xi }+{\mathit{\xi }}^{\prime })·({\mathit{\theta }}_i-{\mathit{\theta }}_j)}{C}_{{\mathit{\theta }}_i,{\mathit{r}}_i}(k_x,\mathit{\xi })}\hfill \\ \hfill & C_{{\mathit{\theta }}_i,{\mathit{r}}_i}(k_x^{\prime },{\mathit{\xi }}^{\prime })C_{{\mathit{\theta }}_j,{\mathit{r}}_j}^{\ast }(k_x,\mathit{\xi })C_{{\mathit{\theta }}_j,{\mathit{r}}_j}^{\ast }(k_x^{\prime },{\mathit{\xi }}^{\prime })\hfill \\ \hfill & ={(n_i{n}_j)}^{-1}({\displaystyle \sum _{\underline{\mathit{k}},\underline{\mathit{l}}}{P}_{3\mathrm{D}}(k)P_{3\mathrm{D}}(l)J_{s_i}(\underline{\mathit{k}},\underline{\mathit{l}})J_{s_j}^{\ast }(\underline{\mathit{k}},\underline{\mathit{l}})}\hfill \\ \hfill & -{\displaystyle \sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}{(k)}^2{J}_{s_i}(\underline{\mathit{k}},\underline{\mathit{k}})J_{s_j}^{\ast }(\underline{\mathit{k}},\underline{\mathit{k}})-\sum _{\underline{\mathit{k}}}{P}_{3\mathrm{D}}{(k)}^2{I}_{s_i}(\underline{\mathit{k}})I_{s_j}(\underline{\mathit{k}})}).\hfill \end{array}$$

This involves the function

$$\begin{array}{c}\hfill J_s(\underline{\mathit{k}},\underline{\mathit{l}})={{\displaystyle \int }}_{\mathit{\theta },|\mathit{r}|=s}{e}^{i\omega (\mathit{\xi }+\mathit{\chi })·\mathit{\theta }}{C}_{\mathit{\theta },\mathit{r}}(\underline{\mathit{k}})C_{\mathit{\theta },\mathit{r}}(\underline{\mathit{l}})\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{r}.\end{array}$$

By developing the expression of C_θ, r and using Parseval’s theorem, we find

$$\begin{array}{c}\hfill J_s(\underline{\mathit{k}},\underline{\mathit{l}})=4\pi {\left(\frac{\omega }{2\pi }\right)}^2{\displaystyle \int }\stackrel{\sim }{\rho }(k_x,\mathit{\xi }+\mathit{\kappa })\stackrel{\sim }{\rho }(l_x,\mathit{\chi }-\mathit{\kappa })(1-J_0(\omega \kappa s))\mathrm{d}\mathit{\kappa }.\end{array}$$(F.2)

Computation of b_D: It corresponds to $\underline{\mathit{k}}=-\underline{\mathit{k}}\mathbf{\prime }=\underline{\mathit{l}}=-\underline{\mathit{l}}\mathbf{\prime }$. It is written

$$\begin{array}{c}\hfill b_D={(n_i{n}_j)}^{-1}{\displaystyle \sum _{\underline{\mathit{k}}}\langle {\left|V_{\underline{\mathit{k}}}\right|}^4\rangle I_{s_i}(\underline{\mathit{k}})I_{s_j}(\underline{\mathit{k}}).}\end{array}$$

Computation of b_F: It corresponds to $\underline{\mathit{k}}=\underline{\mathit{k}}\mathbf{\prime }=-\underline{\mathit{l}}=-\underline{\mathit{l}}\mathbf{\prime }$. It is written

$$\begin{array}{c}\hfill b_F={(n_i{n}_j)}^{-1}{\displaystyle \sum _{\underline{\mathit{k}}}\langle {\left|V_{\underline{\mathit{k}}}\right|}^4\rangle J_{s_i}(\underline{\mathit{k}},\underline{\mathit{k}})J_{s_j}^{\ast }(\underline{\mathit{k}},\underline{\mathit{k}}).}\end{array}$$

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

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}=\langle \mathrm{SF}(s_i)\mathrm{SF}(s_j)\rangle -\langle \mathrm{SF}(s_i)\rangle \langle \mathrm{SF}(s_j)\rangle ,\end{array}$$

and which is written

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}& =\frac{1}{{(2\pi {\mathcal{S}}_{\mathcal{A}})}^2}[2{\displaystyle \sum _{\underline{\mathit{k}},\underline{\mathit{l}}}{P}_{3\mathrm{D}}(k)P_{3\mathrm{D}}(l)J_{s_i}(\underline{\mathit{k}},\underline{\mathit{l}})J_{s_j}^{\ast }(\underline{\mathit{k}},\underline{\mathit{l}})}\hfill \\ \hfill & -{\displaystyle \sum _{\underline{\mathit{k}}}{R}_{\underline{\mathit{k}}}\times \{2I_{s_i}(\underline{\mathit{k}})I_{s_j}(\underline{\mathit{k}})+J_{s_i}(\underline{\mathit{k}},\underline{\mathit{k}})J_{s_j}^{\ast }(\underline{\mathit{k}},\underline{\mathit{k}})\}}].\hfill \end{array}$$(F.3)

The expressions for I and J are given by Eqs. (D.2) and (F.2), respectively. Notably, the second term within brackets vanishes for Rayleigh-distributed coefficients. In principle, finite-sized corrections must also apply, like 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}}={|\hat{\mathcal{A}}|}^2,$ its power spectrum, we obtain $\stackrel{\sim }{\rho }(k_x,\mathit{\xi })=\stackrel{\sim }{ϵ}(k_x)\hat{\mathcal{A}}(\mathit{\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:

$$\begin{array}{cc}& {\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },\mathit{\chi };s)=K{\displaystyle \int }\hat{\mathcal{A}}(\mathit{\xi }+\mathit{\kappa })\hat{\mathcal{A}}(\mathit{\chi }-\mathit{\kappa })(1-J_0(\omega \kappa s))\mathrm{d}\mathit{\kappa }\hfill \\ \hfill & {\mathcal{T}}_{\mathcal{A}}(\mathit{\xi };s)={\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },-\mathit{\xi };s)=K{\displaystyle \int }{P}_{\mathcal{A}}(\mathit{\kappa })(1-J_0(\omega |\mathit{\xi }+\mathit{\kappa }|s))\mathrm{d}\mathit{\kappa }.\hfill \end{array}$$

The normalisation constant is K = (ω/2π)²/𝒮_𝒜 throughout this paper.

Using these functions, the variance is then written

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}& =8{\left(\frac{\omega }{2\pi }\right)}^4{\displaystyle \sum _{\mathit{\xi },\mathit{\chi }}{P}_{2\mathrm{D}}^{\infty }(\xi )P_{2\mathrm{D}}^{\infty }(\chi ){\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },\mathit{\chi };s_i){\mathcal{U}}_{\mathcal{A}}^{\ast }(\mathit{\xi },\mathit{\chi };s_j)}\hfill \\ \hfill & -4{\left(\frac{\omega }{2\pi }\right)}^4{\displaystyle \sum _{\mathit{\xi }}{Q}_{2\mathrm{D}}^{\infty }{(\xi )}^2\{2{\mathcal{T}}_{\mathcal{A}}(\mathit{\xi };s_i){\mathcal{T}}_{\mathcal{A}}(\mathit{\xi };s_j)}\hfill \\ \hfill & +{\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },\mathit{\xi };s_i){\mathcal{U}}_{\mathcal{A}}^{\ast }(\mathit{\xi },\mathit{\xi };s_j)\}.\hfill \end{array}$$(F.4)

where we introduce $P_{2\mathrm{D}}^{\infty }$ the 2D power spectrum of the centroid map over an infinitely extended domain (Appendix C) and

$$\begin{array}{c}\hfill Q_{2\mathrm{D}}^{\infty }{(\xi )}^2={\left(\frac{2\pi }{\omega }\right)}^4{\displaystyle \sum _{k_x}\frac{P_{ϵ}{(k_x)}^2}{F^4}{R}_{k_x,\mathit{\xi }}.}\end{array}$$

If the moduli are Rayleigh-distributed, it is evident that Q_2D = 0 and the expression for the variance (Eq. (F.4)) depends only on the 2D power spectrum of the velocity. As already noted, the terms encapsulating the field-of-view geometry (𝒰_𝒜, 𝒯_𝒜) 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 𝒜 → ∞ then applies and $\hat{\mathcal{A}}$ becomes a strongly peaked function around 0. The above functions are rewritten as

$$\begin{array}{cc}\hfill {\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },\mathit{\chi };s)& \simeq K(1-J_0(\omega \xi s)){\displaystyle \int }\hat{\mathcal{A}}(\mathit{\xi }+\mathit{\kappa })\hat{\mathcal{A}}(\mathit{\chi }-\mathit{\kappa })\mathrm{d}\mathit{\kappa }\hfill \\ \hfill & =K(1-J_0(\omega \xi s))\left\{(\hat{\mathcal{A}}\ast \hat{\mathcal{A}})(\mathit{\xi }+\mathit{\chi })\right\}.\hfill \end{array}$$

The sign indicates the convolution product. Since 𝒜² = 𝒜, we obtain

$$\begin{array}{c}\hfill {\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },\mathit{\chi };s)\simeq K(1-J_0(\omega \xi s)){\left(\frac{2\pi }{\omega }\right)}^2\hat{\mathcal{A}}(\mathit{\xi }+\mathit{\chi }).\end{array}$$

This automatically shows that 𝒯_𝒜(ξ; s)=1 − J₀(ωξs) and 𝒰_𝒜(ξ, ξ; s)=0. Moreover,

$$\begin{array}{cc}& {\displaystyle \sum _{\mathit{\chi }}{P}_{2\mathrm{D}}^{\infty }(\chi ){\mathcal{U}}_{\mathcal{A}}(\mathit{\xi },\mathit{\chi };s_i){\mathcal{U}}_{\mathcal{A}}^{\ast }(\mathit{\xi },\mathit{\chi };s_j)}\hfill \\ \hfill & \simeq {\left(\frac{2\pi }{\omega }\right)}^4{K}^2(1-J_0(\omega \xi s_i))(1-J_0(\omega \xi s_j)){\displaystyle \sum _{\mathit{\chi }}{P}_{\mathcal{A}}(\mathit{\chi }+\mathit{\xi })P_{2\mathrm{D}}^{\infty }(\chi )}\hfill \\ \hfill & ={\left(\frac{2\pi }{\omega }\right)}^2{\mathcal{S}}_{\mathcal{A}}^{-1}{P}_{2\mathrm{D}}^{\infty }(\xi )(1-J_0(\omega \xi s_i))(1-J_0(\omega \xi s_j)).\hfill \end{array}$$

Grouping terms together in Eq. (F.4) leads to Eq. (13).

Appendix G: Fourier transform of the normalised emissivity field ρ (spherical β-model)

We provide useful calculations for the 3D power spectrum of the normalised emissivity ρ in a case of a β-model. We also discuss the limiting case of very extended sources (equivalently, very small fields of view).

G.1. General case

Calculation of the two-dimensional power spectrum involves the calculation of $P_{\rho }={|\stackrel{\sim }{\rho }|}^2$ with ρ(x, θ)=ϵ(x, θ)/F(θ). As noted above (Appendix C), it is compulsory to fill ρ with zeros outside of its domain of definition or outside of the analysis domain. We consider here an arbitrarily large circular analysis domain 𝒜 (of radius R, centred on the source) to perform the following calculations.

Using the expression for F(θ) derived in a previous section (Eq. (B.1)), we obtain

$$\begin{array}{c}\hfill \rho (x,\mathit{\theta })=\frac{1}{r_c{u}_{\beta }}{(1+\frac{{\theta }^2}{r_c^2})}^{3\beta -1/2}{(1+\frac{x^2+{\theta }^2}{r_c^2})}^{-3\beta }·\end{array}$$

The Fourier transform $\stackrel{\sim }{\rho }(k_x,\mathit{\xi })$ is calculated in two steps. The integration over the x axis is performed first and its result is already displayed in Eq. (B.2). Integration over the second axis runs for all θ ∈ 𝒜 and we find

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }(k_x,\mathit{\xi })=\frac{2^{5/2-3\beta }\pi r_c^2}{\mathrm{\Gamma }(3\beta -1/2)}\times {\mathcal{H}}_{(R/r_c);(3\beta -1/2)}(\omega |k_x|r_c,\omega \xi r_c),\end{array}$$

which uses the special integral defined below and represented in Fig. G.1 for n = 3/2 (β = 2/3):

$$\begin{array}{c}\hfill {\mathcal{H}}_{p;n}(u,v)={{\displaystyle \int }}_0^{p}t{J}_0(vt){\mathcal{F}}_n\left(u\sqrt{1+t^2}\right)\mathrm{d}t,\end{array}$$

recalling that ℱ_n(x)=xⁿK_n(x).

G.2. Limit for small analysis domains (R ≪ r_c)

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

$$\begin{array}{c}\hfill P_{\rho }(k_x,\mathit{\xi })\simeq \frac{2^{3-6\beta }}{\mathrm{\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·\end{array}$$

The right-most factor involving the order 1 Bessel function J₁ is the power spectrum of a circular pupil of radius R. Its integral over ξ equals πR²(2π/ω)² and rapidly falls to zero for ω|ξ|≳3.83R⁻¹. This result can easily be generalised to an analysis domain of arbitrary shape, introducing its power spectrum P_𝒜,

$$\begin{array}{c}\hfill P_{\rho }(k_x,\mathit{\xi })\simeq \frac{P_{ϵ}(k_x;\theta =0)}{F{(0)}^2}\times P_{\mathcal{A}}(\mathit{\xi }),\end{array}$$(G.1)

and we naturally recover the limiting case discussed several times in this paper where the emissivity ϵ(x, θ)=ϵ(x) does not depend on the line of sight θ over the analysis domain. Finally, R(≪r_c) can still be very large and therefore $P_{\mathcal{A}}(\mathit{\xi })\to {\left(\frac{2\pi }{\omega }\right)}^2{\mathcal{S}}_{\mathcal{A}}\delta (\mathit{\xi })$: the component of P_ρ in the ξ 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 Figs. 5, 6, 9, and 10 for the two other sets of 100 simulations (Table 1). All parameters remain the same apart from the injection scales, respectively L_inj = 200 kpc (Figs. H.1–H.4) and L_inj = 300 kpc (Figs. H.5–H.8). These comparisons still demonstrate a good match between the simulations and analytic results.

All Tables

All Figures