# Primordial non-Gaussianities with weak lensing: Information on non-linear
scales in the Ulagam full-sky simulations

Dhayaa Anbajagane (  )  Chihway Chang Hayden Lee  Marco Gatti

###### Abstract

Primordial non-Gaussianities (PNGs) are signatures in the density field that
encode particle physics processes from the inflationary epoch. Such signatures
have been extensively studied using the Cosmic Microwave Background, through
constraining the amplitudes, $f^{X}_{\rm NL}$, with future improvements
expected from large-scale structure surveys; specifically, the galaxy
correlation functions. We show that weak lensing fields can be used to achieve
competitive and complementary constraints. This is shown via the Ulagam suite
of N-body simulations, a subset of which evolves primordial fields with four
types of PNGs. We create full-sky lensing maps and estimate the Fisher
information from three summary statistics measured on the maps: the moments,
the cumulative distribution function, and the 3-point correlation function. We
find that the year 10 sample from the Rubin Observatory Legacy Survey of Space
and Time (LSST) can constrain PNGs to $\sigma(f_{\rm NL}^{\rm\,eq})\approx
110$, $\sigma(f_{\rm NL}^{\rm\,or,\,lss})\approx 120$, $\sigma(f_{\rm
NL}^{\rm\,loc})\approx 40$. For the former two, this is better than or
comparable to expected galaxy clustering-based constraints from the Dark
Energy Spectroscopic Instrument (DESI). The PNG information in lensing fields
is on non-linear scales and at low redshifts ($z\lesssim 1.25$), with a clear
origin in the evolution history of massive halos. The constraining power
degrades by $\sim\\!\\!60\%$ under scale cuts of $\gtrsim 20\,{\rm Mpc}$,
showing there is still significant information on scales mostly insensitive to
small-scale systematic effects (e.g., baryons). We publicly release the Ulagam
suite to enable more survey-focused analyses.

##  1 Introduction

An overarching goal of cosmology is to understand the history of the Universe,
both its initial state and its subsequent evolution to the present epoch. The
current paradigm for the generation of the initial density field (i.e. the
initial state) is inflation, a mechanism that generates quantum fluctuations
during an epoch of rapid exponential expansion (see Guth 2004, for a review).
This density field is then evolved under a $\Lambda\rm CDM$ model, where CDM
stands for cold dark matter and $\Lambda$ is the cosmological constant. The
large-scale structure — which is the distribution of matter in the Universe —
depends on both the initial conditions and their subsequent evolution, and is
thus a useful probe for studying the full history of the Universe. Analyses of
this structure, as traced at high redshift by the cosmic microwave background
(CMB) or at low redshift by galaxy surveys, have already shed light on the
properties of our Universe and on the values of the six parameters that make
up the $\Lambda\rm CDM$ model (Planck Collaboration et al. 2016; Asgari et al.
2021; Abbott et al. 2022; More et al. 2023). Other analyses have probed the
physics of the initial conditions, and in particular, have led to constraints
disfavoring certain classes of inflationary models (Planck Collaboration et
al. 2020b).

In the simplest models, only a single field — commonly called the “inflaton”
field — is present during inflation and slowly rolls down the potential,
resulting in simple, weak interactions. In this case, the initial density
field generated is highly Gaussian. Such a field is defined entirely by the
covariance of densities in different parts of the field, which is a Power
spectrum in Fourier space or a 2-point correlation function in real space.
Such functions capture the correlations between any two points in a field (or
two different fields) separated by a given distance. By adding complexity to
the inflation model — such as additional fields that interact with one
another, higher-order interactions within a single field, and so forth — the
inflationary mechanisms gain non-linear terms that then lead to primordial
non-Gaussianities (PNG) in the initial density field. The amplitudes of the
PNGs are captured by the $f_{\rm NL}$ parameters, and directly probe various
energy scales in the theory of inflation (Achúcarro et al. 2022, see their
Figure 1). Given that inflation might have taken place at very high energies
of $10^{14}\,\,{\rm GeV}$ (giga electron volts), which is much larger than
energy scales achievable in terrestrial particle physics experiments, the
parameter $f_{\rm NL}$ probes what has been denoted an “energy frontier” in
cosmology (Achúcarro et al. 2022). Note that $f_{\rm NL}$ only captures the
leading deviation from Gaussianity and corresponds to a bispectrum or three-
point correlation, which defines the correlation between three points in the
field (or three different fields) represented in either Fourier space or real
space, respectively. Such correlations can equivalently be interpreted as
skewness in the field. Higher-order correlations, such as the 4-point function
parameterized by $\tau_{\rm NL}$, can also be non-zero due to inflationary
mechanisms. However, they are not the focus of this work.

The current best constraints on $f_{\rm NL}$ are found in the bispectrum
analysis of the Planck CMB data (Planck Collaboration et al. 2020a) —
$\sigma(f_{\rm NL}^{\rm\,eq})=5$, $\sigma(f_{\rm NL}^{\rm\,eq})=47$,
$\sigma(f_{\rm NL}^{\rm\,or,\,cmb})=24$ — as the spatial correlations of the
observed temperature fluctuations arise from inflationary correlations. These
constraints are primarily limited by cosmic variance, rather than by
measurement noise. Future CMB surveys such as CMB S4 (Abazajian et al. 2019)
can only modestly improve on this result through reduction of the measurement
noise; the Planck data already covers most of the observable sky, thus CMB-S4
will not improve on the cosmic variance limit. However, more significant
improvements are expected from large-scale structure (LSS) galaxy surveys,
using the correlations of galaxy positions as a probe of the inflationary
correlations. These surveys probe a 3D volume, resulting in the number of
countable modes scaling as $N_{k}\propto Vk_{\rm max}^{3}$, while in a 2D CMB
map, the number of modes scales as $N_{k}\propto Ak_{\rm max}^{2}$. Here, $V$
is the 3D survey volume, $A$ is the 2D map area, and $k_{\rm max}$ is the
highest $k$-mode studied in the analysis. As galaxy surveys push to higher
redshift (which increases the survey volume) and higher resolution (which
improves the smallest measured scale), the number of measured modes increases
significantly. Thus, galaxy survey measurements will have superior statistical
power to CMB measurements, and can then be used to constrain inflationary
physics. In addition, galaxy correlations probe a noticeably different set of
length scales than the CMB, and the combination of the two can probe scale-
dependent $f_{\rm NL}$ models (see Section 5.1). Many works have extracted PNG
constraints from galaxy correlation measurements (Mueller et al. 2021; Cabass
et al. 2022b, a; D’Amico et al. 2022; Philcox et al. 2022a).

Another cosmological probe observed by many of the same widefield galaxy
surveys is weak lensing, which is distortions in the shapes of galaxies —
commonly denoted “source galaxies” — due to the foreground structure present
between the observer and the galaxies (Bartelmann & Schneider 2001). The
spatial correlations of these distortions are generated by the foreground
density field and are thus a probe of cosmology. While all observational
constraints on $f_{\rm NL}$ from galaxy surveys have focused on the 3D galaxy
field, none have focused on the weak lensing field (though there exists some
theoretical work in this direction). There are a number of complementary
benefits in using weak lensing as a probe of inflation, such as the lensing
field being a direct, unbiased tracer of the density field insensitive to the
physics of the galaxy–halo connection,111As we will discuss later in Section
5.2, weak lensing is still sensitive to baryonic physics given the latter
impacts the density field that generates the weak lensing signal. However,
this is a distinctly different kind of phenomenon from the galaxy–halo
connection, which concerns itself with the distribution of all galaxies that
can populate a given halo, and thus involves smaller-scale physics. efficient
simulation-based modeling of strongly non-linear scales enabled by less
stringent resolution requirements, and so on (a more detailed discussion is
presented in Section 5.1). These advantages are currently not being utilized
in studies of inflationary physics as lensing-based analyses are yet to be
implemented.

Historically, a limitation in performing such lensing-based studies has been
obtaining a computationally efficient model for $f_{\rm NL}$ signatures in
weak lensing. Given that the weak lensing signal probes a line-of-sight
integral of the density field, a majority of the measurements contain some
contributions from the non-linear density field (e.g., Secco et al. 2022a, see
their Figure 4). Galaxy correlations from $f_{\rm NL}$ signals have been
efficiently modelled using the effective field theory of large-scale structure
(EFTofLSS), which is an analytic approach to modeling the correlations of the
density field and the galaxy field (Baumann et al. 2012; Carrasco et al.
2012). The current calculations of the two-loop power spectrum and one-loop
bispectrum are accurate up to quasi-linear scales and have significant
deviations ($\geq 10\%$) in the non-linear regimes ($k\gtrsim 1h/{\rm Mpc}$);
for example, see Baldauf et al. (2015, their Figure 10) or Sefusatti et al.
(2010, their Figures 2-7). Therefore, it is difficult to employ the EFT
approach to model weak lensing. However, over the past decade, the feasibility
of full, simulation-based modeling has grown significantly and has led to
multiple analyses of the lensing field that are simulation-based (e.g., Fluri
et al. 2018; Fluri et al. 2019; Zürcher et al. 2021; Fluri et al. 2022;
Zürcher et al. 2022), or more often at least simulation-informed (e.g., Secco
et al. 2022a; Amon et al. 2022; Gatti et al. 2022). Thus, with modern
advancements in computing, it is now possible to efficiently and accurately
model the non-linear density field through N-body simulations (see Angulo &
Hahn 2022, for a review), which thereby enables analyses of $f_{\rm NL}$ with
weak lensing.

Previous works have theoretically explored the power of weak lensing in
constraining PNGs, and have found it to be a potentially promising probe
(Marian et al. 2011; Shirasaki et al. 2012; Hilbert et al. 2012). These works
used a modest number of simulations to estimate the signal of a specific type
of $f_{\rm NL}$, called the “local” type (see Section 2) on the lensing field,
and used simple scaling arguments for how constraining power increases with
survey area to approximately estimate the constraints from wide-field lensing
surveys. However, as discussed above, it is now possible to run substantially
larger suites of large-volume high-resolution simulations of the full sky and
estimate the inflationary information in the lensing field. In addition, there
are other types of $f_{\rm NL}$ — beyond the local type — that have valuable
information about inflation and have not been considered in simulation-based
analyses of weak lensing, though some analytical approaches have been
previously employed (Schäfer et al. 2012; Giannantonio et al. 2012).

In this work, we explore the use of the lensing convergence field as a probe
of PNGs. We expand on previous efforts by (i) developing and using a large
simulation suite ($N_{\rm sims}=3600$), which enables better numerical
accuracy of Fisher information estimates and allows a closer match of lensing
survey specifications, (ii) exploring four types of $f_{\rm NL}$, each of
which probes different primordial physics and three of which are being studied
for the first time in the context of simulation-based constraints from weak
lensing on non-linear scales, and (iii) forecasting realistic constraints for
current and upcoming widefield surveys using their expected redshift
distributions, noise amplitudes, survey area etc. This work is organized as
follows: in §2 we describe the suite of simulations developed for this work,
including the types of PNGs we focus on, and also the forward modeling
procedure to simulate the weak lensing observations from each survey. In §3 we
describe the statistics used to summarize the lensing field, which includes
the moments, cumulative distribution functions, and the three-point function.
We present our results in §4, including the Fisher information in different
statistics and surveys, and the physical origin of the PNG signal in lensing.
We discuss the advantages and caveats of lensing-based analyses of inflation
in §5, and then conclude in §6.

##  2 Simulations

The PNGs, by virtue of their impact on the initial density field, can affect
non-linear structure formation and imprint onto any field related to this
structure, such as the galaxy number density fields as in the studies
discussed above. In this work, we are interested in the lensing convergence
field, $\kappa$, which is a line-of-sight integral of the density field,

$\kappa(\mathbf{\hat{n}},z_{s})=\frac{3}{2}\frac{H_{0}^{2}\Omega_{\rm
m}}{c^{2}}\int_{0}^{z_{s}}\\!\\!\\!\delta(\mathbf{\hat{n}},z_{j})\frac{\chi_{j}(\chi_{s}-\chi_{j})}{a(z_{j})\chi_{s}}dz_{j}\frac{d\chi}{dz}\bigg{|}_{z_{j}},$
(2.1)

where $z_{s}$ is the redshift of the “source” plane/galaxies being lensed,
$\mathbf{\hat{n}}$ is the pointing direction on the sky, $\delta$ is the
overdensity field, $\chi$ is the comoving distance from an observer to a given
redshift, $a$ is the scale factor, $H_{0}$ is the Hubble constant,
$\Omega_{\rm m}$ is the matter energy density fraction at $z=0$, and $c$ is
the speed of light. We use the shorthand $\chi(z_{s})\equiv\chi_{s}$ and
$\chi(z_{j})\equiv\chi_{j}$.

We model this convergence — including its dependence on PNGs and cosmology —
using full-sky density maps from N-body simulations. Such simulations are
uniquely suited for modeling these fields in the non-linear regime. The set of
simulations introduced in this work will henceforth be referred to as the
Ulagam suite.222Ulagam is a Tamil word (pronounced “ooh-luh-gum”) that denotes
the World or the Universe, and this naming choice is inspired by the Uchuu
simulations — a set of multi-${\rm Gpc}$, high-resolution simulations — which
were named after the Japanese word for Universe. The simulations are run with
the Pkdgrav3 solver (Potter et al. 2017), which has been used extensively in
modeling the weak lensing field (Fluri et al. 2018; Fluri et al. 2019; Zürcher
et al. 2021; Zürcher et al. 2022; Gatti et al. 2022; Kacprzak et al. 2023;
Gatti et al. 2023). Pkdgrav3 automatically builds lightcones while solving the
gravitational dynamics of the system forward in time, and so our final outputs
are the lightcone shells — i.e. Healpix maps — of the density field at
different redshifts. The simulation box is tiled/repeated as needed to
construct large enough volumes to then build full-sky lightcones to a given
redshift. This repetition will bias any large-scale correlations in the
lightcone, but in this work we only consider scales much smaller than the box
size.

These simulations are run in $L=1\,{\rm Gpc}/h\approx 1.5\,{\rm Gpc}$ boxes,
starting at $z=127$, and with $N=512^{3}$ dark matter particles. The initial
conditions for all runs are obtained from the Quijote suite (Villaescusa-
Navarro et al. 2020) and the Quijote-png extension (Coulton et al. 2022).
Thus, these simulations are lightcone companions of the Quijote simulations
and are specialized for widefield survey analyses. The original Quijote suite
was designed for studying the Fisher information of non-linear structure, as
well as for building emulators sampling different cosmological parameters, but
the available data products provide inadequate redshift resolution for
producing accurate mock lightcones of the lensing/density field. Hence we have
resimulated a subset of these simulations to create accurate full-sky density
and lensing maps. When running simulations with PNGs, these PNGs are included
in the density field of the initial conditions — using the techniques
described below in Section 2.1 — and the field is then evolved via the
fiducial N-body solver with no modifications.

The Ulagam suite contains simulations for computing the derivatives of the
lensing/density field with respect to multiple $w{\rm CDM}$ and PNG
parameters; these are $\Omega_{\rm m}$, $\sigma_{8}$, $w$, and $n_{s}$ for the
$\Lambda{\rm CDM}$ case and different $f_{\rm NL}$ corresponding to four types
— local, equilateral, orthogonal LSS and orthogonal CMB — which are detailed
in Section 2.1. The suite contains 100 simulations per parameter where the
value of that parameter is slightly higher than the fiducial, and another 100
simulations where the value is slightly lower than the fiducial, and these
paired sets are used to compute the derivatives. The fiducial cosmology is
from Planck Collaboration et al. (2016), and the derivatives are computed over
differences of $\Delta\Omega_{\rm m}=0.02$, $\Delta\sigma_{8}=0.03$, $\Delta
w=0.05$, $\Delta n_{\rm s}=0.04$, and $\Delta f_{\rm NL}^{\rm type}=200$,
which are all the same settings as the Quijote and Quijote-PNG suites. The
Ulagam suite also contains $2000$ simulations at the fiducial cosmology which
can be used to compute the covariance matrix for data-vectors. Each all-sky
map can have multiple, completely independent cutouts of the survey
footprints, so in practice, the simulations provide between 6000 to 8000
realizations for the covariance, and also between 300 to 400 realizations for
the derivatives; the lower and upper bound numbers are for LSST and DES
respectively. A summary of the available full-sky runs is listed in Table 1.

