<figure><img src="image_1.png" /></figure>

# Moderate steepening of galaxy cluster dark matter profiles by baryonic pinching

Jesper Sommer-Larsen, Marceau Limousin

► To cite this version:

Jesper Sommer-Larsen, Marceau Limousin. Moderate steepening of galaxy cluster dark matter profiles by baryonic pinching. Monthly Notices of the Royal Astronomical Society, 2010, 408 (4), pp.1998–2007. 10.1111/j.1365-2966.2010.17260.x. hal-01438899

HAL Id: hal-01438899

https://hal.science/hal-01438899

Submitted on 29 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Moderate steepening of galaxy cluster dark matter profiles by baryonic pinching

Jesper Sommer-Larsen<fnref n="1" />,<fnref n="2" />,<fnref n="3" />* and Marceau Limousin<fnref n="1" />,<fnref n="4" />★

Accepted 2010 June 28. Received 2010 June 9; in original form 2009 June 28

## ABSTRACT

To assess the effect of baryonic 'pinching' of galaxy cluster dark matter (DM) haloes, cosmological (Λ cold dark matter) TREESPH simulations of the formation and evolution of two galaxy clusters have been performed, with and without baryons included.

The simulations with baryons invoke star formation, chemical evolution with non-instantaneous recycling, metallicity-dependent radiative cooling, strong starburst-driven galactic superwinds and the effects of a meta-galactic ultraviolet field, including simplified radiative transfer. The two clusters have $T \simeq 3$ and 6 keV, respectively, and, at $z \sim 0$, both host a prominent, central cD galaxy.

Comparing the simulations without and with baryons, it is found for the latter that the inner DM density profiles, $r \lesssim 50–100$ kpc, steepen considerably: $\Delta\alpha \sim 0.5–0.6$, where $-\alpha$ is the logarithmic DM density gradient. This is primarily due to the central stellar cDs becoming very massive, as a consequence of the onset of late-time cooling flows and related star formation. Once these spurious cooling flows have been corrected for, and the cluster gravitational potentials dynamically adjusted, much smaller pinching effects are found: $\Delta\alpha \sim 0.1$. Including the effects of baryonic pinching, central slopes of $\alpha \simeq 1.0$ and 1.2 are found for the DM in the two clusters, interestingly close to recent observational findings for galaxy cluster Abell 1703 based on strong gravitational lensing.

For the simulations with baryons, the inner density profile of DM and cluster gas (intracluster medium) combined is found to be only very marginally steeper than that of the DM, $\Delta\alpha \lesssim 0.05$. However, the total matter inner density profiles are found to be $\Delta\alpha \sim 0.5$ steeper than the inner profiles in the DM-only simulations.

**Key words:** methods: numerical – galaxies: clusters: general – galaxies: evolution – galaxies: formation – cosmology: theory.

## 1 INTRODUCTION

Large N-body cosmological simulations have been carried out for a decade, with the goal of making statistical predictions on dark matter (DM) halo properties. Because of numerical issues, most of these large cosmological simulations contain DM particles only. They all reliably predict that the three-dimensional (3D) density profile $\rho_{DM}(r)$ should fall as $r^{-3}$ at large radii (but within the virial radius). Observations have confirmed these predictions for cluster-sized haloes (e.g. Kneib et al. 2003; Pointecouteau, Arnaud & Pratt 2005; Mandelbaum, Seljak & Hirata 2008; Okabe et al. 2010). This agreement is likely to be connected with the fact that at large radius,

the density profile of a galaxy cluster is DM dominated and the influence of baryons can be neglected. On smaller scales, parametrizing the 3D density profile of the DM using a cuspy profile $\rho_{DM} \propto r^{-\alpha}$; DM-only simulations predict a logarithmic slope $\alpha \sim 1–1.5$ for $r \to 0$. The exact value of the central slope and its universality is debated (Navarro, Frenk & White 1997; Moore et al. 1998; Ghigna et al. 2000; Ricotti 2003; Navarro et al. 2004; Graham et al. 2006; Gao et al. 2008). Theoretical efforts indicate that the inner logarithmic slope of DM haloes should be $\alpha \sim 0.8$ (Austin et al. 2005; Hansen & Stadel 2006). Although clearly of interest, these DM-only studies and their predictions do not help much to make comparison with observations. Indeed, observing the central part (i.e. the inner $\sim 500$ kpc) of a galaxy cluster at any wavelength reveals the presence of baryons (in the forms of stars and X-ray-emitting hot gas). Thus, any attempt to compare observations to simulations in

<footnote n="1">¹Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark</footnote>

<footnote n="2">²Excellence Cluster Universe, Technische Universität München, Boltzmannstrasse 2, 85748 Garching, Germany</footnote>

<footnote n="3">³Marie Kruses Skole, Stavnsholtvej 29-31, DK-3520 Farum, Denmark</footnote>

<footnote n="4">⁴Laboratoire d'Astrophysique de Marseille, UMR 6610, CNRS-Université de Provence, 38 rue Frédéric Joliot-Curie, 13388 Marseille Cedex 13, France</footnote>

<footnote>*E-mail: jslarsen@astro.ku.dk (JS-L); marceau.limousin@oamp.fr (ML)</footnote>
the centre of galaxy clusters has to be made with numerical simulations (or calculations) taking into account the baryonic component and its associated physics.

On the numerical side, efforts are currently being made to include effects of baryons in a realistic way in the simulations (see below). Problems remain, however: for example, due to the so-called ‘overcooling problem’, numerical simulations tend to predict central brightest cluster galaxies which are too blue and too massive relative to observations. Different effects do compete when it comes to the central slope of the density profile: the cooling of gas in the central regions of galaxies and clusters is expected to lead to a more concentrated DM distribution (the so-called adiabatic contraction, see Blumenthal et al. 1986; Gnedin et al. 2004, hereafter G04; Gustafsson, Fairbairn & Sommer-Larsen 2006). On the other hand, dynamical friction heating of massive galaxies against the diffuse cluster DM could, in principle, flattens the slope of the DM density profile (El-Zant, Shlosman & Hoffman 2001; Nipoti et al. 2003, 2004; Ma & Boylan-Kolchin 2004; Johansson, Naab & Ostriker 2009). Also, the properties of the inner part of simulated galaxy clusters (even in DM-only simulations) can depend significantly on initial conditions as demonstrated in Kazantzidis, Magorrian & Moore (2004). Consequently, no coherent picture has yet emerged from N-body simulations when it comes to the shape of the inner density profile of structures.

On the observational side, efforts have been put on probing the central slope $\alpha$ of the underlying DM distribution. These analyzes have led to wide-ranging results, whatever the method used: X-ray (Arabadjis, Bautz & Garmire 2002; Ettori et al. 2002; Lewis, Buote & Stocke 2003; Zappacosta et al. 2006); lensing (Tyson et al. 1998; Smith et al. 2001; Sand, Treu & Ellis 2002; Dahle, Hannestad & Sommer-Larsen 2003; Gavazzi et al. 2003; Sand et al. 2004; Gavazzi 2005; Bradač et al. 2007; Limousin et al. 2008; Sand et al. 2008; Richard et al. 2009; Oguri et al. 2009) or dynamics (Kelson et al. 2002; Biviano & Salucci 2006). This highlights the difficulty of such studies and the possible large scatter in the value of $\alpha$ from one cluster to another.

In summary, one needs to probe observationally and numerically the behaviour of the underlying DM distribution (i.e. after the baryonic component has been separated from the DM component) in the central parts of galaxy clusters. The main difficulties are the following: (i) observationally, to be able to disentangle the baryonic component and the underlying DM distribution; (ii) numerically, to implement the baryonic physics into the simulations and (iii) then to compare both approaches in a consistent way.

