text:
# The radial structure of galaxy groups and clusters

## 1 Introduction

Clusters of galaxies are the largest gravitationally bound structures in the universe, and as such they have often been considered as a canonical data set for cosmological tests. As a first-order approximation for analytical studies, they can be described as spherically symmetric systems in which the baryonic gas is in hydrostatic equilibrium with an underlying cold dark matter (CDM) halo.

Even for this basic model, the radial distributions of both types of matter are far from being well understood. During the last two decades, a great effort has been devoted to investigate the mass distribution of CDM haloes by means of numerical N-body simulations. The systematic study undertaken by Navarro et al. (1997, hereafter NFW) concluded that the simulated density profiles can be fitted by a single two-parameter function, valid from galactic to cluster scales.

Similar results have been found in several independent papers, although there is still some debate about the innermost value of the logarithmic slope of the density profile and its dependence on resolution. For instance, Fukushige & Makino (1997, 2001), Moore et al. (1998, 1999, hereafter MQGSL) or Ghigna et al. (1998, 2000) find a central slope steeper than the NFW profile. On the other hand, Jing & Suto (2000) and Klypin et al. (2001) claim that the actual value might depend on halo mass, merging history and substructure. More recently, Power et al. (2003) point out that the logarithmic slope becomes increasingly shallower inwards, with little sign of approaching an asymptotic value at the resolved radii.

The situation is even less clear for the baryonic component. Most gas in the intracluster medium (ICM) is in the form of a hot diffuse X-ray emitting plasma, where the cooling time (except in the innermost regions) is typically longer than the age of the universe. Adiabatic gasdynamical simulations have therefore been used to study the formation and evolution of galaxy groups and clusters in different cosmologies (e.g. Navarro et al., 1995; Evrard et al., 1996; Bryan & Norman, 1998; Eke et al., 1998). This early work already showed that the gas distribution is also similar in all objects, but more extended than that of the dark matter. Temperature profiles are found to decrease significantly with radius, but an isothermal core (or even increasing temperatures) of variable extent have often been reported.

Recently, many numerical studies have focused on the effects of radiative cooling and non-gravitational heating on the final properties of groups and clusters of galaxies (see e.g. Bialek et al., 2001; Davé et al., 2002; Muanwong et al., 2002; Borgani et al., 2001; Tornatore et al., 2003). Entropy injection at high redshift and/or removal of low-entropy gas modify the radial distribution of gas at small radii, making the ICM even more extended than it would be in the purely gravitational case. Since cool systems are more sensitive to additional heating, non-adiabatic processes have often been invoked to explain the discrepancy between the observed $`L_\mathrm{X}T_\mathrm{X}`$ relation (e.g. Edge & Stewart, 1991) and the self-similar scaling (Kaiser, 1986), as well as the entropy excess detected in galaxy groups (e.g. Ponman et al., 1999).

The extent to which additional physics can influence the cluster structure outside the cooling radius is unfortunately still unclear (see e.g. Lewis et al., 2000; Pearce et al., 2000; Motl et al., 2003). In this work, we are concerned with the theoretical predictions that can be made about the radial structure of galaxy groups and clusters beyond self-similarity arguments, in an attempt to gain some understanding of the purely gravitational case. We will consider that the ICM gas is in thermally-supported hydrostatic equilibrium with the CDM halo, which in turn dominates the mass of the system,

$$\frac{1}{\rho _\mathrm{g}(r)}\frac{\mathrm{d}P(r)}{\mathrm{d}r}=\frac{GM(r)}{r^2}\frac{GM_{\mathrm{dm}}(r)}{r^2}$$
(1)
where $`P=(\rho _\mathrm{g}kT)/(\mu m_\mathrm{p})`$ denotes gas pressure, $`G`$ is Newton’s constant, and $`M`$ is the mass enclosed within radius $`r`$. $`\rho _\mathrm{g}`$ and $`T`$ are the gas density and temperature, $`k`$ is Boltzmann’s constant, $`m_\mathrm{p}`$ is the proton mass and $`\mu `$ is the mean molecular weight of an ionised plasma of primordial composition ($`\mu 0.6`$).

Several approximations are made in equation (1). The first one is spherical symmetry. Second, bulk and random motions within the gas are not included. And third, baryons can contribute a significant fraction of the total mass at large radii. These effects can lead to appreciable departures from equation (1) in systems that are indeed in hydrostatic equilibrium. On the other hand, clusters of galaxies are evolving systems, and such a condition must not necessarily be fulfilled, particularly by dynamically young objects.

Assuming a polytropic equation of state,

$$P(r)\rho _\mathrm{g}^\gamma (r)$$
(2)
a one-to-one correspondence between mass, gas density and temperature can be established (Makino et al., 1998). Note that isothermality is a particular case of a polytropic relation (i.e. $`\gamma =1`$). The physical origin of a polytropic equation of state is still a matter under investigation. A possible justification can be given in terms of the entropy profile, related to the mass accretion history of the system (Voit et al., 2003).

A different approach to the radial structure of galaxy groups and clusters is the the so-called $`\beta `$-model, pioneered by Cavaliere & Fusco-Femiano (1976) and widely used thereafter (see e.g. Rosati et al., 2002, for a recent review). This model assumes that the ICM gas follows a ’universal’ King-like density profile

$$\rho _\mathrm{g}(r)=\rho _0\left[1+(r/r_\mathrm{c})^2\right]^{3\beta /2}$$
(3)
and that its temperature is approximately constant. Then, the observed X-ray surface brightness is given by the expression

$$S_\mathrm{X}(\theta )=S_0\left[1+(\theta /r_\mathrm{c})^2\right]^{3\beta +1/2}$$
(4)
where $`\theta `$ is the projected radius from the observed X-ray centroid and $`S_0`$ is the central surface brightness.

The physical meaning of $`\beta `$ is the ratio between gas temperature and the line-of sight velocity dispersion of the CDM component. The difference between this ratio ($`\beta _{\mathrm{spec}}1`$) and the fits obtained by applying (4) to observed clusters of galaxies ($`\beta _{\mathrm{fit}}0.7`$) has been often referred to as the ’$`\beta `$-discrepancy’ (see e.g. Sarazin, 1986). None the less, it has long been known both from observations (e.g. Neumann & Arnaud, 1999) and simulations (e.g. Bartelmann & Steinmetz, 1996) that the $`\beta `$-model does not provide an optimal description of the gas density, and that the best value of $`\beta _{\mathrm{fit}}`$ tends to increase with the outermost radius used in the fit.

In this paper, we investigate the radial distribution of dark and baryonic matter by means of adiabatic gasdynamical simulations. We assess the validity of the assumptions of hydrostatic equilibrium and a polytropic equation of state, and compare the density and temperature profiles derived from these tenets with a sample of 15 numerical clusters of galaxies.

The accuracy of the $`\beta `$-model has been previously addressed in several numerical studies (e.g. Evrard et al., 1996; Bartelmann & Steinmetz, 1996), focusing mainly on the mass estimates based on the observed X-ray surface brightness. In the present work, we consider four analytical cluster models. Two of them rely on the existence of a ’universal’ CDM density profile, whereas the others assume a $`\beta `$-model for the gas distribution. First, we attempt to obtain a self-consistent fit to our numerical data with each one of these models, noting that baryonic and dark matter profiles are not independent. Then, we try to estimate several cluster profiles from the simulated X-ray surface brightness.

We would like to make special emphasis on the importance of numerical resolution and accurate entropy conservation. We will show that the inner structure of the ICM depends critically on these issues. More specifically, we claim that they are responsible for the systematic differences observed between Eulerian and Lagrangian integration techniques (see e.g. Frenk et al., 1999).

Section 2 describes our numerical experiments. Analytical cluster models are summarised in Section 3. Results are compared in Section 4, and the impact on the analysis and interpretation of current observations is discussed in Section 5. We sum up our main conclusions in Section 6.

## 2 Numerical experiments

We have carried out a series of high-resolution gasdynamical simulations of cluster formation in a flat CDM universe with non-vanishing cosmological constant ($`\mathrm{\Omega }_\mathrm{m}=0.3`$; $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$; $`h=0.7`$; $`\sigma _8=0.9`$; $`\mathrm{\Omega }_\mathrm{b}=0.02h^2`$). A thorough description of the numerical experiments can be found in Ascasibar (2003).

Simulations have been run with the parallel Tree-SPH code Gadget (Springel et al., 2001). We have used a new non-public version of the code, kindly provided by Volker Springel, in which the entropy of SPH particles is explicitly conserved (Springel & Hernquist, 2002). Similar experiments have been accomplished with the Adaptive Refinement Tree (ART) pure N-body code (Kravtsov et al., 1997) from the same initial conditions. In order to test the reliability of our results, one of the clusters has been re-simulated with the gasdynamical version of ART (Kravtsov et al., 2002).

### 2.1 Cluster sample

In a cubic volume of 80 $`h^1`$ Mpc on a side, an unconstrained realization of the power spectrum of density fluctuations corresponding to the $`\mathrm{\Lambda }`$CDM model was generated for a total of $`1024^3`$ Fourier modes. The density field was then resampled to a grid of $`128^3`$ particles, which were displaced from their Lagrangian positions according to the Zeldovich approximation up to $`z=49`$. Their evolution until the present epoch is traced by means of a pure N-body simulation with $`128^3`$ dark matter particles.

Our sample comprises a total number of 15 clusters of galaxies, selected from this preliminary low-resolution experiment. Each object has been re-simulated with higher resolution by means of the multiple mass technique (see Klypin et al., 2001, for details). First, we located the particles in a spherical region around the centre of mass of the $`128^3`$ counterpart at $`z=0`$. Mass resolution is then increased by using smaller masses in the Lagrangian volume depicted by these particles, including the additional small-scale waves from the $`\mathrm{\Lambda }`$CDM power spectrum in the new initial conditions.

We use 3 levels of mass refinement, reaching an effective resolution of $`512^3`$ CDM particles ($`2.96\times 10^8h^1`$ M). Gas has been added in the highest resolved area only. The total number of particles (dark+SPH) in this area is greater than $`2\times 10^6`$ for all clusters. The gravitational softening length was set to $`ϵ=25h^1`$ kpc, depending on number of particles within the virial radius (Power et al., 2003). The minimum smoothing length for SPH was fixed to the same value as $`ϵ`$. In total, 7 independent numerical experiments have been performed, running Gadget on a SGI Origin 3800 parallel computer at Ciemat (Spain), using 32 CPU simultaneously. The average computing time needed to run each simulation was $`8`$ CPU days ($`6\times 10^5`$ s). The same clusters have also been simulated with the N-body version of ART on the Hitachi SVR at the LRZ (Germany).

