Towards mapping turbulence in the intra-cluster medium II. Measurement uncertainties in the estimation of structure functions

X-ray observations of the hot gas ﬁlling the intra-cluster medium (ICM) provide a wealth of information on the dynamics of clusters of galaxies. The global equilibrium of the ICM is believed to be ensured by non-thermal and thermal pressure support sources, among which gas movements and the dissipation of energy through turbulent motions. Accurate mapping of turbulence using X-ray emission lines is challenging due to the lack of spatially resolved spectroscopy. Only future instruments such as the X-ray Integral Field Unit (X-IFU) on Athena will have the spatial and spectral resolution to quantitatively investigate the ICM turbulence over a broad range of spatial scales. Powerful diagnostics for these studies are line shift and the line broadening maps, and the second-order structure function. When estimating these quantities, instruments will be limited by uncertainties of their measurements, and by the sampling variance (also known as cosmic variance) of the observation. Here, we extend the formalism started in our companion Paper I to include the e ﬀ ect of statistical uncertainties of measurements in the estimation of these line diagnostics, in particular for structure functions. We demonstrate that statistics contribute to the total variance through di ﬀ erent terms, which depend on the geometry of the detector, the spatial binning and the nature of the turbulent ﬁeld. These terms are particularly important when probing the small scales of the turbulence. An application of these equations is performed for the X-IFU, using synthetic turbulent velocity maps of a Coma-like cluster. Results are in excellent agreement with the formulas both for the structure function estimation ( ≤ 3%) and its variance ( ≤ 10%). The expressions derived here and in Paper I are generic, and ensure an estimation of the total errors in any X-ray measurement of turbulent structure functions. They also open the way for optimisations in the upcoming instrumentation and in observational strategies.


Introduction
The X-ray emission of clusters of galaxies offers a phenomenal window to observe the thermodynamic and dynamic properties of the hot baryons composing the intra-cluster medium (ICM).The gas trapped in the dark matter potential of these structures holds an untouched fossil record of their formation, giving us a glimpse of the early days of the Universe (see Kravtsov & Borgani 2012;Planelles et al. 2015, for reviews).The first observations of the ICM showed smooth, spherical profiles, well described by β-models (Cavaliere & Fusco-Femiano 1978), suggesting that the gas could be considered in (or close to) hydrostatic equilibrium.Several subsequent X-ray missions have since demonstrated that the ICM is far from homogeneous (Fabian et al. 2006;Vikhlinin et al. 2006;Leccardi & Molendi 2008).Dynamics induced by constant 3D accretion from the medium surrounding the cluster and by merger events throughout their lifetime are strengthened by the role of central active galactic nuclei (AGNs), whose jets, outflows, and bubbles, drive powerful mechanical and radiative motions, stirring the ICM at every spatial scale (Fabian 2012;King & Pounds 2015;Gaspari & Sadowski 2017;Morganti 2017).Other effects present both at small (e.g.galaxy outflows) and large scales (e.g.sloshing, ram-stripping) also create heterogeneities in the gas emission, thereby severely questioning the assumption of hydrostaticity.
Hints of systematic deviations from hydrostatic equilibrium up to a 10 to 20% are indeed found in both state-ofthe-art numerical simulations of the ICM and observational mass measurements (see Pratt et al. 2019, for a review).Other thermal and non-thermal pressure support mechanisms are therefore called upon to compensate the cooling infall of the ICM towards the central parts of the clusters.The identification of the mechanisms responsible for such deviations are crucial to understand the overall equilibrium of clusters (see Werner et al. 2019, for a review), and to have unbiased estimations of their mass, which is key to precision cosmology with clusters.
The dissipation of kinetic energy through either bulk or turbulent motions within clusters is a likely candidate (Gaspari et al. 2018;Voit 2018).In the classical view of the ICM, bulk motions are driven on a full-cluster scale by mechanisms such as mergers or ram-stripping, while turbulence indicates smaller scale motions.Turbulent energy is injected at hundreds of kpc scales and transported through a vortex cascade down to tens of kpc (Donnert et al. 2018).The gas motions induced by the turbulent cascade create a non-thermal pressure support mechanism (Lau et al. 2009), while the subsequent dissipation of the kinetic energy through collisions, small-scale viscosity, and eddies releases heat to the environment, counteracting part of the cooling flows in the cluster core (Zhuravleva et al. 2014).Scale-independent assumptions of the turbulent eddies naturally result in power-law forms of the turbulent power spectrum, characterised by an injection and dissipation scale of the energy, along with the characteristic slope of the spectrum (Kolmogorov 1941a).The determination of the energy injection scale provides information on the dominant energy transport mechanism at the cluster scale, involved for instance in the circulation of metals from the interstellar medium to the ICM (Rebusco et al. 2006).Knowledge of the dissipation length provides instead insight on the viscosity of the ICM and collision mechanisms at small spatial scales (Schekochihin et al. 2009).
Yet, a direct observation of the ICM kinematics remains challenging.Random movements of gas particles related to turbulence create line shifts, induce further broadening of the line, and can add skewness in the projection of the natural line profile along the line-of-sight.The understanding and mapping of these effects therefore require the measurement of centroid, width, and shape of the emission lines with accuracies of a few tens of km/s over the full cluster scale for typical Fe K α lines (∼ 6.4 keV).
Most of the current generation of X-ray instruments cannot provide spatially resolved high-resolution spectroscopy to this level of accuracy.Insight on the ICM kinematics therefore relies on other physical quantities, such as the measurement of bulk motions using cold shock fronts (Markevitch & Vikhlinin 2007), or the investigation of surface brightness, temperature, and density fluctuations in nearby clusters (Churazov et al. 2003(Churazov et al. , 2012;;Zhuravleva et al. 2015).These results are nevertheless insufficient for a definitive understanding of the kinematic pressure support.A ground-breaking step forward was achieved with the soft X-ray spectrometer (SXS) onboard Hitomi (Takahashi et al. 2018).Despite its short lifetime, the SXS mapped for the first time the turbulent velocity of the Perseus cluster, showing a quiescent ICM with velocities ∼ 200 km/s (Hitomi Collaboration 2016, 2018).New results are expected with the X-ray Imaging and Spectroscopy Mission (XRISM, Ishisaki et al. 2018) and its instrument Resolve.However, a spatial mapping of turbulent velocity fields with accuracies of ∼ 10/20 km/s down to a few tens of kpc will require instruments such as the X-ray Integral Field Unit (X-IFU, Barret et al. 2016Barret et al. , 2018) ) on board the future X-ray observatory Athena (Nandra et al. 2013).The X-IFU will provide an unprecedented 2.5 eV spectral resolution below 7 keV with a spatial accuracy of 5 (over a 5 equivalent field-of-view), enabling turbulence measurements through line broadening and deformations of the natural line profile (Ettori et al. 2013).
With the advent of high-resolution X-ray spectroscopy, powerful line diagnostics can be used to investigate turbulent motions.These include the shift and broadening of a spectral line, and the computation of the structure function of the line-of-sight velocity field (Inogamov & Sunyaev 2003), related to the turbulent power spectrum of the ICM (Zhuravleva et al. 2012).Any measurement of these quantities will be limited by statistical uncertainties, linked to the observational set-up, and by the sampling variance (or 'cosmic' variance) of the observation, which refers to the intrinsic variations of a given diagnostic related to the small number of observations of a random process.A theoretical understanding of these effects could provide a significant step forward in our knowledge of the usual line diagnostics used to study the turbulence of the ICM.
An analytical treatment of the cosmic variance is provided in our companion paper I (Clerc et al.,sub.,hereafter CL19), and used here.It allows the fast computation of estimates of the cosmic variance uncertainties, without using iterative Monte-Carlo techniques.In this paper, we extend this approach to include the contribution of statistics to the overall error estimation of the usual line diagnostics (line shift, line broadening and structure function) in the case of turbulence in the optically-thin emitting plasma of clusters of galaxies.Starting from formalism developed in Zhuravleva et al. (2012);ZuHone et al. (2016), we derive in the first part of this paper (Sect.2) the errors associated to the previous line diagnostics, notably on the value and the variance of the estimated structure function.The formulas are generic, and remain valid for any level of statistical error obtained from measurements with an X-ray instrument.A specific application on the future X-IFU instrument is provided, on the basis of synthetic observations (Sect.3) and the comparison of their post-processed outcomes with the prediction from our developed formalism (Sect.4).The implications of these error formulas are then discussed, along with ways to estimate these contributions (Sect.5).Throughout this paper, we assume a Λ-CDM cosmology, with h = 0.72, Ω m = 0.24 and Ω Λ = 0.76.Bold, underlined letters x indicate 3D vectors, bold letters x indicate 2D vectors.In a 3D space mapped by a (x, y, z) orthonormal frame, x is taken as the line-of-sight direction and (y, z) = θ as the plane-of-sky coordinates, • indicates the average operator, || • || 2 the Euclidean norm and X an estimator of the quantity X.Other notations are consistent with paper I.