Sand et al. (2002, 2004, 2008) carried out lensing analyzes aiming to probe the central mass distribution in six galaxy clusters, using the measured velocity dispersion profile of the cD galaxy as an extra constrain. They found central density slopes smaller than 1.

Recently, there have been a growing interest on Abell 1703, a massive $z = 0.28$ (Allen et al. 1992) X-ray cluster with a luminosity $L_X = 8.7 \times 10^{44}$ erg s$^{-1}$ (Böhringer et al. 2000). It contains a large number of gravitational arcs, enabling a very detailed lensing analysis. Moreover, although it displays an intriguing filamentary structure along the north-south direction, it looks rather circular from the geometrical configuration of its multiply-imaged systems, in particular its giant arc, located at large angular separation (~35 arcsec ~147 kpc). It is a unimodal cluster, likely to be relaxed and characterized by a central dominant elliptical cD galaxy. This makes it much easier to interpret the results of the modelling compared to bimodal clusters such as Abell 1689 (Limousin et al. 2007b; Riemer-Sørensen et al. 2009), Abell 2218 (Elíasdóttir et al. 2007), Abell 68 (Richard et al. 2007) or MS 2053.7-0449 (Verdugo,

de Diego & Limousin 2007). In fact, regular relaxed clusters are rare at such redshifts (Smith et al. 2005). Finally, it displays a remarkable lensing configuration, forming a central ‘ring’ composed of four bright images that represent a rare example of lensing by a hyperbolic umbilic catastrophe (Orban de Xivry & Marshall 2009). This ring is located close to the central cD galaxy (4-9 arcsec, corresponding to 17-38 kpc), providing a robust constraint in the very central part of the cluster.

This motivated several lensing works on Abell 1703 aimed to probe the DM density slope of the central region (~20-200 kpc), modelling the DM distribution of Abell 1703 by a generalized Navarro, Frenk & White (NFW) profile, viz.

$$ \rho(r) = \frac{\rho_c \delta_c}{(r/r_s)^{\alpha_{\text{NFW}}} [1 + (r/r_s)^3]^{3-\alpha_{\text{NFW}}}} \quad (1) $$

where $r_s$ is the scale radius. These works (Limousin et al. 2008; Oguri et al. 2009; Richard et al. 2009; Saha & Read 2009) determined the inner slope $\alpha_{\text{NFW}}$ to be consistent with the universal profile ($\alpha_{\text{NFW}} = 1$). The more recent strong lensing analysis by Richard et al. (2009), based on the identification of 16 multiply-imaged systems, of which 10 are spectroscopically confirmed, found $\alpha_{\text{NFW}} = 0.92 \pm 0.04$ (1$\sigma$ confidence level). Note that the scale radius derived by Richard et al. (2009) is found in agreement with that found by Oguri et al. (2009) when combining strong and weak lensing. Since the scale radius needs weak lensing data to be properly constrained, this is important given the degeneracies arising between the scale radius and the inner slope.

Assuming that pure DM haloes are well described by the standard NFW profile with $\alpha_{\text{NFW}} = 1$, the indication of the above results on Abell 1703 is that the baryonic pinching of the DM halo, probably mainly caused by the central dominant cD galaxy, does not substantially affect the DM density slope in the central region of the cluster.

On the theoretical/numerical front, it has only recently been possible to carry out fully cosmological gas-dynamical/N-body simulations of the formation and evolution of galaxy clusters at a sufficient level of numerical resolution and physical sophistication that the cool-out of gas, star formation, chemical evolution and gas inflows and outflows related to individual cluster galaxies can be modelled to a reasonable degree of realism (e.g. Kay et al. 2002; Thomas et al. 2002; Valdarnini 2003; Tornatore et al. 2004; D’Onghia et al. 2005; Romeo, Portinari & Sommer-Larsen 2005; Sommer-Larsen, Romeo & Portinari 2005; Muanwong, Kay & Thomas 2006; Romeo et al. 2006; Saro et al. 2006; Tornatore et al. 2007; Romeo et al. 2008).

Romeo et al. (2005, 2006), D’Onghia et al. (2005) and Sommer-Larsen et al. (2005) presented fully cosmological simulations of galaxy groups and clusters. The TREESPH code used was building on the code used for simulating galaxy formation (e.g. Sommer-Larsen, Götz & Portinari 2003), improved to include modelling of non-instantaneous chemical evolution (Lia, Portinari & Carraro 2002a), metallicity-dependent, atomic radiative cooling, strong supernova and (optionally) active galactic nuclei (AGNs), driven galactic winds and thermal conduction. The two clusters simulated have $z=0$ virial masses<fnref n="1" /> $M_{\text{vir}} \sim 3 \times 10^{14}$ and $1.2 \times 10^{15}$ M$_{\odot}$,

<footnote n="1">¹ In this paper, the virial radius $r_{\text{vir}}$ is the radius corresponding to, at $z=0$, an overdensity 337 times the mean density of the Universe, as appropriate for a top-hat collapse in the adopted ΛCDM cosmology (e.g. Bryan & Norman 1998). The virial mass $M_{\text{vir}}$ of the cluster is the mass enclosed within the region of radius $r_{\text{vir}}$.</footnote>
one approximately the size of the Virgo cluster and the other of the Coma cluster. They were both selected to be fairly relaxed, and both display central prominent cD galaxies at $z=0$ as well as at $z=0.28$. In this paper, we reanalyze the simulations of the two clusters with emphasis on the DM density profiles. We also present resimulations of the clusters aimed at correcting for the effects of spurious, late-time cooling flows. Finally, we present results of new simulations of the *same* two clusters, identical to the previous ones except that they were carried out with DM only. We then compare the results of simulations with and without baryons, to quantify the effects of baryons (especially the cDs) on galaxy cluster DM profiles.

This paper is organized as follows: the code and the simulations are described in Section 2, the results obtained are presented in Section 3 and discussed in Section 4, and, finally, Section 5 constitutes the conclusion.

Throughout this paper, a Λ cold dark matter (ΛCDM) cosmology with $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$ and a Hubble constant $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$ is assumed.

## 2 THE CODE AND SIMULATIONS

The code used for the simulations is a significantly improved version of the TREESPH code we have used previously for galaxy formation simulations (Sommer-Larsen et al. 2003): full details on the code are given in Romeo et al. (2006); here we recall the main improvements over the previous version. (1) In lower resolution regions (which will always be present in cosmological CDM simulations), an improvement in the numerical accuracy of the integration of the basic equations is obtained by solving the entropy equation rather than the thermal energy equation – we have adopted the ‘conservative’ entropy equation solving scheme suggested by Springel & Hernquist (2002). (2) Cold high-density gas is turned into stars in a probabilistic way as described in Sommer-Larsen et al. (2003). In a star formation event, a smoothed particle hydrodynamic (SPH) particle is converted fully into a star particle. Non-instantaneous re-cycling of gas and heavy elements is described through probabilistic ‘decay’ of star particles back to SPH particles as discussed by Lia et al. (2002a). In a decay event, a star particle is converted fully into an SPH particle. (3) Non-instantaneous chemical evolution tracing 10 elements (H, He, C, N, O, Mg, Si, S, Ca and Fe) has been incorporated in the code following Lia et al. (2002a) and Lia, Portinari & Carraro (2002b); the algorithm includes supernovæ of Type II and Type Ia, and mass loss from stars of all masses. Metal diffusion in Lia et al. was included with a diffusion coefficient $\kappa$ derived from models of the expansion of individual supernova remnants. A much more important effect in the present simulations, however, is the redistribution of metals (and gas) by means of starburst-driven ‘galactic superwinds’ (see point 5). This is handled self-consistently by the code, so we set $\kappa = 0$ in the present simulations. (4) Atomic radiative cooling depending both on the metal abundance of the gas and on the meta-galactic ultraviolet (UV) field, modelled after Haardt & Madau (1996), is invoked. We also include a simplified treatment of radiative transfer, by switching off the UV field where the gas becomes optically thick to Lyman limit photons on scales of ~1 kpc. (5) Starburst-driven, galactic superwinds are incorporated in the simulations. This is required to expel metals from the galaxies and get the abundance of the intracluster medium (ICM) to the observed level of about one-third solar in iron. A burst of star formation is modelled in the same way as the ‘early bursts’ of Sommer-Larsen et al. (2003), i.e. by halting cooling in the surrounding gas particles, to mimic the initial heating and subsequent adiabatic expansion phase of the supershell powered by the star-

