# Cosmic shear beyond 2-point statistics: Accounting for galaxy intrinsic alignment with projected tidal fields

Joachim Harnois-Déraps, Nicolas Martinet, Robert Reischke

► To cite this version:

Joachim Harnois-Déraps, Nicolas Martinet, Robert Reischke. Cosmic shear beyond 2-point statistics: Accounting for galaxy intrinsic alignment with projected tidal fields. Monthly Notices of the Royal Astronomical Society, 2022, 509, pp.3868. 10.1093/mnras/stab3222. hal-03313184

HAL Id: hal-03313184

https://hal.science/hal-03313184

Submitted on 14 Apr 2023

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>

# Cosmic shear beyond 2-point statistics: Accounting for galaxy intrinsic alignment with projected tidal fields

Joachim Harnois-Déraps<sup>1,*</sup> Nicolas Martinet<sup>2</sup> and Robert Reischke<sup>3</sup>

Accepted 2021 November 1. Received 2021 November 1; in original form 2021 July 21

## ABSTRACT

Developing analysis pipelines based on statistics beyond two-point functions is critical for extracting a maximal amount of cosmological information from current and upcoming weak lensing surveys. In this paper, we study the impact of the intrinsic alignment of galaxies (IA) on three promising probes measured from aperture mass maps – the lensing peaks, minima, and full PDF. Our 2D IA infusion method converts the light-cone-projected mass sheets into projected tidal tensors, which are then linearly coupled to an intrinsic ellipticity component with a strength controlled by the coupling parameter $A_{IA}$. We validate our method with the $\gamma$-2PCFs statistics, recovering well the linear alignment model of Bridle & King in a full tomographic setting, and for different $A_{IA}$ values. We next use our method to infuse at the galaxy catalogue level a non-linear IA model that includes the density-weighting term introduced in Blazek et al., and compute the impact on the three aperture mass map statistics. We find that large $S/N$ peaks are maximally affected, with deviations reaching 30 per cent (10 per cent) for a Euclid-like (KiDS-like) survey. Modelling the signal in a $w$CDM cosmology universe with $N$-body simulations, we forecast the cosmological bias caused by unmodelled IA for 100 deg$^2$ of Euclid-like data, finding very large offsets in $w_0$ (5-10$\sigma_{stat}$), $\Omega_m$ (4-6$\sigma_{stat}$), and $S_8 \equiv \sigma_8\sqrt{\Omega_m}/0.3$ (~3$\sigma_{stat}$). The method presented in this paper offers a compelling avenue to account for IA in beyond-two-point weak lensing statistics, with a flexibility comparable to that of current $\gamma$-2PCFs IA analytical models.

**Key words:** gravitational lensing: weak – methods: numerical – dark energy – dark matter – large-scale structure of Universe.

## 1 INTRODUCTION

Weak gravitational lensing by large-scale structure, commonly referred to as 'cosmic shear', has provided some of the most stringent constraints on the key cosmological parameters that describe the dark sector of our Universe (Hikage et al. 2019; Hamana et al. 2020; Amon et al. 2021; Asgari et al. 2021). Dedicated Stage-III survey such as the Kilo Degree Survey,<fnref n="1" /> the Dark Energy Survey<fnref n="2" /> and the HyperSuprime Camera Survey<fnref n="3" /> are either complete or nearing completion, and so far mostly agree on the value of $S_8 \equiv \sigma_8\sqrt{\Omega_m}/0.3$, the structure growth parameter that is best measured from lensing. The latter is defined as a combination of the matter density $\Omega_m$ and clumpiness $\sigma_8$, which specifies the amplitude of density fluctuations in spheres of $8h^{-1}$ Mpc.

These cosmological constraints are primarily inferred from measurements of two-point statistics, either the lensing power spectrum or the shear two-point correlation functions ($\gamma$-2PCFs), which are powerful summary statistics that can be analysed with prediction

models now reaching an accuracy of a few per cent (Euclid Collaboration et al. 2019; Nishimichi et al. 2019; Mead et al. 2021). The high accuracy attained by the two-point function analyses come however with a large cost in precision. Indeed, these methods completely overlook the non-Gaussian information contained in the mode coupling and in the phase correlation (Chiang & Coles 2000), which will become increasingly important for the upcoming generation (Stage-IV) of lensing surveys. This has led to the development of a number of alternative measurement techniques, among which the lensing peak count statistics has received a particularly large attention from the community, primarily for the simplicity of its methods and for its effectiveness at capturing additional information (Liu et al. 2015a, b; Kacprzak et al. 2016; Martinet et al. 2018; Shan et al. 2018; Harnois-Déraps et al. 2021). Other promising methods worth mentioning are the lensing PDF (Boyle et al. 2021; Martinet et al. 2021b, sometimes referred to the one-point statistics) and other moments of the convergence map (van Waerbeke et al. 2013; Gatti et al. 2020), lensing minima (Coulton et al. 2020), Minkowski functional (Petri et al. 2015), shear clipping (Giblin et al. 2018), lensing by voids (Davies et al. 2021), deep learning with Convolutional Neural Networks (Fluri et al. 2019), persistent homology analysis and scattering transform of the lensing field (respectively Heydenreich, Brück & Harnois-Déraps 2021; Cheng et al. 2020); these all present an appealing potential at improving the

© The Author(s) 2021.

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, NE1 7RU Newcastle-upon-Tyne, UK</footnote>

<footnote n="2"><sup>2</sup>Aix-Marseille Univ, CNRS, CNES, LAM, Marseille 13013, France</footnote>

<footnote n="3"><sup>3</sup>Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, D-44780 Bochum, Germany</footnote>

<footnote>* E-mail: joachim.harnois-deraps@ncl.ac.uk</footnote>

<footnote n="1">¹kids.strw.leidenuniv.nl</footnote>

<footnote n="2">²www.darkenergysurvey.org</footnote>

<footnote n="3">³hsc.mtk.nao.ac.jp/ssp/</footnote>
parameter constraints. For example, the joint peaks/γ-2PCFs have been shown to increase by a factor two the constraints on S<sub>8</sub>, by a factor three the precision on the dark energy equation-of-state parameter w<sub>0</sub> (Martinet et al. 2021b) and by ~40 per cent the precision on the sum of the neutrino masses (Li et al. 2019; Liu & Madhavacheril 2019) for Stage-IV surveys, compared to the γ-2PCFs alone.

An important feature common to many of these approaches is the absence of analytical models with which to predict the observed signals. Consequently, the cosmology inference must be completely calibrated from mock surveys constructed from suites of numerical weak lensing simulations, such as those described in Dietrich & Hartlap (2010), or the more recent MassiveNuS<sup>4</sup> (Liu et al. 2018) and cosmo-SLICS<sup>5</sup> (Harnois-Déraps, Giblin & Joachimi 2019).

Another key challenge faced by these techniques relates to the handling of systematic uncertainties that are known to exist in the lensing data (see Mandelbaum 2018, for a review), and that must therefore also be accounted for in these alternative measurement methods. In particular, the uncertainties associated to photometric redshifts, shape calibration/shear inference, baryonic feedback mechanism, and intrinsic alignment (IA) of galaxies must be carefully dealt with for the inferred cosmology to be accurate. While the former two of these effects can be included in N-body simulations at the ray-tracing level, the impact of baryon feedback on lensing statistics beyond two-point statistics typically requires to be calibrated on hydrodynamical simulations (Osato, Shirasaki & Yoshida 2015; Castro et al. 2018; Martinet et al. 2021a), although baryonification methods can also be of assistance (Schneider et al. 2019; Weiss et al. 2019; Lu & Haiman 2021).

The impact of IA has been mostly studied in the context of weak lensing two-point statistics, with only a few attempts to propagate the effect on to alternative statistics.<fnref n="6" /> There are in fact a number of theoretical models that attempt to describe the physics of IA, including the Non-Linear tidal Alignment model (NLA hereafter, see Bridle & King 2007), the tidal torquing model (e.g. Catelan, Kamionkowski & Blandford 2001; Hirata & Seljak 2004), a combination of both (the Tidal Alignment and Tidal Torquing model, TATT hereafter, described in Blazek et al. 2019) or halo-based alignment model (Schneider & Bridle 2010; Fortuna et al. 2021). These all produce predictions for two-point statistics, which generally depend on galaxy type, redshift, cosmology, with a few free parameters calibrated on hydrodynamical simulations (see e.g. Zjupa, Schäfer & Hahn 2020; Samuroff, Mandelbaum & Blazek 2021) and observations. For example, Joachimi et al. (2013a) detect and constrain an IA signal in the COSMOS early-type galaxies, but find hardly any signal for late-type galaxies; the WiggleZ blue galaxies are shown in Mandelbaum et al. (2011) to be consistent with no IA; Singh, Mandelbaum & More (2015) find a significant IA signal in the BOSS LOWZ sample; a similar trend was found by Johnston et al. (2019), from the KiDS, SDSS, and GAMA surveys, with no signal detection for the blue galaxies, but a 9σ detection for the red. When interpreted within in the NLA model, they measure an amplitude parameter of $A_{IA} = 3.18^{+0.47}_{-0.46}$. As reported in Samuroff et al. (2021), the different hydrodynamical simulations do not all agree on the physical model that best describes the IA. For example, the tidal torquing model was shown by Zjupa et al. (2020) to be strongly disfavoured over the

NLA model in the *IllustrisTNG* simulation,<fnref n="7" /> while it was observed for high-redshift blue galaxies in the Horizon-AGN simulations<fnref n="8" /> (Codis et al. 2015) and in many other hydrodynamical experiments (see Chisari et al. 2015, and references therein).

There is thus a large uncertainty on the strength of the galaxy alignments, and as a primary contaminant to the cosmic shear signal, it can be modelled and marginalized over from lensing analyses, providing indirect measurements of effective IA parameters. For example, the DES-Y1 reports non-zero values for the NLA model parameters $A_{IA}$ (the amplitude) and $\eta_{IA}$ (the redshift evolution) of $1.3^{+0.5}_{-0.6}$ and $3.7^{+1.0}_{-2.3}$ (Troxel et al. 2018), respectively. Similarly, the KiDS-1000 analysis of Asgari et al. (2021) find values of $A_{IA}$ in the range $0.264^{+0.424}_{-0.337} - 0.973^{+0.292}_{-0.383}$, depending on the choices of two-point statistics, while the HSC cosmic shear analyses report $A_{IA}$ of $0.38 \pm 0.70$ and $0.91^{+0.29}_{-0.38}$, also depending on the choice of statistics (Hikage et al. 2019; Hamana et al. 2020). The recent analysis of DES-Y3 favours slightly negative values of $A_{IA}$, both when analysed with the NLA and the TATT model (Secco et al. 2021).

The effect of IA is significant, and can cancel over 10–100 per cent of the observed cosmic shear signal, depending on the angular scale and redshift samples. If left unmodelled, the DES-Y1 and DES-Y3 inferred cosmologies would be biased by $1\sigma$ (Troxel et al. 2018; Secco et al. 2021), an effect that would worsen to 5–6σ in Stage-IV surveys such as Vera C. Rubin observatory<fnref n="9" /> or the Euclid<fnref n="10" /> and Nancy Grace Roman<fnref n="11" /> space telescopes, as reported in Joachimi et al. (2015, see their fig. 14). In fact, Blazek et al. (2019) have shown that if the true physics of the IA is described by the TATT, but the cosmology inference was assuming the NLA model, an LSST-like survey would still be several $\sigma$ away from the truth in many wCDM parameters.

However, it is still under debate how well the TATT model can describe the alignment process in general (Zjupa et al. 2020), since no alignment signal has been observed so far for blue galaxies. Samuroff et al. (2019) for example find an IA amplitude of $A_{IA} = 2.38^{+0.32}_{-0.31}$ for early-types described by the NLA model, but an amplitude consistent with zero for late-type galaxies which should preferentially align with the torquing mechanism captured by the TATT model according to some hydrodynamical simulation results. With the next generation of surveys the situation will change since the majority of their galaxy sample are expected to be blue galaxies.

IA are expected to be significant on alternative lensing statistics as well and must therefore be harnessed if we are to exploit their increased statistical power and interpret their results correctly. This topic has received very little attention so far, and it is this gap that we intend to fill in this paper. Our approach relies on the infusion of a physically motivated IA signal in mock lensing catalogues based on dark matter-only simulations. This has been explored to limited extend in the literature, either within a Halo Occupation Distribution model based on halo shapes (as in Heymans et al. 2006; Joachimi et al. 2013b), or by reweighting mass shells in the ray-tracing simulations (Fluri et al. 2019; Zürcher et al. 2021), however the accuracy and flexibility of these methods are not satisfactory given the current data, and certainly will not meet the requirements set for the upcoming Stage-IV lensing surveys. In particular, halo-based methods require an accurate measurement of the halo inertia matrix, which in turns requires hundreds of particles for its measurement and therefore

<footnote n="4">⁴columbialensing.org/#massivenus</footnote>

<footnote n="5">⁵slics.roe.ac.uk</footnote>

<footnote n="6">⁶For example, Gatti et al. (2020) modelled the IA in lensing moment via their dependence on the two-point functions.</footnote>

<footnote n="7">⁷www.tng-project.org/about/</footnote>

<footnote n="8">⁸horizon-simulation.org</footnote>

<footnote n="9">⁹lss.org</footnote>

<footnote n="10">¹⁰sci.esa.int/web/euclid</footnote>

<footnote n="11">¹¹roman.gsfc.nasa.gov</footnote>
places a low-mass cut on haloes for which an IA signal can be assigned.

We opted here instead to infuse IA in a manner that is physically consistent with the NLA model, e.g. based on a linear coupling between the intrinsic galaxy shapes and the non-linear projected tidal fields. We further include the first TATT term, a density-weighting component that increases the impact of IA on small scales with respect to the NLA. Our method has a number of advantages over the previous studies. First, an intrinsic shape is given to every galaxy, allowing for easy catalogue-based analyses for any choice of lensing statistics, with the possibility to dis- sect the underlying true IA contribution. We can also correctly combine the ellipticities and compute the reduced shear, which is the true lensing observable. Secondly, the coupling strengths are two free parameters that can be directly associated with the NLA $A_{\text{IA}}$ and the TATT $b_{\text{TA}}$ parameters, such that we can vali- date our infusion methods with two-point theoretical models. This connection further allows for joint analyses between two-point functions and higher order statistics, in which the IA parameters could be coherently varied and marginalized over. Thirdly, our method only requires the projected mass sheets, which are generally stored by default for most weak lensing simulations, meaning that it can be computed straightforwardly at different cosmologies from the public simulations suites mentioned above in order to investigate possible degeneracies. Fourthly, being based on the projected tidal field, it has the flexibility to adapt to many IA models, including those involving higher moments of the tidal field and/or local coupling between the tidal and density fields such as TATT.

