A phenomenological model of galaxy clusters

We present a simple model to describe the dark matter density, the gas density, and the gas temperature profiles of galaxy clusters. Analytical expressions for these quantities are given in terms of only five free parameters with a clear physical meaning: the mass M of the dark matter halo (or the characteristic temperature T_0), the characteristic scale radius a, the cooling radius in units of a (0<alpha<1), the central temperature in units of T_0 (0<t<1), and the asymptotic baryon fraction in units of the cosmic value (f~1). It is shown that our model is able to reproduce the three-dimensional density and temperature profiles inferred from X-ray observations of real clusters within a 20 per cent accuracy over most of the radial range. Some possible applications are briefly discussed.


INTRODUCTION
Observations of galaxy clusters provide invaluable information about the universe and its evolution. The abundance of galaxy clusters as a function of temperature and redshift, their baryon fraction, the X-ray luminosity, the Sunyaev-Zeldovich and lensing effects have often been used as cosmological probes (see e.g. Voit 2005, for a recent review). In particular, they constrain the matter density, the slope and normalization of the power spectrum of primordial fluctuations, and the equation of state of dark energy. Thanks to the advent of high-resolution X-ray observatories, we can now study the structure of the intracluster medium (ICM) with unprecedented accuracy. For the brightest objects, the gas density and temperature profiles can be directly inferred from the observed data, assuming approximate spherical symmetry. However, it is not uncommon that these observations have a limited field of view, and the profiles have to be fitted by some analytical model in order to extrapolate them outwards. In fainter systems, the use of analytical models makes possible a more accurate determination of the cluster parameters. In the most extreme case, the models actually play a crucial role in the very detection of the faintest cluster candidates. It is therefore important to have a simple analytical description of the ICM that captures the most relevant features of the density and temperature profiles of the X-ray emitting gas.
Historically, one of the most popular options is the so-called β-model (Cavaliere & Fusco-Femiano 1976), where ⋆ E-mail: yago@aip.de the gas density is given in terms of three free parameters (the normalization ρ0, a core radius rc, and the exponent β), ρgas(r) = ρ0 The β-model has been widely used over the years to describe the radial structure of the ICM. However, it is well known that equation (1) fails to provide a consistent fit over the whole radial range. Moreover, the gas is assumed to be isothermal, which is certainly not consistent with the observed temperature profiles. Finally, for a gas in hydrostatic equilibrium, the dark matter distribution underlying a βmodel would be which tends to a constant value at the centre. Such a 'cored' density profile is in strong disagreement with the results of numerical simulations, as well as with most recent observational estimates of the mass distribution in galaxy clusters. Polytropic models, where the density and temperature of the gas are related by the effective equation of state have been shown to provide a much better description of the ICM for n ∼ 5. Nevertheless, the temperature profile predicted by these models tends to increase steadily towards the cluster centre, whereas both numerical experiments and observations of real systems show that radiative cooling makes the temperature drop in the central regions of the ICM. On the other hand, the central gas density is often observed to become a power-law rather than a constant density core. This has motivated the introduction of more elaborate analytical models, including additional free parameters in order to describe a wider range of possible behaviours (see e.g. Vikhlinin et al. 2006, hereafter V06, and references therein). In particular, these authors propose a model, characterized by 17 independent free parameters, that accurately reproduces currently available X-ray data. According to the V06 model, the three-dimensional gas density profile would be described by while the gas temperature is modelled as Equations (4) and (5) have great freedom and can provide a good fit to the observed density and temperature profiles, both for the inner and outer regions of the ICM. The main disadvantage, though, is that their 17 parameters are strongly correlated, and thus there are many degeneracies in their best-fitting values. This would be a relatively minor problem (e.g. computational cost) if our only goal was to reproduce the observational data, but it becomes of critical importance when one attempts to extrapolate outside the observed region or when the number of photons from the object under investigation is too low to obtain reliable profiles. In these cases, robustness becomes more important than flexibility, and a model with fewer parameters is preferable in order not to over-fit the available data. A simple, robust model can be extremely helpful, for instance, in cosmological studies, where the observed number counts of galaxy clusters as a function of temperature or luminosity need to be connected with the underlying mass function. Such a model would also be of great interest for multi-wavelength analysis, where data with very different errors are combined in order to recover the three-dimensional structure of the object.
Here we present an analytical model that is able to reproduce the complex behaviour of the gas density and temperature profiles observed in real galaxy clusters by using only five free parameters, all of which have a clear physical meaning. The model and its parameters are described in detail in Section 2. We compare it with observational data in Section 3, where it is shown that our simple model fits all the observable X-ray properties of the ICM, while being consistent with our current knowledge of the structure of dark matter haloes. In Section 4, we consider three possible applications, namely the set up of initial conditions in numerical experiments, the construction of optimal filters for X-ray detection, and the combined analysis of X-ray and Sunyaev-Zel'dovich data. Conclusions are briefly summarized in Section 5.

