text:
# The mass profile of A1413 observed with XMM-Newton: implications for the 𝑀–𝑇 relation.

## 1 Introduction

The simple first-order formation scenario for galaxy clusters, in which they grow through the gravitational infall and subsequent merging of smaller subunits, provides a remarkably good description of the large-scale properties of these objects. Within this hierarchical model, the gas trapped in the potential well of a cluster is heated to the observed X-ray emitting temperatures by the shocks due to the formation process; merger features in the gas distribution are then erased in roughly a sound crossing time ($``$ few Gyr), leaving the gas in hydrostatic equilibrium (HE).

Observation of this gas is a powerful tool for uncovering the physical characteristics and formation history of a cluster. Substructure in X-ray images, combined with optical data, can give clues to the dynamical state (e.g., Buote buote01 (2001)). Direct (temperature maps) and indirect (hardness ratio maps) methods can give an indication of where (and if) interactions and mergers are still occurring (e.g., Markevitch et al. mark99 (1999); Neumann et al. dmnetal01 (2001)). In addition, for clusters in reasonably relaxed state, the assumption of HE and spherical symmetry allow the derivation of the spatial distribution of both the gas and total cluster mass by using the information from the X-ray surface brightness and temperature profiles. This approach, which is of fundamental use in cluster studies, has been shown to give masses which are accurate to about $`\pm 20\%`$ when applied to simulated clusters (e.g., Evrard, Metzler & Navarro emn96 (1996) \[EMN96\]; Schindler schindler96 (1996)).

Numerical simulations based on gravitational collapse are an essential counterpoint to the observations, being as they are ideal scenarios with exactly measurable quantities, thus offering a direct comparison with the real data. A crucial result from these simulations is the suggestion that CDM haloes with masses spanning several orders of magnitude follow a universal density profile independent of halo mass or cosmology (Navarro, Frenk & White nfw97 (1997) \[NFW\]). As the X-ray emitting gas lies in the potential well of the CDM halo, this suggests that many directly measurable cluster properties should display self-similarity. This is observationally testable and indeed, regularity in the local cluster population has been found in previous ROSAT, ASCA and BeppoSAX studies, where the gas density and temperature profiles of hot, relaxed clusters do appear similar when scaled to units of the virial radius111Normally defined from numerical simulations as the radius of fixed density contrast $`\delta =200`$, or $`r_{200}`$ (e.g., Evrard, Metzler & Navarro emn96 (1996)) (Markevitch et al. mark98 (1998); Neumann & Arnaud dmnma99 (1999); Vikhlinin, Forman & Jones vikh99 (1999); Irwin & Bregman irbreg01 (2000); De Grandi & Molendi demol02 (2002); Arnaud, Aghanim & Neumann manadmn02 (2002)). The very existence of these similarities gives strong support to an underlying universality in the dark matter distribution, leading to a pleasing convergence between the observed and simulated properties of galaxy clusters.

However, the temperature profiles in particular have generated much discussion, as rather different profile shapes have been found for similar samples observed by the same satellite (e.g., Markevitch et al. mark98 (1998), White whi00 (2000) \[ASCA\]; Irwin & Bregman irbreg01 (2000), De Grandi & Molendi demol02 (2002) \[BeppoSAX\]). These studies have been hampered somewhat by both PSF issues and sensitivity limits. The former has an inevitable effect on the spatial resolution and is a possible source of systematic uncertainty, the derivation of the profiles being potentially sensitive to the exact correction for the PSF and the detailed modelling of the non-resolved cooling flow component. The latter leads to an inability really to constrain parameters beyond the supposedly isothermal regime, which is expected, from simulations, to extend to $`0.5r_{200}`$. As a direct consequence of this, there are relatively few galaxy clusters for which sufficiently high quality data were available for an accurate determination of the total mass and the corresponding density profile. Furthermore, any systematic uncertainty in the shape of the radial temperature distribution can have a direct effect on the derived mass. For example, the temperature profile obtained by Markevitch et al. (mark98 (1998)) gives mass values that are 1.35 and 0.7 times that derived assuming isothermality at 1 and 6 core radii respectively. As a result, the actual form of the density profile is still a largely untested quantity, at least from an observational point of view.

Clusters can also be used to provide cosmological constraints. For any given cosmology and initial density fluctuation, the mass distribution of virialized objects can be predicted for any given redshift. Constraints on cosmological parameters, $`\sigma _8`$ and $`\mathrm{\Omega }`$, can be found by comparing the predictions with the observed cluster mass function and its evolution (Perrenod perr80 (1980)). For this, however, a great number of accurate observational masses are needed. In the calculation of the observed cluster mass function, the standard way to overcome the paucity of data is to use average cluster temperatures, taking advantage of the tight mass-temperature relation predicted by numerical simulations, where $`MT^{3/2}`$ (e.g., EMN96). While observations have, for hot clusters at least, recovered the slope of this relation, observed masses imply a normalisation consistently lower than found by simulations (e.g., Horner, Mushotzky & Scharf hms99 (1999); Nevalainen, Markevitch & Forman nmf00 (2000); Finoguenov, Reiprich & Böhringer frb01 (2001)). However, these total cluster mass estimates, except in a few cases, required an extrapolation of the data and the level of the discrepancy is sensitive to the assumed temperature profile (e.g. see Horner et al. hms99 (1999), Neumann & Arnaud dmnma99 (1999)).

XMM-Newton and Chandra offer, for the first time, sufficiently good spatial and spectral resolution for self-consistent determinations of global cluster observables such as gas density, temperature and mass profiles. We are now observing clusters with unsurpassed clarity. Chandra, with higher resolution, is the instrument best-suited for the study of cluster cores. In the most recent Chandra study by Allen, Schmidt & Fabian (2001a ), mass-temperature data from 6 clusters are measured up to $`r_{2500}`$, and compared to the reference simulations of EMN96 and Mathiesen & Evrard (me01 (2001)). Once again, a systematic offset of $`40\%`$ is found between the observed and simulated $`MT`$ curves, in the sense that the predicted temperatures are too low for a given mass. XMM-Newton, with its high throughput and large field of view, is the satellite best-matched for the study of the larger scale structure of these objects, and for the determination of essential quantities out to a good fraction of the virial radius. With this capability it is possible to test for other effects, such as potential variations of the normalisation with radius.

In this paper, we use XMM-Newton observations of the relaxed cluster A1413 at $`z=0.143`$ to derive the large scale properties to high resolution, and compare the results to those obtained from both observations and simulations. We address several questions which have been the subject of a large amount of debate in the literature. In particular, we compare our temperature profile with previously derived composite profiles from large samples observed with ASCA and BeppoSAX, and we compare both the form and normalisation of our mass profile with that expected from numerical simulations.

