# Alignments Between Galaxies and the Cosmic Web at $z\sim 1-2$ in the
IllustrisTNG Simulations

Benjamin Zhang Department of Physics and Astronomy, University of Southern
California, Los Angeles, CA 90089, USA<EMAIL_ADDRESS>Khee-Gan Lee Kavli
IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Alex Krolewski Perimeter Institute for Theoretical Physics, 31 Caroline St.
North, Waterloo, ON N2L 2Y5, Canada  Waterloo Centre for Astrophysics,
University of Waterloo, Waterloo, ON N2L 3G1, Canada  Jingjing Shi Kavli IPMU
(WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
Benjamin Horowitz Princeton Department of Astrophysical Sciences, Princeton,
NJ 08544, USA  Robin Kooistra Kavli IPMU (WPI), UTIAS, The University of
Tokyo, Kashiwa, Chiba 277-8583, Japan

###### Abstract

Galaxy formation theories predict that galaxy shapes and angular momenta have
non-random alignments with the cosmic web. This leads to so-called intrinsic
alignment between pairs of galaxies, which is important to quantify as a
nuisance parameter for weak lensing. We study galaxy-cosmic web alignment in
the IllustrisTNG suite of hydrodynamical simulations at redshifts 1 and 2,
finding that alignment trends are consistent with previous studies. However,
we find that the magnitude of the spin alignment signal is $\sim 2.4\times$
weaker than seen in previous studies of the Horizon-AGN simulation, suggesting
that this signal may have significant dependence on subgrid physics. Based on
IllustrisTNG, we then construct mock observational spectroscopic surveys that
can probe shape-cosmic web alignment at $z\sim 1-2$, modeled on the low-$z$
galaxy redshift and IGM tomography surveys on the upcoming Subaru Prime Focus
Spectrograph Galaxy Evolution (PFS GE) survey. However, even over box sizes of
$L=205\,h^{-1}\,\mathrm{Mpc}$, we find that global anisotropies induce a
sample variance in the 2D projected alignment signal that depend on the
projected direction — this induces significant errors in the observed
alignment. We predict a $5.3\sigma$ detection of IllustrisTNG’s shape
alignment signal at $z\sim 1$ from Subaru PFS GE, although a detection would
be challenging at $z\sim 2$. However, a rough rescaling of the relative
alignment signal strengths between the TNG and HorizonAGN simulations suggests
that PFS GE should be able to more easily constrain the latter’s stronger
signal.

Cosmic web (330), Galaxy formation (595), Hydrodynamical simulations (767),
Intergalactic medium (813), Redshift surveys (1378)

††journal: ApJ

tablenum SIXtablenum

##  1 Introduction

In the prevailing paradigm of inflationary Big Bang cosmology, tiny
overdensities in the initial Gaussian density field grow via gravity,
eventually reaching a point at which they decouple from the Hubble flow.
Collapsing under their own gravity, virialized structures referred to as dark
matter halos (Gunn & Gott, 1972) grow via accretion from the surrounding
large-scale structure: smaller halos assemble before merging into larger ones
in a hierarchical fashion (Blumenthal et al., 1984; Wechsler et al., 2002). As
dark matter collapses around halos, on larger scales ($\gtrsim$ Mpc) it forms
a distinctive “cosmic web,” with filaments connecting halos and sheets
surrounding the large empty voids (Zeldovich et al., 1982; Klypin & Shandarin,
1983; Einasto et al., 1984; Davis et al., 1985; Geller & Huchra, 1989; Bond et
al., 1996). The cosmic web is intimately related to the dynamics of dark
matter: gravitational collapse first occurs onto the plane of sheets, then
radially perpendicular to the long axis of filaments, and finally along all
three directions into nodes or halos (Zel’Dovich, 1970). Matter thus moves out
of voids, along filaments and sheets, and ultimately into halos.

Halo shapes and angular momenta are determined by these accretion processes,
linking halo properties to the properties of the surrounding dark matter. On
small ($\sim$ Mpc) scales, the connection between the cosmic web and halo
spins and shapes is in principle a sensitive probe of halo formation and
evolution. This connection also gives rise to so-called intrinsic alignments
between pairs of observed galaxies, a significant systematic in weak lensing
surveys (Kirk et al., 2015; Kiessling et al., 2015) which can bias
cosmological inference if not properly modelled (Bridle & King, 2007; Kirk et
al., 2012; Krause et al., 2016; Blazek et al., 2019).

Simulations find that the major axis of the halo inertia tensor tends to be
preferentially aligned along filaments (Altay et al., 2006; Hahn et al., 2007;
Aragón-Calvo et al., 2007; Zhang et al., 2009; Libeskind et al., 2013; Forero-
Romero et al., 2014) and parallel to the surface of voids (Patiri et al.,
2006; Brunino et al., 2007; Cuesta et al., 2008). Parallel alignment is
generally found across all halo masses, though it is typically weaker at low
masses and stronger at high masses (Hahn et al., 2007; Aragón-Calvo et al.,
2007; Zhang et al., 2009; Libeskind et al., 2013; Forero-Romero et al., 2014;
Chen et al., 2015; Codis et al., 2015; Ganeshaiah Veena et al., 2018; Pandya
et al., 2019).

Theory further predicts that halo spin alignments transition from
preferentially parallel to filaments at low mass, to perpendicular at high
mass. At low mass, halos grow by accretion of matter in the plane
perpendicular to the filament, creating a rotating disk in this plane (Pichon
et al., 2011; Codis et al., 2012). High mass halos grow primarily by mergers,
converting motion along the filament into angular momentum. This leads to a
characteristic transition mass at which halos flip from aligned to anti-
aligned, typically $M_{\rm halo}\sim 10^{12}M_{\odot}$ (Bailin & Steinmetz,
2005; Aragón-Calvo et al., 2007; Hahn et al., 2007; Codis et al., 2012;
Trowland et al., 2013; Aragon-Calvo & Yang, 2014; Ganeshaiah Veena et al.,
2018; Wang & Kang, 2018; López et al., 2019; Ganeshaiah Veena et al., 2021).
The alignment signal appears to depend on redshift: Codis et al. (2012) find
that the transition mass for halo spins declines with time to higher masses.

Of course, strong alignment between the cosmic web and dark matter halos does
not necessarily imply the same degree of alignment between the cosmic web and
baryonic galaxies. However, simulations suggest that during the main epoch of
galaxy formation (at high redshift), the shape of baryonic matter tends to
trace the halo shape. This weakens during later times, but is still
significant at low redshift (Tomassetti et al., 2016). In addition, while
interactions between group centrals and satellites are expected to decrease
alignments between the cosmic web and galaxies, observed satellite fractions
at $z\sim 1-2$ are $\lesssim 0.1$, so that the effects of central-satellite
galaxy alignments are relatively minor at higher redshift (McCracken et al.,
2015; Ishikawa et al., 2017). Looking only at cosmic web-galaxy alignment, we
therefore expect qualitatively similar alignments to the cosmic web-halo case,
which decline in strength at lower redshifts. This is borne out by Codis et
al. (2018), who finds that low-mass galaxies are more strongly aligned at high
redshift. Ganeshaiah Veena et al. (2019) and Bate et al. (2020) also find that
high-mass galaxies show stronger perpendicular alignment at high redshift. As
with halos, galaxy spins undergo a flip from alignment to anti-alignment at
$M_{\star}\sim 10^{10.5}M_{\odot}$ (Dubois et al., 2014; Codis et al., 2015,
2018; Wang et al., 2018a; Kraljic et al., 2020). However, in contrast with
halos, Codis et al. (2018) and Kraljic et al. (2020) do not detect significant
redshift evolution in the galaxy spin transition mass.

Observational evidence for this picture is difficult, owing to the weakness of
the signal and the high completeness and galaxy density required to
reconstruct the cosmic web. A number of low-redshift studies have found
evidence supporting the theoretical picture, in which spiral galaxies spin
parallel to filaments and elliptical galaxies spin perpendicular to filaments
(Tempel et al., 2013; Tempel & Libeskind, 2013; Pahwa et al., 2016; Blue Bird
et al., 2020), although the significance is weak and some studies do not
detect this purported alignment (Krolewski et al., 2019). Similarly, Zhang et
al. (2013) find stronger shape alignments for red galaxies than for blue
galaxies. Meanwhile Welker et al. (2020) have presented the first confident
detection of the mass-dependent spin flip using filaments from GAMA and galaxy
spins from SAMI, and Kraljic et al. (2021) report a similar detection using
MaNGA spins.

The observational requirements for spin and shape alignment measurements are
sufficiently challenging that few studies have pushed beyond $z\sim 0$. The
main difficulty is a lack of a sufficiently high density spectroscopic surveys
to resolve the cosmic web over large volumes, although VIPERS has detected the
cosmic web to $z\sim 0.7$ (Guzzo et al., 2014; Malavasi et al., 2016). At
$z\sim 2$, the technique of Lyman-$\alpha$ forest tomography offers an
alternative method to mapping the density field on Mpc scales (Pichon et al.,
2001; Caucci et al., 2008). As Lyman-$\alpha$ forest spectra are fundamentally
one-dimensional tracers of the underlying density, density reconstruction over
large volumes only requires dense areal sampling, as opposed to the dense
volumetric sampling required for reconstruction from traditional galaxy
surveys (Lee et al., 2014). The CLAMATO Ly$\alpha$ tomographic survey (Lee et
al., 2018; Horowitz et al., 2022) has reconstructed the $z\sim 2$ density
field over a volume of $4.1\text{\times}{10}^{5}$
$h^{-3}\,\mathrm{cMpc}^{3}{}$ in the COSMOS field to a resolution of $2\sim 3$
Mpc. However, with only 810 unique coeval galaxies across existing
spectrographic surveys in the reconstructed volume, it is unable to
conclusively detect alignments (Krolewski et al., 2017a; Momose et al., 2022).
The upcoming Subaru Prime Focus Spectrograph Galaxy Evolution Survey offers an
improvement on both fronts: a high-density spectroscopic survey (Greene et
al., 2022) that can resolve the cosmic web at $z\sim 1-2$ with
$N_{g}\sim$250,000 galaxies, and a Ly$\alpha$ forest tomography survey $z\sim
2.5$ that will also cover $N_{g}\sim$ 20,000 coeval galaxies. In this paper,
we explore the prospects for constraining galaxy alignments with upcoming
Subaru PFS data.

In Section 2, we first explore idealized galaxy alignment mass trends from the
IllustrisTNG cosmological simulation, with an eye towards estimating the
galaxy spin transition mass and constraining the strength of alignment
expected from observational surveys. We describe the IllustrisTNG data used in
Section 2.1, and discuss our results in Section 2.2. In Section 3.1, we
construct a mock Subaru PFS observational survey from the IllustrisTNG
simulation, with appropriate galaxy selections and an accompanying Ly$\alpha$
forest tomography survey. Finally, we present and discuss our results for the
alignment constraints produced by our mock survey in Section 3.2.

Throughout this work, we adopt the Planck Collaboration et al. (2016)
concordance $\Lambda$CDM cosmology ($\Omega_{m}$ = 0.3, $\Omega_{b}$ = 0.047,
$h$ = 0.685, $n_{s}$ = 0.965, $\sigma_{8}$ = 0.8).

##  2 Simulation Results

###  2.1 Data and Methods

####  2.1.1 IllustrisTNG simulations

In this work, we use the IllustrisTNG suite of cosmological
magnetohydrodynamic simulations, which is based on the Arepo adaptive mesh
refinement code. For our main analysis we use the TNG300-1 simulation box,
which has $2500^{3}$ dark matter particles, and a side length
$L=205\,h^{-1}\,\mathrm{cMpc}{}$, for a total volume of
$8.61\text{\times}{10}^{6}$ $h^{-3}\,\mathrm{cMpc}^{3}$ (Nelson et al., 2019;
Naiman et al., 2018; Pillepich et al., 2018; Marinacci et al., 2018; Nelson et
al., 2018; Springel et al., 2018). This offers the largest cosmological volume
and galaxy statistics in the publicly available IllustrisTNG suite, and is
better matched to upcoming spectroscopic surveys than its higher-resolution
counterparts. We will also use the TNG100-1 volume to compare the effect of
simulation grid resolution on our shape and spin alignment trends. This
higher-resolution suite has 18203 dark matter particles with a side length of
$L=75$ $h^{-1}\,\mathrm{cMpc}{}$. Of the $100$ cosmological-redshift snapshots
of the simulations that span $z=127$ to $z=0$, we use the data products at
redshifts $z=2$ and $z=1$ for this work. While this does not exactly match the
expected redshift ranges of $z\sim 2.2-2.7$ and $z\sim 0.7-1.7$ for which the
Subaru PFS Galaxy Evolution Survey will obtain galaxy spectra, we believe the
two selected redshifts are representative and our conclusion should still
hold. The baryon particles in TNG300-1/TNG100-1 have mass resolution of
$7.6\text{\times}{10}^{6}$ M⊙/$9.4\text{\times}{10}^{5}$ M⊙ respectively,
while the dark matter particle masses are $5.9\text{\times}{10}^{7}$
M⊙/$7.5\text{\times}{10}^{6}$ M⊙ respectively. These resolutions are
sufficient to resolve both the large-scale morphologies of galaxies and the
surrounding cosmic web in detail.

The IllustrisTNG suite assumes the same Planck Collaboration et al. (2016)
cosmology we use in this work. To construct the matter density field, we
construct a $512^{3}$ grid mesh and assign each dark matter particle to its
nearest cell.

####  2.1.2 Galaxy properties

In IllustrisTNG, subhalos are identified via the subfind algorithm developed
by Springel et al. (2001). For our analysis of morphological alignment trends
in TNG300-1 and TNG100-1, we use the morphology calculated by Shi et al.
(2021). This morphology sample is calculated for subhalos which have a
cosmological origin (SubhaloFlag $\neq$ 0), and have $M_{R/2}{}>10^{9}$
$h^{-1}\,M_{\odot}$, where $M_{R/2}$ is the stellar mass within twice the
half-stellar-mass radius. However, we note that for the rest of this work, we
use “stellar mass” to refer to total stellar mass to be consistent with
previous works, such as Codis et al. (2018). Our sample includes both galaxies
classified as central galaxies (most massive within each dark matter halo) and
satellite galaxies (all others associated with a given dark matter halo).

To measure galaxy morphology, Shi et al. (2021) adopts the reduced mass
inertia tensor formalism described in Tenneti et al. (2015). In this
formalism, galaxies are modelled as 3D triaxial ellipsoids, with axis lengths
as $a\geq b\geq c$, and corresponding axis vectors as
$\vec{a},\vec{b},\vec{c}$. These quantities are found by calculating the
sorted eigenvalues and eigenvectors respectively of the reduced mass inertia
tensor $\tilde{I}$:

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

where $0\leq i,j\leq 2$ represent coordinate indices, $m_{n}$ is the mass of
the $n^{th}$ particle in the galaxy, and $x_{ni}$ represents the $i^{th}$
coordinate of the $n^{th}$ particle (with respect to the galaxy center of
mass). $r_{n}^{2}$ is defined as:

$r_{n}^{2}=\sum_{i}x^{2}_{ni}$

We henceforth refer to the longest axis $\vec{a}$ as the ‘shape’ direction of
each galaxy in this sample.

For our analysis of spin alignment trends in TNG300-1 and TNG100-1, we will
use the stellar particle-only spin also calculated by Shi et al. (2021) which
we refer to as the ‘spin’ sample. We limit our sample to subhalos with more
than 50 stellar particles, which ensures that they have a well-defined spin
vector following the convention set by Codis et al. (2018).

Table 1 summarizes the number of galaxies in each simulation and redshift
sample, and gives 5% and 95% total stellar mass percentiles for each sample.

Simulation | Redshift | Sample | # galaxies | 5%ile $M_{*}$ ($h^{-1}\,M_{\odot}$) | 95%ile $M_{*}$ ($h^{-1}\,M_{\odot}$)
---|---|---|---|---|---
TNG300-1 | z=2 | Shape | 128,019 | $1.348\,882\text{\times}{10}^{09}$ | $3.827\,169\text{\times}{10}^{10}$
|  | Spin | 256,309 | $3.735\,926\text{\times}{10}^{08}$ | $2.392\,023\text{\times}{10}^{10}$
| z=1 | Shape | 194,189 | $1.384\,419\text{\times}{10}^{09}$ | $4.568\,499\text{\times}{10}^{10}$
|  | Spin | 345,296 | $3.731\,775\text{\times}{10}^{08}$ | $3.260\,130\text{\times}{10}^{10}$
TNG100-1 | z=2 | Shape | 8,530 | $1.409\,582\text{\times}{10}^{09}$ | $3.984\,801\text{\times}{10}^{10}$
|  | Spin | 61,976 | $4.677\,695\text{\times}{10}^{07}$ | $6.630\,453\text{\times}{10}^{09}$
| z=1 | Shape | 13,245 | $1.394\,461\text{\times}{10}^{09}$ | $5.170\,797\text{\times}{10}^{10}$
|  | Spin | 65,641 | $4.503\,330\text{\times}{10}^{07}$ | $1.266\,712\text{\times}{10}^{10}$
Table 1: Number of galaxies in each sample, with 5% and 95% percentiles of the
total stellar mass. Shape samples have reduced inertia tensor morphologies
calculated by Shi et al. (2021), and are selected via a stellar mass
threshold, while spin samples only have angular momentum vectors, and are
selected via a stellar particle number threshold.

####  2.1.3 Deformation tensor and alignment angle

In order to quantify the large-scale dynamics of the matter density field
without directly using peculiar velocities, we adopt the deformation tensor
approach of Krolewski et al. (2017a) and Lee & White (2016), which in turn is
based on the basic formalism described in Hahn et al. (2007) and Forero-Romero
et al. (2009). The deformation tensor is defined as the Hessian of the
gravitational potential $\phi$ of each point in space:

$D_{ij}=\partial x_{i}\partial x_{j}\phi$ (2)

In practice, the deformation tensor is efficiently calculated in Fourier space
directly from the density, as described in Krolewski et al. (2017a).

Diagonalizing the tensor yields a set of 3 eigenvalues and associated
normalized eigenvectors for each point. By convention, we sort the eigenvalues
such that $\lambda_{1}>\lambda_{2}>\lambda_{3}$. Whether the eigenvalues are
indexed in increasing or decreasing order varies from work to work, but
produces the same results with the indices switched accordingly. In the
Zel’dovich approximation, the deformation tensor solely determines the time
evolution of the initial density field (Zel’Dovich, 1970): matter first
collapses along the eigenvector $\boldsymbol{\hat{e}_{1}}{}$ associated with
$\lambda_{1}$, and collapses last along $\boldsymbol{\hat{e}_{3}}{}$. Within
the approximation, a positive/negative eigenvalue indicates that matter is
collapsing/expanding along the axis of the associated eigenvector.

The dark matter density is first smoothed to eliminate the effect of highly
nonlinear fluctuations in density before calculating the deformation tensor.
We apply a Gaussian kernel with standard deviation 2
$h^{-1}\,\mathrm{cMpc}{}$, which is comparable to smoothing scales used by
other shape and spin alignment studies (Codis et al., 2015; Forero-Romero et
al., 2014).

In order to characterize the degree to which a galaxy’s shape or spin vector
is aligned with the underlying matter deformation tensor at the galaxy’s
position, we calculate $\lvert\cos{\theta}\rvert$, where $\theta$ is the angle
between the shape/spin and each of the deformation tensor’s eigenvectors at
the nearest density grid point to the galaxy. We take the absolute value since
the sign convention of the eigenvectors is arbitrary.
$\lvert\cos{\theta}\rvert{}>0.5$ means the two vectors are largely parallel
with each other, while $\lvert\cos{\theta}\rvert{}<0.5$ means they tend to be
perpendicular with each other.

###  2.2 Idealized Simulation Alignments

In this section, we discuss the alignment trends of the simulated galaxy
shapes and spins, with respect to the cosmic web that traces the underlying
matter density field.

####  2.2.1 Alignment of galaxy shapes

(a) $z=1$

(b) $z=2$

Figure 1: Shape-cosmic web alignment signal for all galaxies with radius-
limited stellar mass $M_{*}>10^{10}\,h^{-1}\,M_{\odot}$, as a function of
total stellar mass. For each redshift, the top panel shows the alignment
signal calculated in total stellar mass bins. In the bottom panel, the
alignment mean is calculated as a function of a sliding stellar mass bin.
Error bars are calculated using bootstrap sampling, with error bars/shaded
regions indicating the 16%-84% percentile ranges. The subscript on
$\theta_{s}$ indicates that we are calculating shape alignment.

In Figure 1, we show the 3-dimensional alignment signal between shape vectors
and the cosmic web for the full TNG300-1 galaxy sample from Section 2.1.2, for
$z=1$ and $z=2$ respectively. As shape-$\boldsymbol{\hat{e}_{2}}{}$ and
shape-$\boldsymbol{\hat{e}_{1}}{}$ alignments are largely the inverse of
shape-$\boldsymbol{\hat{e}_{3}}{}$ alignments (i.e. perpendicular when the
latter is parallel), we opt to plot only shape-$\boldsymbol{\hat{e}_{3}}{}$
alignments for the sake of brevity. In the top panel of each figure, we
compare the distribution of alignment angle $\lvert\cos{\theta}\rvert$ in four
stellar mass bins. We refer to a positively-sloped bin trend (i.e. toward an
excess at positive $\lvert\cos{\theta}\rvert$ to the two vectors) as a
positive alignment, and vice versa as a negative alignment (i.e.
perpendicular).

On the bottom panel of each figure, we adopt a “sliding-window” mass binning
approach for characterizing the overall alignment trend as a function of
stellar mass. From all galaxies within a mass range of width $0.5$ in
log10$M_{*}$ space, we compute the mean alignment
$\langle\lvert\cos{\theta}\rvert{}\rangle$. A mean alignment of
$\langle\lvert\cos{\theta}\rvert{}\rangle=0.5$ corresponds to random
orientations and no systematic alignment. In this panel we also show results
from TNG100-1, which has a higher particle/cell mass resolution than TNG300-1.
The trade off, however, is its reduced volume compared to TNG300-1 (753
$h^{-3}\,\mathrm{cMpc}^{3}$ vs. 2053 $h^{-3}\,\mathrm{cMpc}^{3}$, i.e. $\sim
20\times$ smaller) which leads to greater cosmic variance in the simulation
box. This is reflected in the larger bootstrap error bars. The number of
galaxies in the highest-mass bin (smallest galaxy number across all bins) is
186 for TNG100-1, so we believe bootstrap errors should still be fairly
accurate for the high-mass end.

The $\boldsymbol{\hat{e}_{3}}$-shape alignment shows a clear transition mass
at both redshifts: at lower stellar masses there is a flat trend consistent
with no alignment, which changes to an increasing alignment trend with stellar
mass beyond the transition. This is consistent with past studies of dark
matter halo shape alignments with $\boldsymbol{\hat{e}_{3}}$ (Libeskind et
al., 2013; Zhang et al., 2009), which are similar to galaxy shape alignments
due to central galaxies being mostly aligned with their host halo at high
redshifts (Tomassetti et al., 2016). The transition mass (the bin center at
which the bin’s mean alignment is 0.5) is $M_{*}\sim 10^{10.5}$
$h^{-1}\,M_{\odot}$ for $z=2$, and $M_{*}\sim 10^{10.1}$ $h^{-1}\,M_{\odot}$
for $z=1$.

The mass-alignment trends for TNG100-1 are broadly consistent with TNG300-1
(within 1$\sigma$ errorbars for TNG100-1) for $\boldsymbol{\hat{e}_{3}}$ at
both redshifts.

####  2.2.2 Alignment of galaxy spins

(a) $z=1$

(b) $z=2$

Figure 2: Spin-cosmic web alignment signal for all galaxies with more than
$50$ stellar particles, as a function of total stellar mass. For each
redshift, the top panel shows the detailed alignment signal calculated in
total stellar mass bins. The alignment mean is calculated as a function of a
sliding stellar mass bin in the below panel. The error bars are calculated
using bootstrap sampling, with error bars/shaded regions indicating the
16%-84% percentile ranges. The subscript on $\theta_{j}$ indicates that we are
calculating spin alignment.

In Figure 2, we show the alignments between galaxy spins (i.e. angular momenta
of their stellar content) with respect to the $\boldsymbol{\hat{e}_{3}}{}$
cosmic web eigenvector for varying stellar mass at $z=1$ and $z=2$
respectively. As our galaxy sample for spin is larger than that of shape at
the low-mass end due to relaxed selection criteria, we show sliding-window
bins down to a bin centered at $M_{*}=10^{8.5}$ $h^{-1}\,M_{\odot}{}$, in
contrast to the lowest-mass bin of $M_{*}=10^{9}$ $h^{-1}\,M_{\odot}{}$ for
shape alignment.

At $z=2$, we appear to see no trends for spin-$\boldsymbol{\hat{e}_{3}}$
alignment at the >$1\sigma$ level for both TNG100-1 and TNG300-1. This is
somewhat surprising, as previous DM halo alignment studies (Zhang et al.,
2009; Libeskind et al., 2013) and hydrodynamic galaxy alignment studies (Wang
et al., 2018b; Codis et al., 2018, 2015; Kraljic et al., 2020) have observed a
clear “spin-flip” transition at stellar masses of $M_{\text{flip}}\sim
10^{9.5}-10^{10.5}$ $h^{-1}\,M_{\odot}$ akin to the shape alignment transition
we observe at $z=2$.

At $z=1$, the spin-$\boldsymbol{\hat{e}_{3}}$ alignment trend for TNG300-1 has
a significant (>$1\sigma$) transition towards negative alignment at
$M_{*}=10^{10.2}$ $h^{-1}\,M_{\odot}$. This mass-dependent behavior has not
been seen in past simulated galaxy spin alignment studies such as Wang et al.
(2018b), who studied the mass dependence of galaxy spin alignment in the
predecessor Illustris simulation and report a monotonically decreasing
spin-$\boldsymbol{\hat{e}_{3}}$ alignment for both $z=1$ and $z=2$. While this
difference may come from differing simulation physics and differences in
galaxy selection criteria, it may also be partly a result of the broad (order-
of-magnitude in $M_{*}$) and non-overlapping mass bins used by Wang et al.
(2018b), which could obscure finer mass trends for alignment.

For both redshifts, TNG100-1 is again broadly consistent with TNG300-1 (within
1$\sigma$ errorbars for TNG100-1).

####  2.2.3 Comparison with Horizon-AGN results

To investigate the effects of the implemented sub-grid physics models in
cosmological hydro-dynamical simulations on galaxy-cosmic web alignments, in
this subsection, we compare the spin alignments of IllustrisTNG with that from
the Horizon-AGN simulation. The Horizon-AGN simulation (Dubois et al., 2014)
is a hydrodynamic adaptive mesh simulation over a box size of
$L=100\,h^{-1}\,\mathrm{cMpc}$, with a dark matter particle mass resolution of
$8\text{\times}{10}^{7}$ M⊙, which is comparable to the mass resolution of
TNG300-1 ($5.9\text{\times}{10}^{7}$ M⊙). As spin-cosmic web alignments in
Horizon-AGN have been extensively studied (Dubois et al., 2014; Codis et al.,
2015, e.g.), we opt to use it for comparison. In particular, we compare the
spin alignments from the IllustrisTNG $z=1$ with the spin alignments
calculated by Codis et al. (2015) for $z\sim 1.2$.

We impose a minimum stellar mass cut on the galaxy sample for TNG300-1 and
TNG100-1 used in the previous sections of $M_{*}\geq 10^{8.06}$
$h^{-1}\,M_{\odot}$, to match the lower stellar mass limit used in Codis et
al. (2015)’s analysis. The TNG300-1 sample is unchanged due to the minimum
stellar particle count imposed, but the TNG100-1 sample is reduced by
$\sim$37%. We then calculate the spin-cosmic web alignment using the density
field smoothed with a kernel size of $200$ $h^{-1}$ ckpc used in Codis et al.
(2015), instead of the 2 $h^{-1}\,\mathrm{cMpc}$ kernel used in the rest of
our work.

Figure 3: Comparison of Horizon-AGN ideal spin alignments from Figure 3 of
Codis et al. (2015) (in black) with corresponding spin alignments from
TNG300-1 (blue) and TNG100-1 (red). TNG samples are derived from the samples
described in Section 2.1.2, with an additional minimum-mass cut of
$M_{*}>10^{8.06}$ $h^{-1}\,M_{\odot}$. Errorbars for the TNG300-1/TNG100-1
alignments are derived from bootstrap sampling, while errors on Codis et al.
(2015) are Poisson, taken from the original figure.

In Figure 3, we compare TNG spin alignments to those found by Codis et al.
(2015) for the total (non-binned) Horizon-AGN galaxy sample, with
corresponding errorbars (bootstrap for TNG alignments, Poisson for Horizon-
AGN). Horizon-AGN displays a clear, high-amplitude negative trend for
$\boldsymbol{\hat{e}_{3}}$ and positive trends for $\boldsymbol{\hat{e}_{2}}$
and $\boldsymbol{\hat{e}_{1}}$ as expected from theory. Spin alignments with
comparable amplitudes at $z=1$ have also been observed for the SIMBA
hydrodynamic simulation111However, note that the DISPERSE cosmic web finder
(Sousbie, 2011) is used instead of the deformation tensor formalism for that
work. (Kraljic et al., 2020). However, we find that alignment directions for
TNG100-1 in the same mass range are more ambiguous, which can also be seen in
Figure 2a. While TNG300-1 displays alignments in the same sign as theory and
Horizon-AGN, their amplitudes are much lower than Horizon-AGN’s.

To quantify the level of disagreement between Horizon-AGN and TNG, we compute
the ratio of $\langle\lvert\cos{\theta}\rvert{}\rangle$ above null (0.5) of
Horizon-AGN and TNG, $r_{\text{MA}}$:

$r_{\text{MA}}=\frac{\langle\lvert\cos{\theta}\rvert{}\rangle{}_{\text{Horizon-
AGN}}-0.5}{\langle\lvert\cos{\theta}\rvert{}\rangle{}_{\text{TNG}}-0.5}$ (3)

As we expect alignments for $\boldsymbol{\hat{e}_{3}}$ and
$\boldsymbol{\hat{e}_{1}}$ to be negatively/positively trending from theory
for both Horizon-AGN and TNG, we only calculate $r_{\text{MA}}$ for alignments
with those two eigenvectors.

Between Horizon-AGN and TNG300-1, $r_{\text{MA}}$ for
$\boldsymbol{\hat{e}_{3}}$ and $\boldsymbol{\hat{e}_{1}}$ respectively is
$(10.8,6.9)$. Between Horizon-AGN and TNG100-1, $r_{\text{MA}}=(11.5,2.4)$.
Compared to Horizon-AGN’s galaxy sample, TNG300-1 is missing some galaxies on
the low mass end due to mass resolution limits. As we see in Figure 2 (and
expect from theory), lower-mass galaxies have less alignment strength, so we
caution that our TNG300-1 $r_{\text{MA}}$ is likely overestimated. To be
conservative, we therefore take the smallest ratio between
$\boldsymbol{\hat{e}_{3}}$ and $\boldsymbol{\hat{e}_{1}}$ for the
TNG100-1/Horizon-AGN comparison to be the “overall ratio of alignment
strength” between Horizon-AGN and TNG, $r^{*}_{\text{MA}}=2.4$.

While a full analysis of the reasons behind this large alignment discrepancy
is beyond the scope of this work, we speculate that this discrepancy may be
partially caused by differences in feedback prescriptions between the
simulations. Tomassetti et al. (2016) find that the alignment between galaxies
and their surrounding halos has a large dependence on feedback strength. In
particular, stronger feedback increases the strength of galaxy-halo alignment
(and thus, galaxy-cosmic web alignment) at high redshifts, by delaying the
time at which the baryonic component of galaxies compactifies. Observational
measurements of high-redshift alignments could therefore be a very promising
constraint on feedback models for simulations; and, in general, on their
subgrid physics models.

##  3 Mock Subaru PFS Surveys

###  3.1 Methods

In the previous section, we examined the level of alignment between the
galaxies and underlying cosmic web within the IllustrisTNG simulations, which
reflect the underlying ‘truth’ in the simulations. We next turn to the
question of whether these alignments, e.g. as shown in Figures 1, can actually
be detectable in upcoming spectroscopic surveys. This is especially relevant
in the context of the Galaxy Evolution component of the Subaru PFS Strategic
Program survey (Greene et al., 2022), which will for the first time cover
large enough cosmological volumes with sufficient fidelity to reconstruct the
cosmic web at high redshifts. We will construct mock galaxy observations and
cosmic web reconstructions, adopting parameters planned for the PFS Galaxy
Evolution Survey components at $z=1$ and $z=2$.

We first select a sample of galaxies that match the expected numbers of
spectroscopic redshifts by matching their photometric magnitudes. We then
displace the galaxies into redshift space, i.e. convolving their radial
positions with their line-of-sight velocities (Section 3.1.2). For $z=1$, we
infer the reconstructed dark matter density directly from the observed galaxy
distribution (Section 3.1.3), while for $z=2$, we independently reconstruct
the dark matter density using Lyman-alpha tomography (Section 3.1.4). The
redshift-space distortion and both of these reconstruction methods depend on
the viewing angle onto the simulation, which we marginalize over for reasons
detailed in Section 3.1.1. For both redshifts, we calculate the deformation
tensor (the reconstructed cosmic web) from the reconstructed dark matter
density.

Galaxy shapes for PFS galaxies can be estimated from deep high-resolution
images at wavelengths covering their stellar components. Galaxy spins, the
other hand, require resolved spectroscopy (typically with integral field unit
spectrographs, or IFUs) to measure their rotations and hence angular momenta.
While it is in principle possible to obtain IFU data of galaxies at our $z\sim
1-2$ redshifts of interest, the number of galaxies that would be need to be
observed to cover a substantial fraction of the PFS sample needs to
$>100\times$ larger than existing high-redshift IFU samples such as Förster
Schreiber et al. 2019. We therefore do not expect the galaxy spin-cosmic web
alignment at $z=1$ and $z=2$ to be measurable in the foreseeable future, and
will focus solely on shape-cosmic web alignments in the remainder of this
work.

####  3.1.1 Viewing Angle Variance

When observing galaxies in reality, we are only able to observe their shape
projected onto the 2D plane of the sky. While for disc-like galaxies it might
be possible to estimate the inclination from the galaxy’s detailed structure
and therefore reconstruct 3D shape, this has only been attempted with
reasonable effectiveness in the Local Universe (e.g., Haynes & Giovanelli,
1984) — we are certainly unable to resolve comparably detailed galaxy
structure for large galaxy samples at $z\gtrsim 1$. In addition, this method
is not applicable to elliptical-type galaxies, which typically have stronger
alignment with the cosmic web (Pandya et al., 2019). For our mock
observational survey analysis, we therefore project the shape of each galaxy
onto the plane normal to a given viewing angle, which we approximate as the
same across the simulation box222The angular scale subtended by a transverse
comoving distance of $\sim 200\,h^{-1}\,\mathrm{cMpc}$ (the extent of
TNG300-1) at $z\sim 2$ is $\sim 3\deg$, for which small-angle approximations
are still reasonable..

For galaxy shapes, we start with the full galaxy reduced inertia tensor, which
we interpret as an ellipsoid quadratic-form matrix. We project the inertia
tensor ellipsoid onto the plane normal to the viewing angle, and take the new
longest axis of the projected elliptical shadow as our projected shape vector.
We also project all deformation tensor eigenvectors onto the viewing angle
normal plane, and then calculate $\lvert\cos{\theta}\rvert$. As vectors with
independent random directions in 3D reduce to vectors with independent random
directions in 2D when projected onto a plane, we now use the alignment
distribution of independently random 2D vectors as our null distribution. The
distribution of $\cos\theta$ for such vector pairs is the beta distribution
rescaled onto the interval $\cos\theta\in[-1,1]$ (centered at 0) with
distribution parameters $\alpha=\beta=0.5$ [97]. The distribution of
$\lvert\cos{\theta}\rvert$ is the renormalized upper half of that beta
distribution.

However, in preliminary analysis we find that this projected alignment is
significantly dependent on the selected viewing angle to the simulation
volume. In other words, a given sample of galaxies in the simulation could
vary significantly by choosing the viewing angle from which the 2D shape-
cosmic web alignments are computed. When averaged over large numbers of
viewing angles, the alignment signal varies significantly: in some cases, the
mean alignment magnitude is reduced up to a factor of $\sim 2$ compared to the
naive result computed from some of the individually strong viewing angles. We
hereafter refer to this effect as “viewing-angle variance”.

Figure 4: Cartoon illustrating viewing-angle variance on alignments caused by
cosmic variance. For the small simulation volume depicted here with a filament
running through it, galaxies are highly aligned with the filament in 3D. This
extends to the projected alignment when looking “side-on” to the filament
(left panel). However, if we view the volume “face-on” to the filament (right
panel), the projected alignment is nearly random.

We hypothesize that this variation is a result of global anisotropies in the
eigenvector orientations (i.e. dark matter flows) caused by cosmic variance
within the limited simulation volumes. As a toy example, consider a small
($L\sim$ 3-10 $h^{-1}\,\mathrm{cMpc}$ side length) volume depicted in Figure
4, which is taken up by a single large filament passing through, with galaxies
aligned with the filament.

For viewing angles looking “down the barrel” at the filament, projected
eigenvectors would be sensitive to small perturbations in viewing angle, and
vice versa for viewing angles looking length-wise at the filament333A similar
effect is responsible for the systematics discussed in Lamman et al. (2022),
although their application is very different from ours.. This would therefore
induce a high degree of viewing-angle variance in projected alignment and
derived quantities such as mean projected alignment. Given this hypothesis, we
expect the viewing-angle variance in mean projected alignment to scale
inversely with volume size.

We investigate the possibility of this inverse scaling relationship between
volume size and viewing-angle variance using the “true” IllustrisTNG data
discussed in Section 2.1, for $z=1$. To do so, we project the galaxy shape
vectors and deformation tensor eigenvectors onto 64 evenly-spaced viewing
angles into the volume. Next, we divide the $L=205\,h^{-1}\,\mathrm{cMpc}{}$
side-length volume into $N=${1, $2^{3}$, $4^{3}$, $8^{3}$} evenly-sized, non-
overlapping sub-cubes. For each sub-cube, we consider the $64\times N$
projected alignments computed from the $N$ galaxy-eigenvector pairs within the
sub-cube. For each viewing angle individually, we compute the mean projected
alignment across all galaxy-eigenvector pairs. Finally, we compute the
standard deviation of those 64 mean projected alignments, which we take as a
measure of the viewing-angle variance for that sub-cube, $\sigma_{VA}$.

Figure 5: Standard deviation of projected $z=1$ real-space mean alignment
$\sigma_{VA}$ computed from 64 viewing angles of the same sub-volumes, plotted
as a function of subvolume side length. As we divide the TNG300 volume into
non-overlapping subvolumes, we can also plot the error on $\sigma_{VA}$. The
data are consistent with a power-law with fixed exponent $\beta=-1.5$, which
is equivalent to a $\beta=-0.5$ power-law against volume.

Figure 5 displays the results of this analysis. Since for box sizes besides
$L=205^{3}$ $h^{-3}\,\mathrm{cMpc}^{3}{}$ we have multiple sub-cubes, we
display the mean and standard deviation across sub-cubes of $\sigma_{VA}$. To
avoid using a logarithmic scale on the x-axis, we plot $\sigma_{VA}$ against
sub-cube side length instead of sub-cube volume. As expected, $\sigma_{VA}$
increases with decreasing side length/decreasing volume. To determine the
exact relationship, we try fitting a power-law with fixed exponent
$\beta=-1.5$ to the mean $\sigma_{VA}$ (rather than the individual
$\sigma_{VA}$ data points to avoid overweighting small volumes). We find the
data are consistent with this $\beta=-1.5$ power-law against side length,
implying a $\beta=-0.5$ power-law against volume.

The strength of this viewing angle variance even in the TNG300 volume suggests
that even $L=205\,h^{-1}\,\mathrm{cMpc}$ simulation box sizes are affected by
cosmic variance of large-scale filaments as far as galaxy-cosmic web
alignments are concerned. In our mock survey analysis, we therefore
marginalize over viewing angles by adopting a viewing angle-agnostic approach.
Specifically, we generate 64 evenly spaced points on the unit half-sphere
using a Fibonacci lattice (Hannay & Nye, 2004). We take the unit vectors
connecting those points and the origin as our set of observer viewing angles
into the TNG300-1 simulation. For every viewing angle, we construct mock
galaxy observations with the procedure detailed in the rest of the methods
section.

####  3.1.2 Galaxy sample selection

For our mock observational survey analysis only, we subselect the galaxies in
the simulation to match the expected sample from the Subaru Prime Focus
Spectrograph Galaxy Evolution (hereafter PFS GE) survey, which is planned to
take place over $\sim 100$ nights of observation on the Subaru Telescope. At
$z=0.7\sim 1.7$, the possible targets for the PFS GE Survey number $N_{g}\sim$
357,000 galaxies over a total survey volume of $V_{S}\sim$
$32\,467\,532.467\,532\,467$ $h^{-3}\,\mathrm{cMpc}^{3}$, and at $z=2.1\sim
2.5$, $N_{g}\sim$ 100,000 galaxies over a volume of $V_{S}\sim$
$2.7\text{\times}{10}^{7}$ $h^{-3}\,\mathrm{cMpc}^{3}$ (Greene et al., 2022).

As only some rest-frame galaxy magnitudes and not full synthetic galaxy
spectra are provided in the IllustrisTNG release, we proceed to define the
equivalent targets in the simulated IllustrisTNG catalog by performing an
abundance-matching based on their number densities and ranked-ordered
luminosities within a given filter wavelength range. For the $z=1$ redshift
range, the PFS survey target sample are selected based on Hyper Suprime-Cam
$Y$-band limiting magnitude (centered at approximately $\lambda_{C}\approx$
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ wavelength), so we use the rest-
frame Johnson $V$-band ($\lambda_{C}\approx 0.551$
$\text{\,}\mathrm{\SIUnitSymbolMicro m}$) photometric magnitude supplied in
the IllustrisTNG catalog for each galaxy as an approximate proxy for the
observed-frame $Y$-band magnitude. This approach ignores filter bandpass
shapes, and so should be considered approximate; additionally, the intrinsic
magnitudes given by IllustrisTNG do not include the effects of dust
extinction. Using the V-band magnitude, we select the brightest 94,766
galaxies in our original sample, to match the expected number density of
possible PFS $z\sim 1$ survey targets:
$n_{g}=0.011h^{3}\,\mathrm{cMpc}^{-3}{}$.

We perform the same procedure for selecting survey galaxies at the $z=2$
range, which are part of the ‘IGM foreground’ galaxy sample designed to
provide a sample of galaxies within the volume probed by IGM tomography from
background galaxy and quasar Lyman-$\alpha$ forest sightlines. Here, observed-
frame $J$-band at $\lambda_{C}\approx 1.22$
$\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is the limiting magnitude in target
selection, which we approximate by selecting the 31,908 brightest galaxies in
rest-frame $B$-band magnitude ($\lambda_{C}\approx 445$ nm) We match the
selection of the targets to the number density of candidates, which is
$n_{g}=3.7\times 10^{-3}\,h^{3}\,\mathrm{cMpc}^{-3}{}$.

We emphasize that the number of possible survey targets does not equal the
number of galaxies which will have successful redshift measurements. PFS
expects $\sim$70% of observed targets to yield successful redshifts, for a
final $z\sim 1$ data sample of $N_{g}\approx$ 250,000 galaxies. For $z\sim 2$,
22% of possible targets will actually be observed (i.e. the sampling rate),
and of those targeted, 70% are expected to yield successful redshifts. This
translates to a final data sample of $N_{g}\approx$ 15,400 galaxies at $z\sim
2$. We also consider a scenario where the $z=2$ component is allocated double
the baseline observing time within the survey, which increases the sampling
rate to 44%, and the final data sample to 30,800 galaxies. Although the
spectroscopic redshift success rate will depend on a variety of factors in the
real observations, here we will make the simplifying assumption that every
galaxy above the magnitude limit has an equal chance to yield a successful
redshift for our analysis. Due to this assumption, it is valid for us to
sample without replacement the absolute number of successful galaxies PFS
expects from our simulation’s possible target pool, even if the number density
of galaxies we select is higher than the actual number density of successful
galaxies that PFS will observe. We therefore create two samples with 15,400
and 30,800 galaxies for $z=2$ out of the 31,908 possible targets (ensuring
that the former sample is a strict subset of the latter). For the remainder of
this work, we refer to these two samples as Z2_FID (“$z=2$ fiducial”) and
Z2_2X (“$z=2$ with 2$\times$ greater number density”) respectively. The
$N_{g}=$ 250,000 observed galaxies expected at $z=1$ in fact exceeds the total
number of galaxies within TNG300-1 that have the correct luminosities.
Therefore we use all 94,766 possible targets within the TNG300-1 volume and
perform a post-hoc rescaling of the uncertainties. This will be described in
Section 3.1.5. For the remainder of the work, we refer to this $z\sim 1$
dataset as Z1_FID.

####  3.1.3 $z=1$ cosmic web reconstruction

At $z=1$, the expected number density of spectroscopically-confirmed galaxies
from PFS is high enough to enable accurate reconstruction of the underlying 3D
matter density field using the observed galaxy positions as a biased tracer of
density, as in redshift surveys of the Local Universe ($z\lesssim 0.1$,
Kitaura et al. (2009); Wang et al. (2014); Jasche et al. (2015)). For density
reconstruction, we use the TARDIS density reconstruction code (Horowitz et
al., 2019, 2021), which uses maximum-likelihood analysis in conjunction with a
gravitationally evolved forward model of the matter density field. TARDIS
constrains its search space to density configurations that could have evolved
from initial conditions, and has been shown to accurately reconstruct the
cosmic web eigenvectors of the deformation tensor. We use TARDIS’s “galaxy-
only” density reconstruction mode, in which the matter density is computed
from the distribution of observed galaxies in the forward model using fixed
bias parameters in Lagrangian space. In order to fit these bias parameters, we
match PFS’s expected number density by uniformly subsampling without
replacement 70% of our $z=1$ magnitude-limited data sample, which is 66,336
galaxies. With this galaxy subsample, we estimate the bias parameters to be
$b_{1}=1.8$, $b_{2}=0.1$.

For each viewing angle, we take the same galaxy subsample used to fit the
biases, compute the RSD along the viewing angle’s line of sight, and input
into TARDIS. We opt not to select different subsamples for each viewing angle
because the purpose of varying the viewing angle is to quantify the effect of
viewing angle alone on the alignment. For all viewing angles, we reconstruct
the density on a 384-cell length cube, with box length of $L=$205
$h^{-1}\,\mathrm{cMpc}$. We then smooth the reconstructed real-space density
of TARDIS with a 2 $h^{-1}\,\mathrm{cMpc}$ Gaussian kernel, and compute the
deformation tensor for each smoothed, reconstructed volume.

####  3.1.4 $z=2$ cosmic web reconstruction

In order to reconstruct the cosmic web at $z=2$ where only a relatively low
number density of galaxies will be observed, we carry out a joint
reconstruction with Lyman-$\alpha$ forest tomography (also known as IGM
tomography) as the primary tracer of the cosmic web, with galaxies included as
a secondary tracer. In Ly$\alpha$ forest tomography, redshifted hydrogen
Ly$\alpha$ “forest” absorption lines in the spectra of background sources are
used to infer the neutral density of diffuse IGM hydrogen along the line of
sight. As first argued by Lee et al. (2014), with a sufficient transverse
density of background quasars and UV-bright star-forming galaxies, the 3D IGM
density field can be reconstructed to a 2-3 $h^{-1}\,\mathrm{cMpc}$ spatial
resolution. The COSMOS Lyman-alpha Mapping and Tomography Observations
(CLAMATO, Lee et al. 2018) survey has to-date performed Ly$\alpha$ tomographic
reconstruction over a 0.21 deg2 footprint within the COSMOS field,
reconstructing the IGM density at a spatial resolution of 2.5
$h^{-1}\,\mathrm{cMpc}$ over a survey volume of
$V_{S}=$$4.1\text{\times}{10}^{5}$ $h^{-3}\,\mathrm{cMpc}^{3}$ (Horowitz et
al., 2022). The sources used had an average on-sky transverse separation of
2.4 $h^{-1}\,\mathrm{cMpc}$ at $z=2$. PFS aims to perform Ly$\alpha$
tomographic reconstruction using a similar background source separation, but
for a much larger total volume of $V_{S}\sim$ $2.7\text{\times}{10}^{7}$
$h^{-3}\,\mathrm{cMpc}^{3}$ (an ongoing intermediate-sized survey is the LATIS
survey, Newman et al. (2020)). While the median redshift of the PFS IGM
tomography map is expected to be around $\langle z\rangle\approx 2.5$, this
represents only around a 700 Myr difference in cosmic time compared to $z=2$,
so we expect the overall trends we find in the IllustrisTNG $z=2$ snapshot to
qualitatively reflect the PFS survey if not quantitatively exact.

For each viewing angle, we construct a “joint” galaxy and Ly$\alpha$
tomographic survey using mock data derived from the TNG300-1 simulation. For
the galaxy portion, we again match expected galaxy number densities through
random subsampling of our abundance-matched, magnitude-limited galaxies.

The Ly$\alpha$ tomographic portion of our mock survey is composed of a number
of “skewers” that record relative Ly$\alpha$ absorption in a background
object’s spectrum, along the viewing angle’s line of sight. Since IllustrisTNG
is a hydrodynamic simulation, we directly construct mock Ly$\alpha$ absorption
skewers using the fake_spectra code (Bird et al., 2015).

We follow expected PFS survey parameters by generating skewers at random sky
positions with a number density such that the average transverse separation
between skewers is 2.4 $h^{-1}\,\mathrm{cMpc}$. These skewers do not initially
contain instrumental noise. Skewers which contain a damped Ly$\alpha$ absorber
(DLA), i.e. with integrated HI column density of $N_{\mathrm{HI}}\geq$
$2\text{\times}{10}^{20}\text{\,}\mathrm{c}\mathrm{m}^{-2}$ are filtered out
since we expect to be able to remove them in observational data. We then rebin
the skewers to an instrumental spectral resolution of 0.7
$\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ per pixel, and simulate spectral
dispersion by smoothing the skewers with a 1D Gaussian kernel of size 2.7
$\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. This matches the expected
instrumental performance of the PFS spectrograph’s blue arm, which at a range
of 380-650 nm covers the observed Ly$\alpha$ spectra (Sugai et al., 2015).
Next, we add instrumental noise to the skewers using the noise model described
in Horowitz et al. (2019, 2021) and add an offset to the overall transmitted
flux in each skewer to model systematic fitting errors of the skewer’s
background source continuum Krolewski et al. (2017a).

To reconstruct the underlying dark matter density from our mock skewers, we
apply TARDIS jointly to the galaxy and Ly$\alpha$ tomography data (Horowitz et
al., 2021) with the likelihood now calculated from a weighted combination of
the galaxy density tracer and the Ly$\alpha$ absorption tracer. For all
viewing angles, we again reconstruct the density on a 384-cell length cube and
smooth with a 2 $h^{-1}\,\mathrm{cMpc}$ Gaussian kernel before computing the
deformation tensor.

####  3.1.5 Error modeling

In the mock PFS survey, our galaxy redshift-space positions and projected
shapes depend on the viewing angle. Both Z1_FID and Z2_FID density
reconstructions also depend on the viewing angle: RSD is applied to the input
galaxies along the viewing angle line of sight, while for Z2_FID, Lyman-alpha
skewers are also created along each line of sight. In order to quantify the
effects of differing viewing angles, when we compute alignment statistics for
a sample or subsample of $N_{g}$ galaxies, we marginalize over all viewing
angles. For a given eigenvector, we therefore compute alignment statistics on
64$\times N$ galaxy-nearest eigenvector pairs.

Since the actual survey will only be on (approximately) one viewing angle but
over a larger overall volume than that of TNG300, we adopt a modified
bootstrap sampling scheme to quantify the viewing angle variance in addition
to our usual estimate of cosmic variance. For each bootstrap iteration, we
perform the following procedure. First, we randomly select one viewing angle,
and consider only the galaxy-eigenvector pairs which use that viewing angle,
for RSD, density reconstruction, and shape projection (if applicable). We then
perform normal bootstrap sampling on the pairs, and compute a histogram on
$\lvert\cos{\theta}\rvert$ along with
$\langle\lvert\cos{\theta}\rvert{}\rangle$. The 16th and 84th percentiles on
the computed histogram bins and $\langle\lvert\cos{\theta}\rvert{}\rangle$
from 10,000 bootstrap sampling iterations (with a random viewing angle for
each iteration) are then adopted as our error bars.

As mentioned in Section 3.1.2, our pool of 94,766 possible survey targets
within the TNG-300-1 volume is not large enough to sample the $N_{g}$ =
250,000 galaxies expected for the $z=1$ survey planned for Subaru PFS GE. We
therefore perform bootstrap error estimation on all survey targets, and scale
the resulting errorbars by $\sqrt{94,766/250,000}\simeq 0.615$. This is an
approximate scaling assumption based on the Central Limit Theorem’s expression
for sample standard deviation, which scales as 1/$\sqrt{N}$ with a sample size
$N$.

In Section 3.1.1, we found that the variance in projected alignment across
differing viewing angles scales as $\sigma_{VA}\propto V_{S}^{-1/2}$, where
$V_{S}$ is the survey volume. We therefore additionally scale the errorbars of
the projected alignment only by
$\sqrt{205^{3}h^{-3}\,\mathrm{cMpc}^{3}{}/V_{S}}$, where $V_{S}$ is the
expected volume of the survey. For $z=1$,
$V_{S}\sim$32\,467\,532.467\,532\,467$$ $h^{-3}\,\mathrm{cMpc}^{3}{}$, and
$V_{S}\sim$2.7\text{\times}{10}^{7}$$ $h^{-3}\,\mathrm{cMpc}^{3}{}$ for $z=2$.
This yields an error scaling factor of 0.515 and 0.565 respectively for $z=1$
and $z=2$.

###  3.2 Results and Discussion

The first step of validating the results is to test the recovery of the cosmic
web from the density reconstructions described in the previous subsections.
Since this largely repeats similar exercises previously done in Krolewski et
al. (2017b) and Horowitz et al. (2021), we move this discussion to Appendix C.
In all cases, the cosmic web eigenvectors derived from the reconstructed
density fields are in excellent agreement with the underlying ‘truth’, which
sets us up for the subsequent analysis.

Sample | Mass threshold | $N_{gal}$ | 3D $\sqrt{\chi^{2}}$ | Projected $\sqrt{\chi^{2}}$
---|---|---|---|---
Z1_FID | N/A | 250000* | $5.426\,82$ | $5.321\,684\,971\,945\,624$
| $\leq 10^{10.1}h^{-1}\,M_{\odot}{}$ | 122056* | $1.170\,19$ | $3.491\,704\,810\,046\,807$
| $>10^{10.1}h^{-1}\,M_{\odot}{}$ | 127944* | $6.805\,59$ | $4.520\,962\,629\,900\,436$
Z2_FID | N/A | 15400 | $0.944\,623\,508\,444\,108\,1$ | $1.304\,733\,088\,140\,860\,5$
| $\leq 10^{10.5}h^{-1}\,M_{\odot}{}$ | 11264 | $0.372\,337\,404\,368\,316\,1$ | $0.884\,183\,649\,905\,495\,7$
| $>10^{10.5}h^{-1}\,M_{\odot}{}$ | 4136 | $1.056\,814\,897\,993\,778$ | $1.188\,308\,650\,988\,831\,6$
Z2_2X | N/A | 30800 | $1.498\,047\,191\,147\,385\,2$ | $1.471\,300\,765\,505\,908$
| $\leq 10^{10.5}h^{-1}\,M_{\odot}{}$ | 22652 | $0.760\,527\,297\,689\,611\,6$ | $0.946\,938\,045\,832\,180\,2$
| $>10^{10.5}h^{-1}\,M_{\odot}{}$ | 8148 | $1.574\,819\,455\,352\,614\,3$ | $1.428\,487\,144\,096\,319\,2$
Table 2: 3D and projected significances ($\sqrt{\chi^{2}}$) for alignments
between galaxy shapes and cosmic web eigenvectors calculated from the mock
observational surveys, at $z=1$ and $z=2$. We note that $N_{gal}$ is the
number of galaxies that PFS is predicted to observe for each redshift. For
$z=1$, the stated $N_{gal}$ (marked with asterisks) is $\sim$2.6x the number
of galaxies actually used in our analysis, due to volume limitations on
TNG300-1. We therefore apply additional error scaling to compensate, as
detailed in Section 3.1.5.

####  3.2.1 $z=1$ mock observational survey alignment

In Figure 6, we present results for our mock survey alignment analysis for the
Z1_FID galaxy sample at $z=1$. In order to isolate the various observational
effects on any potential measurement, we measure the alignment relative to
both the ‘true’ cosmic web, as well as that derived from the reconstructed
mock density field.

In the left half of each subpanel, we show the mean alignment,
$\langle\lvert\cos{\theta}\rvert{}\rangle$, of the galaxies relative to the
cosmic web eigenvectors computed from the true density field. The right half
of each subpanel shows $\langle\lvert\cos{\theta}\rvert{}\rangle$ computed
from eigenvectors in the TARDIS reconstructed density fields, that were
applied to realistic mock data as detailed in Section 3.1.3. We hereafter
refer to this as the “mock reconstruction” alignment. In both cases, all
measurements are calculated in redshift space.

The top row of panels in the Figure 6 is the
$\langle\lvert\cos{\theta}\rvert{}\rangle$ computed for the 3D galaxy shapes,
while the bottom row is the $\langle\lvert\cos{\theta}\rvert{}\rangle$ from
the projected (2D) galaxy shapes and eigenvectors calculated over different
viewing angle planes. In every subpanel, a dashed black line represents the
$\langle\lvert\cos{\theta}\rvert{}\rangle$ from the null distribution. As
described in Section 3.1.5, we marginalize over all viewing angles to derive
errors for $\langle\lvert\cos{\theta}\rvert{}\rangle$, and additionally scale
our errors to correct for galaxy number counts (for both 3D and projected
alignment) and for viewing-angle variance (for projected alignment only). For
the alignments computed from the true density field, the differing redshift-
space galaxy displacements along different viewing angles are accounted for in
the bootstrap error estimation process; and for the mock reconstruction
alignment, the differing density reconstructions are accounted for as well. In
addition to presenting $\langle\lvert\cos{\theta}\rvert{}\rangle$ for the full
samples, we also calculate $\langle\lvert\cos{\theta}\rvert{}\rangle$ with the
sample split into low and high stellar mass samples. We set the threshold mass
to the transition masses found in Section 2.2.1, i.e. $M_{*}=[3.16\times
10^{10},\;1.26\times 10^{10}]\,h^{-1}\,M_{\odot}$ at $z=[1,2]$ respectively.
For the purpose of assessing observational feasibility, we use the stellar
mass as a rough proxy for brightness. While TNG includes synthetic magnitudes
for galaxies, we opt not to split our sample by magnitude, as we are unable to
incorporate the effects of dust extinction in our viewing-angle agnostic
approach.

In order to quantify the overall detection significance relative to the null
distribution, we adopt a metric based on the reduced chi-squared statistic. We
calculate the covariance of $\langle\lvert\cos{\theta}\rvert{}\rangle$ across
all three eigenvectors using bootstrap sampling, and define our significance
metric as:

$\sqrt{\chi^{2}}=\sqrt{(\vec{m}-\vec{n})^{T}C^{-1}(\vec{m}-\vec{n})}$ (4)

where $\vec{m}$ is the vector of $\langle\lvert\cos{\theta}\rvert{}\rangle$
with the elements populated by the values from the three eigenvectors,
$\vec{n}$ is a constant vector of the null distribution’s mean alignment, and
$C$ is the covariance matrix. For our analysis of Z1_FID, the covariance
matrix is computed using the rescaled 16th and 84th percentile bootstrap
errors computed for the mean alignments. We present mock reconstruction
alignment significances in Table 2.

Figure 6: Z1_FID $\langle\lvert\cos{\theta}\rvert{}\rangle$ summary boxplot
for true density field and mock reconstruction stages, for shape alignments.
3D mean alignments are shown on the top row, while projected mean alignments
are on the bottom. The null distribution’s
$\langle\lvert\cos{\theta}\rvert{}\rangle$ is marked with a dashed line for
both rows. The boxplot whiskers have been modified to display 16th and 84th
percentiles $\langle\lvert\cos{\theta}\rvert{}\rangle$.

In 3D, the alignments seen between the galaxy shapes and cosmic web
eigenvectors still preserve the underlying trends seen in Section 2.2. When
shape alignments are split by the transition mass found in Section 2.2.1, the
low-mass sample has near-null alignment as expected, while the high-mass
alignment is stronger than the total alignment. As TARDIS yields a high
fidelity reconstruction of the cosmic web eigenvectors in the volume, the mock
reconstruction alignment is nearly identical to true density field alignment,
with the only notable difference being a negligible
($\Delta$$\langle\lvert\cos{\theta}\rvert{}\rangle$$\sim$0.001) decrease in
mean alignments.

When the 3D alignments are projected into 2D, the general alignment trends as
a function of stellar mass remain the same as the 3D case. Projection tends to
increase the scatter of $\langle\lvert\cos{\theta}\rvert{}\rangle$ due to the
viewing-angle variance effect described in Section 3.1.1, but the volume
rescaling we apply to the projected $\langle\lvert\cos{\theta}\rvert{}\rangle$
appears to largely counteract the variation increase. We note that the viewing
angle projection also appears to introduce an error asymmetry into the mean
alignments for both $\boldsymbol{\hat{e}_{3}}$ and $\boldsymbol{\hat{e}_{1}}$,
with increased variation in the sign-direction of alignment (positive for
$\boldsymbol{\hat{e}_{3}}$ and negative for $\boldsymbol{\hat{e}_{1}}$). This
may be caused by a bias in eigenvector direction in the IllustrisTNG volume
used, arising from cosmic variance as explored in Section 3.1.1.

The effect of projection on the mock reconstruction alignment can be
quantified in terms of the significance metric $\sqrt{\chi^{2}}$ listed in
Table 2. Due to the volume error scaling, $\sqrt{\chi^{2}}$ is slightly
decreased by a factor of $\sim$0.66 for the high-mass samples, and remains
almost constant when considering the total (low-mass and high-mass) sample.
The 5.3$\sigma$ significance found in the galaxy shape alignments with the
reconstructed cosmic web in the Z1_FID sample, when combined over all stellar
masses, crosses the conventional detection significance threshold of
4$\sigma$, so we anticipate a successful detection in the PFS survey.

####  3.2.2 $z=2$ mock observational survey alignment

Figure 7: Z2_FID $\langle\lvert\cos{\theta}\rvert{}\rangle$ summary boxplot
for true density field and mock reconstruction stages, for shape alignments.
3D mean alignments are shown on the top row, while projected mean alignments
are on the bottom. The null distribution’s
$\langle\lvert\cos{\theta}\rvert{}\rangle$ is marked with a dashed line for
both rows. The boxplot whiskers have been modified to display 16th and 84th
percentiles $\langle\lvert\cos{\theta}\rvert{}\rangle$.

At $z=2$, we display the mean alignment results for only Z2_FID for visual
clarity, though we will later discuss mock reconstruction results for the
increased sample size (Z2_2X). In Figure 7, we show shape alignments for
Z2_FID. We find that even in 3D, the relatively small galaxy sample size
produces significant uncertainties in both samples compared to the Z1_FID
alignments. While the $\langle\lvert\cos{\theta}\rvert{}\rangle$ central
values for $\boldsymbol{\hat{e}_{3}}$ and $\boldsymbol{\hat{e}_{1}}$ are
comparable to those for Z1_FID, the central 68th percentile spread increases
from $\Delta$$\langle\lvert\cos{\theta}\rvert{}\rangle$ $\sim$ 0.002 to
$\Delta$$\langle\lvert\cos{\theta}\rvert{}\rangle$ $\sim$ 0.005. This increase
is also greater for the high-mass sample, as the proportion of high-mass
galaxies is $\sim$ 27% at $z=2$, as opposed to $\sim$ 50% for $z=1$.
Projection further compounds this alignment degradation, with central
$\langle\lvert\cos{\theta}\rvert{}\rangle$ dropping below 1$\sigma$ for
$\boldsymbol{\hat{e}_{3}}$ and $\boldsymbol{\hat{e}_{1}}$ in all cases. As
opposed to Z1_FID where high reconstruction quality has no visible effect on
3D alignment spread, the volume error scaling is unable to compensate for the
comparatively poorer reconstruction quality in Z2_FID. We therefore see a
small increase in the central 68th percentile spread of
$\Delta\langle\lvert\cos{\theta}\rvert{}\rangle{}\sim$ 0.004 for the mock
reconstruction alignment, compared to the true density field alignment. In
fact, the volume error scaling we apply may be detrimentally amplifying
variance in the mean alignment caused by the small sample. When considering
the mock reconstruction significance shown in Table 2, we find that the
projected significance is unexpectedly increased in comparison to the 3D
significance for Z2_FID. However, Z2_2X shows an equal significance for
projected and 3D mock reconstruction alignment, lending credence to the small-
sample amplification effect. For both sample sizes and all mass splits, the
mock reconstruction shape alignment significance is $<$1.5, so we do not
anticipate a significant detection of the IllustrisTNG shape-cosmic web
alignment at $z=2$ by PFS.

####  3.2.3 Horizon-AGN Observability Estimate

Due to the large difference in spin-cosmic web alignment strength between
IllustrisTNG and Horizon-AGN which we pointed out in Section 2.2.3, in this
section we aim to provide a rough estimate on the detectability of the
Horizon-AGN alignment signal.

As the Horizon-AGN mean alignment signal is only published for $z=1$ spin-
cosmic alignment by Codis et al. (2015), we first take the ratio between
Horizon-AGN and IllustrisTNG signal strength at $z=1$ $r_{\text{MA}}$ that we
have found in Section 2.2.3. This is $r^{*}_{\text{MA}}=2.4$. We find in
Appendix A that spin alignments and shape alignments in IllustrisTNG are
highly correlated, as expected. We therefore use the spin-alignment signal
ratio as a very rough guide for the expected detection significance of the
Horizon-AGN _shape_ -cosmic web alignment signal at both $z=1$ and $z=2$, by
multiplying the IllustrisTNG shape “$\langle\lvert\cos{\theta}\rvert{}\rangle$
over null” $\vec{m}-\vec{n}$ by the same ratio. As
$\sqrt{\chi^{2}}{}\propto|\vec{m}-\vec{n}|$, this is equivalent to multiplying
$\sqrt{\chi^{2}}$ by that factor. This assumes that the significance
covariance matrix (see Equation 4) that we have calculated for IllustrisTNG
alignments remains exactly the same if using Horizon-AGN. We believe this to
be an appropriate conservative assumption, as we find in Figure 3 that the
error bars for Horizon-AGN are much smaller than for the corresponding
IllustrisTNG sample (though note that Horizon-AGN uses Poisson errors, while
we use bootstrap errors for TNG).

Following this procedure, we would expect the Horizon-AGN alignment signal to
be detectable at $\sim 3\sigma$ at both redshift bins, with the lowest total
significance of Z2_FID improving from $1.3\sigma$ to $3.1\sigma$ and the
Z1_FID total significance increasing from $5.3\sigma$ to $12.7\sigma$.
However, we note again that this is a very rough estimate, but does suggest
that the Subaru PFS Galaxy Evolution Survey would be able to strongly detect
or rule out the level of galaxy-cosmic web alignment predicted by Horizon-AGN.

####  3.2.4 Possible survey improvements

Assuming that detection significance scales as $\propto\sqrt{N_{gal}}$, we
estimate that in order to constrain $z=2$ IllustrisTNG shape alignments with a
significance of $\sqrt{\chi^{2}}{}=4$, next-generation surveys must be able to
map the $z\sim 2-3$ Ly$\alpha$ forest over a $\sim$10 times larger volume than
fiducial Z2_FID survey we have defined for Subaru PFS, i.e. $V_{S}\sim 3\times
10^{8}\,h^{-3}\,\mathrm{cMpc}^{3}$. Assuming a similar redshift range and
target density as the PFS IGM tomography survey, this implies an survey area
of $\sim 100\,\mathrm{deg}^{2}$, which is likely achievable in the next
generation of large-scale spectroscopic surveys on large telescopes such as
the Maunakea Spectroscopic Explorer (McConnachie et al., 2016) or MegaMapper
(Schlegel et al., 2022) concepts. This also assumes the same magnitude limit
and reconstruction quality as PFS, so this is a conservative estimate as
future surveys will likely be able to improve on both avenues. In addition, as
we find that projection from 3D onto 2D significantly degrades alignment
significance, improved efforts to estimate 3D shape would provide immediate
measured significance improvements.

##  4 Conclusion

In this paper, we analyzed the possible alignments between the large-scale
($\gtrsim$ Mpc) cosmic web and galaxy shapes/spins within the IllustrisTNG
simulation at $z=1$ and $z=2$, with the goal of predicting whether they might
be observable with the upcoming Galaxy Evolution Survey to be carried out on
the Subaru PFS massively-multiplexed fiber spectrograph.

We measure the relative 3D alignment between galaxy shapes and spins to the
underlying cosmic web eigenvectors of the matter density field calculated
using the pseudo-deformational tensor approach. In the TNG300 simulation,
above a transition stellar mass of $M_{*}>10^{10.5}\,M_{\odot}$ we find that a
clear excess alignment between galaxy shapes and the
$\boldsymbol{\hat{e}_{3}}$ eigenvector, which typically traces cosmic web
filaments. This alignment is present at comparable strengths at both $z=1$ and
$z=2$, although the relative lack of high stellar mass galaxies at $z=2$
reduces the statistical significance of our finding at that earlier epoch. The
$M_{*}>10^{10.2}\,M_{\odot}$ galaxy spins show a weak anti-alignment trend
with respect to $\boldsymbol{\hat{e}_{3}}$ at $z=1$, although the trend is
difficult to detect clearly for lower-mass galaxies and at $z=2$. We find that
these alignment trends are especially weak compared with the Horizon-AGN
simulation, with Horizon-AGN alignment being $\sim 5-6\sigma$ stronger than
IllustrisTNG in the most conservative case. This is surprising given that
Horizon-AGN is relatively contemporary with IllustrisTNG, indicating that
galaxy-cosmic web alignment is a sensitive probe of the physics of galaxy
formation and evolution. A detailed analysis of what process leads to this
divergence is a promising future direction. We also compare with the smaller,
but higher-resolution, TNG100 simulations where we generally find consistent
alignment trends, although statistical uncertainties dominate in this smaller
volume.

We proceed to assess the observability of the this alignment signal in the
Galaxy Evolution Survey planned to be carried out with the upcoming Subaru
Prime Focus Spectrograph, which will map the cosmic web with a high-number
density galaxy redshift survey at $0.7<z<1.7$ as well as an IGM tomography
survey of the Lyman-$\alpha$ forest at $2.2<z<2.7$. We project the galaxy
morphologies onto a 2D plane to mimic the observational shapes, but within a
volume the size of the TNG300 box we find that the projected alignment signal
can vary depending on the direction from which the observer is viewing the
data. We find that this viewing angle variance is a form of cosmic variance
scaling with the observed volume, which we attribute to anisotropic large-
scale structures in the simulation volume being viewed from different
directions. This effect suggests even the $L=300$ cMpc TNG300 volume is not
quite large enough to provide a fair sampling of the cosmic web alignments —
future theoretical studies should ideally use larger volumes with $L\gtrsim 1$
cGpc. Equivalently, this additional cosmic variance should be taken into
account in future observational analyses of galaxy alignments with the cosmic
web, although the cosmic volumes mapped by the Subaru PFS surveys will be
$\sim 3-4\times$ the size of the TNG300 volume.

Given the survey volume and number of galaxies expected to be observed in the
lower-redshift ($z\sim 1$) PFS galaxies, we expect the shape alignment signal
in TNG300 to be detectable at $5.3\,\sigma$ significance. At $z\sim 2$, the
signal will be more challenging to detect, rising only to $1.3\,\sigma$
significance above a random distribution. This is based on the current
baseline survey plan for $\sim 15,000$ coeval galaxies within the volume of
the IGM tomographic map, but a possible survey plan that doubles the number of
coeval galaxies would only yield the expected $\sqrt{2}\sim 1.4\times$ boost.

We estimate, however, that with the stronger alignment seen in previous
analyses of the Horizon-AGN simulation, there should a very strong alignment
signal detectable in the PFS $z=1$ samples and even potentially at $z=2$.

However, the main obstacle to a real observational measurement of galaxy-
cosmic web alignments on Subaru PFS is not the galaxy sample sizes nor cosmic
web reconstruction, but in the paucity of resolved near infrared imaging of
galaxy morphologies within the proposed PFS footprint at the time of writing.
The typical effective radii of quiescent galaxies ranges from $R_{e}\sim 1-2$
pkpc over $z\sim 1-2$, which correspond to $\theta\sim 0.2-0.3^{\prime\prime}$
on the sky. Meanwhile, to cover the quiescent stellar distribution redward of
the restframe $\lambda\approx 4000\,\AA$ break for $z>1.5$ galaxies would
require filter coverage at wavelengths beyond $1$ micron wavelengths, i.e.
into the near-infrared. This therefore calls for space-based $H$-band
($\lambda_{C}\sim 1.6\,\mu$m) imaging with the Hubble Space Telescope (HST) to
resolve stellar morphologies of galaxies up to $z\approx 2.7$. Such imaging
currently exists only over 1.4 deg2 in the COSMOS field thanks to the recent
HST-DASH survey, while another aggregate of $\sim 1$ deg2 exists taking into
account all other HST imaging. This currently falls far short of the 12.3 deg2
that the PFS GE survey will cover, although the Roman Space Telescope will
efficiently allow this full footprint to be imaged once it is launched after
2029 — some time after the completion of the Subaru PFS Galaxy Evolution
Survey.

However, once all the data is in place to allow constraints on the galaxy-
cosmic web alignment signals at $z\sim 1-3$, one can expect this quantity to
become an interesting probe of galaxy formation and evolution.

##  5 Acknowledgments

We thank Viraj Pandya for useful discussions that helped initiate this
project. Kavli IPMU was established by World Premier International Research
Center Initiative (WPI), MEXT, Japan. K.G.L. acknowledges support from JSPS
Kakenhi grants JP18H05868 and JP19K14755. B.H. is supported by the AI
Accelerator program of the Schmidt Futures Foundation. 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 authors thank the Yukawa Institute for Theoretical Physics
at Kyoto University. Discussions during the YITP workshop YITP-T-21-06 on
“Galaxy shape statistics and cosmology” were useful to complete this work.

##  Appendix A TNG Shape-Spin Alignment

Figure 8: Cumulative density function of alignment between the spin and three
shape ellipsoid axes of TNG300-1 and TNG100-1 galaxies, at $z=1$. The TNG
galaxy sample is the same as Figure 3, with a minimum-mass cut of
$M_{*}>10^{8.06}$ $h^{-1}\,M_{\odot}$ to match Codis et al. (2015).

In Figure 8, we plot the cumulative density function of alignment between the
spin and three shape (inertia tensor) ellipsoid axes of TNG300-1 and TNG100-1
galaxies at $z=1$, without mock observational effects. We use the same galaxy
sample as Section 2.2.3, which aims to match the Horizon-AGN galaxy sample
used in Codis et al. (2015) by imposing a minimum stellar-mass cut of
$M_{*}>10^{8.06}$ $h^{-1}\,M_{\odot}$.

Comparing to the corresponding CDF plotted in Figure 2 of Codis et al. (2015),
we see close qualitative agreement, with the minor shape axis closely aligned
with the spin axis of galaxies. The slightly weaker degree of alignment
between shape and spin in TNG300-1 compared to TNG100-1 may be due to worse
mass resolution in the former simulation, as it persists across all mass bins.

The high level of alignment between the minor shape axis and spin axis of
galaxies suggests that the much weaker degree of spin alignment seen in TNG
compared to Horizon-AGN extends to shape alignment.

##  Appendix B Alignment Dependence on Smoothing Scale

Figure 9: TNG300-1 $\boldsymbol{\hat{e}_{3}}{}$-shape mean alignment, at
$z=1$, as a function of density smoothing scale. We use the magnitude-selected
Z1_FID galaxy sample, but do not apply any mock observational effects. Errors
are derived from bootstrap resampling (but do not include any viewing-angle
effects).

As cosmic web filaments exist on multiple scales, the degree of shape or spin
alignment observed has a large dependence on the scale (Gaussian kernel size)
at which the density is smoothed before calculating the deformation tensor.
Lee et al. (2020) find that as density smoothing scale increases, the
magnitude of alignments decreases, which is corroborated in our analysis. In
Figure 9, we plot shape-$\boldsymbol{\hat{e}_{3}}$ alignment as a function of
density smoothing scale. We use the Z1_FID galaxy sample outlined in Section
3.1.2, but do not apply any mock observational effects described in the
section.

Our choice of a 2 $h^{-1}\,\mathrm{cMpc}$ smoothing scale is motivated by the
maximum reconstruction resolution we expect from Ly$\alpha$ tomography; as
described in Lee et al. (2014), this is set by the average transverse
Ly$\alpha$ sightline spacing, which for PFS is $\sim 2$
$h^{-1}\,\mathrm{cMpc}$. However, as Figure 9 demonstrates, choosing a
smoothing scale of $\sim 7$ $h^{-1}\,\mathrm{cMpc}$ would have cut the
magnitude of alignment by a factor of two. We therefore recommend that
analyses of alignment from PFS observational data adopt a cosmic web formalism
that takes into account its multi-scale nature: for instance, by using the
NEXUS cosmic web formalism (Cautun et al., 2013, 2014), which extends the
deformation tensor approach we use to multiple scales; or the Multiscale
Morphology Filter (Aragón-Calvo et al., 2007), a similarly multi-scale
algorithm which is based on the deformation tensor formalism.

##  Appendix C Mock Reconstruction Quality

Figure 10: Probability density function of $\lvert\cos{\theta}\rvert$ between
true and TARDIS-reconstructed eigenvectors, for each point in the
reconstructed IllustrisTNG volume. Along with the $z=1$ and joint $z=2$
reconstructions used for the final analysis, we include a Lyman-alpha forest-
only reconstruction at $z=2$ for comparison. Each viewing angle’s
reconstruction is separately plotted, as a semi-transparent line. In the case
of a perfect reconstruction, the PDF would be a delta-function at
$\lvert\cos{\theta}\rvert{}=1$ for each of the panels.

In this appendix, we validate the quality of the TARDIS density
reconstructions on the galaxy redshift-space distribution of the Z1_FID sample
at $z=1$, as well as the joint reconstructions of the galaxy distribution as
well as Ly$\alpha$ forest absorption in the Z2_FID and Z2_2X mock samples at
$z\sim 2$. For this purpose, we compute the eigenvectors of the pseudo-
deformation tensor (as described in Section 2.1.3) for both the reconstructed
volumes and true density field. Figure 10 illustrates the accuracy of the
eigenvector recovery by plotting the PDF of $\lvert\cos{\theta}\rvert$ between
true and reconstructed eigenvectors at each grid point in the volumes, for all
three of our mock observational survey configurations. We also include the
performance of a $z=2$ reconstruction using only Ly$\alpha$ tomography, to
provide a baseline for the effect of incorporating the galaxy sample on
reconstruction quality. In order to display variance in reconstruction quality
across random viewing angles for a particular survey configuration, each
viewing angle’s PDF is plotted separately as a semi-transparent line.

We find that for all survey configurations, $\boldsymbol{\hat{e}_{3}}$ and
$\boldsymbol{\hat{e}_{1}}$ are reconstructed more accurately than
$\boldsymbol{\hat{e}_{2}}$. This is consistent with results reported by
Horowitz et al. (2021). The $z=1$ galaxy-only reconstruction is significantly
more accurate than both $z=2$ joint reconstructions, which are in turn both
more accurate than the Lyman-$\alpha$ only reconstruction. Doubling the galaxy
density for $z=2$ yields a modest improvement, pushing the reconstruction
quality closer to that of $z=1$. Overall reconstruction quality remains
consistent across viewing angles for all survey configurations. However, we
note that this does not imply that for a particular grid point in the volume,
the reconstructed eigenvector will be as consistent. We find that for
individual grid points, the reconstructed-true eigenvector alignment follows
the same distribution as the overall reconstruction quality’s, when comparing
across viewing angles.

Following the cosmic web classification scheme described in Forero-Romero et
al. (2009) and Lee & White (2016) where the cosmic web is divided into
nodes/sheets/filaments/voids using the eigenvalues of the deformation tensor,
we also find that our classification accuracy is consistent with the accuracy
reported by Horowitz et al. (2021).

## References

  * Altay et al. (2006) Altay, G., Colberg, J. M., & Croft, R. A. C. 2006, MNRAS, 370, 1422, doi: 10.1111/j.1365-2966.2006.10555.x
  * Aragón-Calvo et al. (2007) Aragón-Calvo, M. A., Jones, B. J. T., van de Weygaert, R., & van der Hulst, J. M. 2007, A&A, 474, 315, doi: 10.1051/0004-6361:20077880
  * Aragón-Calvo et al. (2007) Aragón-Calvo, M. A., van de Weygaert, R., Jones, B. J. T., & van der Hulst, J. M. 2007, ApJl, 655, L5, doi: 10.1086/511633
  * Aragon-Calvo & Yang (2014) Aragon-Calvo, M. A., & Yang, L. F. 2014, MNRAS, 440, L46, doi: 10.1093MNRASl/slu009
  * Bailin & Steinmetz (2005) Bailin, J., & Steinmetz, M. 2005, ApJ, 627, 647, doi: 10.1086/430397
  * Bate et al. (2020) Bate, J., Chisari, N. E., Codis, S., et al. 2020, MNRAS, 491, 4057, doi: 10.1093/mnras/stz3166
  * Bird et al. (2015) Bird, S., Haehnelt, M., Neeleman, M., et al. 2015, MNRAS, 447, 1834, doi: 10.1093/mnras/stu2542
  * Blazek et al. (2019) Blazek, J. A., MacCrann, N., Troxel, M. A., & Fang, X. 2019, Phys. Rev. D, 100, 103506, doi: 10.1103/PhysRevD.100.103506
  * Blue Bird et al. (2020) Blue Bird, J., Davis, J., Luber, N., et al. 2020, MNRAS, 492, 153, doi: 10.1093/mnras/stz3357
  * Blumenthal et al. (1984) Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517, doi: 10.1038/311517a0
  * Bond et al. (1996) Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603, doi: 10.1038/380603a0
  * Bridle & King (2007) Bridle, S., & King, L. 2007, New Journal of Physics, 9, 444, doi: 10.1088/1367-2630/9/12/444
  * Brunino et al. (2007) Brunino, R., Trujillo, I., Pearce, F. R., & Thomas, P. A. 2007, MNRAS, 375, 184, doi: 10.1111/j.1365-2966.2006.11282.x
  * Caucci et al. (2008) Caucci, S., Colombi, S., Pichon, C., et al. 2008, MNRAS, 386, 211, doi: 10.1111/j.1365-2966.2008.13016.x
  * Cautun et al. (2013) Cautun, M., van de Weygaert, R., & Jones, B. J. T. 2013, MNRAS, 429, 1286, doi: 10.1093/mnras/sts416
  * Cautun et al. (2014) Cautun, M., van de Weygaert, R., Jones, B. J. T., & Frenk, C. S. 2014, MNRAS, 441, 2923, doi: 10.1093MNRAS/stu768
  * Chen et al. (2015) Chen, Y.-C., Ho, S., Tenneti, A., et al. 2015, MNRAS, 454, 3341, doi: 10.1093MNRAS/stv2260
  * Codis et al. (2018) Codis, S., Jindal, A., Chisari, N. E., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 4753, doi: 10.1093/mnras/sty2567
  * Codis et al. (2012) Codis, S., Pichon, C., Devriendt, J., et al. 2012, MNRAS, 427, 3320, doi: 10.1111/j.1365-2966.2012.21636.x
  * Codis et al. (2015) Codis, S., Gavazzi, R., Dubois, Y., et al. 2015, MNRAS, 448, 3391, doi: 10.1093/mnras/stv231
  * Cuesta et al. (2008) Cuesta, A. J., Betancort-Rijo, J. E., Gottlöber, S., et al. 2008, MNRAS, 385, 867, doi: 10.1111/j.1365-2966.2008.12879.x
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371, doi: 10.1086/163168
  * Dubois et al. (2014) Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453, doi: 10.1093/mnras/stu1227
  * Einasto et al. (1984) Einasto, J., Klypin, A. A., Saar, E., & Shandarin, S. F. 1984, MNRAS, 206, 529, doi: 10.1093MNRAS/206.3.529
  * Forero-Romero et al. (2014) Forero-Romero, J. E., Contreras, S., & Padilla, N. 2014, MNRAS, 443, 1090, doi: 10.1093MNRAS/stu1150
  * Forero-Romero et al. (2009) Forero-Romero, J. E., Hoffman, Y., Gottlöber, S., Klypin, A., & Yepes, G. 2009, MNRAS, 396, 1815, doi: 10.1111/j.1365-2966.2009.14885.x
  * Förster Schreiber et al. (2019) Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21, doi: 10.3847/1538-4357/ab0ca2
  * Ganeshaiah Veena et al. (2019) Ganeshaiah Veena, P., Cautun, M., Tempel, E., van de Weygaert, R., & Frenk, C. S. 2019, MNRAS, 487, 1607, doi: 10.1093/mnras/stz1343
  * Ganeshaiah Veena et al. (2021) Ganeshaiah Veena, P., Cautun, M., van de Weygaert, R., Tempel, E., & Frenk, C. S. 2021, MNRAS, 503, 2280, doi: 10.1093/mnras/stab411
  * Ganeshaiah Veena et al. (2018) Ganeshaiah Veena, P., Cautun, M., van de Weygaert, R., et al. 2018, MNRAS, 481, 414, doi: 10.1093/mnras/sty2270
  * Geller & Huchra (1989) Geller, M. J., & Huchra, J. P. 1989, Science, 246, 897, doi: 10.1126/science.246.4932.897
  * Greene et al. (2022) Greene, J., Bezanson, R., Ouchi, M., Silverman, J., & the PFS Galaxy Evolution Working Group. 2022, arXiv e-prints, arXiv:2206.14908.  https://arxiv.org/abs/2206.14908
  * Gunn & Gott (1972) Gunn, J. E., & Gott, J. Richard, I. 1972, ApJ, 176, 1, doi: 10.1086/151605
  * Guzzo et al. (2014) Guzzo, L., Scodeggio, M., Garilli, B., et al. 2014, A&A, 566, A108, doi: 10.1051/0004-6361/201321489
  * Hahn et al. (2007) Hahn, O., Carollo, C. M., Porciani, C., & Dekel, A. 2007, MNRAS, 381, 41, doi: 10.1111/j.1365-2966.2007.12249.x
  * Hahn et al. (2007) Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489, doi: 10.1111/j.1365-2966.2006.11318.x
  * Hannay & Nye (2004) Hannay, J. H., & Nye, J. F. 2004, Journal of Physics A: Mathematical and General, 37, 11591, doi: 10.1088/0305-4470/37/48/005
  * Haynes & Giovanelli (1984) Haynes, M. P., & Giovanelli, R. 1984, AJ, 89, 758, doi: 10.1086/113573
  * Horowitz et al. (2019) Horowitz, B., Lee, K.-G., White, M., Krolewski, A., & Ata, M. 2019, ApJ, 887, 61, doi: 10.3847/1538-4357/ab4d4c
  * Horowitz et al. (2021) Horowitz, B., Zhang, B., Lee, K.-G., & Kooistra, R. 2021, ApJ, 906, 110, doi: 10.3847/1538-4357/abca35
  * Horowitz et al. (2022) Horowitz, B., Lee, K.-G., Ata, M., et al. 2022, ApJS, 263, 27, doi: 10.3847/1538-4365/ac982d
  * Ishikawa et al. (2017) Ishikawa, S., Kashikawa, N., Toshikawa, J., et al. 2017, ApJ, 841, 8, doi: 10.3847/1538-4357/aa6d64
  * Jasche et al. (2015) Jasche, J., Leclercq, F., & Wandelt, B. D. 2015, J. Cosmology Astropart. Phys, 2015, 036, doi: 10.1088/1475-7516/2015/01/036
  * Kiessling et al. (2015) Kiessling, A., Cacciato, M., Joachimi, B., et al. 2015, Space Sci. Rev., 193, 67, doi: 10.1007/s11214-015-0203-6
  * Kirk et al. (2012) Kirk, D., Rassat, A., Host, O., & Bridle, S. 2012, MNRAS, 424, 1647, doi: 10.1111/j.1365-2966.2012.21099.x
  * Kirk et al. (2015) Kirk, D., Brown, M. L., Hoekstra, H., et al. 2015, \ssr, 193, 139, doi: 10.1007/s11214-015-0213-4
  * Kitaura et al. (2009) Kitaura, F. S., Jasche, J., Li, C., et al. 2009, MNRAS, 400, 183, doi: 10.1111/j.1365-2966.2009.15470.x
  * Klypin & Shandarin (1983) Klypin, A. A., & Shandarin, S. F. 1983, MNRAS, 204, 891, doi: 10.1093MNRAS/204.3.891
  * Kraljic et al. (2020) Kraljic, K., Davé, R., & Pichon, C. 2020, MNRAS, 493, 362, doi: 10.1093/mnras/staa250
  * Kraljic et al. (2021) Kraljic, K., Duckworth, C., Tojeiro, R., et al. 2021, MNRAS, 504, 4626, doi: 10.1093/mnras/stab1109
  * Krause et al. (2016) Krause, E., Eifler, T., & Blazek, J. 2016, MNRAS, 456, 207, doi: 10.1093/mnras/stv2615
  * Krolewski et al. (2019) Krolewski, A., Ho, S., Chen, Y.-C., et al. 2019, ApJ, 876, 52, doi: 10.3847/1538-4357/ab1010
  * Krolewski et al. (2017a) Krolewski, A., Lee, K.-G., Lukić, Z., & White, M. 2017a, ApJ, 837, 31, doi: 10.3847/1538-4357/837/1/31
  * Krolewski et al. (2017b) —. 2017b, ApJ, 837, 31, doi: 10.3847/1538-4357/837/1/31
  * Lamman et al. (2022) Lamman, C., Eisenstein, D., Aguilar, J. N., et al. 2022, arXiv e-prints, arXiv:2209.03949.  https://arxiv.org/abs/2209.03949
  * Lee et al. (2020) Lee, J., Libeskind, N. I., & Ryu, S. 2020, ApJ, 898, L27, doi: 10.3847/2041-8213/aba2ee
  * Lee et al. (2014) Lee, K.-G., Hennawi, J. F., White, M., Croft, R. A. C., & Ozbek, M. 2014, ApJ, 788, 49, doi: 10.1088/0004-637X/788/1/49
  * Lee & White (2016) Lee, K.-G., & White, M. 2016, ArXiv e-prints.  http://arxiv.org/abs/1603.04441
  * Lee & White (2016) Lee, K.-G., & White, M. 2016, ApJ, 831, 181, doi: 10.3847/0004-637X/831/2/181
  * Lee et al. (2014) Lee, K.-G., Hennawi, J. F., Stark, C., et al. 2014, ApJ, 795, L12, doi: 10.1088/2041-8205/795/1/L12
  * Lee et al. (2018) Lee, K.-G., Krolewski, A., White, M., et al. 2018, ApJS, 237, 31, doi: 10.3847/1538-4365/aace58
  * Libeskind et al. (2013) Libeskind, N. I., Hoffman, Y., Forero-Romero, J., et al. 2013, MNRAS, 428, 2489, doi: 10.1093/mnras/sts216
  * López et al. (2019) López, P., Merchán, M. E., & Paz, D. J. 2019, MNRAS, 485, 5244, doi: 10.1093/mnras/stz762
  * Malavasi et al. (2016) Malavasi, N., Arnouts, S., Vibert, D., et al. 2016, ArXiv e-prints.  https://arxiv.org/abs/1611.07045
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113, doi: 10.1093/mnras/sty2206
  * McConnachie et al. (2016) McConnachie, A. W., Babusiaux, C., Balogh, M., et al. 2016, ArXiv e-prints.  https://arxiv.org/abs/1606.00060
  * McCracken et al. (2015) McCracken, H. J., Wolk, M., Colombi, S., et al. 2015, MNRAS, 449, 901, doi: 10.1093/mnras/stv305
  * Momose et al. (2022) Momose, R., Lee, K.-G., Horowitz, B., Ata, M., & Kartaltepe, J. S. 2022, arXiv e-prints, arXiv:2212.05984, doi: 10.48550/arXiv.2212.05984
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206, doi: 10.1093/mnras/sty618
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624, doi: 10.1093/mnras/stx3040
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Newman et al. (2020) Newman, A. B., Rudie, G. C., Blanc, G. A., et al. 2020, ApJ, 891, 147, doi: 10.3847/1538-4357/ab75ee
  * Pahwa et al. (2016) Pahwa, I., Libeskind, N. I., Tempel, E., et al. 2016, MNRAS, 457, 695, doi: 10.1093MNRAS/stv2930
  * Pandya et al. (2019) Pandya, V., Primack, J., Behroozi, P., et al. 2019, MNRAS, 488, 5580, doi: 10.1093/mnras/stz2129
  * Patiri et al. (2006) Patiri, S. G., Cuesta, A. J., Prada, F., Betancort-Rijo, J., & Klypin, A. 2006, ApJl, 652, L75, doi: 10.1086/510330
  * Pichon et al. (2011) Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493, doi: 10.1111/j.1365-2966.2011.19640.x
  * Pichon et al. (2001) Pichon, C., Vergely, J. L., Rollinde, E., Colombi, S., & Petitjean, P. 2001, MNRAS, 326, 597, doi: 10.1046/j.1365-8711.2001.04595.x
  * Pillepich et al. (2018) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648, doi: 10.1093/mnras/stx3112
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Schlegel et al. (2022) Schlegel, D. J., Kollmeier, J. A., Aldering, G., et al. 2022, arXiv e-prints, arXiv:2209.04322.  https://arxiv.org/abs/2209.04322
  * Shi et al. (2021) Shi, J., Kurita, T., Takada, M., et al. 2021, J. Cosmology Astropart. Phys, 2021, 030, doi: 10.1088/1475-7516/2021/03/030
  * Sousbie (2011) Sousbie, T. 2011, MNRAS, 414, 350, doi: 10.1111/j.1365-2966.2011.18394.x
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726, doi: 10.1046/j.1365-8711.2001.04912.x
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Sugai et al. (2015) Sugai, H., Tamura, N., Karoji, H., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 035001, doi: 10.1117/1.JATIS.1.3.035001
  * Tempel & Libeskind (2013) Tempel, E., & Libeskind, N. I. 2013, ApJl, 775, L42, doi: 10.1088/2041-8205/775/2/L42
  * Tempel et al. (2013) Tempel, E., Stoica, R. S., & Saar, E. 2013, MNRAS, 428, 1827, doi: 10.1093MNRAS/sts162
  * Tenneti et al. (2015) Tenneti, A., Singh, S., Mandelbaum, R., et al. 2015, MNRAS, 448, 3522, doi: 10.1093/mnras/stv272
  * Tomassetti et al. (2016) Tomassetti, M., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 458, 4477, doi: 10.1093/mnras/stw606
  * Trowland et al. (2013) Trowland, H. E., Lewis, G. F., & Bland-Hawthorn, J. 2013, ApJ, 762, 72, doi: 10.1088/0004-637X/762/2/72
  * Wang et al. (2014) Wang, H., Mo, H. J., Yang, X., Jing, Y. P., & Lin, W. P. 2014, ApJ, 794, 94, doi: 10.1088/0004-637X/794/1/94
  * Wang et al. (2018a) Wang, P., Guo, Q., Kang, X., & Libeskind, N. I. 2018a, ApJ, 866, 138, doi: 10.3847/1538-4357/aae20f
  * Wang et al. (2018b) —. 2018b, ApJ, 866, 138, doi: 10.3847/1538-4357/aae20f
  * Wang & Kang (2018) Wang, P., & Kang, X. 2018, MNRAS, 473, 1562, doi: 10.1093/mnras/stx2466
  * Wechsler et al. (2002) Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52, doi: 10.1086/338765
  * Welker et al. (2020) Welker, C., Bland-Hawthorn, J., van de Sande, J., et al. 2020, MNRAS, 491, 2864, doi: 10.1093/mnras/stz2860
  * whuber (2020) whuber. 2020, Distribution of scalar products of two random unit vectors in $D$ dimensions, Cross Validated
  * Zeldovich et al. (1982) Zeldovich, I. B., Einasto, J., & Shandarin, S. F. 1982, Nature, 300, 407, doi: 10.1038/300407a0
  * Zel’Dovich (1970) Zel’Dovich, Y. B. 1970, A&A, 500, 13
  * Zhang et al. (2009) Zhang, Y., Yang, X., Faltenbacher, A., et al. 2009, ApJ, 706, 747, doi: 10.1088/0004-637X/706/1/747
  * Zhang et al. (2009) Zhang, Y., Yang, X., Faltenbacher, A., et al. 2009, ApJ, 706, 747, doi: 10.1088/0004-637X/706/1/747
  * Zhang et al. (2013) Zhang, Y., Yang, X., Wang, H., et al. 2013, ApJ, 779, 160, doi: 10.1088/0004-637X/779/2/160