Dark matter
Perhaps one of the best known results of cosmological Nbody simulations is that the radial density profiles of dark matter haloes can be well fitted by a relatively simple analytical function with very few parameters (e.g. Navarro et al. 1997). The precise shape of such function, particularly near the centre, is still a matter of heated debate (see e.g. Merritt et al. 2006, for a recent discussion), but there is general agreement in that it should be shallower than isothermal (ρ ∝ r −2 ) as r → 0 while significantly steeper as r → ∞.
In order to model the cluster's dark matter halo, we use a Hernquist (1990) density profile, where M denotes the total mass and a is a characteristic scale length. The cumulative mass inside radius r is given by and the gravitational potential is simply where G is Newton's constant. Analytical expressions for other quantities, such as the velocity dispersion, distribution function or projected surface density, can be found in Hernquist (1990). For the sake of simplicity, we will assume that equations (6) and (7) correspond to the total density and mass, respectively, i.e. the sum of the dark and baryonic components. The dark matter profiles can be trivially obtained by subtracting the contribution of the gas.

Polytropic equation of state
Non-radiative gasdynamical simulations show that, in the absence of additional physics, the ICM of relaxed clusters can be approximately described as a polytropic gas in hydrostatic equilibrium with the gravitational potential created by the dark matter (e.g. Ascasibar et al. 2003). Under these conditions, one can derive analytical expressions for both the gas temperature and density profiles, where T0 and ρ0 correspond to the central values, and n is the effective polytropic index. Hydrostatic equilibrium also imposes the mass-temperature relation where k is the Boltzmann constant, mp denotes the proton mass, and µ ≃ 0.6 is the molecular weight of the gas. In order to obtain a constant baryon fraction at large radii, we set n = 4. , and values of the best-fitting parameters are given in Table 1. Last panels show the fractional deviation between our model and V06 (horizontal dotted lines indicate ±0.05 in logarithmic scale, and a vertical dotted line is drawn at the Chandra detection radius R det ).

Cool core
The polytropic model roughly agrees with observational results (e.g. Markevitch et al. 1998), but the observed temperature profiles often feature a central drop that is not well reproduced by equation (9). Therefore, we introduce a modified temperature where 0 < t < 1 is a free parameter that measures the amount of central cooling with respect to the polytropic solution, and ac reflects the cooling radius below which the effect is important. Physically, one expects ac < a, so we consider the parametrization α = ac/a with 0 < α < 1. Dropping the polytropic assumption, we can substitute this expression into the hydrostatic equilibrium equation and compute the corresponding gas density, where the normalization ρ0 can be expressed in terms of the cosmic baryon fraction, with f ∼ 1.