Throughout this paper we use $`H_0=50`$ km s-1 Mpc-1, and unless otherwise stated,$`\mathrm{\Omega }_\mathrm{m}=1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ ($`q_0=0.5`$). In this cosmology, at the cluster redshift of $`z=0.143`$, one arcminute corresponds to 196 kpc.

## 2 Data analysis and preparation

### 2.1 Observations

A1413 was observed in Guaranteed Time for 29.4 ks during XMM-Newton revolution 182 (2000 December 16). Calibrated event files were provided by the XMM-Newton SOC. The MOS and pn data were obtained with the THIN1 and MEDIUM filters, respectively. For the pn data set, we extracted single events, which correspond to PATTERN 0, while for the MOS data sets PATTERNs 0-12 were selected.

Dedicated blank-sky data sets, which consist of several high-galactic latitude pointings with sources removed (Lumb lumb (2002)), were used as background for the whole of this analysis. These data sets are distributed as calibrated event files which have already been treated with the SAS. We extract the background events using the same PATTERN selection criteria as outlined above. In addition, we transformed the coordinates of the background file such that they were the same as for the A1413 data set. In this way we can ensure that all source/background products come from exactly the same regions of the detector, thus minimising detector variations.

### 2.2 Vignetting correction

The method described in Arnaud et al. (2001b ) was used to correct spectra and surface brightness profiles for vignetting effects. Briefly, this method involves weighting each photon with energy $`(E)`$, detected at position $`(x_j,y_j)`$, by the ratio of the effective area at the detected position $`A_{x_j,y_j}(E)`$ to the central effective area $`A_{0,0}(E)`$.

The background data were treated in the same manner as the source. Note that the background component induced by cosmic rays (see below) is not vignetted, but since source and background observations are treated in the same way, the correction factor is the same and thus cancels.

### 2.3 Background subtraction

The XMM-Newton background, consisting of several components, is both time and energy-dependent, and so subtraction is a subtle process. Furthermore it is essential that the subtraction is done correctly, especially so for extended sources like clusters of galaxies, where background effects begin to play a role at large off-centre distances where the surface brightness declines approximately as $`r^4`$ (e.g., Vikhlinin et al. vikh99 (1999)).

The soft proton background due to solar flares cannot be corrected for in the normal fashion (e.g., spectral subtraction) as it displays extreme temporal and flux variability, causing the spectrum to change rapidly with time. At the moment it can only be removed by excising all frames above a certain count-rate threshold, the main effect of which is to considerably reduce the effective exposure time. For these observations, the $`3\sigma `$ threshold for each camera was calculated using the method described in Appendix A, and all frames not meeting this criterion were rejected. In practice the observation is very clean. Note however that the pn is considerably more sensitive to the flares. The final exposure times were 24163s, 24567s and 10254s for MOS1, MOS2 and pn cameras, respectively. The blank-sky backgrounds were cleaned using the same criteria.

The blank-sky background represents effectively the particle induced background, dominant in the hard X-ray band, which is, both spatially and temporally, relatively constant. Nevertheless, this background is variable at the $`10\%`$-level, and so it is frequently necessary to normalise the background. We normalise these observations using the count rate in the \[10-12\] keV and \[12-14\] keV bands, for MOS and pn respectively, treating each camera separately. We varied the normalisation by $`\pm 10\%`$ to assess any systematic uncertainties.

However, the blank-sky data set does not necessarily represent the cosmic X-ray background (CXB), because this is variable across the sky, especially the soft X–ray component (see Snowden et al. snow97 (1997)). We use the method described in Pratt, Arnaud & Aghanim (pratt (2001)), and Arnaud et al. (maetal02 (2002)) to correct for the difference of the CXB. An annular region external to the cluster emission (between 9’ and 13’ in this case) is used to estimate the local background. The normalised spectrum of the same region of the blank-sky background is then subtracted, giving a difference spectrum, which can then be scaled according to the size of any extraction region and subtracted directly from the source spectra. A similar procedure is applied to subtract the residual CXB component for the surface brightness profile (see Arnaud et al. maetal02 (2002) for details).

In addition to the above, the pn data were corrected for ‘out of time events’, which occur when a photon hits the CCD during the read-out process in the imaging mode.

## 3 Morphology

### 3.1 Image

We show in Figure 1 the MOS1+MOS2 image of the cluster, produced simply by adding the data from each camera without accounting for vignetting. The image is striking: the cluster displays an unmistakably elliptical shape, with a clear brightness enhancement directly to the south, and there are a large number of sources in the field of view.

### 3.2 2D $`\beta `$-model fitting

Motivated by the apparent excess of counts to the south of the cluster (see Figure 1) we fitted the MOS1 +MOS2 image with a 2D $`\beta `$-model in order to quantify the significance of this feature. In fitting the image, we followed closely the procedure described in Neumann & Böhringer (neuboh97 (1997)). Images were extracted in the \[0.3-1.4\] keV band from the weighted MOS event files in pixels of size 3″.4 and added to make a combined MOS image. Since in the case of weighted events Poissonian errors do not apply, errors were correctly calculated from the weights using $`\sigma =\sqrt{\mathrm{\Sigma }_jw_j^2}`$ (see Arnaud et al. 2001b ). An error image was generated for each instrument, and these images were added quadratically. The fitting procedure described below was tested and optimised on simulated data before application to the real data.

The $`\chi ^2`$ test used in the fitting procedure assumes Gaussian statistics, for which the mean is the most likely value. The statistics are actually not Gaussian in the external regions of the field of view, dominated by the background. In these regions the number of photons per pixel is low and follows a Poisson distribution for which the mean is larger than the most likely value. If the image is not smoothed before fitting, there is thus a tendency to underestimate the mean background level, leading to erroneous values for the fitted cluster parameters. The combined MOS image was thus smoothed with a Gaussian of with $`\sigma =10\mathrm{}`$ before fitting. The error image was treated according to the error propagation function for Gaussian filtered images described in Neumann & Böhringer (neuboh97 (1997)). We fix all error pixels with a value of 0 to have a value of 1 before fitting, meaning that we can use $`\chi ^2`$ fitting but are unable to determine confidence parameters on the fit. The data were then fitted with a 2D $`\beta `$-model of the form:

$$S(x,y)=S_0(1+F_1+F_2)^{3\beta +\frac{1}{2}}+B$$
(1)
where

$$F_1=\frac{[\mathrm{cos}\alpha (xx_0)+\mathrm{sin}\alpha (yy_0)]^2}{a_1^2}$$

$$F_2=\frac{[\mathrm{sin}\alpha (xx_0)+\mathrm{cos}\alpha (yy_0)]^2}{a_2^2}$$

Here, $`x_0,y_0`$ is the position of the centre of the cluster; $`x,y`$ are the coordinate positions of each pixel; $`a_1,a_2`$ are the major and minor core radii; $`\alpha `$ is the position angle; and the background is included in the model via the parameter $`B`$.

We fit the image between $`1\mathrm{}.3`$ (see below) and $`13\mathrm{}`$ from the cluster center, excluding obvious point sources. The results of the 2D fit are shown in Table 1. Note that the fitted parameters are slightly dependent on the outer radius and the $`\sigma `$ of the Gauss filter, but the results are always in good agreement with the 1D fit, described below.

In order to quantify the significance of the excess of flux to the south of the cluster, we subtract the 2D $`\beta `$-model from the data and calculate the significance of the residuals using the prescription described in the Appendix of Neumann & Böhringer (neuboh97 (1997)). The excess is an extended source detected at $`>10\sigma `$, as shown in Figure 2.

We extracted the spectrum from a circular region of radius $`1\mathrm{}`$ centred on the excess. This spectrum unfortunately does not contain sufficiently strong line emission for a redshift estimate, so we fitted using a MEKAL model with the same redshift as A1413, absorbed with the galactic column density toward the cluster ($`2.19\times 10^{20}`$ cm-2 from Dickey & Lockman dickey (1990)). We find $`\mathrm{k}T=3.1\mathrm{keV}`$. An overlay of the significance contours on the DSS plate of the image did not reveal any obvious sources associated with the excess, and a hardness ratio map did not reveal any interaction with the main cluster. Our tentative conclusion is that the source is either a foreground or background cluster: deeper optical observations of the region should resolve the issue.

## 4 Gas density profile

### 4.1 Surface brightness profile

For each camera, we generated an azimuthally averaged surface brightness profile for both source and background observations. Weighted events from the corresponding event files were binned into circular annuli centred on the position of the cluster emission peak. We cut out serendipitous sources in the field of view and the southern sub-structure. The background subtraction was performed as described in Sec. 2. We consider the profiles in several energy bands. Due to the contribution of the instrumental Al K line around 1.5 keV, we ignored the \[1.4-2.0\] keV band in both cameras. To maximise the signal to noise (S/N) ratio, particularly in the outer cluster region, we choose to base the following on analysis of the \[0.3-1.4\] keV band.

We checked that the vignetting corrected and background subtracted profiles of the three cameras are consistent: they differ only by a normalisation factor within the error bars. We thus coadd the profiles and bin the resulting profile is the following way. Starting from the central annulus, we re-binned the data in adjacent annuli so that i) at least a $`S/N`$ ratio of $`3\sigma `$ is reached and ii) the width of the bin increases with radius, with $`\mathrm{\Delta }(\theta )>0.1\theta `$. Such a logarithmic radial binning insures a $`S/N`$ ratio in each bin roughly constant in the outer part of the profile, when the background can still be neglected.

The resulting surface brightness profile, $`S(\theta )`$, is shown in Fig. 3. The cluster emission is significantly detected up to $`R_{\mathrm{det}}=8.6^{}`$ or $`1.7`$ Mpc.

### 4.2 Density profile modelling

We fitted $`S(\theta )`$ with various parametric models convolved with the XMM-Newton PSF (Ghizzardi ghizzardi (2001), Griffiths & Saxton griffiths02b (2002)), binned into the same bins as the observed profile.

A single $`\beta `$–model cannot account for the data. When the entire radial range is fitted, the reduced $`\chi ^2`$ is $`13`$; for the best fit slope, $`\beta =0.60`$, and core radius, $`\theta _\mathrm{c}=0.66\mathrm{}`$. An excess of emission is readily apparent in the centre and a lower reduced $`\chi ^2`$ is obtained when excluding the central region from the fit. The reduced $`\chi ^2`$ decreases with increasing cut-out radius until it stabilises for $`R_{\mathrm{cut}}1.3^{}`$. In that case we obtained $`\chi ^2=47`$ for 31 d.o.f, with $`\beta =0.71\pm 0.02`$ and $`\theta _\mathrm{c}=1.30\mathrm{}\pm 0.09\mathrm{}`$. The best fit model is plotted as a dotted line in Fig. 3. The $`\beta `$ value is not surprisingly larger than the value ($`\beta =0.62`$) derived by Cirimele, Nesci & Trevese (cnt97 (1997)) from their global fit to the ROSAT profile, but is in excellent agreement with the value $`\beta =0.70\pm 0.02`$ obtained by Vikhlinin, Forman & Jones (vikh99 (1999)) by fitting the outer cluster region ($`\theta >4.9^{}`$). There is also an excellent agreement between the 1D and 2D $`\beta `$-values.

We note that the last two points ($`7.75\mathrm{}<\theta <8.6\mathrm{}`$) lie significantly below the best fit model (a $`3\sigma `$ effect for the last bin). The cluster flux in the last bin is about 16% of the total background and we cannot totally exclude that this discrepancy is an artifact due to remaining systematic uncertainties in the background subtraction. This is further discussed in Sect. 9.4. These last two points are discarded in the present analysis.

For the mass analysis which follows (Sect. 8) it is convenient to have an analytical description of the gas density radial profile ($`n_\mathrm{H}(r)`$) at all radii. We thus tried several alternative parameterisations, with behaviour at large radii similar to a $`\beta `$–model:

* AB model: A cusped profile similar to the NFW universal density profile:

$$n_\mathrm{H}(r)=A\left(\frac{r}{r_\mathrm{c}}\right)^\alpha \left[1+\left(\frac{r}{r_\mathrm{c}}\right)^2\right]^{\frac{3\beta }{2}+\frac{\alpha }{2}}$$
(2)
where $`\alpha `$ is the slope at small radii.
* BB model: We also introduce a double isothermal $`\beta `$ model (BB), assuming that both the inner and outer gas density profile can be described by a $`\beta `$-model, but with different parameters.

$$\begin{array}{ccccc}r<R_{\mathrm{cut}}\hfill & n_\mathrm{H}(r)\hfill & =\hfill & n_{\mathrm{H},0}\hfill & \left[1+\left(\frac{r}{r_{\mathrm{c},\mathrm{in}}}\right)^2\right]^{\frac{3\beta _{\mathrm{in}}}{2}}\hfill \\ r>R_{\mathrm{cut}}\hfill & n_\mathrm{H}(r)\hfill & =\hfill & N\hfill & \left[1+\left(\frac{r}{r_\mathrm{c}}\right)^2\right]^{\frac{3\beta }{2}}\hfill \end{array}$$
(3)
The boundary between the two regions, $`R_{\mathrm{cut}}`$, is a free parameter of the model and we took care that the density distribution is continious across $`R_{\mathrm{cut}}`$, as well as its gradient (for continuity of the total mass profile, see Eq. 12):

$$N=n_{\mathrm{H},0}\frac{\left[1+\left(\frac{R_{\mathrm{cut}}}{r_{\mathrm{c},\mathrm{in}}}\right)^2\right]^{\frac{3\beta _{\mathrm{in}}}{2}}}{\left[1+\left(\frac{R_{\mathrm{cut}}}{r_\mathrm{c}}\right)^2\right]^{\frac{3\beta }{2}}}$$
(4)
and

$$\beta _{\mathrm{in}}=\beta \frac{1+\left(\frac{r_{\mathrm{c},\mathrm{in}}}{R_{\mathrm{cut}}}\right)^2}{1+\left(\frac{r_\mathrm{c}}{R_{\mathrm{cut}}}\right)^2}$$
(5)
* KBB model: We finally consider a generalisation of the $`\beta `$–model for the inner region, allowing a more centrally peaked gas density profile in the core:

$$\begin{array}{ccccc}r<R_{\mathrm{cut}}\hfill & n_\mathrm{H}(r)\hfill & =\hfill & n_{\mathrm{H},0}\hfill & \left[1+\left(\frac{r}{r_{\mathrm{c},\mathrm{in}}}\right)^{2\xi }\right]^{\frac{3\beta _{\mathrm{in}}}{2\xi }}\hfill \\ r>R_{\mathrm{cut}}\hfill & n_\mathrm{H}(r)\hfill & =\hfill & N\hfill & \left[1+\left(\frac{r}{r_\mathrm{c}}\right)^2\right]^{\frac{3\beta }{2}}\hfill \end{array}$$
(6)
where $`\xi <1`$ and

$$N=n_{\mathrm{H},0}\frac{\left[1+\left(\frac{R_{\mathrm{cut}}}{r_{\mathrm{c},\mathrm{in}}}\right)^{2\xi }\right]^{\frac{3\beta _{\mathrm{in}}}{2\xi }}}{\left[1+\left(\frac{R_{\mathrm{cut}}}{r_\mathrm{c}}\right)^2\right]^{\frac{3\beta }{2}}}$$
(7)
$$\beta _{\mathrm{in}}=\beta \frac{1+\left(\frac{r_{\mathrm{c},\mathrm{in}}}{R_{\mathrm{cut}}}\right)^{2\xi }}{1+\left(\frac{r_\mathrm{c}}{R_{\mathrm{cut}}}\right)^2}$$
(8)
The parameters $`\xi `$ and $`r_{\mathrm{c},\mathrm{in}}`$ are strongly correlated. An arbitrary low value of $`\xi `$ can fit the data if no upper limit is put on $`r_{\mathrm{c},\mathrm{in}}`$. The lower limit on this parameter given in Table 2 is obtained by imposing $`r_{\mathrm{c},\mathrm{in}}<1\mathrm{}`$.

The corresponding surface brightness profile is computed numerically by integration of the emission measure along the line of sight:

$$S_\mathrm{X}(\theta )_r^{\mathrm{}}\frac{n_\mathrm{H}^2(b)}{\sqrt{b^2r^2}}𝑑r^2$$
(9)
where $`r=d_\mathrm{A}\theta `$ and $`d_\mathrm{A}`$ is the angular distance. The emissivity in the considered energy band was estimated using an absorbed isothermal model at the cluster mean temperature (given Sect. 5.1), taking into account the instrument response. In the soft energy band considered, this emissivity is insensitive to the observed temperature gradient (shown Sect. 5.2). Note that the profile beyond $`R_{\mathrm{cut}}`$ obtained for the BB and KBB models is a classical $`\beta `$–model. The inner surface brightness profile for the BB model can be analytically computed using incomplete Beta functions.

The best fit models are plotted in Fig. 3, together with the residuals. The corresponding best fit parameters with errors and $`\chi ^2`$ values are given in Table 2. In all cases, the outer slope, $`\beta `$, is consistently found to be similar to the slope obtained by fitting only the outer part of the profile. We found that the AB model does not provide a particularly good representation of the data: the reduced $`\chi ^2`$ is $`\chi _{\mathrm{red}}^22`$ and the residual profile below $`R_{\mathrm{cut}}`$ clearly indicates that the gas distribution is less peaked than a cusped profile. In other words, the gas distribution possesses a core. The best fit is obtained with the KBB model, but the reduced $`\chi ^2`$, $`\chi _{\mathrm{red}}^21.38`$ is still larger than 1. However, the residuals are small (at the $`3\%`$ level on average) and might be due in part to the observed departure from spherical symmetry. As the BB model is a special case of the KBB model ($`\xi `$ fixed to $`\xi =1`$), we can compare both models using a F-test. The KBB model provides a better fit than the BB model at the $`95\%`$ confidence level, suggesting that the density distribution in the core is indeed more centrally peaked than for a conventional $`\beta `$–model. This KBB model is thus adopted for the remainder of the analysis.

## 5 Spatially resolved spectroscopy

### 5.1 Global spectrum

For each instrument, a global spectrum was extracted from all events lying within 5′.1 from the cluster emission peak. This radial range was chosen to maximise the S/N ratio, allowing us to check in detail the consistency between the three cameras. Each global spectrum was fitted with an absorbed MEKAL model with the redshift fixed at $`z=0.143`$. The normalisation for each instrument was left as an additional free parameter. We excluded the energy bins around the strong fluorescence lines of Ni,Cu & Zn from the pn spectrum fit. These lines, present in the background, are not well subtracted by the procedure described in Sec. 2 because they do not scale perfectly with the continuum of the particle-induced background. In all fits we used the following response matrices: m1\_thin1v9q20t5r6\_all\_15.rsp (MOS1), m2\_thin1v9q20t5r6\_all\_15.rsp (MOS2) and epn\_ef20\_sY9\_medium.rsp (pn).

Fitting the data from all instruments above $`0.3\mathrm{keV}`$, with the absorption fixed at the galactic value of $`N_\mathrm{H}=2.2\times 10^{20}`$ cm-2, we found inconsistent values for the temperature derived with the MOS and pn cameras: $`kT=6.91_{0.23}^{+0.23}`$ keV (MOS1), $`6.33_{0.23}^{+0.23}\mathrm{keV}`$ (MOS2) and $`5.76_{0.19}^{+0.19}\mathrm{keV}`$ keV (pn). A better agreement between the cameras, together with a lower $`\chi ^2`$ value, are obtained if the $`N_\mathrm{H}`$ value is let free, but then the best fit $`N_\mathrm{H}`$ values are significantly lower than the 21 cm value (see Table 3). In other words the data presents an excess at low energy as compared to an isothermal model absorbed with the galactic hydrogen column density. This effect could be due to a true soft excess component (e.g Durret et al. durret (2002)) and/or an artifact due to remaining calibration uncertainties. In particular, it is known that the EPIC-pn and MOS cameras show a relative flux difference which increases with energy above $`4.5\mathrm{keV}`$, resulting in a MOS spectral slope flatter than the pn (Saxton saxton (2002); Griffiths et al. griffiths02a (2002)).

We then performed a systematic study of the effect of imposing various high and low-energy cutoffs for each instrument. The $`N_\mathrm{H}`$ is fixed to the 21 cm value. Having first found that progressive cutting of the high energy channels had a negligible effect on the derived temperatures, we then varied the low energy cutoff, for which the results are shown in Table 3.

This table shows that there is an optimum low-energy cutoff for each instrument, above which no amount of further cutting of the low-energy response will significantly affect the temperature. The temperature stabilises above a certain cutoff point for each instrument, this being $`0.6`$ keV for the MOS cameras and $`1.0`$ keV for the pn camera. The adoption of these low-energy cutoffs has the pleasing effect of bringing the temperatures for each instrument into agreement both with each other and with previous ASCA analysis. The combined MOS +pn global temperature is $`\mathrm{k}T=6.85_{0.16}^{+0.15}\mathrm{keV}`$ ($`90\%`$ confidence for one interesting parameter, $`\chi ^2=1459.6`$ for $`1436`$ dof) in agreement with the results of Ikebe et al. (ikebe02 (2002)), who find $`\mathrm{k}T=6.56_{0.44}^{+0.65}\mathrm{keV}`$ and Matsumoto et al. (matsu01 (2001)), who find $`\mathrm{k}T=6.72\pm 0.26\mathrm{keV}`$, and marginally consistent with the result of White (whi00 (2000)), who finds $`\mathrm{k}T=7.32_{0.24}^{+0.26}\mathrm{keV}`$.

It thus appears that the discrepancies observed by fitting the whole energy range are mostly due to some residual calibration uncertainties in the low-energy response of all instruments and/or a true soft excess. The scientific analysis of such a possible soft excess is beyond the scope of this paper. To minimise these effects, we adopted the low-energy cutoffs derived above for the spatially-resolved analysis discussed below.

### 5.2 Radial temperature profile

We produced a radial temperature profile by excluding sources and extracting spectra in annuli centred on the peak of the X-ray emission. All spectra were binned to $`3\sigma `$ above background level (except the final annulus, which was binned to $`2\sigma `$) to allow the use of Gaussian statistics. We show the fitted spectra for annuli 5 and 9 in Figure 4.

These spectra were fitted using the absorbed MEKAL model described above; we fitted separately the spectra from each instrument as well as making a simultaneous MOS+pn fit, as detailed in Table 6. All the temperatures are consistent within the respective errors. It is also evident from Figure 5 that the form of each profile is similar. As a further test, we fitted the annular spectra with an absorbed MEKAL model with the absorption left as a free parameter. This produced profiles with, again, exactly the same shape, giving us high confidence in the form of the profile which we have derived.

### 5.3 Projection and PSF effects

#### 5.3.1 PSF correction

The PSF of XMM-Newton is a potential cause of concern, especially in the inner regions, where the bin sizes are small. To assess the effect of the PSF, we first calculate a redistribution matrix, $`F(i,j)`$, where $`F(i,j)`$ is the fractional flux in annulus $`i`$ coming from annulus $`j`$. These redistribution factors were derived from our best model of the gas density profile, converted to emission measure profile and convolved with the XMM-Newton PSF. The fractional contribution in each bin of the emission coming from the bin, as well as all inner and outer bins are plotted in Fig. 6. The PSF mostly affects the central regions and, above $`2\mathrm{}`$, the contamination from adjacent bins is less than $`25\%`$.

We have a total of 30 spectra (10 annuli $`\times `$ 3 cameras) to be fitted with a model consisting of 10 MEKAL models (corresponding to the 10 ‘true’ temperatures) absorbed by a common (frozen) absorption. The normalisations of the MEKAL models, for annulus $`i`$, are linked by the factors $`F(i,j)`$, such as to reflect the contribution of each annulus $`j`$ due to the PSF. In practice we ignore any contributions at the less than $`1\%`$ level. Each MEKAL model has 6 parameters, which, together with the absorption, makes 7. If we fit the 30 spectra simultaneously, the model has $`(6\times 10+1)\times 30=1830`$ parameters. XSPEC can only handle 1000 parameters, (even if most of them are frozen), so we have to find a way to reduce their total number. One way is to group the spectra, but for this to work the spectra in each group need a common normalisation.

We typically find differences of $`20\%`$ between the normalisations of MOS and pn annular spectra (MOS1 and MOS2 normalisations are the same to within $`5\%`$). We fit the global spectrum discussed in Sect. 5.1 to find the overall difference in normalisation between MOS and pn. We then checked that the annular fit results were the same when the global difference in normalisation was applied. This being so, we multiplied the pn annular spectra by this factor to bring their normalisations into line with those for the MOS cameras.

We now grouped the MOS and pn spectra of each annulus, giving 10 groups of 3 spectra, leaving us with $`(6\times 10+1)\times 10=610`$ parameters, enabling a simultaneous fit. We froze the metallicity of each MEKAL model at the best-fit value found for each projected annulus. The free parameters in the fit are then the temperature and normalisation of each MEKAL model. The resulting PSF-corrected profile is shown in Figure 7 ($`1\sigma `$ errors). The PSF corrected results are entirely consistent with the projected temperature profile, with systematic differences of about half the $`1\sigma `$ errors in the first 3 bins and smaller beyond. This result is not surprising, since the profile is relatively flat. Consideration of the PSF has a much smaller effect on the temperature profile of A1413 than for the bright cooling flow cluster A1835 at $`z=0.25`$ (Markevitch mark02 (2002), Majerowicz et al. majerow (2002)). In contrast to A1835, A1413 displays neither an extremely steep rise in the central gas density, nor a sharp drop in the temperature towards the center. As a result the contamination of central bins is first smaller: for A1835, more than 1/3 of the observed brightness at any radius is due to PSF scattering at smaller radii (Markevitch mark02 (2002), Majerowicz et al. majerow (2002)) while for A1413 this contamination is already less than $`25\%`$ at $`2\mathrm{}`$, and decreases beyond (Fig. 6). Secondly, the smaller temperature gradient towards the center means that the redistribution biases less the temperature determinations.

#### 5.3.2 Spectral deprojection

Another possible source of error in the derived profile comes from projection effects.

A deprojected temperature profile was produced by first simultaneously fitting the MOS and pn spectra of the outer annulus with a MEKAL model absorbed by a fixed galactic column density. The spectrum of the next annulus inward was then fitted with a two-temperature model with the parameters of one of the models fixed to the best-fitting values derived for the outer spectrum. The normalisation of each fixed model must account for the volume within the outer shell projected along the line of sight toward the next shell inward. Furthermore, as the gas density profile is not flat, the normalisation must also account for this effect. We model the gas density profile using the parameters from the best-fit KBB model described in Sect. 4. The normalisations are then adjusted by the emission weighted volume factors. This process was continued inward, adding a MEKAL model for each annulus, with the parameters of the outer annulus models frozen to their previously determined best-fit values. The abundances of the two outer annuli were frozen to the global value, so for these fits the free parameter is the temperature. For all other annuli both the temperature and abundance were free parameters. The deprojected temperature profile is shown compared to the projected profile in Figure 7. In practice we find very little difference between the projected and deprojected results. The jump in the temperature of the ninth annulus is somewhat an artifact of the fitting process. In this case the software tries to compensate for the contribution of the low temperature found in the tenth (and first-fitted) annulus by putting a higher temperature in the subsequent annular bin. The error on the tenth temperature is large, the contribution of the outer emission in the ninth bin depends on the actual cluster extent and thus the deprojected ninth temperature is probably more uncertain than found in this simple procedure. Note, however, that the projected and deprojected temperatures agree well within the $`1\sigma `$ errors.

In summary, we find that neither a consideration of the PSF or projection effects substantially changes the form of the temperature profile. The profiles obtained by taking into account these effects are consistent with the projected profile, within the $`1\sigma `$ errors. For all subsequent analysis, we thus used the projected profile.

### 5.4 Modelling the temperature profile

We now consider the scaled temperature profile, $`\tau (x)=T(r)/T_\mathrm{X}`$, where $`T_\mathrm{X}`$ is the average temperature and $`x`$ is the scaled radius, normalised to $`r_{200}`$. $`r_{200}`$ is estimated from the average temperature $`T_\mathrm{X}`$ and the $`r_{200}`$$`T`$ relation of EMN96 at the cluster redshift. $`T_\mathrm{X}`$ is estimated by fitting the global spectrum, extracted from the $`[0.5\mathrm{}9\mathrm{}]`$ region, i.e outside the possible cooling flow region (see below). We found $`\mathrm{k}T_\mathrm{X}=6.49\pm 0.15\mathrm{keV}`$ ($`1\sigma `$ error). Note that the temperature profile is determined up to $`0.7r_{200}`$ or $`r_{500}`$.

We then modelled this projected temperature profile with a polytropic model:

$$\tau (x)=\tau _0\left(\frac{n_\mathrm{H}(xr_{200})}{n_{\mathrm{H},0}}\right)^{(\gamma 1)}$$
(10)
where $`n_\mathrm{H}`$ is the gas density profile given by the KBB model (Eq.5), $`\gamma `$ is the polytropic index and $`\tau _0`$ is the normalised temperature at $`x=0`$.

We fitted the profile with all $`\beta `$, $`r_\mathrm{c}`$ and $`\xi `$ parameters fixed at the values best fitting the surface brightness profile (Table 2), so that the free parameters for the fit are $`\tau _0`$ and $`\gamma `$. When the whole radial range is fitted, we find $`\tau _0=1.11\pm 0.03`$ and $`\gamma =1.03\pm 0.01`$, with $`\chi ^2=7.62/8`$. If we then exclude the inner point, we find a better fit with $`\tau _0=1.22\pm 0.03`$ and $`\gamma =1.07\pm 0.01`$, and $`\chi ^2=4.45/7`$. The polytropic fits to the temperature profile are shown in Figure 8.

We note that the fit is considerably better if the central point is excluded from the fit. The resulting polytropic profile rises to a peak in the centre which is not seen in the annular temperature determinations, which may lead us to believe that there is a small cooling flow (CF) at work in the very central regions. This possibility is discussed further in Sect. 6.

We further note that the derived value for the $`\gamma `$-parameter is very close to isothermal and moreover, is not very sensitive to whether the central bin is included in the fit. There is a tantalising hint that the temperature profile may drop further in the very outer regions, but the errors on this last data point are large enough that it is easily compatible with the derived polytropic model.

## 6 Properties of the central gas

Our best-fit polytropic model is an excellent description of the observed temperature profile barring the inner point, which is significantly lower. We examined if this temperature drop could be due to a cooling flow.

The cooling time, the enthalpy of the ICM divided by the energy loss due to X-rays, is calculated using:

$$t_{\mathrm{cool}}=\mathrm{2.9\; 10}^{10}\mathrm{yrs}\sqrt{\frac{\mathrm{k}T_\mathrm{X}}{1\mathrm{keV}}}\left(\frac{n_\mathrm{H}}{10^3\mathrm{cm}^3}\right)^1$$
(11)
from Sarazin (sarazin (1986)). Using the central density derived from the KBB model fit (Table 2), we find $`t_{\mathrm{cool}}2.4\times 10^9`$ yr, or about one quarter of the age of the Universe at the cluster redshift. This suggests that a CF should exist. Furthermore, the cooling radius, defined as the radius where the cooling time is equal to the age of the Universe, is $`r_{\mathrm{cool}}0.6\mathrm{}`$ meaning that any CF should reside in the central bin. This is consistent with the observed temperature drop.

We thus fitted the spectrum of the inner bin with more complicated models:

* The sum of two MEKAL models absorbed by a common column density fixed at the galactic value. When fitting, the second temperature is limited to be less than or equal to the temperature of the main component, and the abundances of the two components are tied together.
* The sum of a MEKAL and a cooling flow model, again with a fixed common absorption. Here the abundance of the CF is tied to that of the thermal spectrum, and the upper temperature for the CF is limited to be the temperature of the thermal gas.

The results are shown in Table 4. Both the two temperature and the MEKAL +CFLOW models are better fits than the single temperature model at the $`>99.99\%`$ level. In addition, the MEKAL+CFLOW model is a better description of the data than the two-temperature model at the $`93\%`$ level.

The secondary temperatures and CF properties are not well constrained. More sophisticated modelling is needed, preferably including RGS data, which is beyond the scope of this paper.

## 7 Abundance

The combined MOS+pn spectral fits also allow us to calculate a radial abundance profile, which can be compared with the profiles derived for cluster samples observed with BeppoSAX (Irwin & Bregman ib01 (2001) (IB01); De Grandi & Molendi dm01 (2001) (DM01)). Fig. 9 shows the profile derived in the conservative radial range where we have information from the Fe K$`\alpha `$ line.

The average abundance is $`0.27\pm 0.03`$, more typical of the value found by DM01 for non-CF clusters ($`0.25\pm 0.01`$) than for CF clusters ($`0.34\pm 0.01`$). There is a clear decrease of abundance with radius. However, the profile exhibits a two-step behaviour, rather than a steadily declining profile, as found by DM01 and IB01. It remains more or less constant (a mean of 0.34, excluding the central bin) up to $`400`$ kpc, dropping to another plateau (with a mean of 0.22) beyond. IB01 found that inside $`0.075r_{200}`$ their CF subsample had an average abundance of $`0.46\pm 0.06`$, while the non-CF subsample had an average of $`0.33\pm 0.04`$. The outer region abundances were $`0.30\pm 0.02`$ and $`0.24\pm 0.03`$ for the CF and non-CF subsamples, respectively. The inner two bins of our observation correspond roughly to $`0.075r_{200}`$. The mean error-weighted value for this region is $`0.34\pm 0.08`$ (3$`\sigma `$ errors, for direct comparison with IB01). Outside this region, the abundance value is $`0.28\pm 0.09`$. Thus when the errors are taken into account, we cannot distinguish between the CF and non-CF subsamples of IB01. However, taken at face value, these figures appear to suggest that the abundance profile shape of A1413 displays characteristics intermediate between the CF and non-CF subsamples. This may be because A1413 appears to host only a modest CF. It must also noted that the abundance profiles of A1795 and A2142, both clusters with strong CF signatures (Peres et al. p98 (1998)), appear to have relatively flat abundance profiles, as shown in IB01 and DM01. It has been suggested that both of these clusters are undergoing (or have undergone) mergers (Oegerle & Hill oh94 (1994); Oegerle et al. ohf95 (1995)), which have presumably not been sufficiently strong to disrupt the CF, but which have effectively mixed the metals and thus flattened the radial profile. A study of cluster abundance profiles as a function of the strength of the CF signature and dynamical state would help to better understand the origin of cluster abundance gradients.

## 8 Total Mass profile

### 8.1 Calculation of the mass profile

The mass profile is calculated under the usual assumptions of hydrostatic equilibrium and spherical symmetry. The integrated mass profile can be calculated from the gas density, $`n_\mathrm{g}`$, and temperature profiles:

$$M(r)=\frac{\mathrm{k}Tr}{\mathrm{G}\mu \mathrm{m}_\mathrm{p}}\left[\frac{d\mathrm{ln}n_\mathrm{g}}{d\mathrm{ln}r}+\frac{d\mathrm{ln}T}{d\mathrm{ln}r}\right]$$
(12)
where G and $`\mathrm{m}_\mathrm{p}`$ are the gravitational constant and proton mass and $`\mu =0.597`$.

If the gas density profile is described by the KBB model (Eq. 5), then the mass profile is described by:

$`r<R_{\mathrm{cut}}M(r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{k}r^2}{\mathrm{G}\mu \mathrm{m}_\mathrm{p}}}\left[{\displaystyle \frac{3\beta _{\mathrm{in}}r^{(2\xi 1)}T(r)}{r^{2\xi }+r_{\mathrm{c},\mathrm{in}}^{2\xi }}}+{\displaystyle \frac{dT}{dr}}\right]`$
$`r>R_{\mathrm{cut}}M(r)`$ $`=`$ $`{\displaystyle \frac{\mathrm{k}r^2}{\mathrm{G}\mu \mathrm{m}_\mathrm{p}}}\left[{\displaystyle \frac{3\beta rT(r)}{r^2+r_\mathrm{c}^2}}+{\displaystyle \frac{dT}{dr}}\right]`$ (13)
where $`r_{\mathrm{c},\mathrm{in}}`$, $`\xi `$, $`r_\mathrm{c}`$ and $`\beta `$ are the parameters of the KBB model, $`\beta _{\mathrm{in}}`$ being linked to them by Eq. 8.

The mass profile is calculated using the Monte Carlo method of Neumann & Böhringer (neubo95 (1995)). The gas density profile parameters are fixed to their best fit values. This method calculates random temperature profiles within the bounds of the observed profile. We calculated 10000 such random temperature profiles, using a window size of 150kpc and a smoothing parameter of $`0.1\mathrm{keV}`$. The final output is the mean mass profile and the corresponding errors for each data point in the input temperature profile. The errors are calculated using the $`90\%`$ errors in the temperature profile and the standard deviation of the mass at any given radius. The resulting errors are calculated to correspond to $`1\sigma `$ errors in the mass profile.

The errors on the mass profile due to the error on the density gradient, $`d\mathrm{ln}n_\mathrm{g}/d\mathrm{ln}r`$, are then calculated. As the parameters $`r_{\mathrm{c},\mathrm{in}}`$, $`\xi `$, $`r_\mathrm{c}`$ and $`\beta `$ are correlated, this error cannot be deduced directly from the errors on these parameters. We used a method similar to the one described in Elbaz et al. (elbaz (1995)). For each considered radius, the surface brightness profile is fitted, considering $`d\mathrm{ln}n_\mathrm{g}/d\mathrm{ln}r`$ estimated at this radius as a free parameter, instead of $`\beta `$. The $`1\sigma `$ error on this parameter is then classically derived from the $`\chi ^2`$ variation, the other parameters (normalisation, $`r_{\mathrm{c},\mathrm{in}}`$, $`r_\mathrm{c}`$ and $`\xi `$) being optimised. Finally the errors due to the density and temperature profiles (derived from the Monte-Carlo method) are added quadratically.

The resulting mass profile, with $`1\sigma `$ error bars, is plotted in Fig. 10.

### 8.2 Factors influencing the mass profile

The temperature profile of A1413 is well determined out to $`0.7r_{200}`$, and shows a gradual decline which is well described by a polytrope of index $`\gamma =1.07`$ beyond the CF region. Assuming such a polytropic description, the mass profile can be calculated analytically from the best fit gas density KBB model. The polytropic mass profile lies well within the errors of the Monte Carlo profile, except for the central point due to the temperature drop observed in the center (full line in Fig. 10). We also note that the derived mass at large radii ($`r>1.3`$ Mpc) lies at the lower range of the Monte Carlo mass. This is due to the drop of temperature in the last radial bin, the best fit temperature being below the polytropic value.

In the classic approach, the gas density is described with an isothermal $`\beta `$-model, in which the temperature profile is assumed to be isothermal and the gas density distribution is parameterised by a $`\beta `$–model. In Fig. 10 (dotted line) we show the mass profile obtained using this approach, with the $`\beta `$–model best fitting the outer part of the cluster and the average cluster temperature outside the CF region, $`\mathrm{k}T_\mathrm{X}`$. Not surprisingly, the mass is greatly underestimated in the centre ($`r<R_{\mathrm{cut}}`$), where the gas density profile is more concentrated (higher gradient) than the extrapolated $`\beta `$–model. If we instead parameterise the gas density using the best fit KBB model, the mass distribution towards the centre is recovered (Fig. 10, dashed line). Beyond $`R_{\mathrm{cut}}`$ the mass profile is slightly steeper than that derived from the true temperature profile, as expected from the observed $`\gamma `$ value, slightly larger than 1. This comparison shows that the temperature gradient has a small but systematic effect on the derived mass profile.

The mass profile is remarkably well constrained: the $`1\sigma `$ error is less than $`\pm 5\%`$ below $`1.4`$ Mpc and rises to $`\pm 18\%`$ at 1.8 Mpc. The temperature logarithmic gradient is much smaller than the density logarithmic gradient ($`7\%`$ for $`\gamma =1.07`$), except in the very outer part, where the temperature gradient is both larger and the constraints are poorer. As a consequence (see Eq. 12), except in this outer region, the error on the mass profile is dominated by the error on the gas density gradient (in the range $`0.5\%3\%`$) and on the average temperature ($`2.3\%`$). For the same reason the mass profile is very robust versus possible systematic errors on the temperature profile. We have shown in Sect. 5.2 that spectral deprojection or PSF correction do not have a significant effect on the form of the profile. One might also ask what effect the ellipticity of the cluster might have on the derived radial quantities. We extracted spectra in elliptical annuli and compared the projected temperature profile with that produced using circular annuli. All temperatures agree within the respective errors, and so we conclude that the cluster ellipticity is also a minor source of error.

### 8.3 Modelling of the mass profile

Navarro, Frenk & White (nfw97 (1997), NFW) performed high resolution N-body simulations which showed that the density profiles of dark matter halos have a universal shape, regardless of halo mass and values of cosmological parameters. The NFW profile is given by:

$$\rho (r)=\frac{\rho _\mathrm{c}(z)\delta _c}{(r/r_\mathrm{s})(1+r/r_\mathrm{s})^2}$$
(14)
where $`\rho (r)`$ is the mass density and $`\rho _\mathrm{c}(z)`$ is the critical density at the observed redshift, which, for a matter dominated $`\mathrm{\Omega }=1,\mathrm{\Lambda }=0`$ Universe is:

$$\rho _\mathrm{c}(z)=\frac{3H_0^2}{8\pi G}(1+z)^3.$$
(15)
The parameters of the model are $`r_\mathrm{s}`$, a scale length and $`\delta _\mathrm{c}`$, a characteristic dimensionless density dependent on the formation epoch of the dark matter halo. $`\delta _\mathrm{c}`$ can be expressed in term of the equivalent concentration parameter, $`c`$:

$$\delta _\mathrm{c}=\frac{200}{3}\frac{c^3}{[ln(1+c)c/(1+c)]},$$
(16)
The radius corresponding to a density contrast of 200 is $`r_{200}=cr_\mathrm{s}`$.

The NFW density profile varies from $`\rho _{\mathrm{NFW}}r^1`$ at small radii to $`\rho _{\mathrm{NFW}}r^3`$ at large radii. As we are fitting the mass profile $`M(r)`$, we use the integrated mass of the NFW profile for the fit (e.g Suto, Sasaki & Makino suto98 (1998)):

$`M(r)`$ $`=`$ $`4\pi \rho _\mathrm{c}(z)\delta _cr_\mathrm{s}^3m(r/r_\mathrm{s})`$ (17)
$`m(x)`$ $`=`$ $`ln(1+x){\displaystyle \frac{x}{1+x}}`$ (18)
More recent, higher resolution simulations by Moore et al. (moore99 (1999), hereafter MQGSL) suggest a profile described by:

$$\rho (r)=\frac{\rho _\mathrm{c}(z)\delta _\mathrm{c}}{(r/r_\mathrm{s})^{3/2}\left[1+\left(r/r_\mathrm{s}\right)^{3/2}\right]},$$
(19)
where

$$\delta _\mathrm{c}=\frac{100c^3}{\mathrm{ln}(1+c^{3/2})}$$
(20)
This is essentially identical to the NFW profile at large radii but is steeper near the centre ($`\rho _mr^{1.5}`$). Again, as we are fitting the mass profile, we use the integrated mass of the MQGSL profile, given by (Suto et al. suto98 (1998)):

$`M(r)`$ $`=`$ $`4\pi \rho _\mathrm{c}(z)\delta _cr_\mathrm{s}^3m(r/r_\mathrm{s})`$ (21)
$`m(x)`$ $`=`$ $`{\displaystyle \frac{2}{3}}\mathrm{ln}(1+x^{3/2})`$ (22)
The derived parameters from each fit are given in Table 5 and the best fit models are compared to the data in Fig. 11. We find that the data are extremely well described by the MQGSL profile across the entire radial range. The NFW profile can also be used to describe the data, but shows a small divergence at small radii ($`\chi ^2=4.54`$ for the first 3 points). The radius corresponding to a density contrast of $`500`$, computed from the data, is indicated by an arrow. There is a slight hint that the measured mass, $`M_{500}=7.7_{0.8}^{+1.2}\times 10^{14}\mathrm{M}_{}`$, is higher than the MQGSL and NFW models ($`6.8\times 10^{14}\mathrm{M}_{}`$ and $`6.5\times 10^{14}\mathrm{M}_{}`$, respectively) around that radius (see Sect. 9.3 for further discussion).

## 9 Discussion

### 9.1 The shape of the temperature profile

It is instructive to compare the projected XMM-Newton temperature profile of A1413 with the composite profiles found for larger cluster samples. The most extensive samples come from ASCA and BeppoSAX data; these are, in order of publication: Markevitch et al. (1998; MFSV98), White (2000; W00) Irwin & Bregman (2000; IB00) and De Grandi & Molendi (2002; DM02). The MFSV98 ASCA-derived profile is sharply decreasing, such that for a typical $`7\mathrm{keV}`$ cluster the temperature drop is characterised by a polytropic index of 1.2-1.3. W00 finds that $`90\%`$ of the cluster profiles in his ASCA sample are consistent with isothermality at the $`3\sigma `$-level. The IB00 BeppoSAX-derived profile extends only out to $`0.3r_{200}`$ and is flat or even slightly increasing. In contrast, the overall DM02 profile, from a larger sample of BeppoSAX observations, is characterised by an isothermal core extending to $`0.2r_{200}`$. Their CF subsample exhibits a temperature drop of a factor of 1.7 between $`0.2r_{200}`$ and $`0.5r_{200}`$. The non-CF clusters exhibit a sharper temperature drop in the outer regions. DM02 suggest that an incorrect treatment of the BeppoSAX strongback may explain the discrepancy between their result and that of IB00.

We have derived the projected temperature profile of A1413 out to $`0.7r_{200}`$ or $`r_{500}`$, in much finer detail than is possible with either ASCA or BeppoSAX. All indications are that A1413 is a relaxed cluster.

Our data are compared to the DM02 and MFSV98 composite profiles in Fig. 12. Although each individual data point is (marginally) consistent with the typical region defined by these composite profiles, there is an obvious systematic difference in shape. The A1413 profile does not decline sharply like the composite profile of MFSV98, or the profile of DM02 beyond $`0.2r_{200}`$. In Sect. 5.4, we show that the polytropic model gives an acceptable fit to the data. Excluding the central bin, the $`\gamma `$ value (1.07) implies an almost isothermal temperature profile, and is not compatible with that found by MFSV98. It is very similar to that found by DM02 for CF clusters, but DM02 reject the polytropic model on the grounds of poorness of fit, which is not surprising given the decline of a factor of 1.7 in the temperature of their composite profile between $`0.2`$ and $`0.5r_{200}`$. We do not see a similar decline, and so a polytropic model is a good fit to these data. On the other hand, their best fit broken line model is a poor fit to our data: we find a $`\chi ^2=22.1/6`$ for their CF best-fit, and the fit is worse for their non-CF relation ($`\chi ^2=27.7/7`$).

It is obvious that, given the extra radial range afforded by these XMM-Newton data, a flat or increasing profile, such as that of W00 or IB00 extrapolated to high radii, does not describe the A1413 data either. We emphasize again that the temperature gradient is modest: the temperature decreases by $`15(20)\%`$ between $`0.1r_{200}`$ and $`0.3(0.5)r_{200}`$. Fitting the temperature profile up to $`0.3r_{200}`$ (i.e., excluding the last three temperature bins, and the inner bin) with a polytrope allows us to compare our profile directly with that of IB00. We find $`\tau _0=1.19\pm 0.03`$, $`\gamma =1.06\pm 0.01`$, consistent with the value derived for the full radial range. This gradient is in agreement with the level of isothermality found by W00 and IB00 in that radial range taking into account their errors, as well as that found by Allen et al. (2001a ) from Chandra data below $`r_{2500}0.3r_{200}`$. Our observation is also consistent with other XMM-Newton observations of nearby clusters, e.g. the slightly decreasing XMM-Newton temperature profile of Coma ($`10\%`$ at $`0.2r_{200}`$, Arnaud et al. 2001a ) and the temperature profile of A1795, measured up to $`0.4r_{200}`$ and found to be flat within $`\pm 10\%`$ beyond the CF region ($`0.1r_{200}`$, Arnaud et al. 2001b ).

### 9.2 Shape of the total mass profile

In Sec. 8.3 we showed that the NFW form can describe the total mass profile of A1413. However a slightly better agreement in the center is obtained with a MQGSL profile, derived from higher resolution simulations.

With Chandra, it is possible to examine the central regions in great detail, at the expense of information at large radii. At present, it is unclear whether the NFW or MQGSL profiles provide the better description of the mass profiles derived from Chandra observations. Allen, Schmidt & Fabian (2001b ) investigate several forms for the mass profile of RXJ1347.5-1145, finding that both the NFW and MQGSL provide an acceptable fit, although the NFW profile is favoured in terms of $`\chi ^2`$. Perhaps the highest resolution examination of a cluster mass profile is that of Hydra A by David et al. (david01 (2001)), who find $`\rho r^{1.3}`$ between 30 and 200 kpc, which is intermediate between the NFW and MQGSL forms. The addition of a mass point from Hα observations leads them to favour the NFW profile, although the result is still consistent with the MQGSL result.

One sticking point is the value of the concentration parameter from the NFW fit by David et al. (david01 (2001)). They find $`c=12`$, which is 3 times that expected for a cluster of the mass of Hydra A. Interestingly, a similar value of $`c`$ was found by Arabadjis, Bautz & Garmire (abg02 (2002)) from a Chandra study of EMSS 1358+6245. On the other hand, the c parameters of Allen et al. (2001a ) are better in agreement with the theoretical predictions. At large radii the NFW and MQGSL profiles coincide and we can compare the c value derived from our NFW fit, $`c=5.4\pm 0.2`$ for a $`M_{200}10^{15}h_{50}^1\mathrm{M}_{}`$ cluster, to numerical simulations. Teyssier (teyssier (2002)) derived c parameters in the range $`4.99.5`$ for 5 clusters in this very mass range. The average $`c`$ parameter derived by Eke et al. (eke (1998)) for hot massive clusters is $`c6`$. It must be noted that a relatively large dispersion on this parameter is expected from numerical simulations, with a $`1\sigma `$ $`\mathrm{\Delta }(\mathrm{log}c)=0.18(50\%)`$ at a given mass (Bullock et al. bullock (2001)). In conclusion, we emphasize the excellent agreement in shape between the mass profile derived for A1413 and the theoretical expectations, all the more remarkable in view of the very small statistical errors on the profile.

### 9.3 Normalisation of the $`M`$$`T`$ relation

We now examine the normalisation of the mass profile. We will classically define $`M_\delta `$, the mass within a given radius $`r_\delta `$, inside which the mean mass density is $`\delta `$ times the critical density, $`\rho _\mathrm{c}(z)`$ at the cluster redshift. For clusters obeying HE and self-similarity, the mass $`M_\delta (T,z)`$, scales with the cluster temperature and redshift as:

$$h(z)M_\delta (T,z)=M_{10}(\delta )T_{10}^{3/2}$$
(23)
where $`h(z)`$ is the Hubble constant normalised to its local value and $`M_{10}(\delta )`$ is the normalisation at density contrast $`\delta `$ (here $`T_{10}`$ is arbitrarily expressed in unit of $`10\mathrm{keV}`$). The above relation is remarkably well verified by adiabatic numerical simulations, down to $`\delta 200`$, for both SCDM and $`\mathrm{\Lambda }\mathrm{CDM}`$ cosmology, with $`M_{10}(\delta )`$ independent of cosmology (e.g Mathiesen math (2001)). The variation of the normalisation $`M_{10}(\delta )`$ with $`\delta `$ reflects the (universal) cluster internal structure and is the same for all clusters for a given density contrast (although some scatter is present in practice).

In Figure 13 we show the scaled mass profile of A1413 defined as:

$`\stackrel{~}{M}(\delta )`$ $`=`$ $`M(r)h(z)\left({\displaystyle \frac{\mathrm{k}T_\mathrm{X}}{10\mathrm{k}\mathrm{e}\mathrm{V}}}\right)^{3/2}`$
$`\mathrm{where}`$
$`\delta `$ $`=`$ $`{\displaystyle \frac{3M(r)}{4\pi r^3\rho _\mathrm{c}(z)}}`$ (24)
with

$`\rho _\mathrm{c}(z)`$ $`=`$ $`h^2(z){\displaystyle \frac{3H_0^2}{8\pi G}}`$ (25)
and $`h(z)=(1+z)^{3/2}`$ for the SCDM cosmology considered here. This profile can be compared with the normalisations $`M_{10}(\delta )`$ derived from numerical simulations, allowing us to check the normalisation of the $`M_\delta `$$`T`$ relation at different density contrasts. The results of EMN96 are indicated with a thick line. We also plot the scaled profiles corresponding to our best-fit MQGSL model (thin line), and to our best fit polytropic model (dashed line). The density contrast is computed self-consistently for each profile using Eq. 24.

Both the data points and the MQGSL and polytropic model curves run parallel to the EMN96 profile down to $`\delta 600`$. This simply reflects the excellent agreement in shape of the A1413 profile with numerical simulations, as outlined above. However there is a very significant offset in normalisation: $`\stackrel{~}{M}(2500)=7.3\pm 0.2\times 10^{14}\mathrm{M}_{}`$, compared with $`M_{10}(2500)=9.95\times 10^{14}\mathrm{M}_{}`$ from EMN96. In other words, the predicted normalisation of the $`M`$$`T`$ relation lies $`36\%`$ higher than the observed value, in excellent agreement with the Chandra finding of Allen et al. (2001a , $`40\%`$ at $`\delta =2500`$).

Below $`\delta =600`$ the observed profile levels off, so that the data points seem to converge towards the EMN96 predictions at small $`\delta `$. At $`\delta =500`$ $`\stackrel{~}{M}(500)=1.8_{0.2}^{+0.3}\times 10^{15}\mathrm{M}_{}`$, only $`20\%`$ lower than the EMN96 value of $`M_{10}(500)=2.2\times 10^{15}\mathrm{M}_{}`$ and actually consistent with it, especially if we also take into account the dispersion observed in the simulations ($`15\%`$). If this effect is real, this would point to a fundamental difference of form in the total mass profile at large radii, i.e the real cluster dark matter density profile drops less steeply than the canonical $`r^3`$ law. The observed discrepancy in the normalisation of the $`M`$$`T`$ relation at high density contrast would thus be due to a flaw in the numerical simulations for the dark matter component.

However, it is more likely that this level off of the observed profile is an artifact due to incomplete virialisation. If there is residual kinetic energy due to infall, the HE equation applied to the observed temperature profile would over-estimate the true mass. Such an effect is observed in the simulations of EMN96, although its expected magnitude is somewhat smaller that we observe, $`10\%`$ at $`\delta =500`$ on average (but a large scatter exists). A further indication that this is the correct explanation comes from a comparison with the polytropic temperature model. Let us assume that the best fit MQGSL model indeed reflects the true total mass distribution. The fact that the mass profile derived from the best fit polytropic model closely follows this profile down to $`\delta =200`$ (see Fig. 13), indicates that this model is a true representation of the thermodynamic state of the gas if it was in HE up to there. The drop in temperature in the last bin ($`\delta >550`$), as compared to this model, would thus be an direct indication of incomplete thermalisation. We also recall that there is a sudden drop in the surface brightness profile at $`\theta =7.8\mathrm{}`$, corresponding to a density contrast of 450 (computed with the MQGSL model). This further supports our interpretation: we might actually be seeing the expected drop of the X-ray brightness beyond the edge of the virialised (and hot) part of the cluster. Finally, incomplete equipartition between the electrons and the ions at the border of the cluster (Chieze, Alimi & Teyssier chieze (1998)) could also contribute to the low (electronic) temperature observed.

The observed scaled mass profile depends on the cosmological model via the function $`h(z)`$, used in the scaling, and the angular distance, $`d_\mathrm{A}`$, used to convert angular radius to physical radius. On the other hand the theoretical normalisation $`M_{10}(\delta )`$ appears to be insensitive to cosmology. We thus examined if a better agreement with the theoretical normalisation is obtained for the currently most favoured $`\mathrm{\Lambda }`$CDM model ($`\mathrm{\Omega }_\mathrm{m}=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$). From Eq. 12, Eq. 24 and Eq. 25, the derived $`\delta `$ value scales as $`(d_\mathrm{A}h(z))^2`$ and $`\stackrel{~}{M}`$ as $`d_\mathrm{A}h(z)`$, with $`h^2(z)=\mathrm{\Omega }_\mathrm{m}(1+\mathrm{z})^3+\mathrm{\Omega }_\mathrm{\Lambda }`$. For a $`\mathrm{\Lambda }`$CDM model (open circles in Fig. 13), as compared to the SCDM model (filled circles), the data points are moved down and left along a line of slope 1/2 in the log-log plane, with $`\delta `$ multiplied by 1.137 and $`\stackrel{~}{M}`$ multiplied by 0.938. The translation is modest and its slope is similar to the slope of the scaled mass profile around $`\delta =1000500`$, so the observed scaled mass profile remains essentially unchanged and the agreement with the theoretical curve is no better.