Table 1 displays the physical properties of our clusters at $`z=0`$. Objects have been sorted (and named) according to their virial mass at the present day. Clusters J2 and K2 are an exception to this rule, since they are two small groups falling into J1 and K1, respectively.

In order to study the effect of mergers and close encounters, we used the Bound Density Maxima galaxy finding algorithm (see e.g. Klypin et al., 1999; Colín et al., 1999). We label as major merger any cluster in which we are able to identify a companion structure inside $`R_{200}`$ whose mass is greater than 0.5 $`M_{200}`$; if the most massive companion is in the range $`[0.10.5M_{200}]`$, the object is classified as a minor merger; otherwise, we consider it a relaxed system in virial equilibrium. The results of this classification scheme are quoted in the second column of Table 1. Clusters named Ji and Ki are relatively close pairs, but they are separate enough ($`23`$ Mpc) not to be considered as mergers.

We find a remarkable amount of dynamical activity in our randomly selected sample. Only 6 objects have been catalogued as relaxed clusters, whereas 5 fall into the category of minor mergers and 4 have been classified as major merging systems. As pointed out by Gottlöber et al. (2001), we find that the typical merging rate in groups is much higher than in the more massive clusters.

Columns 3–7 display the bulk properties of each object. The subscript ’200’ refers to the overdensity with respect to the critical value, $`\rho _c2.8\times 10^{11}h^2`$ M Mpc-3. Although the density contrast predicted in the spherical collapse model is closer to 100 for our $`\mathrm{\Lambda }`$CDM cosmology, we chose this value for consistency with most observational studies. Since most of the mass is concentrated in the central regions, the difference between the two prescriptions is not large.

Throughout this paper, we assume that X-ray emission arises from bremsstrahlung radiation only. The total luminosity of a set of $`N`$ gas particles is therefore computed as

$$L_\mathrm{X}=\lambda \underset{i=1}{\overset{N}{}}m_i\rho _iT_i^{1/2}$$
(5)
where $`\lambda =1.68\times 10^{17}`$ erg s-1 M$`{}_{}{}^{2}{}_{}{}^{}`$ Mpc3 keV-1/2. For equal mass particles, the emission-weighted temperature is given by

$$T_\mathrm{X}=\frac{_{i=1}^N\rho _iT_i^{3/2}}{_{i=1}^N\rho _iT_i^{1/2}}$$
(6)
These formulae have been used to compute the corresponding entries in Table 1, as well as the analytical and numerical profiles of X-ray related quantities.

### 2.2 Comparison with other numerical techniques

Since the late 1980s a variety of techniques have been developed to simulate gasdynamics in a cosmological context. In part inspired by the success of the N-body scheme, the first gasdynamical techniques were based on a particle representation of Lagrangian gas elements using the smoothed particle hydrodynamics (SPH) technique (Lucy, 1977; Gingold & Monaghan, 1977). Soon thereafter, fixed-mesh Eulerian methods were adapted (Cen et al., 1990; Cen, 1992) and, more recently, Eulerian methods with sub-meshing (Bryan & Norman, 1995), deformable moving meshes (Gnedin, 1995; Pen, 1995, 1998) or adaptive mesh refinement (AMR, Bryan et al., 1995; Kravtsov et al., 2002) have been developed, as well as extensions of the SPH technique (Shapiro et al., 1996; Springel & Hernquist, 2002; Serna et al., 2003).

The Santa Barbara cluster comparison project (Frenk et al., 1999) attempted to assess the extent to which existing modelling techniques gave consistent results in a realistic astrophysical application. The formation of a galaxy cluster in a SCDM universe was simulated with 12 independent codes, seven of them based on the SPH scheme and five on Eulerian methods. The properties of the CDM component were encouragingly similar, most discrepancies arising from small differences in timing (there was a merging event at $`z0`$). Less agreement was found in the gas-related quantities, but usually still within $`1020`$ per cent. Only the X-ray luminosity showed a strong dependence on the resolution of the different codes, with a spread as high as a factor of 10.

One of the most remarkable findings was a systematic trend in the temperature profiles obtained for the inner regions ($`r100`$ kpc). Near the centre, SPH codes generate a flat (or even slightly declining inwards) temperature profile, while grid codes produce temperature profiles that are still rising at the resolution limit. The entropy profile in SPH codes decreases continously towards the centre, while grid codes develop an isentropic core at small radii.

Several authors (e.g. Hernquist, 1993; Nelson & Papaloizou, 1993, 1994; Serna et al., 1996; Serna et al., 2003; Springel & Hernquist, 2002) have pointed out that an important shortcoming of conventional SPH formulations is the poor conservation of entropy when a low number of particles is used. We have run several test simulations of the Santa Barbara cluster to study the effects of numerical resolution and entropy conservation. The standard implementation of Gadget (Springel et al., 2001) gives similar results as the SPH-based codes in Frenk et al. (1999), while the temperature and density profiles of the entropy-conserving version (Springel & Hernquist, 2002) are closer to the results of the Eulerian code by Bryan et al. (1995) up to the resolution limit (see Ascasibar, 2003, for details). The innermost radius $`r_{\mathrm{min}}`$ was defined according to the most restrictive of the following criteria:

1. 200 dark matter particles (Klypin et al., 2001).
2. 100 gas particles (Borgani et al., 2002).
3. 3 times the gravitational softening (Power et al., 2003).

Usually, the resolution limit is set by the last two conditions. Excessively low values of the smoothing length can lead to the formation of cold compact groups of SPH particles that decouple from the surrounding medium. Although Eulerian simulations are not affected by this problem, an unphysical temperature drop is often found in the cores of clusters simulated with SPH-based codes at low and medium resolutions (e.g. Eke et al., 1998; Frenk et al., 1999; Mathiesen & Evrard, 2001). When the smoothing length is properly set, a high number of particles is required to notice the effects of numerical entropy loss, because the entropy profile artificially flattens at $`rr_{\mathrm{min}}`$ due to the softened density and temperature (see e.g. Borgani et al., 2002).

The maximum resolution allowed by the Santa Barbara initial conditions is only $`256^3`$ particles, and the comparison between different runs is further complicated by the merger at $`z0`$. Therefore, we decided to test the gasdynamical integration scheme with the most massive of our clusters. The results of the Eulerian code ART for Cluster A (Nagai & Kravtsov, 2002) are compared with both implementations of Gadget in Figure 1.

As expected, CDM density profiles agree within a few percent. Moreover, we find that the mass distribution in the pure N-body ART runs is almost identical ($`\mathrm{\Delta }M/M10`$ per cent) to that obtained in Gadget, once rescaled to account for the slightly different values of $`\mathrm{\Omega }_{\mathrm{dm}}`$ due to the presence of baryons. Only clusters F, I, J and K display bigger discrepancies ($`40`$ per cent) related to offsets in timing (note that all these systems are major mergers).

Differences between both implementations of Gadget are mostly evident in the gas distribution at small radii ($`20h^1`$ kpc). The density profile is steeper in the standard formulation of SPH, whereas the temperature is systematically higher when entropy conservation is enforced. There is an excellent agreement between the gas density found in ART and the entropy-conserving version of Gadget, but we find a lower temperature in the latter throughout most of the cluster.

The shape of the temperature profile, though, is similar in the two codes, and noticeably different from the approximately isothermal structure found in the standard SPH version. This affects the gas entropy near the centre, which shows a ’floor’ value of 30-40 keV cm-2 at $`r20`$ $`h^1`$ kpc in both ART and entropy-conserving Gadget runs, in contrast with most SPH implementations (including the public version of Gadget), which tend to produce a decreasing entropy profile down to the resolution limit (see e.g. Frenk et al., 1999).

The agreement between the Eulerian code ART and the explicit entropy-conserving scheme by Springel & Hernquist (2002) provides encouraging support to this new SPH formulation. Furthermore, it constitutes a sound piece of evidence that the conventional scheme suffers from severe entropy losses, with important consequences on the shape of the profiles in the central regions of groups and clusters of galaxies. We conclude that a careful modelisation of adiabatic processes is an essential requisite prior to the inclusion of additional physics, such as cooling or preheating of the intergalactic medium at high redshift.

## 3 Analytical models

We compare our numerical data with self-consistent analytical cluster models based on the assumptions of hydrostatic equilibrium and a polytropic equation of state. As a first step, we check whether these approximations are consistent with our simulated sample of galaxy groups and clusters.

Once the reliability of the basic assumptions is tested, we focus on two different kinds of models: one is based on the existence of a ’universal’ dark matter density profile (e.g. NFW or MQGSL) and the other relies on a ’universal’ gas density (i.e. $`\beta `$-model). In total, we test four different descriptions of the radial structure of galaxy groups and clusters. Hereafter, we will refer to them as NFW, MQGSL, $`\beta `$-model (BM), and polytropic $`\beta `$-model (PBM).

### 3.1 Basic assumptions

The extent to which equations (1) and (2) hold for our numerical clusters is measured in Figure 2. On the top panel, we plot the ratio between the gravitational term $`GM/r^2`$ and the pressure gradient $`\rho _\mathrm{g}^1\mathrm{d}P/\mathrm{d}r`$. For a gas in thermally-supported hydrostatic equilibrium (i.e. neglecting infall, angular momentum and turbulence), both quantities should be equal to satisfy equation (1). Clusters that have been classified as relaxed on dynamical grounds can be considered in hydrostatic equilibrium with an accuracy better than $`10`$ per cent. For minor mergers, this assumption holds only marginally ($`25`$ per cent). Clusters undergoing a major merger display an even larger scatter, particularly at small radii.

Deviations from hydrostatic equilibrium arise mostly from the contribution of kinetic energy to the gravitational support of the ICM gas. In relaxed clusters, thermal energy dominates the total energy budget. Velocity dispersion increases substantially during merging events, and the gas needs some time to dissipate this extra kinetic energy into heat. Since relaxed systems constitute only $`40`$ per cent of our sample, some caution must be kept in mind whenever equation (1) is applied to a real cluster ensemble. We note as well that infall and random motions are important in relaxed clusters for $`r0.8R_{200}`$.