The simulations have a total of 100 timesteps/shells, with 95 shells between
$0<z<10$, and a high redshift resolution of $\Delta z\approx 0.01-0.05$ in the
latter redshift range; the exact value of $\Delta z$ depends on the shell. The
timesteps in this redshift range are spaced uniformly in proper time, $t$, and
this corresponds to different $z$ and comoving distances depending on the
cosmology. The density shells are then post-processed via Equation 2.1, with
the integral over $z_{j}$ now replaced by a simple discrete sum, to create a
lensing convergence field at each source plane redshift, $z_{s}$. This
technique uses the Born approximation, which computes the total effective
deflection due to lensing but along an undeflected ray path. A more precise
calculation uses full raytracing, which calculates these deflections while
also constantly deflecting/updating the ray path. Petri et al. (2017) found
the Born approximation leads to differences of $\lesssim 5\%$ for the third
moments statistic we will use in Section 4. This effect is important when
requiring a certain absolute accuracy in the simulation predictions, whereas
for estimating the Fisher information — as we do in this work — these
requirements can be relaxed given we only require accuracy in the relative
differences between different simulations (as discussed further below).

Halos are identified in the 3D simulation volume using a friends-of-friends
(FoF) percolation technique with a linking length of $b=0.2$, in units of the
mean inter-particle distance, consistent with the choice made for Quijote. The
halo catalogs contain all identified objects with at least 100 particles,
which corresponds to a minimum mass of $M>10^{14}\,{\rm M}_{\odot}$ across all
redshifts. This is different from the Quijote catalogs, which keep all halos
down to 20 particles. The change was made to reduce the total storage
footprint of the simulation suite.

While the original Quijote suite was run using Gadget3 (last described in
Springel 2005), we use Pkdgrav3 in this work given its specialization and
extensive use in modeling full-sky observables. We verify in Figure 9 that we
can reproduce the results from the Quijote simulations to within expected
accuracy given the difference in the two N-body solvers. Similarly, we have
also performed checks on the choices of particle count in Appendix A. The
numerical requirements for this work are less stringent than a full
simulation-based model as we do not use the simulations for cosmological
inference, but rather for (i) computing derivatives for the Fisher
information, where the relevant quantities are relative differences in the
simulations as we vary cosmological parameters, and for (ii) computing
covariance matrices for the Fisher information, where once again the relevant
quantities are relative differences in simulations across different
realizations. As a result, our requirements for absolute accuracy/calibration
of the simulations are relaxed. Thus, while some results in Appendix A suggest
the simulations’ accuracy can benefit from a higher particle count, the
current suite is still adequate for estimating the Fisher information — which
was the original purpose of the Quijote simulations that the Ulagam suite now
builds on.

Run | $\mathbf{P_{\rm fid}}\pm\Delta P$ | $N_{\rm sim}$
---|---|---
Fiducial | — | 2000
Local PNG, $f_{\rm NL}^{\rm\,loc}$ | $\mathbf{0}\pm 100$ | 100
Equilateral PNG, $f_{\rm NL}^{\rm\,eq}$ | $\mathbf{0}\pm 100$ | 100
LSS Orthogonal PNG, $f_{\rm NL}^{\rm\,or,\,lss}$ | $\mathbf{0}\pm 100$ | 100
CMB Orthogonal PNG, $f_{\rm NL}^{\rm\,or,\,cmb}$ | $\mathbf{0}\pm 100$ | 100
Matter density, $\Omega_{\rm m}$ | $\mathbf{0.3175}\pm 0.01$ | 100
Density fluctuations amplitude. $\sigma_{8}$ | $\mathbf{0.834}\pm 0.015$ | 100
Dark energy EoS $w_{0}$ | $\mathbf{-1}\pm 0.05$ | 100
Spectral index $n_{s}$ | $\mathbf{0.9624}\pm 0.02$ | 100
Table 1: The simulation runs presented in this work. The fiducial simulation
parameters follow those of Planck Collaboration et al. (2016) and are shown in
bold, while the variations to the parameters for calculating derivatives are
shown as the $\pm\Delta P$ values. The cosmological parameter values used in
the fiducial runs are all the bolded values. We always assume a flat cosmology
with $\Omega_{\Lambda}=1-\Omega_{\rm m}$. The parameters not shown above are
$h=0.6711$ and $\Omega_{\rm b}=0.049$.

###  2.1 Initial conditions with primordial non-Gaussianities

Generating initial conditions for a purely Gaussian initial density field is a
well-studied procedure with established numerical recipes (e.g., Crocce et al.
2006). Generating those for a field with PNGs, however, requires careful
transformations of the Gaussian field. The initial conditions of Quijote-Png,
which we use in this work for our PNG simulations, are generated using the
methodology of Scoccimarro et al. (2012). We briefly summarize this process
below for the four different PNGs (see Chen (2010); Achúcarro et al. (2022)
for reviews on inflation-driven PNGs) we consider in this work, as described
in Coulton et al. (2022).

In general, given some Gaussian initial conditions for the gravitational
potential, $\phi(\mathbf{x})$, or its Fourier equivalent $\phi(\mathbf{k})$,
we can generate a field, $\Phi$, with a chosen bispectrum as

$\Phi(\mathbf{k})=\phi(\mathbf{k})+\int f_{\rm
NL}[\delta_{D}]K(\mathbf{k}_{1},\mathbf{k}_{2})\phi(\mathbf{k}_{1})\phi(\mathbf{k}_{2})d^{3}k_{1}d^{3}k_{2}\,,$
(2.2)

where $[\delta_{D}]=\delta_{D}(\mathbf{k}-\mathbf{k}_{1}-\mathbf{k}_{2})$ is
the Dirac delta function enforcing momentum conservation and
$K(\mathbf{k}_{1},\mathbf{k}_{2})$ is a coupling kernel that contains
information about the chosen bispectrum. By choosing different
$K(\mathbf{k}_{1},\mathbf{k}_{2})$, we can generate density fields with
different PNGs.

The bispectrum of the field $\Phi(\mathbf{k})$ defined in Equation 2.2 above
is

$B_{\Phi}=2f_{\rm
NL}K(\mathbf{k}_{1},\mathbf{k}_{2})P_{\Phi,1}P_{\Phi,2}+\text{cyc.}$ (2.3)

Given a particular bispectrum template, one can find the coupling kernel,
$K(\mathbf{k}_{1},\mathbf{k}_{2})$, that transforms the Gaussian field so as
to imprint a chosen bispectrum. In this work, we focus on the same four PNG
templates used in Coulton et al. (2022).

First, the local-type $\boldsymbol{f_{\rm NL}^{\rm\,loc}}$ is the most well-
studied $f_{\rm NL}$ in the context of LSS and of simulation-based studies.
This is largely due to the simplicity in generating the associated initial
conditions, which can be done entirely in real-space. The bispectrum of the
field goes as,

$B^{\rm loc}_{\Phi}(k_{1},k_{2},k_{3})=2f_{\rm
NL}^{\rm\,loc}P_{\Phi}(k_{1})P_{\Phi}(k_{2})+2\,\text{perms.}\,,$ (2.4)

and the field itself can be generated easily by adding the square of the
Gaussian field to the linear term,

$\Phi^{\rm loc}(\mathbf{x})=\phi(\mathbf{x})+f_{\rm
NL}^{\rm\,loc}\big{[}\phi(\mathbf{x})^{2}-\langle\phi(\mathbf{x})^{2}\rangle\big{]}\,,$
(2.5)

where the ensemble average must be subtracted out to enforce that the field,
$\phi(\mathbf{x})^{2}$, has zero mean and is thus a purely perturbative field
that does not alter the mean value of $\Phi^{\rm loc}(\mathbf{x})$. The
bispectrum of this model peaks at “squeezed” configurations, $k_{1}\ll
k_{2},k_{3}$, and can be generated by the presence of a second, light scalar
field during inflation, often called a “curvaton” (Moroi & Takahashi 2001;
Enqvist & Sloth 2002; Lyth & Wands 2002; Sasaki et al. 2006). Such a
bispectrum could also be generated during reheating — a process during which
inflatons decay into the standard model particles — via a fluctuating,
inhomogeneous decay rate (Kofman 2003; Dvali et al. 2004a, b).

Second, the equilateral-type $\boldsymbol{f_{\rm NL}^{\rm\,eq}}$ is a “non-
local” PNG since it cannot be generated as local real-space transforms of the
initial Gaussian field. It was derived in Senatore et al. (2010, see their
Equation 3.1) and has a bispectrum of the form

$\displaystyle B^{\rm eq}_{\Phi}(k_{1},k_{2},k_{3})=6f_{\rm
NL}^{\rm\,eq}\bigg{[}$ $\displaystyle-P_{\Phi}(k_{1})P_{\Phi}(k_{2})+\textit{
2
\text{perm.}}-2\bigg{(}P_{\Phi}(k_{1})P_{\Phi}(k_{2})P_{\Phi}(k_{3})\bigg{)}^{2/3}$
$\displaystyle+P_{\Phi}(k_{1})^{1/3}P_{\Phi}(k_{2})^{2/3}P_{\Phi}(k_{3})+\textit{
5 \text{perm.}}\bigg{]}.$ (2.6)

Such PNGs are generated from inflation models that have “non-canonical”
kinetic terms, and their amplitude peaks in the limit $k_{1}\approx
k_{2}\approx k_{3}$. This template approximates the bispectrum that arises
from leading-derivative cubic interactions in the effective field theory (EFT)
of inflation (Cheung et al. 2008). Prototypical models with a non-canonical
kinetic term and a subluminal sound speed include Dirac-Born-Infeld, or DBI,
inflation (Silverstein & Tong 2004; Alishahiha et al. 2004) or k-inflation
(Armendáriz-Picón et al. 1999). The corresponding real-space expression at the
field-level is

$\displaystyle\Phi^{\rm eq}(\mathbf{x})=\phi+f_{\rm
NL}^{\rm\,eq}\bigg{[}-3\phi^{2}+4\partial^{-1}(\phi\partial\phi)+2\nabla^{-2}(\phi\nabla^{2}\phi)+2\nabla^{-2}(\partial\phi)^{2}\bigg{]}\,.$
(2.7)

Third, the CMB Orthogonal-type $\boldsymbol{f_{\rm NL}^{\rm\,or,\,cmb}}$ is
also a non-local PNG template derived in Senatore et al. (2010, see their
Equation 3.2), which has a shape that is approximately orthogonal to both
local-type and equilateral-type PNG. Together with the equilateral type, the
two shapes cover the parameter space spanned by the two leading-derivative
cubic interactions in the EFT of inflation. The template takes the form

$\displaystyle B^{\rm or,cmb}_{\Phi}(k_{1},k_{2},k_{3})=6f_{\rm
NL}^{\rm\,or,\,cmb}$
$\displaystyle\,\,\bigg{[}-3P_{\Phi}(k_{1})P_{\Phi}(k_{2})+\textit{ 2
\text{perm.}}-8\bigg{(}P_{\Phi}(k_{1})P_{\Phi}(k_{2})P_{\Phi}(k_{3})\bigg{)}^{2/3}$
$\displaystyle+3P_{\Phi}(k_{1})^{1/3}P_{\Phi}(k_{2})^{2/3}P_{\Phi}(k_{3})+\textit{
5 \text{perm.}}\bigg{]}\,.$ (2.8)

which leads to the real-space expression

$\displaystyle\Phi^{\rm or,cmb}(\mathbf{x})=\phi+f_{\rm
NL}^{\rm\,or,\,cmb}\bigg{[}-9\phi^{2}+10\partial^{-1}(\phi\partial\phi)+8\nabla^{-2}(\phi\nabla^{2}\phi)+8\nabla^{-2}(\partial\phi)^{2}\bigg{]}\,.$
(2.9)

Fourth and finally, the LSS Orthogonal-type  $\boldsymbol{f_{\rm
NL}^{\rm\,or,\,lss}}$ is another template derived in Senatore et al. (2010,
see their Appendix B), which is also orthogonal to both local-type and
equilateral-type PNG like $f_{\rm NL}^{\rm\,or,\,cmb}$. When considering the
squeezed limit, it is a better approximation to the true bispectrum shape —
where the true shape is determined by the EFT of inflation — when compared to
the CMB orthogonal type. The bispectrum here is written as

$\displaystyle B^{\rm or}_{\Phi}(k_{1},k_{2},k_{3})=$
$\displaystyle\,\,6f_{\rm
NL}^{\rm\,or,\,lss}\bigg{(}P_{\Phi}(k_{1})P_{\Phi}(k_{2})P_{\Phi}(k_{3})\bigg{)}^{2/3}\bigg{[}\frac{p}{27}\frac{k_{1}^{4}}{k_{2}^{2}k_{3}^{2}}+\textit{2
perms.}$ $\displaystyle-\frac{20p}{27}\frac{k_{1}k_{2}}{k_{3}^{2}}+\textit{2
perms.}-\frac{6p}{27}\frac{k_{1}^{3}}{k_{2}k_{3}^{2}}+\textit{5 perms.}$
$\displaystyle+\frac{15p}{27}\frac{k_{1}^{2}}{k_{3}^{2}}+\textit{5
perms.}+\bigg{(}1+\frac{9p}{27}\bigg{)}\frac{k_{3}^{2}}{k_{1}k_{2}}+\textit{2
perms.}$
$\displaystyle+\bigg{(}1+\frac{15p}{27}\bigg{)}\frac{k_{1}}{k_{3}}+\textit{5
perms.}-\bigg{(}2+\frac{60p}{27}\bigg{)}\bigg{]}\,,$ (2.10)

where the constant $p$ is given by

$p=\frac{27}{-21+\frac{743}{7(20\pi^{2}-193)}}\,.$ (2.11)

The real-space expression for this template is lengthy and so instead of
reproducing it here, we direct readers to Equation A11 of Coulton et al.
(2022).

The “non-local” PNG amplitudes — $f_{\rm NL}^{\rm\,eq},f_{\rm
NL}^{\rm\,or,\,lss},f_{\rm NL}^{\rm\,or,\,cmb}$ — depend on the self-coupling
of the inflationary perturbations. There is a natural theoretical threshold
for these PNG at $f_{\rm NL}\sim 1$. If $f_{\rm NL}$ is much larger than
unity, then the inflationary theory is favored to be strongly coupled, which
in turn disfavors the simplest single-field, slow-roll model. Thus,
constraining these $f_{\rm NL}$ would have profound implications for
understanding the physics of inflation. Given the formalism of Equation 2.2,
one could define other templates as well, each corresponding to a different
bispectrum signature that probes different interactions in the inflationary
field. However, a number of models will have bispectrum shapes that overlap
either/both of the local and equilateral PNG templates. The inclusion of two
more templates in this work — $f_{\rm NL}^{\rm\,or,\,lss}$ and $f_{\rm
NL}^{\rm\,or,\,cmb}$ — further expands the breadth of models we can probe.

Note that all templates above are designed to induce a specific bispectrum in
the initial density field. These templates may induce additional, unintended
corrections to the power spectrum, trispectrum (the Fourier space version of
the 4-point correlation function), and higher-order spectra. Coulton et al.
(2022, see their Figure 1) show that the impact of $f_{\rm NL}$ on the
primordial power spectra is negligible — which is by construction as the
method requires corrections to the power spectrum to be subdominant
(Scoccimarro et al. 2012) — while Jung et al. (2023, see their Appendix A)
also show that the templates generate no unphysical trispectra. In summary,
all the signals we study further below are physical and not artifacts of the
PNG generation procedure.