In summary the good agreement between the mass profile shape and the numerical simulations, measured for the first time in the whole virialised part of the cluster, suggest that the modelling of the Dark Matter component is correct. However, the offset in the normalisation of the $`M`$$`T`$ relation suggests that some physics is lacking in the modeling of the gas. Several groups have studied the effect of non gravitational physics, like pre-heating or cooling, on the $`M`$$`T`$ relation (e.g. Loewenstein loewenstein (2000); Bialek et al. bialek (1999); Tozzi & Norman tozzi (2001); Babul et al. babul (2002); Thomas et al. thomas (2002); Voit et al. voit (2002)). A detailed comparison with numerical simulations, which would require a statistically representative sample, is beyond the scope of this paper. We simply note that pure pre-heating seems to have a small effect of the $`M`$-$`T`$ relation in the high temperature range of A1413 (Loewenstein loewenstein (2000); Tozzi & Norman tozzi (2001); Babul et al. babul (2002)) and that models including cooling (e.g. Thomas et al. thomas (2002)) seem to be more successful. We also emphasize that some care must be exercised when comparing theoretical predictions with these X–ray observations. The magnitude of the effect is not large as compared to the typical difference of $`50\%`$ (Henry henry (2000)) in the normalisation derived by various groups using purely adiabatic simulations. The dispersion of the relation, observed for both simulated ($`20\%`$) and real clusters, require the use of statistically well controlled samples. Finally, a further ambiguity lies in the definition of the temperature. The temperature profile is not exactly isothermal (although our data suggest that the departure is small). Ideally we should compare the data from a given instrument using the spectral temperature estimated from simulation, after full modelling of the plasma emission folded with the instrument response. The study of Mathiesen & Evrard (me01 (2001)) indicates that the spectral temperature could be an underestimate of the mass-weighted temperature, by about $`20\%`$. Note that this effect would worsen the discrepancy observed above. All these systematic effects have to be well controlled, if we want to confirm the departure from the self-similar scaling and identify the physical process responsible for it.

