© J. S. Sanders et al. 2020

1. Introduction

The velocity structure of the intracluster medium (ICM), the dominant baryonic component of galaxy clusters, remains poorly observationally constrained except in the core of the Perseus cluster (Hitomi Collaboration 2016, hereafter H16). Simulations predict that the ICM should contain turbulent or random motions, and bulk flows caused by the merger of other clusters and subcomponents (e.g. Lau et al. 2009; Vazza et al. 2011). In addition, there can be relative bulk motions of a few hundred km s⁻¹ caused by sloshing of the ICM in the potential well, generated by merging substructures (e.g. Ascasibar & Markevitch 2006; ZuHone et al. 2018). The central active galactic nucleus (AGN), via the action of its jets and inflation of relativistic bubbles, also likely generates motions of a few hundred km s⁻¹ (e.g. Brüggen et al. 2005; Heinz et al. 2010).

Measuring the velocities of the ICM is important for several reasons. Turbulent motions provide additional pressure support, particularly at large radii, which affects calculations of hydrostatic equilibrium and cluster mass estimates (e.g. Lau et al. 2009). AGN feedback injects a great deal of energy into the ICM which replaces that which is lost radiatively by X-ray emission (reviewed in Fabian 2012). Feedback has to provide distributed heating, but how this works depends on the microphysics of the ICM and the balance between turbulence and sound waves or shocks. Measuring the associated velocities would help constrain AGN feedback models. Velocity is an excellent probe of the microphysics of the ICM, such as viscosity, which acts to smooth velocity structure (e.g. Reynolds et al. 2005; ZuHone et al. 2018). Simulations predict a close connection between turbulent thermodynamic and velocity power spectra (Gaspari et al. 2014) which should be tested. In addition, measuring velocities will help explain the processes taking place in cluster mergers and the sloshing of gas in cold fronts, which is believed to last for several Gyr (Ascasibar & Markevitch 2006; Roediger et al. 2012).

Hitomi (Takahashi et al. 2016) with its high spectral resolution microcalorimeter Soft X-ray Spectrometer detector (Kelley et al. 2016) directly determined bulk and random motions in the Perseus cluster by measuring the position and width of the Fe-K emission lines (H16). It found a line-of-sight velocity dispersion of 164 ± 10 km s⁻¹ between radii of 30 and 60 kpc and a gradient of 150 km s⁻¹ bulk flow across 60 kpc of the cluster core. Away from the central nucleus and the northwestern cavity, the dispersion drops to ∼100 km s⁻¹ (Hitomi Collaboration 2018a, hereafter H18). Therefore, the Perseus core is not strongly turbulent, despite the obvious impact of the AGN and its jets on the surrounding ICM (e.g. Fabian et al. 2000). Unfortunately, Hitomi was lost and will not be available to observe more clusters, in particular merging galaxy clusters which would have probed different physical processes. The next missions equipped with a high-resolution instrument capable of measuring velocities are likely to be XRISM (the replacement for Hitomi, formerly known as XARM), to be launched in 2022, and Athena in the early 2030s.

Prior to Hitomi, Suzaku obtained velocities in several systems using the Fe-K line thanks to its accurate calibration, despite only having a charge-coupled device (CCD) detector. Tamura et al. (2014) examined multiple pointings, placing limits on relative velocities of 300 km s⁻¹ over scales of 400 kpc in Perseus, while on larger scales below 600 km s⁻¹, with hints of detection of velocity structure at 150−300 km s⁻¹ in the cold front 45 to 90 kpc west of the cluster. Ota & Yoshida (2016) examined several clusters with Suzaku, including the central region of the Perseus cluster. Systematic errors from the Suzaku calibration were around 300 km s⁻¹ and it had a large point-spread function (PSF) of around 1.8 arcmin (half power diameter). Chandra has also been used to measure gas velocities, finding evidence of bulk motions in merging massive clusters (Liu et al. 2016).

Another direct method for measuring velocities is to use the XMM-Newton Reflection Grating Spectrometer (RGS) to limit or measure line widths (e.g. Sanders et al. 2010; Pinto et al. 2015). Because of the slitless nature of the RGS, this is only effective in the central parts of peaked clusters. In several systems, limits below 300 km s⁻¹ were obtained using this technique (Sanders & Fabian 2013), which agrees with the low turbulent velocities later found by Hitomi in Perseus.

There are other indirect methods for measuring motions, including examining resonant scattering (e.g. Werner et al. 2009; Ahoranta et al. 2016; Ogorzalek et al. 2017). Measurements using resonant scattering are currently only possible for low-temperature X-ray emitting gas in elliptical galaxies or the cores of cool-core groups or clusters, except in the Perseus cluster (Hitomi Collaboration 2018b). The power spectrum of density fluctuations has been used to obtain upper limits on the amount of turbulence in some systems (Zhuravleva et al. 2014a). Repeating this analysis on simulated Perseus-like data (Walker et al. 2018a), however, shows that the signal from sloshing gas could contribute strongly outside 60 kpc radius. Turbulence can also be constrained based on arguments about its radial propagation (Bambic et al. 2018).

In this paper, we present a novel technique to calibrate the pn detector of the European Photon Imaging Camera (EPIC) onboard XMM-Newton, allowing us to measure bulk motions in clusters of galaxies. Section 2 describes our calibration methods. Using the calibration and the procedures described in Sect. 3, we analyse observations of the Perseus (Sect. 4) and Coma (Sect. 5) galaxy clusters to study bulk motions. Further systematic errors and future improvements are discussed in Sect. 6. Uncertainties quoted are at the 1σ level. Any astrophysical velocities quoted have a positive sign for material flowing away from the observer. All sky images are shown with north to the top and east to the left.

2. Data calibration

To measure velocities with EPIC-pn we need an excellent characterisation of the energy scale of the instrument. Because the cluster observations we wish to analyse are spread over the lifetime of the mission, we want to study this entire period. The currently-quoted accuracy for the detector energy-scale calibration for imaging modes is better than 12.5 eV (XMM-SOC-CAL-TN-0018). At the energy of Fe-K emission, this is equivalent to roughly 550 km s⁻¹. As we later show, although this quoted accuracy may be an average over the detector at the energy of the Mn-Kα calibration line, it appears to underestimate the energy shifts in some locations of the detector at higher energies.

In our analysis we used the XMM Science Analysis Software (SAS) version 16.1.0. When transforming from raw (PHA) to corrected (PI) event energies, SAS applies several different corrections to the energies of the observed events using its calibration database. This work is done by EPEVENTS, used as part of EPCHAIN, one of the EPIC-pn data-reduction chains. The effects corrected for include charge transfer inefficiency (CTI), operating along the columns of the detector and gain changes, including those obtained from long-term calibration source monitoring, and the effects of detector temperature variation. The standard calibration uses the measured position of the 5.89 keV Mn-Kα and 1.48 keV Al-Kα calibration lines with time. EPEVENTS applies a linear correction to events between these energies. Above the energy of Mn-Kα it does not introduce non-linearity of the energy scale. Owing to the nature of the development of the calibration, SAS applies several corrections in sequence when converting from PHA to PI values. We continue this trend by applying further corrections to the standard PI values. The current correction procedure is only applicable if the astrophysical line to be examined has an energy above the instrumental Mn-Kα, because of the energy scale interpolation used below this energy by SAS. Therefore if using Fe-K to measure redshifts, objects must lie below z ∼ 0.12.

Our first step is a correction for the average gain of the detector during the observation (Sect. 2.8), the second is a correction for the spatial gain variation across the detector with time (Sect. 2.9) and the final correction is to linearise the energy scale as a function of detector position and time (Sect. 2.10).

2.1. The EPIC-pn detector

The EPIC-pn camera (Strüder et al. 2001) contains a monolithic silicon CCD array, made up of 4 × 3 CCDs (Fig. 1). Each of the individual CCDs has 64 columns and 200 rows. The detector is read out along the columns as quadrants, each containing three CCDs. There are a variety of operating modes of the detector, but we are primarily interested in Full Frame (FF) and Extended Full Frame (EFF) modes, which read out the whole detector. In the EFF mode the fraction of out-of-time events (Sect. 2.6) is reduced by using longer frame times, at the expense of being more susceptible to pile-up (more than one photon arriving within a frame time). At an energy of 5.89 keV the detector achieved in its first nine months an average of ∼150 eV full width at half maximum (FWHM) spectral resolution.

Fig. 1. Geometry of the EPIC-pn detector. Shown are the numbers of the CCDs, the field of view exposed to the sky (FoV), the raw coordinates (RAWX and RAWY), the detector coordinates (DETX and DETY) and the approximate location of the “Cu hole”. The position angle (PA) is the angle between the detector and sky coordinates, an example of which is shown. The four quadrants are shaded in different colours.

2.2. Instrumental background emission lines

To calibrate the energy scale of the detector, we measured the position of fluorescent emission lines produced by elements present in the electronics boards behind the detector, and visible in all observations (Freyberg et al. 2001). The dominant emission line is the Cu-Kα line at 8.04 keV (Fig. 2). Weaker surrounding lines include Ni-Kα (7.47 keV), Zn-Kα (8.63 keV) and Cu-Kβ (8.90 keV). The spatial variation of the lines depends on the variation in the composition of the electronics boards. The distribution of Cu is relatively flat, except for a lack of emission in the centre of the detector (Fig. 3), which we call the “Cu hole”. The presence of this hole leads us to exclude the central part of the detector in our astrophysical analysis, as it cannot be calibrated through our procedure, although there is some residual signal from out-of-time events. Ni-Kα and Cu-Kβ/Zn-Kα emission shows spatial variation which appears to be anti-correlated.

Fig. 2. Total stacked spectra of the fluorescent background emission lines for 32.6 Ms of FF (red) and 19.7 Ms EFF (blue) data. Before stacking the individual spectra, the three orders of correction described in Sect. 2 were applied. Only events with PATTERN==0 were included. Au-Lα was not used in our analysis because of its weakness.

Fig. 3. Images of the detector in different fluorescent lines, created from stacked astrophysical observations. Shown are images in Cu-Kα (left; 7.805 to 8.285 keV), Ni-Kα (centre; 7.280 to 7.680 keV) and Cu-Kβ and Zn-Kα (right; 8.455 to 9.075 keV). Approximate background has been subtracted using scaled continuum images on either side of the lines. Some residual astrophysical emission can be seen in the centre of the image. The total exposure is around 64 Ms after cleaning for flares, combining both FF and EFF data and using PATTERN ≤ 4. The image was binned by 64 detector coordinate units and smoothed by a Gaussian of σ = 1 pixel width.

2.3. Datasets used for calibration