###  2.2 Simulating DES and LSST skies

Figure 1: The redshift distribution of the signal in each tomographic bin for
each survey. The DES Y3 distributions are from Myles et al. (2021), and have
been smoothed with a narrow Gaussian kernel for visualization purposes only.
The colored numbers are the mean redshift of each bin. The bottom panel shows
the quantity, $W(z_{j})=\chi_{j}(\chi_{s}-\chi_{j})/({a(z_{j})\chi_{s}})$ as
used in Equation 2.1, with the different colors showing different $z_{s}$. The
lines terminate at the redshift of the source plane, $z_{s}$. LSST has a tail
to high redshifts that is likely unrealistic, but Figure 6 below shows our
results are insensitive to the presence of this tail.

In this work, we are interested in the impact of $f_{\rm NL}$ as observed by a
weak lensing survey. For this, we must construct lensing maps. This can be
done by using the density fields of N-body simulations to create lensing
convergence fields, and then post-processing these fields to match the
observed data. Various aspects of these procedures have been utilized in
existing analyses/forecasts of weak lensing data (Fluri et al. 2019; Zürcher
et al. 2021; Fluri et al. 2022; Gatti et al. 2022; Zürcher et al. 2022; Gatti
et al. 2023; Anbajagane et al. 2023a). The two surveys we focus on are: (i)
the Dark Energy Survey (DES, The Dark Energy Survey Collaboration 2005), which
is an optical imaging survey of 5,000 deg2 of the southern sky, and is
currently the largest precision photometric dataset for cosmology. The Year 3
data products and cosmology results are available (Sevilla-Noarbe et al. 2021;
Abbott et al. 2022), while the legacy Year 6 dataset is not yet available at
the time of publishing this work. (ii) the Rubin Observatory Legacy Survey of
Space and Time (LSST), which is a 14,000 deg2 survey that probes higher
redshifts, and is the successor to current weak lensing surveys. We detail
below the exact steps in our forward modeling procedure to make lensing fields
corresponding to these surveys:

Constructing lensing convergence shells. The main simulation product used in
this work is lightcone shells of the particle counts, which are a set of two-
dimensional, HEALPix maps of projected particle counts at different redshifts.
The particle count in pixel $i$ is trivially converted into the overdensity
field as

$\delta^{i}=N^{i}_{\rm p}/\langle N_{\rm p}\rangle-1,$ (2.12)

where $N^{i}_{\rm p}$ is the number of particles in pixel $i$, and the average
is computed over all pixels in the shell. The density shells can then be
converted into the convergence $\kappa$ using Equation 2.1, after converting
the integral over redshift into a discrete sum over lightcone shells.

Source galaxy redshift distributions. Once we have convergence shells at
different redshifts, we construct the convergence within the different
tomographic bins of each survey. This binned convergence is computed as a
weighted average, where the weights are the source galaxy $n(z)$ of the chosen
bin and survey,

$\kappa^{A}(\mathbf{\hat{n}})=\sum_{j=1}^{N_{\rm
steps}}n^{A}(z_{j})\kappa(\mathbf{\hat{n}},z_{j})\Delta z,$ (2.13)

where $\kappa^{A}$ is the true convergence of a tomographic bin, $A$. The
$n(z)$ is obtained from the following: for DES Year 3 we use the results from
Myles et al. (2021), for DES Year 6, LSST Year 1 and Year 10 we use the same
$n(z)$ utilized in Zhang et al. (2022, see their Table 1). The LSST modeling
in that work follows the baseline analysis choices of The LSST Dark Energy
Science Collaboration et al. (2018, see their Appendix D2.1). The redshift
distributions for DES Y6, and LSST Y1 and Y10 are parameterized as,

$\frac{dN}{dz}\propto
z^{2}\exp\bigg{[}-\bigg{(}\frac{z}{z_{0}}\bigg{)}^{\alpha}\bigg{]},$ (2.14)

with parameters given in Table 2. Once the $n(z)$ of the full survey is
defined, we split it into 4 (5) tomographic bins for DES Y6 (LSST Y1/Y10) of
equal number density. Each bin is then convolved with a Gaussian of width
given by the photometric redshift uncertainty, also quoted in Table 2. The
final $n(z)$ for each survey is shown in Figure 1. The DES Y6 distribution is
non-zero only between $0.2<z<1.3$, following Zhang et al. (2022, see their
Table 1). The LSST distributions are cut at $z<3.5$. The fifth LSST bin
(purple line) has a tail to high redshifts that is likely unrealistic.
However, we will show below (in Figure 6) that this bin has a negligible
impact on our final constraints.

Survey | $(z_{0},\alpha)$ | $n_{\rm gal}/{\rm arcmin}^{2}$ | $\sigma_{z}$
---|---|---|---
DES Y6 | (0.13, 0.78) | 9 | 0.1(1 + z)
LSST Y1 | (0.13, 0.78) | 10 | 0.05(1 + z)
DES Y10 | (0.11, 0.68) | 27 | 0.05(1 + z)
Table 2: The redshift distribution and source galaxy number density assumed
for the upcoming surveys. All numbers are taken from Zhang et al. (2022, see
their Table 1).

Constructing lensing shear shells. Weak lensing surveys do not directly
observe the convergence, $\kappa$, as their main observable is the galaxy
shapes that are tracers of the shear field, $\gamma$.333Formally, the galaxy
ellipticities trace the reduced shear, $e=\gamma/(1-\kappa)$. See Section 5.2
for more details. The shear and convergence field can be transformed between
each other using the Kaiser-Squires (KS) transform (Kaiser & Squires 1993),
implemented in harmonic space as

$\gamma^{\ell m}_{E}+i\gamma^{\ell
m}_{B}=-\sqrt{\frac{(\ell+2)(\ell-1)}{\ell(\ell+1)}}\bigg{(}\kappa^{\ell
m}_{E}+i\kappa^{\ell m}_{B}\bigg{)},$ (2.15)

where $X_{\\{E,B\\}}$ are the E-mode and B-mode (or Q and U polarizations, in
HealPix notation) of the field.

Intrinsic alignments (IA). Even in the absence of weak lensing, the observed
galaxy shapes will have non-zero correlation due to the presence of a tidal
gravitational field aligning the galaxy orientations. This effect is called
intrinsic alignments (Troxel & Ishak 2015; Lamman et al. 2023). Under
perturbation theory, the leading-order contribution of this effect is

$\kappa_{\rm IA}(\mathbf{\hat{n}},z)=-\frac{A_{\rm IA}\rho_{\rm
crit,0}\Omega_{\rm m}}{D_{+}(z,\Omega_{\rm
m})}\bigg{(}\frac{1+z}{1+0.6}\bigg{)}^{\eta_{\rm
IA}}\delta(\mathbf{\hat{n}},z),$ (2.16)

where $\kappa_{\rm IA}(\mathbf{\hat{n}},z)$ is the convergence signal
corresponding to IA, $\delta(\mathbf{\hat{n}},z)$ is the density field,
$A_{\rm IA}$ is the amplitude of the IA effect, $\rho_{\rm crit}$ is the
critical density of the universe at $z=0$, $\Omega_{\rm m}$ is the fraction of
matter energy density at $z=0$, $D_{+}$ is the linear growth function
normalized to $D_{+}(z=0)=1$, and $\eta_{\rm IA}$ is the redshift scaling.
Both $A_{\rm IA}$ and $\eta_{\rm IA}$ are free parameters, and will be
referred to as IA parameters. The convergence field with IA is then simply
$\kappa\rightarrow\kappa+\kappa_{\rm IA}$. This parameterization of IA is
called the non-linear linear alignment (NLA) model, named so because it is the
“linear”, or first-order, IA correction but uses the “non-linear” density
field. We take the first 100 simulations from the fiducial runs to include the
effects of IA, and these simulations will be used to take derivatives of our
data-vectors with respect to IA parameters. In specific, we create four
different variations with $\eta_{\rm IA}^{+}=-1.6$, $\eta_{\rm IA}^{-}=-1.8$
and $A_{\rm IA}^{+}=0.8$, $A_{\rm IA}^{-}=0.6$. These are variations of $0.1$
around the fiducial values of $A_{\rm IA}=0.7$ and $\eta_{\rm IA}=-1.7$ as
chosen in Krause et al. (2021, see their Table 2).444Note that Krause et al.
(2021) use a more sophisticated IA model, called TATT (see Section 5.2), which
has additional terms beyond the NLA model of Equation 2.16. We take the
fiducial values of their $a_{1}$ and $\eta_{1}$ parameters, which correspond
to $A_{\rm IA}$ and $\eta_{\rm IA}$ in Equation 2.16. Through these four runs,
we compute the derivatives of the data-vector with respect to the IA
parameters, and marginalize over these parameters in the analyses to follow.
Other, more sophisticated parameterizations of the IA effect also exist, and
we discuss our IA modeling approach in Section 5.2.

Shape noise. Once the two shear fields are generated, we add the relevant
shape noise in real space. In most cases, the forward-modelled field includes
Gaussian shape noise with a standard deviation given as

$\sigma_{\gamma}=\frac{\sigma_{e}}{\sqrt{n_{\rm gal}A_{\rm pix}}},$ (2.17)

where $n_{\rm gal}$ is the source galaxy number density, and $A_{\rm pix}$ is
the pixel area for a given map resolution. All maps in this work use
$\texttt{NSIDE}=1024$, corresponding to a pixel resolution of
$3.2\hbox{${}^{\prime}$}$. The per-galaxy shape noise is taken to be
$\sigma_{e}=0.26$. This technique is utilized for the DES Y6 and LSST Y1/Y10
fields. For DES Y3 we use a different technique given the weak lensing
catalogs are already available. In this case, we use the public galaxy shape
catalog555https://des.ncsa.illinois.edu/releases/y3a2/Y3key-catalogs and
rotate all galaxy ellipticities to remove any spatial correlation in the
shapes. The rotated shape catalog is used to make a shear field, where each
pixel value is the weighted average of the galaxy ellipticities in that pixel.
This is the same technique used by other works to estimate the DES Y3 shape
noise field (Gatti et al. 2022, 2023; Anbajagane et al. 2023a).

Survey Mask & Mass map construction. The noisy shear field — which is the sum
of the true shear fields and the shape noise fields — is then masked according
to the survey footprint. For DES Y3 and Y6 we use the provided survey mask in
the Y3 data release (Sevilla-Noarbe et al. 2021). For LSST Y1/Y10, we divide
the sky into three equal-area cutouts of roughly 14,000 deg2 each, which is
the expected area coverage for Y10 and $\approx 15\%$ larger than the expected
area for Y1. The noisy shear maps are converted to convergence using Equation
2.15. We only use the resulting E-mode field, $\kappa_{E}$, for our analyses.
The B-mode field is non-zero in this case as the presence of a mask transfers
some fraction of power (and thus cosmological information) from E-modes to
B-modes. The loss of Fisher information due to this leakage can be tested by
including the B-mode maps in the analysis. This inclusion will double the
data-vector size — assuming we extend the data-vector only with replications
of all measurements on the B-mode field and do not also consider cross-
correlations between the E-mode and B-mode fields — which will lead to poorer
numerical convergence of the final constraints. Thus, we do not test the loss
in Fisher information due to this leakage.

To summarize, lensing convergence maps are constructed from the raw particle
number count maps. The $n(z)$ distributions for a given survey are used to
obtain the convergence map in a given tomographic bin. IA effects are added if
desired. This convergence map is converted to shear maps, the relevant shape
noise is added, the relevant survey mask is applied, and then the noisy shear
maps are converted back to a noisy convergence map. The set of procedures
listed above is the standard approach for forward modeling the lensing field
(e.g., Zürcher et al. 2021; Zürcher et al. 2022; Gatti et al. 2022; Anbajagane
et al. 2023a). Thus, our final convergence maps will be an accurate
representation of the survey data.

Some known effects have been left out of our forward modeling procedure above:
mean redshift uncertainties (e.g., Myles et al. 2021), multiplicative bias
(e.g., MacCrann et al. 2022), clustering of source galaxies (e.g., Krause et
al. 2021; Gatti et al. 2020, 2023), and reduced shear (e.g., Krause & Hirata
2010; Gatti et al. 2020), to name a few. This choice has been made for
simplicity, and because these factors are not expected to change the Fisher
information constraints by a notable amount; either because the effect does
not include a nuisance parameter to marginalize over (e.g., reduced shear) or
is an effect with accurate enough calibration such that the nuisance parameter
marginalization has a negligible effect on the final constraints (e.g., mean
redshift, multiplicative bias). Nevertheless, we discuss the effects of these
components in more detail in Section 5.2.

##  3 Higher-order statistics

As mentioned prior, if a cosmological field can be defined by the covariance
between different pixels, then the field’s only degree of freedom is its power
spectrum or 2-point correlation function. Thus the latter are the “optimal”,
i.e. lowest noise, estimators to capture this information and maximize
constraints on any physics that imprints onto this field. In this work, the
PNG signal of interest is inherently non-Gaussian. Therefore, our chosen
summary statistic must extend beyond the 2-point function to capture the
relevant information. This information is often denoted “non-Gaussian” or
“higher-order” information.

In this section, we describe the different summary statistics employed in this
work to capture the non-Gaussian component of the field. These include (i) the
moments, which have been utilized for cosmological constraints in DES Y3
(Gatti et al. 2020, 2023), (ii) the CDFs, which have been tested for DES Y3
data (Anbajagane et al. 2023a), and (iii) the three-point function, which
contains the shape/configuration information whereas the former two only
contain the angle-averaged component (see Section 3.3 for more details). There
exist many more choices for a summary statistic sensitive to non-Gaussian
information, such as the wavelet phase harmonics, scattering transforms, mass-
aperture moments, integrated 3-point functions, skew-spectrum, and more (Allys
et al. 2020; Cheng & Ménard 2021; Halder et al. 2021; Secco et al. 2022b;
Munshi et al. 2022). We choose the moments as they have been applied to data
to extract cosmology, and choose the CDFs (which have been tested on data) and
3-point function as they are theoretically motivated extensions to the moments
approach.

###  3.1 Moments

The moments of the field provide an efficient way to summarize the information
from different orders. They are computed as

$\langle\kappa^{(1)}\kappa^{(2)}\ldots\kappa^{(N)}\rangle(\theta)=\frac{1}{N_{\rm
pix}-1}\sum_{i=1}^{N_{\rm
pix}}\kappa^{(1)}_{i}\kappa^{(2)}_{i}\ldots\kappa^{(N)}_{i}\,,$ (3.1)

where $\kappa^{(j)}$ is the convergence field of the $j^{\rm th}$ tomographic
bin, $N_{\rm pix}$ is the number of pixels in the survey footprint. In all
cases, $\langle\kappa^{(j)}\rangle=0$ is enforced and the scale dependence on
$\theta$ is obtained by making measurements after smoothing the fields with a
harmonic-space tophat filter,

$W_{\ell}(\theta)=2\frac{j_{1}(\ell\theta)}{\ell\theta}\,,$ (3.2)

where $\theta$ is the tophat radius in radians. We use ten bins of $\theta$,
between $3.2\hbox{${}^{\prime}$}<\theta<200\hbox{${}^{\prime}$}$, which
follows the analysis choices of Gatti et al. (2022). All fields are smoothed
by the same $\theta$. Varying this choice so each field in the triplet has a
different smoothing scale can probe additional information in the field. We
prioritize using the same choices as Gatti et al. (2022), who used a single
$\theta$. The quantity $\kappa_{i}$ in Equation 3.1 is the same as the
$\kappa(\mathbf{\hat{n}})$ defined in equation 2.1, but with the continuous
direction $\mathbf{\hat{n}}$ now replaced by a discrete one defined by pixels
in the Healpix map.

