The Open Journal of Astrophysics submitted XXX; accepted YYY

††thanks: ⋆ E-mail<EMAIL_ADDRESS>
# An Empirical Model for Intrinsic Alignments: Insights from Cosmological
Simulations

Nicholas Van Alfen1 Duncan Campbell2 Jonathan Blazek1 C. Danielle Leonard3
Francois Lanusse4 Andrew Hearin5 Rachel Mandelbaum2 the LSST Dark Energy
Science Collaboration  1Department of Physics, Northeastern University,
Boston, MA 02115, USA  2McWilliams Center for Cosmology, Department of
Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA  3School of
Mathematics, Statistics and Physics, Newcastle University, Newcastle upon
Tyne, NE1 7RU, UK  4Université Paris-Saclay, Université Paris Cité, CEA, CNRS,
AIM, 91191, Gif-sur-Yvette, France  5Argonne National Laboratory, Lemont, IL
60439, USA

###### Abstract

We extend current models of the halo occupation distribution (HOD) to include
a flexible, empirical framework for the forward modeling of the intrinsic
alignment (IA) of galaxies. A primary goal of this work is to produce mock
galaxy catalogs for the purpose of validating existing models and methods for
the mitigation of IA in weak lensing measurements. This technique can also be
used to produce new, simulation-based predictions for IA and galaxy
clustering. Our model is probabilistically formulated, and rests upon the
assumption that the orientations of galaxies exhibit a correlation with their
host dark matter (sub)halo orientation or with their position within the halo.
We examine the necessary components and phenomenology of such a model by
considering the alignments between (sub)halos in a cosmological dark matter
only simulation. We then validate this model for a realistic galaxy population
in a set of simulations in the Illustris-TNG suite. We create an HOD mock with
Illustris-like correlations using our method, constraining the associated IA
model parameters, with the $\chi^{2}_{\rm dof}$ between our model’s
correlations and those of Illustris matching as closely as 1.4 and 1.1 for
orientation–position and orientation–orientation correlation functions,
respectively. By modeling the misalignment between galaxies and their host
halo, we show that the 3-dimensional two-point position and orientation
correlation functions of simulated (sub)halos and galaxies can be accurately
reproduced from quasi-linear scales down to $0.1\leavevmode\nobreak\
h^{-1}{\rm Mpc}$. We also find evidence for environmental influence on IA
within a halo. Our publicly-available software provides a key component
enabling efficient determination of Bayesian posteriors on IA model parameters
using observational measurements of galaxy-orientation correlation functions
in the highly nonlinear regime.

###### keywords:

Cosmology, Weak Gravitational Lensing, Intrinsic Alignments

##  1 Introduction

Weak lensing is an important probe for cosmological galaxy imaging surveys,
and will be of major interest for the Rubin Observatory’s Legacy Survey of
Space and Time (LSST, Ivezić et al., 2019) as well as other Stage IV surveys
such as Euclid (Scaramella et al., 2022) and Roman (Akeson et al., 2019).
These new surveys will be much more powerful than those in the past, and will
require correspondingly more precise analyses. The intrinsic alignments (IA)
exhibited by galaxies as a result of their alignment with large-scale
structure and other nearby galaxies can contaminate cosmic shear measurements,
driving the need for accurate IA modeling (Heymans et al., 2006; Joachimi et
al., 2013; Blazek et al., 2011; Krause et al., 2015).

IA has typically been modeled using analytic approaches, usually motivated
through cosmological perturbation theory (e.g. Hirata & Seljak, 2004; Bridle &
King, 2007; Blazek et al., 2015, 2019). However, these approaches are unable
to describe alignments in the fully nonlinear regime. To address this issue,
in recent years, interest in halo-based models of the intrinsic alignment of
galaxies has taken off, with Fortuna et al. (2020) updating previous work by
Schneider & Bridle (2010) to provide a robust theoretical description of
intrinsic alignment within dark matter halos. At the same time, efforts have
emerged to infuse dark-matter-only N-body simulations with mock galaxies which
have correlated alignments. These methods include halo-based modeling (e.g.
Joachimi et al., 2013; Hoffmann et al., 2022), tidal field descriptions (e.g.
Harnois-Déraps et al., 2022), and machine learning approaches (e.g. Jagvaral
et al., 2022). These simulation capabilities are important for validating our
approaches to mitigate IA effects in weak lensing surveys, providing
realistically complex mock data on which to test our analysis pipelines. These
simulations can also be used to directly model IA, including in situations
where analytic modeling is not feasible, such as higher-order or map-based
analyses. Finally, the resulting mock galaxy catalogs can also help us in
understanding and testing galaxy formation and evolution models which could
impact the large-scale alignment of galaxies.

In this work, we provide a new, flexible, and efficient halo-based method for
including IA in mock galaxy catalogs. This method is built on the
halotools111https://github.com/astropy/halotools package for generating mock
galaxy catalogs (Hearin et al., 2017) and can be applied within larger
simulation and analysis frameworks to produce mock galaxy catalogs or for IA
modeling (Van Alfen et al., 2023).

Models of the Halo Occupation Distribution (HOD) typically have components to
capture the number of galaxies within halos (occupation), as well as the
distribution of their positions and velocities within halos (phase space). In
this work, we add two new components to this framework: one component to model
the statistical distribution of galaxy orientations within their halos, and a
second component, the IA strength model, that characterizes the probabilistic
alignment between the galaxy orientation and some reference vector associated
with the parent halo. These new ingredients to HOD modeling are quite flexible
and allow for multivariate dependence on galaxy properties.

This paper is organised as follows: Section 2 presents background on two-point
clustering, giving the forms of the correlation functions used to analyze our
models. Section 3 describes the cosmologies of the simulations used in this
work and discusses how we determine our covariance. In Section 4 we introduce
the Dimroth-Watson distribution that we use to characterize the distribution
of galaxy orientations within their halos, and we describe the probabilistic
models of central and satellite galaxy alignments that will be used throughout
the rest of the paper. In Section 5, we examine the impact of altering central
and satellite alignment strengths. This provides a baseline for understanding
the simplest effect our alignment strengths have on our model. In Section 6,
we use high-resolution N-body simulations to investigate the accuracy with
which subhalo-orientation correlation functions can be captured by our HOD-
type models using radial alignments, as subhalo orientations are often not
available. We use the Illustris-TNG hydrodynamical simulation to test the
flexibility of our model of alignment strength in Section 7. Section 8
summarizes our results and future work.

##  2 Two-Point Correlation Functions and the Halo Occupation Distribution

In this section, we introduce the basic equations of the analytical model of
halo occupation statistics. Historically, these analytic expressions have
provided the principal framework underlying theoretical predictions of two-
point galaxy clustering in the nonlinear regime (e.g., Berlind et al., 2003;
Yang et al., 2003). In the present work, the theoretical prediction pipeline
we introduce is purely simulation-based, and so our Monte Carlo methods for
computing the integrals in this section are quite distinct from the algorithms
implemented in other widely-used halo model libraries such as the Core
Cosmology Library (CCL, Chisari et al., 2019) or HMCode (Mead et al., 2015).
We defer a comparative discussion of these methodologies to Section 8, and we
include the analytical expressions here to elucidate the connection between
simulation-based and analytical pipelines.

###  2.1 Theory: The Halo Model and correlation functions

Here we overview the theoretical equations for the two-point correlation
functions. The functions are given to provide a conceptual background, but our
analysis does not use these formulae to calculate correlation functions.
Instead, we measure the correlations from mock data using the estimators
described in Section 2.2.

Below, we see the theoretical formulation of the position-position
correlations.

$\displaystyle\xi_{\rm gg}(r)=\xi_{\rm gg}^{\rm 1h}(r)+\xi_{\rm gg}^{\rm
2h}(r)$ (1)

$\displaystyle\xi_{\rm gg}^{\rm 1h}(r)=\xi_{\rm cs}^{\rm 1h}(r)+\xi_{\rm
ss}^{\rm 1h}(r)$ (2)

$\displaystyle\xi_{\rm gg}^{\rm 2h}(r)=\xi_{\rm cc}^{\rm 2h}(r)+\xi_{\rm
cs}^{\rm 2h}(r)+\xi_{\rm ss}^{\rm 2h}(r)$ (3)

$\displaystyle\begin{split}\xi_{\rm gg}^{\rm 2h}(r)=\frac{1}{\bar{n}_{\rm
g}^{2}}\int\rm dM_{1}\int\rm dM_{2}&\frac{\rm dn}{\rm dM_{1}}\frac{\rm dn}{\rm
dM_{2}}\xi_{\rm hh}(r|M_{1},M_{2})\\\ &\times\langle N_{\rm
g}|M_{1}\rangle\langle N_{\rm g}|M_{2}\rangle,\end{split}$ (4)

where $\frac{\rm dn}{\rm dM}$ is the halo mass function (the differential
distribution of number density of halos $n$ as a function of halo mass $M$),
$\langle N_{\rm g}|M\rangle$ is the expectation value of the number of
galaxies in a halo of mass $M$, and $\bar{n}_{\rm g}$ is the overall number
density of the population of galaxies being considered. In the above
equations, subscripts indicate the two populations being correlated. The
subscript “g” refers to all galaxies, “h” refers to halo, “c” refers to
central galaxy (the most massive galaxy at/near the center of the host dark
matter halo), and “s” refers to satellite galaxy (a smaller galaxy in a host
dark matter halo besides the central galaxy). The superscripts “1h” and “2h”
denote the scale contributing to the correlations. Namely, “1h” refers to the
1-halo regime (galaxies within the same halo), “2h” refers to the 2-halo
regime (galaxies residing in separate halos), and no superscript indicates
correlations from all scales. In this paper, we use 1 Mpc/$h$ as the
approximate delineation between the 1-halo and 2-halo regimes.

The assumptions underlying Eq. (4) are:

  * •

$r\gg R_{\rm vir}^{\rm max}$

  * •

On very large scales, $\xi_{\rm hh}(r)$ varies slowly across the length scale
of halos

where $R_{\rm vir}^{\rm max}$ refers to the largest virial radius of all halos
being considered for the 2-halo contribution to the correlation.