This paper is structured as follow. We first briefly review in Section 2 the theory of weak lensing and intrinsic alignments, and describe in Section 3 our weak lensing simulations and IA infusion models. We validate in Section 4 our IA-infused mocks against the theoretical NLA predictions at the level of correlation functions, in the context of Stage-III lensing surveys. As a first demonstration of our methodology, we measure from the same galaxy catalogues the impact on peak statistics in Section 5, then explore the dependencies on smoothing scales and coupling strength. We next investigate the cosmological biases caused by IA for different analysis designs in the context of Stage-IV lensing surveys. More precisely, we carry out likelihood analyses based on three aperture mass maps statistics – lensing peaks, lens- ing minima, and lensing PDF – and their combination with the shear correlation functions, measuring in each case the impact on the inferred wCDM parameters $S_8$, $\Omega_m$, and $w_0$. We finally discuss our results and present our conclusions afterwards, in Section 6.

## 2 THEORY AND MODELLING

As mentioned in the introduction, there exist robust predictions for two-point statistics, however analytical modelling is highly limited when it comes to alternative measurement methods. This section briefly reviews the theory behind the modelling of cosmic shear correlation function and of the IA signal.

### 2.1 Cosmic shear 2-point functions

The cosmological information contained in cosmic shear data on large linear scales is well captured by two-point statistics, either in the form of the lensing power spectrum $C_\ell$ or the shear two-point correlation functions $\xi_\pm(\vartheta)$ ($\gamma$-2PCFs hereafter). Both are related to

the matter power spectrum $P_\delta(k, z)$, for which accurate prediction tools such as HALOFIT (Smith et al. 2003; Takahashi et al. 2012), HMCODE (Mead et al. 2021), the DARKEMULATOR (Nishimichi et al. 2019), or the EUCLIDEMULATOR (Euclid Collaboration 2019) exist. This connection largely explains why the two-point statistics stand as a natural choice for compressing and interpreting the weak lensing data. In the Limber approximation, the lensing power spectrum between tomographic bins ‘i’ and ‘j’ is computed as:

$$C_{\ell}^{ij} = \int_{0}^{\chi_{\text{H}}} \frac{q^{i}(\chi) q^{j}(\chi)}{\chi^{2}} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, z(\chi) \right) d\chi, \quad (1)$$

where $\chi_H$ is the comoving distance to the horizon. The lensing kernels $q^i$ depend on the redshift distribution $n(z)$ as:

$$q^{i}(\chi) = \frac{3}{2} \Omega_{\text{m}} \left(\frac{H_0}{c}\right)^2 \frac{\chi}{a(\chi)} \int_{\chi}^{\chi_{\text{H}}} n^{i}(\chi') \frac{dz}{d\chi'} \frac{\chi' - \chi}{\chi'} d\chi', \quad (2)$$

where $c$ is the speed of light, $a$ is the scale factor, and $H_0$ the Hubble parameter, note that $\chi_H = c/H_0$. Predictions for the $\gamma$-2PCFs are computed from equation (1) as:

$$\xi_{\pm}^{ij}(\vartheta) = \frac{1}{2\pi} \int_{0}^{\infty} C_{\ell}^{ij} J_{0/4}(\ell\vartheta) \ell \, d\ell, \quad (3)$$

where $J_{0/4}(x)$ are Bessel functions of the first kind. We adopt in this paper the Takahashi et al. (2012) HALOFIT model when computing $P_\delta(k, z)$ and use the cosmological parameter estimation code COSMOSIS$^{12}$ (Zuntz et al. 2015) to produce our $\xi_\pm^i$ predictions.

### 2.2 Non-linear alignment model

The observed ellipticity of a galaxy $\epsilon_{\text{obs}}$ is a combination of its intrinsic shape $\epsilon_{\text{int}}$ and a cosmic shear signal $\gamma$, the former of which can be further divided in a random component $\epsilon^{\text{ran}}$ and an alignment term $\epsilon^{\text{IA}}$. According to the NLA model,$^{13}$ IA are caused by a linear coupling between galaxy shapes and the non-linear large-scale tidal field at the galaxy position.

In the context of two-point function analyses, these intrinsic shapes contribute to an intrinsic–intrinsic (II) term as well as to an intrinsic-shear coupling (GI) term (Hirata & Seljak 2004), both secondary signals to the true cosmic shear (GG) term. The II and GI terms can be computed from the matter power spectrum as:

$$P_{\text{II}}(k, z) = \left( \frac{A_{\text{IA}} \bar{C}_1 \bar{\rho}(z)}{\bar{D}(z)} \right)^2 a^4(z) P_\delta(k, z) \quad (4)$$

and

$$P_{\text{GI}}(k, z) = - \frac{A_{\text{IA}} \bar{C}_1 \bar{\rho}(z)}{\bar{D}(z)} a^2(z) P_\delta(k, z), \quad (5)$$

which can then be Limber-integrated (as in equations 1 and 3) to compute the secondary signals $C_\ell^{II}$, $C_\ell^{GI}$ and $\xi_\pm^{II}(\vartheta)$, $\xi_\pm^{GI}(\vartheta)$. Here $\bar{D}$ is the ‘rescaled linear growth factor’ defined as $\bar{D} = D(1+z)$, $\bar{\rho}(z)$ is the matter density at redshift $z$ and $\bar{C}_1$ is a constant calibrated in Brown et al. (2002) that takes on the default value of $5 \times 10^{-14} M_\odot^{-1} h^{-2} \text{Mpc}^3$. The amplitude parameter $A_{\text{IA}}$ describes the strength of the tidal coupling and is the main NLA parameter constrained by current cosmic shear surveys. When phrased in terms of intrinsic ellipticities and tidal field $s_{ij}$, the scale factors from the

<footnote>$^{12}$bitbucket.org/joezuntz/cosmosis/wiki/Home</footnote>

<footnote>$^{13}$The term ‘non-linear’ in the model name refers to the use of the non-linear matter power spectrum $P(k)$ in the calculations; the coupling with the tidal field is still linear.</footnote>
previous equations exactly cancel the redshift evolution of the matter density field, leaving:

$$ 
\begin{align} 
\epsilon_1^{\text{IA}} &= -\frac{A_{\text{IA}} \bar{C}_1 \bar{\rho}(z=0)}{D(z)} (s_{xx} - s_{yy}), \\ 
\epsilon_2^{\text{IA}} &= -\frac{2 A_{\text{IA}} \bar{C}_1 \bar{\rho}(z=0)}{D(z)} s_{xy}, 
\end{align} 
\tag{6} $$

where $s_{ij} = \partial_{ij}\phi$ are the Cartesian components of the tidal tensor of the gravitational potential $\phi$. Note that it is now the standard growth factor that appears in the denominator.

As discussed later, the key quantity of interest to cosmic shear analyses is not the tidal tensor itself but its trace-free version, since shape correlations with the density field are subdominant, even though peaks in high density environments are less elliptical on average. Consequently, the model itself has a restricted range of validity: on small scales, higher order couplings to the ellipticity field $\epsilon^{\text{IA}}$ become important, however these are neglected in the NLA model. Only the non-linear evolution of the tidal tensor itself is taken into account, since it does fit observations quite well.

We note that the NLA predicts additional higher order terms and non-zero B-modes (Hirata & Seljak 2004) that we neglect in this analysis. Also, as shown in Vlah, Chisari & Schmidt (2020), the NLA model can be interpreted as the lowest order description of the alignment process of galaxies in the light of an effective field theory description.

## 2.3 Extension: δ-NLA

As mentioned in Blazek, Vlah & Seljak (2015) and Blazek et al. (2019), the intrinsic alignment of galaxies can only be observed at the galaxy positions, which are not randomly distributed on the sky but instead trace the underlying matter density ‘$\delta$’ with some biasing scheme. Accounting for this can be done at the level of theoretical predictions with a density weighting term computed in the abovementioned references from perturbation theory, and at the mock-infusion stage, by imposing (or) not our mock galaxies to trace the underlying matter field. What matters for a good match is that the simulated galaxy catalogues are being analysed with a consistent model. In presence of this sampling effect, the observed ellipticities from equation (6) become

$$ \epsilon_{1/2}^{\text{IA},\delta} = \epsilon_{1/2}^{\text{IA}} \times (1 + b_{\text{TA}}\delta). \tag{7} $$

The term $b_{\text{TA}}$ encodes the coupling strength with the local density field, and is often set to unity but could be allowed to vary. Note that the standard NLA described in last section does not include this term, hence we refer to this improvement as the δ-NLA model. The enhancement is mostly seen on small scales (Blazek et al. 2019), and has also been observed in hydrosimulations (Hilbert et al. 2017). This is fully consistent with the TATT model in which the torque term ($C_2$ or $A_2$, see Blazek et al. 2019; Secco et al. 2021) is nulled, and further reduces to the NLA described in the last section by setting $b_{\text{TA}}$ to zero.

## 3 SIMULATIONS

In this section, we briefly review the back-bone simulated light cones, then describe how we construct weak lensing catalogues infused with IA terms consistent with the NLA and δ-NLA models. We also detail the four weak lensing estimators on which the impact of IA are measured (with results presented in Sections 4 and 5).

### 3.1 Weak lensing light cones

Our simulated cosmic shear data are based on the Scinet Light-Cone Simulations (SLICS hereafter, see Harnois-Déraps et al. 2018) and the cosmo-SLICS (Harnois-Déraps et al. 2019), which are two complementary N-body suites designed for the analysis of Stage-III weak lensing data. The SLICS are a series of 928 fully independent runs at a fixed flat ΛCDM cosmology ($\Omega_m = 0.2905$, $\sigma_8 = 0.826$, $\Omega_b = 0.0474$, $h = 0.6898$, and $n_s = 0.969$), each evolving 1536³ matter particles in a box of side $L_{\text{box}} = 505 h^{-1}$ Mpc. These simulations are specifically tailored for estimating weak lensing covariance matrices, and were central to a number of cosmic shear (e.g. Hildebrandt et al. 2017; Giblin et al. 2018; Martinet et al. 2018; Harnois-Déraps et al. 2020) and combined-probe (Harnois-Déraps et al. 2016, 2017; Brouwer et al. 2018; Joudaki et al. 2018; van Uitert et al. 2018) data analysis. The cosmo-SLICS, in contrast, cover a wide range of wCDM parameter values, more specifically $\Omega_m$, $S_8$, $h$, and $w_0$. They sample 25 points organized in a Latin hypercube in order to minimize the interpolation error, each evolving a pair of N-body runs designed to suppress the sampling variance. One of these nodes, the fiducial model, lies at a ΛCDM cosmology identical to that of the SLICS, albeit with a $\sigma_8$ value that is 1.2 per cent higher.

Every N-body simulations produced sequences of projected density fields chosen such as to fill the past light cones up to $z = 3$ in steps of $L_{\text{box}}/2$, with randomized origins and projection axes. These $10 \times 10 \text{ deg}^2$ mass maps, which we label $\delta_{2D}(\theta)$, are used to produce multiple convergence and shear maps, $\kappa(\theta)$ and $\gamma_{1/2}(\theta)$, from which lensing quantities can be interpolated at a given galaxy redshift and position. We refer the interested reader to Harnois-Déraps et al. (2018) for complete details on how this is achieved.

### 3.2 Weak lensing catalogues

As mentioned in the introduction, we investigate in this paper the impact of IA for two different surveys configurations. We first construct mock data sets that resemble the Fourth Data Release of the Kilo Degree Survey, with shapes noise, galaxy density, and redshifts statistically matching the data properties described in Giblin et al. (2021) and Hildebrandt et al. (2021), respectively. These KiDS-like simulations are split in five tomographic bins, with redshift distributions shown in the upper panel of Fig. 1 and key survey properties listed in Table 1. Contrarily to the public SLICS and cosmo-SLICS mock data, in this work the positions of the mock galaxies are not placed at random, but assigned such as to trace the projected matter density with a linear galaxy bias set to unity. More precisely, the galaxies in each catalogue are ranked according to their redshift, they are next associated with the mass sheet $\delta_{2D}(\theta)$ that is the closest in redshift, then assigned a position by sampling $\delta_{2D}(\theta)$ (as described in Appendix A2 of Harnois-Déraps et al. 2018). This extra level of realism better reproduces the clustering properties seen in the data, but, more importantly for this study, it also up-weights the number of galaxies that live in dense environment and that therefore experience strong tidal fields. Overlooking this aspect would result in an underestimation of the importance of IA, and is the key difference between the NLA and the δ-NLA model (Blazek et al. 2019). These mock data have a relatively low signal-to-noise per tomographic bin, hence we produce 50 IA-infused KiDS-like catalogues from as many cosmo-SLICS light cones, at the fiducial cosmology. They primarily serve to establish the validity of our infusion method with $\gamma$-2PCFs in this paper, but are also ideally suited for including IA modelling in ongoing cosmic shear analyses beyond two-points statistics. As an example, we examine the effect of
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Redshift distribution of the five tomographic bins, measured from the KiDS-like (upper) and Euclid-like (lower) simulations. Note the change in the axes scaling between the two panels arising from different bin widths. In both panels, the $n(z)$ distributions for the different tomographic bins are normalized such that $\sum n(z)dz = 1.0$.</figcaption></figure>

<table><caption>**Table 1.** Properties of our KiDS-like and Euclid-like surveys. For the former survey, the effective number densities $n_{\text{eff}}$ (in gal arcmin$^{-2}$), shape noise per component $\sigma_{\epsilon}$ and redshift distributions listed here match those of the KiDS-1000 data presented in Giblin et al. (2021) and Hildebrandt et al. (2021). The column 'Z$_B$ range' refers to the photometric selection that defines the five KiDS-1000 tomographic bins, while the mean redshift in each bin is listed under $\langle z \rangle$. The specifications of the Euclid-like survey follows those presented in Martinet et al. (2021b), with $n_{\text{eff}} = 6.0$ gal arcmin$^{-2}$ for each tomographic bin and $\sigma_{\epsilon} = 0.26$ per component.</caption><thead><tr><th rowspan="2">tomo</th><th rowspan="2">Z<sub>B</sub> range</th><th colspan="3">KiDS-like</th><th colspan="2">Euclid-like</th></tr><tr><th>n<sub>eff</sub></th><th>σ<sub>ε</sub></th><th>⟨z⟩</th><th>z range</th><th>⟨z⟩</th></tr></thead><tbody><tr><td>bin1</td><td>0.1–0.3</td><td>0.616</td><td>0.270</td><td>0.257</td><td>0.0–0.4676</td><td>0.286</td></tr><tr><td>bin2</td><td>0.3–0.5</td><td>1.182</td><td>0.258</td><td>0.402</td><td>0.4676–0.7194</td><td>0.600</td></tr><tr><td>bin3</td><td>0.5–0.7</td><td>1.854</td><td>0.273</td><td>0.563</td><td>0.7194–0.9625</td><td>0.841</td></tr><tr><td>bin4</td><td>0.7–0.9</td><td>1.259</td><td>0.254</td><td>0.792</td><td>0.9625–1.3319</td><td>1.134</td></tr><tr><td>bin5</td><td>0.9–1.2</td><td>1.311</td><td>0.270</td><td>0.984</td><td>1.3319–3.0</td><td>1.852</td></tr></tbody></table>

IA on the peak count statistics measured from our mock catalogues, which will serve in an upcoming analysis of the actual KiDS DR4 data.

