# SEARCHING FOR SUPERNOVA NEUTRINOS WITH GW MEMORY TRIGGERS

MAINAK MUKHOPADHYAY  Department of Physics, Arizona State University,
Tempe, AZ 85287-1504, USA.

###### Abstract

Anisotropic neutrino emission from a core-collapse supernova (CCSN) causes a
permanent change in the local space-time metric, called the gravitational wave
(GW) memory. Long considered unobservable, this effect will be detectable in
the near future, at deci-Hertz GW interferometers. I will present a novel
idea, where observations of the neutrino GW memory from CCSNe will enable
time-triggered searches of supernova neutrinos at megaton (Mt) scale
detectors, which will open a new avenue to studying supernova neutrinos.

##  1 Introduction

In the current era of multi-messenger astronomy, GW, neutrino, photon, and
cosmic ray observations are combined to extract information about
astrophysical sources and phenomena in the universe [1]. Core-collapse
supernovae (CCSNe) have been one of the prime focuses of multi-messenger
astronomy [2]. Amongst the messengers from CCSNe, neutrinos are the most
dominant ones. Around, $3\times 10^{53}$ ergs ($\sim 99\%$) of the
gravitational binding energy of a star is emitted in neutrinos over a time
scale of $\sim 10$ s. But, the majority of CCSNe have been observed through
optical telescopes and we have observed supernova neutrinos from only one
CCSN, SN1987A [3]. This makes the aspect of detecting supernova neutrinos one
of profound interests and it leads to an obvious question - _why is it that we
have supernova neutrino observations from just one CCSN?_

(a) (a)

(b) (b)

Figure 1: (a) The total number of CCSNe versus the total number of neutrinos
obtained per megaton from those supernovae in 30 years (see text for details).
(b) Schematic diagram to illustrate the idea of memory-triggered supernova
neutrino detection.

The answer to this can be seen in Fig. 1a. We have plotted the total number of
core-collapse supernovae in 30 years versus the total number of neutrinos that
could be obtained per megaton from the supernovae in 30 years. This is done
for mean supernova neutrino energy of $11$ and $15$ MeV. The distance from
earth is marked by markers in both cases. We assume successful supernovae,
that is, neutron-star forming collapses (NSFCs) only. The details of the
parameters and the formulae required to obtain the plot will be discussed
later. But the evident fact is for nearby distances which would lead to very
high statistics, the rate of core-collapse supernova is extremely low. For
large distances, the number of core-collapse supernovae is very large but
because the neutrino flux falls off as $1/D^{2}$, the statistics is extremely
low.

The current situation on observing CCSN neutrinos is - with the upcoming
neutrino detectors, we will manage to have neutrinos observed from nearby
CCSNe ($D\leq 3$ Mpc) with significant statistics [4]. The diffuse supernova
neutrino background (DSNB) neutrinos will provide information from CCSNe from
$\mathcal{O}(1)$ Gpc distances and beyond [5, 6]. This leaves neutrino
observations from CCSNe in the local universe ($D\sim 3-100$ Mpc) to be an
area requiring considerable focus. Ideas proposing multi-megaton neutrino
detectors [7], using GW detectors as triggers [8] exist, but are limited to a
few Mpc. In this work, we present a summary of a novel idea of using the
neutrino GW memory effect to perform time-triggered searches for supernova
neutrinos [9], which would in essence be sensitive to $\mathcal{O}(100)$ Mpc.

A schematic representation of the main idea is shown in Fig. 1b. The core-
collapse supernova leads to an observable GW memory signal in the GW detectors
say at different times, $t_{1},t_{2},t_{3},t_{4},\dots$. We know the time
duration of a supernova neutrino burst, $\Delta t\sim 10$ s. Now, one can go
back to the neutrino detector data and analyze the events within time-windows,
$t_{1}+\Delta t,t_{2}+\Delta t,t_{3}+\Delta t,t_{4}+\Delta t,\dots$. Since,
the background in each of the time-windows is very low, owing to the time
duration of just $\sim 10$ s; the probability of finding a signal event (in
this case a supernova neutrino) is high. This method of reducing the
background events by using a triggered time-window of a short duration allows
collecting a clean sample of signal events over time.

In Sec. 2 we briefly discuss the dynamics of a CCSN and the main neutrino
emission phases. Sec. 3 focuses on the neutrino GW memory effect from a CCSN.
In Sec. 4, we discuss the formalism to calculate the number of memory-
triggered neutrino events possible. We discuss our main results in Sec. 5 and
conclude in Sec. 6.

##  2 Core-collapse supernova neutrinos

This section provides a very short review on the dynamics and neutrino
emission from CCSN [10, 11]. A CCSN is the death of a massive star ($M\geq
8\mathrel{{M_{\odot}}}$). The massive stars go through advanced stages of
nuclear burning until iron starts forming in the core. As iron forms in the
core, it triggers a set of events which finally leads to a pressure loss
followed by the collapse of the stellar core. The inner parts of the core
reach nuclear densities ($\sim 10^{13}$g/cm3), at which matter has been shown
to be incompressible. The infalling matter from the outer layers of the core,
stops abruptly because of this and produces a forward moving shockwave. The
shockwave propagates outwards but it keeps accreting matter, leading it to
stop eventually for a fraction of a second. Following this, it is either re-
energized when neutrinos deposit their energy resulting in a successful
supernova explosion, leading to the formation of a neutron star (hence also
called neutron star forming collapses (NSFCs)) or, sometimes the shock just
dies down resulting in a failed supernova, leading to a black hole being
formed (hence also called black hole forming collapses (BHFCs)). The main
neutrino emission phases from a CCSN are:

  * •

Neutronization burst phase: The sharp peak in the neutrino luminosity right
after the core bounce is known as the neutronization burst. This phase is
mostly dominated by electron neutrinos ($\nu_{e}$) and lasts for $\sim 2$ ms
post-bounce.

  * •

Accretion phase: The accretion phase lasts from $\sim 0.1-1$ s post-bounce.
During this phase, infalling matter keeps accreting on the forward moving
shock. The neutrino emission is nearly thermal and of the order $\sim 10^{52}$
ergs. Radial instabilities develop in this phases because of turbulence,
hydrodynamics instabilities and standing accretion shock instability (SASI).
Ths leads the neutrino emission in this phase being anisotropic, as has been
shown by numerical simulations [12, 13].

  * •

Cooling phase: In case of a NSFC, the proto-neutron star formed, cools by
emitting thermal neutrinos of all flavors characterizing the cooling phase.
This phase lasts $\sim 10$ s. Not much is known about the anisotropy in
neutrino emission in this phase.

##  3 The neutrino gravitational wave memory effect

The gravitational wave memory effect is the _permanent_ distortion of the
local space-time metric caused due to the propagation of a GW. At the linear
level, it is caused by gravitationally unbound systems that have anisotropic
emission of energy in the form of mass or radiation. Unfortunately, although
theoretically this phenomenon has been long known, we are yet to observe it.
This is because, it is a sub-dominant effect thus requires a very powerful
source, some anisotropy111Isotropic emissions will lead to the quadrupole
moment to be zero, resulting in no GW emission., and detectors that would
probe the frequency regimes of interest.

CCSNe have been considered to be an ideal source for GW memory observations.
The energy emitted in neutrinos is $\sim 10^{53}$ ergs, as discussed in Sec.
1, there is some anisotropy ($\sim 10^{-3}-10^{-2}$) in the neutrino emission
and the frequency regimes of interest to us is the low frequency (deci-Hz)
scale222Since the neutrino emission time scale is $\mathcal{O}(10)$ s, the
frequency is $\sim 0.01$ Hz. which will be explored by a bunch of planned
next-generation GW detectors like DECIGO [14], BBO [15]. Recently, the
neutrino GW memory effect was explored [16], in the context of building a new
phenomenological description, highlighting its detectability, and physics
potential in the present context.

Figure 2:  The neutrino luminosity $L_{\nu}(t)$ and anisotropy parameter
$\alpha(t)$ for NSFC (right) and BHFC (left).

The neutrino GW memory strain is given by [17, 16],

$h(r,t)=\frac{2G}{rc^{4}}\int_{-\infty}^{t-r/c}dt^{\prime}L_{\nu}(t^{\prime})\alpha(t^{\prime})\,.$
(1)

where $c$ is the speed of light, $t$ is the time post bounce and $G$ is the
universal gravitational constant. $L_{\nu}(t)$ is the all-flavors neutrino
luminosity and $\alpha(t)$ is the anisotropy parameter. It is evident from Eq.
1, that the GW strain depends on time-integral of the product of the neutrino
luminosity and the anisotropy parameter. We consider two distinct scenarios
and model the neutrino luminosity and anisotropy accordingly (shown in Fig.
2). In particular, the NSFC is based on the Ac3G model and the BHFC is based
on the LAc3G model [16]. In both cases we consider a non-zero anisotropy only
in the accretion phase. For the NSFC, we model the neutrino luminosity in the
accretion phase using a decaying exponential. The anisotropy parameter
consists of three-Gaussian bumps. In case of the BHFC, the neutrino luminosity
in the accretion phase is modeled using a constant. The anisotropy parameter
is broader, that is, it lasts for a longer duration and is modeled using
three-Gaussians combined together.

##  4 Formalism

The formalism mainly relies on, the probability that the GW wave detector sees
the memory signal from a CCSN in some galaxy that is at a distance say $r$
from the earth, the number of supernova neutrinos that would be detectable
from the same CCSN in a neutrino detector and the number of such CCSN possible
at that distance, which will depend on the CCSN rate at that distance. The
probability that the GW detector detects the memory signal depends on the
signal-to-noise ratio (SNR) - which depends on the signal strength, distance
to the source, and characteristics (power spectral noise density) of the
detector. One can then define the probability of detection $P_{det}^{GW}$
given a fixed false alarm probability [18]. The probability of detection
depends on the SNR, which in turn as mentioned depends on the distance to the
source.

For the neutrino detectors, we restrict ourselves to water Cherenkov detectors
(like HyperK [19]). The main channel of detection considered is inverse-beta
decay (IBD)333IBD: $\bar{\nu}_{e}+p\rightarrow e^{+}+n$.. The time-integrated
$\bar{\nu}_{e}$ flux $\Phi(r,E_{\nu})$ is given by the well-known analytical
quasi-thermal spectra [20], where we vary the time-integrated average energy
for the electron antineutrino flavor between [21] $11-15$ MeV for NSFC, and
$15-20$ MeV for BHFC. The predicted number of events in the detector from a
CCSN at distance $r$ is:

$N(r)=\int_{E^{th}_{\nu}}^{E^{max}_{\nu}}N_{p}\eta\sigma(E_{\nu})\Phi(r,E_{\nu})\
dE_{\nu}\,,$ (2)

where $N_{p}$ is the number of target protons, $\eta$ is the detection
efficiency, and $\sigma(E_{\nu})$ is the IBD cross-section. We take an energy
interval $\left[E^{th}_{\nu},E^{max}_{\nu}\right]=[19.3,50]\ \text{MeV}$ to
avoid the spallation background at low energy and the atmospheric neutrino
background at high energy [19].

Below $11$ Mpc the rate of CCSNe is taken from the CCSNe rates in the
individual galaxies [2]. Beyond $11$ Mpc, we use a constant volumetric rate of
CCSN given by444We ignore the redshift $z$ for distances concerning us.,
$R_{SN}=1.5\times 10^{-4}\ \text{Mpc}^{-3}\text{yr}^{-1}$. The cumulative rate
of CCSN can be calculated from this analytically by performing a suitable
volume integral.

##  5 Memory-triggered supernova neutrino search: Results

The total number of memory-triggered neutrino events from all CCSN within a
distance $D$ is given by,

$N^{trig}_{\nu}(D)=\Delta
T\sum_{j,r_{j}<D}R_{j}N(r_{j})P_{det}^{GW}(r_{j})\,,$ (3)

where, $\Delta T$ is the detector running time, $R_{j}$ indicates the CCSN
rate in the galaxy $j$. The other elements in the above formula are discussed
in Sec. 4. We show the main results ($N^{trig}_{\nu}(D)$ as a function of $D$
from Eq. 3) in Fig. 3, where we have considered - a) a population composed of
just NSFCs, and b) a _mixed_ population consisting of $60\%$ NSFCs + $40\%$
BHFCs. The shadings indicate the mean neutrino energy being varied for both
the NSFCs and BHFCs. The summation in Eq. 3 becomes an integral, beyond $11$
Mpc when the analytical rates for CCSN need to used.

Figure 3:  Number of memory-triggered neutrino events from CCSN at distance
$r<D$, as a function of $D$, for a Mt water Cherenkov detector and $\Delta
T=30$ years. The triggered-background is also shown.