$\displaystyle\begin{split}\xi_{\rm gg}^{\rm 1h}(r)=\frac{1}{\bar{n}_{\rm
g}^{2}}\int\rm dM\frac{\rm dn}{\rm dM}\int\rm
d^{3}&x\lambda(\bm{x}|M)\lambda(\bm{x}+\bm{r}|M)\\\ &\times\langle N_{\rm
g}(N_{\rm g}-1)|M\rangle\end{split}$ (5)

In Eq. (5), $\lambda(\bm{x}|M)$ is the satellite profile, i.e., the unit-
normalized spatial distribution of satellite galaxies within a halo of mass
$M$.

Above, we looked at the position-position correlation functions. Here, we
overview the orientation-orientation and position-orientation correlation
functions. Following typical naming conventions, we use the subscripts “e” and
“d” to refer to ellipticity and direction respectively, where direction refers
to the separation between galaxies, which in some contexts we call “position.”
As discussed below, we do not use full shape information in this work, so our
use of “ellipticity” refers in practice to orientation.

$\displaystyle\begin{split}\xi_{\rm ee}^{\rm 2h}(r)=\frac{1}{\bar{n}_{\rm
g}^{2}}\int\rm dM_{1}\int\rm dM_{2}&\frac{\rm dn}{\rm dM_{1}}\frac{\rm dn}{\rm
dM_{2}}\xi_{\rm hh}(r|M_{1},M_{2})\\\ &\times\langle|\bm{e}_{\rm
1}\cdot\bm{e}_{\rm 2}|^{2}\rangle\end{split}$ (6)

In Eq. (6), $\bm{e}$ is the three-dimensional orientation vector of the
galaxy. Similarly, the 1-halo part of the ee correlation function and both
parts of the ed correlation function are

$\displaystyle\begin{split}\xi_{\rm ee}^{\rm 1h}(r)=\frac{1}{\bar{n}_{\rm
g}^{2}}\int\rm dM\frac{\rm dn}{\rm dM}\int\rm
d^{3}&x\lambda(\bm{x}|M)\lambda(\bm{x}+\bm{r}|M)\\\ &\times\langle|\bm{e}_{\rm
1}\cdot\bm{e}_{\rm 2}|^{2}\rangle,\end{split}$ (7)

$\displaystyle\begin{split}\xi_{\rm ed}^{\rm 2h}(r)=\frac{1}{\bar{n}_{\rm
g}^{2}}\int\rm dM_{1}\int\rm dM_{2}&\frac{\rm dn}{\rm dM_{1}}\frac{\rm dn}{\rm
dM_{2}}\xi_{\rm hh}(r|M_{1},M_{2})\\\ &\times\left(\langle|\bm{e}_{\rm
1}\cdot\hat{r}|^{2}\rangle+\langle|\bm{e}_{\rm
2}\cdot\hat{r}|^{2}\rangle\right),\end{split}$ (8)

and

$\displaystyle\begin{split}\xi_{\rm ed}^{\rm 1h}(r)=\frac{1}{\bar{n}_{\rm
g}^{2}}\int\rm dM\frac{\rm dn}{\rm dM}\int\rm
d^{3}&x\lambda(\bm{x}|M)\lambda(\bm{x}+\bm{r}|M)\\\
&\times\left(\langle|\bm{e}_{\rm
1}\cdot\hat{r}|^{2}\rangle+\langle|\bm{e}_{\rm
2}\cdot\hat{r}|^{2}\rangle\right).\end{split}$ (9)

For $\xi_{\rm ee}(r),$ we can decompose the final term in the integrand into
central and satellite contributions:

$\displaystyle\begin{split}\langle|\bm{e}_{\rm 1}\cdot\bm{e}_{\rm
2}|^{2}\rangle&=f_{\rm cen,1}f_{\rm cen,2}\langle|\bm{e}_{\rm
cen,1}\cdot\bm{e}_{\rm cen,2}|^{2}\rangle\\\ &+f_{\rm sat,1}f_{\rm
sat,2}\langle|\bm{e}_{\rm sat,1}\cdot\bm{e}_{\rm sat,2}|^{2}\rangle\\\
&+f_{\rm cen,1}f_{\rm sat,2}\langle|\bm{e}_{\rm cen,1}\cdot\bm{e}_{\rm
sat,2}|^{2}\rangle\\\ &+f_{\rm sat,1}f_{\rm cen,2}\langle|\bm{e}_{\rm
sat,1}\cdot\bm{e}_{\rm cen,2}|^{2}\rangle\end{split}$ (10)

where $f_{\rm cen,\it i}$ and $f_{\rm sat,\it i}$ are the fraction of galaxies
in sample $i$ that are centrals and satellites, so that $f_{\rm cen,\it
i}+f_{\rm sat,\it i}=1.$

For $\xi_{\rm ed}(r),$ this central-satellite decomposition takes the form:

$\displaystyle\begin{split}\langle|\bm{e}_{\rm
i}\cdot\hat{r}|^{2}\rangle=f_{\rm cen,\it i}\langle|\bm{e}_{\rm cen,\it
i}\cdot\hat{r}|^{2}\rangle+f_{\rm sat,\it i}\langle|\bm{e}_{\rm sat,\it
i}\cdot\hat{r}|^{2}\rangle\end{split}$ (11)

where $\bm{r}=r\cdot\hat{r}$ refers to the position vector between two
galaxies.

###  2.2 Two-point Correlation Functions Estimators

Below are the estimators for the position-position, position-orientation, and
orientation-orientation correlation functions that we use to analyze our
model. While earlier we showed the theoretical formulation, here we focus on
the estimators we use on mock galaxy catalogs. As stated earlier, to remain
consistent with the names used in the literature, we have described these
functions using “ellipticity” rather than “orientation.” However, as defined
in Eqs. (6)-(9), and in the estimators below, only the galaxy orientation is
used. Our method, as presented here, does not include full galaxy shape
information, which is the topic of ongoing work. For consistency with the
function names, we sometimes use the word “ellipticity” when referring to
galaxy orientations. We also note that the functions utilized here are
convenient for analysing simulations where three-dimensional position and
orientation information is available. However, for analysis of observed data,
where only projected shape information is present, other correlation functions
are typically used (e.g. Mandelbaum et al., 2006).

The galaxy-galaxy two-point correlation function is given by:

$\xi_{gg}(r)=\frac{DD(r)}{RR(r)}$ (12)

where $DD(r)$ is the number of galaxy pairs separated by $r$ and $RR(r)$ is
the expected number of pairs for a random distribution. We use this estimator
rather than the Landy-Szalay estimator (Landy & Szalay, 1993) even though it
is somewhat sub-optimal in some cases (Singh et al., 2017). This estimator is
faster, and the periodic nature of our box allows us to use analytic randoms,
negating much of the sub-optimality.

The ellipticity-direction (ED) correlation function is defined as:

$\omega(r)=\langle|\hat{e}({\bf x})\cdot\hat{r}|^{2}\rangle-\frac{1}{3}$ (13)

The ellipticity-ellipticity (EE) correlation function is defined as:

$\eta(r)=\langle|\hat{e}({\bf x})\cdot\hat{e}({\bf x}+{\bf
r})|^{2}\rangle-\frac{1}{3}$ (14)

##  3 Simulations and Covariance

We construct our HOD models using halotools, a library for simulation-based
predictions of the galaxy–halo connection. Catalogs of dark matter halos are
the root data product underlying halotools-based predictions, and throughout
the paper we make use of catalogs provided by the library, although we note
that halotools provides full support for user-supplied simulation data. Upon
ingesting a catalog of simulated halos, halotools populates each halo with a
synthetic galaxy population according to an occupation model, distributing
galaxies within each halo according to an assumed phase space model.
Orientations are bestowed upon each synthetic galaxy according to the
alignment model outlined in Section 2 and described in detail in the sections
below.

Throughout the paper, we will distinguish between a host halo, which is a
distinct gravitationally self-bound collection of simulated dark matter
particles, and a subhalo, a gravitationally self-bound object that is located
inside the boundary of a larger host halo. We will use the term (sub)halo
whenever referring to a collection of host halos and subhalos together and/or
when referring to an individual object that may be either a subhalo or a host
halo.

###  3.1 Simulations

Simulation | Particle Mass | $\Omega_{m,0}$ | $\sigma_{8}$ | $n_{s}$ | $h$ | $L_{\rm box}$ | $z$
---|---|---|---|---|---|---|---
| ($M_{\odot}/h$) |  |  |  |  | (Mpc/h) |
Bolshoi-Planck | $\sim 10^{8}$ | 0.30711 | 0.82 | 0.96 | 0.70 | 250 | 0
(Bolplanck) |  |  |  |  |  |  |
Small MultiDark Planck | $\sim 10^{8}$ | 0.307115 | 0.8228 | 0.96 | 0.6777 | 400 | 0
(SMDPL) |  |  |  |  |  |  |
Illustris TNG-300 | $\sim 10^{10}$ | 0.3089 | 0.8159 | 0.9667 | 0.6774 | 205 | 0
Table 1: Cosmology and other parameters for the three simulations described in
Section 3. Additional information about each simulation can be found in said
section.

In this paper, unless otherwise stated, whenever we create a catalog using
halotools, we do so using the available halo catalogs from the Bolshoi-
Planck222https://doi.org/10.17876/cosmosim/bolshoip (BolPlanck) simulation
output at $z=0$, the cosmology of which is given in Table 1 (Klypin et al.,
2011).

While most of the paper uses BolPlanck, Figures 3 and 4 use the Small
MultiDark Planck (SMDPL) simulation with cosmology described in Table 1. The
simulation is evolved by solving for gravitational interactions only using the
L-GADGET-2 code, a version of the publicly available cosmological code
GADGET-2 (Springel, 2005) with a force resolution of $1.5\leavevmode\nobreak\
{\rm kpc}/h$. This simulation belongs to the series of MultiDark simulations
with Planck cosmology. More details for this simulation are described in
Klypin et al. (2016).

For both Bolplanck and SMDPL, (sub)halos are found using the phase-space halo
finder ROCKSTAR (Behroozi et al., 2013a), which uses adaptive, hierarchical
refinement of friends-of-friends groups in six phase-space dimensions and one
time dimension, and tracked over time using the Consistent Trees algorithm
(Behroozi et al., 2013b). As demonstrated in Knebe et al. (2011, 2013), this
results in a very robust tracking of (sub)halos (also see Jiang & van den
Bosch, 2016).