The Nth moment is sensitive to information from the N-point correlation
function; the 2nd moment depends on the 2-point function, while the third
moment depends on the 3-point function. In specific, the moments depend on the
volume-integrated N-point functions of the same order. This means that any
dependence of the correlation on the specific shape of the N-point function is
integrated over when measuring the moments. In this work, we use the 2nd, 3rd,
4th, and 5th moments. Anbajagane et al. (2023a, see their Figure 6) measure
this set of moments on DES Y3 data and find the 5th moment is consistent with
no cosmological signal. However, we include it in this work as the improved
precision of an LSST Y10-like dataset could enable the extraction of
cosmological signals at this order. Including the sixth moment will double the
length of the existing data-vector in return for marginal improvements in the
Fisher information, and so we do not consider it.

The second and third moments of the field have been validated extensively to
determine their robustness as a summary statistic (Gatti et al. 2020) and this
has led to their use in DES Y3 to constrain cosmology parameters (Gatti et al.
2022). This is in contrast to the other two statistics we consider here, which
are not yet validated at the rigor required for extracting lensing-based
constraints.

###  3.2 Cumulative Distribution Function (CDF)

The CDFs are a statistic closely connected to the moments as they are both
different representations of the same distribution of lensing convergence,
$P(\kappa)$. The CDF measurements are defined as,

$\displaystyle{\rm CDF}(\theta,\nu)$
$\displaystyle=P(\kappa^{(1)}>\nu,\kappa^{(2)}>\nu,\ldots\kappa^{(N)}>\nu)$
$\displaystyle=\frac{1}{N_{\rm pix}}\sum_{i=1}^{N_{\rm
pix}}\Theta(\kappa^{(1)}_{i}-\nu)\Theta(\kappa^{(2)}_{i}-\nu)\ldots\Theta(\kappa^{(N)}_{i}-\nu)\,,$
(3.3)

where $\Theta(\kappa-\nu)$ is the Heaviside step function, which takes values
of 1 if $\kappa-\nu>0$ and 0 otherwise, and $\langle\kappa^{(j)}\rangle=0$ is
enforced similarly to the moments measurement. The quantity ${\rm
CDF}(\theta,\nu)$ captures what fraction of the field, smoothed on a scale
$\theta$, is above a threshold $\nu$. The $\nu$ have the same units as
$\kappa$ and are dimensionless quantities. The choice of smoothing scales is
the same as that of the moments,
$3.2\hbox{${}^{\prime}$}<\theta<200\hbox{${}^{\prime}$}$, and for each scale,
we use 7 thresholds $\nu\in\\{-20,-6,-2,0,20,6,20\\}\times 10^{-3}$.
Anbajagane et al. (2023a) determined these thresholds through an approximate
optimization procedure for DES Y3 data. We employ the same for simplicity and
do not explore survey-specific thresholds. While we refer to these
measurements as the CDFs, they are formally the CDF “complement” as the
traditional CDFs are defined as the fraction of a distribution below a certain
threshold, and not above the threshold as we do here. Once again the scale
dependence of the measurement enters through smoothing the field,
$\kappa^{(j)}$, with a tophat filter; the same filter as defined in equation
3.2. We will consider measurements of up to the 3-field CDFs, so
$N\in\\{1,2,3\\}$.

As mentioned prior, the CDFs are intimately connected to the moments. In
particular, for a given threshold $\nu$, the measurement in Equation 3.2 is
sensitive to moments of all orders. In the limit where the CDFs are computed
with arbitrarily high number of thresholds, and the moments are computed to
arbirarily high orders, the two constraints will match.666This is not formally
true for every field as certain distributions — with the log-normal being the
most well-known — cannot be represented using a finite combination of their
moments. In the practical limit with noise, however, this matching is possible
for every field. Banerjee & Abel (2023) formally derive that the CDFs contain
volume-integrals of all N-point correlation functions. As a consequence, the
CDFs do not contain any shape/configuration information, as is the case with
the moments. Naively, this would suggest the moments and CDFs cannot
distinguish contributions from the different $f_{\rm NL}$ types. However, we
will show in Section 4.3 that each $f_{\rm NL}$ has a different scale- and
redshift-dependence that enables distinguishing between them even without the
shape/configuration information.

Anbajagane et al. (2023a) computed the impact of various lensing-based
systematics on the CDF data-vector for DES Y3 data quality, and identified the
relevant/negligible systematic effects. However, a theoretical model of the
CDFs and an end-to-end validation of a CDF inference pipeline are required
before this statistic can be used on data.

###  3.3 3-point correlation function

Both the moments and the CDFs probe information to higher orders but only do
so for the angle-averaged correlations; they contain volume integrals of the
N-point functions rather than the functions themselves. For example, if the
field contains a 3-point function that is non-zero only when the three points
are equidistant (i.e. they form an equilateral triangle) then the moments and
CDFs measurements detailed above would be unable to distinguish this from a
field where a different type of triangle is the sole contributor. As an
alternative, we consider the full 3-point function which includes all of this
shape/configuration dependence. Inflationary signatures are often constrained
using the galaxy bispectrum (Cabass et al. 2022b; D’Amico et al. 2022; Philcox
et al. 2022a) which contains all of this shape information and is particularly
useful in distinguishing between the different types of $f_{\rm NL}$. Since
lensing convergence is a projected integral of the density field, the shape
information will be less useful in distinguishing these types but is still
expected to provide additional information beyond the volume-integrated
statistics considered above.

Estimating the full 3-point function has a naive computational complexity of
$\mathcal{O}(N^{3})$, where $N$ is the number of points being correlated. In
our work, these points correspond to convergence map pixels. Tree-based
calculations, such as TreeCorr (Jarvis et al. 2004), achieve a computational
complexity of $\mathcal{O}(N^{2}\log_{2}N)$. While this is a significant
speedup, it is not adequate for our requirements as we compute this statistic
on $\mathcal{O}(10^{4})$ survey realizations and for 35 (20) triplets between
the tomographic bins of LSST (DES) per realization. Thus, computational speed
is a necessity in allowing a statistic to be used in the simulation-based
model.

An alternative approach explored by Philcox et al. (2022b); Sunseri et al.
(2023) is a Fast Fourier Transform (FFT)-based calculation of the 3-point
function with a complexity of $\mathcal{O}(N\log_{2}N)$, which is the same as
that of tree-based 2-point function estimators. We reproduce below the main
aspects of the computational procedure for completeness and direct readers to
Sunseri et al. (2023, see their Section 2.2) for a detailed derivation.

A 3-point correlation function of a scalar field, henceforth denoted $\zeta$,
is computed by counting triangles formed by different triplets of points.
Thus, the function can be parametrized by the length of two sides, $r_{1}$ and
$r_{2}$, and the angle between them $\phi$. A key choice in Sunseri et al.
(2023) is splitting the radial and angular dependence of $\zeta$ as

$\zeta(r_{1},r_{2},\hat{r}_{1}\cdot\hat{r}_{2})=\sum_{m=0}^{\infty}\zeta_{m}(r_{1},r_{2})e^{im\phi}.$
(3.4)

Note that this expansion pertains to any projected, 2D scalar fields, which in
our case is the convergence field, $\kappa(\vec{x})$. Thus, $\zeta$ is a
projected 3-point function, and $r_{i}$ are projected separations. The radial
coefficients, $\zeta_{m}(r_{1},r_{2})$, can be computed as an area average
over the convergence field $\kappa(\vec{x})$,

$\zeta_{m}(r_{1},r_{2})=\frac{1}{2\pi^{2}}\int\frac{d^{2}\vec{x}}{A}\,\kappa(\vec{x})\,c_{m}(r_{1},\vec{x})\,c^{*}_{m}(r_{2},\vec{x})\,,$
(3.5)

where the Fourier coefficients $c_{m}$ are

$\displaystyle c_{m}(r_{i};\vec{x})\equiv\int
d\phi\;e^{-im\phi}\kappa(\vec{x}+\vec{r}_{i}).$ (3.6)

Thus far, Equations 3.4, 3.5, and 3.6 are written as a function of distances
$r_{1}$ and $r_{2}$. However, our final calculation will involve $\zeta$
computed in different radial and angular bins. We can thereby bin the Fourier
coefficients in circular annuli as,

$c_{m}(r_{i},\vec{x})\to c_{m}^{\rm b}(\vec{x})=\int_{r_{\rm min}^{b}}^{r_{\rm
max}^{b}}\frac{dr_{i}}{A^{b}}r_{i}\int
d\phi\,\,\kappa(\vec{x}+\vec{r}_{i})e^{-im\phi},$ (3.7)

where $A^{b}=\pi[(r_{\rm max}^{b})^{2}-(r_{\rm min}^{b})^{2}]$ is the area
corresponding to the annular bin, and $r_{\rm min}^{b}$ and $r_{\rm max}^{b}$
are the minimum and maximum radius of bin $b$. With this binning, we write the
$\zeta_{m}$ coefficients as

$\displaystyle\zeta_{m}(r_{1},r_{2})\rightarrow\zeta_{m}^{\mathrm{b_{1}}\mathrm{b_{2}}}=\frac{1}{2\pi^{2}}\int\frac{d^{2}\vec{x}}{A}\,\kappa(\vec{x})\,c^{\rm
b_{1}}_{m}(\vec{x})\,c^{\rm b_{2}\,*}_{m}(\vec{x}),$ (3.8)

where $\rm b_{1}$ and $\rm b_{2}$ are bin indices. Equation 3.4 shows that the
3-point function $\zeta$ can be constructed by summing over the coefficients
$\zeta_{m}$ with some angular phase,

$\displaystyle\zeta^{\rm
b_{1},b_{2}}(\phi)=\sum_{m=0}^{\infty}\zeta_{m}^{\mathrm{b_{1}}\mathrm{b_{2}}}e^{im\phi}.$
(3.9)

Equation 3.9 is the final 3-point correlation function used in this work. The
product within the sum is still a complex number, but we only keep the real
component given the 3-point function, $\zeta$, describes correlations in real-
space and has no imaginary component. When building a three-point cross-
correlation between tomographic bins, we still compute Equation 3.8 and choose
which triplet of bins provides the fields $\kappa(\vec{x})$, $\,c^{\rm
b_{1}}_{m}(\vec{x})$ and $c^{\rm b_{2}\,*}_{m}(\vec{x})$. If all fields come
from the same bin then $\zeta$ is the auto-correlation, and if at least one
field comes from a different bin then it is a cross-correlation.

The formalism above utilizes FFTs to compute the coefficients in Equation 3.7.
This inherently assumes the field can be approximated as a flat field. The
wide-field surveys we consider in this work, however, cannot be treated in
such a manner and require spherical harmonics, which account for the curvature
of the maps across the sky. To resolve this difference, we split the survey
footprint into a set of smaller patches — for which the flat-field
approximation is adequate — and compute $c^{\rm b}_{m}(\vec{x})$ separately
for each patch. For each HEALPix pixel, we compute the quantity
$\kappa(\vec{x})\,c^{\rm b_{1}}_{m}(\vec{x})\,c^{\rm b_{2}\;*}_{m}(\vec{x})$,
as defined within the integral of Equation 3.8, and average it across the
survey footprint to obtain $\zeta_{m}^{\rm b_{1},b_{2}}$.

In our implementation, the patches follow a resolution of $\texttt{NSIDE}=4$
which corresponds to a size $15\times 15\deg^{2}$. In practice, each flat
field also includes an additional buffer region (of $8\deg$) around the four
sides, which alleviates edge effects during FFT calculations. We note,
however, that any artifacts induced by the specific choices above (e.g., patch
size, buffer width) do not induce biases in the final Fisher information since
these choices are applied consistently across all measurements in the
analysis. This is also true of future applications to data: as long as the
exact same computational procedure is performed on data — and the simulations
processed to look like data (where the latter is done to build a simulation-
based model) — the inference of any parameters will be unbiased.

The other required choices in measuring this 3-point function are the
tomographic bin combinations for which we compute $\zeta$, as well as the
maximum $m$ mode used in computing Equation 3.9. For the former, we choose all
triplet combinations of the tomographic bins, and this choice generates 20
(35) different combinations for DES (LSST). The maximum $m$ mode is $m_{\rm
max}=5$, which is the same choice made in Sunseri et al. (2023). We do not
explore higher $m_{\rm max}$ given computing limitations. We compute $r_{i}$
in 6 logarithmic bins between
$3.2\hbox{${}^{\prime}$}<r_{i}<200\hbox{${}^{\prime}$}$, and $\phi$ in 6
linear bins between $0<\phi<\pi$. The choice of 6 bins in the former is set by
computational cost. Increasing the number of bins extends the total compute
time as $N_{\rm bin}^{2}$. The choice of 6 $\phi$ bins is because $\zeta$ is
computed using 6 $m$-modes. The quantities $\zeta_{m}$ and $\zeta$ are related
via Fourier modes, and in principle, a fine sampling in real-space is required
to capture the Fourier coefficients (and vice-versa). Thus, it may be
advantageous to continue with the $m$-mode coefficients without converting to
$\phi$ bins. However, we perform the conversion so the final quantity is a
real-space measure, the 3-point function, that is directly related to other
existing measurements of the 3-point function (e.g., Secco et al. 2022b). In
addition, the conversion reduces the data-vector’s memory footprint as
$\zeta_{m}$ are complex coefficients, requiring a number for each of the real
and complex parts, while $\zeta$ can only be real.

This 3-point estimator has not yet been applied to the weak lensing field.
Thus an application of this statistic to data will require further validation.
In principle, the $\zeta$ above is closely related to the 3-point shear
correlation functions studied and tested in Secco et al. (2022b) for DES Y3
data. However, in practice, the computational procedures of the two are
significantly different — for example, Secco et al. (2022b) measure the shear
3-point functions using galaxy shape catalogs, whereas we measure the
convergence 3-point functions using pixelized maps — and so an explicit
validation must still be performed.

The total compute time for this 3-point estimator is roughly 10 times that of
the moments (when computing 2nd to 5th moments). Given these computational
demands, we limit our use of the 3-point function to one specific comparison
test (see Section 4 below) aimed at determining the Fisher information in the
configuration/shape component of the 3-point function. We do not use this
statistic in our fiducial constraints.

##  4 $f_{\rm NL}$ constraints from weak lensing

We now present the Fisher information on PNGs as probed by the weak lensing
fields. In §4.1, we determine the optimal summary statistic between the ones
described in §3, and in §4.2, show the results for different survey datasets
(both current and upcoming). In §4.3, we isolate the signatures of inflation
as seen in the lensing field, and identify the physical processes involved in
generating such signatures. In all results below, the Fisher information is
estimated as

$\textbf{F}_{ij}=\sum_{m,n}\frac{d\widetilde{X}_{m}}{d\theta_{i}}\big{(}\mathcal{C}^{-1}\big{)}_{mn}\frac{d\widetilde{X}_{n}}{d\theta_{j}},$
(4.1)

where $\frac{d\widetilde{X}_{m}}{d\theta_{i}}$ is the mean derivative of point
$m$ in data-vector $X$ with respect to parameter $\theta_{i}$, where the mean
is computed using 300 to 400 independent survey realizations depending on the
survey. $\mathcal{C}^{-1}$ is the inverse of the numerically estimated
covariance matrix and is computed while accounting for the Hartlap factor
(Hartlap et al. 2007),

$\mathcal{C}^{-1}\rightarrow\frac{N_{\rm sims}-N_{\rm data}-2}{N_{\rm
sims}-1}\,\mathcal{C}^{-1}.$ (4.2)