Line centroid and broadening
Two tools to investigate the gas motions projected along the line-of-sight are the line shift, δE, and the line broadening, Σ. Line shift is defined as the difference between the energy of the line in the inertial frame of the observer, E 0 , and the measured value.It can be related to either gas motions along the line-of-sight, or to the cosmological redshift z of the source.By noting E z the energy of the line in the frame of the source, and I l (E) the line profile, 1 δE along a given line-of-sight θ is defined as where F (θ) = I l (E, θ)dE is the integrated flux of the line.Correspondingly, δE can be expressed as a centroid velocity shift C of the projected line-of-sight component of the velocity field as (with c the speed of light): The broadening of a line is the dispersion around its centroid value, and can be expressed similarly using Small turbulent motions create shifts in the corresponding line centroids of the gas particles.Integrated over along the line-of-sight, these result in a broadening of the observed line.The velocity broadening is thus defined by where Σ is the broadening after subtraction of the instrument spectral resolution and other physical broadening effects (e.g.thermal broadening), assumed perfectly known here.The measurement of S provides insight on the velocity distribution along the line-of-sight, making it a tool widely used to study turbulence (Hitomi Collaboration 2018).

The structure function
Another line diagnostic tool for turbulence is the structure function.Its use in turbulence analysis originates from the early studies of turbulent motions in fluid dynamics (Kolmogorov 1941b) which is a particular case of the auto-covariance function of the velocity field, K v .Under the assumption of an isotropic velocity field, the second-order structure function of the 3D field v, SF 2 , can be expressed exclusively as a function of spatial separation s between two points in space, and is related to where we average over all points x ∈ R 3 separated by a distance ||r|| 2 = s.The measurement of SF 2 provides a view of the underlying turbulent velocity power spectrum through a 'modified' second order moment of the velocity field.Although the properties of a turbulent field are not fully characterised by its power spectrum, the properties of structure functions and their simple estimation in fluid dynamics explains the success of this approach in all turbulence-related subjects.Notably, SF 2 can be used to estimate the characteristic lengths of the turbulence (Miniati 2015).More generally, we can define the n th moment of the structure function (n ∈ N) as for a separation ||r|| 2 = s.In the following sections, SF indicates the second-order structure function, and D the first-order structure function, or incremental function.