The essential orientation properties used in this paper (the x, y, and z
values of the unit vector for the major axis of the ellipsoid) are also
generated with the ROCKSTAR halo finder, shown in Appendix B of Rodríguez-
Puebla et al. (2016). For a more detailed look at shape finding, see Appendix
A.

Halos in this catalog are defined to be spherical volumes centered on a local
density peak (SO hereafter), such that the average density inside the sphere
is $\bar{\rho}_{h}(z)=\Delta_{\rm vir}(z)\rho_{\rm m}(z)$. Here $\rho_{\rm
m}(z)=\Omega_{\rm m}(z)\rho_{\rm crit}(z)$, where $\rho_{\rm
crit}(z)=3H(z)2/8\pi G$ is the critical energy density of the Universe, and
$\Delta_{\rm vir}(z)$ is given by a fitting function (Bryan & Norman, 1998):

$\Delta_{\rm
vir}(z)=\left[18\pi^{2}-82\Omega_{\Lambda}(z)-39\Omega_{\Lambda}^{2}(z)\right]\Omega_{\rm
m}^{-1}$ (15)

For the Planck cosmology, $\Delta_{\rm vir}(z)\simeq 360$. The radius of each
such sphere defines the virial radius $R_{\rm vir}$ of the halo, which is
related to the mass of the halo via $M_{\rm vir}=(4/3)\pi R_{\rm
vir}^{3}\bar{\rho}_{h}$. Additionally, subhalos in this catalog are distinct,
self-bound structures whose centers are found within the virial radius of a
more massive host halo.

In Section 7, we compare results from our HOD models (built using the
BolPlanck catalogs described above) to the Next Generation Illustris
Simulations333http://www.tng-project.org (Illustris-TNG), a suite of
hydrodynamical simulations of galaxy formation in cosmological volumes (Nelson
et al., 2017; Pillepich et al., 2017; Naiman et al., 2018; Springel et al.,
2017; Marinacci et al., 2018). Illustris resolves galaxy shapes, therefore the
IA seen in Illustris is modeled directly by hydrodynamic interactions. The
Illustris-TNG runs used in this work are two uniform mass resolution
cosmological volume simulations with side lengths 205 $h^{-1}{\rm Mpc}$, one
“full physics” run, TNG300, including all of the complex physics of galaxy
formation, and a gravity-only counterpart, TNG300-Dark. The initial conditions
of the simulations were set at z = 127 using the Zeldovich approximation. The
adopted cosmological parameters are given in Table 1.

###  3.2 Covariance

Building our HOD model for IA is inherently stochastic, as both the selection
of halos being populated and the alignment of galaxies (discussed below) are
drawn from distributions. When fitting this model to “data” (which in this
work include both measurements of the simulated halos being populated as well
as galaxies in another simulation such as Illustris-TNG), we must account for
the fact that both the data and model have stochasticity. When fitting to the
halo measurements (Section 6), we find that calculating the covariance using
multiple realizations of the galaxy alignments, without changing the halo
occupation, is sufficient to allow a test of the model performance while
avoiding double counting the variance from the occupation step. When fitting
to Illustris-TNG (Section 7), we calculate the covariance from Illustris-TNG
and average over multiple model realizations to remove model variance. Further
details concerning our covariance are given in Appendix D.

##  4 Modeling Alignments

In this paper, we make the assumption that both galaxies and (sub)halos can be
modeled as triaxial homologous ellipsoids (see Appendix A for details). The
orientation of galaxies/halos can then be entirely described by specifying the
vectors associated with the minor, intermediate, and semi-major axes of
ellipsoids. Here we focus on modeling the orientation of a single axis of the
ellipsoid describing a galaxy within its host halo, but we note that the
following framework could be extended to model a three-dimensional galaxy
orientation and shape, albeit with increased complexity.

Our strategy will be to statistically model the misalignment angle between a
vector specifying the orientation of each galaxy’s (sub)halo, and a vector
specifying the orientation of the galaxy, as shown in Figure 1.
Mathematically, if we assume that the vector specifying the orientation of
each galaxy’s host (sub)halo is the unit vector, $\hat{z}$, we can model the
orientation vector of the galaxy by specifying the angular coordinates on the
unit-sphere.

Figure 1: Cartoon of alignment model when subhalos with accurate shapes are
available (top) and without subhalos (bottom). In both cases, the misalignment
angle is determined with respect to some reference axis, which differs for
satellite galaxies depending on the model used, either using subhalo alignment
(top), or radial alignment (bottom).

For this purpose, we utilize the Dimroth-Watson distribution, which specifies
the spherical polar coordinates of a unit vector. We chose the Dimroth-Watson
distribution as it provides a maximum entropy distribution on a sphere while
accounting for the 180 degree symmetry we see in our galaxy orientations
(Watson, 1965). These criteria are difficult to meet with other distributions,
but the Dimroth-Watson distribution naturally fits the role for directional,
axially symmetric systems on a sphere (Sra & Karp, 2013). The probability
distribution for the polar angle, $\theta$, and the azimuthal angle, $\phi$,
are given by

$P(\theta,\phi)=\frac{B(\kappa)}{2\pi}e^{-\kappa\cos^{2}(\theta)}\sin(\theta)\mathrm{d}\theta\mathrm{d}\phi$
(16)

where the normalization factor is given by:

$B(\kappa)=\frac{1}{2}\int_{0}^{1}e^{-\kappa t^{2}}\mathrm{d}t$ (17)

Note that the azimuthal angle, $\phi$, is modeled as a uniform distribution,
so that the angle $\theta$ is then by definition the misalignment angle
between galaxy and host (sub)halo orientation, $\theta_{\rm MA}$. The degree
to which galaxies and (sub)halos align is then controlled by the $\kappa$
parameter.

It is convenient to re-parameterize the strength of alignment $\kappa$ by:

$\mu=\frac{-2\tan^{-1}(\kappa)}{\pi}$ (18)

such that $\mu=1,0,-1$ corresponds to perfect alignment, random alignment, and
perfect anti-alignment, respectively. As illustrated in Figure 2, a positive
value of $\mu$ corresponds to a preferential parallel alignment between a
galaxy and its host (sub)halo. A negative value of $\mu$ corresponds to a
preferential perpendicular alignment between a galaxy and its host halo.

The distribution of misalignment angles is symmetric about $\cos(\theta_{\rm
MA})=0$, in accord with the symmetries we have assumed of the shapes of
galaxies and halos. There is thus no difference in the morphological
orientations between a galaxy rotated by $180^{\circ}$ and its original
orientation. In principle, for disk galaxies whose orientation is determined
by angular momentum, it may be interesting to distinguish between parallel and
anti-parallel angular momentum vectors. Although the difference should not
impact the observable galaxy ellipticities, correlations between the angular
momentum direction and the density field can allow us to probe tidal torque
theory (e.g. Lee, 2019) – we leave exploration of this and related effects to
future work.

Figure 2: Dimroth-Watson distribution model, equation (16), for the
distribution of galaxy–halo orientation misalignment angles. Different color
lines show distributions with different alignment strengths, $\mu$, given by
equation (18).

###  4.1 Galaxy Alignment Models

For all investigations, we align the major axes of central galaxies with the
major axes of their host halos.

For satellite galaxies, we explore two models in this paper for determining
their orientations:

  1. 1.

satellites are oriented relative to that of their host dark matter subhalo,
called subhalo alignment (as seen in the upper image of Figure 1);

  2. 2.

satellites are oriented relative to the host halo centric radial vector,
called radial alignment (illustrated in the lower image of Figure 1). In this
case, there are two different alignment strengths explored.

    1. (a)

Constant alignment strength, where the alignment strength is the same for all
satellite galaxies.

    2. (b)

Distance-dependent alignment strength, where the alignment strength of a
satellite depends on the distance to its central galaxy.

###  4.2 Building an HOD Mock

To build an HOD-based mock catalog, we use the HODModelFactory in halotools.
Central and satellite galaxies have separately-defined occupation models that
determine the number of galaxies of each population present in a given halo.
Additionally, each population has its own model for the intra-halo phase space
distribution. Throughout the paper, central galaxies are placed at the center
of the host halo with the same peculiar velocity as the halo. For satellite
galaxies, we consider two different classes of phase space models. The first
is defined by the NFWPhaseSpace class: the spatial distribution of satellites
is spherically symmetric and has a radial density profile given by the NFW
distribution (Navarro et al., 1997); satellite velocities are isotropic and
exhibit a velocity dispersion profile determined by the Jeans equation. The
second satellite model we consider is given by the SubhaloPhaseSpace class, in
which the position and velocity of each satellite is defined by a subhalo that
has been randomly selected from the simulated halo, preferentially selecting
the most massive subhalos first; when $N_{\rm sat}$ is smaller than the number
of subhalos in the host halo, the lowest-mass subhalos are left unoccupied,
and when $N_{\rm sat}$ exceeds the number of available subhalos, the remaining
satellites are distributed according to the NFWPhaseSpace class. The
SubhaloPhaseSpace class provides a way for the synthetic satellite population
to inherit the complexity of the intra-halo distributions of simulated
subhalos, while at the same time retaining the flexibility of the HOD to
occupy host halos with a variable number of satellites.

##  5 Impact of Galaxy–Halo Misalignment on Orientation Correlation Functions

In this section, we use our halo model of IA to provide a pedagogical
investigation of how galaxy–halo misalignment manifests in orientation
correlation functions. We will separately study effects arising from central
and satellite galaxies, and we will consider the full range of spatial scales,
spanning both 1-halo and 2-halo regimes. Our general strategy for this
investigation is to begin by orienting galaxies to be perfectly aligned with
the reference vector of the (sub)halo; we then leverage the flexibility of our
probabilistic model by programmatically weakening the alignment strength, and
studying how the orientation correlation function changes as a result.

###  5.1 Central and Satellite Galaxy Misalignment

We begin by examining the effect of varying the central galaxy – host halo
alignment strength on the orientation correlation functions. In Figure 3, we
show the ellipticity–direction (ED) on the left and ellipticity–ellipticity
(EE) correlation functions on the right. In both plots, we show subhalos and
host halos with a peak mass greater than $10^{12}\leavevmode\nobreak\
M_{\odot}/h$ as black points with error bars. Errors on these measurements are
estimated using jackknife re-sampling of the simulation box, by splitting into
$4^{3}$ equal volume sub-samples.