For the NSFC population (a), time triggers from DECIGO+ will result in
$N^{trig}_{\nu}\sim 10-30$. For Ult. DECIGO, $N^{trig}_{\nu}\sim 100-300$ is
expected. For the mixed population (b), the prospects increase substantially
for DECIGO+, where it can help detect $\sim 100$ neutrino events due to its
increased distance sensitivity for memory in the BHFC scenario. For Ult.
DECIGO the results do not change significantly from the previous case. Even
DECIGO in its currently planned state, would help trigger and detect $\sim 10$
supernova neutrinos. We also show the background from the triggered-searches
in Fig. 3 (shaded violet area). This background associated with a triggered-
search is, $N^{trig}_{bckg}(D)=N^{trig}_{SN}(D)\lambda\Delta t$, where
$\lambda\simeq 1313$ events/year is the background rate in the detector [19],
$N^{trig}_{SN}(D)=\Delta T\sum_{j,r_{j}<D}R_{j}P_{det}^{GW}(r_{j})$ is the
number of CCSN memory signals observed. It is evident right away, that the
triggered-search background is orders of magnitude lesser than the untriggered
background, thus highlighting the efficiency of triggered-searches.
Furthermore, we see that even the triggered-background dominates the signal
events beyond $\sim 350$ Mpc and hence we restrict our analysis to $4<D<350$
Mpc.

##  6 Conclusion and Discussions

The technique described in this work is a new multi-messenger approach, in
which, the observation of the neutrino GW memory signals enable time-triggered
searches for supernova neutrinos. It could be realized in a few decades with
the upcoming deci-Hz GW interferometers and Mt-scale neutrino detectors. This
will help us collect samples of supernova neutrinos from the local universe
($\sim 10-100$ Mpc) which when compared with galactic SN and DSNB neutrinos
may help us perform population averaged energy and luminosity studies,
understand the distribution of SN populations including NSFCs and BHFCs, and
assist in the joint analysis of CCSNe using EM signals, taking our prospects
of analyzing CCSNe through their various messengers a step forward.

## Acknowledgments

We thank Zidu Lin, Carlos Cardona and Cecilia Lunardini in collaboration with
whom this work was performed. We are grateful to Cecilia Lunardini for
comments and careful reading of the manuscript. MM acknowledges funding from
the National Science Foundation grant number PHY-2012195 and the Fermi
National Accelerator Laboratory (Fermilab) Award No. AWD00035045 during this
work.

## References

## References

  * [1] Kristi Engel, Tiffany Lewis, Marco Stein Muzio, and Tonia M. Venters.  Advancing the Landscape of Multimessenger Science in the Next Decade.  In 2022 Snowmass Summer Study, 3 2022.
  * [2] Ko Nakamura, Shunsaku Horiuchi, Masaomi Tanaka, Kazuhiro Hayama, Tomoya Takiwaki, and Kei Kotake.  Multimessenger signals of long-term core-collapse supernova simulations: synergetic observation strategies.  Mon. Not. Roy. Astron. Soc., 461(3):3296–3313, 2016.
  * [3] Adam Burrows and James M. Lattimer.  Neutrinos from SN 1987A.  Astrophys. J. Lett., 318:L63–L68, 1987.
  * [4] Shin’ichiro Ando, John F. Beacom, and Hasan Yüksel.  Detection of neutrinos from supernovae in nearby galaxies.  Phys. Rev. Lett., 95:171101, Oct 2005.
  * [5] John F. Beacom.  The Diffuse Supernova Neutrino Background.  Ann. Rev. Nucl. Part. Sci., 60:439–462, 2010.
  * [6] Cecilia Lunardini.  Diffuse Neutrino Flux from Supernovae, page 1637.  2017\.
  * [7] Matthew D. Kistler, Hasan Yuksel, Shin’ichiro Ando, John F. Beacom, and Yoichiro Suzuki.  Core-Collapse Astrophysics with a Five-Megaton Neutrino Detector.  Phys. Rev. D, 83:123008, 2011.
  * [8] G. Pagliaroli, F. Vissani, E. Coccia, and W. Fulgione.  Neutrinos from supernovae as a trigger for gravitational wave search.  Phys. Rev. Lett., 103:031102, Jul 2009.
  * [9] Mainak Mukhopadhyay, Zidu Lin, and Cecilia Lunardini.  Memory-triggered supernova neutrino detection.  10 2021.
  * [10] Carlo Giunti and Chung W. Kim.  Fundamentals of Neutrino Physics and Astrophysics.  2007\.
  * [11] Hans-Thomas Janka, K. Langanke, A. Marek, G. Martinez-Pinedo, and B. Mueller.  Theory of Core-Collapse Supernovae.  Phys. Rept., 442:38–74, 2007.
  * [12] Kei Kotake, Wakana Iwakami, Naofumi Ohnishi, and Shoichi Yamada.  Ray-tracing analysis of anisotropic neutrino radiation for estimating gravitational waves in core-collapse supernovae.  The Astrophysical Journal, 704(2):951, 2009.
  * [13] David Vartanyan and Adam Burrows.  Gravitational Waves from Neutrino Emission Asymmetries in Core-collapse Supernovae.  Astrophys. J., 901(2):108, 2020.
  * [14] Seiji Kawamura, Masaki Ando, Naoki Seto, and et.al. Sato, Shuichi.  Current status of space gravitational wave antenna DECIGO and B-DECIGO.  Progress of Theoretical and Experimental Physics, 2021(5):05A105, May 2021.
  * [15] Kent Yagi and Naoki Seto.  Detector configuration of decigo/bbo and identification of cosmological neutron-star binaries.  Phys. Rev. D, 83:044011, Feb 2011.
  * [16] Mainak Mukhopadhyay, Carlos Cardona, and Cecilia Lunardini.  The neutrino gravitational memory from a core collapse supernova: phenomenology and physics potential.  JCAP, 07:055, 2021.
  * [17] E. Mueller and H. T. Janka.  Gravitational radiation from convective instabilities in Type II supernova explosions.  AAP, 317:140–163, January 1997.
  * [18] Piotr Jaranowski and Andrzej Krolak.  Data analysis of gravitational wave signals from spinning neutron stars. 3. Detection statistics and computational requirements.  Phys. Rev. D, 61:062001, 2000.
  * [19] K. Abe et al.  Letter of Intent: The Hyper-Kamiokande Experiment — Detector Design and Physics Potential —.  9 2011.
  * [20] Mathias Th. Keil, Georg G. Raffelt, and Hans-Thomas Janka.  Monte Carlo study of supernova neutrino spectra formation.  Astrophys. J., 590:971–991, 2003.
  * [21] Daniel Kresse, Thomas Ertl, and Hans-Thomas Janka.  Stellar Collapse Diversity and the Diffuse Supernova Neutrino Background.  Astrophys. J., 909(2):169, 2021.


# Power spectrum models for H$\alpha$ emission line galaxies from the Nancy
Grace Roman Space Telescope

Kevin S. McCarthy1, Zhongxu Zhai2,3,4,5, and Yun Wang6
1 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA
91109, USA
2 Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai 200240, China
3 Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai 200240,
China
4 Waterloo Center for Astrophysics, University of Waterloo, Waterloo, ON N2L
3G1, Canada
5 Department of Physics and Astronomy, University of Waterloo, Waterloo, ON
N2L 3G1, Canada
6 IPAC, California Institute of Technology, Pasadena, CA 91125, USA
E-mail<EMAIL_ADDRESS>
###### Abstract

The High Latitude Spectroscopic Survey (HLSS) is the reference baseline
spectroscopic survey for NASA’s Nancy Grace Roman space telescope, measuring
redshifts of $\sim 10$M H$\alpha$ emission line galaxies over a $2000$ deg2
footprint at $z=1-2$. In this work, we use a realistic Roman galaxy mock
catalogue to explore optimal modeling of the measured power spectrum. We
consider two methods for modelling the redshift-space distortions (one with
the canonical Kaiser squashing term $\mathcal{M}_{A}$, and another with a
window function on $\beta$ that selects out the coherent radial infall
pairwise velocities $\mathcal{M}_{B}$), two models for the nonlinear impact of
baryons that smears the BAO signal (one with a fixed ratio between the
smearing scales in the perpendicular $k_{*}^{\perp}$ and parallel
$k_{*}^{\parallel}$ and another where these smearing scales are kept as a free
parameters, P${}_{dw}(k|k_{*})$ and
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$), and two analytical nonlinear
growth corrections (one employing the halo model $\mathcal{F}_{HM}$ and
another formulated from simulated galaxy clustering of a semi-analytical model
$\mathcal{F}_{SAM}$). We find that the best model is
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})*\mathcal{M}_{B}$, which leads
to unbiased measurements of cosmological parameters. We expect the tools that
we have developed to be useful in probing dark energy and testing gravity
using Roman in an accurate and robust manner.

###### keywords:

large-scale structure of Universe – galaxies: statistics – cosmological
parameters

††pagerange: Power spectrum models for H$\alpha$ emission line galaxies from
the Nancy Grace Roman Space Telescope–References

##  1 Introduction

Next generation wide-field galaxy redshift surveys, such as those from the
Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016),
ESA’s Euclid mission (Laureijs et al., 2011), and NASA’s Nancy Grace Roman
Space Telescope (Spergel et al., 2015), will measure the 3D distribution of
galaxies to high precision, collecting more than $10$ times the amount of
redshift data currently available. To garner as much cosmological information
from these surveys as possible, we need accurate modelling of the non-linear
galaxy clustering signal. In this paper, we measure the galaxy clustering
signal in Fourier-spac, using a realistic Roman galaxy mock catalog (Zhai et
al., 2021a)111https://www.ipac.caltech.edu/doi/irsa/10.26131/IRSA546 between
$1.0\leq z\leq 2.0$, and use different methods to modify the linear
predictions of CDM for non-linear evolution of the dark matter halos and
baryonic effects so as to recover the input cosmology.

