# Three-point intrinsic alignments of dark matter halos in the IllustrisTNG
simulation

Susan Pyne,1 Ananth Tenneti1 and Benjamin Joachimi1
1Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
E-mail<EMAIL_ADDRESS>
###### Abstract

We use the IllustrisTNG suite of cosmological simulations to measure intrinsic
alignment (IA) bispectra of dark matter subhalos between redshifts 0 and 1. We
decompose the intrinsic shear field into E- and B-modes and find that the
bispectra $B_{\delta\delta\mathrm{E}}$ and $B_{\delta\mathrm{EE}}$, between
the matter overdensity field, $\delta$, and the E-mode field, are detected
with high significance. We also model the IA bispectra analytically using a
method consistent with the two-point non-linear alignment model. We use this
model and the simulation measurements to infer the intrinsic alignment
amplitude $A_{\mathrm{IA}}$ and find that values of $A_{\mathrm{IA}}$ obtained
from IA power spectra and bispectra agree well at scales up to
$k_{\mathrm{max}}=2\,h\mathrm{Mpc}^{-1}$, for example at $z=1$
$A_{\mathrm{IA}}=2.13\pm$ 0.02 from the cross power spectrum between the
matter overdensity and E-mode fields and $A_{\mathrm{IA}}=2.11\pm$ 0.03 from
$B_{\delta\delta\mathrm{E}}$. This demonstrates that a single physically
motivated model can jointly model two-point and three-point statistics of
intrinsic alignments, thus enabling a cleaner separation between intrinsic
alignments and cosmological weak lensing signals.

###### keywords:

cosmology: theory – large-scale structure of Universe – gravitational lensing:
weak – methods: numerical

††pagerange: Three-point intrinsic alignments of dark matter halos in the
IllustrisTNG simulation–17

##  1 Introduction

Gravitational tidal effects cause dark matter halos and the galaxies within
them to preferentially align with the cosmological large-scale structure. Thus
the shapes and orientations of galaxies and halos may be correlated across
cosmological distances. This phenomenon is known as intrinsic alignment (IA)
(Kirk et al., 2012; Joachimi et al., 2015; Kiessling et al., 2015; Kirk et
al., 2015; Troxel & Ishak, 2015).

Intrinsic alignments have been widely studied, mainly because the IA effects
mimic cosmic shear and thus are an undesirable contaminant of weak lensing
shape measurements. Controlling this systematic uncertainty has been a major
concern for recent weak lensing surveys, such as the Kilo-Degree
Survey111http://kids.strw.leidenuniv.nl/index.php (Joachimi et al., 2021) and
the Dark Energy Survey222https://www.darkenergysurvey.org (Secco et al.,
2022b). For forthcoming surveys such as Euclid333http://sci.esa.int/euclid/
(Laureijs et al., 2011) and the Vera C. Rubin Observatory Legacy Survey of
Space and Time444https://www.lsst.org (LSST), control of systematics will
become even more pressing. These surveys will measure the shapes of billions
of galaxies, reducing statistical uncertainties, so that uncontrolled
systematics will be a limiting issue. Looking further ahead, as understanding
of intrinsic alignments develops it is also possible to see this signal as a
valuable cosmological probe in its own right which future surveys will be able
to exploit (Chisari & Dvorkin, 2013; Taruya & Okumura, 2020). All these
considerations suggest that it is worthwhile to explore intrinsic alignments
as broadly as possible.

So far most theoretical and observational studies of intrinsic alignments have
used two-point statistics; relatively few have taken the further step of
considering three-point measurements. However there is evidence that intrinsic
alignments affect two-point and three-point weak lensing statistics
differently (Huterer et al., 2006; Semboloni et al., 2008; Troxel & Ishak,
2011, 2012; Pyne & Joachimi, 2021). For this reason most work which has
investigated three-point statistics aimed to use them to self-calibrate
intrinsic alignment contamination in weak lensing data.

Three-point IA statistics have already been successfully measured from both
survey data and simulations. On the observational side Semboloni et al. (2010)
measured three-point aperture mass statistics from the Cosmic Evolution Survey
(Scoville et al., 2007) and Fu et al. (2014) built on this work using the
larger Canada France Hawaii Telescope Lensing Survey (CFHTLenS; Heymans et al.
2012). Both these studies showed that using three-point statistics could help
improve constraints on cosmological parameters. Early simulation results were
described by Semboloni et al. (2008) who measured three-point intrinsic
alignment aperture mass statistics from simulations described in Heymans et
al. (2006).

Another strand of work has involved analytical modelling of three-point IA
statistics, building on methods developed for two-point statistics, in
particular the non-linear alignment model (NLA) (Hirata & Seljak, 2004; Bridle
& King, 2007). This postulates that the intrinsic alignment of halos is
related to the tidal gravitational field at an earlier redshift. This
methodology can been extended to three-point statistics in a natural but
empirical way (Troxel & Ishak, 2012; Merkel & Schäfer, 2014; Pyne & Joachimi,
2021). Other intrinsic alignment models, for example the halo-model of
Schneider & Bridle (2010) and tidal alignment models such as Blazek et al.
(2015), have generally only been implemented for two-point statistics. An
exception is the development by Vlah et al. (2020, 2021) of effective field
theory models of galaxy alignments. These authors modelled two-point
statistics at next-to-leading order and three-point statistics at leading
order.

Despite the body of work described above, an analytical model for three-point
intrinsic alignment statistics has never been tested against simulations. In
this work we aim to fill this gap. As a starting point we consider dark matter
subhalos rather than galaxies. Since the intrinsic alignment of galaxies is
related to that of their host halos, our work can be seen as a first step
towards a model for galaxies.

Our approach to measurement from simulations is based on that of Kurita et al.
(2021) who developed a novel method to measure three-dimensional intrinsic
alignment power spectra of dark matter halos which they applied to the high-
resolution DarkQuest suite of simulations (Nishimichi et al., 2019; Miyatake
et al., 2021). This has a resolution of 20483 particles in a periodic cubic
box size of 1 $h^{-1}$ Gpc. Using the same methodology, we measure all
theoretically non-zero intrinsic alignment bispectra for equilateral and
specific isosceles triangles from the IllustrisTNG555https://www.tng-
project.org simulation suite.

Kurita et al. (2021) compared their simulation results with the two-point non-
linear alignment model and concluded that at linear scales the model matches
the simulation results, with a scale-independent intrinsic alignment
amplitude, but that at non-linear scales the model breaks down and the
amplitude is no longer constant. We similarly compare our bispectrum
measurements with an analytical model based on perturbation theory which is in
keeping with the two-point NLA model.

In Sect. 2 we describe the simulation suite and the methods used to measure
dark matter subhalo shapes. Sect. 3 explains how we measure IA power spectra
and bispectra from simulations, and Sect. 4 discusses the NLA model and our
analytical model for three-point IA statistics. In Sect. 5 we present our IA
power spectrum and bispectrum measurements from simulations, and compare these
with analytical results. We summarise and discuss possible further work in
Sect. 6. Appendices discuss other details, including signal-to-noise ratios
and a possible phonomenological modification to the non-linear alignment
model.

##  2 Simulations

###  2.1 Characteristics of the simulations

We use the publicly available cosmological simulation suite IllustrisTNG
(Marinacci et al., 2018; Naiman et al., 2018; Nelson et al., 2018, 2019;
Pillepich et al., 2018; Springel et al., 2018). Each simulation in the suite
uses the moving-mesh code Arepo (Springel, 2010)666https://arepo-code.org and
self-consistently solves for the coupled evolution of dark matter, cosmic gas,
luminous stars and supermassive black holes from the starting redshift $z=127$
to the present day, based on a cosmology consistent with Planck (Ade et al.,
2016):
$\Omega_{\Lambda}=0.6911,\Omega_{\mathrm{m}}=0.3089,\Omega_{\mathrm{b}}=0.0486,\sigma_{8}=0.8159,n_{\mathrm{s}}=0.9667,h=0.6774$.
Specifically we use the IllustrisTNG-300-1 hydrodynamic simulation which has
the largest simulation box size within the suite of simulations available from
the public data release (Nelson et al., 2019). The simulation box volume is
300 Mpc3 (comoving) and evolved with $2\,500^{3}$ dark matter particles and an
equal initial number of gas cells.

The publicly available data includes catalogues of halos identified using the
Friends-of-Friends algorithm (Davis et al., 1985). The subfind algorithm
(Springel et al., 2001) is used to identify substructures within the halos.
This algorithm defines the centre of the subhalo as the location of the most
bound particle, and provides the positions of all dark matter particles in the
halo.

For our analysis we select dark matter subhalos from the simulations at three
redshifts, $z=0.0,0.5,1.0$. The subhalo mass is the total sum of the masses of
all individual particles belonging to a given subhalo as identified by the
subfind algorithm. The maximum subhalo masses in the simulation are $\sim$
$3.4\times 10^{14}$, $7.4\times 10^{14}$ and $1.2\times 10^{15}$
$\mathrm{M}_{\odot}$ at $z=1.0,0.5,0$ respectively. We select subhalos within
the mass range $4\times 10^{10}$ to $10^{14}\;h^{-1}\mathrm{M}_{\odot}$ for
our analysis.

###  2.2 Measurement of 3D halo shapes from simulations

We start by measuring 3D ellipsoidal subhalo shapes. For this we adopt the
widely used method based on the inertia tensor. This is defined as

$\displaystyle I_{ij}=\frac{\sum_{n}m_{n}x_{ni}x_{nj}}{\sum_{n}{m_{n}}}\ ,$
(1)

where $m_{n}$ is the mass of the $n$th particle and $x_{ni},x_{nj}$ are its
position coordinates relative to the centre of the halo. The semi-axes $a,b,c$
of the ellipsoid are obtained from the eigenvalues
$\lambda_{a},\lambda_{b},\lambda_{c}$ of the inertia tensor, with
$a=\sqrt{\lambda_{a}}$ and so on. We set $c\leq b\leq a$ and define axis
ratios $s=c/a$ and $q=b/a$. The eigenvectors of the inertia tensor determine
the orientation of the axes.

In order to obtain well resolved shapes, we choose only subhalos with a
minimum of 1000 dark matter particles, which is consistent with the particle
number threshold adopted in earlier studies, for example Tenneti et al.
(2015). With this threshold the minimum subhalo mass is $\sim$ $4\times
10^{10}h^{-1}\mathrm{M}_{\odot}$ .

