# Universality of dark matter haloes shape over six decades in mass: insights from the Millennium XXL and SBARBINE simulations

Mario Bonamigo, Giulia Despali, Marceau Limousin, Raul Angulo, Carlo Giocoli, Genevieve Soucail

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

## ► To cite this version:

Mario Bonamigo, Giulia Despali, Marceau Limousin, Raul Angulo, Carlo Giocoli, et al.. Universality of dark matter haloes shape over six decades in mass: insights from the Millennium XXL and SBARBINE simulations. Monthly Notices of the Royal Astronomical Society, 2015, 449 (3), pp.3171-3182. 10.1093/mnras/stv417. hal-01434391

HAL Id: hal-01434391

https://hal.science/hal-01434391

Submitted on 15 Sep 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Universality of dark matter haloes shape over six decades in mass: insights from the Millennium XXL and SBARBINE simulations

Mario Bonamigo,<fnref n="1" />* Giulia Despali,<fnref n="2" /> Marceau Limousin,<fnref n="1" /> Raul Angulo,<fnref n="3" /> Carlo Giocoli<fnref n="4" />,<fnref n="5" />,<fnref n="6" /> and Geneviève Soucail<fnref n="7" />,<fnref n="8" />

Accepted 2015 February 23. Received 2015 January 14; in original form 2014 September 29

## ABSTRACT

For the last 30 yr many observational and theoretical evidences have shown that galaxy clusters are not spherical objects, and that their shape is much better described by a triaxial geometry. With the advent of multiwavelength data of increasing quality, triaxial investigations of galaxy clusters is gathering a growing interest from the community, especially in the time of ‘precision cosmology’. In this work, we aim to provide the first statistically significant predictions in the unexplored mass range above $3 \times 10^{14} M_{\odot}h^{-1}$, using haloes from two redshift snapshots ($z = 0$ and $z = 1$) of the Millennium XXL simulation. The size of this cosmological dark matter-only simulation (4.1 Gpc) allows the formation of a statistically significant number of massive cluster scale haloes ($\approx 500$ with $M > 2 \times 10^{15} M_{\odot}h^{-1}$, and $780 000$ with $M > 10^{14} M_{\odot}h^{-1}$). Besides, we aim to extend this investigation to lower masses in order to look for universal predictions across nearly six orders of magnitude in mass, from $10^{10}$ to almost $10^{16} M_{\odot}h^{-1}$. For this purpose we use the SBARBINE simulations, allowing us to model haloes of masses starting from $\approx 10^{10} M_{\odot}h^{-1}$. We use an elliptical overdensity method to select haloes and compute the shapes of the unimodal ones (approximately 50 per cent), while we discard the more unrelaxed. The minor to major and intermediate to major axis ratio distributions are found to be well described by simple universal functional forms that do not depend on cosmology or redshift. Our results extend the findings of Jing & Suto to a higher precision and a wider range of mass. This ‘recipe’ is made available to the community in this paper and in a dedicated web page.

**Key words:** methods: numerical – galaxies: clusters: general – galaxies: haloes – cosmology: theory – dark matter.

## 1 INTRODUCTION

Spectroscopic galaxy redshift surveys and numerical N-body simulations have revealed a large-scale distribution of matter in the Universe featuring a complex network of interconnected filamentary galaxy associations. Vertices, i.e. intersections among the filaments, correspond to the very dense compact nodes within this cosmic web where one can find massive galaxy clusters.

These objects have been first assigned a spherical geometry, being the easiest way to characterize a shape in three dimensions; at the time this fitted the available data well enough. Nowadays, with the advent of multiwavelength data of increasing quality, there is a growing interest from the community to go beyond the spherical assumption, which is inaccurate and misleading. At first, clusters of galaxies have been characterized as spherical objects, a model that fitted well enough the limited data available at the time. Nowadays, with the advent of multiwavelength data of increasing quality, there is a growing interest from the community to go beyond the spherical assumption.

<footnote n="1">¹Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="2">²Dipartimento di Fisica e Astronomia, Università degli Studi di Padova, vicolo dell'Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote n="3">³Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan 1, Planta-2, E-44001 Teruel, Spain</footnote>

<footnote n="4">⁴Dipartimento di Fisica e Astronomia, Alma Mater Studiorum Università di Bologna, viale Berti Pichat, 6/2, I-40127 Bologna, Italy</footnote>

<footnote n="5">⁵INAF - Osservatorio Astronomico di Bologna, via Ranzani 1, I-40127 Bologna, Italy</footnote>

<footnote n="6">⁶INFN - Sezione di Bologna, viale Berti Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote n="7">⁷Université de Toulouse, UPS-OMP, Institut de Recherche en Astrophysique et Planétologie (IRAP), F-31400 Toulouse, France</footnote>

<footnote n="8">⁸CNRS, IRAP, 14 Avenue Edouard Belin, F-31400 Toulouse, France</footnote>

<footnote>* E-mail: mario.bonamigo@lam.fr</footnote>
Many observational evidences have been collected showing that clusters are not spherical objects; in particular from the non-circular projection of various probes; density of cluster galaxies (Carter & Metcalfe 1980; Binggeli 1982); X-Ray surface brightness (Fabricant, Rybicki & Gorenstein 1984; Buote & Canizares 1992, 1996; Kawahara 2010; Lau et al. 2012); Sunyaev Zel'dovich pressure maps (Sayers et al. 2011); strong gravitational lensing (Soucail et al. 1987) and weak gravitational lensing (Evans & Bridle 2009; Oguri et al. 2010, 2012).

Recently, an azimuthal variation of galaxy kinematics has been detected in a stacked sample of 1 743 galaxy clusters from the SDSS (Skielboe et al. 2012). The line-of-sight velocity dispersion of galaxies lying along the major axis of the central galaxy is found to be larger than those along the direction of the minor axis, further evidence supporting the asphericity of galaxy clusters.

On the numerical side, haloes forming in cosmological simulations have been found to be triaxial in shape, with a preference for prolateness over oblateness (Frenk et al. 1988; Dubinski & Carlberg 1991; Warren et al. 1992; Cole & Lacey 1996; Jing & Suto 2002; Bailin & Steinmetz 2005; Hopkins, Bahcall & Bode 2005; Kasun & Evrard 2005; Allgood et al. 2006; Paz et al. 2006; Bett et al. 2007; Muñoz-Cuartas et al. 2011; Gao et al. 2012; Schneider, Frenk & Cole 2012; Despali, Tormen & Sheth 2013). These simulations also predict an evolution of the shape with mass and redshift; low-mass haloes appear more spherical than high-mass haloes, essentially because high-mass haloes have formed later on (Despali, Giocoli & Tormen 2014).

Finally, it can be shown (Doroshkevich 1970) that triaxial collapse is a straightforward prediction of structure growth driven by self-gravity of Gaussian density fluctuations.

These evidences shows that the triaxial framework, though still an approximation, encapsulates halo shapes much more accurately than the spherical counterpart.

Besides, it has been shown that cluster properties (mass, concentration parameter, slope of the inner dark matter density profile, strong lensing cross-section) can differ significantly depending on the shape assumed in the analysis (see, e.g. the discussion in Limousin et al. 2013, regarding Abell 1689); see also Giocoli et al. (2012a,b). Even the galaxy correlation function can be affected by wrong assumptions on the triaxiality of haloes (van Daalen, Angulo & White 2012).

Since the mentioned properties constitute key ingredients of important cosmological tests, this suggests that triaxial modelling is the next milestone in the road map of 'precision cosmology' with galaxy clusters.