COMPARISON WITH REAL CLUSTERS
In order to test whether our simple prescription provides an adequate model of galaxy clusters, we consider the sample of 13 low-redshift, relaxed objects studied by V06. Rather than fitting the raw observational data, we simply attempt to reproduce the three-dimensional density and temperature profiles, described by equations (4) and (5) with values of the 17 free parameters according to tables 2 and 3 of V06.
We account for the errors in the deprojected quantities by assuming a log-normal distribution with σρ = σT = 0.05. For each object, we fit the three-dimensional profiles within the radii 0.2Rmin and R det , defined in V06. We vary T0, t, a, α, and f in 30 logarithmic steps and compute the reduced chi-square as with and The gas density and temperature are evaluated at N = 30 points where the radius ri is also increased logarithmically.  We assume a cosmic baryon fraction Ω b /Ω dm = 0.133 and we relate the electron and proton number densities to the gas density as ne = np ≈ ρgas/mp. Values of the best-fitting parameters are given in Table 1, and the resulting gas density and temperature profiles are plotted for each object in Figure 1, where the last panels show the fractional deviations with respect to the V06 model. For all systems, the discrepancy is of the order of 10 − 20 per cent within the fitted range, consistent with our adopted estimate of the error bars, σρ = σT = 0.05 (see also the reduced χ 2 values in Table 1). Outside this range, the model tends to yield densities and temperatures at small radii that are systematically lower than those measured by V06. Although this might actually be correct for some systems (see the observational data points in V06), it will not be so in others. It is likely that structure on small scales (e.g. cold and shock fronts, the very presence of a central galaxy) breaks down the assumptions of perfect spherical symmetry and hydrostatic equilibrium, making any simple model inadequate to describe the ICM (and maybe even the dark matter potential) in the innermost part of the cluster, and indeed the objects that are worst described by our model display one or more inflection points in their profiles. At large radii, the density profile given by equation (13) seems to be steeper than the results of V06. In order to test the ability of our model to infer the cluster properties from data of poorer quality, restricted to a smaller field of view, we repeated our analysis considering only N = 3 radii between Rmin and 0.5R det . The accuracy of the recovered profiles was similar within the fitted region, but the extrapolation towards large radii was not entirely reliable, with errors of the order of a factor of two or above in the worst cases.
The shape of the likelihood function, or, more precisely, χ 2 (To, t, a, α, f ), is studied in Figure 2. For each cluster, a diamond marks the best-fitting parameters and a contour is drawn at χ 2 = 3. The underlying grayscale maps show the projected χ 2 , marginalized over all clusters and all possible values of the other parameters (i.e. the minimum χ 2 attained at each point by any object). Our results suggest that the likelihood function is well behaved, in the sense that it displays well-defined, unique maxima for each cluster, but it seems to be significantly skewed (i.e. the best fit does not always coincide with the geometrical centre of the contour) given our particular choice of variables. For each object, there is only a relatively mild degeneracy between the best-fitting values of the five free parameters of our model, although a certain mutual dependence is obvious in several pairs (for instance, those involving α, t and f ). Considering the whole sample, we find strong correlations (i.e. scaling relations) between most parameters and the mass of the cluster, or equivalently T0. Although more observational data and/or numerical simulations would be required in order to make a quantitative statistical assessment, it is interesting to note that these correlations could be exploited to reduce the number of free parameters in the model even further. If all parameters could indeed be expressed as a function of T0, our model would effectively have one single free parameter, specifying the overall scale of the cluster.

APPLICATIONS
Having a simple analytical model to describe the threedimensional distribution of gas, temperature and matter in galaxy clusters can be useful in many respects. On the one hand, it may serve as a base to make theoretical predictions (e.g. to study the effect of a given perturbation on the ICM of a relaxed cluster). On the other hand, it also helps to interpret observational data (e.g. to estimate the mass of the cluster, or the temperature profile, with a small number of photons). In this section we briefly discuss three examples of possible applications of a model like the one presented in this paper.