burst; this scheme ensures effective energy coupling and feedback between the bursting star particle and the surrounding gas. The strength of the superwinds is modelled through a free parameter $f_{\text{wind}}$ which determines how large a fraction of the newborn stars partake in such bursting, superwind-driving star formation. We find that in order to get an iron abundance in the ICM comparable to observations, $f_{\text{wind}} \gtrsim 0.5$ and at the same time a fairly top-heavy initial mass function (IMF) has to be used. (6) Thermal conduction was implemented in the code following Cleary & Monaghan (1999).

In this paper, we present results for two simulated clusters, ‘Virgo’ and ‘Coma’. Both systems were chosen to be fairly relaxed (no $\gtrsim 1:2$ merging at $z \lesssim 1$). Virial masses at $z=0$ are approximately $2.0 \times 10^{14}$ and $8.7 \times 10^{14} h^{-1} \text{ M}_\odot$ and X-ray emission weighted temperatures 3.0 and 6.0 keV, respectively. The clusters were selected from a cosmological, DM-only simulation of a flat ΛCDM model, with $\Omega_M = 0.3$, $\Omega_b = 0.036$, $h = 0.7$ and $\sigma_8 = 0.9$ and a box length of $150 h^{-1} \text{ Mpc}$. Mass and force resolution were increased in, and gas particles added to, Lagrangian regions enclosing the clusters. The Virgo cluster was then resimulated using 2.3 million baryonic+DM particles with $m_{\text{gas}} = m_* = 3.1 \times 10^7$ and $m_{\text{DM}} = 2.3 \times 10^8 h^{-1} \text{ M}_\odot$ for the high-resolution gas, star and DM particles. Gravitational (spline) softening lengths of $\epsilon_{\text{gas}} = \epsilon_* = 1.4$ and $\epsilon_{\text{DM}} = 2.7 h^{-1} \text{ kpc}$, respectively, were adopted. The more massive Coma cluster was resimulated using 1.0 million baryonic+DM particles with $m_{\text{gas}} = m_* = 2.5 \times 10^8$ and $m_{\text{DM}} = 1.8 \times 10^9 h^{-1} \text{ M}_\odot$, $\epsilon_{\text{gas}} = \epsilon_* = 2.8$ and $\epsilon_{\text{DM}} = 5.4 h^{-1} \text{ kpc}$. As a resolution test, the Virgo cluster was also simulated at this (eight times lower) mass and (two times lower) force resolution. This simulation will be denoted as ‘VirgoLR’.

For all simulations, gravitational softening lengths were fixed in comoving coordinates till $z=6$, and from then kept constant in physical coordinates at the $z=6$ values.

For the simulation presented in this paper $f_{\text{wind}} = 0.8$, and an Arimoto–Yoshii IMF (of slope $x=-1$, shallower than the Salpeter slope $x=-1.35$) with mass limits [0.1–100] $M_\odot$ was adopted. AGN-driven winds were not invoked. Moreover, the thermal conductivity was set to zero assuming that thermal conduction in the ICM is highly suppressed by magnetic fields (e.g. Ettori & Fabian 2000). We note that, in relation to all relevant quantities presented in this paper, no significant difference is found between simulations with zero thermal conductivity and simulations invoking thermal conduction at one-third of the Spitzer level (Romeo et al. 2005).

### 2.1 Correcting for the effects of late quasi-stationary cooling flows

After a period of major merging at $z \gtrsim 2$ strong, quasi-stationary cooling flows develop at the centres of the clusters despite the strong, supernova-driven energy feedback to the IGM/ICM through galactic superwinds and the use of a fairly top-heavy IMF. As a consequence, stars continue to form steadily at the centres of the cDs ($r \lesssim 10 \text{ kpc}$) at rates of $\sim 100$ and $\sim 500 M_\odot \text{ yr}^{-1}$ at $z=0$ for Virgo and Coma, respectively. As discussed in Sommer-Larsen et al. (2005), this results in the cDs becoming too blue over time, with central B–R colours of $\sim 1.0$ at $z=0$, rather than $1.4–1.5$ observed for real cDs (note though that some star formation is observed in many cDs at the base of the cooling flow, e.g. McNamara 2004). Given the significant gas cool-out rates, the increase in the stellar masses of the cDs since $z \sim 2$ is considerable: for the Virgo cD $M_{*12} = 0.42$, 1.0 and 1.4 in units of $10^{12} M_\odot$, at $z=2, 1, 0$, respectively (inside of 20 kpc). For Coma, the corresponding masses are 0.95, 1.9 and 4.2. In real cDs, the late-time cooling flows will be counteracted, most likely by the energy release related to the accretion of mass on
<table><caption>**Table 1.** Numerical and physical characteristics of the cluster simulations: mass of DM/gas/star particles and the respective gravitational (spline) softening lengths; total number of particles and initial redshift of each run; virial mass, radius and temperature of simulated clusters (last three columns, referring to $z=0$).</caption><thead><tr><th>Run</th><th>m<sub>DM</sub></th><th>m<sub>gas</sub><br>(10<sup>7</sup> M<sub>☉</sub> h<sup>-1</sup>)</th><th>m<sub>*</sub></th><th>ε<sub>DM</sub></th><th>ε<sub>gas</sub><br>(kpc h<sup>-1</sup>)</th><th>ε<sub>*</sub></th><th>N<sub>tot</sub></th><th>z<sub>4</sub></th><th>M<sub>vir</sub><br>(10<sup>14</sup> M<sub>☉</sub>)</th><th>R<sub>vir</sub><br>(Mpc)</th><th>kT<br>(keV)</th></tr></thead><tbody><tr><td>Coma</td><td>180</td><td>25</td><td>25</td><td>5.4</td><td>2.8</td><td>2.8</td><td>950 000</td><td>19</td><td>12.4</td><td>2.9</td><td>6</td></tr><tr><td>Virgo</td><td>23</td><td>3.1</td><td>3.1</td><td>2.7</td><td>1.4</td><td>1.4</td><td>2235 000</td><td>39</td><td>2.8</td><td>1.7</td><td>3</td></tr><tr><td>VirgoLR</td><td>180</td><td>25</td><td>25</td><td>5.4</td><td>2.8</td><td>2.8</td><td>260 000</td><td>19</td><td>2.8</td><td>1.7</td><td>3</td></tr><tr><td>ComaDM</td><td>26</td><td>–</td><td>–</td><td>2.8</td><td>–</td><td>–</td><td>4000 000</td><td>39</td><td>12.4</td><td>2.9</td><td>–</td></tr><tr><td>VirgoDM</td><td>26</td><td>–</td><td>–</td><td>2.8</td><td>–</td><td>–</td><td>1200 000</td><td>39</td><td>2.8</td><td>1.7</td><td>–</td></tr></tbody></table>

