# The mass dependence of dark matter halo alignments with large-scale structure

Davide Piras, Benjamin Joachimi, Björn Malte Schäfer, Mario Bonamigo,

Stefan Hilbert, Edo van Uitert

► To cite this version:

Davide Piras, Benjamin Joachimi, Björn Malte Schäfer, Mario Bonamigo, Stefan Hilbert, et al.. The mass dependence of dark matter halo alignments with large-scale structure. Monthly Notices of the Royal Astronomical Society, 2018, 474 (1), pp.1165-1175. 10.1093/mnras/stx2846. hal-01669619

HAL Id: hal-01669619

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

Submitted on 26 Mar 2025

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

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

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

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

# The mass dependence of dark matter halo alignments with large-scale structure

Davide Piras,<fnref n="1" />,<fnref n="2" />* Benjamin Joachimi,<fnref n="1" /> Björn Malte Schäfer,<fnref n="3" /> Mario Bonamigo,<fnref n="4" />,<fnref n="5" /> Stefan Hilbert<fnref n="6" />,<fnref n="7" /> and Edo van Uitert<fnref n="1" />

Accepted 2017 October 31. Received 2017 October 6; in original form 2017 July 20

## ABSTRACT

Tidal gravitational forces can modify the shape of galaxies and clusters of galaxies, thus correlating their orientation with the surrounding matter density field. We study the dependence of this phenomenon, known as intrinsic alignment (IA), on the mass of the dark matter haloes that host these bright structures, analysing the Millennium and Millennium-XXL N-body simulations. We closely follow the observational approach, measuring the halo position–halo shape alignment and subsequently dividing out the dependence on halo bias. We derive a theoretical scaling of the IA amplitude with mass in a dark matter universe, and predict a power law with slope $\beta_M$ in the range 1/3 to 1/2, depending on mass scale. We find that the simulation data agree with each other and with the theoretical prediction remarkably well over three orders of magnitude in mass, with the joint analysis yielding an estimate of $\beta_M = 0.36_{-0.01}^{+0.01}$. This result does not depend on redshift or on the details of the halo shape measurement. The analysis is repeated on observational data, obtaining a significantly higher value, $\beta_M = 0.56_{-0.05}^{+0.05}$. There are also small but significant deviations from our simple model in the simulation signals at both the high- and low-mass end. We discuss possible reasons for these discrepancies, and argue that they can be attributed to physical processes not captured in the model or in the dark matter-only simulations.

**Key words:** galaxies: haloes – dark matter – large-scale structure of Universe.

## 1 INTRODUCTION

In a statistically isotropic universe, galaxy images have no globally preferred orientation; this, however, does not preclude local correlations of galaxy shapes, with each other and with the large-scale structure. Several mechanisms have been proposed to explain this, such as the accretion of new material along favoured directions, and the effect of gravitational tidal fields from the surrounding dark matter distribution. In this latter picture, in particular, the shape of luminous structures is affected by tidal interactions by the surrounding dark matter halo, and, consequently, galaxies and clusters also align (Kiessling et al. 2015).