### 9.4 Gas distribution

From our determination of the central gas density we have calculated a cooling radius of $`r_{\mathrm{cool}}=0^{}.6`$. Figure 3 shows that the gas density has begun to rise above the $`\beta `$–model fit at about twice this radius. This suggests that the main driver of the central peak in the gas distribution is the cusp in the dark matter profile, not the CF. It is likely that the cusp in the dark matter profile acts as a focus for the gas and the CF. This has implications for CF mass deposition rates $`\dot{M}`$ deduced from excesses to the $`\beta `$–model, in that any application to this cluster would result in a gross overestimate of $`\dot{M}`$.

The ratio of the X-ray gas mass to the total gravitating mass is shown as a function of density contrast $`\delta `$ in Figure 14. The $`f_{\mathrm{gas}}`$ rises slowly with increasing radius, up $`\delta 1000`$. After this point, the $`f_{\mathrm{gas}}`$ either stabilises within the error if we use the total mass data points, or continues increasing beyond if we used the best fit MQGSL total mass profile. A conservative estimate of the gas mass fraction at $`\delta 500`$ is $`f_{\mathrm{gas}}=0.2\pm 0.02`$. This shows that the hot ICM is more extended than the dark matter distribution, as has been found in previous studies (e.g., David, Jones & Forman djf95 (1995)), and is also seen in numerical simulations, and expected from purely dynamical reasons (Chieze, Teyssier & Alimi chieze97 (1997)).

