A&A 413, 415-439 (2004)
DOI: 10.1051/0004-6361:20031512

Spatially resolved X-ray spectroscopy of cooling clusters of galaxies

J. S. Kaastra 1 - T. Tamura 1 - J. R. Peterson 2 - J. A. M. Bleeker 1 - C. Ferrigno 1 - S. M. Kahn 2 - F. B. S. Paerels 2 - R. Piffaretti 3,4 - G. Branduardi-Raymont 5 - H. Böhringer 6


1 - SRON National Institute for Space Research Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands
2 - Department of Physics, Columbia University, 550 West 120th Street, New York, NY 10027, USA
3 - Paul Scherrer Institute, Laboratory for Astrophysics, 5232 Villigen, Switzerland
4 - Institute of Theoretical Physics, University of Zürich, Winterthurerstrasse, 190, 8057 Zürich, Switzerland
5 - Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, UK
6 - Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse, 85748 Garching, Germany

Received 8 January 2003 / Accepted 22 September 2003

Abstract
We present spatially resolved X-ray spectra taken with the EPIC cameras of XMM-Newton of a sample of 17 cooling clusters and three non-cooling clusters for comparison. The deprojected spectra are analyzed with a multi-temperature model, independent of any a priori assumptions about the physics behind the cooling and heating mechanisms. All cooling clusters show a central decrement of the average temperature, most of them of a factor of ${\sim} 2$. Three clusters (Sérsic 159-3, MKW 3s and Hydra A) only show a weak temperature decrement, while two others (A 399 and A 2052) have a very strong temperature decrement. All cooling clusters show a weak pressure gradient in the core. More important, at each radius within the cooling region the gas is not isothermal. The differential emission measure distribution shows a strong peak near the maximum (ambient) temperature, with a steep decline towards lower temperatures, approximately proportional to T3, or alternatively a cut-off at about a quarter to half of the maximum temperature. In general, we find a poor correlation between radio flux of the central galaxy and the temperature decrement of the cooling flow. This is interpreted as evidence that except for a few cases (like the Hydra A cluster) heating by a central AGN is not the most common cause of weak cooling flows. We investigate the role of heat conduction by electrons and find that the theoretically predicted conductivity rates are not high enough to balance radiation losses. The differential emission measure distribution has remarkable similarities with the predictions from coronal magnetic loop models. Also the physical processes involved (radiative cooling, thermal conduction along the loops, gravity) are similar for clusters loops and coronal loops. If coronal loop models apply to clusters, we find that a few hundred loops per scale height should be present. The typical loop sizes deduced from the observed emission measure distribution are consistent with the characteristic magnetic field sizes deduced from Faraday rotation measurements.

Key words: galaxies: clusters: general - cooling flows - X-rays: galaxies: clusters

1 Introduction

The visible mass in clusters of galaxies is dominated by hot diffuse gas. Due to the high temperature of the gas it is visible in the X-ray band. In the central part of the cluster the gas density is high. There the radiative cooling time is shorter than the age of the cluster and the age of the universe. The temperature decreases due to radiative cooling. The corresponding pressure decrease then causes a net inflow towards the center of the cluster. A review of this so-called cooling flow process is given by Fabian (1994). The associated mass deposition rate $\dot{M}$ may reach several thousands  $M_{\hbox{$\odot$ }}$ per year in the strongest cooling flows. In general, the time the gas needs to cool down to very low temperatures is smaller than the flow time towards the center of the cluster. Thus, gas drops out locally from the cooling flow.

In the isobaric cooling flow model (Johnstone et al. 1992) gas cools in pressure equilibrium with its surroundings. The amount of gas at each temperature is controlled by how fast it can cool, i.e. the radiative cooling rate through - mainly - X-ray radiation. The mass deposition rate is not the same at all radii but scales linearly with the radius, $\dot{M}\sim r$.

The first high-resolution X-ray spectra of clusters of galaxies taken with the Reflection Grating Spectrometers (RGS) of XMM-Newton showed indeed the presence of cooler gas in the cores of several clusters. However, the amount of cool gas at lower temperatures was much smaller than predicted by the isobaric cooling flow model (Sérsic 159-3, Kaastra et al. 2001; A 1835, Peterson et al. 2001; A 1795, Tamura et al. 2001a). The second paper also contains a list of possible explanations for this phenomenon.

This discovery has triggered a series of papers offering a broad range of explanations for the apparent failure of the isobaric cooling flow model. Worked-out examples include metallicity inhomogeneities (Fabian et al. 2001), buoyantly rising radio bubbles transporting cool gas outwards (Churazov et al. 2001), halo-in-halo structure (Makishima et al. 2001), turbulent mixing due to rising and falling hot gas bubbles (Quilis et al. 2001), heating by AGN activity (Böhringer et al. 2002), contamination due to non-thermal X-ray emission (McCarthy et al. 2002), heating due to dead radio galaxies (Reynolds et al. 2002), rapid cooling due to mixing with cold gas (Fabian et al. 2002), and heat conduction by electrons (Voigt et al. 2002).

The lack of relatively cool gas has been confirmed by the RGS spectra of other clusters, like A 496 (Tamura et al. 2001b) and Virgo (Sakelliou et al. 2002). The lack of this cool gas is deduced from both RGS data as well as from XMM-Newton EPIC data. These EPIC data have also been reported in the RGS papers mentioned above. In addition, Molendi & Pizzolato (2001) found evidence for an apparent low temperature cut-off in the cooling gas at about 1-2 keV in the EPIC spectra of A 1835, A 1795 and Virgo. Matsushita et al. (2002) report that the EPIC data of Virgo show that the cool gas has a single temperature at each radius.

Also Chandra observations have confirmed the lack of cool gas in clusters. In Hydra A (David et al. 2001) the gas has essentially a single temperature at each radius between 30-200 kpc, with a steadily outwards temperature increase. Only in the innermost 30 kpc there is evidence for multiphase gas, however with a ten times smaller mass deposition rate as derived from the morphologically derived mass accretion rate at 30 kpc.

The superb energy resolution of the RGS allowed a detailed investigation of the temperature structure of the core. Although the RGS has spatial resolution in the cross dispersion direction of the gratings, for most of the available data the statistics are not good enough to map the cool gas in the cross dispersion direction. In order to be able to distinguish between the various theoretical models that have been proposed, it is important to measure the temperature structure for each radius. As shown above, for some clusters a single-phase cooling gas distribution at each radius has been measured (Hydra A, David et al. 2001; Virgo, Matsushita et al. 2002; A 2029, Lewis et al. 2002).

In this paper we study a large sample of clusters and use the spatially resolved XMM-Newton/EPIC spectra to study the temperature structure at each radius. These measurements allow us to distinguish between isothermal, isobaric cooling or cut-off cooling flow models.

The layout of this paper is as follows. In Sect. 2 we give a description of our data analysis procedures. This description is extensive since it also forms the basis of a set of papers on other properties of these clusters that are in preparation (e.g., soft excess emission (Kaastra et al. 2003); abundances (Tamura et al., in preparation). In Sect. 3 we describe the results of our spectral modeling using several models, and in Sect. 4 we compare our results to a range of physical models for the cooling gas. Throughout this paper we use H0 = 50 km s-1 Mpc-1 and q0=0.5.

  
2 Data analysis

2.1 Cluster sample

Our sample consists of 17 clusters with cooling flows, taken from the Guaranteed Time program of the RGS consortium and other available data sets. These cooling flow clusters were selected mainly based upon their suitability for spectroscopic studies with the Reflection Grating Spectrometers (RGS) of XMM-Newton. This selection criterion favors relatively medium distant, compact and cool clusters. Nearby clusters in general have a large angular size, resulting in a degradation of the spectral resolving power of the RGS. In very hot clusters the line emission is weak with respect to the Bremsstrahlung continuum. In addition we added three non cooling clusters as a control case for clusters without cooling gas: the Coma cluster and A 3266, as well as the merging cluster A 754. The cluster sample is listed in Table 1. We list in this and all subsequent tables the clusters in increasing hot gas temperature order.

Table 1: Basic cluster parameters.

Table 2: Observation log.

Temperatures in Table 1 were taken from Ebeling et al. (1998) for NGC 533 and A 1835, from Finoguenov et al. (2001) for MKW 9 and from David et al. (1993) for all other clusters. Unabsorbed X-ray fluxes (0.1-2.4 keV) were taken from Ebeling et al. (1996) for A 426, A 496, A 780, A 1837, and A 4059; from Cruddace et al. (2002) for Sérsic 159-3; from Kriss et al. (1983) for MKW 9; and from Ebeling et al. (1998) for all others. Note that the data presented for A 189 by Ebeling et al. (1998) correspond to NGC 533 (as analyzed in this paper), not to A 189.

2.2 Data extraction

Details about the observations and data processing are given in Table 2. For A 1835 we discarded the MOS1 data (taken in a large window mode). The data were extracted using the standard SAS software, equivalent to version 5.3. Only events with pattern < 13 (MOS) or pattern =0 (pn) were accepted; in addition, we selected flag =0.

  
2.3 Background subtraction

The background spectrum of the XMM-Newton EPIC camera's is described in detail by Lumb et al. (2002). The main components are the cosmic X-ray background, dominating the background at low energies, and a time-variable particle component, dominating the high energy background. Background subtraction is always done using the same detector regions for source and background.

2.3.1 Time-variable particle background

The time-variable particle component is caused by clouds of soft protons. Apart from long quiescent intervals, there are sporadic episodes where this background component is enhanced by orders of magnitude, as well as time intervals with only a low level of flaring.

The contribution of X-ray photons to the full field 10-12 keV band count rate is almost negligible as compared with the background. This energy band is dominated by low energy protons. Therefore it can be used as a monitor for the time-variable background component.

Thus we divided the observation into time intervals of 260 s, and for each time interval H, the number of counts in the 10-12 keV band was determined. During quiescent periods, the average value for H is typically 23-24 for MOS and 31 for pn. We excluded time intervals with high particle background as measured by $H>H_{\rm t}$ where for most observations we took the threshold $H_{\rm t}=35$for MOS and $H_{\rm t}=50$ for pn. This threshold is about 2.5 and 3.5$\sigma$ above the average quiescent level, respectively. Therefore the fraction of time bins rejected due to Poissonian fluctuations in H is very small. For some clusters we had to take slightly higher thresholds in order to retain a minimum of exposure time. The thresholds used as well as the background after filtering are listed also in Table 2.

The above procedure is sufficient to remove the large flares, simply by discarding the time intervals with high 10-12 keV count rates. However after removing these intense large proton events the average background level still varies slightly. The rms variation from cluster to cluster is $\sim $10-20%. This is due to the fact that for some observations the count rate is slightly enhanced but smaller than the cutoff threshold. This holds in particular for the clusters with a higher threshold but also to some extent for the others. These remaining weak flares cannot be removed by time selection. For them, however, the increase in background flux in the 0.2-10 keV raw spectrum is to first order proportional to the increase in 10-12 keV count rate, as we illustrate below.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f1.ps}\end{figure} Figure 1: Background spectra for MOS1. Solid histogram: background spectrum for H between 19-28; squares, triangles and stars: difference of background spectra for H between 29-38, 39-48 and 49-58, respectively, with the 19-28 spectrum.
Open with DEXTER