to supermassive black holes (e.g. Bower et al. 2006; Croton et al. 2006; Dubois et al. 2010; Teyssier et al. 2010). This will lead to a strong suppression of the late-time star formation rate, such that the simulated cDs become too massive. As will be shown in the next section, the gravitational force of this excess mass influences the central DM profile and must be corrected for.

Assuming that the excess mass has been added adiabatically, i.e. on a time-scale large compared to the stellar orbital time-scale in the cD, the ‘correct’ DM distribution can be recovered by slowly removing the excess stellar mass. This is accomplished as follows.

In order to enable a direct comparison to the observational results obtained for Abell 1703, we focus on DM profiles at $z=0.28$, corresponding to $t=10.2$ Gyr. At $t=8.2$ Gyr, stars at $r<20$ kpc from the centre of the cD are identified. For each of the two simulated clusters two additional simulations are carried out – in one the masses of all identified stars formed since $z=1$ are then gradually reduced to zero over a period of 1 Gyr, in the other the masses of all identified stars formed since $z=2$. At the same time radiative cooling as well as star formation is switched off. Subsequently the additional simulations are then continued for 2 Gyr more, still with radiative cooling and star formation switched off. We refer to these additional simulations as Rz1 and Rz2, respectively.

The purpose of these additional simulations is to, without resorting to completely new simulations with AGN feedback, etc., included, determine the cluster DM profiles in the (more realistic) case where effects of late-time cooling flows are counteracted since redshifts of either 1 or 2.

In removing the excess stellar mass, we are not putting the same amount of gas back into the simulation, i.e. mass is formally not conserved. We are justified in doing so, if this gas in a realistic simulations involving AGN feedback will reside well outside of the $r \sim 20$ kpc region. This will most likely be the case – in the model of Springel, Di Matteo & Hernquist (2005), for example, strong winds tend to drive the gas to considerable distances from the galaxy hosting the supermassive black hole.

## 2.2 DM-only simulations

To further complement the existing simulations and determine the role of the inclusion of baryonic physics, we carried out DM-only simulations of the two clusters using the same initial conditions as previously, but not splitting the particles into SPH and DM particles.

The Virgo and Coma cluster were resimulated using 1.2 and 4.0 million DM particles, respectively. These simulations will be denoted as ‘VirgoDM’ and ‘ComaDM’ in the following. The high-resolution DM particles had $m_{DM} = 2.6 \times 10^8 h^{-1} M_{\odot}$ and $\epsilon_{DM} = 2.8 h^{-1}$ kpc.

Numerical and physical parameters for the simulations and the $z=0$ clusters are summarized in Table 1.

# 3 RESULTS

## 3.1 Inner DM halo profiles

To enable a precise determination of the galaxy cluster DM density profiles at $z=0.28$, 21 time frames, with a time spacing of 0.1 Gyr (large enough to enable quasi-random sampling of the central halo profile) and centred on $t=10.2$ Gyr, were co-added. The DM profile power-law index

$$ \alpha(r) = - \frac{d \ln[\rho_{DM}(r)]}{d \ln(r)}, \qquad (2) $$

where $r$ is the cluster-centric radial distance, was then determined by averaging the co-added DM distribution over spherical shells.

In Fig. 1, the resulting $\alpha(r)$ [strictly speaking $-\alpha(r)$, but we shall neglect this distinction in the following) for the DM-only simulations is shown by dotted lines. Moreover, $\alpha(r)$ is shown by thick solid curves for the original simulations including baryons, i.e. simulations including baryons, but without any corrections for late cooling flows and related star formation. Finally, shown by thin solid lines is $\alpha(r)$ for generalized NFW profiles (equation 1) with $\alpha_{NFW} = \alpha(0) = 1.09$ and $0.92$, as determined for Abell 1703 by Limousin et al. (2008) and Richard et al. (2009), respectively, and concentration parameter $c=6$. It is seen from the figure that (a) all the simulations match the NFW profiles in the intermediate and outer parts of the haloes, $r/r_{200} \gtrsim 0.07$, fairly well, (b) the DM-only simulation of Virgo is matched quite well by an $\alpha(0) \simeq 1.1$ NFW profile and that of Coma by an $\alpha(0) \simeq 0.9$ profile and (c) the simulations with baryons result in too large values of $\alpha(r)$ (~1.5–1.6) in the inner parts, $r/r_{200} \lesssim 0.04$ (the results are affected by effects of gravitational softening at $r/r_{200} \lesssim 0.007–0.008$ – see Section 3.2.1 for more detail). In order to assess whether this is due to baryonic ‘over’-pinching, i.e. due to the simulated cDs being unrealistically massive, we compare the simulated and observed Virgo and Coma clusters:

At $z=0$, the cDs in the simulated clusters have masses (inside of $r=50$ kpc – see below) of $7.4 \times 10^{12}$ and $2.2 \times 10^{12} M_{\odot}$ for Coma and Virgo, respectively (at $z=0.28$ the corresponding numbers are $5.1 \times 10^{12}$ and $1.8 \times 10^{12} M_{\odot}$). The real Coma and Virgo clusters both contain two dominant galaxies: NGC 4874 and NGC 4889 in Coma, and M86 and M87 in Virgo. The two dominant galaxies in each cluster will eventually merge into one cD in each. The B-band luminosities of NGC 4874 and NGC 4889 are $1.5 \times 10^{11}$ and $1.8 \times 10^{11} L_{\odot,b}$, respectively (Michard & Andreon 2008). Assuming a B-band mass-to-light ratio (M/L) of about 8 (as found for the cDs in the present simulations based on the Arimoto–Yoshii IMF; the more
<figure><img src="image_3.png" /></figure>

**Figure 1.** Logarithmic density slopes of the spherically averaged cluster DM distributions at $z = 0.28$. Solid red curve shows result for the ‘Virgo’ cluster including the cD uncorrected for effects of late-time cooling flows. Dotted red curve shows the result of the Virgo DM-only simulation. The two red dashed lines show the results of the Rz1 (lower) and Rz2 (upper) resimulations, where effects of spurious late-time cooling flows have been corrected for. The blue curves show the same for the ‘Coma’ cluster. The two black solid curves show generalized NFW models (equation 1) with $\alpha_{NFW} = 0.92$ (upper) and 1.09 (lower), and concentration parameter $c = r_{200}/r_s = 6$. Statistical uncertainties are indicated for the Rz2 Coma resimulation. The uncertainties are similar for the other simulation curves shown in this figure, as well as the ones shown in Figs 5 and 6.

standard Salpeter IMF results in similar M/L), the combined stellar mass of the two galaxies (at $z \simeq 0$) is $2.6 \times 10^{12} M_{\odot}$. Moreover, the simulated Coma cluster is a $T_X \sim 6$ keV cluster, whereas the real Coma cluster is a $T_X \sim 8$ keV cluster. Assuming a standard scaling of mass $\propto T_X^{3/2}$, the cD in a 6 keV cluster should have a stellar mass of $1.7 \times 10^{12} M_{\odot}$. The B-band luminosities of M86 and M87 are $4.2 \times 10^{10}$ and $6.6 \times 10^{10} L_{\odot,b}$, respectively (Bingelli, Sandage & Tammann 1985). The stellar masses of the two galaxies thus add up to about $8.6 \times 10^{11} M_{\odot}$.

The simulated cDs are hence about three to four times more massive than observed, indicating that the simulated DM haloes are too contracted compared to real ones, which at least qualitatively could explain some of the differences seen in Fig. 1.