Galaxy redshift surveys such as the Sloan Digital Sky Survey (SDSS; Eisenstein
et al. 2005; Anderson et al. 2012, 2014; Ross et al. 2015; Beutler et al.
2017), 2dF Galaxy Redshift Survey (2dFGRS; Percival et al. 2001; Cole et al.
2005), WiggleZ (Blake et al., 2011; Kazin et al., 2014), 6dF Galaxy Survey
(6dFGS; Beutler et al. 2011, and Dark Energy Survey (DES; Abbott et al. 2022)
have enabled measurements of the baryon acoustic oscillation (BAO) signal at
various redshifts through the observation of galaxy clustering. Additional
cosmological information can be gained through the measurement of the growth
rate of large-scale structure ($f\sigma_{8}$; Blake et al. 2013; Reid et al.
2014; Alam et al. 2017; Zhai et al. 2019; Lange et al. 2022). These two
statistics (BAO and $f\sigma_{8}$) provide invaluable information about the
evolution of the universe and the nature of cosmic acceleration (Guzzo et al.,
2008; Wang, 2008b), and have so far been consistent with the Planck flat
$\Lambda$CDM model (Planck Collaboration et al., 2016). However, there are
tensions in the results from current data (Reid et al., 2014; Zhai et al.,
2022; Lange et al., 2022; Yuan et al., 2022). Therefore, it is important to
extend the current data set, as well take a closer look at how we model the
non-linear clustering signal.

Galaxy clustering in the quasi-linear regime can be modelled analytically with
linear perturbation theory (Goldstein et al., 2022) or with an effective field
theory (Perko et al., 2016). Although improvements in these analytical methods
are being made such that they can be applicable across smaller scales, the
usual practice for modeling galaxy clustering down into the non-linear regime
is to perform a cosmological N-body simulation, extract a halo catalog, and
assign galaxies to the host halos with an occupation model (Halo Model; Cooray
& Sheth 2002). The challenge with the halo model is that you are limited to
the cosmology of the N-body simulation. If you wish to have cosmological
parameters free you must perform many simulations and employ some technique to
interpolate between the cosmologies, i.e., use an emulator to sample the
response field (see Zhai et al. 2019). Additionally, how you establish the
dark matter halo-galaxy connection can impact cosmological constraints
(Wechsler & Tinker, 2018).

We explore corrections to linear cold dark matter (CDM) theory, as predicated
using the forward modelling code package camb222CAMB.info (Lewis et al., 2000;
Lewis & Bridle, 2002; Howlett et al., 2012), that account for non-linear
baryonic and dark matter halo growth effects. For the smearing out of the BAO
peak due to baryonic effects, we explore two methods to ’de-wiggling’ the BAO
signal, one with a single transition parameter $k_{*}$ (Wang et al., 2013),
that fixes the ratio between the parallel and perpendicular transition scale
at $k^{\parallel}_{*}=[1+f_{g}^{*}(z)]^{-1}k^{\perp}_{*}$ where
$k^{\perp}_{*}=k_{*}$, and another where the parallel and perpendicular
transition parameters are left as independent free variables (Beutler et al.,
2017), $\Sigma_{\perp}=G(z)/k^{\perp}_{*}$ and
$\Sigma_{\parallel}=G(z)/k^{\parallel}_{*}$, respectively. For non-linear
structure growth, we explore two analytical methods, one which incorporates
the halo model (HaloModel or HM; Mead et al. 2015) and one with a fitting
formula to a galaxy semi-analytical model (SAM) (Cole et al. 2005), including
an analysis on incorporating an additional term for higher $k$’s (Sánchez et
al., 2008).

In this work, we use a realistic galaxy mock for the High Latitude
Spectroscopic Survey (HLSS), the reference baseline survey for NASA’s Nancy
Grace Roman space telescope. It covers 2000 deg2 over the wavelength range of
1-1.93$\mu$m (H$\alpha$ redshift range $\sim$1-2), with a depth corresponding
to the H$\alpha$ line flux of $10^{-16}$ergs/s/cm2 (6.5$\sigma$), see Wang et
al. (2022). Our results should be qualitatively applicable to other galaxy
redshift surveys as well.

This paper is organized as follows. Sec.2 presents the methodology. Sec.3
presents results and discussions. We summarize and conclude in Sec.4.

##  2 Methodology

To determine the applicability of the various non-linear corrections to linear
theory mentioned in the Introduction, we will use a realistic Roman galaxy
mock to measure the galaxy power spectrum monopole and quadrupole moments,
$P_{0}(k)$ and $P_{2}(k)$, respectively. This measurement will be performed
following the standard FKP methodology (Feldman et al., 1994) with jackknife
covariance estimations. This will require us to model, in addition to the
physical non-linear effects, the geometry of the survey with window functions
and the redshift-space distortions (RSD). We will then employ a Monte-Carlo
Markov Chain (MCMC) statistical analysis to perform a recovery test of the
input cosmological parameters.

###  2.1 Mock lightcone galaxy catalog

We use the $2000$ deg2 lightcone mock (Zhai et al., 2021a) constructed using
the galacticus (Benson, 2012) semi-analytical model (SAM) and the dark matter
unit (Chuang et al., 2019) simulation, with cosmological parameters
$\Theta=[h,\Omega_{b},\Omega_{m},\Omega_{c},\sigma_{8},n_{s},A_{s}]=[0.6774,0.0462,0.3089,0.2627,0.8147,0.9667,2.06\times
10^{-9}]$, where $A_{s}$ is the amplitude of the primordial matter power
spectrum, which correspond to Planck 2015 cosmology model (Planck
Collaboration et al., 2016). The emission line luminosity of the galaxies is
computed using the cloudy photoionization code (Ferland et al., 2013). A full
description of the technique can be found in Merson et al. (2018).

This galaxy mock simulates the galaxy redshift catalog expected from the Roman
High Latitude Spectroscopic Survey (HLSS), where H$\alpha$ galaxies will be
observed over an area of $\sim 2000$ deg2 mapping out the 3D distribution of
$\sim 10$ million galaxies at $1.0<z<2.0$. We explore this distribution in
redshift slices $z=[(1.0,1.2),(1.2,1.4),(1.4,1.6),(1.6,2.0)]$. Two different
dust models where used to calibrate the SAM to replicate either the H$\alpha$
luminosity function observed in the ground-based narrow-band High-z Emission
Line Survey (HiZELS; Sobral et al. 2012) or the H$\alpha$ number counts
observed in the Hubble Space Telescope (HST) Wide-Field Camera 3 (WFC3)
Infrared Spectroscopic Parallel Survey (WISPS; Mehta et al. 2015). In this
work, we focus on the dust model calibrated to HiZELS. Note that the
observational systematic effects have not been added to the galaxy mock, since
they are not yet quantitatively modeled.

###  2.2 Power Spectrum Measurement

The power spectrum $P(\mathbf{k})$ measures the power of fluctuations in the
matter field with wavelengths $\lambda$ denoted by the wavenumber
$k=2\pi/\lambda$; it is the Fourier transform of the two point correlation
function $\xi(\mathbf{r})$. If we define the cosmological overdensity field as
$\delta(\mathbf{x})\equiv\rho(\mathbf{x})/\bar{\rho}-1$, where $\rho$ is the
matter density and $\bar{\rho}$ is its mean value, then the matter power
spectrum is

$\left<\delta_{k}\delta_{k^{\prime}}\right>\equiv(2\pi)^{3}P(\mathbf{k})\delta^{D}(\mathbf{k}+\mathbf{k}^{{}^{\prime}}),$
(1)

where $\delta^{D}$ is the Dirac delta function and
$\delta_{k}\equiv\int\mathrm{d}^{3}\mathbf{x}\delta(\mathbf{x})\mathrm{exp}(-i\mathbf{k}\cdot\mathbf{x})$
is the Fourier transform of the overdensity $\delta(\mathbf{x})$.

We measure $P(\mathbf{k})$ with an estimation of the density field and a
weighting scheme for the tracers following the FKP approach (Feldman et al.,
1994). This method is explained in more detail in $\S$3.1 $\&$ 3.2 of Gil-
Marín et al. (2015). In short, it amounts to performing a Fourier transform of
the galaxy fluctuation field, which is defined as

$F_{i}(\mathbf{r})\equiv
w_{FKP}(\mathbf{r})\lambda_{i}\left[w_{c}(\mathbf{r})n(\mathbf{r})-\alpha
n_{s}(\mathbf{r})\right],$ (2)

where $n_{s}$ is the number density of a random catalogue,
$\alpha\equiv\sum^{N_{gal}}_{i}w_{c}/N_{s}$, where $N_{s}$ is the number of
randoms and $N_{gal}$ is the number of galaxies, and $\lambda_{i}$ is a
normalization to produce unbiased estimators (for power spectrum $i=2$ and for
bispectrum $i=3$). $w_{FKP}$ and $w_{c}$ are weights encoding information on
observational systematic effects.

For the power spectrum, $i=2$ and when we perform a Fourier transform on
$F_{2}$, and convolve the field at $\mathbf{k}$, we have an estimator to
measure the observation, which relates to the theoretical model of
$P(\mathbf{k}$) in the following way:

$\left<|F_{2}(\mathbf{k})|^{2}\right>=\int\frac{\mathrm{d}^{3}\mathbf{k}^{\prime}}{(2\pi)^{3}}P_{gal}(\mathbf{k}^{\prime})|W_{2}(\mathbf{k}-\mathbf{k}^{\prime})|^{2}+P_{noise},$
(3)

where the left-hand side (LHS) of the equation represents what is measured
from observations and the right-hand side (RHS) is how we will model the
observation. Here, $P_{gal}$ is the theoretical prediction for the galaxy
power spectrum, $W_{2}$ is the window function that accounts for the finite
volume of the survey (to be discussed in more detail in $\S$2.6), and
$P_{noise}$ is the shot noise term that results from the discreetness of the
measurement.

We employ the nbodykit333nbodykit.readthedocs.io python package, which follows
the FKP weighting scheme, to measure the LHS of Eq. (3). In fact, nbodykit
models $P_{noise}$ as well, so we can go ahead and subtract that from the RHS.
The outputs from this measurement will, therefore, have the shot noise removed
but will require our theoretical model to be convolved with a window function.

Figure 1: Power spectrum multipole measurement scaled by wavenumber ($kP_{l}$)
of the Roman HLSS $2000$ deg2 lightcone mock in 4 redshift slices,
$(1.0,1.2),(1.2,1.4),(1.4,1.6),(1.6,2.0)$, denoted by square, triangle,
diamond, and cross, respectively. The monopole moments ($P_{0}$) are shown in
solid lines (with shot noise removed) and the quadrupole moments ($P_{2}$) are
shown in dashed lines.

For this analysis, we model the measurement of the power spectrum monopole
moment for the range $0.02\leq k\leq 0.25$ $h\,$Mpc-1 and the quadrupole
moment for the range $0.04\leq k\leq 0.25$ $h\,$Mpc-1, per Gil-Marín et al.
(2016), with bin-width $\Delta k=0.01$$h\,$Mpc-1 for galaxies with an
H$\alpha$ flux higher than $10^{-16}$ ergs s-1 cm-2. We perform these
measurements for the redshift slices
$z=[(1.0,1.2),(1.2,1.4),(1.4,1.6),(1.6,2.0)]$. These measurements are shown in
Fig. 1 where the monopole already has the shot noise removed.

####  2.2.1 Covariance Matrix

In order to perform the model likelihood analysis, we need a covariance matrix
for the observations. Typically, this is computed using a suite of mock
catalogs with the same galaxy clustering statistics as the observational data
set. For instance, Zhai et al. (2021a) used EZmocks (Chuang et al., 2014) to
construct their covariance matrices.

For convenience, we have chosen to construct our covariance matrix $C$ from
jackknife samples of the lightcone mock for each redshift slice. For the data
vector $\mathbf{P}=\\{P_{0},P_{2}\\}$, this is as follows:

$C(k_{1},k_{2})=\frac{n-1}{n}\sum_{i=1}^{n}(\mathbf{P}(k_{1})_{i}-\bar{\mathbf{P}}(k_{1}))(\mathbf{P}(k_{2})_{i}-\bar{\mathbf{P}}(k_{2})),$
(4)

where $n$ is the number of jackknife subsamples, where some portion of the
survey is removed and the statistics recomputed, and $\bar{\mathbf{P}}(k)$ is
the average of $\mathbf{P}(k)$ over all the subsamples. We chose to divide the
survey footprint into $n=400$ subsamples.

The problem with the jackknife method is that we are subsampling the single
observation that has been made and there is a limit to how often this can be
done. Therefore, our approximate covariance matrix will be noisy and how noisy
depends on how many subsamples we produce. When the inverse of this noisy
matrix is performed, to find the precision matrix $C^{-1}$, our results will
be biased. We can correct for the bias in our inversion with either the
Hartlap correction (Hartlap, J. et al., 2007) or through a Gaussian smoothing
technique (Mandelbaum et al., 2013), where $C$ is convolved with a 2D Gaussian
kernel with width $\sigma$ to remove the noise prior to inversion. We choose
to use the Gaussian smoothing technique with a smoothing kernel width
$\sigma=1.00$. We find, as have others (Mandelbaum et al., 2013; Lange et al.,
2022; Storey-Fisher et al., 2022), that this choice gives consistent result as
the Hartlap correction.

###  2.3 Power Spectrum Modeling

As the universe evolves, prior to recombination and the last scattering of the
cosmic microwave background (CMB), there are acoustic oscillations in the
baryon-photon fluid that couple gravitationally to the dark matter, which
itself is undergoing a mixing of modes across small-scales. These baryon
acoustic oscillations (BAO) become frozen post recombination, while the power
spectrum amplitude grows with redshift $z$ according to the linear growth
factor, $G(z)$. The transfer of power across frequencies, from the primordial
$P(k)$ to the beginning of linear growth, is captured by the linear transfer
function, $T_{L}(k)$, resulting in the linear matter power spectrum,

$P_{L}(k,z|P_{o})=G(z)P_{o}k^{n_{s}}T^{2}_{L}(k).$ (5)

The behavior of the linear matter power spectrum is that only the amplitude of
this signal will change as the universe evolves, the amplitude being a
function of $z$ such that $G(0)=1$ and $P_{o}$ is the amplitude in the present
universe, when $z=0$.

For $T_{L}$ of matter in a universe with baryons, where there are wiggles in
the power spectrum that produce the BAO bump in configuration-space, we
utilize the camb code package to forward model $T_{L}$ from cosmological
parameters, i.e.,
$\Theta=[h,\Omega_{b},\Omega_{m},\Omega_{c},\sigma_{8},n_{s},A_{s}]$. The
shape of the linear transfer function $T_{L}$ will be mostly determined by the
fraction of baryons to matter $\Omega_{b}/\Omega_{m}$. In the extreme case
where there are no baryons ($\Omega_{b}=0.0$), no wiggles appear in $T_{L}$.

####  2.3.1 Nonlinear smearing of BAO by baryons

Eisenstein & Hu (1998) found a simple fitting formula to this ’non-wiggled’
transfer function $T_{nw}$ as

$\begin{split}T_{nw}(q)=\frac{L}{L+Cq^{2}},\\\ L(q)=\mathrm{ln}(2e+1.8q),\\\
C(q)=14.2+\frac{731}{1+62.5q},\\\
q=k(\mathrm{T}_{CMB}/2.7)^{2}/\Omega_{m},\end{split}$ (6)

where $\mathrm{T}_{CMB}=2.726$ is the CMB temperature at $z=0$. Non-linear
baryonic effects have a tendency to broaden the BAO bump in configuration
space, which is equivalent to smearing out small-scale wiggles from the BAO
signal in $P_{L}(k)$. We account for this behavior following the method laid
out in Wang et al. (2013), where a ’de-wiggled’ linear transfer function is
given as

$\begin{split}T_{dw}^{2}(\mathbf{k},z|k_{*})=T_{L}^{2}(k)e^{-g_{\mu}k^{2}/(2k^{2}_{*})}\\\
+T_{nw}^{2}(k)[1-e^{-g_{\mu}k^{2}/(2k^{2}_{*})}],\end{split}$ (7)

with $g_{\mu}$ as

$g_{\mu}=G(z)^{2}[1-\mu^{2}+\mu^{2}(1+f_{g}^{*}(z))^{2}],$ (8)

where $f_{g}^{*}$ is related to the linear growth rate and
$\mu=\mathrm{cos}(\theta)$, $\theta$ being the angle of $k$ relative to the
line-of-sight (LoS). The $e^{-g_{\mu}k^{2}/(2k_{*}^{2})}$ term is derived in
Eisenstein et al. (2007), and captures the effects of nonlinear structure
formation on the signature of acoustic oscillations in the late-time galaxy
distribution (Angulo et al., 2008). The inclusion of the $\mu$-dependence in
the $g_{\mu}$ term, given in Eq. (8), accounts for the additional damping
along the LoS due to redshift-space distortions (RSD).

We can also think of this ‘de-wiggling’ as the addition of a BAO transfer
function to the zero-baryon transfer function, where the BAO damps away at
$k>k_{*}$. We see more clearly if we define
$T_{BAO}^{2}(k)=T_{L}^{2}(k)-T_{nw}^{2}(k)$, which makes Eq. (7) as follows:

$\begin{split}T_{dw}^{2}(\mathbf{k},z|k_{*})=T_{nw}^{2}(k)+T_{BAO}^{2}(k)e^{-g_{\mu}k^{2}/(2k^{2}_{*})}.\end{split}$
(9)

Figure 2: Upper panel: De-wiggled linear transfer function with single
parameter $k_{*}$ perpendicular and parallel to the line-of-sight, $\mu=0$ and
$\mu=1$, respectively. Lower Panel: De-wiggled linear transfer function with
two parameters $(\Sigma_{\perp},\Sigma_{\parallel})$ perpendicular and
parallel to the line-of-sight, $\mu=0$ and $\mu=1$, respectively. The scaling
parameters in each dimension are given in text in their respective panels. For
$T_{dw}(k|k_{*})$, $k_{*}^{\perp}=k_{*}$ and
$k_{*}^{\parallel}=k_{*}^{\perp}/(1+f_{g}^{*}(z))$. For
$T_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$,
$k_{*}^{\perp}=G(z)/\Sigma_{\perp}$ and
$k_{*}^{\parallel}=G(z)/\Sigma_{\parallel}$. These particular choices of
scaling parameters are from the MCMC fit to the observation between
$1.2<z<1.4$ with the lowest $\chi^{2}$.

Although $f_{g}^{*}$ is allowed to be a free parameter, it cannot be measured
from actual data due to parameter degeneracies, so it’s usually fixed, thus
effectively fixing the ratio of damping scales in the perpendicular and
parallel $k$-dimensions. To mitigate this, we could go back to the more basic
equation:

$\begin{split}T_{dw}^{2}(\mathbf{k}|\Sigma_{\perp},\Sigma_{\parallel})=T_{nw}^{2}(k)+T_{BAO}^{2}(k)e^{-[(1-\mu^{2})\Sigma_{\perp}^{2}+\mu^{2}\Sigma_{\parallel}^{2}]k^{2}/2},\end{split}$
(10)

where $G(z)/k^{\perp}_{*}$ is replaced by $\Sigma_{\perp}$ and
$G(z)/k^{\parallel}_{*}$ by $\Sigma_{\parallel}$.

The ‘de-wiggled’ power spectrum is then

$P_{dw}(\mathbf{k},z|P_{o},\Theta)=G(z)P_{o}k^{n_{s}}T^{2}_{dw}(\mathbf{k},z|\Theta),$
(11)

where $\Theta=k_{*}$ when there is a fixed ratio between parallel and
perpendicular dimensions, or $\Theta=\Sigma_{\perp},\Sigma_{\parallel}$ when
each dimension is free. Note that the $k_{*}$ technique effectively fixes
$\Sigma_{\parallel}/\Sigma_{\perp}=1+f_{g}^{*}(z)$. Also, for the
[$\Sigma_{\perp},\Sigma_{\parallel}$] technique, $T_{dw}$ does not explicitly
depend on $z$.

In Fig. 2, we show the behavior of the de-wiggled linear transfer function,
which is to transition from the linear transfer function predicted by CDM
($T_{L}$) that has wiggles present to the non-wiggled linear transfer function
($T_{nw}$) given in Eq. 6 that occurs when there are no baryons in the
universe. The effect is that at some scale denoted by $k_{*}$ or
$[\Sigma_{\perp},\Sigma_{\parallel}]$, $T_{dw}$ transitions from a lower power
with wiggles to a higher power without wiggles, which smears the BAO signal
and restores power for small-scale clustering.

In the upper panel, we see an example of the $T_{dw}(k|k_{*})$ model with
$k_{*}=0.248\,h$/Mpc, and in the lower panel, we see an example of the
$T_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ model with
$\Sigma_{\perp}=2.078\,h^{-1}$Mpc and $\Sigma_{\parallel}=5.909\,h^{-1}$Mpc,
perpendicular and parallel to the line-of-sight when $\mu=0$ and $\mu=1$,
respectively. The main difference to note here is that the transfer functions
when $\mu=0$ and when $\mu=1$ are fixed relative to each other for
$k_{*}$-parameterization and are not for
$[\Sigma_{\perp},\Sigma_{\parallel}]$-parameterization, which allows
$[\Sigma_{\perp},\Sigma_{\parallel}]$ to be more flexible. We will discuss
this more in the Results ($\S$3). The scaling parameters in these dimensions
are given in each panel in text, where $k_{*}^{\perp}=k_{*}$ and
$k_{*}^{\parallel}=k_{*}^{\perp}/(1+f_{g}^{*}(z))$ for $T_{dw}(k|k_{*})$ and
$k_{*}^{\perp}=G(z)/\Sigma_{\perp}$ and
$k_{*}^{\parallel}=G(z)/\Sigma_{\parallel}$ for
$T_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$. These values are from the MCMC
chain with the lowest $\chi^{2}$ in our fit to the data between $1.0<z<1.2$.

####  2.3.2 Nonlinear Matter Growth

The growth of the matter power spectrum post-recombination undergoes further
transfer across frequencies as gravity mixes the small-scale modes. We do not
detail the nonlinear growth according to CDM as a function of $z$ but rather
employ an analytical correction to the linear theory found when the linear
power spectrum was corrected to the non-linear clustering result of N-body
simulations: either the nonlinear matter power spectrum inferred from
clustering according to the halo model (HMcode or $\mathcal{F}_{HM}$; Mead et
al. 2015) or the nonlinear galaxy power spectrum found using a galaxy catalog
constructed with a semi-analytic model (SAM) ($\mathcal{F}_{SAM}$; Cole et al.
2005).

The nonlinear corrections of Mead et al. (2015) uses the halo model directly
in the modelling framework as opposed to the halofits method of (Smith et al.,
2003; Takahashi et al., 2012) that employs an empirical fitting algorithm to
N-body simulation clustering statistics but does not directly use the halo
model equations in their formula. Additionally, Mead et al. (2015) incorporate
parameters to account for the impact on the dark matter from baryonic
feedback. We performed a few minimization tests leaving the baryonic feedback
parameters free and found best fits similar to the default values as set in
camb. Since we also found little dependence on the particular choice of fixed
values we left them fixed at default. Joachimi et al. (2021) found results
similar to these defaults in their analysis of the KiDS-1000 joint galaxy
clustering and weak lensing analysis.

The nonlinear correction to $P_{L}$ according to HMcode for our chosen
cosmology is simply the ratio:

$\mathcal{F}_{HM}(k,z)=\frac{P_{NL}(k,z)}{P_{L}(k,z)}=\frac{T^{2}_{NL}(k,z)}{T^{2}_{L}(k)},$
(12)

where the nonlinear $P_{NL}$ and linear $P_{L}$ power spectrum has been fully
predicted to $z$ with the camb code, the ratio of which becomes the ratio of
the linear transfer functions squared when $G(z)P_{o}k^{n_{s}}$ drop out.
Multiplying $P_{L}$ by $\mathcal{F}_{HM}$ then gives us the nonlinear matter
power spectrum with baryonic feedback impacting the distribution of matter.
Joachimi et al. (2021) found HMcode to be comparable to the state-of-the-art
emulator CosmicEmu (Heitmann et al., 2014).

Figure 3: Upper Panel: Function form of $F_{HM}$ for the center of each
redshift slice observed in this work. Since $F_{HM}$ depends only on the input
cosmology and we keep this fixed, these are the exact functions used in each
model that incorporates $\mathcal{F}_{HM}$ in a particular redshift slice bin.
Lower Panel: Examples of $\mathcal{F}_{SAM}$ with the values of $Q_{1}$ and
$Q_{2}$ used in Eq. (13) given in the legend. These represent the best-fit
values found when modeling the data at each redshift slice with the
$P_{L}*\mathcal{F}_{SAM}$, i.e., blue solid line for best-fit to data between
$1.0<z<1.2$, dashed red line for $1.2<z<1.4$, dash-dot green line for
$1.4<z<1.6$, and dotted black line for $1.6<z<2.0$

Another method we use to incorporate nonlinear growth into our model is
similar to the halofits technique, in that a fitting formula is used to find
the correction to linear theory, but fits to a galaxy catalog rather than the
matter density field via halo clustering. Cole et al. (2005) found that when
comparing linear matter theory to the galaxy catalog output by a SAM (Benson
et al., 2000), where the halo merger trees have been used to paint galaxies
onto the halos, the following function can be used:

$\mathcal{F}_{SAM}(k|Q_{1},Q_{2}(,Q_{3}))=\frac{1+Q_{1}k^{2}}{1+Q_{2}k(+Q_{3}k^{2})},$
(13)

which when multiplied by $P_{L}$ produces the nonlinear galaxy power spectrum.
The $Q_{1}$ term has more impact on small-scales while $Q_{2}$ has more of an
impact on large-scales. The additional term in parenthesis $Q_{3}$ was
proposed by Sánchez et al. (2008) to improve modelling across larger
$k$-values. We will explore the need for $Q_{3}$ in this work.

Cole et al. (2005) compared using $\mathcal{F}_{SAM}$ to predictions of
halofits (see their paper $\S$7.2) for different cosmologies and found it to
be fairly robust with slight systematic deviations on small-scales from the
behaviors of galaxies. A similar result should be expected when comparing
$\mathcal{F}_{SAM}$ to $\mathcal{F}_{HM}$, since HMcode, like halofits, does
not model galaxy bias but rather relies on the linear galaxy bias $b_{g}$ of
our RSD models (see $\S$2.4), which is constant over all $k$’s, whereas
$\mathcal{F}_{SAM}$ corrects the linear matter power spectrum to model the
nonlinear galaxy power spectrum, thereby including nonlinear galaxy biasing in
its functional form.

Since our galaxy mock catalog was built using a SAM (see $\S$2.1), the naive
expectation is that $\mathcal{F}_{SAM}$ should be the correct model versus
using $\mathcal{F}_{HM}$ with $b_{g}$. That being said, we are observing only
galaxies above a particular $H\alpha$ flux cut and between redshifts
$1.0<z<2.0$, therefore it may be that the nonlinear galaxy clustering signal
is not as nuanced as at lower redshifts and either method is applicable. We
should note that $\mathcal{F}_{SAM}$ was not tested for robustness with
regards to changes in $z$.

In Fig. 3, we show examples of these nonlinear prefactors to the linear matter
power spectrum that result in a prediction of the nonlinear matter power
spectrum when using $\mathcal{F}_{HM}$ and the nonlinear galaxy power spectrum
when using $\mathcal{F}_{SAM}$. In the upper panel, we show the function
$\mathcal{F}_{HM}$ at the center of each redshift slice in our observation,
$z=1.1,1.3,1.5$ and $1.8$. Since this function comes from a ratio of the
nonlinear to linear matter power spectrum in camb, where the nonlinear matter
power spectrum has been given according to the HaloModel prescription which
only depends on cosmology, these will be the exact functions used for each
model that incorporates $\mathcal{F}_{HM}$ at each redshift, since we keep the
input cosmology fixed. This is different from the use of $\mathcal{F}_{SAM}$,
examples of which are shown in the lower panel, which does not depend on
cosmology directly but rather the $Q$-values used in Eq. (13), which are kept
as free variables in the modelling fit. Each of these 4 examples are actually
the best fit values found with the $P_{L}*\mathcal{F}_{sam}$ model at each
redshift slice, i.e., $Q_{1}=14.5$ and $Q_{2}=1.1$ when using
$P_{L}*\mathcal{F}_{sam}$ in the redshift slice $1.0<z<1.2$ while $Q_{1}=22.0$
and $Q_{2}=2.6$ for $1.6<z<2.0$. The inclusion of $Q_{3}$ will be discussed in
the Results ($\S\ref{sec:results}$).

####  2.3.3 Model Selection Choices

We have discussed $4$ different ways to model the nonlinear corrections to
linear theory of the matter power spectrum. However, these models are not
mutually exclusive and some combination of these $4$ could be more appropriate
than any single correction. That being said, if the galaxy bias is constant
and the linear approximation is correct, then $\mathcal{F}_{HM}$ and
$\mathcal{F}_{SAM}$ may be equivalent and we do not need to combine them.
Additionally, the two-parameter de-wiggle technique, with
$[\Sigma_{\perp},\Sigma_{\parallel}]$, can reduce to the one-parameter
technique, with $k_{*}$, if $\Sigma_{\parallel}/\Sigma_{\perp}=1+f_{g}^{*}(z)$
is preferred. The resulting interesting combinations are one of the de-wiggled
techniques with one of the non-linear growth functions, i.e.,
$P_{dw}(\mathbf{k}|k_{*})$ with either $\mathcal{F}_{HM}$ or
$\mathcal{F}_{SAM}$ and $P_{dw}(\mathbf{k}|\Sigma_{\perp},\Sigma_{\parallel})$
with either $\mathcal{F}_{HM}$ or $\mathcal{F}_{SAM}$. With each nonlinear
model separately and these combinations we then have $8$ different nonlinear
models (as listed in Table 1 and 2) to analyze.

###  2.4 Redshift Space Distortions

In real-space, where the universe is isotropic, there is no preferred LoS. The
clustering of galaxies is spherical. However, we observe galaxies from
redshift surveys in redshift-space and must account for RSDs, where there are
additional redshifts due to the peculiar motions of galaxies within the Hubble
flow (see Hamilton (1998) for review).

On larger scales, there is a squashing of the two-point correlation function
(2PCF) as galaxies fall into overdensities causing an additional red- or blue-
shift along the LoS. The squashing effect of the RSD is captured by the Kaiser
(Kaiser, 1987) factor $b^{2}(1+\beta\mu^{2})^{2}$ with the anisotropic
parameter $\beta=f_{g}/b$ detailing deviations from sphericity, where $f_{g}$
is the linear growth rate, and $b$ is the tracer bias. The linear growth rate
is the change in the linear growth factor $G$ with scale factor $a$, i.e.,

$f(a)=\frac{d\mathrm{ln}G}{d\mathrm{ln}a},$ (14)

which can be approximated as the cosmic matter density $\Omega_{m}(a)\equiv
8\pi G\rho_{m}(a)/(3H^{2}(z))$ raised to the growth index, $\gamma$ (Wang &
Steinhardt 1998), that is a prediction of the cosmological model ($\Lambda
CDM$ predicts $\gamma\simeq 0.55$, see Lue et al. 2004).

On smaller scales, galaxy peculiar velocities result in elongations of the
2PCF at small perpendicular separations. Since these elongations are always
along the LoS, it appears as if they are always pointing back to the observer
and are so denoted as the Finger-of-God (FoG) effect. The FoG is most often
analytically incorporated to the galaxy clustering by either an exponential
distribution or a Gaussian distribution, the Fourier transform of which
produces either a Lorentz damping or Gaussian damping term (Percival & White,
2009).

We model the RSD in two ways. The first has the standard Kaiser term for
large-scale squashing and an exponential distribution for the small-scale FoG,
which becomes a Lorenzian damping term in Fourier space. We refer to this
method for modeling the RSD as $\mathcal{M}_{A}$ and is:

$\mathcal{M}_{A}=\frac{b_{g}^{2}(1+\beta\mu^{2})^{2}}{1+\frac{1}{2}(k\mu\sigma_{r,v})^{2}}$
(15)

where $\sigma_{r,v}=\sigma_{v}/[H(z)a(z)]$, with $\sigma_{v}$ being the
pairwise velocity dispersion in [km/s] of the galaxies.

The second method for modeling the RSD more accurately model the Kaiser
squashing with the understanding that the peculiar velocities of the galaxy
populations are not strictly radial, or rather, the standard Kaiser term is
derived with a curl-free assumption about the pairwise velocity. However, the
measurement of the RSD will include some contribution from velocities that
have a curl. Zhang et al. (2013) decomposed the peculiar velocity into three
components: an irrotational component correlated(uncorrelated) with the
underlying density field, vδ(vS) and a rotational component vB. They found
that selecting the cosmological information from the irrotational correlated
velocity term vδ is equivalent to applying a window function ($\tilde{W}$) to
the $\beta$ parameter in the Kaiser term. Zheng et al. (2013) found that with
this windowed $\beta$ modification a Gaussian distribution for the FoG was
best able to recover the expected cosmology. Thus we have our second model for
the RSD

$\mathcal{M}_{B}=b_{g}^{2}\left(1+\beta\tilde{W(k,z)}\mu^{2}\right)^{2}\rm{exp}\left[-\frac{(k\mu\sigma_{r,v})^{2}}{2}\right]$
(16)

where

$\tilde{W}(k,z)=\frac{1}{1+\Delta\alpha(z)\Delta^{2}(k,z)}$ (17)

with $\Delta\alpha(z)$ being a free parameter to be determined by
observational data and the dimensionless power spectrum
$\Delta^{2}(k,z)=k^{3}P_{L}/(2\pi^{2})$.

###  2.5 Redshift-Space Multipoles

In redshift-space, where the clustering is no longer isotropic, additional
information can be gained by measuring the higher order multipoles of the
power spectrum

$P_{l}(k,z)=\frac{2l+1}{2}\int^{1}_{-1}P^{s}(\mathbf{k},z)L_{l}(\mu)d\mu,$
(18)

such that the sum of the multipoles produces the total power spectrum, i.e.,

$P^{s}(\mathbf{k},z)=\sum^{\infty}_{l=0,2,4,\cdots}P_{l}(k,z)L_{l}(\mu),$ (19)

where $L_{l}$ is the Legendre polynomial. For the linear power spectrum, the
three non-vanishing multipoles are the monopole ($l=0)$, quadrupole ($l=2$),
and hexadecapole ($l=4$). As discussed in $\S$2.2, we will model the monopole
and quadrupole moments, $P_{0}$ and $P_{2}$, respectively.

###  2.6 Survey window

Now with a theoretical model for the galaxy power spectrum, we simply need to
convolve $P(k)$ with a window function as required by Eq. (3), to account for
the finite survey volume, as follows:

$P^{c}(\mathbf{k})=\int\frac{d^{3}\mathbf{k}^{\prime}}{(2\pi)^{3}}P(\mathbf{k}^{\prime})|W_{2}(\mathbf{k}-\mathbf{k}^{\prime})|^{2}.$
(20)

Note this is the same function we saw on the RHS of Eq. (3). For our purpose,
it is efficient and sufficient to compute this by performing a convolution
with the power spectrum multipoles requiring only 1D FFTs. We will follow the
method detailed in Beutler et al. (2016) that converts the monopole and
quadrupole from Fourier into configuration-space, applies a survey mask, and
then performs an inverse transformation to the corrected correlation functions
to model the window corrected power spectrum multipoles. This technique uses
Hankel transformations to perform the FFTs, going from Fourier to
configuration-space with

$\xi_{l}(s)=\frac{4\pi(-i)^{l}}{(2\pi)^{3}}\int dkk^{2}P_{l}(k)j_{l}(sk),$
(21)

where $j_{l}$ is the spherical Bessel function of order $l$.

Specifically, we use the monopole and quadrupole with appropriate window
functions to find the corrected correlation functions by

$\displaystyle\xi_{0}^{c}$ $\displaystyle=$
$\displaystyle\xi_{0}W_{0}^{2}+\frac{1}{5}\xi_{2}W_{2}^{2}+\cdots$ (22)
$\displaystyle\xi_{2}^{c}$ $\displaystyle=$
$\displaystyle\xi_{0}W_{2}^{2}+\xi_{2}\left[W_{0}^{2}+\frac{2}{7}W_{2}^{2}+\frac{2}{7}W_{4}^{2}\right]+\cdots$
(23)

where the window function multipoles can be derived from the random pair
distribution as

$W_{l}^{2}(s)\propto RR(s,\mu)L_{l}(\mu)$ (24)

with the normalization that $W_{0}^{2}(s\rightarrow 0)=1$.

After applying these window functions to the two-point correlation statistics,
we perform an inverse transform to get back the window function corrected
power spectrum monopole and quadrupole moments, via

$P^{c}_{l}(k)=4\pi(-i)^{l}\int dkk^{2}\xi^{c}_{l}(s)j_{l}(sk).$ (25)

###  2.7 Scaling parameters

In constructing a realistic model for the power spectrum in redshift-space, we
have assumed a fiducial cosmological model (with parameters $\Theta$). In
actual measurements, if the true cosmology is different than the fiducial,
there will be a distortion in the 3D comoving coordinates known as the Alcock-
Paczynski, or AP, effect (Alcock & Paczynski, 1979). To account for this
effect, we introduce two scaling parameters, parallel and perpendicular to the
line-of-sight,

$\alpha_{\parallel}=\frac{H^{fid}(z)r_{s}^{fid}(z_{d})}{H(z)r_{s}(z_{d})},\hskip
14.22636pt\alpha_{\perp}=\frac{D_{A}(z)r_{s}^{fid}(z_{d})}{D_{A}^{fid}(z)r_{s}(z_{d})},$
(26)

where $H$ and $D_{A}$ are the Hubble parameter and angular diameter distance
for the fiducial model, and $r_{s}$ is the sound horizon at the drag epoch,
the superscript $fid$ indicating predictions from the fiducial model. These
scaling terms will distort the parallel and perpendicular components of the
$k$-vector: $k^{\prime}_{\parallel}=k_{\parallel}/\alpha_{\parallel}$ and
$k^{\prime}_{\perp}=k_{\perp}/\alpha_{\perp}$. This becomes a distortion in
the $k$ and $\mu$ via

$\displaystyle
k^{\prime}=\frac{k}{\alpha_{\perp}}\left[1+\mu^{2}\left(\frac{1}{F^{2}}-1\right)\right]^{1/2},$
(27)
$\displaystyle\mu^{\prime}=\frac{\mu}{F}\left[1+\mu^{2}\left(\frac{1}{F^{2}}-1\right)\right]^{-1/2},$
(28)

where $F=\alpha_{\parallel}/\alpha_{\perp}$. This modification will require
the multipole moments be multiplied by
$(\alpha_{\perp}^{2}\alpha_{\parallel})^{-1}$ such that Eq. (18) becomes

$P_{l}(k,z)=\frac{2l+1}{2\alpha_{\perp}^{2}\alpha_{\parallel}}\int^{1}_{-1}P^{s}(\mathbf{k^{\prime}},z)L_{l}(\mu^{\prime})d\mu.$
(29)

By incorporating this AP effect into our power spectrum model, we can very
simply test our cosmological model through the constraints on
$\alpha_{\parallel}$ and $\alpha_{\perp}$, the expectation being that if our
fiducial cosmology is correct we will find $\alpha_{\parallel,\perp}=1.0$.
These terms allow us to constrain cosmology from distortions in the shape of
the BAO signal, detailing how the parallel and perpendicular dimensions are
altered from an incorrect cosmological prior.

Similarly, we can use an alpha term to parameterize deviations from the
fiducial linear growth parameter (Wang et al., 2013),

$f_{g}(z)\sigma_{m}(z)\equiv f_{g}(z)G(z)\sqrt{P_{o}},$ (30)

relative to predictions of the fiducial cosmology as Zhai et al. (2021a)

$\displaystyle\alpha_{g}\equiv\frac{f_{g}(z)\sigma_{m}(z)}{f_{g,fid}(z)\sigma_{m,fid}(z)}=\frac{f_{g}(z)G(z)\sqrt{P_{o}}}{f_{g,fid}(z)G_{fid}(z)\sqrt{P_{o,fid}}}$
(31)
$\displaystyle=\frac{\beta(z)\sqrt{P_{n}(z)}}{f_{g,fid}(z)G_{fid}(z)\sqrt{P_{o,fid}}},$
(32)

where $P_{n}(z)=P_{o}G^{2}(z)b^{2}(z)$. Compared to the widely used
$f_{g}(z)\sigma_{8}(z)$ parametrization, $f_{g}(z)\sigma_{m}(z)$ has the
advantage of being independent of the Hubble constant.

###  2.8 Parameter Constraints

To test our nonlinear models, we will evaluate how well each model is able to
recover the input cosmology through analysis on the constraints of
$\alpha_{\perp}$, $\alpha_{\parallel}$, and $\alpha_{g}$. We expect
$\alpha_{\perp,\parallel,g}=1.0$ within $1\sigma$ if our modeling is
sufficiently accurate.

To find the posterior distribution of the modelling parameters, we will
perform a likelihood analysis with the Markov Chain Monte Carlo (MCMC)
technique, with $\chi^{2}$ as follows:

$\chi^{2}=\sum_{i,j}(\mathbf{P}_{obs,i}-\mathbf{P}_{th,i})C^{-1}_{ij}(\mathbf{P}_{obs,j}-\mathbf{P}_{th,j}),$
(33)

where the index $i$($j$) indicates the data vector at $k_{i}$($k_{j}$),
$\mathbf{P}_{obs}$ is from the measurement of the mock galaxy lightcone while
$\mathbf{P}_{th}$ is the prediction of the data vector from theory, monopole
and quadrupole in Eq. (25) respectively, and $C^{-1}$ is the inverse of the
covariance matrix found in Eq. (4).

All models have the set of parameters
$[\beta,P_{o},\sigma_{v},\alpha_{\perp},\alpha_{\parallel}]$, with priors
$\beta=[0.3,0.85]$, log${}_{10}(P_{o})=[4,8]$, $\sigma_{v}=[10,700]$,
$\alpha_{\perp}=[0.5,1.5]$, and $\alpha_{\parallel}=[0.5,1.5]$. Both FoG
models, $\mathcal{M}_{A}$ and $\mathcal{M}_{B}$ have the parameter
$\sigma_{v}$. If $\mathcal{M}_{B}$ is used, there is an additional parameter
$\Delta\alpha$ with priors $\Delta\alpha=[0,1.6]$. For the single parameter
de-wiggle power spectrum the parameter and bounds $k_{*}=[0.01,1.5]$. For the
two parameter de-wiggle model we have the parameters
$[\Sigma_{\perp},\Sigma_{\parallel}]$ both with bounds $[0.5,30.0]$. The
HaloModel nonlinear correction $\mathcal{F}_{HM}$ does not have any additional
parameters since, as discussed in $\S$2.3.2, we fix the baryon feedback
parameters to the default value in camb. For $\mathcal{F}_{SAM}$ there are
three parameters $[Q_{1},Q_{2},Q_{3}]$ with bounds
$[(0.01,15.0),(2.0,30.0),(0.2,10.0)]$.

Table 1:  Nonlinear $P(k)$ model reduced $\chi^{2}$ values for $k_{max}=0.25h$/Mpc. Bold values are the best for a given redshift while italics indicate second best for a particular $z-$slice. For reasons discussed in $\S$3, we only detail results for the nonlinear models using the $\mathcal{M}_{B}$ FoG model.  Model ($k_{max}=0.25h$/Mpc) | $1.0<z<1.2$ | $1.2<z<1.4$ | $1.4<z<1.6$ | $1.6<z<2.0$ | Sum${}_{\chi^{2}_{r}}$
---|---|---|---|---|---
$P_{dw}(k|k_{*})$ | 0.938 | 0.483 | 0.454 | 1.326 | 3.201
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ | 0.936 | 0.478 | 0.470 | 1.301 | 3.185
$P_{L}*F_{HM}$ | 1.105 | 0.619 | 0.596 | 1.978 | 4.298
$P_{L}*F_{SAM}$ | 1.245 | 0.684 | 0.607 | 1.393 | 3.929
$P_{dw}(k|k_{*})*F_{HM}$ | 1.072 | 0.556 | 0.519 | 1.310 | 3.457
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})*F_{HM}$ | 1.093 | 0.567 | 0.495 | 1.321 | 3.373
$P_{dw}(k|k_{*})*F_{SAM}$ | 1.013 | 0.492 | 0.527 | 1.341 | 3.409
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})*F_{SAM}$ | 1.031 | 0.477 | 0.527 | 1.386 | 3.421
Table 2:  Nonlinear $P(k)$ model Figure-of-Bias (FoB) for $(\alpha_{\perp},\alpha_{\parallel},\alpha_{g})$-parameter space for $k_{max}=0.25h$/Mpc. Bold values are the best for a given redshift while italics indicate second best for a particular $z-$slice. For reasons discussed in $\S$3, we only detail results for the nonlinear models using the $\mathcal{M}_{B}$ FoG model.  Model ($k_{max}=0.25h$/Mpc) | $1.0<z<1.2$ | $1.2<z<1.4$ | $1.4<z<1.6$ | $1.6<z<2.0$ | SumFoB
---|---|---|---|---|---
$P_{dw}(k|k_{*})$ | 1.97 | 0.35 | 0.66 | 1.02 | 4.00
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ | 1.15 | 0.25 | 1.14 | 1.48 | 4.02
$P_{L}*F_{HM}$ | 2.23 | 1.11 | 2.06 | 1.90 | 7.30
$P_{L}*F_{SAM}$ | 1.48 | 0.33 | 1.20 | 1.50 | 4.51
$P_{dw}(k|k_{*})*F_{HM}$ | 1.63 | 0.48 | 1.32 | 0.84 | 4.27
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})*F_{HM}$ | 1.05 | 1.05 | 1.06 | 1.55 | 4.71
$P_{dw}(k|k_{*})*F_{SAM}$ | 1.02 | 0.57 | 0.58 | 1.71 | 3.88
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})*F_{SAM}$ | 0.64 | 0.90 | 1.30 | 1.43 | 4.27