We show in Fig. 1 the full-field MOS1 spectrum of our background field (discussed in the next section) for all 977 bins with H between 19-28, close to the average value of 23.82. The soft excess due to the cosmic X-ray background as well as the instrumental fluorescence lines near 1.5 keV are clearly seen. We also show the difference spectra as compared to this "average'' spectrum for H ranges of 29-38 (average value 31.84, for 249 time bins); 39-48 (average value 42.35, for 37 time bins); and 49-58 (average value 53.17, for 12 time bins). These difference spectra are smoother than the average background spectrum, and they are approximately exponential ${\sim}\exp({-E/E_{\rm c}})$, with $E_{\rm c}$ values of typically 4.2, 4.5 and 5.2 from low to high background values, respectively. Thus, as a first approximation, for small enhancements of the background level, the shape of the background enhancement as a function of energy is more or less fixed while its amplitude scales with the difference of H with the average value of H.

We utilize this behaviour to make a first order correction to the estimated background spectrum. The usual procedure would be to subtract the background spectrum using a blank field observation screened in the same way as the cluster data. Here in order to correct for the weak variability we divide the cluster and blank field observation into subsets with the same 10-12 keV count rate (same H). The subtracted background B(E) was then determined as

\begin{displaymath}B(E) = \sum\limits_{H}^{} B_{H}(E) [f_{\rm c}(H)/f_{\rm b}(H)][t_{\rm c}/ t_{\rm b}],
\end{displaymath} (1)

where BH(E) is the histogram of background events at a given 10-12 keV count rate level H in the blank field, $f_{\rm c}(H)$ and $f_{\rm b}(H)$ are the frequency distributions of H for the cluster and blank field, respectively, and $t_{\rm c}$and $t_{\rm b}$ are the exposure times of the cluster and blank fields. This ensures correct particle background subtraction, under the assumption that for low proton count rates the shape of the proton spectrum does not vary when its flux increases by a small amount. However in order to account for any possible remaining particle background subtraction problems, we have included a systematic uncertainty of at least 10% of the total background in all our fits.

The deep field used in our background subtraction procedure (Lumb et al. 2002) consisted of the background event file provided by the XMM-Newton Science Operations Center, containing a total of 320-410 ks exposure time on 8 deep fields (the exposure time differs slightly for the different instruments). This background event file was filtered the same way as the source file, except that we took the threshold $H_{\rm t}$ slightly higher in order to have some estimate of the background for mildly enhanced soft proton levels (see Table 2).

As stated above, the contribution of X-ray photons to the 10-12 keV band count rate is almost negligible. There are three possible exceptions to this: Coma, Virgo and Perseus.

For Coma, the second brightest cluster in our sample and a very hot cluster, the full-field 10-12 keV cluster emission reaches a level of 10% of the quiescent background count rate in the pn camera only. Owing to the lower sensitivity of the MOS cameras at high energies, the relative cluster contamination in the 10-12 keV band is at most 5% there. Thus, we expect that our method overestimates the time-variable particle contribution in Coma by 5-10%. This does not affect our science results for Coma, however, since (i) the difference is within the systematic background uncertainty used by us; (ii) at lower energies the background is dominated by the constant cosmic X-ray background, and not by the soft proton component; (iii) except for the outermost annulus, the X-ray flux of Coma is a factor of 5-100 times brighter than the subtracted background for all energies below $\sim $6 keV; (iv) the spectral fit is dominated by the higher signal-to-noise part of the spectrum below $\sim $6 keV.

Although Virgo is brighter than Coma (Table 1), its temperature is so much lower that the amount of 10-12 keV photons can be neglected safely with respect to the soft proton background.

Only Perseus has really an enhanced 10-12 keV count rate, as is evident from the large average values of H as listed in Table 2. Here we subtracted first a constant offset from the 10-12 keV light curves, to be consistent with the average background level in the deep fields. Note that Churazov et al. (2003) suggest that for this observation there is a steady enhanced soft proton background. However, Perseus is so bright that in none of our spectra the subtracted background is important. Moreover, our results for the cooling region are most sensitive to Fe-L emission, and this is the spectral range where the background is relatively very low.

2.3.2 Cosmic X-ray background

While at high energies the particle contribution dominates the background spectrum, at low energies the cosmic X-ray background yields the largest contribution to the measured background. The cosmic X-ray background varies from position to position on the sky. This is in particular important at low energies. To get a typical estimate for the sky variation, we took the PSPC count rates for 378 regions at high Galactic latitude as studied by Snowden et al. (2000) and determined the population variation for this component in the R1 (low energy) band; the rms variation is typically $35\pm3$%. In principle, the background estimate for any location on the sky can be improved by taking the measured PSPC count rates into account. However, this requires detailed spectral modeling of the soft X-ray background for each position on the sky that is being studied, which can be rather uncertain due to the poor spectral resolution of the PSPC. Here we have taken a conservative approach by using the average background as contained in the standard EPIC background files, but we included a systematic background error of 35% of the total background below 0.5 keV, 25% between 0.5-0.7 keV, 15% between 0.7-2 keV, and 10% above 10 keV. This last 10% also includes the systematic uncertainty in the subtracted time-variable particle background component.

Note that the above 35% systematic background uncertainty at low energies is a rms estimate. In individual cases there may be a larger deviation from the average high latitude X-ray background.

The differences between the spectral response for pn data taken with the full frame mode (as in our deep fields) and the extended full frame mode (as in some of our clusters) are negligible compared to the systematic uncertainties on the background and effective area that we use in this paper. Of more importance is the difference between data taken with the medium filter (Coma and a few other clusters) and the thin filter (most of our clusters and our deep field). For a power law photon spectrum with photon index 2, the medium filter produces about 10% less counts at 0.2 keV, 6% less at 0.3 keV, and less than 2% above 0.4 keV, as compared to the thin filter. These differences are all well within our adopted systematic uncertainties.

Finally, in all fields the brightest X-ray point sources were removed. This was done by making an image, sampled with a pixel size of 4 $^{\prime\prime}$. At each pixel, the net flux within a square of size 5 pixels (20 $^{\prime\prime}$) was determined and compared with the local background as measured in the surrounding square of 17$~\times~$17 pixels (68$~\times~$68 $^{\prime\prime}$). All sources stronger than 4$\sigma$(MOS) or 4 $\sqrt{2}\sigma$ (pn) were marked. Events within a square of 32$~\times~$32 $^{\prime\prime}$ around the marked sources were discarded. This procedure was only done outside a radius of 2$^\prime $ from the cluster center. Typically, a median of 8 point sources (for MOS1 and MOS2) or 15 (for pn) was discarded this way. In general, we detect more point sources in the high signal to noise data. For those clusters point source rejection is indeed relatively important. Our sensitivity is not seriously affected by this selection, since even for Virgo the discarded detector area is only 5% of the total detector area.

  
2.4 Radial profiles and psf effects

The psf of the X-ray telescopes of XMM-Newton has a rather narrow core (FWHM of order 6 $^{\prime\prime}$), but extended wings. In fact, the shape of the psf closely resembles the shape of a typical cluster radial intensity distribution (see for instance Ghizzardi 2001). Thus, spectra accumulated in a given region of the detector are contaminated by flux from other regions of the cluster, which may have different spectral shapes. In this section we make quantitative estimates of this effect. For simplicity we assume here spherical symmetry for both the source radial profile $\mu (r)$ and the instrumental psf. We determine how much the observed radial intensity profile $\phi (r)$ differs from the true radial profile $\mu (r)$. We denote the encircled energy fraction of the psf by F(r), with by definition F(0)=0 and $F(\infty)=1$. The encircled energy F(r) of X-ray telescope 1 (MOS1) can be described quite accurately with a King profile:

 \begin{displaymath}
F(r) = 1 - \left[ 1 + (r/a)^2 \right]^{-b},
\end{displaymath} (2)

where for photon energies of 1.5 keV a=4.8 $^{\prime\prime}$ and b=0.45, and for 6 keV photons a=3.6 $^{\prime\prime}$ and b=0.40. Unfortunately the distribution (2) has no finite variance. For example, if we cut-off all photons outside a radius of 7.5, 15, 30 and 60$^\prime $, the equivalent 1-dimensional rms width $\sigma$ equals 0.61, 0.90, 1.32 and 1.93$^\prime $, respectively (for an energy of 1.5 keV). This is much larger than the core of the psf, and comparable to the typical core radius of the clusters we are interested in. Thus, the psf effects are potentially large.

To illustrate this, we have calculated using a Monte-Carlo simulation with 107 photons the ratio $\phi(r)/\mu(r)$ for MOS1 (Fig. 2). For $\mu (r)$ we took a $\beta $ profile with $\beta =2/3$ (or equivalently $\mu(r)=\left[1+(x/r_{\rm c})^2\right]^{-1.5}$). For $r_{\rm c}=1$$^\prime $, at small rboth profiles differ by 18%, and at large radii (also at 15$^\prime $ off-axis) there is a 10% excess. Even for $r_{\rm c}=3$$^\prime $, there is a significant effect. This is all due to the strong psf tails.

The good news, however, is that the ratio $\phi/\mu$ does not vary more than 1% (for a core radius of 1$^\prime $) when the photon energy is changed from 0 to 9 keV. This indicates that the psf does not introduce significant spectral distortions for the clusters that we consider in this paper. Thus, we use a correction factor $\phi/\mu$ that is independent of energy, and determine it from a single energy band. Therefore the psf correction does not affect the shape of the spectra, only their normalization.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f2.ps}\end{figure} Figure 2: Ratio of psf-convolved cluster profile $\phi (r)$ over true cluster profile $\mu (r)$ for MOS1 and a photon energy of 1.5 keV. For $\mu (r)$ we have taken a $\beta $ profile with $\beta =2/3$ and core radii as indicated.
Open with DEXTER

  
2.5 Deprojecting spectra

We have demonstrated in the previous section that contamination of spectra in different regions by flux from other regions is an effect that should be taken into account in the analysis of these spectra. We describe in the Appendix how this is formally implemented. However as shown in that Appendix simultaneous fitting of all spectra in different annuli taking all relevant effects (vignetting, psf, projection) into account in a forward folding technique is problematic. Therefore we correct and deproject our spectra as follows.

We start with a template X-ray image of the cluster. We have taken for this a vignetting-corrected Rosat PSPC image as obtained from NASA's skyview facility. The constant background level was determined at large radii and subtracted from the image. In the central parts (typically a few arcmin), where the spatial variations are large, this PSPC image was merged smoothly with a background-subtracted Rosat HRI image. The resulting image was then convolved with the psf of the EPIC camera's and from this the ratio $\phi(r)/\mu(r)$ cf. Sect. 2.4 was determined.