To understand the differences seen in Fig. 1 more quantitatively, we now turn to the resimulations aimed at correcting for the effects of the late cooling flows. In Fig. 1, we show by thick dashed lines $\alpha(r)$ for the simulations with baryons, where the mass of all stars in the central cD, formed since either $z_f = 1$ or 2, has been adiabatically reduced to zero – we shall denote these two sets of resimulations by Rz1 and Rz2 in the following. For the inner parts of Coma, a large effect on $\alpha(r)$ results already from removing central stars formed since $z_f = 1$. This is related to the fact that the mass of the cD more than doubles since $z = 1$ in the original simulation including baryons. For Virgo a strong effect is seen only between the original simulation and $z_f = 2$ correction. This is also reasonable, since the mass of the cD more than doubles going from $z = 2$ to 1, but then only increases by additional 40 per cent going to redshift 0. In order to assess which of the Rz1 and Rz2 resimulations, if any, result in cDs which are realistic, we now derive the properties of the resulting cDs.

## 3.2 Properties of the simulated cD galaxies

### 3.2.1 cD stellar masses

In Fig. 2, the cumulative stellar mass for the Coma and Virgo cDs from the $z_f = 1$ and 2 adiabatically corrected simulations is shown, respectively. It is seen that at $r = 50$ kpc, corresponding to approximately the $25B \text{ mag arcsec}^{-2}$ (optical) isophote (see below), for $z_f = 1$ the cumulative stellar masses are $3.1$ and $1.3 \times 10^{12} M_{\odot}$ for Coma and Virgo, respectively. For $z_f = 2$, the corresponding numbers are $2.1$ and $0.74 \times 10^{12} M_{\odot}$. Comparison to the observational estimates above, indicates that the $z_f = 2$ cDs provide a better match to reality than the $z_f = 1$ ones.

An alternative comparison can be made to the cD of Abell 1703, for which Limousin et al. (2008) estimate a stellar mass of $1.3 \pm 0.3 \times 10^{12} M_{\odot}$ inside of $r = 30$ kpc. The temperature of the Abell 1703 ICM is uncertain, but (also uncertain) estimates of the virial mass of Abell 1703 indicate that the temperature is likely at least as large as for the simulated Coma cluster (Riemer-Sørensen, private communication, Oguri et al. 2009). For the Coma Rz1 and Rz2 cDs, at $z = 0.28$, we find stellar masses inside of $r = 30$ kpc of $2.2 \times 10^{12}$ and $1.2 \times 10^{12} M_{\odot}$. Again, the stellar mass appears realistic for the Rz2 Coma cD, but too large for the Rz1 cD.

### 3.2.2 cD effective radii

Next, we determine the effective radii for the cDs at $z = 0$. To this end, the projected B-band surface brightness profiles are determined for the cDs (see Romeo et al. 2005 for details on calculating broad-band photometric properties of the star particles): first, a B-band surface brightness map is constructed by projecting the star-particle distribution within a cube of side-length $L = 200$ kpc, centred on the cD, along the three cardinal directions and then superposing the three maps (larger values of $L$ result in similar results). Secondly, a radial surface profile is constructed by azimuthally averaging over the combined image. Thirdly, the radial profile obtained is divided by a factor of 3. At $R \simeq 50$ kpc, the B-band surface brightness for both cDs drops to $25 \text{ mag arcsec}^{-2}$ (Sommer-Larsen et al. 2005), which is taken to be the radial extent of the cD itself. $R_{\text{eff}}$ is then determined as the projected radius enclosing half of the cD B-band

<figure><img src="image_4.png" /></figure>

**Figure 2.** For resimulations Rz1 (solid curves) and Rz2 (dashed curves), at $z = 0$, the cD cumulative stellar masses are shown (Virgo: red curves; Coma: blue curves).
luminosity. For the Coma and Virgo cDs in the Rz1 simulations, $R_{\text{eff}} = 8$ and 7 kpc is found. The corresponding values for the Rz2 simulations are $R_{\text{eff}} = 16$ and 11 kpc. For a large (observational) sample of Sloan Digital Sky Survey (SDSS) galaxies (at $z \sim 0$), Shen et al. (2003) find for early-type galaxies of stellar masses 3.1 and $1.3 \times 10^{12} M_{\odot}$ median effective radii of $\sim 25$ and 16 kpc, and for stellar masses of 2.1 and $0.74 \times 10^{12} M_{\odot}$, $R_{\text{eff}} \sim 20$ and 13 kpc. It follows that the Rz1 cDs have too small effective radii for their mass, whereas the Rz2 cDs have reasonable effective radii compared to observations.

### 3.2.3 cD metal abundances and $\alpha$/Fe ratios

In Fig. 3, for cDs in the Rz1 and Rz2 resimulations, the projected average stellar oxygen abundance as a function of projected radius $R$ is shown ([O/H] for a particular bin of projected radius is calculated on the basis of the metallicity information for the individual star particles comprising the bin, as described in Sommer-Larsen et al. 2005 and Romeo et al. 2005). In these $\alpha$-element enhanced systems (see below), [O/H] can be taken as a measure of the total abundance [Z/H]. It is seen that the oxygen abundance is supersolar all the way to 50 kpc projected radius. For the Rz1 simulations, the integrated (global) cD abundances are [O/H] = 0.59 for both clusters. For the Rz2 simulations, the integrated abundances are [O/H] = 0.40 and 0.26 for Coma and Virgo, respectively. Following Thomas, Maraston & Bender (2005), observed early-type galaxies of stellar masses 3.1 and $1.3 \times 10^{12} M_{\odot}$ have median $\langle [Z/H] \rangle = 0.38$ and 0.33, respectively, and for stellar masses of 2.1 and $0.74 \times 10^{12} M_{\odot}$, $\langle [Z/H] \rangle = 0.36$ and 0.30. So the Rz1 cDs have too large abundances compared to what is observed, whereas the Rz2 cDs provide a good match to observations. Taking [O/Fe] as a proxy for $\alpha$/Fe, the results for the Rz1 and Rz2 simulations are similar, [O/Fe] $\sim 0.3-0.35$. These values are within the range observed for large early-type galaxies (Thomas et al. 2005).

### 3.2.4 cD $U-V$ colours

For $z \sim 0$ early-type galaxies, the rest-frame $U-V$ broad-band colour allows sampling of the 4000 Å break and therefore is par-

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stellar oxygen abundance as a function of projected radial distance from the centre of the cD for resimulations Rz1 and Rz2, at $z = 0$. The curves are labelled as in Fig. 2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** $U-V$ colour as a function of projected radial distance from the centre of the cD for resimulations Rz1 and Rz2, at $z = 0$. The curves are labelled as in Fig. 2.</figcaption></figure>

ticularly sensitive to age and metallicity variations of the stellar populations. In order to check whether the simulated cDs fall on the observational ‘red sequence’ for early-type galaxies, we hence calculate the $U-V$ colours for the cDs. The $U-V$ colours are calculated on the basis of the abundances and ages of the star particles comprising the galaxies, as discussed in Romeo et al. (2005). In Fig. 4, for cDs in the Rz1 and Rz2 resimulations, the $U-V$ colour as a function of projected radius $R$ is shown. In the inner parts of the galaxies, all the cDs display a negative colour gradient, qualitatively consistent with what is observed for large early-type galaxies. As the mean stellar age is approximately constant with projected radius $R$, the reason for the gradient in colour is the metallicity gradient seen in Fig. 3. This is also the generally accepted explanation for the observed colour gradients in early-type systems (e.g. Tamura et al. 2000). For the Rz1 cDs, integrated (global) $U-V$ colours of 1.31 and 1.27 are found for Coma and Virgo, respectively. For the Rz2 cDs, the corresponding numbers are $U-V = 1.29$ and 1.28. The two Rz1 cDs have $M_V = -24.6$ and $-23.7$ for Coma and Virgo, respectively. For the Rz2 cDs, the corresponding numbers are $-23.8$ and $-23.0$. Comparing to the observational $U-V$ versus $M_V$ sequence (e.g. Romeo et al. 2008), early-type galaxies of such absolute $V$ magnitudes have a median $U-V \sim 1.4$. We hence find that the Rz1 and Rz2 cDs have very similar $U-V$ colours and that these are slightly bluer than the observed average. This is, at least partly, due to the large, somewhat bluer stellar envelope surrounding the cDs (cf. Fig. 4).