Our calibration analysis made use of two different sets of data. The first was a set of observations of astrophysical objects taken from the XMM archive (Table A.1). We obtained as many long public observations using EPIC-pn full-frame modes as possible, manually filtering them by examining the preprocessed images to avoid those with bright X-ray point sources and optical sources. The observations were chosen to cover the full mission lifetime (from 2000 to 2019), although proprietary periods prevented us from analysing some recent data. These astrophysical datasets contain the background Cu-Kα, Ni-Kα, Zn-Kα and Cu-Kβ fluorescent lines. The observations were split into FF and EFF sets for separate analysis. Those with a cleaned exposure of less than 10 ks were ignored. The cumulative exposures as a function of revolution for the two modes are shown in Fig. 4, where the XMM revolution is around 47.9 h, or two days.

Fig. 4. Cumulative FF (red) and EFF (blue) cleaned exposure of the astrophysical datasets used in the calibration as a function of time (in XMM revolutions).

For the energy scale correction (Sect. 2.10) we also made use of detector calibration observations (Table A.2). We ignored those datasets where the cleaned exposure was less than 10 ks (the filtering is described in Sect. 2.5). EPIC-pn includes a ⁵⁵Fe radioactive calibration source which can be viewed by the detector, illuminating much of the field of view with a spectrum of Al-Kα, Mn-Kα and Mn-Kβ emission lines, with energies of 1.48, 5.89 and 6.49 keV, respectively. The half life of the source is 2.7 yr, so it has become much fainter over the mission lifetime. The illumination of the detector by the source is also not flat, with some CCDs having poor coverage.

2.4. Spectral models

When fitting the spectral lines, we used the line shape measurements of Hölzer et al. (1997) for Mn-Kα, Ni-Kα, Cu-Kα, Mn-Kβ and Cu-Kβ. We implemented an XSPEC model for each line, summing the provided Lorentzian fit sub-components. The shift in the energy of a line is parametrised by a redshift, measuring the fractional energy shift (not a real physical redshift). For the weaker Zn-Kα line we made a model containing two Lorentzians with the energies taken from Bearden (1967), assuming a flux ratio of 0.52 for the ratios of the Kα₂ and Kα₁ lines, similar to the other Kα lines in Hölzer et al. (1997), and choosing line widths of 3.3 and 2.4 eV, respectively.

2.5. Selected events

In our spectral analysis, we only considered single-pixel events (PATTERN==0), which have the best energy resolution. When examining real cluster emission, we only examined events with FLAG==0, which discards bad pixels and neighbouring regions, and the outer part of the detector which is not illuminated by the sky. When examining the fluorescent lines and calibration source, we instead filtered with (FLAG & 0xfa097d)==0, which also includes events outside the sky field of view, enabling us to improve our understanding of the detector near its corners.

To reduce the effect of the flaring component of the XMM background, we applied a uniform maximum count rate threshold of 1.0 ct s⁻¹ in 100s bins, in the 10 to 15 keV band with FLAG==0 and PATTERN==0.

When extracting spectra, we binned them by the standard factor of five PI channels (i.e. 5 eV bins). Slightly improved uncertainties on measured line energies could have been obtained using a factor of one, although this would have meant we could not have used the standard response matrices.

2.6. Out-of-time events

If a photon hits the detector while it is being read out, it will be recorded with a wrong RAWY position, producing an out-of-time (OoT) event. For bright sources these events can be seen as a stripe along the RAWY direction. In FF mode, the fraction of events which are OoT is 6.3 per cent, while it is 2.3 per cent in EFF mode.

In our analysis we accounted for OoT events with an OoT event file, created by running EPCHAIN a second time with the parameter “withoutoftime = yes”. This option randomly shifts events along the RAWY axis, before the standard gain and CTI corrections are applied. We extracted OoT spectra using the same extraction regions as the sources. These spectra were used as background spectra when spectral fitting, scaling the exposure times up by the appropriate fractions shown above. If we added the source spectra together, these background spectra were similarly added. The effect of correcting for OoT events in our analysis was small; the galaxy cluster astrophysical results differed by less than their statistical uncertainties.

2.7. Response matrix

When fitting spectra we used a response matrix with fine energy bins (0.3 eV) created with RMFGEN. This energy bin size is sufficient to resolve velocities down to ∼15 km s⁻¹. With larger binning, line energies tend to be discretised onto the energy grid and error bars are often asymmetric. For the stacked and calibration observations, the response matrix was generated using an early observation, 0114120101, and for small RAWY coordinates, giving the highest achievable energy resolution. A variable Gaussian smoothing component was used when fitting, to account for the change in energy resolution over the detector position and over time. EPIC-pn response matrices only differ between spatial positions and observations in their spectral resolution. Fitting for the energy resolution significantly reduces the storage and computational time which would be required to make individual responses.

The ancillary response file used for fitting the calibration observations and fluorescent line data was generated from the same observation, disabling the effect of the mirrors and the detector filter (using parameters modelfiltertrans = false and modeleffarea = false), which do not affect the background lines.

2.8. First order gain correction (mean Cu-Kα)

The first order energy correction we made was to measure the average “redshift”, z, of the Cu-Kα line in each observation and correct it back to zero. The total Cu-Kα spectrum was extracted from each observation, including events outside the sky region, but excluding the central Cu hole. The hole was defined as a polygon with DETX and DETY coordinate pairs of 7601, 4390, −6209, 5255, −4531, 5978, −2361, 6383, −480, 6181, 995, 5573, 2094, 4994, 3133, 3905, 3136, −6116, 1747, −7447, 300, −8286, −1232, −8459, −3431, −8459, −4646, −8228, −5688, −7794, −6874, −6926, −7616, −6052.

We fitted the spectrum between 7.00 and 9.25 keV by a model with Cu-Kα, Cu-Kβ, Ni-Kα, Zn-Kα and powerlaw components. The redshift of each line component was allowed to vary independently. We applied smoothing by a common variable-width Gaussian component, to allow for changes in energy resolution. The spectra were fitted by minimising the C-statistic in XSPEC.

Figure 5 shows the Cu-Kα redshifts of the astrophysical observations as a function of time. It can be seen that there is a systematic offset between FF and EFF observations, which widens with time. At the 8 keV energy of Cu-Kα, the maximum shift is around 20 eV for FF and 35 eV for EFF modes.

Fig. 5. Best-fitting Cu-Kα redshifts for individual astrophysical observations as a function of time (in XMM revolutions). A 0.1 per cent correction is 8 eV at the energy of Cu-Kα.

We corrected observations for their average Cu-Kα gain shift by scaling the PI values of the observation by 1 + z. In detail, we take the integer PI value in an event file, add a uniformly-distributed random number between zero and one, scale by 1 + z, then round downwards back to an integer value.

2.9. Second order gain correction (stacked Cu-Kα)

After correcting for the average Cu-Kα redshift in each observation, we measured the Cu-Kα redshift as a function of detector position and time. We stacked the spectra from the astrophysical observations in time bins of 500 revolutions (overlapping, as the bins start in 250 revolution intervals), examining the FF and EFF modes separately. As the detector is read out in columns, we stacked the spectra for each CCD in bins of 1×20 raw detector pixels. We fitted each spectrum between 7.7 and 8.4 keV, by the Cu-Kα model with variable redshift (z) and smoothing, plus a powerlaw.

Figure 6 shows the redshift in six non-overlapping revolution intervals as a function of position, for the FF and EFF modes. Excluding the central Cu hole region, in FF mode the position of the Cu-Kα line is shifted between around −0.2 and +0.8 per cent (−15 to +65 eV), with shifts most positive furthest away from the readout, and most negative in the centre of the CCDs. In EFF mode there is a different distribution in shift over the detector, with a range of around −0.5 to +0.6 per cent (−40 to +50 eV), with the most positive shifts in the centre of each CCD, and the most negative shifts furthest from the readout. The difference between the FF and EFF maps is fairly constant in time.

Fig. 6. Redshift of the Cu-Kα line (in per cent) as a function of detector position over the mission lifetime in FF (left) and EFF (right) modes, after correcting for the average value. The range shown in the top-left of each panel is the XMM revolution. The shift has been measured in 1 × 20 pixel regions and a σ = 2 pixel Gaussian smoothing applied for display. Detector regions with a large uncertainty in a combined error map have been masked out.

The ranges of the shift appear to increase with time, but change fairly gradually, allowing us to correct for this effect. To correct the gain for an observation, we took the average Cu-Kα line shift, when measured from the uncorrected observation (as in Sect. 2.8) and the correction map with the central revolution closest in value with the observation revolution. Each PI value was corrected by the combined shifts, looking up the correction on the map according to CCD number and pixel coordinate.

We tested that the PI scaling procedure works to correct the stacked Cu-Kα line energy by correcting the input event files and repeating the procedure in this section. We found that the produced maps are centred around zero shift, with a standard deviation close to the average statistical uncertainty on each line measurement.

Another test is to look for position-dependent energy changes which are not corrected by the average Cu-Kα for an observation and the time-dependent average detector maps we have produced here. This may occur, for example, if there are random fluctuations in the gain of each CCD within single observations. To test for this, we extracted and fitted the spectra from individual CCDs for each observation which went into the stacking analysis, after applying our correction for the average Cu-Kα and the time-dependent detector maps in this section, excluding the central Cu-Kα hole. The distribution of residual Cu-Kα redshifts was well modelled by a Gaussian distribution, after excluding points with large error bars. We modelled the intrinsic distribution of the points before measurement uncertainties for each CCD as a Gaussian. Figure 7 shows the obtained Gaussian centre of the distribution and the intrinsic distribution width (σ) for each CCD in the two different detector modes.

Fig. 7. Intrinsic Gaussian centre and width (σ) of the distribution of Cu-Kα redshifts in each CCD for the observations which went into the stacking analysis, after correcting each observation for the average Cu-Kα value and the time-dependent maps (Fig. 6). This analysis excludes observations with statistical uncertainties on the Cu-Kα redshift of greater than 0.1 per cent. Results for FF and EFF modes are shown separately.

This analysis shows no evidence of strong variation within observations of the Cu-Kα line energy. We see a systematically low Cu-Kα energy after correction, but this is only of the order of 0.005 per cent. This may be because when we corrected for the average Cu-Kα energy for an observation, we fitted the whole fluorescent line complex energy range, but when we examined spatial bins, we only fitted for the Cu-Kα line. The shape of the underlying continuum may affect the Cu-Kα line energy slightly. The corrected Cu-Kα values have a distribution width of around 0.01 per cent, which is very good for our purposes (∼30 km s⁻¹).

We also tested a modified version of the 1st and 2nd order corrections, where the initial 1st order correction was done to each CCD, rather than to the detector as a whole. The changes to the galaxy cluster astrophysical results, using these corrections instead of those we present, were much smaller than the statistical uncertainties.

2.10. Energy scale correction

After applying the correction for the observation-specific Cu-Kα position and our correction maps (the gain corrections), it became apparent that the energies of other lines required further correction. The data preferred Zn-Kα and Cu-Kβ line redshifts greater than zero, and Ni-Kα line redshifts less than zero (Fig. 8). We examined calibration observations, also finding that the Mn-Kα and β lines were at too high energies.