Estimators and value: the influence of finite statistics
The measurement of a velocity shift or a velocity broadening is related to a choice of the line-of-sight.Similarly, the definition of the structure function is related to a spatial average, such that an exact value can only be accessed either by averaging over a large number of spatial data points, or -if we assume ergodicity -by averaging over the same area for a large number of realisation of the turbulent velocity field.When observing astrophysical sources only a finite number of pointings and a limited exposure time are possible.Hence, it is essential to distinguish, for any line diagnostic, between the true value and its estimation.

Definitions and estimators of the structure functions
In the rest of this study, we assume ergodicity and isotropy of the turbulence processes.For a given point in space x ∈ R 3 with a speed v in the referential of the observer, we only consider the velocity component along the line-of-sight, that is v(x)•e x = v(x), with no loss of generality due to isotropy.
In astrophysical observations, velocities can only be measured in the 2D space of the detector.Per pixel, the result will be the projection of the line-of-sight component of the velocity field modulated by emissivity effects.We define the subset S s ⊂ R 4 , which contains all the doublets in the plane (x, y) separated by exactly s in the sense of the Euclidean norm.By convention, S 0 = ∅.We also define Ss as the 'halved' subset, not counting for x ↔ y permutations.The cardinal of Ss is noted N p (s) and represents the number of evaluations of the spatial average at a separation s.Under these assumptions, the previous estimators can be transposed to their 2D projected equivalents, SF and D, using the centroid shift C: where • Ss is the average over the data points in Ss .In practice, for any pixel (or centre of a region of pixels) in an observation, SF and D are computed using the following estimators: where C is the estimator of the centroid shifts (i.e.actual measurement per pixel).In real data sets, only one, or a few realisations of these quantities will be computed.To determine whether these estimators are biased, one has to compute their expected value (in the statistical sense) and compare it to the real value.X-ray observations of the centroid shift and line broadening will be affected by sources of statistical and systematic errors, such that in every pixel or region x, C(x) = C(x)+δC(x) stat +δC(x) syst and S(x) = S(x)+δ S(x) stat + δ S(x) syst .The former is related to the exposure time of the observation, the latter to uncertainties in the calibration (energy scale and energy redistribution function) or in the fit.We assume no systematic error is present in the observations.The statistical error on each point is represented as a random variable, normally distributed and centred.We define σ stat,C and σ stat, S the standard deviation of the statistical error for C and S respectively (not necessarily equal).
This distribution is considered spatially independent (i.e.valid on any subset of pixels).We further assume that the bulk motion is perfectly known, such that its contribution can be systematically subtracted from the measurements.Any turbulent velocity field is thus considered as centred, with a dispersion σ turb .

Expected value of the velocity shift and broadening
The estimator of the centroid shift along a given line-ofsight is obtained directly from the measurements.The expected value of the velocity shift for a pixel (or region) x over multiple realisations of the same random process is simply (for a centred field): The corresponding variance is (Appendix A.1): where Var[C] is the intrinsic variance of the centroid shift of the projected line-of-sight component of the velocity field over different random observations, which is affected by emissivity (see CL19).
Similarly, the estimator of the line broadening is simply the measured broadening of the line.After subtraction of the instrumental effects and other physical effects, the estimator may be affected by the statistics of the measurements such that over multiple realisations of the velocity field (see Appendix A.1): where E[ S2 (x)] and Var[ S2 (x)] are respectively the expected value and the intrinsic variance of the line broadening of the projected line-of-sight component of the velocity field (see CL19).Statistics induce an additional broadening, which adds to the intrinsic variance through statistical terms and cross products.

Expected value of the structure function
In the case of the structure function an ampler analytical approach is needed to quantify the effect of limited statistics in the measurements.For many observations of the same random turbulent process, the expected value of SF(s) is: The development (see Appendix A) yields the biased expected value of SF(s) shown by ZuHone et al. (2016): Equation 16 -valid throughout this paper -shows that the measurement of the SF using a statistically inaccurate measurement of the turbulent velocity is systematically biased, regardless of the number of points used to derive the structure function.

Variance of the structure function
It is important to determine whether the variance of the structure function is also affected by systematic biases.The accurate knowledge of the errors is crucial to understand the measurements and to distinguish between SF-related quantities (e.g.injection or dissipation scales).The same approach as in Sect.2.3.3 is thereby extended to the variance of the estimator: Under the previous assumptions, the variance of the estimator is given by (see Appendix A for development): (18) where N nei (s) is the number of neighbours at a distance s of any given point (see Appendix A.2 for a mathematical definition) and Var indicates the variance of the quantity over multiple observations of the same random process.
In the absence of statistical error (i.e.σ stat,C = 0), we recover the intrinsic variance of the structure function, which can be determined using the approach presented in CL19.With statistics, we distinguish between three different terms.Each can be interpreted as follows: (1) Velocity field fluctuation term: The first term is related to the variance of the incremental function.This term provides a sense of the velocity fluctuations over the observational filed-of-view (FoV).If fluctuations are small (i.e.nearby pixels have similar velocities), the effect of a statistical error will be small.On the contrary, when pixel-to-pixel fluctuations are large, statistics will affect the computation of the estimator for a given observation of the turbulent velocity power spectrum.This term will therefore be small for dissipation scales larger than the pixel (or binned region) size, or large otherwise.(2) Structure function fluctuation term: This second term can be related to the uncertainty with which the structure function is computed when using turbulent velocities affected by statistical errors.Its value follows closely the shape of the 'true' structure function (i.e.not positively biased, as in Sect.2.3.3) and is therefore negligible at low spatial separations (where SF is small), but increases with s.This term becomes negligible when a large number of pairs is used to estimate the structure function.
(3) Statistical fluctuations term: This term is the sheer contribution of the statistics to the overall variance of the structure function estimator.It appears from the covariance terms of the velocity field and is associated to the topology of the detector through the neighbour term N nei (s).Its contribution is most important at low spatial separations, where velocities are similar, and large statistical errors may introduce biases in the structure function estimation.
Alternatively, the first two terms can be written under a single term (1)+(2), related to the intrinsic nature of the line-of-sight component of the velocity field (both scale with σ 2 stat,C ).Formally, Var[D] is linked to SF (see Appendix A) and both show similar properties, notably at large spatial scales.Whenever N p is large, terms (2) and (3), related to the number of regions used to evaluate SF, go to zero.The first term however, is intrinsically linked to the number of observations of the turbulent process, and remains nonzero, even for large N p .It can be interpreted as a cross product between the cosmic variance of the field and the statistics, which can only be determined through multiple observations of the random process.For this reason, in the case of a Gaussian field, we made the choice to separate it from the sheer contribution of SF (2).A verification of these formulas on simple test cases (e.g.constant velocity fields, Gaussian fields) yields excellent results.