### 3.2.5 cD central physical stellar densities

From the results presented above it is possible to estimate the stellar core density of the cDs, defined as

$$ \rho_{c,*} = \frac{M_{*}(r < R_{\text{eff}})}{4\pi/3R_{\text{eff}}^3}. \qquad (3) $$

For the Rz1 cDs, we find $\rho_{c,*} = 6$ and $4 \times 10^8 M_{\odot} \text{ kpc}^{-3}$ for Coma and Virgo, respectively. For the Rz2 cDs, the corresponding numbers are 5 and $4 \times 10^7 M_{\odot} \text{ kpc}^{-3}$. Extrapolating the observational estimates of $\rho_{c,*}$, given in van Dokkum et al. (2008) on the basis of SDSS data, to galaxy stellar masses of $\sim 10^{12} M_{\odot}$, the Rz2
cD core densities appear the most reasonable, although no strong conclusions can be made on the basis of the available data.

## 3.3 How much are galaxy cluster DM haloes pinched by the central cDs?

Based on the results given in the previous subsection, it is clear that the cDs formed in the Rz2 resimulations give a much better match to reality than the ones formed in the Rz1 simulations, and in fact, in general, meet all observational constraints. We shall hence in the following focus on the Rz2 resimulations, and the comparison of these to the DM-only simulations.

### 3.3.1 Resolution limitations due to gravity softening

The Coma DM-only simulation has $\epsilon_{DM} = 2.8 h^{-1} \text{kpc}$, and the gravity force of the central DM particles is hence purely Newtonian (i.e. unsoftened) at $r \simeq 8 \text{kpc}$. At $z = 0.28$, $r_{200}$ for the Coma cluster is 1820 kpc, which means that the simulation result is affected by resolution effects inside of $r/r_{200} \simeq 0.004$. The Coma simulation with baryons has $\epsilon_{DM} = 5.4 h^{-1} \text{kpc}$, so although the star and gas particle softening lengths are considerably smaller, the resulting DM profile may be affected by gravity softening inside of $r/r_{200} \simeq 0.008$. The Virgo simulations with and without baryons have $\epsilon_{DM} = 2.7$ and $2.8 h^{-1} \text{kpc}$, respectively. As $r_{200}$ for this cluster at $z = 0.28$ is 1105 kpc, the simulation results may be affected by gravity softening inside of $r/r_{200} \simeq 0.007$. Finally, the VirgoLR simulations have $\epsilon_{DM} = 5.4$ and $5.6 h^{-1} \text{kpc}$, respectively, resulting in an inner radial resolution limit of $r/r_{200} \simeq 0.014$.

In the following we shall only discuss simulation results outside of the above gravity softening resolution limits.

### 3.3.2 Steepening of the inner dark matter profiles

As can be seen from Fig. 1, the DM density profile obtained in the Coma DM-only simulation can be fairly well described by a generalized NFW profile of $\alpha_{NFW} = \alpha(0) \sim 0.9$. Comparing it to what is obtained for the Rz2 resimulation of Coma the steepening of the inner DM profile, caused by the cD, is quite moderate, $\Delta\alpha \sim 0.1$. For the Virgo cluster, the DM-only profile is also fairly well described by a generalized NFW profile, this time of $\alpha_{NFW} = \alpha(0) \sim 1.1$. Again, the Rz2 resimulation only leads to a very moderate steepening, $\Delta\alpha \sim 0.1$, relative to the DM-only simulation. So for realistic cDs, the steepening of the DM profile in galaxy clusters is very moderate, $\Delta\alpha \sim 0.1$. This is the most important result presented in this paper.

### 3.3.3 Total matter density profiles

As gravitational lensing probes the total matter (surface) density, in order to deduce DM profiles observationally, one has to do a DM versus baryonic matter decomposition. This obviously adds some uncertainty (statistical as well as possibly systematic) to the estimate of the DM density profile. To bypass this, it is clearly also of interest to calculate total matter density profiles from the simulations. These can then more directly be compared to observational finding, based on lensing. In Fig. 5, the result for the Rz2 simulations is shown. Both the Virgo and the Coma profile are characterized by a rather sharp break at $r/r_{200} \simeq 0.05-0.06$. Inside of this the profiles are rather flat, with $\alpha(r) \sim \alpha_{NFW;DMonly} + 0.5$. The profiles are not well described by generalized NFW models, as exemplified in the figure.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Logarithmic density slopes for the Rz2 resimulations. Results for total matter densities are shown by thick, short-dashed lines (red: Virgo; blue: Coma). The thin long-dashed lines show the result for the DM (as in Fig. 1). The two black solid curves show generalized NFW models (equation 1) with $\alpha_{NFW} = 1.0$ and 1.3, and concentration parameters $c = r_{200}/r_s = 4$ and 6, respectively.</figcaption></figure>

### 3.3.4 Dark matter + gas density profiles

In performing the DM versus baryonic matter decomposition mentioned above, one typically assumes that the baryonic component is dominated by the stellar mass of the cD, for which the mass distribution can be fairly well modelled. It is hence important to show that neglecting the hot gas (ICM) component, for which the mass distribution is much more uncertain, does not significantly bias the estimate of the DM distribution. In Fig. 6, for the Rz2 simulations, $\alpha(r)$ for the DM alone, as well as for DM and hot ICM gas combined, is shown. As can be seen, the effect on $\alpha(r)$ of including the ICM is very modest, $\Delta\alpha \lesssim 0.05$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** For the Rz2 resimulations, $\alpha(r)$ for the DM alone (thick short-dashed lines), and for DM and ICM gas combined (thin long-dashed lines) is shown. Results for Coma are shown by blue curves and for Virgo by red curves. The two black curves are the NFW models shown in Fig. 1.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of the Virgo and VirgoLR simulations, at $z = 0$. The red curves display the results for Virgo, labelled in the same way as in Fig. 1. The black curves show the corresponding results for VirgoLR. The statistical uncertainty for the VirgoLR data is $\Delta\alpha \sim 0.05$; for Virgo about a factor of $\sqrt{8}$ less. The light blue curve shows a generalized NFW models (equation 1) with $\alpha_{\text{NFW}} = 1.09$ and concentration parameter $c = r_{200}/r_s = 6$.</figcaption></figure>

## 3.4 A numerical resolution test

In Fig. 7 $\alpha(r)$ for the sets of Virgo simulations as well as VirgoLR simulations is shown. Statistical uncertainties on the VirgoLR results are $\Delta\alpha \sim 0.05$, and on the Virgo results about a factor of $\sqrt{8}$ less. Note that the VirgoLR results are affected by gravity softening at about $r/r_{200} \simeq 0.014$. Outside of this region, there is reasonably good agreement between the two sets of simulations, indicating that the results presented in Fig. 1 are not seriously affected by numerical resolution limitations outside of $r/r_{200} \simeq 0.007-0.008$.

## 4 DISCUSSION

