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

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 their amplitudes, fX 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 σ(f NL eq) ≈ 110, σ(f NL or, lss) ≈ 120, σ(f NL loc) ≈ 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 ≲ 1.25), with a clear origin in the evolution history of massive halos. The constraining power degrades by ∼60% under scale cuts of ≳ 20 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 ΛCDM model, where CDM stands for cold dark matter and Λ is the cosmological constant.The large-scale structurewhich 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 Λ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 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 GeV (giga electron volts), which is much larger than energy scales achievable in terrestrial particle physics experiments, the parameter f NL probes what has been denoted an "energy frontier" in cosmology (Achúcarro et al. 2022).Note that f 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 τ 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 NL are found in the bispectrum analysis of the Planck CMB data (Planck Collaboration et al. 2020a) -σ(f eq NL ) = 5, σ(f eq NL ) = 47, σ(f or, cmb NL ) = 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 ∝ V k 3 max , while in a 2D CMB map, the number of modes scales as N k ∝ Ak 2 max .Here, V is the 3D survey volume, A is the 2D map area, and k 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 scaledependent f 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 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,1 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 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 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 (≥ 10%) in the non-linear regimes (k ≳ 1h/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. 2018Fluri 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 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 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 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 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 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.

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, κ, which is a line-of-sight integral of the density field, where z s is the redshift of the "source" plane/galaxies being lensed, n is the pointing direction on the sky, δ is the overdensity field, χ is the comoving distance from an observer to a given redshift, a is the scale factor, H 0 is the Hubble constant, Ω m is the matter energy density fraction at z = 0, and c is the speed of light.We use the shorthand χ(z s ) ≡ χ s and χ(z j ) ≡ χ 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. 2 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. 2021Zü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 Gpc/h ≈ 1.5 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 wCDM and PNG parameters; these are Ω m , σ 8 , w, and n s for the ΛCDM case and different f 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 ∆Ω m = 0.02, ∆σ 8 = 0.03, ∆w = 0.05, ∆n s = 0.04, and ∆f type NL = 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 ∆z ≈ 0.01−0.05 in the latter redshift range; the exact value of ∆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 ≲ 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.objects with at least 100 particles, which corresponds to a minimum mass of M > 10 14 M ⊙ 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.

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, ϕ(x), or its Fourier equivalent ϕ(k), we can generate a field, Φ, with a chosen bispectrum as is the Dirac delta function enforcing momentum conservation and K(k 1 , k 2 ) is a coupling kernel that contains information about the chosen bispectrum.By choosing different K(k 1 , k 2 ), we can generate density fields with different PNGs.
The bispectrum of the field Φ(k) defined in Equation 2.2 above is Given a particular bispectrum template, one can find the coupling kernel, K(k 1 , 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 f loc NL is the most well-studied f 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, and the field itself can be generated easily by adding the square of the Gaussian field to the linear term, where the ensemble average must be subtracted out to enforce that the field, ϕ(x) 2 , has zero mean and is thus a purely perturbative field that does not alter the mean value of Φ loc (x).
The bispectrum of this model peaks at "squeezed" configurations, k 1 ≪ 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 f eq NL 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 Such PNGs are generated from inflation models that have "non-canonical" kinetic terms, and their amplitude peaks in the limit k 1 ≈ k 2 ≈ 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 (2.7) Third, the CMB Orthogonal-type f or, cmb NL 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 which leads to the real-space expression (2.9) Fourth and finally, the LSS Orthogonal-type f or, lss NL 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 or, cmb NL .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 where the constant p is given by . (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 eq NL , f or, lss NL , f or, cmb NL -depend on the selfcoupling of the inflationary perturbations.There is a natural theoretical threshold for these PNG at f NL ∼ 1.If f 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 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 or, lss NL and f or, cmb NL -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 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.

Simulating DES and LSST skies
In this work, we are interested in the impact of f 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 Nbody 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 deg 2 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 deg 2 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 where N i 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 κ 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,   DES Y6, and LSST Y1 and Y10 are parameterized as, 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.
Constructing lensing shear shells.Weak lensing surveys do not directly observe the convergence, κ, as their main observable is the galaxy shapes that are tracers of the shear field, γ. 3 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 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 where κ IA (n, z) is the convergence signal corresponding to IA, δ(n, z) is the density field, A IA is the amplitude of the IA effect, ρ crit is the critical density of the universe at z = 0, Ω m is the fraction of matter energy density at z = 0, D + is the linear growth function normalized to D + (z = 0) = 1, and η IA is the redshift scaling.Both A IA and η IA are free parameters, and will be referred to as IA parameters.The convergence field with IA is then simply κ → κ + κ 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 η + IA = −1.6,η − IA = −1.8 and A + IA = 0.8, A − IA = 0.6.These are variations of 0.1 around the fiducial values of A IA = 0.7 and η IA = −1.7 as chosen in Krause et al. (2021, see their Table 2).4Through 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 where n gal is the source galaxy number density, and A pix is the pixel area for a given map resolution.All maps in this work use NSIDE = 1024, corresponding to a pixel resolution of 3.2 ′ .The per-galaxy shape noise is taken to be σ 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 catalog5 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(Gatti et al. , 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 deg 2 each, which is the expected area coverage for Y10 and ≈ 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, κ 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 crosscorrelations 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. 2021Zü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. 2020Gatti et al. , 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.

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(Gatti et al. , 2023)), (ii) the CDFs, which have been tested for DES Y3 data (Anbajagane et al. 2023a), and (iii) the threepoint 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 3point 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.

Moments
The moments of the field provide an efficient way to summarize the information from different orders.They are computed as where κ (j) is the convergence field of the j th tomographic bin, N pix is the number of pixels in the survey footprint.In all cases, ⟨κ (j) ⟩ = 0 is enforced and the scale dependence on θ is obtained by making measurements after smoothing the fields with a harmonic-space tophat filter, where θ is the tophat radius in radians.We use ten bins of θ, between 3.2 ′ < θ < 200 ′ , which follows the analysis choices of Gatti et al. (2022).All fields are smoothed by the same θ.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 θ.The quantity κ i in Equation 3.1 is the same as the κ(n) defined in equation 2.1, but with the continuous direction 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 3point 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.

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 (κ).The CDF measurements are defined as, where Θ(κ − ν) is the Heaviside step function, which takes values of 1 if κ − ν > 0 and 0 otherwise, and ⟨κ (j) ⟩ = 0 is enforced similarly to the moments measurement.The quantity CDF(θ, ν) captures what fraction of the field, smoothed on a scale θ, is above a threshold ν.
The ν have the same units as κ and are dimensionless quantities.The choice of smoothing scales is the same as that of the moments, 3.2 ′ < θ < 200 ′ , and for each scale, we use 7 thresholds ν ∈ {−20, −6, −2, 0, 20, 6, 20} × 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, κ (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 ∈ {1, 2, 3}.
As mentioned prior, the CDFs are intimately connected to the moments.In particular, for a given threshold ν, the measurement in Equation 3.3 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. 6anerjee & 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 NL types.However, we will show in Section 4.3 that each f 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-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 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 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 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 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 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 ζ, 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 ϕ.A key choice in Sunseri et al. (2023) is splitting the radial and angular dependence of ζ as (3.4) Note that this expansion pertains to any projected, 2D scalar fields, which in our case is the convergence field, κ(⃗ x).Thus, ζ is a projected 3-point function, and r i are projected separations.The radial coefficients, ζ m (r 1 , r 2 ), can be computed as an area average over the convergence field κ(⃗ x), where the Fourier coefficients c m are c m (r i ; ⃗ x) ≡ dϕ e −imϕ κ(⃗ x + ⃗ 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 ζ computed in different radial and angular bins.We can thereby bin the Fourier coefficients in circular annuli as, where is the area corresponding to the annular bin, and r b min and r b max are the minimum and maximum radius of bin b.With this binning, we write the ζ m coefficients as where b 1 and b 2 are bin indices.Equation 3.4 shows that the 3-point function ζ can be constructed by summing over the coefficients ζ m with some angular phase, (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, ζ, 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 κ(⃗ x), c b 1 m (⃗ x) and c b 2 * m (⃗ x).If all fields come from the same bin then ζ 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 flatfield approximation is adequate -and compute c b m (⃗ x) separately for each patch.For each HEALPix pixel, we compute the quantity In our implementation, the patches follow a resolution of NSIDE = 4 which corresponds to a size 15 × 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 ζ, 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 max = 5, which is the same choice made in Sunseri et al. (2023).We do not explore higher m max given computing limitations.We compute r i in 6 logarithmic bins between 3.2 ′ < r i < 200 ′ , and ϕ in 6 linear bins between 0 < ϕ < π.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 2 bin .The choice of 6 ϕ bins is because ζ is computed using 6 m-modes.The quantities ζ m and ζ 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 ϕ 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 ζ m are complex coefficients, requiring a number for each of the real and complex parts, while ζ 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 ζ 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.

f 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 where d Xm dθ i is the mean derivative of point m in data-vector X with respect to parameter θ i , where the mean is computed using 300 to 400 independent survey realizations depending on the survey.C −1 is the inverse of the numerically estimated covariance matrix and is computed while accounting for the Hartlap factor (Hartlap et al. 2007), We verify in Section B that the covariance is well-converged for all summary statistics.The Hartlap factor is ≳ 0.95 for the 2nd and 3rd moments, the fiducial data-vector used in this work, and is ≳ 0.55 for the full 3-point function.Note that for the analysis using the 3point 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.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 NL , and the rows progressively step from constraints varying just f NL , to marginalizing over σ 8 and Ω m , and finally to also marginalizing over the intrinsic alignment parameters, η IA and A IA .The constraints for f loc NL and f eq NL , when varying just f 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 or, cmb NL .The 2nd moment constraints are above the range of many panels.

Dependence on summary statistic
Figure 2 shows the Fisher information in different summary statistics for an LSST Y10 survey.The first row shows constraints when varying just f NL (with y-labels denoting the specific type being varied).In the second row we marginalize over σ 8 and Ω 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 eq NL and f loc NL , 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 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 NL , the constraints from the 2nd moments alone are significantly weaker  The purple bands in most panels overlap completely with the yellow bands and are not seen.The 3point 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.
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 or, lss NL (f or, cmb
The relevance of the higher-order information improves significantly when extending the parameter space to perform a marginalized analysis of f 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 eq NL , worse for f loc NL , and comparable for f or, lss NL and f eq NL .This is generally consistent with the behaviors of the moments and CDFs found in the wCDM analysis of Anbajagane et al. (2023a).Further marginalizing over IA parameters results in minimal degradation of the f 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 ′ < θ < 200 ′ , as opposed to the 10 bins used in the main analysis. 7To 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 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 NL , and increases to 50% improvements when marginalizing over cosmology and/or IA.The exception is f or, cmb NL 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 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  practical challenges in robust use of the three-point function while performing simulationbased 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.2022), and for f loc NL we take the numbers from DESI Collaboration et al. (2016).Lensing measurements at θ min = 3.2 ′ (θ min = 20 ′ ) correspond to power spectra at k ∈ [1, 2.5] h/Mpc (k ∈ [0.3, 0.6] h/Mpc), depending on the redshift bin.These are higher than the current k max = 0.2h/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.

Fisher information in DES and LSST
The key goal of this work is extracting the Fisher information on f 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 eq NL and f or, lss NL , 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 loc NL .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 NL explored: For f loc NL , galaxy correlations have significantly more information than weak lensing.This is expected given the signature of f loc NL in galaxy correlation functions is a scaledependent galaxy bias in the 2-point function, where the bias increases towards large scales as b ∝ 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 loc NL 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 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 loc NL constraints.The analysis in this work does not compute the Fisher information in the galaxy correlations and is unable to test this.
For f eq NL , 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 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 f or, lss NL and f or, cmb NL , the constraints from weak lensing are about 1.5 to 4 times broader than those from galaxy clustering.We only compare f or, lss NL as the BOSS/DESI analyses do not measure f or, cmb NL .The DES constraints for the former f 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 NL amplitudes for Stage IV surveys.
The sound speed, c s , and the EFT parameter c 3 , capture the two leadingderivative 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 eq NL and f or, lss NL as where the expressions are taken from the conversions presented in Equation 5of Cabass et al. (2022b).Since the constraints allow the value c s = 1 (and thus 1/(c   Cabass et al. (2022b) in marginalizing over 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 eq NL and f or, lss NL , and not c s and 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 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 ≥ 0.013, and the CMB analysis of Planck Collaboration et al. (2020a), which finds c s ≥ 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 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 NL given their different redshift and scale-dependence (see Section 4.3).
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 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 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 eq NL , 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 gal = 10 (close to the DES Y6 fiducial value) and LSST Y10 n gal = 27 (the LSST Y10 fiducial value) is 20-30%.
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 eq NL , the constraints from a DES Y6 survey with n gal = 24 are 40% better than those of a fiducial survey with n gal = 9.At fixed source galaxy number density, the LSST Y10 constraints are ≈ 45% better than those of DES Y6.We also compute constraints for an LSST Y10-like survey but with n 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 eq NL and f or, lss NL .Using simple scaling arguments for the dependence of measurement noise on galaxy number density and sky area, the Fisher information is roughly proportional to √ n gal and f sky , where f sky is the fraction of the full sky covered by the survey.The exact scaling with n 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 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.

Physical origin of f 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 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 loc NL 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 θ > 10 ′ or θ > 20 ′ ; such cuts are employed in the moments-based weak lensing analyses of Gatti et al. (2020Gatti et al. ( , 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 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,8 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 max .Any cross-correlations between bins below z max with those above z 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 NL but the increased amplitude of weak lensing signals for highredshift observations -result in the sensitivity of f NL from weak lensing peaking at a redshift of z max ≈ 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 max ≈ 1.25 to explicitly confirm this.Note, however, that this discussion above concerns signatures (and thus constraints) of only f NL .We have verified that in analyses that also marginalize over Ω m , σ 8 , and the IA parameters, the higher redshift bins still add information (though, the choice z max ≈ 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 NL .
Halo mass function.We can then also directly explore the variation in the halo mass function (HMF) as we change f NL , as this will then imprint on the weak lensing field.Figure 7 shows the fractional change in the HMF for ∆f NL = 200, for all four types of f 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 M ⊙ .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 NL .Figure 7 presents a clear impact of f NL on halo counts towards the massive halo end.The sign of the fractional change is positive for f eq NL and f loc NL , meaning the abundance of massive halos increases with f NL , and this is expected as for positive f NL the initial conditions have a larger skewness and have more high-density peaks.At z = 0, the HMF for lower mass halos, M ≈ 10 14 M ⊙ , is reduced with increasing f loc NL and f eq NL .This sign flip is seen more prominently in both orthogonal-type f NL , though the relation is reversed as high f NL implies a lower halo count.For f or, cmb NL , the change is factors of 2-3 smaller than for all the other types.For both orthogonal-type f NL we once again find that the sign of the change is flipped at redshift z = 0, for halos with M ≈ 10 14 M ⊙ , as their abundance increases with f 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 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 or, lss NL reduces the amplitude of the power spectra on small scales, while increasing f or, cmb NL 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 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 mean > z 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 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 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.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 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 NL on the halo field is the scale-dependent bias found in the local-type f NL (Dalal et al. 2008).This bias goes as b ∝ f loc NL /k 2 and is a diverging signal for k → 0. Other types of f 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, where C ℓ is the power at different multipole ℓ.The choice of using C hm ℓ over C hh ℓ removes the impact of shot noise in our estimate of b ℓ .redshifts and for the four types of f 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 loc NL , where the bias of the halo sample grows at large scales (low ℓ).For f loc NL 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 or, cmb NL , which is expected to have a scaling of b ∝ 1/k on linear scales (Schmidt & Kamionkowski 2010). 9Finally, the equilateral type shows the halo bias effect is nearly scale independent for ℓ < 200; meaning, f eq NL simply alters the linear bias.For ℓ > 200, there is a mild scale-dependent, redshift-dependent behavior, though the amplitude of this variation is still low.The bias for f or, lss NL shows similar features to that of f eq NL , with a scale-independent bias below ℓ < 200, and a mildly dependent one above it.The behaviors of the bias on large scales, for all f NL types, are consistent with the perturbation theory predictions of Schmidt & Kamionkowski (2010).
All scale-dependent signatures of the different f NL are completely unused if we choose the halo counts as our statistic.On the other hand, the weak lensing field (or any large-9 The f or, cmb NL template is an approximation of the full EFT template (Senatore et al. 2010).The b ∝ 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 or, cmb NL 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 or, cmb NL template and the true, EFT template are inconsequential.The sample has a fixed lower mass cut of M > 10 14 M ⊙ at every redshift, and this is set by the resolution limit of the simulation.The bias depends on f 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 NL = ±100 use the same random seed for the initial conditions.
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 NL even on scales of θ > 10 ′ , where this angular scale corresponds to 5 to 20 Mpc across different redshifts, as such scales probe the "two-halo" term and the signatures of f NL imprinted into this term.
Figure 8 also shows that the impact of f 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 M ⊙ at z = 2 is a significantly rarer structure than a halo of M > 10 14 M ⊙ 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 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.

Discussion
Having shown that weak lensing can provide usable constraints on f 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.

Advantages of weak lensing as a probe of f 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 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. is therefore justified to assume weak lensing provides substantial/comparable information on f 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 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 sim ∼ 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.10These 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 NL .However, following the usage of weak lensing in large surveys, we can infer that of equal importanceif 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 NL from spectroscopic galaxy surveys utilize a one-loop bispectrum model, which has 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, k > 1Mpc −1 .Secco et al. (2022a, see their Figure 4) show that a lensing measurement at angular scale θ probes a range of wavenumbers k, represented heuristically as where w(k) is the weight of each mode, quantifying the sensitivity of measurement X(θ) to this mode, P (k) is the convergence power spectrum, and ξ κ (θ) 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 θ = 3.2 ′ correspond to mean wavenumbers of 1 < k [h/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 max(w(k)) -we find 3.5 < k [h/Mpc] < 6, depending on the tomographic bin.Even under our conservative angular scale cut of θ min = 20 ′ , we find the maximum contributing wavenumber is 0.  2022) use up to k ≲ 0.2 h/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 θ min = 20 ′ , 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).

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 γ 1,2 .However, this is an approximation as the galaxy shapes actually trace, e obs 1,2 = γ 1,2 /(1 − κ), where κ is the convergence field.In the limit of κ ≪ 1, the expression can be Taylor-expanded to e obs 1,2 = γ 1,2 (1 + κ + κ 2 /2 + . ..).The assumption of e obs 1,2 ≈ γ 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 ∝ (1 + qκ), where q = 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(Gatti et al. , 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 Ω m and σ 8 ) leads to LSST Y10-based constraints that are still comparable to DESI for multiple types of f 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 Nbody 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,b).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;Anbajagane et al. 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/Mpc (Schneider et al. 2019;Giri & Schneider 2021) and for some higher-order statistics down to θ = 1 ′ 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.

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 NL parameter.We consider four types of f NL -f loc NL which arises from multi-field inflation, f eq NL from a strong self-coupling of the inflaton field, f or, lss NL and f or, cmb NL from the same physics as f eq NL 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 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 loc NL and f eq NL .However, the latter are clearly better for f or, lss NL and f or, cmb NL , and for analyses of all f NL that marginalize over cosmology (Ω m and σ 8 ) and intrinsic alignment parameters (Figure 2).
• The shape/configuration information in the lensing field adds significantly to f 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 eq NL , LSST Y10 (DES Y6) is competitive/better than DESI (BOSS).For f or, lss NL , 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 θ > 20 ′ (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 ≈ 1.25.Including source galaxies beyond this redshift helps modestly with the constraining power when varying only f NL .However, the high redshift data is valuable in marginalized f NL constraints, as it helps in constraining the cosmological and IA parameters that are marginalized (Figure 6).
• The impact of f NL on the lensing field, including the range of scales that probe f NL , can be understood through the impact of f 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 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 eq NL (and possibly f or, lss NL ) -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 eq NL and f or, lss NL 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 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 twopoint 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.There is clear agreement between our simulations and the quijote runs.The specific ordering of the 10 most massive halos varies slightly.The distribution of 1 + δ is within 1% for δ > 1, and grows to 10% for δ ≈ 0. The power spectra agree to 1% at z = 0 and grow to 2 − 3% at higher redshifts, which is within the expected deviations due to differences in N-body solvers (Schneider et al. 2016) range of δ values, and the difference grows to 10% only for δ ≈ 0, near the tails of the distribution.The bottom right panel shows the comparison of the power spectra for the two runs at different redshifts.We once again find that the deviations are minimal; they are within 1% at z = 0 and grow to 2 − 3% at higher redshift.Such differences, at our resolution level of 512 3 particles, are consistent with expectations from comparison studies of the Pkdgrav3 and Gadget3 (Schneider et al. 2016).
Power Spectra.We then validate the simulations against a number of theoretical models, starting with the power spectrum of the density field.Figure 10 shows the fractional deviation between the Ulagam suite (from the average of five realizations) alongside four different theoretical predictions.Two are simply the average of five lower/higher resolutions runs from the Ulagam suite, while the other two are predictions from Halofit (Takahashi et al. 2012) andEuclid-Emu (Euclid Collaboration et al. 2019).Both models take the perturbation theory result from the Class boltzmann code (Lesgourgues 2011) and modify it to model the non-linear regime.The agreement with all models (other than the lower resolution run) is within 5% up to k = 0.3 Mpc −1 at z = 2, and up to k = 1 Mpc −1 at z = 0.While the deviations with Euclid-Emu increase towards high k, improving the simulation resolution to 1024 3 particles results in significantly better agreement.The Euclid-Emu model was built using Pkdgrav3 simulations, but run in a larger volume of L = 1.9 Gpc and with 2048 NL , Ω m , σ 8 , A IA and η IA .The dashed contours show constraints over all five paramters, while solid contours are constraints varying f NL and cosmology parameters alone.There is only marginal differences between the two.The 1D posterior for f eq NL also shows the unmarginalized constraint in the dotted lines.
.13) where κ 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 fromMyles et al. (2021), for DES Year 6, LSST Year 1 and Year 10 we use the same n(z) utilized inZhang et al. (2022, see their as defined within the integral of Equation 3.8, and average it across the survey footprint to obtain ζ b 1 ,b 2 m .

Figure 3 .
Figure 3. Similar to Figure2but 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 3point 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 Figure2as 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 4 .
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 eqNL are comparable to the existing BOSS constraints.The LSST Y10 constraints are comparable or potentially better than DESI for f eq NL , a factor of 2 broader for f or, lss NL , and a factor of 8 broader for f loc NL .The DES Y3 constraints are above the range of the plots for most panels.

Figure 5 .
Figure5.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 gal ).The fiducial values are n 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 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 Figure6).For f eq NL , 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 gal = 10 (close to the DES Y6 fiducial value) and LSST Y10 n gal = 27 (the LSST Y10 fiducial value) is 20-30%.

Figure 8 Figure 7 .
Figure 8 presents the fractional change in the halo bias as a function of ℓ, for four

Figure 8 .
Figure 8.The fractional change in halo bias due to f NL , shown as a function of multipole, ℓ, for the halo samples of different redshift bins.The variation is computed as ln b h(f NL = 100) − ln b h (f NL = −100).The sample has a fixed lower mass cut of M > 10 14 M ⊙ at every redshift, and this is set by the resolution limit of the simulation.The bias depends on f 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 NL = ±100 use the same random seed for the initial conditions.

Figure 9 .
Figure 9.The density fields of a Ulagam and Quijote simulation with the same initial condition.The slices have 500 Mpc/h thickness.The density distribution is shown in the top right and the power spectra in the bottom right.White circles show the location of the 10 most massive halos in each slice.There is clear agreement between our simulations and the quijote runs.The specific ordering of the 10 most massive halos varies slightly.The distribution of 1 + δ is within 1% for δ > 1, and grows to 10% for δ ≈ 0. The power spectra agree to 1% at z = 0 and grow to 2 − 3% at higher redshifts, which is within the expected deviations due to differences in N-body solvers(Schneider et al. 2016)

Figure 15 .
Figure15.The full Fisher constraints for f eq NL , Ω m , σ 8 , A IA and η IA .The dashed contours show constraints over all five paramters, while solid contours are constraints varying f NL and cosmology parameters alone.There is only marginal differences between the two.The 1D posterior for f eq NL also shows the unmarginalized constraint in the dotted lines.
The halo catalogs contain all identified

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 ±∆P values.The cosmological parameter values used in the fiducial runs are all the bolded values.We always assume a flat cosmology with Ω Λ = 1 − Ω m .The parameters not shown above are h = 0.6711 and Ω b = 0.049.
Myles et al. (2021)t distribution of the signal in each tomographic bin for each survey.The DES Y3 distributions are fromMyles 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 ) = χ j (χ s − χ j )/(a(z j )χ s ) as used in Equation2.1,withthedifferent 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 Figure6below shows our results are insensitive to the presence of this tail.

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

Table 2 .
Zhang et al. (2022,bution and source galaxy number density assumed for the upcoming surveys.All numbers are taken fromZhang et al. (2022,see their Table 1).

Table 3 .
The Fisher information presented in Figure4for the full range of scales (top) and a conservative scale cut (middle).The numbers in square brackets are constraints after marginalizing over Ω Philcox et al. (2022a)arameters.The "BOSS (fix)" results, analyzed at fixed cosmology, are from D'Amico et al. (2022) while the "BOSS (vary)" results from also varying σ 8 are fromPhilcox et al. (2022a).The expected constraints from DESI for f eq NL and f or, lss NL are obtained by rescaling the corresponding BOSS results by 1.6 as mentioned in D'Amico et al. (

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 Ω m , σ 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.can only constrain the combination (c −2 s − 1) c 3 rather than c 3 .Equations 4.3 and 4.4 require joint constraints on f eq NL and f or, lss NL , whereas we have thus far only varied one f 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,