Concerning the equation of state, the most straightforward way to check whether the ICM follows a polytropic relation is to compute the effective value of the polytropic index from the spherically-averaged density and temperature profiles,

$$\gamma (r)=1+\frac{\mathrm{d}\mathrm{log}[T(r)]}{\mathrm{d}\mathrm{log}[\rho _\mathrm{g}(r)]}$$
(7)
This quantity is shown on the bottom panel of Figure 2. Albeit some scatter, gas in our relaxed clusters and minor mergers is consistent with a polytropic equation of state with $`\gamma 1.18`$ (dashed line). An isothermal profile ($`\gamma =1`$) can be confidently ruled out from these data. A density-temperature histogram of the gas particles of Cluster A is plotted in Figure 3. We see that a polytropic equation of state reflects not only the mean behaviour of the ICM gas, but also describes the actual ratio between the density and temperature of individual mass elements. Although only one of our clusters is shown in Figure 3, the plot looks similar for any other object.

Gas in small haloes is typically much colder than the intracluster medium. Cold gas clumps corresponding to infalling galaxies can be distinctly appreciated in Figure 3, but their mass is too low to alter the effective polytropic index of Cluster A. However, the profiles become very noisy beyond the virial radius, as well as in major mergers. A constant $`\gamma `$ is not an accurate approximation in these cases, where the amount of substructure within a spherical shell is large.

### 3.2 Cluster models based on a ’universal’ CDM density profile

As pointed out in Makino et al. (1998) and Suto et al. (1998), the gas density and temperature profiles can be derived analytically from the CDM mass distribution under the assumptions of hydrostatic equilibrium and a polytropic equation of state. In the present work, we will consider the functional forms proposed by NFW,

$$\rho _{\mathrm{dm}}^{\mathrm{NFW}}(r)=\frac{\rho _\mathrm{s}}{(r/r_\mathrm{s})(1+r/r_\mathrm{s})^2}$$
(8)
and MQGSL,

$$\rho _{\mathrm{dm}}^{\mathrm{MQGSL}}(r)=\frac{\rho _\mathrm{m}}{(r/r_\mathrm{m})^{3/2}\left[1+(r/r_\mathrm{m})^{3/2}\right]}$$
(9)
Gas density and temperature profiles can be computed by substituting these mass distributions into equation (1) and combining it with equation (2). Following the notation of Suto et al. (1998), the gas temperature is given by

$$\frac{T(r)}{T_0}=1B_\mathrm{i}f_\mathrm{i}(r/r_\mathrm{i})$$
(10)
where the subscript ’i’ can be either ’s’ or ’m’ for NFW and MQGSL models,

$$f_\mathrm{s}(x)=1\frac{\mathrm{ln}(1+x)}{x}$$
(11)
and111Note the typos in formulas (15) and (16) of Suto et al. (1998).

$`f_\mathrm{m}(x)`$ $`=`$ $`{\displaystyle \frac{\pi }{3\sqrt{3}}}+{\displaystyle \frac{2}{\sqrt{3}}}\mathrm{arctan}\left({\displaystyle \frac{2x^{1/2}1}{\sqrt{3}}}\right)`$ (12)
$`+{\displaystyle \frac{1}{3}}\mathrm{ln}\left({\displaystyle \frac{x+2x^{1/2}+1}{xx^{1/2}+1}}\right){\displaystyle \frac{2\mathrm{ln}(1+x^{3/2})}{3x}}`$
In both cases,

$$B_i=\frac{4\pi G\mu m_\mathrm{p}(\gamma 1)\rho _\mathrm{i}r_\mathrm{i}^2}{\gamma kT_0}$$
(13)
relates the central temperature $`T_0`$ to the underlying dark matter distribution. The gas density can be easily obtained from the polytropic relation

$$\frac{\rho _\mathrm{g}(r)}{\rho _0}=\left[\frac{T(r)}{T_0}\right]^{1/(\gamma 1)}$$
(14)
According to this simple picture, clusters of galaxies can be described as a function of 5 parameters. A characteristic density and radius define the properties of the CDM halo, whereas for the ICM we also need to specify the gas density and temperature at $`r=0`$, as well as the polytropic index $`\gamma `$.

However, equation (13) relates the values of these variables. The number of free parameters can be further reduced by imposing additional constraints. For instance, Komatsu & Seljak (2001) choose $`T_0`$ and $`\gamma `$ by enforcing approximately constant baryon fraction at large radii. In our case, we just impose vanishing density at infinity, i.e. $`B_\mathrm{s}=1`$ and $`B_\mathrm{m}=3\sqrt{3}/(4\pi )`$, and we compute the central gas density from the condition that the baryon fraction never exceeds the cosmic value,

$$\frac{\rho _\mathrm{g}(r_{\mathrm{max}})}{\rho _{\mathrm{dm}}(r_{\mathrm{max}})}=\frac{\mathrm{\Omega }_\mathrm{b}}{\mathrm{\Omega }_{\mathrm{dm}}}$$
(15)
where $`r_{\mathrm{max}}`$ is the radius at which the predicted baryon fraction reaches its maximum ($`2.71r_\mathrm{s}`$ or $`1.34r_\mathrm{m}`$). With this prescription,

$$\rho _0^{\mathrm{NFW}}1.51\frac{\mathrm{\Omega }_\mathrm{b}}{\mathrm{\Omega }_{\mathrm{dm}}}\rho _\mathrm{s},\rho _0^{\mathrm{MQGSL}}23.08\frac{\mathrm{\Omega }_\mathrm{b}}{\mathrm{\Omega }_{\mathrm{dm}}}\rho _\mathrm{m}$$
(16)
In principle, there is no reason why the baryon fraction cannot exceed the cosmic value somewhere within the cluster. Yet, our numerical results (see Section 4.3 below) indicate that the gas to dark matter ratio increases monotonically (i.e. $`r_{\mathrm{max}}=\mathrm{}`$) up to $`\mathrm{\Omega }_\mathrm{b}/\mathrm{\Omega }_{\mathrm{dm}}`$. Since the analytical estimate obtained for NFW and MQGSL models does not tend to an asymptotic value, we decided to set the normalisation at the maximum. The value of $`\rho _0`$ varies by $`25`$ per cent for any choice between $`r_{\mathrm{max}}`$ and $`R_{200}`$, but inaccuracies in the assumption of thermally-supported hydrostatic equilibrium might bias a prescription based on large radii.

The last parameter in this family of models is the polytropic index $`\gamma `$. According to our results (see Figure 2), we set it to $`\gamma =1.18`$. Although this value is consistent with recent observations (e.g. Markevitch et al., 1998; Sanderson et al., 2003), it would be interesting to investigate its physical origin, as well as whether it can be related to other parameters such as halo mass or concentration. However, such a study requires a much larger number of objects, over a broader mass range, in order to be statistically significant.

### 3.3 $`\beta `$-models

Most observational studies do not explicitly rely on the assumption of a given CDM density profile. Instead, they fit the observed X-ray surface brightness with equation (4), and then compute the dark and baryonic profiles by applying hydrostatic equilibrium and a polytropic (often isothermal) equation of state.

To a great extent, the so-called ’$`\beta `$-discrepancy’ is due to the fact that expression (3) does not provide a good fit to the gas distribution throughout the whole cluster (see e.g. Navarro et al., 1995; Bartelmann & Steinmetz, 1996). Mathiesen & Evrard (2001) argue that line emission from cold gas in the outskirts of the cluster is entirely responsible for the observed value $`\beta 2/3`$.

Polytropic $`\beta `$-models have also been proposed to account for the presence of a temperature gradient, but equation (3) is still used to model the gas density (e.g. Ettori, 2000). In the general polytropic case, the surface brightness

$$S_\mathrm{X}(\theta )=S_0\left[1+(\theta /r_\mathrm{c})^2\right]^{1/23\beta (1+\frac{\gamma 1}{4})}$$
(17)
is similar to the conventional form (4) for any reasonable value of the polytropic index ($`1<\gamma <5/3`$). Therfore, the gas density inferred from a fit to the observed X-ray profile has virtually the same shape in both the isothermal (BM) and polytropic (PBM) $`\beta `$-models. The central surface brightness is given by

$$S_0=\sqrt{\pi }\rho _0^2r_\mathrm{c}\lambda T_0^{1/2}\frac{\mathrm{\Gamma }[3\beta (1+\frac{\gamma 1}{4})1/2]}{\mathrm{\Gamma }[3\beta (1+\frac{\gamma 1}{4})]}$$
(18)
Hydrostatic equilibrium leads to an underlying mass distribution

$$M(r)=\frac{3\beta kT_0\gamma r_\mathrm{c}}{G\mu m_\mathrm{p}}(r/r_\mathrm{c})^3\left[1+(r/r_\mathrm{c})^2\right]^{3\beta (\gamma 1)/21}$$
(19)
which features a finite central density

$$\rho _{\mathrm{dm}}(0)\rho _\beta =\frac{9\beta kT_0\gamma }{4\pi G\mu m_\mathrm{p}r_\mathrm{c}^2}$$
(20)
Therefore, the structure of the dark matter halo in the $`\beta `$-model is related to gas temperature in a way analogous to the models described above. We can follow the same procedure to normalise the central gas density, imposing that the maximum baryon fraction equals the cosmic value. We will consider a standard $`\beta `$-model ($`\gamma =1`$) with $`\beta =2/3`$ and a polytropic form with $`\gamma =1.18`$ and $`\beta =1`$. For these values, the maximum baryon fraction occurs at $`r_{\mathrm{max}}^{\mathrm{BM}}\mathrm{}`$ and $`r_{\mathrm{max}}^{\mathrm{PBM}}4.18r_\mathrm{c}`$, respectively, and hence

$$\rho _0^{\mathrm{BM}}=\frac{1}{3}\frac{\mathrm{\Omega }_\mathrm{b}}{\mathrm{\Omega }_{\mathrm{dm}}}\rho _\beta ,\rho _0^{\mathrm{PBM}}0.39\frac{\mathrm{\Omega }_\mathrm{b}}{\mathrm{\Omega }_{\mathrm{dm}}}\rho _\beta $$
(21)
We quote in Table 2 the distributions of dark and baryonic matter expected in our four analytical cluster models. Cumulative gas mass and projected quantities in NFW and MQGSL must be integrated numerically. Relations between the characteristic parameters are summarised in Table 3.