Generation of synthetic turbulent velocity fields
Future micro-calorimeter instruments such as the X-IFU will provide the required spatially resolved high-resolution spectroscopy to measure line shifts or line broadening, thus setting constraints on the turbulent velocity fields of nearby clusters (Roncarelli et al. 2018).The typical measurable scales of turbulence span from the size of the FoV -or the mapped area in case of multiple pointings -to the size of the pixel.A study of the turbulent cascade over different spatial scales with X-ray instruments is therefore limited to nearby clusters, where kpc to Mpc scales are accessible with arcmin-like FoV and arcsec-like pixels (provided a sufficient angular resolution).
The assessment of turbulence in these objects, however, will be hindered by cosmic variance and further degraded by the limited statistics of the observations.We provide in the rest of this paper an application of the previous equations (Sect.2) in the case of synthetic X-IFU observations to validate these formulas and demonstrate the capabilities of the instrument.The generation of the turbulent velocity fields is inspired from ZuHone et al. ( 2016) (hereby Z16).Simulations are based on the official E2E simulator SIXTE (Wilms et al. 2014), and performed similarly to Cucchetti et al. (2018) (hereby C18).

Emission profile and turbulent power spectrum
Forecasted targets to investigate the ICM turbulence with the X-IFU are local and massive clusters, such as the Perseus and Coma clusters.Hereafter, we consider a Comalike cluster, as in Z16, assuming an emission profile described by a β-model: where n e,0 is the electron density at the core and r c the 'core' radius of the cluster.In the rest of this paper we assume n e,0 = 3 × 10 −3 cm −3 , r c = 400 kpc, β = 2/3, and n e = 1.2 n H .We consider observations of the core of the cluster, where emissivity varies but slightly (∼ 2%) over an X-IFU FoV (i.e. 5 in equivalent diameter).At the redshift of Coma (z 0 = 0.023), 1 kpc corresponds to 2.21 on the sky.For simplicity, we assume an isothermal cluster at k B T = 7 keV, with a constant metallicity Z = 0.7 Z (see Ettori et al. 2015, abundances are given with respect to Anders & Grevesse 1989).The x direction remains the line-of-sight of the observations and we chose the centre of the cluster as the origin of a (x, y, z) orthonormal frame.For a point r, its 3D wave-vector is k = (k x , k y , k z ), with (k y , k z ) = ξ.Each gas particle in the ICM is simulated with a velocity v(r) along the line-of-sight.A full description of turbulence requires hydrodynamical treatments (Gaspari & Churazov 2013;Gaspari et al. 2014).We simplify here this approach (also for computational reasons) by assuming that turbulence follows an isotropic Kolmogorov 3D power spectrum where ṽ is the 3D Fourier transform of the velocity field along line-of-sight, C n is a normalisation factor of the power spectrum and k = ||k|| 2 .We note k dis , k inj the dissipation and injection scale respectively, and α the turbulent powerslope (Figure 1 -Left).