Numerical simulations
Our model can be used, for instance, to set up the initial conditions in idealized numerical experiments (and, in fact, it has already been used for this purpose in Ascasibar & Markevitch 2006). In a few words, the procedure to generate a synthetic cluster is as follows: 1 First, the radius of each particle is obtained by generating a uniform random number µ between 0 and 1 and inverting the appropriate (gas or dark matter) mass profile, i.e. solving for M (r) = µMx, where Mx denotes the corresponding total mass. Angular coordinates φ and cos θ are uniform random numbers in the range [0, 2π] and [−1, 1], respectively.
Gas temperature and density are given by equations (12) and (13) as a function of radius. When modelling a composite system containing several objects, the density of each particle is computed as ρ = ρ(ri), where ri is the distance to object i, and temperatures are set according to T = ρ(ri)T (ri)/ρ. Velocities are given by v = ρ(ri)v i /ρ, with v i being the centre of mass velocity of the i-th object.
For collisionless dark matter particles, each object is completely independent from the others. Velocities with respect to the relevant centre of mass are assigned from the probability distribution where Φ(r) denotes the gravitational potential, and the distribution function f (E) is computed by using Eddington's formula (Binney & Tremaine 1987), The probability distribution (18) is sampled by means of the von Neumann rejection algorithm (see e.g. Press et al. 1992). We generate a tentative velocity v, uniformly distributed between 0 and vmax = −2Φ(r), and an auxiliary random number p between 0 and v 2 max f (vmax). The velocity is accepted only if p < v 2 f (v). Otherwise, two new random numbers are generated for v and p until a value is finally accepted for the velocity of the particle. As for positions, the angular coordinates of the velocity are obtained from two uniform random variables, 0 < φ < 2π and −1 < cos θ < 1.

Optimal filtering
A major issue in the construction of galaxy cluster samples is the optimization of the detection algorithm. One of the most successful approaches is based on the use of wavelets (see e.g. Rosati et al. 1995;Lazzati et al. 1999). Simple bases, like the Mexican Hat wavelet, are optimal when the underlying signal is Gaussian and the background has a k −2 power spectrum, but neither condition is met by cluster Xray data. The underlying signal would be better described by our model, and the background can be modelled as a random Poisson variable whose normalization varies from pointing to pointing. In the case of the Sunyaev-Zel'dovich effect, some authors have improved the detection algorithms by adopting optimal filters instead of wavelets and assuming that the gas is described by a β-model (Herranz et al. 2002). As mentioned in the introduction, the β-model does not capture the complex behaviour of the gas density and 1 Computer code is available upon request. temperature profiles; moreover, it fails to provide a good description of the ICM in the outer regions, to which the Sunyaev-Zel'dovich effect is also sensitive. Therefore, an optimal filter based on our simple model could help to improve the detection rate in both the X-ray and millimeter bands.
For this purpose, simplicity is a major concern, since the shape and scale of the filter depend on the values of the free parameters. Our model provides an accuracy comparable to that of V06, at the expense of only one extra parameter with respect to the β-model. As discussed above, the correlations observed in Figure 2 suggest that it might be possible to define an optimal filter in terms of only one single scale parameter. Although exploring in detail the performance of such a filter is well beyond the scope of the present work, we show in Figure 3 an example of its application to cluster detection in X-ray data. A simulated cluster (based on one of the 13 models) has been added to real XMM data (see left panel of figure 3). After subtracting the brightest sources, we build an optimal filter based on the model of the cluster and the XMM background. The result is shown on the right panel of figure 3, where the simulated object can be seen in the center of the image, as well as another cluster candidate to its left.