Fitting the DM density profile of Abell 1703 by a generalized NFW profile, different authors find $\alpha(0)$ between 0.8 and 1.2 (Limousin et al. 2008; Richard et al. 2009; Oguri et al. 2009; Saha & Read 2009). Although we in this paper are analyzing simulations of just two clusters, and although other observational values of $\alpha(0)$ have previously been reported in the literature, it is obviously interesting that we for the two clusters simulated, including the effects of the central cD on the DM halo, find $\alpha(0) \simeq 1.0-1.2$.

Given that $\Lambda$CDM haloes undoubtedly have a range of inner slopes, and given that we are probing only two haloes, an even more interesting result is that for both clusters, the steepening of the inner slopes is both moderate and consistent, $\Delta\alpha \sim 0.1$. For Milky Way sized galaxies, as well as smaller galaxies, Gustafsson et al. (2006) found much larger effects of baryonic pinching of the DM haloes, $\Delta\alpha \sim 0.6$. To understand the reason for this difference, we show in Fig. 8, for the Rz2 resimulations at $z = 0.28$, the cumulative masses of DM, stars and stars+gas. As can be seen, the cumulative masses of baryons and DM are equal at $r_{eq} \simeq 7$ and 9 kpc for Virgo and Coma, respectively. In both cases, this is well inside of the effective radius, i.e. the characteristic scale of the galaxy. In contrast, Gustafsson et al. (2006) found $r_{eq} \sim 10$ kpc for the individual (field) galaxies considered. These galaxies have characteristic scales of $\sim 2-4$ kpc, so the chief difference is

<figure><img src="image_10.png" /><figcaption>**Figure 8.** For the Rz2 runs, at $z = 0.28$, cumulative mass of DM (solid curves; Coma: thick, Virgo: thin), stars (dotted curves; Coma: blue, Virgo: red) and stars+gas (dashed curves; Coma: blue, Virgo: red) is shown.</figcaption></figure>

that the cD galaxies considered in this paper are much less baryon dominated than normal (substantially smaller) field galaxies. In addition, one should note that the present simulations can only probe the increase in DM inner profile steepening outside of $\sim 8$ and 16 kpc for Virgo and Coma, respectively (whereas the innermost, baryon-dominated regions are about twice as well resolved due to the shorter gravitational softening lengths of the star and gas particles).

It is also worth noting that G04 found a considerably larger steepening of galaxy cluster inner DM profiles due to baryonic pinching, on the basis of their simulations invoking baryonic physics. This, however (as acknowledged by these authors), is due to the spurious effects of ‘over-cooling’, discussed and corrected for in the present work. G04 in fact carried out a test simulation for which radiative cooling was switched off at $z = 2$. For this experiment, they found that the inner DM profile was still considerably steeper than what was found for the corresponding adiabatic (no radiative cooling) simulation, seemingly in contradiction to our findings in this paper. This, however, is due to a difference in the physical processes invoked in the simulations: for a cluster only slightly larger than our Virgo cluster, G04 find for their ‘standard’ simulation (radiative cooling at all $z$) a cumulative stellar mass at $r = 20$ kpc of about $7 \times 10^{12} M_{\odot}$ (assuming that the bulk of the baryonic mass inside of $r = 20$ kpc is stellar, as is the case for our simulations) at $z = 0$. For their Rz2-type simulation, the corresponding cumulative stellar mass is about $3 \times 10^{12} M_{\odot}$, which is actually more than the $z = 0$ stellar mass inside of $r = 50$ kpc in our ‘standard’ Virgo simulation (cooling at all $z$). Given the findings in this paper (e.g. Fig. 1), it is hence not surprising that G04 find considerable steepening of their inner DM profile, even for their Rz2 simulation. The reason for the difference between the cD masses obtained in our simulations and those of G04 must be ascribed to differences in the physics invoked. Most likely, the significantly lower cD masses obtained in our simulations are a consequence of the strong starburst-driven galactic superwinds invoked here.

An underlying assumption in the Rz1 and Rz2 resimulations is that the rate of change of cD gravitational potential, caused by stars either being formed *in situ* or being added through accretion since $z_f = 1$ or 2, is slow (adiabatic), compared to the dynamical rate.
Although this is true for the stars forming *in situ* at the base of the cooling flows, in the Coma simulation ~1:3 and 1:4 mergers take place at $z \simeq 1.4$ and 0.7, respectively, and in the Virgo simulation ~1:4 mergers take place, at $z \simeq 1.5$ and 1.0, respectively. In such mergers, the adiabatic approximation is not fulfilled, and more realistic simulations, including effects of supermassive black hole growth and AGN feedback, would have to be undertaken to check whether the results we find from the adiabatic stellar mass removal resimulations are in fact fully correct. However, the fact that the Rz2 cDs appear to be realistic in terms of stellar mass, linear extent, metallicity, etc., can be considered very strong evidence that the inner DM profile steepening found for the Rz2 simulations is entirely reasonable. It is obvious that the primary drivers of the baryonic pinching are (i) cD stellar mass, and, to a lesser extent, (ii) the cD linear extent.

## 5 CONCLUSION AND OUTLOOK

We have in this paper presented results of cosmological (ΛCDM) TREESPH simulations of the formation and evolution of two galaxy clusters. In order to quantify the effect of baryonic 'pinching' of galaxy cluster DM haloes, we have performed simulations of the two clusters, with and without baryons included (starting from otherwise identical initial conditions).

The hydrodynamical simulations include a number of important physical processes, such as star formation, chemical evolution with non-instantaneous recycling, metallicity-dependent radiative cooling, strong starburst-driven galactic superwinds and the effects of a meta-galactic UV field, including simplified radiative transfer. The two clusters have $T \simeq 3$ and 6 keV, respectively, and both host at $z \sim 0$ a prominent, central cD galaxy.

Comparing the simulations without and with baryons, it is found that for the latter, the inner DM density profiles, $r \lesssim 50-100$ kpc, steepen considerably: $\Delta\alpha \sim 0.5-0.6$, where $-\alpha$ is the logarithmic DM density gradient. This is primarily due to the central stellar cDs becoming very massive, as a consequence of the onset of late-time cooling flows and related star formation. Once these spurious cooling flows have been corrected for, and the cluster gravitational potentials dynamically adjusted, much smaller pinching effects are found: $\Delta\alpha \sim 0.1$.

The main result of this paper is hence that our work strongly supports the notion that the steepening of the inner DM density profile, due to baryonic pinching, is very moderate. This finding is at stark variance with findings for individual field galaxies. The main reason for this is that galaxy clusters, even in the presence of central dominant cD galaxies, are DM dominated outside of ~7–8 kpc. As the effective radii of the two cDs considered are 11 and 16 kpc, this result indicates that cD galaxies are dark matter dominated at the effective radius.

Including the effects of baryonic pinching, central slopes of $\alpha \simeq 1.0$ and 1.2 are found for the DM in the two clusters, interestingly close to recent observational findings for Abell 1703 based on gravitational lensing studies.

For the simulations with baryons, the inner density profile of DM+ICM combined is found to be only very marginally steeper than that of the DM. The total matter inner density profiles, however, are found to be $\Delta\alpha \sim 0.5$ steeper than the inner profiles in the DM-only simulations. The total matter density profiles are not well described by generalized NFW models.

Given that we have presented results for just two simulated clusters in this paper, it is clearly of importance to increase the sample size. On the observational side, strong lensing constraints for more

unimodal clusters should be obtained. We plan to present such results in forthcoming papers.

## ACKNOWLEDGMENTS

