A look to the inside of haloes: a characterisation of the halo shape as a function of overdensity in the Planck cosmology
Abstract
In this paper we study the triaxial properties of dark matter haloes of a wide range of masses extracted from a set of cosmological body simulations. We measure the shape at different distances from the halo centre (characterised by different overdensity thresholds), both in three and in two dimensions. We discuss how halo triaxiality increases with mass, redshift and distance from the halo centre. We also examine how the orientation of the different ellipsoids are aligned with each other and what is the gradient in internal shapes for halos with different virial configurations. Our findings highlight that the internal part of the halo retains memory of the violent formation process keeping the major axis oriented toward the preferential direction of the infalling material while the outer part becomes rounder due to continuous isotropic merging events. This effect is clearly evident in high mass haloes  which formed more recently  while it is more blurred in low mass haloes. We present simple distributions that may be used as priors for various mass reconstruction algorithms, operating in different wavelengths, in order to recover a more complex and realistic dark matter distribution of isolated and relaxed systems.
keywords:
galaxies: halos  cosmology: theory  dark matter  methods: numerical1 Introduction
Different wide field surveys, observing at various wavelengths, are revealing that most of the matter density content of our Universe is in form of collisionless particles (Amara et al., 2012; Guzzo et al., 2014; Covone et al., 2014). These particles do not emit radiation and interact only gravitationally with the surrounding density field: they are generally termed Dark Matter. Following the standard scenario of structure formation, dark matter drives the structure evolution processes: systems up to protogalactic scales form as consequence of gravitational collapse and then merge together, along the cosmic time, forming the more massive ones (White & Rees, 1978; Tormen, 1998; Springel et al., 2001b; Tormen et al., 2004). Galaxy clusters sit at the top of this hierarchical pyramid being the most massive and late forming virialized structures of our Universe (Frenk et al., 1990; Borgani & Kravtsov, 2011; Angulo et al., 2012).
Various studies of time evolving isolated perturbations seeding in the dark matter density field have lead many scientist to the development of the spherical collapse model (White & Silk, 1979; Press & Schechter, 1974; Pace et al., 2010). A density perturbation grows with the expansion of the Universe in concentric shells and, if it is dense enough, will pull away from the background expansion, and will collapse after reaching a maximum size characterised by null kinetic energy. The collapse happens when the perturbation exceeds the predicted critical value by the spherical collapse model forming a so called dark matter halo (Bond et al., 1991; Eke et al., 1996; Bryan & Norman, 1998). It is interesting to notice that density perturbations from which haloes form are not independent with each other nor isolated, but during the expansion and collapse perturbations are typically pulled, stretched and sheared by the surrounding density field (Doroshkevich, 1970; Despali et al., 2013). All these also translate in a mass dependence of the collapse threshold (Eisenstein & Loeb, 1995; Sheth & Tormen, 1999; Sheth et al., 2001; Sheth & Tormen, 2002) and in the formation of haloes that are generally triaxial, more in particular prolate (Jing & Suto, 2002; Despali et al., 2014; Bonamigo et al., 2015; Vega et al., 2016). Thus, the standard spherical modelling of the dark matter, stars, intergalactic and the intracluster medium is only a rough approximation and in particular both theory and observations agree on the general picture that haloes in which galaxies and clusters live are very well approximated by a triaxial ellipsoid (Morandi et al., 2011, 2012; Sereno & Zitrin, 2012; Limousin et al., 2013; Groener & Goldberg, 2014). The study of the asphericity of galaxy clusters is growing, in light of the analyses performed on different numerical simulations during the last years (Angrick & Bartelmann, 2010; Rossi et al., 2011a, b; Despali et al., 2014; Vega et al., 2016).
In this paper we aim to discuss the dependence of the halo shape on the distance from the centre. Other previous works (Allgood et al., 2006; Bailin & Steinmetz, 2005) have measured shapes at different fractions of the virial radius. For reasons that are linked to the different halo mass definitions, , or (Despali et al., 2016) and that we will better discuss later in the text, in this work we choose to present our results in term of different overdensity thresholds: we define halos as triaxial regions enclosing a desired multiple of the critical density of the Universe and, for each halo identified using , we measure the triaxial shape at other four overdensity thresholds, multiples of the critical density of the Universe (200, 500, 1000 and 2000). The choice of these values is also motivated by the fact that, typically, the Xray community adopts a mass definition that is associated to the region enclosing 500 times the critical density, while weak lensing and dynamical analyses usually prefer 200 times the critical density. On the other side, strong lensing researchers make use of the very central region of clusters or galaxies, where critical lines emerge; in this case we can refer to the region enclosing at least 1000 times the critical density of the Universe (Broadhurst et al., 2005b; Coe et al., 2012; Zitrin et al., 2011). In light of the various observational analyses, it is important to underline that a study of the degree of alignment of the mass density distribution at different distances from the centre can help us to understand dark matter and baryonic properties in which the various physical processes are taking place. Hopkins et al. (2005) have studied the shape properties of haloes extracted from a lightcone up to redshift constructed from a largescale highresolution body simulation. They discuss that the mean halo ellipticity increases with redshift as and with the cluster mass, as also found by Despali et al. (2014); Bonamigo et al. (2015). On the other hand, for Hopkins et al. (2005) the cluster ellipticity decreases with radius in disagreement with other results (Bailin & Steinmetz, 2005; Hayashi et al., 2007; VeraCiro et al., 2011; Velliscig et al., 2015). In particular VeraCiro et al. (2011), studying the body haloes from the Aquarius simulation (Springel et al., 2008), found that the evolution in halo shape correlates well with the distribution of the infalling material: prolate configurations arise when haloes are fed through narrow filaments whereas triaxial/oblate configurations result as the accretion turns more isotropic at later times. The geometrical properties of haloes at different epochs are not lost: haloes retain memory of their structure at earlier times. This memory is imprinted in their presentday shape trends with radius, which change from typically prolate in the inner (earlier collapsed) regions to a triaxial in the outskirts – corresponding to the shells that have collapsed last and are now at the virial radius.
In this work we will present a study of the shape properties of haloes extracted from DM only simulations. We underline that our results do not account for the presence of baryons – mainly influencing the most internal shells – and could eventually be adapted to their presence using precalibrated analytical recipes.
The paper is organised as follows: Section 2 describes the numerical simulations and the halo catalogues; in Section 3 we discuss how we selected relaxed and regular halos; our results are presented in Sections 4 and 5, which show respectively the distributions derived in three or two dimensions; Section 6 is dedicated to summarise our results and draw our conclusions.
2 The Numerical Simulations
2.1 Le SBARBINE simulations
Le SBARBINE simulations are six cosmological simulations which have been run in Padova using the publicly available code GADGET2 (Springel, 2005); these are part of a series of new simulations which have been presented in a previous work (Despali et al., 2016). The adopted cosmology follows the recent Planck results (Planck Collaboration et al., 2014), in particular we have set: , , and . The initial power spectrum was generated with the code CAMB (Lewis et al., 2000) and initial conditions were produced perturbing a glass distribution with NGenIC (http://www.mpagarching.mpg.de/gadget). They all follow collisionless particles in a periodic box of variable length (we refer the reader to Table 1 for more details).
name  box [Mpc ]  []  soft [kpc ]  

Ada  62.5  124  1.5  2264847  36561  
Bice  125  99  3  2750411  44883  
Cloe  250  99  6  3300880  54467  
Dora  500  99  12  3997898  58237  
Emma  1000  99  24  4739379  38636  
Flora  2000  99  48  5046663  5298 
2.2 Halo catalogues
At each stored snapshot, we identified the dark matter haloes using the Ellipsoidal Overdensity algorithm, as described in Despali et al. (2013, 2014); Despali et al. (2016) and Bonamigo et al. (2015). This algorithm identifies ellipsoidal haloes in numerical simulations: it works similarly to the more common Spherical Overdensity criterion (Lacey & Cole, 1994; Tormen et al., 2004; Planelles & Quilis, 2010; Knebe et al., 2011, 2013), with the difference that the halo shape is refined using an iterative procedure to find the best triaxial ellipsoid that follows the mass density distribution  instead of forcing a spherical shape. For example, at the present time we used the virial overdensity as a density threshold for the main halo catalogues (Eke et al., 1996; Bryan & Norman, 1998). Then, we identified the haloes at other four overdensity thresholds, corresponding to 200, 500, 1000 and 2000 (as in Despali et al. (2016)); each run has been made independently, so that the resulting shape and direction of each shell is not influenced by the virial value.
We calculate halo shapes using eigenvalues of the mass tensor, defined as:
(1) 
where is the position vector of the the particle and and are the tensor indexes. By diagonalising we obtain the eigenvalues and eigenvector: the axes of the ellipsoid () are then defined as the square roots of the eigenvalues.
3 Halo selection
Since the purpose of this work is to provide reliable prediction of the halo shapes as a function of radius – and overdensities, we decided to restrict our halo catalogue in order to exclude irregular, merging or highly unrelaxed haloes. First of all we applied one of the common criteria to define relaxed haloes (Neto et al., 2007; Macciò et al., 2008)  method 1: we calculated the distance between the position of the minimum of potential and the centre of mass of the halo; then, we maintain only systems in which this difference is less than of the corresponding halo virial radius. As seen in Bonamigo et al. (2015), this criterion is able to exclude most of the irregular haloes – meaning those that cannot be reliably described by one single triaxial ellipsoid – and their fraction increases with the mass, due to the fact that high mass haloes, forming later (Giocoli et al., 2007; Zhao et al., 2009; Giocoli et al., 2012b), are still in a merging phase. As a second criterion  method 2, we calculated the total energy of haloes – as a sum of the kinetic and potential energies of the constituting particles – and discarded those with positive energy, getting rid of some other irregular systems. In Figure 1 (top and bottom panels refer to and , respectively) we illustrate the effect of these two selection criteria on the halo catalogue showing the percentage of irregular haloes detected (and so excluded) by each method and by their combination (requiring that at least one of the two methods excludes the halo).
Nevertheless, after this first selection, we noticed that some irregular haloes were still present in our catalogues, as for example the two projections of the halo displayed on top panels of Figure 2. This system was chosen randomly between those exhibiting extremely low virial axial ratios (), who were not excluded by the first selections: it is clearly an unrelaxed halo, being composed by multiple mass clumps and in a merging phase. It survived through the previous selection because (even if it may not be clear from the figure due to the projection effect and the colour combination) its main clump is still more massive than the other, keeping the centre of mass near the minimum of potential. From the figure, we can notice how the shape of the ellipsoid (shown in blue) is very different from the overall virial shape (black dots): while , the axial ratio in the inner ellipsoid is , in contrast with the common finding that inner parts are more elongated than the outer ones (see Section 4.1). For comparison, the two bottom panels, show two projections of a relaxed halo with similar mass: its axial ratios do not change dramatically as a function of radius and the ellipsoid enclosing is clearly larger and more massive than in the previous case.
In order to capture systems like the one displayed in the top panel of the previous figure, we suggest and use a third selection criterion (method 3), based on the discreet measurement of the density profile of the halo computed on elliptical shell of defined overdensities. In particular, we measure the mass fraction contained in each overdensity ellipsoid, with respect to the total virial mass. Figure 3 shows the mass enclosed by each overdensity threshold, as a function of the overdensity with respect to the critical one. The points show the median values in each mass bin: the mass decreases with overdensity, with a mass dependent slope because smaller haloes are more concentrated than the more massive ones (Giocoli et al., 2012; Meneghetti et al., 2014). In those haloes which are in a merging phase, the central clump can be expected to be less massive than in relaxed haloes, since a significant part of the mass still resides in the infalling clumps. Thus, this median density profile can be used as a selection criterion, in particular to characterise a relaxed sample in addition to the method 1 and method 2: we exclude also all haloes for which the mass fraction always lies in the lower quartile ( of the distribution) for all the four discrete overdensity ellipsoid. We remind the reader that this last selection criterion is very analogous to the substructure mass fraction method adopted by Neto et al. (2007) and to the residuals from a NFW fit used in Macciò et al. (2007) to characterise relaxed and unrelaxed haloes. The fraction of haloes excluded only by method 3 is shown by the magenta histogram in Figure 1: they are approximately of the whole halo sample both at and . Adding method 3 as an exclusion criterion, we are able to discard some irregular haloes that are not identified by the first two methods: the final cut in the halo catalogue is shown by the black histogram. The relaxed halo in the bottom panels of Figure 2 satisfy all our selection criteria, thus proving to be relaxed and regular – in all considered overdensities, while the top one is successfully excluded from the whole sample by method 3.
4 3D shape as a function of overdensity
In this section we analyse how the threedimensional shape of relaxed haloes changes as a function of the overdensity, and characterise the variation as a function of mass and redshift. This work is complementary to what has been done in Despali et al. (2014) and Despali et al. (2016); by presenting our results in terms of overdensity instead of radius, we want to provide useful distributions for observational studies in general and for strong lensing parametric mass modelling in particular, as these distributions can be inserted as priors in algorithms for the representation of the position of lensed multiple image systems ( for example Lenstool (Jullo et al., 2007) ).
We start by showing how the overall distribution of shapes depends on overdensity (for different masses and redshifts) and then we explore the conditional distribution of shapes, binning in the virial axial ratio: in this way we can give realistic prediction of the change in shape within individual halos going from the outside to the inside. Finally, we address the misalignment between the ellipsoids and show how it depends on mass and an the virial shape properties.
A general colour/style code is used through the paper: different overdensities are represented by various colours (from black for the virial case to red for the innermost one, going through green, blue and magenta) and results are shown by solid lines, while at we chose dashed lines – we mention that for better display our data and results we do not show the case that lays in the middle between and .
4.1 General distributions
First, we confirm that halo shapes are not selfsimilar as a function of distance from the centre (Allgood et al., 2006; VeraCiro et al., 2011; Jing & Suto, 2002; Bailin & Steinmetz, 2005). Figure 4 shows how the axial ratios varies as a function of the virial mass for our five overdensity thresholds. We remind the reader that the three dimensional ellipticity is defined as
(2) 
(with ) and it is equal to zero for a spherical system. Apart from the well known dependence on mass (Allgood et al., 2006; Despali et al., 2014; Bonamigo et al., 2015), we notice how the dependence on overdensity is almost the same and regular for all relaxed masses, with inner ellipsoids being more triaxial than the outermost virial one, shown in black. The same behaviour can be observed at (dashed lines): at this time, for a given mass bin, haloes are generally more triaxial, but the inner ordering is unaltered. Note that the black and green dashed lines (virial and ) almost coincide, since at this redshift the two overdensities are very close to each other (Despali et al., 2016). The analogous distributions for unrelaxed haloes at are shown in Appendix A, proving why unrelaxed haloes cannot be easily described by simple relations, and a regular trend is absent. This highlight the fact that the morphological properties of galaxies and clusters at different radii may be used to infer the state of relaxation (Donahue et al., 2016) and that our relations for relaxed haloes do not hold for some of the recently observed clusters (in particular for five out of six Frontier Fields clusters), which appear unrelaxed and present multiple components.
In Figure 5 we present the shape distributions in more details: we consider five mass bins, centred on masses from to and outline the cumulative distributions of the shapes – axial ratios and ellipticities – for different overdensities. From the figure we notice that all trends are very regular both in mass and overdensities. Comparing the distributions at and we notice that at higher redshifts haloes of same mass are more triaxial, once the mass is fixed – the case would lay in the middle between and not shown here avoiding to overcrowd the panels. As described in Despali et al. (2014) the redshift dependence can be removed comparing haloes possessing the same peakheight .
(0.732,0.878,0.052)  (0.720,0.882,0.054)  (0.693,0.874,0.060)  (0.669,0.863,0.066)  (0.646,0.844,0.071)  
(0.687,0.851,0.061)  (0.670,0.851,0.066)  (0.637,0.832,0.074)  (0.611,0.810,0.081)  (0,585,0.782,0.087)  
(0.629,0.803,0.076)  (0,612,0.798,0.080)  (0.577,0.767,0.090)  (0.549,0.735,0.099)  (0.506,0.701,0.106)  
(0.546,0.710,0.100)  (0.537,0.713,0.102)  (0.506,0.678,0.113)  (0.481,0.646,0.121)  (0.465,0.620,0.128)  
(0.475,0.624,0.124)  (0.471,0.631,0.126)  (0.449,0.610,0.133)  (0.430,0.587,0.141)  (0.420,0.566,0.145) 
4.2 Conditional distributions
As already mentioned, strong lensing measurements focus on the very central region of galaxy cluster and so are able to probe the shape of the inner highdensity shells (Meneghetti et al., 2016). At the same time, theoretical prediction of mass and shape at the virial (or ) radius are often used in their analysis and the shape in the central parts is assumed to be selfsimilar to the virial one. Since we know that this is not the case and that the difference may not be the same for all masses, we think that it is important to link the shapes at various distances with the virial measurement. For this reason, in this Section we present conditional distributions, binning our data in virial shape.
Figure 6 shows the conditional distribution of the axial ratios, for different mass bins and redshifts and , using solid and dashed line styles, respectively. In order to characterise more precisely the shape variation inside individual haloes, we bin the axial ratio distributions using the computed values at the virial overdensity: each colour shows the median distribution of the axial ratios for haloes in a certain bin of virial axial ratio. Each bin in or has a width of . Thus, the virial axial ratios of haloes represented by the yellow lines are contained in the interval [, ], the green ones in [,] and so on and so forth. The lowest axial ratios [,] are represented in purple colour. Objects with even lower axial ratios are excluded by our selection criteria, since extreme elongations often coincide with haloes in a merging phase, as already discussed above in the text. A new feature, that was not visible in the previous figures, emerges from the conditional distributions: the variation of shape with overdensity – or radius – depends on their outer shape. While triaxial haloes with axial ratios are selfsimilar in their inner parts, the more spherical ones present a greater shape variation, becoming considerably more triaxial inside. This effect is present for all masses, even if it may be caused by different phenomena: in general, it has been argued that the outer parts of haloes become rounder due to the interactions with the surrounding density field, which take place after their formation, while the inner parts maintain the original triaxiality due to the collapse process. This is true in particular for low mass haloes, which formed earlier and are more influenced by the surrounding tidal field or by encounters with more massive structures. For highmass haloes, apart from the few with high formation redshifts, the physical explanation of this dependence may be different: it is possible that, while matter is still in an accretion phase from many directions onto a halo  as they live at the intersection of filaments, the centre already collapsed in a well defined triaxial object.
We underline that these distributions are also useful to generate mock mass density distributions of dark matter in galaxies and clusters producing for example more realistic lens models (Giocoli et al., 2012a). Moreover, using the virial shape as a prior in the analysis of strong lensing clusters and assuming that the halo shape is selfsimilar at all radii can introduce a bias in the calculation, as this is true only for a subset of objects and the innermost part of the halo can be much more elongated than the outside (Meneghetti et al., 2016). Looking at Figure 6, we also notice that the axial ratio of the central parts of the haloes tend to converge to similar values in all cases: this means that retrieving the true virial shape from observations, using only the inner shape as an information, can be risky. For this reason, the shape informations must be combined with those on mass, concentration and other halo properties in order to derive meaningful estimates.
4.3 Misalignment of the different ellipsoids
It has been previously shown that shapes measured at different radii within the same halo are not perfectly aligned with each other. Having independently measured the triaxial shapes of four inner ellipsoids, we can easily compare their relative orientations – with respect to a predefined direction of the three dimensional ellipsoid – in a way to better understand how on average the different ellipsoids are misaligned with each other. We remind the reader that our shape measurements include all the components of haloes and do not discriminate between the main smooth component and the substructures, as has been done in other previous works (as for example VeraCiro et al., 2011). In Figure 7 we show the median misalignment angle between the four inner ellipsoids with respect to the virial one – that we view as reference; we considered the misalignment between the two longest () and the two shortest axes () of the 3D mass ellipsoid – this will give us also an idea of the deformation of the triaxial mass ellipsoid. The measurements are divided in five mass bins, represented by different point types. From the figure we notice that while the median misalignment at is around 10 degrees only, it becomes larger when going toward the halo centre. We stress that there is a considerable dispersion in the data, which increases for low mass haloes: to give an idea about this, the dashed (dotted) lines show the 25% and 75% quartiles associated with the highest (lowest) mass bin  thus (). The median misalignment appears to depend on mass: in particular, for cluster size haloes, which formed very recently – or are in their formation phase, the various ellipsoids are aligned with each other within 10 degrees  probably to the direction of compression of the gravitational collapse, while low mass haloes – that formed typically at higher redshifts (Lacey & Cole, 1993; Giocoli et al., 2007) – present greater variations, again due to the interactions with the surrounding field and their evolution after the formation time (more evident is the halo to halo variation) – also being less gravitationally strong they tend to be more influenced by the surrounding density field. As shown in Despali et al. (2014), the ellipticity anticorrelates with the formation redshift and so more elongated haloes formed more recently. Also, this mass dependence may hide a geometrical shape dependence: low mass halos are rounder and thus the axes direction at the virial radius may be less defined than in a very elongated system. This seems to be supported also by the results of Figure 8, where we display the conditional distributions of the misalignment of the different ellipsoids enclosing various overdensities. The data are divided according to the virial axial ratio , as in Figure 6. From this figure it appears more clear that for very triaxial haloes all the ellipsoids are well aligned, probably due to the phase of collapse or a recent formation, while this is not true for rounder haloes. This figure shows all the masses together, as we noticed that the conditional distribution are very similar for different mass bins, reinforcing our framework and our interpretation.
Together with the results on conditional axial ratios presented at the beginning of this section, these distribution can be used to produce selfconsistent mock mass density distribution of realistic triaxial and perturbed haloes.
5 Projected 2D shape as a function of overdensity
After the discussion presented in the previous section about three dimensional shapes, we proceed analysing the shapes in two dimensions (2D), which can be more directly related with observed quantities projected on the plane of the sky and that have not been modelled in previous works. We present the 2D results with figures similar to the 3D ones discussed in the previous sections: cumulative distributions for different overdensities and mass bins and conditional distribution of the axial ratios.
We project each halo along three random directions in particular along the three axes of the coordinate system of the simulation (,,) – considering each projection a random measure of the 2D shape of the halo ellipsoid. Since we already have a relatively large number of haloes we do not consider necessary to project each object along different possible random line of sights. We then look at the distributions of halo shapes and orientations for different overdensities and masses, as done and discussed for the 3D case. In 2D we calculate the ellipticity as
(3) 
As a general result, 2Ddistributions maintain the same properties and ordering of the 3D ones, but they become shallower due to projection effects. Even the extreme cases, such as very elongated shapes, are blurred by being projected in random directions. This general effect can be seen in Figure 9 where the contours show the point density of the relationship between the axial ratios (left) and ellipticities (right) as measured in 2D versus the 3D ones. The black points show the median of the distribution at a fixed 3D value.
In Figure 10 we present the cumulative probability distribution of axial ratios (top panel) and ellipticities (bottom panel) for haloes of different masses and at redshift (solid curves) and (dashed curves). The ordering due to the different overdensity definitions is not altered when one looks at two dimensional quantities (Figure 11), even if the lowaxial ratios tail is reduced in comparison with the three dimensional quantities, as presented in Figure 5.
Same considerations also hold for the conditional distributions as presented in Figure 12, where we show the 2D axis ratios at various overdensities binning the haloes in term of the 2D axis ratios at the virial definition, and in Figure 13 which focuses on the misalignment.
6 Summary and Conclusions
The aim of this work is to provide simple and clear estimates of how the shape of relaxed haloes changes as a function of overdensity, mass and redshift. We looked at five different ellipsoids, enclosing various overdensities, and measured their shape in three and two dimensions. The main results of our work can be summarised as follows:

conditional distributions: the rate at which the shape varies through the halo depends on that at the virial radius; very triaxial haloes show a similar shape at all the overdensity ellipsoids and they are quite well aligned with each other, while for rounder haloes the inner ellipsoids are, proportionally, both more misaligned and more triaxial than the virial one.

2D projections: we calculated projected shapes by taking three different projections for each halo; the conclusions coming from projected quantities are similar to the 3D ones, even though the differences between halo ellipticities and orientations are shallower due to projection effects.
Our findings are consistent with the standard picture of structure formation, in which the central part of haloes may maintain its original triaxiality longer than the outskirts which are subjected to stronger interactions with the surrounding field; also, haloes formed recently will be still aligned with the direction of the last merger or of the filament along which matter accreted onto the halo, and so their whole shape will probably be well aligned. The distributions presented in this work may be used as priors for mass reconstruction algorithms working in different wavelengths, in order to recover a more realistic triaxial matter distribution of galaxies and clusters. In this framework, it is important to keep in mind that haloes are not perfectly selfsimilar and that to reconstruct the virial properties of the cluster dark matter halo from the strong lensing a multiparametric approach is needed.
Acknowledgments
CG thanks CNES for financial support. ML acknowledges the Centre National de la Recherche Scientifique (CNRS) for its support. This work was performed using facilities offered by CeSAM (Centre de donneS Astrophysique de Marseille (http://lam.oamp.fr/cesam/). This work was granted access to the HPC resources of AixMarseille Universite financed by the project Equip@Meso (ANR10EQPX2901) of the pro gram Investissements d’Avenir supervised by the Agence Nationale pour la Recherche (ANR). This work was car ried out with support of the OCEVU Labex (ANR 11 LABX0060) and the A*MIDEX pro ject (ANR11IDEX 000102) funded by the Investissements d’Avenir French government program managed by the ANR. We also thank the anonymous referee for her/his useful comments that helped to improve the presentation of our results.
Appendix A Unrelaxed haloes
In the analysis preformed in our paper we have chosen to discard unrelaxed and irregular haloes from our sample, because they tend to introduce more scatter in the distributions, not easily to explain and model in the same way for all the objects. In the various cases, the analyses should take into account the specific features of each system and of the field surrounding it. Figure 14 shows the analogous of Figure 4 for unrelaxed haloes, showing clearly why they cannot be modelled together with the more relaxed ones; in particular, the ellipticity trend is completely reversed, leading to more triaxial shapes in the outer ellipsoids (instead of the inner ones), probably due to the presence of multiple components or infalling material along a preferential direction at the observing time.
References
 Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
 Amara et al. (2012) Amara A., Lilly S., Kovač K., Rhodes J., Massey R., Zamorani G., Carollo et al. 2012, MNRAS, 424, 553
 Angrick & Bartelmann (2010) Angrick C., Bartelmann M., 2010, A&A, 518, A38
 Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins a., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
 Bailin & Steinmetz (2005) Bailin J., Steinmetz M., 2005, ApJ, 627, 647
 Bonamigo et al. (2015) Bonamigo M., Despali G., Limousin M., Angulo R., Giocoli C., Soucail G., 2015, MNRAS, 449, 3171
 Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
 Borgani & Kravtsov (2011) Borgani S., Kravtsov A., 2011, Advanced Science Letters, 4, 204
 Broadhurst et al. (2005b) Broadhurst T., Benítez N., Coe D., Sharon K., Zekser K., White R., Ford H., et al. B., 2005b, ApJ, 621, 53
 Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
 Coe et al. (2012) Coe D., Umetsu K., Zitrin A., Donahue M., Medezinski E., Postman M., Carrasco M., Anguita T., et al. 2012, ApJ, 757, 22
 Covone et al. (2014) Covone G., Sereno M., Kilbinger M., Cardone V. F., 2014, ApJ, 784, L25
 Despali et al. (2016) Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486
 Despali et al. (2014) Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208
 Despali et al. (2013) Despali G., Tormen G., Sheth R. K., 2013, MNRAS, 431, 1143
 Donahue et al. (2016) Donahue M., Ettori S., Rasia E., Sayers J., Zitrin A., Meneghetti M., Voit G. M., Golwala S., Czakon N., Yepes G., Baldi A., Koekemoer A., Postman M., 2016, ApJ, 819, 36
 Doroshkevich (1970) Doroshkevich A. G., 1970, Astrophysics, 6, 320
 Eisenstein & Loeb (1995) Eisenstein D. J., Loeb A., 1995, ApJ, 439, 520
 Eke et al. (1996) Eke V. R., Cole S., Frenk C. S., 1996, MNRAS, 282, 263
 Frenk et al. (1990) Frenk C. S., White S. D. M., Efstathiou G., Davis M., 1990, ApJ, 351, 10
 Giocoli et al. (2012a) Giocoli C., Meneghetti M., Bartelmann M., Moscardini L., Boldrin M., 2012a, MNRAS, 421, 3343
 Giocoli et al. (2012) Giocoli C., Meneghetti M., Ettori S., Moscardini L., 2012, MNRAS, 426, 1558
 Giocoli et al. (2007) Giocoli C., Moreno J., Sheth R. K., Tormen G., 2007, MNRAS, 376, 977
 Giocoli et al. (2012b) Giocoli C., Tormen G., Sheth R. K., 2012b, MNRAS, 422, 185
 Groener & Goldberg (2014) Groener A. M., Goldberg D. M., 2014, ApJ, 795, 153
 Guzzo et al. (2014) Guzzo L., Scodeggio M., Garilli B., Granett B. R., Fritz A., Abbas U., Adami C., Arnouts et al. 2014, A&A, 566, A108
 Hayashi et al. (2007) Hayashi E., Navarro J. F., Springel V., 2007, MNRAS, 377, 50
 Hopkins et al. (2005) Hopkins P. F., Bahcall N. A., Bode P., 2005, ApJ, 618, 1
 Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
 Jullo et al. (2007) Jullo E., Kneib J.P., Limousin M., Elíasdóttir Á., Marshall P. J., Verdugo T., 2007, New Journal of Physics, 9, 447
 Knebe et al. (2011) Knebe A., Knollmann S. R., Muldrew S. I., Pearce F. R., AragonCalvo M. A., Ascasibar Y., Behroozi P. S., Ceverino D., et al. 2011, MNRAS, 415, 2293
 Knebe et al. (2013) Knebe A., Pearce F. R., Lux H., Ascasibar Y., Behroozi P., Casado J., Moran C. C., Diemand J., et al. 2013, MNRAS, 435, 1618
 Lacey & Cole (1993) Lacey C., Cole S., 1993, MNRAS, 262, 627
 Lacey & Cole (1994) Lacey C., Cole S., 1994, MNRAS, 271, 676
 Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
 Limousin et al. (2013) Limousin M., Morandi A., Sereno M., Meneghetti M., Ettori S., Bartelmann M., Verdugo T., 2013, Space Sci.Rev.
 Macciò et al. (2008) Macciò A. V., Dutton A. A., van den Bosch F. C., 2008, MNRAS, 391, 1940
 Macciò et al. (2007) Macciò A. V., Dutton A. A., van den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55
 Meneghetti et al. (2016) Meneghetti M., Natarajan P., Coe D., Contini E., De Lucia G., Giocoli C., Acebron A., Borgani S., et al. 2016, ArXiv eprints: 1606.04548
 Meneghetti et al. (2014) Meneghetti M., Rasia E., Vega J., Merten J., Postman M., Yepes G., Sembolini F., Donahue M., Ettori S., Umetsu K., Balestra I., et al. 2014, ApJ, 797, 34
 Morandi et al. (2012) Morandi A., Limousin M., Sayers J., Golwala S. R., Czakon N. G., Pierpaoli E., Jullo E., Richard J., Ameglio S., 2012, MNRAS, 425, 2069
 Morandi et al. (2011) Morandi A., Pedersen K., Limousin M., 2011, ApJ, 729, 37
 Neto et al. (2007) Neto A. F., Gao L., Bett P., Cole S., Navarro J. F., Frenk C. S., White S. D. M., Springel V., Jenkins A., 2007, MNRAS, 381, 1450
 Pace et al. (2010) Pace F., Waizmann J.C., Bartelmann M., 2010, MNRAS, 406, 1865
 Planck Collaboration et al. (2014) Planck Collaboration Ade P. A. R., Aghanim N., Alves M. I. R., ArmitageCaplan C., Arnaud M., Ashdown M., AtrioBarandela F., Aumont J., Aussel H., et al. 2014, A&A, 571, A1
 Planelles & Quilis (2010) Planelles S., Quilis V., 2010, A&A, 519, A94
 Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
 Rossi et al. (2011a) Rossi G., Sheth R. K., Tormen G., 2011a, MNRAS, 416, 248
 Rossi et al. (2011b) Rossi G., Sheth R. K., Tormen G., 2011b, MNRAS, 416, 248
 Sereno & Zitrin (2012) Sereno M., Zitrin A., 2012, MNRAS, 419, 3280
 Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
 Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
 Sheth & Tormen (2002) Sheth R. K., Tormen G., 2002, MNRAS, 329, 61
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Springel et al. (2008) Springel V., Wang J., Vogelsberger M., Ludlow A., Jenkins A., Helmi A., Navarro J. F., Frenk C. S., White S. D. M., 2008, MNRAS, 391, 1685
 Springel et al. (2001b) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001b, MNRAS, 328, 726
 Tormen (1998) Tormen G., 1998, MNRAS, 297, 648
 Tormen et al. (2004) Tormen G., Moscardini L., Yoshida N., 2004, MNRAS, 350, 1397
 Vega et al. (2016) Vega J., Yepes G., Gottlöber S., 2016, ArXiv eprints
 Velliscig et al. (2015) Velliscig M., Cacciato M., Schaye J., Hoekstra H., Bower R. G., Crain R. A., van Daalen M. P., Furlong M., McCarthy I. G., Schaller M., Theuns T., 2015, MNRAS, 454, 3328
 VeraCiro et al. (2011) VeraCiro C. A., Sales L. V., Helmi A., Frenk C. S., Navarro J. F., Springel V., Vogelsberger M., White S. D. M., 2011, MNRAS, 416, 1377
 White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
 White & Silk (1979) White S. D. M., Silk J., 1979, ApJ, 231, 1
 Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Bnörner G., 2009, ApJ, 707, 354
 Zitrin et al. (2011) Zitrin A., Broadhurst T., Barkana R., Rephaeli Y., Benítez N., 2011, MNRAS, 410, 1939