# Intrinsic alignments of bulges and discs

Yesukhei Jagvaral1, Sukhdeep Singh1, Rachel Mandelbaum1
1Department of Physics, McWilliams Center for Cosmology, Carnegie Mellon
University, Pittsburgh, PA 15213, USA E-mail<EMAIL_ADDRESS>https://orcid.org/0000-0001-7068-7037 0000-0003-2271-1527

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Galaxies exhibit coherent alignments with local structure in the Universe.
This effect, called Intrinsic Alignments (IA), is an important contributor to
the systematic uncertainties for wide-field weak lensing surveys. On
cosmological distance scales, intrinsic shape alignments have been observed in
red galaxies, which are usually bulge-dominated; while blue galaxies, which
are mostly disc-dominated, exhibit shape alignments consistent with a null
detection. However, disc-dominated galaxies typically consist of two prominent
structures: disc and bulge. Since the bulge component has similar properties
as elliptical galaxies and is thought to have formed in a similar fashion,
naturally one could ask whether the bulge components exhibit similar
alignments as ellipticals? In this paper, we investigate how different
components of galaxies exhibit IA in the TNG100-1 cosmological hydrodynamical
simulation, as well as the dependence of IA on the fraction of stars in
rotation-dominated structures at $z=0$. The measurements were controlled for
mass differences between the samples. We find that the bulges exhibit
significantly higher IA signals, with a nonlinear alignment model amplitude of
$A_{I}=2.98^{+0.36}_{-0.37}$ compared to the amplitude for the galaxies as a
whole (both components), $A_{I}=1.13^{+0.37}_{-0.35}$. The results for bulges
are statistically consistent with those for elliptical galaxies, which have
$A_{I}=3.47^{+0.57}_{-0.57}$. These results highlight the importance of
studying galaxy dynamics in order to understand galaxy alignments and their
cosmological implications.

###### keywords:

methods: numerical – cosmology: theory – galaxies: statistics – galaxies:
kinematics and dynamics – galaxies: structure – gravitational lensing: weak

††pubyear: 2021††pagerange: Intrinsic alignments of bulges and discs–Intrinsic
alignments of bulges and discs

##  1 Introduction

Deflection of light by matter inhomogeneities between the emitter and the
observer produces an effect known as gravitational lensing which results in
distorted images of the light source (Dodelson, 2017). The lensing is dubbed
strong when multiple images of the same source, or in rare cases Einstein
rings, are produced. However, the images of most objects in the Universe
experience minute distortions en route to the observer in a phenomenon known
as weak gravitational lensing. The statistical correlations of the lensed
shapes produced by coherent structure along the line of sight, known as cosmic
shear, is used as a direct probe of large-scale structure, and provides
cosmological information on dark energy (Kilbinger, 2015). Upcoming weak
lensing surveys such as the Vera C. Rubin Observatory Legacy Survey of Space
and Time (LSST)111 https://www.lsst.org/ , Roman Space Telescope222
https://roman.gsfc.nasa.gov/ High Latitude Survey (HLS) and Euclid333
https://www.euclid-ec.org/ should provide unprecedented precision and
constrain cosmological parameters to the percent level (Weinberg et al.,
2013).

However, weak lensing measurements suffer from a number of systematic biases
that must be corrected for precisely in order to unlock this cosmological
constraining power. One important contributor is intrinsic alignments (IA),
the tendency of galaxies to coherently align with the density field and
produce correlations among galaxy shapes (see, e.g., Troxel & Ishak, 2015, for
a review). This effect can masquerade as a weak lensing signal and bias the
constraints on cosmological parameters such as the equation of state of dark
energy. Thus, for precision cosmological measurements this effect has to be
taken into account either by mitigation or marginalization (Mandelbaum, 2018).
On the flip side, IA are of interest because they can give us insight into
galaxy formation and evolution, and into the large-scale structure, since IA
is produced by the gravitational interplay between galaxies and the underlying
large-scale structure.

Early analytic models for galaxy intrinsic alignments included the linear
alignment model (Catelan et al., 2001; Hirata & Seljak, 2004). Later
extensions included non-linear contributions (Bridle & King, 2007; Blazek et
al., 2019). These models can effectively describe the large-scale alignment
behavior of elliptical galaxies; however for intermediate and small scales,
these models tend to underestimate the alignment. In order to address this
problem, halo models of IA (Schneider & Bridle, 2010; Fortuna et al., 2021a)
were developed to describe the small scale effects. These analytic models
usually depend on assumptions relating the galaxy and the host dark matter
halo orientations which may or may not be valid.

On another frontier, N-body cosmological simulations containing only dark
matter (DM) have also been used to study IA by “painting” galaxies onto the DM
haloes. However, the predicted IA signals proved to be entirely dependent on
the model used (Heymans et al., 2006). Recently, with the advent of large-
volume hydrodynamic simulations, direct study of large-scale IA within
simulations has been enabled. There have been numerous studies of IA within
the Illustris, Illustris-TNG, MassiveBlack-II, Horizon-AGN and EAGLE
simulations (e.g., Tenneti et al., 2014; Velliscig et al., 2015b; Tenneti et
al., 2016; Chisari et al., 2015; Samuroff et al., 2020) with a goal of
developing accurate models to incorporate into future survey analysis
pipelines and to provide insight into how galaxy formation/evolution processes
influence the alignment signals. While such IA measurements will still depend
on the subgrid galaxy formation and evolution physics models used in the
hydrodynamical simulations, compared with ‘galaxy painting’ models, hydro
simulations are more complex and try to capture more aspects of the baryonic
physics.

Not all galaxies exhibit similar IA; many studies show that large-scale IA
depend on the galaxy location within its dark matter halo (satellite vs
central), color (red vs blue), and luminosity (Kirk et al., 2015). In
observational studies of large-scale IA, galaxy color typically serves as a
proxy for morphology and the samples are split by color (Kiessling et al.,
2015). Red galaxies (which are usually elliptical) are dynamically dispersion-
dominated and are believed to exhibit coherent aligned with the large-scale
tidal field. This model would explain the measured IA effects for red galaxies
(Mandelbaum et al., 2006; Hirata et al., 2007; Okumura et al., 2009; Joachimi
et al., 2011; Singh et al., 2015). In contrast, the dynamics of blue galaxies
(which are usually disc-dominated) are determined by their angular momentum,
and are influenced by the small scale tidal fields via torque, leading to a
correlation in their orientations (Hirata & Seljak, 2004). As a result, the
common understanding is that blue galaxies tend to have spin alignments with
the cosmic web (Dubois et al., 2014). In observational studies, large-scale IA
of galaxy shapes has not been detected for blue galaxies so far (Mandelbaum et
al., 2011; Samuroff et al., 2019; Johnston et al., 2019).

In simulations, all of the phase space data are available, making dynamical
morphological classification of galaxies possible. Tenneti et al. (2016)
studied dynamically classified disc and elliptical galaxies in MBII and
Illlustris, and found a similar trend that early-type galaxies are more
strongly aligned than late-types. Later, Chisari et al. (2015) came to the
conclusion that spheroidal galaxies tend to be aligned radially towards over-
densities and other spheroidals, whereas disc galaxies tend to be oriented
tangentially around spheroidals in 3D for the Horizon AGN simulation.
Additionally, Shao et al. (2016) studied IA in EAGLE and cosmo-OWLS
simulations and reported that dispersion-dominated (spheroidal) centrals have
a stronger alignment than rotation-dominated (disc) centrals with both their
dark matter halo and the distribution of their satellites.

Most disc galaxies consist of two prominent structures: bulge and disc.
Naively, one may assume that bulges are scaled down ellipticals. However,
there are some differences. Galactic bulges are further broken down into
classical bulges and disky-bulges (pseudobulges). Ellipticals and classical
bulges are hypothesized to be the products of major galaxy mergers. On the
other hand, pseudobulges are hypothesized to be the product of secular
evolution of galactic discs (Laurikainen et al., 2016). Also, Gadotti (2009)
concludes that classical bulges and ellipticals follow offset mass–size
relations, suggesting that high-mass bulges cannot be considered as high-mass
ellipticals that happen to be surrounded by a disc. Despite the differences
between ellipticals and bulges in general, they have similar stellar
kinematics and internal properties such as low gas content, low star formation
rate and red color. Since bulges are thought to form and evolve in a similar
fashion to elliptical galaxies and share the same properties, we investigate
how different components of galaxies exhibit IA in a hydrodynamical
simulation. Although disc galaxy samples have tended to produce null large-
scale IA signals (with substantial errors in some cases), it does not mean
future surveys can safely ignore them. This paper aims to provide a better
understanding of the morphological dependence of IA.

Also, there is some evidence that a simple split of red and blue galaxies is
not the full story. For example, Johnston et al. (2019) reports that the
modelling of red/blue-split galaxy alignments may be insufficient to describe
samples with variable central/satellite galaxy fractions; Georgiou et al.
(2019) came to a similar conclusion. Furthermore, S0 galaxies, which appear
red in color, usually have disc-like morphology (van der Kruit & Freeman,
2011); also, there are low mass blue elliptical galaxies. Both of these
factors complicate the color-morphology relation. IA is a dynamical effect and
has no direct physical relationship with galaxy color and chemical
composition, with any correlations between IA and color presumably being
inherited via the pronounced color-morphology correlation. Therefore, our
study focusing on IAs of dynamically classified galaxies will provide better
insight. Additionally, since these intrinsic alignment mechanisms are
inevitably tied with the environment and the formation of galaxies in their
host halo, IA also provides valuable insight into the large-scale structure
and galaxy formation/evolution.

In this work, we do not consider additional structures such as bars, rings,
stellar halos and spiral arms. We use projected shape correlations of galaxies
that are directly comparable and relatable to observed galaxy shape
correlations, instead of IA estimators that rely on a knowledge of the 3D
orientation of the galaxies (e.g., Tenneti et al., 2014; Velliscig et al.,
2015b).