Fig. 8. Demonstration of the necessity of an energy scale correction. Shown in black are the total spectra of calibration observations between revolutions of 1000 and 1500 (<7 keV) and the total spectra of calibration observations plus stacked observations in the same time period (>7 keV), outside the Cu hole region. The count axis for data >7 keV has an arbitrary scaling. Top panel: models and data, while the bottom panel shows their ratio. There are two “6-line” models. One (blue) assumes that the lines are at their standard energies. The other (pink) has a linear compression model for the line energies about Cu-Kα where δE/E = 0.0013 keV−1. The poorly-modelled line wings do not affect the best-fitting compression factors.

Fig. 9. Energy scale correction factor for each CCD as a function of time. Positive values mean that events have too low an energy at high energies and too high an energy at low energies. This plot shows the results of fitting the “6-line” model to stacked calibration observations and observational datasets after gain correction, excluding the Cu hole. The lines show linear plus constant fits to the results excluding the first revolution bin.

We found that a linear model, where observed line positions are compressed in energy by some factor proportional to their energy difference from the Cu-Kα line, modelled these shifts well. Unfortunately the correction factor has both spatial and time dependence. Figure 9 shows this energy scale correction factor as a function of time for each CCD. These results were obtained by fitting a spectral model made of Mn-Kα, Mn-Kβ, Ni-Kα, Cu-Kα, Zn-Kα and Cu-Kβ lines (the “6-line” model). The redshifts of the lines were parametrised by an energy correction factor multiplied by the difference in the expected energy of each line to the Cu-Kα line energy (−2.150, −1.555, −0.570, 0.590 and 0.865 keV, respectively). The Cu-Kα line position was free in these fits. We also included powerlaw components (each with a free slope) between energies of 5.2 and 7.1 keV, 7.1 and 8.0 keV, and 8.0 and 9.5 keV to account for continuum and broad line components. The line components were convolved by a Gaussian broadening component parametrised by its energy width. As the calibration source becomes weaker with time and has spatial variation, we fitted both stacked calibration observations (below 7 keV) and stacked calibration plus real observations (above 7 keV), in bins of revolution. We fitted the spectrum between 5.825 and 6.245 keV, 6.375 and 7.000 keV, 7.170 and 7.795 keV, and 7.905 and 9.200 keV (as in Fig. 8). The energy ranges chosen were designed to exclude the asymmetries seen in the Mn-Kα and Mn-Kβ line shapes by excluding the broader lower energy part of the line below the peak. We used the same response and ancillary response as previously (Sect. 2.7) and minimised the C-statistic in the spectral fits. We did not find any difference when comparing FF and EFF mode data, so we combined both in this analysis.

Figure 9 shows that there is an increasing energy scale compression (PI values are compressed around the Cu-Kα line) in each CCD with time. There appears to be an offset at zero revolutions, and in some CCDs a jump is seen at around revolutions of 250. The plot includes a constant plus linear with time model, fitted excluding datasets below a revolution of 250 and by minimising the χ². To avoid the fits being strongly biased by the points with small error bars (particularly those early in the mission when the calibration source was strongest), we increased the error bars on the data points by adding 10⁻⁴ keV⁻¹ in quadrature.

We do not know the physical origin of this non-linearity in the PI values from the detector. It is affected by the ageing of the detector unit in space, although it was apparent from launch. By examining spectra in different spatial bins, we found that differences were clearer in the RAWX axis of the CCDs (i.e. groups of columns were more different). There is also some variation in the RAWY axis, with the steepness of the slope in time being stronger towards large RAWY values.

To remove this effect from our astrophysical data and calibration observations, we fitted stacked spectra in bins of 250 revolutions, and 4 × 2 spatial bins in each CCD (using inclusive pixel ranges in RAWX of 0–16, 17–32, 33–48 and 49–64 and in RAWY 0–110 and 111–200). As above, the spectra were fitted (taking both stacked calibration and observational data) with the “6-line” model and the energy scale correction factor was obtained. We fitted the energy correction factors as a function of time for each spatial region with a linear model in time as in Fig. 9, excluding the first time bin.

To correct the PI values for an event file, we compute the correction factors as a function of position given its revolution number. For revolutions below 250, we just use the results for that time bin, but for the others use the linear relation fitted to all the time bins. We choose a correction factor given by which RAWX bin the event is in, but use linear interpolation of the correction factors in RAWY.

To test how well the energy scale correction worked, we corrected our input event files and stacked them in revolution bins, including calibration observations (for the Mn-Kα and β lines) and standard observations (for Ni-Kα, Zn-Kα and Cu-Kβ). Figure 10 shows the maps of the residual energy correction factor, after all our corrections. The corrections are well-scattered around zero.

Fig. 10. Maps of the residual energy correction factor in stacked observations, after applying our correction to this in the individual observations. The units are per cent per keV. The panels show revolution bins of 0–500, 500–1000, 1000–1500, 1500–2000, 2000–2500 and 2500–3000 (left to right and top to bottom). The uncertainty of these values is not the same across the detector or in time, so we use 8 and 16 pixel binning. The statistical uncertainty is much larger in the bottom-left area because of the non-uniformity of the calibration source. Pixels where the uncertainty is larger than 10−3 keV−1 are not shown.

There may be some variation in the correction factor between observations. We checked for this by fitting the spectra extracted from CCDs in the corrected calibration observations for the residual correction factor. Assuming the distribution can be represented as a Gaussian, we measured the width beyond what is expected from the statistical uncertainties. Figure 11 shows the centre and width of the distribution for each CCD. We obtain centres close to zero, with offsets less than about 10⁻⁴ keV⁻¹ and a width around 1.5 × 10⁻⁴ keV⁻¹. We excluded observations 0060740801 and 0108860301 in this analysis, as they showed very anomalous values of the correction factor. Including these observations in this analysis makes little difference to the results, but we suspect there is some problem with these individual observations, a small fraction of the sample.

Fig. 11. Distribution of residual correction factors (centre and width) obtained for each CCD, after correcting the calibration observations and refitting them. These results primarily use the Mn-Kα and β lines. Also shown are the distribution of correction factors obtained in each quadrant (three CCDs) for the individual astrophysical observations after correction, obtained from the Ni-Kα, Zn-Kα and Cu-Kβ lines.

We also examined the distribution of the residual correction factor for individual astrophysical observations. This is much harder to measure as the Ni-Kα, Zn-Kα and Cu-Kβ lines are weaker and closer to the Cu-Kα. We therefore measure the residual factor in CCD quadrants (Fig. 1) rather than in CCDs, and limit the analysis to those observations where it can be measured in each quadrant to better than 2 × 10⁻³ keV⁻¹. This signal is difficult to measure, as the individual measurements have an average statistical error of 1.1 × 10⁻³ keV⁻¹ and the statistical width is an order of magnitude smaller. Using the same Gaussian distribution assumption, we find consistent values of the width parameter. The centre parameter appears a bit larger, on average around 7 × 10⁻⁵ keV⁻¹, perhaps indicating that the linear correction model becomes inadequate at this high level of precision, or there is slightly inaccurate modelling of the details of the lines.

If we instead examine whole observations rather than separate quadrants, we decrease the measurement errors by a factor of two. In this case we obtain an average residual correction factor of (7 ± 1) × 10⁻⁵ keV⁻¹ and a distribution width of keV⁻¹.

In summary, these tests show that we have accurately calibrated the energy scale of the detector in different spatial regions. Summing the residual correction factor offset and distribution width for the calibration observations (10⁻⁴ keV⁻¹ and 1.5 × 10⁻⁴ keV⁻¹, respectively), our energy scale appears accurate to better than 2.5 × 10⁻⁴ keV⁻¹, or 150 km s⁻¹ (rounding up) at the energies of cluster Fe-K emission. When we simultaneously fit the spectra from observations at different epochs and with different detector positions, we may obtain better precision, as some of the calibration differences will cancel out.

The overall accuracy of the correction for a single observation is determined by the precision to which the Cu-Kα line energy can be measured. In archival astrophysical observations (Table A.1), for an exposure time of T_exp, the Cu-Kα redshift uncertainties have a functional form of 10⁻⁴ A(T_exp/10 ks)^−1/2, where A is between 1.6 and 2.5. This individual uncertainty is significantly smaller than the statistical uncertainty on astrophysical redshift for a single observation. This systematic error will be reduced if multiple observations are combined, as it is random and Gaussian. Therefore, this systematic error is unimportant, except perhaps for extremely bright regions, such as the core of the Perseus cluster, which we cannot examine here owing to a lack of offset observations.

We applied further checks on the calibration procedure using the Perseus and Coma datasets in Appendix A.

3. Data analysis

Described here are the data analysis procedures applied to the observations of the Perseus and Coma clusters (observation details are listed in Table 1). Firstly, we processed each observation data file with the standard EPCHAIN SAS tool. To filter bad time intervals from flares, we applied the same 1.0 ct s⁻¹ rate threshold as described in Sect. 2.5. The event files were energy corrected using the three steps described in Sect. 2, where we applied an overall gain correction using the observation Cu-Kα redshift, the spatial and time-dependent gain map correction and finally the energy scale correction. When analysing spectra, we used only single-pixel events (PATTERN==0) to be consistent with our calibration analysis. We also filtered using FLAG==0 to avoid regions close to CCD edges or bad pixels.

When spectral modelling within XSPEC, we used the APEC spectral code (Smith et al. 2001) to model cluster emission, absorbed with a PHABS (Balucinska-Church & McCammon 1992) component to account for Galactic photoelectric absorption. We used APEC version 3.0.9, which includes corrections derived from the Hitomi observations (Hitomi Collaboration 2018c) and assumed the Solar abundance ratios of Lodders & Palme (2009) used in the Hitomi analysis. In the fits the temperature, metallicity, redshift and normalisation were allowed to be free. The spectra were fitted between 4.00 and 9.25 keV, avoiding low energies where our energy corrections are most likely incorrect, but allowing us to derive a temperature using the continuum. We assumed a Galactic column density of 1.4 × 10²¹ cm⁻² for Perseus and 8.54 × 10¹⁹ cm⁻² for Coma (Kalberla et al. 2005), although the absorbing component has only a weak effect in the energy range examined. As background components we included Cu-Kα, Ni-Kα, Zn-Kα and Cu-Kβ lines, and a powerlaw component with its photon index fixed at 0.136 (an average taken from our archival observations). The redshifts of the background lines were tied together and allowed to vary. Although these background redshifts are not astrophysically interesting as they should be zero, they allowed us to test our gain correction. The normalisations of the lines and powerlaw were free to vary. For the background components we used an ancillary response file which did not include the effect of the filter or mirrors. We minimised the C-statistic when fitting the spectra.

We used the response matrices generated for each observation and spatial region. Spatial weighting maps with the number of counts between 3.00 and 7.00 keV were used when making the responses. We used the same 0.3 eV energy bin size in the response matrices as for our calibration analysis. No additional broadening was used in these fits, however. The ancillary response files were also calculated using the same spatial weighting maps.