Also shown in Figure 3 are color-coded lines, showing $\omega^{\rm g}_{\rm
gg}$ and $\eta^{\rm g}_{\rm gg}$ (position-orientation and orientation-
orientation correlations for all galaxies) for model galaxies in the same
sample of (sub)halos with varying levels of central-host halo major axis
alignments. The satellite–subhalo alignment strength is kept fixed at full
strength, $\mu_{\rm sat}=1$. Because there is still significant uncertainty in
the EE measurement in a $400\leavevmode\nobreak\ {\rm Mpc}$ box, we use a
fitting function to smooth the results (see Appendix B for details). For
maximum alignment, $\mu_{\rm cen}=1$, the correlation functions reproduce the
(sub)halo measurement by construction.

Figure 3: The effect of degrading the central galaxy-host halo alignment
strength on the ED (left) and EE (right) orientation correlation functions.
The black points with error bars (which are often smaller than the points) are
measurements made directly on (sub)halos. The colored lines are for model
galaxies with varying levels of alignment strength between central galaxies
and their host halo, from random alignments (purple) to perfect alignments
(pink). In each model, satellite galaxies take on the same orientation as
their subhalo. On the right, a fitting function is used to smooth the results,
letting us focus on the overall trend from lowering the satellite alignments
which tends to be noisy for ED correlations. We also see a downward bump
around the 1-halo to 2-halo scale transition in the right panel. This suggests
to us an effect related to the satellite–satellite correlations, with the
sudden upturn going into the1-halo regime being related to satellite pairs
within a halo becoming more significant. This effect also only becomes
apparent once the central alignment strength has become weak enough that the
satellite alignments become more important. Figure 4: The effect of degrading
the satellite galaxy-subhalo alignment strength on the ED (left) and EE
(right) orientation correlation functions. The black points with error bars
are measurements made directly on (sub)halos. The colored lines are for model
galaxies with varying levels of alignment strength between satellite galaxies
and their subhalo, from random alignments (purple) to perfect alignments
(pink). In each model, central galaxies take on the same orientation as their
host halo. On the right, a fitting function is used to smooth the results.

Next we show the effect of varying the satellite–subhalo alignment strength on
the orientation correlation functions. Similar to previous section, in Figure
4 we show the ED on the left and EE on the right, with correlation functions
for subhalos and host halos given as black points with error bars. In Figure
4, the color-coded lines correspond to varying levels of satellite–subhalo
alignments, while the central–halo alignment strength is kept fixed at full
strength, $\mu_{\rm cen}=1$. Again, when $\mu_{\rm sat}=1$, the (sub)halo
orientation correlation functions are reproduced by construction.

Broadly speaking, we can see that central galaxy alignments produce a much
larger effect on the correlation functions relative to satellites, especially
at larger scales. The satellite correlations in Figure 4 show noticeable
effects at smaller scales, but become much less important on scales larger
than $\sim 1$ Mpc/$h$. Meanwhile, Figure 3 shows that changing the alignment
strength of centrals has a noticeable effect at all scales.

##  6 Modeling Satellite Orientations without Subhalo Orientations

Our analyses in the previous sections have relied upon the SubhaloPhaseSpace
model of HOD satellite galaxies, in which the orientation of a synthetic
satellite was variably-correlated with the major axis of its parent subhalo.
But relying on subhalo catalogs as the core data product of the model
introduces a significant source of systematic error arising from uncertainty
in subhalo-finding and numerical/resolution issues on small scales (e.g., van
den Bosch et al., 2018; Campbell et al., 2018). The resolution requirements on
subhalo finding are steep (Mansfield & Avestruz, 2021), and so making
theoretical predictions with subhalo-based models leads to costly demands on
the Gpc-sized simulations that are required to analyze cosmological surveys.
These considerations motivate an alternative formulation of satellite
orientations that is based on purely synthetic satellites, as in standard HOD
type models.

In this section we carry out a targeted study of the flexibility of IA models
of satellites that do not rely upon subhalos. Our primary assumption will be
that the orientation of a satellite galaxy tends to align radially towards the
center of the host halo. We examine this assumption in the simplest alignment
model in §6.1, and in subsequent subsections we proceed to incorporate
additional complexity and study the impact on orientation correlation
functions. In order to assess the flexibility required by an HOD-type model of
satellite orientations, we begin by focusing on recovering the orientation
correlation functions exhibited by (sub)halos at $z=0$ in the gravity-only
Bolshoi-Planck simulation. Following this preliminary investigation, in
Section 7, we will proceed to apply our IA model to full-physics
hydrodynamical simulations (described in Section 3).

###  6.1 Radial Alignment Modeling

Constant radial misalignment uses a single, constant alignment strength for
all galaxies, giving them all the same Dimroth-Watson distribution for their
radial alignments. A distance-dependent alignment strength changes the shape
of the Dimroth-Watson distribution for each satellite galaxy depending on its
distance from its central galaxy. We choose a form of this relationship given
by a power law – equation (19). Because $\mu\in[-1,1]$, this equation is
limited to that range as well, assigning -1 to all $\mu<-1$, and 1 to all
$\mu>1$.

$\mu=a\,(r/r_{\textrm{vir}})^{\gamma}$ (19)

Since $\mu$ refers to how strongly galaxies tend to align with respect to a
given reference vector, the $\mu$ values in this section refer to how strongly
satellite galaxies align with the radial vector between the satellite and the
center of its host halo. In previous sections, we used (sub)halo major axes as
reference vectors, so matching the halo correlation in that context would mean
perfectly aligning with the halo major axes. Here, however, matching the halo
correlations is less straightforward. We would only expect a value of $\mu=1$
in this case if subhalo shapes are perfectly aligned with their radial
vectors.

In any given mock galaxy catalog that we generate, there are roughly three
times as many central galaxies as there are satellite galaxies (i.e., the
satellite fraction is roughly 25% for the galaxy populations we consider), so
looking at the correlation functions between all galaxies dilutes the effect
of the satellites, and is heavily influenced by the central correlations. In
order to compare the effects of these satellite alignments strengths, we look
at $\omega^{\rm g}_{\rm cs}$ (the correlation functions for the central galaxy
position and satellite galaxy orientation) rather than the full ED correlation
function for all galaxies to all galaxies. In this way, we can better see the
effect of only the satellite orientations. All fits discussed in this section
are fits to central position, satellite orientation correlation functions. We
do note that in the 2-halo regime, it is expected that the overall alignment
correlation functions are dominated by the impact of central-central alignment
only (Schneider & Bridle, 2010; Fortuna et al., 2020); this is an important
caveat to bear in mind in the following discussions in regards to the impact
of our model’s ability or inability to fit central-satellite correlations well
in the 2-halo regime.

As shown in central plot of Figure 5, and discussed in more detail in Section
6.4, the distance-dependent and constant radial alignment strengths perform
comparably. The correlation function from the distance-dependent radial
alignment may match that of its halos slightly better, but adds an extra
degree of complexity. We also see that the correlation functions match much
better in the 1-halo regime (here, defined as anything within 1 Mpc/$h$) than
overall.

Figure 5: Left: $\xi_{\rm cs}$ correlation function (central galaxy position
correlated with satellite galaxy position). Galaxy correlation functions are
shown in the solid red and dashed blue lines, with the corresponding halo
correlation function shown in black dots. Here, all three share the same $\xi$
(position–position) correlation function as the only difference in the model
is the alignment. Middle: $\omega_{\rm cs}$ correlation function (central
galaxy position and satellite galaxy orientation). Here, we see a slight
difference between $\omega^{\rm g}_{\rm cs}$ (the galaxy correlation
functions) and $\omega^{\rm h}_{\rm cs}$ (the corresponding halo correlation
function) with the distance-dependent alignment (solid red line) performing
slightly better at small scales than the constant alignment strength (shown as
a dashed blue line). Right: $\eta_{\rm cs}$ correlation function (central
galaxy orientation and satelleite galaxy position). In all panels (and in
figures below), the shaded regions show the 1-sigma confidence interval of
each model.

While we haven’t used the subhalo orientations to build our model, mimicking a
scenario where they are not available, we do have them and can use them for
comparison purposes now. Comparing to the subhalo correlation functions, we
see in Figure 5 that the distance-dependent alignment strength produces
correlation functions that match within $20\%$ of the halo correlation
functions overall, and within $6\%$ in the 1-halo regime. Meanwhile, the
constant alignment strength produces results that match within $20\%$ overall,
and within $8\%$ in the 1-halo regime. While there may be a slight preference
for a distance-dependent radial alignment, the two methods perform comparably.
Certainly, constant radial alignment provides a simpler model and may be
sufficient for many needs.

###  6.2 Fitting to the 1-halo and 2-halo Regimes

We have compared the galaxy correlation functions to those of their respective
halos for both the full range, as well for the 1-halo regime alone, generally
finding that our models fared better on small scales. Because of this, we
restrict ourselves to the 1-halo regime when finding the best-fit parameters
for the constant and distance-dependent radial alignment models. In this
section, we further justify that decision by showing the limits of our model
when considering satellite alignments. For both the constant radial alignment
strength and the distance-dependent radial alignment strength cases, we ran
MCMC chains and compared the galaxy correlation functions to the halo
correlation functions in the 1-halo regime. To run these chains, and any time
we use MCMC, we use the emcee python package (Foreman-Mackey et al., 2013).

Our attempts to fit to the correlation function across all scales led to a
wide spread of $a$ and $\gamma$ values for equation (19), producing extreme
values for $\mu$, and resulting in poor fits and no clear best-fit parameter
values. Our attempts to fit to large-scale correlations alone (i.e., the
2-halo regime) resulted in similar issues. Meanwhile, when restricting
attention to the small-scale correlations alone, our fitter enjoyed higher-
accuracy success. These results, and the fact that we underestimate the 2-halo
correlations when fitting to the 1-halo regime only, imply that there is extra
environmental information needed to reproduce correlation functions between
central position and satellite orientation in the 2-halo regime, and that a
simple radial alignment does not fully capture subhalo orientations. This
situation is not surprising, since our model only directly includes
information from the 1-halo regime (e.g. satellite positions and (sub)halo
orientations).