As an improvement on Eq. (1) we use the reduced inertia tensor (Tenneti et
al., 2015) which gives more weight to particles which are closer to the centre
of a subhalo, thus avoiding potential problems with defining its outer edge.
The reduced inertia tensor is defined as

$\displaystyle\tilde{I}_{ij}=\frac{\sum_{n}m_{n}\frac{x_{ni}x_{nj}}{r_{n}^{2}}}{\sum_{n}{m_{n}}}\
,$ (2)

where

$\displaystyle r_{n}^{2}=\sum_{i=1}^{3}x_{ni}^{2}\ .$ (3)

Rather than taking the ‘raw’ axis ratios defined above we use the iterative
approach described in Tenneti et al. (2015) to recover the shape of an
isodensity surface (Schneider et al., 2012). In this method the principal axes
of the ellipsoids are rescaled iteratively while the enclosed volume is kept
constant. After each rescaling, particles outside the ellipsoidal volume are
discarded. The process is repeated until the fractional change in the axis
ratios is below a predefined limit, in our case 1 per cent.

##  3 Intrinsic alignment spectra from simulations

###  3.1 Ellipticity and tidal shear

The ellipticity, $\epsilon$, of a subhalo shape is a spin-2 quantity (it is
invariant under rotations of integer multiples of $\pi$) which can be
parameterized in several ways in terms of the shape and orientation of the
subhalo. All parameterizations are essentially equivalent so we can make a
choice which suits the problem at hand. In many cases it is convenient to
express the ellipticity as $\epsilon=\epsilon_{+}+\mathrm{i}\epsilon_{\times}$
where $\epsilon_{+}$ represents stretching along a defined axis and
$\epsilon_{\times}$ represent stretching along an axis at 45° to this. It is
always possible to find such a decomposition (Stebbins, 1996).

We can then express the two components of the ellipticity as

$\displaystyle\epsilon_{+}$ $\displaystyle=\bigg{(}\frac{a-b}{a+b}\bigg{)}\cos
2\theta$ (4) $\displaystyle\epsilon_{\times}$
$\displaystyle=\bigg{(}\frac{a-b}{a+b}\bigg{)}\sin 2\theta\ ,$ (5)

where $a$ and $b$ are the semi-major and semi-minor axes of the ellipse. Here
we take the line of sight to be the $z$-axis and consider ellipses projected
on to the plane perpendicular to this (Kurita et al., 2021; Shi et al., 2021).
The angle $\theta$ is an arbitrary choice which we take to be the position
angle of the major axis with respect to the $x$-axis in the $x$-$y$ plane, as
determined by the relevant eigenvector. Eqs.(4) and (5) clearly satisfy the
requirements for a spin-2 quantity. Moreover, under a parity transformation
$\epsilon_{+}$ is unchanged but
$\epsilon_{\times}\rightarrow-\epsilon_{\times}$ (Schneider et al., 2002).

This definition of $\epsilon$, in terms of $(a-b)/(a+b)$, is commonly used for
weak lensing shear because it provides an unbiased estimator of the shear and
does not depend on the ellipticity distribution of source galaxies (Seitz &
Schneider, 1997; Viola et al., 2014). Even though we are considering intrinsic
ellipticity rather than cosmic shear we adopt this definition for consistency.
It has the added advantage that the tidal shear
$\gamma=\gamma_{+}+\mathrm{i}\gamma_{\times}$ can be assumed to be directly
proportional to the ellipticity.

Other studies of intrinsic alignments, for example Blazek et al. (2015) and
Shi et al. (2021), used an alternative definition of ellipticity which
replaces $(a-b)/(a+b)$ with $(a^{2}-b^{2})/(a^{2}+b^{2})$ in Eqs. (4) and (5).
This version is more readily comparable with observations but does not provide
an unbiased shear estimator (Schneider & Seitz, 1995). Moreover rather than
the shear being directly related to the ellipticity, an extra responsivity
factor $\mathcal{R}=1-\epsilon_{\mathrm{rms}}^{2}$ is required (Bernstein &
Jarvis, 2002), where
$\epsilon_{\mathrm{rms}}=\sqrt{\langle\epsilon_{+}^{2}\rangle}=\sqrt{\langle\epsilon_{\times}^{2}\rangle}$.
The responsivity measures the average response of $\epsilon_{+,\times}$ to
$\gamma_{+,\times}$ so that
$\gamma_{+,\times}=\epsilon_{+,\times}/2\mathcal{R}$. With our definition of
the ellipticity we do not need to consider the responsivity and can assume
that $\epsilon$ directly traces the tidal shear field $\bm{\gamma}$. Thus
$\gamma_{+}$ and $\gamma_{\times}$ are given by Eqs. (4) and (5).

###  3.2 Decomposition into E- and B-modes

For studies of intrinsic alignments it is convenient to go a step further and
decompose the shear field into a curl-free (E-mode) component
$\gamma_{\mathrm{E}}$ and a gradient-free (B-mode) component
$\gamma_{\mathrm{B}}$ (Kamionkowski et al., 1998; Crittenden et al., 2002). We
define these by the equations

$\displaystyle\nabla^{2}\gamma_{\mathrm{E}}(\bm{x})$
$\displaystyle=(\partial_{x}\partial_{x}-\partial_{y}\partial_{y})\gamma_{+}(\bm{x})+2\partial_{x}\partial_{y}\gamma_{\times}(\bm{x})$
(6) $\displaystyle\nabla^{2}\gamma_{\mathrm{B}}(\bm{x})$
$\displaystyle=(\partial_{x}\partial_{x}-\partial_{y}\partial_{y})\gamma_{\times}(\bm{x})-2\partial_{x}\partial_{y}\gamma_{+}(\bm{x})\
,$ (7)

where $\bm{x}$ is the configuration space position and
$\partial_{x}\partial_{x}\equiv\partial^{2}/\partial x^{2}$, and so on.

This decomposition takes a simpler form in Fourier space. We define Fourier
space coordinates to be $\bm{k}=(k_{x},k_{y},k_{z})$ and choose the
$k_{z}$-axis to be along the line of sight.

The derivatives in Eqs. (6) and (7) change to multiplicative factors so we get

$\displaystyle k_{xy}^{2}\tilde{\gamma}_{\mathrm{E}}(\bm{k})$
$\displaystyle=(k_{x}^{2}-k_{y}^{2})\tilde{\gamma}_{+}(\bm{k})+2k_{x}k_{y}\tilde{\gamma}_{\times}(\bm{k})$
(8) $\displaystyle k_{xy}^{2}\tilde{\gamma}_{\mathrm{B}}(\bm{k})$
$\displaystyle=(k_{x}^{2}-k_{y}^{2})\tilde{\gamma}_{\times}(\bm{k})-2k_{x}k_{y}\tilde{\gamma}_{+}(\bm{k})\
,$ (9)

where $k_{xy}^{2}=k_{x}^{2}+k_{y}^{2}$ and $\tilde{\gamma}_{+}$ and
$\tilde{\gamma}_{\times}$ are the Fourier transforms of Eqs. (4) and (5).
Alternatively these expressions can be written in terms of the angle
$\phi=\tan^{-1}(k_{x}/k_{y})$ giving (Kurita et al., 2021; Shi et al., 2021)

$\displaystyle\tilde{\gamma}_{\mathrm{E}}(\bm{k})$
$\displaystyle=\tilde{\gamma}_{+}(\bm{k})\cos{2\phi}+\tilde{\gamma}_{\times}(\bm{k})\sin{2\phi}$
(10) $\displaystyle\tilde{\gamma}_{\mathrm{B}}(\bm{k})$
$\displaystyle=\tilde{\gamma}_{\times}(\bm{k})\cos{2\phi}-\tilde{\gamma}_{+}(\bm{k})\sin{2\phi}\
.$ (11)

Even though $\tilde{\gamma}_{\mathrm{E}}$ and $\tilde{\gamma}_{\mathrm{B}}$
are measured in the plane perpendicular to the line of sight, the wavevector
$\bm{k}$ is three-dimensional:
$\bm{k}=k(\sqrt{1-\mu^{2}}\cos{\phi},\sqrt{1-\mu^{2}}\sin{\phi},\mu)$, where
$\cos^{-1}(\mu)$ is the angle between $\bm{k}$ and the $k_{z}$-axis (Blazek et
al., 2015; Kurita et al., 2021).

###  3.3 Intrinsic alignment power spectra and bispectra

Eqs. (10) and (11) lead directly to the intrinsic alignment power spectrum
between the Fourier transforms of the E-mode shear,
$\tilde{\gamma}_{\mathrm{E}}$, and matter density contrast, $\tilde{\delta}$,
and also the auto power spectra of $\tilde{\gamma}_{\mathrm{E}}$ and
$\tilde{\gamma}_{\mathrm{B}}$:

$\displaystyle\langle\tilde{\delta}(\bm{k})\tilde{\gamma}_{\mathrm{E}}(\bm{k}^{\prime})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}+\bm{k}^{\prime})P_{\delta\mathrm{E}}(\bm{k})$
(12)
$\displaystyle\langle\tilde{\gamma}_{\mathrm{E}}(\bm{k})\tilde{\gamma}_{\mathrm{E}}(\bm{k^{\prime}})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}+\bm{k}^{\prime})P_{\mathrm{E}\mathrm{E}}(\bm{k})$
(13)
$\displaystyle\langle\tilde{\gamma}_{\mathrm{B}}(\bm{k})\tilde{\gamma}_{\mathrm{B}}(\bm{k^{\prime}})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}+\bm{k}^{\prime})P_{\mathrm{B}\mathrm{B}}(\bm{k})\
,$ (14)

where $\langle\rangle$ denotes the ensemble average and
$\delta_{\mathrm{D}}^{3}$ is the three-dimensional Dirac delta function. From
parity considerations these are the only possible non-zero intrinsic shear
power spectra (Stebbins, 1996; Kamionkowski et al., 1998; Crittenden et al.,
2002; Schneider et al., 2002). Power spectra involving a single B-mode shear
will switch sign if $\bm{k}\rightarrow-\bm{k}$ which is physically impossible
unless the spectra are zero. Theoretically the B-mode auto power spectrum also
vanishes to first order but in practice it may be non-zero because of Poisson
shot noise due to the finite sampling of the halo positions, which is also
present in the E-mode auto power spectrum. (See Blazek et al. (2019) and
Kurita et al. (2021) for more detailed discussions of shape noise).

The formalism of Eqs. (12) to (14) can be extended to intrinsic alignment
bispectra. Again, any bispectrum which includes an odd number of B-modes can
be expected to be zero by parity arguments so there are five possible non-zero
bispectra:

$\displaystyle\langle\tilde{\delta}(\bm{k}_{1})\tilde{\delta}(\bm{k}_{2})\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{3})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\delta\delta\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
(15)
$\displaystyle\langle\tilde{\delta}(\bm{k}_{1})\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{2})\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{3})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\delta\mathrm{E}\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
(16)
$\displaystyle\langle\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{1})\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{2})\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{3})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\mathrm{E}\mathrm{E}\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
(17)
$\displaystyle\langle\tilde{\delta}(\bm{k}_{1})\tilde{\gamma}_{\mathrm{B}}(\bm{k}_{2})\tilde{\gamma}_{\mathrm{B}}(\bm{k}_{3})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\delta\mathrm{B}\mathrm{B}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
(18)
$\displaystyle\langle\tilde{\gamma}_{\mathrm{E}}(\bm{k}_{1})\tilde{\gamma}_{\mathrm{B}}(\bm{k}_{2})\tilde{\gamma}_{\mathrm{B}}(\bm{k}_{3})\rangle$
$\displaystyle=(2\mathrm{\pi})^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\mathrm{E}\mathrm{B}\mathrm{B}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\
.$ (19)

###  3.4 Methodology for measuring power spectra and bispectra from
simulations

To measure the power spectra and bispectra we follow the methodology in Kurita
et al. (2021). We refer the reader to Kurita et al. (2021) and Shi et al.
(2021) for more detailed descriptions of the measurement methodology. To
measure power spectra we use the publicly available package
nbodykit777https://github.com/bccp/nbodykit (Hand et al., 2018) which provides
a wide range of tools to analyse cosmological simulations. For bispectra we
adapt bskit888https://github.com/sjforeman/bskit, developed by Foreman et al.
(2020). This package is based on nbodykit together with the Fast Fourier
Transform-based bispectrum measurement algorithm presented in Tomlinson et al.
(2019). In both cases we incorporate Eqs. (10) and (11) into the existing code
in order to measure the IA spectra.

In measuring from simulations we align all axes in Fourier space with the
simulation box sides. We then define a regular grid of size 5123 within the
box and assign the subhalo shape measurements to the grid using the cloud-in-
cell interpolation algorithm. We investigated the alternative triangle shaped
cloud assignment method but found it made little difference to the results. We
confine our bispectrum measurements to equilateral triangles and to a
representative isosceles configuration whose sides have magnitudes in the
ratio 2:2:1.

To quantify uncertainty in the estimates, we divide the simulation box into
$3^{3}$ subboxes, and estimate standard errors using jackknife sampling,
excluding each subbox in turn.

In all subsequent results we show IA auto power spectra with Poisson shot
noise subtracted, but without any allowance for non-linear effects due to
intrinsic alignments of shapes. Kurita et al. (2021) explored the latter in
detail but found it to be only 5-10 per cent of Poisson shot noise in their
halo sample from the Dark Quest simulations. Thus we measure shot noise in the
IA power spectra as $\epsilon_{\mathrm{rms}}^{2}/n_{\mathrm{eff}}$ where
$n_{\mathrm{eff}}=n_{\mathrm{h}}/L_{\mathrm{box}}^{3}$ is the effective number
density within the simulation box with side length $L_{\mathrm{box}}$ and
$n_{\mathrm{h}}$ is the number of subhalos in the box.

The measured shear ‘fields’ are weighted by the number density of halos.
Density weighting is important because halos are biased tracers of the matter
density field and ellipticity/shear measurements can only be made at the
positions where halos exist.

##  4 Analytical modelling

Two-point IA statistics are commonly modelled by the linear alignment model
(Hirata & Seljak, 2004). This model assumes that the intrinsic ellipticity of
a halo is linearly related to the local quadrupole of the gravitational
potential at the redshift at which the halo formed. Thus in Fourier space we
can write

$\displaystyle\tilde{\gamma}_{({+/\times})}(\bm{k},z)$
$\displaystyle=-A_{\mathrm{IA}}f_{({+/\times})}\frac{C_{1}\Omega_{\mathrm{m}}\rho_{\mathrm{cr}}}{D(z)}\tilde{\delta}(\bm{k},z)\
,$ (20)

where $\Omega_{\mathrm{m}}$ is the total matter density parameter,
$\rho_{\mathrm{cr}}$ is the critical density at the present day and $D(z)$ is
the linear growth factor. The functions $f_{(+/\times)}$ are defined as

$\displaystyle f_{+}$ $\displaystyle=(1-\mu^{2})\cos 2\phi\ ,$ (21)
$\displaystyle f_{\times}$ $\displaystyle=(1-\mu^{2})\sin 2\phi\ ,$ (22)

where, as before, $\cos^{-1}(\mu)$ is the angle between $\bm{k}$ and the
$k_{z}$-axis and $\phi=\tan^{-1}(k_{x}/k_{y})$. The parameter $C_{1}$ in Eq.
(20) is a normalisation factor which in principle can be determined from
observations or simulations. The amplitude $A_{\mathrm{IA}}$ quantifies the
magnitude of the intrinsic alignment effect. This quantity, commonly used in
cosmological inference, is what we are particularly interested in.

Substituting from Eq. (20) into Eqs. (10) and (11) we have

$\displaystyle\tilde{\gamma}_{\mathrm{E}}(\bm{k})$
$\displaystyle=f_{\mathrm{IA}}\tilde{\delta}(\bm{k})\ ,$ (23)
$\displaystyle\tilde{\gamma}_{\mathrm{B}}(\bm{k})$ $\displaystyle=0\ ,$ (24)

where we define $f_{\mathrm{IA}}$ as

$\displaystyle f_{\mathrm{IA}}$
$\displaystyle=-A_{\mathrm{IA}}(1-\mu^{2})\frac{C_{1}\Omega_{\mathrm{m}}\rho_{\mathrm{cr}}}{D(z)}\
.$ (25)

From Eq. (23) we can obtain the three-dimensional E-mode intrinsic alignment
power spectra

$\displaystyle P_{\delta\mathrm{E}}(k)$
$\displaystyle=f_{\mathrm{IA}}P_{\mathrm{NL}}(k)\ ,$ (26) $\displaystyle
P_{\mathrm{E}\mathrm{E}}(k)$
$\displaystyle=f_{\mathrm{IA}}^{2}P_{\mathrm{NL}}(k)\ ,$ (27)

where we have used the non-linear matter power spectrum, $P_{\mathrm{NL}}(k)$,
as suggested by Bridle & King (2007). This modification, known as the non-
linear alignment model, has been found empirically to improve the fit of the
model at non-linear scales. Given the definition of $f_{\mathrm{IA}}$, it is
clear that $P_{\delta\mathrm{E}}(k)$ is negative. Theoretically the model
predicts that B-mode power spectra are zero to first order, although as
discussed previously $P_{\mathrm{BB}}$ may be non-zero in practice due to
shape noise.

Note that Eqs. (26) and (27) also depend on redshift but for brevity we have
omitted the $z$-dependence, here and in all subsequent related equations.

###  4.1 Extension of linear alignment model to bispectra

To develop a similar analytical model for three-point statistics we need to
relate the intrinsic alignment bispectra to the non-linear matter bispectrum.
It is most straightforward to use a fitting function for the matter bispectrum
based on tree-level perturbation theory, for example those given in
Scoccimarro & Couchman (2001) and Gil-Marín et al. (2012). These have the
generic form