We next assemble a series of Euclid-like mock galaxy catalogues specifically designed for assessing the impact of IA on Stage-IV cosmic shear measurements, from the summary statistics all the way to the cosmological inference stage. For this purpose, we exploit the 928 SLICS and the 26 × 10 cosmo-SLICS light-cones introduced in Martinet et al. (2021b) to estimate the covariance matrix and to model the wCDM cosmological dependence of our estimators, respectively.<fnref n="14" /> We further construct a series of IA-infused mocks following the same method as for the KiDS-like catalogues, except

that the redshift distribution is now given by:

$$n(z) = A \frac{z^a + z^{ab}}{z^b + c}, \quad (8)$$

with $A = 1.7865$, $a = 0.4710$, $b = 5.1843$, $c = 0.7259$, and a galaxy density of $n_{\text{gal}} = 30$ gal arcmin$^{-2}$. The catalogues are next divided evenly across five tomographic bins, for a density of $n_{\text{gal}} = 6.0$ gal arcmin$^{-2}$ per bin. The shape noise is set to $\sigma_{\epsilon} = 0.26$ per component (see Table 1 for a summary of these properties). The resulting $n(z)$ is presented in the lower panel of Fig. 1. With a lower level of shape noise and a significantly higher galaxy density compared to the KiDS-like catalogues, Martinet et al. (2021a) find that 10 light-cones per cosmology is enough for their measurement (an analysis similar to ours but in the context of baryonic feedback). For all Euclid-like catalogues, five shape noise realizations are produced and averaged over in the end, to better capture the non-Gaussian signal.

### 3.3 Infusion of intrinsic alignments

The next step consists in computing the projected tidal fields within our simulations, which can then be coupled to an intrinsic galaxy shape according to our IA models. Given a 3D mass overdensity field $\delta(\mathbf{x})$ and a smoothing scale $\sigma_{\text{G}}$, the trace-free tidal tensor $s_{ij}(\mathbf{x})$ can be computed as (Catelan et al. 2001):

$$\tilde{s}_{ij}(\mathbf{k}) = \left[ \frac{k_i k_j}{k^2} - \frac{1}{3} \right] \tilde{\delta}(\mathbf{k}) \mathcal{G}(\sigma_{\text{G}}), \quad (9)$$

where tilde symbols $\tilde{\cdot}$, $\tilde{\cdot}^*$ denote Fourier transformed quantities, the indices $(i, j)$ label the components of the Cartesian wave vector $\mathbf{k}^T = (k_x, k_y, k_z)$, and $k^2 = k_x^2 + k_y^2 + k_z^2$. The smoothing kernel $\mathcal{G}(\sigma_{\text{G}})$ is typically a 3D Gaussian function described by a single parameter $\sigma_{\text{G}}$, which determines what physical scales are allowed to affect the IA term in our model. We show in Appendix A (see equation A3) that the components of the projected tidal shear $s_{ij,2D}(\theta) \equiv \sum_z s_{ij}(\mathbf{x})$ can be computed from the projected density $\delta_{2D}(\theta)$ as:

$$\tilde{s}_{ij,2D}(\mathbf{k}_{\perp}) = \sum_{z} \tilde{s}_{ij}(\mathbf{k}) = 2\pi \left[ \frac{k_i k_j}{k^2} - \frac{1}{3} \right] \tilde{\delta}_{2D}(\mathbf{k}_{\perp}) \mathcal{G}_{2D}(\sigma_{\text{G}}), \quad (10)$$

where $\mathbf{k}_{\perp}$ represents the two Fourier components perpendicular to the line of sight, the indices $(i, j)$ hereafter refer to either the x or the y component, and $\mathcal{G}_{2D}$ is now a 2D smoothing kernel. An important aspect of our model is that neither the functional form of $\mathcal{G}_{2D}(\sigma_{\text{G}})$ nor the value of $\sigma_{\text{G}}$ are specified within the NLA model. Smoothing amounts to selecting physical scales that do not contribute to the alignment of galaxies, which is still a debatable quantity, but is also introduced for numerical stability. Blazek et al. (2015) argues that $1.0 \ h^{-1} \ \text{Mpc}$ could be a reasonable fiducial value, being larger than the typical halo size, but recognizes that one-halo terms are also required to better match the observations. In their later work, Blazek et al. (2019) do not include smoothing at all, and neither do the KiDS-1000 nor DES-Y1 analyses based on the NLA model (Troxel et al. 2018; Asgari et al. 2021). In this work, we used a 2D Gaussian filter and calibrated $\sigma_{\text{G}}$ empirically to $\sigma_{\text{G}} = 0.1 \ h^{-1} \ \text{Mpc}$ (see Section 4.2), however these two choices are arbitrary and could possibly be better optimized in the future; we also explore $\sigma_{\text{G}} = 0.5 \ h^{-1} \ \text{Mpc}$ later on. We note that the smoothing scale is degenerate with the resolution of the simulation itself, and that one should use caution when smoothing on scales that approach the resolution limits.

We employ a numerical technique worth mentioning here: the Fourier transforms involved in computing equation (10) are computed from the full (projected) periodic boxes of the simulations, then interpolated on the light cones. We find that tidal field computed

<footnote n="14">¹⁴The galaxy positions were assigned at random in these wCDM and covariance simulations, and as opposed to the density-weighted positions used in the IA-infused mocks. Since we measure the relative impact of IA relative to the density-weighted catalogues with $A_{1\text{A}}$ set to 0.0, this systematic difference should not impact the inferred cosmology.</footnote>
directly from the 10 × 10 deg² light cones suffer from important large scales features, largely caused by the non-periodic boundary conditions, which our method avoids.

The projected tidal field maps $s_{11}$, $s_{22}$, and $s_{12}$ are constructed for each of the 18 mass sheets in every light cones and interpolated at the position of every galaxy. Fig. 2 illustrates this process for a small zoomed-in patch at $z \sim 0$, starting from a $\delta_{2D}(\theta)$ map (upper large panel), computing the tidal field components (middle four panels) and comparing the results with the cosmic shear signal generated by the same mass distribution (bottom two panels, see equation A5). The tidal field maps clearly reproduce the cosmic shear maps, and the minus sign in front of both terms in equation (6) causes the IA to undo some of the lensing signal.

Given a value of $A_{\text{IA}}$ and following equation (6), we couple the tidal field $s_{ij,2D}$ with the complex intrinsic ellipticities $\epsilon^{\text{IA}}$, from which we compute observed ellipticities as:

$$ \epsilon^{\text{obs}} = \frac{\epsilon^{\text{int}} + g}{1 + \epsilon^{\text{int,*}} g}, \quad \text{with} \quad \epsilon^{\text{int}} = \frac{\epsilon^{\text{IA}} + \epsilon^{\text{ran}}}{1 + \epsilon^{\text{IA,*}} \epsilon^{\text{ran}}}. \qquad (11) $$

In the above expressions, the denominators guarantee that no ellipticity component exceeds unity. The complex spin-2 reduced shear $g \equiv (\gamma_1 + i\gamma_2)/(1 + \kappa)$ is computed from the shear $(\gamma_{1/2})$ and convergence $(\kappa)$ maps, interpolated at the galaxy positions and redshifts. The complex spin-2 random orientation term $\epsilon^{\text{ran}}$ is drawn from two Gaussians (one per component) with their standard deviations provided by the shape noise level detailed in Table 1. We further constrain the random ellipticity to satisfy $|\epsilon^{\text{ran}}| \le 1.0$.

We recall that since our simulated galaxy catalogues trace the dark matter density, our default IA-infusion method is consistent with the $\delta$-NLA model (see Section 2.3). Analytical predictions for the two-point functions are more involved in this case (see Blazek et al. 2019) and not yet available on the public COSMOSIS release. We therefore validate our methods with the standard NLA first, which we infuse by ‘correcting’ the $\delta$-NLA measurement with equation (7), i.e. by replacing $\epsilon^{\text{IA}}$ by $\epsilon^{\text{IA}}/(1 + b_{\text{IA}}\delta)$ in our simulations. Once again, these are less realistic but better suited for validation of the measured $\gamma$-2PCFs against a theoretical model; after this is established, we use the $\delta$-NLA catalogues as our fiducial infusion model.

To be clear, the current version makes no differentiation between galaxy colours or type, and instead treats the full sample as a single population for which we measure an effective alignment signal, similar to the colour-free incarnation of the NLA model (see Samuroff et al. 2019, for an example with a red/blue split).

## 3.4 Data vector

We first set out to validate our IA-infusion method with a two-point analysis, in which we compare the $\gamma$-2PCFs extracted from the simulations with the analytical NLA predictions described in equation (3) and Section 2.2. The measurements are carried out with the public codes TREECORR (Jarvis, Bernstein & Jain 2004) and ATHENA (Kilbinger, Bonnett & Coupon 2014), which both compute, for each tomographic bin combination $(\alpha, \beta)$, the estimator $\tilde{\xi}_{\pm}^{\alpha\beta}(\vartheta)$ as:

$$ \tilde{\xi}_{\pm}^{\alpha\beta}(\vartheta) = \frac{\sum_{ab} W_a W_b \left[ \epsilon_{a,t}^{\alpha}(\theta_a) \epsilon_{b,t}^{\beta}(\theta_b) \pm \epsilon_{a,\times}^{\alpha}(\theta_a) \epsilon_{b,\times}^{\beta}(\theta_b) \right] \Delta\vartheta_{ab}}{\sum_{ab} W_a W_b}. \qquad (12) $$

The lensing weights $W_{a/b}$ are set to 1.0 throughout this paper, and $\epsilon_{a,t/\times}^{\alpha}(\theta_a)$ refers to the tangential/cross-component of the observed

<figure><img src="image_4.png" /><figcaption>Figure 2. (top:) Projected mass overdensity, measured at $z = 0.04$ at our fiducial cosmology, subtending 2.6 deg on the side. (second and third rows:) Different components of the tidal field tensor computed from the overdensity map shown on the top (see equation 10). (lower row:) The two lensing shear components produced by the same mass distribution, assuming a galaxy source plane at $z_s = 0.08$, i.e. just behind.</figcaption></figure>
ellipticity from galaxy 'a' with respect to the centre of the line that connects galaxies a and b. The binning operator $\Delta\theta_{ab}$ is set to unity if the angular distance between the galaxy pair falls inside the $\theta$ bin, and to zero otherwise; the sums run over all galaxy pairs in the selected tomographic samples. The measurements on the KiDS-like samples are carried with TREECORR in 9 logarithmically spaced $\theta$-bins using min/max separation angles of 0.5 and 475.5 arcmin, and a bin-slope parameter set to 0.01. In contrast, the Euclid-like measurements are performed with ATHENA in 10 logarithmically spaced $\theta$-bins spanning the range [0.1 - 300] arcmin, with a tree opening angle set to the ATHENA recommended value of 1.7 deg.

We also explore the effect of IA on three non-Gaussian statistics constructed from aperture mass maps $M_{ap}(\theta)$ (Schneider 1996). Closely following the measurement methods described in Martinet et al. (2021b, a), these maps are constructed directly from the galaxy catalogues by stacking the tangential component of ellipticities around pixels at position $\theta$, weighted by an aperture filter $Q(\theta, \theta_{ap}, x_c)$:

$$M_{ap}(\theta) = \frac{1}{n_{gal}(\theta) \sum_a W_a} \sum_a \epsilon_{a,t}(\theta, \theta_a) Q(|\theta - \theta_a|, \theta_{ap}, x_c) W_a. \quad (13)$$

$n_{gal}(\theta)$ is the galaxy density in the filter centred at $\theta$, $\theta_a$ is the position of galaxy a, and the tangential ellipticity with respect to the aperture centre is computed as $\epsilon_{a,t}(\theta, \theta_a) = -[\epsilon_1(\theta_a) \cos(2\phi(\theta, \theta_a)) + \epsilon_2(\theta_a) \sin(2\phi(\theta, \theta_a))]$, where $\phi(\theta, \theta_a)$ is the angle between both coordinates. We choose the same aperture filter as the previous two references:

$$Q(x) = \frac{\tanh(x/x_c)}{x/x_c} [1 + \exp(6 - 150x) + \exp(-47 + 50x)]^{-1}, \quad (14)$$

which is optimized for detecting haloes following an NFW profile (Schirmer et al. 2007). Here again, $x = \theta/\theta_{ap}$, $\theta$ is the distance to the filter centre, $x_c = 0.15$ and $\theta_{ap}$ is set to 12.5 arcmin for Stage III mocks, and to 10 arcmin for Stage IV mocks, as this parameter is also optimized against shape noise.

The local variance in these maps is evaluated from the magnitude of the ellipticities:

$$\sigma_{ap}^2(\theta) = \frac{1}{2n_{gal}^2(\theta) (\sum_a W_a)^2} \sum_a |\epsilon_a|^2 Q^2(|\theta - \theta_a|, \theta_{ap}, x_c) W_a^2, \quad (15)$$

allowing us to construct signal-to-noise maps $S/N(\theta) = M_{ap}(\theta)/\sigma_{ap}(\theta)$ in which we identify and count, in bins of $S/N$:

(i) peaks ($N_{peaks}$), defined as pixels with $S/N$ values larger than their eight neighbours,

(ii) pixel values ($N_{pix}$), in other words the full PDF,<fnref n="15" />

(iii) minima ($N_{min}$), defined as pixels with $S/N$ values smaller than their eight neighbours.

These make up the three aperture map-based (non-Gaussian) data vectors investigated in this paper; we focus our attention on the peak statistics at first and come back to minima and PDF later on. Closely following Martinet et al. (2021b), these measurements are conducted on each of the five tomographic bins, but also on joint catalogues constructed from the union between multiple bins, i.e. $1\times2, 1\times3... 4\times5$, since these 'cross-bins' contain a significant amount

of additional information that is not captured by the 'single-bins tomography' configuration (referred to as the 'auto-bins'). For the Euclid-like measurements we additionally include the union of more than two bins, e.g. $1\times2\times3, 1\times2\times4... 1\times2\times3\times4\times5$. Finally, it is shown in Martinet et al. (2021b) that equations (13) and (15) can be evaluated with fast convolutions on ellipticity grids as opposed to performing the sum over the exact galaxy positions in the aperture. This comes at a negligible cost in accuracy, and we exploit this speed-up approach for the analysis of our Euclid-like catalogues. Note that in this case we directly re-use the measurements from the SLICS and the cosmology simulations presented in Martinet et al. (2021b) to estimate the covariance matrix and model the cosmology dependence, which by design have the same survey properties as the Stage-IV IA-infused mocks introduced in this paper.

# 4 RESULTS: VALIDATION WITH $\gamma$-2PCFS

This section presents a comparison between the $\gamma$-2PCFs measurements from our NLA-infused KiDS-like catalogues and the analytical predictions.

## 4.1 GG, GI, and II