Assuming that the gas density profile follows the dark matter profile at large radii (within a factor of 2 between $`r_{200}/2`$ and $`2r_{200}`$) and a polytropic equation of state, Komatsu & Seljak (komatsu (2001)) derived an analytical solution of the gas in HE in an NFW profile. They obtained $`\gamma =1.15+0.01(c0.5)`$, or $`\gamma =1.14`$ for the concentration parameter $`c=5.4`$ derived for A1413, and a X-ray outer slope of $`\beta 0.65`$ for a typical 6.5 keV cluster (their Fig. 14). The first assumption is roughly verified for our best fit profiles and not surprisingly, our derived values are in fair agreement with their model, although we obtain a slightly but significantly lower value of $`\gamma =1.07\pm 0.01`$ and a slightly higher $`\beta `$ value.

The value at $`\delta =600`$, corresponding roughly to the virialised part of the cluster, $`f_{\mathrm{gas}}0.2`$, can be used to calculate the total mass density in the Universe, following the arguments that assume that the properties of clusters constitute a fair sample of those of the Universe as a whole (e.g., White et al. whi93 (1993)). Assuming that the luminous baryonic mass in galaxies in A1413 is approximately one-fifth of the X-ray ICM mass (e.g., White et al. whi93 (1993)), and neglecting other possible sources of baryonic dark matter, $`\mathrm{\Omega }_m=(\mathrm{\Omega }_b/1.2f_{\mathrm{gas}})`$, where $`\mathrm{\Omega }_b`$ is the mean baryon density. For $`\mathrm{\Omega }_bh_{50}^2=0.0820`$ (O’Meara et al. ome (2001), $`h_{50}`$ is the Hubble constant in units of 50 km s-1 Mpc-1), we obtain $`\mathrm{\Omega }_m=0.34h_{50}^{0.5}`$.