The use of a PSPC image as the template for the true cluster image does not introduce much bias in our analysis. This can be seen as follows. For a Gaussian psf $F(r) = 1-\exp\left[ -r^2/2\sigma^2\right]$ the ratio $\phi(r)/\mu(r)$ can be developed in a series as:

 \begin{displaymath}
\phi(r)/\mu(r) = 1 + 0.5\sigma^2 \left[ \mu^{\prime\prime}(r)
+ \mu^{\prime}(r)/r
\right]/\mu(r) + O\left(\sigma^4\right).
\end{displaymath} (3)

Therefore the correction factor is of order $(\sigma / a)^2$ with a the characteristic scale height of $\mu (r)$ (and in the limit of $\sigma\ll a$). For a non-Gaussian psf similar expressions can be obtained that can be cast in the shape of Eq. (3) with an effective value for $\sigma$ that depends upon the details of the psf. The above (3) holds both for convolution of any radial profile with the Rosat ( $\sigma_{\rm R}$) and XMM-Newton ( $\sigma_{\rm X}$) psf. In both cases, $(\sigma / a)^2$ is of the order of ${\sim} 10$%. Our estimate for $\phi(r)/\mu(r)$ is $\mu (r)$ convolved with $\sigma_{\rm X}$ and $\sigma_{\rm R}$, dividided by $\mu (r)$ convolved with $\sigma_{\rm R}$. The true value for $\phi(r)/\mu(r)$ is $\mu (r)$ convolved with $\sigma_{\rm X}$, dividided by $\mu (r)$. It is easy to show that the difference between both expressions is of order $(\sigma_{\rm R} / a)^2 (\sigma_{\rm X} / a)^2$, i.e. on the one percent level. Therefore there is no need to consider higher resolution template images.

Then for the EPIC data we make a radial, background subtracted profile for each energy band. These radial profiles are corrected for exposure (CCD gaps, removed hot pixels, removed point sources etc.). We then apply the psf correction factor $\mu(r)/\phi(r)$, and finally the (energy and position dependent) vignetting correction in order get the exposure, psf and vignetting corrected radial profile for each energy band. Vignetting is both due to the telescope and in case of the MOS cameras also due to the Reflection Grating Array.

From these radial profiles we then construct the deprojected count rates in various shells, assuming spherical symmetry. It can be shown in general that if the spatial emissivity is $\epsilon(r)$, the projected surface brightness $\mu (r)$, and the integrated emissivity y(r) is defined by

\begin{displaymath}y(r) \equiv \int\limits_0^{r} 4\pi s^2 \epsilon(s) {\rm d}s,
\end{displaymath} (4)

then y(r) is related to $\mu$ by:

 \begin{displaymath}
y(r) = y(\infty) - \int\limits_r^{\infty}
2\pi t \mu(t) g(t/r) {\rm d}t,
\end{displaymath} (5)

where

\begin{displaymath}g(x) \equiv \left\{
\begin{array}{ll}
0 & \qquad\mbox{if $x ...
...}\right)} &
\qquad\mbox{if $x \geq 1$ }.
\end{array} \right.
\end{displaymath} (6)

The integration kernel $2\pi t \mu(t){\rm d}t$ is proportional to the number of observed counts in the projected annulus centered around t with width ${\rm d}t$. The difference $Y_{12}\equiv y(r_2)-y(r_1)$ is the number of emitted counts in a spherical shell between r1 and r2. Hence, Eq. (5) gives the (deprojected) number of emitted counts in a spherical shell as a linear combination of the observed number of counts in the (projected) surface brightness profile. We sample the surface brightness $\mu$in circular annuli with thickness $\Delta$ and inner and outer radii $(i-1)\Delta$ and $i\Delta$, and take Ni the number of observed counts in these projected annuli. It follows that

 \begin{displaymath}