We begin in § 2 by describing the theory and the IA model we have used. In § 3
we describe the simulation suite and the analysis methods that we use to
quantify IA, including galaxy shape estimators, two-point statistics and
misalignment angles. Next, in § 4 we present the measured shapes and discuss
the various misalignment angles we have investigated, and interpret the
measured two-point correlation functions. Finally, in § 5 we discuss and
conclude our paper.

##  2 Theory

Weak lensing measurements typically utilize the correlation function of the
ellipticity measurements of many galaxies in order to achieve the statistical
power to probe the lensing-induced galaxy shape correlations, compared to the
more dominant shape noise (Dodelson, 2017). The galaxy shears can be modeled
as a random term plus two coherent shears. The two coherent shears can be
written as $\gamma=\gamma^{G}+\gamma^{I}$, where $\gamma^{G}$ is due to
gravitational lensing and $\gamma^{I}$ is due to intrinsic alignments. Then,
the two-point correlation function is

$\langle\gamma\gamma\rangle=\langle\gamma^{G}\gamma^{G}\rangle+\langle\gamma^{G}\gamma^{I}\rangle+\langle\gamma^{I}\gamma^{G}\rangle+\langle\gamma^{I}\gamma^{I}\rangle$
(1)

The first term, $\langle\gamma^{G}\gamma^{G}\rangle$ is the desired signal in
weak lensing surveys. For completely random shape orientations in the absence
of lensing, the terms $\langle\gamma^{G}\gamma^{I}\rangle$ and
$\langle\gamma^{I}\gamma^{I}\rangle$, average to zero. In reality, due to
coherent alignments of galaxy shapes with the large-scale structure, these two
terms acquire non-zero values, thus contaminating the weak lensing signal.
Intuitively, one can think of the $\langle\gamma^{G}\gamma^{I}\rangle$ and the
$\langle\gamma^{I}\gamma^{G}\rangle$ terms as describing a pair of galaxies,
where one in the foreground (via its tidal field) lenses the one in the
background whilst the same tidal field (of the foreground galaxy) also affects
its own shape and that of nearby galaxies. In contrast, the
$\langle\gamma^{I}\gamma^{I}\rangle$ term reflects a scenario where two nearby
galaxies are affected by the same tidal field (Hirata & Seljak, 2004; Troxel &
Ishak, 2015; Singh & Mandelbaum, 2016).

###  2.1  The nonlinear alignment (NLA) model

In this paper, we will interpret the intrinsic alignment two-point correlation
functions measured in simulations in terms of the nonlinear alignment model.
This subsection describes how we calculate the theoretical predictions for
this model. We will start with the linear alignment (LA) model, and then the
non-linear alignment (NLA) model that incorporates the non-linear power
spectrum will be briefly discussed. The LA model assumes that IA are
determined by the tidal field at the time of galaxy formation (Catelan et al.,
2001). Consequently we can write the intrinsic shear in terms of the
primordial potential $\phi_{p}$:

$\gamma^{I}=(\gamma^{I}_{+},\gamma^{I}_{-})=-\frac{C_{1}}{4\pi
G}(\partial^{2}_{x}-\partial^{2}_{y},\partial_{x}\partial_{y})\phi_{p}$ (2)

where $C_{1}$ is the alignment amplitude and $(x,y)$ forms an arbitrary
coordinate system in the plane of the sky. Here our sign convention is the
same as in Singh et al. (2015). Next, assuming a linear galaxy bias between
matter overdensities and galaxy densities $\delta_{g}=b\delta_{m}$, the power
spectrum of galaxy-shape and shape-shape correlations are as follows (Hirata &
Seljak, 2004):

$\displaystyle P_{g+}(\bm{k})$
$\displaystyle=A_{I}b\frac{C_{1}\rho_{\text{crit}}\Omega_{m}}{D(z)}\frac{k_{x}^{2}-k_{y}^{2}}{k^{2}}P_{\delta}^{\text{lin}}(\bm{k})$
(3) $\displaystyle P_{++}(\bm{k})$
$\displaystyle=\left(A_{I}\frac{C_{1}\rho_{\text{crit}}\Omega_{m}}{D(z)}\frac{k_{x}^{2}-k_{y}^{2}}{k^{2}}\right)^{2}P_{\delta}^{\text{lin}}(\bm{k})$
(4) $\displaystyle P_{g\times}(\bm{k})$
$\displaystyle=A_{I}b\frac{C_{1}\rho_{\text{crit}}\Omega_{m}}{D(z)}\frac{k_{x}k_{y}}{k^{2}}P_{\delta}^{\text{lin}}(\bm{k}).$
(5)

Here $P_{\delta}^{\text{lin}}$ is the linear matter power spectrum, while
$P_{g+}$ is the cross-power spectrum between the galaxy density field and the
shear component along the line joining the galaxy pair. $P_{++}$ is the shape-
shape correlation of the shear component along the line joining the galaxy
pair. $P_{g\times}$ is the cross-power spectrum between the galaxy density
field and the shear component at 45∘ from the line joining the galaxy pair.
Following the convention in the literature (Joachimi et al., 2011; Singh et
al., 2015), we fix $C_{1}\rho_{\text{crit}}=0.0134$ and use the $A_{I}$
parameter to quantify the IA amplitude. The NLA model incorporates the non-
linear matter power spectrum to partially describe the non-linear regime
(Bridle & King, 2007), so the non-linear matter power spectrum
($P_{\delta}^{\text{nl}}(\bm{k})$) replaces the linear power spectrum
($P_{\delta}^{\text{lin}}(\bm{k})$) in Eqs. (3)–(5).

###  2.2 Modeling the real-space observables

To measure intrinsic alignments in data and simulations, we typically use
correlation functions, which can be obtained by fourier transforming the power
spectra as

$\xi_{ab}(\bm{r})=\int d^{3}\bm{k}P_{ab}(\bm{k})e^{i\bm{k}\cdot\bm{r}}$ (6)

where the indices $a,b\in(g,+)$ represent the two different fields being
correlated (in an auto correlation, $a=b$), and $r$ is the 3D separation for a
pair of galaxies. Two different samples of galaxies are cross-correlated in
order for one sample to trace the biased matter density ($g$) and the other to
trace the intrinsic shear ($+$). We further integrate the correlation
functions along the line-of-sight to obtain the 2D projected correlation
functions as (see e.g. Singh et al., 2015, for more details):

$w_{ab}(r_{\mathrm{p}})=\int_{-\Pi_{\mathrm{max}}}^{\Pi_{\mathrm{max}}}\mathrm{d}\Pi\,\xi_{ab}(r_{\mathrm{p}},\Pi).$
(7)

$\Pi$ is the line-of-sight separation and $r_{\mathrm{p}}$ is the projected
separation for a pair of galaxies, with $r^{2}=r_{p}^{2}+\Pi^{2}$.

For faster and numerically stable calculations of the theoretical model, we
follow the method of Singh et al. (2021), where the projected correlation
function from the previous equation is written as

$w_{ab}(r_{p})=\sum_{\ell=0}^{2}2\int_{0}^{\Pi_{\text{max}}}d\Pi\,\xi_{ab,2\ell}(r)L_{2\ell}\left(\frac{\Pi}{r}\right).$
(8)

where $L_{2\ell}$ are the Legendre polynomials of order $2\ell$, and
$\xi_{ab,2\ell}(r)$ are the multipoles of the correlation function obtained
from the power spectra as

$\xi_{ab,2\ell}(r)=(-1)^{\ell}\alpha_{2\ell}(\beta_{a},\beta_{b})\frac{1}{2\pi^{2}}\int
dkk^{2}P_{ab}(k)j_{2\ell}(kr).$ (9)

Here $j_{2\ell}$ denotes the spherical Bessel functions of order $2\ell$,
$P_{ab}$ refers to the cross power spectra as defined in Eq. (5) and the
galaxy power spectrum for a sample with galaxy bias $b$ is
$P^{\mathrm{nl}}_{gg}=b^{2}P^{\mathrm{nl}}_{\delta}$. We use the FFTlog
implementation in the MCFIT package (Li et al., 2019) to obtain
$\xi_{ab,2\ell}$. The prefactors $\alpha_{2\ell}(\beta_{a},\beta_{b})$ are as
follows:

$\displaystyle\begin{split}\alpha_{0}(\beta_{a},\beta_{b})&=1+\frac{1}{3}(\beta_{a}+\beta_{b})+\frac{1}{5}\beta_{a}\beta_{b}\\\
\alpha_{2}(\beta_{a},\beta_{b})&=\frac{2}{3}(\beta_{a}+\beta_{b})+\frac{4}{7}\beta_{a}\beta_{b}\\\
\alpha_{4}(\beta_{a},\beta_{b})&=\frac{8}{35}\beta_{a}\beta_{b}\end{split}$
(10)