## 10 Conclusions

The main conclusions of this work may be summarised as follows:

1. We have reported XMM-Newton observations of A1413, a relaxed galaxy cluster at z=0.143.
2. In a 2D $`\beta `$ model analysis, we detect substructure to the south which does not appear to be interacting with the main cluster.
3. Excluding the data from this region and all obvious point sources, we have measured the gas density and temperature profiles up to $`r_{500}`$ (corresponding to a density contrast $`\delta 500`$, with respect to the critical density at the redshift of the cluster). With the assumptions of HE and spherical symmetry, we have calculated the mass profile out to the same distance.
4. The gas density profile is well described with a $`\beta `$-model beyond $`250h_{50}^1\mathrm{kpc}`$. We further parameterise the inner regions with a modified version of the $`\beta `$-model (the KBB, Eq. 5), which allows a more centrally peaked gas distribution.
5. The temperature profile (excluding the inner point) is well described by a polytropic model with $`\gamma =1.07\pm 0.01`$. The decline is modest: a decrease of $`20\%`$ between $`0.1r_{200}`$ and $`0.5r_{200}`$.
6. The mass profile, derived from the HE equation, is determined with an accuracy of about $`\pm 5\%`$ up to $`r_{600}`$ and $`\pm 18\%`$ at $`r_{500}`$. It can be remarkably well described by a Moore et al. (moore99 (1999)) profile with a scale radius $`r_\mathrm{s}=845\pm 43`$ kpc and concentration parameter $`c=2.6\pm 0.1`$. An NFW profile also gives an acceptable fit but describes less well the central regions. The $`c`$ values we find are in good agreement with those expected from numerical simulations for a cluster of this mass. The Dark Matter modelling in these simulations is thus strongly supported by the excellent agreement between observed and simulated profiles.
7. Beyond $`r_{600}`$, the observed temperature and derived mass profiles begin to depart systematically from, respectively, the polytropic description and the Moore et al. (moore99 (1999)) profile. There is also a sudden drop of the surface brightness profile at $`r_{450}`$. This suggests that the gas in these regions may not be in HE, and we may thus be seeing the outer edge of the virialized parts of this cluster.
8. The offset in the normalisation of the $`M_\delta T`$ relation, with respect to the simulations of Evrard et al. (1996) is now confirmed to be $`40\%`$ across the entire radial range up to $`r_{500}`$ (i.e., in the virialised part of the cluster).
9. The gas distribution is peaked primarily as a result of the cusp in the dark matter profile. The gas mass fraction increases with increasing radius, to reach $`0.2`$ at $`r_{500}`$.

