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

HAL
open science

# On the shape of dark matter haloes from MultiDark Planck simulations

Jesús Vega-Ferrero, Gustavo Yepes, Stefan Gottlöber

► To cite this version:

Jesús Vega-Ferrero, Gustavo Yepes, Stefan Gottlöber. On the shape of dark matter haloes from MultiDark Planck simulations. Monthly Notices of the Royal Astronomical Society, 2017, 467 (3), pp.3226-3238. 10.1093/mnras/stx282. hal-01554257

HAL Id: hal-01554257

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

Submitted on 27 Aug 2025

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

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

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

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

# On the shape of dark matter haloes from MultiDark Planck simulations

Jesús Vega-Ferrero,<fnref n="1" />,<fnref n="2" />* Gustavo Yepes<fnref n="1" />,<fnref n="3" /> and Stefan Gottlöber<fnref n="4" />

Accepted 2017 January 31. Received 2017 January 16; in original form 2016 March 31

## ABSTRACT

The halo shape plays a central role in determining important observational properties of the haloes such as mass, concentration and lensing cross-sections. The triaxiality of lensing galaxy clusters has a substantial impact on the distribution of the largest Einstein radii, while weak lensing techniques are sensitive to the intrinsic halo ellipticity. In this work, we provide scaling relations for the shapes of dark matter haloes as a function of mass (peak height) and redshift over more than four orders of magnitude in halo masses, namely from $10^{11.5}$ to $10^{15.8} h^{-1} M_{\odot}$. We have analysed four dark matter only simulations from the MultiDark cosmological simulation suite with more than 56 billion particles within boxes of 4.0, 2.5, 1.0 and $0.4 h^{-1}$ Gpc size assuming Planck cosmology. The dark matter haloes have been identified in the simulations using the ROCKSTAR halo finder, which also determines the axis ratios in terms of the diagonalization of the inertia tensor. In order to infer the shape for a hypothetical halo of a given mass at a given redshift, we provide fitting functions to the minor-to-major and intermediate-to-major axis ratios as a function of the peak height.

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

## 1 INTRODUCTION

According to the current paradigm of structure formation isolated galaxies, groups and clusters of galaxies are hosted by dark matter haloes of increasing mass. Quantum fluctuations generated during an early, inflationary epoch of the Universe become classical density fluctuations: the seed of today's observed large-scale structures. In an expanding Universe dominated by weakly interacting dark matter, these fluctuations grow with time due to the gravitational instability that leads to the formation of self-gravitating dark matter haloes. The first haloes emerge from small-scale fluctuations and progressively merge to form later generations of more massive haloes, where groups and clusters of galaxies are hosted. These haloes can be associated with the peaks in the primordial Gaussian density field of dark matter, and their mass distributions can be derived in terms of the statistics of Gaussian random fields and the spherical collapse model of halo formation (Press & Schechter 1974).

Contrary to the model assumption, the mass accretion on to haloes comes along a preferential direction (mostly along a filament) and tends to be clumpy, i.e. massive haloes form by merging with smaller haloes. For these reasons, one should not expect that dark matter haloes show a spherical shape if their relaxation times are longer

than the time between mergers or their accretion events occur from a preferred direction. For the most massive haloes that host clusters of galaxies, observational studies based on optical and density maps (Carter & Metcalfe 1980; Binggeli 1982), X-ray data (Fabricant, Rybicki & Gorenstein 1984; Buote & Canizares 1992, 1996; Kawahara 2010; Lau, Nagai & Nelson 2013), Sunyaev–Zel'dovich pressure maps (Sayers et al. 2011), strong and weak gravitational lensing methods (Soucail et al. 1987; Evans & Bridle 2009; Oguri et al. 2010, 2012) indicate that cluster-size dark matter haloes are not spherical at all.

Additionally, the study of dark matter haloes formed in cosmological simulations has shown that their shapes are well described by a triaxial model (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). Clusters of galaxies seem to be less spherical towards their centres (Warren et al. 1992; Jing & Suto 2002; Bailin & Steinmetz 2005; Allgood et al. 2006; Schneider et al. 2012), probably due to interactions with other clusters or groups of galaxies. Besides, triaxiality of clusters increases with increasing both the mass and the redshift of the cluster (Jing & Suto 2002; Allgood et al. 2006; Muñoz-Cuartas et al. 2011), essentially because they are affected by the direction of the last major merger or the presence of filaments around them, and because

<footnote n="1">¹Departamento de Física Teórica, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="2">²LERMA, CNRS UMR 8112, Observatoire de Paris, 61 Avenue de l'Observatoire, F-75014 Paris, France</footnote>

<footnote n="3">³Astro-UAM, UAM, Unidad Asociada CSIC, E-28006 Madrid, Spain</footnote>

<footnote n="4">⁴Leibniz-Institut für Astrophysik, D-14482 Potsdam, Germany</footnote>

<footnote>* E-mail: astrovega@gmail.com</footnote>
massive haloes were formed later on (Despali, Giocoli & Tor-
men 2014). This scenario has been confirmed by the recent work
of Bonamigo et al. (2015) on the shape of simulated dark mat-
ter haloes extracted from the Millennium XXL simulation (Angulo
et al. 2012).

It is well known that the triaxiality of dark matter haloes is par-
ticularly important for the study of clusters of galaxies. In fact, the
triaxiality of dark matter haloes is supposed to play a central role in
several important properties: cluster mass and concentration, inner
slope of the dark matter density profile and strong lensing cross-
sections (Bartelmann et al. 1998; Meneghetti et al. 2000, 2001;
Oguri et al. 2005; Giocoli et al. 2012a,b; Wojtak 2013). In par-
ticular, the triaxiality of lensing galaxy clusters has a substantial
impact on the distribution of the largest Einstein radii (Waizmann,
Redlich & Bartelmann 2012; Waizmann et al. 2014; Vega-Ferrero
et al., in preparation). Moreover, semi-analytic models designed to
derive the gravitational lensing properties in galaxy clusters gen-
erally described the mass distribution by a triaxial model as pro-
posed by Jing & Suto (2002). The results presented by Jing & Suto
(2002) were obtained from numerical simulations that hardly con-
tained any halo above 10¹⁴ h⁻<fnref n="1" /> M☉. Therefore, arc statistics from
galaxy clusters inferred by these models (see Giocoli et al. 2012a;
Redlich et al. 2012, for some examples) rely on extrapolations from
lower mass haloes. In order to study the impact of the triaxiality,
Waizmann et al. (2012) introduced a cut-off in the distribution of
axis ratios of Jing & Suto (2002) to remove extreme axis ratios,
finding that the impact of the tail of the axis ratio distribution on
the largest Einstein radii distribution is substantial. However, due
to the limited knowledge of the statistics of haloes with extremely
small axis ratios, it is not possible to clearly define a proper choice
of the cut-off (if present) until the triaxiality distributions of large
halo samples (up to cluster masses) are derived from large volume
numerical simulations.

Apart from some analyses over few massive clusters (Muñoz-
Cuartas et al. 2011; Schneider et al. 2012; Despali et al. 2014;
Groener & Goldberg 2014), the work of Bonamigo et al. (2015)
is the only one that presents statistically significant predictions on
the shapes of haloes with masses above 3 × 10¹⁴ h⁻<fnref n="1" /> M☉. The
MultiDark suite of cosmological simulations (Klypin et al. 2016)
– in particular the Huge-MultiDark (HMD) simulation, containing
almost 70 billion particles in a 4 h⁻<fnref n="1" /> Gpc box volume – allows us
to study in detail the shape of simulated dark matter haloes in the
cluster mass range with unprecedented statistics and in the most up-
to-date cosmological parameters consistent with the recent results
from the Planck Collaboration.

Moreover, weak gravitational lensing (see Bartelmann & Schneider 2001, for more details) has been successfully used in a wide variety of applications: from the mass calibration of galaxy clusters (Ford et al. 2015; Hoekstra et al. 2015; Kettula et al. 2015; Applegate et al. 2016; Medezinski et al. 2016; Mo et al. 2016) to the measurement of the ellipticity of dark matter haloes (Natarajan & Refregier 2000; Evans & Bridle 2009; Oguri et al. 2010, 2012; Schrabback et al. 2015). Since weak lensing is sensitive to halo ellipticity (or triaxiality), a detailed study of the intrinsic shapes of galaxy-scale simulated haloes is essential for present and future weak lensing studies, especially in the time of ‘precision cosmology’.

In this paper, we aim at characterizing the shape of simulated
dark matter haloes over more than four orders of magnitude in halo
masses. To this end,