There were multiple observations for each spatial region we examined. We tried two techniques for fitting the spectra. The first, with which we show the results for the individual spectral regions, was to do a joint simultaneous fit of the spectra. In this analysis we only included a spectrum for an observation if there were more than 2000 counts between 4.00 and 9.25 keV, to be able to sufficiently constrain the velocity. We allowed the normalisations of each component to be separate in each observation, but forced metallicities, temperatures and APEC redshifts to be the same. The background line energies were allowed to vary between observations. The second technique was to add the spectra from the different observations together and to create average weighted response and ancillary response files, using the number of counts between 4.00 and 9.25 keV as a weighting factor. The total spectrum was fitted, reducing the number of free parameters, improving the speed of the fitting procedure and making it easier for XSPEC to find the best-fitting parameters. We used this technique when fitting to produce the spectral maps of the clusters. As we will show, there is generally good agreement between the “joint” and “total” fitted results.

We also repeated the fits for the total spectra with the SPEX spectral code (Kaastra et al. 1996) instead of APEC, using version 3.04 with version 3.04.00 of its atomic code and tables. As the SPEX plasma model was designed to be used within SPEX only, we tabulated the lines and continuum produced by the model at 0.02 dex temperature intervals and exported it into XSPEC in APEC table format¹. This allowed us to use exactly the same modelling procedure for the two codes. As we show later, there is very good agreement between the APEC and SPEX results.

4. The Perseus cluster

The Perseus cluster, Abell 426, is the X-ray brightest galaxy cluster in the sky. There is clear evidence of feedback in the centre of the cluster, where the lobes of the central AGN, 3C 84 are displacing the ICM (Böhringer et al. 1993; Fabian et al. 2000). The AGN is also generating weak shocks in the ICM. Ripples, which may be sound waves generated by the AGN activity, are observed (Fabian et al. 2003). If sound waves, these ripples could combat a significant fraction of the radiative cooling taking place in a distributed fashion (Sanders & Fabian 2007). Alternatively, the surface brightness fluctuations in the cluster, if interpreted as turbulence, could energetically offset the radiative cooling taking place (Zhuravleva et al. 2014a), although sloshing of the gas in the potential well could be a significant contribution to the density variations with increasing radius (Walker et al. 2018a).

There is a well-known east-west asymmetry in the X-ray surface brightness (Churazov et al. 2003). The cluster is also extended in this direction on the largest scales. These edges in surface brightness and temperature are seen out to 700 kpc radius (Simionescu et al. 2012; Urban et al. 2014; Walker et al. 2018b). An explanation for these structures is that the gas is sloshing in the potential well of the system (Markevitch et al. 2001; Churazov et al. 2003; Walker et al. 2017), caused by a disturbance to the potential because of a minor-merger, as is also seen in simulations of clusters (e.g. Ascasibar & Markevitch 2006). In Perseus there is a chain of galaxies in the western side of the cluster which may be responsible for the sloshing (Churazov et al. 2003).

Perseus was the main target for observation by Hitomi (H16), providing an excellent opportunity to compare our technique with high resolution spectroscopy.

4.1. Fitting spectra from regions

In Perseus we made 21 non-overlapping regions by hand (Fig. 12 left panel) containing similar numbers of counts in the Fe-K complex (around 1000, but increasing to 4500 in the centre; Fig. 12 top-right panel), and following the surface brightness contours, such as those generated by cold fronts. We excluded bright point sources by hand, as automatic source detection was confused by the strong spatial variation of the source. Unfortunately, with the archival pointings the main AGN feedback region and main Hitomi pointing is within the Cu hole, so we cannot examine this area. However, the offset pointing of Hitomi mostly overlaps with our data, so we created an additional region for comparison with this (Region 22), which overlaps with our other regions 1 and 2. Region 22 is very similar to region 6 in H18, although the PSF of XMM is much smaller than Hitomi so our result will be less effected by flux from neighbouring regions.

Fig. 12. Regions examined in Perseus, chosen to cover areas with similar numbers of counts in the Fe-K complex. (Left panel) Regions plotted on the fractional difference in 0.5 to 2.0 keV surface brightness from the average at each radius. Region 22 overlaps with other regions as it covers the outer Hitomi pointing. (Top-right panel) Fe-K count map, showing the number of counts in each 3.2 arcsec pixel in the Fe-K complex (6.50 to 6.90 keV, restframe), after subtracting neighbouring scaled continuum images (6.06 to 6.43 and 6.94 to 7.16 keV, restframe). Gaussian smoothing of σ = 4 pixels was applied. (Bottom-right panel) Exposure map, plotting the total Fe-K exposure (log s) excluding the Cu hole and accounting for vignetting.

Figure 13 (left panel) shows the redshifts and velocities (relative to the H16 average velocity) we obtained for the regions in Perseus. Plotted are the results from joint simultaneous fits with the APEC model. In Fig. 13 (right panels) we also show the best-fitting temperatures and metallicities in each region. We list the results of the spectral fits in Table 2.

Fig. 13. Results for the Perseus regions. Shown are the best-fitting values and statistical uncertainties for the joint fit using the APEC model. (Left panel) Redshifts and velocities. The average shown is the weighted average for the non-overlapping regions 1 to 21. The right axis shows the redshifts as velocities relative to the H16 result (z = 0.01756). The NGC 1275 velocity of z = 0.017284 is from H18. Region 22 was chosen to overlap with the outer Hitomi pointing, for which we also show their result. (Right panels) Temperatures and metallicities. Also shown are the values for the outer Hitomi pointing from Hitomi Collaboration (2018d) and Simionescu et al. (2019).

For our region 22 which overlaps with the outer Hitomi pointing, we obtained a redshift of 0.0172 ± 0.0004 (all uncertainties quoted only include statistical errors), while a value of 0.0171 ± 0.0001 was found using Hitomi (H18), with a model for the PSF. The difference between these is within the statistical uncertainties, even if systematic effects are not included. We note, however, that there is a small portion of the Hitomi pointing region we cannot examine and that the PSF of Hitomi could blur out some signal over the cluster despite the applied correction. We can also compare our temperature value for that region ( keV) with the Hitomi PSF-corrected continuum (5.11 ± 0.05 keV) and line temperatures (5.00 ± 0.16 keV) from Hitomi Collaboration (2018d). Our temperature value is 6 to 8 per cent greater, although again our region does not overlap exactly with the outer Hitomi pointing, which also includes cooler gas at smaller radii. We obtain a metallicity (dominated by the Fe lines) of 0.61 ± 0.02 Z_⊙, while Simionescu et al. (2019) obtained a consistent Fe metallicity of 0.57 ± 0.03 Z_⊙ with Hitomi.

The average value of the redshift in the non-overlapping regions 1 to 21 is 0.01767 ± 0.00015, agreeing with the value obtained by H16 for the central region (0.01756). The fit of a constant value for the regions is poor (reduced χ² = 37.3/20 = 1.86). If we assume the intrinsic bulk velocities come from a Gaussian distribution and obtain the model parameters, we find a central redshift of 0.01778 ± 0.00025 and a distribution width of 0.00071 ± 0.00030 or 214 ± 85 km s⁻¹.

The offset between the average and Gaussian centre is probably caused by the higher velocities in regions 14–18, which are 200–300 kpc from the nucleus to the east, appearing to have systematically larger velocities than those on the other side (regions 11–13). The eastern regions are on average 480 ± 210 km s⁻¹ higher than the H16 average result, while those to the west are consistent at 20 ± 190 km s⁻¹ relative to Hitomi. Regions 14 to 18 have systematically higher surface brightnesses at that radius (Fig. 12 left panel). The cluster is also cooler towards the eastern side (Fig. 13 top-right panel). To see the significance of the signal in this region, we show an example spectrum for region 16 (Fig. 14). We show the residuals for a model with the best-fitting redshift and those using the reference cluster Hitomi value. If the cluster redshift is used, there are clear positive and negative residuals either side of the He-like Fe XXV line complex, indicating a line shift. In Fig. 15 we show the change in fit statistic if the redshift is varied from the best-fitting value and the spectra refitted, for selected spatial regions. The curves are smooth and only have a single minimum, and only become asymmetric far from the best-fitting redshifts, showing that the best-fitting values and statistical uncertainties are robust.

Fig. 14. Example fit of the total spectrum for region 16 in Perseus, rebinned for display to have a signal to noise ratio of 10 in each spectral bin. The vertical line is at an energy of 6.55 keV. (Top panel) Fitted data and best-fitting cluster model (fitted redshift), showing its cluster and background components. Also shown is the cluster component with the redshift changed to the Hitomi cluster value. (Centre panel) Residuals of the total best-fitting model. (Bottom panel) Residuals of the total model with the APEC redshift changed to the Hitomi value.

Fig. 15. Change in fit statistic as a function of redshift for selected regions in the Perseus cluster, for the total spectral fits.

If we model the distribution of velocities in only the central region (regions 1 to 10) as a Gaussian we find a central redshift of 0.01739 ± 0.00021, while the distribution width is small (28 to 173 km s⁻¹; at the 1σ level).

There are further possible systematic errors, which we discuss in Sect. 6.1, although they appear unimportant compared to the energy calibration and statistical uncertainties.

4.2. Spectral maps

We can see the velocity, temperature and metallicity structure better by making a smooth map of the cluster. To do this we extracted and fitted spectra from regions over the cluster, moving in grid with a spacing of 0.5 arcmin. To make the regions better follow the steeply-peaked surface brightness profile of the cluster, we used elliptical regions with a 3:1 axis ratio, rotating them so that the longest axis lay tangentially to the vector to the central AGN. We adaptively changed the radii of the ellipses to have ∼700 counts in the Fe-K complex after continuum subtraction, to have similar uncertainties on the velocities. Point sources were excluded from the regions. The map was created by fitting the total spectra rather than by a joint fit. We used the same response matrix for all the pixels, but created weighted-average ancillary responses for each ellipse, created from the individual ones for each observation.

Figure 16 shows the resulting maps of the velocity, temperature and metallicity. For comparison we also show an X-ray image of the cluster from XMM, and the fractional difference of the X-ray image to the average at each radius from the nucleus. The temperature plot shows that there is a cool spiral of gas going from the west of the cluster to the east via the north. Its outer edge varies from ∼170 to 260 kpc radius. The outer edge of the cool spiral can be seen in the surface brightness image, although there is another region of higher surface brightness which lies at larger (∼330 kpc) radius to the east. There is a spiral metallicity structure associated with the cooler spiral. To the east of the core, the low temperature spiral gas has low metallicity, with the hotter material at smaller and larger radii having higher metallicities.

Fig. 16. Perseus cluster maps. Shown is the 0.5 to 2 keV X-ray surface brightness and contours (log cts s−1 per 1.6 arcsec pixel; top left), fractional difference of surface brightness to radial average (top right), temperature (keV; centre left), metallicity (Z⊙; centre right), velocity relative to the cluster value in H16 (km s−1; bottom left) and 1σ statistical uncertainty on velocity (km s−1; bottom right). Maps were created by a moving 3:1 elliptical region (rotated to lie tangentially to the nucleus) containing ∼700 Fe-K counts (see example ellipses in velocity error plot). The pixels are not statistically independent.