In this paper, we aim to characterize the shape of numerically simulated clusters, described within a triaxial framework. Apart from the three Euler angles, a triaxial geometry is characterized by three axes ($a < b < c$), hence two axial ratios: minor to major ($s = a/c$ in the following) and intermediate to major ($q = b/c$).

Shape of triaxial haloes have been investigated theoretically in a number of works which aim to characterize the dependence of shapes on mass, redshift, radius and so on. Most of the works agree on the fact that massive haloes are on average more elongated than low-mass haloes (Jing & Suto 2002; Allgood et al. 2006; Muñoz-Cuartas et al. 2011; Despali et al. 2013, 2014), since they form at later times and thus still retain memory of their original shape; which is influenced by the direction of the surrounding filaments or of the last major merger. Moreover, shapes depend also on redshift with haloes of all masses having on average smaller axial ratios at higher $z$; even though, the rank in mass is maintained at all times (Muñoz-Cuartas et al. 2011; Despali et al. 2014). Other works

have investigated halo shapes as a function of radius, measuring axial ratios of shells at different distances from the centre and the alignment between shells (Warren et al. 1992; Jing & Suto 2002; Bailin & Steinmetz 2005; Allgood et al. 2006; Schneider et al. 2012); haloes are more elongated in the central regions, while the outskirts are more rounded, probably due to interactions with the surrounding environment. Obviously the available number of haloes increased in parallel with computational resources; the analysis of Jing & Suto (2002) was based on simulations with 512³ particles in a 100 Mpc $h^{-1}$ box, which contained hardly any halo above $10^{14} M_{\odot} h^{-1}$ and some higher resolution runs which provided only 12 haloes with more than $10^6$ particles. On the other hand more recent works, i.e. Schneider et al. (2012), have been able to analyse larger data sets like the Millennium I and II simulations (Springel 2005; Boylan-Kolchin et al. 2009). The mass range between $10^{12}$ and $10^{14} M_{\odot} h^{-1}$ has been widely explored in all these works, while only recently small haloes down to $10^{10} M_{\odot} h^{-1}$ (Muñoz-Cuartas et al. 2011; Schneider et al. 2012) and some massive haloes of $10^{15} M_{\odot} h^{-1}$ (Despali et al. 2014) have been included in this kind of analysis. So far, no statistically significant predictions are available above $3 \times 10^{14} h^{-1} M_{\odot}$ and we rely on extrapolations from lower mass haloes when it comes to predict the shapes of massive galaxy clusters. With about 300 billion particles and a box size of 3 Gpc $h^{-1}$, the Millennium XXL (MXXL) simulation (Angulo et al. 2012) fills the range of high masses and explore the properties of cluster size haloes.

In this work our aims are twofold:

(i) using cluster scale haloes ($M > 10^{14} M_{\odot} h^{-1}$) from the MXXL simulation, we aim to provide predictions for the shape of massive clusters.

(ii) Then, we extend the mass range by considering haloes from the SBARBINE simulations, applying similar methods in order to investigate the shapes of haloes and provide predictions over five decades in mass, from $\sim 3 \times 10^{10}$ to $\sim 4 \times 10^{15} M_{\odot} h^{-1}$.

This paper is organized as follows: in Section 2, we present the simulations and the methodology used to extract haloes and measure their shapes. In Section 3, we present our results for the massive cluster scale haloes, then in Section 4 we extend our analysis to a broader mass range. In Section 5, we compare our findings with previous works. We discuss our results and conclude in Section 6.

## 2 HALO CATALOGUE

We have derived the shape of galaxy clusters from the MXXL simulation (Angulo et al. 2012). To generalize our analysis to lower masses, we used a new set of simulations (Despali et al., in preparation), which extended the mass range to more than five orders of magnitudes. From both simulations, we have analysed haloes from two redshifts: $z=0$ and 1. The main features of the simulations are described in the following sections and summarized in Table 1.

### 2.1 MXXL simulation

With a box side of 3 Gpc $h^{-1}$ (4.1 Gpc), this simulation was especially tailored to study massive haloes which can be only found in very large volumes, because of their nature of extremely rare objects and due to the dampening of large fluctuation modes in smaller boxes. The $6720^3 \approx 3 \times 10^{11}$ dark matter particles have a mass of $6.174 \times 10^9 M_{\odot} h^{-1}$; the Plummer-equivalent softening length is $\epsilon = 13.7$ kpc. For reasons of consistency with the previous Millennium runs, the adopted $\Lambda$CDM cosmology is the *Wilkinson*
<table><caption>**Table 1.** Main features of the simulations used in this work. The last two columns report the total number of haloes with more than 1000 particles ($N_h$) and the corresponding fraction of 'regular' haloes ($N_{reg}$), at redshift $z = 0$.</caption><thead><tr><th></th><th>Box (Mpc h<sup>-1</sup>)</th><th>z<sub>i</sub></th><th>m<sub>p</sub>(M<sub>☉</sub> h<sup>-1</sup>)</th><th>Soft (kpc h<sup>-1</sup>)</th><th>N<sub>h</sub>(z = 0)</th><th>N<sub>reg</sub>(z = 0)</th></tr></thead><tbody><tr><td>Ada</td><td>62.5</td><td>129</td><td>1.94 × 10<sup>7</sup></td><td>1.5</td><td>39 445</td><td>28 005</td></tr><tr><td>Bice</td><td>125</td><td>99</td><td>1.55 × 10<sup>8</sup></td><td>3</td><td>49 100</td><td>32 107</td></tr><tr><td>Dora</td><td>500</td><td>99</td><td>9.92 × 10<sup>9</sup></td><td>12</td><td>66 300</td><td>33 970</td></tr><tr><td>Emma</td><td>1000</td><td>99</td><td>7.94 × 10<sup>10</sup></td><td>24</td><td>46 665</td><td>20 696</td></tr><tr><td>Flora</td><td>2000</td><td>99</td><td>6.35 × 10<sup>11</sup></td><td>48</td><td>7754</td><td>2997</td></tr><tr><td>MXXL</td><td>3000</td><td>63</td><td>6.17 × 10<sup>9</sup></td><td>18.8</td><td>937 755</td><td>568 477</td></tr></tbody></table>

*Microwave Anisotropy Probe* one: total matter density $\Omega_m = 0.25$, baryons density $\Omega_b = 0.045$, cosmological constant $\Omega_\Lambda = 0.75$, power spectrum normalization $\sigma_8 = 0.9$ and dimensionless Hubble parameter $h = 0.73$. The main properties of the simulation are summarized in Table 1.

Due to the huge number of haloes in the simulation (almost 68 millions at redshift 0), we restricted the analysis to only a random subsample; for each logarithmic mass bin of size 0.2 [mass inside a spherical overdensity (SO) of $200\Omega_{crit}$] we chose either $10^5$ random objects or all, for the higher masses where the number of haloes in the bin is lower. This cut happens at a logarithmic mass of about 14.4 and 14.0 $M_{\odot} h^{-1}$ for redshifts 0 and 1, respectively. We have then re-identified haloes at redshift $z = 1$ and 0 using an ellipsoidal halo finder, which will be described in Section 2.3.

## 2.2 LE SBARBINE simulations