(i) we have analysed four high-resolution dark matter only sim-
ulations within boxes of 4.0, 2.5, 1.0 and 0.4 h⁻<fnref n="1" /> Gpc;

(ii) we examine the shapes of haloes with more than 3000 par-
ticles and provide predictions within the halo mass range: 3 ×
10¹<fnref n="1" /> ≲ M ≲ 8 × 10¹⁵ h⁻<fnref n="1" /> M☉;

(iii) we investigate the radial dependence of the halo shapes by
measuring the axis ratios: at virial radius (defined at an overdensity
Δ<sub>vir</sub>ρ<sub>b</sub>) and at an overdensity of 500ρ<sub>c</sub>, along with the alignment of
the major axis at these overdensities;

(iv) additionally, we examine the probability distribution of the
minima (minor-to-major) axis ratio.

This paper is organized as follows. In Section 2, we describe the
set of simulations. In Section 3, we describe the different methods
used to measure the shape of dark matter haloes and present our
results for the four simulated boxes. In Section 4, we compare our
findings for the minor-to-major axis ratios as a function of halo
mass with previous studies. Finally, in Section 5, we summarize our
results and conclusions.

## 2 HALO CATALOGUES

We derived the shape of simulated dark matter haloes over more
than four orders of magnitude in halo masses. To this end, we anal-
ysed four dark matter only simulations from the MultiDark cos-
mological simulation suite. This set of simulations consist on four
cosmological volumes of 4.0, 2.5, 1.0 and 0.4 h⁻<fnref n="1" /> Gpc cube sides,
which are labelled here as Huge-MultiDark (HMD), Big-MultiDark
(BMD), New-MultiDark (NMD) and Small-MultiDark (SMD), respec-
tively. The largest simulation volume (i.e. the HMD) contains 4096³
particles, while the other simulations are performed with 3840³ par-
ticles. All the simulations were carried out with L-GADGET-2 code
(Springel 2005) and were performed assuming the same Λ cold
dark matter (ΛCDM) cosmological model with *Planck* cosmologi-
cal parameters (Planck Collaboration XVI 2014): (Ω<sub>M</sub>, Ω<sub>b</sub>, Ω<sub>Λ</sub>, σ<sub>8</sub>, n<sub>s</sub>, h) = (0.307, 0.048, 0.693, 0.829, 0.96, 0.678). For these reasons,
the MultiDark suite constitutes a coherent and homogeneous data
set of cosmological simulations. The particle masses along with the
parameters for each simulation are shown in Table 1. All the infor-
mation concerning these simulations is available at the MultiDark
data base<fnref n="1" /> and the CosmoSim data base.<fnref n="2" />

Dark matter haloes are identified using the **ROCKSTAR** halo finder (Behroozi, Wechsler & Wu 2013). From the **ROCKSTAR** halo catalogues, we select all the distinct haloes at four different redshifts, $z = (0.00, 0.33, 0.66, 1.00)$. Among other parameters, the **ROCKSTAR** halo finder provides virial masses and the length of the minor, intermediate and major axis. Thereafter, we use the definition of the virial mass as the mass inside the radius encompassing a given density contrast $\Delta_{vir}\rho_b$ and 500$\rho_c$, where $\rho_b$ and $\rho_c$ are the mean and the critical density of the Universe, respectively. The overdensity $\Delta_{vir}$ is computed using the approximation given by Bryan & Norman (1998). For the given cosmology, $\Delta_{vir} \approx 333$ times the mean density of the Universe at $z=0$.

To ensure a good resolution we limited our analysis to distinct haloes with more than 3000 particles within $R_{vir}$. A sample build in this way also constitutes a complete volume-limited sample of the HMD, BMD, NMD and SMD box volumes. In Fig. 1, we show the halo mass function for the different simulations with the *Planck* cosmology. The Tinker et al. (2008) mass function provides an excellent fit to $z=0$ results and slightly underestimates the mass function at $z=1$.

<footnote n="1">¹ http://www.multidark.org</footnote>

<footnote n="2">² https://www.cosmosim.org</footnote>
<table><caption>**Table 1.** Numerical parameters for the simulations. The columns correspond to: the simulation identifier; the size of the simulated box in $h^{-1}$ Gpc; the number of particles and the mass per simulation particle $m_p$ in units $h^{-1}$ M$_\odot$; the Plummer equivalent gravitational softening length $\epsilon$ in units of physical $h^{-1}$ kpc; the number of haloes at $z=0$; the number of relaxed haloes at $z=0$ with more than 3000 particles; the fraction of relaxed haloes at $z=0$.</caption><thead><tr><th>Simulation</th><th>Box</th><th>Particles</th><th>$m_p$</th><th>$\epsilon$</th><th>$N_h$ ($z=0$)</th><th>$N_{rel}$ ($z=0$)</th><th>$f_{rel}$ ($z=0$)</th></tr></thead><tbody><tr><td>SMD</td><td>0.4</td><td>3840<sup>3</sup></td><td>9.6 × 10<sup>7</sup></td><td>1.5</td><td>814 977</td><td>638 414</td><td>0.78</td></tr><tr><td>NMD</td><td>1.0</td><td>3840<sup>3</sup></td><td>1.5 × 10<sup>9</sup></td><td>5</td><td>1004 645</td><td>720 799</td><td>0.72</td></tr><tr><td>BMD</td><td>2.5</td><td>3840<sup>3</sup></td><td>2.4 × 10<sup>10</sup></td><td>10</td><td>692 647</td><td>432 815</td><td>0.62</td></tr><tr><td>HMD</td><td>4.0</td><td>4096<sup>3</sup></td><td>7.9 × 10<sup>10</sup></td><td>25</td><td>387 849</td><td>219 567</td><td>0.57</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Halo mass function in simulations with the Planck cosmological parameters. Different points correspond to the different simulations at $z=0$ and 1. Error bars are computed considering that the error in the number of haloes is $\Delta n = \sqrt{n}$. The full curves correspond to the Tinker et al. (2008) mass function. It gives an excellent fit to $z=0$ results and slightly underestimates the mass function at $z=1$.</figcaption></figure>

Additionally, the modelling of non-spherical objects is not simple and can lead to uncertain results. Since highly elongated haloes are mostly due to unrelaxed or merging haloes, we divide our sample in relaxed and unrelaxed haloes. Relaxed haloes are defined by these three criteria (Klypin et al. 2016):

$$r_{\text{off}} < 0.07, \quad \lambda < 0.07, \quad 2K/|W| < 1.5, \qquad (1)$$

where $\lambda$ is the spin parameter as introduced by Peebles (1969), $2K/|W|$ is the virial parameter and $r_{\text{off}} = |r_{\text{cm}} - r_{\phi}|/R_{\text{vir}}$ is defined as the offset between the centre of mass ($r_{\text{cm}}$) and the centre of the potential ($r_{\phi}$) in units of the virial radius. The fraction of relaxed haloes as a function of the halo mass for the different simulations are shown in Fig. 2. At redshift $z=0$, this fraction drops from 80 per cent down to 50 per cent in the mass range considered in this work. The merger rate of massive haloes is usually larger and more violent compared to lower mass haloes, which explains the negative slope in the relaxed fraction. At $z=1$, this fraction is smaller in about 25 per cent than that at $z=0$ within the whole mass range. The fraction of relaxed haloes at $z=0$ for each simulation are shown in Table 1.

## 3 HALO SHAPES

The common procedure to determine the shape of haloes is to model them as ellipsoids. The length of the minor, intermediate

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Fraction of relaxed haloes as a function of the halo mass. Different points correspond to the median relaxed fraction for the different simulated boxes at $z=0$ and 1. Vertical error bars are computed considering that the error in the number of haloes is $\Delta n = \sqrt{n}$, while horizontal error bars enclosed the Q1 and Q3 quartiles. The solid and dashed lines indicate the linear best fit to the data points for $z=0$ and 1, respectively.</figcaption></figure>

and major axis ($c \le b \le a$) of the ellipsoid can be computed by diagonalizing the inertia tensor. In this work, we determine the axis ratios using two forms of the inertia tensor that are commonly used in the literature: the unweighted inertia tensor,

$$I_{ij} = \sum_{n} x_{i,n}x_{j,n}; \qquad (2)$$

and the weighted (or reduced) inertia tensor,

$$\tilde{I}_{ij} = \sum_{n} \frac{x_{i,n}x_{j,n}}{r_{n}^{2}}, \qquad (3)$$

where

