Impact of baryon physics on dark matter structures: a detailed simulation study of halo density profiles

The back-reaction of baryons on the dark matter halo density profile is of great interest, not least because it is an important systematic uncertainty when attempting to detect the dark matter. Here, we draw on a large suite of high resolution cosmological hydrodynamical simulations, to systematically investigate this process and its dependence on the baryonic physics associated with galaxy formation. The inclusion of baryons results in significantly more concentrated density profiles if radiative cooling is efficient and feedback is weak. The dark matter halo concentration can in that case increase by as much as 30 (10) per cent on galaxy (cluster) scales. The most significant effects occur in galaxies at high redshift, where there is a strong anti-correlation between the baryon fraction in the halo centre and the inner slope of both the total and the dark matter density profiles. If feedback is weak, isothermal inner profiles form, in agreement with observations of massive, early-type galaxies. However, we find that AGN feedback, or extremely efficient feedback from massive stars, is necessary to match observed stellar fractions in groups and clusters, as well as to keep the maximum circular velocity similar to the virial velocity as observed for disk galaxies. These strong feedback models reduce the baryon fraction in galaxies by a factor of 3 relative to the case with no feedback. The AGN is even capable of reducing the baryon fraction by a factor of 2 in the inner region of group and cluster haloes. This in turn results in inner density profiles which are typically shallower than isothermal and the halo concentrations tend to be lower than in the absence of baryons.