$\displaystyle B_{\delta\delta\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=2\,[F_{2}^{\mathrm{eff}}(\bm{k}_{1},\bm{k}_{2})P_{\mathrm{NL}}(k_{1})P_{\mathrm{NL}}(k_{2})+2\
\text{perms.}]\ ,$ (28)

where $F_{2}^{\mathrm{eff}}(\bm{k}_{1},\bm{k}_{2})$ is a modification of the
standard perturbation theory kernel (Bernardeau et al., 2002).

This leads directly to expressions for the intrinsic alignment bispectra which
are in the spirit of the two-point non-linear alignment model (Pyne &
Joachimi, 2021):

$\displaystyle B_{\delta\delta\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=2\,\Big{[}f_{\mathrm{IA}}^{2}F_{2}^{\mathrm{eff}}(\bm{k}_{1},\bm{k}_{2})P_{\mathrm{NL}}(k_{1})P_{\mathrm{NL}}(k_{2})$
(29) $\displaystyle\hskip
11.38092pt+f_{\mathrm{IA}}F_{2}^{\mathrm{eff}}(\bm{k}_{2},\bm{k}_{3})P_{\mathrm{NL}}(k_{2})P_{\mathrm{NL}}(k_{3})$
$\displaystyle\hskip
11.38092pt+f_{\mathrm{IA}}F_{2}^{\mathrm{eff}}(\bm{k}_{3},\bm{k}_{1})P_{\mathrm{NL}}(k_{3})P_{\mathrm{NL}}(k_{1})\Big{]}\
,$ $\displaystyle
B_{\delta\mathrm{E}\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=2\,\Big{[}f_{\mathrm{IA}}^{3}F_{2}^{\mathrm{eff}}(\bm{k}_{1},\bm{k}_{2})P_{\mathrm{NL}}(k_{1})P_{\mathrm{NL}}(k_{2})$
(30) $\displaystyle\hskip
11.38092pt+f_{\mathrm{IA}}^{2}F_{2}^{\mathrm{eff}}(\bm{k}_{2},\bm{k}_{3})P_{\mathrm{NL}}(k_{2})P_{\mathrm{NL}}(k_{3})$
$\displaystyle\hskip
11.38092pt+f_{\mathrm{IA}}^{3}F_{2}^{\mathrm{eff}}(\bm{k}_{3},\bm{k}_{1})P_{\mathrm{NL}}(k_{3})P_{\mathrm{NL}}(k_{1})\Big{]}\
,$ $\displaystyle
B_{\mathrm{E}\mathrm{E}\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=f_{\mathrm{IA}}^{4}B_{\delta\delta\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\
.$ (31)

The bispectrum $B_{\mathrm{E}\mathrm{E}\mathrm{E}}$ is positive but the signs
of $B_{\delta\delta\mathrm{E}}$ and $B_{\delta\mathrm{E}\mathrm{E}}$ depend on
the triangle configuration.

For equilateral triangles Eqs. (29) to (31) reduce to

$\displaystyle B_{\delta\delta\mathrm{E}}$
$\displaystyle=\frac{1}{3}\Big{[}f_{\mathrm{IA}}^{2}+2f_{\mathrm{IA}}\Big{]}B_{\delta\delta\delta}\
,$ (32) $\displaystyle B_{\delta\mathrm{E}\mathrm{E}}$
$\displaystyle=\frac{1}{3}\Big{[}2f_{\mathrm{IA}}^{3}+f_{\mathrm{IA}}^{2}\Big{]}B_{\delta\delta\delta}\
,$ (33) $\displaystyle B_{\mathrm{E}\mathrm{E}\mathrm{E}}$
$\displaystyle=f_{\mathrm{IA}}^{4}B_{\delta\delta\delta}\ ,$ (34)

where we have omitted the $\bm{k}$ arguments for brevity.

In this work we use the fitting formula from Gil-Marín et al. (2012) in Eqs.
(29) to (31) with the non-linear matter power spectrum estimated by the
fitting formula in Takahashi et al. (2012). We note, however, that the
formulas of both Scoccimarro & Couchman (2001) and Gil-Marín et al. (2012) are
known to have deficiencies: they are fitted over limited $k$ scales and have
been found to be inaccurate for squeezed triangles (Namikawa et al., 2019).
These issues were explored in detail by Takahashi et al. (2020) who developed
a new formula, Bihalofit, based on the halo model. It is more accurate than
the perturbation theory based fitting functions over a wider range of
wavenumbers and redshifts but does not lend itself to use in our model.
However, as we show in Appendix A, for scales up to at least $k\approx
3\,h\mathrm{Mpc}^{-1}$ there is good agreement between the matter bispectrum
estimated from Gil-Marín et al. (2012) and from Takahashi et al. (2020) and
between these results and our measurements from simulations.

##  5 Results

We present our results as follows. First, in Sect. 5.1, we discuss the halo
ellipticity distributions which underlie our intrinsic alignment measurements.
In Sect. 5.2 we present our measurements of intrinsic alignment power spectra
and bispectra from simulations, and discuss consistency with previous power
spectrum results and between our two-point and three-point measurements. In
the rest of Sect. 5 we use these simulation measurements to validate our
analytical models. In Sects. 5.3 and 5.4 we focus on IA power spectra, using
the NLA model to estimate intrinsic alignment amplitudes and explore their
mass dependence. These power spectrum results establish our methodology and
provide context for our principal results from bispectra based on our three-
point analytical IA model. Results for equilateral triangles are given in
Sect. 5.5 and compared to power spectrum results in 5.6. Finally, in Sect. 5.7
we discuss non-equilateral IA bispectra.

###  5.1 Distribution of halo ellipticities by mass

Previous measurements of halo intrinsic alignments from simulations, for
example Jing (2002), Lee et al. (2008), Xia et al. (2017), Piras et al. (2018)
and Kurita et al. (2021), have found that intrinsic alignment increases with
increasing halo mass and to a lesser extent with redshift. Similar trends have
been noted from measurements of galaxy shapes from simulations (Tenneti et
al., 2014, 2015), and in survey data for galaxies (Joachimi et al., 2011;
Joachimi et al., 2013; Singh et al., 2015) and clusters (Van Uitert &
Joachimi, 2017). In view of this we split our measured subhalos into four mass
bins, with each bin spanning one order of magnitude from $10^{10}$ to
$10^{14}\ h^{-1}\mathrm{M}_{\odot}$.

Figure 1 shows the distribution of ellipticity, defined as
$\epsilon=\sqrt{\epsilon_{+}^{2}+\epsilon_{\times}^{2}}$, in each of the four
bins at $z=0.5$. (Similar distributions are found for other redshifts.) Dotted
vertical lines show the median ellipticity in each bin. This confirms previous
authors’ findings: higher-mass halos are more elliptical. It is noticeable
that the two lowest-mass bins have similar ellipticities, which are lower than
those of the two high-mass bins.

Figure 1: Distribution of dark matter subhalos by projected ellipticity
$\epsilon=\sqrt{\epsilon_{+}^{2}+\epsilon_{\times}^{2}}$ in four mass bins at
$z=0.5$. Vertical dotted lines indicate the median ellipticity in each bin.

###  5.2 IA power spectra and bispectra measured from simulations

Figure 2 shows our measured intrinsic alignment power spectra,
$P_{\delta\mathrm{E}}$, $P_{\mathrm{E}\mathrm{E}}$ and
$P_{\mathrm{B}\mathrm{B}}$, defined by Eqs. (12) to (14), at three redshifts,
with the non-linear matter power spectrum $P_{\delta\delta}$ shown for
comparison. In this and all similar figures we show the absolute value of the
spectra. $P_{\mathrm{B}\mathrm{B}}$ is essentially equal to the shot noise at
all but the largest scales considered here. The non-vanishing cross-power
spectrum $P_{\delta\mathrm{E}}$ confirms that subhalo shapes are correlated
with the matter overdensity field at all scales. The power spectra show little
variation across redshifts. These results are consistent with IA power
spectrum measurements in Kurita et al. (2021). (See their Figs. 2 and 4).

Figure 2:  Intrinsic alignment power spectra
$P_{\delta\mathrm{E}},P_{\mathrm{E}\mathrm{E}}$ and $P_{\mathrm{B}\mathrm{B}}$
measured from simulations. The matter power spectrum $P_{\delta\delta}$ is
also shown for reference. Shot noise has been subtracted from the auto power
spectra. Shaded areas are 68% confidence intervals.

Figures 3 and 4 show similar results for measured E-mode and B-mode bispectra
respectively (Eqs. (15) to (19)), for equilateral triangles and for isosceles
triangles with sides in the ratio 2:2:1. The B-mode bispectra have very low
signal-to-noise ratios (see Appendix B) and we do not consider them further in
this work. However the E-mode bispectra, in particular
$B_{\delta\delta\mathrm{E}}$, have relatively strong signals, especially for
isosceles configurations.

There is no convenient benchmark with which to compare our bispectrum results.
The only previous measurements of three-point intrinsic alignments from
simulations were by Semboloni et al. (2008). Unfortunately it is difficult to
compare our results with theirs for a number of reasons: they worked in
configuration space and measured aperture mass statistics, and they focused on
the magnitudes of intrinsic alignments relative to lensing signals, rather
than on the strength of the IA signal itself.

Figure 3:  Intrinsic alignment bispectra
$B_{\delta\delta\mathrm{E}},B_{\delta\mathrm{E}\mathrm{E}}$ and
$B_{\mathrm{E}\mathrm{E}\mathrm{E}}$ measured from simulations. Top:
Equilateral triangles. Bottom: Isosceles triangles with sides in the ratio
2:2:1. The bispectrum $B_{\delta\delta\delta}$ is also shown for reference.
Shaded areas are 68% confidence intervals.  Figure 4:  Intrinsic alignment
bispectra $B_{\delta\delta\mathrm{B}}$ and
$B_{\mathrm{E}\mathrm{B}\mathrm{B}}$ measured from simulations. Top:
Equilateral triangles. Bottom: Isosceles triangles with sides in the ratio
2:2:1. The bispectrum $B_{\delta\delta\delta}$ is also shown for reference.
Shaded areas are 68% confidence intervals.

Figure 5:  Intrinsic alignment spectra in four mass bins. Top:
$P_{\delta\mathrm{E}}$, Centre: $P_{\mathrm{EE}}$, Bottom:
$B^{\mathrm{equi}}_{\delta\delta\mathrm{E}}$.

In Fig. 5 we split the IA power spectra $P_{\delta\mathrm{E}}$ and
$P_{\mathrm{EE}}$ and the bispectrum $B_{\delta\delta\mathrm{E}}$ between mass
bins. In all three cases the magnitudes of the spectra increase with mass.
Again we note that the spectra in the two lowest-mass bins are similar to each
other and generally smaller than those in higher-mass bins. The results for
$P_{\mathrm{EE}}$ are noisy and the relationship with mass is not as
pronounced, especially for $z=0$.

###  5.3 Intrinsic alignment amplitude from power spectra

Having shown that E-mode spectra are detectable from the simulations, we now
estimate the intrinsic alignment amplitude $A_{\mathrm{IA}}$ in Eq. (20) from
simulation measurements. We first consider power spectra. From Eqs. (26) and
(27) we can obtain approximate estimates of $f_{\mathrm{IA}}$, and hence
$A_{\mathrm{IA}}$, from the ratio
$P_{\delta\mathrm{E}}(k)/P_{\delta\delta}(k)$ or from
$\sqrt{P_{\mathrm{E}\mathrm{E}}(k)/P_{\delta\delta}(k)}$. Alternatively,
following Kurita et al. (2021), we can use least-squares minimization to fit
$A_{\mathrm{IA}}$ from these equations. To do this we find the value of the
parameter $\hat{A}_{\mathrm{IA}}$ which minimizes $\chi^{2}$ given by

$\displaystyle\chi^{2}$ $\displaystyle={\sum}_{k}\
{\frac{[R(k)-F(\hat{A}_{\mathrm{IA}})]^{2}}{\sigma^{2}_{R}(k)}}\ .$ (35)

Here, in the first case, based on $P_{\delta\mathrm{E}}$,
$R(k)=P_{\delta\mathrm{E}}(k)/P_{\delta\delta}(k)$ and

$\displaystyle F(\hat{A}_{\mathrm{IA}})$
$\displaystyle=\int_{0}^{1}\,(1-\mu^{2})\
\mathrm{d}\mu\,c(z)\hat{A}_{\mathrm{IA}}=\frac{2}{3}c(z)\hat{A}_{\mathrm{IA}}\
,$ (36)

where $c(z)=C_{1}\Omega_{\mathrm{m}}\rho_{\mathrm{cr}}/D(z)$ (see Eq. (25))
and $\sigma^{2}_{R}(k)$ is the variance of $R(k)$ calculated using jackknife
sampling from 27 simulation subboxes. To calculate $F(\hat{A}_{\mathrm{IA}})$
we obtain the growth factor from the cosmological parameter estimation code
CosmoSIS999https://bitbucket.org/joezuntz/cosmosis/wiki/Home (Zuntz et al.,
2015). We use the value of $C_{1}$ derived by Bridle & King (2007) which is
$5\times 10^{-14}\,h^{-2}\mathrm{M}_{\odot}^{-1}\mathrm{Mpc}^{3}$, leading to
$C_{1}\rho_{\mathrm{cr}}=0.0134$ (Joachimi et al., 2011). All cosmological
parameters values are identical to those used in the simulations.

In the second case, based on $P_{\mathrm{EE}}$,
$R(k)=\sqrt{P_{\mathrm{EE}}(k)/P_{\delta\delta}(k)}$ and

$\displaystyle F(\hat{A}_{\mathrm{IA}})$
$\displaystyle=\int_{0}^{1}\,(1-\mu^{2})^{2}\
\mathrm{d}\mu\,c(z)\hat{A}_{\mathrm{IA}}=\frac{8}{15}c(z)\hat{A}_{\mathrm{IA}}\
,$ (37)

with $c(z)$ and $\sigma^{2}_{R}(k)$ defined as before.

Figure 6 compares the values of $A_{\mathrm{IA}}$ obtained from ratios of
power spectra to those obtained using Eq. (35). Results are for
$k_{\mathrm{min}}=0.1\,h\mathrm{Mpc}^{-1}$ and varying $k_{\mathrm{max}}$.
There is good agreement for all values of $k_{\mathrm{max}}$ for estimates
from $P_{\delta\mathrm{E}}$ and for $k_{\mathrm{max}}\lesssim
0.5\,h\mathrm{Mpc}^{-1}$ for estimates from $P_{\mathrm{E}\mathrm{E}}$. This
validates the NLA model. The results for $P_{\mathrm{E}\mathrm{E}}$ are
consistently higher than those for $P_{\delta\mathrm{E}}$, suggesting that
this discrepancy may be due to an unexplained systematic trend. However, as
discussed in Appendix B, the signal-to-noise ratio for
$P_{\mathrm{E}\mathrm{E}}$ is low and thus the results are less robust than
those for $P_{\delta\mathrm{E}}$.

Figure 6:  The intrinsic alignment amplitude $A_{\mathrm{IA}}$ as a function
of $k_{\mathrm{max}}$ estimated approximately as the ratio between the IA
power spectra and the matter power spectrum and also fitted using Eq. (35)
over the range $[0.1,k_{\mathrm{max}}]\ h\mathrm{Mpc}^{-1}$.

###  5.4 Dependence of IA amplitude on mass

Since the intrinsic alignment power spectra depend strongly on subhalo mass,
we also expect $A_{\mathrm{IA}}$ to depend on mass. We therefore obtain
estimates of $A_{\mathrm{IA}}$ from $P_{\delta\mathrm{E}}$ and
$P_{\mathrm{E}\mathrm{E}}$ for each of our four mass bins. We assume that the
same linear alignment model applies for all halo masses, with any mass-
dependence being absorbed into the estimated amplitude. Alternative ways in
which the mass-dependence has been incorporated include inserting a mass-
dependent halo bias term into the linear alignment model (Xia et al., 2017)
and using a virial argument to derive a scaling by mass (Piras et al., 2018).

We perform the $\chi^{2}$ minimization using Eq. (35) for each of the four
mass bins and for three redshifts, and also vary the maximum $k$ value used.
The results are shown in Fig. 7. The values of $A_{\mathrm{IA}}$ obtained from
$P_{\delta\mathrm{E}}$ are consistent with Fig. 6 in Kurita et al. (2021),
given slightly different mass bins and $k$ and $z$ ranges.

Figure 7:  The intrinsic alignment amplitude $A_{\mathrm{IA}}$ estimated from
$P_{\delta\mathrm{E}}$, $P_{\mathrm{E}\mathrm{E}}$ and
$B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ for four subhalo mass ranges, as
a function of $k_{\mathrm{max}}$. Typical uncertainties in the fits are given
in Table 1.

We also fit a power law of the form

$\displaystyle A_{\mathrm{IA}}$ $\displaystyle\propto M_{\mathrm{h}}^{\beta}\
,$ (38)

where $M_{\mathrm{h}}$ is the mean halo mass per bin. We estimate this
relationship using both $P_{\delta\mathrm{E}}$ and $P_{\mathrm{E}\mathrm{E}}$.
The results are shown in Fig. 8. At all redshifts the two estimates are
reasonably consistent for $k_{\mathrm{max}}\leq 2\ h\,\mathrm{Mpc}^{-1}$, but
show some variability for higher $k_{\mathrm{max}}$.

We can compare our estimated power spectrum slopes in Fig. 8 with those found
by Piras et al. (2018) from the Millennium simulation. These authors obtained
$\beta\approx 0.35$ for $z=0.46$, with a slightly lower-mass sample of halos
($10^{11.36}<M_{\mathrm{h}}<10^{13.36}\ h^{-1}\mathrm{M}_{\odot}$). Our value
$\beta\approx 0.43$ for $k_{\mathrm{max}}\approx 2\ h\,\mathrm{Mpc}^{-1}$ at
$z=0.5$ is consistent with this. We note, however, that Piras et al. (2018)
repeated their analysis using observational data and found $\beta\approx
0.56$. They suggested several possible reasons for this difference, including
the effect of baryons on halo shapes and the relative strength of the stellar
and dark matter signals in different mass bins.

Figure 8:  Slope $\beta$ from a power law relationship $A_{\mathrm{IA}}\propto
M_{\mathrm{h}}^{\beta}$, estimated from $P_{\delta\mathrm{E}}$,
$P_{\mathrm{E}\mathrm{E}}$ and $B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$.
Typical uncertainties in the fits are given in Table 1.

###  5.5 Intrinsic alignment amplitude from equilateral bispectra

We now turn to our main aim: to investigate whether our three-point analytical
model is consistent with the two-point NLA model over the non-linear scales of
interest. We start with the equilateral bispectrum
$B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ given by Eq. (32). We later
expand this to non-equilateral versions of $B_{\delta\delta\mathrm{E}}$, but
do not consider other IA bispectra because we find these have insufficient
signal to provide robust estimates of the amplitude $A_{\mathrm{IA}}$.

We again use Eq. (35) to fit a parameter $\hat{A}_{\mathrm{IA}}$ but in this
case $R(k)=B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}/B_{\delta\delta\delta}$.
Also, using Eq. (32), we now have

$\displaystyle F(\hat{A}_{\mathrm{IA}})$
$\displaystyle=\int_{0}^{1}\,(1-\mu^{2})\
\mathrm{d}\mu\,\frac{[f_{\mathrm{IA}}^{2}+2f_{\mathrm{IA}}]}{3}$ (39)
$\displaystyle=\frac{2}{9}[(c(z)\hat{A}_{\mathrm{IA}})^{2}+2c(z)\hat{A}_{\mathrm{IA}}]\
,$ (40)

where $c(z)$ and $\hat{A}_{\mathrm{IA}}$ are defined as in Sect. 5.3. In
contrast to the power spectra, there is no simple relationship between the
ratio $B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}/B_{\delta\delta\delta}$ and
the intrinsic alignment amplitude so we cannot produce a plot similar to Fig.
6.

As shown in Fig. 5, $B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ increases
with mass, albeit not as consistently as $P_{\delta\mathrm{E}}$, so we again
obtain estimates of $A_{\mathrm{IA}}$ for each mass bin, and fit a power law
to this. Figures 7 and 8 show the resulting estimates of $A_{\mathrm{IA}}$ and
of the slope $\beta$. Bispectrum measurements for small $k$ values are noisy
so we were only able to fit for $k\gtrsim 0.7\ h\,\mathrm{Mpc}^{-1}$ in
contrast to $k\gtrsim 0.1\ h\,\mathrm{Mpc}^{-1}$ for the power spectrum.

###  5.6 Comparison of results from power spectra and equilateral bispectra

Table 1 summarises our main results. It shows estimates of $A_{\mathrm{IA}}$
derived from $P_{\delta\mathrm{E}}$, $P_{\mathrm{E}\mathrm{E}}$ and
$B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ for each mass bin. Also in this
table are the estimated power-law slopes $\beta$ from Eq. (38) and the mean
value of $A_{\mathrm{IA}}$, weighted by the number of subhalos per mass bin.
For illustrative purposes we choose $k_{\mathrm{max}}=2\ h\,\mathrm{Mpc}^{-1}$
in this table, well within non-linear scales, since our models appear reliable
up to this value. This range also covers the typical fit range for cosmic
shear studies. The overall picture would be qualitatively similar with a
different $k_{\mathrm{max}}$.

The IA amplitudes obtained from $P_{\delta\mathrm{E}}$ and
$B_{\delta\delta\mathrm{E}}$ are completely consistent, underlining the
validity of the analytical modelling up to at least $k_{\mathrm{max}}=2\
h\,\mathrm{Mpc}^{-1}$ and confirming that our modelling of three-point IA
statistics is consistent with the two-point non-linear alignment model. As
discussed earlier, results from $P_{\mathrm{E}\mathrm{E}}$ agree only
qualitatively with the other results. The power-law relationship which we
obtain using $B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ is broadly
consistent with the power spectrum results but is rather lower in most cases
and more dependent on the value of $k_{\mathrm{max}}$.

Table 1: Estimates of the intrinsic alignment amplitude $A_{\mathrm{IA}}$ and the power-law exponent $\beta$ in Eq. (38) for subhalos in four mass ranges, with 68% confidence intervals, from the power spectra $P_{\delta\mathrm{E}}$ and $P_{\mathrm{E}\mathrm{E}}$ and the bispectrum $B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ with $k_{\mathrm{max}}=2\ h\,\mathrm{Mpc}^{-1}.$ Also shown are the mean $A_{\mathrm{IA}}$ for all masses, calculated as the average over mass bins weighted by the number of halos per bin. For both power spectra $k_{\mathrm{min}}\approx 0.1\ h\,\mathrm{Mpc}^{-1}$ and for the bispectrum $k_{\mathrm{min}}\approx 0.7\ h\,\mathrm{Mpc}^{-1}$.  |  | Subhalo mass range $[h^{-1}\mathrm{M}_{\odot}]$ |  |
---|---|---|---|---
| Redshift | $10^{10}-10^{11}$ | $10^{11}-10^{12}$ | $10^{12}-10^{13}$ | $10^{13}-10^{14}$ | Slope $\beta$ | Mean $A_{\mathrm{IA}}$
From power spectrum $P_{\delta\mathrm{E}}$ | 0.0 | 0.74$\pm$0.01 | 0.86$\pm$ 0.01 | 2.40$\pm$ 0.03 | 11.13$\pm$0.12 |  0.392$\pm$0.002 | 0.91$\pm$0.01
| 0.5 | 1.02$\pm$0.01 | 1.31$\pm$0.01 | 4.26$\pm$0.04 | 19.38$\pm$0.14 |  0.432$\pm$0.001 | 1.36$\pm$0.01
| 1.0 | 1.59$\pm$0.01 | 2.10$\pm$0.02 | 6.68$\pm$0.06 | 28.87$\pm$0.21 |  0.427$\pm$0.001 | 2.13$\pm$0.02
From power spectrum $P_{\mathrm{E}\mathrm{E}}$ | 0.0 | 1.40$\pm$0.02 | 1.61$\pm$ 0.02 | 6.84$\pm$ 0.09 | 45.08$\pm$0.81 |  0.507$\pm$0.002 | 1.90$\pm$0.02
| 0.5 | 1.96$\pm$0.02 | 2.25$\pm$0.02 | 9.14$\pm$0.09 | 63.43$\pm$1.01 |  0.507$\pm$0.001 | 2.64$\pm$0.02
| 1.0 | 2.69$\pm$0.02 | 3.21$\pm$0.02 | 13.54$\pm$0.15 | 112.53$\pm$1.47 |  0.543$\pm$0.001 | 3.81$\pm$0.03
From bispectrum $B_{\delta\delta\mathrm{E}}$ | 0.0 | 0.94$\pm$0.01 | 0.99$\pm$0.01 | 2.33$\pm$0.03 | 9.74$\pm$0.16 |  0.332$\pm$0.003 | 1.06$\pm$0.01
| 0.5 | 1.07$\pm$0.01 | 1.25$\pm$0.01 | 3.79$\pm$0.05 | 17.59$\pm$0.23 |  0.407$\pm$0.002 | 1.34$\pm$0.01
| 1.0 | 1.58$\pm$0.02 | 2.07$\pm$0.02 | 5.94$\pm$0.09 | 33.88$\pm$0.39 |  0.446$\pm$0.002 | 2.11$\pm$0.03

###  5.7 Analytical predictions for isosceles bispectra

It is worth exploring whether our analytical model can also be applied to non-
equilateral bispectra since these can contain more information than
equilateral configurations (Appendix B). In this case it is not possible to
obtain simplified expressions similar to those in Eqs. (32) to (34) because
the perturbation theory kernels do not cancel out. Instead we consider whether
the intrinsic alignment amplitudes predicted from
$B^{\mathrm{equi}}_{\delta\delta\mathrm{E}}$ are also valid for non-
equilateral bispectra. For illustration we again consider isosceles triangles
with one side with magnitude $k$ and two sides with magnitude $2k$.

We take the estimates of $A_{\mathrm{IA}}$ obtained from
$B^{\mathrm{equi}}_{\delta\delta\mathrm{E}}$ in the final column of Table 1,
noting that these are very similar to the estimates obtained from
$P_{\delta\mathrm{E}}$. We calculate $f_{\mathrm{IA}}$ from Eq. (25) using the
same assumptions as in Sect. 5.3, obtain the non-linear matter power spectrum
from Halofit (Takahashi et al., 2012), and insert the resulting values into
Eqs. (29) to (31). In Fig. 9 we compare these predicted bispectra with those
measured from simulations, finding a reasonable fit across scales and
redshifts. This confirms that our empirical IA bispectrum model works for both
equilateral and isosceles triangle configurations.

Figure 9:  Intrinsic alignment bispectra for isosceles triangles with one side
of magnitude $k$ and two sides of magnitude $2k$. Solid lines: Measured from
simulations. Dashed lines: Calculated from Eqs. (29) to (31) with best-fitting
values of $A_{\mathrm{IA}}$ estimated from
$B^{\mathrm{equi}}_{\delta\delta\mathrm{E}}$.

###  5.8 Decorrelation of matter overdensity and E-mode fields

In Appendix C we discuss a possible phenomenological change to our two-point
and three-point analytical models which adjusts them for the correlation
between the matter overdensity field and the E-mode ‘field’. This is not
central to our main results but may be worth considering for future intrinsic
alignment modelling.

##  6 Summary and discussion

We have measured intrinsic alignment bispectra of dark matter subhalos from
the IllustrisTNG300-1 cosmological simulation suite, building on the power
spectrum methodology developed by Kurita et al. (2021). We also measured the
IA power spectra $P_{\delta\mathrm{E}}$ and $P_{\mathrm{E}\mathrm{E}}$ and
confirmed that they are consistent with results obtained by Kurita et al.
(2021) using the DarkQuest simulation suite.

At all redshifts the cumulative signal-to-noise ratios of the IA bispectra
were well below those we obtained for power spectra – for example, 10-15 per
cent at $k_{\mathrm{max}}\approx 4h\,\mathrm{Mpc}^{-1}$ compared with around
40 per cent for power spectra. The bispectrum
$B_{\mathrm{E}\mathrm{E}\mathrm{E}}$ and measured bispectra involving B-modes
have very low signal-to-noise ratios and provide no useful information.
However we found that the E-mode bispectra $B_{\delta\delta\mathrm{E}}$ and
$B_{\delta\mathrm{E}\mathrm{E}}$ do have useful information content. Signal-
to-noise ratios for the non-equilateral triangles which we studied are notably
higher than for equilateral triangles, suggesting that the common
simplification of using only equilateral triangles in bispectrum analyses is
sub-optimal (Yankelevich et al., 2022).

All the IA power spectra and bispectra we studied showed similar strong
relationships with subhalo mass. This can be traced back to the relationship
between halo ellipticity and mass. As discussed in Sect. 5.1 this mass
dependence has been noted by others in both simulations and observations. Its
origin is debated. Smith & Watts (2005) attributed it to higher mass halos
forming later so they have had less time to virialise and therefore retain
more memory of the tidal fields at the time they formed. More recently Xia et
al. (2017) postulated that the IA strength depends independently on both halo
formation time and mass, whereas Piras et al. (2018) suggested that higher-
mass halos experience stronger tidal fluctuations.

We used the standard non-linear alignment model to estimate the intrinsic
alignment amplitude from both $P_{\delta\mathrm{E}}$ and
$P_{\mathrm{E}\mathrm{E}}$, and found that our estimates are consistent with
corresponding results obtained by Kurita et al. (2021). This validated our
two-point modelling.

For IA bispectra we used the analytical model from Pyne & Joachimi (2021),
which is in the spirit of the NLA model. From this we again estimated the
intrinsic alignment amplitude, this time from equilateral bispectra. We found
that the best-fitting amplitudes $A_{\mathrm{IA}}$ obtained using
$B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$ were extremely close to those
obtained from $P_{\delta\mathrm{E}}$. It is not possible to use the same
methodology to estimate $A_{\mathrm{IA}}$ from non-equilateral triangles but
we showed that the predicted $A_{\mathrm{IA}}$ from
$B^{\mathrm{equi}}_{\delta\delta\mathrm{E}}$ produced an acceptable fit to
simulation measurements of IA bispectra of specific isosceles triangles.

We fitted power-law relationships between the estimated $A_{\mathrm{IA}}$ and
subhalo mass, obtaining almost identical relationships from
$P_{\delta\mathrm{E}}$ and $B_{\delta\delta\mathrm{E}}^{\mathrm{equi}}$. These
also agreed with the relationship found by Piras et al. (2018) using power
spectrum measurements from the Millennium simulation. It is interesting that
all our estimated IA amplitudes are similar for the two lowest-mass bins and
larger for the high-mass bins. This hints at a possible broken power-law
relationship with mass similar to the luminosity relationship found by Fortuna
et al. (2021) for a sample of luminous red galaxies from KiDS-1000 (Kuijken et
al., 2019).

Our bispectrum measurements are an advance on the early three-point IA
measurements from simulations reported in Semboloni et al. (2008). We have
been able to take advantage of the improved resolution and methodology of the
IllustrisTNG simulation suite, and have also produced Fourier space
measurements for the first time. It is difficult to compare the two sets of
results since Semboloni et al. (2008) measured aperture mass statistics and
also focused on the magnitude of the IA effect relative to cosmic shear,
whereas we considered the IA amplitude only. However, Pyne & Joachimi (2021)
estimated ratios between IA and lensing signals using the same analytical
models as in this paper. The results were not completely consistent with
Semboloni et al. (2008) although both studies found that the two-point and
three-point ratios differed significantly. The present work has shown that our
analytical models agree well with measurements from IllustrisTNG which
validates the results reported in Pyne & Joachimi (2021).

The fits between the measured and modelled IA power spectra and bispectra are
not perfect. In Appendix C we suggest a possible phenomenological modification
of the analytical models based on the correlation between the matter
overdensity and E-mode fields measured in simulations. This may be worth
considering for future IA modelling.

Overall our results demonstrate that a single physically-motivated analytical
approach can be applied to both two-point and three-point IA statistics,
enabling a cleaner separation between IA and weak lensing signals. This opens
up the prospect of using such a model in joint power spectrum–bispectrum
analysis. Pyne & Joachimi (2021) showed that such analysis can allow self-
calibration to mitigate intrinsic alignment contamination of weak lensing data
in forthcoming surveys. This is particularly pertinent in the light of
advances in the measurement of three-point statistics from survey data. Secco
et al. (2022a) recently reported high signal-to-noise detections of three-
point shear correlations and aperture mass statistics in the first three years
of data from the Dark Energy Survey. These measurements lay the foundations
for joint two- and three-point cosmological analyses which will of course need
tight control of systematics such as intrinsic alignments.

This work has considered only dark matter halos and so has only limited
application to observational data. In future work we plan to build on the
measurement techniques and modelling in this paper to confirm that the three-
point approach can be extended to galaxies and to investigate the impact of
galaxy morphology and environment.

## Acknowledgements

We thank the IllustrisTNG team for making their simulation data publicly
available. SP thanks Harry Johnston and Anik Halder for helpful discussions.
This work was partially enabled by funding from the UCL Cosmoparticle
Initiative.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ade et al. (2016) Ade P. A., et al., 2016, Astronomy & Astrophysics, 594, A13
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztanaga E., Scoccimarro R., 2002, Physics reports, 367, 1
  * Bernstein & Jarvis (2002) Bernstein G., Jarvis M., 2002, The Astronomical Journal, 123, 583
  * Blazek et al. (2015) Blazek J., Vlah Z., Seljak U., 2015, Journal of Cosmology and Astroparticle Physics, 2015, 015
  * Blazek et al. (2019) Blazek J. A., MacCrann N., Troxel M., Fang X., 2019, Physical Review D, 100, 103506
  * Bridle & King (2007) Bridle S., King L., 2007, New Journal of Physics, 9, 444
  * Chisari & Dvorkin (2013) Chisari N. E., Dvorkin C., 2013, Journal of Cosmology and Astroparticle Physics, 2013, 029
  * Crittenden et al. (2002) Crittenden R. G., Natarajan P., Pen U.-L., Theuns T., 2002, The Astrophysical Journal, 568, 20
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D., 1985, The Astrophysical Journal, 292, 371
  * Foreman et al. (2020) Foreman S., Coulton W., Villaescusa-Navarro F., Barreira A., 2020, Monthly Notices of the Royal Astronomical Society, 498, 2887
  * Fortuna et al. (2021) Fortuna M. C., et al., 2021, Astronomy & Astrophysics, 654, A76
  * Fu et al. (2014) Fu L., et al., 2014, Monthly Notices of the Royal Astronomical Society, 441, 2725
  * Gil-Marín et al. (2012) Gil-Marín H., Wagner C., Fragkoudi F., Jimenez R., Verde L., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 047
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, The Astronomical Journal, 156, 160
  * Heymans et al. (2006) Heymans C., White M., Heavens A., Vale C., Van Waerbeke L., 2006, Monthly Notices of the Royal Astronomical Society, 371, 750
  * Heymans et al. (2012) Heymans C., et al., 2012, Monthly Notices of the Royal Astronomical Society, 427, 146
  * Hirata & Seljak (2004) Hirata C. M., Seljak U., 2004, Physical Review D, 70, 063526
  * Huterer et al. (2006) Huterer D., Takada M., Bernstein G., Jain B., 2006, Monthly Notices of the Royal Astronomical Society, 366, 101
  * Jing (2002) Jing Y., 2002, Monthly Notices of the Royal Astronomical Society, 335, L89
  * Joachimi et al. (2011) Joachimi B., Mandelbaum R., Abdalla F., Bridle S., 2011, Astronomy & Astrophysics, 527, A26
  * Joachimi et al. (2013) Joachimi B., Semboloni E., Bett P., Hartlap J., Hilbert S., Hoekstra H., Schneider P., Schrabback T., 2013, Monthly Notices of the Royal Astronomical Society, 431, 477
  * Joachimi et al. (2015) Joachimi B., et al., 2015, Space Science Reviews, 193, 1
  * Joachimi et al. (2021) Joachimi B., et al., 2021, Astronomy & Astrophysics, 646, A129
  * Kamionkowski et al. (1998) Kamionkowski M., Babul A., Cress C. M., Refregier A., 1998, Monthly Notices of the Royal Astronomical Society, 301, 1064
  * Kiessling et al. (2015) Kiessling A., et al., 2015, Space Science Reviews, 193, 67
  * Kirk et al. (2012) Kirk D., Rassat A., Host O., Bridle S., 2012, Monthly Notices of the Royal Astronomical Society, 424, 1647
  * Kirk et al. (2015) Kirk D., et al., 2015, Space Science Reviews, 193, 139
  * Kuijken et al. (2019) Kuijken K., et al., 2019, Astronomy & Astrophysics, 625, A2
  * Kurita et al. (2021) Kurita T., Takada M., Nishimichi T., Takahashi R., Osato K., Kobayashi Y., 2021, Monthly Notices of the Royal Astronomical Society, 501, 833
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv preprint arXiv:1110.3193
  * Lee et al. (2008) Lee J., Springel V., Pen U.-L., Lemson G., 2008, Monthly Notices of the Royal Astronomical Society, 389, 1266
  * Marinacci et al. (2018) Marinacci F., et al., 2018, Monthly Notices of the Royal Astronomical Society, 480, 5113
  * Merkel & Schäfer (2014) Merkel P. M., Schäfer B. M., 2014, Monthly Notices of the Royal Astronomical Society, 445, 2918
  * Miyatake et al. (2021) Miyatake H., et al., 2021, arXiv preprint arXiv:2101.00113
  * Naiman et al. (2018) Naiman J. P., et al., 2018, Monthly Notices of the Royal Astronomical Society, 477, 1206
  * Namikawa et al. (2019) Namikawa T., Bose B., Bouchet F. R., Takahashi R., Taruya A., 2019, Physical Review D, 99, 063511
  * Nelson et al. (2018) Nelson D., et al., 2018, Monthly Notices of the Royal Astronomical Society, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 1
  * Nishimichi et al. (2019) Nishimichi T., et al., 2019, The Astrophysical Journal, 884, 29
  * Pillepich et al. (2018) Pillepich A., et al., 2018, Monthly Notices of the Royal Astronomical Society, 475, 648
  * Piras et al. (2018) Piras D., Joachimi B., Schäfer B. M., Bonamigo M., Hilbert S., van Uitert E., 2018, Monthly Notices of the Royal Astronomical Society, 474, 1165
  * Pyne & Joachimi (2021) Pyne S., Joachimi B., 2021, Monthly Notices of the Royal Astronomical Society, 503, 2300
  * Schneider & Bridle (2010) Schneider M. D., Bridle S., 2010, Monthly Notices of the Royal Astronomical Society, 402, 2127
  * Schneider & Seitz (1995) Schneider P., Seitz C., 1995, Astronomy & Astrophysics, 294, 411
  * Schneider et al. (2002) Schneider P., van Waerbeke L., Mellier Y., 2002, Astronomy & Astrophysics, 389, 729
  * Schneider et al. (2012) Schneider M. D., Frenk C. S., Cole S., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 030
  * Scoccimarro & Couchman (2001) Scoccimarro R., Couchman H., 2001, Monthly Notices of the Royal Astronomical Society, 325, 1312
  * Scoville et al. (2007) Scoville N., et al., 2007, The Astrophysical Journal Supplement Series, 172, 1
  * Secco et al. (2022a) Secco L., et al., 2022a, arXiv preprint arXiv:2201.05227
  * Secco et al. (2022b) Secco L., et al., 2022b, Physical Review D, 105, 023515
  * Seitz & Schneider (1997) Seitz C., Schneider P., 1997, Astronomy and Astrophysics, 318, 687
  * Semboloni et al. (2008) Semboloni E., Heymans C., Van Waerbeke L., Schneider P., 2008, Monthly Notices of the Royal Astronomical Society, 388, 991
  * Semboloni et al. (2010) Semboloni E., Schrabback T., van Waerbeke L., Vafaei S., Hartlap J., Hilbert S., 2010, Monthly Notices of the Royal Astronomical Society, 410, 143
  * Shi et al. (2021) Shi J., Kurita T., Takada M., Osato K., Kobayashi Y., Nishimichi T., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 030
  * Singh et al. (2015) Singh S., Mandelbaum R., More S., 2015, Monthly Notices of the Royal Astronomical Society, 450, 2195
  * Smith & Watts (2005) Smith R. E., Watts P., 2005, Monthly Notices of the Royal Astronomical Society, 360, 203
  * Springel (2010) Springel V., 2010, Monthly Notices of the Royal Astronomical Society, 401, 791
  * Springel et al. (2001) Springel V., White S. D., Tormen G., Kauffmann G., 2001, Monthly Notices of the Royal Astronomical Society, 328, 726
  * Springel et al. (2018) Springel V., et al., 2018, Monthly Notices of the Royal Astronomical Society, 475, 676
  * Stebbins (1996) Stebbins A., 1996, arXiv preprint astro-ph/9609149
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, The Astrophysical Journal, 761, 152
  * Takahashi et al. (2020) Takahashi R., Nishimichi T., Namikawa T., Taruya A., Kayo I., Osato K., Kobayashi Y., Shirasaki M., 2020, The Astrophysical Journal, 895, 113
  * Taruya & Okumura (2020) Taruya A., Okumura T., 2020, The Astrophysical Journal Letters, 891, L42
  * Tenneti et al. (2014) Tenneti A., Mandelbaum R., Di Matteo T., Feng Y., Khandai N., 2014, Monthly Notices of the Royal Astronomical Society, 441, 470
  * Tenneti et al. (2015) Tenneti A., Singh S., Mandelbaum R., Matteo T. D., Feng Y., Khandai N., 2015, Monthly Notices of the Royal Astronomical Society, 448, 3522
  * Tomlinson et al. (2019) Tomlinson J., Jeong D., Kim J., 2019, The Astronomical Journal, 158, 116
  * Troxel & Ishak (2011) Troxel M. A., Ishak M., 2011, Monthly Notices of the Royal Astronomical Society, 419, 1804
  * Troxel & Ishak (2012) Troxel M. A., Ishak M., 2012, Monthly Notices of the Royal Astronomical Society, 423, 1663
  * Troxel & Ishak (2015) Troxel M., Ishak M., 2015, Physics Reports, 558, 1
  * Van Uitert & Joachimi (2017) Van Uitert E., Joachimi B., 2017, Monthly Notices of the Royal Astronomical Society, 468, 4502
  * Viola et al. (2014) Viola M., Kitching T., Joachimi B., 2014, Monthly Notices of the Royal Astronomical Society, 439, 1909
  * Vlah et al. (2020) Vlah Z., Chisari N. E., Schmidt F., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 025
  * Vlah et al. (2021) Vlah Z., Chisari N. E., Schmidt F., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 061
  * Xia et al. (2017) Xia Q., Kang X., Wang P., Luo Y., Yang X., Jing Y., Wang H., Mo H., 2017, The Astrophysical Journal, 848, 22
  * Yankelevich et al. (2022) Yankelevich V., McCarthy I. G., Kwan J., Stafford S. G., Liu J., 2022, arXiv preprint arXiv:2202.07680
  * Zuntz et al. (2015) Zuntz J., et al., 2015, Astronomy and Computing, 12, 45

##  Appendix A Comparison between simulations and analytical calculations of
the matter bispectrum

To confirm the suitability of the perturbation theory-based formula from Gil-
Marín et al. (2012) for our purposes, Fig. 10 compares it with our matter
bispectrum measurements at $z=1$ from the IllustrisTNG300-1 simulations. Also
shown are more accurate analytical estimates using Bihalofit (Takahashi et
al., 2020) and tree-level perturbation theory bispectra (Bernardeau et al.,
2002). We show results for equilateral triangles and for isosceles triangles
with sides in the ratio 2:2:1. All the separate results for the non-linear
bispectrum are consistent at the scales which we are interested in, with the
simulation measurements and analytical estimates from Gil-Marín et al. (2012)
only diverging significantly for $k\gtrsim 3h\,\mathrm{Mpc}^{-1}$.

Figure 10:  Comparison between the matter bispectrum at $z=1$ measured from
the IllustrisTNG300-1 simulation, calculated using the fitting formula from
Gil-Marín et al. (2012), and calculated using the more accurate Bihalofit
fitting formula (Takahashi et al., 2020). Also shown are the bispectra
calculated from tree-level perturbation theory (Bernardeau et al., 2002).
Left: Equilateral triangles, Right: Isosceles triangles with sides in the
ratio 2:2:1.

##  Appendix B Signal-to-noise ratios

We approximate the cumulative signal-to-noise ratio (SNR) of a power spectrum
$P(k)$ as

$\displaystyle\Bigg{(}\frac{\mathrm{S}}{\mathrm{N}}\Bigg{)}^{2}$
$\displaystyle=\sum_{k_{i}=k_{\mathrm{min}}}^{k_{\mathrm{max}}}\frac{[P(k_{i})]^{2}}{\sigma^{2}_{i}}\
,$ (41)

where $\sigma^{2}_{i}$ is the variance in the $i$th bin, measured using
jackknife sampling, and we consider only diagonal terms of the covariance.
Similar definitions apply to the IA power spectra and bispectra. Equation (41)
is a simplification which somewhat overestimates the signal-to-noise ratio.
Nevertheless it allows a useful comparison between the information content of
different spectra.

Figure 11 shows, as a function of $k_{\mathrm{max}}$, the cumulative SNR for
the matter power spectrum and also for E-mode IA power spectra. Figure 12
shows similar information for the E-mode intrinsic alignment bispectra, for
equilateral triangles and for isosceles triangles with sides in the ratio
2:2:1.

The main messages to take from these figures are that the SNR for
$P_{\mathrm{E}\mathrm{E}}$ is much weaker than for $P_{\delta\mathrm{E}}$; the
power spectra contain much more information than the bispectra; and the
intrinsic alignment spectra $B_{\delta\mathrm{E}\mathrm{E}}$ and
$B_{\mathrm{E}\mathrm{E}\mathrm{E}}$ contain very little signal in the
IllustrisTNG300-1 volume. For this reason we use only
$B_{\delta\delta\mathrm{E}}$ in most parts of this work. The
$P_{\delta\mathrm{E}}$ SNR appears surprisingly strong in comparison to that
for $P_{\delta\delta}$ but this could be due to the simplifications in our
calculations. We also confirmed that in each case the cumulative SNRs for
individual mass bins are typically within 80 per centof the SNRs for the whole
sample.

Figure 11:  Cumulative signal-to-noise ratios for the matter power spectrum
and E-mode intrinsic alignment power spectra measured from simulations.
Figure 12:  Cumulative signal-to-noise ratios for the matter bispectrum and
E-mode intrinsic alignment bispectra measured from simulations. Top:
Equilateral triangles. Bottom: Isosceles triangles with sides in ratio 2:2:1.

##  Appendix C Phenomenological adjustment to analytical models

In Fig. 9 we compared IA isosceles bispectra measurements from simulations
with analytical predictions and showed that there is a close but not perfect
match between the two estimates. Figures 13 and 14 respectively show analogous
results for intrinsic alignment power spectra and equilateral bispectra. Again
the fits are close but not exact.

Figure 13:  Intrinsic alignment power spectra $P_{\delta\mathrm{E}}$ and
$P_{\mathrm{E}\mathrm{E}}$ estimated from simulations (solid lines) and
calculated from the non-linear alignment model with bestt-fitting values of
$A_{\mathrm{IA}}$ estimated from $P_{\delta\mathrm{E}}$ (dashed lines).
Figure 14:  Intrinsic alignment bispectra for equilateral triangles. Solid
lines: Measured from simulations. Dashed lines: Calculated from Eqs. (29) to
(31) with best-fitting values of $A_{\mathrm{IA}}$ estimated from
$B^{\mathrm{equi}}_{\delta\delta\mathrm{E}}$.

Looking more closely at Fig. 13, the simulations and analytical calculations
match well for $P_{\mathrm{E}\mathrm{E}}$; the discrepancy is only in
$P_{\delta\mathrm{E}}$ which correlates the matter overdensity and E-mode
fields (the latter is not a true field but we treat it as such). From Fig. 9
and Fig. 14 we see that the fit is not perfect for any of the bispectra. The
measured $B_{\mathrm{E}\mathrm{E}\mathrm{E}}$ (isosceles and equilateral) are
noisy and in this case part of the discrepancy may be due simply to the
uncertainty in the measurements. We therefore consider here how to improve the
fits for the spectra which depend on both the matter density and E-mode
fields. We tackle this by considering the correlation between the measured
power spectra of the two fields. The correlation coefficient $r(k)$ between
the power spectra can be defined as (Kurita et al., 2021)

$\displaystyle r^{2}(k)$
$\displaystyle=\frac{6}{5}\frac{P_{\delta\mathrm{E}}(k)^{2}}{P_{\delta\delta}(k)P_{\mathrm{E}\mathrm{E}}(k)}\
.$ (42)

The normalising factor 6/5 arises from integration over $\mu$ of the terms in
$(1-\mu^{2})$ in Eqs. (26) and (27):

$\displaystyle\frac{\big{[}\int_{0}^{1}(1-\mu^{2})\,\mathrm{d}\mu\big{]}^{2}}{\int_{0}^{1}(1-\mu^{2})^{2}\,\mathrm{d}\mu}=\frac{5}{6}\
.$ (43)

With this normalisation, $r(k)$ should be equal to unity if the non-linear
alignment model is valid. Our calculated correlation coefficients at three
redshifts are shown in blue in Fig. 15. The coefficients asymptote to
approximately one as $k\to 0$ but reduce to a lower, redshift-dependent value
for $k\gtrsim 1$, apart from $z=0$ where the correlation coefficient is very
noisy at large values of $k$.

To model the correlation coefficients at each redshift we fit a generalised
logistic curve of the form

$\displaystyle r(k)$
$\displaystyle=r_{0}+\frac{r_{1}-r_{0}}{\big{[}1+a\exp{(-b(k-k_{0})\big{]}}^{1/\nu}}\
.$ (44)

This has six free parameters: $r_{0}$ and $r_{1}$ are the lower and upper
asymptotes of the curve, $k_{0}$ determines the $k$ value at which the curve
starts to decline, and the remaining parameters, $a$, $b$ and $\nu$, control
the shape of the curve. In particular, $b$ controls the rate at which the
curve decreases as $k$ increases. We normalise the fitted curves to be equal
to one at $k_{\mathrm{min}}$. For $z=0$, because of the noisy data at larger
$k$, we set $r(k)=r(1.0)$ for $k\geq 1.0\,h\mathrm{Mpc}^{-1}$. The fitted
curves are shown in green in Fig. 15.

It is difficult to fit a single model which takes redshift into account since
we only have data for three redshifts. Instead we fit the model separately for
each redshift, leading to the best-fitting parameter values shown in Table 2.
The values for $z=0.5$ and $z=1$ are similar so it is plausible that single
model could in fact be constructed for both these redshifts.

Table 2: Best-fitting parameter values and normalisation factor for the modelled correlation coefficients given by Eq. (44). Redshift | $r_{0}$ | $r_{1}$ | $a$ | $b$ | $\nu$ | $k_{0}$ | Norm
---|---|---|---|---|---|---|---
0.0 | 0.6 | 0.95 | 1.03 | 614 | 57.8 | 0.19 | 0.95
0.5 | 0.43 | 1.02 | 1.72 | 12.5 | 6.00 | 0.16 | 1.00
1.0 | 0.43 | 1.06 | 2.20 | 19.9 | 13.0 | 0.08 | 1.04
Figure 15: Blue: Correlation coefficients $r(k)$ between the E-mode field and
the matter density field as defined by Eq. (42), for three redshifts. Green:
Logistic curve fitted to the calculated correlation coefficients using Eq.
(44) with parameter values from Table 2.

We now introduce the correlation coefficient into the two-point analytical
model by rewriting Eq. (26) as

$\displaystyle P_{\delta\mathrm{E}}(k)$
$\displaystyle=r(k)f_{\mathrm{IA}}P_{\mathrm{NL}}(k)\ ,$ (45)

for each redshift. Thus the model now includes the simulation-based
correlation between the fields. Similarly we incorporate correlation
coefficients into the three-point model given by Eqs. (29) and (30). So, since
$P_{\mathrm{EE}}$ is unchanged, we have

$\displaystyle B_{\delta\delta\mathrm{E}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=2\,\Big{[}f_{\mathrm{IA}}^{2}F_{2}^{\mathrm{eff}}(\bm{k}_{1},\bm{k}_{2})P_{\mathrm{NL}}(k_{1})P_{\mathrm{NL}}(k_{2})$
$\displaystyle\hskip
11.38092pt+r(k_{2})f_{\mathrm{IA}}F_{2}^{\mathrm{eff}}(\bm{k}_{2},\bm{k}_{3})P_{\mathrm{NL}}(k_{2})P_{\mathrm{NL}}(k_{3})$
$\displaystyle\hskip
11.38092pt+r(k_{3})f_{\mathrm{IA}}F_{2}^{\mathrm{eff}}(\bm{k}_{3},\bm{k}_{1})P_{\mathrm{NL}}(k_{3})P_{\mathrm{NL}}(k_{1})\Big{]}\
,$ (46)

and similarly for $B_{\delta\mathrm{EE}}$.

Figure 16 compares the unadjusted and adjusted analytical results to the
simulation measurements for $P_{\delta\mathrm{E}}$. In this case, by
construction, the phenomenological modification virtually eliminates the
discrepancy between the simulations and analytical calculations.

Figure 17 shows similar results for the bispectra $B_{\delta\delta\mathrm{E}}$
and $B_{\delta\mathrm{E}\mathrm{E}}$ for equilateral and isosceles triangles.
Here the benefits of introducing the correlation coefficient are less clear,
especially for $B_{\delta\mathrm{E}\mathrm{E}}$ where neither version of the
analytical model fits the simulations very well, especially at higher
redshifts. Nevertheless we conclude that it may be worth considering the
decorrelation between the density field and the intrinsic shape in future two-
point and three-point IA models.

Figure 16:  Intrinsic alignment power spectrum $P_{\delta\mathrm{E}}$ measured
from simulations and calculated from the non-linear alignment model with and
without adjustment by the correlation coefficient $r(k)$ given by Eq. (42).

Figure 17:  Intrinsic alignment bispectra $B_{\delta\delta\mathrm{E}}$ and
$B_{\delta\mathrm{E}\mathrm{E}}$ estimated from simulations and calculated
analytically with and without adjustment by the correlation coefficient $r(k)$
given by Eq. (42). Top: Equilateral triangles, Bottom: Isosceles triangles
with sides in ratio 2:2:1.