# Density profile of dark matter haloes and galaxies in the horizon-agn simulation: the impact of AGN feedback

S. Peirani, Y. Dubois, M. Volonteri, J. Devriendt, K. Bundy, J. Silk, C.

Pichon, S. Kaviraj, R. Gavazzi, M. Habouzit

► To cite this version:

S. Peirani, Y. Dubois, M. Volonteri, J. Devriendt, K. Bundy, et al., Density profile of dark matter haloes and galaxies in the horizon-agn simulation: the impact of AGN feedback. Monthly Notices of the Royal Astronomical Society, 2017, 472 (2), pp.2153-2169. 10.1093/mnras/stx2099. hal-01645203

HAL Id: hal-01645203

https://hal.science/hal-01645203v1

Submitted on 29 Aug 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Density profile of dark matter haloes and galaxies in the HORIZON–AGN simulation: the impact of AGN feedback

Sébastien Peirani,<sup>1,2,3,4</sup> Yohan Dubois,<sup>2*</sup> Marta Volonteri,<sup>2</sup> Julien Devriendt,<sup>5,6</sup>

Kevin Bundy,<sup>3</sup> Joe Silk,<sup>2,5,7</sup> Christophe Pichon,<sup>2,8</sup> Sugata Kaviraj,<sup>9</sup>

Raphaël Gavazzi<sup>2</sup> and Mélanie Habouzit<sup>2,10</sup>

Accepted 2017 August 11. Received 2017 August 11; in original form 2016 November 22

## ABSTRACT

Using a suite of three large cosmological hydrodynamical simulations, HORIZON–AGN, HORIZON–NOAGN (no AGN feedback) and HORIZON–DM (no baryons), we investigate how a typical sub-grid model for AGN feedback affects the evolution of the inner density profiles of massive dark matter haloes and galaxies. Based on direct object-to-object comparisons, we find that the integrated inner mass and density slope differences between objects formed in these three simulations (hereafter, $H_{AGN}$, $H_{noAGN}$ and $H_{DM}$) significantly evolve with time. More specifically, at high redshift ($z \sim 5$), the mean central density profiles of $H_{AGN}$ and $H_{noAGN}$ dark matter haloes tend to be much steeper than their $H_{DM}$ counterparts owing to the rapidly growing baryonic component and ensuing adiabatic contraction. By $z \sim 1.5$, these mean halo density profiles in $H_{AGN}$ have flattened, pummelled by powerful AGN activity (‘quasar mode’): the integrated inner mass difference gaps with $H_{noAGN}$ haloes have widened, and those with $H_{DM}$ haloes have narrowed. Fast forward 9.5 billion years, down to $z = 0$, and the trend reverses: $H_{AGN}$ halo mean density profiles drift back to a more cusped shape as AGN feedback efficiency dwindles (‘radio mode’), and the gaps in integrated central mass difference with $H_{noAGN}$ and $H_{DM}$ close and broaden, respectively. On the galaxy side, the story differs noticeably. Averaged stellar profile central densities and inner slopes are monotonically reduced by AGN activity as a function of cosmic time, resulting in better agreement with local observations.

**Key words:** methods: numerical – galaxies: evolution – galaxies: haloes – galaxies: jets – dark matter.

## 1 INTRODUCTION

The inner structure of dark matter (DM) haloes represents one of the most important constraints on cosmological and galaxy formation models. Within the cold dark matter (CDM) paradigm, earlier works ignoring the effects of gas dissipation and star formation processes, have suggested that DM haloes have spherically averaged density profiles that can be well fitted by simple analytical function such as

the NFW profile (Navarro, Frenk & White 1996a; Navarro, Frenk & White 1997) in which the inner slope tends to −1 or could even be as steep as −1.5 (e.g. Moore et al. 1998; Jing & Suto 2000). Later numerical works favour the Einasto profile with a finite central density though this profile is still cuspy and close to the NFW one (Merritt et al. 2006; Stadel et al. 2009; Navarro et al. 2010).