With the purpose of comparing different data sets and extending the available mass range, we use (from Section 4 on) the results from five cosmological simulations which have been run in Padova using the publicly available code GADGET-2 (Springel 2005); these are part of a series of new simulations which will be presented in a subsequent work ('LE SBARBINE' simulations, Despali et al., in preparation). The adopted cosmology follows the recent Planck results (Planck Collaboration XVI 2014): $\Omega_m = 0.307$, $\Omega_\Lambda = 0.693$, $\sigma_8 = 0.829$ and $h = 0.677$. The initial power spectrum was generated with the code CAMB (Lewis, Challinor & Lasenby 2008) and initial conditions were produced perturbing a glass distribution with N-GenIC (http://www.mpa-garching.mpg.de/gadget). They all follow $1024^3$ particles in a periodic box of variable length. Table 1 shows some of the main characteristics of these simulations. Haloes were identified using a SO algorithm (Tormen, Moscardini & Yoshida 2004; Giocoli, Tormen & van den Bosch 2008) and then the best-fitting ellipsoid was found using an ellipsoidal overdensity method, already presented in Despali et al. (2013, 2014) and similar to the one used on the MXXL haloes and described in the next section; the two codes produce equivalent results. We selected only haloes with more than 1000 particles to ensure a good resolution and to have a good comparison with the haloes of the MXXL simulation.

## 2.3 Ellipsoidal halo finder

It is known that Friends Of Friends finders tend to connect together multiple virialized haloes via thin bridges of particles (Jing & Fang 1994); thus, to characterize halo shapes more precisely, we used a second halo finder that iteratively selects particles inside an ellipsoid

and then uses their mass distribution to compute the ellipsoid for the next step in the iteration. A more detailed description of the 'ellipsoidal halo finder' in general and of the iterative procedure can be found in Despali et al. (2013).

We start with a traditional SO algorithm which selects particles inside a sphere of given overdensity, namely the value from the spherical collapse model at $z = 0$: $\Delta_{vir} = 359.7$ times the background density (Eke, Cole & Frenk 1996), and centred in the particle with lowest potential (most bound particle). We then compute the mass tensor<fnref n="1" />

$$M_{\omega\beta} = \sum_{i=1}^{N_V} \frac{m_i r_{i,\alpha} r_{i,\beta}}{M_{\text{TOT}}} \quad (1)$$

of the particles inside the virial radius of the sphere of mass $M_{TOT}$, where $r_i$ is the distance of the $i$th particle, of mass $m_i$, from the most bound particle. The tensor's eigenvectors give the direction of the ellipsoid that approximate the mass distribution, while the square roots of the eigenvalues are proportional to the axes length ($c > b > a$).

Having derived the triaxial distribution of dark matter for the SO, we use it to select particles inside an ellipsoid. This technique has been already adopted in the literature (Allgood et al. 2006; Schneider et al. 2012; Despali et al. 2013); however, different authors use different criteria to define the ellipsoid. We select particles inside an ellipsoid, centred in the most bound particle, that encloses an overdensity equal to the virial one, as provided by the spherical collapse model $\Delta_{vir}$; we do not fix the mass, the volume or the major axis to be equal to the spherical values, as has been done in previous works. We then recompute the mass tensor with the new subset and we iterate this procedure until both the ratios of minor to major axis $s = a/c$ and intermediate to major axis $q = b/c$ converge within a 0.5 per cent of error. This method allows us to adopt a more general description while being still close to theoretical predictions and is the simplest possible extension of the SO, which actually becomes just the first step in our iteration. The mass difference between the spherical and the ellipsoidal identifications goes from 2 per cent at $10^{13} M_{\odot} h^{-1}$ to 5 per cent for very massive haloes of $5.5 \times 10^{15} M_{\odot} h^{-1}$. On the other hand the change in the measured shapes is about 30 per cent and cannot be ignored, as shown also by Despali et al. (2013). Ellipsoidal masses are in general larger than the spherical ones, since a triaxial shape is expected to follow the actual distribution of matter better than a sphere, and so follow the overdense regions, adding more particles. The difference in the following results between using the mass inside a sphere or an ellipsoid is negligible.

<footnote n="1">¹ Not to be confused with the inertia tensor (Bett et al. 2007).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Density distribution (colour scale) of dark matter particles inside a 10 Mpc h⁻¹ side cube centred in two different haloes and the respective computed ellipsoids (red) that approximate the mass distribution of the halo. The halo shown on the left-hand panel has a virial mass of 5.29 × 10¹⁴ M☉ h⁻¹, the one on the right has a mass of 6.90 × 10¹⁴ M☉ h⁻¹. These represent two families of objects: a ‘regular’ haloes (left) and a perturbed one (right), due to the large amount of substructures the latter has to be discarded, as it cannot be well described by a triaxial approximation.</figcaption></figure>

## 2.4 Halo selection

Fig. 2 shows the mass function of all haloes (upper panel). Data from redshift $z=0$ and $z=1$ are indicated by red squares and blue circles, respectively. As previously explained, we have analysed only a random sample of the entire halo catalogue of the MXXL simulation; this is causing the flattening at the mass bins which have more than 10⁵ objects in the entire box ($\approx 10^{14} M☉ h^{-1}$). To avoid any resolution effect, we have kept only haloes with at least 1000 particles within the ellipsoid (vertical dashed line) in both simulations.

For our analysis, we cleaned the halo catalogue from many unrelaxed systems. As an example, Fig. 1 shows the density distribution of dark matter of two haloes (colour scale) and, in red, the computed ellipsoid which encloses an overdensity of $\Delta_{vir}$. The object on the left has a virial mass of 5.29 × 10¹⁴ M☉ h⁻¹ and represent a relaxed halo; the mass of the one on the right is 6.90 × 10¹⁴ M☉ h⁻¹ and it is clearly multimodal. The ellipsoid seems to capture quite well the overall three-dimensional matter distribution of the relaxed halo; though it fails, as expected, with the perturbed one. This system is highly asymmetrical and lacks of a well-defined centre, therefore it cannot be described using a single triaxial model. Since triaxial multimass modelling is beyond the purpose of this paper we decided to not consider for our analysis multimodal haloes, like the one on the right-hand panel.

In general, defining relaxed and unrelaxed haloes is not a trivial task; there is no absolute definition and the limiting criteria depend on the model that must be tested. In our case, where we are interested in haloes shapes, we would like at least to be able to describe them with a single ellipsoid. In particular, we need to well define their centres as well as their symmetry, up to a certain degree.

A possible way to discriminate this kind of objects is to look at the offset between centre of mass and geometrical centre of the ellipsoid, which is one of the commonly used criteria for the selection of relaxed haloes. While the latter is associated with the minimum of the potential (most bound particle) of the most massive substructure (Springel et al. 2001), the other represents the centre of mass of the whole particle distribution. This means that, if a significant number of massive substructures is present and per-

turbs the whole halo distribution, there will be an offset between the centre of the ellipsoid and the centre of mass. We decided to select only haloes for which the offset is less than 5 per cent of their virial radius,

$$ \frac{|\bar{x}_{MBP} - \bar{x}_{cm}|}{R_{vir}} < 0.05. \qquad (2) $$

The lower panel of Fig. 2 shows the percentage of cleaned haloes as a function of the mass. As expected, the number of perturbed haloes increases with the mass, due to more massive haloes being assembled recently (Giocoli et al. 2007). At high redshift (blue circles), the percentage of ‘regular’ haloes was lower and more constant with mass, than at $z=0$ (red squares). For cluster masses the percentage is roughly 50 per cent.

Generally, ‘relaxed’ haloes are selected using both this and two other criteria: the amount of mass in substructures and the ratio of kinetic to potential energy as measurements of the dynamical state of a halo (Meneghetti et al. 2014). For this reason, we choose to call our cleaned sample of haloes ‘regular’ and not relaxed haloes. However, as can be seen in Neto et al. (2007), the selection in the centre offset is responsible for the majority of the rejected haloes; this means that our selection is still able to eliminate the most unrelaxed and irregular objects. Ludlow et al. (2012) used a similar selection ($N_{200} > 5000$ and spherically defined haloes) and found different results: the fraction of objects with an offset less than 5 per cent is 0.536, while, combining all the three relaxation criteria, the fraction of selected haloes is 0.285. In comparison with their work, our selection is still able to capture approximately 65 per cent of all perturbed haloes.

Thus, the choice of the criteria to distinguish between relaxed and unrelaxed haloes is still different in different works. Since we are interested in the overall shape of haloes, we decided to use only the centre offset as a selection criterion, since it is able to exclude very irregular haloes which could not be well fitted by an ellipsoid; we believe that adding the other two criteria would not change our results more than a few percent. Moreover, our choice is motivated also by the fact that we do not want to restrict our analysis to a very limited and regular sample, since our future plans include a
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Mass function and halo selection of the MXXL sample at redshift $z=0$ (red squares) and $z=1$ (blue circles) obtained with the ellipsoidal overdensity. The vertical dashed line indicates the mass of a halo with 1000 particles. Top panel: points show the mass function of the whole selected halo catalogue. The cut at low masses is clearly visible. Bottom panel: points show the percentage of regular haloes in each mass bin (i.e. objects with a centre offset smaller than 5 per cent of their virial radius).</figcaption></figure>

comparison with observational results, which are far from being homogeneous and regular.

Another used criterion is the rms of the fit adopting an NFW profile (Macciò et al. 2007; Muñoz-Cuartas et al. 2011) as reference. As we are only interested in the overall shape of a halo, this does not play an important role; it is also interesting to mention that the halo profile may vary for the simple NFW predictions (Navarro et al. 2004; Prada et al. 2012; Klypin et al. 2014). This can be seen, for example, in fig. 6 of Macciò et al. (2007), where the authors compare the dependence of shape on mass for different selection criteria. A selection only on centre offset (as done here) correspond to what the authors called GOOD and NOISY halo samples (offset less than 4 per cent) and differ only by the goodness of the NFW fit. The two curves in the plot are almost identical, while there is a noticeable difference when compared with the BAD and UGLY sample (offset larger than 4 per cent).

These methods for selecting regular haloes have also been applied to the five SBARBINE simulations, obtaining a catalogue equivalent to the MXXL haloes. The selected number haloes for both simulations at redshift $z=0$ is shown in Table 1.

# 3 TRIAXIAL SHAPES OF MASSIVE GALAXY CLUSTERS FROM MXXL

## 3.1 MXXL results

In this first analysis we are mostly interested in the clusters mass range, therefore we will use only a portion of the available MXXL

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Probability distribution functions – differential and cumulative – of $s = a/c$. The distributions for the entire haloes population is shown in grey (and with dashed lines), while the red (solid) ones refer to the cleaned population.</figcaption></figure>

<table><caption>**Table 2.** Number of haloes in each logarithmic mass bin(in log (M/M$_{⊙}$h)) and percentage of regular haloes for redshifts $z=0$ and 1.</caption><thead><tr><th rowspan="2">log (M)[M<sub>☉</sub> h<sup>-1</sup>]</th><th colspan="3">z = 0</th><th colspan="2">z = 1</th></tr><tr><th>N<sub>h</sub></th><th>N<sub>reg</sub>/N<sub>h</sub></th><th>N<sub>h</sub></th><th>N<sub>reg</sub>/N<sub>h</sub></th><th>N<sub>h</sub></th><th>N<sub>reg</sub>/N<sub>h</sub></th></tr></thead><tbody><tr><td>14.0–14.2</td><td>57 759</td><td>58.56 per cent</td><td>30 823</td><td>41.19 per cent</td><td></td><td></td></tr><tr><td>14.2–14.4</td><td>56 083</td><td>56.61 per cent</td><td>13 271</td><td>39.11 per cent</td><td></td><td></td></tr><tr><td>14.4–14.6</td><td>42 951</td><td>53.52 per cent</td><td>3914</td><td>38.24 per cent</td><td></td><td></td></tr><tr><td>14.6–14.8</td><td>20 715</td><td>50.60 per cent</td><td>919</td><td>39.39 per cent</td><td></td><td></td></tr><tr><td>14.8–15.0</td><td>7823</td><td>48.50 per cent</td><td>134</td><td>36.81 per cent</td><td></td><td></td></tr><tr><td>15.0–15.2</td><td>2305</td><td>46.46 per cent</td><td>6</td><td>19.35 per cent</td><td></td><td></td></tr><tr><td>15.2–15.4</td><td>523</td><td>45.84 per cent</td><td></td><td></td><td></td><td></td></tr><tr><td>15.4–15.6</td><td>84</td><td>46.15 per cent</td><td></td><td></td><td></td><td></td></tr></tbody></table>

data. By taking the ratio of minor to major axis $s = a/c$ we can measure the degree of triaxiality of a halo; the closer $s$ is to 0, the less spherical the object is. If we combine this information with the value of the intermediate to major axis ratio $q = b/c$, we can infer how much prolate or oblate the halo is. In Fig. 3, the distribution of $s$ is shown for the entire halo catalogue (dashed grey curves), and for the regular one (solid red curves). The filled histograms represent the differential distributions, while the curves are cumulative distributions of the two different samples. In the original population there is a noticeable bump at low $s$ which corresponds to highly aspherical objects; clearly this is the case of unrelaxed or merging clusters. As it can be seen in the red histogram, the selection criteria we adopted have helped to remove this unwanted feature, since modelling them is beyond the goal of this work. We have divided our sample in eight logarithmic mass bins, from $10^{14}$ to $3.98 \times 10^{15} M_{☉} h^{-1}$. Table 2 reports the total number of haloes $N_h$ and the percentage of regular ones $N_{rel}/N_h$ for each mass bin for both redshifts of the MXXL. As expected, the number of clusters at high redshift is lower and we do not have any halo in the highest mass bins. As noted before, the percentage of regular haloes is higher at low masses, which formed earlier and thus had more time to reach an equilibrium state.

It has already been established (Jing & Suto 2002; Allgood et al. 2006; Bett et al. 2007; Schneider et al. 2012) that the axis ratio $s$ depends on the mass of the halo, however this dependence has not been tested at the high masses available in large simulation boxes such as the MXXL. Fig. 4 shows the distributions of $s$ for different mass bins in our sample – only five mass bins of Table 2,
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Probability distribution function of $s = a/c$ binned in mass using a fixed bin of $0.2 M_{\odot} h^{-1}$ for both redshifts. For clarity, we show the results for only five of the mass bins reported in Table 2.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Conditional distributions $p(q|s)$, with $q = b/c$. Different colours represent the distributions for six bins in $s$; solid and dashed histogram shows data from redshift 0 and 1, respectively. The vertical dashed lines of corresponding colour show the median value of $s$ for each bin.</figcaption></figure>

to avoid an overcrowded plot; as halo mass increases, the median value of the axis ratio becomes smaller, that is, the halo is less spherical. This effect is barely visible at redshift $z = 1$. Moreover the dispersion in $s$ is larger in the lower bins. It is also noticeable that the distributions are not symmetric, particularly they are skewed to low values of the axis ratio. To fully describe the shape of haloes, we need also the conditional probability distribution function $p(q|s)$, which is the distribution of $q$ for a given value of $s$. Fig. 5 shows the conditional distributions obtained for six bins in $s$: solid histogram for $z = 0$ and dashed for $z = 1$. The two redshifts are almost indistinguishable, which hints at the universality of the conditional distribution that will be discussed later on. For any interval, the median value of $b/c$ is fairly close to the median of $a/c$ (dashed vertical lines); although still fully triaxial, haloes tend to be prolate rather than oblate. For example, in the case of a ‘disc-like’ object, all the distributions would have been prominently shifted to values close to unity, because, in this case, $b \approx c$ independently of the minor axis $a$.

## 3.2 Minor to major axis ratio distribution: functional form

We aim to obtain a functional form to describe the axial ratio distributions at different masses. Due to the low statistic, Jing & Suto (2002) were not able to fully resolve the shape of the distribution and therefore assumed a Gaussian distribution. On the other hand, Schneider et al. (2012) claimed to be able to fit all the masses with a single beta distribution, although, even after a rescaling of $s$, they mention some residual mass dependence. Thanks to the large statistic in the MXXL simulation we are able to reconstruct the distributions with greater detail, even at large masses. Moreover, we are only interested in clusters, so we do not need the same level of generalization of the previous authors (see Section 4 for broader analysis). These two conditions allow us to simplify the analysis and obtain a better fit of the axial ratio distributions.

As shown by various authors (Press & Schechter 1974; Bond et al. 1991; Lacey & Cole 1993) the mass function written as a function of peak height $v = \delta_c(z)/\sigma(M)$ does not depend on redshift nor on cosmology (see appendix A for the details on how to compute $v$). It can be understood as follows: $\delta_c(z)$ is the critical overdensity of the spherical collapse model (initial density required for a fluctuation to collapse at redshift $z$), it increases with $z$; $\sigma(M)$ is the variance in the initial density field smoothed on a scale of a uniform sphere of mass $M$ and is higher for small masses. Then, since at high-redshift haloes were less massive, the dependences on time of the two quantities compensate with each other. For example, $v(M_*, z) = 1$ at every redshift, and $v > 1$ always represent a halo with a mass larger than the typical haloes collapsing at that time, even though the exact value of $M_*$ changes with redshift.

Fig. 6 shows the logarithm of $s$ versus the logarithm of peak height ($\approx$mass) for the selected haloes. Medians of $\log(s)$ for the two redshifts are shown in red squares and blue circles; the redshift dependence seen in Fig. 4 has disappeared completely. As already shown by Despali et al. (2014, fig. 5), the universality of haloes properties seems to extend also to the shape when using $v$ instead of mass. The change of variable allows us to provide results that are independent of the redshift and valid for different cosmologies. This idea was already in the original Jing & Suto (2002) paper, as the mass was given in units of $M_*$, but the use of $v$ is more

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (logarithmic) Distribution of $s$ as function of peak height: the black boxes and whiskers represent the quartiles and 1.5 the quartiles range of the combined distributions. The horizontal error shows the different bins, while the green solid line is the linear fit to the medians. Red squares and blue circles are redshift 0 and 1 subsamples. As reference, the corresponding mass for MXXL cosmology at $z = 0$ is shown in the top axis.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7. Distribution of the scaled axial ratio $\tilde{s}$ for masses shown in Table 2. It can be easily seen that the distributions at all masses are well represented by an unique fitting function.</figcaption></figure>

general and gives a more direct connection to the theory of structure formation. As a result, we can safely treat the two data sets as a single population, shown by the box and whiskers plot for a given $v$ bin (horizontal error bars). This plot confirms the previously mentioned trend: more massive haloes (higher $v$) are more aspherical. We have looked for a linear relation between $v$ and axis ratio in the log-log space (green line in Fig. 6) and obtained an inclination $a = -0.257 \pm 0.01$ and an intercept $b = -0.219 \pm 0.005$. The intercept is the logarithm of the median axis ratio at $M_*$: $\tilde{s}(M_*) = 10^b = 0.604 \pm 0.006$, which however does not enter directly in the following relations. The relation translate to a rescaling similar to the one adopted by previous authors: $\log(s) = a \log(v) + b \Rightarrow \tilde{s} = 10^b = 10^{\log s - a \log(v)} = s v^{-a}$. Therefore

$$ \tilde{s} = s v^{0.257} \qquad (3) $$

as $v$ takes care of any time and cosmology dependence, this rescaling is valid also for different redshifts and cosmologies. As Fig. 7 shows, distributions of the rescaled axis ratios (coloured histograms) are nearly indistinguishable from each other, meaning that we have eliminated all the dependence on the mass, in contrast with the findings of Schneider et al. (2012). Moreover, we were not able to fit the histogram of $\tilde{s}$ with a beta distribution. As it can be seen in Fig. 7, the distributions are non-zero at values greater than $\tilde{s} = 1$; this does not mean that there are haloes with axis ratio greater than 1: $\tilde{s}$ is not a physical quantity, this effect is due to the rescaling. Nevertheless, one can argue that $\tilde{s}$ represents the physical axis ratio at $v = 1$ ($M = 5.8 \times 10^{12} M_{\odot} h^{-1}$); still, this rescaling has been obtained only for $M > 10^{14} M_{\odot} h^{-1}$, leaving the unscaled axis ratio well within the physically meaningful boundaries. We have chosen to fit the minor to major axis ratio using a lognormal distribution

$$ p(x, \mu, \sigma) = \frac{1}{x\sqrt{2\pi}\sigma} \exp\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right), \qquad (4) $$