We show in Fig. 3 the different components of the $\gamma$-2PCF, namely the GG, GI, and the II terms, measured from the KiDS-like catalogues. The fiducial infusion model adopted here assumes a smoothing scale of $0.1 \ h^{-1} \ M_{\odot}$, a value that we justify in the next section and which corresponds to a mass of $4.8 \times 10^9 \ h^{-1} \ M_{\odot}$ at redshift zero. The shape noise is switched off in this measurement, easing the comparison with the theoretical model. As shown in Harnois-Déraps et al. (2019, see fig. 6 therein), the GG signal (black squares) matches the theoretical predictions to a few per cent at most angular scales, with some power loss at small physical scales due to limitations in the mass resolutions of the N-body runs. This is mostly visible in the $\xi_{-}$ panel for $\theta < 3$ arcmin, and hence these scales should be recalibrated or excluded when entering a cosmological inference based on real data. The GI and II terms (blue triangles and green circles, respectively) also show a similar overall agreement, undershooting the theory line only for $\theta < 10$ arcmin in $\xi_{-}$.

The II and GI terms, the later of which dominates the contamination within the LNA and $\delta$-NLA models, are accurately captured.<fnref n="16" />

While the II contribution increases the overall signal, the GI contribution typically removes power since the intrinsic alignment of galaxies and cosmic shear are anticorrelated. Moreover, the relative importance of the IA contribution to the cosmic shear signal rapidly decreases towards higher redshift, as expected from the NLA model: the II term becomes weaker there, while at the same time the GG contribution increases by about an order of magnitude. Therefore, the II signal is only important for the lowest redshift bins in autocorrelation, while the GI contribution is the largest in the cross-redshift bins, when correlating a mix of very high and very low redshift galaxies.

## 4.2 $\xi_{\pm}^{\text{IA}}$ versus $\xi_{\pm}^{\text{noIA}}$

Fig. 4 shows the theoretical predictions for the fractional difference between the $\gamma$-2PCFs with and without IA – the different types of black lines represent the NLA model with $A_{\text{IA}} = 1.0, 1.5,$

<footnote n="16">¹⁶This statement will need to be revisited for more complex scenarios, since it was shown in Secco et al. (2021) that the II term can be dominant in the TATT model.</footnote>

<footnote n="15">¹⁵These PDF are sometimes referred to the 'lensing PDF' or the 1D distribution in the literature.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Different components of the γ-2PCFs $\xi_+^{ij}$ (upper right-hand panels) and $\xi_-^{ij}$ (lower left-hand panel), for different tomographic bin combinations (i, j). The black squares show the measurements from 50 noise-free light cones without IA, while the blue triangles and green circles show the −GI and II terms measured from the IA-infused simulations, respectively. The $A_{IA}$ parameter is set to 1.5 in the simulations, while the dashed, solid, and dotted line show the NLA predictions for $A_{IA} = 1.0, 1.5$, and $2.0$, respectively. The tidal fields are smoothed with a Gaussian kernel with a smoothing scale of $\sigma_G = 0.1~h^{-1}~\text{Mpc}$. In this figure, the error bars indicate the error on the mean; these are difficult to distinguish here given their small sizes, but become more apparent later on when presenting ratios.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Fractional difference between the shear 2PCFs with and without the IA signal, $\xi_{IA}^+/\xi_{\pm}^{\text{noIA}} = 1$, in the NLA model. The theory lines show the predictions with three values of the $A_{IA}$ parameter indicated in the legend, while the measurements are for $A_{IA} = 1.5$ with different smoothing scales: the magenta squares are our default model, with $\sigma_G = 0.1~h^{-1}~\text{Mpc}$, while the other coloured curves show $\sigma_G = 0.5$ (green), 0.25 (blue), and $0.0~h^{-1}~\text{Mpc}$ (purple). The error bars show the error on the mean measured from 50 light cones. Note the different scaling of the y-axis. As for Fig. 3, no shape noise is included in this figure, to better exhibit the IA signature. The thin horizontal line represents the $A_{IA} = 0.0$ case.</figcaption></figure>
and 2.0. We also overplot in Fig. 4 the measurements from our simulations assuming $A_{IA} = 1.5$. The magenta squares show the fiducial smoothing scale, while the results obtained with different values of $\sigma_G$ are presented with the other coloured symbols. There is an excellent agreement with the theory line at large angular separation, with some deviations seen in the highly non-linear scales, notably for $\theta < 10$ arcmin in $\xi_-$, and $\theta < 3$ arcmin in $\xi_+$. These are the same scales for which the *II* and *GI* were found to be inaccurate in the previous section, due to limits in the mass resolution of the simulations. As $\sigma_G$ grows, the impact of IA at small-scale is increasingly erased, causing the measurements to reconnect with the no-IA case as $\theta \to 0.0$ arcmin, which is mostly visible in $\xi_-$ for the $\sigma_G = 0.5 \ h^{-1}\ \text{Mpc}$ case (green symbols). The agreement with theory is the best for the case without smoothing (purple lines in Fig. 4), however for physical reasons it is unlikely that the tidal fields on scales smaller than 100 kpc could influence the galaxy orientations, which themselves extend over tens of kpc. We have selected $\sigma_G = 0.1 \ h^{-1}\ \text{Mpc}$ for our fiducial model as it traces well the underlying theory predictions for most scales of interest, a value that could be revisited in the future. There is a visible vertical offset between the simulations and the model for bins 1-1 and 2-2, which is due to a slight underestimation of the *II* term in the simulations, as reported above. This is expected to have only a small impact on the cosmological inference, given that the signal is maximal at higher redshifts, and that only two out of the 15 tomographic configurations are affected. The source of this underestimation is unclear, it is likely due to the large projection length (257.5 $h^{-1}\ \text{Mpc}$) along the redshift direction involved in the construction of our tidal field maps. Indeed, we could possibly have obtained a better match had we computed the tidal fields from the full particle distribution instead or finer redshift sampling, however these are not available to us and hence we leave this verification for future work.

It should also be mentioned that the smoothing prescription is a physical criterion that does not appear explicitly in the theoretical NLA model, aside from optional cuts in k-modes that are allowed to contribute to the $P_{II}(k)$ and $P_{GI}(k)$ spectra, as in Blazek et al. (2015). We do not include such selections here.

### 4.2.1 Varying $A_{IA}$

The final validation step consists in verifying that our infusion model correctly responds to changes in the coupling strength between the intrinsic galaxy shapes and the tidal fields, as controlled by the $A_{IA}$ parameter. This is tested by consistently modifying the parameter value in the model (equation 6) and at the infusion stage (equations 4 and 5). We present our results for $A_{IA} = 1.0, 1.5$, and $2.0$ in Fig. 5, where we observe that the scaling of the signal with $A_{IA}$ is as expected, with some deviations mostly in the autotomographic bins, and a good overall agreement. Achieving an accurate $A_{IA}$ scaling is an important milestone: we can hereafter vary this parameter in the simulated lensing catalogues and use equation (11) to construct new mock observations, from which we can interpret cosmic shear data with a flexible IA sector, for any weak lensing estimator. This is therefore compatible with joint [γ-2PCFs; Peaks] analyses as in Martinet et al. (2018) and Harnois-Déraps et al. (2020), but now augmented with a new joint $A_{IA}$ marginalization capacity.

Of course, this presupposes that the infusion model adequately describes the physical IA at play in the real Universe, which is still highly uncertain: according to some recent studies (e.g. Blazek et al. 2015; Troxel et al. 2018; Fortuna et al. 2021) the standard

NLA model is disfavoured over more complex models, whereas Secco et al. (2021) find that although consistent, the TATT model is unnecessarily flexible for the analysis of the DES-Y3 data and in fact degrades the cosmological constraints compared to the NLA case. If required, additional features can be infused in equation (6) including redshift scaling (e.g. with the $\eta_{IA}$ parameter, see equation 15 of Asgari et al. 2021), or even completely different physical models that can be constructed from the tidal tensor such as the TATT. One of such extensions has been introduced in Section 2.3, the δ-NLA model, and results in significantly different IA signatures at small scales, shown by the green lines in Fig. 5. This has been computed by fixing the smoothing length to $0.1 \ h^{-1} \ \text{Mpc}$ in both models, a choice that might be revisited in the future, but nevertheless the deviations from the NLA are clear and visible at all scales.

## 5 RESULTS: IMPACT ON $M_{\text{AP}}$ STATISTICS

### 5.1 Peak count + IA: Impact on Stage-III surveys

Having validated our IA infusion model with the γ-2PCFs in the last section, we now turn our investigation towards the aperture mass statistics defined in Section 3.4. With the shape noise turned on hereafter, we construct $S/N$ maps and measure the lensing peak function, $N_{\text{peaks}}$, in the KiDS-like catalogues. We carry out this measurement for all tomographic bins and their pair-wise combinations, with and without the IA. Following Harnois-Déraps et al. (2020), we use $S/N$ bins sizes of 1/3, here over the wider range of $[-3 \le S/N \le 4]$. We consider hereafter the δ-NLA model with $\sigma_G = 0.1 \ h^{-1} \ \text{Mpc}$ as our fiducial case, being more realistic than the standard NLA, but discuss some NLA results as well. The multiple panels in Fig. 6 illustrate the different contributions to the observed aperture mass map. The top panels show the $S/N$ maps from tomographic bins 1, 5, and $1 \cup 5$; the middle panels map out the aperture number count $N_{\text{ap}}$, obtained by setting the ellipticities to unity in equation (13); the bottom panels present the IA contamination map, $M_{\text{ap}}^{\text{IA}}$, present in the observed $S/N$ shown in the top panels. We can compute $M_{\text{ap}}^{\text{IA}}$ either from the difference between the measured $M_{\text{ap}}$ with and without IA, or by passing the pure IA signal to (equation 13) instead of the observed ellipticities. Note the strong correlation between $N_{\text{ap}}$ and $M_{\text{ap}}^{\text{IA}}$, expected from the coupling between the IA ellipticities, the tidal field, and the mass overdensity. The relative impact of IA on the observed $M_{\text{ap}}$ is only a few per cent over most of the pixels and is therefore not visible by eye. As for two-point correlations, the IA contribution is higher in low redshift slices (bin 1 compared to bin 5), which both contributes in the combination of slices (bin $1 \cup 5$). The effect on the aperture mass statistics further depends on the noise level, which is the lowest for combined bins.

When counting peaks with and without the IA however, we notice a stronger effect on the high $S/N$ peaks, as shown in Fig. 7. This suppression is particularly strong in the cross-tomographic bins, causing in some cases a 10 per cent suppression of peaks with $S/N = 4.0$. This can be understood by the fact that high $S/N$ peaks are found in regions of high foreground density, where substantial alignment of galaxies is produced, leading to a large contribution from a GI-like term.<fnref n="17" /> There is also a noticeable effect on peaks

<footnote n="17">¹⁷The expressions *GG*, *II* and *GI* strictly apply only to two-point functions, however we can generalize this concept for other types of measurements. Here, we use the expression *GI*-like to describe the contribution to the total peak count statistics coming from cosmic shear signal (*G*) of high-redshift galaxies, with the intrinsic signal (*I*) from the low-redshift ones.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Similar to Fig. 4, but here the strength of the $A_{IA}$ parameter is varied in the mocks, taking values of 2.0 (magenta dotted lines), 1.5 (magenta symbols with solid lines), and 1.0 (magenta dashed lines) – the same as those assumed for the three NLA theory lines (shown in black). The smoothing scale is fixed to $\sigma_G = 0.1 \, h^{-1} \, \text{Mpc}$. Again, no shape noise is included here, and the error is about the mean, measured from 50 light cones. The green line shows the results from the infused $\delta$-NLA model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fractional difference between the number of lensing peaks with and without the IA signal, measured from the KiDS-like simulations. Here again, the measurements are obtained after smoothing the tidal field with a smoothing scale of $0.1 \, h^{-1} \, \text{Mpc}$. The error bars about the yellow symbols show the error on the mean, measured from 50 light cones with shape noise, for $A_{IA} = 1.0$ (green), 1.5 (yellow), and 2.0 (red).</figcaption></figure>

<caption>**Figure 6.** (upper panels:) $S/N$ maps, constructed from different tomographic bins in one of the KiDS-1000-like simulations, including intrinsic alignments in the galaxy orientations computed from the tidal fields shown in Fig. 2, for a smoothing length of $0.1 \, h^{-1} \, \text{Mpc}$ and $A_{IA} = 1.5$. (middle panels:) Aperture number count map, constructed from the same galaxy samples. (bottom panels:) Contamination from the intrinsic alignment of the galaxies (shown in the middle panels) to the aperture mass maps (shown in the upper panels). Note that the upper and lower panels do not show the same quantity, whence the differences in the colour scale.</caption>

with large negative $S/N$ values, which is not surprising given that these are highly correlated with the high $S/N$ peaks (Martinet et al. 2018). Similarly to the $\gamma$-2PCFs, the cross-redshift bins are more severely affected by IA whereas the autotomographic bins are relatively immune; this could prove helpful in a data analysis for choosing elements of the measurement vector that are less affected

by uncontrolled IA systematics (as done in Harnois-Déraps et al. 2020, and further explored in Section 5). Although this conclusion certainly depends on the alignment model in place, one would still expect the cross-redshift bins to be most affected because of the $GI$-like terms. For instance, impact from the NLA model is shown by the black lines in Fig. 7, which exhibits a milder effect than the $\delta$-NLA
model, consistent with the γ-PCFs. Alternatively, if, for example, a non-perturbative IA model (e.g. Fortuna et al. 2021) was to be considered instead of the current δ-NLA, much larger alignments around haloes could occur; the same applies also for higher order operators of the tidal field, however, we leave this investigation for future work.

Changing the smoothing scale from 0.1 to 0.5 $h^{-1}$ Mpc has almost no effect of the results presented in Fig. 7. This is likely caused by the additional smoothing inherent to aperture map statistics (see equation 13), which suppresses features smaller than the filter scale. These two smoothing scales correspond to angular scales of less than 2 arcmin (except for at the very lowest redshifts, where no lensing signal originate anyway) and are therefore not well distinguishable in the aperture mass maps. Note that this conclusion will not hold for smaller aperture filter sizes however, for which we expect to recover stronger IA effects for smaller $\sigma_G$.

When it comes to understanding what happens to individual peaks in presence of IA, Kacprzak et al. (2016) suggest that the suppression of large peaks arises from foreground satellite galaxies in-falling on to massive foreground clusters, but being incorrectly included in the background source galaxy sample. Inspired by Schneider & Bridle (2010), their model assumes that the radial alignment of these cluster-member satellites would tend to undo the tangential cosmic shear signal imparted by the said cluster on to the true background galaxies. This hypothesis is however not supported by the observational analysis of Sifón et al. (2015) who measured negligible radial alignment of satellite galaxies around galaxy clusters. In addition, Fig. 7 reveals that IA affect the peak function even in absence of cluster member contamination, since this is not explicitly included in our simulations. To be precise, our method introduces the equivalent of a small local contamination, since all simulated galaxies that populate a given mass sheet receive the same IA contribution. This, however, would not explain the IA effect observed in cross-tomographic bins, which are well separated in redshift; the presence of a GI-like term is a more accurate explanation.