We are now in a position directly to confront simulations with observations. The results are encouraging (the obvious validity of the modelling of the Dark Matter distribution at large scale) but many questions remain. How peaked are dark matter profiles? What is the relationship between between central dark matter cusps and CFs? Why are some studies finding unrealistic values of the concentration parameter? What is the source of the discrepancy in the M-T relation?

The statistical errors on the observed quantities are now small enough so that we can determine in detail the intrinsic dispersion in cluster properties and systematic discrepancies with the classical self-similar model. To answer the above and other questions, a statistical sample of cluster properties would be of great help, preferably using Chandra to probe the central regions and XMM-Newton to determine properties at great distances from the cluster centre. Confrontation with numerical simulation is essential. The full range of observations, correlations between cluster properties, and detailed internal gas structure should be derived taking into account that they are viewed through a given instrument, so that we are able truly to compare like with like.

###### Acknowledgements.

We thank the referee, J. Irwin, for insightful comments and suggestions which have improved the paper. We thank T. Ponman for providing the script used in the deprojection of the temperature profile, and D. Neumann for providing the Monte-Carlo code to derive the mass profile. We thank R. Teyssier and A. Refregier for useful discussions, and S. De Grandi for providing the BeppoSAX results shown in Fig. 12. We thank A.Evrard, M.Markevitch and T.Ponman for their comments on the manuscript. The present 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).

## Appendix A Cleaning of data for soft proton solar flares

Various methods have been used to remove background flares from XMM-Newton and Chandra data sets. For XMM-Newton observations, perhaps the most simple is direct visual inspection of the binned high-energy \[10-12\] keV light curve over the whole field of view, the adoption of a threshold level, and the exclusion of any intervals above the selected threshold (see e.g., Arnaud et al. 2001). For Chandra data, an $`n\sigma `$ clipping is used, where the average count rate in the \[3 \- 6\] keV band over the field of view is calculated, thresholds are set depending on this value, and the light curve cleaned for any intervals where the thresholds are not met.

It is known that the XMM-Newton quiescent background level is variable by $`\pm 10\%`$, and so it is not possible to set a rigid threshold level for flare rejection because of the risk of losing good data. The threshold level should ideally be dependent on the quiescent rate of the observation in question.

For the A1413 observation, we extracted the MOS \[10-12\] keV light curves in the field of view in 104s bins (chosen as this is an integral multiple of the frame readout time). Similar light curves were made for the pn, but in the \[12-14\] keV band. We then made a histogram of each light curve and fitted this histogram with a Poisson distribution

$$y=\frac{\lambda ^xe^\lambda }{x!}$$
(26)
where the mean of the distribution, $`\lambda `$, is the free parameter of the fit. Following Poisson statistics, the error on the mean, $`\sigma =\sqrt{\lambda }`$. We found that the Poisson distribution was, without exception, a better fit than a Gaussian distribution. We then defined thresholds at the $`\pm 3\sigma `$ level, and rejected any time intervals outside these thresholds. We show the Poisson fit and the original and cleaned light curves for A1413, which is a quiet observation, and for comparison, the observation of MKW9 (Neumann et al., in preparation), which shows several solar flares, in Figures 15 and Figures 16.

This method is effective at finding the quiescent periods, even for data strongly affected by flares. As such it is not prone to the overestimation of the mean rate, the main problem with the $`n\sigma `$ clipping method.

## Appendix B Table


format: text

text:
# Secondary Neutrinos from Tau Neutrino Interactions in Earth

## Abstract

The energy dependence of “secondary” neutrinos from the process $`(\nu _\tau \tau \overline{\nu }_\mu \overline{\mu })`$ for two input tau neutrino fluxes ($`F_\nu ^0E_\nu ^1`$ and $`E_\nu ^2`$), assumed to have been produced via neutrino oscillations from extragalactic sources, is evaluated to assess the impact of secondary neutrinos on upward muon rates in a km3 detector. We show that the secondary fluxes are considerably suppressed for the steeper flux, and even for fluxes $`E_\nu ^1`$, the secondary flux will be difficult to observe experimentally.

Evidence of neutrino oscillations from measurements of the atmospheric $`\nu _e`$ and $`\nu _\mu `$ fluxes leads one to the conclusion that $`\nu _\mu \nu _\tau `$ over distances characterized by the radius of the Earth for neutrino energies $`1`$ GeV . The results of a two-flavor analysis yield $`\mathrm{\Delta }m^2=2.5\times 10^3`$ eV2 and bi-maximal mixing, $`\mathrm{sin}^22\theta =1`$. Three flavor analyses are consistent with this result .

Atmospheric neutrinos come from cosmic ray interactions with air nuclei, yielding hadrons, especially mesons, which decay to neutrinos. Neutrinos may also be produced at the sources of cosmic rays, where energetic protons interact with nucleons and photons at the source. The pions and kaons produced at these distant astrophysical sources are parents of neutrino fluxes. Active galactic nuclei and gamma ray bursters are two proposed astrophysical neutrino sources . Other possible sources include exotic particle annihilations . Bi-maximal mixing, in the context of astrophysical neutrino sources, results in flavor oscillation from the source ratio of fluxes: $`\nu _e:\nu _\mu :\nu _\tau =1:2:0`$ to flux ratios at the Earth of $`\nu _e:\nu _\mu :\nu _\tau =1:1:1`$, independent of the neutrino energy, given that astronomical distances are so large .

Halzen and Saltzberg pointed out in Ref.  that high energy $`\nu _\tau `$ flux attenuation in the Earth differs from $`\nu _\mu `$ flux attenuation due to the fact that the $`\tau `$ produced in charged-current (CC) interactions with nucleons decays before it loses energy. For each $`\nu _\tau `$ lost in CC interactions, another $`\nu _\tau `$ appears following each $`\tau `$ decay, albeit at a lower energy. Detailed evaluations of $`\nu _\tau `$ flux attenuation in the Earth appear in Refs. . Depending on the incident flux, the $`\nu _\tau `$ flux shows a degree of “pile-up” as neutrinos of sufficiently high energy interact in the Earth and yield neutrinos at lower energies.

In a recent paper, Beacom, Crotty and Kolb  have suggested that in addition to a pile-up of tau neutrinos, the signal of astrophysical tau neutrinos will be enhanced by the appearance of “secondary” neutrinos. These secondary neutrinos come from purely leptonic decays of $`\tau `$’s. The idea is that while $`\nu _{\mathrm{}}`$ ($`\mathrm{}=e,\mu `$) fluxes starting, for example, at nadir angle 0 are extinguished for sufficiently high energies, they will be regenerated by the $`B=0.18`$ branching fraction for $`\tau \nu _\tau \mathrm{}\overline{\nu }_{\mathrm{}}`$, the $`\tau `$ being produced by $`\nu _\tau \tau `$ CC interactions. The flux of $`\overline{\nu }_{\mathrm{}}`$ is less likely to be extinguished due to the shorter path-length through the Earth (the $`\nu _\tau `$ already had to travel its “interaction distance,”) and its lower energy due to the combined energy loss in the CC process and the decay of the $`\tau `$.

The evaluations in Ref.  were for mono-energetic neutrinos. In this paper, we consider a few energy dependences for the $`\nu _\tau `$ fluxes. We show that the secondary fluxes are considerably suppressed for steeply falling fluxes and even for fluxes $`E_\nu ^1`$, the secondary flux will be difficult to observe experimentally.

Neutrino attenuation in the Earth is governed by a coupled set of partial differential equations. To illustrate, we write the coupled equations for the $`\nu _\tau `$ flux. For energy dependent flux $`F_{\nu _\tau }`$ (neutrinos/(cm2s sr GeV)),

$`{\displaystyle \frac{F_{\nu _\tau }(E,X)}{X}}={\displaystyle \frac{F_{\nu _\tau }(E,X)}{_\nu ^{int}(E)}}`$ (1)
$`+{\displaystyle _E^{\mathrm{}}}𝑑E_y[G^{\nu _\tau \nu _\tau }(E,E_y,X)+G^{\tau \nu _\tau }(E,E_y,X)]`$ (2)
$`{\displaystyle \frac{F_\tau (E,X)}{X}}={\displaystyle \frac{F_\tau (E,X)}{_\tau ^{int}}}{\displaystyle \frac{F_\tau (E,X)}{_\tau ^{dec}}}`$ (3)
$`+{\displaystyle _E^{\mathrm{}}}𝑑E_y[G^{\tau \tau }(E,E_y,X)+G^{\nu _\tau \tau }(E,E_y,X)].`$ (4)
Here, $`_\nu ^{int}=1/N_A\sigma _{\nu N}`$  and similarly for the $`\tau `$ interaction length, and $`_\tau ^{dec}=\gamma c\tau \rho `$ for density $`\rho `$ and Lorentz factor $`\gamma =E_\tau /m_\tau `$. The quantity $`X`$ is the column depth (in g/cm2), and for example,