We verify in Section B that the covariance is well-converged for all summary
statistics. The Hartlap factor is $\gtrsim 0.95$ for the 2nd and 3rd moments,
the fiducial data-vector used in this work, and is $\gtrsim 0.55$ for the full
3-point function. Note that for the analysis using the 3-point function, we
generate additional “pseudo-independent” realizations. These additional
realizations are made by randomly rotating the original convergence maps, and
adding a completely independent noise realization to them. For every
independent realization, we make roughly two pseudo-independent ones and have
a total of $16,000$ realizations. This increase in the number of realizations
is required given the large length of the 3-point data-vector (when using no
scale cuts). The set of $16,000$ realizations is used to compute the
covariance of all data-vectors only in the analysis comparing the
configuration/shape information (Figure 3). We emphasize that all other
results in this work do not use any pseudo-independent realizations and only
work with fully independent ones.

In all presentations of the Fisher information, we show both the raw estimate
as well as one degraded by 40%. The latter is used as a potentially
pessimistic estimate of the Fisher information, and the specific choice of 40%
is because this degradation would correspond to a survey with half the
expected survey volume. As a result, all Fisher constraints below are shown as
bands, with the lower (upper) limit corresponding to the raw (pessimistic)
estimate.

###  4.1 Dependence on summary statistic

Figure 2: The Fisher information, as a function of the minimum angular scale
of the data-vector, for different statistics measured on an LSST Y10-like
survey. The lower (upper) bound is the Fisher (pessimistic Fisher)
information, where the degradation factor of 40% for the pessimistic case
approximates constraints from half the expected survey volume. The columns
correspond to different $f_{\rm NL}$, and the rows progressively step from
constraints varying just $f_{\rm NL}$, to marginalizing over $\sigma_{8}$ and
$\Omega_{\rm m}$, and finally to also marginalizing over the intrinsic
alignment parameters, $\eta_{\rm IA}$ and $A_{\rm IA}$. The constraints for
$f_{\rm NL}^{\rm\,loc}$ and $f_{\rm NL}^{\rm\,eq}$, when varying just $f_{\rm
NL}$, are consistent across all statistics. However, once other parameters are
included in the analysis, the statistics sensitive to non-Gaussianities are
clearly more favorable. The combination of the 2nd and 3rd moments does fairly
similarly to all other non-Gaussian statistics considered, with the exception
of analyses of $f_{\rm NL}^{\rm\,or,\,cmb}$. The 2nd moment constraints are
above the range of many panels. Figure 3: Similar to Figure 2 but for the
moments and the full 3-point correlation function. The purple bands in most
panels overlap completely with the yellow bands and are not seen. The 3-point
function constraints are considerably better than the 3rd moment constraints,
indicating the configuration information found in the former (and missing in
the latter) is significant. Note that the constraints for the moments are
weaker here than those in Figure 2 as all summary statistics are measured
within 6 radial bins rather than the fiducial choice of 10. This is motivated
by computing limitations for the 3-point function, and is discussed further in
the text. The 3rd moment constraints are above the range of many panels.

Figure 2 shows the Fisher information in different summary statistics for an
LSST Y10 survey. The first row shows constraints when varying just $f_{\rm
NL}$ (with y-labels denoting the specific type being varied). In the second
row we marginalize over $\sigma_{8}$ and $\Omega_{\rm m}$, and in the third
row we also marginalize over the IA parameters. All results are shown as a
function of minimum angular scale of the data-vector, and all discussions on
scale-dependence are in Section 4.3. For $f_{\rm NL}^{\rm\,eq}$ and $f_{\rm
NL}^{\rm\,loc}$, in the unmarginalized case, all statistics considered here
are roughly equivalent. In particular, the 2nd moments alone are an adequate
statistic, and including the higher order moments adds only marginally to the
Fisher information. On non-linear scales, the late-time matter power spectra
contain strong signatures from these two types of $f_{\rm NL}$ (Coulton et al.
2022, see Figure 6) and can thus constrain these parameters on their own.
Higher-order spectra, such as the matter bispectrum and trispectrum, will
contain additional information but are also more noise-dominated than the
power spectra. Hence, their contribution to the total Fisher information can
be minuscule compared to the contribution from the power spectra. This
behavior is seen in Figure 2, where the 2nd moments dominate the constraining
power over the 3rd, 4th and 5th moments. For both orthogonal-type $f_{\rm
NL}$, the constraints from the 2nd moments alone are significantly weaker than
combinations with higher-order moments. This is expected as Coulton et al.
(2022, see Figure 1) shows the non-linear power spectrum has only minimal
(negligible) changes when varying $f_{\rm NL}^{\rm\,or,\,lss}$ ($f_{\rm
NL}^{\rm\,or,\,cmb}$).

The relevance of the higher-order information improves significantly when
extending the parameter space to perform a marginalized analysis of $f_{\rm
NL}$. The changes in the power spectrum due to PNGs have some overlap with
changes due to gravitational evolution. Thus, marginalizing over cosmology
(which is marginalizing over gravitational evolution) results in a strong
degradation in the Fisher information of PNGs as probed by the 2nd moments.
Including the 3rd moment significantly improves the constraints. The
degeneracy-breaking from combining 2nd-order and 3rd-order information has
been explored extensively in the literature (e.g., Gatti et al. 2020; Zürcher
et al. 2021; Gatti et al. 2022; Zürcher et al. 2022; Anbajagane et al. 2023a).
They have also been explicitly shown for the PNGs we study here (Coulton et
al. 2022). Including the 4th and 5th moments improves the constraints
slightly. The CDFs are generally similar to the combinations of 2nd and 3rd
moments, and are better/worse depending on the exact analysis being performed:
they are better in constraining $f_{\rm NL}^{\rm\,eq}$, worse for $f_{\rm
NL}^{\rm\,loc}$, and comparable for $f_{\rm NL}^{\rm\,or,\,lss}$ and $f_{\rm
NL}^{\rm\,eq}$. This is generally consistent with the behaviors of the moments
and CDFs found in the $w{\rm CDM}$ analysis of Anbajagane et al. (2023a).
Further marginalizing over IA parameters results in minimal degradation of the
$f_{\rm NL}$ constraints.

Importance of configuration/shape information. All the statistics above have
utilized the angle-averaged information in the field. These statistics involve
volume integrals over the correlation functions of the field and have no
sensitivity to the shape of the correlations. In Section 3.3, we discussed the
potential information in the configuration/shape information of the field’s
spatial correlations. Here, we explicitly check this by computing the full
3-point correlation function and comparing its constraints to those from the
3rd moments, which have no shape information. Due to computing limitations, we
only measure the 3-point functions in 6 radial bins between
$3.2\hbox{${}^{\prime}$}<\theta<200\hbox{${}^{\prime}$}$, as opposed to the 10
bins used in the main analysis.777The calculation complexity goes as $N_{\rm
r,bin}^{2}$ so using 10 radial bins (instead of 6) extends the compute time by
a factor of 3. To perform a fair comparison between statistics, we also
remeasure the moments in 6 radial bins as well. However, we will show that
this reduction in bins reduces the Fisher information on PNGs as probed by the
moments. Therefore, we use the measurements of moments with 6 radial bins only
for the comparison in this section and use the fiducial measurements with 10
bins for all other analyses. The moments here continue to be measured directly
on the spherical sky, and do not use the patch-by-patch flat-sky approach used
for the 3-point function estimator.

Figure 3 shows the constraints from the 3rd moments and the 3-point function,
as well as from combinations with the 2nd moments and from the combination of
all three. The Fisher information for all $f_{\rm NL}$ types is higher in the
full 3-point function in comparison to the 3rd moment. This increase in
information is still notable when varying only $f_{\rm NL}$, and increases to
$50\%$ improvements when marginalizing over cosmology and/or IA. The exception
is $f_{\rm NL}^{\rm\,or,\,cmb}$ where the constraints are improved by a factor
of 2 or more in all settings. The combination of 2nd moment and 3-point
function does better than the latter alone. Adding the 3rd moment to this
combination leads to no improvement; the purple and yellow bands are atop each
other for most panels. This emphasizes that the 3-point function contains all
the information in the 3rd moments, as expected. This behavior is consistent
regardless of the set of parameters being varied in the analysis.

We have verified that the marginalized $f_{\rm NL}$ constraints from using the
3-point function (either on its own or combined with other statistics) are
numerically converged to within 10-20%, indicating that the constraints in
Figure 3 are potentially overestimated by 10-20%. This non-convergence arises
due to noise in the numerically estimated derivative. It is not related to the
covariance, as we have verified the constraints change by $<1\%$ when changing
the number of realizations used in estimating the covariance. The numerical
convergence-based overestimate of 10-20% is lower than the 50% improvement
mentioned above. Therefore, it is still likely that the configuration/shape
measurement leads to an increase in the Fisher information of PNGs. A robust
estimate of this increase will require better numerical convergence in the
estimates of the derivatives.

The poorer numerical convergence of constraints from the 3-point function,
compared to the percent-to-sub-percent convergence of constraints from the
CDFs and the moments, can be improved by employing more independent
realizations to estimate the derivatives. Note that the convergence issues
have occured even after we consciously reduced the size of the data-vector by
limiting the radial binning to 6 bins instead of the fiducial 10 bins.
Comparisons of the moments-based constraints in Figure 2 and Figure 3 also
show that this change in binning leads to significantly degraded constraints.
This highlights the practical challenges in robust use of the three-point
function while performing simulation-based modeling. One could compress the
datavector to reduce its size and thus, the number of simulations needed to
estimate the covariance. Studies on data have used SVD decompositions (e.g.,
Zürcher et al. 2022) or the MOPED compression (e.g., Gatti et al. 2022). Given
our goal of performing a true comparison between the 3-point function and the
3rd moments, we do not employ any compression. Though, we note the MOPED
compression, in principle, should not degrade the Fisher constraints. In our
case, this is complicated by the fact that the MOPED compression depends on
our simulation-based derivatives of the data-vectors, and the noise in these
derivative estimates will lead to poorer compression. We do not explore this
further.

Given the discussions above on the benefits and detriments in practical
implementations of each statistic, we henceforth use the combination of the
2nd and 3rd moments as the fiducial statistic for this work and for all the
Fisher information analyses below.

###  4.2 Fisher information in DES and LSST

Figure 4: The Fisher information in the 2nd and 3rd moments of the lensing field, shown as a function of the minimum angular scale of data-vector, for four different survey configurations. The lower (upper) bound is the Fisher (pessimistic Fisher) information, where the degradation factor of 40% for the pessimistic case approximates constraints from half the expected survey volume. The existing constraints from BOSS (D’Amico et al. 2022) are shown as the black line and the potential constraint from DESI (D’Amico et al. 2022; DESI Collaboration et al. 2016) in red. The presented lensing and galaxy constraints are both consistent in fixing cosmological parameters. The DES Y6 constraints for $f_{\rm NL}^{\rm\,eq}$ are comparable to the existing BOSS constraints. The LSST Y10 constraints are comparable or potentially better than DESI for $f_{\rm NL}^{\rm\,eq}$, a factor of 2 broader for $f_{\rm NL}^{\rm\,or,\,lss}$, and a factor of 8 broader for $f_{\rm NL}^{\rm\,loc}$. The DES Y3 constraints are above the range of the plots for most panels. Survey | $\sigma(f_{\rm NL}^{\rm\,eq})$ | $\sigma(f_{\rm NL}^{\rm\,loc})$ | $\sigma(f_{\rm NL}^{\rm\,or,\,lss})$ | $\sigma(f_{\rm NL}^{\rm\,or,\,cmb})$
---|---|---|---|---
Fiducial Fisher information, $\theta_{\rm min}>3.2\hbox{${}^{\prime}$}$
DES Y3 | 334 [984] | 125 [315] | 562 [850] | 1136 [1456]
DES Y6 | 187 [575] | 70 [186] | 300 [524] | 648 [905]
LSST Y1 | 136 [295] | 50 [90] | 169 [234] | 337 [457]
LSST Y10 | 109 [186] | 39 [55] | 116 [142] | 220 [278]
$\theta_{\rm min}>20\hbox{${}^{\prime}$}$
DES Y3 | 530 [1279] | 203 [397] | 745[946] | 1451 [1809]
DES Y6 | 355 [833] | 137 [262] | 471 [618] | 993 [1355]
LSST Y1 | 200 [437] | 75 [124] | 246 [280] | 463 [555]
LSST Y10 | 168 [338] | 63 [92] | 187 [195] | 343 [381]
Galaxy correlation functions
BOSS (fix) | 212 | 29 | 72 | —
BOSS (vary) | 350 | — | — | —
DESI (fix) | 133 | — | 45 | —
DESI (vary) | 220 | 5 | — | —
Table 3: The Fisher information presented in Figure 4 for the full range of
scales (top) and a conservative scale cut (middle). The numbers in square
brackets are constraints after marginalizing over $\Omega_{\rm m}$,
$\sigma_{8}$, and two IA parameters. The “BOSS (fix)” results, analyzed at
fixed cosmology, are from D’Amico et al. (2022) while the “BOSS (vary)”
results from also varying $\sigma_{8}$ are from Philcox et al. (2022a). The
expected constraints from DESI for $f_{\rm NL}^{\rm\,eq}$ and $f_{\rm
NL}^{\rm\,or,\,lss}$ are obtained by rescaling the corresponding BOSS results
by 1.6 as mentioned in D’Amico et al. (2022), and for $f_{\rm NL}^{\rm\,loc}$
we take the numbers from DESI Collaboration et al. (2016). Lensing
measurements at $\theta_{\rm min}=3.2\hbox{${}^{\prime}$}$ ($\theta_{\rm
min}=20\hbox{${}^{\prime}$}$) correspond to power spectra at
$k\in[1,2.5]\,\,h/{\rm Mpc}$ ($k\in[0.3,0.6]\,\,h/{\rm Mpc}$), depending on
the redshift bin. These are higher than the current $k_{\rm max}=0.2h/{\rm
Mpc}$ limit of galaxy correlation analyses. The difference is more substantial
if we consider the maximum scale (and not the average scale) that contributes
to the lensing measurement. These scale differences are discussed further in
Section 5.1.

The key goal of this work is extracting the Fisher information on $f_{\rm NL}$
from weak lensing, and comparing it to the information in galaxy correlation
functions from Baryon Oscillation Spectroscopic Survey (BOSS) or expected from
DESI. Such comparisons motivate/inform the utility of weak lensing in PNG
analyses. We study two versions of the DES and LSST surveys, and compare their
constraints with current constraints from BOSS and expected ones from DESI.
The BOSS constraints come from the analysis of D’Amico et al. (2022), while
the expected DESI numbers are taken from (i) D’Amico et al. (2022) for $f_{\rm
NL}^{\rm\,eq}$ and $f_{\rm NL}^{\rm\,or,\,lss}$, as they compute the expected
constraints to be $1.6$ times higher than their BOSS constraints, and (ii)
DESI Collaboration et al. (2016) for $f_{\rm NL}^{\rm\,loc}$. Note that other
works (Cabass et al. 2022a; Philcox et al. 2022a) show significantly different
results from D’Amico et al. (2022), primarily due to the choice of which
parameters to marginalize over; the latter fixes cosmology while the former
vary cosmology and other nuisance parameters. We will present results for
constraint with and without such marginalization (Table 3).

Figure 4 presents the Fisher information for weak lensing surveys as a
function of the minimum angular scale used in the analysis. We detail our
findings below for each of the four different $f_{\rm NL}$ explored:

For $\boldsymbol{f_{\rm NL}^{\rm\,loc}}$, galaxy correlations have
significantly more information than weak lensing. This is expected given the
signature of $f_{\rm NL}^{\rm\,loc}$ in galaxy correlation functions is a
scale-dependent galaxy bias in the 2-point function, where the bias increases
towards large scales as $b\propto k^{-2}$ with $k$ being the Fourier
wavenumber (Dalal et al. 2008). Thus, this effect has a large (diverging)
amplitude towards large scales and is more easily distinguished from
gravitational evolution. The Fisher information in LSST (DES) is factors of 3
to 8 times lower than the information in DESI (BOSS). Thus, the benefit from
including weak lensing in analyses of $f_{\rm NL}^{\rm\,loc}$ is unlikely to
be from the weak lensing-only constraints and would instead be from
constraints of galaxy bias parameters via the cross-correlation of lensing and
galaxies. These galaxy bias parameters, numbering $\mathcal{O}(10)$ in the
latest models (D’Amico et al. 2022; Cabass et al. 2022b, a; Philcox et al.
2022a), are required in the modeling of the correlation functions and cannot
be known a priori. A more detailed discussion of this aspect can be found in
Section 5.1. It is also possible that the weak lensing-only constraints enable
degeneracy breaking that does benefit the final $f_{\rm NL}^{\rm\,loc}$
constraints. The analysis in this work does not compute the Fisher information
in the galaxy correlations and is unable to test this.

For $\boldsymbol{f_{\rm NL}^{\rm\,eq}}$, weak lensing provides constraints
that are competitive with galaxy clustering. The Fisher information in LSST
(DES) is similar to, and potentially better than, DESI (BOSS). Note that the
difference in constraining power between DES Y3 and DES Y6 is about 60%. The
Y3 analysis uses the actual Y3 noise fields and redshift
distribution/uncertainties, while the Y6 results use expected noise fields and
redshift distributions. The improvement in DES Y6 over DES Y3 is estimated to
be (i) a 50% increase in galaxy sample size, which implies a 25% increase in
constraints, and (ii) an increase in the highest redshifts measured, where the
amplitude of the lensing signal grows with redshift and this is expected to
cause significant improvement in $f_{\rm NL}$ constraints (see Figure 6). We
also find that a DES Y6 survey with the DES Y3 number density still performs
better than the DES Y3 fiducial survey (see Figure 5) which further signifies
the importance of the higher maximum redshift in DES Y6 compared to DES Y3.
These all signify that the 60% improvement is reasonable.

For $\boldsymbol{f_{\rm NL}^{\rm\,or,\,lss}}$ and $\boldsymbol{f_{\rm
NL}^{\rm\,or,\,cmb}}$, the constraints from weak lensing are about 1.5 to 4
times broader than those from galaxy clustering. We only compare $f_{\rm
NL}^{\rm\,or,\,lss}$ as the BOSS/DESI analyses do not measure $f_{\rm
NL}^{\rm\,or,\,cmb}$. The DES constraints for the former $f_{\rm NL}$ type are
more than a factor of 3 wider than the existing BOSS constraints, while the
LSST constraints are a factor of 1.5 to 2 wider than those of DESI. Therefore,
lensing can still provide valuable information in constraining these $f_{\rm
NL}$ amplitudes for Stage IV surveys.

The sound speed, $\boldsymbol{c_{s}}$, and the EFT parameter
$\boldsymbol{\widetilde{c}_{3}}$, capture the two leading-derivative cubic
interactions of the EFT of inflation (e.g., Cabass et al. 2022b, see their
Equation 1). These parameters can be directly inferred from $f_{\rm
NL}^{\rm\,eq}$ and $f_{\rm NL}^{\rm\,or,\,lss}$ as

$\displaystyle c_{s}$ $\displaystyle=\bigg{[}1+\frac{324}{85}\bigg{(}1.3f_{\rm
NL}^{\rm\,eq}+8.97f_{\rm NL}^{\rm\,or,\,lss}\bigg{)}\bigg{]}^{-1/2},$ (4.3)
$\displaystyle(c_{s}^{-2}-1)\widetilde{c}_{3}$
$\displaystyle=\frac{243}{10}\bigg{(}-0.293f_{\rm NL}^{\rm\,eq}-7.71f_{\rm
NL}^{\rm\,or,\,lss}\bigg{)}-\frac{3}{2}(c_{s}^{2}-1),$ (4.4)

where the expressions are taken from the conversions presented in Equation 5
of Cabass et al. (2022b). Since the constraints allow the value $c_{s}=1$ (and
thus $1/(c_{s}^{2}-1)$ is undefined) one can only constrain the combination
$(c_{s}^{-2}-1)\widetilde{c}_{3}$ rather than $\widetilde{c}_{3}$. Equations
4.3 and 4.4 require joint constraints on $f_{\rm NL}^{\rm\,eq}$ and $f_{\rm
NL}^{\rm\,or,\,lss}$, whereas we have thus far only varied one $f_{\rm NL}$ at
a time. Upon doing the joint analysis, we obtain the following $c_{s}$
constraints for DES Y6 and LSST Y10, depending on whether or not we
marginalize over cosmology and IA parameters,

$\displaystyle c_{s}^{\rm
DES\,Y6}\gtrsim\begin{cases}0.021,&\text{unmargin.}\\\
0.011,&\text{margin.}\\\ \end{cases}\hskip 20.0pt[\text{lower}\,\,95\%]$ (4.5)
$\displaystyle c_{s}^{\rm
LSST\,Y10}\gtrsim\begin{cases}0.028,&\text{unmargin.}\\\
0.016,&\text{margin.}\\\ \end{cases}\hskip 20.0pt[\text{lower}\,\,95\%]$ (4.6)

where we have followed Planck Collaboration et al. (2020a); Cabass et al.
(2022b) in marginalizing over $\widetilde{c}_{3}$ and presenting constraints
on just $c_{s}$. We stress that the results above are only rough estimates.
This is because our analysis varies $f_{\rm NL}^{\rm\,eq}$ and $f_{\rm
NL}^{\rm\,or,\,lss}$, and not $c_{s}$ and $\widetilde{c}_{3}$ directly. Not
all of the parameter space spanned by the former leads to well-defined
quantities in the latter (Cabass et al. 2022b). A more robust estimate would
be obtained by running simulations that vary $c_{s}$ and $\widetilde{c}_{3}$
directly. We do not have such simulations so do not do this. If we consider
the marginalized case, bounds for both DES Y6 and LSST Y10 are similar to
those from the BOSS analysis of Cabass et al. (2022b), which finds $c_{s}\geq
0.013$, and the CMB analysis of Planck Collaboration et al. (2020a), which
finds $c_{s}\geq 0.021$. Both are 95% confidence lower bounds corresponding to
our marginalized case above. Given the approximate nature of our estimates
above (as the simulations do not directly vary $c_{s}$ and
$\widetilde{c}_{3}$) we do not interpret the comparison with more detail. Note
that even though the moments integrate over shape information via the volume
integral, they can still distinguish — and thus jointly constraint — different
$f_{\rm NL}$ given their different redshift and scale-dependence (see Section
4.3).

Figure 5: The Fisher information for DES Y6 and LSST Y10 measured using the 2nd and 3rd moments, as a function of source galaxy number density ($n_{\rm gal}$). The fiducial values are $n_{\rm gal}=9\,(27)$ for DES Y6 (LSST Y10), and are denoted by the vertical dotted line. This tests the improvement in constraints from source galaxy number density improvements alone. At fixed $n_{\rm gal}$, the difference between DES and LSST shows the difference in constraints primarily from survey area (with some impact from increasing redshift limits in LSST, see Figure 6). For $f_{\rm NL}^{\rm\,eq}$, where weak lensing is a promising probe, the LSST survey area improves constraints by 30% over DES Y6. The difference between constraints of LSST Y10 at $n_{\rm gal}=10$ (close to the DES Y6 fiducial value) and LSST Y10 $n_{\rm gal}=27$ (the LSST Y10 fiducial value) is 20-30%. Mom. Order | $\sigma(f_{\rm NL}^{\rm\,eq})$ | $\sigma(f_{\rm NL}^{\rm\,loc})$ | $\sigma(f_{\rm NL}^{\rm\,or,\,lss})$ | $\sigma(f_{\rm NL}^{\rm\,or,\,cmb})$
---|---|---|---|---
$N=2$ | 67 [96] | 25 [33] | 75 [95] | 164 [196]
$N=3$ | 99 [184] | 30 [92] | 82 [178] | 128 [180]
$N=4$ | 126 [235] | 39 [173] | 85 [171] | 121 [155]
$N=5$ | 152 [246] | 52 [197] | 97 [165] | 121 [140]
$N\leq 3$ | 51 [66] | 16 [17] | 39 [43] | 68 [75]
$N\leq 4$ | 40 [48] | 13 [13] | 28 [32] | 47 [55]
$N\leq 5$ | 20 [24] | 6 [7] | 13 [16] | 22 [25]
Table 4: Constraints from different orders of the true convergence field (i.e.
noiseless) for LSST Y10-like source galaxy redshift bins/distributions. We
show results both for the individual moments and for combinations that
successively include higher-order moments. Constraints from including up to
the 5th moment improve on those of the 2nd moment alone by factors of 2 to 6.
The square brackets denote constraints after marginalizing over $\Omega_{\rm
m}$, $\sigma_{8}$, and the IA parameters. Constraints from the combination of
moments are degraded by 10% to 20% upon marginalization, whereas those from
the individual moments are impacted more.

Variation with shape noise. We then consider variations of the DES Y6 and LSST
Y10 surveys where the shape noise amplitude is artificially
increased/decreased compared to the fiducial value. This tests the dependence
of the constraints on the noise alone. Therefore, for this test, the source
galaxy redshift distributions, the survey footprints, and all other aspects
are still fixed to fiducial values for each survey. Figure 5 shows the Fisher
information in each survey as a function of source galaxy number density,
which is inversely related to the noise level as shown in Equation 2.17.
Focusing on $f_{\rm NL}^{\rm\,eq}$, the constraints from a DES Y6 survey with
$n_{\rm gal}=24$ are 40% better than those of a fiducial survey with $n_{\rm
gal}=9$. At fixed source galaxy number density, the LSST Y10 constraints are
$\approx 45\%$ better than those of DES Y6. We also compute constraints for an
LSST Y10-like survey but with $n_{\rm gal}=50$, which corresponds to the
number density of a potential weak lensing sample from the Roman space
telescope (Eifler et al. 2021). The constraints, in this case, improve only by
$10\%$ for $f_{\rm NL}^{\rm\,eq}$ and $f_{\rm NL}^{\rm\,or,\,lss}$. Using
simple scaling arguments for the dependence of measurement noise on galaxy
number density and sky area, the Fisher information is roughly proportional to
$\sqrt{n_{\rm gal}}$ and $\sqrt{f_{\rm sky}}$, where $f_{\rm sky}$ is the
fraction of the full sky covered by the survey. The exact scaling with $n_{\rm
gal}$ can deviate from expectations depending on the scale-dependence of the
signal in question.

Interplay of different orders. It is also useful to understand the “true”
information at different orders of the true convergence field. Table 4 shows
this, as probed by the moments. In specific, we extract the information
content in the individual moments, and in their combinations. For the 4th and
5th moments, we only considered the connected piece that is obtained by
subtracting out different combinations of 2nd moments or 2nd and 3rd moments,
respectively. For an auto-correlation with a single field, the expression is

$\langle\kappa^{4}\rangle^{\rm
conn}=\langle\kappa^{4}\rangle-4\langle\kappa^{2}\rangle\langle\kappa^{2}\rangle,$
(4.7) $\langle\kappa^{5}\rangle^{\rm
conn}=\langle\kappa^{5}\rangle-10\langle\kappa^{2}\rangle\langle\kappa^{3}\rangle.$
(4.8)

This modification accounts for the fact that a field with a 2nd and 3rd moment
will automatically have a 4th and 5th moment. The subtraction removes such
contributions and extracts the “connected” 4th and 5th moment. All numbers in
Table 4 correspond to only connected information. The table shows that
including up to the 5th moment improves the constraints from the 2nd moment-
only case by factors of 2 to 6, highlighting the significant information from
inflation contained in the higher-order moments. We have verified that all
numbers are converged to within $1\%$. Note that Figure 2 has already
identified that in the practical limit, the 2nd and 3rd moments contain almost
all of the information. Table 4 shows the impact of the higher-order
information that has been lost due to the larger amplitude of noise (relative
to the 2nd moment case) in these higher-order moments.

###  4.3 Physical origin of $f_{\rm NL}$ signatures in the weak lensing field

The discussions above have thus far established there is valuable information
in the weak lensing field on PNG signatures. It is therefore imperative to
identify how these signatures imprint into this field and into the data-
vectors we study. Such identification will further focus our efforts on
mitigating the relevant systematics and/or on selecting a data-vector that is
more tuned for inflationary signatures.

Scale dependence. Figure 4 already shows that the more non-linear scales
contain the most information, and Figure 2 and 3 confirm this is true for all
summary statistics we discuss here and for both the marginalized and
unmarginalized constraints. Such behavior is expected as varying $f_{\rm NL}$
changes the tails of the initial density distribution, which then changes the
abundance of massive halos and in turn alters the non-linear regime of the
density and lensing fields. Previous works studying $f_{\rm NL}^{\rm\,loc}$
have also identified the abundance of massive halos as the key observable
difference in the lensing field (Dalal et al. 2008; Shirasaki et al. 2012;
Marian et al. 2011; Hilbert et al. 2012). Note, however, that the constraints
in Figure 4 and Table 3 (particularly for LSST) are relatively similar even
under scale cuts of $\theta>10\hbox{${}^{\prime}$}$ or
$\theta>20\hbox{${}^{\prime}$}$; such cuts are employed in the moments-based
weak lensing analyses of Gatti et al. (2020, 2022) with DES Y3 data, and
mitigate the impact of all lensing-based systematics in DES Y3. These angular
scales correspond to physical scales of $10$ to $30\,{\rm Mpc}$ (comoving)
depending on the redshift, which are much larger than the virial radius of
massive halos and thus correspond to the local, large-scale halo environment.
Our constraints after such scale cuts are still comparable/competitive with
BOSS/DESI,888In principle, the scale cuts for LSST may need to be larger than
in DES. This is because scale cuts are designed to minimize the impact on the
$\chi^{2}$ of a given datavector. Since the $\chi^{2}$ increases with improved
precision, this implies that for a given systematic with a fixed amplitude, an
LSST data-vector will require more scales to be cut compared to a DES data-
vector. and this suggests the analysis’ sensitivity to baryon effects is
manageable. We discuss this further in Section 5.2.

Redshift dependence. Figure 6 presents the Fisher information for the
different surveys, but now limits the maximum redshift of the bins used in the
analysis. In each estimate, we use all tomographic bins with a mean redshift
below $z_{\rm max}$. Any cross-correlations between bins below $z_{\rm max}$
with those above $z_{\rm max}$ are also not considered. Given the signatures
in the lensing field arise from modifications to the halo mass function, the
strongest signatures would imprint at low redshift where the structure is most
non-linear. However, the weak lensing amplitude depends directly on the amount
of matter that lenses the source galaxies. Therefore, high redshift source
galaxies contain a larger lensing signal since they are lensed by more
foreground structure. These two opposing effects — the sensitivity of the low-
redshift non-linear density field to $f_{\rm NL}$ but the increased amplitude
of weak lensing signals for high-redshift observations — result in the
sensitivity of $f_{\rm NL}$ from weak lensing peaking at a redshift of $z_{\rm
max}\approx 1.25$. This behavior is seen clearly in Figure 6 for LSST Y1 and
Y10, while for DES we see similar qualitative trends but do not have
tomographic bins with average redshifts beyond $z_{\rm max}\approx 1.25$ to
explicitly confirm this. Note, however, that this discussion above concerns
signatures (and thus constraints) of only $f_{\rm NL}$. We have verified that
in analyses that also marginalize over $\Omega_{\rm m}$, $\sigma_{8}$, and the
IA parameters, the higher redshift bins still add information (though, the
choice $z_{\rm max}\approx 1.25$ still approximately maximizes constraints).
This improvement is because the high redshift bins, by virtue of their larger
lensing signal, have significant information on cosmology and IA parameters
and thus improve the marginalized constraints on $f_{\rm NL}$.