## 4 Results

We have shown that our clusters are, to a fair extent, in hydrostatic equilibrium, and that a polytropic equation of state with $`\gamma =1.18`$ provides a good description of the hot diffuse component of the ICM. Under these conditions, models of the radial structure of galaxy groups and clusters can be derived from phenomenological approaches to either the gas or dark matter density. In this Section, we compare these simple, self-consistent models with the spherically-averaged distributions of mass, gas density and temperature found in our sample of simulated galaxy clusters.

For NFW and MQGSL, we fit the cumulative dark mass of each halo, while for the $`\beta `$-models the gas mass is used instead. We compute the numerical profiles in 26 logarithmic bins between $`0.05R_{200}`$ and $`R_{200}`$. We set the lower cut-off well above our resolution limit ($`0.01R_{200}`$) not only to avoid numerical effects but also to investigate how accurate are the extrapolations of the analytical cluster models towards $`r=0`$.

The parameter space is explored by letting $`r_{200}`$ and $`c_\mathrm{i}r_{200}/r_\mathrm{i}`$ vary uniformly in the range $`0.25<r_{200}<1.25h^1`$ Mpc and $`0.1<c_\mathrm{i}<16`$ (50 for BM). For the gas mass, the radius $`r_{10}`$ (enclosing a mean density of $`10\rho _\mathrm{c}`$) is used instead of $`r_{200}`$. Both quantities are roughly equivalent for the baryon fraction assumed in this work.

Results of the $`\chi ^2`$ minimisation are plotted in Figure 4. In most cases, the quality of the fit is only slightly better for NFW and MQGSL models ($`\sqrt{\chi ^2/(dof)}0.050.1`$) than for $`\beta `$-models ($`0.050.15`$). Best-fitting values of $`r_{200}`$ are close to the actual $`R_{200}`$ (quoted in Table 1) obtained directly from the simulated mass profile. The values of $`r_{10}`$ in BM ($`\beta =2/3`$) are less well correlated with $`R_{200}`$ than those obtained for PBM ($`\beta =1`$).

Although the full sample is consistent with constant dark matter concentration over the restricted mass range probed, we notice a well-defined bimodal behaviour: while relaxed clusters follow the usual relation, approximated by $`c=186M_{200}^{0.1}`$ for the NFW profile (Burkert & Silk, 1999), mergers are dramatically offset down, particularly low-mass systems. In the $`\beta `$-models, we find that the best-fitting gas concentration even *increases* with cluster mass. In several systems, the core radius inferred by the BM is only slightly larger than the innermost radial bin included in the fit.

An important fact is that dark and baryonic matter distributions are expected to be strongly correlated, by virtue of the hydrostatic equilibrium equation and the constraint of a cosmic baryon fraction. All analytical models of galaxy groups and clusters studied in the present work have only two free parameters. According to our fitting procedure, gas density and temperature are genuine predictions (i.e. not fits) of the models based on a dark matter profile. Conversely, the same is true for the CDM density and gas temperature in the $`\beta `$-models.

### 4.1 Mass distribution

The radial density profiles of our dark matter haloes are shown on the left panel of Figure 5, rescaled by their best-fitting characteristic densities and radii. We plot the expected mass distribution in each model (see Table 2) as a solid line. On the right panel, we plot the accuracy of the analytical estimates, defined as

$$\frac{\mathrm{\Delta }\rho _{\mathrm{dm}}}{\rho _{\mathrm{dm}}}\frac{\rho _{\mathrm{est}}(r)\rho _{\mathrm{num}}(r)}{\rho _{\mathrm{num}}(r)}$$
(22)
where $`\rho _{\mathrm{est}}`$ and $`\rho _{\mathrm{num}}`$ are the estimated and numerical profiles, respectively.

Not surprisingly, both NFW and MQGSL offer an excellent fit to the dark matter distribution, but the core-like predictions of the $`\beta `$-models represent a very poor approximation to the numerical profiles. The conventional $`\beta `$-model has a much better accuracy throughout the fitted range, but the extrapolation to larger radii is more reliable in the polytropic version, thanks to its higher value of $`\beta `$. 222In the isothermal case, though, dark matter density is insensitive to the actual value of this parameter (see Table 2).

Concerning the controversy about the inner slope of the density profile, it is apparent from Figure 5 that MQGSL formula is slightly more accurate than NFW near the centre, albeit the scatter of individual profiles around the average is somewhat higher. This issue is addressed in more detail in Figure 6, where we plot the logarithmic slopes of the CDM density and cumulative mass profiles, computed directly from the simulation data

$$\alpha _\rho =\frac{\mathrm{d}\mathrm{log}[\rho _{\mathrm{dm}}(r)]}{\mathrm{d}\mathrm{log}(r)},\alpha _\mathrm{M}=\frac{\mathrm{d}\mathrm{log}[M_{\mathrm{dm}}(r)]}{\mathrm{d}\mathrm{log}(r)}$$
(23)
In general terms, both NFW and MQGSL formulae are consistent with our results up to the resolution limit. When we divide our sample according to the dynamical state of each cluster, we find a well defined trend in the sense that relaxed haloes display shallower central slopes (close to NFW) than merging systems (better described by MQGSL). As an extreme case, major mergers feature pure power-law mass distributions for more than one decade in radius, whose slopes are even steeper than the value proposed by MQGSL. In agreement with Power et al. (2003), we do not find any sign of an asymptotic slope for either relaxed clusters or minor mergers. More resolution is clearly required before achieving a firm conclusion on this matter.

### 4.2 Gas temperature

Mass determinations from X-ray emission in clusters usually assume that relaxed clusters (i.e. morphologically symmetric) are isothermal. However, this assumption has been questioned recently by Markevitch et al. (1998), who found evidence of a decreasing temperature profile for nearby clusters observed with ASCA. This result has been confirmed for cool clusters by Finoguenov et al. (2001). On the contrary, White (2000) and Irwin & Bregman (2000), using data from BeppoSAX and ASCA satellites, do not find any decrease of the temperature in a large collection of clusters. More recently, the analysis of BeppoSAX observations accomplished by De Grandi & Molendi (2002) concluded that temperature profiles of galaxy clusters can be described by an isothermal core followed by a rapid decrease.

Numerical simulations of galaxy clusters support a significant decrease in the ICM temperature from the central regions to the virial radius (e.g. Frenk et al., 1999). Recent results from high resolution AMR gasdynamical simulations seem to indicate that the temperature profile has an universal form (Loken et al., 2002). These authors proposed a simple formula to fit the projected emission-weighted temperature of their clusters:

$$T_\mathrm{p}^\mathrm{X}(\theta )=T_0(1+\theta /a_\mathrm{x})^\delta $$
(24)
where $`T_0`$ is the central temperature and $`a_\mathrm{x}`$ a core radius. Loken et al. (2002) quote best-fitting values of these parameters $`a_\mathrm{x}=r_{\mathrm{vir}}/1.5`$ and $`\delta =1.6`$.

The average mass-weighted temperature profile of our sample of galaxy clusters (excluding major mergers) is shown on the left panel of Figure 7. The normalisation $`T_0`$ has been computed analytically from the fit to the dark matter (NFW, MQGSL) or gas (BM, PBM) distribution of each cluster. As can be seen on the right panel, the theoretical prediction is accurate within $`15`$ per cent for all prescriptions, with the obvious exception of the isothermal $`\beta `$-model. NFW and MQGSL models are virtually indistinguishable.

In order to compare with observations, as well as with the numerical work of Loken et al. (2002), we plot in Figure 8 the average projected X-ray temperature profile. We computed the emission-weighted temperature according to equation (6) in cylindrical shells of 3 $`h^1`$ Mpc length, oriented along the main axes of the simulation box. The final profile of each cluster is given by the average over the three orthogonal projections, and the normalisation comes from the emission-weighted temperature, $`T_{200}^X`$, within $`R_{200}`$ (both quantities are quoted in Table 1).

Our results are in excellent agreement with the observational estimate of Markevitch et al. (1998), who claim that the emission-weighted temperature is well represented by a polytropic $`\beta `$-model (indeed, our data are almost identical to their best-fitting profile). However, our profiles are only marginally consistent with the results reported by De Grandi & Molendi (2002). In particular, we do not find any evidence of the temperature flattening observed by these authors at $`r0.2R_{200}`$.

The projected X-ray temperature of our clusters rises uninterruptedly until the innermost radial bin, in agreement with the results of Eulerian codes, such as ART (see Section 2.2) or that employed by Loken et al. (2002). Equation (24) fits well our numerical data, but we obtain rather different values of the core radius and the asymptotic exponent, $`a_\mathrm{x}=r_{200}/4.5`$ and $`\delta =0.7`$. Since both simulated temperature profiles are consistent within the error bars, it seems that the best-fitting values of $`a_\mathrm{x}`$ and $`\delta `$ are rather sensitive to the details of the data.

Never the less, we would like to stress that a decreasing temperature profile is a robust result of numerical simulations regardless of the integration scheme (either AMR or SPH), as long as entropy conservation is enforced in the latter. Furthermore, such behaviour is strongly supported by most recent observational measurements. The presence of an isothermal core in the observed galaxy clusters could reflect the need for additional non-adiabatic physics in the simulations.

### 4.3 Gas density

Finally, we compare the simulated gas density profile with the distributions expected in our four analytical cluster models. It is evident from Figure 9 that the discrepancies are remarkably larger than for the dark matter density or ICM temperature profiles.

It is somewhat surprising that the $`\beta `$-models show the worst inconsistencies with the data, since they are specifically intended to fit the gas distribution. When $`\beta =2/3`$, the asymptotic gas density falls too slowly, as $`r^2`$. However, the central regions cannot be properly fitted if $`\beta =1`$ is assumed. The final outcome is that the best-fitting value of $`\beta `$ increases from one value to the other as larger radii are considered in the fit.

Models based on a ’universal’ dark matter profile provide a fair description of the ICM density. None the less, we caution that in many cases the simulated gas density can differ substantially from the analytical prediction, giving rise to the huge scatter ($`40`$ per cent) observed in Figure 9. On one hand, the CDM distribution at small radii depends on the dynamical state of the halo. Near the virial radius, even relaxed clusters deviate noticeably from thermally-supported hydrostatic equilibrium (see Figure 2), and we also observe that the dark matter density is slightly higher than indicated by both NFW and MQGSL formulae (see Figure 5). More importantly, the polytropic index is not exactly the same for all clusters, which strongly alters the shape of the inferred gas density profile.