Y_{12} = \sum_{i}^{~}
\left\{
g\left( {\frac{(i+1/2)\Del...
...
- g\left( {\frac{(i+1/2)\Delta}{r_2}} \right)
\right\}
N_i.
\end{displaymath} (7)

Evaluating Eq. (7) for all annuli and energy bins, the deprojected spectra for each shell are obtained directly. It is also straightforward to evaluate the statistical error bars on the spectrum. In the choice of the shell boundaries r1 and r2 we have searched a compromise between spatial resolution and signal to noise ratio. Therefore we choose the boundary separations not smaller than 0.5$^\prime $ (for smaller separations, in particular in the core, psf effects cause strong overlap between the neighbouring shells). The shell boundaries used in this work are listed in Table 3. These 9 shells cover approximately the entire field of view of the XMM-Newton telescopes.

Table 3: Boundaries between the shells. For each shell nr., the lower and upper radius (in arcmin) are listed.

A possible complication arises for those clusters that have X-ray emission outside the field of view of the XMM-Newton telescopes. We have not attempted to estimate a correction for this in Eq. (7). For all our clusters, the surface brightness drops rapidly towards large radii. Therefore the projected contribution from shells outside the field of view to the X-ray flux of the innermost shells with cooling gas is negligible. In the outer parts of the clusters, the shape of the X-ray spectrum does not vary rapidly with radius, so neglecting the out-of-field emission only gives a small bias to the normalisation of the spectra of the outermost shells, not to their shape. But as noted, several of our clusters have no significant or only a small flux outside 16$^\prime $ radius. However, for safety we do not consider here the spectra of the outermost annulus 9 (12$^\prime $-16$^\prime $), although the counts in this region have been taken into account in the deprojection for the inner shells.

Finally, we note that since all correction factors (vignetting, exposure, psf, projection) are taken into account in the deprojected spectra, we can use the nominal, on-axis response matrices for the spectral analysis. For these we took the v9q20 precanned matrices from the XMM-Newton calibration database (version April/May 2002).

2.6 Projected spectra

Deprojecting spectra is absolutely necessary for studying the innermost, cooling parts of the cluster. For proper modeling, any local hot gas in the core must be separated from projected hot gas from the outer parts. This is evident if one considers for example the strong abundance gradients observed in several clusters (Kaastra et al. 2001; Tamura et al. 2001a,b). However, the deprojection increases slightly the noise in the spectra, because also counts from the outermost annuli are taken into account in arriving at the deprojected spectrum. For studies of the outer parts of the cluster (for example the mass profile at large radii) deprojection is less urgent due to the slow variations of the spectral parameters with radius in those regions.

In the present work we do not use the non-deprojected spectra, however in a subsequent paper (Tamura et al. 2003, on abundances in clusters) we use extensively the non-deprojected spectra; since in that paper we will not repeat extensively all these analysis issues, we give them here for completeness.

When we use these non-deprojected spectra, we follow a more straightforward approach. We start with the template image (a Rosat PSPC image), apply the vignetting correction, convolve it explicitly with the instrumental psf and apply the exposure corrections. We then determine for each annulus the ratio of the number of PSPC counts in this corrected image to the number of PSPC counts in the original template image, and this determines the relative effective area correction (arf) that we apply to the standard response matrix. This is essentially the same method as used in the analysis of BeppoSAX MECS data of clusters by D'Acri et al. (1998).

2.7 Spectral binning

The spectra for each annulus were accumulated in 15 eV bins. These bins were then rebinned further onto a grid with a spacing of about 1/3 FWHM for all three detectors. We used the same energy grid for all three detectors. The data were binned then further by a factor between 2-12 as listed in Table 4. This was done in order to enhance the sensitivity in parts of the spectrum where there is only continuum or with very weak lines. It enhances the S/N ratio for the weak spectra of the outer parts. Spectral fitting was restricted to the 0.2-10 keV range. The energy range below 0.2 keV is currently too poorly calibrated to be useful for spectral analysis, while above 10 keV most of our spectra lack sufficient flux.

Table 4: Spectral binning. Number of bins of 1/3 FWHM taken together in intervals E1-E2, in the rest frame of the cluster.

For the spectral analysis, we used the SPEX package (Kaastra et al. 2002). For the interstellar absorption, we used the Morrison & McCammon (1983) model, and for the thermal plasma emission we used the collisional ionisation model (CIE) as available in SPEX. Systematic errors, both as a fraction of the source spectrum and the background spectrum, were added according to the prescriptions of the next section. The spectra of the three camera's were fit simultaneously, but in our figures we added them for clarity.

2.8 Systematic uncertainties

There are remaining systematic uncertainties in both the effective area to be used in the spectral response matrices, as well as in the subtracted background level. Based upon a study of several calibration sources by us, we have put conservatively systematic errors to the spectrum of 10% of the source flux below 0.3 keV and above 2 keV, and 5% in the 0.3-2 keV range.

In addition we included a systematic background error of 35% of the total background below 0.5 keV, 25% between 0.5-0.7 keV, 15% between 0.7-2 keV, and 10% above 10 keV. This was based upon our discussion in Sect. 2.3. Both systematic errors were added in quadrature to the binned spectra.

Finally, there is still a systematic difference in normalisation between the three EPIC camera's. Therefore in our spectral fits we allowed the relative normalisation of MOS2 and pn with respect to MOS1 to vary. In our fits, we found no systematic dependency of these renormalisations as a function of radius. For our full cluster sample, the average MOS2/MOS1 renormalisation factor was $1.011\pm 0.014$, while the pn/MOS1 renormalisation factor was $0.898\pm 0.013$. All spectral normalisations in the present work are related to MOS1.

  
3 Spectral modeling

3.1 Spectral model

There are several ways to model the spectrum of cooling gas in clusters of galaxies. These include multi-temperature fitting, parameter fitting of isobaric cooling flow models, both with and without the inclusion of an additional hot gas component, etc. The key ingredient in all these models is a prescription of the differential emission measure distribution ${\rm d}Y(T)/{\rm
d}T$, where the emission measure Y is the volume integral of the product $n_{\rm e}n_{\rm H}$, with $n_{\rm e}$ the electron density and $n_{\rm H}$ the hydrogen density.

In the present paper we attempt to do the spectral fitting as much as possible independently of any assumed physical model for ${\rm d}Y(T)/{\rm
d}T$. Based upon our spectral fits we make a numerical representation of ${\rm d}Y(T)/{\rm
d}T$. This numerical representation can be compared afterwards with any model for the cooling gas. Therefore, a kind of multi-temperature fitting is the preferred way to proceed.

In practice, the best temperature diagnostics for clusters of galaxies are obtained from the exponential turnover of the Bremsstrahlung continuum, in particular for the hottest gas, and the line emission from mainly the Fe-L and Fe-K complexes. In particular the Fe-L complex is a good temperature indicator for the cooler gas in the cluster. The line power emitted by any ion peaks at a characteristic temperature, and has a typical FWHM of the order of a factor of two in temperature. At each temperature, there are usually only a few ions with large emissivity. This combined with the factor of two FWHM for each ion makes it useless to sample ${\rm d}Y(T)/{\rm
d}T$ on a temperature grid that is finer than $\sim $0.15 in $\log T$. We performed simulations and found that a temperature binning of 0.15 in $\log T$ is sufficient to reconstruct any cluster spectrum, as compared to an arbitrary continuous or discrete emission measure distribution. However, with a spacing of 0.15 in $\log T$ we sometimes found strong correlations between the emission measure values of neighbouring bins. The reason is that the isothermal spectrum of a given temperature component can be approximated by a linear combination of the spectra of its neighbouring temperature bins. When the observed spectrum is noisy, it is hard to see the difference, and as a result strong correlations between the components arise. We found that with a spacing of 0.30 in $\log T$ these correlations disappear, while the temperature resolution is still sufficient to reconstruct ${\rm d}Y(T)/{\rm
d}T$ completely.

Therefore we have used a multi-temperature CIE model with temperatures T1, T2, etc., emission measures Y1, Y2, etc., constrained by the condition that Ti+1 = 0.5 Ti. We leave the temperature of the hottest component, T1, a free parameter, in order not to distribute the emission measure over two bins in case the spectrum is isothermal. This multi-temperature CIE model was redshifted with the cosmological redshift of each cluster and absorbed by Galactic absorption. Redshifts and column densities are given in Table 1. We have not found evidence for significant absorption above the Galactic values; this is confirmed by the higher resolution RGS data for the same clusters (Peterson et al. 2003); but even if there were additional absorption, this would not affect our multi-temperature fitting too much, since the strongest effect of absorption is the depth of the oxygen edge, and this edge energy is well below the energy of most of the strong Fe-L lines.

In the present paper we focus upon the temperature structure; as we explained the Fe-L complex is a good temperature indicator. Therefore we left the abundance of Fe free. Since K-shell lines from Ne and Mg blend with Fe-L lines, we also left the Ne and Mg abundance free. Both O and Si produce relatively strong spectral lines, and therefore we also left their abundances free. The other elements produce weaker lines that do not influence the emission measure estimate; since we are not interested in abundances in this paper, we couple the abundances of C and N to O, and those of Na, Al, S, Ar, Ca and Ni to Si. We assumed that all CIE components have the same abundance. The abundances are studied in greater detail by Tamura et al. (2003).

In order to keep the fitting procedure stable and reproducible, we started each time with a single temperature component with fixed temperature and abundances. In subsequent steps the number of free parameters is increased, and each time a spectral fit is made, starting with the most sensitive parameters. These parameters are in order: the temperature of the hot component, its iron abundance, its other abundances, the renormalisation between the different instruments. Then one by one the cooler temperature components were added. As a last step, the error bars for the relevant parameters were determined.

3.2 Single-temperature fitting

To facilitate comparison with previous work, we first present the results for single temperature fitting (Table 5). In general, we give temperatures and densities for the innermost 8 shells, except when the data become too noisy in the outermost shell, or in the innermost shell whenever there is a known strong X-ray emitting AGN in the center (Virgo, Perseus) or the data are noisy due to low surface brightness in combination with small extraction radius (A 3266). In Coma, we omit the innermost shell since it is strongly contaminated by NGC 4874, which we adopted as the center of our annuli. Since we study cooling flows in this paper, we also restrict the radial range to cooling times less than 100 Gyr.

The temperature profiles that we derive for the single-temperature fits are in good agreement with what we or others derived before from measurements with other satellites or from the same XMM-Newton data (M 87: Böhringer et al. 2001; Matsushita et al. 2002; Sérsic 159-3: Kaastra et al. 2001; A 496: Tamura et al. 2001b; A 1795: Tamura et al. 2001a; Arnaud et al. 2001; A 1835: Majerowicz et al. 2002).

Although the formal statistical error bars on the temperature and density are often very small, this does not prove that the spectra at each radius are strictly isothermal. The small error bars merely indicate the high statistical quality of our spectra. In fact, in particular in the inner parts of many clusters the $\chi^2$ of the fits are statistically not acceptable, contrary to the multi-temperature fits of Sect. 3.3 below.

Table 5: Spectral fit results for a single-temperature model. The # indicates the shell number (cf. Table 3).

  
3.3 Multi-temperature fitting

We list our results obtained using the multi-temperature fitting described in the previous section in Tables 67. We only list parameters for shells with a cooling time less than 30 Gyr, as determined from our single-temperature fits of the previous section. For most clusters the density outside that radius is low and as a result the X-ray flux is low, resulting in a spectrum that is often too noisy for multi-temperature fitting.

We give the root-mean-square error bars on all emission measures ( $\Delta\chi^2=2$). We averaged the positive and negative error bars: in most cases the differences between these errors is not very large. In those cases where the best-fit emission measure is zero, we only list the upper error bar: our model does not allow for negative emission measures for obvious reasons.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f3.ps}\end{figure} Figure 3: Differential emission measure distribution for A 262 in the innermost four shells. Line width decreases with radius; symbols: triangles: 0-0.5$^\prime $; squares: 0.5-1$^\prime $; stars: 1-2$^\prime $; circles: 2-3$^\prime $. The dashed line approximates the slope of the isobaric cooling flow model.
Open with DEXTER

In Fig. 3 we plot as a typical example the differential emission measure distribution for A 262 in the innermost four shells (where the radiative cooling time is shorter than the Hubble time). It is seen that the spectra are not strictly isothermal: there is a significant emission measure contribution at T2=0.5T1, although the emission measure drops rapidly towards lower temperatures. Note also that T1 increases with radius. Finally, the lowest temperature in an outer shell extends lower than the highest temperature in an inner shell. So we are not just seeing a temperature gradient smoothed in a given shell. The properties described in this paragraph for A 262 are typical for all the cooling clusters in our sample.

Table 6: Spectral fit results for the multi-temperature model. The temperature is the temperature of the hottest component, T1; the emission measures Y1 to Y5 have been renormalised to 1 for Y1. The # indicates the shell number. A "C" in the last column indicates consistency with an isobaric cooling flow spectrum, an "I" consistency with an isothermal spectrum.

Table 7: Spectral fit results for the multi-temperature model (continued). Same notation as Table 6.

In most cases (except for the innermost annuli of 2A 0335+096 and A 2052) two temperature components are sufficient to describe the spectrum, i.e. in those cases Y3 is within its error bars consistent with zero. This does not necessarily mean that Y3=0 in a strict mathematical sense, the error bars are such that in most cases we cannot exclude the presence of emission at T=0.25 T1 at a level below 1% of the emission measure Y1 of the hottest gas.

We have also investigated whether these emission measure distributions are consistent with the (cut-off) isobaric cooling flow model. The differential emission measure distribution ${\rm d}Y(T)/{\rm
d}T$ for the isobaric cooling flow model can be written as (see for example Fabian 1994):

 \begin{displaymath}
D(T)\equiv
{\rm d}Y(T)/{\rm d}T = {5\dot{M}k \over 2\mu m_{\rm H} \Lambda(T)},
\end{displaymath} (8)

where $\dot{M}$ is the mass deposition rate, k is Boltzmann's constant, $\mu$the mean molecular weight (0.618 for a plasma with 0.5 times solar abundances), $m_{\rm H}$ is the mass of a hydrogen atom, and $\Lambda(T)$ is the cooling function. The function D(T) for a 0.5 times solar abundances plasma is shown in Fig. 4. For any given cluster of galaxies and shell, the function D(T) should be put to zero naturally above a temperature $T_{\max}$, corresponding to the hottest gas at that location. It is seen from Fig. 4 that for the relevant part of D(T) for most of our data (between 0.3-3 keV) this function does not vary by more than 20% around a mean value of $1.3\times 10^{70}$ m $^{-3}~M_{\hbox{$\odot$ }}^{-1}$ yr. Hence, since in our modeling we have taken temperature steps of a factor of 2, the emission measures Yi for each temperature component Ti should scale approximately as $Y_i\sim T_i$. In several cases the deduced emission measure distribution decreases more steeply towards lower temperatures than proportional to T(Fig. 3), or may be even cut-off. Therefore we modeled the emission measures of Tables 67 with the function D(T), but cut-off to zero outside the interval ( $T_{\min},T_{\max}$), where $T_{\min}$ is a cut-off temperature. For the standard isobaric cooling flow model, $T_{\min}$ is zero.
  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f4.ps}\end{figure} Figure 4: Differential emission measure distribution for the isobaric cooling flow model, for 0.5 times solar abundances.
Open with DEXTER

We cannot assume that the division boundaries between the temperature components are exactly spaced by factors of two in temperature, for an arbitrary emission measure model. The division boundaries should be chosen in such a way that the emission measure weighted temperature for each bin equals the temperature assigned to it in our spectral fitting procedure. We adopted an iterative procedure where for a given set of parameters ( $T_{\min},T_{\max}$) and the emission measure distribution D(T) the bin boundaries were adjusted until the correct weighted average temperatures were determined. Then in a least squares fit the normalisation $\dot{M}$ was determined.

In our grid search we have put an artificial lower limit of 0.1 keV for $T_{\min}$; XMM-Newton is not sensitive to lower temperatures.

In a number of cases we find that the spectrum is consistent with $T_{\min}=T_{\max}$, i.e. the spectrum is formally indistinguishable from an isothermal plasma. This is the case for truly isothermal plasmas or spectra with poorer statistics, both conditions occurring in general in the outer parts of clusters. We indicate those cases with the letter I in the last column of Tables 67.

In other cases we find that $T_{\min}$ is less than 0.1 keV, hence the spectrum is formally consistent with the standard (non-cut off) isobaric cooling flow model, as far as the sensitivity range of XMM-Newton is concerned. Again, this is the case for truly isobaric cooling plasmas, or for spectra with poorer statistics. We indicate those cases with the letter C in the last column of Tables 67.

There are several instances where the statistics do not allow us to distinguish the full isobaric cooling flow model (C) from the purely isothermal model (I). These are the shells for which we list both a C and I in the last column of Table 6. In all other cases we have either a good lower limit to $T_{\min}$ or upper limit to $T_{\max}$.

Finally, in a few rare cases (Virgo shells 6-8, Perseus shells 2-6) the emission measure for T5 has formally a significant non-zero value, despite the fact that the table indicates that the spectrum is consistent with an isothermal model ("I"). This interpretation is chosen due to the fact that Y2-Y4 are zero within the error bars, and therefore the best fit to the isobaric cooling flow model finds $T_{\min}=T_{\max}$. Note however that the flux associated with T5 is very small, for example <0.4% of the 0.2-10 keV flux of shell 6 for Virgo, due to both the relatively small emission measure and the low temperature of $\sim $0.2 keV. These emission measures are well within the systematic uncertainties and we are only able to "detect" them thanks to the high statistical quality of the Virgo spectrum. In yet another case (Sérsic 159-3) the enhancement of Y5 for shells 2-4 despite the negligible values for Y3 and Y4 is due to the strong soft X-ray excess in this cluster (see Kaastra et al. 2003).

  
3.4 Alternative spectral models

The observed differential emission measure distribution in several of our clusters shows a rather steep temperature gradient downwards from the maximum temperature, as our multi-temperature fitting shows. Here we try to parameterise this behaviour, by fitting the spectra with a DEM model of the following shape:

 \begin{displaymath}
{ {\rm d}Y \over {\rm d}T } = \left\{
\begin{array}{ll}
c T^...
...
0 & \qquad\mbox{if $T \geq T_{\max}$ }.
\end{array} \right.
\end{displaymath} (9)

For $\alpha\rightarrow\infty$, we obtain the isothermal model, for large $\alpha $ a steep temperature decline is recovered while for $\alpha=0$ the emission measure distribution is flat. Note that Peterson et al. (2003) use a similar parameterisation but then for the differential luminosity distribution. In practice, we have implemented the model (9) by using the integrated emission measure $Y_{{\rm tot}}$instead of c for the normalisation, and instead of $\alpha $ its inverse $p=1/\alpha$, so that we can test isothermality by taking p=0. The emission measure distribution for the model is binned to bins with logarithmic steps of 0.10 in $\log T$, and for each bin the spectrum is evaluated at the emission measure averaged temperature and with the integrated emission measure for the relevant bin (this is needed since for large $\alpha $ the emission measure weighted temperature is very close to the upper temperature limit of the bin, and not to the bin centroid). This model is now implemented as the wdem model in the SPEX package. The results of our spectral fits are displayed in Table 9.

While multi-temperature fitting (Sect. 3.3) has the advantage of being completely model independent, the present wdem model has the advantage that the estimated uncertainties on the derived parameters, in particular $T_{\max}$ and $\alpha $ (and hence the degree of non-isothermality) are much smaller than the equivalent parameters for the multi-temperature model, provided of course that the wdem model applies. This is due to the fact that if the spectrum shows weak emission above $T_{\max}$ (either true or due to noise), the multi-temperature model may increase $T_{\max}$ to accommodate for this high temperature tail by adjusting all values of the emission measures accordingly, yielding almost the same best-fit $\chi^2$, but with larger error bars on the derived emission measures. A good example is the second shell of MKW 3s, where the multi-temperature model gives $\chi^2=242$ for 233 degrees of freedom, with $T_{\max}$ in the range between 3.8-10 keV , and with a 100% error on the emission measure of the hottest component. The wdem model has $\chi^2=244$ for 237 degrees of freedom (i.e. almost the same goodness of fit), but $T_{\max}$ is now well constrained: 4.34$~\pm~$0.38 keV, and with $1/\alpha=0.36^{+0.13}_{-0.17}$, the spectrum is non-isothermal at the 2$\sigma$confidence level. Therefore the wdem model offers a more sensitive test for isothermality than the multi-temperature fitting, and this is why in several cases the wdem model is able to detect the emission measure of the cooler gas and correspondingly finds a significant $1/\alpha>0$ (i.e., non-isothermality), despite the fact that the multi-temperature fit is consistent with an isothermal model, as indicated by an "I" in Tables 67. Note however that the emission measure distribution from the multi-temperature fits in the cases indicated above is not only consistent with an isothermal spectrum but in general also with a broad range of non-isothermal models, including the solution from the wdem model.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f5.ps}\end{figure} Figure 5: Spectrum of A 2052 in the inner three shells (with cooling gas). The spectra of all three EPIC camera's have been added. The spectra of the 0.5-1.0$^\prime $ and 1-2$^\prime $ shells have been multiplied by factors of 5 and 25, respectively. The spectra are shown as energy times counts/s/keV. Thin solid histograms: fits with a single temperature model; thick solid histograms: fits with the wdem model described in the text.
Open with DEXTER

Figure 5 shows an example of the fit improvement with the wdem model as compared to a single temperature model (the best fit model for the multi temperature model of Sect. 3.3 is very similar to the fit with the wdem model). The single temperature model in the innermost shells fails in two ways. First, it under-estimates the Bremsstrahlung continuum at high energies. This is due to the fact that the statistical quality of the spectrum near the Fe-L complex is much better than at high energies, due to a larger number of photons and a higher effective area. Therefore, the temperature of the hottest gas is dominated by the line emission from the most highly ionised ions in the Fe-L complex. Secondly, the single temperature model lacks the emission from the lowly ionised Fe-L ions (with predominantly line emission in the 0.7-0.9 keV range). This emission is included both in the multi-temperature fitting and in the wdem model. For larger radii, the spectrum becomes close to isothermal and there all three spectral models agree (see the 1-2$^\prime $ shell in Fig. 5).


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f6.ps}\end{figure} Figure 6: Steepness of the differential emission measure distribution versus cooling time for four clusters as indicated on the plot.
Open with DEXTER