Examining the velocity structure, there appears to be a clear correlation between the temperature and velocity maps. These correlations can also be seen in radial profiles of the maps (Fig. 17). The cool spiral has lower velocities, while the material lying outside it to the east and south has larger velocities. There is possibly higher velocity material to the south-west, although that location is close to the edge of the field of view and the elliptical regions are rather large, smoothing out small scale structure. We do not see clear velocity structure within the cool spiral. These maps suggest that some of the regions in Fig. 12 could contain velocity structure within them which is being averaged out. Better regions could be defined from the velocity map itself, but this is likely to bias the resulting velocities towards the extremes.

Fig. 17. Radial profiles of temperature, metallicity and velocity pixels from the Perseus maps in Fig. 16 coloured by position angle. Error bars are only shown for every 10 points. Radii are from the position of 3C 84. Nearby points are not statistically independent.

4.3. Perseus discussion

Other than Hitomi, previous measurements of the velocity structure in Perseus include those of Dupke & Bregman (2001), who claimed a 1000 km s⁻¹ velocity difference on 30 arcmin scales on an axis inclined with a position angle of 135°, 45° to the major east-west axis of the cluster. This is on larger scales than we probe here, so we cannot make comparisons.

Tamura et al. (2014) examined multiple pointing positions in the Perseus cluster and in the core using Suzaku. On larger scales they limited motions to <600 km s⁻¹. Our results are consistent with this on limit. They found a hint of bulk motions around −150 to −300 km s⁻¹ around 45 to 90 kpc to the west of the cluster centre. This is a high surface brightness region inside the inner western cold front and includes our regions 1, 2 and 22. Their result is roughly consistent with our results, with regions 22 and 1 being lower than the Hitomi average by ∼220 km s⁻¹, although the lower surface brightness region 2 is greater than Hitomi by ∼50 km s⁻¹. More detailed comparisons are difficult as their results are measured relative to each other and their reference region is inside the Cu-hole in our data. Ota & Yoshida (2016) placed limits on the velocities in Perseus of 860 km s⁻¹ at the 90 per cent level, which is consistent with our results.

We found evidence of a bulk motions of ∼480 ± 210 km s⁻¹ away from the observer, at a radius 250 kpc to the east of the Perseus cluster core. At this radius the gas velocities are unlikely to be affected by the central AGN. High spatial resolution maps (Fig. 16) suggest that there is a good correlation between the structure of the cold front to the east, as seen in temperature map, and the gas velocity. The material at smaller radius with lower temperature appears to have a low negative velocity, while the hotter gas outside this radius has a high positive velocity. This cooler material has a spiral-like twisted structure. Such structures are seen in simulations of ICM sloshing in cluster potential wells. Figure 18 compares a simulation of a Perseus-like cluster with bulk-sloshing motions and turbulence initiated by an off-centre encounter with a subcluster, with our Perseus maps. The surface brightness residuals, temperature and velocity have similar morphologies to what we see in Perseus, and are indicative of sloshing, although the simulations are not an exact match to what is observed.

Fig. 18. Comparison of a sloshing simulation (top row) and data for the Perseus cluster (Fig. 16; bottom row). The columns show surface brightness residuals (left), velocity (centre) and temperature (right). The simulation is of a Perseus-like cluster which has bulk sloshing motions and turbulence initiated by an off-centre encounter with a subcluster of 1/5 its mass, from ZuHone et al. (2010, 2018).

Unusually, the higher surface brightness, cooler material has a lower metallicity than the hotter material at larger radius and at smaller radius, which is somewhat contrary to simulations of sloshing of gas in potential wells (e.g. Roediger et al. 2012), where high metallicities are correlated with high surface brightnesses, and as seen in other clusters (e.g. Centaurus; Sanders et al. 2016). There is also an unusual surface brightness structure at the temperature interface. Moving across the front, the X-ray surface brightness is higher than the radial average, but drops temporarily at the interface between the hotter and colder gas. Usually across sloshing structures there are an alternate series of cooler, metal rich relatively-brighter, and hotter, metal poor and relatively-fainter regions. The lack of a typical alternating structure in Perseus may indicate that the structure is not simple sloshing. It is possible that there are additional structures induced by a merging subcluster, as is indicated by the east-west chain of galaxies (Churazov et al. 2003).

ZuHone et al. (2018) simulated the emission-weighted velocity field in Perseus from sloshing motions, projecting it along different lines of sight. The range of line shifts found in the simulations was between ∼ − 250 and 300 km s⁻¹. We find velocities which are a bit larger than this (500 ± 200 km s⁻¹), although the measurement uncertainties and ∼150 km s⁻¹ systematic uncertainties could bring these into agreement. Alternatively, if there is an undergoing merger, then this could also be responsible for the increased velocities, although the excess chain of galaxies lies to the western and not the eastern side.

If the gas were sloshing in the potential well, it might be expected to see a corresponding negative velocity in the cold front on the other side of the core at smaller radius (ZuHone et al. 2018). The average velocity for regions 1 to 3 on the western side of the core where the X-ray surface brightness is high, is −150 ± 90 km s⁻¹. This magnitude of this velocity is much lower than our detection to the east.

Assuming a Gaussian distribution, the width of the velocities for all the non-overlapping regions is km s⁻¹. If we examine the distribution of velocities in the inner regions (1 to 10), the width is lower ( km s⁻¹) and close to the 100 km s⁻¹ Doppler width of the Fe-K line as observed by Hitomi regions where there is no AGN disturbance (H18), although the spread of bulk velocities and the line width are not necessarily measuring the same set of motions.

Examining the velocity map structure inside the cold front to the west shows little evidence of structure. This western cold front appears much more typical, with cooler, more metal rich and higher surface brightness material inside the front and the opposite outside. There is little evidence of a change in velocity across this edge to the west, but there is a possible high velocity region to the south outside the front, seen in region 14.

5. The Coma cluster

The Coma cluster, Abell 1656, contains two dominant galaxies, which are NGC 4874 and NGC 4889. The two galaxies have a large velocity difference of ∼700 km s⁻¹ between them (Fitchett & Webster 1987), likely caused by a major merger in the cluster. The flat central X-ray surface brightness supports the idea that the cluster is unrelaxed (Briel et al. 1992; White et al. 1993). The cluster hosts a giant radio halo (Willson 1970), indicating it is a merging system. In addition, an analysis of the velocity distribution of the galaxies in Coma shows that there are subgroups associated with the two central galaxies (Adami et al. 2005) and up to 15 further subgroups in the system. A weak lensing analysis (Okabe et al. 2010) also suggests multiple substructures within the cluster, two of which are associated with the central galaxies. It is unclear which of the two central galaxies is the newest arrival in the cluster, although the velocity structure of the Coma galaxies appears to better match NGC 4874 (Adami et al. 2005).

Previous XMM observations have shown that the group NGC 4839 to south of the cluster (Neumann et al. 2001) appears to be merging, as indicated by a temperature enhancement between the group and the cluster. There is also a radio relic in the direction of the merging group (Jaffe & Rudnick 1979; Giovannini et al. 1985).

In addition, there is another group to the east seen by its X-ray emission (Vikhlinin et al. 1997), associated with NGC 4911 and NGC 4921. This system has an estimated mass of 5−8 × 10¹² M_⊙ (Neumann et al. 2003). The outer edge of its X-ray emission lies at the edge of the radio halo (Simionescu et al. 2013).

There are also a number of cooler denser arms of X-ray emitting gas that extend into the core of the cluster, which may be stripped material from merging subclusters (Sanders et al. 2013), some of which appear connected to the NGC 4911/4921 structure.

5.1. Results

Owing to the lower total number of counts in the Fe-K complex in Coma compared to Perseus, we split the cluster into 11 regions, with ten in the centre and one in the merging NGC 4839 group (Fig. 19 left panel). The group region also contains gas ahead of the group in the cluster, as there were not enough counts to make a pure extraction region, and it is therefore not a clean measurement of its properties. In this system we are able to get complete coverage of the central region (Fig. 19 right panels), because of the large number of separate pointings with different positions. We arrange our central regions in a 3 × 3 grid, with the central box split into two to separate the two central galaxies.

Fig. 19. Regions examined in Coma. (Left panel) The regions plotted on an exposure-corrected surface brightness image in the 0.5 to 2.0 keV band. (Top-right panel) Image in the Fe-K band, after subtracting continuum, in 10−3 counts per 1.6 arcsec pixel. (Bottom-right panel) Fe-K exposure map in log s, including the effect of vignetting and the Cu hole.

Figure 20 (left panel) shows our obtained velocities for the different regions, compared to the velocities of the central galaxies and merging group. Plotted are the velocities of the two central galaxies, NGC 4889 and NGC 4874, and modelled values for the NGC 4839 and NGC 4911 groups using a substructure analysis (Adami et al. 2005). The temperatures and metallicities we obtain for the regions are shown in Fig. 20 (right panels). We list the results of the spectral fits in Table 3.

Fig. 20. Results for the Coma regions. Shown are the best-fitting values for the joint fits using the APEC model. (Left panel) Redshifts and velocities. The right axis shows the redshifts as velocities relative to the redshift of NGC 4889. (Right panels) Temperatures and metallicities.

The results show that we obtain gas velocities which are very similar to the optical velocities of the two central galaxies, appearing to cluster around these values. The regions in the centre, to the west, to the south-west and to the south have velocities consistent with NGC 4889, while those to the north-west, north, north-east, east and south-east appear to better match NGC 4874. The reduced-χ² for a single velocity value for regions 1 to 10 is poor (1.7 = 15/9). If we model the intrinsic velocity distribution of the jointly-fitted results for these regions by a Gaussian, we obtain a central redshift of 0.02250 ± 0.00045 and a width of (100−440 km s⁻¹), although a Gaussian distribution appears to be a poor model given the data.

We can see the structure in more detail by mapping the parameters on smaller spatial scales. We moved a 4 arcmin radius circular region across the cluster in increments of 1 arcmin. The spectra were added in each region to make total spectra and fitted using the same model as previously. Similarly to Perseus, we used an average response matrix for all regions but created separate weighted ancillary responses for each region. Figure 21 shows the best-fitting temperature, metallicity and velocity for the nearest overlapping position, when the data are sufficient to constrain the parameters to the level given. For the velocity parameter we also show the statistical uncertainty for each position. Statistical uncertainties on the other quantities, such as temperature, increase in the outer parts, because of the declining surface brightness. We also show images of the cluster to the same scale, and a map generated from the fractional difference of the average surface brightness at each radius. Values in the maps are not statistically independent, unless they are separated by more than 8 arcmin.

Fig. 21. Coma cluster maps. Shown is the 0.5 to 2.0 keV X-ray surface brightness with contours (log cts s−1 per 1.6 arcsec pixel; top left), fractional difference of surface brightness to radial average (top right), temperature (keV; centre left), metallicity (Z⊙; centre right), velocity relative to NGC 4889 (km s−1; bottom left) and 1σ statistical uncertainty on velocity (km s−1; bottom right). Spectral results were obtained from spectra extracted from a sliding 4 arcmin radius extraction region (circle) over a grid. Regions where the uncertainty on the temperature is >1 keV, metallicity >0.07 Z⊙ or velocity >600 km s−1 are not shown. The position of the four central lensing-map peaks of Okabe et al. (2010) are indicated by squares.