Inaccuracies in the estimated dark matter and/or gas profiles lead to a significant mismatch between the analytical and numerical baryon fractions. The local ratio of gas to CDM density is plotted in Figure 10 as a function of radius, rescaled by the cosmic value. None of the analytical cluster models is able to fit satisfactorily the simulated profile. The canonical $`\beta `$-model is the only scenario in which the baryon fraction reaches an asymptotic value, but this is due to the equally wrong gas and dark matter densities. In the other models, the baryon fraction reaches a maximum around three times the characteristic radius, but then it drops rather steeply.

We used this maximum to find out the normalisation of the ICM density, $`\rho _0`$. However, the steep asymptotic slope of the gas density profile hints that more realistic models have to be considered in order to reproduce the radial structure near the virial radius. On the other hand, such a modelisation (e.g. including departures from spherical symmetry, infall and turbulent motions) would introduce additional free parameters that would complicate the physical interpretation of the results.

## 5 Observable consequences

Analytical models of galaxy groups and clusters are a key ingredient in the interpretation of observational data, since they relate the X-ray emission of the ICM to the underlying gas and dark matter distributions. In particular, X-ray observations are usually restricted to the central regions, and these models allow the extrapolation of the spherically-averaged profiles up to the virial radius. In some occasions, they are also extrapolated inwards, in order to correct for the presence of a central cooling flow. Quite often, the virial radius itself is derived by assuming a particular model.

Therefore, we would like to investigate in some detail to what extent the assumed radial structure can bias the conclusions drawn from observational studies. We construct X-ray surface brightness profiles in the same manner as we did for the projected emission-weighted temperature (i.e. assuming pure bremsstrahlung radiation and averaging over the three axes). Then, we fit the simulated profiles with our four analytical models, varying the projected characteristic radius, $`\theta _\mathrm{i}`$, and the central surface brightness, $`S_0`$. Results are shown in Figure 11.

Only data between 0.05 and 0.5 $`R_{200}`$ were used. The lower cut-off is similar to region excised in observational studies to avoid the cooling flow region. The upper scale is similar to that attained by X-ray data. Furthermore, we find that the emission at larger radii is significantly enhanced by small mass subhaloes, which manifest as peaks in the surface brightness. We discard these features by fitting only those bins that lead to a monotonically decreasing profile. As pointed out by Mathiesen & Evrard (2001), the inclusion of emission lines would still increase the contribution of substructure to the total X-ray luminosity.

The dependence of the fit on the choice of the outer radius $`R_{\mathrm{out}}`$ is plotted in Figure 12. NFW and MQGSL models are more stable than $`\beta `$-models, and their best-fitting $`\theta _\mathrm{i}`$ are much closer to the actual characteristic radius $`r_\mathrm{i}`$, computed from the mass distribution between 0.05 $`R_{200}`$ and $`R_{200}`$ (see Section 4). The conventional $`\beta `$-model offers a poor fit to our simulated surface brightness at large radii, due to the shallow slope associated to $`\beta =2/3`$. For several objects, the best fit is obtained with the lowest possible $`\theta _\mathrm{c}`$ (i.e. a pure power law, $`S_\mathrm{X}r^{1.5}`$). Since we assumed $`\beta =1`$ for the PBM, a better agreement is found with the numerical results, although the best-fitting $`\theta _\mathrm{c}`$ is systematically higher than the core radius inferred from the gas mass.

We assume that the characteristic radii of our clusters in each model are equal to their best-fitting projected equivalents, $`\theta _\mathrm{i}`$. To obtain the characteristic CDM densities, as well as the central gas densities and temperatures, we must solve a system of equations, which in the $`\beta `$-model is given by (18), (20) and (21). For the models based on a ’universal’ dark matter density profile, the surface brightness is integrated numerically, yielding the central values quoted in Table 3. They must be combined with equations (13) and (16), substituting the appropriate value of $`B_\mathrm{i}`$.

A last interesting point to note is that we are not using any information about the X-ray temperature. In doing so, we could skip one of our constraints, such as the value of the polytropic index or the normalisation of the central gas density with respect to the CDM component.

### 5.1 X-ray surface brightness

The first question we would like to address is how well can our analytical cluster models fit the X-ray surface brightness. We show in Figure 13 that the numerical profiles are fairly described by NFW or MQGSL, although the accuracy of the fits degrades dramatically at large radii. Never the less, we find that the difference between the analytical and simulated profiles is of the order of $`3040`$ per cent up to the virial radius. Some part of it is caused by the low gas density predicted by these models at large radii. The rest is due to the emission of substructure beyond 0.5 $`R_{200}`$.

Consistently with observational results, the isothermal $`\beta `$-model with $`\beta =2/3`$ provides a fairly good fit to the X-ray surface brightness only for $`r0.3R_{200}`$. As was shown for the gas density, the outer parts are better approximated by $`\beta 1`$. When we let $`\beta `$ vary as a free parameter, the best fit to our data is obtained by $`\beta 0.8`$. However, this is a ’compromise solution’ that fails to properly describe the ICM at very small or large radii.

### 5.2 Mass estimates

One of the most important applications of hydrostatic equilibrium and a polytropic equation of state is the estimate of the cumulative mass profile of both dark and baryonic components, given the observed X-ray emission. We plot in Figure 14 the CDM mass distribution inferred for our cluster sample (excluding major mergers), scaled by the characteristic densities and radii derived from the fit to the surface brightness profile.

The accuracy of the mass estimate is similar in all models ($`20`$ per cent errors), but the shape of the mass distribution is substantially more accurate when a NFW or MQGSL ’universal’ CDM density profile is assumed. The isothermal $`\beta `$-model gives an acceptable estimate of the dark matter mass for $`0.1<r/R_{200}<1`$ because of the small values inferred for the core radius. In the PBM, this quantity is substantially higher and thus a ’core’ in the dark matter distribution can be noticed for $`r0.3R_{200}`$.

Cumulative gas mass estimates are plotted in Figure 15. In this case, the small values of $`r_\mathrm{c}`$ in the BM lead to extremely high central densities, whereas the PBM (with $`\beta =1`$) predicts an excessively large core. Again, a profile closer to the numerical results can be obtained by treating $`\beta `$ as a free parameter. The estimates based on a ’universal’ CDM profile are only slightly more accurate. The scatter, as well as the overall shape, are greatly improved by considering individual values of the polytropic index.

Regarding the total baryon fraction, we can see in Figure 16 that the $`\beta `$-models can lead to severe overestimates of this quantity at small radii ($`r0.10.3R_{200}`$). Near the virial radius, though, uncertainties are of the order of 10 per cent. NFW and MQGSL models are accurate within $`20`$ per cent for $`r0.2R_{200}`$ and $`50`$ percent as we move closer to the centre.

Another (sometimes more important) source of uncertainty comes from the estimation of the virial radius itself. Furthermore, $`R_{200}`$ is not a fixed multiple of the characteristic radius in any analytical model of galaxy clusters, and hence the cumulative baryon fraction at a fixed overdensity is extremely sensitive to halo concentration, particularly at low radii. Assuming a NFW model, the baryon fraction within $`R_{200}`$ can change by less than 10 per cent from $`c=4`$ to $`c=9`$, but at 0.1 $`R_{200}`$, the ratio of gas to dark matter mass varies from .5 to .75 times the cosmic value.

### 5.3 Entropy

Entropy plays a fundamental role in clusters of galaxies because convection acts as an entropy-sorting device, moving low entropy material to the cluster core and high entropy material to the outskirts. Gas density and temperature in hydrostatic and convective equilibrium are just manifestations of the underlying entropy distribution. The claim of an entropy floor in small groups (Ponman et al., 1999) has often been regarded as evidence of preheating of the ICM at high redshift. Using a much larger sample, Ponman et al. (2003) find that the observed profiles appear to be approximately self-similar apart from a normalisation constant. They obtain a best-fit behaviour $`S(0.1R_{200})T^{0.65}`$ for systems spanning a temperature range from 1 to 10 keV.

The entropy profile has recently been studied by Borgani et al. (2001) and Finoguenov et al. (2003) by means of numerical simulations that include cooling and preheating. With gravitational heating alone, these authors find a power-law entropy profile, consistent with other results based on the standard implementation of SPH (see e.g. Frenk et al., 1999). As was discussed in Section 2.2, the situation changes dramatically when entropy conservation is explicitly enforced, either in the SPH algorithm or in Eulerian schemes. Shock heating becomes substantially more efficient, leading to a much higher entropy near the cluster centre, as well as stripping low-entropy gas from the infalling galaxies that otherwise would survive for several crossing times (Borgani et al., 2001).

We plot our entropy profiles in Figure 17. As expected for a polytropic gas in hydrostatic equilibrium, they are not pure power laws for any reasonable choice of the ’universal’ gas or dark matter distribution. Furthermore, the shape of the entropy profile does not depend systematically on the mass or temperature of the object, in agreement with recent observations (Ponman et al., 2003). NFW or MQGSL models provide a much better estimate of the entropy distribution than any version of the $`\beta `$-model. However, the low gas densities predicted at large radii yield a slope at $`rR_{200}`$ significantly steeper than the numerical data.

Given our limited temperature coverage, the normalisation of the entropy profiles is consistent within the error bars with both the self-similar scaling ($`ST`$) as well as with the observed trend, $`ST^{0.65}`$. In any case, we would like to stress that we do not expect to find the self-similar scaling albeit we are neglecting non-gravitational processes. The central entropy in *any* model scales as $`T_0\rho _0^{2/3}`$, where $`\rho _0\rho _\mathrm{i}`$ and $`T_0\rho _\mathrm{i}r_\mathrm{i}^2`$. Recovering a normalisation exactly proportional to the central temperature would require that the characteristic CDM density was independent on halo mass (i.e. constant concentration).

As in the baryon fraction, we argue that the characteristic density and radius can be more physically meaningful than the mass and radius at a given overdensity. Depending on concentration, the entropy at some fraction of $`R_{200}`$ can vary as much as a factor of 3. This can have important consequences on the $`ST`$ relation, since we have shown that merging systems are systematically much less concentrated than relaxed haloes of the same mass. A lower concentration implies a lower central gas density, but also a lower temperature and a smaller value of $`R_{200}/r_\mathrm{i}`$. Therefore, the net effect depends on the details of the mass-concentration relation.