We could continue to add parameters and find a model that fits the clustering
data perfectly, but this does not mean we have found the best model according
to Occam’s Razor. There are many ways to penalize a model for having more
parameters, such as the Bayesian or Akaike Information Criterion and the
reduced chi-squared value, for instance. We choose to evaluate the reduced
chi-square, which is

$\chi^{2}_{r}=\chi^{2}/\nu,$ (34)

where the degrees-of-freedom $\nu=n-m-1$, with $n$ and $m$ being the number of
data points fit and the number of modelling parameters, respectively.

We should then select the model with the lowest $\chi^{2}_{r}$ value and call
that our best model. However, fitting the data more correctly does not mean we
have the most physically relevant model. For instance, a particular model
could chase the data points with smaller errors and miss the data with larger
errors, which in this case means fitting the small-scale clustering data but
missing the large-scale clustering where most of the constraining power on
cosmology resides. We will then use $\chi^{2}_{r}$ as a validation test of the
model, i.e., the model is valid if $\chi^{2}_{r}\sim 1$, but not use it to
determine which model is best as this quantity ignores correlations between
parameters, and different models can have different dependencies on the same
parameter.

Since we are more interested in how well a model is able to recover the input
cosmology through analysis of the posteriors on $\alpha_{\perp,\parallel,g}$,
any model that has $\alpha_{\perp,\parallel,g}$ within $1\sigma$ of unity is a
good model. To determine the best model(s) we will quantify this with a
Figure-of-Bias (FoB) value (Eggemeier et al., 2021), which is closely related
to how Figure-of-Merit (FoM) is defined (Wang, 2008a). The FoB for the three
alpha parameters is as follows:

$FoB=\left[\sum_{i,j}(\bar{\theta}_{i}-\theta_{fid,i})S^{-1}_{ij}(\bar{\theta}_{j}-\theta_{fid,j})\right]^{1/2},$
(35)

where
$\bar{\theta}=[\bar{\alpha}_{\perp},\bar{\alpha}_{\parallel},\bar{\alpha}_{g}]$
is the average of the posterior distributions of the $\alpha$-terms from the
MCMC chains, $\theta_{fid}=[1,1,1]$, and $S_{ij}$ is the covariance matrix of
the $[\alpha_{\perp},\alpha_{\parallel},\alpha_{g}]$ posterior distribution.
Generally speaking, smaller FoB indicates better accuracy in recovering the
input cosmological model. However, it is possible that two posterior
distributions score the same FoB if one has a mean closer to unity with a
smaller $\sigma$ while the other is further from unity with a larger $\sigma$.
We should then think of the FoB as an indication of how well the truth falls
within $1\sigma$ of the resulting posterior distribution.

Table 3:  Nonlinear P(k) model reduced $\chi^{2}$ values with $k_{max}=0.3h$/Mpc.  Model ($k_{max}=0.3h$/Mpc) | $1.0<z<1.2$ | $1.2<z<1.4$ | $1.4<z<1.6$ | $1.6<z<2.0$ | Sum${}_{\chi^{2}_{r}}$
---|---|---|---|---|---
$P_{dw}(k|k_{*})$ | 0.889 | 0.457 | 0.453 | 1.157 | 2.956
$P_{dw}(k|k_{*})*F_{SAM}$ | 0.982 | 0.460 | 0.532 | 1.241 | 3.215
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ | 0.874 | 0.457 | 0.465 | 1.128 | 2.924
Table 4:  Nonlinear P(k) model Figure-of-Bias (FoB) for $(\alpha_{\perp},\alpha_{\parallel},\alpha_{g})$-parameter space with $k_{max}=0.3h$/Mpc.  Model ($k_{max}=0.3h$/Mpc) | $1.0<z<1.2$ | $1.2<z<1.4$ | $1.4<z<1.6$ | $1.6<z<2.0$ | SumFoB
---|---|---|---|---|---
$P_{dw}(k|k_{*})$ | 1.87 | 0.09 | 1.30 | 0.94 | 4.20
$P_{dw}(k|k_{*})*F_{SAM}$ | 1.54 | 0.34 | 1.45 | 1.09 | 4.42
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ | 0.92 | 0.54 | 1.13 | 1.10 | 3.69