In Fig. 6 we plot $1/\alpha $ versus cooling time for some of our clusters. For cooling times larger than 10-20 Gyr, the error bars on $1/\alpha $ are in general large, and include the value 0, which indicates that the spectra are consistent with an isothermal model in the outer parts of the cluster. The slightly enhanced values of $1/\alpha $ in a few rare cases at these locations are probably due to azimuthal temperature variations or other irregularities at larger distance from the core of the cluster. In the cooling region (cooling times less than 10-20 Gyr), none of our clusters shows a significant variation of $\alpha $ versus cooling time (and hence radius). Only in three cases (NGC 533, A 262 and A 2052) the innermost shell shows an enhanced value of $1/\alpha $ (see Fig. 6 for the case of A 2052). Perhaps in these cases an additional soft spectral component related to the central galaxy is present.

The average value for $1/\alpha $ for our clusters is listed in Table 10. For the three clusters mentioned above we omitted the central shell from the average. We also list the temperature of the hot gas as determined in the region where the cooling time is between 10-100 Gyr. These temperatures are in general somewhat higher than the temperatures derived from the single temperature fits in the same region. This is due to the fact that in the wdem model used here, a certain range of temperatures is allowed and hence the maximum temperature is higher than the average temperature. Figure 7 shows the data graphically.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f7.ps}\end{figure} Figure 7: Steepness of the differential emission measure distribution versus maximum temperature. The diamonds represent the data points with error bars. The dashed line indicates the weighted average value of $1/\alpha $, the solid line the best fit power law scaling as described in the text.
Open with DEXTER

The weighted average value for $1/\alpha $ for all clusters is 0.32$~\pm~$0.02, implying $\alpha=3.1\pm 0.2$. Thus, at each radius in the cooling region the emission measure distribution decreases rapidly, proportional to T3.1. There is a weak indication that $1/\alpha $ is slightly larger for the hotter clusters, see Fig. 7. Fitting formally a power law yields $\alpha=(7.1\pm 1.2) T^{-0.61~\pm~ 0.13}$, where T is in keV. However, this fit is driven by only a few cool or hot clusters, so more data is needed to confirm such a correlation.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f8.ps}\end{figure} Figure 8: Scaled temperature profile versus cooling time. Temperatures have been normalised to the temperature $T_{\rm h}$ of the hot gas outside the cooling region. Solid lines, triangles: clusters cooler than 3 keV; dashed lines, squares: clusters with temperatures from 3-4 keV; dash-dotted lines, stars: clusters with temperatures from 4-6 keV; dotted lines, circles: clusters with temperatures above 6 keV. The non-cooling clusters Coma, A754 and A 3266 have been omitted from the plot.
Open with DEXTER

  
4 Discussion

4.1 Present results

The models with a single temperature at each radius (Table 5) do not always provide the best spectral fits to the inner cooling parts in our clusters, although in some cases the single temperature model gives a $\chi^2$value that is formally acceptable. However, in almost all cases, the fits with the multi-temperature model (Tables 67) or wdem model (Table 9) have much lower $\chi^2$ values, and often show a significant deviation from isothermality (as expressed for example in the ratio Y2/Y1 and $\alpha $, respectively), although the emission measure distribution is dominated by emission from the hottest gas at each radius. As a typical example, the second shell of A 262 has $\chi^2=241$ for 237 degrees of freedom when fitted with the single temperature model, which is formally a good fit. But the multi-temperature model decreases $\chi^2$ to 220 (for 233 degrees of freedom), and gives a 3$\sigma$ detection of a cooler component ( $Y_2/Y_1=0.06\pm 0.02$).

However the single temperature fits show that gas cooling occurs in the center of these clusters. In all clusters where the cooling time at any radius is less than 10-15 Gyr we see a statistical significant temperature drop within that radius. Consistent with this, the non-cooling clusters A 3266, Coma and A 754 (cooling time always larger than 20 Gyr) show no central temperature drop.

  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f9.ps}\end{figure} Figure 9: Scaled temperature profile versus scaled radius. Temperatures have been normalised as indicated using the temperature $T_{\rm h}$ of the hot gas outside the cooling region as well as the central temperature $T_{\rm c}$. Triangles: clusters cooler than 3 keV; squares: clusters with temperatures from 3-4 keV; stars: clusters with temperatures from 4-6 keV; circles: clusters with temperatures above 6 keV. The non-cooling clusters Coma, A754 and A 3266 have been omitted from the plot. The solid line is the scaling curve Eq. (10) with $\mu =2$.
Open with DEXTER

However, there is no universal scaling law with cooling time (or gas density). This is illustrated in Fig. 8. Despite short cooling times, in Sérsic 159-3 and MKW 3s there is only a very weak temperature gradient. Contrary to this, in A 496 the temperature gradient is very steep.

We have also scaled our temperature profiles using the following expression:

 \begin{displaymath}
T(r) = T_{\rm c} + (T_{\rm h}-T_{\rm c}) \frac{(r/r_{\rm c})^{\mu}}
{1+(r/r_{\rm c})^{\mu}},
\end{displaymath} (10)

which is a reasonable universal fit to the temperature profile of cooling clusters (Allen et al. 2001). In addition, the choice of this profile is motivated by the analysis of Voigt et al. (2002). In order to reduce the number of parameters we set $T_{\rm c}$ equal to the temperature of the central bin and use $\mu =2$, as Allen et al. (2001) find $\mu=1.9 \pm 0.4$ when fitting the scaled temperature profiles of 6 cooling clusters.

The scaled temperature profile are shown in Fig. 9. Data outside a radius of 5$r_{\rm c}$ are not shown, as for some clusters the temperature starts decreasing again, therefore violating Eq. (10). All clusters obey the scaling curve (10) remarkably well. The cooling radius that we find correlates well with $T_{\rm h}$ ( $r_{\rm c} = (4.1\pm 0.8) T_{\rm h}^{1.84\pm 0.14}$, with $r_{\rm c}$in kpc and $T_{\rm h}$ in keV). Since the virial radius scales also with $T_{\rm h}^{0.5}$, the ratio of $r_{\rm c}$ to the virial radius is not constant, contrary to the findings of Allen et al. (2001). Note however that Allen et al. have a smaller sample with a more limited temperature range (hot, distant clusters).

We have investigated correlations of the cooling gas with other parameters. In Table 11 we list the average temperature $T_{\rm h}$ of the hot gas outside the cooling radius, the relative temperature decrement $(T_{\rm h}-T_{\rm c})/T_{\rm h}$, as well as the cooling radius $r_{{\rm cool}}$, defined here as the radius where the cooling time $t_{{\rm cool}}$ equals 15 Gyr. Over a wide range in cooling time (typically for cooling times between 5-50 Gyr), the cooling radius is found to scale as $r_{{\rm cool}}\sim t_{{\rm cool}}^a$. We also list the best-fit parameter a in Table 11.