which corresponds to the probability distribution function of a variable which is normally distributed in the logarithmic space. The parameters of the fitted function are the following:

$$ \begin{aligned} \mu &= -0.49 \\ \sigma &= 0.20; \end{aligned} \qquad (5) $$

they can be converted to more familiar quantities

$$ \begin{aligned} \text{median} &= e^{\mu} = 0.61, \\ \text{std} &= \sqrt{(e^{\sigma^2} - 1)e^{2\mu + \sigma^2}} = 0.13. \end{aligned} \qquad (6) $$

In this framework, for a simple analysis, one can just use the scaled median value $\tilde{s} = e^{\mu} = 0.61$ with asymmetric quartiles at 0.53 and 0.70; then use equation (3) to obtain the physical value. On the other hand, it is possible to use the fit to obtain the whole distribution for a given mass. For example, to use it as a prior distribution of the minor to major axis ratio, one draws a value $x$ from a normal (Gaussian) distribution with mean $\mu = -0.49$ and standard deviation $\sigma = 0.20$, the scaled axis ratio is then $e^x$ (or directly extract $\tilde{s}$ from a lognormal distribution); inverting the rescaling relation one can obtain the axis ratio of the halo at a given peak height, which can be subsequently converted in mass for a given cosmology at a given redshift.