$$r_n = \sqrt{x_n^2 + y_n^2/q^2 + z_n^2/s^2} \qquad (4)$$

is the elliptical distance in the eigenvector coordinate system from the halo centre to the nth particle. The square root of the eigenvalues obtained from the diagonalization of both expressions provides the axis ratios ($c \le b \le a$), while the corresponding eigenvectors determine their orientations. Halo centres are identified as the densest point in phase space, which normally corresponds to the densest point in position space.

The minor-to-major and intermediate-to-major axis ratios are defined as $s \equiv c/a$ and $q \equiv b/a$, respectively. Hereafter, we denote the virial axis ratios measured at an overdensity $\Delta_{\text{vir}}\rho_b$ (i.e. $R_{\text{vir}}$) as $s$ and $q$. Additionally, the axis ratios measured at an overdensity of $500\rho_c$ ($R_{500}$) are denoted as $s_{500}$ and $q_{500}$.
<figure><img src="image_6.png" /><figcaption>Figure 3. Left-hand panel: distribution of virial axis ratios for HMD relaxed haloes ($z=0$) obtained by different methods. Central panel: distribution of virial axis ratios in the plane $s_w$ versus $s_{\text{sph}}$. The black circles and error bars are the running means and $1\sigma$ levels in eight bins in $s_{\text{sph}}$. The grey contours overlaid to the image show the intensity levels corresponding to 1, 10, 50 and 90 per cent of the probability peak. The dashed black line indicates the result of the fit with a power law to the data points ($s_w \simeq s_{\text{sph}}^{1.86}$), while red dotted line shows the findings of Bailin & Steinmetz (2005). The solid black line corresponds to the linear best fit to the data points. Right-hand panel: distribution of virial axis ratios in the plane $s_w$ versus $s_{\text{un}}$. The black circles and error bars are the running means and $1\sigma$ levels in eight bins in $s_w$. The solid black line corresponds to the linear best fit to the data points, while the dashed black line shows the results of Allgood et al. (2006) for all haloes (both relaxed and unrelaxed ones).</figcaption></figure>

The procedure used to recover the shape begins by determining the inertia tensor with $s=1$ and $q=1$, including all the particles inside a sphere of a given overdensity. Subsequently, new values for $s$ and $q$ are computed as the volume analysis is deformed along the eigenvectors (proportionally to the eigenvalues), while keeping the major axis equal to the original spherical radius. After the deformations of the original spherical region, the inertia tensor is calculated again using the new values of $s$ and $q$, and including only the particles within the new ellipsoidal region. This iterative process can be repeated until the variance in both axis ratios, $s$ and $q$, is less than a given tolerance or up to a maximum number of iterations (see Allgood et al. 2006, for more details). In ROCKSTAR, the halo centres are not recalculated within the iterative procedure. In this work, we derive the axis ratios using three different methods:

(i) a weighted inertia tensor in a spherical window (i.e. only one iteration) with axis ratios labelled as $s_{\text{sph}}$;

(ii) an iterative weighted inertia tensor (10 iterations) with axis ratios labelled as $s_w$;

(iii) an iterative unweighted or reduced inertia tensor (10 iterations) with axis ratios labelled as $s_{\text{un}}$.

In order to compare between the different methods, we computed the axis ratios considering the three methods described above for relaxed haloes only in the HMD simulation at $z=0.00$. In Fig. 3, we show the distribution of the different virial axis ratios. For a direct comparison, we also show a scatter plot for the axis ratios $s_w$ versus $s_{\text{sph}}$. We perform a power-law fit to the data finding $s_w = s_{\text{sph}}^{\gamma}$ with $\gamma = 1.86 \pm 0.13$ and $R^2 = 0.969$ (where $R^2$ is the coefficient of determination). A similar result was found in Bailin & Steinmetz (2005) with $\alpha = \sqrt{3} \approx 1.73$. However, our data points are best fit by a linear relation

$$s_w = (-0.13 \pm 0.02) + (0.89 \pm 0.02)s_{\text{sph}}, \qquad (5)$$

with higher $R^2 = 0.996$ and $\sigma = 0.02$. We also examined the difference in the determination of the virial axis ratios derived with the weighted ($s_w$) and the unweighted ($s_{\text{un}}$) iterative inertia tensor in Fig. 3 (right-hand panel). Data points are best fit by a linear relation as well

$$s_w = (0.096 \pm 0.002) + (0.755 \pm 0.004)s_{\text{un}}, \qquad (6)$$

with a dispersion of $\sigma = 0.004$. The weighted axis ratios are $\sim 15$ per cent lower than those obtained with the unweighted inertia tensor.

Although not plotted for clarity, we do not find any significant dependence on the halo mass in the above relations. For comparison, we also show the results obtained by Allgood et al. (2006) with a slope of about 0.87. Although the slope is different than our findings ($\sim 0.75$), we recall that Allgood et al. (2006) included the whole sample in the analysis (both relaxed and unrelaxed haloes). Besides, they measured the axis ratios at an inner radius ($R \sim 0.3R_{\text{vir}}$), while in the right-hand panel of Fig. 3, we determine the axis ratios at $R = R_{\text{vir}}$, which can lead to differences in the measurement of the axis ratios, particularly for the unweighted inertia tensor.

Since the weighted inertia tensor weights the contribution from each particle in the sum by the distance to the particle squared, it is less sensitive to large substructure in the outer regions of the volume analysed. Moreover, the unweighted inertia tensor is more biased than the weighted inertia tensor due to particles at large radii. For these reasons, hereafter, we show the axis ratios derived from the weighted inertia tensor (i.e. $s = s_w$) for the four simulated volumes at $z = (0.00, 0.33, 0.66, 1.00)$.

### 3.1 Minor-to-major axis ratio

As already proposed by many authors (Jing & Suto 2002; Allgood et al. 2006; Bett et al. 2007; Schneider et al. 2012), the axis ratios depend on the mass and the redshift of the halo; however, this dependence needs to be tested for the most massive haloes in large simulation boxes. Recently, Bonamigo et al. (2015) proposed an interesting procedure to account for the redshift dependence of the halo shapes in terms of the peak height $v = \delta_c/\sigma(M)$, where $\delta_c$ is the critical overdensity of the spherical collapse model and $\sigma(M)$ is the variance in the initial density field smoothed on a linear scale of a uniform sphere of mass $M$. Within this framework, the characteristic mass ($M_*$) is defined as the mass giving $v(z, M_*) = 1$ and, therefore, dark matter haloes with $v > 1$ will have masses larger than the typical haloes collapsing at a given redshift. See Appendix A for the details on how to compute the peak height $v$ for a halo of a given mass $M$.

In Fig. 4, we show the axis ratio $s$ distributions as a function of the peak height $v$ of the halo for the four simulated box volumes (only relaxed haloes). Since the merger rate of massive haloes is larger, we expect them to show smaller axis ratio than low-mass haloes (with median $s \approx 0.4$), as already pointed out by Jing &
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Virial axis ratio *s* as a function of the peak height in log–log space for the four simulated boxes. Different points correspond to the median values of *s* in each bin considered. The dotted lines enclose the Q1 and Q3 quartiles for the distribution of *s* at *z* = 0.</figcaption></figure>

Suto (2002), Allgood et al. (2006), Muñoz-Cuartas et al. (2011) and Despali et al. (2014). Moreover, haloes at higher redshift are also more elongated, in agreement with the findings of Muñoz-Cuartas et al. (2011) and Despali et al. (2014). At the same time, axis ratios for a given halo mass are slightly smaller as the simulated box volume increases (mostly at *z* = 1). This effect can be explained because the halo shape is noisier (and therefore more likely to be aspherical) for lower particle counts. Nevertheless, these differences are less than 20 per cent of the interquartile range of the relation (i.e. IQR = Q3 − Q1).

It is also important to notice that using the same procedure as in Bonamigo et al. (2015), we do not completely remove the redshift dependence on the axis ratio distributions (Fig. 4). For this reason, we proposed a modified scaling relation between the peak height *v* and the axis ratio *s* in log–log space, which explicitly includes the redshift dependence as function of the linear growth rate *D(z)* as follows:

$$ \log(s) = b + a \log(v'), \quad (7) $$

where we have defined $v' = vD(z)^c$. Since $D(0) = 1$ by definition, this scaling relation is identical to the one presented by Bonamigo et al. (2015) only at $z = 0$ (or in the case of $c = 0$). Therefore, the redshift dependence of the scaling relation $s-v'$ can be assessed considering that

$$ v'(M, z) = v(M, z)D(z)^c = v_0(M)D(z)^{c-1}. \quad (8) $$

In Fig. 5, we show the median values for the axis ratio *s* as a function of *v'* for the relaxed haloes in the four simulated boxes. We binned the haloes according to their *v* in log–log space and considering $\Delta\log v = 0.06$. The error bars in the figure correspond to the 25 and 75 per cent probabilities of *v'* in each bin. The best fit by equation (7) to the data points is obtained with $ac = 0.08 \pm 0.01$, which reveals a stronger dependence with redshift than found by Bonamigo et al. (2015). The standard deviation of the residuals of the fit is $\sigma = 0.008$, defined as

$$ \sigma = \frac{1}{N_{\text{dof}}} \sum_{i} R_{i}^{2}, \quad (9) $$

where $R_i$ is the residual measured in the *i*th bin and $N_{\text{dof}}$ is the number of degrees of freedom (i.e. the number of bins minus the

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Virial axis ratio *s* as function of the peak height (expressed as *v'* ) in log–log space for relaxed haloes in the four simulated boxes. Different points correspond to the median values of the distribution at different redshifts, while error bars enclose the Q1 and Q3 quartiles of the peak height in each bin considered. Dashed lines show the Q1 and Q3 quartiles around the median values of *s* in each bin of *v*. The solid black line shows the linear fit in log–log space to the data points using equation (7), while the shade grey region corresponds to the standard deviation of the residuals (equation 9).</figcaption></figure>

number of free parameters in the fit). The procedure used to perform the fit is a general-purpose optimization based on Nelder & Mead (1965) algorithm.

The best-fitting parameters are shown in Table 2 for relaxed haloes, unrelaxed haloes and the whole sample. We also show the results for the axis ratio measured at $\Delta = 500\rho_c$ as function of the corresponding peak height ($v_{500}$) for a halo of a given mass $M_{500}$ at redshift *z*.

Given the linear dependence (in log scale) of the axis ratio as a function of the peak height ($v'$), we derived the scaled axis ratios in a similar way as performed by previous authors: $\log(s) = a \log(v') + b \Rightarrow \tilde{s} = 10^b = 10^{\log s - a \log(v')} = s (v')^{-a}$. The black line in Fig. 5 corresponds to the best fit to the median axis ratios for relaxed haloes with slope $a = -0.295$. The *y*-intercept point is the logarithm of the median axis ratio at $v' = 1$, which corresponds to $M = M_*$ at $z = 0$. Thus, the fit yields to $\tilde{s}(z = 0, M_*) = 10^b \approx 0.59$.

The scaled axis ratio distribution $\tilde{s}$ of relaxed haloes shows a fairly universal distribution almost independently on the mass and redshift of the haloes for the four simulated boxes (Fig. 6). The probability distribution function of the scaled axis ratio ($\tilde{s}$) is well fitted by a lognormal distribution as proposed by Bonamigo et al. (2015):

$$ p(\tilde{s}, \mu, \sigma) = \frac{1}{\tilde{s}\sqrt{2\pi}\sigma} \exp\left(-\frac{(\ln\tilde{s} - \mu)^2}{2\sigma^2}\right). \quad (10) $$

The parameters of the fitted function are $\mu = -0.52$ and $\sigma = 0.19$, which can be converted to the median and the standard deviation of the Gaussian distribution:

$$ \text{median} = e^{\mu} = 0.59, $$

$$ \text{std} = \sqrt{(e^{\sigma^2} - 1)e^{2\mu + \sigma^2}} = 0.11. \quad (11) $$

The results for the scaled axis ratio of relaxed haloes are in good agreement with those obtained by Bonamigo et al. (2015), which found a median of 0.61 and a standard deviation of 0.13.
<table><caption>**Table 2.** Best-fitting parameters to the scaling relation between $v$ and axis ratio (equation 7). First column indicates the selection criteria (rlx = relaxed haloes, unr = unrelaxed haloes, all = full sample); second column shows the definition of virialization adopted; columns 3, 4 and 5 correspond to the b, a and c parameters; column 6 shows the product of columns 4 and 5; last column shows the standard deviation of the residuals.</caption><thead><tr><th>Select.</th><th>Δ</th><th>b</th><th>a</th><th>c</th><th>ac</th><th>σ</th></tr></thead><tbody><tr><td>rlx</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td><td>-0.226 ± 0.002</td><td>-0.295 ± 0.005</td><td>-0.27 ± 0.05</td><td>0.08 ± 0.01</td><td>0.008</td></tr><tr><td>unr</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td><td>-0.333 ± 0.002</td><td>-0.227 ± 0.004</td><td>-0.11 ± 0.04</td><td>0.03 ± 0.01</td><td>0.005</td></tr><tr><td>all</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td><td>-0.246 ± 0.002</td><td>-0.313 ± 0.005</td><td>-0.31 ± 0.05</td><td>0.09 ± 0.02</td><td>0.008</td></tr><tr><td>rlx</td><td>500ρ<sub>c</sub></td><td>-0.270 ± 0.003</td><td>-0.303 ± 0.009</td><td>-0.13 ± 0.07</td><td>0.04 ± 0.02</td><td>0.01</td></tr><tr><td>unr</td><td>500ρ<sub>c</sub></td><td>-0.301 ± 0.003</td><td>-0.329 ± 0.009</td><td>-0.11 ± 0.07</td><td>0.03 ± 0.02</td><td>0.01</td></tr><tr><td>all</td><td>500ρ<sub>c</sub></td><td>-0.279 ± 0.003</td><td>-0.310 ± 0.007</td><td>-0.21 ± 0.06</td><td>0.07 ± 0.02</td><td>0.01</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Distribution of the scaled axis ratio $\tilde{s}$ in different mass bins for relaxed haloes in the four simulated boxes. Different line types indicate the different redshifts. The black line corresponds to the best fit to a lognormal distribution (equation 10).</figcaption></figure>

In the case of axis ratios measured at 500$\rho_c$ (see Table 2), the best-fitting parameters show a lower normalization than for $\Delta_{vir}\rho_b$. In particular, for relaxed haloes $\tilde{s}_{500}(z=0, M_*) = 10^b \approx 0.54$. The redshift dependence is also lower for the axis ratios measured at 500$\rho_c$ with $c = -0.13$. The parameters of the lognormal distribution for $\tilde{s}_{500}$ are $\mu = -0.62$ and $\sigma = 0.19$ (median = 0.54 and std = 0.10, respectively).

Therefore, we can now study the distribution of axis ratios for any halo of a given mass (i.e. peak height) and redshift. First we generate one random value of the scaled axis from the lognormal distributions, either for $\tilde{s}$ or $\tilde{s}_{500}$. Then, by inverting the rescaling relation, $s = \tilde{s}v'^a$, it is possible to infer the axis ratio of the halo given its peak height $v(z, M) = v'D(z)^{-c}$. Then it can then be converted in mass for a given redshift of the halo and for a particular cosmology (see Appendix A).

## 3.2 Intermediate-to-major axis ratio

As already proposed by several authors (Schneider et al. 2012; Bonamigo et al. 2015), the minor-to-major axis ratio is better described by the rescaled quantity, $\tilde{q} = (q - s)/(1 - s)$, that eliminates the limited interval at which $q$ is defined by construction ($q \ge s$ and $q < 1$). Fig 7 shows the relation between the rescaled axis ratio $\tilde{q}$ and $s$ for relaxed haloes. We computed the probability distributions of $\tilde{q}$ in bins of different $s$. We have fitted the proba-

bility distribution of $\tilde{q}$ in each bin of $s$ to a $\beta$-function (Bonamigo et al. 2015):

$$p(\tilde{q}, \alpha, \beta) = \frac{1}{B(\alpha, \beta)} \tilde{q}^{\alpha-1} (1 - \tilde{q})^{\beta-1}, \qquad (12)$$

with two shape parameters $\alpha$ and $\beta$ and the normalization factor $1/B(\alpha, \beta)$. The mean value of the $\beta$ distribution is defined as $\mu = 1/(1 + \beta/\alpha)$. From the fitting procedure, we derived $\mu$ and $\beta$ parameters for each bin in $s$. Fig. 8 shows the dependence of these parameters on the axis ratio $s$. The solid black lines in the figure correspond to the best fit of these two parameters as follows:

$$\mu(s) = -0.002 + 0.661s, \qquad \beta(s) = 1.170s^{-1.865}. \qquad (13)$$

These results are in fairly good agreement with those reported by Bonamigo et al. (2015). Nevertheless we find a slightly higher mean value of the $\beta$ distribution for low values of $s$. Given these two equations and $p(\tilde{q}|s)$, it is possible to infer the distribution of $\tilde{q}$ for a given $s$ along with its scatter.

Additionally, we have repeated this calculation considering the axis ratios at an overdensity $500\rho_c$ (denoted as $s_{500}$ and $\tilde{q}_{500}$). The probability distributions of $\tilde{q}_{500}$ in different bins of $s_{500}$ are also well described by a $\beta$ distribution with parameters $\mu_{500}$ and $\beta_{500}$. The best fit of these two parameters are

$$\mu_{500}(s_{500}) = -0.014 + 0.676s_{500}, \qquad \beta_{500}(s_{500}) = 1.159s_{500}^{-1.792}. \qquad (14)$$

## 3.3 Triaxiality and major axis alignment

The triaxiality of an ellipsoid can be expressed as (Franx, Illingworth & de Zeeuw 1991)

$$t = \frac{a^2 - b^2}{a^2 - c^2} = \frac{1 - q^2}{1 - s^2}. \qquad (15)$$

Depending on the value of the parameter $t$, an ellipsoid is considered as oblate if $0 < t < 1/3$; triaxial if $1/3 < t < 2/3$ and prolate if $2/3 < t < 1$. In Fig. 9, we show the fraction of haloes classified in this way as a function of the halo mass (only relaxed haloes). Haloes show preferentially a prolate shape with fraction larger than the 45 per cent for haloes with masses above $10^{12} h^{-1} M_{\odot}$ at $z = 0.00$. At higher redshifts, the fraction of prolate haloes increase, while the fraction of triaxial and oblate decrease. The more massive the haloes are, the larger the fraction of prolate haloes. Consequently, the opposite is true for the triaxial and oblate fraction. In particular, for haloes at $z = 0.00$ with masses larger than $10^{14} h^{-1} M_{\odot}$ are dominated by a prolate shape (with fraction above 70 per cent), while only a few massive haloes are found to be triaxial or oblate in shape (with fraction below 20 and 5 per cent, respectively). At higher redshift, the fraction of prolate massive haloes are even larger,
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left-hand panel: distribution of $\tilde{q} = (q - s)/(1 - s)$ as a function of the virial axis ratio ($s$) for relaxed haloes. Points indicate the median values of $s$ and $\tilde{q}$, while solid (dashed) error bars represent the 25–75 per cent (1–99 per cent) probabilities of the distribution. Right-hand panel: probability distributions of $\tilde{q}$ for different values of $s$ (histograms) and the fitting function to a $\beta$ distribution (lines) for relaxed haloes. Both panels include all the redshifts analysed in this work.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Best-fitting parameters of equation (12) for different values of the virial axis ratio ($s$) for relaxed haloes (solid points). Solid lines show the best fit to the data points given by equation (13), while dashed lines correspond to the standard deviation of the residuals. For comparison, the results of Bonamigo et al. (2015) are shown as dotted lines.</figcaption></figure>

with values of about 80 per cent. This is in support of the idea that haloes with peak height $v > 1$ are expected to be undergoing a higher rate of merging processes than haloes with $v < 1$, and that the mergers happen along preferred directions (Knebe et al. 2004; Faltenbacher et al. 2005; Zentner et al. 2005).

Besides, the fraction of prolate haloes at inner radius ($R_{vir}$) are larger on average than those measured at $R_{vir}$, while the fractions of triaxial and oblate haloes are slightly lower towards the centre, independently of the halo mass. About 5–10 per cent of the haloes classified as triaxial attending to the triaxiality measured at $R_{vir}$, would be classified as prolate when looking at inner radius. In the case of oblate fraction these differences are below the 2 per cent.

Since we measured the axis ratios at two different overdensities, we examine the angle between the major axis at $\Delta_{vir}\rho_b$ and $500\rho_c$, define here as $\theta$. For the given cosmology, the radius measured at

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Fraction of oblate (red lines), triaxial (green lines) and prolate (blue lines) haloes as a function of mass for relaxed haloes according to their triaxiality parameter $t$ (equation 15). Thicker lines show to the triaxiality parameter measured at $R_{vir}$, while thinner lines show to the same parameter measured at $R_{500}$. Different line styles (solid, dashed, dashed-dotted, dotted) correspond to the different redshift analysed in this work.</figcaption></figure>

$500\rho_c$, denoted as $R_{500}$, roughly corresponds to $0.5R_{vir}$. In Fig. 10, we show the alignment of the major axis as a function of the axis ratio ($s$). We find that the major axis aligns pretty well for prolate haloes, with 99(75) per cent of them below $\theta < 40^{\circ}(11^{\circ})$ (Fig. 10). Since the most massive haloes are preferentially prolate (Fig. 9), they appear to be well aligned, but this would be simply because less massive haloes are more spherical and thus the direction of the major axis is less accurate than in the case of cluster-sized haloes. On the other hand, triaxial and oblate haloes show larger misalignments, with $\theta \gtrsim 10^{\circ}$ (on average). The angle between the major axis increases to $\theta \simeq 12^{\circ}-13^{\circ}$ (on average) for haloes $t \lesssim 0.2$. We also find a very small fraction of haloes with misalignments close to $\theta = 90^{\circ}$ and low triaxiality. These large misalignments can be due to less massive haloes with larger $q$ (thus $b \approx a$), for which the direction of the major axis is arbitrary within a plane. Finally,
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Distribution of the major axis alignment ($\theta$, in degrees) as a function of the triaxiality parameter $t$ for the whole sample at $z = 0.00$. Solid and dashed black lines correspond to the median values of the distribution for the whole and relaxed samples, respectively. Dark grey shaded region indicates the Q1 and Q3 quartiles around the median values of $\theta$, while grey shaded region indicates the 1 and 99 per cent probabilities of the distribution of $\theta$ for the whole sample only.</figcaption></figure>

relaxed haloes show, on average, slightly lower misalignments with respect to the whole sample (less than 2° differences in the case of prolate haloes).

## 3.4 Distribution of the minima axis ratio $s$

As already mentioned, triaxiality of dark matter haloes has a strong impact on the lensing cross-sections. Since massive haloes are more prolate than less massive ones, this effect is even stronger in the case of clusters of galaxies. In particular, Waizmann et al. (2012) found that the impact of the tail of the axis ratio distribution on the largest Einstein radii distribution is substantial. However, they also argued that due to the limited knowledge of the statistics of extremely small axis ratios, it is not possible to clearly define a proper choice of the cut-off (if present) until the triaxiality distributions of large halo samples (up to cluster masses) are studied in numerical simulations. The results presented in this work provide scaling relations to characterize the halo triaxiality that are derived from a huge sample of haloes (including thousands of galaxy clusters).

Since new semi-analytic lensing models could be implemented with our prescriptions on the shape of dark matter haloes, we consider important to check that extremely elongated haloes do not propagate through the sampling procedure. For this reason and in order to clarify whether a cut-off in the minor-to-major axis ratio distribution is needed or not, we present in this section the distribution of the minima of the minor-to-major axis ratio. In this analysis, we include all the haloes (both relaxed and unrelaxed) in the four simulated boxes at the four redshifts analysed. The cumulative distribution of $s$, denoted as $P(s < S)$, is shown in Fig. 11. As expected, unrelaxed haloes populate the lower tail of the distribution of axis ratios. The probability of finding a halo with $s < 0.1$ is $P(s < 0.1) = 10^{-3.20} \approx 6 \times 10^{-4}$, while for relaxed haloes this probability goes down to $P(s < 0.1) \approx 10^{-5.62}$ (about two haloes per million). In particular, we find only 13 relaxed haloes with $s < 0.1$ from a total of almost 5.5 million relaxed haloes (all redshifts included). The relaxed halo with the minimum axis ratio in our sample

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Cumulative distribution of the minima axis ratio, $P(s \le S)$, for the whole sample (solid black line) and for relaxed haloes (dashed black line), including all the redshifts analysed in this work. Red solid and dashed lines correspond to the best fitted to the whole and relaxed distributions (equations 16 and 17, respectively).</figcaption></figure>

corresponds to $s = 0.07$. For the whole sample (composed of more than 8.6 million haloes), we find 5430 haloes with $s < 0.1$.

Given the linear shape of the distribution in log-log space, we give in this section explicit relations for $P(s < S)$ as a function of $S$ that are valid for $S \lesssim 0.4$ for the whole sample,

$$ \log[P(s < S)] = (1.08 \pm 0.03) + (4.26 \pm 0.03)\log(S), \quad (16) $$

and for relaxed haloes,

$$ \log[P(s < S)] = (2.38 \pm 0.08) + (8.58 \pm 0.12)\log(S). \quad (17) $$

Although not plotted in Fig. 11, we do not find any significant dependence with redshift on the $P(s < S)$ distribution.

# 4 COMPARISON WITH PREVIOUS RESULTS

In this section, we present a comparison of the mass-shape results presented in this work with previous estimations from literature. In Fig. 12, we show the median minor-to-major axis ratios ($s$) as a function of the halo mass from the four simulated boxes at $z = 0$. We display separately the axis ratio for relaxed haloes (black filled circles) and for the whole sample (blue filled diamonds) obtained with the weighted inertia tensor, $s_w$. The dashed black line indicates the scaling relation for the unweighted inertia tensor using equation (6). In order to show the scatter of the mass-shape relation, we also plot the IQR corresponding to the median axis ratios for the relaxed sample and using the weighted inertia tensor (shaded grey region). For clarity we do not include that value for the rest of the scaling relation, although their IQR is very similar to that shown in the figure. Results from other authors are shown with thinner lines, the solid part of the curves represents the mass limits studied in each case.

As it can be seen in Fig. 12, there is an overall agreement in the minor-to-major axis ratio ($s$) as a function of halo mass, with more massive haloes showing less spherical shapes. The scatter in the axis ratio relations, which is constrained within ~10 per cent, is mainly due to the differences in the method of measuring the shapes (Section 3) and the definition of halo masses and the axis radius (i.e. the overdensity at which it is measured). For a reliable
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of the mass-shape scaling relation with previous studies at $z=0$. Black and dark-blue lines show the results of this work using the weighted inertia tensor at $R_{vir}$ for relaxed haloes and the whole sample, respectively. The grey shaded region indicates the IQR associated with the black line. The dashed black line shows the axis ratio derived using the unweighted inertia tensor at $R_{vir}$ (only for hMD simulation), while the dashed dark-blue line corresponds to the axis ratio measured at $0.3R_{vir}$ for the whole sample (weighted inertia tensor). Black filled circles and blue filled diamonds correspond to the median axis ratio derived at $R_{vir}$ using the weighted inertia tensor for relaxed haloes and the whole sample, respectively. Thin solid lines show the results from other studies in the mass range analysed in each case, while extrapolation outside the mass range is denoted by thin dashed lines.</figcaption></figure>

comparison, we have converted the results from other authors to redshift $z=0$ and the cosmology assumed in this work. Therefore, the differences shown in Fig. 12 cannot be attributed to differences in the cosmology assumed in each work. Nevertheless, we would like to remark that the MultiDark suite constitutes a coherent and homogeneous data set of cosmological simulations (performed with L-GADGET-2 code for identical cosmological parameters), for which all the haloes are identified using the ROCKSTAR halo finder. In particular, the HMD simulation ($4h^{-1}$ Gpc cube side) allows us to study in detail the shape of simulated dark matter haloes in the cluster mass range with unprecedented statistics and resolution (more than $8 \times 10^5$ haloes at $z=0$ with more than 3000 particles per halo). As we described in Section 2, given that massive haloes are more aspherical, larger simulated volumes show (on average) larger unrelaxed fraction than smaller ones (~40 per cent for the HMD simulation, while ~20 per cent for the SMD simulation). The fraction of unrelaxed haloes found in the simulations have a significant impact on the mass-shape scaling relations when all the haloes are included in the analysis (i.e. not selected by their relaxation state), as we have shown in Section 3.1. On the other hand, small simulated volumes will not be strongly affected by this selection criteria because of the small fraction of unrelaxed haloes found in them. In order to facilitate the comparison we list the information and details of previous studies on the shapes of dark matter haloes in Table 3.

The definition of the halo masses can be obviously a source of the scatter in the mass-shape scaling relations. In this work we used the definition given by the spherical collapse model (i.e. the mass inside the radius encompassing a given density contrast), then the procedure used to compute the axis ratio iteratively keeps the major axis equal to the original sphere radius. According to this

description the ellipsoid is, by definition, inside the sphere and therefore its mass is always smaller than the mass of the initial sphere. Nevertheless, Muñoz-Cuartas et al. (2011) checked that these differences are below 20 per cent. Alternatively, in Bonamigo et al. (2015) and Despali et al. (2014) the halo mass is defined as the mass within the ellipsoid that encloses a given overdensity. Therefore, the best-fitting ellipsoid always contains a mass larger or equal to that of the original sphere. Despali et al. (2013) found that elliptical masses are on average about 3 per cent larger than spherical ones for cluster-size haloes and 5 per cent larger for lower mass haloes. Taking all this into account, if spherical masses showed in Fig. 12 would be converted into elliptical masses, the results slightly shifted to the right (i.e. towards higher masses), leading to a better agreement with Bonamigo et al. (2015).

Additionally, since ellipsoid used in this work to compute the axis ratio is inside the original sphere, we measure an inner shell of the halo mass distribution, which we know may be on average more aspherical towards the halo centre. For the same reason, the method used in Bonamigo et al. (2015) to determine the axis ratios could lead to more spherical shapes as it extends the computation to radii beyond the original sphere.

Given that the peak height is defined in terms of the critical overdensity of the spherical collapse model ($\delta_c$), we consider more appropriate to derive the peak height from the spherical halo mass. Besides, the mass function of dark matter haloes is generally described according to the spherical overdensity criteria, so that the virial halo masses obtained from the mass function are spherical rather than elliptical. For all these reasons, we consider the spherical definition of the virial halo mass more useful for semi-analytical lensing studies that infer the lensing
<table><caption>**Table 3.** Details of the simulations and the methods used to derive the mass-shape scaling relations from previous studies (shown in Fig. 12). First column indicates the referenced paper; second column shows the ΛCDM cosmological model considered; third column shows the size limits of the simulated boxes analysed (in $h^{-1}$ Gpc); column 4 corresponds mass limits per particle, $m_p$ (in units $10^{10} h^{-1}$ M$_\odot$); column 5 shows the minimum number of particles per halo ($N_p$); column 6 indicates the selection criteria (rlx = relaxed haloes, all = full sample); column 7 indicates the inertia tensor used (w = weighted, unw = unweighted); column 8 corresponds to the overdensity at which the axis ratios are measured.</caption><thead><tr><th>Reference</th><th>Cosmo.<sup>a</sup></th><th>Box</th><th>m<sub>p</sub></th><th>N<sub>p</sub></th><th>Select.</th><th>Tensor</th><th>Δ</th></tr></thead><tbody><tr><td>This work</td><td>Planck</td><td>[0.4–4.0]</td><td>[0.01–7.9]</td><td>3000</td><td>rlx, all</td><td>w, unw</td><td>Δ<sub>vir</sub>ρ<sub>b</sub>, 500ρ<sub>c</sub></td></tr><tr><td>Bonamigo et al. (2015)</td><td>WMAP1, Planck</td><td>[0.06–3.0]</td><td>[0.002–0.6]</td><td>1000</td><td>rlx</td><td>unw</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td></tr><tr><td>Despali et al. (2014)</td><td>Planck</td><td>[0.1–2.0]</td><td>[0.06–63]</td><td>200</td><td>all</td><td>unw</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td></tr><tr><td>Schneider et al. (2012)</td><td>WMAP3</td><td>[0.1–0.5]</td><td>[< 10<sup>-3</sup>–0.09]</td><td>500</td><td>rlx</td><td>w</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td></tr><tr><td>Muñoz-Cuartas et al. (2011)</td><td>WMAP5</td><td>[0.01–0.15]</td><td>[0.001–1.13]</td><td>4000</td><td>rlx</td><td>unw</td><td>Δ<sub>vir</sub>ρ<sub>b</sub></td></tr><tr><td>Allgood et al. (2006)</td><td>ΛCDM</td><td>[0.08–0.2]</td><td>[0.03–4.0]</td><td>7000</td><td>all</td><td>w</td><td>Δ<sub>vir</sub>ρ<sub>b</sub><sup>b</sup></td></tr><tr><td>Jing & Suto (2002)<sup>c</sup></td><td>ΛCDM</td><td>0.1</td><td>0.06</td><td>10 000</td><td>all</td><td>unw</td><td>2500ρ<sub>c</sub></td></tr></tbody></table>

properties of dark matter haloes sampled from a particular mass function.

The various definitions of the halo centre can also introduced some scatter in the mass-shape scaling relations. In this work, we define the halo centres as the densest point in phase space, which are not recalculated within the iterative procedure. Therefore, if there is an offset between the centre of mass and the density peak, this contributes to the inferred shape. Since unrelaxed clusters tend to have larger offsets, this effect is stronger for studies including all haloes in their analysis. However, we do not expect the mass-shape scaling relations to be very sensitive to the choice of the halo centres if they are defined as the densest point or the most bounded particle. In other words, even in the case of unrelaxed haloes, the position of the densest point (or most bounded particle) will remain almost constant, while the centre of mass will change in position.

Jing & Suto (2002) measured the axis ratios using particles within the isodensity surface of $2500\rho_c$, which roughly corresponds to $0.3R_{vir}$. Therefore, for a better comparison, these results should be compared with our results for $s_{500}$ (dashed blue line). There is a good agreement on average for low-mass haloes, but given that our relation is steeper, more massive haloes are less aspherical than found by Jing & Suto (2002). It should be notice that Jing & Suto (2002) derived the shape properties of dark matter haloes extracted from five realizations of $100 h^{-1}$ Mpc box size and from 12 haloes with $10^6$ particles within the virial radius, and therefore they studied a narrower mass range.

In the case of Allgood et al. (2006), where analytical expressions on the radial dependence of the axis ratios are explicitly given, we convert their axis ratio relation from $0.3R_{vir}$ to $R_{vir}$. Although we find a steeper slope than Allgood et al. (2006), the results are in agreement within errors for low-mass haloes when unrelaxed haloes are also included in the analysis (solid blue line).

Both Allgood et al. (2006) and Schneider et al. (2012) used the method based on the weighted inertia tensor. In particular, Schneider et al. (2012) derived the axis ratios measured at the virial radius, which allows a meaningful comparison with our results for relaxed haloes (solid black line). On the other hand, Muñoz-Cuartas et al. (2011) and Bonamigo et al. (2015) studied the shapes using the unweighted inertia tensor. Although the results are in good agreement when comparing with our results for relaxed haloes using the unweighted inertia tensor (dashed black line), we obtain a lower normalization. In the case of Muñoz-Cuartas et al. (2011), shapes are

obtained from a small sample of haloes extracted from small simulated volumes, less than 12 000 with masses above $5 \times 10^{12} h^{-1}$ M$_\odot$ (i.e. $N_p > 500$). For the mass-shape scaling relation, the authors restricted the analysis to an even smaller sample, including only haloes with $N_p > 4000$.

On the other hand, for massive clusters we find axis ratios lower than those presented by Bonamigo et al. (2015) in about 5–8 per cent. It should be noticed that most of the haloes analysed in Bonamigo et al. (2015) are extracted from the Millennium XXL simulation, which is performed with 1-year Wilkinson Microwave Anisotropy Probe (WMAP1) cosmological parameters ($\Omega_M = 0.25$, $\Omega_b = 0.045$, $\Omega_\Lambda = 0.75$, $\sigma_8 = 0.9$, $h = 0.73$). Nevertheless, the authors scaled their results in order to account for the differences in the cosmological model. Despite the fact of the different mass resolution and box volume in the simulations, the main difference between the results of this work and the results presented by Bonamigo et al. (2015) is the definition of the halo mass, the method to determine the axis ratios (weighted and unweighted, respectively) and the criteria considered to select relaxed haloes. As we mentioned above, since massive haloes are more aspherical than lower mass haloes, we expect the differences between the spherical and the elliptical masses would slightly shift our mass-shape scaling relations to the right (i.e. towards higher mass), and hence reducing the observed scatter. Although we use the scaling relation given by equation (6) to derive the unweighted shape-mass scaling relation (dashed black line), there is a better agreement with Bonamigo et al. (2015). It is important to notice the different criteria adopted to select relaxed haloes: Bonamigo et al. (2015) fixed the offset between the most bounded particle and the centre of mass of the particles enclosed by the ellipsoid to be less than 5 per cent of their virial radius, while in this work, we assume a more flexible value of 7 per cent. Therefore, the relaxed sample of haloes analysed in Bonamigo et al. (2015) included less distorted haloes than those incorporated in our analysis. For that reason, we expect that the shape-mass relation of Bonamigo et al. (2015) lies above our results (i.e. higher values of $s$) since unrelaxed haloes tend to be more elongated. In particular, this effect is more important for more massive haloes given that unrelaxed fraction are higher for higher halo masses, as it can be clearly seen in Fig. 2.

Finally, the findings of Despali et al. (2014) are obtained from a set of simulations performed with almost the same cosmological parameters as the simulations used in this work and based on the

<footnote><sup>a</sup>Cosmological parameters ($\Omega_M$, $\Omega_\Lambda$, $h$, $\sigma_8$): WMAP1 (0.25, 0.75, 0.73, 0.9); WMAP3 (0.25, 0.75, 0.73, 0.9); WMAP5 (0.258, 0.742, 0.72, 0.796); Planck (0.307, 0.693, 0.678, 0.829); ΛCDM (0.3, 0.7, 0.7, 0.9).</footnote>

<footnote><sup>b</sup>The halo shapes were determined at five different fraction of their virial radius in order to examine the radial dependence of $s$ with halo mass. We use this radial-shape relation to convert from $0.3R_{vir}$ to $R_{vir}$.</footnote>

<footnote><sup>c</sup>The analysis also includes 12 high-resolution haloes with more than $10^6$ particles within the virial radius.</footnote>
unweighted inertia tensor. We find a very good agreement with their predictions for relaxed haloes, even though they included in the analysis all the haloes regardless of their state of relaxation. Nevertheless, given that their simulated volumes are smaller than the ones analysed in this work (and therefore their fraction of unrelaxed haloes are also smaller), we do not expect a strong impact on the distribution of the median axis ratios when including the unrelaxed haloes. Again in this case, the scatter in the scaling relation can be due to the definition of the halo mass (either spherical or elliptical) and the used of the unweighted inertia tensor method to measure the axis ratios.

It is clear that there is a general agreement between the studies performed over a large range in halo mass, such as the results presented in this work and in Despali et al. (2014) and Bonamigo et al. (2015). Besides, these mass-shape scaling relations are obtained in terms of the peak height $v$ that leads to steeper relations than those found in rest of the studies included in the comparison (which are derived in terms of the characteristic mass, $M_*$). Despite the different simulations and the various methodologies used to determine the halo shapes in the mentioned studies, the mass-shape relations show a similar slope but with a different ($\lesssim 10$ per cent) normalization.

## 5 CONCLUSIONS

We have presented a detailed analysis on the triaxiality of dark matter haloes in four different computational volumes from the MultiDark Planck simulation suite (Klypin et al. 2016). In particular, the HMD simulation allows us to characterize the shapes of cluster mass haloes with unprecedented statistics (almost 400 000 haloes with masses above $2 \times 10^{14} h^{-1} M_{\odot}$). Using the three simulations with smaller box volumes, we extend our analysis to lower masses down to $3 \times 10^{11} h^{-1} M_{\odot}$. The main results of our analysis can be summarized as follows.

(i) Dark matter haloes are well described by a triaxial model. In particular, due to interactions with other clusters or groups of galaxies, more massive haloes are less spherical than less massive ones (with a tendency of being prolate). Besides, massive haloes are also less spherical at high redshift essentially because they are affected by the direction of the last major merger and/or by the presence of filaments around them, and because massive haloes were formed later on.

(ii) Since the modelling of non-spherical objects is not simple and highly elongated haloes are mostly due to unrelaxed, we analyse separately relaxed and unrelaxed haloes based on the criteria proposed by Klypin et al. (2016).

(iii) The iterative weighted and unweighted inertia tensors give similar results, while the inertia tensor in a spherical window is considerably higher than the formers.

(iv) Minor-to-major axis ratio can be expressed by a universal function in terms of the peak height $v$ (equation 7). The peak height is related to the halo mass and incorporates the dependence on the cosmological parameters.

(v) Rescaled minor-to-major axis ratio obtained from this functional form is well described by a lognormal distribution within the whole mass range considered.

(vi) Dark matter haloes are less spherical at inner radius, with axis ratios measured at $R_{500}$ lower than those measured at $R_{vir}$ on about 9 per cent.

(vii) The conditional intermediate-to-major axis ratio distribution ($p(q|s)$) is well fitted by a beta distribution, which depends only on the minor-to-major axis ratio (but not on the halo mass).

(viii) For $z=0$ haloes with masses larger than $10^{14} h^{-1} M_{\odot}$ are dominated by prolate shapes (with fraction above 70 per cent), while only a few massive haloes are found to be triaxial or oblate in shape (with fraction below 20 and 5 per cent, respectively). At a higher redshift, the fraction of prolate massive haloes are even larger, with values of about 80 per cent.

(ix) The angle between the major axis measured at $R_{500}$ and $R_{vir}$ is below $30^{\circ}$ for the 75 per cent of the triaxial and prolate haloes analysed, meaning that the direction of the major axis at $0.5R_{vir}$ aligns pretty well with the direction of the major axis at $R_{vir}$.

(x) Extremely low axis ratio are mostly due to massive and unrelaxed haloes; overall, the probability of finding a halo with $s < 0.1$ is $P(s < 0.1) = 10^{-3.20} \approx 6 \times 10^{-4}$. When the analysis is restricted only to relaxed haloes these probabilities decrease down to $P(s < 0.1) \approx 10^{-5.62}$ (about two haloes per million).

Finally, in Section 4, we have compared our results with previous studies. Although there is an overall agreement with previous studies, we find a $\sim 10$ per cent scatter mainly due to the different methods used to determine the axis ratios (i.e. inertia tensor, the radius at which the axis are measured) and the definition of the halo mass (spherical or elliptical). Additionally, in studies that include unrelaxed haloes, the fraction of unrelaxed haloes found in the simulated volumes can considerably affect the mass-shape scaling relation, leading to a lower normalization value of the axis ratios. Nevertheless, the best agreement is found for studies in which a large halo mass range is considered.

Using a suite of high-resolution simulations, the shape of dark matter haloes has been accurately calculated over five orders of magnitude in halo mass. These precise estimates are in general agreement with previous work done for smaller mass ranges with differences due to the different estimation of the shapes.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE; http://www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at the Leibniz Supercomputing Centre (LRZ; http://www.lrz.de) in Munich. We would like to thank Peter Behroozi for his helpful comments on this paper and for making the ROCKSTAR halo finder software publicly available. We would like to acknowledge the anonymous referee for detailed comments that have contributed to improve the quality of this paper. The MultiDark data base was developed in cooperation with the Spanish MultiDark Consolidator Project CSD2009-00064. The CosmoSim data base used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). JV-F acknowledges financial support from MINECO (Spain) under the FPI grant Ref. BES-2010-039705 and travel grants Ref. EEBB-I-12-05072 and Ref. EEBB-I-14-08162. GY acknowledges financial support from MINECO (Spain) under research grants AYA2012-31101 and AYA2015-63810-P (MINECO/FEDER, EU).

## 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., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046

Applegate D. E. et al., 2016, MNRAS, 457, 1522

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

Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

Bartelmann M., Huss A., Colberg J. M., Jenkins A., Pearce F. R., 1998, A&A, 330, 1

Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109

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

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

Bonamigo M., Despali G., Limousin M., Angulo R., Giocoli C., Soucail G., 2015, MNRAS, 449, 3171

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

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

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

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

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

Faltenbacher A., Allgood B., Gottlöber S., Yepes G., Hoffman Y., 2005, MNRAS, 362, 1099

Ford J., Van Waerbeke L., Milkeraitis M., Laigle C., Hildebrandt H., Erben T., Heymans C., Hoekstra H., 2015, MNRAS, 447, 1304

Franx M., Illingworth G., de Zeeuw T., 1991, ApJ, 383, 112

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., Meneghetti M., Bartelmann M., Moscardini L., Boldrin M., 2012a, MNRAS, 421, 3343

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

Groener A. M., Goldberg D. M., 2014, ApJ, 795, 153

Hoekstra H., Herbonnet R., Muzzin A., Babul A., Mahdavi A., Viola M., Cacciato M., 2015, MNRAS, 449, 685

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

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

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

Kawahara H., 2010, ApJ, 719, 1926

Kettula K. et al., 2015, MNRAS, 451, 1460

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

Knebe A., Gill S. P. D., Gibson B. K., Lewis G. F., Ibata R. A., Dopita M. A., 2004, ApJ, 603, 7

Lau E. T., Nagai D., Nelson K., 2013, ApJ, 777, 151

Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., Navarro J. F., Cole S., Frenk C. S., 2016, MNRAS, 460, 1214

Medezinski E., Umetsu K., Okabe N., Nonino M., Molnar S., Massey R., Dupke R., Merten J., 2016, ApJ, 817, 24

Meneghetti M., Bolzonella M., Bartelmann M., Moscardini L., Tormen G., 2000, MNRAS, 314, 338

Meneghetti M., Yoshida N., Bartelmann M., Moscardini L., Springel V., Tormen G., White S. D. M., 2001, MNRAS, 325, 435

Mo W. et al., 2016, ApJ, 818, L25

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

Murray S. G., Power C., Robotham A. S. G., 2013, Astron. Comput., 3, 23

Natarajan P., Refregier A., 2000, ApJ, 538, L113

Nelder J. A., Mead R., 1965, Comput. J., 7, 308

Oguri M., Takada M., Umetsu K., Broadhurst T., 2005, ApJ, 632, 841

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

Peebles P. J. E., 1969, ApJ, 155, 393

Planck Collaboration XVI, 2014, A&A, 571, A16

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

Redlich M., Bartelmann M., Waizmann J.-C., Fedeli C., 2012, A&A, 547, A66

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., 5, 30

Schrabback T. et al., 2015, MNRAS, 454, 1432

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

Springel V., 2005, MNRAS, 364, 1105

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Waizmann J.-C., Redlich M., Bartelmann M., 2012, A&A, 547, A67

Waizmann J.-C., Redlich M., Meneghetti M., Bartelmann M., 2014, A&A, 565, A28

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

Wojtak R., 2013, A&A, 559, A89

Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H., 2005, ApJ, 624, 505

## APPENDIX A: DENSITY PEAK HEIGHT

The density peak height $v$ for a virialized halo with mass $M$ at redshift $z$ for a given cosmological model is given by

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

where $\delta_c(z)$ is the critical overdensity of the spherical collapse model (i.e. the initial density required for a fluctuation to collapse at redshift $z$) and $\sigma(M)$ 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$).

The critical overdensity can be rewritten as the collapse overdensity at present time ($\delta_c(z=0) = \delta_c \approx 1.69$) rescaled to a given time as follows: $\delta_c(z) = \delta_c/D(z)$, with $D(z)$ being the linear growth rate of a density fluctuation normalized to unity at the present time.

For the computation of $D(z)$ we use the approximation given by Carroll, Press & Turner (1992), which can be expressed as $D(z) \propto g(z)/(1+z)$ with

$$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 (A2)$$

The variance in the initial density field, $\sigma(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 (A3)$$

where $\tilde{W}$ is the Fourier transform of the top-hat window function. For a given cosmology, both the power spectrum ($P(k)$) and the variance in the initial density field ($\sigma(M)$) can be computed using a software like HMF_CALC (Murray, Power & Robotham 2013).

On one hand, the overdensity $\delta_c$, which depends only on redshift and not on the mass, increases with redshift due to $D(z)$; while, on the other hand, the $\sigma(M)$ depends on the mass but not on redshift (since $\sigma(M) = \sigma(M, z=0)$ by definition). Therefore, to compute the peak height $v$ for a halo of mass $M$ only the linear growth rate $D(z)$ and the initial power spectrum $P(k)$ are needed.

Additionally, we proposed a relation between the halo mass $M$ and the peak height $v$, for the cosmological parameters given in Section 2. At present time, the relation between $v(z=0) = v_0$ and the halo mass (in log-log space) can be expressed as a third-order polynomial:

$$\log v_0 = n_0 + n_1 x + n_2 x^2 + n_3 x^3, \qquad (A4)$$
<figure><img src="image_16.png" /></figure>

**Figure A1.** Peak height as a function of the virial halo mass for the Planck cosmology. Solid lines correspond with the results obtained in this work, while dashed lines show the findings of Ludlow et al. (2016). Blue and red lines indicate the corresponding redshift $z = 0$ and 1, respectively. The dashed green lines correspond to the best fit with equation (A4). The vertical lines indicate the characteristic masses at each redshift given by $v(z, M_*) = 1$.

with $x = \log M$. As can be clearly seen in Fig. A1, this relation gives an excellent fit within the whole mass range studied here with best-fitting parameters: $n_0 = -2.63 \pm 0.04$, $n_1 = 0.504 \pm 0.009$, $n_2 = -0.0434 \pm 0.0007$ and $n_3 = 0.00158 \pm 0.00002$. This result is in good agreement with the functional form proposed by Ludlow et al. (2016) for a similar set of cosmological parameters.

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