$$G^{\nu _\tau \nu _\tau }(E,E_y,X)=\left[\frac{F_\nu (E_y,X)}{_\nu ^{int}}\right]\frac{dn^{NC}}{dE}(E_y,E).$$
(5)
The cross section normalized energy distribution of neutrinos with incident energy $`E_y`$ and final energy $`E`$ is represented by $`dn^{NC}/dE`$ in Eq. (2).

There are similar equations for $`\overline{\nu }_\tau `$ and $`\overline{\tau }`$ fluxes. In our discussion of flux ratios in the introduction, we didn’t distinguish between neutrinos and anti-neutrinos. Models predict that $`F_\nu =F_{\overline{\nu }}`$ to a good approximation. For our discussion in this section, we distinguish between particles and antiparticles since the secondary neutrinos are actually anti-neutrinos ($`\overline{\nu }_{\mathrm{}}`$) for incident $`\nu _\tau `$. One must keep in mind that there are equal incident neutrino and antineutrino fluxes of all three flavors in the context of neutrino oscillations. In the event rates evaluated below, we sum both neutrino and antineutrino contributions.

For nadir angle 0, $`_\nu ^{int}`$ equals the Earth’s diameter when $`E_\nu 40`$ TeV. Tau neutrino pile-up and associated secondary anti-neutrino production is relevant at this energy and higher. At very high energies, we have shown that the effect of attenuation for a surface detector viewing at zero nadir angle is significant even for $`\nu _\tau `$. For example, for a $`1/E_\nu `$ flux (with a smooth cutoff at $`E_\nu =10^8`$ GeV), the attenuated flux of $`\nu _\tau `$ at $`E_\nu =10^6`$ GeV is only 6% of the incident flux at that same energy . The $`\nu _e`$ and $`\nu _\mu `$ fluxes are about 1% of the incident flux at the same energy. Steeper fluxes have an even more significant attenuation and thus will have lower relative detection rates at high energies, whether it be $`\nu _e,\nu _\mu `$ or $`\nu _\tau `$ incident fluxes.

Guided by the falling fluxes and the increased attenuation, we confine our attention to the energy range of $`E_\nu =10^310^8`$ GeV. For energies below $`10^8`$ GeV, $`_\tau ^{int}>_\tau ^{dec}`$. Starting at $`E10^8`$ GeV, $`G^{\tau \tau }`$ becomes important , however, we neglect it here and confine our attention to lower energies. We also neglect the contribution of $`_\tau ^{int}`$ in Eq. (1). With the fluxes considered below, the error in this approximation should be small for rates evaluated with minimum energies of $`10^410^5`$ GeV. With these approximations, the set of differential equations is simplified and solved  using a modification of the iterative method detailed by Naumov and Perrone in Ref. . For $`\nu _\mu `$ fluxes, $`G^{\mu \mu }`$ from electromagnetic muon energy loss effectively eliminates any return of $`\nu _\mu `$ from CC interactions followed by $`\mu `$ decay in the energy range of interest.

The $`\tau `$ flux solution to Eq. (1) is responsible for generating the secondary neutrino flux. In the energy range of interest, one can write the differential equation for the $`\overline{\nu }_{\mathrm{}}`$ flux as:

$`{\displaystyle \frac{F_{\overline{\nu }_{\mathrm{}}}(E,X)}{X}}={\displaystyle \frac{F_{\overline{\nu }_{\mathrm{}}}(E,X)}{_{\overline{\nu }}^{int}(E)}}`$ (6)
$`+{\displaystyle _E^{\mathrm{}}}𝑑E_y[G^{\overline{\nu }_{\mathrm{}}\overline{\nu }_{\mathrm{}}}(E,E_y,X)+G^{\tau \overline{\nu }_{\mathrm{}}}(E,E_y,X)].`$ (7)
In what follows, we set $`G^{\nu \nu }=G^{\overline{\nu }\overline{\nu }}`$. At sufficiently high energies, neutrino and antineutrino interaction rates are equal because the cross sections are dominated by the sea quark distributions. The energy at which the interaction length equals the column depth increases with nadir angle, so the approximation is best at larger nadir angles. Because the pile-up comes from higher energy neutrino interactions, this approximation is not unreasonable, as discussed below.

A second approximation to obtain the secondary flux is to take:

$$G^{\tau \overline{\nu }_{\mathrm{}}}BG^{\tau \nu _\tau },$$
(9)
for $`B=0.18`$, the branching fraction for $`\tau \mathrm{}`$. In fact, the $`\overline{\nu }_{\mathrm{}}`$ spectrum from $`\tau `$ decay is a little softer than the $`\nu _\tau `$ spectrum, so this is an approximation that will slightly overestimate the secondary flux.

With these two approximations, one finds that the combination $`\mathrm{\Delta }(E,X)=F_{\nu _\tau }(E,X)F_{\overline{\nu }_{\mathrm{}}}(E,X)/B`$ satisfies the same transport equation as $`F_{\nu _\mu }(E,X)`$. At $`X=0`$, $`\mathrm{\Delta }(E,0)=F_{\nu _\tau }(E,0)=F_{\nu _\mu }(E,0)`$, so we can write

$$F_{\overline{\nu }_{\mathrm{}}}(E,X)B\left(F_{\nu _\tau }(E,X)F_{\nu _\mu }(E,X)\right).$$
(10)
Eq. (5) will be used in what follows to approximate the secondary antineutrino flux for two different incident spectra: $`E_\nu ^1`$ and $`E_\nu ^2`$.

In Fig. 1, we show the ratio of the attenuated flux to the incident flux at nadir angle $`\theta =0,30`$ and 60 degrees for

$$F_\nu ^0=F_\nu (E,X=0)=N_1/E1/(1+E/E_0)^2$$
(11)
where $`E_0=10^8`$ GeV and $`N_1`$ is a normalization factor, and for $`F_\nu ^01/E^2`$. The dashed curve is the $`\nu _\tau `$ flux, the dotted curve shows the attenuation of the other two neutrino species. The solid curve, the result of Eq. (5), is the secondary $`\overline{\nu }_{\mathrm{}}`$ flux, were $`\mathrm{}=e`$ or $`\mu `$. Except for the smallest nadir angles for the $`E_\nu ^1`$ spectrum, the secondary antineutrino flux is a small correction to the primary attenuated electron neutrino or muon neutrino flux. Even at nadir angle zero, the secondary flux is negligible compared to the transmitted primary flux for the $`1/E_\nu ^2`$ spectrum.

One should note that even though the $`\nu _\tau `$ flux dominates the primary and secondary $`\nu _\mu `$ fluxes, it does not dominate the contributions to the muon event rate because of the branching fraction $`B=0.18`$ of $`\tau \mu `$ together with the effect of energy loss as the $`\nu _\tau `$ converts to a $`\tau `$ which then decays to a $`\mu `$.

As a quantitative illustration of the implications of the secondary neutrino flux, we consider the muon event rate from the $`F_{\nu +\overline{\nu }}^0N_1/E_\nu `$ case (Eq. (6)) with $`N_1=10^{13}`$/(GeV cm2 s sr) for each neutrino flavor. The normalization factor $`N_1`$ is chosen to be in line with the Waxman and Bahcall gamma ray burster flux of Ref. , in which $`N_1=4\times 10^{13}`$ (in the same units) for the sum of all neutrino species and for $`E_\nu <10^5`$ GeV. The event rates shown below are for an underground detector of 1 km2 effective area, for example, the proposed IceCube detector . Details of the calculation appear in Ref. . Our integrals over neutrino energies were performed from the minimum neutrino energy $`E_\mu >10^4`$ or $`10^5`$ GeV, up to a maximum neutrino energy of $`10^8`$ GeV.

In Fig. 2, we include the contributions from $`\nu _\mu \mu `$ and $`\nu _\tau \tau \mu `$ (shown with the dashed lines) and in addition, the corresponding antineutrino induced antimuons from $`\nu _\tau \tau \overline{\nu }_\mu \overline{\mu }`$ (solid line). Anti-neutrino induced muons and antimuons are also included. We note that for the $`\nu _\mu \mu `$ rates, doing the neutrinos and antineutrinos separately results in at most a $`10\%`$

correction to the event rates at the energies shown here. This small correction is an indication that our approximations to obtain the secondary flux of neutrinos are not unreasonable.

The secondary neutrino contribution to the muon event rate has its largest relative contribution at nadir angle zero, with an enhancement over the $`\nu _\mu \mu `$ plus $`\nu _\tau \tau \mu `$ rate of 50-60% for the $`1/E_\nu `$ flux. At this angle, the ratio of rates of the secondary contribution to the muon event rate relative to the event rate of muons from tau decays is quite large, about 1.5 for $`E_\mu >10^4`$ GeV and 2.6 for $`E_\mu >10^5`$ GeV. Unfortunately, this is where the event rate is smallest and statistics are low. By a nadir angle of $`60^{}`$ (1 rad), where the event rate is roughly a factor of 10-20 larger, depending on the minimum muon energy, the enhancement in the overall muon rate is about $`25\%`$. At this nadir angle, the secondary $`\nu `$ produced muons are equal to the tau decay muon rate for $`E_\mu >10^5`$ GeV. The crossover occurs at $`\theta 0.7`$ rad for $`E_\mu >10^4`$ GeV. At the larger nadir angles, the $`\nu _\mu \mu `$ contribution to the muon event rate is dominant.

The normalization of the isotropic $`1/E_\nu `$ flux has been guided by the Waxman-Bahcall gamma ray burster flux of Ref. . With $`E_0=10^8`$ GeV in Eq. (6), our flux violates the Waxman-Bahcall bound  above $`E_\nu 10^6`$ GeV, so the event rates in Fig. 2 may be optimistic. If the normalization $`N_1`$ is reasonable, then the secondary contribution to the muon rate will be difficult to observe. Low statistics will make it hard to have a meaningful comparison between the small and large nadir angle rates. Compounding the problem is that one does not expect to know the input flux energy dependence or normalization exactly.

For more rapidly falling fluxes, the contribution of secondary neutrinos to the event rates is smaller. The attenuated fluxes shown in Fig. 1 for $`1/E_\nu ^2`$ yield secondary enhancements which are quite small. The secondary muon rate is only about 10-15% of the primary $`\nu _\mu \mu `$ rate at nadir angle zero. Typical theoretical neutrino fluxes have spectra that lie somewhere between the $`1/E_\nu `$ and $`1/E_\nu ^2`$ cases in this energy range .

In summary, the energy dependence of the incident tau neutrino flux is crucial in evaluations of the implications of $`\nu _\tau `$ interactions to regenerate $`\nu _\tau `$ and secondary $`\overline{\nu }_\mu `$ and $`\overline{\nu }_e`$. Neutrino fluxes are attenuated due to their passage through the Earth, even the tau neutrinos, thus moderating tau neutrino contributions to secondary neutrino fluxes. Our evaluation has relied on approximations including Eq. (4), setting $`_\tau ^{int}=0`$ and $`G^{\nu \nu }=G^{\overline{\nu }\overline{\nu }}`$. If muon rates are determined to be large, even near nadir angle zero, then a more detailed evaluation of the secondary flux may be in order. With our current theoretical expectations for flux normalizations, however, the secondary neutrinos coming from $`\tau `$ decays will be difficult to observe experimentally, as they contribute significantly to a muon excess only at small nadir angles where the fluxes are already strongly attenuated and for spectra like $`1/E_\nu `$.

We thank J. Beacom for discussions. The work of S.I.D. has been supported in part by NSF Grant 0070998. The work of I.S. has been supported in part by the DOE under contracts DE-FG02-95ER40906 and DE-FG03-93ER40792. The work of M.H.R. has been supported in part by NSF Grant No. PHY-9802403 and DOE under contract FG02-91ER40664. M.H.R. thanks the Fermilab Theory Group for its hospitality.