Another difference with the Kacprzak et al. (2016) model is that while the fractional effect is the largest on large positive and negative $S/N$ peaks, we find that the absolute effect is larger for peaks with $S/N \sim [0.0 - 1.0]$, as seen in Fig. 8. The bottom panel of this figure shows the difference $N_{\text{peaks}}^{\text{IA}} - N_{\text{peaks}}^{\text{noIA}}$ instead of the ratio, for the auto-tomographic bin 2 and the combination $2\cup5$. The impact is negligible in the former but strong in the latter, and demonstrates that IA tend to increase the number of peaks in the range $-0.5 < S/N < 2$ and suppress those outside that range. The ‘bubble plot’ in the upper panel shows the peak positions found in the $2\cup5$ catalogue with (red) and without IA (black), where the size of each bubble scales with the peak’s $S/N$ value. We observe a clear correlation between the large red and black circles, indicating that both cases identify common foreground overdensities with high significance. However the largest black circles generally have their red counterparts slightly smaller in size, corresponding to a reduction in $S/N$ observed in the peak function. Moreover their positions are slightly shifted, but not as much as for the smaller peaks, for which there are many red circles with no black counterpart and vice versa. Because of their increased number, the low-$S/N$ peaks are in an absolute sense maximally affected by IA, as already seen in the lower panel of Fig. 8, but relatively speaking, these only affect a small fraction of the peaks. We further note that although this detailed inspection is presented here only for the bin combination $2\cup5$, others cross-bins show similar but milder features, while in the autotomographic cases the relative effect of IA is consistent with noise, as shown in the diagonal panels of Fig. 7.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** (upper:) Position and height of the peaks measured in the tomographic bin combination $2\cup5$, with $A_{\text{IA}} = 2.0$ (red) and 0.0 (black), inside a 1.0 deg$^2$ patch. The four concentric circles indicate the sizes for $S/N = 1, 2, 3$ and 4. (lower:) Difference between the peak function with and without IA as measured from our 50 KiDS-like catalogues, for the tomographic bin 2 and the combination $2\cup5$, shown for $A_{\text{IA}} = 1.0$ (green), 1.5 (yellow) and 2.0 (orange).</figcaption></figure>

We can conclude from this investigation that the impact of IA is non-negligible in the peak count analyses of Stage-III lensing survey as it suppresses systematically the high-$S/N$ end of data vector, which is also highly sensitive to cosmology (Martinet et al. 2018). With their cluster-member contamination model, Kacprzak et al. (2016) reports a small 1.5 per cent shift in $S_8$ due to IA in a non-tomographic setting, which was also adopted in the KiDS-450 peak count analysis of Martinet et al. (2018). The tomographic DES-Y1 analysis of Harnois-Déraps et al. (2020) modelled the IA with a halo occupation distribution model, aligning the central galaxies with the shape of the host dark matter haloes, without any mis-alignment scatter. They find a sub-per cent impact on $S_8$, but their model includes IA only in the bins 1, 2, and $1\cup2$ (out of four bins), albeit with a stronger relative effect on the peak function. In light of the IA infusion model investigated in this paper, it is clear that all cross-tomographic bins are affected by the IA, which suggests that the three analyses abovementioned could have underestimated the impact of IA. The two earlier analyses were carried in a non-tomographic setting and therefore are as
vulnerable to IA contamination, as we demonstrate in the next section.

Comparing the amplitude of the impact with γ-2PCFs, we expect that completely neglecting the IA in a peak count case is less catastrophic: whereas the IA affects ξ_± at all angular bins by sometimes 50–100 per cent (Fig. 5), only the largest peaks are affected here, and by no more than 5–10 per cent (Fig. 7). We show in the next section how this affects the inferred cosmological constraints for the upcoming generation of lensing surveys.

## 5.2 Peak count + IA: Impact on Stage-IV surveys

With highly improved statistics, Stage-IV lensing surveys will be increasingly sensitive to systematics and secondary signals. In this section, we aim to forecast the impact of IA on a cosmic shear measurement based on 100 square degrees of Euclid-like data, analysed with the peak count statistics, γ-2PCFs, and in the joint case. As detailed in Section 3.2, we use the cosmo-SLICS catalogues to model the cosmology dependence of the signal and the 928 SLICS catalogues to estimate the joint covariance matrix. The likelihood is modelled as a multivariate t-distribution:

$$ \mathcal{L}(\pi|\mathbf{d}) \propto \frac{N_{\text{sim}}}{2} \ln \left[ 1 + \chi^2/(N_{\text{sim}} - 1) \right], \quad \text{with} \\ \chi^2 = \sum [\mathbf{x}(\pi) - \mathbf{d}]^\mathrm{T} \mathrm{Cov}^{-1} [\mathbf{x}(\pi) - \mathbf{d}], \quad (16) $$

where d and x(π) are, respectively, the measurement and model vectors, π is the cosmological point at which the model is evaluated, and Nsim = 928. As explained in Sellentin & Heavens (2016), this choice of likelihood accurately captures the residual noise caused by the inversion of a covariance matrix estimated from a finite number of independent simulations.<fnref n="18" /> With at least hundreds of peaks in each of our S/N bins, adopting this form is justified. The likelihood is sampled within a parameter range that is determined by the cosmo-SLICS training set:

$$ \begin{aligned} 0.10 &\le \Omega_m \le 0.55 \\ 0.60 &\le S_8 \le 0.90 \\ -2.0 &\le w_0 \le -0.5 \\ 0.60 &\le h \le 0.82 \end{aligned} \quad (17) $$

and is evaluated by interpolating the measurements at the cosmo-SLICS nodes with radial basis functions. We follow the approach developed for hydrodynamical simulations in Martinet et al. (2021a) and infuse the bias measured in the IA catalogues to the measurement data vector. This ensures that differences between the mocks used for the wCDM model and those used for IA infusion (for instance, the different sampling of galaxy positions) are not biasing our model, nor our estimate of the resulting bias on cosmology. The interpolation is performed in a two-step approach on a regular grid of 40 points for each parameter, the second step refining the parameter step to the hyperspace where the likelihood is non-zero.

The relative effect of IA on peak statistics is measured from our Euclid-like simulations in eight bins spanning the range [−2.5 ≤ S/N ≤ 5.5], and shown in the upper panel of Fig. 9. The trends seen in the KiDS-like data are similar and amplified here, with elements in the cross-terms bins affected by up to 40 per cent (see e.g. the bin 1U3). Once again, IA are the strongest at low redshift and in cross-tomographic bins. Overplotted on these measurements are the effect

<figure><img src="image_10.png" /><figcaption>**Figure 9.** (upper:) Same as Fig. 7, but measured from the 10 Euclid-like simulations, each averaged over 5 noise realizations. The thin black lines indicate the impact from the baryonic feedback reported in Martinet et al. (2021a), which generally affects the signal more severely than the IA except at the lowest redshifts. The different coloured lines show the effect of varying AIA. (middle and lower:) Same as upper panel, but for lensing PDF (Npix) and minima (Nmin), respectively. A smoothing of 0.1 h⁻¹ Mpc is used here, however these curves are mostly unchanged if employing σG = 0.5 h⁻¹ Mpc instead. Note the change in the y-axis scaling for the lowest redshift bin.</figcaption></figure>

<footnote n="18">¹⁸The t-distribution likelihood reduces to a multivariate Gaussian likelihood in the limit of infinite number of simulations.</footnote>
<table><caption>**Table 2.** Different survey configurations investigated with Stage-IV cosmology inference pipeline. For each of these, we consider three aperture mass statistics $N_\alpha(S/N)$, namely peaks ($\alpha$ = peaks), minima ($\alpha$ = min), and lensing PDF ($\alpha$ = pix).</caption><thead><tr><th rowspan="2">case</th><th colspan="4">N<sub>α</sub>(S/N)</th><th colspan="4">γ-2PCFs</th><th rowspan="2">N<sub>d</sub></th></tr><tr><th>no-tomo</th><th>auto</th><th>pairs</th><th>all</th><th>no-tomo</th><th>auto</th><th>cross</th></tr></thead><tbody><tr><td>(i)</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>8</td></tr><tr><td>(ii)</td><td>-</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>40</td></tr><tr><td>(iii)</td><td>-</td><td>√</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>120</td></tr><tr><td>(iv)</td><td>-</td><td>√</td><td>√</td><td>√</td><td>-</td><td>-</td><td>-</td><td>200</td></tr><tr><td>(v)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>√</td><td>-</td><td>-</td><td>16</td></tr><tr><td>(vi)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>√</td><td>-</td><td>80</td></tr><tr><td>(vii)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>√</td><td>√</td><td>240</td></tr><tr><td>(viii)</td><td>√</td><td>-</td><td>-</td><td>-</td><td>√</td><td>-</td><td>-</td><td>24</td></tr><tr><td>(ix)</td><td>-</td><td>√</td><td>√</td><td>-</td><td>-</td><td>√</td><td>√</td><td>360</td></tr><tr><td>(x)</td><td>-</td><td>√</td><td>√</td><td>√</td><td>-</td><td>√</td><td>√</td><td>440</td></tr></tbody></table>

of baryonic feedback measured from the Magneticum simulations in Martinet et al. (2021a), which unlike IA is approximately constant for all redshift bins. Since the cosmological information extracted by lensing is the highest at high redshift, we could expect peak count statistics from deep surveys to be more affected by baryonic feedback than IA, however this conclusion could of course vary when considering models different from the δ-NLA and magneticum.

Fig. 9 also presents the results for lensing PDF ($N_{pix}$) and minima ($N_{min}$) measured from the same simulations and $S/N(\theta)$ maps, with the exception that the former uses nine bins in the range $[-4 \le S/N \le 5]$, and the latter eight bins with $[-5 \le S/N \le 3]$, as in Martinet et al. (2021b). While the lensing PDF responds to IA in a manner almost identical to the peaks, the lensing minima seems to be more immune, with at most a 10 per cent effect even in the lowest redshift bins. We could therefore expect that lensing minima would incur a smaller cosmological bias from unmodelled IA, since these are the strongest in high-density regions, which are mostly avoided by this statistics. We show next that although this reasoning has some logic to it, it does not always hold in practice.

Combining the data vectors, the wCDM models and the covariance matrices, we now forecast Stage-IV parameter constraints and compare the marginalized a posteriori distributions after running the likelihood pipeline on data, with and without the alignments. Any observed difference can be interpreted as a bias on the inferred cosmology caused by the inclusion of IA. We explore different analysis strategies including aperture mass statistics and γ-2PCFs in a variety of tomographic choices, for the purpose of identifying the configurations that are least affected. As in Martinet et al. (2021b), we reject the two largest 9-bins in $\xi_+$ to protect us from residual finite simulation box effects, and further remove the two smallest 9-bins in $\xi_-$ as they are subject to uncertain non-linear physics. For each aperture mass statistics $N_\alpha(S/N)$, with $\alpha$ = 'peaks', 'min' or 'pix', we consider four distinct cases summarized in Table 2 (cases i-iv). We further compare these results to the γ-2PCFs with cases (v-vii), and the joint probes [γ-2PCFs; $N_\alpha(S/N)$] in cases (viii-x).

The results for cases (i), (ii), (iii), (iv), and (vii) with the peak count statistics are shown in Fig. 10. As first found in Martinet et al. (2021b), tomographic peak count analyses can improve by a factor of three the precision on $w_0$, and by a factor of 2.0 that of $S_8$, in absence of systematics. This can be seen by comparing the size of the marginalized posterior distribution in the different panels of Fig. 10, more precisely the unbiased cases, and from comparing the marginal constraints reported in Table 3 in the 'no IA' cases. We expect that the inclusion of additional systematics (e.g. photometric

redshifts, baryons) will further deteriorate both γ-2PCFs and peaks constraints, but this can be achieved without biasing the inferred value if modelled correctly.

We draw the attention on the difference between the biased (green contours) and the unbiased (blue contours) cases, which illustrate the impact IA could have on the inferred cosmology if left unmodelled, assuming $A_{IA} = 1.5$. There is a clear shift towards lower $S_8$ and $\Omega_m$ values in most cases, as reported in Table 3. However, these are mild for peak statistics, in comparison with the bias experienced by the γ-2PCFs ($-1.0$ per cent and $-6.3$ per cent, versus $-8.8$ per cent and $-31.1$ per cent, for $S_8$ and $\Omega_m$ respectively). It is encouraging to see that except for the full tomographic case, the dark energy parameter $w_0$ is less affected by the IA. Whether this secondary spike observed in the full tomography is a real IA-$w_0$ degeneracy or noise caused by our relatively coarse cosmological sampling remains to be demonstrated in future work based on the next generation of wCDM weak lensing simulations, however it is consistent with the shift seen in the γ-2PCFs (case vii). Blazek et al. (2019) finds similar results, i.e. that unmodelled IA tends to lower $S_8$ by as much as $\Delta S_8 = 0.1$, and the dark energy equation of state by $\Delta w_0 = 0.5$. Note that their setup also include a redshift evolution, namely $w(a) = w_0 + w_a(1 - a)$, and $w_a$ is shifted by $7\sigma$ to $-3.0$, compared to the input value of 0.0.

We note from Fig. 10 that none of the projected posteriors are prior dominated, however the residual noise observed in the shape is caused by the sparse cosmological sampling and could be improved with more cosmo-SLICS nodes.

The results for all cases and all estimators are also summarized in Table 3 (see Appendix B for an equivalent table measured from different $A_{IA}$ values). In the full tomographic joint [γ-2PCFs; $N_{peaks}$] analysis explored in this paper (case x), an unmodelled δ-NLA IA signal would bias $S_8$ by $\sim 7$ per cent, $\Omega_m$ by $\sim 40$ per cent, and $w_0$ by $\sim 53$ per cent, driven by both the γ-2PCFs and peaks. In the joint [γ-2PCFs; $N_{min}$] analysis, $S_8$, $\Omega_m$, and $w_0$ could be biased by 9.6 per cent, 43.7 per cent, and 52.8 per cent, respectively, and similar results are obtained from the [γ-2PCFs; $N_{pix}$] case (with 4.5 per cent, 41.1 per cent, and 57.0 per cent). Contrary to our expectations, lensing minima alone show a similar bias in $S_8$ compared to peaks, but the lensing PDF alone are relatively immune, with a shift of less than 1 per cent. Most non-Gaussian statistics show a secondary peak develop in $w_0$, which can bias the measurement beyond 50 per cent. The effect on $\Omega_m$ are the smallest for the lensing PDF, reaching 23.8 per cent (compared to 34.1 per cent for minima and 38.9 per cent for peaks). As the reported biases correspond to the shift in the maximum of the likelihood between the IA and no IA cases, part of the bias could be due to the sparsity of the cosmological parameter space, although it is dense enough to get the right order of magnitude. We also observe that the inclusion of IA affects the size of the marginalized errors, as the best-fitting values lie in a different region in parameter space, therefore changing the relative importance of signal and shape noise.