There is interesting structure in the maps, with the temperature, metallicity and velocity providing complementary information. The gradient in temperature across the core of the cluster has been seen using ASCA (Donnelly et al. 1999), XMM (Arnaud et al. 2001) and Chandra (Sanders et al. 2013). The boundary between the hot and cold emission has a similar shape to the surface brightness ratio image. The hotter gas appears to lie outside of the northern surface brightness edge. Using Chandra we found cool arms of dense gas which are hinted at with this analysis.

There is a region extending north-south of high metallicity gas which lies between the two central galaxies. As the individual regions showed, the gas to north and east of the cluster has a higher redshift than that to the south and west. The north-south higher metallicity gas appears to have lower velocity than the surrounding material. There is also a weaker extension of high metallicity gas from NGC 4889 to NGC 4874 and beyond, which also has lower velocity.

To the south-east and south-west of the two central galaxies are regions with around 400 km s⁻¹ lower velocities, while to the north the velocity jumps up by around 700 km s⁻¹, in line with the measurements for the individual regions.

5.2. Power spectrum of density fluctuations

Recent theoretical progress has shown that the velocity field in galaxy clusters can be constrained using the power spectrum of density fluctuations (Schuecker et al. 2004; Churazov et al. 2012; Gaspari & Churazov 2013). In the presence of random, isotropic gas motions, the amplitude of density fluctuations is linearly related to the velocity dispersion of turbulent motions (Gaspari et al. 2014; Zhuravleva et al. 2014b), provided other contributions to surface brightness fluctuations are not important (Walker et al. 2018a). The velocity dispersion indirectly recovered through this method can be compared with the direct measurements we present.

We extracted count images in the 0.5–2.0 keV band from all available XMM pointings. We used the XMMSAS task EEXPMAP to compute an effective exposure map of each observation, taking vignetting and quantum efficiency into account. To model the particle background, we extracted count images in the same energy band as for our observations from a large set of filter-wheel-closed (FWC) data. We rescaled the FWC images to match the count rates measured in the unexposed corners of the three EPIC detectors. We then combined the data from all observations and from the three EPIC cameras to create a mosaic image of the cluster, and we applied the same procedure to the exposure and background maps. Point sources were detected using the XMMSAS task EWAVELET and excised from the image. We also masked the region surrounding the infalling galaxy group around NGC 4839 to avoid inducing additional fluctuations in the main cluster.

To determine the amplitude of density fluctuations in Coma, we followed the procedure outlined in Eckert et al. (2017). Namely, we determined the ellipticity and the ellipse angle of the observed emission using a principal component analysis, and we fitted the large-scale gas distribution with an elliptical beta model, which was found to describe the gas distribution accurately. We then constructed a model image by convolving the best-fitting model with the exposure map and adding the model particle background. Finally, we used the modified Δ-variance method of Arévalo et al. (2012) to reconstruct the power spectrum from the residual image in a region of 1.5 Mpc radius. For more details on the reconstruction of the power spectrum we refer to Arévalo et al. (2012).

In Fig. 22 we show the amplitude of surface-brightness fluctuations as a function of the wave number k = 1/ℓ. On the same figure we show the three-dimensional amplitude of density fluctuations A_3D = δρ/ρ, which is recovered after correcting for projection effects using the window function W(k) (see Eq. 11 of Churazov et al. 2012). We calculated the window function numerically by simulating fluctuations on top of the best-fitting three-dimensional model and computing the ratio between 3D and projected power spectra.

Fig. 22. Amplitude of gas density fluctuations in the Coma cluster as a function of the wave number k = 1/ℓ. The curves show the amplitude of projected surface brightness fluctuations (red) and the amplitude of 3D density fluctuations (blue). The uncertainties are indicated by the shaded areas.

The Coma power spectrum recovered here agrees well with the measurements of Churazov et al. (2012) over the common range of scales. The fractional amplitude of density fluctuations peaks at 1/k ≈ 1 Mpc. This agrees with the results of Khatri & Gaspari (2016), albeit with an overall lower normalisation. The difference with the latter analysis probably lies in the treatment of the NGC 4839 group, which the authors did not mask out from their residual image. We measure a maximum 3D amplitude A_max(δρ/ρ) = 0.092 ± 0.008. In case the detected fluctuations can be entirely ascribed to random gas motions, the maximum amplitude is linearly related to the 1D Mach number ℳ_1D = σ_v/c_s as ℳ_1D ≈ 2.3A_max(δρ/ρ) = 0.22 ± 0.02 (Gaspari et al. 2014). Given the average sound speed in the medium c_s = (γkT/μm_p)^1/2 ≈ 1480 km s⁻¹, surface-brightness fluctuations in the Coma cluster imply a velocity dispersion σ_1D = 310 ± 25 km s⁻¹, which matches well with the dispersion of the line-of-sight bulk velocities we present.

5.3. Coma discussion

The Coma cluster is a complex merger environment. In addition to the two-subclusters which presumably hosted the two central galaxies, NGC 4889 and NGC 4874, there is the well-known NGC 4839 group to the south and a group containing NGC 4911/4921 to the east.

Previous studies have suggested that the NGC 4839 group is currently merging with the cluster and has not yet passed through the cluster centre (Neumann et al. 2001). Evidence of this includes the high temperature region between the group and cluster, and its morphology, with the most of the group galaxies on the cluster-side of the group, and a tail of X-ray emission on the other side. However, it has been recently argued that a post-merger scenario better matches the observed features when compared to simulations (Lyskova et al. 2019). We find that the gas velocity of the group (including some material between the group and cluster) is consistent with the optical velocity of the group. The gas in the group should be slowed down as it interacts with the cluster by around 540 km s⁻¹ (Neumann et al. 2001), while the galaxies are not affected. We do not find evidence of this, although our measurement uncertainties are large, and the spectrum also includes some cluster emission.

Examining the individual regions in the Coma cluster, we find the X-ray gas appears to trace the optical velocities of the two central galaxies. The central cluster regions suggest velocities which are lower than that of the NGC 4911/4921 group. A map of the gas velocity (Fig. 21 bottom-left panel), shows that there is lower velocity material extending from NGC 4889 to the south-east and to the south-west of the cluster core, while the material to the north and east of the core better matches the NGC 4874 velocity. This material with velocities similar to NGC 4889 also has higher metallicity and lower temperature. An alternative origin for the high velocity gas is that it is material stripped from NGC 4911/4921 group, although the cool arms seen in the cluster core appear to be lower velocity material.

The higher metallicity north-south filament might be stripped material or a residual core from a merging subcluster, possibly the one which contained NGC 4889, considering the similarity in velocity. In addition, there are the arms of cool material extending towards the east which are likely to be stripped from the NGC 4911 subcluster (Neumann et al. 2003; Sanders et al. 2013). These arms make a U-shaped structure, however, with NGC 4889 at the tip of one of the arms and the group at the bend.

Adami et al. (2005) found that the galaxies within their cluster-core subcomponent have a velocity distribution peaking between the velocities of the two central galaxies. In contrast, Gerhard et al. (2007) found that the main component of the intracluster planetary nebulae within the cluster core have a peak velocity of ∼6500 km s⁻¹ (z ∼ 0.022), close to that of NGC 4889, rather than NGC 4874. Our velocity for the gas in the central region (regions 5 and 6) is similar to NGC 4889. However, if we take the whole examined cluster centre (regions 1–10), the modelled Gaussian distribution has a centre of 0.0225, or a weighted average of 0.0227; both of these values are close to the average velocity of the two central galaxies.

Biviano et al. (1996) also examined the velocity distribution of the galaxies in the cluster, finding a velocity gradient in the faint galaxies, decreasing from the north-east to south-west. This is similar to the morphology of the gas velocity seen in our maps, although the galaxy velocity gradient covers a wider spatial area. This velocity distribution is reminiscent of rotation, although the connection between structures in the velocity map and the other X-ray properties argues against this simple interpretation.

Okabe et al. (2010) found several weak-lensing subclumps within Coma, the inner four of which are indicated in Fig. 21. Two of these clumps are clearly associated with the central galaxies. Andrade-Santos et al. (2013) modelled the X-ray emission using the inner three mass peaks under the assumptions of slow motions and no gas, and found a good agreement with the ICM distribution. This model predicts that the south-east subclump lies around 400 kpc from the cluster mid-plane, in contrast with the two central galaxies, which would lie at ∼120 and ∼10 kpc, for NGC 4889 and 4874, respectively. Examining our maps, it is possible that the subclumps may be associated with some of the features we see. For example, the south-east subclump is coincident with a region of lower velocity, lower temperature and high metallicity. It is, however, unclear without further work whether subclumps could be entirely responsible for the structures.

If we take the line-of-sight velocity width of the bulk motions (V_los), we can estimate the ratio of energy density in bulk motions to the thermal energy density, using (Werner et al. 2009), where μ is the mean particle mass, m_p is the proton mass and kT is the temperature (∼8.5 keV). This measure excludes any turbulent component along the line of sight. Taking our Gaussian width of 100−440 km s⁻¹, implies that f_bulk < 0.15, at the 1σ confidence level. This estimate, however, assumes a Gaussian distribution of the gas velocity, which may not be the case. In addition, simulations show that the assumption of a symmetric velocity distribution lead to overestimates of the pressure support from motions (Vazza et al. 2018).

6. Discussion

6.1. Further possible systematic uncertainties

There are possible systematic uncertainties which may be present in addition to the energy calibration uncertainty. As stated in Sect. 3, we fitted the spectra jointly (producing the results above) or fitted the total spectra. Figure 23 shows the difference between the two by plotting bars between the two sets of results. The variation in redshift is small (much smaller than the statistical error bars), except for region 8 in Coma.

Fig. 23. How further systematic effects change the obtained Perseus (left panel) and Coma (right panel) redshifts. For each region, the ends of the leftmost (red) bars are the results for the joint and for the total spectral fits, using the APEC model. The ends of the centre (blue) bars are the results using the APEC and using the SPEX spectral models, when fitting the total spectra. The right (green) points show the effect of using a continuum for the SPEX model which is for a temperature which is 10 per cent greater or less than the line temperature, fitting the total spectra.

Fig. 24. Total spectra and APEC fits for the Perseus cluster regions. Top panels: spectra (blue), total best-fitting model (orange), foreground components (green) and background components (purple) in units of 10−1 cts−1 keV−1. Bottom panels: ratio between the spectrum and model.

Fig. 25. Total spectra and APEC fits for the Coma cluster regions. Top panels: spectra, total fit, foreground fit and background fit in units of 10−1 cts−1 keV−1. Bottom panels: ratio between the spectrum and model.