This phenomenon is known as intrinsic alignment (henceforth IA; see Joachimi et al. 2015; Kiessling et al. 2015; Kirk et al. 2015;

Troxel & Ishak 2015 for recent reviews, Schmidt, Chisari & Dvorkin 2015 for its impact in the early-Universe cosmology and L'Huillier et al. 2017 for its role in the study of different modified gravity and dark energy models): besides carrying information about galaxy formation, IA has an impact on the measurements of the cosmological weak lensing effect, the induced correlations of galaxy image shapes due to gravitational lensing by the mass distribution along the line of sight. IA must therefore be taken into consideration when analysing lensing surveys, such as LSST (LSST Science Collaboration 2009) and Euclid (Laureijs et al. 2011). Heavens, Refregier & Heymans (2000) and Croft & Metzler (2000) first predicted the non-negligible contamination of the weak lensing signal due to correlations in the intrinsic shapes of galaxies, and a number of works afterwards (e.g. Heymans et al. 2006; Semboloni et al. 2008; Kirk et al. 2012) confirmed that this effect must be accounted for in order not to bias the cosmological inference.

<footnote n="1">¹Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="2">²Dipartimento di Fisica ‘G. Galilei’, Università di Padova, via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote n="3">³Astronomisches Recheninstitut, Zentrum für Astronomie der Universität Heidelberg, Philosophenweg 12, D-69120 Heidelberg, Germany</footnote>

<footnote n="4">⁴Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark</footnote>

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

<footnote n="6">⁶Exzellenzcluster Universe, Boltzmannstr 2, D-85748 Garching, Germany</footnote>

<footnote n="7">⁷Ludwig-Maximilians-Universität, Universitäts-Sternwarte, Scheinerstr 1, D-81679 München, Germany</footnote>

<footnote>* E-mail: d.piras@ucl.ac.uk</footnote>
The IA signal has been measured in N-body simulations (Heymans et al. 2006; Kuhlen, Diemand & Madau 2007; Lee et al. 2008; Schneider, Frenk & Cole 2012; Joachimi et al. 2013a,b), in hydrodynamical simulations (Chisari et al. 2015; Codis et al. 2015; Velliscig et al. 2015b; Tenneti, Mandelbaum & Di Matteo 2016; Hilbert et al. 2017) and in imaging data (Hao et al. 2011; Joachimi et al. 2011; Li et al. 2013; Singh, Mandelbaum & More 2015; van Uitert & Joachimi 2017). Results from hydrodynamical simulations and observations, in particular, claimed that massive red galaxies point towards matter overdensities, while blue galaxies do not show any clear sign of alignment (Hirata et al. 2007; Mandelbaum et al. 2011).

The amplitude of the IA signal has been found to increase with mass in observations. Hao et al. (2011), for example, studied a large sample of galaxy clusters from the SDSS DR7 finding a dependence of the alignment with the mass of the brightest cluster galaxy, and an increasing IA amplitude with luminosity, or the corresponding mass, has been identified for luminous red galaxies (LRGs, Joachimi et al. 2011; Singh et al. 2015) and for galaxy clusters (van Uitert & Joachimi 2017).

Simulations agree with this picture: using a 512³-particle N-body simulation, Jing (2002), for example, found that the alignment increases with mass over three orders of magnitude, up to about 10¹³ M☉. Moreover, Lee et al. (2008) used data from the Millennium simulation, claiming stronger correlations with higher mass over two mass bins around 10¹² M☉, and Joachimi et al. (2013b), using the same set of data, found an increasing trend over the same mass range only for early-type galaxies (whose shapes were assumed to follow those of the underlying haloes), while no dependence on luminosity or particular trend with mass or luminosity for late-type galaxies (whose orientation was determined by the halo spin) was identified.

These results motivate us to search for a universal relation between the alignment strength and the mass of dark matter haloes. This could lend support to an IA mechanism that successfully explains the trends; moreover, we need to reduce the degrees of freedom in modelling the IA signal to obtain tighter cosmology constraints from lensing observations. To achieve this, we delve into the dependence of the amplitude of the IA signal on the mass of the halo, which is postulated to be the main driver of the shape and orientation of the hosted galaxies (e.g. Joachimi et al. 2015). We first derive the expected scaling from the theory, and then test our predictions using data from two N-body simulations, mimicking the observational approach in order to be able to straightforwardly compare our results with real data.

We give details about our theoretical model and derive the expected scaling of the IA with halo mass in the tidal alignment paradigm in Section 2. We then present the simulations that we use (Section 3.1), and how we define the shapes of the dark matter haloes they contain (Section 3.2). In Section 4.1, we explain how we measure the IA signal, and then (Section 4.2) we describe our mass-dependent IA model. We finally show our results and compare them with our theoretical predictions (Section 5.1) and real data (Section 5.2).

## 2 THEORETICAL BACKGROUND

The physical picture of tidal interaction of a self-gravitating system that is in virial equilibrium, with a velocity dispersion $\sigma^2$, such as an elliptical galaxy or a relaxed galaxy cluster, is a distortion of the system's gravitational potential through tidal gravitational forces. The particles of the system remain in virial equilibrium and fill

up the distorted potential along an isocentre of the gravitational potential, which results in a change in the shape of the system. This shape modification reflects the magnitude and the orientation of the tidal gravitational fields, and the magnitude of the change in shape depends on how tightly the system is bound.

In isolated virialized systems the Jeans equation applies (see e.g. Binney & Tremaine 2008):

$$ \frac{1}{\rho} \frac{\partial}{\partial r} (\rho \sigma^2) + \frac{2}{r} \beta_{\text{aniso}} \sigma^2 = - \frac{\partial \Phi}{\partial r}, \quad (1) $$

with $\rho$ the particle density, $r$ the distance, $\Phi$ the gravitational potential and $\beta_{\text{aniso}}$ the anisotropy parameter, which we set to $\beta_{\text{aniso}} = 0$ because we aim to derive only the scaling behaviour of the alignment amplitude. In the case of vanishing anisotropy and a constant velocity dispersion, the Jeans equation can be solved to yield an exponential dependence $\rho \propto \exp(-\Phi/\sigma^2)$ between the density and the gravitational potential.

Gravitational tidal fields generated by the ambient large-scale structure distort the gravitational potential $\Phi$, and we will work in the limit that the distortion is well described by a second-order Taylor-expansion of the potential relative to the centre of the galaxy at $r_0 = 0$. Therefore, the potential relevant for the motion of test particles is given by

$$ \Phi(\mathbf{r}) \rightarrow \Phi(\mathbf{r}) + \frac{1}{2} \sum_{i,j=1}^{3} \frac{\partial^2 \Phi(r_0)}{\partial r_i \partial r_j} r_i r_j, \quad (2) $$

with $i, j$ indicating the spatial dimensions. Consequently, the density of particles changes according to

$$ \rho = \rho_0 \left( 1 - \frac{1}{2\sigma^2} \sum_{i,j=1}^{3} \frac{\partial^2 \Phi(r_0)}{\partial r_i \partial r_j} r_i r_j \right), \quad (3) $$

with $\rho_0 \propto \exp(-\Phi(r)/\sigma^2)$, under the assumption of a weak distortion such that a Taylor-expansion of the exponential to first order is sufficient.

The projected ellipticity is calculated through the tensor of second brightness moments (Bartelmann & Schneider 2001), which we define as

$$ q_{ij} = \int d^2 x \rho_0(\mathbf{x}) x_i x_j, \quad (4) $$

where the integral is calculated over the plane of the projected sky, with coordinates $x_i, x_j$ and $i, j \in \{1, 2\}$. If we consider the distortion of the density described in equation (3), we obtain

$$ \begin{aligned} \tilde{q}_{ij} &= \int d^2 x \rho_0(\mathbf{x}) \left( 1 - \sum_{a,b=1}^{2} \frac{\partial^2 \Phi(\mathbf{x}_0)}{2\sigma^2} x_a x_b \right) x_i x_j \\ &\approx \int d^2 x \rho_0(\mathbf{x}) x_i x_j - \sum_{a,b=1}^{2} \frac{\partial^2 \Phi}{2\sigma^2} \int d^2 x \rho_0(\mathbf{x}) x_a x_b x_i x_j \\ &\equiv q_{ij} + \psi_{ij}, \end{aligned} \quad (5) $$

where we switch to a 2D sum since we assume that most relevant are the components of the tidal shear perpendicular to the line of sight, and where the equivalence is not exact due to the approximation that $\partial^2 \Phi(\mathbf{x}_0)/\sigma^2$ is constant and can be thus taken out of the integral. From this latter equation, we can see that because of the tidal shear fields the second moments of the brightness distribution get a correction $\psi_{ij}$. This term is small in the limit of weak tidal fields, which is characterized by $R^2 \partial^2 \Phi/\sigma^2 \ll 1$ with R being the size of the halo and providing a bound for $r$.
If we then define the complex ellipticity as in Bartelmann & Schneider (2001):

$$ \varepsilon = \frac{q_{11} - q_{22}}{q_{11} + q_{22} + Q} + 2i \frac{q_{12}}{q_{11} + q_{22} + Q}, \qquad (6) $$

with $Q = 2\sqrt{q_{11}q_{22} - q_{12}^2}$, and consider the correction obtained in equation (5), we can write

$$ \tilde{\varepsilon} = \varepsilon + \psi_{\varepsilon}, \qquad (7) $$

with $\varepsilon$ the unperturbed shape as in equation (6), which we assume to be randomly oriented, and

$$ \psi_{\varepsilon} = \frac{\psi_{11} - \psi_{22}}{q_{11} + q_{22} + Q} + 2i \frac{\psi_{12}}{q_{11} + q_{22} + Q}, \qquad (8) $$

assuming that $\psi_{ij}$ is a small correction.

Therefore, a measurement of the ellipticity of a tidally distorted halo through the second moments of the projected density yields a proportionality to $R^2\partial^2\Phi/\sigma^2$, in accordance with the linear alignment model for elliptical galaxies, where tidal interaction imprints a shape distortion that causes correlations between shapes of neighbouring objects, and the tidally induced ellipticity is proportional to the magnitude of the tidal fields (Hirata & Seljak 2004; Blazek, McQuinn & Seljak 2011; Blazek, Vlah & Seljak 2015). It is worth pointing out that the tidal effect on galaxies is determined by the velocity scale $\sigma^2$, whereas in the case of gravitational lensing the tidal shear is measured in units of the velocity scale $c^2$.

In the case of virialized systems it is possible to relate the velocity dispersion $\sigma^2$ with the size of the object R: the virial relationship $\sigma^2 = GM/R$ assumes a proportionality between the specific kinetic energy, $\sigma^2/2$, and the magnitude of the specific potential energy, $GM/R$. With the scaling $M \propto R^3$ one expects $\sigma^2 \propto M^{2/3}$, such that $R^2/\sigma^2$ is constant. Therefore, any scaling of the ellipticity with mass is entirely due to the dependence of tidal gravitational fields on the mass scale, and more massive systems are subjected to stronger tidal interactions because of the stronger fluctuations in the tidal fields that they are experiencing.

The variance of tidal shear fields can be inferred from the variance of the matter density by the Poisson equation, $\Delta\Phi = 3\Omega_m/(2\chi_H^2)\delta$, with $\Omega_m$ the matter density, $\chi_H = c/H_0$ the Hubble-distance and $\delta$ the density field. Computing the tidal shear fields $\partial^2\Phi$ shows that they must have the same fluctuation statistics as $\delta$: in Fourier-space, the solution to the Poisson equation is $\Phi \propto \delta/k^2$ with the wave vector $k$, and the tidal shear fields become $k_i k_j \Phi \propto k_i k_j / k^2 \delta$. Therefore, the power spectrum $P_{\partial^2\Phi}(k)$ of the tidal shear fields is proportional to the power spectrum $P_\delta(k)$ of the density fluctuations.

Consequently, one can derive the variance of the tidal shear fields from the variance of the density fluctuations, i.e. from the standard cold dark matter (CDM) power spectrum $P_\delta(k)$. Doing so, one can relate a mass scale M to the wave vector k by requiring that the mass M should be contained in a sphere of radius R,

$$ M = 4\pi\Delta\rho_{\text{crit}}R^3/3 = 32\pi^4\Delta\rho_{\text{crit}}/3k^3, \qquad (9) $$

with $k = 2\pi/R$, $\rho_{\text{crit}} = 3H_0^2/8\pi G$ the critical density and $\Delta = 200$ an overdensity factor.<fnref n="1" /> This defines a scale k in the power spectrum which is proportional to $M^{-1/3}$. This implies that on galaxy and cluster scales, where the CDM-spectrum scales $\propto k^\gamma$, with

$\gamma \in [-3, -2]$, one obtains for the standard deviation of the tidal shear field a behaviour $\propto M^{\beta_M}$, with $\beta_M \in [1/3, 1/2]$; in particular, for cluster-sized objects ($R = 0.5-1.5$ Mpc), where the non-linear matter power spectrum is proportional to $k^{-2.2}$ (Blas, Lesgourgues & Tram 2011), our prediction is $\beta_M \simeq 0.36$.

We can take this approach a step further: given the premises above, for the amplitude of the IA $A_{IA}$ we can write

$$ A_{IA} \propto \partial^2\Phi \propto \sqrt{P_{\partial^2\Phi}(k)} \propto \sqrt{P_\delta(k)} \propto \sqrt{P_\delta(M)}, \qquad (10) $$

where in the last step we assume that equation (9) holds and provides the link between the wave vector k and the mass M. This approach overcomes the previous approximation that the power spectrum follows a power law, and provides a single-parameter model with only a free amplitude; we describe how we test this model in Section 4.2.

We emphasize that we are primarily interested in the scaling between ellipticity and tidal shear, which is entirely due to the scale- (or mass-) dependence of the tidal shear field and not due to the internal dynamics of the halo. Predicting the dimensionless constant of proportionality would require many more assumptions in relation to the internal structure of the halo and the size of the system for restricting an otherwise diverging integration.

Comparing our work with Catelan, Kamionkowski & Blandford (2001), we would like to point out that our virial argument provides indications that the ellipticity of an aligned galaxy is proportional to the tidal shear field and that this proportionality does not depend on mass or redshift. Any scale-dependence of alignments is due to the scale- (or mass-) dependence of the tidal shear fields themselves.

We also emphasize that our model assumes a spherically symmetric, self-gravitating halo in virial equilibrium, with an isotropic velocity distribution with constant dispersion. A consistent solution for such a system would be the singular isothermal sphere with a density profile $\rho \propto 1/r^2$, for which we need to assume a cut-off radius in order to obtain finite values for the quadrupole moments, similarly to Camelio & Lombardi (2015). Our virial argument, on the other hand, does not differentiate between the dark matter and stellar components and, in particular, does not assume a potential of the unperturbed galaxy. Gravitational tidal fields would, due to equivalence, act on both components, which should be in virial equilibrium to each other, in the same way: it would not be possible to keep the potential generated by the dark matter component fixed and expose the stellar component to this potential. Furthermore, the potential of the system needs to be consistent with the particle distribution in configuration- and velocity-space, and cannot be chosen independently. These are two points in which we hope to improve the investigation by Camelio & Lombardi (2015).

# 3 DATA

## 3.1 Simulations

In this work, we consider haloes from two different simulations:

(i) the Millennium simulation (MS), first presented in Springel et al. (2005), which uses $2160^3$ dark matter particles of mass $m_p^{\text{MS}} = 8.6 \times 10^8 \text{ M}_\odot h^{-1}$ enclosed in a $500 \text{ Mpc} h^{-1}$-side box. In particular, we consider 2 of its 64 snapshots, i.e. the one at $z = 0$ (snapshot 63), for our baseline, and the one at $z \simeq 0.46$ (snapshot 49), to study a potential redshift dependence of our results. Dark matter haloes are identified as in Joachimi et al. (2013a) and references therein: a simple 'friends-of-friends' group-finder

<footnote n="1">¹ Here, the overdensity is relative to the critical density, but note we also consider overdensities relative to the underlying mean matter density. In those cases, we just rescale the mass by a factor of $\Omega_m$, which is the ratio between the mean matter density and the critical density.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Histogram showing the number density distribution of the mass of the haloes in the two simulations. Each logarithmic mass bin has a size of 0.1 dex. In this case, the mass of the halo is defined as $M_{200c}$, the mass within a region where the density exceeds 200 times the critical density; more details about the chosen mass ranges are presented in Section 4.2. In the range where the selected bins overlap, the trend of the density agrees for the two simulations.</figcaption></figure>

(Davis et al. 1985) is run first to spot virialized structures, followed by the SUBFIND algorithm (Springel et al. 2001, 2005) to identify sub-haloes, some of which are then treated as separate haloes if they are only temporarily close to the halo. We only consider haloes with a minimum number of particles $N_P = 300$ (Bett et al. 2007).

(ii) The Millennium-XXL simulation (MXXL), which samples 6720³ dark matter particles of mass $m_P^{MXXL} = 6.174 \times 10^9 M_{\odot} h^{-1}$ confined in a cubic region of 3000 Mpc $h^{-1}$ on a side (Angulo et al. 2012). In this case, we consider only one snapshot, at $z=0$. Haloes are selected using an ellipsoidal overdensity algorithm, as described in Despali, Tormen & Sheth (2013) and Bonamigo et al. (2015): a traditional spherical overdensity algorithm (Lacey & Cole 1994) gives an initial estimate of the true shape and orientation of the halo, which is then improved by building up an ellipsoid using the previously selected particles.

We define the mass of the objects in the catalogues as the mass within a region where the density is 200 times the critical density at the redshift corresponding to the respective snapshot ($M_{200c}$). Note that for the MS we first convert the halo mass from $M_{\text{Dhalo}}$, as defined in Jiang et al. (2014), to $M_{200c}$ using the median line in Jiang et al. (2014, fig. 2); this transformation is necessary in order to have a consistent definition of the mass of the haloes in the simulations, but its impact on our results is negligible. The number density distribution of the haloes used in our analysis is shown in Fig. 1.

In the simulations the same set of cosmological parameters is adopted, namely they both assume a spatially flat $\Lambda$CDM universe with the total matter density $\Omega_m = \Omega_b + \Omega_{\text{dm}} = 0.25$, where $\Omega_b = 0.045$ indicates the baryon density parameter and $\Omega_{\text{dm}} = 0.205$ represents the dark matter density parameter, a cosmological constant $\Omega_{\Lambda} = 1 - \Omega_m = 0.75$, the dimensionless Hubble parameter $h = 0.73$, the scalar spectral index $n_s = 1$ and the density variance in spheres of radius $8 \text{ Mpc } h^{-1}$, $\sigma_8 = 0.9$.

## 3.2 Halo shapes

We define the simple inertia tensor,<fnref n="2" /> whose eigenvalues and eigenvectors describe the shape of the halo, as

$$M_{\mu\nu} \propto \sum_{i=1}^{N_P} r_{i,\mu} r_{i,\nu}, \qquad (11)$$

where $N_P$ is the total number of particles within the halo, $\mu, \nu \in \{1, 2, 3\}$, and $r_i$ is the vector that indicates the position of the $i$th particle with respect to the centre of the halo, i.e. the location of the gravitational potential minimum, in the reference frame of the simulation box. For the MS only, we also consider a reduced inertia tensor, which is defined as (Pereira, Bryan & Gill 2008)

$$M_{\mu\nu}^{\text{red}} \propto \sum_{i=1}^{N_P} \frac{r_{i,\mu} r_{i,\nu}}{r_i^2}, \qquad (12)$$

with $r_i^2$ the square of the 3D distance of the $i$th particle from the centre of the halo. The reduced inertia tensor is more weighted towards the centre of the halo, and may yield a more reliable approximation of the shape of the galaxy at its centre (Joachimi et al. 2013b; Chisari et al. 2015).

The eigenvectors and eigenvalues define an ellipsoid, which we project on to one of the faces of the simulation box along the z-axis: the resulting ellipse is the projected shape of the halo. We proceed as in Joachimi et al. (2013a) to define the ellipticity $\epsilon$ of the objects, adopting their procedure for early-type galaxies. We denote the eigenvalues of the inertia tensor as $\omega_{\mu}$, and the three eigenvectors as $s_{\mu} = \{s_{x,\mu}, s_{y,\mu}, s_{z,\mu}\}^{\top}$, $\mu \in \{1, 2, 3\}$, where the coordinates refer to the Cartesian system of the simulation box. The projected ellipse is given by the points $x$ which satisfy $x^{\top} W^{-1} x = 1$, defining a symmetric tensor

$$W^{-1} = \sum_{\mu=1}^{3} \frac{s_{\perp,\mu} s_{\perp,\mu}^{\top}}{\omega_{\mu}^{2}} - \frac{\kappa \kappa^{\top}}{\alpha^{2}}, \qquad (13)$$

with $s_{\perp,\mu} = \{s_{x,\mu}, s_{y,\mu}\}^{\top}$,

$$\kappa = \sum_{\mu=1}^{3} \frac{s_{\perp,\mu} s_{\perp,\mu}^{\top}}{\omega_{\mu}^{2}}, \qquad (14)$$

and

$$\alpha^2 = \sum_{\mu=1}^{3} \left( \frac{s_{\perp,\mu}}{\omega_{\mu}} \right)^2. \qquad (15)$$

We compute the two Cartesian components of the ellipticity (Bartelmann & Schneider 2001)

$$\epsilon_1 = \frac{W_{11} - W_{22}}{W_{11} + W_{22} + 2\sqrt{\det W}}, \qquad (16)$$

$$\epsilon_2 = \frac{2W_{12}}{W_{11} + W_{22} + 2\sqrt{\det W}}, \qquad (17)$$

<footnote n="2">² MS and MXXL use two different tensor definitions to describe the shape, but they result in the same halo ellipticity (see also Bett et al. 2007 for further details).</footnote>
<table><caption>**Table 1.** Mass range, number of haloes $N_h$, mean ellipticity $\langle\epsilon\rangle$ and ellipticity dispersion $\sigma_\epsilon$ for each of the four bins of both the MS and MXXL. The masses are defined as $M_{200c}$, the mass within a region where the density exceeds 200 times the critical density.</caption><thead><tr><th></th><th>Mass range (M<sub>☉</sub> h<sup>-1</sup>)</th><th>N<sub>h</sub></th><th>⟨ε⟩</th><th>σ<sub>ε</sub></th></tr></thead><tbody><tr><td rowspan="4">Millennium simulation</td><td>10<sup>11.36</sup>–10<sup>11.86</sup></td><td>916 109</td><td>0.13</td><td>0.15</td></tr><tr><td>10<sup>11.86</sup>–10<sup>12.36</sup></td><td>328 364</td><td>0.13</td><td>0.15</td></tr><tr><td>10<sup>12.36</sup>–10<sup>12.86</sup></td><td>113 917</td><td>0.14</td><td>0.16</td></tr><tr><td>10<sup>12.86</sup>–10<sup>13.36</sup></td><td>37 511</td><td>0.16</td><td>0.18</td></tr><tr><td rowspan="4">Millennium-XXL simulation</td><td>10<sup>13</sup>–10<sup>13.5</sup></td><td>5450 501</td><td>0.19</td><td>0.22</td></tr><tr><td>10<sup>13.5</sup>–10<sup>14</sup></td><td>1618 707</td><td>0.21</td><td>0.24</td></tr><tr><td>10<sup>14</sup>–10<sup>14.5</sup></td><td>370 911</td><td>0.23</td><td>0.26</td></tr><tr><td>10<sup>14.5</sup>–10<sup>15</sup></td><td>48 714</td><td>0.25</td><td>0.28</td></tr></tbody></table>

which we then translate in the radial (+) component, following the IA sign convention:<fnref n="3" />

$$ \epsilon_+ = \epsilon_1 \cos(2\varphi) + \epsilon_2 \sin(2\varphi), \quad (18) $$

where $\varphi$ is the polar angle of the line that connects a halo pair with respect to the x-axis, in the reference frame of the simulation box. We show how we use $\epsilon_+$ to measure the correlation between halo shapes in the next section.

As a sanity check, we calculate the intrinsic ellipticity dispersion, $\sigma_\epsilon = \sqrt{\sum_{i=1}^{N_h} (\epsilon_1^2 + \epsilon_2^2)/N_h$ for several mass bins, with $N_h$ the number of haloes (details about the bins are reported in Section 4.2 and in Table 1). The magnitude of the dispersion is in good agreement with that of early-type galaxies (Joachimi et al. 2013a) and clusters (0.13–0.19 for the cluster samples used in van Uitert & Joachimi 2017). The dispersion increases with mass, as found in other simulations (e.g. Despali, Giocoli & Tormen 2014; Schrabback et al. 2015), with a small excess for the MXXL samples due to the lower sampling of the halo shapes.

# 4 METHODOLOGY

## 4.1 Measurements

To measure the alignment of every dark matter halo with other haloes, we mimic the standard analysis that is usually performed with observations (e.g. van Uitert & Joachimi 2017), i.e. we use the halo catalogue as the tracer of the underlying density field. Instead of fitting physical models to the alignment signal and dividing out the galaxy/cluster bias dependence afterwards, though, we take a shortcut: we choose one large bin in the comoving transverse separation $R_p$, and proceed as follows to remove the bias factor dependence.

We first define an estimator as a function of $R_p$ and the line-of-sight distance $\Pi$:

$$ \hat{\xi}_{g+}(R_p, \Pi) = \frac{S_+ D}{DD}, \quad (19) $$

where $S_+D$ represents the raw correlation between halo shapes ($\epsilon_+$) and the density sample, and DD the number of halo shape–density pairs. Per default, the halo shape sample is also used as the density

sample. We then integrate along the line of sight to obtain the total projected IA signal:

$$ \hat{w}_{g+}(R_p) = \int_{-\Pi_{\text{max}}}^{\Pi_{\text{max}}} d\Pi \, \hat{\xi}_{g+}(R_p, \Pi). \quad (20) $$

Throughout this work, we adopt $\Pi_{\text{max}} = 60 \, \text{Mpc} \, h^{-1}$, a value large enough not to miss part of the signal, but small enough not to pick up too much noise. We describe the IA signal by simplifying the model in equation (5) of van Uitert & Joachimi (2017), namely we assume

$$ w_{g+}(R_p, M) = A_{\text{IA}}(M) b_h(M) w_{\delta+}^{\text{model}}(R_p), \quad (21) $$

with $A_{\text{IA}}(M)$ the amplitude of the IA signal, $b_h(M)$ the halo bias and $w_{\delta+}^{\text{model}}(R_p)$ a function in which we include the dependence on $R_p$. In the tidal alignment paradigm, $w_{\delta+}^{\text{model}}$ is independent of halo mass, since it is fully determined by the properties of the dark matter distribution, assuming that any mass dependence of the response of a halo shape to the tidal gravitational field is captured by $A_{\text{IA}}(M)$.

We evaluate the expression in equation (21) in the interval which covers $10 \, \text{Mpc} \, h^{-1} < R_p < 20 \, \text{Mpc} \, h^{-1}$, denoted by $R_p^*$, to remove the dependence on $R_p$; in other words, we define the halo-pair weighted average in $R_p^*$ as

$$ w_{g+}^*(M) = \langle w_{g+}(R_p, M) \rangle_{R_p^*} \quad (22) $$

We choose this particular interval for two reasons. First, $10 \, \text{Mpc} \, h^{-1}$ is well above the minimum scale usually adopted in observational papers, below which the bias becomes non-linear (Tasitsiomi et al. 2004; Mandelbaum et al. 2006). Secondly, we gain little from extending the range beyond $20 \, \text{Mpc} \, h^{-1}$, as the signal-to-noise ratio of the IA signal is small at larger scales. As a sanity check, we repeat our analysis in an extended interval, which covers $6 \, \text{Mpc} \, h^{-1} < R_p < 30 \, \text{Mpc} \, h^{-1}$, without significantly changing our results. In a very broad bin the effective radial weighting of halo pairs could vary significantly across the mass range considered; hence we prefer the narrower bin as our default.

In order to constrain $b_h(M)$ in equation (21), we use the LS (Landy & Szalay 1993) estimator to calculate the clustering signal:

$$ \hat{\xi}_{gg}(R_p, \Pi) = \frac{DD - 2DR + RR}{RR}, \quad (23) $$

where DD represents the number of halo pairs, DR the number of halo-random point pairs and RR the number of random point pairs. To measure DR and RR, we generate random catalogues that contain objects uniformly distributed between the minimum and maximum value of the x, y and z coordinates of each simulation sub-box.<fnref n="4" /> These catalogues normally are three times denser, but in some cases, when a sub-box encloses very few objects, we switch to random catalogues which are 10 times denser. We always use equation (23) since we re-normalize the estimators according to the sample size. We then integrate along the line of sight to obtain the total projected clustering signal:

$$ \hat{w}_{gg}(R_p) = \int_{-\Pi_{\text{max}}}^{\Pi_{\text{max}}} d\Pi \, \hat{\xi}_{gg}(R_p, \Pi). \quad (24) $$

We describe the clustering signal with a simple model:

$$ w_{gg}(R_p, M) = b_h^2(M) w_{\delta\delta}^{\text{model}}(R_p) - C_{\text{IC}}, \quad (25) $$

with $w_{\delta\delta}^{\text{model}}(R_p)$ a function in which we include the dependence on $R_p$ (van Uitert & Joachimi 2017, equation 9), and $C_{\text{IC}}$ the integral

<footnote n="3">³ Our definition of $\epsilon_+$, as commonly done in IA works, has an opposite sign with respect to, for example, Bartelmann & Schneider (2001).</footnote>

<footnote n="4">⁴ For further information about the sub-boxes, see Section 4.2.</footnote>
constraint, which accounts for the offset due to the restricted area of the simulation box. We assess that this correction is negligible for the MXXL by measuring that the clustering signals are unaltered if we normalize the random catalogues with respect to the whole simulation box rather than the sub-boxes. We estimate the integral constraint for the MS by selecting an identical mass bin in both simulations ($10^{13} M_{\odot} h^{-1} - 10^{13.5} M_{\odot} h^{-1}$), and calculating $w_{gg}$ as a function of $R_p$. Assuming that the MXXL signal has a negligible $C_{\text{IC}}$, we find the optimal $C_{\text{IC}}$ for the MS signal via weighted least-squares fitting of the MS $w_{gg}$ to the MXXL $w_{gg}$. We then deduce the individual integral constraints for each MS mass bin by assuming that they scale with $b_h^2(M)$. Note that, since we are dealing with dark matter clustering, $w_{\delta\delta}^{\text{model}}(R_p)$ does not depend on the mass of the halo. Again, we average the previous expression in $R_p^*$, obtaining

$$w_{gg}^*(M) = \langle w_{gg}(R_p, M) \rangle_{R_p}^* \quad (26)$$

In observational analyses, the halo bias factor is usually measured from a fit to $w_{gg}$ and then included in the modelling of $w_{g+}$; equivalently, to remove the mass dependence of the halo bias $b_h(M)$, we define

$$r_{g+}(M) = \frac{w_{g+}^*(M)}{\sqrt{w_{gg}^*(M) + C_{\text{IC}}}} = \frac{A_{\text{IA}}(M) \langle w_{\delta+}^{\text{model}}(R_p) \rangle_{R_p}^*}{\sqrt{\langle w_{\delta\delta}^{\text{model}}(R_p) \rangle_{R_p}^*}} \propto A_{\text{IA}}(M), \quad (27)$$

where we assume that the clustering signal $w_{gg}^*(M)$ is positive (see Section 4.2 for further discussion). We stress that, under our assumptions, $r_{g+}$ depends only on the mass of the halo $M$.

## 4.2 Modelling

The goal of this paper is to study the halo mass dependence of the IA amplitude $A_{\text{IA}}(M)$, which we constrain by fitting $r_{g+}(M)$.

We select the haloes from the catalogues described in Section 3.1 in $n_M = 4$ logarithmic mass bins, each extending over two orders of magnitude, between $10^{1.36}$ and $10^{1.36} M_{\odot} h^{-1}$ for the MS and between $10^{13}$ and $10^{15} M_{\odot} h^{-1}$ for the MXXL (see Table 1). We split them in $N_{\text{sub}} = 3^3 = 27$ sub-boxes based on their positions inside the cube of the respective simulation, and calculate $w_{g+}^*$ and $w_{gg}^*$ for each sub-sample by replacing the integrals in equations (20) and (24) with a sum over 20 line-of-sight bins, each $2\Pi_{\text{max}}/20 = 6$ Mpc $h^{-1}$ wide, defining the line of sight along the z-axis.

We estimate the covariance matrix on the mean of the sub-boxes:

$$C_{\mu\nu} = \frac{1}{N_{\text{sub}}(N_{\text{sub}} - 1)} \sum_{j=1}^{N_{\text{sub}}} (d_{j,\mu} - \bar{d}_{\mu})(d_{j,\nu} - \bar{d}_{\nu}), \quad (28)$$

with $\mu, \nu \in \{1, \dots, n_M\}$, $\bar{d}_{\mu} = \frac{1}{N_{\text{sub}}} \sum_{j=1}^{N_{\text{sub}}} d_{j,\mu}$ and $d_{j,\mu} = r_{g+}(M)$ for each sub-box and each mass bin, as defined in equation (27). We then invert the covariance matrix and correct the bias on the inverse to obtain an unbiased estimate of the precision matrix, given by

$$C_{\text{unbiased}}^{-1} = \frac{N_{\text{sub}} - n_M - 2}{N_{\text{sub}} - 1} C^{-1}, \quad (29)$$

where $N_{\text{sub}} > n_M + 2$ clearly holds (Kaufman 1967; Hartlap, Simon & Schneider 2007; Taylor, Joachimi & Kitching 2013). We assess that this procedure does not introduce any systematic errors in the covariance assigned to our measurements by calculating our signals in $N_{\text{sub}}$ MS-sized boxes selected from the MXXL, and then observing that the results do not significantly differ from the MS signals.

The choice of $n_M$ and $N_{\text{sub}}$ is constrained by several factors: first of all, if $N_{\text{sub}}$ is too large, the single values of $w_{gg}^*$ (and $w_{g+}^*$)

tend to fluctuate around the mean, thus increasing the error bar and sometimes dropping below 0, which is unacceptable for our choice of $r_{g+}(M)$; see equation (27). Furthermore, we want $n_M$ to be large enough to be capable of displaying the trend of the signals along the whole mass range chosen. Finally, we need to take $n_M \ll N_{\text{sub}}$ to avoid divergences related to the fact that we estimate the covariance from a finite number of samples (Taylor et al. 2013).

As a first approach, which we will denote as ‘power-law’ approach, we adopt the following model for $r_{g+}(M)$:

$$r_{g+}(M) = A \left( \frac{M}{M_p} \right)^{\beta_M}, \quad (30)$$

with A a generic amplitude which we will treat as a nuisance parameter, $M_p = 10^{13.5} M_{\odot}/h_{70}$ a pivot mass and $\beta_M$ a free power-law index, which we intend to compare with the value predicted in Section 2. We perform a likelihood analysis over the data to infer the posteriors of A and $\beta_M$. According to Bayes’ theorem, if d is the vector of the data and p the vector of the parameters,

$$P(p|d) \propto P(d|p) P(p) \propto e^{-\frac{1}{2} \chi^2} P(p), \quad (31)$$

with $P(p|d)$ the posterior probability, $P(d|p)$ the likelihood function, $P(p)$ the prior probability and $\chi^2 = (d - m)^T C^{-1} (d - m)$, with m the vector of the model and C$^{-1}$ the precision matrix, the inverse of the covariance matrix C. We assume uninformative flat priors in the fit with ranges $\log_{10} \frac{A}{(M_{\text{pc}}h^{-1})^{1/2}} \in [-1.3; -0.6]$ and $\beta_M \in [0.2; 0.7]$.

As a second approach, which we will denote as ‘power spectrum’ approach, we use the CLASS (Cosmic Linear Anisotropy Solving System, Blas et al. 2011) algorithm to generate a non-linear matter power spectrum using the simulation cosmology. We then use equation (9) to convert the wavenumber to the corresponding mass M, and calculate the square root of the values of the power spectrum $P_\delta(M)$, as in equation (10). We adopt the following model for $r_{g+}(M)$:

$$r_{g+}(M) = A_{\text{PS}} \sqrt{P_\delta(M)}, \quad (32)$$

with $A_{\text{PS}}$ a generic amplitude which we will treat as a nuisance parameter.

To be able to directly compare our results with those in van Uitert & Joachimi (2017, fig. 7), we further convert the halo mass definition from $M_{200c}$ to $M_{200m}$, defined as the mass enclosed in a region inside of which the density is 200 times the mean density at the redshift corresponding to the respective snapshot.

# 5 RESULTS AND DISCUSSION

## 5.1 Simulation data

The trend of $w_{g+}$ and $w_{gg}$ with $R_p$ at $z=0$ for the four mass bins for each simulation is shown in Fig. 2. The points shown in Fig. 2 are the arithmetic mean of the $N_{\text{sub}}$ values for each mass bin, while the error bars are the standard deviation of the mean values. The overall behaviour of $w_{g+}$ and $w_{gg}$ agrees with previous works (Joachimi et al. 2011; van Uitert & Joachimi 2017); in particular, we clearly detect positive IA in all samples, implying that the projected ellipticities of dark matter haloes tend to point towards the position of other haloes. Also, it is worth noting that both signals increase with increasing mass.

We then study the dependence of $w_{g+}^*, w_{gg}^*$ and $r_{g+}$ on the mass of the halo. As one can see from Fig. 3, despite the use of two different halo finders, the MS and the MXXL follow the same
<figure><img src="image_4.png" /><figcaption>(a) $w_{g+}$ and $w_{gg}$ for the Millennium simulation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) $w_{g+}$ and $w_{gg}$ for the Millennium-XXL simulation.</figcaption></figure>

**Figure 2.** The IA signal $w_{g+}$ and the clustering signal $w_{gg}$ as a function of the comoving transverse separation $R_p$ at $z=0$ for (a) the MS and (b) the MXXL. The pink lines indicate the $10 < R_p/\text{Mpc} h^{-1} < 20$ interval, the range used to model the signal. The mass ranges are displayed considering $M_{200m}$ as the mass of the halo. In the graph, points are slightly horizontally shifted, so that they do not overlap; negative values are displayed in absolute value with open symbols of the same colour. An increasing trend with mass is clear in each panel separately, and comparing the two panels as well.

trend, and yield consistent results in the small mass range where they overlap; furthermore, all three $w_{g+}^*$, $w_{gg}^*$ and $r_{g+}$ increase with increasing mass. In Fig. 3 we also include, for the MS only, two more results: grey diamonds represent the signal from the objects at redshift $z=0.46$, while open black diamonds represent the signal from the objects at $z=0$ obtained using the reduced inertia tensor (rit, as in equation 12) to measure the shapes of the haloes. We find that the use of the reduced inertia tensor leads to lower alignment signals, and that these signals increase with increasing redshift, as found in Joachimi et al. (2013b). Importantly, we do not see any substantial deviations from the general trend with mass for these alternative measurements, suggesting that the details about redshift and halo shape definition are absorbed in the nuisance parameter A.

We proceed by showing the results of the likelihood analysis of the ‘power-law’ approach described in Section 4.2: Fig. 4(a) shows the results from the individual data sets and from the joint analysis of the two simulations, obtained by multiplying the likelihood functions and assuming the same flat priors on the parameters. In

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

**Figure 3.** The IA signal $w_{g+}^*$, the clustering signal $w_{gg}^*$ and $r_{g+}$ as defined in equation (27) as a function of halo mass $M_{200m}$ for the MS and the MXXL, calculated as a weighted average in the interval which covers $10 \text{ Mpc} h^{-1} < R_p < 20 \text{ Mpc} h^{-1}$. The label *sit* stands for simple inertia tensor, while *rit* means reduced inertia tensor; note that for the MXXL data only the simple inertia tensor is available, as mentioned in Section 3.2. The points are not placed at the midpoint of the bin, but at the value corresponding to the arithmetic mean of the mass of the objects. The red line and the green line represent the best-fitting lines for the MS and MXXL likelihood analyses, respectively, which are drawn using the parameters reported in Table 2, while the purple dashed line represents the best-fitting line for the ‘power spectrum’ model. With both approaches, we exclude the highest mass bin from the fit. Points showing the results from the *rit* choice are horizontally shifted by a small amount, so that they do not overlap with the corresponding *sit* dots.

the MXXL and in the joint fit, we exclude the highest mass point, since it clearly does not follow our model predictions. We test and discuss possible reasons for this deviation below.

The most stringent bounds come from the MXXL, while the MS yields larger errors on the parameters, albeit consistent with the results of the MXXL within $3\sigma$. The joint analysis returns a value for the slope compatible with $\beta_M = 1/3$, which agrees with the predictions made in Section 2 for a DM-only universe; in particular, the tightest constraints are obtained with cluster-sized objects, for which we predicted $\beta_M \simeq 0.36$, in agreement with our MXXL high-mass results ($\beta_M = 0.38^{+0.01}_{-0.01}$). Moreover, we find that neither the inertia tensor definition nor the chosen redshift for our default analysis have significant impact on our conclusions for $\beta_M$. We report all the best-fitting values, together with their respective errors and reduced $\chi^2$, in Table 2.

In Fig. 3, we also show the result of our ‘power spectrum’ approach: we choose $\Delta = 200 \Omega_m$ in equation (9) to be consistent with the definition of the mass, and in the fit we exclude the highest mass result. Our best-fitting line corresponds to $\log_{10}(A_{\text{PS}}/\text{Mpc} h^{-1}) = -1.54$ and $\chi^2/\text{d.o.f.} = 6.16$. In this case, the large value of the reduced $\chi^2$ is largely driven by the two lowest mass points.

Although our simple model impressively reproduces the general mass trend, the bad goodness of the fit in most cases indicates that there are deviations which are significant beyond the very small statistical error bars of our simulation analysis. To test whether
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Posterior on the fit parameters of the IA model, obtained for (a) the MS, the MXXL, joint MS and MXXL and (b) real data. Note that the ranges of the prefactors are quite different. The bottom-left graph in the top set of panels shows the contour lines and the 2D grey-scale posterior for the real data, while the bottom-left graph in the lower set of panels shows the contour lines of the 2D posteriors for all the simulations (single and joint), but the 2D grey-scale posterior for the joint analysis only. All other sub-panels show the marginalized 1D posterior normalized to a peak amplitude of 1. Contour lines enclose the 68 per cent and 95 per cent confidence intervals, crosses and vertical solid lines indicate the best-fitting values, while dashed lines represent the $1\sigma$ confidence interval. We note that the MS returns larger error bars, but the results are consistent for the two catalogues, while real data yield a value of the slope which is incompatible with the one from the joint analysis. The different degeneracy for the parameters in the two simulations is due to the choice of the pivot mass $M_p = 10^{13.5} M_{\odot}/h_{70}$: for the MS, it is at the upper end of the mass scale, while for the MXXL it is at the lower end. The mean values and 68 per cent confidence intervals of A, $A_r$ and $\beta_M$ are listed in Table 2.</figcaption></figure>

these deviations could originate from our treatment of the halo bias, we repeat our measurements taking the fourth and the first mass bin for the MS and the MXXL, respectively, as our density tracers for all samples in each simulation. We re-calculate $w_{g+}^*$ as the cross-correlation between the new density tracer and the shape sample in the respective mass bins, and derive $r_{g+}$ using $w_{gg}^*$ of the new density tracer. We present the results in Fig. 5. We do not observe any significant changes to the trends at the low- and high-mass ends and hence conclude that our removal of halo bias is robust. The deviation from our theoretical prediction at low masses appears even slightly more significant. It also suggests that the integral constraint

is accounted for correctly, and that the deviations are driven by the alignment rather than the clustering of haloes.

The decreased alignment strength in the $10^{14.5}-10^{15} M_{\odot}$ bin could be caused by a substantial fraction of haloes still in the process of formation at $z=0$, violating the assumption of virial equilibrium in our model. At the low-mass end, instead, alignments are stronger than predicted by the model, which suggests that a different, or additional, alignment mechanism is at play. While our MS halo sample does not contain satellites, defined here as gravitationally bound sub-haloes, many haloes of Milky Way mass and less could be infalling along filaments on to larger haloes, a
<table><caption>**Table 2.** Mean and 68 per cent confidence interval of the power-law fit parameters of the IA model from the likelihood analysis over the MS, the MXXL, their joint contribution, the MS at $z = 0.46$, the MS using the reduced inertia tensor (rit) and real data. We exclude the highest mass point in the MXXL only and joint fits. Note that the values from the snapshot at different redshift and from the reduced inertia tensor assumption are compatible with the outcomes of the MS only. A discussion about the reasons why the reduced $\chi^2$ values obtained considering the MXXL (and, consequently, the joint analysis) significantly differ from unity is presented in the text.</caption><thead><tr><th></th><th>MS only</th><th>MXXL only</th><th>Joint</th><th>MS, z = 0.46</th><th>MS, rit</th></tr></thead><tbody><tr><td>β<sub>M</sub></td><td>0.31<sup>+0.03</sup><sub>-0.03</sub></td><td>0.38<sup>+0.01</sup><sub>-0.01</sub></td><td>0.36<sup>+0.01</sup><sub>-0.01</sub></td><td>0.35<sup>+0.03</sup><sub>-0.03</sub></td><td>0.29<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>log<sub>10</sub>(A [Mpc h<sup>-1</sup>]<sup>1/2</sup>)</td><td>-0.96<sup>+0.04</sup><sub>-0.05</sub></td><td>-0.947<sup>+0.003</sup><sub>-0.003</sub></td><td>-0.942<sup>+0.003</sup><sub>-0.003</sub></td><td>-0.72<sup>+0.04</sup><sub>-0.04</sub></td><td>-1.21<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>χ<sup>2</sup>/d.o.f.</td><td>0.27</td><td>7.71</td><td>4.80</td><td>2.27</td><td>0.70</td></tr><tr><td colspan="6" style="text-align: center;">Real data</td></tr><tr><td>β<sub>M</sub></td><td>0.56<sup>+0.05</sup><sub>-0.05</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>log<sub>10</sub>A<sub>r</sub></td><td>0.61<sup>+0.03</sup><sub>-0.04</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>χ<sup>2</sup>/d.o.f.</td><td>1.68</td><td></td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Ratio of $r_{g+}$ as a function of halo mass for two sets of signals, with respect to the results of our main analysis shown in Fig. 3 as black and blue diamonds. Open diamonds are derived from the signals obtained with constant density tracers, which we measure to assess our treatment of the halo bias. Filled diamonds are derived from the signals obtained with all ellipticities rescaled to the mean in each sample (see equations 33 and 34), and demonstrate that the intrinsic halo ellipticity affects the mass scaling of the alignment signal. We neglect any correlations between the $r_{g+}$ measurements in estimating the error on the ratio.</figcaption></figure>

process that is distinct from our modelling ansatz (Forero-Romero, Contreras & Padilla 2014, and references therein). We assumed that the alignment is a perturbative effect to a largely randomly oriented intrinsic ellipticity $\epsilon$ (see equation 7), which is unlikely to hold in an infall/filament alignment scenario. Indeed, we see evidence that the intrinsic halo ellipticity affects the mass scaling of the alignment signal, violating a tenet of our model. We demonstrate this by normalizing the components of the ellipticity of each halo to

$$ \epsilon'_{1} = \frac{\epsilon_{1}}{\sqrt{\epsilon_{1}^{2} + \epsilon_{2}^{2}}} = \cos 2\phi, \qquad (33) $$

$$ \epsilon'_{2} = \frac{\epsilon_{2}}{\sqrt{\epsilon_{1}^{2} + \epsilon_{2}^{2}}} = \sin 2\phi, \qquad (34) $$

where $\phi$ is the polar angle of the major axis of the ellipse in the Cartesian coordinate system of the simulation box. Near-spherical haloes with $|\epsilon| < 0.01$ are excluded from further analysis. We then re-measure $w_{g+}^*$ and $r_{g+}$, which are multiplied by the mean ellipticity of the haloes in each bin, to make the signals directly compatible with the original results. As can be seen in Fig. 5, the

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Amplitude of the IA model as a function of halo mass from observations. The results from Joachimi et al. (2011) and Singh et al. (2015) are for LRGs, while van Uitert & Joachimi (2017) measured the signal for galaxy clusters. The solid purple line shows the best-fitting line from our likelihood analysis. The exact values of the parameters are shown in Table 2.</figcaption></figure>

mean-ellipticity $r_{g+}$ is generally about 15 per cent below the original signals, which implies that haloes with above-average ellipticity drive the alignment. Significant deviations arise for the second and third mass bin of the MS, confirming a breakdown of our model on the scales of galaxies.

## 5.2 Observation data

We repeat the likelihood analysis for the collection of observational data sets shown in van Uitert & Joachimi (2017, fig. 7): we consider all 21 data points, which include literature results for LOWZ LRGs from Singh et al. (2015), MegaZ-LRGs and SDSS LRGs from Joachimi et al. (2011), as well as results for the clusters contained in the redMaPPer catalogue version 6.3 from van Uitert & Joachimi.

We neglect the error bars on the mass, which are smaller than the errors on $A_{1A}$ and whose impact is subdominant (van Uitert & Joachimi), and treat all the data as independent, as in van Uitert & Joachimi. In this way, the covariance matrix is diagonal. We show the points, together with the best-fitting power law from our analysis, in Fig. 6.

We adopt the model of equation (30) but with a different prefactor $A_r$, which has an altered meaning and is now dimensionless,
thus making it impossible to directly compare its value to the one obtained with simulation data. We also assume a different range for the flat prior in the fit for this new parameter, namely $log_{10}A_r \in [0.4; 0.9]$. The outcomes of our analysis are shown in Fig. 4(b) and in Table 2: we observe that the value of the reduced chi-square for this latter analysis can be improved to 1.36 by excluding the high-redshift SDSS results (filled-blue diamonds) without affecting the value of the slope in a significant way.

The incompatibility between the values of the slope $\beta_M$ between simulation and real data is significant. This discrepancy is likely to be attributed to an additional mass dependence of the response of a galaxy ellipticity to the ellipticity of its host halo or the local tidal field, and, more generally, to the fact that, while we observe luminous matter, the simulations and the theory model only consider dark matter. Moreover, the effective scale at which ellipticities are measured in observational data changes from the central part of galaxies (since weak lensing techniques are employed) at the low-mass end, to the distribution of satellites in clusters that trace the halo out to the virial radius. We also note that a bias in $\beta_M$ could be determined from data if, with increasing halo mass, there is also a trend in redshift, which however is not the case of the data collection of Fig. 6.

# 6 CONCLUSIONS

In this work, we studied the dependence of the IA amplitude on the mass of dark matter haloes, using data from the MS and MXXL N-body simulations. We derived the IA amplitude scaling with mass in the tidal alignment paradigm for a dark matter-only universe. Our analytical estimate assumed a virialized system which is weakly perturbed by ambient tidal shear fields. In this model, the magnitude of the tidal distortion of a halo is determined entirely by the amount of fluctuations of the tidal fields on the mass scale of the halo. The model predicts a scaling with halo mass $\propto M^{\beta_M}$, with $\beta_M \in [1/3, 1/2]$, if we approximate the matter power spectrum with a power law.

We mimicked the observational approach to measure the halo shape-position alignments, and we performed a Bayesian analysis on the mass dependence of the alignments with a simple power-law model. We found that the results from the two simulation data sets agree very well in the mass range covered by both. Our model predicts the overall trend of the mass scaling of halo alignments remarkably well, but we observed significant deviations from the predictions at masses below $\sim 10^{13} M_{\odot}$ as well as in the highest mass bin near $10^{15} M_{\odot}$. We demonstrated that those are not caused by our analysis choices, nor by changes of the slope of the matter power spectrum as a function of scale. Rather, we found evidence that the physical processes driving the halo alignments in these mass ranges violate the assumptions of our simple model.

The joint analysis of the MS and MXXL yields $\beta_M = 0.36^{+0.01}_{-0.01}$, and for cluster-sized scales, for which we predicted $\beta_M = 0.36$, we found $\beta_M = 0.38^{+0.01}_{-0.01}$, however generally with bad goodness of fit due to the aforementioned deviations. Furthermore, there is no obvious dependence of the mass scaling on redshift or on the definition of the inertia tensor which describes the shape of the halo.

We repeated our statistical analysis using observational data, inferring a value of $\beta_M = 0.56^{+0.05}_{-0.05}$, which is not compatible with the simulation results or our model prediction. We argue that the incompatibility is attributed to the fact that simulations consider a dark

matter-only universe, while we observe luminous matter: hydrodynamical simulations suggest in fact that the presence of baryons could significantly affect the shape of the host halo, especially in the inner regions (Kiessling et al. 2015, and references therein).

For example, Bailin et al. (2005) found that baryons cause haloes to be closer to spherical but that their impact diminishes at larger radii (and thus with larger mass values), in agreement with Kazantzidis et al. (2004), and more recently Tenneti et al. (2014) and Velliscig et al. (2015a) discovered that galaxies are more misaligned with the hosting halo at lower masses. Also, Tenneti et al. (2015) compared the projected density-shape correlation function in three mass bins covering four orders of magnitude between a hydrodynamical and a dark matter-only simulation, showing that for the lowest mass bin the stellar matter signal is more than halved with respect to the dark matter signal, while at higher mass the fractional difference between the two does not exceed 40 per cent. Finally, the effective scale at which observations evaluate ellipticities moves outwards with increasing mass. All these trends contribute to a steepening of the mass scaling compared to dark matter halo alignments.

More quantitatively, Okumura & Jing (2009) showed the presence of a significant amount of misalignment between the orientations of LRGs and their dark matter haloes. They estimated a typical misalignment angle of $\sigma_{\theta} = 34.9^{+1.9}_{-2.1}$ deg, which causes the galaxy-shape correlation function to be about half the halo shape correlation function on all scales. We used this result to rescale the MXXL $r_g$+ values according to their mass, retrieving the median cosine of the misalignment angle for the highest mass bins from Velliscig et al. (2015a, fig. 8). With this approach, we obtain $\beta_M = 0.66^{+0.01}_{-0.01}$, a significantly higher value with respect to what we found using N-body simulations, and compatible within $2\sigma$ with the observational outcome.

These results hint that this discrepancy could be addressed by measuring the slope $\beta_M$ considering a hydrodynamical simulation large enough to contain clusters, which could then account for the additional effects of baryons and gas.

After completing this analysis, another study of the mass scaling of IAs in N-body simulations appeared. Xia et al. (2017) also develop an analytic model based on the tidal alignment paradigm, incorporating a mass dependence solely via the halo bias that enters by employing the halo power spectrum in lieu of the matter power spectrum. The authors find good agreement between predicted and measured halo ellipticity-ellipticity correlation functions. It will be interesting to compare the two ansätze in future work.

# ACKNOWLEDGEMENTS

We thank Raul Angulo for helpful information about the MXXL catalogue. DP acknowledges support by an Erasmus+ traineeship grant. BJ acknowledges support by an STFC Ernest Rutherford Fellowship, grant reference ST/J004421/1. MB acknowledges 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 program managed by the ANR. SH acknowledges support by the DFG cluster of excellence ‘Origin and Structure of the Universe’ (www.universe-cluster.de). EvU acknowledges support from an STFC Ernest Rutherford Research Grant, grant reference ST/L00285X/1. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program ‘Investissements d’Avenir’ supervised by the Agence Nationale pour la Recherche (ANR).
REFERENCES

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

Bailin J. et al., 2005, ApJ, 627, L17

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

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

Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, NJ

Blas D., Lesgourgues J., Tram T., 2011, J. Cosmol. Astropart. Phys., 7, 034

Blazek J., McQuinn M., Seljak U., 2011, J. Cosmol. Astropart. Phys., 5, 010

Blazek J., Vlah Z., Seljak U., 2015, J. Cosmol. Astropart. Phys., 8, 015

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

Camelio G., Lombardi M., 2015, A&A, 575, A113

Catelan P., Kamionkowski M., Blandford R. D., 2001, MNRAS, 320, L7

Chisari N. et al., 2015, MNRAS, 454, 2736

Codis S. et al., 2015, MNRAS, 448, 3391

Croft R. A. C., Metzler C. A., 2000, ApJ, 545, 561

Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371

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

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

Forero-Romero J. E., Contreras S., Padilla N., 2014, MNRAS, 443, 1090

Hao J., Kubo J. M., Feldmann R., Annis J., Johnston D. E., Lin H., McKay T. A., 2011, ApJ, 740, 39

Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399

Heavens A., Refregier A., Heymans C., 2000, MNRAS, 319, 649

Heymans C., White M., Heavens A., Vale C., Van Waerbeke L., 2006, MNRAS, 371, 750

Hilbert S., Xu D., Schneider P., Springel V., Vogelsberger M., Hernquist L., 2017, MNRAS, 468, 790

Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526

Hirata C. M., Mandelbaum R., Ishak M., Seljak U., Nichol R., Pimbblet K. A., Ross N. P., Wake D., 2007, MNRAS, 381, 1197

Jiang L., Helly J. C., Cole S., Frenk C. S., 2014, MNRAS, 440, 2115

Jing Y. P., 2002, MNRAS, 335, L89

Joachimi B., Mandelbaum R., Abdalla F. B., Bridle S. L., 2011, A&A, 527, A26

Joachimi B., Semboloni E., Bett P. E., Hartlap J., Hilbert S., Hoekstra H., Schneider P., Schrabback T., 2013a, MNRAS, 431, 477

Joachimi B., Semboloni E., Hilbert S., Bett P. E., Hartlap J., Hoekstra H., Schneider P., 2013b, MNRAS, 436, 819

Joachimi B. et al., 2015, Space Sci. Rev., 193, 1

Kaufman G. M., 1967, Report No. 6710. Center for Operations Research and Econometrics, Catholic University of Louvain, Heverlee, Belgium

Kazantzidis S., Kravtsov A. V., Zentner A. R., Allgood B., Nagai D., Moore B., 2004, ApJ, 611, L73

Kiessling A. et al., 2015, Space Sci. Rev., 193, 67

Kirk D., Rassat A., Host O., Bridle S., 2012, MNRAS, 424, 1647

Kirk D. et al., 2015, Space Sci. Rev., 193, 139

Kuhlen M., Diemand J., Madau P., 2007, ApJ, 671, 1135

L'Huillier B., Winther H. A., Mota D. F., Park C., Kim J., 2017, MNRAS, 468, 3174

Lacey C., Cole S., 1994, MNRAS, 271, 676

Landy S. D., Szalay A. S., 1993, ApJ, 412, 64

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lee J., Springel V., Pen U.-L., Lemson G., 2008, MNRAS, 389, 1266

Li C., Jing Y. P., Faltenbacher A., Wang J., 2013, ApJ, 770, L12

LSST Science Collaboration, 2009, preprint (arXiv:0912.0201)

Mandelbaum R., Hirata C. M., Ishak M., Seljak U., Brinkmann J., 2006, MNRAS, 367, 611

Mandelbaum R. et al., 2011, MNRAS, 410, 844

Okumura T., Jing Y. P., 2009, ApJ, 694, L83

Pereira M. J., Bryan G. L., Gill S. P. D., 2008, ApJ, 672, 825

Schmidt F., Chisari N. E., Dvorkin C., 2015, J. Cosmol. Astropart. Phys., 10, 032

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

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

Semboloni E., Heymans C., van Waerbeke L., Schneider P., 2008, MNRAS, 388, 991

Singh S., Mandelbaum R., More S., 2015, MNRAS, 450, 2195

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

Springel V. et al., 2005, Nature, 435, 629

Tasitsiomi A., Kravtsov A. V., Wechsler R. H., Primack J. R., 2004, ApJ, 614, 533

Taylor A., Joachimi B., Kitching T., 2013, MNRAS, 432, 1928

Tenneti A., Mandelbaum R., Di Matteo T., Feng Y., Khandai N., 2014, MNRAS, 441, 470

Tenneti A., Mandelbaum R., Di Matteo T., Kiessling A., Khandai N., 2015, MNRAS, 453, 469

Tenneti A., Mandelbaum R., Di Matteo T., 2016, MNRAS, 462, 2668

Troxel M. A., Ishak M., 2015, Phys. Rep., 558, 1

van Uitert E., Joachimi B., 2017, MNRAS, 468, 4502

Velliscig M. et al., 2015a, MNRAS, 453, 721

Velliscig M. et al., 2015b, MNRAS, 454, 3328

Xia Q., Kang X., Wang P., Luo Y., Yang X., Jing Y., Wang H., Mo H., 2017, ApJ, 848, 22

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