In Figure 6, we see a direct comparison of $\omega^{{}^{g}}_{cs}$ (the central
galaxy position, satellite galaxy orientation correlation function), as
measured directly from simulated (sub)halos, and two different HOD models: one
with alignment strengths found by fitting to the 1-halo regime only (labeled
as the 1-halo fit model), and one defined by fitting to the 2-halo regime only
(which we will call the 2-halo fit model). Overall, the 1-halo fit model
performs better, and is able to match $\omega^{\rm h}_{\rm cs}$ at all scales
to better than $20\%$, while the 2-halo fit model is only accurate at the
$\sim 40\%$ level. Neither fit suggests that our model can accurately capture
the full range of spatial scales, but looking at the 1-halo and 2-halo regimes
separately, we see the strengths and weaknesses of each model in more detail.
Of course, each model presents comparable better accuracy over the range of
scales used in the fit, but the 1-halo fit model generally performs much
better. In the 1-halo range of scales, the 1-halo fit model achieves $\sim
6\%$ accuracy while the 2-halo fit model only matches at the $60\%$ level. In
the 2-halo regime, the 1-halo fit model is only accurate at the $30\%$ level,
but even the 2-halo fit model performs no better than $20\%$.

For a clearer picture, the $\chi^{2}$ per degree of freedom ($\chi^{2}_{\rm
dof}$) (describing how well the galaxy correlation functions match those of
their halos) of the 1-halo fit model is $5.2$, $8.6$, and $3.8$ for the full,
1-halo, and 2-halo regimes respectively, while the 2-halo fit model had values
of $26160$, $55792$, and $10.8$. At this point, the only clear meaning of the
2-halo fit $\chi^{2}_{\rm dof}$ values is that fitting this model to the
2-halo regime does not work without extra information not yet captured by our
model. Surprisingly, the 2-halo regime has better $\chi^{2}_{\rm dof}$ in both
the 1-halo fit model and the 2-halo fit model, though we note that the 1-halo
regime has noticeably smaller errors, giving it a larger $\chi^{2}_{\rm dof}$
even with closer fractional agreement. We remind the reader that at this
point, we still have made no attempt to refine the radial alignment model, but
rather we have thus far only illustrated why we have chosen to fit our
parameters to the 1-halo regime alone. In Section 6.4, we turn attention to
fitting distance-dependent radial alignment parameters and refining the
goodness of fit.

Figure 6: Comparison of the central-satellite position-orientation correlation
function for halos, $\omega^{\rm h}_{\rm cs}$ (black points), with that of the
galaxies, $\omega^{\rm g}_{\rm cs}$ (colored lines). Two galaxy series are
shown, each one having been generated by an HOD model where the galaxy
alignment strengths were fit using MCMC. The series where the parameters were
fit to the 1-halo part of the halo correlation functions fits remarkably
closely in the 1-halo regime (agreeing within $6\%$), as expected, while
failing to accurately represent the 2-halo regime (differing by $30\%$).
Meanwhile, the HOD with alignment strengths obtained by fitting to the 2-halo
regime fail to model the 1-halo signal (differing by $60\%$), and don’t follow
the 2-halo correlation function well (only agreeing within $20\%$) despite
having been fit to that area.

We conclude this section by noting that our results imply a need for an
additional modeling ingredient that captures the effect of environment beyond
what is contained in the radial vector, similar to such ingredients that are
now commonly used to capture assembly bias in HODs (e.g., Hearin et al., 2016;
Yuan et al., 2018). The nature of this correction may require appeal to some
measure of large-scale density smoothed on some scale, and/or tidal
information. We acknowledge this limitation, however, as demonstrated here,
satellite alignments based on radial vectors are a simple, flexible, and
effective model when simulated subhalos are unavailable. The development of a
model with additional environmental dependence is beyond our current scope and
is left as a topic for future work.

###  6.3 Satellite Anisotropy

In §6.1-6.2, we restricted consideration to HOD-style models of satellite
alignments in which the spatial positions of satellites within their halos
were spherically symmetric. However, we know that satellite anisotropy affects
correlations (Samuroff et al., 2020), and for IA models of satellites based on
the radial vector, anisotropy in satellite locations has an influence on the
orientation correlation functions. And so in this section we study the effect
of satellite anisotropy on measurements of IA.

We model the anisotropic number density profile of satellites in halos as a
triaxial NFW profile following Jing & Suto (2002) and Schneider et al. (2012).

$\rho(R)=\frac{\rho_{c}}{\frac{R}{R_{s}}\left(1+\frac{R}{R_{s}}\right)^{2}}$
(20)

where the relation between the Cartesian and elliptical coordinates is given
by:

$\displaystyle x$
$\displaystyle=r\sin(\theta)\cos(\phi)=R\left(\frac{a}{c}\right)^{-\beta}\sin(\Theta)\cos(\Phi)$
$\displaystyle y$
$\displaystyle=r\sin(\theta)\sin(\phi)=R\left(\frac{b}{c}\right)^{\beta}\sin(\Theta)\sin(\Phi)$
$\displaystyle z$ $\displaystyle=r\cos(\theta)=R\cos(\Theta)$ (21)

where $a$, $b$, and $c$ are the principle axis lengths of the host halo
normalized such that $a=1$. The anisotropy parameter, $\beta$, controls the
magnitude of the anisotropy; when $\beta=0$, satellites follow a spherically
symmetric NFW profile; when $\beta=1$, satellites follow a triaxial NFW
profile with axis ratios equal to that of the underlying dark matter halo;
when $\beta>1$, the satellite distribution exhibits a greater degree of
anisotropy than the underlying halo.

We compare three levels of anisotropy in the spatial distribution of
satellites within their halos:

  * •

A subhalo position model, which places satellites at the locations of
subhalos.

  * •

A semi-isotropic model, which initially places galaxies using subhalo
positions, but then randomly rotates each host halo’s satellite population
about the halo major axis. This model therefore preserves correlations
involving the major axis but erases other anisotropies in the satellite
distribution.

  * •