### 3.3 Intermediate to major axis ratio distribution: functional form

Once we are able to describe $s$ as a function of mass we can look at the correlation between the two axial ratios. For construction, $q$ is always greater (or equal) than $s$; also it is always less than 1. These limits have the effect of distorting the distribution of intermediate to major axis ratio in a way that depends directly on $s$. To avoid this problem, we use the rescaled quantity $\tilde{q} = (q - s)/(1 - s)$ instead of the simple axial ratio (Schneider et al. 2012), eliminating the issues of a limited interval; the correlation between the rescaled second axial ratio and $s$ can be seen in the left-hand panel of Fig. 8, where medians (red error bars) and quartiles (box and whiskers plot) for different values of the first axis ratio are shown. We have divided $\tilde{q}$ in bins of different $s$ and extracted the distributions $p(\tilde{q}|s)$ (right-hand panel of Fig. 8). From both plots, it is quite evident that $\tilde{q}$ strongly depends on the first axial ratio, with higher values at higher $s$, which is in agreement with haloes that tend to be prolate. Moreover the scatter is larger at higher $s$, though this is mostly due to the rescaling which extends the allowed interval of $\tilde{q}$.

Because of the strong correlation between $\tilde{q}$ and $s$, we cannot just give $\tilde{q}$ as a function of mass to obtain the second axis ratio distribution for a given mass, we have to describe $p(\tilde{q}|s)$ and then get the first axis ratio from its distribution at different masses (as shown in Section 3.2). Given the large differences in the shapes of the distributions of $\tilde{q}$ at a given $s$, the rescaling needed to reduce them to a single one needs to be much more complex than the one adopted in the last section. Therefore, we fit each single histogram with a different beta distribution, which has the following analytical expression:

$$ p(x, \alpha, \beta) = \frac{1}{B(\alpha, \beta)} x^{\alpha-1} (1-x)^{\beta-1}. \qquad (7) $$

This function has two shape parameters $\alpha$ and $\beta$; the factor $1/B(\alpha, \beta)$ is a normalization constant that can be computed by requiring that the integral of the probability distribution function is equal to unity.

From the fitting procedure, we obtained a pair of parameters for each bin in $s$; however, $\alpha$ has a complicated dependence on the first axial ratio (almost constant with an average value of $\alpha = 2.15$), while the mean value of the beta distributions $\mu = 1/(1 + \beta/\alpha)$ follows a linear relation. Fig. 9 shows the dependence of the mean $\mu$ (red squares on left-hand panel) and $\beta$ parameter (blue circles on right-hand panel) of the fitted beta functions on the first axial ratio
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distribution of $\tilde{q} = (q - s)/(1 - s)$ as function of $s$: the black boxes and whiskers represent the quartiles and 1.5 the quartiles range, respectively. The horizontal red error bars represent the bin inside which the medians have been computed. Right: distributions of $\tilde{q}$ for different values of $s$ (histograms) and fitting function resulting from the model presented in the section (curves).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Parameters of the fitted beta functions. Red is the mean of the distribution, in blue the second parameter $\beta$.</figcaption></figure>

$s$. The coloured lines in each respective panel show a fit of these two parameters:

$$ 
\begin{aligned}
\mu(s) &= 0.633s^{-0.007} \\
\beta(s) &= 1.389s^{-1.685}.
\end{aligned} 
\qquad (8) $$

These two equations give us a functional form of $p(\tilde{q}|s)$: starting from a value of $s$, one can retrieve the mean $\mu$ and $\beta$ from which the other parameter can be computed $\alpha = \beta/(1/\mu + 1)$. This gives what is needed to reconstruct the distribution of $\tilde{q}$ of a given $s$ and the scatter, if needed. The final step is to revert the change of coordinates and compute the physical axial ratio $q$.

## 4 EXPLODING THE MASS RANGE TO FIVE ORDERS OF MAGNITUDE

The next step of our work is to explode the recipes for dark matter halo shapes to lower masses; in the following sections we describe how to generalize the axial ratio distribution to a wider mass range. To do so, we combined the MXXL data with the SBARBINE simulations, a set of cosmological simulations that will allow us to study the shape of dark matter haloes from $3 \times 10^{10}$ to $6 \times 10^{15} M_{\odot} h^{-1}$.

As before, we express the mass dependence in terms of peak height $v$. By doing this, it is possible to treat homogeneously data from different redshifts and cosmologies, such as the SBARBINE and the MXXL simulations.

### 4.1 Axis ratio distribution: minor to major

On left-hand panel of Fig. 10, the logarithm of the minor to major axial ratio $s$ is shown as a function of the logarithm of $v$. As before, horizontal error bars represent the interval in $v$ and the box and whiskers are the quartiles and 1.5 the quartiles range for the combined sample, while coloured points are medians of individual catalogues. Again, there is no difference in the medians between redshifts, neither between the single simulations. It can be seen that $s$ has a nearly linear dependence on $\log(v)$, with a hint of flattening at both high and low masses.

For each bin, we extracted the probability distribution function of $\log(s)$ (right-hand panel of Fig. 10). The resulting curves exhibit an interesting pattern: high and low $v$ histograms are mirrored with respect to a central symmetric distribution which corresponds to $v \approx 1.21 (M \approx M_*)$. The rescaling adopted in Section 3.2 for cluster-size haloes does not compensate this large variation in the form of the distributions and it is not able to remove entirely the mass dependence. Instead of using a different rescaling relation to obtain a single probability distribution function (pdf), we decided to follow the same recipe we used for the second axial ratio; first of all we separately fit each distribution and then we relate the resulting parameters to the binning quantity. This is shown in Fig. 11, where we fit the mean (left-hand panel) and $\beta$ parameter (right-hand panel) of the Beta distributions we derived by fitting the histograms of the right-hand panel of Fig. 10. In order to keep the procedure simple we fit with a linear relation both $\mu$ and $\log \beta$,

$$ 
\begin{aligned}
\mu(v) &= -0.322 \log v + 0.620 \\
\log(\beta(v)) &= 0.560 \log v + 0.836.
\end{aligned} 
\qquad (9) $$

As before, the dependence of $\alpha$ is difficult to describe and it is almost constant with a value of about 11.21.

Using this fits we are now able to approximate the probability distribution function of the first axial ratio with a Beta function with parameters $\alpha = \beta/(1/\mu - 1)$ and $\beta$, over a range in mass of almost
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left: distribution of *s* as function of peak height for all the haloes selected from both redshifts the two simulations; the black boxes and whiskers represent the quartiles and 1.5 the quartiles range, respectively, computed within the bins shown by the horizontal error bars. The coloured points represent the medians for individual redshifts for the two simulations. Right: differential distribution of *s* for 8 bin in *v* (histograms) and the respective approximating functions obtained as shown in section 4.1 (curves).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Parameters of the fitted beta functions. Red is the mean of the distribution, in blue the second parameter $\beta$.</figcaption></figure>

<caption>**Figure 12.** Axis ratio *q* as function of *s* for different masses, represented by the points of different colours. Since there is no residual mass dependence in the conditional distribution, we get the same result as in the MXXL with all the simulations, confirming that this relation is universal.</caption>

six orders of magnitudes. Moreover the use of $v$ allows us to extend these results to different cosmologies and different redshifts.

## 4.2 Axis ratio distribution: intermediate to major

Finally, to fully describe a triaxial halo of a given mass the intermediate to major axis ratio has to be parametrized. As Fig. 12 shows, the relation between *q* and *s* at redshift *z* = 0 does not depend on the mass; the curves of different colours represent different mass bins and still trace the same relation. The fact that all the mass dependence is already inside *s*, allows us to use for $p(q|s)$ the same functional form of Section 3.3, independently of the mass we choose. The same applies to different redshifts (not shown here, but see Fig. 5 for a limited comparison), with the relation between

the two quantities being indistinguishable from the one in Fig. 12. Moreover, this independence of the conditional distribution from both mass and redshift is in agreement with the theoretical predictions from Rossi, Sheth & Tormen (2011).

# 5 COMPARISON WITH PREVIOUS WORKS