The $`ST`$ relation for our sample of numerical clusters is shown in Figure 18. The entropy, $`S=kT[\rho _\mathrm{g}/(\mu m_\mathrm{p})]^{2/3}`$, is evaluated at 0.1 $`R_{200}`$ and plotted as a function of the X-ray temperature of the object (see Table 1). Systems that have been catalogued as mergers according to our substructure criterium display higher entropies than relaxed clusters. Since mergers are more common on group scales, the $`ST`$ relation becomes considerably flattened when these systems are taken into account.

## 6 Conclusions

In this paper, we have considered four self-consistent analytical models of galaxy clusters, based on the hypotheses that the hot ICM gas is in hydrostatic equilibrium with the dark matter halo and that it follows a polytropic equation of state. Two of our models assume NFW and MQGSL formulae to describe the CDM density profile, whereas the other two assume a $`\beta `$-model for the gas distribution. One is an isothermal version with $`\beta =2/3`$ (BM) and the other is a polytropic model with $`\gamma =1.18`$ and $`\beta =1`$ (PBM).

We performed a set of high-resolution adiabatic gasdynamical simulations to assess the accuracy of the basic approximations (i.e. hydrostatic equilibrium and a polytropic relation). Then, we compared the radial distributions of gas and dark matter expected in each analytical model with the results for a sample of 15 numerical groups and clusters between 1 and 3 keV.

Our main conclusions can be summarised as follows:

1. We find additional evidence of non-physical entropy losses in the standard implementation of the SPH algorithm. This effect can be of critical importance in the computation of the temperature and entropy profiles. The new formulation proposed by Springel & Hernquist (2002) appears to be a promising alternative to overcome this problem.
2. Thermally-supported hydrostatic equilibrium is a valid approximation for objects classified as ’relaxed’ on dynamical grounds, up to $`0.8R_{200}`$. We find evidence of non-thermal support in merging systems.
3. A polytropic equation of state with $`\gamma 1.18`$ provides a fairly good description of the ICM gas. Substructure induces deviations from this relation.
4. The density profile of our clusters can be well fitted by either NFW or MQGSL formulae. $`\beta `$-models predict constant dark matter density near the centre, in conflict with the numerical results. Major mergers feature a power-law density profile with slope $`2\alpha 1`$.
5. Gas temperature declines by a factor of $`23`$ from the centre to the virial radius. Our projected X-ray temperature profiles can be fitted by the ’universal’ form proposed by Loken et al. (2002), derived from Eulerian gasdynamical simulations, although we obtain different values for the free parameters. Purely adiabatic physics cannot account for a large isothermal core.
6. $`\beta `$-models fail to reproduce the gas density profile. While a fit to the inner part yields $`\beta 2/3`$, larger values are required at large radii. NFW and MQGSL models are able to predict the average distribution, but variations in the polytropic index lead to a very large scatter. The asymptotic slope of the gas density profile is too steep compared to the simulation data.
7. Models based on a ’universal’ CDM profile are able to estimate the ICM properties within $`3040`$ per cent errors, when applied to the simulated X-ray surface brightness. Assuming a $`\beta `$-model yields similar estimates for $`r0.1R_{200}`$, but the shape of the inferred profiles at smaller radii can be severely misleading.
8. Contrary to conventional SPH estimates, the entropy distribution is not a pure power law. All our analytical cluster models predict an increasing slope with radius, consistent with our numerical data. The entropy profile is entirely determined by the characteristic density and radius of the CDM halo, and the normalisation depends on the mass-concentration relation. In particular, it does not necessarily scale linearly with gas temperature.

In general terms, we claim that the radial structure of galaxy groups and clusters can be understood in terms of hydrostatic equilibrium and a polytropic equation of state for the gas, at least when only gravitational processes and adiabatic gasdynamics are considered.

Although $`\beta `$-models can sometimes provide a fair estimate of several cluster properties, they fail to provide an overall description that matches all the numerical data. Models based on a ’universal’ CDM density profile have a similar or better accuracy, and they are far more reliable in a qualitative sense (i.e. concerning the shape of the profiles). We argue that these models are preferable in order to analyse and interpret observational data.

## Acknowledgements

We thank Ben Moore for a useful discussion on the CDM density profiles. We also thank Volker Springel for providing the entropy conserving version of Gadget, and Andrey Kravtsov for providing the ART data in Figure 1. This work has been partially supported by the MCyT (Spain) under project number AYA-0973, by the Acciones Integradas Hispano-Alemanas HA2000-0026 and by Deutscher Akademischer Austauschdienst DAAD (Germany). We thank the Ciemat and CEPBA (Spain) and LRZ (Germany) for allowing us to use their supercomputers to perform the simulations reported in this paper.

## Appendix A Cluster Atlas

In the following image, we display our 15 simulated clusters, projected along the three major axes. Colours indicate the dark matter surface density along the line of sight. Contour plots represent the X-ray emission in conventional units.

Quite interestingly, it can be appreciated that many images of merging systems (both minor and major) appear to be spherically symmetric due to projection effects. Usually, a lot more substructure is present in the CDM density (which roughly corresponds to the optical emission) than in the hot ICM gas (X-ray).

Indeed, we find that the gaseous component of infalling galaxies is very efficiently stripped, while dark matter haloes are able to survive during many orbits around the cluster centre 333 A visual impression of this phenomenon can be obtained from a series of animations that can be accessed at http://pollux.ft.uam.es/gustavo/VIDEOS/LCDM80/clustervideo.html . We expect these haloes to retain some of their original gas when cooling is taken into account, probably undergoing an intense star formation burst. However, all gas in a warm-hot reservoir is lost during the first orbit, quenching any subsequent star formation activity. Numerical experiments that include self-consistent cooling, star formation and feedback are needed in order to make a quantitative assessment of this scenario.


format: text

text:
# Magnetophoresis in the Rubinstein – Duke model

## I Introduction

The magnetophoresis problems is a member of the class of reptation problems in which a long polymer is driven through a gel. The reptative motion of the polymer can be succesfully modeled by a lattice version using sections of the polymer (the reptons) as the mobile units, which hop stochastically from cell to cell. The driving field is incorporated as a bias in the hopping rates, favoring the motion in the field direction. A simple and adequate model is the Rubistein - Duke model RD , which represents the polymer as a chain of $`N`$ independently moving reptons, with the restriction that the integrity of the chain is preserved. The reptons trace out a connected string of cells in space, each cell containing at least one repton. The cells, which can be multiply occupied, carry the extra reptons as units of stored length. In order to preserve the integrity of the chain, only those reptons which are located in cells with stored length can hop.

Usually one considers the case of polyelectrolytes in which the reptons are uniformly charged. Thus the driving field pulls equally on all reptons and the bias is the same for the hopping rates all along the chain. A practical situation where this occurs is in DNA electrophoresis viovy . The DNA molecules, being acid, get charged in solution and when they are placed in a gel subject to an external electric field they perform a biased reptative motion along the field direction. Electrophoresis is a technique of great importance in molecular biology and sequence analysis, as it allows to separate DNA strands according to their length viovy .

In nature the charge distribution is of course not always uniform. The extreme alternative is the case where all reptons are neutral except one end repton, which is charged. A possible realization of such a situation is a magnetic bead, attached at one end of the polymer, which is driven by a magnetic field. Therefore this case can be referred to as the magnetophoresis (MP) problem bark96 . It is the subject of this paper. As we will frequently compare our findings with the more common case of the uniformly charged polymer, we will, for briefness, refer to the latter as the electrophoresis (EP) problem and to the present one as the MP problem, although the distinction between the two is not electromagnetic, but only in the forces exerted on the reptons (see Fig. 1).

The MP problem, within the framework of the Rubinstein–Duke model, has so far been studied by means of Monte Carlo simulations bark96 , and the calculations were mostly restricted to the drift velocity as function of the applied field and chain length. In this paper we analyze the MP problem by means of density-matrix renormalization-group (DMRG) techniques, which allow us to perform a detailed analysis of both global quantities as drift and curvilinear velocities and diffusion coefficient, but also on local average shapes of the polymer.

The dynamics of the reptating chain is governed by the Master Equation which we put in the form

$$\frac{P(𝐲,t)}{t}=\underset{𝐲^{}}{}H(𝐲,𝐲^{})P(𝐲^{},t)$$
(1)
Here $`𝐲`$ stands for the set of links $`y_1,\mathrm{},y_{N1}`$, where $`y_j`$ measures the distance between the repton $`j`$ and $`j+1`$, along the direction of the applied field. In our lattice representation the $`y_j`$ can take the values $`\pm 1`$ and $`0`$. The value $`y_j=0`$ corresponds to the case that the reptons $`j`$ and $`j+1`$ occupy the same cell. Thus each zero is a unit of stored length. The non–zero values represent the cases where $`j+1`$ occupies a cell “higher” (1) or “lower” (-1) than $`j`$. Higher and lower refer to a position in the direction of the field. $`𝐲`$ represents a complete configuration of the chain (see Fig. 1). $`P(𝐲,t)`$ is the probability distribution of the configuration at time $`t`$ and the matrix $`H(𝐲,𝐲^{})`$ contains the gain and loss transitions from $`𝐲^{}`$ to $`𝐲`$. The bias in the hopping rate is contained in the matrix elements of $`H`$. Generally we have for the bias factor

$$B_j=\mathrm{exp}(aq_jE/k_BT)$$
(2)
with $`q_j`$ the charge of repton $`j`$, $`E`$ the driving field, $`a`$ the distance between adjacent cells (measured along the field direction) and $`k_BT`$ the standard combination of Boltzmann’s constant and the absolute temperature. In the MP problem all $`q_j=0`$ except for $`j=N`$. We put

$$B_NB=\mathrm{exp}(\epsilon /2)$$
(3)
and use $`\epsilon `$ as the parameter for the driving field. The other parameter of the model is the number of reptons N.