##  3 Results and Discussion

###  3.1 Comparison of different models

Evaluation of $\chi^{2}_{r}$ and FoB have shown that the FoG model
$\mathcal{M}_{B}$ is consistently better than $\mathcal{M}_{A}$. Zhai et al.
(2021a) also found that using $\mathcal{M}_{B}$, i.e., placing a window
function on $\beta$ to select the pairwise velocities proper for the Kaiser
formulation as suggested in Zhang et al. (2013), was generally more accurate.
Specifically, we found that $\mathcal{M}_{A}$ would consistently produce
constraints on $\alpha_{\parallel}$ lower than found when using
$\mathcal{M}_{B}$, indicating the bias induced by the choice of theoretical
template. For the $\mathcal{F}_{SAM}$ nonlinear correction, there is an
additional term $Q_{3}$ introduced by Sánchez et al. (2008) to produce better
fits at larger $k$’s. We found that including this parameter did not improve
our constraints, with $Q_{3}$ having a flat posterior distribution. Therefore,
when we talk about $\mathcal{F}_{SAM}$ we refer to the original formulation by
Cole et al. (2005) without the $Q_{3}$ term. We present the $\chi^{2}_{r}$ and
FoB for the $8$ non-linear modelling combinations all with $\mathcal{M}_{B}$
as the FoG model in Tables 1 and 2. These constraints are all produced from
measurements with $k_{max}=0.25\,h$/Mpc.

In Table 1, we detail the $\chi^{2}_{r}$ values for each model in each
redshift-slice with the lowest $\chi^{2}_{r}$ in bold and the second lowest in
italics. In the rightmost column we sum each redshift slice and again show the
first and second best-fits in bold and italics. When employing only a de-
wiggled model, using two damping parameters for the de-wiggled power spectrum,
notated as $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$, is better than
having a single parameter, notated as $P_{dw}(k|k_{*})$, with a fixed ratio
between the parallel and perpendicular directions, such that
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ is equivalent to
$P_{dw}(k|k_{*})$ if $\Sigma_{\parallel}=(1+f_{g}(z))*\Sigma_{\perp}$ and
$\Sigma_{\perp}=k_{*}^{-1}$. If we only consider the nonlinear correction
functions multiplied by the linear power spectrum ($P_{L}*\mathcal{F}_{HM}$ vs
P${}_{L}*\mathcal{F}_{SAM}$), both fit the clustering data similarly. We find
that incorporating the de-wiggled power spectrum, either one or two free
parameters, with a nonlinear correction function is always better than if we
use only a nonlinear correction function. That being said, these
$\chi^{2}_{r}$ values suggest that all models are ‘good fits’, or rather, the
differences are not significant enough to rule out any model, with all
$\chi^{2}_{r}$ being close to unity. This is likely due to the fact that at
$z>1$, matter density fluctuations are close to linear at
$k_{max}=0.25\,h$/Mpc.

Evaluating the FoB results in Table 2 we find
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$ performs the best. Second best recovery of
the input cosmology is accomplished by $P_{dw}(k|k_{*})$ while third is
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$. That being said, we are
interested in the best and most unique models. Considering that
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$ is an extension of $P_{dw}(k|k_{*})$, it
wouldn’t make much sense to consider P$(k|k_{*})$ as being the second best
model rather than a model that performs slightly worse than the best model, or
rather adding $\mathcal{F}_{SAM}$ to the $k_{*}$ de-wiggle power spectrum is
better than using the $k_{*}$ alone. If we consider the de-wiggled model with
$[\Sigma_{\perp},\Sigma_{\parallel}]$, we see that using this alone is better
than adding either $\mathcal{F}_{HM}$ or $\mathcal{F}_{SAM}$. This suggests
that the $k_{*}$ model is missing something in the behavior of the observed
clustering and, therefore, needs an additional modification, with
$\mathcal{F}_{SAM}$ being better than $\mathcal{F}_{HM}$, whereas
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ does not benefit from these
additions, suggesting that whatever $\mathcal{F}$ adds to the model is already
captured by the behavior of $[\Sigma_{\perp},\Sigma_{\parallel}]$.

###  3.2 Impact from $k_{max}$

When we extend the observation to $k_{max}=0.3\,h$/Mpc for these three best
models, see Table’s 3 and 4, we find that
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ now performs the best in both
$\chi^{2}_{r}$ and FoB while $P_{dw}(k|k_{*})$ is second best in both
statistics. We also re-evaluate the inclusion of $Q_{3}$ in
$P_{L}*\mathcal{F}_{SAM}$ and find that the posterior is no long flat for
$P_{L}*\mathcal{F}_{SAM}$, preferring the upper end of the bounds, but that
the improvement on the constraints is not significant and not as beneficial as
using $\mathcal{F}_{SAM}$ without $Q_{3}$ but with $k_{*}$. So if we are going
to add one additional parameter to $\mathcal{F}_{SAM}$ we would select $k_{*}$
over $Q_{3}$. If we include both, the posterior on $Q_{3}$ is again flat.