Figure 6: The Fisher information in the 2nd and 3rd moments of the lensing
field, shown as a function of the maximum redshift of the data-vector, for
four different survey configurations. The maximum redshift is enforced by
removing all tomographic bins with a mean redshift of $z_{\rm mean}>z_{\rm
max}$. The existing constraints from BOSS are shown as the black line and the
potential constraint from DESI are in red. The Fisher information saturates
near $z_{\rm max}=1.25$, due to opposing effects of the lensing amplitude
growing towards high redshift and non-linear evolution growing towards low
redshift. Note that we only vary the $f_{\rm NL}$ parameter and do not
marginalize over cosmology and IA parameters. This choice is consistent with
the BOSS/DESI estimates shown above, which fix the cosmological parameters.
Figure 7: The fractional differences in the halo mass function (HMF) — which
is the counts of halos as a function of mass — as we vary $f_{\rm NL}$. The
colors show results from different redshifts, and the panels show different
$f_{\rm NL}$ types. The variation is taken as $\ln{\rm HMF}(f_{\rm
NL}=100)-\ln{\rm HMF}(f_{\rm NL}=-100)$, which gives us the fractional
deviation. The cosmic variance term is suppressed as the simulations with
$f_{\rm NL}=\pm 100$ use the same random seed for the initial conditions.
Figure 8: The fractional change in halo bias due to $f_{\rm NL}$, shown as a
function of multipole, $\ell$, for the halo samples of different redshift
bins. The variation is computed as $\ln b_{h}(f_{\rm NL}=100)-\ln b_{h}(f_{\rm
NL}=-100)$. The sample has a fixed lower mass cut of $M>10^{14}\,{\rm
M}_{\odot}$ at every redshift, and this is set by the resolution limit of the
simulation. The bias depends on $f_{\rm NL}$ over a wide range of scales,
including both linear and non-linear scales. The cosmic variance term is
suppressed as the simulations with $f_{\rm NL}=\pm 100$ use the same random
seed for the initial conditions.

Halo mass function. We can then also directly explore the variation in the
halo mass function (HMF) as we change $f_{\rm NL}$, as this will then imprint
on the weak lensing field. Figure 7 shows the fractional change in the HMF for
$\Delta f_{\rm NL}=200$, for all four types of $f_{\rm NL}$, and for five
different redshifts. Given the particle resolution of the Ulagam suite and the
requirement of at least 100 particles per halo, our halo catalogs are limited
to $M>10^{14}\,{\rm M}_{\odot}$. While the weak lensing signal is also
sensitive to masses below this scale, studying the HMF for such masses is
still informative in understanding the qualitative impact of $f_{\rm NL}$.
Figure 7 presents a clear impact of $f_{\rm NL}$ on halo counts towards the
massive halo end. The sign of the fractional change is positive for $f_{\rm
NL}^{\rm\,eq}$ and $f_{\rm NL}^{\rm\,loc}$, meaning the abundance of massive
halos increases with $f_{\rm NL}$, and this is expected as for positive
$f_{\rm NL}$ the initial conditions have a larger skewness and have more high-
density peaks. At $z=0$, the HMF for lower mass halos, $M\approx 10^{14}\,{\rm
M}_{\odot}$, is reduced with increasing $f_{\rm NL}^{\rm\,loc}$ and $f_{\rm
NL}^{\rm\,eq}$. This sign flip is seen more prominently in both orthogonal-
type $f_{\rm NL}$, though the relation is reversed as high $f_{\rm NL}$
implies a lower halo count. For $f_{\rm NL}^{\rm\,or,\,cmb}$, the change is
factors of 2-3 smaller than for all the other types. For both orthogonal-type
$f_{\rm NL}$ we once again find that the sign of the change is flipped at
redshift $z=0$, for halos with $M\approx 10^{14}\,{\rm M}_{\odot}$, as their
abundance increases with $f_{\rm NL}$ now. The redshift- and mass-dependence
of the sign flip are consistent with findings in Jung et al. (2023). This has
also been analytically derived in LoVerde et al. (2008, see their Section 4.2)
and is due to the fact that for a fixed matter energy density, increasing
$f_{\rm NL}$ causes more massive halos to form at the expense of less matter
available to form lower mass systems. Coulton et al. (2022, see their Figure
1) also show that the increasing $f_{\rm NL}^{\rm\,or,\,lss}$ reduces the
amplitude of the power spectra on small scales, while increasing $f_{\rm
NL}^{\rm\,or,\,cmb}$ has a much more mild effect that is nearly non-existent
at $z=0$. While these results are for the power spectrum, they can be
translated to signatures in the HMF given halo formation defines the structure
of the power spectrum on small scales. Thus, more/less massive halos imply
more/less power on small scales.

Halo bias. Given the above discussion on the HMF, a natural conclusion is that
the optimal statistic for inflationary constraints from non-linear structure
is the HMF itself, where the latter can be measured through the counts of
galaxy clusters (e.g., Abbott et al. 2020; Costanzi et al. 2021). While
practical considerations motivate weak lensing, in comparison to the HMF, as a
simpler observable to robustly analyze, there are theoretical motivations as
well. The signature from inflationary models associated with these $f_{\rm
NL}$ types is not solely in the number of high-density peaks of the initial
conditions, but also in the way the peaks are spatially clustered. One of the
first, well-known signatures of $f_{\rm NL}$ on the halo field is the scale-
dependent bias found in the local-type $f_{\rm NL}$ (Dalal et al. 2008). This
bias goes as $b\propto f_{\rm NL}^{\rm\,loc}/k^{2}$ and is a diverging signal
for $k\rightarrow 0$. Other types of $f_{\rm NL}$ can also have scale-
dependent biases (Schmidt & Kamionkowski 2010). This bias can imprint on the
small-scale density/lensing field in the “two-halo” term/regime, which is
comprised of the signal from neighboring halos and therefore depends on the
halo clustering. We can compute the scale-dependent halo bias in the Ulagam
suite as the ratio of the measured halo-matter angular power spectrum and the
matter-matter angular power spectrum,

$b_{\ell}=C_{\ell}^{hm}/C_{\ell}^{mm}\,,$ (4.9)

where $C_{\ell}$ is the power at different multipole $\ell$. The choice of
using $C_{\ell}^{hm}$ over $C_{\ell}^{hh}$ removes the impact of shot noise in
our estimate of $b_{\ell}$.

Figure 8 presents the fractional change in the halo bias as a function of
$\ell$, for four redshifts and for the four types of $f_{\rm NL}$. Since this
quantity is estimated on mock full-sky maps, and not 3D real-space fields, we
do not compute/show the $z=0$ trends. We reproduce the result of Dalal et al.
(2008) for $f_{\rm NL}^{\rm\,loc}$, where the bias of the halo sample grows at
large scales (low $\ell$). For $f_{\rm NL}^{\rm\,loc}$ we also see a change in
the sign of the bias on small scales to negative values. Similar but weaker
features are seen in $f_{\rm NL}^{\rm\,or,\,cmb}$, which is expected to have a
scaling of $b\propto 1/k$ on linear scales (Schmidt & Kamionkowski
2010).999The $f_{\rm NL}^{\rm\,or,\,cmb}$ template is an approximation of the
full EFT template (Senatore et al. 2010). The $b\propto 1/k$ scaling of the
former template is considered an artifact of this approximation as the scaling
is not found in the latter. In our work, we use $f_{\rm NL}^{\rm\,or,\,cmb}$
only to broaden the range of non-Gaussian templates whose signatures we study
with weak lensing (and we do not use it to constrain any EFT parameters), in
which case differences between the approximate, $f_{\rm NL}^{\rm\,or,\,cmb}$
template and the true, EFT template are inconsequential. Finally, the
equilateral type shows the halo bias effect is nearly scale independent for
$\ell<200$; meaning, $f_{\rm NL}^{\rm\,eq}$ simply alters the linear bias. For
$\ell>200$, there is a mild scale-dependent, redshift-dependent behavior,
though the amplitude of this variation is still low. The bias for $f_{\rm
NL}^{\rm\,or,\,lss}$ shows similar features to that of $f_{\rm NL}^{\rm\,eq}$,
with a scale-independent bias below $\ell<200$, and a mildly dependent one
above it. The behaviors of the bias on large scales, for all $f_{\rm NL}$
types, are consistent with the perturbation theory predictions of Schmidt &
Kamionkowski (2010).

All scale-dependent signatures of the different $f_{\rm NL}$ are completely
unused if we choose the halo counts as our statistic. On the other hand, the
weak lensing field (or any large-scale structure field in general) is
sensitive both to the number of massive halos and to the halos’ spatial
clustering properties. Thus, the lensing field utilizes more available
information than the halo counts alone. This is also consistent with the
lensing field’s sensitivity to $f_{\rm NL}$ even on scales of
$\theta>10\hbox{${}^{\prime}$}$, where this angular scale corresponds to $5$
to $20\,\,{\rm Mpc}$ across different redshifts, as such scales probe the
“two-halo” term and the signatures of $f_{\rm NL}$ imprinted into this term.

Figure 8 also shows that the impact of $f_{\rm NL}$ on the halo bias grows
with redshift. This, however, is an artifact of selection effects induced by a
common halo mass cut across all redshifts. A halo of $M>10^{14}\,{\rm
M}_{\odot}$ at $z=2$ is a significantly rarer structure than a halo of
$M>10^{14}\,{\rm M}_{\odot}$ at $z=1$. Thus, our fixed mass cut is selecting
rarer structures at higher redshift, and therefore the halo bias of the high
redshift samples will be larger.

While it may appear contradictory that we discuss the halo bias in a weak
lensing field — which we describe above as an unbiased, direct tracer of the
density field — this is still consistent with our previous statement that the
lensing field is insensitive to the galaxy-halo connection. On small scales
the clustering of matter can be represented as the combination of three halo
properties: their spatial clustering, their number density, and their density
profiles. This is denoted the “halo model” approach (Cooray & Sheth 2002) and
utilizes (among other quantities) the halo bias. The galaxy bias and the
galaxy-halo connection do not appear in the halo model prediction for the
density/lensing field. We decompose the signatures of $f_{\rm NL}$ in the weak
lensing field into signatures in the halo mass function and the halo bias, as
this can be a more intuitive picture for understanding the physical origin and
scale-dependence of the weak lensing signatures. Both the HMF and the halo
bias features discussed above are statistically significant even if we account
for cosmic variance, where the latter is generally factors of 3 to 5 larger
than the uncertainties shown in Figure 7 and 8.

##  5 Discussion

Having shown that weak lensing can provide usable constraints on $f_{\rm NL}$,
we now discuss in §5.1 its unique advantage as a probe, and in §5.2 the
potential modeling challenges, including existing methods to alleviate them.

###  5.1 Advantages of weak lensing as a probe of $f_{\rm NL}$

Weak lensing has a number of advantages as a cosmological probe, which have
aided its development as a leading probe for constraining cosmological
parameters (e.g., Asgari et al. 2021; Abbott et al. 2022; More et al. 2023).
Here, we highlight some of these advantages that are specific to the $f_{\rm
NL}$ signatures we study.

Unbiased, direct tracer of the density field. As mentioned before, the primary
advantage of weak lensing as a cosmological probe is that it is a direct
tracer of the density field. A vast majority of cosmological signals imprint
directly into the correlations of the density field. Historically, the spatial
correlations of galaxy fields are a more popular probe as the measurement
signal-to-noise is high. However, such analysis requires either
marginalization, or prior knowledge, of the galaxy bias parameters. These
parameters are required to translate the correlations of the galaxy field,
which is the key observable, into those of the density field, which contains
the physical signatures of interest. In many cosmological analyses using both
lensing and galaxies, the former provides a significant fraction of total
constraining power — even though it is a lower signal-to-noise measurement
than the latter; the difference in DES Y3 is 50% (Rodríguez-Monroy et al.
2022; Secco et al. 2022a; Amon et al. 2022) — as it does not need to
marginalize over the unknown galaxy bias. It is therefore justified to assume
weak lensing provides substantial/comparable information on $f_{\rm NL}$ when
compared to the information from galaxy correlations. The results of Section
4.2 confirm this to be the case.

Ease of simulation-based modeling. Any simulations used to model the weak
lensing field are defined by two scales: the volume of the survey which sets
the simulation size, and the smallest scale in the analysis of the lensing
field which sets the simulation resolution. For example, if the analysis does
not use scales below $10\,{\rm Mpc}$, the simulation can simply be run to be
accurate only above those scales. This enables the production of large suites
of simulations with coarse resolution that can still be used to infer
cosmological constraints. For simulation-based modeling of galaxy fields,
however, there is an additional scale in the size of halos/galaxies. Modeling
the galaxy field starts from the accurate simulating of halos, proceeded by
the assignment of galaxies to these halos (using some form of the galaxy–halo
connection). Thus, the smallest scale that must be resolved in the simulation
is a few times smaller than the smallest halo size. The lack of such a
limitation for weak lensing has led to multiple large suites of simulations,
some with $N_{\rm sim}\sim\mathcal{O}(10^{3}-10^{4})$, being developed and
used in weak lensing analysis (Zürcher et al. 2021; Zürcher et al. 2022; Gatti
et al. 2023; Kacprzak et al. 2023). This then directly enables the use of non-
linear scales in the lensing measurements.101010If the chosen non-linear
scales are sufficiently small, then the weak lensing fields will enforce the
same small-scale requirement as the galaxy field for simulation-based
modeling. These non-linear scales can currently be modelled only through
simulations as analytic approaches are inaccurate in this regime. There are,
however, ongoing efforts for hybrid approaches that combine the ideas of
EFTofLSS with the non-linear predictions of N-body simulations and can thereby
extend the range of usable scales (Modi et al. 2020; Kokron et al. 2022;
Banerjee et al. 2022).

Constraining galaxy bias parameters with cross-correlations. The constraints
in Figure 4 motivate the use of weak lensing-only data as a probe of $f_{\rm
NL}$. However, following the usage of weak lensing in large surveys, we can
infer that of equal importance — if not more — is the usage of the lensing-
galaxy cross-correlation. This correlation constrains, or “self-calibrates”,
the galaxy bias parameters and thus enables/improves the constraints from the
galaxy clustering measurements. The latest analyses of $f_{\rm NL}$ from
spectroscopic galaxy surveys utilize a one-loop bispectrum model, which has
$\mathcal{O}(10)$ bias parameters compared to the single linear galaxy bias
parameter used in most analyses of photometric surveys. Thus, the potential
improvement in constraints, due to the self-calibration of bias parameters via
the lensing-galaxy cross-correlation, will be greater in this scenario.
Philcox et al. (2022a, see their Section 7) also identified that the
uncertainty in these bias parameters as the biggest limitation in improving
the theoretical modeling of the perturbation theory approach. There are also
indications that some assumptions and/or choices of priors for the bias
parameters need to be relaxed (e.g., Barreira 2022; Brinch Holm et al. 2023),
in which case self-calibration via lensing-galaxy cross-correlations can
provide an even larger benefit. Exploring this cross-correlation requires a
common modeling framework, and the hybrid-EFT approach mentioned above is a
potential method forward.

Accessing smaller scales, $\boldsymbol{k>1{\rm Mpc}^{-1}}$. Secco et al.
(2022a, see their Figure 4) show that a lensing measurement at angular scale
$\theta$ probes a range of wavenumbers $k$, represented heuristically as

$\xi_{\kappa}(\theta)=\int_{0}^{\infty}P_{\kappa}(k)w(k)\,,$ (5.1)