We have chosen this “hamiltonian form” of the Master Equation in order to stress the formal correspondence with a quantum mechanical model governed by a hamiltonian matrix. The DMRG method exploits this analogy and indeed its success in one-dimensional quantum problems carries over to reptation problems paper1 ; paper2 ; paes02 . We are interested in the stationary state of the probability distribution. In the quantum language this corresponds to finding the right eigenvector of $`H`$ belonging to the eigenvalue zero. The adaptation of the DMRG method to the MP problem is straightforward and the data presented in this paper are obtained by the DMRG method. An important difference with the Master Equation is that in quantum mechanical problems the hamiltonian is hermitian, whereas in the reptation problem the matrix $`H`$ is non–hermitian, due to the influence of the driving field. This restricts the applicability of the DMRG–method to moderately long chains and/or small driving fields.

The physics of the EP and MP problems is qualitatively different. As illustration consider the weak field case ($`\epsilon `$ small), where the Nernst–Einstein relation $`v=FD`$ relates the drift velocity $`v`$ to the total applied force $`F`$ and the diffusion coefficient $`D`$. The force $`F`$, equals $`N\epsilon `$ in the EP case, since one pulls at each repton. It is well known that $`v`$ scales as $`v\epsilon /N`$, yielding for the diffusion the non–trivial result $`DN^2`$. In the MP problem $`D`$ will be essentially the same, as hopping is limited by the availability of stored length. In both cases the motion can be considered as diffusion of stored length. Since $`F=\epsilon `$ in the MP problem, we expect the drift velocity to scale as $`v\epsilon N^2`$, a feature which is born out by our calculations.

In the Sec. II we discuss some moment equations derived from the Master Equation which are more helpful then in the EP problem in analyzing the drift– and curvilinear velocity. They are expressed in terms of the probabilities $`n_j^k`$ that the link $`j`$ has the value $`y_j=k`$. The sum of the $`n_j^k`$ adds up to 1

$$n_j^0+n_j^++n_j^{}=1$$
(4)
So it suffices to consider the two quantities $`n_j^0`$ and $`m_j`$ defined by

$$n_j^0=\mathrm{\hspace{0.17em}1}y_j^2m_j=y_j=n_j^+n_j^{}$$
(5)
$`n_j^0`$ can be called the local density of stored length. $`m_j`$ is a measure for the local orientation and will be referred to as the profile of the chain.

In the Sections III and IV we present our data for the velocities and the profiles. The analysis is most transparent in the strong field limit where we can make an ansatz which almost perfectly represents the data. In section IV we discuss the behavior of the profile for weak and strong pulling fields.

## II Moments of the Master Equation

The DMRG method deals with the whole probability distribution $`P(𝐲)`$. In the MP problem it is fruitful to consider moments of the Master Equation. One set of moments is obtained by multiplying Eq. (1) with $`y_j`$ and then summing over all $`𝐲`$. This leads to $`N1`$ relations which can be seen as an expression of the fact that the drift velocity $`v`$ across all the $`N1`$ links of the chain must be the same on the average. An even more useful set of relations is obtained by multiplying the Master Equation by $`y_j^2`$ and summing over all $`𝐲`$. The resulting $`N1`$ relations are an expression of the fact that the curvilinear velocity $`J`$ is the same across all links. These relations obtain the form

$$J=n_{j1}^0n_j^01<j<N$$
(6)
which can be viewed as the familiar law that the current $`J`$ equals minus the gradient of the density of the stored length. In addition to Eq. (6) one has two relations footnote  concerning the traffic in and out both ends of the chain. They read

$`J`$ $`=`$ $`13n_1^0,`$
$`J`$ $`=`$ $`(n_{N1}^0n_{N1}^{})B+(n_{N1}^0n_{N1}^+)B^1`$ (7)
The expression for the drift velocity involves correlations between neighboring links

$$v=(1y_{j1}^2)y_j(1y_j^2)y_{j1}$$
(8)
and it is therefore not as informative as Eq. (6). As one sees Eq. (7) involves only averages over the first (last) link. This holds also for the expressions for the drift velocity in terms of the averages of the first and last link.

$`v`$ $`=`$ $`m_1,`$
$`v`$ $`=`$ $`(n_{N1}^0+n_{N1}^{})B(n_{N1}^0+n_{N1}^+)B^1`$ (9)
These equations have been derived by Barkema and Schütz bark96  using balance arguments.

Equation (6) is a powerful relation since it allows to express the density of stored length in terms of the curvilinear velocity

$$n_j^0=n_1^0(j1)J$$
(10)
showing that the density profile is linear in the position of the cell. In particular (10) implies a relation between the densities of the first and last cell

$$n_{N1}^0=n_1^0(N2)J$$
(11)
We have used the linearity of the density $`n_j^0`$ as a check of the numerical calculations.

Counting the number of unknowns ($`v`$, $`J`$, $`n_1^0`$, $`n_{N1}^0`$, $`m_1`$, $`m_{N1}`$) and the number of equations (7), (9) and (11) we see that we have one more unknown than equations. This situation is similar to the EP problem. There the expression for the curvilinear velocity does not obtain the simple form (6), due to the bias on the internal reptons. So one misses relation (11). On the other hand $`J=0`$ for EP, due to the symmetry of the polymer on exchanging head and tail. Thus in both cases the moment equations are not sufficient to determine the velocities. Higher moments do not lead to additional information since again higher order correlations appear.

## III Global quantities

We discuss first the behavior of global quantities as the drift and curvilinear velocities and the diffusion constant.

### III.1 The Weak Field Limit

In the weak field limit the polymer assumes mostly a random configuration and all the densities $`n_j^k`$ are close to 1/3. The overal behavior of the drift velocity $`v`$ as function of $`\epsilon `$ and $`N`$ is given in Fig. 2. Note that the drift velocity becomes proportional to $`\epsilon `$ for small $`\epsilon `$, as expected on the basis of the Nernst – Einstein relation discussed above. For stronger fields, the velocity saturates to a finite values, as discussed in the next paragraph. A similar dependence of $`v`$ on $`\epsilon `$ and $`N`$ was observed in the Monte Carlo study of Ref. bark96 . The limiting behavior for small $`\epsilon `$ and large $`N`$ is thus:

$$v(N)\epsilon D(N)\frac{\epsilon }{N^2},$$
(12)
with $`D(N)`$ the zero field diffusion coefficient. The scaling behavior of $`D(N)`$ as function of the length $`N`$ has been studied quite intensively wido91 ; bark97 ; bark98 ; fris00 . Reptation theory predicts that $`D(N)1/N^2`$, while conflicting results appeared on experimental measurements, for which both $`1/N^2`$ and, more recently, $`1/N^{2.3}`$ lodge  have been reported.

A detailed study of the scaling of $`D(N)`$, within the Rubinstein–Duke model by means of DMRG method was recently performed paper1 ; paper2 , for various end-point stretching rates. In that case the diffusion coefficient was calculated from the limiting value of the drift velocity for $`\epsilon 0`$, with the field acting on all reptons (the EP problem). Here we repeat the same analysis only for a single case (using a stretching rate $`d=1`$ following the definition of $`d`$ of Refs. paper1 ; paper2 ). The advantage of calculating $`D(N)`$ with a small field acting only on an end repton is that the DMRG procedure is much more stable in this case and one can compute longer chains. This is due to the fact that in the MP problem non-hermiticity is restricted only to the repton where the field is applied. As mentioned in the introduction non-hermiticity hampers the efficiency of the DMRG method.

In order to calculate the diffusion coefficient from the Nernst-Einstein relation $`D=lim_{\epsilon 0}v/\epsilon `$ in practice, we used a small field ($`\epsilon =10^3`$) and checked explicitely that results do not change for smaller fields. The scaling behavior of the diffusion coefficient is expected to be:

$$D(N)N^2=A+\frac{A^{}}{\sqrt{N}}+\mathrm{}$$
(13)
with $`A`$ and $`A^{}`$ some constants. The form of the subleading correction to $`D(N)`$ has been debated for a while bark ; prah  and recent DMRG results suggest that it is of the type $`1/\sqrt{N}`$ paper1 , supporting Eq. (13). The coefficient was determined exactly vLK ; prah ; Al-L : $`A=1/3`$.

To analyze the scaling behavior of $`D(N)`$ it is most convenient to use the logarithmic derivative of the DMRG data:

$$\alpha (N)=\frac{\mathrm{ln}\left[D(N)\right]\mathrm{ln}\left[D(N+2)\right]}{\mathrm{ln}N\mathrm{ln}(N+2)},$$
(14)
which is shown for $`N=9`$, $`11`$, …, $`51`$ in Fig. 3. Plugging Eq. (13) in Eq. (14) one finds for the effective exponent $`\alpha (N)=2+A^{}/(2A\sqrt{N})`$, a behavior which is accurately reproduced by our numerical data of Fig. 3. The present results corroborate previous claims paper1  about the scaling form of $`D(N)`$.

### III.2 The Strong Field Limit

In the strong field limit the polymer assumes an oriented configuration, with the ’+’ links dominating at the pulled end. At the other end we still have a substantial amount of links ’0’, since the polymer can only move by the diffusion of stored length from the tail to the head. Eliminating $`n_1^0`$ from Eq. (10) with the use of Eq. (7) we get

$$n_{N1}^0=\frac{1}{3}(1KJ),$$
(15)
where $`K=3N5`$. In order that $`n_{N1}^0`$ stays finite for $`N\mathrm{}`$, the curvilinear velocity must vanish as

$$JK^1N\mathrm{}$$
(16)
As one sees from (15) this limiting value is not sufficient to determine the limiting value of density $`n_{N1}^0`$, which is sensitive to the corrections to (16). For the strong field limit it is useful to relate the drift velocity to the curvilinear velocity. With (7) and (8) we get

$$v=J+2(n_{N1}^{}Bn_{N1}^0B^1),$$
(17)
Now, if the polymer is fully stretched, $`v`$ and $`J`$ become the same. In Fig. 4 we have plotted the difference $`vJ`$ as calculated by DMRG for various fields and chain lengths. We note that it is small for all values of $`\epsilon `$ and $`N`$, in particular for strong fields and that this tendency is enforced for long polymers. That it also is small in the small field limit is a consequence of the fact that both quantities vanish in that limit. In order for the difference to vanish we must have

$$n_{N1}^0B^2n_{N1}^{}$$
(18)
Now we may use this relation as the 6th relation, which enables us to make all the desired quantities explicit functions of $`\epsilon `$ and $`N`$. We find for instance