Multi-wavelength 3D deprojection
There is already a number of galaxy clusters that have been observed in different wavebands, and this number is expected to increase dramatically during the forthcoming years. Since different observations probe the gas or matter distribution in different ways, a combined analysis provides extremely tight constraints on the three-dimensional properties of the cluster. For instance, the X-ray emissivity of the plasma, the magnitude of the thermal Sunyaev-Zel'dovich effect at millimeter bands and the convergence or shear of gravitational lensing (measured from optical data) are proportional to the integral along the line of sight of the gas density squared, the gas pressure, and the mass density, respectively. Taking into account all these observables, it is possible to reconstruct the three-dimensional structure of clusters with an extraordinary accuracy. Here we illustrate the problem with a simple example, involving two observations of the same object (an X-ray image and a Sunyaev- Zel'dovich image). Both datasets were simulated by adding uniform random noise to one realization of our model, and likelihoods on the α − a plane were derived in each case. For the sake of simplicity, the remaining three parameters have been fixed to the values used to generate the data. As can be seen in Figure 4, each data set is slightly biased with respect to the fiducial model (asterisk). When both results are combined together, the best-fitting α and a are much closer to their true values.
An interesting possibility, further along this direction, would be to recover the triaxial structure of the dark matter halo. Although the model we described here is spherically symmetric, it can be trivially modified to account for a triaxial gas distribution. The shape of the gravitational potential and the density of the dark matter halo can then be obtained from the hydrostatic equilibrium condition. Since the procedure involves second derivatives, it is extremely sensitive to the details of the gas distribution. It is important that the model captures these details accurately, while at the same time it must be simple enough to filter the high-frequency noise, which otherwise would dominate the final result.

CONCLUSIONS
We have presented an analytical model of galaxy clusters based on spherical symmetry and hydrostatic equilibrium. Our model is fully specified by equations (6), (11), (12), (13), and (14), and its five parameters have a well-defined physical interpretation: M is the mass of the system, a is a characteristic scale length, 0 < α < 1 is the cooling radius in units of a, 0 < t < 1 is the central temperature in units of T0, and f ∼ 1 is the asymptotic baryon fraction in units of the cosmic value.
Its main advantage with respect to previous work is the reduced number of free parameters. Besides a more straight-forward interpretation of the results, describing the observations in terms of a few parameters is also very convenient from a computational point of view, since (depending on the details of the algorithm) the complexity of the fitting procedure scales roughly exponentially with the dimensionality of the parameter space. An additional advantage are the smaller degeneracies, which make both the inferred profiles and the best-fitting values of the parameters more robust. This, in turn, makes possible to investigate correlations (scaling relations) that reduce the effective number of free parameters. Our results suggest that, in fact, it might be possible that the whole structure of a galaxy cluster could be completely determined by just one single characteristic scale. This very interesting possibility needs to be investigated further with more data. If confirmed, it would provide a crucial improvement to our current ability to model -and understand -the internal structure of galaxy clusters.
Comparing our model to the three-dimensional density and temperature profiles obtained from high-resolution Xray observations, we find that it provides an accurate description of the intracluster medium. In particular, the fractional deviations with respect to the more elaborate model proposed by V06 are comparable to the measurement errors (of the order of 10 − 20 per cent). The systems that feature high values of the reduced χ 2 usually display a complex behaviour that may be associated to departures from perfect hydrostatic equilibrium, but a much more detailed study would be required in order to test this hypothesis.
We have discussed a few possible applications of the model presented here. First, it provides a simple description of galaxy clusters to set up the initial conditions in idealized numerical experiments. Second, it can be helpful in the detection and characterization of observed systems, specially when the signal-to-noise ratio is not extremely large. Finally, a simple model of the matter and gas distribution makes possible to reconstruct the three-dimensional structure of the cluster from multiwavelength data. The construction of an optimal filter for cluster detection based on our model, as well as the combination of X-ray and Sunyaev-Zel'dovich information, will be the subject of future work.

ACKNOWLEDGMENTS
This work has been funded by the Spanish Ministerio de Educación y Ciencia, under project AYA2006-06266. YA would like to thank A. Vikhlinin for useful discussions, and the Instituto de Física de Cantabria for their kind hospitality. JMD benefits from a Ramón y Cajal contract from the Ministerio de Educación y Ciencia.