Considering both the fits with $k_{max}=0.25\,h$/Mpc and $k_{max}=0.3\,h$/Mpc,
it would seem that $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ and either
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$ (without $Q_{3}$) or $P_{dw}(k|k_{*})$ are
the best models for this set of simulated Roman HLSS data, all with
$\mathcal{M}_{B}$ as the FoG model. The level of uncertainties on
($\alpha_{\perp}$, $\alpha_{\parallel}$, $\alpha_{g}$) for these three models
are approximately at the ($2.45$, $3.66$, $8.28$), ($2.75$, $3.87$, $9.04$),
and ($2.67$, $4.24$, $8.28$) percent-level, for $P_{dw}(k|k_{*})$,
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$, and
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$, respectively, when
$k_{max}=0.25\,h$/Mpc and at the ($2.43$, $3.54$, $8.21$), ($2.92$, $3.54$,
$9.00$), and ($2.65$, $4.08$, $8.71$) percent-level when $k_{max}=0.3\,h$/Mpc.
These values are averages over all redshift slices.

Figure 4: The best-fit clustering results from the MCMC chains are shown for
the three best models, $P_{dw}(k|k_{*})$, $P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$,
and $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$, for each redshift slice:
$1.0<z<1.2$, $1.2<z<1.4$, $1.4<z<1.6$, $1.6<z<2.0$. All constraints are with
FoG model $\mathcal{M}_{B}$ and with $k_{max}=0.3$. Figure 5: Constraints on
$\alpha_{\perp}$, $\alpha_{\parallel}$, and $\alpha_{g}$ parameters are shown
for the three best models, $P_{dw}(k|k_{*})$,
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$, and
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$, for each redshift slice:
$1.0<z<1.2$, $1.2<z<1.4$, $1.4<z<1.6$, $1.6<z<2.0$. All constraints are with
FoG model $\mathcal{M}_{B}$ and with $k_{max}=0.3$.

We show the best-fit clustering in Fig. 4 and the $\alpha_{\perp,\parallel,g}$
constraints for these best models for the different redshift slices when
$k_{max}=0.3\,h^{-1}$Mpc in Fig. 5. From these figures we can see that there
is not much difference between these three best models. What is clear is that
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$ does not always perform the same as the
other models. For instance, in Fig. 4 $P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$
tends to oscillate more in the quadrupole moment, $P_{2}$, and in Fig. 5 it is
more degenerate in the first redshift slice than the others, upper left panel.
Comparing $P_{dw}(k|k_{*})$ to $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$,
both have essentially the same behavior when modelling the monopole and
quadrupole moments over all redshift slices. This is represented in the
$\alpha_{\perp,\parallel,g}$ constraints as well. However, the additional free
parameter in $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ does add some
additional degeneracy and inflates the contours.

###  3.3 Comparison with earlier works

One of the key goals of this paper is to evaluate the uncertainty of BAO and
RSD signals from Roman. The simulation and modeling is similar to Zhai et al.
2021a, but with more thorough investigation of the galaxy power spectrum
template, thus we are able to recover the input cosmological model more
accurately and consistently. When the same model is used, we find consistency
of the $\alpha$ parameters as expected. However, we should note one difference
from the estimate of the covariance matrix. Zhai et al. (2021a) adopt an
EZmock approach and produce thousands of approximate mocks. In this work, we
use the data set itself to estimate the covariance matrix with a jackknife
subsampling method.

The constraints found here are less tight than results from the Fisher matrix
forecast of Wang et al. (2022), indicating that we can further tighten the
constraints by improving the modeling approach. The lack of significant
improvement in constraints by going to smaller scales (from
$k_{max}=0.25\,h$/Mpc to $k_{max}=0.3\,h$/Mpc) supports this. One possibility
is to use a suite of tailored galaxy mocks, e.g., from BAM (Balaguera-
Antolínez et al. 2019), to compute the covariance matrix, as jackknife
approach is known to overestimate the measurement uncertainties. Although a
detailed comparison of different methods is beyond the scope of this paper, we
should be aware of the possible impact from the modeling given the high
statistical precision expected from Roman HLSS. Examples for the test of
covariance matrix along this direction can be found in Mohammad & Percival
(2022) and Percival et al. (2022).

###  3.4 Modeling BAO Damping

The fact that $\mathcal{F}_{SAM}$ always performs better that
$\mathcal{F}_{HM}$ was expected given that the nonlinear correction of the
HaloModel was constructed to produce the nonlinear power spectrum for matter
not for the galaxies. We, therefore, had to include a galaxy bias term $b_{g}$
in our $FoG$ models. For $\mathcal{F}_{SAM}$, we do not need this term since
the nonlinear correction for this technique was constructed to produce the
nonlinear power spectrum of the galaxies. The reason this performs better is
likely because the linear galaxy bias $b_{g}$ that we use is a constant value
across all scales and it is known that on small-scales the galaxy bias becomes
a function of separation. It can be assumed that $\mathcal{F}_{SAM}$ has
incorporated into its functional form a nonlinear galaxy bias.

Figure 6: Comparison between $\Sigma_{\perp}$ and $\Sigma_{\parallel}$ from
constraints for P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ model for
different redshift slices and the prediction for the relationship
$\Sigma_{\parallel}=(1+f(z))*\Sigma_{\perp}$, which is the assumption made in
the P${}_{dw}(k|k_{*})$ model, and also with the relationship
$\Sigma_{\parallel}=(b_{g}(z)+f(z))*\Sigma_{\perp}$, where $b_{g}(z)$ is
estimated with Eq. 36. Solid markers are for constraints with $k_{max}=0.25$
and the slightly offset empty markers are for fits with $k_{max}=0.3$.

We should then ask, why are P${}_{dw}(k|k_{*})*\mathcal{F}_{SAM}$ and
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ comparable models? To answer
this we shall first consider the physical differences between the de-wiggled
models, P${}_{dw}(k|k_{*})$ and
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$. Mathematically,
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ reduces to P${}_{dw}(k|k_{*})$
when $\Sigma_{\parallel}=[1+f_{g}(z)]\Sigma_{\perp}$. This particular
relationship is derived in Eisenstein et al. (2007), where they explore the
nonlinear degradation of the BAO signal and propose what we call the
P${}_{dw}(k|k_{*})$ model. The inspiration for this model comes from the
consideration that structure formation and bulk flows on $6\,h^{-1}$Mpc scales
will smear out the BAO signal at separations of $100\,h^{-1}$Mpc, which
broadens the BAO peak. Since it is the wiggles in Fourier-space that produce
the BAO peak in configuration-space, hence the reason for the term Baryon
Acoustic Oscillations, this smearing of the bump in configuration-space
equates to fewer wiggles in Fourier-space. This behavior is captured in the
first portion of the right-hand side of Eq. (7), where $T_{L}$ is modulated by
the smearing Gaussian term. The second part of the right-hand side of this
equation, where the non-wiggled $T_{nw}$ is included, works to restore power
to the small-scales from small-scale clustering where there are no wiggles.

Eisenstein et al. (2007) derived this model using the Zel’dovich approximation
of the Lagrangian displacement (Zel’dovich, 1970), where the displacement
along the line-of-sight in redshift-space is increased by a factor of $f_{g}$,
the linear growth rate given in Eq. (14), hence the relationship
$\Sigma_{\parallel}=(1+f_{g}(z))*\Sigma_{\perp}$. The reason for this term is
the same for the Kaiser squashing term $b_{g}(1+\beta\mu^{2})$, which reduces
to $(1+f_{g})$ when considering only the parallel direction where $\mu=1$ and
the RSD of matter overdensity when $b_{g}=1$ while $\beta=f/b_{g}$. Since we
are using galaxies as tracers of the matter field we should allow cases when
$b_{g}\neq 1$, and instead have the relationship
$\Sigma_{\parallel}=(b_{g}(z)+f_{g}(z))*\Sigma_{\perp}$.

In Fig. 6, the constraints on $\Sigma_{\perp}$ (squares) and
$\Sigma_{\parallel}$ (diamonds) are given from the
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ modelling of the simulated
observation at the different redshift slices when $k_{max}=0.25\,h$/Mpc (solid
markers) and when $k_{max}=0.3\,h\,$/Mpc (offset empty markers). In this plot
we also give the predication for $\Sigma_{\parallel}$ if the relationship
$\Sigma_{\parallel}=(1+f(z))*\Sigma_{\perp}$ is adopted, the assumption in the
$k_{*}$ model, and we see in the triangles that these values fall lower than
the constraints on $\Sigma_{\parallel}$ when left free in the MCMC. However,
these are only lower by about $1\sigma$, so may not indicate a total break
down of the $k_{*}$ model, hence the good performance reported previously. If
we think that it might just be the wrong ratio given the details of the
derivation of $k_{*}$, perhaps the relationship
$\Sigma_{\parallel}=(b_{g}(z)+f(z))*\Sigma_{\perp}$, where $b_{g}(z)$ is
estimated with the functional form

$b_{g}(z)=0.88z+0.49$ (36)

derived in Zhai et al. (2021b), is more appropriate. We show this ratio’s
prediction in the circles and find that they tend to be more in-line with
$\Sigma_{\parallel}$ measured. This would then give us a physical motivation
to either alter the ratio in the $k_{*}$ model or leave $k_{*}^{\perp}$ and
$k_{*}^{\parallel}$ free as is done in the
P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ model. Since we can not always
be sure what the proper functional form of $b_{g}(z)$ is, this one (Eq. 36)
was specifically derived for the galaxy mock used as the simulated
observation, it would be more robust to leave the perpendicular and parallel
damping parameters as free variables.

##  4 Conclusion

In this work, we have explored possible nonlinear corrections to the linear
predication of cold dark matter theory, in order to accurately model the
galaxy clustering signal in Fourier-space for the Roman High Latitude
Spectroscopic Survey (Wang et al., 2022). Roman is expected to observe $\sim
10M$ H$\alpha$ emission line galaxies between redshifts $1.0-2.0$, with an
H$\alpha$ flux $>10^{-16}$ [ergs/s/cm2] at a signal-to-noise ratio of
$6.5\sigma$. Redshifts from [OIII] will also be obtained by Roman at higher
redshifts, but we focus on the Hα survey only in this work.

We use a lightcone galaxy mock created using the semi-analytical model (SAM)
galacticus to paint galaxies onto a cosmological N-body simulation. The method
requires a tuning of the dust model to match observational data. The mock we
use adopts the dust model calibrated to the HiZELS H$\alpha$ luminosity
function. We make a cut in the H$\alpha$ flux at $>10^{-16}$ [ergs/s/cm2] in
redshift slices $z=[(1.0,1.2),(1.2,1.4),(1.4,1.6),(1.6,2.0)]$. We then measure
the power spectrum monopole and quadrupole multipole moments following FKP
methodology.

We utilize the code package camb to predict the linear power spectrum from
cosmological parameters for the redshift of interest, shift the model into
redshift-space with the inclusion of redshift-space distortions, i.e., Kaiser
squashing and Finger-of-God effect, and apply a window function to account for
the geometry of our simulated survey. We then assess the appropriate non-
linear corrections to this model, since although the universe becomes more and
more linear at higher redshifts, a purely linear model is not sufficient at
$>k_{max}=0.25\,h$/Mpc between $z=1.0-2.0$. Therefore, it is important that we
consider carefully how we make nonlinear corrections to the linear theory
since the statistical power available to constrain cosmological models on
large-scales is limited and if we want to improve our understanding of dark
energy via observations of the BAO signal and the growth-rate of large-scale
structure, we must push to smaller scales.

We evaluate two methods to account for the nonlinear evolution of the baryons
that smear out the BAO signal, a single parameter method P${}_{dw}(k|k_{*})$
that fixes the ratio between the smearing scale in the perpendicular and
parallel dimension and P${}_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ that
keeps the smearing scale free in each dimension. These two ‘de-wiggled’
methods also account for growth of nonlinear structure on small-scales. We
also evaluate two methods to account for the nonlinear growth of the density
field, one that emulates the behavior of the halos seen in N-body simulations
with the halo model as its analytical form $\mathcal{F}_{HM}$, and another
that emulates the behavior of the galaxies seen in a SAM created with N-body
simulations $\mathcal{F}_{SAM}$. We explore each of these independently and
also combinations of a de-wiggled power spectrum with a growth prefactor
resulting in $8$ total methods to account for nonlinearities in the galaxy
clustering signal.

To determine which model(s) are best able to reproduce the observed clustering
signal, we perform a recovery test, keeping the input cosmology fixed and
producing constraints on the AP-effect parameters $\alpha_{\perp}$ and
$\alpha_{\parallel}$ and a similar parameter for the linear growth parameter
$\alpha_{g}$ through application of an MCMC technique. If we have a correct
model then $\alpha_{i}=1$ ($i=\perp,\parallel,g$), since the observed
cosmology is the same as the fiducial cosmology. We use two quality
statistics, the reduced chi-square ($\chi^{2}_{\nu}$) and the Figure-of-Bias
(FoB), to determine which models are valid and most appropriate, respectively.
Since investigations of the two different ways to account for the Finger-of-
God effect show that $\mathcal{M}_{B}$, where a window function is added as a
prefactor to $\beta$ to select out the radial correlated infall pairwise
velocities, is better able to recover the input cosmologies Zhai et al.
(2021a), we have focused on the results using this model, $\mathcal{M}_{B}$
(see Eq.(16)).