$`J(B)=v(B)={\displaystyle \frac{B^42B^2+1}{K(B^4+B^2+1)+3B^3}},`$ (19)
$`n_{N1}^0(B)={\displaystyle \frac{1+B/K}{B^2+1+3B/K+1/B^2}}.`$ (20)
This explicit field dependence is compared to the data in Fig. 5 and Fig. 6. The agreement is excellent in both cases. Note also that Eq. (19) is consistent with Eq. (16) and that it provides the proportionality coefficient.

However, the crossover formulae (19)-(20) do not describe the subtile dependencies in the limit of small fields $`B=\mathrm{exp}(\epsilon /2)1`$. In this limit the drift velocity vanishes as $`v\epsilon `$, while one observes from (19), that the curvilinear velocity vanishes as $`J\epsilon ^2`$. For this reason, in the limit $`\epsilon 0`$, the crossover formula (19) predicts $`v\epsilon ^2/(3N5)`$, in disagreement with the correct scaling behavior of Eq. (12). The strong field limit does not suffer from this problem. We note, for instance, that in the limit $`N\mathrm{}`$ the saturation value of the velocity for $`B\mathrm{}`$ is in agreement with the exact expression given in Ref. bark96 : $`v=1/(3N5)`$.

## IV Profiles

Next we discuss some profiles, i.e. the local orientation $`m_iy_i`$ as function of the segment position along the chain. We consider $`N1`$ segments, thus $`N`$-reptons with the charged one at head position $`N`$. Fig. 7 shows a plot of $`m_i/ϵ`$ as function of the scaled variable $`(i1)/(N2)`$ for chains of various lengths and at fixed field $`\epsilon =0.001`$. This profile corresponds to the linear regime where the drift velocity scales as $`v\epsilon `$. The notable feature is a symmetry between head and tail with respect to the center of the chain, although the magnetophoresis problem is clearly asymmetric. This symmetry can be shown next  to be strict in the weak field limit. It disappears for stronger values of the field as Fig. 8 shows, where profiles are plotted for $`ϵ=1`$ and various lengths N.

In order to analyze the data further we also plot the individual probabilities $`n_i^k`$ for having a $`+,0`$ or $``$ at the site $`i`$ of the chain. For small values of $`ϵ`$ (not shown here) the curves are all near $`1/3`$, with a slight excess of $`+`$ links at the head and a depletion of $``$ links. The densities of $`+`$ and $``$ links are monotonically increasing and decreasing functions of the position $`i`$. For intermediate fields $`ϵ=1`$, the densities are more interesting and in Fig. 9 we plot the values of $`n_i^0`$, $`n_i^+`$ and $`n_i^{}`$ for $`N=51`$. The linear behavior for $`n_i^0`$ is consistent with Eq. (10). The curve for $`n_i^+`$ is monotonically increasing, but that for $`n_i^{}`$ is not monotonically decreasing.

The qualitative behavior of the orientation profile can be understood by considering the “origin” of the non-zero links ($`y_i=\pm 1`$) as has been introduced by Barkema and Newman bark97 . In the MP problem more links are created at the pulled head than at the tail. They stream gradually down to the tail. We can keep track for every link $`y_i=\pm 1`$, whether it is formed at the head or at the tail. After sufficient time the chain is divided into two zones: a head zone and a tail zone. They are separated by a small intermediate region with zeroes (we do not follow the origin of the zeroes). The zones remain separated because the $`y_i=\pm 1`$, created at the head cannot cross the $`y_i=\pm 1`$ created at the tail. The division between the two zones fluctuates in time and occasionally the tail zone disappears, while very rarely (particularly at large fields) the head zone vanishes. The larger the force on the head, the larger the asymmetry between the head and tail zones. We supplement these speculations by making an assumption on the ratios

$$r_i(j)=p_i^+(j)/p_i^{}(j)$$
(21)
where $`p_i^\pm (j)`$ is the probability of finding a $`\pm `$ at $`i`$ when the division is at $`j`$. We put

$`r_i(j)=r_h`$ $`\mathrm{for}j<i`$ (22)
$`r_i(j)=1`$ $`\mathrm{for}ji`$ (23)
Since in the tail zone there is no distinction between $`+`$ and $``$ we have set the ratio equal to 1. The idea, underlying this assumption, is that the $`+`$ and $``$ links are interlocked. So while moving in their zone their ratio can not change.

At position $`i`$ the average number of nonzero links ($`+`$ and $``$) is equal to $`1n_i^0`$. We introduce $`f_i`$ as the fraction of such nonzero links which are in the head zone. One can express the densities $`n_i^\pm `$ in terms of $`f_i`$ as

$$\{\begin{array}{ccc}\hfill n_i^+& =& \left(f_i\frac{r_h}{r_h+1}+(1f_i)\frac{1}{2}\right)[1n_i^0]\hfill \\ \hfill n_i^{}& =& \left(f_i\frac{1}{r_h+1}+(1f_i)\frac{1}{2}\right)[1n_i^0]\hfill \end{array}$$
(24)
In both equations the terms proportional to $`f_i`$ are the contributions that $`i`$ is in the head zone while the terms proportional to $`1f_i`$ are the contributions from the case in which $`i`$ is in the tail zone. We see that in the magnetophoresis problem the situation simplifies, since we do not have to worry about the tail zone. It drops out when we consider the profile

$$m_i=n_i^+n_i^{}=\frac{r_h1}{r_h+1}f_i[1n_i^0]$$
(25)
Thus the profile $`m_i`$ is, apart from the known factor $`1n_i^0`$, directly related to the fraction $`f_i`$. The latter has a simpler interpretation. It starts out at $`i=1`$ with a value nearly zero, since the head zone will only seldomly extend over the whole chain. It ends at $`i=N1`$ at a value very close to 1, since the tail zone will hardly ever extend over the whole zone. We can use this fact to tie the ratio $`r_h`$ to the end point values, discussed earlier, by considering (25) for $`i=N1`$

$$y_{N1}=\frac{r_h1}{r_h+1}[1n_{N1}^0]$$
(26)
and solving for $`r_h`$. It leads to

$$y_i=y_{N1}f_j\frac{1n_i^0}{1n_{N1}^0}$$
(27)
This form contains only $`f_i`$ as unknown. We can draw some conclusions from Eqs. (25) and (27) for weak fields, as well as for long chains at stronger fields.

### IV.1 Weak Fields

For $`ϵ0`$ we may put

$$r_h=1+a_hϵ$$
(28)
The function $`1n_i^0`$ will approach the limit $`2/3`$, so (25) becomes

$$y_i=ϵ\frac{a_h}{3}f_i$$
(29)
For the zero field limit of the profile we can take the zero field limit of $`f_i`$. It has the property that head and tail become equivalent or

$$f_i=1f_{N+1i}\mathrm{or}f_i+f_{N+1i}=1$$
(30)
The zero field limit of $`f_i`$ has been determined in bark97  by Monte Carlo simulations. We note that (30) is consistent with the mentioned next  symmetry in $`y_i`$. One should have $`a_h=2`$ in order that the profile becomes $`2ϵ/3`$ at the head, as is observed (see Fig. 7). This is perfectly in agreement with the value $`r_h=B^41+2\epsilon `$ for small $`\epsilon `$. Combining (29) and the first Eq. (9) we find that $`f_1`$ is a measure for the drift velocity $`v`$. According to (12) $`f_1`$ should vanish as $`1/3N^2`$. This result has been derived in bark97 .

Another feature of Fig. 7 seems to be the collapse of the data on a single curve. Further data on longer chains show that the flattening–off at the ends of the chain shrinks with the size of the system and that the slope in the middle slowly decreases. This is another manifestation of the slow approach towards the asymptotic behavior next  for large $`N`$.

Note that if the division between the head and tail region were located with equal probability on all sites of the chain then one would have simply $`f_i=i/N`$, which from Eq. (29) implies a linear profile. The profile of Fig. 7 is linear only at the center of the chain, while it strongly deviates from linearity close to the edges. This implies that the probability of finding the division between the head and tail regions is flat in the center of the chain and drops off at the chain edges.

### IV.2 Long Chains and stronger Fields

In this case the head zone will be dominant beyond a certain point (i.e. $`f_i=1`$ for $`i>i_0`$ ) in the chain, thus the division between the head and tail zones is expected to become localized close to the end of the chain which is not pulled. The curves in Fig. 9 convincingly show this behavior. It is interesting to note that when $`f_i1`$, Eq. (24) becomes:

$$\{\begin{array}{ccc}\hfill n_i^+& =& [1n_i^0]\frac{r_h}{r_h+1}\hfill \\ \hfill n_i^{}& =& [1n_i^0]\frac{1}{r_h+1}\hfill \end{array}$$
(31)
It immediately implies that both $`n_i^+`$ and $`n_i^{}`$ are linearly increasing functions of $`i`$ in the head zone, being $`n_i^0`$ a linearly decreasing function of $`i`$ (see Eq. (10)). This explains the monotonic increase of $`n_i^{}`$ close to the pulled end shown in the inset of Fig. 9. Note that using Eq. (31) one can estimate $`r_h`$ from the ratio of the slopes of $`n_i^+`$ and $`n_i^{}`$ in the head region. We find a ratio $`r_hB^4`$ in agreement with our crossover formulae.

## V Discussion

We have presented a series of numerical and analytical results for the MP problem in the Rubinstein - Duke model, where a single reptating polymer is pulled by a constant driving field applied to one polymer end. We have shown that the numerical data for the drift- and curvilinear velocities can be quite well reproduced by simple interpolating formulas following from the assumption that both velocities are equal in the limit of long chains. Indeed the measured differences are small which shows that the polymer is fairly stretched by the pulling force.

We studied also local quantities, as the profiles, which provide information on the shape of the reptating chain. These are quite well understood using a representation in which the polymer is divided into a head and tail region, with different ratios of $`+`$ and $``$ links. At small fields the division between the two regions meanders trough the whole chain, and the probability of finding it close to the edges drops off. At strong fields the division gets localized close to the free end of the chain. Moreover some profiles show an unexpected non-monotonic behavior which has a simple interpretation in the interface picture. The precise shape of the profiles at weak fields close to the polymer edges, both at finite $`N`$ and in the asymptotic limit $`N\mathrm{}`$, will be discussed in details elsewhere next .

Acknowledgements: This work has been supported by the Polish Science Committee (KBN) under grant in years 2003-2005. A.D. acknowledges a grant from the Royal Netherlands Academy of Science (KNAW) enabling him to stay at the Leiden University where part of this investigation were carried out.