The $\beta$ factors model the redshift space anisotropy of the correlation
function ($\beta=f/b$ in redshift space distortion measurements with $f$ being
the rate of growth of structure. Since we use the simulated data in real space
rather than redshift space, $\beta_{g}=0$. The IA model, on the other hand,
has an additional anisotropy from the use of projected shapes (Singh &
Mandelbaum, 2016), such that $\beta_{+}=-1$.

These projected two point functions $w_{g+}$ and $w_{++}$ are widely used in
observational IA studies. Together with $w_{gg}$, these equations can be used
to model the IA in the simulations and derive IA model parameters that can be
directly compared with observational studies.

##  3 The Simulation and Analysis Methods

Here we describe the simulated data and explain the methods we use to measure
simulated galaxy and DM halo shapes, misalignment angles, and two point
statistics, and to decompose/classify galaxies dynamically.

###  3.1 Simulated data

Here, we succinctly introduce the IllustrisTNG simulation used in this work
(for more information, please refer to Nelson et al., 2018; Pillepich et al.,
2018b; Springel et al., 2018; Naiman et al., 2018; Marinacci et al., 2018;
Nelson et al., 2019). The IllustrisTNG100-1 is a large-volume hydrodynamical
simulation with a box side length of 75 Mpc/h. The simulation was run using
the moving-mesh code Arepo (Springel, 2010) and has $2\times 1820^{3}$
resolution elements with a gravitational softening length of 0.7 kpc/h for
dark matter and star particles, respectively. The dark matter particle mass is
$7.46\times 10^{6}M_{\odot}$ and star particle masses are variable. The model
in the simulations for galaxy formation and evolution includes radiative gas
cooling and heating; star formation in the ISM; stellar evolution with metal
enrichment from supernovae; stellar, AGN and blackhole feedback; formation and
accretion of supermassive blackholes (Pillepich et al., 2018a). The DM halos
within the simulation were cataloged using friends-of-friends (FoF) methods
(Davis et al., 1985), and the subhalos were cataloged using the SUBFIND
algorithm (Springel et al., 2001). The simulation suite includes 100 snapshots
at different redshifts; we use the latest snapshot at $z=0$ for our analysis.

###  3.2 Shapes of Halos and Galaxies

To measure the shapes of galaxies and DM halos we utilize the mass quadrupole
moments (often incorrectly referred to as the inertia tensor). We use three
different definitions of these moments – ${I}_{ij}$ (simple), $\hat{I}_{ij}$
(reduced) and $\widetilde{I}_{ij}$ (reduced iterative) – with the first two
defined as:

${I}_{ij}=\frac{\sum_{n}m_{n}{r_{ni}r_{nj}}}{\sum_{n}m_{n}}$ (11)
$\hat{I}_{ij}=\frac{\sum_{n}\frac{m_{n}}{r_{n}^{2}}r_{ni}r_{nj}}{\sum_{n}\frac{m_{n}}{r_{n}^{2}}}.$
(12)

Here the summation index $n$ runs over all particles of a given type in a
given galaxy, where $m_{n}$ is the mass of the $n^{\rm th}$ particle and

$r_{n}^{2}=\sum_{i=1}^{3}r_{n,i}^{2}=x_{n}^{2}+y_{n}^{2}+z_{n}^{2}$ (13)

is the distance between the centre of mass and the $n^{\rm th}$ particle, with
$i$ indexing the three spatial directions.

The reduced mass quadrupole moment in Eq. (12) upweights particles that are
closer to the centre of the halo/galaxy, thus downweighting the loosely bound
particles in the outer regions. However, this specific moment definition
imposes a spherical symmetry on the halo/galaxy, and produces rounded shapes
compared to the simple mass quadrupole moment and compared to the true shape
of the halo/glaxy.

To devise an alternate estimator that upweights particles at small separations
without producing a bias, an iterative procedure is used to obtain
$\widetilde{I}_{ij}$: after the initial calculation of the $\hat{I}$, the
particles are rotated so that the three unit eigenvectors of the mass
quadrupole moment, defined as
$\mathbf{s}_{\mu}=\\{s_{x,\mu},s_{y,\mu},s_{z,\mu}\\}^{\tau}$ and
$\mu\in\\{a,b,c\\}$, are aligned with the x-, y-, z-axes, respectively. The
half-lengths of the principal axes of the ellipsoid are given by
$a\propto\sqrt{\omega_{a}}$, $b\propto\sqrt{\omega_{b}}$, and
$c\propto\sqrt{\omega_{c}}$, such that $a\geq b\geq c$ and
$\omega_{a},\omega_{b},\omega_{c}$ are the eigenvalues of the mass quadrupole
moment. Then, the radial distance, $r_{n}^{2}$ in Eq. (13), is replaced with
the elliptical radial distance:

$\tilde{r}_{n}^{2}=\left(\frac{x_{n}}{a}\right)^{2}+\left(\frac{y_{n}}{b}\right)^{2}+\left(\frac{z_{n}}{c}\right)^{2}.$
(14)

Then, $\widetilde{I}_{ij}$ is calculated using $\tilde{r}_{n}^{2}$ as:

$\widetilde{I}_{ij}=\frac{\sum_{n}\frac{m_{n}}{\tilde{r}_{n}^{2}}r_{n,i}r_{n,j}}{\sum_{n}\frac{m_{n}}{\tilde{r}_{n}^{2}}}.$
(15)

The iteration continues until the change in eigenvalues to the next step is
less than 1 percent (Warren et al., 1992).

To predict the projected alignment signals, we need to use the 3D mass
quadrupole moments to define 2D projected shapes. Following Joachimi et al.
(2013), we can obtain the projected 2D ellipse by solving
$\bm{x}^{\tau}{\mathbf{W}}^{-1}\bm{x}=1$, where

${\mathbf{W}}^{-1}=\sum_{\mu=1}^{3}\frac{\bm{s}_{\perp,\mu}\mathbf{\bm{missing}}{s}_{\perp,\mu}^{\tau}}{\omega_{\mu}^{2}}-\frac{\bm{k}\bm{k}^{\tau}}{\alpha^{2}}\;,$
(16)

and

$\mathbf{k}=\sum_{\mu=1}^{3}\frac{s_{\parallel,\mu}\bm{s}_{\perp,\mu}}{\omega_{\mu}^{2}}\;~{}~{}~{}\mbox{and}~{}~{}\alpha^{2}=\sum_{\mu=1}^{3}\left(\frac{s_{\parallel,\mu}}{\omega_{\mu}}\right)^{2}\;.$
(17)

Here, $\mathbf{s}_{\perp,\mu}=\\{s_{x,\mu},s_{y,\mu}\\}^{\tau}$ are the
eigenvectors projected along the projection axis (for which we arbitrarily
choose the z-axis of the 3D simulation box).

Then, the two components of the galaxy ellipticity can be expressed in terms
of the symmetric tensor ${\mathbf{W}}$

$(e_{1},e_{2})=\frac{(W_{xx}-W_{yy},2W_{xy})}{W_{xx}+W_{yy}+2\sqrt{\mathrm{det}\mathbf{W}}}.$
(18)

For the special case that the $\bm{s}_{c}$ lies perfectly along the projection
axis, the absolute value of the ellipticity is $|e|=(a-b)/(a+b)$. In terms of
the projected simulation box, the $x,y$ directions correspond to the positive
and negative direction of $e_{1}$ (since we projected along the $z$
direction).

###  3.3 Misalignment angle

DM halos are usually modeled by assuming the collisionless cold dark matter
particles gravitationally collapse in an ellipsoidal shape (Sheth et al.,
2001). However, the galaxies themselves exhibit a variety of shapes ranging
from very oblate thin discs to ellipticals. This shape variety among galaxy
populations is hypothesized to be a consequence of angular momentum
redistribution during galaxy formation and evolution (van der Kruit & Freeman,
2011). Therefore, the study of how galaxy shapes and angular momenta align
with those of their host DM halo will provide insight into the astrophysical
processes that guide galaxy evolution. The orientation of a galaxy relative to
its host DM halo can be described by three misalignment angles defined as

$\theta_{\mu\lambda}=\cos^{-1}(\bm{s}^{\mathrm{dm}}_{\mu}\cdot\bm{s}^{\mathrm{g}}_{\lambda})$
(19)

where $\lambda,\mu=a,b,c$, and
$\bm{s}^{\mathrm{dm}}_{\mu}(\bm{s}^{\mathrm{g}}_{\mu})$ is the unit vector
defining one of the principal axes of the DM halo (galaxy). In addition, we
calculate the misalignment angle of the minor axis
$\mathbf{s}^{\mathrm{g}}_{c}$ with the angular momentum of the DM halo and of
the galaxy:

$\theta_{\bm{L}c}=\cos^{-1}(\bm{L}^{\mathrm{k}}\cdot\bm{s}^{\mathrm{g}}_{c})$
(20)

where

$\bm{L}^{\mathrm{k}}=\sum_{n=1}^{N}m^{\mathrm{k}}_{n}\bm{x}_{n}\crossproduct\bm{v}_{n}.$
(21)

Here $\mathrm{k}\in\text{(dm, gal)}$, $m^{\mathrm{k}}_{n}$ is the mass of
particle $n$, $\bm{x}_{n}$ is its position, and $\bm{v}_{n}$ is its velocity
relative to the origin of the subhalo/galaxy. The origin of the
subhalo/galaxyis taken as the position of the most bound particle in the
subhalo/galaxy.

###  3.4 Two-point correlation function estimators

We measure the galaxy two-point correlation functions using the standard
Landy-Szalay estimator (Landy & Szalay, 1993),

$\xi_{gg}(r_{\rm p},\Pi)=\frac{DD-DR-DR+RR}{RR},$ (22)

where $DD$, $RR$ and $DR$ are weighted counts of galaxy-galaxy, random-random
and galaxy-random pairs, binned based on their perpendicular and line-of-sight
separation, $r_{p}$ and $\Pi$ 444All of the two-point statistic were measured
using the HALOTOOLS package v0.7 (Hearin et al., 2017) and the supporting
halotools_ia package ..

The cross correlation function of galaxy positions and intrinsic
ellipticities, $\xi_{g+}(r_{\mathrm{p}},\Pi)$, can be similarly measured using
a modified Landy-Szalay estimator (Mandelbaum et al., 2011) as a function of
$r_{\mathrm{p}}$ and $\Pi$:

$\xi_{g+}(r_{\mathrm{p}},\Pi)=\frac{S_{+}D-S_{+}R}{RR}.$ (23)

Similarly, the shape-shape correlation can be estimated as

$\xi_{++}(r_{\mathrm{p}},\Pi)=\frac{S_{+}S_{+}}{RR}.$ (24)

Here

$S_{+}D\equiv\frac{1}{2}\sum_{\alpha\neq\beta}w_{\alpha}w_{\beta}\,e_{+}(\beta|\alpha),$
(25)
$S_{+}S_{+}\equiv\frac{1}{4}\sum_{\alpha\neq\beta}w_{\alpha}w_{\beta}\,e_{+}(\alpha|\beta)e_{+}(\beta|\alpha),$
(26)

represent the shape correlations, where $e_{+}(\beta|\alpha)$ is the $+$
component of the ellipticity of galaxy $\beta$ (from the shape sample)
measured relative to the direction of galaxy $\alpha$ (from the density tracer
sample) and $w_{\alpha}$ ($w_{\beta})$ are weights associated with galaxy
$\alpha$ ($\beta$).

Finally, the projected two-point correlation functions (Eq. 7) are
approximated as sums over the line-of-sight separation ($\Pi$) bins:

$w_{ab}(r_{\mathrm{p}})=\sum_{-\Pi_{\mathrm{max}}}^{\Pi_{\mathrm{max}}}\Delta\Pi\,\xi_{ab}(r_{\mathrm{p}},\Pi),$
(27)

where $a,b\in(g,+)$.

Additionally, in both the model and the measurement we employ a
$\Pi_{\text{max}}$ value of 20 Mpc, and the model has a $k_{\text{min}}$
cutoff of $\pi/L_{\text{box}}=0.4$ h/Mpc. These choices are made in order to
exclude the modes that are not captured by the finite size of the simulation
box.

We fit the three correlation functions ($w_{gg}$, $w_{g+}$, and $w_{++}$)
jointly to the NLA model, with the galaxy bias $b$ and the intrinsic alignment
amplitude AI as the fit parameters. The same set of density tracers, including
all galaxies above $\log_{10}(M_{*}/M_{\odot})=10$, were used for all
correlation functions and fits. For our fits, we have used an analytic
estimate of the covariance matrix derived from the dominant Gaussian
contribution; under this assumption, the covariance matrix includes a noise
term and cosmic variance (for more details see Singh et al. 2017; Samuroff et
al. 2020).

###  3.5 2D kinematic decomposition model

In this subsection, we describe the dynamical model that probabilistically
assigns each particle to either the bulge or the disc component of the galaxy
(for more details, see Jagvaral et al. 2021). The model identifies the two
galaxy components through two physically-motivated assumptions:

  1. 1.

Disc stars’ angular momentum is approximately aligned with the total angular
momentum of the galaxy, while the orientation of bulge stars’ angular momentum
is randomly distributed.

  2. 2.

Disc stars’ orbits are approximately circular, while bulge stars’ orbits are
elongated or circular.

To implement a dynamical decomposition based on those two principles, we
define the following two parameters of interest:

  * •

$j_{\text{r}}\equiv\frac{j_{\text{star}}}{j_{\text{circ}}(r)}$, where
$\bm{j}_{\text{star}}$ is the angular momentum of a single star particle and
$j_{\text{star}}$ is its magnitude;
$j_{\text{circ}}(r)=r\,v_{\text{circ}}(r)=r\,\sqrt{\frac{GM(r)}{r}}$ is the
expected angular momentum for a circular orbit at the same position as that
star, where $M(r)$ is the total mass (across all types of particles – stars,
gas, dark matter) contained within that radius. This parameter indicates the
type of orbit taken by the star particle. Stars on circular orbits will have
$j_{\text{r}}\sim 1$, while those on elliptical orbits can have $j_{\text{r}}$
either above or below 1.

  * •

$\cos\alpha$ is the cosine of the angle between the angular momentum vector of
the star particle and the total angular momentum of the galaxy. A
concentration of particles at $\cos\alpha\sim 1$ signals that the galaxy
contains a disc structure, since particles preferentially have their angular
momentum aligned with that of the galaxy overall. Spread within the two
angular parameters $\cos\alpha$ and $\cos\beta$ signifies the disordered
motion among bulge stars.

To build this 2D model, we assume bulge stars should exhibit a flat
distribution in $\cos(\alpha)$, so their distribution in this plane should be
solely dependent on $j_{r}$. On the other hand, we assume the distribution of
disc stars depends on both parameters. Hence, we consider the following model
for the probability distribution of star particles:

$p_{\text{star}}(j_{\text{r}},\cos\alpha)\equiv(1-f^{\text{disc}})\,p_{\text{bulge}}(j_{\text{r}},\cos\alpha)+f^{\text{disc}}\,p_{\text{disc}}(j_{\text{r}},\cos\alpha).$
(28)

Here $p_{\text{bulge}}$ and $p_{\text{disc}}$ are the probability
distributions (both normalized to integrate to 1) that a star at a given point
in this 2D space belongs to the bulge or the disc.

The 2D model is built deterministically, after which we can generate Monte
Carlo realizations of the model, assigning star particles to the bulge or to
the disc. In this work, we will focus on a single realization of the Monte
Carlo simulation; the probabilistic nature of the Monte Carlo method was
explored in Jagvaral et al. (2021).

###  3.6 Sample selection

We employ a minimum stellar mass threshold of $\log_{10}(M_{*}/M_{\odot})=10$
for all galaxies, using their stellar mass from the SUBFIND catalog. A
previous study (Jagvaral et al., 2021) motivates the mass cut on the shape
sample, where we concluded that the disc fractions from TNG100-1 agree well
with observational values down to $\log_{10}(M_{*}/M_{\odot})=10$, below which
we believe the resolution of the simulation may be affecting the results.
Also, we have checked that using a lower mass on the density tracers, to
include more galaxies, did not change our $A_{I}$ amplitudes.

Based on the model discussed in the previous section, we classify galaxies
into three bins in $f_{\mathrm{disc}}$, each containing the same number of
galaxies, where $f_{\mathrm{disc}}$ is the fraction of stellar mass in the
disc for a given galaxy. For the rest of this paper, we will refer to the bin
with the lowest mean $\langle f_{\mathrm{disc}}\rangle=0.37$ as the Elliptical
sample, the bin with the highest $\langle f_{\mathrm{disc}}\rangle=0.83$ as
the Pure Disc sample, and the middle bin with $\langle
f_{\mathrm{disc}}\rangle=0.68$ as the two-component Disc+Bulge sample.
Further, the calculations for the Disc+Bulge sample are carried out separately
for the two components of each galaxy: the Bulge (only) sample and the Disc
(only) sample.

Often galaxies are split into two groups based on morphology, but in this
study we chose to split the sample into three bins in order to gain more
insight into how IA depends on $f_{\mathrm{disc}}$. We chose the sample names
for simplicity, but should emphasize that the galaxies in the Pure Disc sample
contain a non-negligible fraction of dispersion-dominated stars. Conversely,
the galaxies in the Elliptical sample contain a non-negligible fraction of
rotation-dominated stars. These non-negligible but small fractions do not lend
themselves to making reliable shape measurements, so we do not study them as
separate structures.

Since IA is known to depend on both the stellar and total mass (e.g., Tenneti
et al., 2015), we have controlled for differences in the mass distributions in
our samples by implementing a mass-dependent weight for each galaxy. We binned
the galaxies in each sample by their total subhalo masses. Then we took the
ratio of the histogram bin heights for the Pure disc sample to the histogram
bin heights for the given sample (because the Pure disc sample has the
narrowest distribution). These ratios were used as weights when calculating
the two-point functions using the weighted estimators in Sec. 3.4.

Also, IA is known to depend on the satellite fraction in a given sample
(Johnston et al., 2019). In our samples, the satellite fractions were within 5
per cent of each other,so we did not control for this; the satellite fractions
were 0.62, 0.63, 0.58 for the Elliptical, Disc+Bulge and Pure Disc samples,
respectively. Lastly, we have checked the total host (parent) halo mass
distributions of all samples after controlling for total subhalo mass and the
mass distribution of the host halos is very similar, suggesting that the
samples are experiencing very similar environments.

##  4 Results

In this section, we will first describe the measured shapes of the classified
galaxies and galaxy components. The misalignments of galaxy shapes and angular
momenta with respect to those of the host DM halos are examined next. Third,
we measure the intrinsic alignment signals of the various samples. For large-
scales we compare them with the predictions of the NLA model, while for small
scales we quantify the alignment amplitude by fitting a simple power law
model.

###  4.1 Shapes of galaxies and components

Our first test is to compare the galaxy shapes measured using the three
different mass quadrupole moment definitions described in Sec. 3.2. Fig. 1
shows the distribution of 3D shapes, defined as $s=c/a$ and $q=b/a$, obtained
from the three different methods of mass quadrupole moment: simple, reduced,
reduced iterative. We will use it to explore the shape distributions of
morphologically classified samples. The reduced method produces the most round
shapes, by imposing a spherically symmetric down-weighting of particles in the
galaxy outskirts, thus biasing the measurement. By comparing the axis ratios
for the reduced method (middle column) and reduced iterative methods (right
column), we can see that the latter are systematically less round, because the
iterative procedure with an elliptical weight function avoids the biases of
the reduced method. The simple method shows overlap among the morphologically
separated samples in the plots comparing axis ratios for elliptical, pure
disc, and disc$+$bulge samples, while the reduced and iterative reduced
methods show similar overlap between the Disc+Bulge and Elliptical samples. In
contrast, for the Disc (only) and the Bulge (only) samples the contours do not
overlap, implying these are morphologically distinct structures within a
common set of galaxies (at least when the outer particles are downweighted),
however, there is still some overlap for the simple method at $\langle
s\rangle\sim 0.5$. This suggests that the disc structures may be thicker than
expected from observations (Sandage et al., 1970), since for the simple method
particles are not weighted by their distance.

Interestingly, the shapes obtained by the reduced method correlate moderately
(Pearson-$r$ coefficient of +0.51) rather than strongly with the simple method
shapes, which may suggest that for a significant number of galaxies the outer
region is quite different than the inner region. The shapes from the reduced
and the iterative reduced methods are strongly correlated, with a Pearson-$r$
coefficient of 0.9. Even when corrected for the rounding of the reduced method
via iteration, the shape measurements still show moderate correlation with the
simple one. Therefore, we will ignore the reduced method and focus on the
shapes obtained by the simple and the reduced iterative methods.

Next, in Fig. 2 we investigate how the average axis ratio changes with mass.
All samples show relatively flat trends with mass. For the most part, when
comparing $\langle s\rangle$ and $\langle q\rangle$, $\langle s\rangle$
appears to be lower than $\langle q\rangle$ by 20-40 per cent for a given
sample, which shows that the galaxies in this simulation are generally oblate.
As expected for $\langle s\rangle$, the Disc (only) sample is the thinnest
(most oblate), followed by the Disc+Bulge sample. Ellipticals and Bulge (only)
are at a higher value, $\sim$0.7. Also, looking at plot showing dependence of
$\langle s\rangle$ on mass, we see that the Pure Disc and Disc (only) samples
are thinner by $\sim$ 20-25 percent compared to the Ellipticals and the Bulges
(only), which implies that the measured shapes of disc structures are
relatively thick.

Figure 1: Comparison of shapes obtained from the three different mass
quadrupole moment definitions: simple, reduced, reduced iterative ordered by
columns; and rows are split by different samples as labeled in the legend,
where the middle row shows three distinct samples and the bottom row shows the
results for the two components of the galaxies in the Disc+Bulge sample. The
simple method generally produces a broader range of shapes compared with the
other two methods, stretching from $s=0.2$–$0.9$ and $q=0.6$–$1$. The reduced
method produces the most round shapes, as illustrated by the fact that the
population in all samples is closer to the top-right edge of the plots than
for the other two methods. The reduced iterative method corrects the over-
rounding effect and produces a wider range of shapes than the reduced method.
Comparing the Disc+Bulge sample from row 2 to row 3 (where it is split into
components), the population separates out into a bimodal distribution.
Additionally, for row 2, we can see that the populations with higher disc
fraction are to the left, whereas populations with lower disc fraction are to
the right.  Figure 2: Comparison of average axis ratios as a function of
stellar mass. For $\langle s\rangle$, the Disc (only) populations as expected
are the lowest, followed by Disc+Bulge, and Bulge (only) and Ellipticals are
higher. For $\langle q\rangle$, all samples are relatively flat between 0.8
and 0.9.  Figure 3: Distributions of the projected ellipticity magnitude $|e|$
for different simulated samples and observed values from the HST COSMOS
samples from (Joachimi et al., 2013). Top panel: Ellipticity distribution
obtained using the simple mass quadrupole moment. Bottom panel: Ellipticity
distribution obtained using the reduced iterative mass quadrupole moment.
Comparing the different samples, the bulges exhibit the roundest ellipticity
distribution, with a peak around 0.09. Samples with significant discs exhibit
the broadest distribution. The Elliptical and the Bulge (only) samples have
similar distributions as the COSMOS Elliptical sample. However, samples with
significant discs do not agree with the COSMOS Disc sample. For the bottom
panel, Disc (only) and Disc+Bulge samples both show broad distributions that
are closer to the the COSMOS Disc sample than the Disc (only) and Disc+Bulge
samples from the top panel.

Next, we compare the 2D projected shapes, described at the end of Section 3.2,
with observed values as seen in Fig. 3. Here, we show the ellipticity
distributions of our samples and the HST COSMOS samples from Joachimi et al.
(2013). These 2D (projected) shape distributions follow the same trend as the
3D shape distributions where Elliptical and Bulge (only) samples display round
shapes and Bulge+Disc, Disc (only) and Pure Disc  samples display elongated
shapes. Interestingly, the projected shapes obtained by the simple method show
all of the samples being close to round, since the distributions have peaks at
around $|e|=0.1$ but with a tail to higher values. Also, the Elliptical and
the Disc+Bulge samples show very similar distribution, as was seen in the 3D
case. Additionally, when compared with the COSMOS data, the Elliptical and the
Bulge (only) samples have similar projected shape distributions as the COSMOS
Elliptical sample; on the contrary, the Disc-related samples do not agree with
COSMOS Disc sample, the latter of which has a very flat ellipticity
distribution.

The reduced iterative 2D shape distributions in the bottom panel of Fig. 3
exhibit more variety between samples compared to those for the simple method.
The Bulge (only) sample exhibits the roundest ellipticity distribution, with a
peak around 0.09 followed by a steep drop-off, followed by the Elliptical
sample. The Disc (only) and Disc+Bulge samples exhibit very similar, broad
distributions that extend to $|e|\sim 0.6$; these two distributions are closer
to the distribution for the COSMOS Disc sample than what was obtained using
the simple method. Additionally, the COSMOS Elliptical data is midway between
the TNG-100 Bulge (only) and the Elliptical samples. The disagreement between
the observed and simulated shapes of the disc samples may be due to the
simulation producing thick discs, as was noted in discussion of the 3D shape
distributions.

###  4.2 Misalignment of galaxies and components with their host DM halos

In order to characterize the relative orientation between the galaxy and the
DM halo we focus on the relative angles between their axes. For rotation-
dominated structures, the minor axis should strongly coincide with the total
angular momentum of the system. Since these objects are often circular in
shape in the plane of the galaxy, there is also some ambiguity in determining
their major and the intermediate axis, causing a degeneracy in these axes. In
Fig. 4 we show the distribution of various misalignment angles of different
galaxy types and components and in Table 1 we tabulate the median, the 16th
and the 84th percentiles for the top-left, middle-right, bottom-right panels
in Fig. 4. In all panels, the grey dotted lines denote the expected
distribution for purely random orientations; as shown, all galaxy types and
components have alignments that are not random with respect to the halo shape,
and the galaxy and halo angular momentum. In the plots showing the alignment
angles between the galaxy axis and halo axis in the top row, Ellipticals show
stronger alignment compared with the other samples, consistent with the
physical picture that elliptical galaxies tend to be stretched out along the
direction of the tidal fields. Interestingly, the Bulge (only) sample follows
a very similar trend as the Ellipticals suggesting that bulges also tend to be
stretched out by the tidal fields. Comparing the alignments of the
intermediate and major axes, we see that the alignments are less significant
for the intermediate axes for the Ellipticals and Bulges (only), however for
Pure Disc and Disc (only) samples the curves look very similar, since there is
some ambiguity in determining the major and intermediate axes for oblate
objects. For all samples, the minor axes of the galaxy and DM halo shape
exhibit the strongest alignment compared to the other two axes, again with
Ellipticals and Bulges (only) being the most strongly aligned (with median
$\theta=20.34,21.77$ respectively), followed by the disc-related structures.

Next, we inspect how the galaxy axis aligns with respect to the halo angular
momentum vector as shown in the second row of Fig. 4. As expected, the galaxy
major and intermediate axes exhibit a tendency to anti-align with the halo
angular momentum vector (i.e., systematically less aligned than a random
distribution) for all samples. The galaxy minor axis, as expected, has a
strong alignment with the halo angular momentum for the Pure Disc , Disc+Bulge
and Disc (only) cases (with median $\theta=15.99,21.13,18.97$, respectively) ,
suggesting a strong torquing mechanism between the DM halos and the resident
galaxies. In contrast, the Bulge (only) and Elliptical samples show weaker
alignments with the halo angular momentum (with median $\theta=36.59,28.87$,
respectively), implying a weaker torquing mechanism from the DM halo to the
resident galaxies.

Lastly, the bottom row of Fig. 4 shows that there is a strong alignment
between the minor axis of disc-like samples and the galaxy’s total angular
momentum, with Pure Disc , Disc+Bulge and Disc (only) showing low median
$\theta$ of 5.51, 9.29, 6.09, respectively . While the Bulge (only) and
Elliptical samples exhibit lower alignment with higher median $\theta$ of
29.41 and 17.16, respectively. Comparing how the galaxy shape axes are aligned
with the halo and galaxy angular momentum vectors, the galaxy shape minor axis
is more aligned, and the major and the intermediate shape axes are more anti-
aligned with the galaxy angular momentum vector. Also, for these angular
momentum vector alignments with the galaxy axes, the ambiguity in determining
the major and intermediate axes is not relevant, since we are primarily
interested in the alignment of the minor axes with the angular momentum
vectors.

Overall, all samples exhibit a preferred alignment of their galaxy and DM halo
shape axes, with the Ellipticals and Bulge (only) samples showing a stronger
alignment compared to the other samples. Also, all samples tend to align their
galaxy minor axes with the galaxy angular momentum and to a lesser degree with
the DM halo angular momentum, with the Pure Disc, Disc+Bulge and Disc (only)
samples exhibiting stronger alignment than the Elliptical and Bulge (only)
samples.

Further, to test the dependence of the alignments with the total angular
momentum vectors shown in Fig. 4 on the angular momentum magnitude, we explore
the relationship between the alignment angle $\theta_{\bm{L}c}$ with the
magnitude of the angular momentum vector $|\bm{L}|$ of the galaxy and the DM
halo in Fig. 5. When comparing different samples, the Disc+Bulge, Disc (only)
and Pure Disc samples show a high concentration of probability density near 0
degrees, indicating strong alignment, whereas the Elliptical and the Bulge
(only) samples have a smooth distribution with no highly concentrated areas.
All panels in this plot show that a low angular momentum magnitude $|\bm{L}|$
leads to significant misalignment, consistent with the physical picture that
significant angular momentum is required in order for it to drive galaxy
alignments.

Galaxy-halo axis misalignment has been studied in various cosmological
hydrodynamical simulations: in MassiveBlack-II and Illustris (Tenneti et al.,
2016), in Horizon-AGN (Chisari et al., 2017), and in EAGLE (Velliscig et al.,
2015a). Both Tenneti et al. (2016) and Chisari et al. (2017) found that
elliptical galaxies in the simulations have stronger alignments with their
host halos than do disc galaxies, similar to what we have found in TNG100. In
contrast, Velliscig et al. (2015a) reports that elliptical galaxies are more
misaligned with their host halos compared to disc galaxies. We note that each
of these studies used slightly different dynamical methods to identify the
galaxy morphology, and the studies did not control for intrinsic alignment
trends with mass, though they explored trends in mass by breaking their
samples into mass bins. Furthermore, this study and those cited above find
that galaxy shapes exhibit coherent alignments with respect to host halo
shapes, and that the elliptical versus disc galaxy distinction is necessary
for accurate models of misalignment angles. As we will see in Sec. 4.3, galaxy
samples that are more aligned with their host DM halos exhibit higher large-
scale galaxy alignment correlations, which we can quantify using the NLA model
amplitude.

Figure 4: Distribution of misalignment angles of various shape axes and total angular momenta. [G = galaxy, $\mathrm{\mu}$ = axis $\in(a,b,c)$, H = DM halo, L = total angular momentum, e.g. G$\mathrm{\mu}$-HL means misalignment angle between the galaxy axis and the DM halo angular momentum]. The grey dotted line in all panels represents random alignment. The minor axis of the galaxy for the most part aligns well with the minor axis of the DM halo and the two angular momentum vectors, especially for the Disc+Bulge and Disc (only) samples. Conversely, the major and intermediate axes are mostly perpendicular to the two $L$ vectors. The galaxy and DM halo major and intermediate axes shows sign of preferred alignment compared to a random distribution for all samples.  | Elliptical | Disc+Bulge | Bulge (only) | Disc (only) | Pure Disc
---|---|---|---|---|---
Gc-Hc (top-left in Fig. 4) | $20.34^{+49.96}_{-8.10}$ | $26.79^{+57.78}_{-9.99}$ | $21.77^{+50.81}_{-8.50}$ | $31.92^{+63.68}_{-12.36}$ | $32.97^{+66.11}_{-12.80}$
Gc-HL (middle-right in Fig. 4) | $28.87^{+64.01}_{-9.69}$ | $21.13^{+53.95}_{-7.27}$ | $36.59^{+70.68}_{-13.23}$ | $18.97^{+47.42}_{-6.77}$ | $15.99^{+39.39}_{-6.42}$
Gc-GL (bottom-right in Fig. 4) | $17.16^{+51.77}_{-5.74}$ | $9.29^{+32.04}_{-3.23}$ | $29.41^{+68.79}_{-9.94}$ | $6.09^{+18.83}_{-2.27}$ | $5.51^{+17.60}_{-1.96}$
Table 1:  Median and 16th and 84th percentiles of misalignment angles of
galaxy-halo major axes, and galaxy minor axes against total angular momenta of
galaxy and halo. This Table is intended to supplement Fig. 4 by providing
quantitative information on the plotted distributions.  Figure 5: Top row:
Distribution of the magnitude of the total angular momentum of galaxies vs.
the misalignment angle with the minor axis of the galaxy, where $\mathbf{L}$
is in units of $\mathrm{Mpc\,km}\sqrt{a}/s\,M_{\odot}$. Bottom row: The same,
for the magnitude of the total angular momentum of halos. The color bar
indicates the number of galaxies in each bin, while the red line indicates the
median for each panel at a given angle. When contrasting the top and the
bottom row, the galaxy angular momentum is more aligned with the minor axis of
the galaxy shape than the halo’s angular momentum for each sample. For the
bottom row, the Disc+Bulge, Disc (only) and the Pure Disc samples show high
concentration of population near 0 degrees, whereas the Elliptical and the
Bulge (only) samples show smooth distribution with no highly concentrated
areas. Overall, the decreasing trend indicated by the red line shows that the
alignment is not dominated by low magnitude angular momentum vectors.

###  4.3 Alignment measurements from two-point functions

In this section we present and discuss the results from the joint fits to the
NLA model that were performed using the $w_{gg}$, $w_{g+}$, $w_{++}$ data
vectors. We found that $w_{g+}$ for the reduced iterative method is consistent
with zero, similar to the findings of a previous study done by Chisari et al.
(2015) using the reduced method. This may be due to particles in the outer
regions of galaxies being affected more by the tidal forces from the large-
scale structure compared to the particles in the inner regions (Singh &
Mandelbaum, 2016; Chisari et al., 2017). Hence in this section we focus on the
results obtained using the simple mass quadrupole moment. Since the NLA model
is known to break down at small scales ($r_{p}$ less than a few Mpc), our fits
to the NLA model only include $r_{p}>4$ Mpc (Bridle & King, 2007).

The measured values and the best-fitting curves for $w_{g+}$, $w_{++}$ are
shown in Fig. 6. We omitted $w_{gg}$ for the density tracer sample since it
was the same in all the fits. In general, the NLA model fits the measured
values from the simulation well on large scales. As in previous studies
(Mandelbaum et al., 2006; Singh et al., 2015), the $w_{++}$ curves are too
noisy to carry out meaningful model fits and were consistent with zero; we
only include them in the joint fits for the NLA model at large-scales rather
than trying to model them separately. Also, from the extrapolation of the
best-fitting curves, one can see that the NLA model severely underestimates
$w_{g+}$ at small scales, confirming our previous assumption that the model
will fail on smaller scales. In all of the samples the galaxy bias parameter
from the joint fits was around 1; the bias parameter being the same across
samples is expected, since we have used the same density tracer sample in all
cases.

In row 1 of Table 2 we show the best fitting parameters of the NLA model for
our mass-controlled samples. These results are visualized in Fig. 7 with
purple points showing $A_{I}$ versus $f_{\text{disc}}$. As expected, the
Elliptical galaxies have the highest $A_{I}$ of $3.47^{+0.57}_{-0.57}$,
whereas the Pure Disc population has an $A_{I}$ value consistent with zero. In
addition, when the two-component Disc+Bulge galaxies are split into their
separate components, the Bulge (only) sample exhibits a high $A_{I}$ signal
comparable to that of Ellipticals galaxies, at $2.98^{+0.36}_{-0.37}$. In
contrast, the Disc (only) component has a lower $A_{I}$ value of
$1.02^{+0.34}_{-0.36}$, still different than that for the Pure Disc sample,
corresponding to a marginal detection $2.5\sigma$ from zero. The two-component
galaxies as a whole have a weak signal at $A_{I}=1.13^{+0.37}_{-0.35}$ that is
consistent with that for the Disc (only) sample; this is understandable since
the shape of those galaxies is dominated by the disc structure when using the
simple shape estimator. In Fig. 7 we see a clear downward trend of $A_{I}$
with $f_{\text{disc}}$, and with the Bulge (Only) sample showing consistent
alignment as the Elliptical sample.

Compared with previous IA studies using hydrodynamical simulations, our
results for the whole non-mass controlled sample agree with Samuroff et al.
(2019), who report $A_{I}=1.71^{+0.17}_{-0.17}$ and
$b_{g}=1.11^{+0.07}_{-0.07}$, though it should be noted that they used the
larger volume TNG300 and a lower mass-cut of
$\log_{10}(M_{*}/M_{\odot})\approx 9.37$ on their sample.

We have investigated how the intrinsic alignment results for mass-controlled
samples compare with the full (i.e. non mass-controlled samples) samples. We
did the same analysis for the full samples; row 2 of Table 2 shows the
resulting NLA fit parameters. As expected given the higher average mass, the
full samples shows higher $A_{I}$ value by $\sim$5-15% compared to the mass
controlled samples. Nonetheless, the trends with $f_{\text{disc}}$ are very
similar in both cases as illustrated by Fig. 7: both the mass-controlled and
full samples exhibit a decreasing trend in $A_{I}$ with increasing
$f_{\text{disc}}$.

In observational tests of intrinsic alignments, the galaxy samples are usually
split by color, since morphological information is typically unavailable and
color serves as a proxy for morphology. Hence, in this study we also measure
the alignments for samples split by color to provide better comparison with
observational studies. We split the full sample in the $g-r$ color versus $r$
magnitude plane using the same method as in Samuroff et al. (2020). The
results of NLA fits to the alignment signals are summarized in Table 3. For
the Red sample we present both the mass-controlled and the full sample;
however the Blue sample was not mass controlled because of its very narrow
mass range. Here, we should note that the simulation produces a slightly
higher number of red disc galaxies compared to observations, which complicates
the interpretation of these results (Nelson et al., 2018; Jagvaral et al.,
2021). The Blue sample has an IA signal that is consistent with zero, similar
to previous studies in real data such as Johnston et al. (2019); Mandelbaum et
al. (2011). The Red sample has alignment amplitude
$A_{I}=2.79^{+0.50}_{-0.35}$, consistent with the $A_{I}=3.63^{+0.79}_{-0.79}$
of the GAMA Z1R sample in (Johnston et al., 2019). The two samples have
comparable mean luminosities of $\langle L\rangle/L_{0}=0.45$ and $0.50$,
respectively, where $L_{0}$ is the pivot luminosity corresponding to
$M_{r}=-22$. One interesting feature of this result is that the
morphologically-selected samples, Ellipticals and Bulge (only), both exhibit a
higher $A_{I}$ signal compared to the color-selected Red sample. The mass-
controlled Elliptical sample shows almost thrice as strong IA signal compared
to the mass-controlled Red sample. Compared with observational results from
Fortuna et al. (2021b), their dense red sample with average luminosity of
$\langle L\rangle/L_{0}=0.38$ has the same average luminosity as our
Elliptical sample, and has $A_{I}=3.69^{+0.66}_{-0.65}$, which agrees with
both our mass-controlled and full elliptical sample. Also their L1 sample,
which is a subset of their luminous red galaxy sample, with $\langle
L\rangle/L_{0}=0.46$ and $A_{I}=1.80^{+0.96}_{-0.95}$, has consistent
alignment amplitude with both the mass-controlled and the full Red samples
within the error bar.

Figure 6: Measured quantities $w_{g+}$ (top row) and $w_{++}$ (bottom row) in blue for 6 different dynamically classified/decomposed mass-controlled samples as indicated in the column titles. Data on scales below $r_{p}<4$ Mpc (indicated by the orange dashed vertical line) were excluded from the NLA fits. For each sample (i.e., for the data in a given column), $w_{g+}$, $w_{++}$ along with density tracer $w_{gg}$ ($w_{gg}$ curve not shown here, since it was the same for all samples) were fit simultaneously to the NLA model using a Markov Chain Monte Carlo method in the range $r_{p}<4$ Mpc; the best-fitting parameters are shown in row 1 of Table 2. In order to quantify the small scale alignments, we fit a power-law function (Eq. 29) for $r_{p}<1$ Mpc (indicated by the red dashed vertical line). The power-law fit parameters are presented in Table 4. The best fitting curves for the NLA model are shown in orange and for the power-law model in red. We show the extrapolation of the best-fitting NLA model to small scales, to justify our choice of $r_{p}$ cut.  Figure 7: Dependence of the alignment strength parameter $A_{I}$ on $f_{\text{disc}}$. The purple points are the main mass-controlled samples presented in the top row of Table 2. The orange points are the non-mass-controlled samples (i.e., the full samples from the simulation, shown in the second row of Table 2). In both cases, the alignment strength exhibits a decreasing trend with $f_{\text{disc}}$. The purple points are slightly below the orange points, indicating that there is some dependence on mass (with higher mass resulting in greater alignments), however the overall trend with $f_{\text{disc}}$ is the same. Additionally we have plotted the color-split samples with the diamond-shaped points; for the Red sample, the half-filled diamond represents the mass-controlled sample and the full diamond represents the full sample. Purple points were arbitrarily shifted by 1% to the left for visual clarity.  | All | Elliptical | Disc+Bulge | Bulge (only) | Disc (only) | Pure Disc
---|---|---|---|---|---|---
mass controlled $A_{I}$ | $1.64^{+0.28}_{-0.28}$ | $3.47^{+0.57}_{-0.57}$ | $1.13^{+0.37}_{-0.35}$ | $2.98^{+0.36}_{-0.37}$ | $1.02^{+0.34}_{-0.36}$ | $-0.34^{+0.35}_{-0.36}$
non-mass controlled $A_{I}$ | $1.99^{+0.19}_{-0.19}$ | $4.03^{+0.36}_{-0.34}$ | $1.31^{+0.22}_{-0.20}$ | $3.24^{+0.22}_{-0.25}$ | $1.09^{+0.27}_{-0.29}$ | $-0.34^{+0.35}_{-0.36}$
Table 2: Best fitting parameters of the NLA model using $w_{gg}$, $w_{g+}$,
$w_{++}$ data points simultaneously, for $r_{p}>4$ Mpc. As expected, bulge-
dominated galaxies exhibit a higher $A_{I}$. In addition, when disc-dominated
galaxies are split by their components, Bulge (only) shows high $A_{I}$ signal
comparable to that of Elliptical galaxies, followed by Disc+Bulge and Disc
(only) samples. The Pure Disc sample shows an alignment signal consistent with
zero. For the non-mass controlled samples (i.e., full samples without any
downselection to match the mass distributions), the alignment signals are
typically stronger by 10-20% compared to the mass-controlled samples, but the
trends with morphological subsample are very similar to those for the mass-
controlled samples, see Fig. 7 for direct illustration.

In § 4.2, we saw that the Ellipticals and Bulge (only) samples were more
aligned with their host DM halo compared to the Disc (only) and Pure Disc
samples. That finding, taken together with the results of this section,
suggests that the local alignment of galaxy and DM halo orientations
correlates with the large-scale alignment of galaxies (Heymans et al., 2006).
In order to quantify the small-scale intrinsic alignment amplitude, we employ
an empirical power-law model of the form:

$w_{g+}(r_{p})=w^{0}_{g+}\Bigg{(}\frac{r_{p}}{1\,h^{-1}\mathrm{Mpc}}\Bigg{)}^{-1},$
(29)

fitting the measured data below 1 Mpc. In Table 4 we present the best fitting
parameters of the power law model. All the fits for $w_{g+}$ were done
independently, as opposed to the NLA model where $w_{gg}$, $w_{g+}$, and
$w_{++}$ were fit jointly, for $r_{p}<$1 Mpc. Comparing with other work, we
note that Mandelbaum et al. (2006) (which is from SDSS) and Chisari et al.
(2015) (which is from the Horizon-AGN simulation) report values of
$w_{g+}^{0}=0.10^{+0.07}_{-0.07}$ and $0.13^{+0.03}_{-0.03}$, respectively,
for their whole sample. Both of these are about a factor of 7 smaller than
what we have measured; however, they had their power-law index as an
additional fit parameter, unlike ours, which was fixed at -1, and they report
$-0.59^{+0.65}_{-2.30}$ and $-0.75^{+0.25}_{-0.25}$ for their power-law index.
Additionally, we expect the amplitude of Mandelbaum et al. (2006) and Chisari
et al. (2015) to be lower since they used lower mass cuts of
$\log_{10}(M_{*}/M_{\odot})\sim 9$ and $\log_{10}(M_{*}/M_{\odot})\sim 8$,
respectively. Also, non-linear bias was not modeled, so it is hard to explain
the disagreement quantitatively. Thus, our fixed power-index of -1 is within
the range of both of these studies, despite the strong disagreement in the
amplitude $w_{g+}^{0}$.

| Red (full sample) | Blue | Red (mass controlled)
---|---|---|---
$A_{I}$ | $2.79^{+0.42}_{-0.43}$ | $0.13^{+0.35}_{-0.34}$ | $0.95^{+0.41}_{-0.45}$
Table 3: Best fitting parameters of NLA model when we split the whole sample by color rather than morphology. The two color bins were split in the $g-r$ color versus $r$ magnitude plane. The full Red sample has an alignment amplitude $A_{I}$ of $2.79^{+0.42}_{-0.43}$, whereas the Blue sample has $A_{I}$ consistent with zero. When mass-controlled, the Red sample’s alignment strength decreases by more than a factor of two. The mean $f_{\text{disc}}$ values were 0.63, 0.61 and 0.59 for the Blue, mass-controlled Red and full Red samples, respectively. The Blue sample was not mass controlled, since the Blue sample’s mass range was very narrow.  | All | Elliptical | Disc+Bulge | Bulge (only) | Disc (only) | Pure Disc
---|---|---|---|---|---|---
$w_{g+}^{0}$ | $0.79^{+0.07}_{-0.06}$ | $1.15^{+0.11}_{-0.11}$ | $0.67^{+0.16}_{-0.15}$ | $1.25^{+0.11}_{-0.11}$ | $0.59^{+0.17}_{-0.16}$ | $0.60^{+0.15}_{-0.15}$
Table 4: Best fitting amplitude parameter $w_{g+}^{0}$ from the power law
model fits to $w_{g+}$ for $r_{p}<1$ Mpc. These amplitudes show a very similar
trend as the NLA fit amplitudes on large scales: the Elliptical and Bulge
(only) samples have higher amplitudes, followed by the Disc+Bulge and Pure
Disc samples, with the Disc (only) sample having the weakest amplitudes.

##  5 Conclusions

In this work, we have investigated the intrinsic alignments of dynamically
decomposed/classified galaxies using the TNG100 hydrodynamical simulation from
the IllustrisTNG simulation suite. As a first test, we measured the 2D and 3D
shapes of simulated galaxies divided into samples based on their morphologies.
We found that the distributions of the absolute value of the ellipticities for
the Elliptical and Bulge (only) samples compare well with those from real
COSMOS data. However, the shapes of the disc populations differ from those of
the COSMOS dataset, overall showing more round shapes with the absolute
magnitude of the ellipticity distribution always being lower than 0.6 compared
to that for the COSMOS dataset, which stretches far beyond 0.6.

Second, we have investigated how the shapes of the different galaxy components
and types align with the shapes of their host halos and with their total
angular momentum vectors. Our investigation shows that the Elliptical and
Bulge (only) samples are more strongly aligned with their host DM halos
compared to the disc-related samples, in agreement with previous studies.
Additionally, when measuring the galaxy shape alignments with the total
angular momentum vectors of the DM halo and of the galaxy, all sample show
preferred alignment, with the angular momentum-dominated Pure Disc, Disc+Bulge
and Disc (only) samples showing a higher degree of alignment compared to the
Bulge (only) and the Elliptical samples.

Third, we have measured two-point correlation functions and fit them using the
NLA model for scales $r_{p}>4$ Mpc and a power-law model for scales $<1$ Mpc.
As expected, the Elliptical sample exhibits the highest value of
$A_{I}=3.47^{+0.57}_{-0.57}$, closely followed by the Bulge (only) sample with
$A_{I}=2.98^{+0.36}_{-0.37}$. These two samples have statistically consistent
alignments. In contrast, the Disc+Bulge and Disc (only) samples exhibit a
relatively low alignment of $A_{I}=1.13^{+0.37}_{-0.35}$ and
$A_{I}=1.02^{+0.34}_{-0.36}$,respectively. Lastly, the Pure Disc sample shows
an alignment consistent with zero. In order to account for potential trends in
intrinsic alignments with mass, we have implemented mass-dependent weights in
our two-point function calculations. Therefore, our results cannot be
explained by different host halo masses for morphologically-divided samples,
but rather must be explained by different alignments at fixed mass.

We have also considered samples divided based on their rest-frame color. The
full unweighted Red sample has a value of $A_{I}=2.79^{+0.42}_{-0.43}$, which
is smaller than the results for the Elliptical and Bulge (only) samples and
even smaller when mass-controlled at $A_{I}=0.95^{+0.41}_{-0.45}$. Though it
should be noted that the simulation produces higher number of red galaxies
compared to observations which makes interpretability of the color classified
results complicated. In contrast, the Blue sample (which was not mass
controlled due to its very narrow mass range) has an IA signal consistent with
zero.

Modeling and accounting for intrinsic alignments will be essential for the
next generation of cosmological surveys such as LSST, Roman and Euclid in
order to achieve unbiased cosmological parameter estimates given their
unprecedented statistical precision. In this study, we had tried to provide
better understanding of IA based on the stellar dynamics of galaxies. We have
shown that IA decreases as angular momentum starts to dominate the stellar
dynamics: large-scale IA becomes undetectable in our sample when angular
momentum almost completely dominates the stellar dynamics of galaxies.
Dispersion-dominated systems like Ellipticals and Bulges exhibit a
preferential alignment towards matter densities and are generally more
strongly aligned with the shape of their host DM halos. This finding that the
Elliptical and Bulge samples both have a strong IA signal on cosmological
scales provides further evidence for the Elliptical-Bulge likeness (or the
hypothesis that bulges are just scaled down ellipticals). In future work, it
would be useful to pursue a similar study in higher redshifts, in order to
track the evolution of IA in time, akin to the analysis of Bhowmick et al.
(2020), but with morphological separation. Another important aspect is
baryonic physics effect on IA; it would be valuable in the future to perform a
similar analysis on the gas structure of galaxies and investigate the effects
of star formation and active galactic nuclei on IA, similar to Tenneti et al.
(2016); Samuroff et al. (2020); these papers studied IA in different
simulations with different baryonic physics implementations.

Our findings indicate that intrinsic alignments are more complex than is
implied by the conventional Red-Blue (Elliptical-Disc) separation, with the
underlying stellar dynamics playing a significant role in determining the
alignments. For current and future weak lensing studies, IA models that
connect to the underlying disc fraction and relate to individual galaxy
components may improve the robustness and flexibility of the IA models.

## Data Availability

The data used in this paper is publicly available. The IllustrisTNG data can
be obtained through the website at https://www.tng-project.org/data/. The
catalog data with morphological decompositions of galaxies is available at
https://github.com/McWilliamsCenter/gal_decomp_paper

## Acknowledgements

We thank Simon Samuroff, Duncan Campbell, Markus Michael Rau, Scott Dodelson
for useful discussion that informed the direction of this work. We, also thank
Benjamin Joachimi for providing us the COSMOS ellipticity dataset. This work
was supported in part by the National Science Foundation, NSF AST-1716131 and
by a grant from the Simons Foundation (Simons Investigator in Astrophysics,
Award ID 620789). SS is supported by a McWilliams postdoctoral fellowship at
Carnegie Mellon University.

## References

  * Bhowmick et al. (2020) Bhowmick A. K., Chen Y., Tenneti A., Di Matteo T., Mandelbaum R., 2020, MNRAS, 491, 4116
  * Blazek et al. (2019) Blazek J. A., MacCrann N., Troxel M. A., Fang X., 2019, Phys. Rev. D, 100, 103506
  * Bridle & King (2007) Bridle S., King L., 2007, New Journal of Physics, 9, 444
  * Catelan et al. (2001) Catelan P., Kamionkowski M., Blandford R. D., 2001, MNRAS, 320, L7
  * Chisari et al. (2015) Chisari N., et al., 2015, MNRAS, 454, 2736
  * Chisari et al. (2017) Chisari N. E., et al., 2017, Mon. Not. Roy. Astron. Soc., 472, 1163
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Dodelson (2017) Dodelson S., 2017, Gravitational Lensing
  * Dubois et al. (2014) Dubois Y., et al., 2014, MNRAS, 444, 1453
  * Fortuna et al. (2021a) Fortuna M. C., Hoekstra H., Joachimi B., Johnston H., Chisari N. E., Georgiou C., Mahony C., 2021a, MNRAS, 501, 2983
  * Fortuna et al. (2021b) Fortuna M. C., et al., 2021b, A&A, 654, A76
  * Gadotti (2009) Gadotti D. A., 2009, MNRAS, 393, 1531
  * Georgiou et al. (2019) Georgiou C., et al., 2019, A&A, 628, A31
  * Hearin et al. (2017) Hearin A. P., et al., 2017, AJ, 154, 190
  * Heymans et al. (2006) Heymans C., White M., Heavens A., Vale C., van Waerbeke L., 2006, MNRAS, 371, 750
  * Hirata & Seljak (2004) Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526
  * Hirata et al. (2007) Hirata C. M., Mandelbaum R., Ishak M., Seljak U., Nichol R., Pimbblet K. A., Ross N. P., Wake D., 2007, MNRAS, 381, 1197
  * Jagvaral et al. (2021) Jagvaral Y., Campbell D., Mandelbaum R., Rau M. M., 2021, arXiv e-prints, p. arXiv:2105.02237
  * Joachimi et al. (2011) Joachimi B., Mandelbaum R., Abdalla F. B., Bridle S. L., 2011, A&A, 527, A26
  * Joachimi et al. (2013) Joachimi B., Semboloni E., Bett P. E., Hartlap J., Hilbert S., Hoekstra H., Schneider P., Schrabback T., 2013, MNRAS, 431, 477
  * Johnston et al. (2019) Johnston H., et al., 2019, A&A, 624, A30
  * Kiessling et al. (2015) Kiessling A., et al., 2015, Space Sci. Rev., 193, 67
  * Kilbinger (2015) Kilbinger M., 2015, Reports on Progress in Physics, 78, 086901
  * Kirk et al. (2015) Kirk D., et al., 2015, Space Sci. Rev., 193, 139
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, ApJ, 412, 64
  * Laurikainen et al. (2016) Laurikainen E., Peletier R., Gadotti D., 2016, Galactic Bulges
  * Li et al. (2019) Li Y., Singh S., Yu B., Feng Y., Seljak U., 2019, J. Cosmology Astropart. Phys., 2019, 016
  * Mandelbaum (2018) Mandelbaum R., 2018, ARA&A, 56, 393
  * Mandelbaum et al. (2006) Mandelbaum R., Hirata C. M., Ishak M., Seljak U., Brinkmann J., 2006, MNRAS, 367, 611
  * Mandelbaum et al. (2011) Mandelbaum R., et al., 2011, MNRAS, 410, 844
  * Marinacci et al. (2018) Marinacci F., et al., 2018, Mon. Not. Roy. Astron. Soc., 480, 5113
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Nelson et al. (2018) Nelson D., et al., 2018, Mon. Not. Roy. Astron. Soc., 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2
  * Okumura et al. (2009) Okumura T., Jing Y. P., Li C., 2009, ApJ, 694, 214
  * Pillepich et al. (2018a) Pillepich A., et al., 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich A., et al., 2018b, MNRAS, 475, 648
  * Samuroff et al. (2019) Samuroff S., et al., 2019, MNRAS, 489, 5453
  * Samuroff et al. (2020) Samuroff S., Mandelbaum R., Blazek J., 2020, arXiv e-prints, p. arXiv:2009.10735
  * Sandage et al. (1970) Sandage A., Freeman K. C., Stokes N. R., 1970, ApJ, 160, 831
  * Schneider & Bridle (2010) Schneider M. D., Bridle S., 2010, MNRAS, 402, 2127
  * Shao et al. (2016) Shao S., Cautun M., Frenk C. S., Gao L., Crain R. A., Schaller M., Schaye J., Theuns T., 2016, MNRAS, 460, 3772
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
  * Singh & Mandelbaum (2016) Singh S., Mandelbaum R., 2016, MNRAS, 457, 2301
  * Singh et al. (2015) Singh S., Mandelbaum R., More S., 2015, MNRAS, 450, 2195
  * Singh et al. (2017) Singh S., Mandelbaum R., Seljak U., Slosar A., Vazquez Gonzalez J., 2017, Monthly Notices of the Royal Astronomical Society, 471, 3827
  * Singh et al. (2021) Singh S., Yu B., Seljak U., 2021, MNRAS, 501, 4167
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2018) Springel V., et al., 2018, Mon. Not. Roy. Astron. Soc., 475, 676
  * Tenneti et al. (2014) Tenneti A., Mandelbaum R., Di Matteo T., Feng Y., Khandai N., 2014, MNRAS, 441, 470
  * Tenneti et al. (2015) Tenneti A., Singh S., Mandelbaum R., di Matteo T., Feng Y., Khandai N., 2015, MNRAS, 448, 3522
  * Tenneti et al. (2016) Tenneti A., Mandelbaum R., Di Matteo T., 2016, MNRAS, 462, 2668
  * Troxel & Ishak (2015) Troxel M. A., Ishak M., 2015, Phys. Rep., 558, 1
  * Velliscig et al. (2015a) Velliscig M., et al., 2015a, MNRAS, 453, 721
  * Velliscig et al. (2015b) Velliscig M., et al., 2015b, MNRAS, 454, 3328
  * Warren et al. (1992) Warren M. S., Quinn P. J., Salmon J. K., Zurek W. H., 1992, ApJ, 399, 405
  * Weinberg et al. (2013) Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87
  * van der Kruit & Freeman (2011) van der Kruit P. C., Freeman K. C., 2011, ARA&A, 49, 301