Comparing now the bias ($\Delta$) to the precision ($\delta$) on the parameters in Table 3, we note that the shifts are ranging from the 0–15$\sigma$, the worst cases being the error on $w_0$ in the joint full tomographic settings, where the precision is the highest. We notice that the posterior distribution about this parameter are in this case rather noisy, developing a secondary peak, which inflates significantly this shift. Nevertheless, all parameters are significantly affected by IA, with levels that vary with the choice of estimator and tomographic bin, and no case is fully protected. The average biases on $S_8$, $w_0$, and $\Omega_m$ are $2.2\sigma_{stat}$, $2.8\sigma_{stat}$, and $4.2\sigma_{stat}$, the latter being systematically severely affected.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Bias on the *w*CDM cosmological parameters from analysing 100 deg² of IA-infused *Euclid*-like data with models calibrated on simulations without IA. The upper four panels are for peak without tomography (upper left-hand panel), with auto-tomographic bins (upper right-hand panel), with autotomographic bins + bin pairs (middle left-hand panel) and with all bin combinations (middle right-hand panel). The full tomographic 2PCFs analysis (lower left-hand panel) and joint analysis (lower right-hand panel) are significantly more affected by IA in comparison, especially in *w*₀. The red lines and '+' indicate the input cosmology while the grey '+' symbols indicate the *cosmo*-SLICS training nodes.</figcaption></figure>

Note that scaling these results to a full *Euclid* area of 15 000 deg² would improve the precision ‘δ’ by a factor of √15 000/100 ~ 12.2 compared to the current estimates, while leaving the IA-induced bias unchanged. In that case, a 0.3σ shift would become a 3.6σ bias, which already is catastrophic for the inference.

To be highlighted, the autotomographic analysis of the lensing minima is mostly immune to the IA, with no shift recorded in $S_8$ nor $w_0$, and a 4 per cent shift in $\Omega_m$. Another robust case is the autotomographic peaks case, where $\Delta S_8$ is about a tenth of $\delta S_8$, $\Delta w_0$ ~ 0.05σ, but $\Delta\Omega_m$ ~ 0.3σ already. These measurements are already
<table><caption>**Table 3.** Bias on the inferred cosmology caused by the intrinsic alignments of galaxies, if left unaccounted for, for the different cases listed in Table 2. Parameter shifts are reported as ‘Δ’, and can be contrasted with the precision expected for a 100 deg² of Euclid-like data, reported as $δ ≡ σ_{stat}$. Numbers in parenthesis show the shifts Δ in units of statistical errors δ, and the percentage of the marginalized error with respect to the input parameter value. An $A_{IA} = 1.5$ and a smoothing scale of $0.1 h^{-1}$ Mpc are assumed here. The cases labelled as ‘M20’ are extracted from table 1 in Martinet et al. (2021b).</caption><thead><tr><th>case</th><th>estimator(s)</th><th>ΔS<sub>8</sub></th><th>Δw<sub>0</sub></th><th>ΔΩ<sub>m</sub></th><th>δS<sub>8</sub></th><th>δw<sub>0</sub></th><th>δΩ<sub>m</sub></th></tr></thead><tbody><tr><td>(i)</td><td>N<sub>peaks</sub>, no tomo</td><td>-0.052 (-1.4)</td><td>0.038 (0.1)</td><td>-0.007 (-0.1)</td><td>0.035 (4.2 per cent)</td><td>0.284 (28.4 per cent)</td><td>0.047 (16.2 per cent)</td></tr><tr><td>(ii)</td><td>N<sub>peaks</sub>, tomo auto</td><td>-0.003 (-0.1)</td><td>-0.009 (-0.0)</td><td>-0.01 (-0.3)</td><td>0.027 (3.3 per cent)</td><td>0.186 (18.6 per cent)</td><td>0.033 (11.3 per cent)</td></tr><tr><td>(iii)</td><td>N<sub>peaks</sub>, tomo auto + pairs</td><td>-0.042 (-1.2)</td><td>-0.606 (-1.7)</td><td>-0.103 (-2.7)</td><td>0.034 (4.1 per cent)</td><td>0.358 (35.8 per cent)</td><td>0.038 (13.3 per cent)</td></tr><tr><td>(iv)</td><td>N<sub>peaks</sub>, tomo full</td><td>-0.058 (-2.9)</td><td>-0.524 (-5.9)</td><td>-0.113 (-6.6)</td><td>0.02 (2.5 per cent)</td><td>0.089 (8.9 per cent)</td><td>0.017 (5.7 per cent)</td></tr><tr><td>(M20)</td><td>N<sub>peaks</sub>, tomo full, no IA</td><td>-</td><td>-</td><td>-</td><td>0.022 (2.7 per cent)</td><td>0.130 (13 per cent)</td><td>0.024 (8 per cent)</td></tr><tr><td>(v)</td><td>γ-2PCFs, no tomo</td><td>-0.042 (-1.0)</td><td>0.000 (0.0)</td><td>-0.069 (-1.1)</td><td>0.042 (5.1 per cent)</td><td>0.307 (30.7 per cent)</td><td>0.062 (21.2 per cent)</td></tr><tr><td>(vi)</td><td>γ-2PCFs, tomo auto</td><td>-0.006 (-0.2)</td><td>-0.037 (-0.1)</td><td>0.005 (0.1)</td><td>0.026 (3.2 per cent)</td><td>0.253 (25.3 per cent)</td><td>0.044 (15.1 per cent)</td></tr><tr><td>(vii)</td><td>γ-2PCFs, tomo auto + pairs</td><td>-0.083 (-4.4)</td><td>-0.878 (-15.1)</td><td>-0.124 (-11.6)</td><td>0.019 (2.3 per cent)</td><td>0.058 (5.8 per cent)</td><td>0.011 (3.9 per cent)</td></tr><tr><td>(M20)</td><td>γ-2PCFs, tomo auto + pairs, no IA</td><td>-</td><td>-</td><td>-</td><td>0.023 (2.8 per cent)</td><td>0.190 (19 per cent)</td><td>0.034 (12 per cent)</td></tr><tr><td>(viii)</td><td>joint, no tomo</td><td>-0.087 (-3.2)</td><td>-0.423 (-1.4)</td><td>-0.029 (0.6)</td><td>0.027 (3.2 per cent)</td><td>0.294 (29.4 per cent)</td><td>0.046 (15.8 per cent)</td></tr><tr><td>(ix)</td><td>joint, tomo auto + pairs</td><td>-0.067 (-4.5)</td><td>-0.553 (-5.2)</td><td>-0.118 (-9.8)</td><td>0.015 (1.9 per cent)</td><td>0.106 (10.6 per cent)</td><td>0.012 (4.0 per cent)</td></tr><tr><td>(x)</td><td>joint, tomo full</td><td>-0.057 (-4.1)</td><td>-0.527 (-6.0)</td><td>-0.116 (-11.6)</td><td>0.014 (1.7 per cent)</td><td>0.088 (8.8 per cent)</td><td>0.010 (3.3 per cent)</td></tr><tr><td>(M20)</td><td>joint, tomo full, no IA</td><td>-</td><td>-</td><td>-</td><td>0.015 (1.8 per cent)</td><td>0.090 (9 per cent)</td><td>0.017 (6 per cent)</td></tr><tr><td>(i)</td><td>N<sub>min</sub>, no tomo</td><td>-0.051 (-1.3)</td><td>0.000 (0.0)</td><td>0.000 (0.0)</td><td>0.040 (4.8 per cent)</td><td>0.273 (27.3 per cent)</td><td>0.092 (31.7 per cent)</td></tr><tr><td>(ii)</td><td>N<sub>min</sub>, tomo auto</td><td>0.000 (0.0)</td><td>0.000 (0.0)</td><td>0.012 (0.2)</td><td>0.029 (3.6 per cent)</td><td>0.195 (19.5 per cent)</td><td>0.049 (16.8 per cent)</td></tr><tr><td>(iii)</td><td>N<sub>min</sub>, tomo auto + pairs</td><td>-0.041 (-1.3)</td><td>-0.582 (-1.7)</td><td>-0.086 (-2.2)</td><td>0.032 (3.9 per cent)</td><td>0.347 (34.7 per cent)</td><td>0.04 (13.7 per cent)</td></tr><tr><td>(iv)</td><td>N<sub>min</sub>, tomo full</td><td>-0.056 (-2.5)</td><td>-0.578 (-10.7)</td><td>-0.099 (-3.8)</td><td>0.022 (2.7 per cent)</td><td>0.054 (5.4 per cent)</td><td>0.026 (9.1 per cent)</td></tr><tr><td>(M20)</td><td>N<sub>min</sub>, tomo full, no IA</td><td>-</td><td>-</td><td>-</td><td>0.024 (2.9 per cent)</td><td>0.120 (12 per cent)</td><td>0.027 (9 per cent)</td></tr><tr><td>(viii)</td><td>joint, no tomo</td><td>-0.088 (-3.1)</td><td>-0.037 (-0.2)</td><td>-0.024 (-0.5)</td><td>0.028 (3.4 per cent)</td><td>0.233 (23.3 per cent)</td><td>0.052 (17.8 per cent)</td></tr><tr><td>(ix)</td><td>joint, tomo auto + pairs</td><td>-0.069 (-3.8)</td><td>-0.506 (-7.1)</td><td>-0.129 (-10.8)</td><td>0.018 (2.2 per cent)</td><td>0.071 (7.1 per cent)</td><td>0.012 (4.3 per cent)</td></tr><tr><td>(x)</td><td>joint, tomo full</td><td>-0.079 (-5.6)</td><td>-0.528 (-6.0)</td><td>-0.127 (-11.5)</td><td>0.014 (1.7 per cent)</td><td>0.088 (8.8 per cent)</td><td>0.011 (3.7 per cent)</td></tr><tr><td>(M20)</td><td>joint, tomo full, no IA</td><td>-</td><td>-</td><td>-</td><td>0.014 (1.7 per cent)</td><td>0.080 (8 per cent)</td><td>0.017 (6 per cent)</td></tr><tr><td>(i)</td><td>N<sub>pix</sub>, no tomo</td><td>-0.063 (-2.1)</td><td>0.000 (0.0)</td><td>0.000 (0.0)</td><td>0.030 (3.6 per cent)</td><td>0.139 (13.9 per cent)</td><td>0.035 (12.1 per cent)</td></tr><tr><td>(ii)</td><td>N<sub>pix</sub>, tomo auto</td><td>-0.000 (0.0)</td><td>0.036 (0.2)</td><td>-0.015 (-0.5)</td><td>0.026 (3.2 per cent)</td><td>0.217 (21.7 per cent)</td><td>0.028 (9.8 per cent)</td></tr><tr><td>(iii)</td><td>N<sub>pix</sub>, tomo autopairs</td><td>0.016 (0.6)</td><td>0.039 (0.5)</td><td>-0.056 (-1.75)</td><td>0.029 (3.6 per cent)</td><td>0.075 (7.5 per cent)</td><td>0.032 (11.2 per cent)</td></tr><tr><td>(iv)</td><td>N<sub>pix</sub>, tomo full</td><td>0.008 (0.4)</td><td>0.052 (1.9)</td><td>-0.069 (-3.1)</td><td>0.021 (2.5 per cent)</td><td>0.028 (2.8 per cent)</td><td>0.022 (7.5 per cent)</td></tr><tr><td>(M20)</td><td>N<sub>pix</sub>, tomo full, no IA</td><td>-</td><td>-</td><td>-</td><td>0.021 (2.6 per cent)</td><td>0.100 (10 per cent)</td><td>0.019 (7 per cent)</td></tr><tr><td>(viii)</td><td>joint, no tomo</td><td>-0.081 (-3.1)</td><td>0.028 (0.23)</td><td>-0.003 (-0.1)</td><td>0.026 (3.2 per cent)</td><td>0.118 (11.8 per cent)</td><td>0.031 (10.7 per cent)</td></tr><tr><td>(ix)</td><td>joint, tomo auto + pairs</td><td>-0.060 (-4.3)</td><td>-0.596 (-3.7)</td><td>-0.093 (-8.5)</td><td>0.014 (1.6 per cent)</td><td>0.162 (16.2 per cent)</td><td>0.011 (3.7 per cent)</td></tr><tr><td>(x)</td><td>joint, tomo full</td><td>-0.037 (-3.4)</td><td>-0.570 (-6.0)</td><td>-0.119 (-13.2)</td><td>0.011 (1.3 per cent)</td><td>0.095 (9.5 per cent)</td><td>0.009 (2.9 per cent)</td></tr><tr><td>(M20)</td><td>joint, tomo full, no IA</td><td>-</td><td>-</td><td>-</td><td>0.013 (1.5 per cent)</td><td>0.060 (6 per cent)</td><td>0.015 (5 per cent)</td></tr></tbody></table>

competitive with the γ-2PCFs projected constraints, but nowhere near the combined-probe results, especially that of $[N_{pix}, \gamma-2PCFs]$, which can reach precision of 0.06 on $w_0$, 0.013 on $S_8$, and 0.015 on $\Omega_m$, but on which unmodelled IA would shift the results by 6.0, 3.4, and $13.2\sigma$, respectively. In light of these, it is clear that IA forward modelling will be critical for these alternative statistics.

When comparing the global impact of an unmodelled δ-NLA IA signal versus the baryons feedback presented in Martinet et al. (2021b, see their table 2), we notice that IA affects the γ-2PCF more severely for all parameters, with $(\Delta S_8, \Delta w_0, \Delta \Omega_m) = (-0.083, -0.878, -0.124)$ for IA compared to $(-0.034, 0.037, -0.035)$ found from analysing the Magneticum simulation with an identical survey setup. Indeed, IA affects all angular scales in $\xi_{\pm}(9)$, with an important redshift modulation, causing a large shift in all parameters. The story is similar for all three aperture mass statistics, where we observe that in the full tomographic case, baryons have a smaller impact on $S_8$ (e.g. $\Delta S_8 = -0.024$ versus $-0.058$ for IA, as measured with peak statistics), on $\Delta w_0$ (0.035 versus $-0.524$), and on $\Delta \Omega_m$ ($-0.005$ versus $-0.113$). As for the γ-2PCFs, this can be explained by the fact that while baryon feedback has a mild redshift dependence, IA

exhibits a more complex structure across the different tomographic bin combinations that is partly degenerate with the growth of structure, directly affected by both $w_0$ and $\Omega_m$ but less so by $S_8$. This picture is also supported by the fact that we find similar biases for baryons and IA in the non-tomographic cases.

We repeated the measurement of these biases when changing $\sigma_G$ from $0.1$ to $0.5 h^{-1}$ Mpc and found that the biases are reduced on average for $\Delta S_8$ and $\Delta w_0$, while $\Delta \Omega_m$ are slightly larger, compared to our fiducial case. The residual noise in the observed marginalized likelihood prevents us from deriving a meaningful scaling relation between the cosmological biases and the smoothing length, and hence future analyses will likely need to include this as a free parameter in the IA model.

## 5.3 Discussion

The methods presented in this paper offer an efficient and accurate avenue to model and mitigate over the uncertain effect and strength of the IA, at the simulation catalogue level. Weak lensing analyses based on statistics beyond two point functions, which are generally relying
on simulation-based inference methods, can now be augmented with a series of dedicated IA-infused mocks from which the secondary signal can be studied for different IA model (NLA, δ-NLA, TATT) as a function of the corresponding parameter values ($A_{\text{IA}}$, $b_{\text{IA}}$, ...).