We have searched for correlations of the relative temperature decrement $(T_{\rm h}-T_{\rm c})/T_{\rm h}$ with various quantities such as the cooling radius $r_{{\rm cool}}$, the hot gas temperature $T_{\rm h}$, the average temperature gradient $T_{\rm h}/r_{{\rm cool}}$, the central cooling time (hence central density), but none of them are significant. There is only a weak anti-correlation between cooling radius and relative temperature decrement, in the sense that the larger cooling regions tend to have relatively weaker cooling flows. However, this is only significant at the 10% significance level, and is perhaps a selection effect in the sense that the larger clusters tend to have lower densities hence longer cooling times. Another explanation is that the larger clusters have a relatively inefficient heat conduction (proportional to ${\rm d}T/{\rm d}r$). Therefore it is more difficult to balance radiative losses by heat conduction in the more extended clusters. The role of heat conduction is investigated in more detail in Sect. 4.4.

We note that most clusters have a rather narrow range of temperature decrement, between 0.35-0.60. The only exceptions are A 399 and A 2052, both with a high temperature decrement around 0.72, and Sérsic 159-3, MKW 3s and Hydra A, all three with a low temperature decrement of about 0.2. These temperature decrements are in excellent agreement with the results from the high-resolution RGS spectra of these clusters (Peterson et al. 2003).

The most important finding of our paper is however that the spectra of the cooling region are not isothermal. They are not isothermal in two aspects: first, the maximum temperature at each radius decreases towards the center; secondly, at each radius the gas has a range of temperatures, the emission measure distribution at each radius showing a steep decline or cut-off towards lower temperatures.

Table 8: Constraints for the isobaric cooling flow model. The lower limit, best-fit value and upper limit are given for both $T_{\min}$ (in keV) and $T_{\min}/T_{\max}$. Only shells with a cooling time less than 30 Gyr (as determined from the single temperature fits) are listed. No entries are listed whenever the data cannot constrain the ratio $T_{\min}/T_{\max}$ (i.e., any value between 0-1 is statistically allowed). The # indicates the shell number.

Table 9: Spectral fit results with the wdem model.

The decline of the maximum temperature of the gas towards the center of the cluster in all our cooling clusters (Tables 67) excludes a class of models where the core consists of a mixture of hot isothermal gas and significantly cooler gas, with a position-dependent filling factor of the cool gas (for example Fukazawa et al. 1994; Ikebe et al. 1999).

The fact that even at a single radius the spectrum is not isothermal contradicts any model for a spherical symmetric, single phase cooling flow model. We will elaborate on this last model in the next section. Here we note that the non-isothermality is not caused by the width of our extraction annuli, because the temperature gradients across our annuli are not large enough for that.

Most of our clusters have a small but significant emission measure down to about 0.5$T_{\max}$ (Tables 67). Our fits to this differential emission measure distribution with a cut-off isobaric cooling flow model (Sect. 3.3, Table 8) show a range of values for the ratio $T_{\min}/T_{\max}$ for the cases with good statistics. As we show in this paper, the isobaric cooling flow model does not apply to our clusters but nevertheless the shape of its emission measure distribution is a convenient parameterisation of the shape of the actual emission measure distribution. This temperature ratio $T_{\min}/T_{\max}$ varies between 0.7-0.9 for a cluster like Hydra A, 0.46 for Coma and A 262, down to 0.30 for 2A 0335+096 and 0.25 for Perseus. Alternatively, the minimum temperature $T_{\min}$ can be as low as 0.7 keV for NGC 533 and A 262 but at least twice as large for other clusters.

There is no clear correlation between this cut-off temperature (both in absolute and relative terms) and other cluster parameters. This excludes for example inhomogeneous metallicity models, as shown already by Peterson et al. (2003) from the RGS data of these clusters.

Our fits to the cut-off isobaric cooling flow model have relatively large error bars. This is due to the fact that our results are based upon a two-stage fitting procedure: first determine the differential emission measure distribution, and from that derive the best-fit cut-off isobaric cooling flow parameters. In particular the temperature of the hottest component T1sometimes has a rather large uncertainty despite the good statistical quality of the spectrum. This is due to the fact that by increasing T1 its emission measure Y1 should decrease, but because the cooler temperature components also shift their temperatures, a rearrangement of the emission measure distribution takes place resulting in almost the same spectrum and differential emission measure distribution, but as a result with relatively large error bars. This implies that our test for consistency with an isothermal model as indicated by an "I" in Tables 67 is not the most sensitive test.

With respect to that our fits with the emission measure distribution (9) is more robust. It shows that in general the emission measure distribution of the cooling gas has a rather constant shape, both as a function of radius as well as from cluster to cluster (Sect. 3.4). We come back to this in our discussion on magnetic loop models (Sect. 4.5).

It is impossible to discuss in this paper all models proposed to solve the cooling flow problem in relation to our current findings. In the next sections we discuss a few of the most promising scenario's. But we hope that the model independent emission measure distributions derived in this paper will be useful to test any alternative cooling flow model.

4.2 The isobaric cooling flow model

It has been argued that a simple stationary flow dominated by radiative cooling is inconsistent with the density profiles of clusters (Fabian et al. 1984; Thomas et al. 1987; Johnstone et al. 1992). If radiative cooling were the only important energetic process, we would expect that the emission measure distribution would resemble that of (8) at each radius. Clearly, the high resolution RGS observations are inconsistent with that distribution integrated over the entire cluster (Kaastra et al. 2001; Peterson et al. 2001; Tamura et al. 2001a,b; Peterson et al. 2003).

Detailed measurements of the temperature distribution at each radius have been relatively difficult until observations with XMM-Newton EPIC and Chandra. Many recent studies (David et al. 2001; Molendi & Pizzolato 2001; Johnstone et al. 2002; Matsushita et al. 2002; Schmidt et al. 2002; Sanders & Fabian 2002; Buote 2002) have argued to varying degrees that the temperature distribution at each radius is isothermal or very close to isothermal at each radius. Our measurements are in broad agreement with the average temperatures found by these studies. We argue, however, that a single-phase temperature distribution at each radius is unlikely to be the case as demonstrated by the better spectral fits with our multi-temperature or wdem models. Instead, either a relatively steep or narrow differential emission measure distribution might be required. The exact form of the differential distribution is difficult to determine at the current spectral resolution. It might be possible that the parameterization of the global temperature distribution offered in Peterson et al. (2003) partly reflects a temperature gradient as in (10) but also represents a steep intrinsic differential emission distribution like that modelled in Sect. 3.4.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f10.ps}\end{figure} Figure 10: Pressure versus temperature for 16 clusters. Only results for shells where the cooling time is less than 15 Gyr are shown. The center of the cluster is the leftmost data point for each cluster. Pressure and temperature are derived from the isothermal model fits (Table 5).
Open with DEXTER

Table 10: Weighted average temperature $T_{\max}$ of the hot gas in the region with cooling time between 10-100 Gyr, and weighted average value for $1/\alpha $ for the full cluster (average is dominated by the cooling cluster core).