Considering the $\chi^{2}_{r}$ statistic for measurements with
$k_{max}=0.25\,h$/Mpc, all models are valid with values close to unity.
Considering the FoB statistic for measurements with $k_{max}=0.25\,h$/Mpc, we
find that the three best models are $P_{dw}(k|k_{*})$,
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$, and
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$. If we want to suggest the two
best most unique models, we should choose either $P_{dw}(k|k_{*})$ or
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$, since one is a combination of the other.
Extending $k_{max}$ to $0.3\,h$/Mpc for these three best models we find that
it is $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ that performs the best
with $P_{dw}(k|k_{*})$ now second best. The difference between
$P_{dw}(k|k_{*})$ and $P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ is that
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$ reduces to $P_{dw}(k|k_{*})$
when $\Sigma_{\parallel}=(1+f_{g}(z))*\Sigma_{\perp}$.

From the posterior distribution of [$\Sigma_{\perp},\Sigma_{\parallel}$] we
find that the relationship is closer to
$\Sigma_{\parallel}=(b_{g}(z)+f_{g}(z))*\Sigma_{\perp}$, see Fig.6. That being
said, we cannot always a priori know the functional form of the galaxy bias
$b_{g}$, and so would suggest that keeping these two parameters free is a
better choice. The level of uncertainties on ($\alpha_{\perp}$,
$\alpha_{\parallel}$, $\alpha_{g}$) for these three models are approximately
at the ($2.45$, $3.66$, $8.28$), ($2.75$, $3.87$, $9.04$), and ($2.67$,
$4.24$, $8.28$) percent-level, for $P_{dw}(k|k_{*})$,
$P_{dw}(k|k_{*})*\mathcal{F}_{SAM}$, and
$P_{dw}(k|\Sigma_{\perp},\Sigma_{\parallel})$, respectively, when
$k_{max}=0.25\,h$/Mpc and at the ($2.43$, $3.54$, $8.21$), ($2.92$, $3.54$,
$9.00$), and ($2.65$, $4.08$, $8.71$) percent-level when $k_{max}=0.3\,h$/Mpc
The addition of $\mathcal{F}_{SAM}$ gives the $k_{*}$ model sufficient
flexibility to compensate for the assumed fixed ratio between the dampling
parameters but does not hold true as we push to smaller scales. This probably
suggests that the functional form of $\mathcal{F}_{SAM}$ is not entirely
physical and extrapolating to other forms of galaxy bias might render this
nonlinear correction invalid.

Our current analysis is a purely cosmological study with only a few parameters
modeling galaxy formation physics. At large scales, this may be sufficient and
the net impact can be described by a single galaxy bias parameter. However, as
we go to smaller scales, the impact due to galaxy formation is more
complicated and becomes non-negligible (McCarthy et al. 2019). The Galacticus
SAM mock that we have used provides a useful framework for such exploration
given its large volume and galaxy property parameters.

Our work is an extended analysis using the Roman SAM galaxy mock for the large
scale structure analysis based on two point statistics. It is anticipated that
constraining power will be significantly enhanced by adding higher order
statistics such as galaxy bispectrum. This is our ongoing work, and will be
presented elsewhere (McCarthy et. al., in prep).

The High Latitude Wide Area Spectroscopic Survey (HLWASS) that Roman will
execute will be determined in an open community process. The HLSS studied here
has served as the reference baseline in Roman mission development. We expect
that the tools that we have developed in this work, and continue to develop,
will be useful in probing dark energy and testing gravity using Roman HLWASS
data in an accurate and robust manner.

## Acknowledgments

KSM is supported by the NASA Postdoctoral Program. This work was carried out,
in part, by IPAC at the California Institute of Technology, and was sponsored
by NASA.

Software: Python, Matplotlib (Hunter, 2007), NumPy (van der Walt et al.,
2011), SciPy (Jones et al., 01), emcee (Foreman-Mackey et al., 2013), camb
(Lewis et al., 2000), nbodykit (Hand et al., 2018)

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Abbott et al. (2022) Abbott T. M. C., et al., 2022, Phys. Rev. D, 105, 043512
  * Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617
  * Alcock & Paczynski (1979) Alcock C., Paczynski B., 1979, Nature, 281, 358
  * Anderson et al. (2012) Anderson L., et al., 2012, MNRAS, 427, 3435
  * Anderson et al. (2014) Anderson L., et al., 2014, MNRAS, 441, 24
  * Angulo et al. (2008) Angulo R. E., Baugh C. M., Frenk C. S., Lacey C. G., 2008, MNRAS, 383, 755
  * Balaguera-Antolínez et al. (2019) Balaguera-Antolínez A., Kitaura F.-S., Pellejero-Ibáñez M., Zhao C., Abel T., 2019, MNRAS, 483, L58
  * Benson (2012) Benson A. J., 2012, New Astronomy, 17, 175
  * Benson et al. (2000) Benson A. J., Baugh C. M., Cole S., Frenk C. S., Lacey C. G., 2000, MNRAS, 316, 107
  * Beutler et al. (2011) Beutler F., et al., 2011, MNRAS, 416, 3017
  * Beutler et al. (2016) Beutler F., et al., 2016, Monthly Notices of the Royal Astronomical Society, 466, 2242
  * Beutler et al. (2017) Beutler F., et al., 2017, MNRAS, 464, 3409
  * Blake et al. (2011) Blake C., et al., 2011, MNRAS, 418, 1707
  * Blake et al. (2013) Blake C., et al., 2013, MNRAS, 436, 3089
  * Chuang et al. (2014) Chuang C.-H., Kitaura F.-S., Prada F., Zhao C., Yepes G., 2014, Monthly Notices of the Royal Astronomical Society, 446, 2621
  * Chuang et al. (2019) Chuang C.-H., et al., 2019, Monthly Notices of the Royal Astronomical Society, 487, 48
  * Cole et al. (2005) Cole S., et al., 2005, MNRAS, 362, 505
  * Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * Eggemeier et al. (2021) Eggemeier A., Scoccimarro R., Smith R. E., Crocce M., Pezzotta A., Sánchez A. G., 2021, Phys. Rev. D, 103, 123550
  * Eisenstein & Hu (1998) Eisenstein D. J., Hu W., 1998, The Astrophysical Journal, 496, 605
  * Eisenstein et al. (2005) Eisenstein D. J., et al., 2005, ApJ, 633, 560
  * Eisenstein et al. (2007) Eisenstein D. J., Seo H.-J., White M., 2007, ApJ, 664, 660
  * Feldman et al. (1994) Feldman H. A., Kaiser N., Peacock J. A., 1994, ApJ, 426, 23
  * Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gil-Marín et al. (2015) Gil-Marín H., Noreña J., Verde L., Percival W. J., Wagner C., Manera M., Schneider D. P., 2015, MNRAS, 451, 539
  * Gil-Marín et al. (2016) Gil-Marín H., et al., 2016, MNRAS, 460, 4188
  * Goldstein et al. (2022) Goldstein S., Pandey S., Slosar A., Blazek J., Jain B., LSST Dark Energy Science Collaboration 2022, Phys. Rev. D, 105, 123518
  * Guzzo et al. (2008) Guzzo L., et al., 2008, Nature, 451, 541
  * Hamilton (1998) Hamilton A. J. S., 1998, in , Astrophysics and Space Science Library. Springer Netherlands, pp 185–275, doi:10.1007/978-94-011-4960-0˙17, https://doi.org/10.1007%2F978-94-011-4960-0_17
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160
  * Hartlap, J. et al. (2007) Hartlap, J. Simon, P. Schneider, P. 2007, A&A, 464, 399
  * Heitmann et al. (2014) Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2014, ApJ, 780, 111
  * Howlett et al. (2012) Howlett C., Lewis A., Hall A., Challinor A., 2012, J. Cosmology Astropart. Phys., 2012, 027
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Joachimi et al. (2021) Joachimi B., et al., 2021, A&A, 646, A129
  * Jones et al. (01 ) Jones E., Oliphant T., Peterson P., et al., 2001–, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Kaiser (1987) Kaiser N., 1987, Monthly Notices of the Royal Astronomical Society, 227, 1
  * Kazin et al. (2014) Kazin E. A., et al., 2014, MNRAS, 441, 3524
  * Lange et al. (2022) Lange J. U., Hearin A. P., Leauthaud A., van den Bosch F. C., Guo H., DeRose J., 2022, MNRAS, 509, 1779
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Lue et al. (2004) Lue A., Scoccimarro R., Starkman G., 2004, Physical Review D, 69
  * Mandelbaum et al. (2013) Mandelbaum R., Slosar A., Baldauf T., Seljak U., Hirata C. M., Nakajima R., Reyes R., Smith R. E., 2013, Monthly Notices of the Royal Astronomical Society, 432, 1544
  * McCarthy et al. (2019) McCarthy K. S., Zheng Z., Guo H., 2019, MNRAS, 487, 2424
  * Mead et al. (2015) Mead A. J., Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015, Monthly Notices of the Royal Astronomical Society, 454, 1958
  * Mehta et al. (2015) Mehta V., et al., 2015, The Astrophysical Journal, 811, 141
  * Merson et al. (2018) Merson A., Wang Y., Benson A., Faisst A., Masters D., Kiessling A., Rhodes J., 2018, MNRAS, 474, 177
  * Mohammad & Percival (2022) Mohammad F. G., Percival W. J., 2022, MNRAS, 514, 1289
  * Percival & White (2009) Percival W. J., White M., 2009, Monthly Notices of the Royal Astronomical Society, 393, 297
  * Percival et al. (2001) Percival W. J., et al., 2001, MNRAS, 327, 1297
  * Percival et al. (2022) Percival W. J., Friedrich O., Sellentin E., Heavens A., 2022, MNRAS, 510, 3207
  * Perko et al. (2016) Perko A., Senatore L., Jennings E., Wechsler R. H., 2016, arXiv e-prints, p. arXiv:1610.09321
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Reid et al. (2014) Reid B. A., Seo H.-J., Leauthaud A., Tinker J. L., White M., 2014, MNRAS, 444, 476
  * Ross et al. (2015) Ross A. J., Samushia L., Howlett C., Percival W. J., Burden A., Manera M., 2015, MNRAS, 449, 835
  * Sánchez et al. (2008) Sánchez A. G., Baugh C. M., Angulo R. E., 2008, MNRAS, 390, 1470
  * Smith et al. (2003) Smith R. E., et al., 2003, Monthly Notices of the Royal Astronomical Society, 341, 1311
  * Sobral et al. (2012) Sobral D., Smail I., Best P. N., Geach J. E., Matsuda Y., Stott J. P., Cirasuolo M., Kurk J., 2012, Monthly Notices of the Royal Astronomical Society, 428, 1128
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Storey-Fisher et al. (2022) Storey-Fisher K., Tinker J., Zhai Z., DeRose J., Wechsler R. H., Banerjee A., 2022, arXiv e-prints, p. arXiv:2210.03203
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, The Astrophysical Journal, 761, 152
  * Wang (2008a) Wang Y., 2008a, Phys. Rev. D, 77, 123525
  * Wang (2008b) Wang Y., 2008b, Journal of Cosmology and Astroparticle Physics, 2008, 021
  * Wang & Steinhardt (1998) Wang L., Steinhardt P. J., 1998, ApJ, 508, 483
  * Wang et al. (2013) Wang Y., Chuang C.-H., Hirata C. M., 2013, Monthly Notices of the Royal Astronomical Society, 430, 2446
  * Wang et al. (2022) Wang Y., et al., 2022, ApJ, 928, 1
  * Wechsler & Tinker (2018) Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435
  * Yuan et al. (2022) Yuan S., Garrison L. H., Eisenstein D. J., Wechsler R. H., 2022, MNRAS, 515, 871
  * Zel’dovich (1970) Zel’dovich Y. B., 1970, A&A, 5, 84
  * Zhai et al. (2019) Zhai Z., et al., 2019, ApJ, 874, 95
  * Zhai et al. (2021a) Zhai Z., Chuang C.-H., Wang Y., Benson A., Yepes G., 2021a, MNRAS, 501, 3490
  * Zhai et al. (2021b) Zhai Z., Wang Y., Benson A., Chuang C.-H., Yepes G., 2021b, MNRAS, 505, 2784
  * Zhai et al. (2022) Zhai Z., et al., 2022, arXiv e-prints, p. arXiv:2203.08999
  * Zhang et al. (2013) Zhang P., Pan J., Zheng Y., 2013, Phys. Rev. D, 87, 063526
  * Zheng et al. (2013) Zheng Y., Zhang P., Jing Y., Lin W., Pan J., 2013, Phys. Rev. D, 88, 103510
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science Engineering, 13, 22