Systematic errors may also be caused by the choice of spectral model. Inaccuracies in modelling may lead to shifts in the convolved line shape. We plot the results for the APEC and SPEX models when fitting the total spectra in Fig. 23. There is little difference between the results using the two models, suggesting this effect is not important. However, deviations from a thermal plasma could theoretically affect our results.

Another effect which may give systematic errors is if the temperatures obtained by spectral fitting are systematically too large or small. This may occur, for example, if instrument calibration inaccuracies lead to a too flat or too steep spectrum, changing the best-fitting temperature. Schellenberger et al. (2015) found that the temperatures of hot galaxy clusters were different in Chandra and XMM data, probably caused by calibration differences. However they found that this effect is relatively small in hard energy bands like we have used, with offsets of less than 10 per cent for high-temperature systems.

To test for errors caused by incorrect temperatures we created modified versions of our SPEX table model, where for a particular line temperature we calculated the model continuum using a temperature which was 10 per cent greater or smaller than this. We refitted the spectra with these two models, showing the difference between them and the standard model in Fig. 23. The effect of continuum inaccuracies appears small for these data. We found that the obtained temperatures changed by around 3 per cent from the temperatures using the standard SPEX model, indicating that the obtained temperature was dominated by the spectral line shape rather than the continuum.

Another potential systematic is that we previously assumed a slope for the powerlaw background component which was fixed. To check the effect of this assumption, we refitted the spectra using powerlaws which differed in their index by 0.2. The changes in the obtained redshifts was very small, with a maximum redshift difference of 10⁻⁴, and an average of 2 × 10⁻⁵.

Our analysis does not apply a heliocentric correction for the motion of the XMM spacecraft in the solar system to the energy scale. This correction, however, is relatively small, as the maximum heliocentric velocity is ∼30 km s⁻¹.

The fluorescent X-ray emission line originates below the detector, unlike the astrophysical and calibration lines which come from the top. The detector could have a different response to these two X-ray sources. Figure 11 suggests that this effect is not important, as we see consistent averages and widths of residual correction factor between the astrophysical observations (in quadrants) and calibration observations (for each CCD).

These results demonstrate that the above systematic effects are small compared with the measurement and energy calibration uncertainty.

6.2. Future improvements to the method

There are several possible future improvements that could be made to map velocities in extended objects with EPIC-pn. Firstly, it should be possible to calibrate the energy of events made of doubles, where the signal is split between two pixels. Including doubles would approximately increase the number of X-ray events by 80 per cent, although they would have poorer spectral resolution (by up to 10 per cent for Fe-K) than the single-pixel events.

We have not attempted to model the Al-Kα lines around 1.5 keV from the calibration source. This would probably not be directly useful for analysing Fe-K from clusters, given the large difference in energy. However, monitoring the shift of this line may provide more information about the change in energy scale we see at larger energies, allowing us to model its effects better. The Al-Kα is also a detector background line. We also could examine the fluorescent detector lines from Ti, V, Cr, Au and Mo, although these are much weaker than the other lines we used. The Mo-Kα line is at 17.5 keV and so could provide an extreme reference value for the energy scale. Ti, V and Mo is also visible in the Cu hole, so would allow some sort of calibration for the central part of the detector.

Furthermore, it may be possible to calibrate the Cu hole region using calibration observations. Calibrating the central region would allow the analysis of a many archival observations where the cluster is placed in the centre. If we corrected the overall gain of the detector using Cu-Kα as we currently do, we could fit for gain and the energy correction factor in the centre using Mn-Kα and β in calibration observations. This technique could be tested in other regions of the detector where we have the other reference lines. It may, however, be difficult to constrain these parameters in later times in the mission, when the calibration source has weakened.

One could redo the EPIC-pn energy calibration from scratch, working from raw PHA values to build up a complete model of the energy scale and avoiding the multiple correction factors currently employed. This would probably be a very time-consuming project, however.

Observations where the target is offset from the centre of the detector best make use of our data analysis method. We have obtained data for the Virgo and the Centaurus cluster using targeted observations of this nature and plan to publish them in future works.

Although XRISM with its high-spectral resolution X-ray spectrometer should launch in a few years, XMM EPIC-pn has the advantage of a much larger field of view, allowing larger fractions of a cluster to be probed without many pointings. In addition, the better angular resolution of XMM enables point sources to be more easily excluded from the analysis and for different structures to be separated.

7. Conclusions

We describe and demonstrate a novel technique to use background X-ray lines seen in the spectra of the XMM EPIC-pn detector, calibrating the absolute energy scale of the outer part of the detector to better than 150 km s⁻¹ at Fe-K. Using this technique we map the bulk velocity distribution of the ICM over a large fraction of the central Mpc of two nearby clusters of galaxies, Perseus and Coma.

We are able to independently confirm our calibration procedure using a a 65 kpc-square region in Perseus, where we find a velocity consistent with results using an overlapping pointing by Hitomi. There is no evidence of significant motions over much of the cluster, and our average velocity measurement is close to the Hitomi central pointing. However, to the east of the cluster core we observe a shift in velocity of 480 ± 210 km s⁻¹ (statistical). This is spatially coincident with a cold front, seen as an edge in surface brightness and temperature and is direct evidence of gas sloshing in the potential well of the cluster. Excluding this region, we find that the width (σ) of the bulk velocity distribution is 20−150 km s⁻¹, if modelled with a Gaussian distribution. This width is consistent to the Hitomi line width measurements in the core of the cluster, excluding the regions containing strong feedback.

In the Coma cluster, we find that the velocity of the gas is close to the optical velocities of the two central galaxies. The ICM to the north and east of the core has a velocity close to that of NGC 4874, while that to the south and south-west has a lower velocity similar to NGC 4889. The NGC 4911/4921 group, which may have passed through the centre of the cluster, has a consistent optical velocity with some of the higher velocity regions to the north and east. Cooler and denser filaments which were thought be stripped from the subgroup, appear, however, to have a substantially lower velocity close to that of NGC 4889. The velocity of material located near the merging subgroup NGC 4839 is consistent with its optical velocity, although the measurement uncertainties are large.

Our results highlight the differences between a relatively relaxed galaxy cluster, such as Perseus, and the product of an ongoing merger, such as Coma. In the relaxed system, where Hitomi previously determined a low amount of turbulence, we see a narrow range in velocity, except for some evidence of the ICM sloshing in the potential well at the level of a few hundred km s⁻¹. In the merging system, where one would expect high levels of turbulence, we see a large ∼1000 km s⁻¹ range in velocity, with significant structure in our maps.

Acknowledgments

JSS thanks M. Freyberg for helpful discussions. HRR is supported by an STFC Rutherford Fellowship and an Anne McLaren Fellowship. CP is supported by European Space Agency (ESA) Research Fellowships. ACF acknowledges ERC Advanced Grant 340442. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.

References

Appendix A: Additional consistency checks

As there are multiple observations for each object, we can check for consistency between them in different regions. For these comparisons we employed the same spectral modelling as used in the astrophysical cluster analysis.

There are two observations of the Perseus cluster, 0085110101 and 0305780101, which are aimed close to the central nucleus (with an offset of 10s of arcsec) and have similar roll angles. We extracted spectra from five different regions. These include a total comparison region made of an annulus centred on the coordinate 03:19:56.9, +41:31:35 between radii of 6.70 and 12.75 arcmin, designed to exclude the Cu-hole and include most of the detector. We also made four further regions consisting of this annulus split into the different detector quadrants. We fitted these spectra in each region with the APEC code plus background model, allowing independent redshifts and normalisations for the observations, but jointly fitting the temperature and metallicity.

Figure A.1 shows the Cu-Kα redshifts (top panel) and astrophysical redshifts (bottom panel) for each of these different regions. The corrected Cu-Kα redshifts for each of the observations are consistent with zero. The astrophysical redshifts for each of the observations should also be consistent in each of the regions.

Fig. A.1. Comparison of redshift between the two central observations of Perseus, 0085110101 and 0305780101, for different detector regions. (Top panel) Comparison of Cu-Kα redshifts after correction, with an estimated systematic of 75 km s−1. (Bottom panel) Comparison of astrophysical redshifts, with two sets of error bars showing the statistical uncertainty on the redshift after correcting for the Cu-Kα best-fitting redshift, and the error bars extended by the uncertainty on the Cu-Kα redshift. Shown is the systematic of 150 km s−1.

The results agree well between the observations, except for a ∼2σ difference for CCDs 4–6, ignoring the systematic uncertainty. The original disagreement may be caused by the short 16 ks exposure time for CCDs 4–6 (the exposure for CCDs 1–3 is 31 ks). Including the systematic uncertainties on the two data points makes the difference insigificant.

For the Coma cluster we cannot make exact consistency tests as the observations have a different pointing positions, so any differences may be caused by astrophysical signal. However, despite this, we compared the redshifts for each observation in the extraction regions 1–10 (Fig. 19) and looked for any additional evidence of variation within the region. In this analysis, we jointly fitted the temperature and abundance for each spatial region, but allowed the redshifts and normalisations to be free for each observation. Figure A.2 (left side) shows the obtained Cu-Kα redshifts after correction. The horizontal axis shows the observation, indexed by its order in the Coma observations listed in Table 1. As during the joint fits of astrophysical spectra, we ignored those observations for which there were less than 2000 counts between 4.00 and 9.25 keV in the extraction region. Similarly, Fig. A.2 (right side) shows the astrophysical redshifts for the regions and observations.

Fig. A.2. Cu-Kα (left side) and astrophysical (right side) redshifts for regions 1–10 in the Coma cluster, as a function of observation. The dataset number is the number of the observations in Table 1. Systematic uncertainties of 75 and 150 km s−1, respectively, are shown as a grey shaded region.

Both the Cu-Kα and astrophysical redshifts appear consistent for each region. We note again that for the astrophysical results, there may be real astrophysical differences in these datapoints. There are no obvious trends observed with observation number, however. To quantify the scatter for the astrophysical redshifts, we modelled the additional width of points beyond the uncertainties for each dataset in a region by a Gaussian model. We increased the statistical uncertainties on the astrophysical redshift input datapoints by the uncertainties on the Cu-Kα redshift and by a systematic of 150 km s⁻¹, adding in quadrature. Figure A.3 shows the posterior probability distribution on the Gaussian σ, obtained from a Markov Chain Monte Carlo analysis assuming a Gaussian distribution. There appears to be no significant evidence of a width beyond the statistical uncertainties in any region. In regions 2, 5 and 9 there is very marginal evidence of additional spread, which may be because of the velocity structure seen in the maps (Fig. 21).

Fig. A.3. Posterior probability density distributions for the additional velocitity width (σ) between observations, for regions 1–10 in the Coma cluster.

All Tables

All Figures

Fig. 1. Geometry of the EPIC-pn detector. Shown are the numbers of the CCDs, the field of view exposed to the sky (FoV), the raw coordinates (RAWX and RAWY), the detector coordinates (DETX and DETY) and the approximate location of the “Cu hole”. The position angle (PA) is the angle between the detector and sky coordinates, an example of which is shown. The four quadrants are shaded in different colours.