INTRODUCTION
Increasingly powerful computers, combined with ever more efficient N-body codes, have permitted accurate numerical tests of structure formation in a Cold Dark Matter (CDM) universe (e.g. Klypin et al. 1996;Navarro et al. 1997; Bullock et al. 2001;Diemand et al. 2008;Gao et al. 2008). It is now an established prediction that CDM haloes form a cuspy mass distribution that is close to a 'universal profile', independent of halo mass (Navarro et al. 1996;Navarro et al. 1997, hereafter NFW). In detail, however, the haloes are not strictly self-similar (Navarro et al. 2010). Foremost amongst these results has been the existence of power law relationships in the phase-space density of haloes across several decades in radius (Taylor & Navarro 2001) and the halo concentrationmass relation over five decades in mass (Bullock et al. 2001;Neto et al. 2007;Macciò et al. 2008;Duffy et al. 2008). High resolution N-body simulations have also shown that DM haloes are not smooth, but contain self-bound substructures that have survived the accretion process (e.g. Klypin et al. 1999). More recently, it has been shown that substructures themselves contain substructures, a pattern that shows no signs of abating (Stadel et al. 2009;. These power laws and embedded substructures are the product of an approximately scale-free system that interacts only via gravity and, for scales large enough that gravity is the only dynamically significant force, can be viewed as an accurate approximation to structure formation. On smaller scales, other processes significantly modify the structure of the baryons, which can subsequently affect the DM in the densest regions. In particular, it has long been known that a galaxy forms because of the ability of its progenitor gas to cool efficiently via emission of radiation (e.g. Hoyle 1953). Such a process introduces a characteristic physical scale into the system and thus breaks the self-similarity of structures (White & Rees 1978). The net effect is that the cooled gas and stars then cluster on smaller scales than the DM, thereby deepening the potential well of the system and causing the DM to orbit closer to the centre. This will modify a number of N -body simulation predictions on small scales, in a way that depends on the efficiency of galaxy formation.
This modification of the DM distribution is usually parametrised in the form of adiabatic invariants, hence the term 'adiabatic contraction' that is used to describe the overall effect, first considered by Eggen et al. (1962) and Zeldovich et al. (1980). The adiabatic parameter considered in the first generation of adiabatic contraction models for collapsing haloes (Blumenthal et al. 1986) was the product of radius with the mass internal to this radius, M (< r)r. Gnedin et al. (2004), henceforth G04, found that such a simple scheme, only strictly applicable for particles on circular orbits, did not provide a good description of what was happening in their hydrodynamic simulations of groups and clusters at z ≈ 0 and galaxies at z > 3. Better agreement was found when they modified the parameter to use the orbit-averaged radius,r, that is a power-law function of radius. The work of G04 was extended to galaxy scales at low redshift by Gustafsson et al. (2006), who modified the slope and normalisation of this power law distribution to provide a better description of their results at this lower mass range. The compression of the DM by the infalling baryons has also been extensively studied by purely analytic methods (e.g. Sellwood & McGaugh 2005;Cardone & Sereno 2005;Klar & Mücket 2008).
In addition to the contraction of the DM halo, the concentration of baryonic mass at its centre also leads to a more spherical structure (Kazantzidis et al. 2004). This effect is due to the more spherical gravitational potential formed by the baryons and it is this that induces an adiabatic decrease in the eccentricity of the orbits of the DM and stars, rather than through the destruction of box orbits . Additionally, in the presence of a gaseous disc, in-falling subhaloes with low orbital inclinations with respect to the disc will experience significant disruption due to the high local density of the baryons. This will lead to the accumulation of stars, and DM, in the same plane, forming both a thin stellar and a DM disc .
With the addition of baryons to the simulation, one can expect further effects due to the transfer of angular momentum from infalling satellites (e.g. Debattista et al. 2008;Romano-Díaz et al. 2008). However, the influence of baryons on these smaller, less well-resolved objects, are by no means certain. For example, the substructure may exhibit greater resistance to tidal disruption due to its deeper potential well in the presence of baryons, while at the same time the substructure will typically suffer from increased dynamical friction with the main halo (Jiang et al. 2008). Recent work by Macciò et al. (2006) found a factor of two increase in the number of surviving satellites in a galaxysized halo, in a hydrodynamical simulation relative to a DMonly simulation. This increased survival rate in the inner regions (Weinberg et al. 2008) may enable satellite infall to partly counteract the contraction of the DM halo, as these objects can efficiently transfer angular momentum to the inner parts of the DM halo, instead of being tidally disrupted at larger radii (for recent considerations of this effect see Pedrosa et al. 2009, Pedrosa et al. 2010and Abadi et al. 2009). The erasure of the central DM cusp in galaxies by infalling satellites was proposed by El-Zant et al. (2001) and for the case of clusters by El-Zant et al. (2004).
The effect of the baryons on the inner DM density profile is of particular interest, as it is currently a significant source of uncertainty when making predictions for DM detection experiments. For example, the expected γ-ray signal from self-annihilation of potential DM candidate particles , which may become detectable in the near future using γ-ray observatories such as Fermi 1 , is proportional to ρ 2 DM . We further note that the total mass density profile is also of interest, as it is this quantity that is determined from strong lensing studies and is relevant for direct comparison with such observations (for a review of the methodologies and uses of strong lensing, see Kochanek 2006).
The incorporation of baryons in a simulation is a challenging theoretical undertaking as the relevant scales of the physical processes are rarely resolved in a cosmological simulation and approximations are therefore necessary. The method by which stellar feedback, for example, is incorporated can have a large influence on the resulting baryonic distribution of the galaxy (e.g. . For this reason, we have attempted to probe these effects in a systematic way by examining the DM halo density profiles from galaxies, groups and clusters, drawn from an extensive series of cosmological simulations with that incorporate a wide variety of prescriptions for baryonic processes. These form a subset of the larger, overall simulation series known as the OverWhelmingly Large Simulations project (OWLS; Schaye et al. 2010).
Our work has a number of novel features that allow us to extend recent studies that have probed the dynamical response of the DM halo profile to the presence of baryons. For example, while Gustafsson et al. (2006) examined a small number of galaxy-sized haloes at z = 0, we provide results on both group and cluster mass scales, with significantly better statistics. G04 examined the DM profiles in a similar mass range to our own, albeit limited to higher redshift for galaxy scales, and with similar mass resolution. They demonstrate the significant impact that metal enrichment and stellar feedback have on the overall gas distribution, although they did not vary the models systematically. The recent work of Pedrosa et al. (2009) investigated the effects of physics on galaxy scales by simulating a single halo. Here we provide a more statistically-robust set of results, due to the large number of haloes we have resolved. (We note that a direct comparison with their z = 0 galaxy is unfortunately not yet possible and must await further simulations). In summary, the large dynamic range in mass, from dwarf galaxies to clusters, across a wide variety of physics implementations and with high statistical significance, is a significant advance in the study of baryons in DM haloes.
The rest of the paper is organised as follows. In Section 2 we describe the range of simulations used in this study, along with the differences between implementations of gas physics. In Section 3 we investigate the distribution of the baryons within these haloes as a function of the physics implementation. In Section 4 we compare our haloes with both strong gravitational lensing observations at high redshift and stellar mass fractions at the present day. We discuss the implications of our findings for structure formation in the CDM model, with the proviso that certain simulations can fulfil some observations but none can be used at all mass scales. The analysis of the haloes formed in these simulations, via their density profiles, is discussed in Section 5. We then attempt to parametrise the DM density profile with existing theoretical profiles, and the total matter profile by comparing non-parametric measures of halo concentrations, in Section 6. In Section 7 we test whether we can mimic the influence of the baryons by making use of the adiabatic contraction models of Blumenthal et al. (1986) and G04. Finally, we summarise our work in Section 8.

SIMULATIONS
This work is based on OWLS; a series of high-resolution simulations of cosmological volumes with differing subgrid physics implementations (Schaye et al. 2010). These simulations were run using a modified version of gadget-3, itself a modified version of the publicly-available gadget-2 code (Springel 2005). gadget calculates gravitational forces using the Particle-Mesh algorithm (e.g. Klypin & Shandarin 1983) on large scales, supplemented by the hierarchical tree method (Barnes & Hut 1986) on smaller scales. Hydrodynamic forces are computed using the Smooth Particle Hydrodynamics (SPH) formalism (Monaghan & Lattanzio 1985;Springel & Hernquist 2002).
The production simulations available in OWLS model cubes of comoving length 25 and 100 h −1 Mpc with 512 3 gas and 512 3 DM particles. For all runs, glass-like cosmological initial conditions were generated at z = 127 using the Zel'dovich approximation and a transfer function generated using cmbfast (v. 4.1; Seljak & Zaldarriaga 1996). Cosmological parameters were set to the best-fit values from the 3rd year Wilkinson Microwave Anisotropy Probe data (Spergel et al. 2007), henceforth known as WMAP3. Specifically, the values for [Ωm, Ω b , ΩΛ, h, σ8, ns] were set to [0.238, 0.0418, 0.762, 0.73, 0.74, 0.95]. The primordial mass fraction of He was set to 0.248. All runs that used the same box size had identical initial conditions.
We have analysed a subset of 6 OWLS runs for each box size: a DM-only run (hereafter labelled DMONLY) and 5 additional runs that also follow the baryonic component with hydrodynamics, gas cooling and star formation. For DMONLY, the particle masses are m = 7.7 × 10 6 and 4.9 × 10 8 h −1 M⊙ for the 25 and 100 h −1 Mpc runs, respectively. For the runs with baryons, each particle was split in two in the initial conditions, with the mass shared according to the universal baryon fraction, f univ b = Ω b /Ωm = 0.176, such that the DM (gas) particle mass is 6.3 (1.4) × 10 6 and 4.1 (0.86) × 10 8 h −1 M⊙ for the 25 and 100 h −1 Mpc runs, respectively. The 25 h −1 Mpc simulations were run to z = 2 while the 100 h −1 Mpc simulations were run all the way to z = 0 (our results will focus on data from these redshifts). At early times the softening length was held constant in comoving co-ordinates at 1/25 of the initial mean inter-particle spacing. Below z = 2.91 the softening length was held fixed in proper units (thus, the z = 0 values were 0.5 h −1 kpc and 2 h −1 kpc respectively).

Baryonic physics
For our purposes two main aspects of the baryonic physics have been varied: radiative cooling (with and without metal lines) and feedback (from supernovae and accreting supermassive black holes). The modification of the baryon distribution will be primarily determined by these two competing effects. Within OWLS other physics prescriptions have been varied, such as the choice of stellar initial mass function for example, however these will at most have a secondary impact on the baryon distribution. We use 'PrimC' and 'ZC' to label two different approaches to cooling, denoting the absence-or inclusion of metals to the cooling rate respectively. For the supernova feedback, we use 'NFB', 'WFB' and 'SFB' to identify whether stellar feedback is absent, weak or strong. In addition, we have analysed a simulation that includes feedback due to accreting supermassive black holes, associated with Active Galactic Nuclei; this run is labelled 'AGN'. Table 1 lists the runs with baryonic physics, along with their original OWLS names (Schaye et al. 2010), and highlights their key differences. We only provide brief descriptions here, for further details please see Schaye et al. (2010) and the references given below.
In all the simulations with baryons, radiative cooling (and heating) rates are implemented element-by-element using tables generated with the cloudy radiative transfer code (Ferland et al. 1998), following the method described in Wiersma et al. (2009). We assume collisional equilibrium before reionisation (z > 9) and photoionisaton equilibrium afterwards, in the presence of an evolving UV/X-ray background (Haardt & Madau 2001). In the 'PrimC' runs cooling rates were computed assuming primordial abundances. In those labelled 'ZC' we also track line emission from nine metals: C, N, O, Ne, Mg, Si, S, Ca and Fe. Both types of simulation also include cooling via free-free Bremsstrahlung emission and Compton cooling due to interactions between Table 1. A list of all simulations, the corresponding names in the OWLS project (Schaye et al. 2010)  When gas is sufficiently dense it is expected to be multiphase and unstable to star formation (Schaye 2004). Our simulations lack both the resolution and the physics to model the cold interstellar gas phase and we therefore impose an effective equation of state, P ∝ ρ γ eff , for densities that exceed our star formation threshold of nH > 0.1 cm −3 . The normalisation is fixed such that P/k = 1.08 × 10 3 K cm −3 at the star formation density threshold. The index is set to γ eff = 4/3, which has the advantage that both the Jeans mass and the ratio of the Jeans length to the SPH smoothing length are independent of density, thus suppressing spurious fragmentation .  showed how the observed Kennicutt-Schmidt law,Σ * = A(Σg/1 M⊙ pc −2 ) n can be analytically converted into a pressure law, which can be implemented directly into the simulations. This has the advantage that the parameters are observables and that the simulations reproduce the input star formation law irrespective of the assumed equation of state. We use the  method, setting A = 1.515× 10 −4 M⊙ yr −1 kpc −2 and n = 1.4 (Kennicutt 1998, note that we renormalised the observed relation for a Chabrier IMF).
These star particles are assumed to be simple stellar populations with an initial mass function given by Chabrier (2003). The metal abundances of the stellar particles are inherited from their parent gas particles. The stellar evolution, and the subsequent release of metals from the star particle, depends on this metallicity (Portinari et al. 1998;Marigo 2001;Thielemann et al. 2003). The delayed release of 11 individual elements (the ones used for the calculation of the cooling rates) by massive stars, from Type Ia and Type II supernovae, as well as AGB stars, is tracked. Every time step they are shared amongst neighbouring gas particles according to the SPH interpolation scheme. For further details please see Wiersma et al. (2009).
Stars formed in the simulations labelled 'WFB' and 'SFB' inject energy into local gas particles. This energy is in kinetic form, i.e. nearby gas particles are kicked away from the stars. On average, each newly formed star particle kicks η times its own mass, where η is the dimensionless mass loading parameter, by adding a randomly oriented velocity vw to the velocity of each kicked particle. The simulations labelled 'WFB' use η = 2 and vw = 600 km s −1 , which corresponds to forty per cent of the SN energy for our initial mass function. Further details of the feedback prescription can be found in Dalla . Note that although we have labelled the supernova feedback in this simulation "weak", it is in fact strong compared to many prescriptions in the literature.
The runs labelled 'SFB' implement wind velocity and mass-loading parameters that depend on the local density of the gas, vw = 600 km s −1 nH/0.1 cm −3 1/6 and η = 2 vw/600 km s −1 −2 , such that in the densest regions winds are launched with the largest speeds and smallest massloadings. For gas that follows the imposed effective equation of state, this implies a wind speed that is proportional to the effective local sound speed, vw ∝ c s,eff ∝ (P/ρ) 1/2 . The result is that winds in the 'SFB' model remove gas from higher mass systems more efficiently than in the 'WFB' model. Note that the WFB and SFB models use the same amount of SN energy. The difference in the efficiency thus results purely from energy is distributed between mass-loading and wind velocity.
Finally, in the simulation labelled 'AGN' supermassive black holes (BHs) are grown, with subsequent feedback, at the centres of all massive haloes according to the methodology of , which is a substantially modified version of that introduced by . Seed BHs of mass m seed = 9 × 10 4 M⊙ (i.e. 10 −3 mg in the 100 h −1 Mpc box, where mg is the initial mass of the gas particles), are placed into every DM halo more massive than 4 × 10 10 M⊙ (i.e. 10 2 DM particles in the 100 h −1 Mpc box). Haloes are identified by regularly running a Friendsof-Friends (FOF) group finder on-the-fly during the simulation. After forming, BHs grow by two processes: accretion of ambient gas and mergers. Gas accretion occurs at the minimum of the Eddington rate and the Bondi-Hoyle (1944) rate, where the latter is multiplied by (nH/10 −1 cm −3 ) 2 for starforming gas (i.e. nH > 10 −1 cm −3 ) to compensate for the fact that our effective equation of state strongly underestimates the accretion rate if a cold gas phase is present. Gas accretion increases the BH masses asṁBH = (1 − ǫr)ṁaccr, where ǫr = 0.1 is the assumed radiative efficiency. We assume that 15 per cent of the radiated energy (and thus 1.5 per cent of the accreted rest mass energy) is coupled thermally to the surrounding medium.  showed that this model reproduces the redshift zero cosmic BH density as well as the observed relations between BH mass and both the stellar mass and the central stellar velocity dispersion.

Halo definitions and density profiles
The principal quantity that is extracted from each simulation is the spherically-averaged halo density profile, dom-inated in all but the central regions by the dark matter. Haloes were first identified using the FOF algorithm, that links dark matter particles using a dimensionless linkinglength, b = 0.2 (Davis et al. 1985). We then used the subfind algorithm (Springel et al. 2001;Dolag et al. 2009) to decompose the FOF group into separate, self-bound substructures, including the main halo itself. Finally, a sphere was placed at the centre of each halo, where the centre was identified with the location of the minimum potential particle in the main halo. The spheres were then grown until a specified mean internal density contrast was reached.
A halo thus consists of all mass, M∆, within the radius, R∆, for which where the mean internal density is ∆ times the critical density, ρcrit(z) = 3H(z) 2 /8πG. If M∆ = Mvir and R∆ = Rvir, then one can compute ∆ from the spherical top-hat collapse model. In this case, we adopt the fitting formula from Bryan & Norman (1998), that depends on both cosmology and redshift where x = Ωm(z) − 1, and H(z) is the usual Hubble parameter for a flat universe In our cosmology at z = 0 (2) ∆ = 92.5 (168.5). We will also consider overdensities that are integer multiples of the critical density; ∆ = 500 and 2500, with radii and masses denoted in the standard way, e.g. R500 and M500 for ∆ = 500. Following Neto et al. (2007), only haloes with more than 10 4 DM particles within Rvir are initially selected. Density profiles are then computed using a similar method to our previous work on DM-only haloes (Duffy et al. 2008). To briefly summarise, we take 32 uniform logarithmic shells of width ∆ log 10 (r) = −0.078, in the range −2.5 log 10 (r/Rvir) 0. We then sum the mass within each bin for the three components (gas, stars and DM) and divide by the volume of the shell. We fit model density profiles by minimising the difference of the logarithmic densities between halo and model, assuming equal weighting (this was shown by Neto et al. 2007 to be a self-consistent weighting scheme).
A key issue is the range of radii over which we can trust that the profiles have converged numerically. To check this, we have performed a convergence study, the results of which are presented in Appendix A. In summary, we show that the minimum radius proposed by Power et al. (2003), henceforth P03, is appropriate for our simulations, even for the ZC WFB model (P03 defined their resolution limit on the two-body dynamical relaxation timescale in the context of DM-only simulations). All haloes used in this work satisfy the criterion RP03 < 0.05 Rvir, so we adopt the latter value as the minimum radius for our density profiles.
A selection of well-resolved haloes are singled out for investigation of the slopes of their inner profiles (Section 5.3).
For these haloes, we place even higher demands on the particle number and ensure that haloes satisfy RP03 < 0.025 Rvir and use the latter value as the minimum radius in this instance. In practice this criterion is applied to the DMONLY simulation to define a lower (virial) mass limit: at z = 0 this is 3 × 10 13 h −1 M⊙ and at z = 2 it is 5 × 10 11 h −1 M⊙. The other runs are then matched 2 to the DMONLY haloes that make this cut. All haloes used to study the inner profile have at least 6 × 10 4 DM particles within Rvir.
Between the two simulation sets, we have selected a total of 552 (204) haloes at z = 0 (2) in DMONLY and of these, 67 (32) meet the more stringent resolution demands for the inner profile study. The difference in halo numbers in the different gas physics simulations is less than 10 per cent, primarily due to low-mass haloes having fewer than 10 4 DM particles and hence missing the initial cut. In the subsequent discussion we will use the terms; Dwarf Galaxy, Galaxy, Group, and Cluster, to correspond to Mvir[h −1 M⊙] in the fixed ranges [< 10 12 ,10 12 − 10 13 ,10 13 − 10 14 ,> 10 14 ] respectively.

HALO BARYON DISTRIBUTION
The ability of baryons to cool radiatively and thereby cluster on smaller scales than the DM is a recurring theme in this work, leading to, potentially, numerous modifications of the standard N -body predictions. The exact distribution of the baryons in the halo will be sensitive to the physics implementation included and hence we should first quantify this key quantity. To that end we have plotted, in Fig. 1, the baryon fraction, f b , inside the virial radius against that within r/Rvir = 0.05. We show the median value and the quartile scatter of f b , within mass bins at the different redshifts. Only haloes that meet our more stringent inner profile criteria are shown.
A number of interesting features are apparent from the figure. Firstly, we see that for Groups and Clusters at z = 0 (bottom panel) the baryon fraction within the virial radius is close to the universal value in all runs except ZC WFB AGN. The other runs are approximately consistent with, though slightly higher than, the prediction from non-radiative gas simulations (Crain et al. 2007). The deep gravitational potential wells and long cooling times in these objects, lead to little gas being expelled from the halo. Only in ZC WFB AGN, which includes AGN feedback and for which the feedback is therefore expected to be stronger on group and cluster scales, does a significant amount (∼ 10−40 per cent) of gas get expelled. It should be noted, however, that a significant fraction of this gas will have been ejected from smaller systems at higher redshift.
The inner baryon fractions, f b (r/Rvir < 0.05), typically decrease as the strength of the feedback increases. Note that ZC WFB AGN has an inner baryon fraction that is lower than the universal value. Runs with weak or no feedback have values that are more than twice as high in the run Figure 1. The panels show the baryon fraction within the virial radius as a function of the baryon fraction within r < 0.05 R vir , at z = 2 (top panel) and z = 0 (bottom panel). The relative values of these quantities gives a useful overview of the effect of cooling and feedback on the baryon distribution in haloes at different redshifts and virial masses. Each point corresponds to median values for haloes within a given mass bin, indicated by the symbol size, while the error bars indicate the quartile scatter. The solid vertical and horizontal lines correspond to f univ b . The dashed horizontal line is from the non-radiative gas simulations of Crain et al. (2007). Galaxies are more sensitive to the various feedback and gas cooling schemes than groups and clusters. For groups and clusters the sensitivity to the feedback and cooling is mostly limited to the inner regions but these processes affect galaxies all the way out to the virial radius. AGN feedback reduces the baryon fractions most strongly.
with AGN feedback. The inclusion of metal-line cooling increases the inner baryon fraction, as can be seen by comparing ZC WFB with PrimC WFB.
The situation for Galaxies at z = 2 (top panel) is very different. Only the model with no supernova feedback, PrimC NFB, has a median baryon fraction within the virial radius that is close to f univ b ; even the weak feedback model is capable of expelling some gas. Note that PrimC NFB has a baryon fraction that is higher than predicted from nonradiative simulations, suggesting that gas cooling in this model has compensated for the heating and expansion of the gas due to energy transfer from the dark matter when the halo collapsed. In Galaxy mass haloes, the gravitational potential is sufficiently low that SN feedback is able to unbind gas from the halo. As a result, there is a strong positive Figure 2. To test the realism of our models we have compared the median stellar mass fractions within R 500 , in units of the universal baryon fraction, as a function of total mass within R 500 , at z = 0. The coloured points represent the various physics prescriptions, apart from the single black data point at low halo mass which is from Conroy et al. (2007) for a sample of isolated L * galaxies. The vertical and horizontal error bars represent, respectively, the quartile spread and mass range within each bin. The solid and dashed curves illustrate observational determinations from Giodini et al. (2009) and Lin et al. (2004), with the outer curves representing the uncertainty in the normalisation of their fits. It is clear that if the gas is allowed to cool by metal-line emission, it must be prevented from forming stars by either a supernova model that targets dense regions of the halo, ZC SFB, or by feedback from AGN, ZC WFB AGN. We note that although estimates of stellar masses are typically uncertain by a factor of 2-3 (e.g. Küpcü Yoldaş et al. 2007;Longhetti & Saracco 2009) the observations still favour the strong feedback prescriptions. correlation between f b (r/Rvir < 1) and f b (r/Rvir < 0.05). The galaxies in both ZC WFB and PrimC NFB have inner regions that are baryon-dominated (f b > 0.5).

COMPARISON WITH OBSERVATIONS
In the previous section we compared simulations run with a number of widely used physical prescriptions, and demonstrated that a large range of inner baryon fractions are potentially possible. We now determine which of these models, if any, are compatible with observations. We consider two probes that offer complementary and contrasting constraints on the baryon physics. First, in Section 4.1, we present a comparison of the predicted redshift zero stellar mass fractions to those observed in Galaxies, Groups and Clusters. Secondly, in Section 4.2, we compare predicted total density profiles to those inferred for a high redshift galaxy sample through strong lensing measurements.

Stellar fractions
The stellar fraction by halo mass is well constrained observationally therefore a comparison of median stellar fraction (within R500) to observation allows us to quickly determine the relative realism of our models. In Fig. 2 we show the stellar fraction, in units of M500, for haloes at z = 0 as a function of mass for our full range of physics models.
To compare with the observations, we have taken the group and cluster samples of Lin et al. (2004) and Giodini et al. (2009). Note that the best fitting power law results we take from these works do not contain the contribution from the diffuse intracluster light, and can hence be viewed as lower limits. This additional contribution likely ranges between 11 and 22 per cent (e.g. Zibetti et al. 2005;Krick & Bernstein 2007). Following Balogh et al. (2008), we convert the observed stellar luminosities from Lin et al. (2004) to masses, assuming the best-fit stellar mass-to-light ratio, M/LK = 0.9. Our simulations adopted a Chabrier (2003) IMF, whereas Balogh et al. (2008) used the Kroupa (2001) IMF; we expect this difference to be unimportant, however, as they are very similar over the relevant range of stellar masses. For Giodini et al. (2009) we have adopted the best fitting stellar fraction for the X-ray COSMOS groups/poor clusters only. Also, we reduced the stellar mass fraction by 30 per cent (Longhetti & Saracco 2009) to account for their use of a Salpeter (1955) IMF. Additionally, we include a result on Galaxy scales at z = from Conroy et al. (2007) who determined stellar masses from the spectroscopic SDSS value-added catalogue (Blanton et al. 2005), along with halo mass estimates utilising the measured velocity dispersions of the satellites (we have scaled their halo mass, from M200 to M500, assuming an NFW profile with concentration given by Duffy et al. 2008). We further assume that the stars are significantly more concentrated than the DM and that modifying the mass cut from R200 to R500 will have a negligible effect on the stellar mass. We note that the result of Conroy et al. (2007) is in agreement with Mandelbaum et al. (2006) who made use of weak lensing to determine halo masses.
The trend seen in our simulations, is a stellar fraction that decreases gently as M500 increases from 10 13 to 10 14 M⊙. This is in accord with the semi-analytic model results in Balogh et al. (2008). Simulations with weak or no feedback contain stellar masses that are 2-3 times higher than observed, PrimC WFB is an exception. This difference is expected as the dominant cooling channel at the virial temperatures of haloes with mass ≈ 10 13 M⊙ is believed to be from metal-line emission and hence 'ZC' will have enhanced cooling rates over 'PrimC'. The simulations with stronger feedback, ZC SFB and ZC WFB AGN, are in better agreement with the observations.
We can see that simulations which have relatively efficient supernova feedback in the dense regions of haloes, or include an additional heating source in the form of an AGN, predict more realistic stellar fractions in present day Groups and Clusters with M500 10 13 M⊙. If we include the full mass range at low redshift, we can conclude that AGN feedback is necessary to prevent stellar mass building up in Galaxies, Groups and Clusters. This is in good agreement with McCarthy et al. (2009), who found that model ZC WFB AGN predicts group K-band luminosities and gas fractions that are in excellent agreement with observations. However, in the following section we will see that at high redshift the baryon-dominated inner halo of the schemes with weak or no feedback is required to match strong lensing results, leading to an intriguing conflict which we will discuss further in Section 8. . Inner slope of the total mass density profile versus central baryon fraction. Only our results for Galaxy haloes at z = 2 are shown, to compare with the observational constraint on lensing galaxies (Koopmans et al. 2006; the vertical size of the, lower, maroon box indicates the intrinsic 1σ Gaussian spread across all redshifts while the horizontal size of the box is the 68 per cent confidence interval of the quoted stellar fractions, estimated by bootstrapping. Note that for the sample of early type galaxies the stellar fraction is essentially the baryon fraction). The different symbol sizes denote different mass ranges. Symbol type and colour are used to distinguish different simulations. The black, hatched region indicates the quartile spread of the DMONLY simulation. As is clear from this figure, only simulations with high central baryon fractions reproduce the observed steep inner density profiles of high redshift Galaxies. This is in contrast with the preferred simulation schemes at higher masses and lower redshift (Fig. 2).

Observed inner profile slopes
Gravitational lensing of light by an intervening galaxy enables us to probe the inner mass profile of the lens galaxy's halo. The slope of the inner mass profile for a sample of lens galaxies at z < 1 was presented by Koopmans et al. (2006). Surprisingly, the inner slope is strongly constrained to be close to isothermal (β ≡ d ln ρ d ln r ≈ −2) with no evidence for evolution. The region where the slope is measured (the Einstein radius is typically around 3 h −1 kpc) is comparable to the scales accessible in our highest-resolution simulations at z = 2, allowing us to compare the two results. Fig. 3 shows the inner slope of the total mass density profile versus central baryon fraction for our Galaxy haloes at z = 2. When feedback is weak or absent, the inner slope is close to the isothermal value (β = −2); when feedback is strong, the slope is flatter and close to the DMONLY values (β ∼ −1.4; shown as the black, hatched region, Section 5.3). Comparing these results with the observations of Koopmans et al. (2006), shown as the lower maroon box in the figure, only the simulations with weak or no feedback produce similar values to the observations. The isothermal profile at z = 2 is also seen in Romano-Díaz et al. (2008), who found that for their simulations the influx of subhaloes at late times, z 1, acted to flatten the profile.
This conclusion is apparently at odds with what can be drawn from the observed stellar fractions. On the one hand, strong feedback is required to keep cooling under control, such that the observed stellar mass fractions in groups and clusters at low redshift are reproduced. On the other hand, Figure 4. We plot the mean density profile of a number of species, normalised by the average virial mass and radius, and scaled by r 2 to reduce the dynamic range. Here we can see the significant dynamical influence of the baryons in this particular object. This is a galaxy-sized object at z = 2 taken from ZC WFB, averaged over 20 relaxed haloes, from the simulation with comoving length 25 h −1 Mpc. The yellow triangles are the total matter profile; the black squares are the DM profile; the blue pentagons indicate the gas profile and the green circles denote the stellar profile. Error bars illustrate the 68 per cent confidence limits within each radial shell, estimated by bootstrap analysis. Vertical lines indicate the region where the inner slope is measured (0.025 < r/R vir < 0.05). The arrow represents R P03 , taken to be the resolution limit. The solid red curve is the best-fit NFW profile to the DM, over the outer region (r/R vir > 0.05). The DM mass was divided by (1 − f b (r < R vir )) such that the total mass of the DM component would equal M vir . With the halo mass known the NFW profile has only one free parameter. The NFW curve has had this factor removed for this plot. The legend contains the mean virial mass and radius of the haloes, M vir and R vir respectively, the bestfit NFW concentration, c vir , and inner profile slope, β. Baryons strongly influence both the total and the DM density profiles within 0.1 R vir . the feedback has to be weak (or absent), to generate the steep inner profiles observed in galaxies at low and intermediate redshift (out to z ∼ 1).
One possible reason for this dichotomy, is that the simulated galaxies are at higher redshift than the lensing observations 3 . However, as shown by Dehnen (2005), the density profile of a collisionless system does not steepen during mergers, so the lensing galaxies must have been isothermal at higher redshift, where it is more likely that significant gas condensation could have occurred. Secondly, we cannot rule out selection effects in the observations that may have a bias towards steeper density profiles. This could be important as isothermal profiles do exist in the ZC WFB AGN simulation, albeit significantly fewer in number than ZC WFB creates. A final possibility is that the simulations themselves do not yet accurately model the high-redshift growth of the brightest galaxies. One would therefore require a feedback prescription that was less effective at high redshift (allowing more gas to accumulate and steepen the central profile). 3 We use the highest resolution simulation at z = 2 due to the close match between typical Einstein radius of galaxies in Koopmans et al. (2006), ≈ 3 h −1 kpc, and our fitting range ≈ 2 − 4.5 h −1 kpc.
The feedback must then rapidly expel the gas before it can form stars, faster than the dynamical timescale of the halo, such that the DM retains an isothermal profile.

HALO DENSITY PROFILES
To first get a general idea of the effects of baryons on the DM halo, we show an example of the halo density profile and its components in Fig. 4, where we have plotted the total mass profile, as well as the individual gas, stellar and DM components. Here, the halo is an average over 20 relaxed Galaxy haloes from ZC WFB at z = 2, and is shown in dimensionless form: (r/Rvir) 2 Φ(r/Rvir), where Φ(r/Rvir) ≡ ρ(r/Rvir)R 3 vir /Mvir. As can be seen in the figure, the baryons are more centrally concentrated than the DM. The stellar component has a steeper profile than the DM at all radii, whereas the gas is steeper in the inner region (r/Rvir < 0.05) and slightly flatter in the outer region (r/Rvir > 0.1). The inner region of the DM profile has been significantly affected by the baryons. To highlight this change, we also show the best-fit Navarro, Frenk & White (Navarro et al. 1997; hereafter NFW) profile (solid, red curve), fit to data with 0.05 r/Rvir 1. The NFW profile, which is a good approximation to the equivalent DMONLY profiles, is of the form where δc is a characteristic density contrast and rs is a scale radius. We fit for only one parameter, rs, as we use the previously-measured virial radius and mass for each halo to define δc.
The simulated DM profiles are significantly steeper in the inner regions than the NFW profile. This is expected: the high central baryon fraction in this run (which, as we can see, is mainly in stars) has pulled the DM in towards the centre. Measurement of the inner slope of the DM profile (β ≈ −2) confirms that the inner profile is isothermal. It is therefore clear that the effects of baryonic physics (radiative cooling in this case) can have a profound impact on the inner DM density profile.

Galaxies and groups at high redshift
We present mean DM density profiles for the z = 2 data in Fig. 5. These consist of 32 haloes from the 25 h −1 Mpc simulation at Dwarf Galaxy and Galaxy mass scales (top panels) and an additional 14 Group scale objects (bottom panels) drawn from the larger 100 h −1 Mpc simulation. Note that the group scale objects do not satisfy our stringent criterion for the inner profile study: the black arrow, depicting the P03 resolution limit, is at r/Rvir ≈ 0.04.
We have grouped the simulations together according to the physics prescriptions that we are testing. The first set is shown in the left panels and corresponds to the introduction of metal-line cooling and weak stellar feedback. The second group, in the right panels, tests the types of feedback, namely stellar (both weak and strong) and AGN, all with the same cooling prescription. The DM density profiles have been divided by (1 − f univ b ) for comparison with the DMONLY profiles. Figure 5. In order to quantify the response of the DM halo to the varying physics prescriptions, we plot the mean DM halo density profile of haloes in fixed mass bins, M vir = 5 × 10 11 − 5 × 10 12 h −1 M ⊙ , the mean mass of which corresponds to Galaxies (top panel) and M vir = 10 13 − 5 × 10 13 h −1 M ⊙ (the latter limit is to include the largest object) which are Groups (bottom) at z = 2 (drawn from the 25 h −1 Mpc simulations). Units are normalised by the virial mass and radius from the equivalent haloes in the DMONLY simulation. The DM density profiles for the models including baryons have been divided by (1 − f univ b ) to facilitate comparison with DMONLY. The left panels show results for the runs with weak or no feedback and the right panels for different feedback schemes. The vertical lines denote r/R vir = 0.025 and 0.05, the region where the inner profile slope is estimated. The vertical arrows denote the P03 resolution limit; in the case of the Groups, the inner profile slope cannot be reliably measured at z = 2 (the other haloes have a resolution limit within the innermost line). Error bars are bootstrap estimates of the 68 per cent confidence limits about the mean density within each bin. The legend contains the mean virial mass of the haloes, M vir , the best-fit NFW concentration, c vir , and inner profile slope, β. As expected, the largest differences occur in the inner regions, where the density is highest. In the left panels, showing results for runs in which feedback effects are weak or absent, there is a significant steepening of the density profile towards higher densities on smaller scales across all mass ranges. The curves then more-or-less converge with DMONLY at radii larger than r/Rvir = 0.2. Given that PrimC WFB and PrimC NFB predict nearly identical DM profiles, but differ strongly from DMONLY, it is clear that cooling plays a crucial role in determining the magnitude of the back-reaction on the DM. Indeed, comparing PrimC WFB with ZC WFB, we see that including metal-line cooling results in a steeper profile for Galaxies. This difference is reduced in more massive systems for which the virial temperatures exceed the regime in which metalline cooling is efficient (e.g. Wiersma et al. 2009).
In the right panels of Fig. 5 we see that the stronger feedback schemes (ZC SFB and ZC WFB AGN) have DM density profiles that are closer to those from DMONLY. These runs have smaller baryon fractions, thus the overall effect of the cooling has been reduced.

Groups and clusters at the present day
For Groups and Clusters at z = 0 we repeat our previous investigation in Fig. 6. Again, the DM profiles in the simulations with baryons diverge from the DMONLY results in the inner regions, although the effect is not as dramatic as it was for the less massive systems shown in Fig. 5. This is because the typical cooling times are longer in these systems. Furthermore, the difference between the runs with weak and strong stellar feedback (ZC WFB and ZC SFB; right panels) is smaller, because the strong stellar feedback is less effective in the more massive Groups and Clusters than it was in the lower-mass systems at high redshift. The ZC WFB AGN model (with the lowest central baryon fractions) produces almost identical profiles to DMONLY in Clusters, but a slightly lower profile in Groups.
Another interesting effect apparent in Fig. 6 for the Clusters at z = 0 concerns the difference between the DM profiles of the no and weak feedback schemes (left panels). The model showing the highest central DM density is PrimC WFB, while PrimC NFB and ZC WFB are statistically indistinguishable. At first glance this is unexpected, as the latter two models have higher central baryon fractions (see Fig. 1, bottom panel). A similar result was found by Pedrosa et al. (2009) for their galaxy simulations at z = 0. We summarise their explanation of the effect within the context of our simulations. When weak stellar feedback is included (going from PrimC NFB to PrimC WFB), a certain amount of the gas will be expelled from satellite galaxies but not from the main (group or cluster) halo itself. As a result, the satellite haloes will be less bound and suffer more tidal disruption as their orbit decays due to dynamical friction. The result is that less mass is transferred to the centre of the halo. When metal enrichment is added to the simulation (going from PrimC WFB to ZC WFB), the cooling time of the satellite gas becomes shorter (due to metal-line emission) and as a result, the satellite is able to hold on to more of its gas. This reduces the effect the feedback has on the evolution of the satellite halo itself (and thus more mass can be transferred to the centre). The most important con-sequence of these effects is where the angular momentum of the satellites gets transferred. In the case of PrimC WFB, less of the angular momentum is transferred to the inner region than in PrimC NFB, for example. As a result, the inner profile is denser in the former case, even though the overall baryonic mass is smaller in the centre of the halo.

Inner profiles
It is clear that the main driver of the change in the (inner) DM profile is the condensation of gas to smaller scales than the DM, thereby increasing both the local baryon fraction and the DM concentration in the now-deepened potential well. If we assume that DM is pulled inwards by the condensing baryons, then we should expect some relation between the DM density profile and the baryon fraction within the inner region. We test this in Fig. 7 by plotting the median inner power law slope of the DM profile, βDM, as a function of the median f b (r/Rvir 0.05) value. The quartile spread in β values for the DMONLY run is also illustrated, as the hatched region. We remind the reader that in the inner profile study we utilised a well-resolved subsample of the full halo catalogue, such that RP03 0.025Rvir.
For the z = 0 Group and Cluster haloes (bottom panel), f b varies by more than a factor of 3 in the inner region of the haloes and yet the resultant inner slope stays within 20 per cent of the DMONLY value. There is a tentative steepening of the inner profile with increasing baryon fraction, but the steepest profile found in the simulation (PrimC WFB) does not have the largest central baryon fraction, as already discussed in Section 5.2.
Our z = 2 Galaxy sample (top panel in Fig. 7) does demonstrate a clear and significant trend of steepening inner DM density profile with increasing central baryon fraction. The lowest baryon fractions are found in the ZC WFB AGN simulation for which the DM haloes are indistinguishable from those in DMONLY. In models PrimC NFB and ZC WFB the baryon-dominated central regions generate nearly isothermal (β = −2) inner DM density profiles.

HALO CONCENTRATIONS
For an NFW halo of a given mass, the DM density profile is specified entirely by one parameter, the concentration. In this section, we measure and compare NFW concentrations for the DM profiles of haloes in our simulations (Section 6.1). We then consider the total mass profile of the halo by utilising simple, non-parametric measures of the concentration, based on ratios of density contrasts at different scales (Section 6.2) and the velocity profile of the halo, in particular the maximum velocity, in Section 6.3. All haloes with a P03 convergence radius, RP03 < 0.05Rvir, are considered in this work on halo concentrations. This sets an effective mass limit of Mvir > 8 × 10 10 h −1 M⊙ at z = 2 and Mvir > 5 × 10 12 h −1 M⊙ at z = 0. The sample of haloes includes over 500 objects at z = 0 and ∼ 20 (∼ 150) systems at z = 2 in the simulation volumes with comoving length 100 (25) h −1 Mpc.

DM profile: NFW concentrations
A well established result from N -body simulations (Bullock et al. 2001) is that the NFW concentration of the DM halo is anti-correlated with its mass (for the latest results see Neto et al. 2007;Duffy et al. 2008;Gao et al. 2008 andMacciò et al. 2008) and takes on a power-law form where Bvir is close to −0.1 when fit to data over nearly five orders of magnitude in mass. However, it is not clear how much this trend, which is primarily driven by the formation time of the halo, is modified by the presence of baryons. Observations of X-ray luminous groups and clusters (Buote et al. 2007;Schmidt & Allen 2007) suggest a steeper dependence of concentration on mass than the DM only simulations predict, as was pointed out by Duffy et al. (2008). This is primarily due to observed groups having ≈ 30 per cent higher concentrations than the simulated objects (the concentrations of clusters were in good agreement with the simulations if a subsample of dynamically-relaxed haloes was used in the comparison). It is therefore important to check whether the inclusion of baryons can bring theory and observations into agreement. We will present concentrations relative to the equivalent values for the DMONLY model. Note that our simulations assume the WMAP3 cosmology, which has an 8 per cent lower value of σ8 than the WMAP5 value (0.74 versus 0.796) assumed by Duffy et al. (2008) and leads to somewhat smaller concentrations at fixed mass 4 . As described in Section 2.2 we fit density profiles over the range 0.05 r/Rvir 1 We first assess the goodness-of-fit of the NFW function when baryons are included. To do this, we compute the usual quantity where N bins is the number of bins in our profile and ρsim and ρNFW are the densities from the simulation and the bestfit NFW profile respectively. For DMONLY we find that σ fit ≈ 0.02. For the runs with baryons the goodness-of-fit is similar (typically within 10 per cent) for Group and Cluster haloes. For smaller objects the difference is more pronounced for the simulations with high central baryon fractions, for which σ fit increases by around a factor 2. Shown in Fig. 8 are the NFW concentrations of the DM haloes in the runs with baryons, relative to the corresponding values from DMONLY, as a function of halo virial mass. For the Group and Cluster haloes at z = 0 (bottom row), the only simulation to show substantial (> 10 per cent) deviations from the DMONLY values, is the simulation that includes AGN feedback, ZC WFB AGN. For this model the NFW concentrations of Mvir ∼ 10 13 h −1 M⊙ haloes are about 15 per cent lower. For the run with strong stellar feedback, ZC SFB, the decrease is about 10 per cent. In both cases the expulsion of gas has caused the DM to expand relative to the DMONLY case (e.g. Hills 1980).
The effect of baryons on the NFW concentrations of z = 2 Dwarf Galaxy, Galaxy and Group haloes is shown in the top row of Fig. 8. As before, the differences are more dramatic for these lower mass, higher redshift objects. In the runs without strong feedback the concentration increases, as expected. The increase is typically 10-20 per cent for Groups, but can be as large as 50 per cent (when supernova feedback is absent) for Dwarf Galaxies. This dramatic increase is similar in magnitude to that found by Romano-Díaz et al. (2009) for the concentration of a Mvir ≈ 10 12 h −1 M⊙ halo. As was the case at z = 0, in runs with effective feedback the presence of baryons makes little difference to the concentration of the DM; the maximum effect being a decrease of ∼ 15 per cent for Galaxies in ZC WFB AGN.
We checked how the results change when only relaxed haloes are selected (as defined in Duffy et al. 2008). As was found in previous work (e.g. Duffy et al. 2008), the average concentration of the haloes increases, but the power-law relation between concentration and mass remains the same within the errors.
We have also checked explicitly that the concentration of a halo increases with its central baryon fraction. This is indeed the case for all simulations and mass scales, except for the Group and Cluster haloes at z = 0, in the PrimC WFB model. As discussed in the previous section, this run has haloes on these mass scales with anomalously high central DM densities (as compared with the run with primordial cooling and no feedback, PrimC NFB).
In Duffy et al. (2008) we demonstrated that the inferred NFW concentrations of groups observed in X-rays at z = 0 are 30 per cent more concentrated than predicted from DM only simulations. It was suggested that the inclusion of baryons would alleviate this discrepancy through the contraction of the DM halo. As is clear from Fig. 8, however, the largest increase of any physics scheme is still less than 10 per cent. Moreover, for those schemes which reproduce the ob- Figure 9. Here we plot the median ratio of the spherical overdensity radii R 500 and R 2500 as a function of halo mass, M 500 in this case, at z = 2 (top) and 0 (bottom). This ratio is a useful non-parametric measure of the concentration of a halo. The vertical error bars are 68% confidence limits about the medians, estimated by bootstrapping the samples within each bin, and the horizontal values represent the mass range of the bin. Note that higher values of the ratio indicate lower concentrations. All simulations show a positive trend with mass; this is indicative of the decreasing amount of baryons able to cool and condense as the virial temperature of the halo increases. The offset in the normalisation between simulation schemes is larger than the scatter, allowing the use of this ratio as a nonparametric concentration proxy.
served z = 0 stellar fractions (shown in Fig. 2), the strength of feedback is such that we actually reduce the concentration of the Groups relative to the DMONLY simulation. As an additional issue, the trend of the concentration ratio with mass is positive in the strong feedback schemes which may further increase the disagreement with observations. Clearly, the inclusion of baryons does not resolve the problem. It is thus important to check whether observational biases or selection effects can account for the mismatch between theory and observation.

Total mass profile: Non-parametric concentrations
Although the NFW profile is a reasonable approximation to the DM profile at large radii (r/Rvir > 0.05, as shown in Fig. 4), the measurement of its concentration parameter requires the shape of the profile to be constrained over a range of scales, with accurate removal of the baryonic mass profile. A simpler, and thus more easily achievable method from an observational point of view, is to consider the entire halo, DM and baryon components together, and measure the mass/radius ratio of a halo at two different spatial scales. In Fig. 9 we plot the ratio of two radii that are commonly used by observers (e.g. in observations of X-ray groups and clusters), R500/R2500 5 , as a function of halo mass. All runs demonstrate a positive, albeit weak, dependence on mass with significant run-to-run variations in the normalisation. (Note that runs with higher central baryon fractions will typically have lower R500/R2500 ratios, because the value of R2500 grows as a result of the increased central mass.) The deviation from DMONLY is at the sub-25 per cent level and much smaller if the feedback is strong. The differences between the models are qualitatively similar to those for the NFW concentrations discussed in the previous section. Note, however, that contrary to the NFW concentrations, the non-parametric total matter concentrations are never significantly reduced relative to the DMONLY case.

Concentration measures: circular velocity
A key, and relatively easily observed measure for the structure of a halo is the circular velocity, vc(r) = (GM (< r)/r) 1/2 , in particular its maximal value, vmax. Like R500/R2500, it can be more robustly determined than the NFW concentration. The creation of realistic velocity profiles for galaxies has, however, been a long standing issue in simulations within the CDM paradigm. For N -body only simulations, the maximum velocity is well approximated by the analytic solution to the NFW profile vmax = vc(r ≈ 2.17rs) (Cole & Lacey 1996;Navarro et al. 1996), with the final result that vmax ≈ Vvir = (GMvir/Rvir) 1/2 . Observations find that the maximum velocity in the disk is similar to the virial velocity (e.g. Dutton et al. 2005). Typically, however, when simulations include baryons, either explicitly or through adiabatic contraction models, the haloes will have a maximum velocity that is a factor of two higher than the halo virial velocity (e.g. Navarro & Steinmetz 2000;Dutton et al. 2007;Pedrosa et al. 2010). This increase in the velocity ratio is a consequence of the contraction of the halo in the presence of a significant central baryon component, with a large velocity ratio indicative of a high NFW concentration.
Recently, Pedrosa et al. (2010) investigated the circular velocity profile in a resimulation of a single high resolution galaxy for various physics implementations. They found that there was a positive correlation between the ratio vmax/V200 and the 'sharpness' of the DM density profile 6 , indicative of the contraction of the halo in the presence of baryons. Furthermore, effective feedback was necessary to obtain a realistic, i.e. low, velocity ratio in their galaxy. In addition to galaxies we have investigated a larger mass range, groups and clusters. We also lend statistical weight to conclusions concerning the importance of the various physics implementations by examining a large sample of haloes. A more extensive investigation of the rotation properties of the OWLS 5 Typical values for these radii in terms of R vir are ≈ 0.5 (0.6) and 0.2 (0.3) for R 500 and R 2500 at z = 0 (2) respectively. The values will change by 10 per cent dependent on the simulation physics due to the baryonic back-reaction. galaxies at z = 2, in addition to vmax, is being performed in a separate study (Sales et al., in prep).
In Fig. 10 we show the maximum circular velocity, in units of the virial velocity, as a function of halo mass at z = 2 (0) in the top (bottom) panel for various implementations of the cooling and feedback prescriptions. The most striking result is the good agreement between the strong feedback schemes (and DMONLY) with vmax ≈ Vvir at all masses and redshifts, and the significant offset for the other physics implementations. This divergence is reduced with increasing mass at all redshifts due to the strong anti-correlation of the velocity ratio with halo mass in the weak/no feedback runs. This is indicative of the reduction in gas cooling efficiency as the halo mass, and hence the virial temperature, increases.
At z = 2 (top panel) there is a dramatic divergence in the velocity profile below 10 13 h −1 M⊙ between the runs with primordial cooling, due to the supernova feedback. When metal-line cooling is included, ZC WFB, the divergence with the no feedback model is reduced; the overall effect of the enhanced cooling is therefore to counteract the gas removal efforts of the supernovae.
At z = 0 (bottom panel of Fig. 10) we once again find the counterintuitive result that the model with weak feedback, PrimC WFB, is more concentrated in clusters (but not in groups) than the model without feedback PrimC NFB. We also see the familiar effect of the metal-line cooling overcoming the weak supernova feedback when comparing PrimC NFB and ZC WFB models.
It is clear that strong feedback is necessary if one wishes to limit the effect of the baryons in increasing the maximum circular velocity of the halo above the virial velocity.

ADIABATIC CONTRACTION
Having demonstrated that baryons can significantly influence the DM density profile, we will now assess the degree to which this modification can be modelled as a secular adiabatic contraction of the DM halo. We will test the models of Blumenthal et al. (1986), henceforth B86, and G04, and will make use of the publicly-available code contra (G04).
The B86 model for adiabatic contraction assumes that the DM halo is spherically-symmetric and that the particles are on circular orbits so that we can think of the halo as a series of shells that can contract but do not cross. Assuming that the baryons initially trace the DM halo density profile and then fall slowly (i.e., such that the mass internal to radius r changes slowly compared to the orbital period at r) towards the centre as they cool, we can compute the response of the dark matter shells. In that case the dark matter particles conserve their angular momentum and hence rvc ∝ [rM (r)] 1/2 is conserved, where M (r) is the total mass internal to r. We thus have Figure 10. We plot the maximum circular velocity of the halo, relative to the circular velocity at the virial radius, R vir , as a function of virial mass at z = 2 (0) in the top (bottom) panel. This measure is sensitive to the degree of contraction of the halo in the presence of baryons. Each point corresponds to the median value for haloes within a given mass bin. Vertical error bars are bootstrap estimates of the 68% confidence interval for each bin and the horizontal bars are the bin widths. A number of familiar halo responses is easily seen in this figure, with galaxies more sensitive to the presence of the baryons than higher mass systems. The very high maximum velocity in PrimC NFB indicates the importance of feedback in forming observed galaxies, whose maximum circular velocity in the disk is similar to the virial velocity (Dutton et al. 2005). Furthermore, the sharp divergence between PrimC NFB and PrimC WFB at a critical mass of 10 13 h −1 M ⊙ indicates the mass range at which the weak feedback model becomes ineffective. Finally, the metal-line cooling has the overall effect of counteracting the supernova feedback in Galaxies and Groups at z = 2, seen by the agreement between PrimC NFB and ZC WFB. The anti-correlation seen at z = 0, for the simulations PrimC NFB and ZC WFB, is indicative of the rising halo virial temperature restricting the efficiency of gas cooling. The strong feedback schemes are in close agreement to the DMONLY simulation predictions at all masses and redshifts, and are therefore necessary to recover observed values.
where M dm,i (ri) is the initial DM profile, M b,f (r f ) is the final baryon profile, and r f is the final radius of the DM shell that was initially (i.e. before the baryons contracted) at ri. Note that we made use of the equality M dm,f (r f ) = M dm,i (ri), which holds because DM shells do not cross. Assuming that we know both M b,f (r f ) and M dm,i (ri), we can solve for ri (and thus M dm,i (ri) = M dm (r f )) as a function of r f . In the following we will make the standard assumption that M dm,i (ri) is given by the density profile (reduced by (1 − f univ b )) of the corresponding DMONLY halo. G04 found that this method did not provide a satisfactory description of their numerical simulations. Because the DM particles are typically not on circular orbits, they suggested replacing M (r)r by M (r)r, wherer is the orbitaveraged radius, which they found can be approximated as where A = 0.85 and w = 0.8 (c.f. the B86 model which has A = w = 1). Gustafsson et al. (2006) extended this work and showed that the values of A and w that best describe the difference between simulations with and without baryons depend on halo mass and the baryonic physics implementation. Furthermore, they showed that the best fitting values of A and w generally differ substantially from the values that provide a good fit to the actual mean radiusr of the DM particle orbits. This suggests that while the introduction of the two free parameters A and w enables better fits, the underlying model does not capture the relevant physics.
Here, we extend the work by Gustafsson et al. (2006) to a larger range of mass scales and gas physics models.

Best-fit G04 parameters
We first compute the best-fitting parameters A and w for a subset of our simulated haloes with baryons to see if our simulations prefer a specific combination and, if so, how this compares with those found in previous work. We focus on 3 models: the runs with no feedback (PrimC NFB), weak stellar feedback (ZC WFB) and strong stellar feedback (ZC SFB). For each simulation we matched the haloes to the corresponding objects in DMONLY, producing an average density profile at each redshift. At z = 2 we consider haloes with Mvir = [5 − 50] × 10 11 h −1 M⊙, while at z = 0 we study the range Mvir = [3 − 60] × 10 13 h −1 M⊙. In Fig. 11 we show the distribution of σ fit values in the A − w plane, for each of the three simulation models at z = 2 (top panels) and z = 0 (bottom panels). We calculate σ fit using equation (7), replacing ρNFW with the appropriate adiabatically contracted density profile, over the range 0.025 r/Rvir 1.
The best-fit values of A and w depend strongly on both halo mass and on the baryonic physics. There is significant degeneracy between the parameters A and w, higher values of A can be compensated by higher values of w. The parameters suggested by G04 7 and the original model by B86 work best for simulation ZC WFB, but even here they differ substantially from our best-fit values. These findings are in good agreement with Gustafsson et al. (2006).

Predicted DM density profiles
The DM density profiles predicted by the adiabatic contraction models are shown in Fig. 12. When the parameter values suggested by G04 are used, the G04 model predicts similar profiles to B86, which do not describe the contracted DM profiles well. The models typically underestimate the DM density for r 10 −1 Rvir and more so for the simulations with stronger feedback. If, on the other hand, we use the best-fit values of A and w for each simulation and halo sample then the predictions of the G04 method agree much better with the simulated profiles, but even in that case one would obtain a closer match to the actual density profile by neglecting adiabatic contraction for r > 0.1Rvir.
It is perhaps not too surprising that the models for adiabatic contraction do not describe the simulations well. The assumption that the baryons initially trace the DM halo profiles is clearly violated in hierarchical models. Haloes are built by mergers of smaller progenitors, and cooling and feedback have already redistributed the baryons in these objects. Rather than contracting slowly as the gas cools, a large fraction of the baryons simply fall in cold (Kay et al. 2000;Kereš et al. 2005). Moreover, in the stronger feedback models a substantial fraction of the baryons is ejected. Tissera et al. (2009) recently demonstrated that the contraction is manifestly not adiabatic. They find that the pseudo-phase space density relation is strongly modified when baryons were added to high resolution DM only simulations.
Models for adiabatic contraction are required when full hydrodynamic simulations are not available. Unfortunately, it is not possible to predict what values of A and w to use for the G04 model without a much better understanding of the baryonic physics. Moreover, even if the physics were known, we would need to simulate many haloes because the bestfit values of the parameters depend on both halo mass and redshift.
Our results suggest that it is better to ignore adiabatic contraction for r > 0.1 Rvir, but that the use of adiabatic contraction models such as those by B86 and G04 does typically represent an improvement at smaller radii, provided the feedback is moderate or weak.

CONCLUSIONS
Our main aim in this work was to investigate the response of the DM halo to the presence of baryons at a variety of masses and redshifts. We utilised a series of high-resolution simulations within a cosmological volume, with a number of different prescriptions for the sub-grid physics, to probe the effect of baryons on the DM distribution of haloes. Our results centred on galaxy-scale haloes at z = 2 and groups and clusters at z = 0. We were particularly interested to discern the effect of the baryons when going from the situation in which radiative cooling dominates (leading to a high central concentration of baryons in the form of stars and cold gas) to one where feedback dominates (reducing the central galaxy mass and expelling gas from the halo). As we showed Figure 12. We test adiabatic contraction models by comparing the predicted mean DM density profiles, in average virial units of the matched DMONLY haloes, for the simulations PrimC NFB (left), ZC WFB (middle) and ZC SFB (right) respectively. Also shown are the mean DMONLY profiles (which have been reduced by (1 − f univ b ) for comparison). Overlaid are the predictions from the adiabatic contraction models of B86 (triple-dot-dashed), G04 using their default parameter values (scaled to our definition of the virial radius) (dashed), and G04 with our best-fit parameter values (solid) as determined separately for each model, mass range, and redshift in Fig. 11 along with the goodness-of-fit measure σ fit . Additionally, in the top legend we give the mean virial mass of the haloes, M vir , the best-fit NFW concentration, c vir , and inner profile slope, β. The top (bottom) row shows results for haloes at z = 2 (z = 0). No one parameter combination can reproduce the range of DM haloes, indicating that the slowly cooling gas picture, on which adiabatic contraction models are based, is not sufficient to model the behaviour of the DM in a live simulation. In fact, ignoring adiabatic contraction altogether gives the best results for r > 0.1 R vir .
in Fig. 1, when comparing central and global baryon fractions, our simulations with strong feedback (ZC WFB AGN and ZC SFB) reduce the baryon fraction f b , in comparison with PrimC NFB (no feedback), by factors of 2-3 for Galaxy haloes at z = 2. In Group and Clusters the AGN can remove nearly half of the baryons from the inner region of the halo, at z = 0.
By comparing with observed stellar fractions in low redshift Groups and Clusters of galaxies (Fig. 2) we found that the simulations with a high baryon fraction (Fig. 1) also predict stellar fractions significantly larger than observed. The simulation with inefficient gas cooling and stellar feedback, PrimC WFB, and the strong feedback models, ZC WFB AGN and ZC SFB, are broadly in agreement with the observed stellar fractions in z = 0 objects of mass Mvir ≈ 10 14 h −1 M⊙. However, observed maximum star formation efficiencies of order 10% -20% are only reproduced by the inclusion of AGN feedback.
However, these same strong feedback simulations are in disagreement with the constraints inferred from combined gravitational lensing and stellar dynamics analyses of the inner total mass density profile of massive, early-type galax-ies (Fig. 3). In this case the efficient feedback prevents the steepening of the density profile that is necessary to reproduce the observed isothermal profiles. Instead, the simulations with high baryon fractions are in closer agreement. A more detailed comparison between the observations and simulations is warranted, especially with regards to the biasing of lensing observations to steeper density profiles.
An enhanced baryon fraction in the inner halo is expected to contract the DM distribution, as a response to the deeper potential well of the system. This effect was clearly seen, from Dwarf Galaxy to Cluster scales and at both low (z = 0) and high (z = 2) redshifts, and is summarised in Fig. 7, where haloes with larger central baryon fractions develop steeper central profiles (especially the Galaxy haloes at high redshift). To quantify the contraction effect, we fit NFW profiles and compared them to a simulation with no baryons. Variations in the concentration are typically around 20 per cent or less, although the concentration of high-redshift dwarf galaxy haloes can be as much as 50 per cent higher when feedback effects are ignored. Strong feedback produces a mild decrease in the concentration of a halo through the removal of a significant amount of gas from the halo, in a similar manner to what was found by Koyama et al. (2008).
We also investigated non-parametric measures of concentration. We found that the ratio R500/R2500 and the maximum circular velocity are both useful indicators of the degree of contraction of the total mass profile. Efficient feedback is required to redistribute the mass in the halo such that the maximum circular velocity is similar to the virial velocity, as is found observationally for disk galaxies.
An interesting counterexample to the rule of increasing baryonic condensation leading to more concentrated DM haloes was witnessed. In our low redshift groups and clusters, the dark matter is denser in the simulation with weak stellar feedback (PrimC WFB) than in the simulation with no feedback (PrimC NFB), even though the central baryon fractions are lower in the latter case. A similar result was found by Pedrosa et al. (2009), who analysed simulations of galaxies with and without feedback. They argued that the contraction of the DM halo is slowed by the infall of satellites (due to dynamical friction). If the satellites are less-bound, they may be tidally disrupted before they sink to the centre and their effect on halting the contraction is thus reduced.
We compared our results to the adiabatic contraction models of Blumenthal et al. (1986) and the revised model of Gnedin et al. (2004). We found that adiabatic contraction only improves the fits in the inner parts of the halo. Outside 0.1 Rvir one would do better not to use any prescription at all. Within 0.1 Rvir the former model is unable to reproduce the back-reaction on the DM. The latter model can do a reasonably good job, but only if its two parameters are allowed to vary with the baryonic physics, halo mass, and redshift, in agreement with the findings of Gustafsson et al. (2006). This is of particular importance when feedback effects are strong, as required to reproduce the cooled baryon fractions in low redshift groups and clusters. This therefore removes the predictive power of the model and we caution its use, particularly if a detailed prediction for the structure of a DM halo is required.
If one wishes to match the stellar fractions at low redshift, then models with strong feedback are required to suppress star formation. The total mass profiles in such simulations are very similar to the DMONLY case, with vmax ≈ Vvir, as observed. Intriguingly, the NFW concentration in these systems is actually reduced relative to the same halo in the DMONLY run. The finding reported in Duffy et al. (2008) that observed groups appear to be significantly more concentrated than simulated haloes, is therefore still unexplained. The hypothesis that the inclusion of baryons would resolve the discrepancy between theory and observation has been shown to be wrong. Worse, the disagreement actually grows larger if one utilises strong feedback physics schemes that can reproduce the observed stellar fractions in these systems.

393, 99
Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574 Zeldovich Y. B., Klypin A. A., Khlopov M. Y., Chechetkin V. M., 1980, Sov. J. Nucl. Phys., 31, 664 Zibetti S., White S. D. M., Schneider D. P., Brinkmann J., 2005 APPENDIX A: RESOLUTION TESTS To assess the effects of resolution on our results, we compare the 512 3 DMONLY and ZC WFB simulations with lower resolution versions (256 3 and 128 3 , with softening lengths increased by factors of 2 and 4 respectively). We initially select haloes with 10 3 DM particles from the lowest resolution simulation. Matching haloes in the higher resolution simulations are found by first identifying high-resolution haloes that lie within the virial radius, R 128 pot , of the low resolution halo. We consider a halo to match if the candidate object fulfils the following criteria This final selection removes ∼ 10 per cent of the haloes from the sample, these objects are in the process of merging and hence would create artifacts in the averaged density profiles. We therefore perform resolution tests on 57 (27) matched haloes from ZC WFB at z = 0 (2) and 51 (27) haloes from DMONLY at z = 0 (2). Note that this matching is different to the scheme employed in the rest of this work which was based on linking identical DM particles between different physics simulations. Fig. A1 compares average total mass density profiles for a group and a cluster-scale halo at z = 0 in the three runs. In the DMONLY case, resolution effects cause the density to be underestimated, whereas the density is overestimated for ZC WFB. At z = 2, the results from (lower mass) dwarf and galaxy-scale haloes are similar, although the density is also underestimated in the under-resolved regions in ZC WFB, (Fig. A2). Power et al. (2003) defined a convergence radius, such that the mean two-body relaxation time of the particles within that radius is of the order a Hubble time. The equivalent radius is shown for each of our simulations in the two figures, as a solid vertical arrow. Although this radius was originally applied to collisionless N -body simulations (like DMONLY), we can see that it also provides a satisfactory indication of where the density profile becomes numerically resolved if baryons are included. This is especially visible in the lower sub-panels, where the fractional deviations of the two lower-resolution profiles from our standard-resolution (512 3 ) results are shown. Figure A1. Comparison of mean DM density profiles (and their residuals) for Group and Cluster haloes (top and bottom panels respectively), drawn from the 100 h −1 Mpc box at z = 0, for runs with different resolutions. Results from DMONLY are shown in the left column, while results from ZC WFB are shown in the right column. The colours indicate the resolution, from 128 3 (red), to 256 3 (blue) to 512 3 (black). Most softening lengths are sufficiently small to fall outside of the plotted area but, where visible, softening scales are indicated with vertical lines with colour and line-style indicating simulation resolutions, described in the legend. The arrows indicate the P03 convergence radius for each simulation, according to colour. The vertical lines denote 2.5 and 5 per cent of the virial radius, corresponding to the scale within which the inner density profile slope is measured. The legend contains the mean virial mass of the haloes, M vir . Figure A2. As in Fig. A1, but for Dwarf Galaxy (top panels) and Galaxy (bottom panels) haloes at z = 2, drawn from the 25 h −1 Mpc box.