An isotropic model, which places galaxies following a spherically symmetric
NFW profile (as in §6.1-6.2.

We stress that in all three models sketched above, the only difference between
the models is the level of anisotropy in the intra-halo positions of
satellites. For the model of orientations explored in this section, here we
assume that centrals are perfectly aligned with the halo shape, and for
satellites we adopt a radial alignment model with constant strength of $0.85$.

As shown in Figure 7, both the subhalo position and semi-isotropic models
allow the galaxy correlation functions to match those of the halos within
$10\%$ for the ED correlation functions and within $30\%$ for the EE
correlation function. Meanwhile, models with isotropic satellites predict much
weaker correlations, especially in the 1-halo regime, resulting in a $30\%$
and a $60\%$ difference from the subhalo ED and EE correlation functions,
respectively. The results illustrated in Figure 7 demonstrate that the level
of anisotropy of satellite positions within the halo comprises an important
degree of freedom for orientation correlation functions and that the major
axis of the anisotropy contains most of the relevant information.

Figure 7: The effect of subhalo anisotropy on EE and ED. Here we compare the
correlations of all of the halos themselves (black points) to those of all of
the galaxies (black lines), where the isotropy of galaxy placement is the only
difference between the three galaxy correlation functions shown. The galaxy
correlations shown are the average of fifty iterations (assigning position and
alignment fifty times). While all three models align galaxies with respect to
the radial vector, the subhalo position model uses subhalo positions to place
galaxies while the isotropic model places the galaxies isotropically, and the
semi-isotropic model uses the subhalo positions then rotates them around the
halo major axis to preserve the radial distribution. Here we see that using
the subhalo positions, even if we rotate them, produces a signal remarkably
close to that of the subhalos. Placing galaxies isotropically appears to lose
some information required to match the halo correlation functions with the
parameters given.

###  6.4 Capturing Subhalo IA in Gravity-only Simulations

The previous sections have dealt only with the position-orientation
correlation functions; in this section we consider additional orientation
correlation functions to study the flexibility of our satellite IA models.
Figure 5 shows how three different correlation functions are impacted by
satellite IA models with either constant or radially-dependent alignment
strength. Each prediction was computed from 50 realizations of the model, and
in all cases galaxy positions are placed at the center of (sub)halos. We
remind the reader that the satellite fraction is $\sim 25\%$ in the predicted
galaxy population, and so our plots show central-to-satellite correlation
functions of each type: the central position to satellite position correlation
function, $\xi,$ the central position to satellite orientation correlation
function, $\omega$, and the central orientation to satellite orientation
correlation function, $\eta$. We plot central-satellite cross-correlation
functions to isolate the effect of satellite alignments on the IA signal.

The orientation correlation function $\eta$ is better fit in the 1-halo
regime, again indicating that additional environmental information beyond
simple radial alignment is needed in order to capture physically plausible
satellite effects on IA in the 2-halo regime. Fitting to the 1-halo regime, we
found the best constant radial alignment strength to be $0.826$, and the best
distance-dependent alignment strength, following equation (19), to have
$\text{a}=0.804$ and $\gamma=-0.04$.

To derive these values, we can directly fit the Dimroth-Watson distribution to
the PDF of misalignment angles in the simulation. The results of this fitting
process are shown in Figure 8, which displays the best-fitting Dimroth-Watson
curve to the distribution of misalignment angles of simulated subhalos. To
study the distance-dependence of alignment strength, we bin the satellite
galaxies by $r/r_{vir}$, and perform separate fits on each bin. In the right
panel of Figure 8, we show a power law fit to the resulting relationship
between $\mu$ and $r/r_{vir}$.

Figure 8: Left: Differential distribution of (sub)halo radial misalignment
angles for extant (sub)halos with a peak mass greater than
$10^{12}\leavevmode\nobreak\ h^{-1}\leavevmode\nobreak\ h^{-1}M{\odot}$. The
best fit Watson distribution is shown as a red line with an alignment
strength, $\mu=0.83$. The shaded region shows the variation in the alignment
strength for 90% of subhalos. Right: Dependence of the alignment strength for
subhalos as a function of radial position, scaled by the virial radius of the
host halo. The red line is a clipped power law fit.

At first glance, both the model with distance-dependent alignment strength and
that with the constant alignment strength appear to give equally valid results
for the $\omega$ correlation function, as they match within $20\%$ overall and
within $6\%$ and $8\%$, respectively, in the 1-halo regime. However, there are
noticeable differences when the full $\chi^{2}_{\rm dof}$ is considered.
Looking at the $\omega$ (ED) correlation function, the model with distance-
dependent alignment strength has a $\chi^{2}_{\rm dof}$ of $4.3$ overall,
$3.9$ in the 1-halo regime, and $4.0$ in the 2-halo regime while the model
with constant alignment strength shows a $\chi^{2}_{\rm dof}$ of $6.3$, $8.2$,
and $4.3$ for the overall, 1-halo, and 2-halo regimes respectively. Similarly,
if we look at the $\eta$ (EE) correlation functions, we see that the distance-
dependent alignment strength matches the halo correlation functions within
$10\%$ overall and within $13\%$ in the 1-halo regime, while the constant
alignment strength model only matches within $24\%$ overall and $15\%$ in the
1-halo regime. From this, we see that there does seem to be minor improvement
from using a distance-dependent alignment model, especially when we limit our
investigation to the 1-halo regime where our model works best. However, we
should note the small value for $\gamma$ found above ($\gamma=-0.04$),
indicating this power law dependence is very close to a constant radial
alignment.

In the 2-halo regime, the radial alignment model fails to capture the
environmental dependence of subhalo orientations. In order to capture
physically plausible satellite contributions to the IA signal on large scales,
improving the physical realism beyond a simple radial alignment model will be
necessary. We will explore such modeling improvements in future work.

##  7 Capturing IA with Realistic Complexity: testing the model with Illustris

In the previous sections, we studied how models of satellite IA are able to
capture the orientation correlation functions of subhalos in gravity-only
N-body simulations. We now turn attention in this section to the IA signal
exhibited by galaxies in the IllustrisTNG hydrodynamical simulation. In so
doing, we will study IA models of both central and satellite galaxies acting
together in concert. For central IA, we use the shape of the parent halo as
the reference vector with respect to which central galaxy orientations are
(mis)aligned. For satellite IA, we consider two different types of model: one
based on the orientation of subhalos, the second HOD-type model based on
radial alignments of constant strength. The free parameters of our model
modulate the strength of central (satellite) galaxies with respect to the
parent (sub)halo, and for all our results we use emcee to fit the parameters
regulating the alignment strength.

Prior to fitting the galaxy alignment strengths, we first adjust the
occupation model parameters. We used the Zheng07Cens and Zheng07Sats for the
central and satellite galaxy occupation models available in halotools as
described in Zheng et al. (2007). Table 3 shows the occupation parameters
found by fitting the central and satellite occupation functions used by
Zheng07Cens and Zheng07Sats to the mean occupations seen in Illustris.
Specifically, the parameters were fit on the region where the central
occupation transitions from zero central galaxies, to one central galaxy. For
simplicity, in our fits $\sigma_{\log{M}}$ and $\alpha$ were held fixed at
$0.26$ and $1.0$ respectively and not allowed to vary.

The corner plots in Figure 9 display the posteriors from our fits to the
orientation correlation functions in Illustris. As indicated in the legend, we
show results for fits to three separate galaxy samples, each defined by a
different stellar mass cut. Sample 1 corresponds to $\log{(M_{*})}\geq 10.5$,
sample 2 to $\log{(M_{*})}\geq 10$, and sample 3 to $\log{(M_{*})}\geq 9.5$.
Two broad trends are apparent in Figure 9. First, the central galaxies have a
noticeably larger alignment strength than the satellite galaxies. Second, as
the stellar mass threshold increases, so does the alignment strength of both
populations. Both of the above trends can be more clearly seen in Figure 10,
which shows how the best-fitting alignment strength $\mu$ changes with the
stellar mass cut.

(a) Sample One

(b) Sample Two

(c) Sample Three

Figure 9: The parameter values found from MCMC for central alignment strength
($\mu_{rmcen}$) and satellite alignment strength ($\mu_{\rm sat}$) to let HOD
models best fit Illustris correlation functions using a subhalo alignment
model for the satellite galaxies. Each sample corresponds to a different mass
cutoff for the halos included in the HOD model. These are the same mass
cutoffs seen in Figure 11, namely $\log{M_{*}}\geq 10.5$ for sample one,
$\log{M_{*}}\geq 10.0$ for sample two, and $\log{M_{*}}\geq 9.5$ for sample
three. In the models used to generate these corner plots, we used a central
alignment for central galaxies and subhalo alignment for satellite galaxies.
We did the same process using constant radial alignment for satellite
galaxies, yielding similar corner plots not shown here. As a general trend, we
see the alignment strengths increase with higher mass cutoffs. Left: Best
parameters found were $\mu_{\rm cen}=0.802$ and $\mu_{\rm sat}=0.449$. Middle:
Best parameters found were $\mu_{\rm cen}=0.683$ and $\mu_{\rm sat}=0.274$.
Right: Best parameters found were $\mu_{\rm cen}=0.528$ and $\mu_{\rm
sat}=0.112$.

As seen in Figure 11, our best-fitting HOD model generates correlation
functions that agree reasonably well with Illustris-TNG. It is by now well
established that an HOD model more complex than the one we use here is
required for precision fits to the spatial clustering of hydro-simulated
galaxies (e.g., Beltz-Mohrmann et al., 2020; Hadzhiyska et al., 2020). For the
$\xi$ correlations functions, the HOD models for samples one, two, and three
produce correlations within $14.0\%$, $9.5\%$, and $12.0\%$ respectively to
their corresponding counterpart from the Illustris data. This level of
agreement is a sufficient baseline for purposes of broadly assessing the
accuracy of our alignment strength models.

The center and right panels of Figure 11 display the level of success of our
best-fitting IA model to capture the position-orientation and orientation-
orientation correlation functions, respectively. Table 2 shows the
$\chi^{2}_{\rm dof}$ values for each of these fits. For the satellite IA, we
find that a model based on radial alignment performs comparably well to a
model based on subhalo alignment, which is encouraging for the prospects of
HOD-type approaches to capturing orientation correlation functions.

Satellite | $\omega$ | $\eta$ | Central | Satellite
---|---|---|---|---
Alignment |  |  | Alignment | Alignment
Type |  |  | Strength | Strength
Radial (Sample One) | 1.9 | 1.3 | 0.780 | 0.335
Radial (Sample Two) | 2.1 | 1.2 | 0.685 | 0.139
Radial (Sample Three) | 1.4 | 0.5 | 0.545 | 0.030
Subhalo (Sample One) | 1.8 | 1.2 | 0.802 | 0.449
Subhalo (Sample Two) | 1.7 | 1.1 | 0.683 | 0.274
Subhalo (Sample Three) | 1.8 | 0.7 | 0.528 | 0.112
Table 2: $\chi^{2}_{\rm dof}$ values for the central-satellite $\omega$
(position-orientation) and $\eta$ (orientation-orientation) correlation
functions comparing the functions from our HOD models with those from
Illustris. Also included are the alignment strengths that best fit our HOD
model to Illustris-TNG. As before, these are the same mass cutoffs seen in
Figure 11, namely $\log{M_{*}}\geq 10.5$ for sample one, $\log{M_{*}}\geq
10.0$ for sample two, and $\log{M_{*}}\geq 9.5$ for sample three. Figure 10:
Alignment using the central alignment and subhalo alignment models for the
central and satellite galaxies respectively. The three points shown for each
type of galaxy correspond to samples one, two, and three. The samples are
separated based on the minimum mass threshold for the galaxies included in the
sample, shown on the x-axis. Figure 11: Left: Two-point correlation functions
offset by 1-dex for clarity Middle: ellipticity-direction (ED) correlation
function Right: ellipticity-ellipticity (EE) correlation function. In each
panel, the points with error bars are measurements made on the Illustris
TNG300-1 simulation with error bars estimated using jackknife re-sampling of
the box. The lines with shaded regions are halo model predictions made by
populating a DMO simulation with mock galaxies where the shaded region shows
the variation from random realizations of the model. The three colors are for
three stellar mass threshold samples.

##  8 Discussion & Conclusions

In this work, we have studied forward modeling techniques that enable
predictions of large-scale structure that include intrinsic alignments in
galaxy shape. Our methodology is fully simulation-based: when predicting
galaxy correlation functions, we first create a synthetic universe of galaxies
by populating simulated halos, and then we compute summary statistics with the
same point estimators used in the corresponding observational measurements.
Our prediction pipeline is computationally efficient, and it has the practical
capability to run likelihood analyses in the deeply nonlinear regime with MCMC
sampling algorithms. We have demonstrated the effectiveness of our model by
reproducing IA signals exhibited by subhalos in gravity-only simulations, as
well as by galaxies in hydrodynamical simulations such as Illustris-TNG.

Our simulation-based methodology is flexible, and it is able to leverage
information from subhalos when available in a high-resolution simulation, or
alternatively to populate catalogs of host halos in lower-resolution
simulations. When subhalo information is available in the underlying
simulation, we can generate more physically realistic predictions by using the
orientation of subhalos as the reference vector for satellite galaxies. When
using simulations that do not include subhalo orientations (as is often the
case due to the computational demands of survey-scale simulations), a radial
satellite alignment model (in which satellites preferentially align with
respect to the radial vector between central and satellite) can approximate
reasonably well the IA signal exhibited by hydro-simulated galaxies in the
1-halo regime.

In detail, we find that a distance-dependent alignment strength, following
equation (19), may provide a slightly closer match to the 1-halo term signal
seen in the correlation functions of either gravity-only subhalos or hydro-
simulated galaxies. However, even a satellite model with a simple constant
radial alignment strength produces a reasonable correlation function on small
scales; looking specifically at $\omega^{\rm g}_{\rm cs}$, we see that the
distance-dependent alignment strength with the radial alignment model gives a
$\chi^{2}_{\rm dof}$ of 4.3 overall, while the constant alignment strength
gives 6.3.

As shown in Figure 5, no form of radial alignment model is able to faithfully
reproduce the satellite correlation functions in the 2-halo regime; this
implies the need for the development of a new model of galaxy alignments that
more flexibly incorporates environmental information, analogous to ongoing
efforts to capture assembly-biased halo occupation statistics (e.g., Contreras
et al., 2021; Yuan et al., 2022; Lange et al., 2023). Reference vectors based
in part on the tidal field are a promising way to incorporate such
environmental effects (Harnois-Déraps et al., 2022).

We find that the anisotropic spatial distribution of satellites within their
host halos plays an important role in the strength of orientation correlation
functions. Figure 7 shows that models based on the anisotropic positions of
subhalos substantially improve upon the realism of the predictions. Meanwhile,
models based on a spherically symmetric NFW profile perform poorly in the
1-halo regime. These results imply the need for IA models with improved
sophistication in the intra-halo phase space distributions of satellites,
echoing the well-established need for such ingredients to capture ordinary
two-point clustering (e.g., Orsi & Angulo, 2018; Hadzhiyska et al., 2022).

One of our main goals was to test how well the IA signal in Illustris-TNG can
be mimicked with a simple HOD-type model of galaxy orientations. In our tests,
we used MCMC to fit our models of central and satellite alignment strength to
the orientation correlation functions of Illustris-TNG galaxies, repeating
this analysis separately for samples defined by three different stellar mass
cuts. Figure 11 summarizes the results of this analysis, and Table 2 records
the alignment strengths of our best-fitting models. Subhalo-based models are
particularly effective in their ability to generate a mock galaxy catalog with
correlation functions similar to Illustris-TNG, attaining $\chi^{2}_{\rm dof}$
close to unity for every sample and correlation function considered. Satellite
models based on constant radial alignment perform comparably well to the
subhalo models, despite their relative simplicity. In future work, we will
generalize the HOD-type methodology used here to incorporate joint
correlations between IA strength and galaxy brightness, color, and redshift as
seen in simulations (Zjupa et al., 2022) and observations (Samuroff et al.,
2023). We will also work towards adding more sophisticated alignment models
that take into account more detailed sets of halo properties and other
information, including redshift, in order to produce more robust results (Xu
et al., 2023).

In a full forward model of galaxy correlation functions $\xi,\omega,$ and
$\eta$, additional ingredients for the probabilistic distribution of galaxy
shapes will be needed in addition to the alignment distribution model studied
here. Our mock-population framework lays the groundwork for carrying out
Bayesian inference with simulation-based forward models of galaxy shape
correlation functions. Our simulation-based approach simplifies the effort
required to model systematic effects that can be challenging to capture
analytically, such as galaxy assembly bias, halo exclusion in the transition
between 1- and 2-halo regimes (van den Bosch et al., 2013; García & Rozo,
2019), and beyond-linear halo bias (Mahony et al., 2022; Mead & Verde, 2021).
Forward models based on mock-population methods also extend naturally to
predict higher-order summary statistics, which can provide additional
information not contained in the two-point functions studied here (Harnois-
Déraps et al., 2021). With the additional effort outlined above, our models
could be used to derive posteriors on the true intrinsic alignment strength of
galaxies from observations on highly nonlinear scales, and to supply priors on
the IA nuisance parameters used in cosmological analyses of larger scales.

## Acknowledgements

This paper has undergone internal review in the LSST Dark Energy Science
Collaboration. The internal reviewers were Joachim Harnois-Déraps, Eve Kovacs,
and Contance Mahony; we thank them for helpful comments and discussion.

Author contributions: NVA wrote code, conducted the primary analysis, and
wrote the main body of the text. DC conceived of and initiated the project and
contributed significantly to the initial text. JB participated in initial
project planning, contributed significant guidance and mentorship during the
code creation and analysis stages, and authored portions of the text. CDL
contributed significant guidance and mentorship during the code creation and
analysis stages, as well as writing and revising portions of the draft. FL
contributed guidance on code, model, and statistical aspects of the analysis.
APH contributed to code development of the alignment model and summary
statistics computations, and helped develop and revise the draft. RM provided
significant mentorship and guidance during the project’s initial stages and
provided input on the paper’s near-complete version.

NVA and JB are partially supported by NSF Award AST-2206563. RM was supported
in part by a grant from the Simons Foundation (Simons Investigator in
Astrophysics, Award ID 620789) and in part by the US Department of Energy
grant DE-SC0010118.

This work was completed in part using the Discovery cluster, supported by
Northeastern University’s Research Computing team. This research used
resources of the National Energy Research Scientific Computing Center (NERSC),
a U.S. Department of Energy Office of Science User Facility located at
Lawrence Berkeley National Laboratory, operated under Contract No. DE-
AC02-05CH11231.

The DESC acknowledges ongoing support from the Institut National de Physique
Nucléaire et de Physique des Particules in France; the Science & Technology
Facilities Council in the United Kingdom; and the Department of Energy, the
National Science Foundation, and the LSST Corporation in the United States.
DESC uses resources of the IN2P3 Computing Center (CC-IN2P3–Lyon/Villeurbanne
- France) funded by the Centre National de la Recherche Scientifique; the
National Energy Research Scientific Computing Center, a DOE Office of Science
User Facility supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, funded
by UK BEIS National E-infrastructure capital grants; and the UK particle
physics grid, supported by the GridPP Collaboration. This work was performed
in part under DOE Contract DE-AC02-76SF00515.

This work made use of Astropy:444http://www.astropy.org a community-developed
core Python package and an ecosystem of tools and resources for astronomy
(Astropy Collaboration et al., 2013, 2018, 2022).

##  Appendix A Galaxy & Halo Shapes

We model (sub)halos and galaxies as 3-dimensional ellipsoids. The shape and
orientation of (sub)halos/galaxies may then be characterized by calculating
the reduced inertia tensor for the particle distribution. We define the
reduced inertia tensor as:

$\tilde{\bf I}_{ij}=\frac{\sum m_{n}\frac{x_{ni}x_{nj}}{r_{n}^{2}}}{\sum
m_{n}}$ (22)

where

$r_{n}^{2}=\sum x_{ni}^{2}$ (23)

is the distance between the center of mass and the $n^{\rm th}$ particle in
the system. The reduced inertia tensor applies more weight to particles that
are near the center of the object, reducing the sensitivity to loosely bound
particles present in the outer regions of the object.

The triaxial shape of the (sub)halo/galaxy is specified by the eigenvalues of
the inertia tensor, $\lambda_{a}>\lambda_{b}>\lambda_{c}$, where the half-
lengths of the principle axis of the ellipsoid are given by
$a=\sqrt{\lambda_{a}}$, $b=\sqrt{\lambda_{b}}$, $c=\sqrt{\lambda_{c}}$,
($a\geq b\geq c$). The orientation of this ellipsoid is then specified by the
eigenvectors, $\hat{e}_{a}$, $\hat{e}_{b}$, and $\hat{e}_{c}$.

The radial weighting applied in equation (22) tends to make
(sub)halos/galaxies appear more spherical. To alleviate this, we iteratively
calculate the reduced inertia tensor for (sub)halos/galaxies. After the
initial calculation of the $\tilde{\bf I}$, the particle distribution is
rotated so that $\hat{e}_{a}$, $\hat{e}_{b}$, and $\hat{e}_{c}$ are aligned
with the x-, y-, x-axis. The radial distance, $r_{n}^{2}$ in equation (22), is
replaced with the elliptical radial distance:

$r_{n}^{2}=\left(\frac{x_{n}}{a}\right)^{2}+\left(\frac{y_{n}}{b}\right)^{2}+\left(\frac{z_{n}}{c}\right)^{2}$
(24)

This process is repeated until the eigenvalues change by less than 1% between
iterations.

We implicitly model (sub)halos/galaxies as homologous triaxial ellipsoids.
Note the caveats here. halos contain substructure. galaxies may contain
multiple components. galaxies are observed to have twisting isophotes. By
modeling galaxies as homologous ellipsoids, we can specify a single shape and
orientation, vastly simplifying the modeling.

##  Appendix B Fitting Functions

In Figures 3 and 4, we fit a smooth function to our measurements of simulated
(subhalo) EE, $\eta(r)$, measurements. We primarily do this to reduce the
noise in the measurement at any given radius in order to make the figures more
clear. We present our fitting function here for completeness.

$\eta(r)=A_{1}\exp\left(-(r/R_{1})^{\gamma}\right)+\frac{A_{2}\exp\left(-(R_{2}/r)\right)}{1.0+(r/R_{2})^{\alpha}}$
(25)

$\displaystyle\omega(r)=$ $\displaystyle
A_{1}\exp\left(-(r/R_{1})^{\gamma}\right)+\frac{A_{2}\exp\left(-(R_{2}/r)\right)}{1.0+(r/R_{2})^{\alpha}}$
(26) $\displaystyle+\frac{A_{3}}{1.0+(r/R_{3})^{\beta}}$ (27)

##  Appendix C Illustris-TNG HOD

As shown in §7, the alignment correlation functions simultaneously depend on
the misalignment distributions as well as the underlying halo occupation
statistics. Thus in order to test the efficacy of our misalignment modeling
using Illustris-TNG, it’s important to first establish an accurate baseline
HOD model. In Table 3, we see the occupation model parameters used. These
values were found by fitting the occupation parameters in the Zheng07Cens and
Zheng07Sats occupation models to fit the occupation distribution of our HOD
models with that of Illustris. Specifically, we fit the parameters in the
region where the central occupation transitions from zero to one. Note that
the values $\sigma_{\log{M}}$ and $\alpha$ have been fixed for simplicity,
only allowing the other three parameters to vary.

$\log(M_{\rm thresh})$ | $\log(M_{\rm min})$ | $\sigma_{\log{M}}$ | $\log(M_{0})$ | $\log(M_{1})$ | $\alpha$
---|---|---|---|---|---
9.0 | 11.37 | 0.26 | 11.55 | 12.35 | 1.0
9.5 | 11.61 | 0.26 | 11.8 | 12.6 | 1.0
10.0 | 11.93 | 0.26 | 12.05 | 12.85 | 1.0
10.5 | 12.54 | 0.26 | 12.68 | 13.48 | 1.0
Table 3: HOD parameters for TNG300-1

##  Appendix D Covariance

Due to the probabilistic nature of populating an HOD model and aligning with a
Dimroth-Watson distribution (discussed later in Section 4), there is
stochasticity in both the data vector and the model. We have two different
ways in which we select our covariance, depending on our data vector.

In Section 6, we fit our satellite galaxy alignment strengths to find an
alignment strength for the radial alignment such that the central position,
satellite orientation correlation functions for the galaxies match as closely
as possible to those of the halos (keeping central galaxies perfectly aligned
with their host halos for simplicity). However, because there is an element of
stochasticity both in the population of the dark matter halos with galaxies,
and in the realignment of the galaxies, both our data vector and our model
have covariance. To account for this, we look at the contribution to the
covariance of repopulating halos versus that of assigning galaxy alignments.
As seen in Figure 12, we find that most of the covariance comes from aligning
galaxies. As such, we use a fixed seed for the galaxy population step in our
MCMC, meaning the same dark matter halos are populated on each run, keeping
the galaxy alignment step the only stochastic process. For our covariance, we
use the covariance of $\omega^{\rm g}_{\rm cs}$ (the central galaxy position,
satellite galaxy orientation correlation function) from 100 instances of
realigning galaxies at that seed, for a given alignment strength (see Section
6 for more detail). We use this covariance in the MCMC, after which a new
alignment strength is determined, and the process repeats (generating a new
covariance from 100 new iterations at the new alignment strength) until the
alignment strength resulting from the MCMC does not change appreciably.

Figure 12: The diagonals of the covariance matrices generated from measuring
the ed correlations of 100 instances of mock galaxy catalogs. The curve for
total covariance comes from catalogs where we performed a full repopulation
(i.e. no random seed was used and each time, halos were selected and populated
randomly before assigning galaxy orientation). The orientation curve was
generated in a similar fashion, but instead of allowing a full repopulation,
we used a fixed seed for the halo population (i.e. the same halos were chosen
each time with the same galaxy populations) meaning the only source of
variance was the galaxy alignments. The final curve shows the difference
between these two, Total - Orientation, representing the covariance that comes
from just the repopulation of the halos, ignoring the galaxy alignments. Here
we see that the majority of the covariance comes from realigning the galaxies.

In Section 7, rather than fitting to the halo orientations, we found both
central and satellite galaxy alignment strengths by fitting $\omega^{\rm
g}_{\rm gg}$ (the gg subscript now referring to all galaxies) correlation
functions to those of Illustris-TNG300. In this case, since we had a fixed
data vector, we used jackknife resampling of the Illustris data, splitting the
simulation box into 5x5x5 jackknife regions to determine our covariance. It is
important to note that we multiplied this jackknife covariance by the Hartlap
factor, shown in equation (28) (taken from equation 17 in Hartlap et al.
(2007))

$\hat{C}^{-1}=\frac{n-p-2}{n-1}\hat{C}^{-1}_{*},$ (28)

where $\hat{C}^{-1}$ is the adjusted inverse covariance, $\hat{C}^{-1}_{*}$ is
the inverse covariance directly from jackknife resampling, $n$ is the number
of independent observations (i.e. the number of jackknife regions), and $p$ is
the number of elements in the data vector (i.e. the number of bins). This
helps approximately debias our estimation of the original inverse jackknife
covariance.

#### Data availability

We used the Bolshoi-Planck catalog and Multidark, both available through
halotools and online at https://www.cosmosim.org/. The Illustris-TNG data we
used is publicly available at https://www.tng-project.org/.

## References

  * Akeson et al. (2019) Akeson R., et al., 2019, The Wide Field Infrared Survey Telescope: 100 Hubbles for the 2020s (arXiv:1902.05569)
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Behroozi et al. (2013a) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109
  * Behroozi et al. (2013b) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013b, ApJ, 763, 18
  * Beltz-Mohrmann et al. (2020) Beltz-Mohrmann G. D., Berlind A. A., Szewciw A. O., 2020, MNRAS, 491, 5771
  * Berlind et al. (2003) Berlind A. A., et al., 2003, ApJ, 593, 1
  * Blazek et al. (2011) Blazek J., McQuinn M., Seljak U., 2011, Journal of Cosmology and Astroparticle Physics, 2011, 010
  * 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. A., Fang X., 2019, Phys. Rev. D, 100, 103506
  * Bridle & King (2007) Bridle S., King L., 2007, New Journal of Physics, 9, 444
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Campbell et al. (2018) Campbell D., van den Bosch F. C., Padmanabhan N., Mao Y.-Y., Zentner A. R., Lange J. U., Jiang F., Villarreal A. S., 2018, MNRAS, 477, 359
  * Chisari et al. (2019) Chisari N. E., et al., 2019, ApJS, 242, 2
  * Contreras et al. (2021) Contreras S., Angulo R. E., Zennaro M., 2021, MNRAS, 504, 5205
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fortuna et al. (2020) Fortuna M. C., Hoekstra H., Joachimi B., Johnston H., Chisari N. E., Georgiou C., Mahony C., 2020, Monthly Notices of the Royal Astronomical Society, 501, 2983
  * García & Rozo (2019) García R., Rozo E., 2019, MNRAS, 489, 4170
  * Hadzhiyska et al. (2020) Hadzhiyska B., Bose S., Eisenstein D., Hernquist L., Spergel D. N., 2020, MNRAS, 493, 5506
  * Hadzhiyska et al. (2022) Hadzhiyska B., et al., 2022, arXiv e-prints, p. arXiv:2210.10068
  * Harnois-Déraps et al. (2021) Harnois-Déraps J., Martinet N., Castro T., Dolag K., Giblin B., Heymans C., Hildebrandt H., Xia Q., 2021, MNRAS, 506, 1623
  * Harnois-Déraps et al. (2022) Harnois-Déraps J., Martinet N., Reischke R., 2022, MNRAS, 509, 3868
  * Hartlap et al. (2007) Hartlap Simon, P. Schneider, P. 2007, A&A, 464, 399
  * Hearin et al. (2016) Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D., Tollerud E., 2016, MNRAS, 460, 2552
  * Hearin et al. (2017) Hearin A. P., et al., 2017, The Astronomical Journal, 154, 190
  * 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
  * Hirata & Seljak (2004) Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526
  * Hoffmann et al. (2022) Hoffmann K., et al., 2022, Phys. Rev. D, 106, 123510
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Jagvaral et al. (2022) Jagvaral Y., Lanusse F., Singh S., Mandelbaum R., Ravanbakhsh S., Campbell D., 2022, Monthly Notices of the Royal Astronomical Society, 516, 2406
  * Jiang & van den Bosch (2016) Jiang F., van den Bosch F. C., 2016, MNRAS, 458, 2848
  * Jing & Suto (2002) Jing Y. P., Suto Y., 2002, The Astrophysical Journal, 574, 538
  * Joachimi et al. (2013) Joachimi B., Semboloni E., Hilbert S., Bett P. E., Hartlap J., Hoekstra H., Schneider P., 2013, Monthly Notices of the Royal Astronomical Society, 436, 819
  * Klypin et al. (2011) Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Knebe et al. (2011) Knebe A., et al., 2011, MNRAS, 415, 2293
  * Knebe et al. (2013) Knebe A., et al., 2013, MNRAS, 435, 1618
  * Krause et al. (2015) Krause E., Eifler T., Blazek J., 2015, Monthly Notices of the Royal Astronomical Society, 456, 207
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, ApJ, 412, 64
  * Lange et al. (2023) Lange J. U., Hearin A. P., Leauthaud A., van den Bosch F. C., Xhakaj E., Guo H., Wechsler R. H., DeRose J., 2023, MNRAS, 520, 5373
  * Lee (2019) Lee J., 2019, ApJ, 872, 37
  * Mahony et al. (2022) Mahony C., et al., 2022, Monthly Notices of the Royal Astronomical Society, 515, 2612
  * Mandelbaum et al. (2006) Mandelbaum R., Hirata C. M., Ishak M., Seljak U., Brinkmann J., 2006, MNRAS, 367, 611
  * Mansfield & Avestruz (2021) Mansfield P., Avestruz C., 2021, MNRAS, 500, 3309
  * Marinacci et al. (2018) Marinacci F., et al., 2018, Monthly Notices of the Royal Astronomical Society, 480, 5113
  * Mead & Verde (2021) Mead A. J., Verde L., 2021, Monthly Notices of the Royal Astronomical Society, 503, 3095
  * Mead et al. (2015) Mead A. J., Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015, MNRAS, 454, 1958
  * Naiman et al. (2018) Naiman J. P., et al., 2018, Monthly Notices of the Royal Astronomical Society, 477, 1206
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Nelson et al. (2017) Nelson D., et al., 2017, Monthly Notices of the Royal Astronomical Society, 475, 624
  * Orsi & Angulo (2018) Orsi Á. A., Angulo R. E., 2018, MNRAS, 475, 2530
  * Pillepich et al. (2017) Pillepich A., et al., 2017, Monthly Notices of the Royal Astronomical Society, 475, 648
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla A., Behroozi P., Primack J., Klypin A., Lee C., Hellinger D., 2016, Monthly Notices of the Royal Astronomical Society, 462, 893
  * Samuroff et al. (2020) Samuroff S., Mandelbaum R., Di Matteo T., 2020, MNRAS, 491, 5330
  * Samuroff et al. (2023) Samuroff S., et al., 2023, MNRAS,
  * Scaramella et al. (2022) Scaramella R., et al., 2022, A&A, 662, A112
  * Schneider & Bridle (2010) Schneider M. D., Bridle S., 2010, Monthly Notices of the Royal Astronomical Society, 402, 2127
  * Schneider et al. (2012) Schneider M. D., Frenk C. S., Cole S., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 030
  * 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
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2017) Springel V., et al., 2017, Monthly Notices of the Royal Astronomical Society, 475, 676
  * Sra & Karp (2013) Sra S., Karp D., 2013, Journal of Multivariate Analysis, 114, 256
  * Van Alfen et al. (2023) Van Alfen N., Campbell D., Hearin A., 2023, The Journal of Open Source Software, in preparation
  * Watson (1965) Watson G. S., 1965, Biometrika, 52, 193
  * Xu et al. (2023) Xu K., Jing Y. P., Zhao D., 2023, Towards Physical Understanding of Galaxy-Halo Alignment (arXiv:2307.12334)
  * Yang et al. (2003) Yang X., Mo H. J., van den Bosch F. C., 2003, MNRAS, 339, 1057
  * Yuan et al. (2018) Yuan S., Eisenstein D. J., Garrison L. H., 2018, MNRAS, 478, 2019
  * Yuan et al. (2022) Yuan S., Garrison L. H., Hadzhiyska B., Bose S., Eisenstein D. J., 2022, MNRAS, 510, 3301
  * Zheng et al. (2007) Zheng Z., Coil A. L., Zehavi I., 2007, The Astrophysical Journal, 667, 760
  * Zjupa et al. (2022) Zjupa J., Schäfer B. M., Hahn O., 2022, MNRAS, 514, 2049
  * van den Bosch et al. (2013) van den Bosch F. C., More S., Cacciato M., Mo H., Yang X., 2013, MNRAS, 430, 725
  * van den Bosch et al. (2018) van den Bosch F. C., Ogiya G., Hahn O., Burkert A., 2018, MNRAS, 474, 3043