Table 11: Characteristic temperatures and cooling times. $T_{\rm h}$ is the temperature outside the cooling region, $T_{\rm c}$ the temperature in the center of the cluster as determined from a single temperature fit; $r_{\rm c}$ is the characteristic radius of the temperature profile (see Eq. (10); $r_{{\rm cool}}$ is the cooling radius for a cooling time of 15 Gyr and a is the power law index of the scaling of $r_{{\rm cool}}$ with cooling time (see text for details).

In any case the present data rule out a model with a local isobaric cooling flow (8) at each radius for most of our clusters. The few cases where we cannot rule it out formally are those cases where the spectra are relatively noisy. This is evident from the steep emission measure distribution or cut-off emission measure distribution that we find in the present work.

Finally, our spectral modeling confirms that the cluster cores are not isobaric: for most of our clusters, the pressure drops by half an order of magnitude from the cluster center to the cooling radius $r_{{\rm cool}}$(Fig. 10).

4.3 Heating by AGN

Table 12: Radio fluxes of the central galaxies in clusters. S1.4 is the flux (mJy) at 1.4 GHz, $L_{\nu }$ the intrinsic luminosity at the same frequency, $\alpha $ is the radio spectral index ( $L_{\nu }\sim \nu ^{\alpha }$). $L_{\rm X}$ is the total 0.2-10 keV luminosity within the cooling radius $r_{{\rm cool}}$ as given in Table 11.

It has been proposed that heating by a central AGN could balance, at least for a part, the radiative losses in a cluster due to cooling (e.g. Böhringer et al. 2002). In order to test this hypothesis on our large sample, we have collected the radio fluxes of the central galaxies in each cluster (Table 12). The fluxes are all evaluated at 1.4 GHz (20 cm) and most of them were taken from the NRAO VLA Sky Survey (Condon et al. 1998) or other surveys available in the literature.

We have searched for a correlation of the relative temperature decrement $(T_{\rm h}-T_{\rm c})/T_{\rm h}$ with the radio luminosity of the central galaxy, but there is no significant correlation (Spearman rank correlation coefficient -0.28, but formally significant only at the 1.1$\sigma$ level). Only the most radio luminous cluster in our sample (Hydra A) is also the cluster with the smallest temperature decrement, despite the short central cooling time of 1.6 Gyr.


  \begin{figure}
\par\includegraphics[angle=-90,width=8.8cm,clip]{h4230f11.ps}\end{figure} Figure 11: Correlation between the 0.2-10 keV X-ray luminosity within the cooling radius $r_{{\rm cool}}$ as given in Table 11 and $\nu L(\nu )$ at 1.4 GHz for the central radio galaxy. Upper limits are indicated by an arrow.
Open with DEXTER

However there is a correlation between the radio luminosity of the central AGN with the X-ray luminosity within the cooling radius (Fig. 11). Interestingly, such a correlation was not found in a similar analysis based upon Chandra data (Voigt & Fabian 2003), however they used luminous cooling clusters (X-ray luminosity within the cooling radius >1037 W), while we also include less luminous cooling clusters. This correlation could be consistent with a heating by the central AGN scenario, but a correlation alone is insufficient proof for this. First, the correlation may be due to some selection effect or common "hidden" variable to which both $L_{\rm X}$ and $L_{{\rm AGN}}$ correlate. Furthermore, radio flux alone is not a proper measure of the time averaged power of an AGN. The relation between total power and observed radio luminosity is not perfect, due to intrinsic differences in radio loudness and orientation (beaming) effects, as well as the fact that radio and X-rays look at different time scales.

  
4.4 Heat conduction by electrons

Thermal conduction by electrons might play an important role in cooling clusters. It is possible that some heat influx from the hot outer region of the cluster to the center, totally or only partially balances radiative losses. We have shown in Sect. 3.3, however, that the ICM is not locally isothermal at each radius but has a range of temperatures with a steep differential emission measure. Although the cause of this particular emission measure distribution is unknown, it is unlikely that thermal conduction would be effective on large (Mpc) scales from the outer regions of clusters if conduction is incapable in removing temperature inhomogeneities on smaller (kpc) scales. This seems unlikely since as the conduction coefficient is increased, thermal instability is removed first at the largest wave numbers (Field 1965).

We can also derive conduction coefficients ignoring the above argument and assuming that heat conduction by electrons balances radiative losses, i.e.:

 \begin{displaymath}
\int_{{\rm V}}n^2 \Lambda(T){\rm d}V=\int_{{\rm S}} \kappa (\nabla T){\rm d}S,
\end{displaymath} (11)

where n, $\Lambda$, T and $\kappa $ are the gas density, cooling function, temperature and thermal conductivity, respectively and V and S are the volume and surface area of the X-ray emitting region, respectively. Assuming spherical symmetry and using data binned into shells, the conductivity coefficients at the boundary of each shell can be computed. Using this procedure, Voigt et al. (2002) investigated thermal conduction in two cooling clusters: A 2199 and A 1835. We extent this analysis to our sample, including the three non cooling clusters as a control case. The temperature gradient between shells is evaluated by fitting the temperature profile with the function (10). We exclude A 1837 from the analysis because its temperature profile cannot be fitted with Eq. (10), since only two bins are within the cooling radius.


  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{h4230f12.ps}\end{figure} Figure 12: The conduction coefficients $\kappa $ required for heat conduction to balance radiation losses as a function of temperature. The solid line is the Spitzer conductivity, the dashed lines are one fifth and one hundredth of the Spitzer conductivity, respectively. Filled diamonds represent the values for the three clusters without cooling flow (Coma, A 3266 and A 754), while open diamonds represent the values for the clusters with a shallow temperature profile (MKW 3s, Sérsic 159-3 and Hydra A). Filled stars represent the remaining clusters. Values for the same cluster are joined by a line and values are given only for bins with cooling times less than $\sim $30 Gyr.
Open with DEXTER

Our estimates of the conduction coefficients $\kappa $ are shown in Fig. 12, in which we do not plot error bars in order to facilitate the reading (the 1$\sigma$ errors for $\kappa $ are $\sim $20-90% depending mainly on the accuracy of the temperature estimates).

The estimated conduction coefficients must be compared to theoretical calculations to see whether heat conduction from the outer regions can totally balance radiative losses. In this paper, we do not discuss the possibility of a partial balance of heat conduction and radiative losses. For a highly ionised plasma such as the ICM, the maximum rate is expected to be the Spitzer conductivity $\kappa_{{\rm s}}$ (Spitzer 1962). In the presence of a homogeneous magnetic field, the conductivity is the Spitzer rate only along the field, but severely decreased in the transverse direction. Thermal conduction in a tangled magnetic field has been studied by Chandran & Cowley (1998). They conclude that in cooling flows, thermal conductivity is below the Spitzer level by a factor of order 102 to 103, depending on the assumed field tangling scale. Conductivity is less severely decreased if the magnetic field behaves chaotically over a wide range of scales: Narayan & Medvedev (2001) estimate that in this case conductivity is only a factor $\sim $5 below the Spitzer rate (given the uncertainties this factor lies in the range $\sim $2.5-10). Therefore, the Spitzer conductivity $\kappa_{{\rm s}}$, $\kappa_{{\rm s}}/5$ and $\kappa_{{\rm s}}/100$ are shown in Fig. 12 for comparison.

For the three non cooling clusters we find very low values for the conductivity coefficients: the absence of significant cooling allows even inefficient heat conduction to remove temperature inhomogeneities. For the cooling clusters we see that clusters with shallow temperature profile show a very different trend and relatively high conductivity because of the small temperature gradients. In addition, in most of the clusters conductivity increases towards the center (lowest temperature) and as we approach the radius at which the temperature reaches its maximum (this can not be clearly seen in Fig. 12 for every cluster because for some the cooling time at this radius is larger than $\sim $30 Gyr and consequently not shown). This trend is due to the fact that our temperature profile model has vanishing gradient both at large and small radii. Most important is the fact that we find conductivities larger that one fifth of the Spitzer rate in almost all the cooling clusters. Since the thermal conductivity in the ICM is estimated to be below the Spitzer rate by at least a factor of 5, we conclude that heat conduction alone is insufficient to balance radiative losses in cooling clusters.

Although we demonstrated that heat conduction on a global scale cannot maintain the large scale temperature gradients in cooling clusters, heat conduction may still play an important role on smaller spatial scales. In the next section we consider a class of magnetic loop models as a possible mechanism in which heat conduction on smaller scales may occur.

  
4.5 Magnetic loop models

4.5.1 Cluster magnetic fields

It has been argued in the past by several groups that magnetic fields may play an important role in the dynamics of clusters. For example, Clarke et al. (2001) show in a study of rotation measures for a large sample of X-ray bright Abell clusters that the magnetic field occupies at least 95% of the area of the non-cooling innermost 750 kpc of these clusters. This suggests that enhanced magnetic field levels permeate the clusters with a high filling factor, with magnetic structure down to the level of 10 kpc or even smaller. The average magnetic field is of the order of 1 nT, and the corresponding magnetic energy density is a few percent of the thermal energy density of the gas. In the centers of cooling clusters, magnetic fields up to 4 nT are found (Taylor et al. 2002).

One of the most important aspects of magnetic fields in clusters of galaxies is that heat conduction in order to balance radiative losses can only occur along the magnetic field lines. Thus, the effectiveness of thermal conduction depends strongly upon the magnetic topology. For example, Norman & Meiksin (1996) have considered models where relatively cold loops reconnect with hot loops that are connected to the thermal reservoir of the cluster. This process would effectively slow down the cooling. A problem with this model is that hot loops must extend down from the outer parts of the cluster to the region where the cool loops reside. In that case, we would still expect - apart from the cooling gas - to see a hot component throughout the cluster volume with an almost constant temperature. In all cooling clusters we find however that the temperature of the hottest gas component decreases significantly towards the center of the cluster.

4.5.2 Coronal loop models

Here we consider a different class of magnetic loops, motivated by a recent paper by Aschwanden & Schrijver (2002, AS hereafter). These authors created an analytical model for solar and stellar coronal loops. The differential emission measure distribution of their models resembles the differential emission measure distributions that we obtained for our clusters: a distribution that is strongly peaked towards the maximum temperature $T_{\rm a}$at the top of the loops.

Their model considers stationary loops (i.e. no net gas flows) for which mass and momentum conservation is taken into account. In the energy balance equation, radiative losses, heat conduction and position-dependent heating are taken into account, as well as the effects of gravity. For the geometry of the loops they took a semicircular shape with height $h(s)=h_1\sin(s/h_1)$ where sis the coordinate along the loop. The loop half length is L and the loop crosses the photosphere at s=s0. The loop cross section a(s) expands by an expansion factor $\Gamma $ from the photosphere to the loop apex. The heating process is not prescribed in the AS models, it is only parameterised with a form $E_{\rm H}(s) = E_0 \exp{(-s/s_{\rm H})}$, with E0 the heating rate (W m-3) at the foot points of the loop and $s_{\rm H}$ a scale height. For $s_{\rm H}=0$, all heating takes place at the loop foot points, while for $s_{\rm H} \gg L$ the heating occurs uniformly across the loop. The hydrodynamic equations are solved with as boundary condition the requirement that the conductive flux at the foot points and top of the loop vanishes. The fact that ${\rm d}T/{\rm d}s=0$ at the apex of the loops causes the sharp peak of the differential emission measure just below the apex temperature $T_{\rm a}$. In general, the solutions have a temperature decrease and density increase from apex to foot points, with only a relatively small pressure increase. The loops also obey scaling laws, which are to lowest order similar to the scaling laws derived first by Rosner et al. (1978), namely $n_0 L\sim T_{\rm a}^2$, with n0 the gas density at the loop apex. At fixed $T_{\rm a}$, the loops have therefore similar column density ( ${\sim} n_0 L$), but the shorter loops have higher density and therefore are relatively very efficient X-ray emitters (emissivity scales with $n_0^2
\Lambda(T)$).

4.5.3 Application of the coronal loop models to clusters

We expect the same physical processes that play a role in coronal loops to play a role in magnetic loops in clusters, be it on a different scale. The centers of clusters contain highly tangled magnetic fields, hence heat conduction takes preferentially place along the magnetic field lines of the loop. Also the densities are high enough that radiative cooling occurs. Loops can be heated by a variety of processes, whatever their origin. In fact the lack of large amounts of cool gas in the cluster cores provides evidence for (effective) heating, either direct or indirect by mixing processes. The precise heating process is not important for our present discussion. Whatever the heating process in clusters is, the heat will predominantly be transported along the magnetic field lines and the heat will finally be lost due to the observed radiative cooling. These are exactly the same conditions that are the basis for the solar coronal loop models.

It is harder to identify the foot points of magnetic loops in clusters. These foot points can be associated with active galactic nuclei, member galaxies, or just enhancements in the highly tangled intracluster field.

Table 13: Emission measure distributions for the Aschwanden & Schrijver (2002) models. $T_{\rm a}$ (keV) is the maximum temperature at the loop apex. $\Gamma $ is the loop cross section ratio between the loop apex and the photosphere. L is the loop half length (kpc), and $s_{\rm H}$ the heating scale height. M13 is the cluster mass in $10^{13}~M_{\hbox{$\odot$ }}$ within the radius R (kpc) where the loop is located. Y1-Y5 are the emission measures Yi (in m-3 of the temperature components 1-5. The corresponding temperatures are Ti = T1 / 2i-1. The density at the loop apex is $n_{\rm a}$(given in units of 104 m-3) and the average density over the loop volume is $\overline {n}$. The mass of an individual loop is Ml (in $10^{10}~M_{\hbox{$\odot$ }}$). Finally, N is the maximum number of loops per scale height, defined by $N\equiv 4\pi R^3/V$ with V the volume of an individual loop. All calculations are done for loops with a cross-sectional radius at the foot points (r0) of 10 kpc. For other values of r0, the following scaling applies: $Y_i\sim r_0^2$, $M_l\sim r_0^2$ and $N\sim r_0^{-2}$.

We used the model of AS and applied it for typical cluster conditions. Typical median masses of $8\times 10^{12}~M_{\hbox{$\odot$ }}$ at 100 kpc and $4\times
10^{13}~M_{\hbox{$\odot$ }}$ at 200 kpc were used (derived from the work of Ettori et al. 2002). We took the differential emission measure distribution from the AS models, and binned it into bins with a factor of two difference in temperature, just as we did for our clusters. The highest temperature bin was centered around the emission measure weighted loop temperature, which is close below the maximum temperature $T_{\rm a}$ at the apex of the loop. Table 13 lists a few typical cases of the expected emission measure distribution.

The table shows that the shape of the emission measure distribution (as characterised for example by the ratio Y2/Y1) depends strongly upon the maximum temperature $T_{\rm a}$ as well as the loop expansion factor $\Gamma $: steep differential emission measure distributions are obtained for large $T_{\rm a}$ or large $\Gamma $. The loop length L is not important for the shape of the emission measure distribution as long as the loops are smaller than the scale height for the gravity of the cluster (set by R). The heating scale height $s_{\rm H}$ is only important when $s_{\rm H}$ becomes smaller than the loop length L, i.e. when the heating occurs preferentially near the foot points of the loop instead of uniformly over the volume. The results do not depend strongly upon the central mass M13 and the position in the cluster R (i.e. the gravity of the cluster is relatively unimportant, at least for sufficiently small loops). It is also seen that for a fixed temperature $T_{\rm a}$, different combinations of $\Gamma $ and L can produce the same emission measure distribution.

An important difference between various models yielding the same emission measure distribution shape is the absolute value of the emission measure as well as the total mass associated with the loops. For some models, the total mass of the loops assuming unity filling factor (NMl) exceeds the gravitating mass within that radius. This occurs when L or $\Gamma $ becomes too small. These unphysical models are to be rejected, hence in practice useful lower limits for $\Gamma $ or L can be obtained.

The ratio Y2/Y1 for the inner parts of our cluster sample (defined here as the part where the cooling time is less than 10 Gyr) has a median value of 0.11, with individual cases typically between 0.00 and 0.40. Virgo, which has one of the best measured spectra, has Y2/Y1 distributed between 0.05-0.13.

For large heating scale heights $s_{\rm H}/L=10$ (almost uniform heating as a function of loop radius), this ratio can be obtained for loop expansion factors $\Gamma $of the order 10, with a typical spread of half a decade. In this scenario, an increase of Y2/Y1 towards the center of a cluster, where the ambient temperature is lower, is expected if all other parameters (including $\Gamma $) remain constant. This is indeed observed in some of our clusters (NGC 533, A 262 and A 2052, see Sect. 3.4). If a uniform value of $\Gamma $for each cluster would apply, then we would expect that the hotter clusters have smaller Y2/Y1, i.e. are closer to isothermal. This is not consistent with the observations, which show a tendency for the hotter clusters to have smaller power law indices $\alpha $ (see Fig. 6). This then would imply that the hotter clusters should have either smaller loop expansion factors $\Gamma $, maybe even down to $\Gamma=1$, or larger loop length L. A larger loop length seems more likely since this keeps the total mass of the loops within acceptable limits.

Another possibility is that the heating scale height $s_{\rm H}$ is substantially smaller than the semi loop length L. In that case, the observed values for Y2/Y1 can be explained by a $\Gamma $ that is a few times smaller than in the case of uniform heating. The increase of Y2/Y1 for hotter clusters could also be explained if the hotter clusters have loops that are more uniformly heated.

Table 14: Application of the coronal loop models of Aschwanden & Schrijver (2002) to the innermost three shells of 2A 0335+096. M13 is the total cluster mass in $10^{13}~M_{\hbox{$\odot$ }}$ within the radius R (kpc), for which we take the average radius of the shell. $\overline {\rho _l}/\rho _{{\rm tot}}$ is the average mass density of the loops divided by the total mass density of the cluster (including dark matter etc.). The total mass density profile was taken from Ettori et al. (2002). In this table, f is the volume filling factor of the loops.

The number of loops can be quite large. As an example we have worked out here the innermost three shells of 2A 0335+096. The cooling time in these shells is less than the age of the universe, and all three shells have significant multi-temperature structure. From the emission measure distribution of Table 6 we obtain the loop parameters as listed in Table 14. For unity volume filling factor f, the loops have a typical length of 0.25-0.50 times the radius R where the loops are located. In the outermost shell this is in fair agreement with the magnetic field reversal length as obtained from numerical simulations of non-cooling regions (e.g., Dolag et al. 2002), in the innermost shell it is significantly smaller. Since the mass density of the loops compared to the total gravitating mass density ( $\overline {\rho _l}/\rho _{{\rm tot}}$) is relatively large for f=1, a smaller volume filling factor is more likely in the innermost two shells, giving even a much smaller loop length L. A smaller filling factor of the hot gas in the cores of several clusters is also indicated by observations of cavities in clusters with the Rosat HRI or Chandra (for example Perseus, Böhringer et al. 1993; Fabian et al. 2000). However, the number of loops needed is large, considering also that for smaller loop length the loop radius will be also smaller than our adopted 10 kpc. The large number of loops as well as their small size then imply a highly tangled magnetic field structure in the central part of the cluster. A large number of loops also causes a relatively smooth surface brightness, consistent with detailed X-ray mapping of several clusters.

Some caution is needed, however. In the context of the loop model, it is not clear how much emission is due to the magnetised loops and how much is due to emission from gas in between. If a substantial fraction of the cluster emission is from very hot gas in between the loops (for example in cavities), then in the measured ratio Y2/Y1 the emission measure Y1 should first be corrected for this hot gas. However, in that case there would still be the problem why we see so little cooling gas in our clusters.

Other points of caution are the different abundances in the cluster as compared to the solar corona, and hence a somewhat different cooling function, and the precise physical conditions at the foot points of the loops.

5 Conclusions

We have determined in our work the spatially resolved differential emission measure distribution of sample of 17 cooling clusters, plus three non-cooling clusters for comparison.

Most of the cooling clusters show a central temperature decrement of a factor of $\sim $2. A few clusters (Sérsic 159-3, MKW 3s and Hydra A) do not show such a strong temperature decrement. In the last case this may be caused by the luminous AGN at the center of the cluster. Two clusters have a very strong temperature decrement (A 399 and A 2052), the reason for this is not clear. We also find a correlation between the X-ray luminosity within the cooling radius with the radio luminosity of the central AGN.

It is clear from our modeling that the isobaric cooling flow model does not work in all our clusters in regions with cooling time less than about 15 Gyr. It does not work on a global scale, since there is a weak pressure gradient in the core and the spectra at each radius are not isothermal. It also does not work on a local scale, because there is a steep power law like emission measure distribution at each radius. This is much steeper than is expected from the isobaric cooling flow model. At each radius the differential emission measure distribution is approximately proportional to $T^{3.1\pm0.2}$.

This emission measure distribution looks remarkably similar to the predicted emission measure distribution of coronal loops. Applying coronal loop models to our emission measure distribution shows that moderate loop expansion factors of order 10 are required (i.e., the loop tops are about a factor of 3 larger than the foot points). Assuming a base loop radius of order 10 kpc (the typical size of a galaxy, or the characteristic size of the magnetic field as derived from Faraday rotation measurements), at least a few hundred loops per shell should be present.

Acknowledgements
This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). The Space Research Organisation of the Netherlands (SRON) is supported financially by NWO, The Netherlands Organisation for Scientific Research. MSSL is supported financially by the UK Particle Physics and Astronomy Research Council. We thank an anonymous referee for making valuable comments on the manuscript.

Appendix A: Determination of spectra and response matrices

We follow the usual method of accumulating spectra in concentric annuli around the cluster center. The accumulation of background-subtracted spectra is easy. Background subtraction is done using the same detector regions for source and background. It is more difficult to obtain the proper response matrices. They are created as follows. We start with a map of the intensity distribution $\mu(x,y)$ of the cluster on the sky. This intensity map is then multiplied by the effective vignetting v(x,y) of the telescope and detector (the off-axis relative effective area). This vignetted map is then convolved with the instrumental point-spread function (psf) that can be expressed as $p(x^\prime-x,y^\prime-y;x,y)$ giving the probability that a photon emitted at coordinates (x,y) is observed at position $(x^\prime,y^\prime)$. In our numerical estimates, we have taken this position dependence of the psf into account. The resulting map is then multiplied by the exposure map $e(x^\prime,y^\prime)$, which gives the fraction of time that each pixel was "on" (excluded hot pixels, gaps between CCD's and point sources are accounted for this way). The resulting predicted intensity distribution $f(x^\prime,y^\prime)$ is then given by

 \begin{displaymath}
f(x^\prime,y^\prime) = e(x^\prime,y^\prime) \int {\rm d}x
\int {\rm d}y
~\mu(x,y) v(x,y)
p(x^\prime-x,y^\prime-y;x,y).
\end{displaymath} (A.1)

The above formalism holds for any energy bin Ek. We divide the cluster into a set of annuli with outer radii ri. The emitted flux $\phi_i$ in annulus i is then given by

 \begin{displaymath}%
\phi_i \equiv
\int\limits_{{0}}^{{2\pi}} {\rm d}\theta
\int\limits_{{r_{i-1}}}^{{r_i}} r {\rm d}r~
\mu(x,y),
\end{displaymath} (A.2)

while the observed flux ci (counts) is given by

 \begin{displaymath}%
c_i \equiv
\int\limits_{{0}}^{{2\pi}} {\rm d}\theta
\int\limits_{{r_{i-1}}}^{{r_i}} r{\rm d}r~
f(x,y).
\end{displaymath} (A.3)

Combining (A.3) with (A.2) and (A.1), we formally write

 \begin{displaymath}
c_i \equiv \sum\limits_{j}^{} R_{ij} \phi_j
\end{displaymath} (A.4)

as the expected count rate in annulus i due to emission in all the annuli j. Again, (A.4) can be evaluated for any photon energy Ek. Because the matrix Rij is based upon a convolution of the source image $\mu$ with the instrumental psf p, the count spectrum in each annulus depends upon the emitted spectrum in all other annuli (i.e., the matrix Rij is non-diagonal).

Our initial spectral fitting method therefore consisted of making model spectra in each annulus, and fitting all count spectra, coupled through (A.4), simultaneously. We also used a variant where we considered a separate spectral model in each of a set of concentric, three-dimensional shells, projecting those on the sky and then convolving with the psf.

As an illustration we present in Table A.1 the predicted correlations between the annuli for the MOS1 data of A 1795. The magnitude of the redistribution due to the instrumental psf as well as the vignetting effects are clearly visible.

Table A.1: Distribution of the observed flux, expressed as a percentage of the original flux emitted in a given annulus. Calculation applies to the MOS1 data of A 1795. The observed flux is corrected for vignetting, psf and exposure map. The column o labels the observed annulus number, the row ethe emitting annulus.

Unfortunately this method appeared not to be very stable. The reason is the significant overlap of neighboring shells. If, due to statistical fluctuations, the best-fit normalisation of the spectrum in one shell happens to be large, this can be compensated by decreasing the normalisation of the neighbouring shells. These then affect their neighbours, etc. Apart from this effect, the fitting procedure and in particular the error search consumes considerable computing time, because the complex spectra of all shells need to be fitted simultaneously, leading quickly to 50-100 free parameters of the fit, depending on the sophistication of the spectral model. And due to the strong correlations between the parameters, error search takes a very long time.

Therefore we have developed a different approach. We deproject spectra directly and take account of the correction factors for the psf as discussed in Sect. 2.4. We outline our method in Sect. 2.5.

References



Copyright ESO 2003