Normalisation of the power spectrum
Given the cluster emission profile, the velocity measured by the instrument will be convolved with the power spectrum of the cluster emission.As shown in Z16, if we note ∝ n e n H Λ(T, Z) the X-ray volume emissivity, the emission-measure weighted projection of the line-of-sight component of the velocity field along a given line-of-sight θ for a Gaussian or Voigt line is simply As the emission is isothermal and isometallic, only depends on the squared density of the cluster.By calling the normalised X-ray emissivity, the velocity dispersion S along a specific line-of-sight is given by The expected value of S2 along θ is related to the turbulent power spectrum (Zhuravleva et al. 2012, CL19) by: where P θ ρ is the 1D power spectrum of the normalised Xray emissivity ρ (Equation 22) for a fixed θ (depends on the selected line-of-sight, see Figure 1  The normalisation C n can be computed through numerical integration of Equation 24 for θ = 0.In our study, the computation is performed using a uniform 8192 × 8192 × 8192 grid of a length k = 0.50 kpc −1 , which corresponds to the spatial scale of pixel.To accurately compute the small scales in the centre of the cluster (notably for P θ ρ ), we take k min = (50r c ) −1 kpc −1 (Figure 1 -Right).This approach yields excellent results with respect to the purely analytical approach (possible with a β-model) with accuracies better than 0.1%.

Generation of the turbulent velocity field
One realisation of the turbulent velocity field is generated under the previous assumptions using the 3D turbulent power spectrum.We operate with a uniform area on the sky of 8.5 × 8.5 (i.e.230 kpc × 230 kpc) and 1.84 Mpc along the line-of-sight.The θ plane size is chosen to include more than 1.5 times the X-IFU FoV to avoid finite box-size effects of the simulation, while the grid is extended over the line-of-sight (8 times larger) to account more accurately for projection effects and ensure a smoother cut-off of the emissivity at the edges of the grid.For each run, we take a 2048 × 256 × 256 mesh, with a step size of /2 = 2.14 (0.97 kpc), which corresponds to the half-width of the X-IFU pixel to avoid aliasing (Shannon criterion), and offers a good computational speed compromise.
Each grid point is given a turbulent velocity in Fourier space ṽ(k) = |V k |e iψ , with |V k | the modulus, and ψ the phase, assumed without spatial correlation.As in Z16, we use a Rayleigh distribution of parameter Σ V k = P 3D (k)/2 for the modulus, and a spatially independent uniform distribution of the phase ψ.The corresponding probability distribution function being Without loss of generality, nor influence on the power spectrum, phases are computed exclusively on the lower triangular matrix of the velocity, and transposed to the upper triangular matrix to obtain a Hermitian velocity grid.Once computed in Fourier space, the velocity grid is transposed into real space to obtain the line-of-sight component of the velocity v(r) in each point in the grid.Given the large arrays the inverse Fourier transform is performed through 2DECOMP&FFT2 (Li & Laizet 2010), which is memory-optimised for large matrices.An example of the emission-measure weighted projection of a simulated velocity grid v(r) is provided Figure 2 (Top left).

Simulation set-up
The generated velocity grids can be used to validate the previous formulas in the case of X-IFU observations.They were therefore used as inputs to perform E2E simulations of our toy model Coma cluster.

Particle emission model
Similarly to C18, each gas particle is associated with a grid point and an element of volume, and assumed to emit isotropically.Since particle volumes are four times smaller than the X-IFU pixel area, and given the Athena telescope required angular resolution of 5 half-equivalent width, we consider them as point-like individual sources on the sky.The emission spectrum for each particle is assumed to follow an unabsorbed thermally broadened plasma emission, modelled through XSPEC using wabs*apec with a constant temperature of 7 keV and a metallicity of 0.7 Z .As in C18, the wabs absorption model (Morrison & McCammon 1983) is preferred for computational efficiency.A column density value of 0.03 × 10 22 cm −2 is used, and represents a typical high Galactic latitude value of the column density seen over the sky (Kalberla et al. 2005).
The turbulent motions of the gas are included by converting the line-of-sight component of the velocity field for each grid point into an additional redshift, which shifts the line by Doppler effect.No excess broadening is considered locally due for instance to microscopic turbulent motions.The corresponding total redshift z i of the i th cell is computed using the classical redshift composition: where z v,i is induced by the velocity of the cell, Finally, the normalisation N i of each emission spectra is provided for each cell through the apec normalisation where D A is the angular distance to the Coma cluster and V cell the volume of each cell (constant for our uniform grid).
The photon generation and the simulation follow the same process as in C18, with the same configuration of the X-IFU instrument through SIXTE (xifupipeline).Since the main objective here is the estimation of turbulent velocities and the velocity power spectrum, no background nor cross-talk are included in the simulation to avoid introducing additional instrumental systematics.However, as shown in C18, an accurate knowledge of the background components should not bias the following results (especially for redshift measurements).Further, since the observations are focused on the centre of a bright Coma-like cluster, the astrophysical and instrumental background levels are expected to be sub-dominant with respect to the cluster emission.A typical exposure time of 100 ks is considered for any of the following synthetic observations.

Post-processing of the data
Due to limited statistics of the observation, pixels are binned into regions with a signal-to-noise (S/N) ratio S/N = 200 (∼ 40 000 counts per region) to reduce the statistical uncertainty on the measurements.To do so, we Theoretical structure function (k inj =200 kpc -1 , k dis =10 kpc -1 ) Theoretical structure function (k inj =150 kpc -1 , k dis =20 kpc -1 ) Theoretical structure function (k inj =200 kpc -1 , k dis =20 kpc -1 ) Uncorrected data points selected a spatial binning using an adapted Voronoï tessellation of the plane 3 (Cappellari & Copin 2003), which provides ∼ 150 regions over the X-IFU FoV (∼ 20 pixels per region, Figure 2 -Bottom).This choice is motivated by the will to remain generic in our approach (this binning can be applied to any detector) and to provide round-shaped regions of constant S/N ratio, which ensure a faster convergence of the cosmic variance computation presented in CL19 than square regions.The spectrum from each region is fitted using the input XSPEC model with an additional broadening component (bapec) to account for the effect of the turbulent velocities.A simultaneous fit of all the free parameters (temperature, abundance, redshift, velocity broadening, and norm) is performed.Results from the fits are excellent.No bias is visible on both parameters, and the statistical error distribution of the measured velocity shift, δC stat , is consistent with a centred Gaussian, of standard deviation σ stat,C = 34 km/s (Figure 2 -Top right).This value of the statistical error was confirmed in every run, when using the same exposure and binning procedure (at constant β-model input).An example of binned input velocity map and the recovered output is provided Figure 2 (Bottom).

General approach
The previous E2E simulations are used as a test case to verify the formulas derived in Sect. 2. To do so, a computation of the structure function and its variance over a very large FoV would be required.This implies however to simulate large spatial grids with a refined mesh, which is rapidly computationally cumbersome (memory-and time-3 https://www-astro.physics.ox.ac.uk/~mxc/software/#binning wise).We take advantage of the ergodicity assumption of the turbulence to simulate many independent velocity fields over the previous 2048 × 256 × 256 mesh, and average over these iterations to derive an estimation of the structure function.In practice, for a given choice of the turbulent velocity power spectrum (i.e. a subset of α, k inj and k dis ), we create 100 different velocity fields, which are then observed using the previous E2E pipeline assuming a 100 ks exposure time.We thus obtain -for one choice of P 3D -100 independent synthetic velocity maps of the same random process.
In the case of the Coma cluster, we use as minimal dissipation scale of the turbulent velocity power spectrum 10 kpc, which represents a good compromise between current observational expectations (Gaspari & Churazov 2013) and future capabilities of the X-IFU (an X-IFU pixel size represents ∼ 2 kpc at the Coma cluster redshift).Given the binned regions of our map (∼ 5 pixels of diameter) any scale larger than 10 kpc should be resolved.We know that large injection scales, a few hundred kpc up to 1 Mpc, can occur due to merger event or subgroup accretion (e.g.Khatri & Gaspari 2016).As we aim here at verifying our analytical formulation of the statistical error, we consider injection scales ≤ 200 kpc, which provide a likely description of the injection scale and keep the velocity field size within range of our available computation power.The slope of the turbulent power spectrum is fixed at α = −11/3.

Structure function estimation
For a turbulent velocity power spectrum, SF can be estimated for each recovered velocity field through the previous formulas, using as separation s the distances between the centres of binned regions W. For a Voronoï tessellation, the centres are found by taking the weighted barycentre with respect to the number of counts in each region W. Given the homogeneous emission of the cluster toy model over the FoV, these points coincide in most of the cases with the geometrical centre of W (see the green crosses on Figure 2 -Bottom right for an example).
Though the emission profile of the cluster is the same, the non-constant region shape provided by the Voronoï tessellation creates slightly different spatial bins from one observation to the other.Hence, separations between regions do not follow a discrete mesh in each of the 100 observations of a given turbulent power spectrum.To compare the structure functions between the runs, we estimate SF on an a priori grid of spatial separations, equal for each iteration, and with a step size of ∼ 5 kpc, which is approximately the equivalent radius of a Voronoï region.For instance, the regions with distances between 10 and 15 kpc will be considered in the same bin to compute the value of SF in each run.The expected value of SF in the bin is then recovered by averaging over the 100 observations.The 'true' value of spatial separation in the bin is taken by averaging all the real distances contained in the bin.
Examples of the estimated structure function for a given P 3D are shown Figure 3, along with the corresponding 1σ deviation of each separation bin.The theoretical structure function associated with each 3D power spectrum is recovered through numerical integration of the analytical formula in a cylindrical frame (r, ϕ, x) along the line-of-sight x by (Zhuravleva et al. 2012): ) where J 0 is the Bessel function of the first kind and θ eff a fixed 'effective' radius to compute the 1D power spectrum P θ ρ .As the emissivity is not constant (see an example in Figure 1 -Right), the structure function also varies depending on the regions considered to compute it (except in annular regions as the emission satisfies a spherical symmetry here).However, since is slowly-varying over the detector FoV in this case, only minor changes of the SF are expected (≤ 5% over the FoV).A good approximation of the observed structure function can be obtained by evaluating the previous formula on the annular radius corresponding to half of the equivalent radius of the detector.For the X-IFU this corresponds to θ eff ∼ 34 kpc (1.24 ).
As expected, the uncorrected values of the structure function are positively biased due to the statistical uncertainties in the measurements .This bias can be corrected by subtracting the variance of the statistical error σ 2 stat,C (see Equation 16).Even with this correction, discrepancies of ∼ 5% remain, which can be related to binning effects (see CL19).After subsequent correction (Figure 3 -Right), the average value of the structure function matches with the analytical structure function (≤ 3% in average on the relative difference between the simulations and the computed structure function over all separations).Remaining sources of error may be related to the number of iteration used here to recover the expected value of SF (100).Numerical effects related to the integration, the choice of θ eff and the binning may also create deviations.

Structure function variance estimation
Similarly, the variance of the structure function over the runs can be compared to Equation 18.The geometrical terms, N nei (s) and N p (s), are derived from the binning map.D can be computed through D, unbiased for a centred Gaussian statistical noise.Its variance however, is biased such that (see Appendix A.5): Var[D(s)] is therefore estimated through Var[D(s)] over the 100 simulations, and corrected of its bias.The intrinsic cosmic variance is obtained using the formulas provided in CL19.To do so, a specific circular FoV of the instrument and a pixel size (or bin size) are needed.We consider here that bins are well described by disks of the same diameter as the Voronoï regions.Similarly the hexagonal detector is approximated by an equivalent disk.Three different options were considered for its radius R: -'Equivalent' radius R eq .If S A is the total area of the detector, R eq = S A /π = 67.6 kpc = 149.4 .-Radius of the inscribed circle, R in = 63.9 kpc = 141.2 .
Figure 4 (Left) shows the comparison for different values of R, and the case without statistical corrections.For large separations, the statistical terms presented in Equation ( 18) have little effects, but must be accounted for when considering smaller separations.Also, despite the slight differences between the considered radii, changes in the analytical values of the variance of a factor 2 are observed when statistics are included.Simulation points are comprised between the R = R ci and R = R in curves, and all three curves show a good agreement within error bars.Deviations between the simulated data and predicted errors is lower than 20% for all separations (i.e. less than 10% in standard deviation) with R eq providing the best results (10% in variance, hence ≤ 5% in standard deviation).
These curves accurately recover the shape of the expected variance, but show a consistent deviation at large separations, which can be related to two distinct effects.On the one hand, these separations are sampled only a handful of times within a single X-IFU pointing (i.e. one or two regions per iteration are separated by a detector diameter), thus creating a simulation-related sample variance in the data.On the other hand, the circular FoV approximation reaches a limit for separations of the same order (or higher) than R eq .Smaller deviations could also be caused by boxsize effects of the turbulent velocity grid.
Other factors can explain the remaining deviations between the theoretical curves and the simulations: -The computation of the analytical error formulas involve complex numerical integrals, which can account for errors in the variance around 1 to 2%.would however render the computational time unreasonable (∼ 1 week of computational time for 250 iterations for a given P 3D with our current set-up).-The approximation of the X-IFU FoV with a disk.A better description of the detector geometry in the computation of the intrinsic variance component could reduce this component, especially for separations larger than the detector equivalent diameter.-The bin shape used for the numerical computation, which assumes a uniform tessellation of the detector.
Voronoï binning was chosen to represent a more generic case than simple square regions.It does not, however, verify a uniform tessellation.Test runs performed using square groups of pixels show similar results in the comparison, suggesting little impact of this particular effect on our previous results.
This analysis was extended to other turbulent velocity power spectra, with k inj within 100 and 200 kpc, and k dis from 10 kpc to 30 kpc, showing similar results.

Validation of error formulas and relative contribution
Figure 4 (Right) shows the contribution of the sources of error for a given P 3D .Among the statistical error contributions, computed in each bin, the velocity field fluctuation term dominates the other two and shows a monotonous trend.The structure function fluctuation term is lower than the previous contributor, but becomes comparable at high s, where the number of regions used to compute SF is small.Finally, as expected, the purely statistical term is the smallest and can in most cases be neglected.Its contribution is mostly on small separations and minimal on average sepa-rations, where N p is highest.We note however that its value becomes comparable to the cosmic variance for small s.
For σ stat,C = 34 km/s, used here, the intrinsic cosmic variance term dominates over most of the separations (s ≥ 30/40 kpc).This result holds for all the P 3D tested.As suggested by Figure 4 (Left), the correction included by the statistical terms in the estimation of the total variance is mainly visible at small separations, validating the formulas out to ∼ 50 kpc.Testing the previous formulas for high s requires to decrease the contribution of the cosmic variance or to increase the other statistical contributors to enhance error terms for high s.All things being equal, these terms scale with σ 2 stat,C or σ 4 stat,C , and become dominant when σ stat,C ≥ 100 km/s.This can thus be achieved using larger pointings with a similar binning and exposure (to reduce the cosmic variance), or by artificially increasing the statistical error by for instance considering shallower exposures.The latter was considered, by reducing the previous runs to 20 ks exposures, and provides a good agreement with the theoretical formulas (see Figure 5).

Practical estimation of errors
As shown in CL19, an estimation of the cosmic variance can be obtained numerically with several approximations on the detector geometry.The statistical terms however are intrinsically related to the observational set-up.We provide here some solutions to estimate the quantities involved in Equations 16 and 18 for a specific instrumental configuration.
The values of N p (s) and N nei (s) are the simplest to derive, as they are related to the binning and geometry of the detector, and can be determined analytically with high accuracy.σ stat,C will be a direct output of the observation.However, an a priori estimation of the statistical error of an observation, which could be used to forecast structure function errors before an observation (e.g. to optimise the exposure strategy) is more challenging.Crude estimations can be derived using a simple Poissonian approach using the flux of the source and the exposure time.Provided E2E simulations of the instrument become sufficiently representative, a promising solution could be to derive σ stat,C numerically (a single pointing is required).
The variance of D is by far the most challenging term to estimate, as it requires multiple observations.In our simulations, the ergodicity assumption simplifies these computations, as accurate estimations of the previous terms can be obtained by averaging over a large number of iterations.This approach does not hold in flight, as multiple realisations of the same turbulent field are unlikely to be met (statistically speaking), even assuming a self-similar behaviour of the turbulence between clusters.Even if single pointings of a larger mapping of the same object are performed, an accurate computation of the variance cannot be achieved unless the emission profile does not vary across the FoV, which only concerns a handful of clusters (at a zeroth order approximation).A solution is once again to use numerical simulations.An idea is to find analytical formulas, similarly to CL19, and use them to obtain a fast numerical estimation of Var[D(s)] for a given power spectrum assumption.As this remains to be investigated, dedicated MC simulations remain at the moment the other possible solution.Although relying once again on an iterative method, these simulations do not require a full E2E approach, and can be performed without any prior on the instrument except its geometry (implying large gains in computational time).

Towards optimising observation strategies
Under the current assumption of 100 ks pointings, we expect statistical errors on X-IFU observations to be low and well constrained.Hence, over a single pointing, an accurate knowledge of the intrinsic cosmic variance will be sufficient to provide an error estimation of the SF at large separations (Figure 4), and thereby to investigate turbulence at large scales (i.e.injection scale).Statistical terms can in this case be neglected in the error computation.On the contrary, when measuring the structure function for small separations (i.e. to probe dissipation scales), a good estimation of all the statistical terms is paramount to ensure a proper error description.Depending on the science objectives of the X-IFU (large-or small-scale turbulence investigations), the observational strategy can then be optimised.
For a constant exposure time, several options can be explored.Deep pointings would provide a constant value of the cosmic variance, but significantly reduce σ stat,C , which may be interesting to explore small separations.On the contrary, since statistical errors are negligible at high s, multiple shallower pointings could be used to explore larger separations (notably larger than the detector) while reducing the cosmic variance.An optimal point where all errors are comparable across the separations could also be used.The formulas derived here and in CL19 demonstrate that accurate error estimations of the SF can be found by numerical integration.This approach is complementary to studies of turbulence limited through an iterative approach (e.g. 100 simulations of a large grids and associated errors), which can now be reduced to a handful of simulations.

Conclusion
With improvements in high-resolution spatially resolved Xray spectroscopy, measurements of line shifts, line broadening and structure functions will provide new insight on the turbulence at play within the ICM.In this paper, we addressed the challenge of computing these diagnostics and estimating their errors, related to both cosmic variance and measurement uncertainties.Specifically, this work extends the approach started in our companion paper I, which derives a formulation for the cosmic variance, and adds the contribution of finite statistics in the observations.All the formulas presented here should thus be coupled to those from Clerc et al. (2019, sub.) on the intrinsic cosmic variance.
We found that all the estimators, notably those of the structure function and its variance over various observations, are biased by the statistics of the measurements.For the variance of the structure function, these biases can be divided into multiple contributors, which add to the intrinsic cosmic variance: a purely statistical contributor, which scales with the squared variance of the statistical uncertainty (σ 4 stat,C ), and a contribution related to the intrinsic nature of the velocity field, which scales with σ 2 stat,C and depends on the spatial binning, the detector geometry, and the scales of the turbulence (injection and dissipation).We divided the latter into two terms in the case of Gaussian error field, a structure function fluctuation term related to the turbulent velocity power spectrum, and a velocity field fluctuation term, related to the scales of the turbulence within the FoV.The equations derived in Sect. 2 are generic and valid for any instrument with measurement uncertainties.
A specific application to the X-IFU instrument was performed to demonstrate the validity of these formulas in the framework of the future mission Athena.For this test, turbulent velocity fields of a toy model Coma-like cluster were generated for different underlying turbulent power spectra (here for different injection and dissipation scales of a Kolmogorov power spectrum), and used as inputs for synthetic observations with the X-IFU E2E simulator SIXTE.When such observations are averaged over a large number of different realisations of the velocity field (assuming ergodicity) and corrected of their biases, our results show excellent agreement with the analytical values, with relative errors below 3% over the spatial scales investigated.A comparison of the variance of the estimated structure function using the analytical error formulas presented here and in Clerc et al. (2019, sub.) also provides accurate result (better than 10% in variance, hence ≤ 5% in standard deviation).
Results presented here demonstrate that we can provide accurate estimations of the total variance of the structure function.For typical X-IFU observations of 100 ks, statistical terms in the structure function errors can be neglected for large spatial separations (s ≥ 70 kpc), but are required to investigate smaller separations (s ≤ 30 kpc).Depending on the science case, efforts can therefore be directed into reducing one or several of the error terms specifically, or into optimising the observational strategy (exposure and spatial map) and the spatial binning (region size and shape).A dedicated analysis of these optimisations, taking advantage of the fast computations enabled by the formalism proposed here, will be discussed in a forthcoming study.Measurements remain all the same a challenging objective, especially when systematics or the physics of the ICM (e.g.AGN, shocks) are considered.New results will thus also rely on alternative diagnostics, such as line non-Gaussianity, achievable through the spectral resolution of the X-IFU.
-Left).As in Z16, the normalisation of the 3D power spectrum is chosen to satisfyE[ S2 (0)] = (Mc sound ) 2, where E[ S2 (0)] is the expected velocity broadening of the line-of-sight component of the velocity field at the centre of the cluster, M the Mach number and c sound the sound celerity in the ICM.For our simulations, we use M = 0.3 and c sound = 1460 km/s (Z16).