The prediction of a peaked inner density profile seems to be seriously challenged by observations. In particular, dwarfs and low surface-brightness galaxies suggest a much shallower profile with a relatively flat slope ($\alpha \ge -0.5$ with $\rho(r) \sim r^{\alpha}$, Palunas & Williams 2000; Salucci & Burkert 2000; de Blok et al. 2001;

<footnote n="1"><sup>1</sup>Observatoire de la Côte d'Azur, Université Côte d'Azur, CNRS, Laboratoire Lagrange, F-06300 Nice, France</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris (UMR 7095: CNRS and UPMC), 98 bis Bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="4"><sup>4</sup>Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan</footnote>

<footnote n="5"><sup>5</sup>Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="6"><sup>6</sup>Observatoire de Lyon, UMR 5574, 9 avenue Charles André, Saint Genis Laval F-69561, France</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, The Johns Hopkins University, Homewood Campus, Baltimore, MD 21218, USA</footnote>

<footnote n="8"><sup>8</sup>Korea Institute of Advanced Studies (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote n="9"><sup>9</sup>Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield, Herts AL10 9AB, UK</footnote>

<footnote n="10"><sup>10</sup>Center for Computational Astrophysics CCA, 160 5th Ave, New York 10010, USA</footnote>

<footnote>* E-mail: sebastien.peirani@oca.eu (SP); dubois@iap.fr (YD)</footnote>
Swaters et al. 2003; Gentile et al. 2004; Spekkens, Giovanelli & Haynes 2005; Goerdt et al. 2006; Walter et al. 2008; de Blok et al. 2008; Oh et al. 2011; Walker & Peñarrubia 2011). On the other extreme of the halo mass scale, galaxy clusters tend to have a central cusp, as evidenced by studies combining strong and weak lensing (e.g. Limousin et al. 2007; Leonard et al. 2007; Umetsu, Takada & Broadhurst 2007), but shallower slopes than the NFW one can also be found within the inner ≈5 kpc (Sand et al. 2004, 2008; Newman et al. 2009, 2011, 2013; Richtler et al. 2011). Note that a recent study found a cusp slope over 5–350 kpc of $\alpha = -1.62$, again from combined weak- and strong-lensing analysis of the complete Cluster Lensing And Supernova Survey with Hubble cluster sample (Zitrin et al. 2015). At the intermediate halo mass scales, it is worth mentioning that Oldham & Auger (2016) have recently demonstrated the presence of a core at the centre of the dark halo of the massive elliptical galaxy M87, by combining stellar kinematics in the central regions with the dynamics of 612 globular clusters.

This so-called cusp-core problem could imply that the CDM paradigm needs to be revised to account for DM self-interaction (Yoshida et al. 2000; Burkert 2000; Kochanek & White 2000; Spergel & Steinhardt 2000; Davé et al. 2001; Ahn & Shapiro 2005; Vogelsberger et al. 2014b; Elbert et al. 2015; Lin & Loeb 2016), a warmer DM candidate (Colín, Avila-Reese & Valenzuela 2000; Bode, Ostriker & Turok 2001; Lovell et al. 2012) or an ultralight axion-particle (e.g. Marsh & Pop 2015), and/or a more exotic coupling between DM and dark energy (e.g. Macciò et al. 2015). Alternatively, it could simply mean that baryonic processes play a key role in shaping the inner parts of haloes and galaxies. Indeed, several astrophysical solutions have been proposed to reconcile observations with theoretical predictions. Stellar feedback could produce rapid variations of the gravitational potential through substantial gas mass outflows from the central region. This would flatten the inner density profile of the DM halo (Navarro, Eke & Frenk 1996b; Gnedin & Zhao 2002; Read & Gilmore 2005; Mashchenko, Couchman & Wadsley 2006; Mashchenko, Wadsley & Couchman 2008; Ogiya & Mori 2011, 2014; Governato et al. 2012; Pontzen & Governato 2012; Macciò et al. 2012; Teyssier et al. 2013; Ónorbe et al. 2015; Chan et al. 2015; El-Zant, Freundlich & Combes 2016; Del Popolo & Pace 2016). On the other hand, DM can also be gravitationally ‘heated’ by baryons through dynamical friction caused either by self-gravitating gas clouds orbiting near the centre of the galaxy (El-Zant, Shlosman & Hoffman 2001; El-Zant et al. 2004; Jardel & Sellwood 2009; Lackner & Ostriker 2010; Cole, Dehnen & Wilkinson 2011; Del Popolo & Pace 2016) by the presence of a stellar bar (Weinberg & Katz 2002; Holley-Bockelmann, Weinberg & Katz 2005; Sellwood 2008), by the radiation recoil from coalescing black holes (BHs, Merritt et al. 2004), or by processes which transfer of angular momentum from baryonic to DM (Tonini, Lapi & Salucci 2006; Del Popolo 2009, 2012, 2014).

On larger mass scales, numerical simulations from Peirani, Kay & Silk (2008, see also Duffy et al. 2010; Dubois et al. 2010; Teyssier et al. 2011; Martizzi et al. 2012; Martizzi, Teyssier & Moore 2013, Ragone-Figueroa & Granato 2011; Ragone-Figueroa, Granato & Abadi 2012; Ragone-Figueroa et al. 2013) have argued that active galactic nuclei (AGNs) feedback plays a similar role to that of stellar feedback in smaller systems: it can heat/expel large amounts of gas from the central regions of galaxy groups and clusters. A fraction of this gas then cools and returns to the centre, generating repeated cycles of significant inflows/outflows which in turn cause rapid fluctuations of the gravitational potential, steepening/flattening out the inner DM halo and galaxy stellar density profiles. Such a mechanism is commonly invoked to explain

the substantial body of observational evidence that the majority of massive elliptical or cD galaxies exhibit very shallow slopes in their inner (≈1 kpc) stellar surface-brightness profiles (Kormendy 1999; Quillen, Bower & Stritzinger 2000; Laine et al. 2003; Graham 2004; Trujillo et al. 2004; Lauer et al. 2005; Ferrarese et al. 2006; Côté et al. 2007; Kormendy et al. 2009; Graham 2013). A related phenomenon is the formation of cores within the central 100 pc of massive ellipticals. These are believed to be formed dynamically, by the scouring effect of possibly stalled supermassive black hole (SMBH) binaries at ~10 pc separation (see, for instance Faber et al. 1997; Thomas et al. 2014, and references therein). Larger cores are also found, for more widely separated SMBH pairs, of up to ~500 pc extent (e.g. Mazzalay et al. 2016, and references therein).

In this paper, we aim to extend previous theoretical work on the role played by AGN feedback, using a statistically representative sample spanning a comprehensive range of DM halo and galaxies masses and looking at the evolution of their inner density profiles throughout a considerable fraction of the age of the Universe. This sample is extracted from our state-of-the-art hydrodynamical cosmological simulation HORIZON-AGN (Dubois et al. 2014) which includes gas cooling, star formation, stellar and AGN feedback, and that we compare with two other simulations HORIZON-NOAGN (no AGN feedback) and HORIZON-DM (no baryons) stripped down in terms of modelled physical processes but featuring identical initial conditions. These simulations have been used to highlight the role of AGN feedback (Volonteri et al. 2016) in regulating the baryon content of massive galaxies (Kaviraj et al. 2017; Beckmann et al. 2017) and their morphological transformations (Welker et al. 2017; Dubois et al. 2016).

The paper is organized as follows. Section 2 briefly summarizes the numerical modelling upon which this work is based (simulations and post-processing). Sections 3 and 4 present our main results concerning the evolution of the inner density profiles of massive DM haloes and galaxies, respectively. Finally, we put forward and discuss our conclusions in Section 5.

# 2 NUMERICAL MODELLING

## 2.1 The three simulations $H_{AGN}$, $H_{noAGN}$ and $H_{DM}$

In this paper, we analyse and compare two large cosmological hydrodynamical simulations, HORIZON-AGN ($H_{AGN}$), HORIZON-NOAGN ($H_{noAGN}$) and one DM only cosmological simulation HORIZON-DM ($H_{DM}$). HORIZON-AGN is already described in Dubois et al. (2014), so we only summarize here its main features. We adopt a standard $\Lambda$CDM cosmology with total matter density $\Omega_m = 0.272$, dark energy density $\Omega_\Lambda = 0.728$, amplitude of the matter power spectrum $\sigma_8 = 0.81$, baryon density $\Omega_b = 0.045$, Hubble constant $H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $n_s = 0.967$ compatible with the WMAP-7. The size of the simulated volume is $L_{\text{box}} = 100 \text{ h}^{-1} \text{ Mpc}$ on a side, and it contains $1024^3$ DM particles, which results in a DM mass resolution of $M_{\text{DM, res}} = 8.27 \times 10^7 \text{ M}_\odot$. The simulation is run with the RAMSES code (Teyssier 2002), and the initially uniform grid is adaptively refined down to $\Delta x = 1$ proper kiloparsecs (kpc) at all times. Refinement is triggered in a quasi-Lagrangian manner: if the number of DM particles becomes greater than 8, or the total baryonic mass reaches 8 times the initial DM mass resolution in a cell.

Gas can radiatively cool down to $10^4 \text{ K}$ through H and He collisions with a contribution from metals using rates tabulated by Sutherland & Dopita (1993). Heating from a uniform UV background takes place after redshift $z_{\text{reion}} = 10$ following Haardt &
Madau (1996). The star formation process is modelled using a Schmidt law: $\dot{\rho}_* = \epsilon_* \rho / t_{ff}$ for gas number density above $n_0 = 0.1 \text{ H cm}^{-3}$, where $\dot{\rho}_*$ is the star formation rate density, $\epsilon_* = 0.02$ the constant star formation efficiency and $t_{ff}$ the local free-fall time of the gas. The stellar-mass resolution is $M_{*,res} = 2 \times 10^6 \text{ M}_\odot$. Feedback from stellar winds, supernovae type Ia and type II are also taken into account for mass, energy and metal release (Kimm et al. 2012).

BH formation is also included, and BHs accrete gas at a Bondi-capped-at-Eddington rate and coalesce when they form a tight enough binary. They also release energy in a quasar (heating) or radio (kinetic jet) mode when the accretion rate is above (below) one per cent of Eddington, with efficiencies tuned to match the BH-galaxy scaling relations (see Dubois et al. 2012 for detail). The presence of both quasar and radio modes is supported by recent observations. In particular, using Mapping Nearby Galaxies at APO data (Bundy et al. 2015), Cheung et al. (2016) report the presence of bi-symmetric emission features in the centre of quiescent galaxies of mass around $2 \times 10^{10} \text{ M}_\odot$ from which they infer the presence of centrally driven winds. On top of the fact that such 'red geysers' galaxies seem to be very common at this mass scale (Bundy et al. in preparation), the energy released by their SMBHs is capable of driving the observed winds and displays a mechanical content sufficient to suppress star formation. It is therefore very likely that such kinetic winds (radio mode) play a crucial role in galaxy formation and should be taken into account in numerical models (see for instance, Weinberger et al. 2017).

HORIZON-NOAGN and HORIZON-DM were performed using the same set of initial conditions and sub-grid modelling of physical processes but with no BH formation (and therefore no AGN feedback) and baryons, respectively.

## 2.2 Dark matter halo and galaxy catalogues

DM haloes are identified using the ADAPTAHOP (sub)halo finder (Aubert, Pichon & Colombi 2004; Tweed et al. 2009). In our different catalogues, host haloes and subhaloes are studied separately. Since we are particularly interested in studying the very inner part of DM haloes, a robust definition of their centre is critical. In general, the position of the most bound particle yields an accurate estimate, especially in case of HORIZON-NOAGN and HORIZON-DM as haloes in these simulations have cuspy inner profiles. However, as far as HORIZON-AGN is concerned, as we will see in the remainder of this paper, the DM profiles of haloes hosting large galaxies can be flatter. Therefore, for some (rare) object in this latter simulation, the centre can be associated with a sub-structure which is offset from the 'real' centre of mass and could lead to the attribution of a spurious core. To circumvent this issue, we use a shrinking sphere approach (Power et al. 2003) whereby starting from the virial radius, we recursively identify the centre of mass within spheres 10 per cent smaller in linear size at each iteration. We stop the procedure once the sphere reaches a 2 kpc radius and identify the centre of the halo with its densest particle. Twenty neighbours are used to compute the local density. Only structures with an average density larger than 200 times the average matter density and containing more than 100 particles become part of the (sub)halo catalogue.

Galaxies are also identified with ADAPTAHOP but applied to star particles. In this case, we use the most bound particle as the definition of the centre of mass and only galaxies with more than 50 star particles are part of the final catalogue.

Finally, we build the merger trees of all the DM haloes and galaxies for each simulation, using TREEMAKER (Tweed et al. 2009)

to link 52 outputs equally spaced in time between $z \sim 5.8$ and $z = 0$, i.e. with a time resolution of $\sim 250$ Myr.

These procedures yield, for the HORIZON-AGN run at redshift zero, about 250 000 haloes and subhaloes with a mass greater than $10^{10} \text{ M}_\odot$ and 125 000 galaxies with a mass greater than $10^8 \text{ M}_\odot$.

## 2.3 Matching dark matter haloes and galaxies

Since we start from the same initial conditions, each DM particle possesses an identity which is identical in any of the three simulations. Thus, if 75 per cent or more of the particles of any given halo in the HORIZON-AGN run also belong to a halo identified in the HORIZON-NOAGN or HORIZON-DM runs, we initially assume that these haloes are twins. However, if the mass ratio of the matching pair is greater than 10 (or lower than 0.1), we exclude it from our comparison sample. This last step is rendered necessary because sub-structures can get their particles stripped by the host halo at different times and with different intensities in the three simulations. As a result, a subhalo could become twinned with a much more massive host halo if its equivalent subhalo in the other simulation has already been destroyed (or has become too small to be detected). In general, we are able to match more than 85 per cent of DM objects at any redshift by applying these two criteria.

We cannot implement the same procedure for galaxies, since a given stellar particle is not necessarily created at the exact same moment in the very same galaxy in HORIZON-AGN and HORIZON-NOAGN. Moreover, the total number of stellar particles will also differ as it depends on the impact AGN feedback has on the star formation process. Therefore, instead of relying upon a common identity of the stellar particles they are composed of, to directly match objects between runs, we first couple each galaxy to a host DM halo in their parent simulation. We determine these galaxy-halo pairs by picking the most massive galaxy whose centre is located within a sphere of radius equal to 5 per cent of the virial radius of its host halo. Galaxy twins between runs are then determined through the matching of their host halo as previously described.

To illustrate the typical efficiency of such a procedure, at $z=0$, we are able to match about 68 per cent of H\_AGN galaxies with a mass greater $10^{10} \text{ M}_\odot$ to H\_noAGN counterparts. This fraction might seem a bit low, but our matching algorithm requires three steps to establish the link (galaxy to host halo, host halo to host halo twin, host halo twin to galaxy twin), with a number of objects dropping out of the sample at each step. Moreover, identifying a galaxy with its host halo can be challenging, especially in dense environments where interactions are more common. So in order to avoid complex situations especially during merger processes, we chose restrictive parameters. In particular, we picked the most massive galaxy whose centre is located within a sphere of radius equal to 5 per cent of the virial radius. Increasing this latter value to 10 per cent for instance would increase the number of matched galaxies (see Chisari et al. 2017). However, during very complex merger processes especially in dense environment where galaxy centre can have a big offset with respect to the centre of its host halo, it is not guaranteed to select the good galaxy. Nevertheless, we have checked that relaxing the quite stringent criteria adopted in this work improves the matching fraction, but since it also increases the number of false matches and does not alter any of our conclusions, we prefer to restrict ourselves to the more conservative sample defined in this section.
<figure><img src="image_3.png" /><figcaption>Figure 1. The number of DM haloes matched between the three simulations for our four distinct fixed mass intervals namely ~5 × 10¹¹ M☉ (dotted line), ~10¹² M☉ (dotted line), ~5 × 10¹² M☉ (dashed line) and ≥10¹³ M☉ (solid line). We derived statistics only when 10 objects can be identified at a specific redshift.</figcaption></figure>

# 3 DARK MATTER HALO DENSITY PROFILES

Let us first study the evolution of the inner density profiles of DM haloes in the three simulations.

## 3.1 Definitions

Our investigation starts with the evolution of the DM component. The questions we want to address here are twofold: (1) Does AGN feedback noticeably modify the inner density profiles of DM haloes? and (2) How does the difference between the inner density profiles of HAGN and HnoAGN (or HDM) haloes evolve? To do so, we use a systematic object-to-object comparison between the HORIZON-AGN, HORIZON-NOAGN and HORIZON-DM simulations based on the matching procedure described in the previous section. Note that when calculating the density profiles of haloes in the HDM run, we rescale the mass of the DM particles by a factor ($\Omega_m - \Omega_b$)/$\Omega_m$ to make it identical to the mass of DM particles in the baryonic runs.

In the following, we split our HAGN DM halo sample into four different mass intervals: 5 × 10¹¹(±10 per cent) M☉, 10¹²(±10 per cent) M☉, 5 × 10¹²(±10 per cent) M☉ and ≥10¹³ M☉, which we match to their HnoAGN and HDM counterparts. For each of these four mass bins, we then compute the mean density profiles (binned in spherical shells equally spaced in log r), i.e. $\rho_{AGN}(r)$, $\rho_{noAGN}(r)$ and $\rho_{DM}(r)$ at every redshift. This allows us to consider the evolution of density profiles at fixed halo mass. In addition, we have also considered the evolution of the density profiles of the progenitors of haloes within these mass bins at z = 0. Thanks to our large simulated volume, each sub-sample consists in general of 10 objects or more in the most massive bin and thousands of objects in the least massive one. However, when this is not the case (at high redshift for the most massive objects: see Figs 1 or 2), we lower the redshift until a minimum of 10 haloes of that mass can be identified and an average density profile computed.

In order to study the (relative) evolutions of the mean inner density profiles of DM haloes extracted from the three different simulations, we use two estimators. The first one is the mass-weighted density slope within $r_1$ and $r_2$ introduced by Dutton & Treu (2014):

$$ \gamma' = \frac{1}{M(r)} \int_{r_1}^{r_2} \gamma(x) 4\pi r^2 \rho(x) \, dx, \qquad (1) $$

where $\gamma \equiv d \log \rho/d \log r$ is the local logarithmic slope of the density profile $\rho$ and M the local mass. Using a discrete representation of the density profiles, we measure $\gamma(r)$ and $M(r)$ for each spherical shell centred on position r. $\gamma'$ is estimated in the range $[r_1 - r_2] = [1.0 - 5]$ kpc.

Secondly, in order to quantitatively study the evolution of the gap between the different density profiles, we compute the quantity $A_{DM}$ and $A_{noAGN}$ defined by:

$$ A_{DM|noAGN} = \int_{\log r_1}^{\log r_2} \frac{\rho_{DM|noAGN}(r) - \rho_{AGN}(r)}{\rho_{AGN}(r)} d \log r. \qquad (2) $$

$A_{DM}$ and $A_{noAGN}$ give an estimation of the gap between the mean profiles of HnoAGN and HDM haloes with respect to the mean profiles of HAGN haloes. In the following, $A_{DM}$ and $A_{noAGN}$ have been estimated in the range $[r_1 - r_2] = [1 - 10]$ kpc.

Note that the minimal value adopted for $r_1$ is dictated by our spatial resolution, which does not allow us to reach smaller scales. Spatial resolution is also the main reason why we decide not to use a fixed fraction of the virial radius to define $r_1$ (and $r_2$): we would be restricted to selecting (at least) $r_1 = 0.01 \times r_{vir}$ at high redshift because haloes are more compact and we cannot probe their inner kpc. In turn, this would then translate into $r_1 \ge 10$ kpc for galaxy cluster size haloes at z = 0, which somewhat contradicts our purpose to examine their inner density profile to the best of our ability. For these reasons, we adopt instead a simple fixed physical scale for $r_1$ and $r_2$: we probe the inner 1–5 ($\gamma'$) or 1–10 ($A_{DM|noAGN}$) kpc of DM haloes for any redshift and halo mass. The different choice of value for $r_2$ used for $\gamma'$ and $A_{DM|noAGN}$ illustrates the fact that the results presented in this paper are robust to variations of up to a factor 2 in the values of $r_1$ and $r_2$.

## 3.2 Visual inspection: a three phase’s scenario?

Fig. 2 shows the averaged spherical density profiles $\rho_{AGN}$, $\rho_{noAGN}$ and $\rho_{DM}$ derived for haloes pertaining to our different mass sub-samples and at three different redshifts. At a given time, they clearly appear different from one other especially in the central region ($r \lesssim 20$ kpc). Conversely, they seem indistinguishable at large radii, suggesting that the presence of baryons and/or AGN feedback induces noticeable effects only on small scales.

Now, if one analyses the evolution of the mean HDM, HnoAGN and HAGN halo density profiles separately, clear trends can be noticed. First, as expected from numerous previous studies, the HDM density profiles are always centrally cuspy: their inner slopes are consistent with an NFW profile (as shown in the next section). Note also that at fixed halo mass, the HDM density profiles are more extended, less concentrated, at present times than at high redshift. This well-known result is mainly due to the fact that low-redshift haloes have undergone more (major) mergers than their high-redshift analogues, and these mergers tend to diffuse material at larger radii (see, for instance Klypin et al. 2016). As a consequence, for a fixed mass interval, the density in the halo inner region progressively decreases as more and more mass is distributed at large radii. The HnoAGN halo density profiles are always much steeper and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The evolution of the mean density profiles of DM haloes extracted from $H_{\text{AGN}}$ (green lines), $H_{\text{noAGN}}$ (red lines) and $H_{\text{DM}}$ (blue lines). We focus on four distinct fixed mass intervals: ~5 × 10¹¹ M☉ (first column), ~10¹² M☉ (second column), ~5 × 10¹² M☉ (third column) and ≥10¹³ M☉ (fourth column). Three different epochs are also considered: high redshift (first line), intermediate redshift (second line) and z = 0 (third line). For indicative purposes only, the two vertical dashed lines at r = 1 and 5 kpc represent, respectively, the simulation grid size and a recommended resolution limit following Power et al. (2003). The error bars correspond to the dispersion. These plots suggest that $H_{\text{noAGN}}$ haloes have always very dense and cuspy central regions. On the contrary, AGN feedback tends to flatten the profiles especially at intermediate redshifts (z ~ [1.6–2.7]) whereas a ‘cusp regeneration’ is observed at z = 0. Finally, we also show at z = 0, the mean density profiles of DM haloes of similar mass range extracted from the Eagle simulation (Schaller et al. 2015a). The latter results suggest there are some slight differences especially for lower mass haloes (5 × 10¹¹ M☉), but quite consistent results for massive ones.</figcaption></figure>

have higher central values compared to those in the other simulations. For a fixed mass interval, no significant variation in profile shape as a function of redshift is observed. Finally, the $H_{\text{AGN}}$ density profiles present a more complex evolution. At high redshift (z ≥ 4), the halo density profiles of the two first mass bins (i.e. 5 × 10¹¹ and 10¹² M☉) are quite steep and close to their $H_{\text{noAGN}}$ matches. At intermediate redshifts (z ~ [2.7 – 1.6]), they appear to have flattened and can even have lower central density values (≤5 kpc) than to their $H_{\text{DM}}$ counterparts. Finally, at present times, the density profiles of $H_{\text{AGN}}$ haloes steepen and approach again those of their $H_{\text{noAGN}}$ twins. In view of these results, the evolution of $H_{\text{AGN}}$ halo density profiles therefore seems to follow three distinct phases.

This scenario is confirmed when studying the evolution of the progenitors of DM haloes of mass 5 × 10¹² M☉ at z = 0. The variations of their density profiles are shown at the same three redshifts in Fig. 3. $H_{\text{AGN}}$ haloes clearly present a mean density profile that is shallower at z = 1.6 than at z = 5 or 0. This strongly suggests the existence of successive phases of expansion and contraction of the mean inner DM halo density profile. Note that at z = 5, $\rho_{\text{AGN}}$ and

$\rho_{\text{noAGN}}$ are almost indistinguishable (red and green dotted lines on Fig. 3). Since the density profiles of proto-DM haloes (i.e. at very high redshifts) of the three simulations are expected to be identical (as galaxy formation has not happened yet), this indicates that $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ haloes density profiles have undergone a high identical first phase of condensation before z ~ 5. It is worth mentioning that similar trends are obtained when following the progenitors of the other halo mass sub-samples. Finally, as we will see in more detail in Section 4, the situation differs for the evolution of galaxy profiles, in the sense that $H_{\text{AGN}}$ stellar density profiles remain quite shallow at low redshift.

### 3.3 Quantitative trends

The three-redshift snapshots presented in Figs 2 and 3 provide a useful qualitative impression of the general evolution of the mean density profiles of $H_{\text{DM}}$, $H_{\text{noAGN}}$ and $H_{\text{AGN}}$ haloes. In order to derive more quantitative estimates, we now look at the evolution of the mass-weighted slope of their inner density profile, $\gamma'_{\text{dm}}$, as defined by equation (1).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Upper panel: the evolution of the mean density profiles of the progenitors of H<sub>AGN</sub> DM haloes of mass ~5 × 10<sup>12</sup> M<sub>☉</sub> at z = 0 (green lines). Three different redshifts have been considered: z = 5 (dotted line), z = 1.6 (dashed line) and z = 0 (solid line). Results for the H<sub>noAGN</sub> counterpart profiles are shown in red colours. For indicative purposes only, the two vertical dashed lines at r = 1 and 5 kpc represent, respectively, the simulation grid size and a recommended lower resolution limit following Power et al. (2003). For clarity, we do not show either the H<sub>DM</sub> counterpart profiles or all the dispersion (which are similar to those in Fig. 2). As for Fig. 2, the mean density profile of DM haloes is flatter at intermediate redshift (z ~ 1.6) and steeper at high and low z. The lower panel shows the corresponding situation for the associated galaxies. In this case, when AGN feedback is included, the profiles of galaxies progressively flatten all the way to z = 0.</figcaption></figure>

Fig. 4 shows the evolution of $\gamma'_{dm}$ for H<sub>AGN</sub> DM haloes with a fixed mass of $10^{12}$ M<sub>☉</sub> (dotted lines), $5 \times 10^{12}$ M<sub>☉</sub> (dashed lines) and $>10^{13}$ M<sub>☉</sub> (solid lines). Results for their H<sub>noAGN</sub> and H<sub>DM</sub> counterparts are also presented in the figure. First, as far as H<sub>DM</sub> haloes are concerned, values of $\gamma'_{dm}$ are very close to -1, consistent with NFW profile expectations, although we note that $\gamma'_{dm}$ is slightly increasing with time by about 0.3. As haloes become more

and more extended at low redshift, the range of [1–5] kpc is probing a relatively ‘deeper’ region in terms of fraction of the virial radius compared to that probed for haloes of the same mass at higher redshifts. Second, as expected from adiabatic contraction considerations (e.g. Blumenthal et al. 1986), the inner density profiles of H<sub>noAGN</sub> are always very steep with slopes that are close to -2. In this case a slight decrease, of a similar amplitude to the increase noted for H<sub>AGN</sub> haloes, is observed in the evolution of $\gamma'_{dm}$. Finally, H<sub>AGN</sub> haloes with a fixed mass of $10^{12}$ M<sub>☉</sub> clearly exhibit shallower inner density slopes at z ~ [1.5–1.0] than at higher or lower redshift, thus confirming the more qualitative results extracted from Fig. 2. Similar trends are observed for lower haloes masses (i.e. $5 \times 10^{11}$ M<sub>☉</sub>) but not shown in Fig. 4 for sake of clarity. More massive haloes, with masses $5 \times 10^{12}$ M<sub>☉</sub> and $>10^{13}$ M<sub>☉</sub>, generally take more time to reach their final mass, so that our analysis of their mean density profiles can only start from redshift z = 2.7 and 2.3, respectively (see Fig. 2). Due to the non-negligible past AGN activity, they have already experienced by that time, one can notice that H<sub>AGN</sub> haloes of these mass sub-samples already feature profiles that are quite flat and can even have mean slopes lower than their H<sub>DM</sub> counterparts. On the contrary, at low redshift $\rho_{AGN}$ profiles steepen considerably, with $\gamma'$ reaching values well below -1, and close to -1.5 by z = 0, thereby confirming the visual trends highlighted in Fig. 2. It is also worth mentioning that at any given redshift, more massive haloes tend to have flatter inner density profiles than low-mass ones, which confirms that AGN feedback impact increases with halo mass. To give explicit numbers, at z = 0, density profiles of haloes with masses $>10^{13}$, $5 \times 10^{12}$, $10^{12}$ and $5 \times 10^{11}$ M<sub>☉</sub> boast an inner slope value of ~-1.3, ~-1.5, ~-1.8 and ~-2.1, respectively (the latter is not shown in Fig. 4). Note that the range of [1–5] kpc is close to the simulation grid size (i.e. 1 kpc) where DM density profiles might not fully converge. This does not mean that the trends observed in Fig. 4 are necessarily wrong but one also need to consider slightly higher values in order to get a clear diagnostic. Following Power et al. (2003), a lower limit value of 5 kpc is recommended for the studied halo mass range, though their analysis concerns pure DM simulations. Thus, in the following, we also study the evolution of DM/stellar density profiles in the range of [5–10] kpc for indicative purposes only. To this effect, we also show in Fig. 4 the evolutions of $\gamma'_{dm}$ but measured this time at [5–10] kpc and one can notice that similar evolutions and trends are obtained, though the variations from H<sub>AGN</sub> haloes are less pronounced, and therefore same overall conclusions can be drawn.

Fig. 5 shows the evolution of $\gamma'_{dm}$ estimated either at [1–5] or [5–10] kpc for the progenitors of DM haloes that have a mass of $5 \times 10^{12}$ M<sub>☉</sub> either at z = 1.0 or 0.0. Focusing first on the highest redshifts (z > 2), H<sub>AGN</sub> and H<sub>noAGN</sub> profiles are clearly much steeper than H<sub>DM</sub> ones. As previously mentioned, this is mainly due to the galaxy formation process: radiative cooling and subsequent star formation lead to a steepening of the DM density profile of DM haloes due to adiabatic contraction (Blumenthal et al. 1986; Gnedin et al. 2004), in agreement with previous works (e.g. Gustafsson, Fairbairn & Sommer-Larsen 2006; Romano-Díaz et al. 2008 Abadi et al. 2010; Pedrosa, Tissera & Scannapieco 2010; Tissera et al. 2010; Duffy et al. 2010). Note that the mean slopes of H<sub>AGN</sub> and H<sub>noAGN</sub> halo density profiles are very close at z ~ 5 suggesting that AGN feedback has not had an important impact yet by this redshift. Fast forwarding to z ~ 1.6, whilst the mean profile of H<sub>noAGN</sub> haloes remains very steep with a slope close to -2, that of their H<sub>AGN</sub> counterparts is progressively flattened by AGN feedback, in good agreement with the results reported in Peirani et al. (2008) and Martizzi et al. (2013). Finally, from z ~ 1.6 to present times, $\rho_{AGN}$
<figure><img src="image_6.png" /></figure>

**Figure 4.** Time evolution of the mean mass-weighted density slope $\gamma'_{dm}$ estimated either within [1–5] kpc (upper panel) or [5–10] kpc (lower panel) of DM haloes with masses $10^{12} M_{\odot}$ (dotted lines), $5 \times 10^{12} M_{\odot}$ (dashed lines) and $>10^{13} M_{\odot}$ (solid lines). Slopes measured for HORIZON–AGN haloes and matched HORIZON–NOAGN and HORIZON–DM counterparts are coloured in green, red and blue, respectively. The three arrows indicate the times when a new refinement level is added in the simulations. For clarity, we do not show the corresponding evolution for the $5 \times 10^{11} M_{\odot}$ mass interval as it is similar to those derived for haloes of mass $10^{12} M_{\odot}$. Typical dispersions are indicated by vertical error bars at $z=0$.

becomes slightly steeper as the mean inner slope of $H_{\text{AGN}}$ haloes progressively decreases by 0.1. The two panels of Fig. 5 explain the origin of the difference in the slope of $H_{\text{AGN}}$ halo density profiles at fixed halo mass for intermediate and low redshifts displayed in Fig. 4: $H_{\text{AGN}}$ haloes tend to have steeper inner profiles at low

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

**Figure 5.** Time evolution of the mass-weighted density slope $\gamma'_{dm}$ of progenitors of DM haloes with a mass of $5 \times 10^{12} M_{\odot}$ at $z=1$ (upper panel) and $z=0.0$ (lower panel). The mean slope has been estimated either in the range of [1–5] kpc (solid lines) or [5–10] kpc (dotted lines). Results from HORIZON–AGN, HORIZON–NOAGN and HORIZON–DM simulations are coloured in green, red and blue, respectively. The three arrows indicate the times when a new refinement level is added in the simulations. The error bars correspond to the standard deviations. To facilitate comparison, filled red circles in each panel indicate the values of $\gamma'_{dm}$ when $H_{\text{AGN}}$ DM haloes reach a mass of $5 \times 10^{12} M_{\odot}$. In spite of belonging to the same mass sub-sample, $H_{\text{AGN}}$ haloes at $z=0.0$ tend to have steeper profiles than those at $z=1$.

redshift because they have undergone a phase of ‘cusp regeneration’. For completeness sake, we recall that additional refinement levels are triggered at $z \sim 5, 1.5$ and $0.26$. They induce a sudden increase in star formation and BH accretion at these different epochs as the gas can collapse to higher densities. These abrupt changes slightly affect the evolution of $\gamma'_{dm}$ especially in the case $H_{\text{noAGN}}$
haloes, as shown on Fig. 5. Finally, it is worth mentioning that no significant difference is seen in the evolutions of $\gamma'_{dm}$ when using the ranges of [1–5] kpc or [5–10] kpc.

Hence, the evolution of the inner part of $\rho_{AGN}$ does exhibit three successive phases which are quantitatively differentiable. A first condensation phase, from very high redshifts to $z \sim 5$ where dissipative galaxy formation processes and consequent adiabatic contraction of DM haloes dominate. This stage is then followed by a flattening phase, driven by important AGN activity until $z \sim 1.6$. Finally, a second condensation phase or cusp regeneration occurs from $z \sim 1.6$ down to $z = 0$ as AGN activity slowly subsides. To understand this latter phase, one can study the evolution of AGN activity in the relevant DM halo sub-samples. Indeed, an analysis of the evolution of gas mass accretion on to BHs $\dot{M}_{BH} \equiv dM_{BH}/dt$ for progenitors of DM haloes taken from each mass sub-sample and displayed in Fig. 6 reveals that for any given halo mass sample, $\dot{M}_{BH}$ and therefore AGN activity progressively decrease with time. In spite of the noticeable bumps induced by the addition of new refinement levels which artificially boost accretion, $\dot{M}_{BH}$ is one order of magnitude lower at $z = 0$ than at $z \sim 1.4$. The cusp regeneration phase measured in the evolution of $\gamma'_{dm}$ is most probably related to the fading of AGN activity at lower redshifts. As emphasized by Peirani et al. (2008) and Martizzi et al. (2013), the flattening of DM halo inner density profiles is due to a repetitive cycle of rapid gas expansions driven by AGN feedback and slower contractions as gas cools and falls back, which increases the velocity dispersion of the DM particles. If AGN activity progressively dwindles, this mechanism becomes less efficient at counterbalancing DM adiabatic contraction and thus keeping the DM density profiles flat. Fig. 6 also shows the evolution of the Eddington ratio $\chi \equiv \dot{M}_{BH}/\dot{M}_{Edd}$ (where $\dot{M}_{Edd}$ is the Eddington accretion rate) for the same haloes. These results suggest that the radio mode tends to be the dominant mode below $z \lesssim 2$. A similar behaviour is observed when considering haloes of fixed mass, rather than progenitors. Such a redshift dependence of the dominant mode is also in agreement with Volonteri et al. (2016) who have studied in detail the cosmic evolution of BHs in the HORIZON-AGN simulation. These authors found that the evolution of the luminosity function (LF) in HORIZON-AGN is consistent with the existing observational determination of the bolometric LF (Hopkins, Richards & Hernquist 2007; Shankar, Weinberg & Miralda-Escudé 2009; Ueda et al. 2014) provided there exists a transition at $z \sim 2$ between quasar- and radio-dominated mode (see their fig. 13). Recall that in the HORIZON-AGN simulation, the explosive quasar mode consists in an isotropic injection of thermal energy into the surrounding gas while at low accretion rates, the more quiescent radio mode deposits AGN feedback energy into a bipolar outflow (see Dubois et al. 2010 for a complete description). It is therefore not very surprising that such a transition also affects the evolution of inner DM halo profiles, since the radio mode of AGN feedback mainly prevents hot gas from cooling and eject little gas from the galaxy (Beckmann et al. 2017).

The cusp regeneration phase was already observed by Di Cintio et al. (2014) and Tollet et al. (2016) in hydrodynamical simulations focusing on lower mass haloes (i.e. $M_{halo} < 10^{12} M_{\odot}$). These authors claimed that the density slope of DM haloes mainly depends on the ratio between stellar mass and total halo mass $M_{star}/M_{halo}$, with large ratios ($M_{star}/M_{halo} > 0.01$) corresponding to contracted profiles. Fig. 7 shows the variation of the slope of DM halo density profiles, estimated this time within 1–2 per cent of the virial radius, as a function of the ratio $M_{star}/M_{halo}$ for $H_{AGN}$ haloes of mass $\ge 10^{13} M_{\odot}$ and matching $H_{noAGN}$ haloes. We find that in general, our haloes have values of $M_{star}/M_{halo}$ greater than 0.01 and

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Upper panel: evolution of the median mass accretion on to BHs $\dot{M}_{BH} \equiv dM_{BH}/dt$ for progenitors of DM haloes of mass $5 \times 10^{11} M_{\odot}$ (black line), $10^{12} M_{\odot}$ (green line), $5 \times 10^{12} M_{\odot}$ (red line) and $\ge 10^{13} M_{\odot}$ (blue line) at $z = 0$. Lower panel: median Eddington ratio $\chi \equiv \dot{M}_{BH}/\dot{M}_{Edd}$ where $\dot{M}_{Edd}$ is the Eddington accretion rate. Arrows indicate when an additional level of refinement is added in the simulation. For a given halo mass, the mass accretion on to BHs and therefore the AGN activity is much lower at low redshifts where the radio mode dominates.</figcaption></figure>

slopes steeper than $-1$, indicative of contraction. At $z = 0$, $H_{AGN}$ haloes have steeper density profiles and higher $M_{star}/M_{halo}$ values compared to haloes of the same mass at $z = 1$. Moreover, in the absence of AGN feedback, higher values of $M_{star}/M_{halo} \gtrsim 0.05$ are obtained which are associated with steeper profiles ($\gamma'_{dm} \sim -2$). This suggests that $\gamma'_{dm}$ is strongly correlated with $M_{star}/M_{halo}$, thus corroborating the previous findings of Di Cintio et al. (2014) and Tollet et al. (2016). However, it is worth mentioning that their
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Upper panel: variation of the mass-weighted density slope $\gamma'_{dm}$ of DM profiles within [1–2] per cent of the virial radius as a function of the ratio $M_{star}/M_{halo}$ and haloes mass for H$_{AGN}$ haloes of mass $\ge 10^{13}$ M$_{\odot}$ (green crosses) and matching H$_{noAGN}$ haloes (red crosses) identified at $z=1$ (light colours) and $z=0$ (dark colours). DM haloes with a higher $M_{star}/M_{halo}$ ratio tend to have steeper inner profiles. Lower panel: variation of $\gamma'_{dm}$ as a function of the halo mass for the same objects. In each panel, we also show the slope of the internal structure of DM haloes for groups-scale lenses (cyan squares; Newman et al. 2015) and cluster-scale lenses (blue squares; Newman et al. 2013). Finally, we also plot results from Eagle clusters (Schaller et al. 2015b) in the lower panel. Here, most massive haloes tend to have a flatter profiles and theoretical predictions seem to agree remarkably well with observations.</figcaption></figure>

simulations did not include AGN feedback and therefore are difficult to extrapolate to high-mass haloes. None the less, we remark that our results from H$_{noAGN}$ haloes are in fair agreement with the predictions of their fitting function for $M_{star}/M_{halo} > 0.05$.

Although a systematic comparison with observational data will be presented in a companion paper (Peirani et al. in preparation), it is still very instructive to make some comparison at this stage

in order to understand whether the modelled processes are in the right efficiency ballpark or not. In Fig. 7, we also show the variation of slope of the internal structure of DM haloes $\gamma'_{dm}$ derived from several groups-scale lenses with $M_{200} \sim 10^{14}$ M$_{\odot}$ at $\langle z \rangle \sim 0.36$ (Newman et al. 2015) and cluster-scale lenses with $M_{200} \sim 10^{15}$ M$_{\odot}$ at $z \sim [0.2-0.3]$ (Newman et al. 2013). Although our samples are dominated by DM haloes with a mass lower than $\sim 10^{14}$ M$_{\odot}$, the theoretical trends obtained from HORIZON–AGN are quite consistent with observational expectations namely DM haloes with a higher $M_{star}/M_{halo}$ ratio tend to have steeper inner profiles. Such results are confirmed when studying the variation of $\gamma'_{dm}$ as a function of halo mass displayed in the lower panel of Fig. 7. In this case, more massive DM haloes tend to have flatter density profiles. These figures also suggest that the inclusion of AGN feedback leads to a much better agreement with observational values and trends which is one of the main conclusion of our forthcoming study. Note also that it is quite promising that our theoretical predictions seem to be in good agreement with the Eagle simulation ones (Schaller et al. 2015b). In this regard, we only show results from the Eagle clusters, but a similar nice agreement is also obtained for lower mass range haloes (see, for instance fig. 4 of Schaller et al. 2015b).

## 3.4 Evolution of the gap between $\rho_{AGN}$, $\rho_{noAGN}$ and $\rho_{DM}$

In this section, we study the evolution of $A_{DM}$ and $A_{noAGN}$ defined by equation (2) which monitors how the relative gaps between $\rho_{AGN}$, $\rho_{noAGN}$ and $\rho_{DM}$ evolve. The variations of $A_{DM}$ and $A_{noAGN}$ are thus expected to provide complementary information to the evolution of the slope of DM halo density profiles previously studied.

### 3.4.1 $H_{AGN}$ versus $H_{DM}$

Let us start with the relative evolution of $\rho_{DM}$ and $\rho_{AGN}$. Fig. 8 shows the variations of $A_{DM}$ for the haloes of our four considered mass intervals at $z=0$ (top panel) and their progenitors (bottom panel). Three successive phases can clearly be identified again, particularly when looking at the bottom panel of Fig. 8. At very high redshift, $A_{DM}$ tends to 0 as galaxy formation has yet to significantly affect the inner structure of DM haloes. Then, from very high redshift to $z \sim 3$, $A_{DM}$ becomes progressively more and more negative due to rapid galaxy formation and subsequent DM adiabatic contraction. From $z \sim 3$ to [1.6–1.2], $A_{DM}$ is increasing whilst retaining its negative value. This suggests again that AGN activity is reducing the central density (and the inner density slope) of H$_{AGN}$ haloes which results in reducing the gap with respect to the matching H$_{DM}$ haloes density profiles. Finally, from $z \sim [1.6-1.2]$ to $z=0$, $A_{DM}$ decreases again, always remaining negative. The gap between the inner density profiles of H$_{AGN}$ and H$_{DM}$ haloes thus progressively becomes more and more important, since this time the central density (and the inner density slope) of H$_{AGN}$ haloes is increasing (phase of cusp regeneration). Note that the second phase (i.e. flattening of H$_{AGN}$ DM halo density profiles) occurs sooner and lasts longer for more massive objects. Furthermore, the third phase occurs later for more massive objects. This reinforces our conclusions based on the evolution of inner profile density slopes that AGN feedback impacts more durably and significantly the most massive haloes.

A similar conclusion can be drawn when studying the time evolution of $A_{DM}$ for haloes with a fixed mass, independent of cosmic time, as displayed in the top panel of Fig. 8. $A_{DM}$ is increasing at high redshifts ($5 > z > 2$) for the first three lower mass intervals while keeping negative values. These black, green and red curves
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Time evolution of $A_{\text{DM}}$ which measures the gap between $H_{\text{AGN}}$ and matched $H_{\text{DM}}$ DM halo density profiles in the range of [1–10] kpc. The upper panel considers haloes within fixed mass intervals, independently of cosmic time, while the lower one follows the evolution of $A_{\text{DM}}$ for the progenitors of haloes of mass ~5 × 10<sup>11</sup>, ~10<sup>12</sup>, ~5 × 10<sup>12</sup> and ≥10<sup>13</sup> M<sub>☉</sub> at z = 0. The vertical dashed lines indicate when additional levels of refinement are triggered. The shaded areas represent the error on the mean. Note that similar trends are obtained when estimating $A_{\text{DM}}$ in the range of [1–5] or [5–10] kpc.</figcaption></figure>

reach a maximal value at redshifts around 1.6, 1.6 and 1.2, respectively: AGN activity has efficiently reduced the central density as well as flattened the profiles of $H_{\text{AGN}}$ haloes, leading to the reduction of the gap between inner density profiles of $H_{\text{AGN}}$ and matching $H_{\text{DM}}$ haloes. Note that for these three mass intervals, $A_{\text{DM}}$ is already

negative at high redshift due to DM adiabatic contraction induced by galaxy formation having already taken place in such massive haloes. Then, after each peak, $A_{\text{DM}}$ monotonically decreases until $z=0$ (phase of cusp regeneration). For the most massive haloes ($\ge 10^{13}$ M<sub>☉</sub>), $A_{\text{DM}}$ is mainly decreasing but remains positive until $z \sim 1$ indicating that $H_{\text{AGN}}$ haloes of this mass at higher redshifts have lower central density (and also flatter density profiles) than their $H_{\text{DM}}$ counterparts. This also reinforces that AGN activity has a stronger impact on the inner DM density profiles of the most massive haloes.

### 3.4.2 $H_{\text{AGN}}$ versus $H_{\text{noAGN}}$

AGN feedback manifests itself through two effects. First, it reduces star formation rates in massive galaxies and therefore the amount of stellar material present in their central parts. Second, as we have previously discussed, it flattens the DM host halo density profile, especially at intermediate redshifts. As a result, the central DM profiles of $H_{\text{AGN}}$ haloes should always be less steep and have lower central values when compared with their $H_{\text{noAGN}}$ halo counterparts. In other words, $A_{\text{noAGN}}$ must always be positive.

Fig. 9 shows the variation of $A_{\text{noAGN}}$ for both progenitors and haloes within fixed mass intervals and confirms again the trends previously pointed out. In the top panel, for the three less massive halo mass intervals, the gap between the DM halo density profiles in $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ first increases until it reaches a maximum value for $z$ close to 1.2. It then slightly decreases as the redshift progresses towards $z=0$. For the most massive haloes, with masses $\ge 10^{13}$ M<sub>☉</sub>, the gap is essentially always increasing until $z = 0.8$ after which it remains constant until $z = 0$. This is consistent with our previous result: the effect of AGN feedback (quasar mode) is more efficient at high redshift, leading to a rapid increase of the gap between the $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ DM halo profiles. Then, in a longer and later phase between 1.2 $\ge z \ge 0$, AGN activity diminishes, and the gap narrows slightly.

The bottom panel of Fig. 9 presents the variations of $A_{\text{noAGN}}$ focusing on the evolution of the DM halo profiles of the progenitors of haloes with a mass of $5 \times 10^{11}$, $10^{12}$, $5 \times 10^{12}$ and $\ge 10^{13}$ M<sub>☉</sub> at $z=0$. Here again, at very high redshift, $A_{\text{noAGN}}$ is close to 0 for all considered mass intervals, since AGN feedback has not yet kicked in. The gap between $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ haloes then increases until $z \sim [1-0.6]$ and stays roughly constant afterwards. It is worth mentioning again that adding an extra refinement level at $z \sim 1.5$ and 0.25 causes spurious artefacts to appear in the otherwise rather smooth evolution of $A_{\text{noAGN}}$. This corresponds to a sudden better force resolution and enhanced gas condensation which has a more dramatic effect in HORIZON-NOAGN than in HORIZON-AGN, where it is somewhat compensated by a simultaneous rise in AGN activity.

# 4 STELLAR DENSITY PROFILES

## 4.1 General trends

In the following, we compare the stellar density profiles of galaxies in the same way as we previously did for the DM component: we use a systematic (AGN) galaxy to (noAGN) galaxy object-to-object comparison using our matching procedure. We also focus our analysis on four ($H_{\text{AGN}}$) galaxy stellar-mass intervals: $10^9 (\pm 10 \text{ per cent})$, $10^{10} (\pm 10 \text{ per cent})$, $10^{11} (\pm 10 \text{ per cent})$ and $\ge 5 \times 10^{11}$ M<sub>☉</sub>. Fig. 10 indicates the number of galaxies identified in each mass bin and for each considered redshift.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Time evolution of $A_{\text{noAGN}}$, which measures the gap between $H_{\text{AGN}}$ and matching $H_{\text{noAGN}}$ DM halo density profiles in the range of [1–10] kpc. Upper and lower panels consider haloes within fixed mass intervals across cosmic time and the progenitors of haloes with masses of ~5 × 10¹¹, ~10¹², ~5 × 10¹² and ≥10¹³ $M_{\odot}$ at $z=0$ respectively. The vertical dashed lines indicate when additional levels of refinement are triggered in the simulation. Shaded areas represent the error on the mean. Here again, similar trends are obtained when estimating $A_{\text{noAGN}}$ in the range of [1–5] or [5–10] kpc.</figcaption></figure>

Fig. 11 displays spherically averaged galaxy stellar-mass density profiles, $\rho_{*,\text{AGN}}$ and matching $\rho_{*,\text{noAGN}}$, derived for each of the four stellar-mass bins and at three different redshifts. We again indicate in each panel, the simulation grid size (1 kpc) and a lower resolution limit (5 kpc) recommended by Power et al. (2003). For the lowest

<figure><img src="image_12.png" /><figcaption>**Figure 10.** The number of galaxies matched between the two hydrodynamical simulations for our four distinct fixed mass intervals namely ~10⁹ $M_{\odot}$ (dot-dashed line), ~10¹⁰ $M_{\odot}$ (dotted line), ~10¹¹ $M_{\odot}$ (dashed line) and ≥5 × 10¹¹ $M_{\odot}$ (solid line). We derived statistics only when 10 objects can be identified at a specific redshift.</figcaption></figure>

mass interval (i.e. 10⁹ $M_{\odot}$), we notice no significant difference between $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ galaxy density profiles: this is not very surprising as AGN feedback is not thought to be effective in such low-mass objects. For the three most massive mass intervals, clear gaps between $\rho_{*,\text{AGN}}$ and $\rho_{*,\text{noAGN}}$ can be seen in the central regions, and up to a radius of 50 kpc for the most massive galaxies. Moreover, the gap between inner density profiles is all the more important than the galaxies are massive, and increases slightly with time. Similarly to what we observe for the DM component, stellar density profiles with and without AGN feedback are similar at large radii. However, in contrast to what happens for DM haloes, $H_{\text{AGN}}$ galaxy density profiles remain quite flat at low redshift, as already suggested by Fig. 3. Such a behaviour is also more consistent with observations (e.g. Kormendy et al. 2009).

Note that, at a given mass, galaxies are in general more extended at low redshifts because they experience more mergers that tend to spread material at large radii, which increase their effective radius over time (Khochfar & Silk 2006; Bournaud, Jog & Combes 2007; Naab, Johansson & Ostriker 2009; Peirani et al. 2010; Oser et al. 2010, 2012; Shankar et al. 2013; Welker et al. 2017; Rodriguez-Gomez et al. 2016). This effect is more pronounced with AGN feedback since the *in situ* star formation is regulated by AGN activity, at the benefit of the accreted stellar mass in the overall stellar-mass budget (Dubois et al. 2013, 2016).

In order to derive a more quantitative evolution, we first study the evolution of the galaxy mass-weighted density slope $\gamma'_*$ measured within [1–5] and [5–10] kpc for $H_{\text{AGN}}$ galaxies within a fixed mass bin, independent of cosmic time, and their matching $H_{\text{noAGN}}$ counterparts. Those evolutions are displayed in Fig. 12. For galaxies with masses of 10¹⁰, 10¹¹ and >5 × 10¹¹ $M_{\odot}$ $H_{\text{AGN}}$, we see clear differences when AGN is included or not. Stellar density profiles of $H_{\text{AGN}}$ galaxies always display shallower inner slopes than their $H_{\text{noAGN}}$ counterparts. Therefore, as was measured for the DM
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the mean stellar density profiles of galaxies extracted from H<sub>AGN</sub> (green lines) and H<sub>noAGN</sub> (matching, red lines). We divide our galaxy sample on four fixed mass intervals: ~10<sup>9</sup> M<sub>☉</sub> (first column), ~10<sup>10</sup> M<sub>☉</sub> (second column), ~10<sup>11</sup> M<sub>☉</sub> (third column) and ≥5 × 10<sup>11</sup> M<sub>☉</sub> (fourth column). As for the DM haloes, three different epochs are considered: high redshift (first row), z = 0.7 (second row) and z = 0 (third row). The two vertical dashed lines at r = 1 and 5 kpc are also displayed to indicate, respectively, the simulation grid size and a recommended resolution limit following Power et al. (2003). The error bars correspond to the dispersion. The discrepancy between density profiles is all the more important than galaxies are massive. Moreover, the mean H<sub>AGN</sub> stellar density profiles remain quite flat at low redshift.</figcaption></figure>

component, AGN feedback tends to flatten the inner stellar density profiles of massive galaxies. We also note that H<sub>AGN</sub> mean stellar density slopes increase rapidly at high redshift (z > 1.5) and then stall at lower z, as a consequence of the evolution of BHs growth and AGN activity history reported in Fig. 6. Finally, the most massive H<sub>AGN</sub> galaxies tend to have a flatter inner profiles at any given redshift than their lower mass equivalents, which suggest again that AGN feedback plays a more important role in shaping more massive objects. Note that no particular difference are obtained when estimating γ<sub>*</sub> within [1–5] or [5–10] kpc.

We derived in Fig. 13, the evolution of γ' (estimated within either at [1–5] or [5–10] kpc) for the progenitors of galaxies of mass 10<sup>11</sup> M<sub>☉</sub> at z = 0. At high redshift (z > 5), since AGN activity has not yet picked up, the density slopes of H<sub>AGN</sub> and H<sub>noAGN</sub> galaxies are very similar. Then, from z ~ 5 to 1, AGN feedback strongly flattens the density profiles of H<sub>AGN</sub> galaxies. At lower redshifts, the inner stellar density profile slopes remain almost constant or increase slightly which confirms trends seen in Fig. 3. It is worth mentioning that similar results are also obtained for progenitors of galaxies with a mass of ≥5 × 10<sup>11</sup> M<sub>☉</sub> at z = 0. Note again that the evolutions of γ' estimating within [1–5] or [5–10] kpc are very similar.

Finally, Fig. 14 shows the time evolution of $A^*_{\text{noAGN}}$ for galaxies within a mass interval fixed throughout cosmic time (right-hand panels). We use the same definition of $A_{\text{noAGN}}$ (but for galaxies) given by equation (2) and consider this time $r_1 = 1$ kpc and $r_2 = 5$ kpc. As expected, for 10<sup>9</sup> M<sub>☉</sub> galaxies, $A_{\text{noAGN}}$ stays constant and equal to 0. On the contrary, for 10<sup>10</sup> M<sub>☉</sub> galaxies, $A^*_{\text{noAGN}}$ is always positive and slightly decreasing below z ~ 2. Although AGN activity is relatively weak in galaxies pertaining to this mass range, it still affects their stellar density profiles all the way to the present time. More massive galaxies, with masses of 10<sup>11</sup> M<sub>☉</sub> feature a rapid increase in $A^*_{\text{noAGN}}$ until z ~ 1 and this quantity then remains roughly constant between z = 1 and 0. Finally, $A^*_{\text{noAGN}}$ is monotonously increasing for the most massive galaxies (≥5 × 10<sup>11</sup> M<sub>☉</sub>). This is explained by two main reasons. First, AGN activity rapidly flattens H<sub>AGN</sub> galaxy density profiles at high redshift, and second, the mass of H<sub>noAGN</sub> galaxies is still noticeably increasing down to z = 0, as there is no upper limit in this mass range. Fig. 14 also presents the evolution of the progenitors of these galaxies (left-hand panels). In this case, $A^*_{\text{noAGN}}$ is always increasing (albeit more rapidly for more massive galaxies) which means that the gap between H<sub>AGN</sub> and H<sub>noAGN</sub> galaxy density profiles is continuously increasing. Note
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of the galaxy mass-weighted density slope $\gamma_{*}'$ estimated within [1–5] kpc (upper panel) and [5–10] kpc (lower panel). We show results for the three more massive $H_{\text{AGN}}$ galaxy mass intervals (green colour) and matched $H_{\text{noAGN}}$ galaxies (red colour). The three arrows indicate the times when an additional refinement level is added to the simulations. Typical standard deviations are represented by vertical error bars at $z=0$. AGN feedback significantly and monotonically flattens the inner stellar density profiles of galaxies.</figcaption></figure>

again that the additional refinement levels have a limited, but noticeable impact on the evolution of $A_{\text{noAGN}}^*$. The extra star formation spuriously induced at these epochs increases the central stellar mass, especially in $H_{\text{noAGN}}$ galaxies, which causes the more pronounced increases seen in the evolution of $A_{\text{noAGN}}^*$. However these numerical effects do not affect our conclusions.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Time evolution of the mass-weighted density slope $\gamma_{*}'$ of the progenitors of galaxies with masses $10^{11} M_{\odot}$ at $z=0$ estimating within [1–5] kpc (solid lines) or [5–10] kpc (dotted line). Results from $H_{\text{AGN}}$ galaxies and their matching $H_{\text{noAGN}}$ counterparts are displayed in green and red colours, respectively. The three arrows represent epochs where an additional refinement level is added in the simulations. Error bars indicate the standard deviations. $H_{\text{AGN}}$ galaxy stellar density profiles are rapidly flattened until $z \sim 1$. Without AGN feedback, galaxy density profiles always remain relatively steep.</figcaption></figure>

## 4.2 Matching or no matching, what is the difference?

All along the paper so far, we have only carried out object-to-object comparisons to ensure that we study the properties and evolution of the same objects in the different simulations. However, one can also be interested in comparing the evolution of the properties of objects of similar mass between the three simulations, which is not enforced by the matching strategy. In this section we examine what the difference is between these two approaches.

As far as the DM component is concerned, the answer is none, because AGN feedback does not significantly affect the virial mass of DM haloes. This is demonstrated in Fig. 15 where we plot the mean density profiles of $H_{\text{noAGN}}$ DM haloes either at $z=1$ or 0, either obtained through matching or simply considering the population in the same mass interval: the two are indistinguishable from one another.

On the contrary, stellar density profiles of galaxies prove to be very different. Indeed, for any given stellar-mass interval of $H_{\text{AGN}}$ galaxies, the $H_{\text{noAGN}}$ counterparts are in general much more massive. These galaxies display stellar density values that are much higher in the central parts but similar at large radii to the $H_{\text{AGN}}$ ones. Therefore, when comparing $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ galaxies of the same mass, Fig. 16 clearly shows that these latter still present higher central stellar density values. However, in order to compensate for the extra mass enclosed within these regions, they also exhibit lower densities at large radii. Thus, for a galaxy of a given mass, AGN feedback leads to density profiles that are more centrally flat, but also more extended.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Time evolution of $A_{\text{noAGN}}^*$ which measures the gap between $H_{\text{AGN}}$ and $H_{\text{noAGN}}$ galaxy density profiles in the range of [1–5] kpc. The left-hand panels consider galaxies within fixed mass intervals throughout cosmic time: $10^9$, $10^{10}$ and $10^{11} M_{\odot}$ (upper panel) and $\ge 5 \times 10^{11} M_{\odot}$ (lower panel). In the right-hand panels, we follow the evolution of the progenitors of galaxies that have a ($H_{\text{AGN}}$) mass of $\sim 10^{10}$, $10^{11}$ and $\ge 5 \times 10^{11} M_{\odot}$ at $z=0$. Shaded areas represent the error on the mean. Vertical dashed lines indicate when an additional level of refinement is introduced in the simulations. Note that we have multiplied results obtained for $10^{10} M_{\odot}$ galaxies by a factor 3 for clarity. It is worth mentioning that similar evolutions are obtained when estimating $A_{\text{noAGN}}^*$ in the range of [1–10] or [5–10] kpc. Only the amplitudes will of course change.</figcaption></figure>

## 5 DISCUSSION AND CONCLUSIONS

By comparing results from two state-of-the-art hydrodynamical cosmological simulations whose only difference is the presence/absence of AGN feedback, and one cosmological simulation without baryons which otherwise shares the same initial conditions, we have explored the impact of AGN feedback on the evolution of the inner density profiles of massive DM haloes and galaxies. We focused on DM haloes and galaxies with a mass greater than $5 \times 10^{11}$ and $10^9 M_{\odot}$, respectively. Since the resolution limit of the simulations is 1 kpc (physical), we have only investigated the (relative) variations of halo and galaxy density profiles within a few kpc from their centre (i.e. [1–5] or [1–10] kpc). Our findings can be summarized as follows:

(i) When AGN feedback is included, the mean inner density profiles $\rho_{\text{AGN}}$ of massive DM haloes undergo successive phases of contraction (steepening) and expansion (flattening). From very high redshift to $z \sim 3$, $\rho_{\text{AGN}}$ becomes steeper than $\rho_{\text{DM}}$ due to adiabatic contraction induced by early galaxy formation in the centre of the host DM haloes. From $z \sim 3$ down to $z \sim 1.6$, $\rho_{\text{AGN}}$ is noticeably flattened by AGN activity which is high (quasar mode). From $z \sim 1.6$ to the present time, $\rho_{\text{AGN}}$ steepens again ('cusp regeneration') as AGN activity considerably reduces.

(ii) The gaps between $\rho_{\text{AGN}}$, $\rho_{\text{noAGN}}$ and $\rho_{\text{DM}}$ are also evolving with time. At high redshift, $\rho_{\text{AGN}}$ and $\rho_{\text{noAGN}}$ tend to be much steeper than $\rho_{\text{DM}}$ due to rapid galaxy formation. Until $z \sim 1.6$, the flattening of $\rho_{\text{AGN}}$ tends to increase the gap with $\rho_{\text{noAGN}}$ and conversely decrease the gap with $\rho_{\text{DM}}$. Finally after $z \sim 1.6$, the

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Mean density profiles of DM haloes with a mass of $5 \times 10^{12} M_{\odot}$ at $z=1$ (upper panel) and $z=0$ (lower panel) and extracted from the HORIZON-AGN (green) or HORIZON-NOAGN (blue) simulations. The DM density profile of matched haloes is plotted in red. The error bars correspond to the standard deviations. Since AGN feedback does not affect the virial mass of DM haloes, the red and blue profiles sample the same population of haloes and are therefore indistinguishable.</figcaption></figure>

phase of 'cusp regeneration' leads to both a slight reduction and increase of the gap with respect to $\rho_{\text{noAGN}}$ and $\rho_{\text{DM}}$, respectively.

(iii) AGN feedback noticeably reduces the central density in massive galaxies and efficiently flattens their inner profiles, which lead to trends more compatible with the observations of massive elliptical or cD galaxies, that exhibit very shallow slopes in the stellar surface-brightness profiles within small radii ($\approx 1$ kpc, Kormendy 1999; Quillen et al. 2000; Laine et al. 2003; Graham 2004; Trujillo et al. 2004; Lauer et al. 2005; Ferrarese et al. 2006; Côté et al. 2007; Kormendy et al. 2009; Graham 2013). In contrast to the DM component, galaxy inner density profiles remain quite flat at low redshifts.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Mean stellar density profiles of galaxies with a mass of $10^{11}$ M$_\odot$ at $z=1$ (upper panel) and $z=0$ (lower panel) and extracted from the HORIZON–AGN (green) or HORIZON–NOAGN (blue) simulations. The stellar density profile of the matched galaxies is plotted in red. The error bars correspond to the standard deviations. Contrary to the DM component, the presence of AGN feedback significantly affects the virial mass of the galaxies. In this case, the red and blue profiles are not derived from the same population of galaxies and are therefore very different.</figcaption></figure>

(iv) At any given redshift, more massive DM haloes or galaxies have in general flatter central density profiles than their less massive counterparts. The impact of AGN feedback in the flattening of DM haloes/galaxies density profiles is all the more important than the objects are massive.

(v) Without AGN feedback, the inner density profiles of DM haloes and galaxies are always very steep.

This study clearly demonstrates that the inner density profiles of DM haloes and galaxies are very sensitive to sub-grid physics and more specifically to AGN feedback. For instance, our model predicts a three-phase scenario in the evolution of DM density profiles which is intimately associated with the strength of AGN activity via the contribution of the different accretion modes on the

central BH. Indeed, as advocated by Peirani et al. (2008), repetitive cycles of gas expansion by AGN feedback and gas cooling are requested to efficiently flatten the DM profiles. If the AGN activity progressively decreases, this proposed mechanism becomes less efficient at counterbalancing the DM adiabatic contraction and at keeping the DM density profiles flat. Furthermore, we also found that the radio mode tends to be the dominant mode after $z \sim 2$, which renders the expansion phase of the gas more difficult. In parallel, AGN feedback is also expected to regulate the star formation in the objects studied here (see Beckmann et al. 2017 for detail). In this regard, Kaviraj et al. (2017) have recently studied the reproduction of quantities in HORIZON–AGN that trace the aggregate stellar-mass growth of galaxies over cosmic time namely the luminosity and stellar-mass functions, the star formation main sequence, rest-frame UV-optical-near-infrared colours and the cosmic star formation history. They found that HORIZON–AGN successfully captures the evolutionary trends of observed galaxies over the lifetime of the Universe, making it an excellent tool for studying the processes that drive galaxy evolution. However, it is worth mentioning that although the galaxy stellar-mass functions at different redshift are consistent at the high-mass end, HORIZON–AGN tend to overshoots the low-mass end ($< 5 \times 10^{10}$ M$_\odot$). Similar conclusions were obtained by Welker et al. (2017) by comparing galaxy stellar-mass functions with observation of Cosmic Assembly Near-infrared Deep Extragalactic Legacy survey - Ultra Deep Survey and the Great Observatories Origins Deep Survey. Future accurate observations will allow us to confirm these theoretical predictions, but, above all, will help constraining AGN models. One complication may arise when trying to probe the distribution of the DM component. If one assumes that the stellar-mass distribution will trace that of the DM then high-resolution observations of large samples of massive galaxies across a large-redshift range are requested. This could be done using *Euclid* in the redshift range $z < 1$ and James Webb Space Telescope in the redshift range $1 < z < 4$.

Basic comparisons with other theoretical works have also been considered in this work. For instance, the mean density profiles of DM haloes derived from the Eagle simulation (Schaye et al. 2015) how noticeable differences especially for lower mass haloes ($5 \times 10^{11}$ M$_\odot$) but more consistent results for the massive ones. Moreover, we found that the variations of the slope of the internal structure of DM haloes $\gamma'_{dm}$ are also in nice agreement with those derived from Eagle data in spite different hydrodynamic solvers and AGN model implementations. Finally, our finding regarding the cusp regeneration phase is in good agreement results derived by Di Cintio et al. (2014) and Tollet et al. (2016) in hydrodynamical simulations focusing on lower mass haloes (i.e. $M_{halo} < 10^{12}$ M$_\odot$). Further relevant comparison, with the Illustris simulation (Vogelsberger et al. 2014a; Genel et al. 2014) for instance, will be done in a companion paper (Peirani et al. in preparation) focusing on the total density slope of massive early-type galaxies.

We mainly focused on the inner part of DM haloes density profiles where clear differences can been seen between the different simulations. Conversely, for a given mass sample, the density profiles of H$_{AGN}$ haloes and their matching H$_{noAGN}$ and H$_{DM}$ haloes converge and become identical at distances $r > 10–20$ kpc, which suggests that these outer regions are not affected by baryons and AGN feedback. However, previous works indicate that AGN feedback is also expected to produce effect on larger scales. For instance, AGN feedback is essential to produce massive galaxies that resemble ellipticals. Without BH activity, massive galaxies are disc-like with kinematics dominated by rotational support (see Dubois et al. 2013, 2016). At even larger scales, Suto et al. (2017) have recently
examined the asphericity of galaxy clusters using the projected axis ratios of X-ray surface brightness, star and DM distributions of the most massive haloes ($M_{200} > 5 \times 10^{13} M_{\odot}$) extracted from HORIZON-AGN, HORIZON-NOAGN and HORIZON-DM. They found that the baryonic physics and especially AGN feedback can significantly affect the asphericity of DM distribution even beyond the central region, approximately up to the half of the virial radius. AGN feedback seems therefore an indispensable ingredient for the formation of massive DM haloes and galaxies.

In a companion paper, we will investigate the effect of AGN feedback on the total (DM + stars) density slope of massive early-type galaxies estimated at the effective radius.

## ACKNOWLEDGEMENTS

We warmly thank the referee for an insightful review that considerably improved the quality of the original manuscript. We also warmly thank Y. Suto and T. Nishimichi for interesting discussions and M. Schaller for providing relevant Eagle simulation data. SP acknowledges support from the Japan Society for the Promotion of Science (JSPS long-term invitation fellowship). This work was granted access to the High Performance Computing resources of Centre Informatique National de l'Enseignement Supérieur under the allocations 2013047012, 2014047012 and 2015047012 made by Grand Équipement National de Calcul Intensif and has made use of the Horizon cluster hosted by the Institut d'Astrophysique de Paris on which the simulation was post-processed. This work was carried out within the framework of the Horizon project (http://www.projet-horizon.fr) and is partially supported by the grants ANR-13-BS05-0005 of the French Agence Nationale de la Recherche. The research of JD is supported by Adrian Beecroft and STFC.

## REFERENCES

* Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435

* Ahn K., Shapiro P. R., 2005, MNRAS, 363, 1092

* Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

* Beckmann R. S. et al., 2017, preprint (arXiv:1701.07838)

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

* Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93

* Bournaud F., Jog C. J., Combes F., 2007, A&A, 476, 1179

* Bundy K. et al., 2015, ApJ, 798, 7

* Burkert A., 2000, ApJ, 534, L143

* Chan T. K., Kereš D., Ónorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 454, 2981

* Cheung E. et al., 2016, Nature, 533, 504

* Chisari N. E. et al., 2017, preprint (arXiv:1702.03913)

* Cole D. R., Dehnen W., Wilkinson M. I., 2011, MNRAS, 416, 1118

* Colín P., Avila-Reese V., Valenzuela O., 2000, ApJ, 542, 622

* Côté P. et al., 2007, ApJ, 671, 1456

* Davé R., Spergel D. N., Steinhardt P. J., Wandelt B. D., 2001, ApJ, 547, 574

* de Blok W. J. G., McGaugh S. S., Bosma A., Rubin V. C., 2001, ApJ, 552, L23

* de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S.-H., Kennicutt R. C., Jr, 2008, AJ, 136, 2648

* Del Popolo A., 2009, ApJ, 698, 2093

* Del Popolo A., 2012, MNRAS, 419, 971

* Del Popolo A., 2014, J. Cosmol. Astropart. Phys., 7, 019

* Del Popolo A., Pace F., 2016, Ap&SS, 361, 162

* Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415

* Dubois Y., Devriendt J., Slyz A., Teyssier R., 2010, MNRAS, 409, 985

* Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

* Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

* Dubois Y. et al., 2014, MNRAS, 444, 1453

* Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

* Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161

* Dutton A. A., Treu T., 2014, MNRAS, 438, 3594

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

* El-Zant A. A., Hoffman Y., Primack J., Combes F., Shlosman I., 2004, ApJ, 607, L75

* El-Zant A. A., Freundlich J., Combes F., 2016, MNRAS, 461, 1745

* Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Ónorbe J., Peter A. H. G., 2015, MNRAS, 453, 29

* Faber S. M. et al., 1997, AJ, 114, 1771

* Ferrarese L. et al., 2006, ApJS, 164, 334

* Genel S. et al., 2014, MNRAS, 445, 175

* Gentile G., Salucci P., Klein U., Vergani D., Kalberla P., 2004, MNRAS, 351, 903

* Gnedin O. Y., Zhao H., 2002, MNRAS, 333, 299

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

* Goerdt T., Moore B., Read J. I., Stadel J., Zemp M., 2006, MNRAS, 368, 1073

* Governato F. et al., 2012, MNRAS, 422, 1231

* Graham A. W., 2004, ApJ, 613, L33

* Graham A. W., 2013, Planets, Stars and Stellar Systems. Volume 6: Extra-galactic Astronomy and Cosmology. Springer, Dordrecht, p. 91

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

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

* Holley-Bockelmann K., Weinberg M., Katz N., 2005, MNRAS, 363, 991

* Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731

* Jardel J. R., Sellwood J. A., 2009, ApJ, 691, 1300

* Jing Y. P., Suto Y., 2000, ApJ, 529, L69

* Kaviraj S. et al., 2017, MNRAS, 467, 4739

* Khochfar S., Silk J., 2006, ApJ, 648, L21

* Kimm T. et al., 2012, MNRAS, 425, L96

* Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

* Kochanek C. S., White M., 2000, ApJ, 543, 514

* Kormendy J., 1999, Galaxy Dynamics – A Rutgers Symposium, Vol. 182. Astron. Soc. Pac., San Francisco

* Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216

* Lackner C. N., Ostriker J. P., 2010, ApJ, 712, 88

* Laine S., van der Marel R. P., Lauer T. R., Postman M., O'Dea C. P., Owen F. N., 2003, AJ, 125, 478

* Lauer T. R. et al., 2005, AJ, 129, 2138

* Leonard A., Goldberg D. M., Haaga J. L., Massey R., 2007, ApJ, 666, 51

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

* Lin H. W., Loeb A., 2016, J. Cosmol. Astropart. Phys., 3, 009

* Lovell M. R. et al., 2012, MNRAS, 420, 2318

* Macciò A. V., Stinson G., Brook C. B., Wadsley J., Couchman H. M. P., Shen S., Gibson B. K., Quinn T., 2012, ApJ, 744, L9

* Macciò A. V., Mainini R., Penzo C., Bonometto S. A., 2015, MNRAS, 453, 1371

* Marsh D. J. E., Pop A.-R., 2015, MNRAS, 451, 2479

* Martizzi D., Teyssier R., Moore B., Wentz T., 2012, MNRAS, 422, 3081

* Martizzi D., Teyssier R., Moore B., 2013, MNRAS, 432, 1947

* Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539

* Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174

* Mazzalay X., Thomas J., Saglia R. P., Wegner G. A., Bender R., Erwin P., Fabricius M. H., Rusli S. P., 2016, MNRAS, 462, 2847

* Merritt D., Milosavljević M., Favata M., Hughes S. A., Holz D. E., 2004, ApJ, 607, L9

* Merritt D., Graham A. W., Moore B., Diemand J., Terzić B., 2006, AJ, 132, 2685

* Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ, 499, L5
Naab T., Johansson P. H., Ostriker J. P., 2009, ApJ, 699, L178

Navarro J. F., Frenk C. S., White S. D. M., 1996a, ApJ, 462, 563

Navarro J. F., Eke V. R., Frenk C. S., 1996b, MNRAS, 283, L72

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

Navarro J. F. et al., 2010, MNRAS, 402, 21

Newman A. B., Treu T., Ellis R. S., Sand D. J., Richard J., Marshall
P. J., Capak P., Miyazaki S., 2009, ApJ, 706, 1078

Newman A. B., Treu T., Ellis R. S., Sand D. J., 2011, ApJ, 728, L39

Newman A. B., Treu T., Ellis R. S., Sand D. J., 2013, ApJ, 765, 25

Newman A. B., Ellis R. S., Treu T., 2015, ApJ, 814, 26

Ogiya G., Mori M., 2011, ApJ, 736, L2

Ogiya G., Mori M., 2014, ApJ, 793, 46

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt R. C., Jr, 2011, AJ, 141, 193

Oldham L. J., Auger M. W., 2016, MNRAS, 457, 421

Oñorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Kere D., Faucher-Giguère C.-A., Quataert E., Murray N., 2015, MNRAS, 454, 2092

Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312

Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, ApJ, 744, 63

Palunas P., Williams T. B., 2000, AJ, 120, 2884

Pedrosa S., Tissera P. B., Scannapieco C., 2010, MNRAS, 402, 776

Peirani S., Kay S., Silk J., 2008, A&A, 479, 123

Peirani S., Crockett R. M., Geen S., Khochfar S., Kaviraj S., Silk J., 2010, MNRAS, 405, 2327

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel
V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Quillen A. C., Bower G. A., Stritzinger M., 2000, ApJS, 128, 85

Ragone-Figueroa C., Granato G. L., 2011, MNRAS, 414, 3690

Ragone-Figueroa C., Granato G. L., Abadi M. G., 2012, MNRAS, 423, 3243

Ragone-Figueroa C., Granato G. L., Murante G., Borgani S., Cui W., 2013, MNRAS, 436, 1750

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Richtler T., Salinas R., Misgeld I., Hilker M., Hau G. K. T., Romanowsky
A. J., Schuberth Y., Spolaor M., 2011, A&A, 531, A119

Rodriguez-Gomez V. et al., 2016, MNRAS, 458, 2371

Romano-Díaz E., Shlosman I., Hoffman Y., Heller C., 2008, ApJ, 685, L105

Salucci P., Burkert A., 2000, ApJL, 537, L9

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

Schaller M. et al., 2015a, MNRAS, 451, 1247

Schaller M. et al., 2015b, MNRAS, 452, 343

Schaye J. et al., 2015, MNRAS, 446, 521

Sellwood J. A., 2008, ApJ, 679, 379

Shankar F., Weinberg D. H., Miralda-Escudé J., 2009, ApJ, 690, 20

Shankar F., Marulli F., Bernardi M., Mei S., Meert A., Vikram V., 2013, MNRAS, 428, 109

Spekkens K., Giovanelli R., Haynes M. P., 2005, AJ, 129, 2119

Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760

Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, 398, L21

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Suto D., Peirani S., Dubois Y., Kitayama T., Nishimichi T., Sasaki S., Suto Y., 2017, PASJ, 69, 14

Swaters R. A., Madore B. F., van den Bosch F. C., Balcells M., 2003, ApJ, 583, 732

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L., 2011, MNRAS, 414, 195

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Thomas J., Saglia R. P., Bender R., Erwin P., Fabricius M., 2014, ApJ, 782, 39

Tissera P. B., White S. D. M., Pedrosa S., Scannapieco C., 2010, MNRAS, 406, 922

Tollet E. et al., 2016, MNRAS, 456, 3542

Tonini C., Lapi A., Salucci P., 2006, ApJ, 649, 591

Trujillo I., Erwin P., Asensio Ramos A., Graham A. W., 2004, AJ, 127, 1917

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104

Umetsu K., Takada M., Broadhurst T., 2007, Mod. Phys. Lett. A, 22, 2099

Vogelsberger M. et al., 2014a, MNRAS, 444, 1518

Vogelsberger M., Zavala J., Simpson C., Jenkins A., 2014b, MNRAS, 444, 3684

Volonteri M., Dubois Y., Pichon C., Devriendt J., 2016, MNRAS, 460, 2979

Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20

Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt R. C., Jr, Thornley M. D., Leroy A., 2008, AJ, 136, 2563

Weinberg M. D., Katz N., 2002, ApJ, 580, 627

Weinberger R. et al., 2017, MNRAS, 465, 3291

Welker C., Dubois Y., Devriendt J., Pichon C., Kaviraj S., Peirani S., 2017, MNRAS, 465, 1241

Yoshida N., Springel V., White S. D. M., Tormen G., 2000, ApJ, 535, L103

Zitrin A. et al., 2015, ApJ, 801, 44

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