We have compared our results with measurement of axis ratios from other authors (Fig. 13). The data from both redshifts of the MXXL and SBARBINE simulations are shown with red squares, the median result from the analysis on cluster masses (Section 3.2) is the blue solid line and the green solid line is from the combined
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison between previous works (dashed lines) and the results of this paper (solid lines). Red squares represent the data from both redshifts of the MXXL and the SBARBINE simulations, converted to redshift $z = 0$ for the Millennium cosmology. The blue solid line is the model for clusters shown in Section 3.2; the green solid line is the fit for the entire mass interval from section 4.1. The dotted parts of the curves show the mass ranges outside where the relations have been derived from.</figcaption></figure>

data sets (Section 4.1). Results from other authors are shown with dashed lines in the mass range where their analysis was carried out and with dotted lines when extrapolated beyond it. Moreover all data and predictions have been converted to redshift $z = 0$ for the Millennium cosmology, when necessary. As it can be seen, there is a general agreement in the dependence of $s$ on the mass, with more massive haloes being less spherical. Although there seems to be a scatter of about 15 per cent, this is due more to the differences in the method of measuring shapes (different finders, radius, cleaning procedure), than an error on the measurement. It must be noticed that instead of the spherical mass, we used the mass within the ellipsoid for consistency reasons; yet, this does not substantially alter the findings presented here.

The most important difference comes from the radius at which the shape is measured. Jing & Suto (2002, blue dashed line) used particles of the isodensity surface corresponding to 2500$\delta_c$, roughly at a radius of $0.3R_{vir}$; this analysis is different from all the following authors, as it reflects the shape of an ellipsoidal shell, and not of all the mass inside the ellipsoid. Their mass range $6 \times 10^{12}–10^{14} M_{\odot} h^{-1}$ was also quite small compared to later analysis.

Studying a larger mass interval, $6 \times 10^{11}–3 \times 10^{14} M_{\odot} h^{-1}$, Allgood et al. (2006, yellow dashed line) derived axis ratios of particles distribution inside $0.3R_{vir}$ diagonalizing the normalized mass tensor (weighted by the distance from the centre); because of this their measure reflects the shape at an even closer radius.

On the other hand, Schneider et al. (2012, black dashed line) extended the analysis up to the virial radius, nevertheless the use of the normalized tensor prevents a meaningful comparison with our results.

All of these results are lower than what we derived, which can be explained by the fact that the shapes were measured at inner radii, where the particle distribution is supposed to be more elongate. However, if we restrict the comparison to works that used particles within the virial radius the agreement becomes much more strong. This is the case of Muñoz-Cuartas et al. (2011, magenta dashed

line), who studied shapes with an ellipsoidal overdensity algorithm similar to the one adopted in this work; their results agree with ours much more than any other work.

Finally, using a different type of halo finder, Bett et al. (2007, red dashed line) measured $s$ for a set of particles that represent all the bounded particles of a halo without assuming any particular shape; the finder also clean the sample for irregular haloes. The agreement with our results is another indication that the adopted selection criteria are justified and ellipsoids are a good approximation for regular haloes.

The other difference can arise from the cleaning of the sample; the green dashed line show the prediction from Despali et al. (2014), which is obtained from all haloes, regardless of their state of relaxation. As expected the values are lower, since the more unrelaxed haloes are typically irregular and so they appear more elongated with lower axial ratios. The difference is greater for less massive haloes.

# 6 SUMMARY AND CONCLUSIONS

We have studied the triaxiality of dark matter haloes from the MXXL simulation, which enabled us to characterize the shape of haloes with extremely good statistic in the galaxy clusters mass range, from $10^{14}$ to $4 \times 10^{15} M_{\odot} h^{-1}$. Using the SBARBINE simulations, we have extended our analysis to lower masses down to $3 \times 10^{10} M_{\odot} h^{-1}$, thus increased the mass range by more than five orders of magnitude. The main results of our analysis are the following:

(i) dark matter haloes are triaxial with a tendency of being prolate and in particular more massive objects are less spherical; as shown in Fig. 3 very unrelaxed haloes have the effect of artificially increasing the axis ratios and cannot be described by this simple ellipsoidal model, which is unimodal by construction;
(ii) for clusters, the distribution of the rescaled minor to major axis ratio is well described by a lognormal distribution, in contrast to previous extrapolations from lower masses that found a simple Gaussian fit;

(iii) over the whole examined mass range, *s* can be approximated by a beta distribution that depends only on the peak height *v*;

(iv) the conditional intermediate to major axis ratio distribution *p(q|s)* can also be described by a beta distribution that depends only on the first axis ratio and not on the mass, thus the same approach can be used for both clusters and the whole mass range of haloes;

(v) overall, the pdf of the shape of a dark matter halo is given by one single parameter *v*, related to its mass, that incorporates the dependence on redshift and cosmology. This goes in support of methods that allows us to change the cosmology of a numerical simulation (Angulo & White 2010), as within good approximation most of halo properties depend only on *v*.

In the recipe that we provide, a halo shape is determined only by its mass and can be changed to different cosmologies and redshifts. Depending on the level of precision desired, it is possible to choose different approximations,

(i) for a simpler analysis that is focused on the entire mass range, Section 4.1 presents a single method that can be applied to masses from $10^{10}$ up to $10^{16} M_{\odot} h^{-1}$. If restricted to masses lower than $10^{14} M_{\odot} h^{-1}$, this is actually a very accurate description of haloes shapes;

(ii) if the interest is only on clusters shapes, then Section 3.2 gives a more precise model;

(iii) finally, it is possible to combine the two description and just use the most suitable one given the mass of the halo, although losing the universality of the description.

A simple implementation of this model can be found on a dedicated website.<fnref n="2" />

In Section 5, we have compared our results with previous findings. There is a general agreement with previous works within a 15 per cent scatter that is due to the different methods used and especially to the radius at which the shape is measured. However, the picture is clear; dark matter haloes are triaxial objects and this effect is more prominent in clusters where the spherical model is quite far from being able to realistically represent the matter distribution.

## ACKNOWLEDGEMENTS

This work has been carried out thanks to the support of the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX* project (ANR-11-IDEX-0001-02) funded by the 'Investissements d'Avenir' French government programme managed by the ANR. ML acknowledges the Centre National de la Recherche Scientifique (CNRS) for its support. This study also benefited from the facilities offered by CeSAM (Centre de donnéeS Astrophysiques de Marseille http://lam.oamp.fr/cesam/). The 3D visualization was performed with GLNEMO2 program, developed by Jean-Charles LAMBERT at CeSAM - (http://projets.lam.fr/projects/glnemo2). GD has been partially financed by the Strategic Research Project AACSE (Algorithms and Architectures for Computational Science and Engineering) of the University of Padova. CG's research is part of the project GLENCO, funded under the European Seventh Framework Programme, Ideas, Grant Agreement no. 259349. CG and RA thank

LAM for supporting their visits during which part of this work has been done. GD and CG thank the whole cosmology group of the University of Padova with whom the SBARBINE simulations were designed; in particular, we thank Giuseppe Tormen for providing the computational resources and Giacomo Baso for running *Ada*.

## REFERENCES

* Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781

* Angulo R. E., White S. D. M., 2010, MNRAS, 405, 143

* Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046

* Bailin J., Steinmetz M., 2005, ApJ, 627, 647

* Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215

* Binggeli B., 1982, A&A, 107, 338

* Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440

* Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150

* Buote D. A., Canizares C. R., 1992, ApJ, 400, 385