where $w(k)$ is the weight of each mode, quantifying the sensitivity of
measurement $X(\theta)$ to this mode, $P(k)$ is the convergence power
spectrum, and $\xi_{\kappa}(\theta)$ is the convergence 2-point function. We
use the same method of e.g., Secco et al. (2022a), introduced in Tegmark &
Zaldarriaga (2002), to estimate $w(k)$ for the shear 2-point correlation
functions — which corresponds to the 2nd moments measured in our work — while
accounting for the redshift distribution of the different tomographic bins. We
will focus on LSST Y10 but note that the results for other surveys are
similar. The measurements at the minimum scale of
$\theta=3.2\hbox{${}^{\prime}$}$ correspond to mean wavenumbers of
$1<k\,\,[h/{\rm Mpc}]<2.5$. This is computed as the weighted average of $k$,
with weights $w(k)$. If we instead consider the maximum contributing
wavenumber — defined here as the maximum scale with a weight, $w(k)$, that is
at least 10% of the maximum weight $\texttt{max}(w(k))$ — we find
$3.5<k\,\,[h/{\rm Mpc}]<6$, depending on the tomographic bin. Even under our
conservative angular scale cut of $\theta_{\rm min}=20\hbox{${}^{\prime}$}$,
we find the maximum contributing wavenumber is $0.5<k\,\,[h/{\rm Mpc}]<2.6$.
The CMB analyses of Planck Collaboration et al. (2020a) probe $\ell_{\rm
max}=2500$, which for the CMB redshift of $z=1100$ corresponds to scales of
$k\lesssim 0.02\,\,h/{\rm Mpc}$. The galaxy correlation function analyses of
Cabass et al. (2022a, b); Philcox et al. (2022a); D’Amico et al. (2022) use up
to $k\lesssim 0.2\,\,h/{\rm Mpc}$, which is set by the current accuracy of the
EFT model above this scale. Note that our maximum scale sensitivity will also
be limited by modeling choices. However, for angular scale cut of $\theta_{\rm
min}=20\hbox{${}^{\prime}$}$, which is consistent with the choices of DES Y3
(Gatti et al. 2022), the model is accurate and the measurement still accesses
smaller scales than those of the CMB and galaxy correlations. We discuss in
Section 5.2 the modeling challenges in accessing even smaller scales. The
sensitivity of lensing measurements to higher wavenumber than the other probes
provides the opportunity to study scale-dependent PNGs, especially when
combined with other, small-scale measurements from the CMB (Emami et al. 2015,
see their Figure 1). Such scale-dependence is directly connected to
inflationary interactions, such as a change in the sound speed $c_{s}$ over
time (e.g., Planck Collaboration et al. 2020a, see their Section 2.4.2).

###  5.2 Modeling challenges

Simulation-based modeling of weak lensing fields has already been utilized in
current surveys to perform precision cosmology (Fluri et al. 2019; Fluri et
al. 2022; Zürcher et al. 2022) and many more have used simulations for certain
aspects of the modeling pipeline (e.g., Secco et al. 2022a; Amon et al. 2022;
Gatti et al. 2022). However, simulation-based modeling will face new
challenges in future surveys, where improved measurement precision requires
improved modeling techniques. We detail some of these upcoming challenges
below:

Higher order shear. In Section 2.2, we discuss that the observable of a weak
lensing survey are the shear fields $\gamma_{1,2}$. However, this is an
approximation as the galaxy shapes actually trace, $e^{\rm
obs}_{1,2}=\gamma_{1,2}/(1-\kappa)$, where $\kappa$ is the convergence field.
In the limit of $\kappa\ll 1$, the expression can be Taylor-expanded to
$e^{\rm obs}_{1,2}=\gamma_{1,2}(1+\kappa+\kappa^{2}/2+\ldots)$. The assumption
of $e^{\rm obs}_{1,2}\approx\gamma_{1,2}$, which we have used in this work, is
the “reduced shear” approximation. It has been explicitly verified to be a
negligible effect for multiple statistics in DES Y3 data (Krause & Hirata
2010; Gatti et al. 2020; Anbajagane et al. 2023a). However, simple scaling
arguments suggest it will be a statistically significant effect for LSST Y10
precision. Of similar impact is the magnification effect, where more source
galaxies are observed in directions with more foreground structure (larger
convergence). The impact is modelled as $\propto(1+q\kappa)$, where
$q=\mathcal{O}(1)$. Thus, the magnification impact is similar to that of the
reduced shear above, which implies it will also be a notable effect for LSST
Y10. Yet another higher-order shear effect is source clustering, which is the
correlation of source galaxy positions with the foreground convergence field.
This effect arises because source galaxies and lensing convergence are both
tracers of the underlying density field, and has been observed in DES Y3 data
with different statistics (Gatti et al. 2023; Anbajagane et al. 2023a). For
the 2nd and 3rd moments, its impact can be greatly minimized (Gatti et al.
2022, 2023). In general, all these effects can be included in (simulation-
based) forward modeling approaches at low computational cost.

Intrinsic Alignments. We have already shown the impact of intrinsic alignments
in Figure 2. This test utilized a specific parameterization of IA, with
theoretically motivated but fixed parameter values. We do not have
observational constraints as the weak lensing data from DES Y3 does not
identify an IA signal (Secco et al. 2022a; Amon et al. 2022). Table 3 shows
that marginalizing over IA (in addition to $\Omega_{\rm m}$ and $\sigma_{8}$)
leads to LSST Y10-based constraints that are still comparable to DESI for
multiple types of $f_{\rm NL}$. The IA approach we have used, NLA, can be
thought of as similar to the hybrid EFT approach, where the underlying
framework is that of perturbation theory while the non-linearities in the
density field are modelled through N-body simulations. While it is possible to
add higher-order terms through the “Tidal-alignemnt tidal-torquing” (TATT,
Blazek et al. 2019) model, the data is currently not precise enough to show a
preference for TATT over NLA. The NLA model has been used extensively for
various lensing-related analyses (e.g., Secco et al. 2022a; Amon et al. 2022;
Gatti et al. 2022). The requirements for the IA modeling in an LSST Y10
dataset are unclear. If the current, fiducial IA model continues to prove
adequate, then we find the IA modeling is not an issue.

Other lensing nuisance parameters. The full analysis of the lensing 2-point
correlations also includes marginalizations over other “nuisance” parameters
that alleviate any systematics-driven biases. These parameters include the
mean redshift of the galaxies in each tomographic bin and a multiplicative
bias in the estimated shapes of galaxies per bin. The latter is a measurement
bias arising primarily from the blending of source galaxies in the image. In
current surveys, marginalizing over these parameters leads to negligible
impact on the total constraints, in comparison to marginalizing over IA. We
can infer this will be more true if we also marginalize over cosmology as we
do in this work. It is highly likely that this behavior continues to be the
case for LSST, under our current understanding of IA. The effects associated
with these nuisance parameters can be easily included in the simulation
modeling, and this has already been utilized in multiple analyses (e.g.,
Zürcher et al. 2021; Fluri et al. 2022).

Baryon Imprints. A significant systematic in all analyses related to the
density field is the impact of baryonic evolution (e.g., Chisari et al. 2018,
see their Figure 6). All existing simulation-based models (of either weak
lensing or galaxy correlations) employ N-body simulations. While it is
possible to use hydrodynamic simulations with galaxy formation models to
perform the modeling, such models require many assumptions on the nature of
galaxy formation. The assumptions required are often on processes like gas
cooling and AGN (Active Galactic Nuclei) feedback, which are the dominant
physical processes behind alterations of the density distribution in and
around halos (Blumenthal et al. 1986; Gnedin et al. 2004; Duffy et al. 2010;
Anbajagane et al. 2022a; Shao et al. 2022). Given the range of possible,
allowed assumptions, the simulations manifest a variety of halo property
behaviors. Comparative studies show the predictions agree in the overall
trends but differ in specific details (e.g., Anbajagane et al. 2020; Lim et
al. 2021; Lee et al. 2022; Cui et al. 2022; Stiskalek et al. 2022; Anbajagane
et al. 2022a; Anbajagane et al. 2022b). Studies on the thermodynamic
properties of gas also find differences between the measurements from data and
the predictions from these hydrodynamic simulations (e.g., Hill et al. 2018;
Amodeo et al. 2021; Pandey et al. 2022; Anbajagane et al. 2022c, 2023b).

An alternative approach to modeling baryons is “baryonification” (Schneider et
al. 2019), which is a flexible, halo-based model that alters the density field
in an N-body simulation to include the baryon imprints. This technique
provides a higher-level, approximate galaxy formation model that depends only
on “macro” properties like the halo baryon fraction, the baryon density
profiles, dark matter density profile etc. The technique has already been
utilized in previous analyses of widefield surveys (Fluri et al. 2022; Chen et
al. 2023; Aricò et al. 2023). The model flexibility/utility has been shown for
the power spectrum/2-point functions — which are directly related to the 2nd
moments we use — up to $k=10\,\,h/{\rm Mpc}$ (Schneider et al. 2019; Giri &
Schneider 2021) and for some higher-order statistics down to
$\theta=1\hbox{${}^{\prime}$}$ scales (Lee et al. 2023), but not for the 3rd
moments used in this work. Thus, additional validation work is required to
apply this model to the data-vector we employ here. Such baryon correction
models will become increasingly necessary for LSST data as simple scale cuts
to remove “baryon contaminated” measurements, akin to those used in DES Y3,
will remove a larger portion of non-linear scales given the increased
precision in the LSST data. While Figure 4 suggests the LSST Y10 constraints
after scale cuts will still be comparable to DESI, these constraints would be
much improved by including such non-linear scales and marginalizing over
baryon evolution instead.

##  6 Conclusion

We explore the weak lensing field as a potential probe of primordial non-
Gaussianities (PNGs), where the amplitude of PNGs is denoted by the $f_{\rm
NL}$ parameter. We consider four types of $f_{\rm NL}$ — $f_{\rm
NL}^{\rm\,loc}$ which arises from multi-field inflation, $f_{\rm
NL}^{\rm\,eq}$ from a strong self-coupling of the inflaton field, $f_{\rm
NL}^{\rm\,or,\,lss}$ and $f_{\rm NL}^{\rm\,or,\,cmb}$ from the same physics as
$f_{\rm NL}^{\rm\,eq}$ but corresponding to different interactions (see
Section 2 for more details) — and run N-body simulations to extract the Fisher
information in DES-like and LSST-like surveys for each type. The Ulagam suite
of simulations allows us to forward model wide-field surveys and use physical
scales that probe deep into the non-linear regime. Our findings are summarized
as follows:

  * •

When varying just $f_{\rm NL}$, the two-point correlation function — as traced
by the 2nd moment of the field — provides constraints comparable to those from
higher-order statistics for $f_{\rm NL}^{\rm\,loc}$ and $f_{\rm
NL}^{\rm\,eq}$. However, the latter are clearly better for $f_{\rm
NL}^{\rm\,or,\,lss}$ and $f_{\rm NL}^{\rm\,or,\,cmb}$, and for analyses of all
$f_{\rm NL}$ that marginalize over cosmology ($\Omega_{\rm m}$ and
$\sigma_{8}$) and intrinsic alignment parameters (Figure 2).

  * •

The shape/configuration information in the lensing field adds significantly to
$f_{\rm NL}$ constraints. The full 3-point correlation function is
significantly more constraining than the 3rd moment, which integrates over the
shape/configurations (Figure 3). A computationally inexpensive implementation
is still challenging.

  * •

Using the combination of 2nd and 3rd moments as our fiducial statistic, we
find the weak lensing-based constraints can be comparable or potentially
better than galaxy clustering-based constraints from spectroscopic surveys.
For $f_{\rm NL}^{\rm\,eq}$, LSST Y10 (DES Y6) is competitive/better than DESI
(BOSS). For $f_{\rm NL}^{\rm\,or,\,lss}$, the LSST Y10 constraints are a
factor of 1.5 to 2 broader than the DESI constraints (Figure 4).

  * •

The LSST constraints are still comparable to DESI (within factor of 1.5) even
with conservative scale cuts of $\theta>20\hbox{${}^{\prime}$}$ (Table 3). At
such scales, all systematics associated with weak lensing — as seen in DES Y3
— are know to be alleviated.

  * •

The redshift dependence of the signal peaks for source galaxies of $z\approx
1.25$. Including source galaxies beyond this redshift helps modestly with the
constraining power when varying only $f_{\rm NL}$. However, the high redshift
data is valuable in marginalized $f_{\rm NL}$ constraints, as it helps in
constraining the cosmological and IA parameters that are marginalized (Figure
6).

  * •

The impact of $f_{\rm NL}$ on the lensing field, including the range of scales
that probe $f_{\rm NL}$, can be understood through the impact of $f_{\rm NL}$
on the halo mass function and the halo bias (Figure 7, 8).

Our results show that weak lensing on its own is a useful probe for analyses
of $f_{\rm NL}$, and enhances the search for scale-dependent PNGs. Including
cross-correlations between the weak lensing and galaxy fields can result in
more significant benefits (see discussion in Section 5.1). While a large part
of our discussion has centered around Stage IV surveys, such as LSST and DESI,
there is also potential for a DES Y6 analysis — particularly of $f_{\rm
NL}^{\rm\,eq}$ (and possibly $f_{\rm NL}^{\rm\,or,\,lss}$) — as it is
comparable to current constraints from BOSS. This would serve as a pathfinder
analysis building towards performing such an analysis with LSST. Improving
constraints on $f_{\rm NL}^{\rm\,eq}$ and $f_{\rm NL}^{\rm\,or,\,lss}$ will
help probe the energy scales of inflation and explore one of the energy
frontiers in cosmology (Achúcarro et al. 2022).

We also show that the configuration information in the 3-point correlations is
significant. However, the extraction of such information from the traditional
weak lensing measurement is not ideal given the weak lensing field is a
projected integral of the density field. Methods exist to reconstruct a
“three-dimensional mass map”, where the lensing fields in different
tomographic bins are used to reduce the contribution from the
foreground/background structure from each bin (Simon et al. 2009). The final
maps are noisier but are better suited for extracting the shape/configuration
information, which can improve the $f_{\rm NL}$ constraints (Figure 3).

Finally, while our work focuses on inflation, and primordial non-Gaussianities
in particular, we emphasize that the arguments made above (particularly in
Section 5) can easily extend to any early universe physics that affects the
initial conditions. Such physics has thus far been constrained primarily using
the galaxy clustering probe, and often with just the two-point function. Weak
lensing was not considered a competitive probe of such physics given the
limitations in analytic modeling the lensing field, and the reduction in
signal amplitude (compared to the amplitude in the 3D density field) due to
the projection integral. However, the computing advances of recent times allow
for simulation-based modeling to replace the purely analytic approach, and
consequently use non-linear scales that are beyond the reach of current
analytical modeling approaches. Weak lensing as a probe of new physics must be
revisited in light of this shift. In upcoming work, our simulation-based
modeling approach will be extended to explore the use of weak lensing in
constraining a multitude of different early Universe physics.

## Acknowledgements

DA is supported by NSF grant No. 2108168. CC is supported by the Henry Luce
Foundation and DOE grant DE-SC0021949. HL is supported by the Kavli Institute
for Cosmological Physics at the University of Chicago through an endowment
from the Kavli Foundation and its founder Fred Kavli.

We thank Will Coulton, Cyrille Doux, Daniel Eisenstein, Giulio Fabbian, Doug
Finkbeiner, Josh Frieman, Wayne Hu, Bhuv Jain, Austin Joyce, Nick Kokron, and
Lucas Secco for helpful conversations during various stages of this work. We
also thank Francisco Villaescusa-Navarro for graciously hosting the Ulagam
simulation data and for advice on our public data release. All analysis in
this work was enabled greatly by the following software: Pandas (McKinney
2011), NumPy (Van der Walt et al. 2011), SciPy (Virtanen et al. 2020), and
Matplotlib (Hunter 2007). We have also used the Astrophysics Data Service