We have benefited considerably from discussions with Peter Johanson, Sune Toft and Signe Riemer-Sørensen. The comments by the anonymous referee considerably helped in improving the presentation of this paper. We gratefully acknowledge abundant access to the computing facilities provided by the Danish Centre for Scientific Computing. This work was supported by the DFG Cluster of Excellence 'Origin and Structure of the Universe'. The Dark Cosmology Centre is funded by the Danish National Research Foundation. ML acknowledges the Centre National d'Etudes Spatiales for their support.

## REFERENCES

* Allen S. W. et al., 1992, MNRAS, 259, 67

* Arabadjis J. S., Bautz M. W., Garmire G. P., 2002, ApJ, 572, 66

* Austin C. G., Williams L. L. R., Barnes E. I., Babul A., Dalcanton J. J., 2005, ApJ, 634, 756

* Bingelli B., Sandage A., Tammann G. A., 1985, AJ, 90, 1681

* Biviano A., Salucci P., 2006, A&A, 452, 75

* Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27

* Bower R. G., Benson A. J., Malbon R., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

* Bradač M. et al., 2007, preprint (arXiv e-prints, 0711.4850)

* Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

* Böhringer H. et al., 2000, ApJS, 129, 435

* Cleary P. W., Monaghan J. J., 1999, J. Comput. Phys., 148, 227

* Croton D. J. et al., 2006, MNRAS, 367, 864

* Dahle H., Hannestad S., Sommer-Larsen J., 2003, ApJ, 588, L73

* D'Onghia E., Sommer-Larsen J., Romeo A. D., Burkert A., Pedersen K., Portinari L., Rasmussen J., 2005, ApJ, 630, L109

* Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, in press (arXiv:1004.1851)

* El-Zant A., Shlosman I., Hoffman Y., 2001, ApJ, 560, 636

* Elíasdóttir, Á. et al., 2007, preprint (arXiv:0710.5636)

* Ettori S., Fabian A. C., 2000, MNRAS, 317, 57

* Ettori S., Fabian A. C., Allen S. W., Johnstone R. M., 2002, MNRAS, 331, 635

* Gao L. et al., 2008, MNRAS, 387, 536

* Gavazzi R., 2005, A&A, 443, 793

* Gavazzi R., Fort B., Mellier Y., Pelló R., Dantel-Fort M., 2003, A&A, 403, 11

* Ghigna S. et al., 2000, ApJ, 544, 616

* Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16 (G04)

* Graham A. W., Merritt D., Moore B., Diemand J., Terzic B., 2006, AJ, 132, 2701

* Gustafsson M., Fairbairn M., Sommer-Larsen J., 2006, Phys. Rev. D., 74, 123522

* Haardt F., Madau P., 1996, ApJ, 461, 20

* Hansen S. H., Stadel J., 2006, J. Cosmol. Astropart. Phys., 5, 14

* Johansson P., Naab T., Ostriker J. P., 2009, ApJ, 697, L38

* Kay S. T., Pearce F. R., Frenk C. S., Jenkins A., 2002, MNRAS, 330, 113

* Kazantzidis S., Magorrian J., Moore B., 2004, ApJ, 601, 37

* Kelson D. D. et al., 2002, ApJ, 576, 720

* Kneib J. et al., 2003, ApJ, 598, 804

* Lewis A. D., Buote D. A., Stocke J. T., 2003, ApJ, 586, 135

* Lia C., Portinari L., Carraro G., 2002a, MNRAS, 330, 821

* Lia C., Portinari L., Carraro G., 2002b, MNRAS, 335, 864

* Limousin M. et al., 2007b, ApJ, 668, 643

* Limousin M. et al., 2008, A&A, 489, L23

* Ma C.-P., Boylan-Kolchin M., 2004, Phys. Rev. Lett., 93, 021301
Mandelbaum R., Seljak U., Hirata C. M., 2008, J. Cosmol. Astropart. Phys., 8, 6

McNamara B. R., 2004, in Ruprecht T. H., Kempner J. C., Soker N., eds, The Riddle of Cooling Flows in Galaxies and Clusters of Galaxies

Michard R., Andreon S., 2008, A&A, 490, 923

Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ, 499, L5

Muanwong O., Kay S. T., Thomas P. A., 2006, ApJ, 649, 640

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Navarro J. F. et al., 2004, MNRAS, 349, 1039

Nipoti C., Stiavelli M., Ciotti L., Treu T., Rosati P., 2003, MNRAS, 344, 748

Nipoti C., Treu T., Ciotti L., Stiavelli M., 2004, MNRAS, 355, 1119

Oguri M. et al., 2009, ApJ, 699, 1038

Okabe N., Takada M., Umetsu K., Futamase T., Smith G. P., 2010, PASJ, 62, 811

Orban de Xivry G., Marshall P., 2009, MNRAS, 399, 2

Pointecouteau E., Arnaud M., Pratt G. W., 2005, A&A, 435, 1

Richard J. et al., 2007, ApJ, 662, 781

Richard J., Pei L., Limousin M., Jullo E., Kneib J., 2009, A&A, 498, 37

Ricotti M., 2003, MNRAS, 344, 1237

Riemer-Sørensen S., Paraficz D., Ferreira D. D. M., Pedersen K., Limousin M., Dahle H., 2009, ApJ, 693, 1570

Romeo A. D., Portinari L., Sommer-Larsen J., 2005, MNRAS, 361, 983

Romeo A. D., Sommer-Larsen J., Portinari L., Antonuccio-Delogu V., 2006, MNRAS, 371, 548

Romeo A. D., Napolitano N. R., Covone G., Sommer-Larsen J., Antonuccio-Delogu V., Capaccioli M., 2008, MNRAS, 389, 13

Sand D. J., Treu T., Ellis R. S., 2002, ApJ, 574, L129

Sand D. J., Treu T., Smith G. P., Ellis R. S., 2004, ApJ, 604, 88

Sand D. J., Treu T., Ellis R. S., Smith G. P., Kneib J.-P., 2008, ApJ, 674, 711

Saha P., Read J. I., 2009, ApJ, 690, 154

Saro A., Borgani S., Tornatore L., Dolag K., Murante G., Biviano A., Calura F., Charlot S., 2006, MNRAS, 373, 397

Shen S. et al., 2003, MNRAS, 343, 978

Smith G. P., Kneib J.-P., Ebeling H., Czoske O., Smail I., 2001, ApJ, 552, 493

Smith G. P. et al., 2005, MNRAS, 359, 417

Sommer-Larsen J., Götz M., Portinari L., 2003, ApJ, 596, 46

Sommer-Larsen J., Romeo A. D., Portinari L., 2005, MNRAS, 357, 478

Springel V., Hernquist L., 2002, MNRAS, 333, 649

Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776

Tamura N., Kobayashi C., Arimoto N., Kodama T., Ohta K., 2000, AJ, 119, 2134

Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2010, MNRAS, submitted (arXiv:1003.4744)

Thomas P. A., Muanwong O., Kay S. T., Liddle A. R., 2002, MNRAS, 330, L48

Thomas D., Maraston C., Bender R., 2005, ApJ, 621, 673

Tornatore L., Borgani S., Matteucci F., Recchi S., Tozzi P., 2004, MNRAS, 349, L19

Tornatore L., Borgani S., Dolag K., Matteucci F., 2007, MNRAS, 382, 1050

Tyson J. A., Kochanski G. P., dell'Antonio I. P., 1998, ApJ, 498, L107

Valdarnini R., 2003, MNRAS, 339, 1117

van Dokkum P. G. et al., 2008, ApJ, 677, L5

Verdugo T., de Diego J. A., Limousin M., 2007, ApJ, 664, 702

Zappacosta L. et al., 2006, ApJ, 650, 777

This paper has been typeset from a *TEX/LATEX* file prepared by the author.