From this point onwards, the simplest approach would be to cut out the data elements that are causing the bias in cosmology. For example, the likelihood analysis presented in the last section could be repeated for different selections of $S/N$ bin and choice of tomographic setup, until the bias $\Delta$ falls well under the precision $\delta$ for every cosmological parameters we wish to measure. This is similar in spirit with the approach adopted by the DES-Y1 cosmic shear analysis of Troxel et al. (2018), where conservative<fnref n="19" /> small-scale cuts are applied such as to guarantee that the effect of baryons and non-linear physics had a controlled and minor impact on the inferred cosmology. This is not adequate for $\gamma$-2PCFs here since most elements are affected by IA, but could be useful for the aperture mass statistics investigated in this paper, for which the effect is localized on to a small number of $S/N$ bins.

The elements of the data vector mostly contaminated by IA are also highly sensitive to cosmology, since most of the contamination arises from the cross-correlation of the signal of interest and the systematic via the *G*1-like contributions. It is thus more optimal to attempt some form of modelling in the future instead. There exists a few options to achieve this. First, one could measure the shift in cosmology and in error bars caused by different realistic IA models separately from the likelihood sampling (e.g. using the results from Table 3, then correct the inferred cosmology and inflate error bars accordingly. This approach was adopted by Martinet et al. (2018) for the treatment of baryon feedback on peak statistics in the analysis of the KiDS-450 data. The risk here is that there might be interactions between the IA effect and any additional parameters sampled in the full likelihood analysis but excluded from the calculation of Table 3 (e.g. shape calibration, redshift uncertainty, etc.).

To avoid this, it is better to model the IA at the level of the measurement data vector and apply a correction factor to the data, as done in e.g. Kacprzak et al. (2016) and in Harnois-Déraps et al. (2021). Repeating for different values of $A_{\text{IA}}$ would allow us to estimate the impact on the cosmology inference and the associated systematic uncertainty, to be combined with the statistical error. The main drawback from this method is that the correction factor is computed at a fixed cosmology. This effectively assumes a weak cosmological dependence of the IA signal relative to the underlying signal, an assumption that has yet to be tested and that could also lead to important biases.

It would be preferable to infuse IA in a series of mocks that cover the full parameter space volume inside of which the likelihood is estimated. This forward modelling approach is more computationally demanding as it requires the calculation of projected tidal fields at every cosmological training node, however the gain in accuracy is guaranteed since it correctly captures the full cosmological dependence of the secondary signal. This is equivalent to the treatment of IA in the $\gamma$-2PCFs pipeline, for which the model (see equations 4 and 5) are re-evaluated alongside the *GG* signal (equation 3) from the matter power spectrum at each step of the likelihood sampling. This will result in a slight increased uncertainty on the cosmological parameters of interest, as found in Zürcher et al. (2021) with an alternative IA-infusion method.

In fact, a number of results presented in this paper can also be obtained from the techniques introduced in Fluri et al. (2019) and Zürcher et al. (2021), which is also based on projected mass sheets. The key improvements over their method is that we work at the level of galaxy catalogues, as opposed to lensing maps, and are therefore closer to the data produced by the lensing surveys. Each galaxy is assigned a random orientation, a cosmic shear signal, and an intrinsic alignment extracted from a flexible tidal field-based model (only the NLA model is implemented in Zürcher et al. 2021). Therefore, combined with the convergence and the local density fields, our method can construct reduced shear values – the lensing observables – with more complex IA models, for example by including higher order operators in the tidal fields.

One of the key open questions is whether the accuracy of our infusion method is high enough for the Stage-IV lensing surveys. As presented in Blazek et al. (2019), next generation data analyses will require flexibility and accuracy in the modelling of the IA sector, without which the inferred cosmological parameters could incur catastrophic biases. In particular, they have demonstrated that an LSST-like data set infused with the TATT IA model but analysed with the NLA would result in $S_8$ and $w_0$ values that are significantly too low. Adding a redshift dependence to the NLA model reduces these biases down to ~1$\sigma$ shifts, however only a full TATT model produces unbiased results. The methods presented in this paper are designed to enable this flexibility, however there are residual discrepancies with the input model that will need to be improved upon, as shown in Fig. 4. It is not entirely clear at this stage what causes this (we suspect a projection effect), nor how these deviations will impact the inferred cosmology after marginalization. For this a full likelihood analysis will need to be carried out. This calls for further tests to be done in the near future, notably on the convergence of the IA signal as a function of the projection thickness.

We present in this paper our implementation of the NLA and $\delta$-NLA models, however additional physics is likely to be required to analyse upcoming data, in particular the tidal torquing term and a red/blue split dependence. Whereas the former is within reach with the *cosmo*-SLICS simulations, the later could be achieved by applying our methods on external light-cone simulations where these galaxy properties exist, for example in the Millenium,<fnref n="20" /> MICE,<fnref n="21" /> *Euclid* Flagship (Potter, Stadel & Teyssier 2017), or the LSST DC2²² simulations. One could construct different samples based on the colour information (or any other existing galaxy property that could correlate with the IA) and assign an intrinsic alignment per galaxy from projected tidal field, with a coupling that is allowed to vary for different samples. This hybrid scheme would allow for more detailed IA models that hopefully better describe the physics observed in spectroscopic data, but without the need to resolve the halo shapes. This is a particularly powerful advantage, since measuring halo shapes through the inertia matrix can only be achieved for haloes resolved with hundreds of particles, introducing a low-mass cut on the haloes for which galaxies have a reliable IA assigned.

We note that having access to the halo information has an additional advantage, in that one could study separately the impact of IA on peaks that correspond to collapsed regions and those that are due to line-of-sight projection, as in White, van Waerbeke & Mackey (2002). Additionally, our IA-infusion method can be combined with other IA mitigation techniques, including exploiting the correlation

<footnote n="19">¹⁹The $\xi_+$ and $\xi_-$ data vectors in Troxel et al. (2018) are truncated for $9 < [4-6]$ arcmin and $< [40-90]$ arcmin, respectively.</footnote>

<footnote n="20">²⁰wwwmpa.mpa-garching.mpg.de/galform/virgo/millennium/</footnote>

<footnote n="21">²¹maia.ice.cat/mice</footnote>

<footnote n="22">²²portal.nersc.gov/project/lsst/cosmoDC2/.README.html</footnote>
length (Catelan et al. 2001), nulling (Huterer & White 2005; Joachimi & Schneider 2010), self-calibration (Bernstein & Jain 2004; Bernstein 2009; Yao et al. 2020; Sánchez et al. 2021), and using galaxy colour information (Tugendhat, Reischke & Schaefer 2020). So far these techniques have been only applied to the standard $\gamma$-2PCFs and it would be essential to study how they perform in non-Gaussian statistics. In fact, this whole approach could be phrased as a search for the (non-Gaussian) lensing estimator that is the best compromise between capturing a maximum of cosmological information, while being minimally sensitive to a suite of IA models.

Although all calculations presented in this paper are carried out in the flat-sky limit, it is straightforward to generalize to full sky, replacing the Fourier transforms with spherical harmonic transforms.

Finally, this approach could easily be extended to simulations with massive neutrinos, modified gravity or a hydrodynamical sector, such as to investigate possible degeneracies between the impact of distinct physical phenomena (i.e. IA, baryons, neutrinos, and gravity) on higher order lensing statistics.

# 6 CONCLUSIONS

Intrinsic alignments of galaxies cause a secondary signal in cosmic shear measurements that, if unaccounted for, will create catastrophic biases (>$10\sigma$) in the cosmology inferred from the next generation of lensing surveys. The physical mechanism of this effect is still under debate, however many models proposed in the literature seem to agree relatively well with current data (Blazek et al. 2019; Johnston et al. 2019; Fortuna et al. 2021). While these are primarily aimed for two-point statistics, we show in this paper that different IA models can also be infused at the catalogue level, which can subsequently serve to model the IA on statistics beyond two-points.

We construct 2D tidal field maps with the projected mass sheets extracted from the cosmo-SLICS wCDM simulations, and couple these maps with lensing galaxy catalogues in order to assign an intrinsic alignment to individual objects. Our current model includes three free-parameters, namely the smoothing scale $\sigma_G$, the tidal coupling parameter $A_{\text{IA}}$, and the density-weighting parameter $b_{\text{TA}}$. We validate our infusion model by comparing two-point statistics against predictions from the NLA model (including the GG, GI, II terms) and recover good match for most tomographic bins. The II term is less well modelled, likely due to its increased sensitivity to projection effects.

We measure the impact of the $\delta$-NLA IA model on three different aperture mass map statistics – peaks, minima, and lensing PDF – and find that region of large (positive or negative) $S/N$ are more affected, especially when combining galaxies with different tomographic bins. Of these three probes, the lensing minima are more immune, primarily due to their lesser sensitivity to regions of high density (and therefore of high IA). When propagated on to a full cosmological inference, we find that the $\gamma$-2PCFs suffer the largest biases in presence of unmodelled IA. Table 3 summarizes the cosmological biases incurred in presence of unmodelled IA, for various tomographic choices. Most notably, and for all probes, the inferred $w_0$ is shifted to more negative values, while $S_8$ and $\Omega_m$ are lower compared to the input truth.

Some sort of IA modelling will be required for all probes, however the effect is more localized than for $\gamma$-2PCFs, for which all angular scales are affected. In general, the study presented here offers an avenue to identify parts of the non-Gaussian data vector that are the least affected by IA, given a specified IA model. We intend to proceed with a dedicated investigation of this in an upcoming work,

notably for Minkowski functionals, voids, and clipped shear. The method we develop here can be easily extended to marginalize over $A_{\text{IA}}$ and $b_{\text{TA}}$ for any weak lensing estimator, simply by multiplying the '$\epsilon_{1/2}^\text{IA}$' catalogue entries by a real number and combining with the cosmic shear $\gamma_{1/2}$ columns. This opens up the possibility to jointly sample over the uncertain IA parameters in a joint [$\gamma$-2PCFs; $N_{\text{peaks}}(S/\mathcal{N})$] weak lensing analysis, improving constraints on both the cosmological and the IA parameters.

# ACKNOWLEDGEMENTS

We would like to thank Simon Samuroff, for sharing his density-to-tidal field conversion code, Jonathan Blazek for his insights on the TATT model and on the interpretation of our measurements, and Ludovic van Waerbeke and François Lanusse for useful discussions. We are thankful to the referee for pointing out an error in our equations (equations 4-5-6), which after correction resulted in an improved agreement with the theoretical NLA model. JHD acknowledges support from an Ernest Rutherford Fellowship (Science and Technology Facilities Council, project reference ST/S004858/1). Computations were carried on the Cuillin cluster at the Royal Observatory of Edinburgh, which is managed by Eric Tittley, and on the Mardec cluster supported by the OCEVU Labex (ANR-11-LABX-0060) and the Excellence Initiative of Aix-Marseille University - A*MIDEX, part of the French ‘Investissements d’Avenir’ program. RR is supported by the European Research Council (Grant No. 770935).

All authors contributed to the development and writing of this paper. JHD led the analysis; NM conducted the cosmology inference, while RR co-developed the IA infusion formalism.

# DATA AVAILABILITY

The SLICS numerical simulations can be found at http://slics.roe.ac.uk/, while the cosmo-SLICS can be made available upon request.

# REFERENCES

* Amon A. et al., 2021, preprint (arXiv:2105.13543)

* Asgari M. et al., 2021, A&A, 645, A104

* Bernstein G. M., 2009, ApJ, 695, 652

* Bernstein G., Jain B., 2004, ApJ, 600, 17

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

* Blazek J. A., MacCrann N., Troxel M. A., Fang X., 2019, Phys. Rev. D, 100, 103506

* Boyle A., Uhlemann C., Friedrich O., Barthelemy A., Codis S., Bernardeau F., Giocoli C., Baldi M., 2021, MNRAS, 505, 2886

* Bridle S., King L., 2007, New J. Phys., 9, 444

* Brouwer M. M. et al., 2018, MNRAS, 481, 5189

* Brown M. L., Taylor A. N., Hambly N. C., Dye S., 2002, MNRAS, 333, 501

* Castro T., Quartin M., Giocoli C., Borgani S., Dolag K., 2018, MNRAS, 478, 1305

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

* Cheng S., Ting Y.-S., Ménard B., Bruna J., 2020, MNRAS, 499, 5902

* Chiang L.-Y., Coles P., 2000, MNRAS, 311, 809

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

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

* Coulton W. R., Liu J., McCarthy I. G., Osato K., 2020, MNRAS, 495, 2531

* Davies C. T., Cautun M., Giblin B., Li B., Harnois-Déraps J., Cai Y.-C., 2021, MNRAS, 507, 2267

* Dietrich J. P., Hartlap J., 2010, MNRAS, 402, 1049

* Euclid Collaboration, 2019, MNRAS, 484, 5509

* Fluri J., Kacprzak T., Lucchi A., Refregier A., Amara A., Hofmann T., Schneider A., 2019, Phys. Rev. D, 100, 063514
Fortuna M. C., Hoekstra H., Joachimi B., Johnston H., Chisari N. E., Georgiou C., Mahony C., 2021, MNRAS, 501, 2983

Gatti M. et al., 2020, MNRAS, 498, 4060

Giblin B. et al., 2018, MNRAS, 480, 5529

Giblin B. et al., 2021, A&A, 645, A105

Hamana T. et al., 2020, PASJ, 72, 16

Harnois-Déraps J. et al., 2016, MNRAS, 460, 434

Harnois-Déraps J. et al., 2017, MNRAS, 471, 1619

Harnois-Déraps J. et al., 2018, MNRAS, 481, 1337

Harnois-Déraps J., Giblin B., Joachimi B., 2019, A&A, 631, A160

Harnois-Déraps J., Martinet N., Castro T., Dolag K., Giblin B., Heymans C., Hildebrandt H., Xia Q., 2021, MNRAS, 506, 1623

Heydenreich S., Brück B., Harnois-Déraps J., 2021, A&A, 648, A74

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

Hikage C. et al., 2019, PASJ, 71, 43

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

Hildebrandt H. et al., 2017, MNRAS, 465, 1454

Hildebrandt H. et al., 2021, A&A, 647, A124

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

Huterer D., White M., 2005, Phys. Rev. D, 72, 043002

Jarvis M., Bernstein G., Jain B., 2004, MNRAS, 352, 338

Joachimi B., Schneider P., 2010, A&A, 517, A4

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

Johnston H. et al., 2019, A&A, 624, A30

Joudaki S. et al., 2018, MNRAS, 474, 4894

Kacprzak T. et al., 2016, MNRAS, 463, 3653

Kaiser N., Squires G., 1993, ApJ, 404, 441

Kilbinger M., Bonnett C., Coupon J., 2014, Astrophysics Source Code Library, record ascl:1402.026

Li Z., Liu J., Zorrilla Matilla J. M., Coulton W. R., 2019, Phys. Rev. D, 99, 063527

Liu J., Madhavacheril M. S., 2019, Phys. Rev. D, 99, 083508

Liu J., Petri A., Haiman Z., Hui L., Kratochvil J. M., May M., 2015a, Phys. Rev. D, 91, 063507

Liu X. et al., 2015b, MNRAS, 450, 2888

Liu J., Bird S., Zorrilla Matilla J. M., Hill J. C., Haiman Z., Madhavacheril M. S., Petri A., Spergel D. N., 2018, J. Cosmol. Astropart. Phys., 2018, 049

Lu T., Haiman Z., 2021, MNRAS, 506, 3406

Mandelbaum R., 2018, ARA&A, 56, 393

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

Martinet N. et al., 2018, MNRAS, 474, 712

Martinet N., Castro T., Harnois-Déraps J., Jullo E., Giocoli C., Dolag K., 2021a, A&A, 648, A115

Martinet N., Harnois-Déraps J., Jullo E., Schneider P., 2021b, A&A, 646, A62

Mead A., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401

Nishimichi T. et al., 2019, ApJ, 884, 29

Osato K., Shirasaki M., Yoshida N., 2015, ApJ, 806, 186

Petri A., Liu J., Haiman Z., May M., Hui L., Kratochvil J. M., 2015, Phys. Rev. D, 91, 103511

Potter D., Stadel J., Teyssier R., 2017, Comput. Astrophys. Cosmol., 4, 2

Samuroff S. et al., 2019, MNRAS, 489, 5453

Samuroff S., Mandelbaum R., Blazek J., 2021, MNRAS, 508, 637

Sánchez C. et al., 2021, preprint (arXiv:2105.13542)

Schirmer M., Erben T., Hetterscheidt M., Schneider P., 2007, A&A, 462, 875

Schneider P., 1996, MNRAS, 283, 837

Schneider M. D., Bridle S., 2010, MNRAS, 402, 2127

Schneider A., Teyssier R., Stadel J., Chisari N. E., Le Brun A. M. C., Amara A., Refregier A., 2019, J. Cosmol. Astropart. Phys., 2019, 020

Secco L. F. et al., 2021, preprint (arXiv:2105.13544)

Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132

Shan H. et al., 2018, MNRAS, 474, 1116

Sifón C., Hoekstra H., Cacciato M., Viola M., Köhlinger F., van der Burg R. F. J., Sand D. J., Graham M. L., 2015, A&A, 575, A48

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

Smith R. E. et al., 2003, MNRAS, 341, 1311

Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152

Troxel M. A. et al., 2018, Phys. Rev. D, 98, 043528

Tugendhat T. M., Reischke R., Schaefer B. M., 2020, 2020, 494, 2969

van Uitert E. et al., 2018, MNRAS, 476, 4662

van Waerbeke L., et al., 2013, MNRAS, 433, 3373

Vlah Z., Chisari N. E., Schmidt F., 2020, J. Cosmol. Astropart. Phys., 2020, 025

Weiss A. J., Schneider A., Sgier R., Kacprzak T., Amara A., Refregier A., 2019, J. Cosmol. Astropart. Phys., 2019, 011

White M., van Waerbeke L., Mackey J., 2002, ApJ, 575, 640

Yao J., Pedersen E. M., Ishak M., Zhang P., Agashe A., Xu H., Shan H., 2020, MNRAS, 495, 3900

Zjupa J., Schäfer B. M., Hahn O., 2020, preprint (arXiv:2010.07951)

Zuntz J. et al., 2015, Astron. Comput., 12, 45

Zürcher D., Fluri J., Sgier R., Kacprzak T., Refregier A., 2021, J. Cosmol. Astropart. Phys., 2021, 028

## APPENDIX A: PROJECTED TIDAL FIELD

We derive in this Appendix the equation for the projected tidal field that is used in the main text (equation 10). Starting from equation (9), the 3D field projected along the z-axis can be computed as:

$$ 
\begin{aligned}
s_{2D}(\mathbf{x}_{\perp}) &= \sum_{z} s_{ij}(\mathbf{x}) \\
&= \sum_{z} \int d^3\mathbf{k} \exp[i\mathbf{k}\mathbf{x}] \left[ \frac{k_i k_j}{k^2} - \frac{1}{3} \right] \tilde{\delta}(\mathbf{k}) \mathcal{G}(\sigma_G) \\
&= \sum_{z} \int d^3\mathbf{k} \exp[i\mathbf{k}x_x + i\mathbf{y}k_y + izk_z] \\
&\qquad \times \left[ \frac{k_i k_j}{k^2} - \frac{1}{3} \right] \tilde{\delta}(\mathbf{k}) \mathcal{G}(\sigma_G).
\end{aligned}
\quad (A1) $$

The sum over $z$ affects only the third term of the exponential function, hence using the identity $\int_{-a}^{a} \exp[i\mathbf{k}z]dz = 2a(\sin k_z a)(k_z a) = 2\pi \delta_D(k_z)$ in the large $a$ limit, we can write:

$$ s_{2D}(\mathbf{x}_{\perp}) = 2\pi \int d^3\mathbf{k}_{\perp} \exp[i\mathbf{k}_{\perp}\mathbf{x}_{\perp} + i\mathbf{y}\mathbf{k}_{\perp}] \times \left[ \frac{k_i k_j}{k^2} - \frac{1}{3} \right] \tilde{\delta}(\mathbf{k}_{\perp}, k_z = 0) \mathcal{G}_{2D}(\sigma_G). \quad (A2) $$

Note that the Dirac delta function, $\delta_D(k_z)$, ensures that $k^2$ appearing in the denominator can be computed with $k_x$ and $k_y$ only, receiving no contribution from $k_z$. The quantity $\tilde{\delta}(\mathbf{k}_{\perp}, k_z = 0)$ can therefore be identified as the $x-y$ Fourier transform of the z-projected 3D density field, $\delta_{2D}(\mathbf{x}_{\perp})$. As a result, the projected tidal field, at position $\mathbf{x}_{\perp}$, can be evaluated in Fourier space by solving:

$$ \tilde{s}_{2D}(\mathbf{k}_{\perp}) = \sum_{z} \tilde{s}_{ij}(\mathbf{k}) = 2\pi \left[ \frac{k_i k_j}{k^2} - \frac{1}{3} \right] \tilde{\delta}_{2D}(\mathbf{k}) \mathcal{G}_{2D}(\sigma_G), \quad (A3) $$

for $ij = (xx, yy, xy)$, then inverse Fourier transformed to produce the three tidal maps, $s_{xx}$, $s_{yy}$, and $s_{xy}$.

The prescription to assign an intrinsic alignment based on the projected tidal field, described in equation (6), involves the combinations $(s_{xx} - s_{yy})$ and $s_{xy}$, which therefore correspond to:
$$ \begin{aligned} \tilde{\epsilon}_1^{1A}(\mathbf{k}_{\perp}) &\propto \left(\frac{k_x^2 - k_y^2}{k^2}\right) \tilde{\delta}_{2D}(\mathbf{k}_{\perp}) \mathcal{G}_{2D}(\sigma_G), \\ \tilde{\epsilon}_2^{1A}(\mathbf{k}_{\perp}) &\propto \left(\frac{k_x k_y}{k^2}\right) \tilde{\delta}_{2D}(\mathbf{k}_{\perp}) \mathcal{G}_{2D}(\sigma_G) \end{aligned} \quad (A4) $$

Aside from the smoothing kernel, these are the same filters that are used for converting convergence maps to shear maps under the Kaiser & Squires (1993, KS hereafter) inversion:

$$ \tilde{\gamma}_1(\ell) = \left( \frac{k_x^2 - k_y^2}{k^2} \right) \tilde{\kappa}(\ell), \qquad \tilde{\gamma}_2(\ell) = \left( \frac{k_x k_y}{k^2} \right) \tilde{\kappa}(\ell) \quad (A5) $$

meaning that one can linearly combine the mass sheets with the correct coefficients and obtain intrinsic ellipticities from a normal KS inversion. This fact is exploited in Fluri et al. (2019) to generate pure 1A convergence maps that are then combined with the cosmic shear convergence maps to create a contamination.

## APPENDIX B: STAGE-IV BIASES FOR DIFFERENT $A_{IA}$ AND $\sigma_G$ VALUES

The main results in Section 5 are reported cosmological biases for a coupling parameter set to $A_{IA} = 1.5$. In this Appendix, we extend these results to other analysis choices. We first show the marginalized constraints from the full tomographic peak statistics for $A_{IA} = 1.0$ and 2.0 in Fig. B1, the comparison between two smoothing scales $\sigma_G$ (0.1 and $0.5 \ h^{-1} \ \text{Mpc}$) in Fig. B2, and the results from the three aperture mass map statistics in Fig. B3. As expected, a stronger coupling with the tidal fields and a smaller smoothing length both enhance the cosmological biases. The marginalized posterior are tabulated in Table B1 for different estimators and two coupling strengths, again assuming the $\delta$-NLA model. The choice of smoothing scale has an important effect on all three parameters, and needs to be physically motivated, and possibly included as a free parameter in the model.

<figure><img src="image_12.png" /><figcaption>Figure B1. Marginalized constraints from the full tomographic peak count statistics, with $A_{IA} = 1.0$ (left-hand panel) and 2.0 (right-hand panel).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure B2. Marginalized constraints from the full tomographic peak count statistics, with $\sigma_G = 0.1 \ h^{-1} \ \text{Mpc}$ (left-hand panel) and $\sigma_G = 0.5 \ h^{-1} \ \text{Mpc}$ (right-hand panel).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure B3.** Marginalized constraints from the full tomographic minima count statistics (left-hand panel) and lensing PDF (right-hand panel), with $A_{\text{IA}} = 1.5$.</figcaption></figure>
<table><caption>**Table B1.** Bias on the inferred cosmology caused by the intrinsic alignments of galaxies, if left unaccounted for, for the different cases listed in Table 2. Parameter shifts are reported for two different values of $A_{1A}$. A smoothing scale of $0.1 h^{-1}$ Mpc is assumed here. Numbers in parenthesis show the percentage of the marginalized error with respect to the input parameter value.</caption><thead><tr><th rowspan="2">case</th><th rowspan="2">estimator(s)</th><th rowspan="2">ΔS<sub>8</sub></th><th colspan="2">A<sub>1A</sub> = 1.0</th><th colspan="2">A<sub>1A</sub> = 2.0</th></tr><tr><th>Δw<sub>0</sub></th><th>ΔΩ<sub>m</sub></th><th>Δw<sub>0</sub></th><th>ΔΩ<sub>m</sub></th></tr></thead><tbody><tr><td>(i)</td><td>N<sub>peaks</sub>, no tomo</td><td>-0.037 (-4.6 per cent)</td><td>0.038 (3.8 per cent)</td><td>-0.007 (-2.2 per cent)</td><td>-0.066 (-8.1 per cent)</td><td>0.038 (3.8 per cent)</td><td>-0.007 (-2.2 per cent)</td></tr><tr><td>(ii)</td><td>N<sub>peaks</sub>, tomo auto</td><td>-0.003 (-0.3 per cent)</td><td>-0.009 (-0.9 per cent)</td><td>-0.01 (-3.4 per cent)</td><td>-0.006 (-0.7 PER CENT)</td><td>-0.009 (-0.9 per cent)</td><td>-0.015 (-5.2 per cent)</td></tr><tr><td>(iii)</td><td>N<sub>peaks</sub>, tomo auto + pairs</td><td>-0.011 (-1.3 per cent)</td><td>0.009 (0.9 per cent)</td><td>-0.028 (-9.6 per cent)</td><td>-0.049 (-6.0 per cent)</td><td>-0.6 (-60.0 per cent)</td><td>-0.113 (-39.0 per cent)</td></tr><tr><td>(iv)</td><td>N<sub>peaks</sub>, tomo full</td><td>-0.049 (-5.9 per cent)</td><td>-0.545 (-54.5 per cent)</td><td>-0.095 (-32.6 per cent)</td><td>-0.06 (-7.3 per cent)</td><td>-0.512 (-51.2 per cent)</td><td>-0.125 (-43.0 per cent)</td></tr><tr><td>(i)</td><td>N<sub>min</sub>, no tomo</td><td>-0.037 (-4.4 per cent)</td><td>0.0 (0.0 per cent)</td><td>0.0 (0.0 per cent)</td><td>-0.073 (-8.9 per cent)</td><td>0.0 (0.0 per cent)</td><td>0.0 (0.0 per cent)</td></tr><tr><td>(ii)</td><td>N<sub>min</sub>, tomo auto</td><td>0.0 (0.0 per cent)</td><td>0.0 (0.0 per cent)</td><td>0.0 (0.0 per cent)</td><td>0.0 (0.0 per cent)</td><td>0.0 (0.0 per cent)</td><td>0.012 (4.0 per cent)</td></tr><tr><td>(iii)</td><td>N<sub>min</sub>, tomo auto + pairs</td><td>-0.02 (-2.4 per cent)</td><td>-0.285 (-28.5 per cent)</td><td>-0.064 (-22.1 per cent)</td><td>-0.051 (-6.2 per cent)</td><td>-0.57 (-57.0 per cent)</td><td>-0.103 (-35.4 per cent)</td></tr><tr><td>(iv)</td><td>N<sub>min</sub>, tomo full</td><td>-0.051 (-6.2 per cent)</td><td>-0.582 (-58.2 per cent)</td><td>-0.084 (-29.0 per cent)</td><td>-0.064 (-7.7 per cent)</td><td>-0.564 (-56.4 per cent)</td><td>-0.116 (-39.8 per cent)</td></tr><tr><td>(i)</td><td>N<sub>pix</sub>, no tomo</td><td>-0.046 (-5.6 per cent)</td><td>0.037 (3.7 per cent)</td><td>0.0 (0.0 per cent)</td><td>-0.078 (-9.5 per cent)</td><td>0.027 (2.7 per cent)</td><td>0.0 (0.0 per cent)</td></tr><tr><td>(ii)</td><td>N<sub>pix</sub>, tomo auto</td><td>-0.0 (-0.0 per cent)</td><td>0.036 (3.6 per cent)</td><td>-0.015 (-5.3 per cent)</td><td>0.003 (-0.3 per cent)</td><td>0.036 (3.6 per cent)</td><td>-0.004 (-1.4 per cent)</td></tr><tr><td>(iii)</td><td>N<sub>pix</sub>, tomo autopairs</td><td>0.018 (2.2 per cent)</td><td>0.039 (3.9 per cent)</td><td>-0.023 (-7.9 per cent)</td><td>0.014 (1.7 per cent)</td><td>0.039 (3.9 per cent)</td><td>-0.057 (-19.8 per cent)</td></tr><tr><td>(iv)</td><td>N<sub>pix</sub>, tomo full</td><td>-0.006 (-0.8 per cent)</td><td>0.059 (5.9 per cent)</td><td>-0.066 (-22.7 per cent)</td><td>-0.025 (-3.0 per cent)</td><td>0.04 (4.0 per cent)</td><td>-0.109 (-37.7 per cent)</td></tr></tbody></table>

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