* Buote D. A., Canizares C. R., 1996, ApJ, 457, 565

* Carroll S. M., Press W. H., Turner E. L., 1992, ARA&A, 30, 499

* Carter D., Metcalfe N., 1980, MNRAS, 191, 325

* Cole S., Lacey C., 1996, MNRAS, 281, 716

* Despali G., Tormen G., Sheth R. K., 2013, MNRAS, 431, 1143

* Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208

* Doroshkevich A. G., 1970, Astrophysics, 6, 320

* Dubinski J., Carlberg R. G., 1991, ApJ, 378, 496

* Eke V. R., Cole S., Frenk C. S., 1996, MNRAS, 282, 263

* Evans A. K. D., Bridle S., 2009, ApJ, 695, 1446

* Fabricant D., Rybicki G., Gorenstein P., 1984, ApJ, 286, 186

* Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507

* Gao L., Navarro J. F., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 425, 2169

* Giocoli C., Moreno J., Sheth R. K., Tormen G., 2007, MNRAS, 376, 977

* Giocoli C., Tormen G., van den Bosch F. C., 2008, MNRAS, 386, 2135

* Giocoli C., Meneghetti M., Bartelmann M., Moscardini L., Boldrin M., 2012a, MNRAS, 421, 3343

* Giocoli C., Meneghetti M., Ettori S., Moscardini L., 2012b, MNRAS, 426, 1558

* Hopkins P. F., Bahcall N. A., Bode P., 2005, ApJ, 618, 1

* Jing Y.-P., Fang L.-Z., 1994, ApJ, 432, 438

* Jing Y., Suto Y., 2002, ApJ, 574, 538

* Kasun S. F., Evrard A. E., 2005, ApJ, 629, 781

* Kawahara H., 2010, ApJ, 719, 1926

* Klypin A., Yepes G., Gottlober S., Prada F., Hess S., 2014, preprint (arXiv:1411.4001)

* Lacey C., Cole S., 1993, MNRAS, 262, 627

* Lau E. T., Nagai D., Kravtsov A. V., Vikhlinin A., Zentner A. R., 2012, ApJ, 755, 116

* Lewis A., Challinor A., Lasenby A., 2008, ApJ, 2, 2

* Limousin M., Morandi A., Sereno M., Meneghetti M., Ettori S., Bartelmann M., Verdugo T., 2013, Space Sci. Rev., 177, 155

* Ludlow A. D., Navarro J. F., Li M., Angulo R. E., Boylan-Kolchin M., Bett P. E., 2012, MNRAS, 427, 1322

* Macciò A. V., Dutton A. A., van den Bosch F. C., Moore B., Potter D., Stadel J., 2007, MNRAS, 378, 55

* Meneghetti M. et al., 2014, ApJ, 797, 34

* Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584

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

* Neto A. F. et al., 2007, MNRAS, 381, 1450

* Oguri M., Takada M., Okabe N., Smith G. P., 2010, MNRAS, 405, 2215

* Oguri M., Bayliss M. B., Dahle H., Sharon K., Gladders M. D., Natarajan P., Hennawi J. F., Koester B. P., 2012, MNRAS, 420, 3213

* Paz D. J., Lambas D. G., Padilla N., Merchán M., 2006, MNRAS, 366, 1503

<footnote n="2">² http://wiki.lam.fr/triaxial</footnote>
Planck Collaboration XVI, 2014, A&A, 571, A16

Prada F., Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012, MNRAS, 423, 3018

Press W. H., Schechter P., 1974, ApJ, 187, 425

Rossi G., Sheth R. K., Tormen G., 2011, MNRAS, 416, 248

Sayers J., Golwala S. R., Ameglio S., Pierpaoli E., 2011, ApJ, 728, 39

Schneider M. D., Frenk C. S., Cole S., 2012, J. Cosmol. Astropart. Phys., 2012, 030

Skielboe A., Wojtak R., Pedersen K., Rozo E., Rykoff E. S., 2012, ApJ, 758L, 16

Soucail G., Fort B., Mellier Y., Picat J. P., 1987, A&A, 172, L14

Springel V., 2005, MNRAS, 364, 1105

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726

Tormen G., Moscardini L., Yoshida N., 2004, MNRAS, 350, 1397

van Daalen M. P., Angulo R. E., White S. D. M., 2012, MNRAS, 424, 2954

Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405

## APPENDIX A: DENSITY PEAK HEIGHT

In this appendix, we describe step by step how to compute density peak height $v$ for a virialized halo with mass $M$ at redshift $z$ for a given cosmological model. Its definition is the following:

$$v \equiv \frac{\delta_c(z)}{\sigma(M)}, \qquad (A1)$$

where $\delta_c(z)$ is the critical overdensity of the spherical collapse model, the initial density required for a fluctuation to collapse at redshift $z$. This in turn can be expressed as the collapse overdensity at redshift $z=0$ rescaled to a given time: $\delta_c(z) = \delta_c/D(z)$, with $D(z)$ being the linear growth rate of a density fluctuation normalized to unity at $z=0$. The overdensity $\delta_c$ depends only on redshift and not on the mass; on the other hand, the denominator $\sigma(M)$, depends on the mass but not on redshift. It is the variance in the initial density field smoothed on a linear scale $R$, which corresponds to the radius of a uniform sphere of mass $M$. Therefore, only the linear growth rate $D(z)$ and the initial power spectrum $P(k)$ are needed.

From the linear perturbation theory, it is possible to compute $D(z)$

$$D(z) \propto H(t) \int_{0}^{t} \frac{dt'}{a^2(t')H^2(t')}, \qquad (A2)$$

which has to be solved numerically. Fortunately, there is an approximated solution (Carroll, Press & Turner 1992) that can be expressed as $D(z) \propto g(z)/(1+z)$, where

$$g(z) = \frac{5/2 \Omega_m(z)}{\Omega_m^{4/7} - \Omega_\Lambda(z) + [1 + \Omega_m(z)/2][1 + \Omega_\Lambda(z)/70]} \qquad (A3)$$

Additionally, the collapse overdensity has an extremely weak dependence on cosmology: $\delta_c \approx 1.686[\Omega_m(t_c)]^{0.0055}$; for realistic cosmologies this can be approximated to $\delta_c \approx 1.69$. Therefore, at $z=0$ the collapse overdensity is $\delta_c$ and it increases with redshift, due to $D(z)$.

The other quantity required, the variance $\sigma^2(M)$, is defined from the power spectrum as

$$\sigma^2(M) = \frac{1}{2\pi^2} \int_0^\infty P(k) \tilde{W}^2(kR) k^2 dk; \qquad (A4)$$

where $\tilde{W}$ is the Fourier transform of a window function. Typically, $W$ is a Top Hat (sphere) in the coordinates space, so that its Fourier transform $\tilde{W}$ is

$$\tilde{W}(kR) = 3 \frac{\sin(kR) - kR \cos(kR)}{(kR)^3}; \qquad (A5)$$

with the radius $R$ given by $M = \rho_b 4\pi/3 R^3$. The power spectrum $P(k)$ of the density fluctuations is the main input; given a set of cosmological parameters it can be computed from a software like CAMB (Lewis et al. 2008). As it is function of initial conditions only, $\sigma(M)$ needs to be computed only once for a given cosmology; all the redshift dependence is inside $D(z)$.

Finally, for a halo of mass $M$, using equation (A4) it is possible to compute $\sigma(M)$ and combine it with the value of $D(z)$ from equation (A3) to obtain the correct density peak height $v$.

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