In the text

	Fig. 2. Total stacked spectra of the fluorescent background emission lines for 32.6 Ms of FF (red) and 19.7 Ms EFF (blue) data. Before stacking the individual spectra, the three orders of correction described in Sect. 2 were applied. Only events with PATTERN==0 were included. Au-Lα was not used in our analysis because of its weakness.
In the text

Fig. 3. Images of the detector in different fluorescent lines, created from stacked astrophysical observations. Shown are images in Cu-Kα (left; 7.805 to 8.285 keV), Ni-Kα (centre; 7.280 to 7.680 keV) and Cu-Kβ and Zn-Kα (right; 8.455 to 9.075 keV). Approximate background has been subtracted using scaled continuum images on either side of the lines. Some residual astrophysical emission can be seen in the centre of the image. The total exposure is around 64 Ms after cleaning for flares, combining both FF and EFF data and using PATTERN ≤ 4. The image was binned by 64 detector coordinate units and smoothed by a Gaussian of σ = 1 pixel width.

In the text

	Fig. 4. Cumulative FF (red) and EFF (blue) cleaned exposure of the astrophysical datasets used in the calibration as a function of time (in XMM revolutions).
In the text

	Fig. 5. Best-fitting Cu-Kα redshifts for individual astrophysical observations as a function of time (in XMM revolutions). A 0.1 per cent correction is 8 eV at the energy of Cu-Kα.
In the text

Fig. 6. Redshift of the Cu-Kα line (in per cent) as a function of detector position over the mission lifetime in FF (left) and EFF (right) modes, after correcting for the average value. The range shown in the top-left of each panel is the XMM revolution. The shift has been measured in 1 × 20 pixel regions and a σ = 2 pixel Gaussian smoothing applied for display. Detector regions with a large uncertainty in a combined error map have been masked out.

In the text

Fig. 7. Intrinsic Gaussian centre and width (σ) of the distribution of Cu-Kα redshifts in each CCD for the observations which went into the stacking analysis, after correcting each observation for the average Cu-Kα value and the time-dependent maps (Fig. 6). This analysis excludes observations with statistical uncertainties on the Cu-Kα redshift of greater than 0.1 per cent. Results for FF and EFF modes are shown separately.

In the text

Fig. 8. Demonstration of the necessity of an energy scale correction. Shown in black are the total spectra of calibration observations between revolutions of 1000 and 1500 (<7 keV) and the total spectra of calibration observations plus stacked observations in the same time period (>7 keV), outside the Cu hole region. The count axis for data >7 keV has an arbitrary scaling. Top panel: models and data, while the bottom panel shows their ratio. There are two “6-line” models. One (blue) assumes that the lines are at their standard energies. The other (pink) has a linear compression model for the line energies about Cu-Kα where δE/E = 0.0013 keV−1. The poorly-modelled line wings do not affect the best-fitting compression factors.

In the text

Fig. 9. Energy scale correction factor for each CCD as a function of time. Positive values mean that events have too low an energy at high energies and too high an energy at low energies. This plot shows the results of fitting the “6-line” model to stacked calibration observations and observational datasets after gain correction, excluding the Cu hole. The lines show linear plus constant fits to the results excluding the first revolution bin.

In the text

Fig. 10. Maps of the residual energy correction factor in stacked observations, after applying our correction to this in the individual observations. The units are per cent per keV. The panels show revolution bins of 0–500, 500–1000, 1000–1500, 1500–2000, 2000–2500 and 2500–3000 (left to right and top to bottom). The uncertainty of these values is not the same across the detector or in time, so we use 8 and 16 pixel binning. The statistical uncertainty is much larger in the bottom-left area because of the non-uniformity of the calibration source. Pixels where the uncertainty is larger than 10−3 keV−1 are not shown.

In the text

Fig. 11. Distribution of residual correction factors (centre and width) obtained for each CCD, after correcting the calibration observations and refitting them. These results primarily use the Mn-Kα and β lines. Also shown are the distribution of correction factors obtained in each quadrant (three CCDs) for the individual astrophysical observations after correction, obtained from the Ni-Kα, Zn-Kα and Cu-Kβ lines.

In the text

Fig. 12. Regions examined in Perseus, chosen to cover areas with similar numbers of counts in the Fe-K complex. (Left panel) Regions plotted on the fractional difference in 0.5 to 2.0 keV surface brightness from the average at each radius. Region 22 overlaps with other regions as it covers the outer Hitomi pointing. (Top-right panel) Fe-K count map, showing the number of counts in each 3.2 arcsec pixel in the Fe-K complex (6.50 to 6.90 keV, restframe), after subtracting neighbouring scaled continuum images (6.06 to 6.43 and 6.94 to 7.16 keV, restframe). Gaussian smoothing of σ = 4 pixels was applied. (Bottom-right panel) Exposure map, plotting the total Fe-K exposure (log s) excluding the Cu hole and accounting for vignetting.

In the text

Fig. 13. Results for the Perseus regions. Shown are the best-fitting values and statistical uncertainties for the joint fit using the APEC model. (Left panel) Redshifts and velocities. The average shown is the weighted average for the non-overlapping regions 1 to 21. The right axis shows the redshifts as velocities relative to the H16 result (z = 0.01756). The NGC 1275 velocity of z = 0.017284 is from H18. Region 22 was chosen to overlap with the outer Hitomi pointing, for which we also show their result. (Right panels) Temperatures and metallicities. Also shown are the values for the outer Hitomi pointing from Hitomi Collaboration (2018d) and Simionescu et al. (2019).

In the text

Fig. 14. Example fit of the total spectrum for region 16 in Perseus, rebinned for display to have a signal to noise ratio of 10 in each spectral bin. The vertical line is at an energy of 6.55 keV. (Top panel) Fitted data and best-fitting cluster model (fitted redshift), showing its cluster and background components. Also shown is the cluster component with the redshift changed to the Hitomi cluster value. (Centre panel) Residuals of the total best-fitting model. (Bottom panel) Residuals of the total model with the APEC redshift changed to the Hitomi value.

In the text

	Fig. 15. Change in fit statistic as a function of redshift for selected regions in the Perseus cluster, for the total spectral fits.
In the text

Fig. 16. Perseus cluster maps. Shown is the 0.5 to 2 keV X-ray surface brightness and contours (log cts s−1 per 1.6 arcsec pixel; top left), fractional difference of surface brightness to radial average (top right), temperature (keV; centre left), metallicity (Z⊙; centre right), velocity relative to the cluster value in H16 (km s−1; bottom left) and 1σ statistical uncertainty on velocity (km s−1; bottom right). Maps were created by a moving 3:1 elliptical region (rotated to lie tangentially to the nucleus) containing ∼700 Fe-K counts (see example ellipses in velocity error plot). The pixels are not statistically independent.

In the text

	Fig. 17. Radial profiles of temperature, metallicity and velocity pixels from the Perseus maps in Fig. 16 coloured by position angle. Error bars are only shown for every 10 points. Radii are from the position of 3C 84. Nearby points are not statistically independent.
In the text

Fig. 18. Comparison of a sloshing simulation (top row) and data for the Perseus cluster (Fig. 16; bottom row). The columns show surface brightness residuals (left), velocity (centre) and temperature (right). The simulation is of a Perseus-like cluster which has bulk sloshing motions and turbulence initiated by an off-centre encounter with a subcluster of 1/5 its mass, from ZuHone et al. (2010, 2018).

In the text

	Fig. 19. Regions examined in Coma. (Left panel) The regions plotted on an exposure-corrected surface brightness image in the 0.5 to 2.0 keV band. (Top-right panel) Image in the Fe-K band, after subtracting continuum, in 10−3 counts per 1.6 arcsec pixel. (Bottom-right panel) Fe-K exposure map in log s, including the effect of vignetting and the Cu hole.
In the text

	Fig. 20. Results for the Coma regions. Shown are the best-fitting values for the joint fits using the APEC model. (Left panel) Redshifts and velocities. The right axis shows the redshifts as velocities relative to the redshift of NGC 4889. (Right panels) Temperatures and metallicities.
In the text

Fig. 21. Coma cluster maps. Shown is the 0.5 to 2.0 keV X-ray surface brightness with contours (log cts s−1 per 1.6 arcsec pixel; top left), fractional difference of surface brightness to radial average (top right), temperature (keV; centre left), metallicity (Z⊙; centre right), velocity relative to NGC 4889 (km s−1; bottom left) and 1σ statistical uncertainty on velocity (km s−1; bottom right). Spectral results were obtained from spectra extracted from a sliding 4 arcmin radius extraction region (circle) over a grid. Regions where the uncertainty on the temperature is >1 keV, metallicity >0.07 Z⊙ or velocity >600 km s−1 are not shown. The position of the four central lensing-map peaks of Okabe et al. (2010) are indicated by squares.

In the text

	Fig. 22. Amplitude of gas density fluctuations in the Coma cluster as a function of the wave number k = 1/ℓ. The curves show the amplitude of projected surface brightness fluctuations (red) and the amplitude of 3D density fluctuations (blue). The uncertainties are indicated by the shaded areas.
In the text

Fig. 23. How further systematic effects change the obtained Perseus (left panel) and Coma (right panel) redshifts. For each region, the ends of the leftmost (red) bars are the results for the joint and for the total spectral fits, using the APEC model. The ends of the centre (blue) bars are the results using the APEC and using the SPEX spectral models, when fitting the total spectra. The right (green) points show the effect of using a continuum for the SPEX model which is for a temperature which is 10 per cent greater or less than the line temperature, fitting the total spectra.

In the text

	Fig. 24. Total spectra and APEC fits for the Perseus cluster regions. Top panels: spectra (blue), total best-fitting model (orange), foreground components (green) and background components (purple) in units of 10−1 cts−1 keV−1. Bottom panels: ratio between the spectrum and model.
In the text

	Fig. 25. Total spectra and APEC fits for the Coma cluster regions. Top panels: spectra, total fit, foreground fit and background fit in units of 10−1 cts−1 keV−1. Bottom panels: ratio between the spectrum and model.
In the text

Fig. A.1. Comparison of redshift between the two central observations of Perseus, 0085110101 and 0305780101, for different detector regions. (Top panel) Comparison of Cu-Kα redshifts after correction, with an estimated systematic of 75 km s−1. (Bottom panel) Comparison of astrophysical redshifts, with two sets of error bars showing the statistical uncertainty on the redshift after correcting for the Cu-Kα best-fitting redshift, and the error bars extended by the uncertainty on the Cu-Kα redshift. Shown is the systematic of 150 km s−1.

In the text

	Fig. A.2. Cu-Kα (left side) and astrophysical (right side) redshifts for regions 1–10 in the Coma cluster, as a function of observation. The dataset number is the number of the observations in Table 1. Systematic uncertainties of 75 and 150 km s−1, respectively, are shown as a grey shaded region.
In the text

	Fig. A.3. Posterior probability density distributions for the additional velocitity width (σ) between observations, for regions 1–10 in the Coma cluster.
In the text