Fig. 1 .
Fig. 1.Example of power spectra used in the simulations.(Left) Turbulent power spectra used in the simulation for different injections and dissipation scales.(Right) Normalised spectrum of X-ray emissivity of the Coma model at the centre used to compute Cn (||θ||2 = 0 kpc, red) and off-axis (||θ||2 = 250 kpc, blue).

Fig. 2 .
Fig. 2. Example of simulated velocity fields.(Top left) Example of an emission-measure weighted projection of the simulated line-of-sight component of the turbulent velocity field.In this case kinj = 1/150 kpc −1 and k dis = 1/20 kpc −1 .The shape and the extent of the X-IFU FoV is shown as a white dashed line (Top right) Absolute error distribution between the recovered line-of-sight velocity in one of the simulations and the corresponding input emission-measure-weighted velocity.The statistical error follows a centred Gaussian distribution where the Gaussian best fit (in red) is found for µstat = 0.2 km/s and σstat,C = 34 km/s.(Bottom left) Example of a synthetic X-IFU observation of bulk motion for kinj = 1/200 kpc −1 and k dis = 1/10 kpc −1 and (Bottom right) corresponding emission-measure weighted input map (binned).The small green crosses indicate the centres of the Voronoï regions.

Fig. 3 .
Fig.3.Estimated structure function (km 2 /s 2 ) averaged over 100 different observations of a velocity field generated with the same underlying turbulent velocity power spectrum as a function of the separation s (kpc).(Left) Raw structure function recovered for kinj = 1/200 kpc −1 and k dis = 1/20 kpc −1 .Different theoretical structure functions are also shown, the one associated to the run is given by red solid line.(Right) Same as left panel but the data points are corrected for the statistical bias and the binning projection effects (see CL19 for more information on the latter).Error bars indicate the ±1σ deviation within the 100 iterations.

Fig. 4 .
Fig. 4. (Left) Structure function variance (km 4 /s 4) averaged over 100 observations of a velocity field generated with the same underlying turbulent power spectrum (kinj = 1/200 kpc −1 , k dis = 1/20 kpc −1 ) and corresponding ±1σ error bars as a function of the separation s (kpc).The comparison to the theoretical models derived from the simulation and the formulas by CL19 is shown for a circular FoV of Rin (blue solid), Req (red dashed) or Rci (green dash-dotted).The sheer contribution of the cosmic variance without statistical terms for R = Req is given in light purple.(Right).Error contributions in the total variance for the previous case for R = Req when σstat,C = 34 km/s.The data points for the statistical and structure function fluctuation term are derived from analytical formulas, while the velocity field fluctuation term is computed using the 100 iterations (see text).

Fig. 5 .
Fig. 5. Same as Figure 4 for 20 ks observations using the same binning map.Only the case R = Req is shown.